polyhedron.h

another source code for reduced basis

  ITail[i] = new IDList();
  IHead[i]->next = ITail[i];
  ITail[i]->back = IHead[i];
  PHead[i] = new PolarList();
  PTail[i] = new PolarList();
  PHead[i]->next = PTail[i];
  PTail[i]->back = PHead[i];
  FHead[i] = new IDList();
  FTail[i] = new IDList();
  FHead[i]->next = FTail[i];
  FTail[i]->back = FHead[i];
  
  VHead[i] = new IDList();
  VTail[i] = new IDList();
  VHead[i]->next = VTail[i];
  VTail[i]->back = VHead[i];
  
}
  void Polyhedron::SetBoundaryLines(){
    int i=0;numboundary=0;
    if(neighborI!=NULL && neighborF!=NULL && boundary!=NULL)
      for(i=0;i<numberV;i++){
	if(((neighborI[i]) == neighborF[i]) && 
	   (neighborF[i] !=0) && 
	   (neighborI[i] !=0)){
	  boundary[i] = 0;
	}else{
	  boundary[i] = 1;
	  numboundary++;
	}
      }
    //printf("number of boundary points = %d\n",numboundary);
  }
  void Polyhedron::setFace(int i,int di,int dj,int dk){
  Face[i] = new int[3];
  Face[i][0] = di;
  Face[i][1] = dj;
  Face[i][2] = dk;
  /* One */
  neighborF[di]++;
  
  IDtool->AppendISort(dj,IHead[di],ITail[di],di,neighborI);
  IDtool->AppendISort(dk,IHead[di],ITail[di],di,neighborI);
  IDtool->AppendVF(dj,VTail[di]);
  IDtool->AppendVF(dk,VTail[di]);
  /* Two */
  neighborF[dj]++;
  
  IDtool->AppendISort(di,IHead[dj],ITail[dj],dj,neighborI);
  IDtool->AppendISort(dk,IHead[dj],ITail[dj],dj,neighborI);
  IDtool->AppendVF(dk,VTail[dj]);
  IDtool->AppendVF(di,VTail[dj]);
  /* Three */
  neighborF[dk]++;
  
  IDtool->AppendISort(di,IHead[dk],ITail[dk],dk,neighborI);
  IDtool->AppendISort(dj,IHead[dk],ITail[dk],dk,neighborI);
  IDtool->AppendVF(di,VTail[dk]);
  IDtool->AppendVF(dj,VTail[dk]);
  }


  void Polyhedron::setPolarMapNaturalB(){

  int i,j,k;
  IDList *now=NULL;
  IDList *now2=NULL;
  double dx,dy;
  PolarList *nowp=NULL;
  PolarList *nowp2=NULL;
  int nowID,nextID;
  double angle = 0.0;
  double theta = 0.0;
  int bid0=0;
  int bid1=0;
   int sw0=0;
  int cn0=0;
  // make Geodesic Polar Map 
  for(i=0;i<numberV;i++){
    dy = 0.0;
    theta = 0.0;
    now = VHead[i];
    while(next(now)!=VTail[i]){
      now = next(now);
      PT->makeVector(bc[0],point[i],point[now->ID]);
      PT->makeVector(bc[1],point[i],point[next(now)->ID]);
      angle = acos((PT->InnerProduct(bc[0],bc[1])/(PT->Point3dSize(bc[0])*PT->Point3dSize(bc[1]))));
      theta += angle;
      now = next(now);
    }
    
    if(boundary[i]!=1){
      now = IHead[i];
      now = next(now);
      nowID = now->ID;
      dx = PT->Distance(point[i],point[now->ID]);
      dy = 0.0;
      IDtool->AppendPolarI(nowID,PTail[i],dx,dy);
      cn0 = 1;
      while(cn0<neighborI[i]){
	now = VHead[i];
	while(next(now)!=VTail[i]){
	  now = next(now);
	  if(now->ID==nowID){
	    nextID = next(now)->ID;
	    break;
	  }
	  now = next(now);
	}
	
	PT->makeVector(bc[0],point[i],point[nowID]);
	PT->makeVector(bc[1],point[i],point[nextID]);
	dx = PT->Point3dSize(bc[1]);
	angle = acos((PT->InnerProduct(bc[0],bc[1])/(PT->Point3dSize(bc[0])*PT->Point3dSize(bc[1]))));
	dy += (2.0*PI*angle)/theta;
	IDtool->AppendPolarI(nextID,PTail[i],dx,dy);
	nowID = nextID;
	cn0++;
      }
    }else{
       int thcheck=0;
      now = IHead[i];
      while(next(now)!=ITail[i]){
	now = next(now);
	if(neighborI[now->ID]==2&&boundary[now->ID]==1){
	  thcheck=1;
	  nowID = now->ID;break;
	}
	
      }
      if(thcheck==1){sw0=0;
	now = VHead[i];
	while(next(now)!=VTail[i]){
	  now = next(now);
	  if(now->ID==nowID){
	    sw0=1;
	    break;
	  }
	  now = next(now);
	}
	if(sw0==0){
	  cn0 = 1;
	
	  while(cn0<neighborI[i]){
	    now = VHead[i];
	    while(next(now)!=VTail[i]){
	      now = next(now);
	      if(next(now)->ID==nowID){
		nextID = now->ID;
	      
		break;
	      }
	      now = next(now);
	    }
	    nowID = nextID;
	    cn0++;
	  }
	}
      }else{

      
	now = VHead[i];
	while(next(now)!=VTail[i]){
	  now = next(now);
	  if(boundary[now->ID]==1){
	    nowID = now->ID;
	    break;
	  }
	  now = next(now);
	}
      }
      
	dx = PT->Distance(point[i],point[nowID]);
	dy = 0.0;
	IDtool->AppendPolarI(nowID,PTail[i],dx,dy);
	cn0 = 1;
	
	while(cn0<neighborI[i]){
	  now = VHead[i];
	  while(next(now)!=VTail[i]){
	    now = next(now);
	    if(now->ID==nowID){
	      nextID = next(now)->ID;
	      
	      break;
	    }
	    now = next(now);
	  }
	  
	  PT->makeVector(bc[0],point[i],point[nowID]);
	  PT->makeVector(bc[1],point[i],point[nextID]);
	  dx = PT->Point3dSize(bc[1]);
	  angle = acos((PT->InnerProduct(bc[0],bc[1])/(PT->Point3dSize(bc[0])*PT->Point3dSize(bc[1]))));
	  
	  dy += (PI*angle)/theta;
	  IDtool->AppendPolarI(nextID,PTail[i],dx,dy);
	  
	  nowID = nextID;
	  cn0++;
	}
           
      
      
      
      
      
    }
  }
}

  void Polyhedron::setPolarMap(){
    
    int i,j,k;
  IDList *now=NULL;
  IDList *now2=NULL;
  double dx,dy;
  PolarList *nowp=NULL;
  PolarList *nowp2=NULL;
  int nowID,nextID;
  double angle = 0.0;
  double theta = 0.0;
  int cn0=0;
  // make Geodesic Polar Map 
  for(i=0;i<numberV;i++){
    if(boundary[i]!=1){
      
      dy = 0.0;
      theta = 0.0;
      now = VHead[i];
      while(next(now)!=VTail[i]){
	now = next(now);
	PT->makeVector(bc[0],point[i],point[now->ID]);
	PT->makeVector(bc[1],point[i],point[next(now)->ID]);
	angle = acos((PT->InnerProduct(bc[0],bc[1])/(PT->Point3dSize(bc[0])*PT->Point3dSize(bc[1]))));
	theta += angle;
	now = next(now);
      }
      
      
      now = IHead[i];
      now = next(now);
      nowID = now->ID;
      dx = PT->Distance(point[i],point[now->ID]);
      dy = 0.0;
      IDtool->AppendPolarI(nowID,PTail[i],dx,dy);
      cn0 = 1;
      while(cn0<neighborI[i]){
	now = VHead[i];
	while(next(now)!=VTail[i]){
	  now = next(now);
	  if(now->ID==nowID){
	    nextID = next(now)->ID;
	    break;
	  }
	  now = next(now);
	}
	
	PT->makeVector(bc[0],point[i],point[nowID]);
	PT->makeVector(bc[1],point[i],point[nextID]);
	dx = PT->Point3dSize(bc[1]);
	angle = acos((PT->InnerProduct(bc[0],bc[1])/(PT->Point3dSize(bc[0])*PT->Point3dSize(bc[1]))));
	dy += (2.0*PI*angle)/theta;
	IDtool->AppendPolarI(nextID,PTail[i],dx,dy);
	nowID = nextID;
	cn0++;
      }
    }
  }
  
  
  
  
}

