polyhedron.h~

another source code for reduced basis

	 mybcg->ija[dlk]=nowp->ID+numberV+1;
       }
    }
    mybcg->ija[i+numberV+1+1]=dlk+1;
  }
  
  for(i=0;i<numberV;i++){
    UaXY[i+1] = pU[i];
    UaXY[i+numberV+1] = pV[i];
  }
  mybcg->linbcg(((unsigned long)(2*(numberV))),vecb,UaXY,1,error,itenum,&iter,&linerr);
  
  for(i=0;i<numberV;i++){
    if(boundary[i]!=1){
      pU[i] = UaXY[i+1];
      pV[i] = UaXY[i+numberV+1];
    }
  }
  
  // Re-solving linear system
  double initialstrech=0.0;
  double currentstrech=0.0;
  Point2d **prevU = new Point2d* [numberV];
  for(i=0;i<numberV;i++){
    prevU[i] = new Point2d(0.0,0.0); 
  }
  initialstrech = getCurrentE();
  int kk = 0;
  
  for(kk=0;kk<itenum;kk++){
    
    setSigmaZero();
    for(i=0;i<numberV;i++){
      
      if(boundary[i]!=1){
	sigsum[i]=0.0;
	nowp = PHead[i];
	while(next(nowp)!=PTail[i]){
	  nowp = next(nowp);
	 
	  nowp->lambda /= sigma[nowp->ID];
	  
	 
	  
	  sigsum[i] += nowp->lambda;
	}
	
	mybcg->sa[i+1] = 1.0;
	mybcg->sa[i+1+numberV] = 1.0;
	
      }else{
      
	mybcg->sa[i+1] = 1.0;
	mybcg->sa[i+1+numberV] = 1.0;
      }
    }
    dlk=2*(numberV)+1;
    
    for(i=0;i<numberV;i++){
      if(boundary[i]!=1){
	nowp = PHead[i];
	
	while(next(nowp)!=PTail[i]){
	  nowp = next(nowp);
	  ++dlk;
	  mybcg->sa[dlk] = -nowp->lambda/sigsum[i];
	}
      }
    }
    for(i=0;i<numberV;i++){
      if(boundary[i]!=1){
	nowp = PHead[i];
	while(next(nowp)!=PTail[i]){
	  nowp = next(nowp);
	  ++dlk;
	  mybcg->sa[dlk] = -nowp->lambda/sigsum[i];
	}
      }
    }
     
    for(i=0;i<numberV;i++){
      UaXY[i+1] = pU[i];
      UaXY[i+numberV+1] = pV[i];
    }
    mybcg->linbcg(((unsigned long)(2*(numberV))),vecb,UaXY,1,error,itenum,&iter,&linerr); 
    for(i=0;i<numberV;i++){
      if(boundary[i]!=1){
	prevU[i]->x = pU[i];
	prevU[i]->y = pV[i];
	
	pU[i] = UaXY[i+1];
	pV[i] = UaXY[i+numberV+1];
      }
    }
    currentstrech = getCurrentE();
    //printf("currentstrech = %lf\n",currentstrech);
    if(initialstrech<currentstrech){
      for(i=0;i<numberV;i++){
	if(boundary[i]!=1){
	  pU[i] = prevU[i]->x;
	  pV[i] = prevU[i]->y;
	}
      }
      break;
    }else{
      initialstrech = currentstrech;
    }
  }
  
  level = kk;
 
 
  printf("STL2 error = %lf\n",getCurrentE());
  delete mybcg;
  for(i=0;i<numberV;i++){
    
    if(boundary[i]!=1){
    IDtool->CleanNeighborPolar(PHead[i],PTail[i]);
      
    PHead[i] = new PolarList();
    PTail[i] = new PolarList();
    PHead[i]->next = PTail[i];
    PTail[i]->back = PHead[i];
    }
  }
  for(i=0;i<numberV;i++){
    delete prevU[i];
  }

 

  delete [] prevU;
  delete [] vecb;
  delete [] UaXY;
  delete [] sigsum;

}
 
 // UOpt: Optimal low stretch parameterization (slow and smooth)
  void Polyhedron::ParametrizationSmoothOptimal(int itenum,double error){
    double *UaXY = new double[2*(numberV)+1];
    double *vecb = new double[2*(numberV)+1];
  
  int i;
  IDList *now;
  IDList *now2;
  PolarList *nowp;
  level=0;
  
  int nonzero=(numberV);
  for(i=0;i<numberV;i++){
    vecb[i+1]=0.0;
    if(boundary[i]!=1){
      now = IHead[i];
      while(next(now)!=ITail[i]){
	now = next(now);
      nonzero++;
      }
    }
  }
  int iter=0;
  double linerr=0.0;
  double weight=0.0;
  
  PCBCGSolver *mybcg = new PCBCGSolver(2*nonzero);
  double *sigsum = new double[numberV];
  if(weighttype==0){
    setFloaterC();
  }else if(weighttype==1){
    setLaplaceC();
  }else if(weighttype==2){
    setEckHC();
  }else if(weighttype==3){
    setDesbrunC();
  }else if(weighttype==4){
    setMVCC();
  }else{
    setFloaterC();
  }
    
  SortIndexP();
  
  
  for(i=0;i<numberV;i++){
    
    if(boundary[i]!=1){
      
      
      
      mybcg->sa[i+1] = 1.0;
      mybcg->sa[i+1+numberV] = 1.0;
      vecb[i+1] = 0.0;
      vecb[i+1+numberV] = 0.0;
    }else{
      mybcg->sa[i+1] = 1.0;
      vecb[i+1] = pU[i];
      mybcg->sa[i+1+numberV] = 1.0;
      vecb[i+1+numberV] = pV[i];
    }
  }
  mybcg->ija[1] = 2*(numberV)+2;
  int dlk=2*(numberV)+1;
  
  for(i=0;i<numberV;i++){
    if(boundary[i]!=1){
      nowp = PHead[i];
      
      while(next(nowp)!=PTail[i]){
	nowp = next(nowp);
	
	
	++dlk;
	mybcg->sa[dlk] = -nowp->lambda;
	mybcg->ija[dlk]=nowp->ID+1;
      }
    }
    mybcg->ija[i+1+1]=dlk+1;
  }
  for(i=0;i<numberV;i++){
    if(boundary[i]!=1){
       nowp = PHead[i];
       while(next(nowp)!=PTail[i]){
	 nowp = next(nowp);
	 ++dlk;
	 mybcg->sa[dlk] = -nowp->lambda;
	 mybcg->ija[dlk]=nowp->ID+numberV+1;
       }
    }
    mybcg->ija[i+numberV+1+1]=dlk+1;
  }
  
  for(i=0;i<numberV;i++){
    UaXY[i+1] = pU[i];
    UaXY[i+numberV+1] = pV[i];
  }
  mybcg->linbcg(((unsigned long)(2*(numberV))),vecb,UaXY,1,error,itenum,&iter,&linerr);
  
  for(i=0;i<numberV;i++){
    if(boundary[i]!=1){
      pU[i] = UaXY[i+1];
      pV[i] = UaXY[i+numberV+1];
    }
  }
  
  // Re-solving linear system
  double initialstrech=0.0;
  double currentstrech=0.0;
  Point2d **prevU = new Point2d* [numberV];
  for(i=0;i<numberV;i++){
    prevU[i] = new Point2d(0.0,0.0); 
  }
  initialstrech = getCurrentE();
  int kk = 0;
  
  for(kk=0;kk<itenum;kk++){
    setSigmaZero();
    for(i=0;i<numberV;i++){
      
      if(boundary[i]!=1){
	sigsum[i]=0.0;
	nowp = PHead[i];
	while(next(nowp)!=PTail[i]){
	  nowp = next(nowp);
	  
	  nowp->lambda /= (0.5*(sigma[nowp->ID]+sigma[i]));
	  
	  
	  
	  sigsum[i] += nowp->lambda;
	}
	
	mybcg->sa[i+1] = 1.0;
	mybcg->sa[i+1+numberV] = 1.0;
	
      }else{
      
	mybcg->sa[i+1] = 1.0;
	mybcg->sa[i+1+numberV] = 1.0;
      }
    }
    dlk=2*(numberV)+1;
    
    for(i=0;i<numberV;i++){
      if(boundary[i]!=1){
	nowp = PHead[i];
	
	while(next(nowp)!=PTail[i]){
	  nowp = next(nowp);
	  ++dlk;
	  mybcg->sa[dlk] = -nowp->lambda/sigsum[i];
	}
      }
    }
    for(i=0;i<numberV;i++){
      if(boundary[i]!=1){
	nowp = PHead[i];
	while(next(nowp)!=PTail[i]){
	  nowp = next(nowp);
	  ++dlk;
	  mybcg->sa[dlk] = -nowp->lambda/sigsum[i];
	}
      }
    }
    
    for(i=0;i<numberV;i++){
      UaXY[i+1] = pU[i];
      UaXY[i+numberV+1] = pV[i];
    }
    mybcg->linbcg(((unsigned long)(2*(numberV))),vecb,UaXY,1,error,itenum,&iter,&linerr);
    for(i=0;i<numberV;i++){
      if(boundary[i]!=1){
	prevU[i]->x = pU[i];
	prevU[i]->y = pV[i];
	
	pU[i] = UaXY[i+1];
	pV[i] = UaXY[i+numberV+1];
      }
    }
    currentstrech = getCurrentE();
    //printf("currentstrech = %lf\n",currentstrech);
    if(initialstrech<currentstrech){
      for(i=0;i<numberV;i++){
	if(boundary[i]!=1){
	  pU[i] = prevU[i]->x;
	  pV[i] = prevU[i]->y;
	}
      }
      break;
    }else{
      initialstrech = currentstrech;
    }
  }
  
  level = kk;
  
  
  printf("STL2 error = %lf\n",getCurrentE());
  delete mybcg;
  for(i=0;i<numberV;i++){
    
    if(boundary[i]!=1){
    IDtool->CleanNeighborPolar(PHead[i],PTail[i]);
      
    PHead[i] = new PolarList();
    PTail[i] = new PolarList();
    PHead[i]->next = PTail[i];
    PTail[i]->back = PHead[i];
    }
  }
  for(i=0;i<numberV;i++){
    delete prevU[i];
  }

  

  delete [] prevU;
  delete [] vecb;
  delete [] UaXY;
  delete [] sigsum;

}
 
 
  // compute natural boundary 
  void Polyhedron::setNaturalB(int itenum,double error){
  double *UaXY = new double[2*(numberV)+1];
  double *vecb = new double[2*(numberV)+1];
  
  int i;
  IDList *now;
  IDList *now2;

  int myite = 2000;
  
  if(numberV > 20000)myite = 5000;
  if(numberV > 50000)myite = 10000;
  
  int Bfix0=0;
  int Bfix1=0;
  for(i=0;i<numberV;i++){
    if(boundary[i]==1){
      Bfix0 = i;
      break;
    }
  }
  
  for(i=0;i<numberV;i++){
    if(boundary[i]==1&&neighborI[i]==2){
      Bfix0 = i;
      break;
    }
  }
  

  double dmax=0.0;
  double ddis=0.0;
  int firstcheck=0;
  double dx,dy;
  for(i=0;i<numberV;i++){
    if(boundary[i]==1){
      dx = pU[i]-pU[Bfix0];
      dy = pV[i]-pV[Bfix0];
      
      if(firstcheck==0){
	firstcheck=1;
	dmax = sqrt((dx*dx+dy*dy));
	Bfix1 = i;
      }else{
	ddis = sqrt((dx*dx+dy*dy));
	if(ddis>dmax){
	  dmax = ddis;
	  Bfix1 = i;
	}
      }    
    }
  }

  
  int nonzero=(numberV);
  for(i=0;i<numberV;i++){
    vecb[i+1]=0.0;
    if(i!=Bfix0&&i!=Bfix1){
      now = IHead[i];
      while(next(now)!=ITail[i]){
	now = next(now);
	nonzero++;
      }
      if(boundary[i]==1)nonzero += 2;
    }
  }
  
  int bbID0=0;
  int bbID1=0;
  
  

  int iter=0;
  double linerr=0.0;
  double weight=0.0;
  double *coef = new double[4];
  Point2d *v0 = new Point2d(0.0,0.0);
  Point2d *v1 = new Point2d(0.0,0.0);
  Point2d *v2 = new Point2d(0.0,0.0);
  Point2d *v3 = new Point2d(0.0,0.0);
  PCBCGSolver *mybcg = new PCBCGSolver(2*nonzero);
  
  int nextID=0;
  int backID=0;

  double conformalterm;
  double *sigsum = new double[numberV];
  
  setPolarMapNaturalB();
  SortIndexPNaturalB();
  



  double ancos1=0.0;
  PolarList *nowp=NULL;
  double dl1,dl2;
  for(i=0;i<numberV;i++){
    nowp = PHead[i];
    while(next(nowp)!=PTail[i]){
      nowp = next(nowp);
      nowp->lambda = 0.0;
      now = VHead[i];
      while(next(now)!=VTail[i]){
	now = next(now);
	
	
	if(nowp->ID==now->ID){
	  PT->makeVector(bc[0],point[next(now)->ID],point[i]);
	  PT->makeVector(bc[1],point[next(now)->ID],point[now->ID]);
	  ancos1 = PT->InnerProduct(bc[0],bc[1])/(PT->Point3dSize(bc[0])*PT->Point3dSize(bc[1]));
	  nowp->lambda += (ancos1/sqrt(1.0-ancos1*ancos1));
	  
	}
	if(nowp->ID==next(now)->ID){
	  PT->makeVector(bc[0],point[now->ID],point[i]);
	  PT->makeVector(bc[1],point[now->ID],point[next(now)->ID]);
	  ancos1 = PT->InnerProduct(bc[0],bc[1])/(PT->Point3dSize(bc[0])*PT->Point3dSize(bc[1]));
	  nowp->lambda += (ancos1/sqrt(1.0-ancos1*ancos1));
	}
	now = next(now);
      }
    }
  }
  



  int cn=0;
  for(i=0;i<numberV;i++){
    if(i==Bfix0||i==Bfix1){
      mybcg->sa[i+1] = 1.0;
      vecb[i+1] = pU[i];
      mybcg->sa[i+1+numberV] = 1.0;
      vecb[i+1+numberV] = pV[i];
    }else{
      sigsum[i]=0.0;
      nowp = PHead[i];
      while(next(nowp)!=PTail[i]){
	nowp = next(nowp);
	sigsum[i] += nowp->lambda;	
      }
      mybcg->sa[i+1] = -1.0;
      mybcg->sa[i+1+numberV] = -1.0;
      vecb[i+1] = 0.0;
      vecb[i+1+numberV] = 0.0;
    }
  
  }
  mybcg->ija[1] = 2*(numberV)+2;
  int dlk=2*(numberV)+1;
  
  int BID0=0;
  int BID1=0;
  
  for(i=0;i<numberV;i++){
    if(i!=Bfix0&&i!=Bfix1){
      
      nowp = PHead[i];
      while(next(nowp)!=PTail[i]){
	nowp = next(nowp);
	++dlk;
	mybcg->sa[dlk] = nowp->lambda/sigsum[i];
	mybcg->ija[dlk]=nowp->ID+1;
      }
      if(boundary[i]==1){
	
	nowp = PHead[i];
	 while(next(nowp)!=PTail[i]){
	   nowp = next(nowp);
	   if(boundary[nowp->ID]==1&&nowp->theta==0.0){
	     BID0 = nowp->ID;
	     break;
	   }
	 }
	 nowp = PHead[i];
	 while(next(nowp)!=PTail[i]){
	   nowp = next(nowp);
	   if(boundary[nowp->ID]==1&&nowp->theta!=0.0){
	     
	     BID1 = nowp->ID;
	     break;
	   }
	 }
	
	 
	 if(BID0<BID1){
	   ++dlk;
	   mybcg->sa[dlk] = 1.0/sigsum[i];
	   mybcg->ija[dlk]=BID0+numberV+1;
	   
	   ++dlk;
	   mybcg->sa[dlk] = -1.0/sigsum[i];
	   mybcg->ija[dlk]=BID1+numberV+1;
	 }else{
	   ++dlk;
	   mybcg->sa[dlk] = -1.0/sigsum[i];
	   mybcg->ija[dlk]=BID1+numberV+1;
	   
	   ++dlk;
	   mybcg->sa[dlk] = 1.0/sigsum[i];
	   mybcg->ija[dlk]=BID0+numberV+1;
	 }
      }
      
      
      
    }    
    
    mybcg->ija[i+1+1]=dlk+1;
  }
  for(i=0;i<numberV;i++){
    if(i!=Bfix0&&i!=Bfix1){
      
    if(boundary[i]==1){