void Polyhedron::setAreaMap3D(){
  int i,j;
  IDList *now=NULL;
  sumarea3D = 0.0;
  for(i=0;i<numberF;i++){
    PT->makeVector(bc[0],point[Face[i][0]],point[Face[i][1]]);
    PT->makeVector(bc[1],point[Face[i][0]],point[Face[i][2]]);
    PT->CrossVector(bc[2],bc[1],bc[0]);
    areaMap3D[i] = PT->Point3dSize(bc[2])/2.0;
    sumarea3D += areaMap3D[i];
  }
  
  if(boundarytype==0||boundarytype==2){
    constsumarea3D = sqrt(1.0/sumarea3D);
  }else if(boundarytype==1){
    constsumarea3D = sqrt(((0.5*0.5*PI)/sumarea3D));
  }
 
}




// U1 low stretch parameterization
 void Polyhedron::ParametrizationSingle(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;
   
   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,(1000+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 

  
  
  // local stretch calculation
  setSigmaZero();
  
  // Re-solving linear system according to local stretches

  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){
      pU[i] = UaXY[i+1];
      pV[i] = UaXY[i+numberV+1];
    }
  }
    
  
  printf("STL2 error = %lf\n",getCurrentE());
  
  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];
    }
  }

  delete mybcg;
  delete [] vecb;
  delete [] UaXY;
  delete [] sigsum;
  
  
  
 }
  // UOpt: Optimal low stretch parameterization
  void Polyhedron::ParametrizationOptimal(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];
  }