polyhedron.h

another source code for reduced basis

	   mybcg->sa[dlk] = -1.0/sigsum[i];
	   mybcg->ija[dlk]=BID0+1;
	   ++dlk;
	   mybcg->sa[dlk] = 1.0/sigsum[i];
	   mybcg->ija[dlk]=BID1+1;
	 }else{
	   ++dlk;
	   mybcg->sa[dlk] = 1.0/sigsum[i];
	   mybcg->ija[dlk]=BID1+1;
	   ++dlk;
	   mybcg->sa[dlk] = -1.0/sigsum[i];
	   mybcg->ija[dlk]=BID0+1;
	 }
      }     
    nowp = PHead[i];
    while(next(nowp)!=PTail[i]){
      nowp = next(nowp);
      ++dlk;
      mybcg->sa[dlk] = nowp->lambda/sigsum[i];
      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,(myite),&iter,&linerr);
  
  for(i=0;i<numberV;i++){
    if(i!=Bfix0&&i!=Bfix1){
      pU[i] = UaXY[i+1];
      pV[i] = UaXY[i+numberV+1];
    }
  }
  for(i=0;i<numberV;i++){
    
    
    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 [] coef;
  delete v0;
  delete v1;
  delete v2;
  delete v3;
  delete [] sigsum;

}
void Polyhedron::setEckHC(){
  PolarList *nowp;
  int i;
  int nextID=0;
  int backID=0;
  double ddweight=0.0;
  for(i=0;i<numberV;i++){ 
    if(boundary[i]!=1){
      nowp =PHead[i];
      ddweight=0.0;
      while(next(nowp)!=PTail[i]){
	nowp = next(nowp);
	if(nowp->next!=PTail[i]){
	  nextID = nowp->next->ID;
	}else{
	  nextID = PHead[i]->next->ID;
	}
	if(nowp->back!=PHead[i]){
	  backID = nowp->back->ID;
	}else{
	  backID = PTail[i]->back->ID;
	}
	nowp->lambda = getHarmonicTerm(i,backID,nowp->ID,nextID);
	ddweight += nowp->lambda;
      }
      
      if(ddweight!=0.0){
	nowp =PHead[i];
	while(next(nowp)!=PTail[i]){
	  nowp = next(nowp);
	  nowp->lambda /= ddweight;
	}
      }
    }
  }
}
double Polyhedron::getHarmonicTerm(int i,int backID,int nowID,int nextID){
  double ancos1,ancos2;
  PT->makeVector(bc[0],point[i],point[backID]);
  PT->makeVector(bc[1],point[i],point[nowID]);
  PT->makeVector(bc[2],point[i],point[nextID]);
  PT->CrossVector(bc[3],bc[0],bc[1]);

  double area1 = 0.5*PT->Point3dSize(bc[3]);if(area1==0.0)area1=1.0;
  PT->CrossVector(bc[3],bc[1],bc[2]);
  double area2 = 0.5*PT->Point3dSize(bc[3]);if(area2==0.0)area2=1.0;
  
  double lenik1 = PT->Point3dSize(bc[0]);
  double lenjk1 = PT->Distance(point[backID],point[nowID]);
  double lenij = PT->Point3dSize(bc[1]);
  double lenik2 = PT->Point3dSize(bc[2]);
  double lenjk2 = PT->Distance(point[nextID],point[nowID]);
  lenij *= lenij;
  
  ancos1 = (lenik1*lenik1 +lenjk1*lenjk1 - lenij)/area1; 
  
  ancos2 = (lenik2*lenik2 +lenjk2*lenjk2 - lenij)/area2;
    

  return (ancos1+ancos2);

}



void Polyhedron::setDesbrunC(){

  PolarList *nowp;
  int i;
  int nextID=0;
  int backID=0;
  double ddweight=0.0;
  for(i=0;i<numberV;i++){ 
    if(boundary[i]!=1){
      nowp =PHead[i];
      ddweight=0.0;
      while(next(nowp)!=PTail[i]){
	nowp = next(nowp);
	if(nowp->next!=PTail[i]){
	  nextID = nowp->next->ID;
	}else{
	  nextID = PHead[i]->next->ID;
	}
	if(nowp->back!=PHead[i]){
	  backID = nowp->back->ID;
	}else{
	  backID = PTail[i]->back->ID;
	}
	nowp->lambda = ((1.0-intrinsiclambda)*getChiTerm(i,backID,nowp->ID,nextID)+intrinsiclambda*getAuthalicTerm(i,backID,nowp->ID,nextID));
	ddweight += nowp->lambda;
      }
      
      if(ddweight!=0.0){
	nowp =PHead[i];
	while(next(nowp)!=PTail[i]){
	  nowp = next(nowp);
	  nowp->lambda /= ddweight;
	}
      }
    }
  }
}
double Polyhedron::getAuthalicTerm(int i,int backID,int nowID,int nextID){
  double ancos1,ancos2;
  PT->makeVector(bc[0],point[backID],point[i]);
  PT->makeVector(bc[1],point[backID],point[nowID]);
  PT->makeVector(bc[2],point[nextID],point[nowID]);
  PT->makeVector(bc[3],point[nextID],point[i]);


  
  ancos1 = PT->InnerProduct(bc[0],bc[1])/(PT->Point3dSize(bc[0])*PT->Point3dSize(bc[1]));
  
  ancos2 = PT->InnerProduct(bc[2],bc[3])/(PT->Point3dSize(bc[2])*PT->Point3dSize(bc[3]));
  
  return (ancos1/sqrt(1.0-ancos1*ancos1) +ancos2/sqrt(1.0-ancos2*ancos2));

}
double Polyhedron::getChiTerm(int i,int backID,int nowID,int nextID){
  double ancos1,ancos2,len;
  PT->makeVector(bc[0],point[nowID],point[i]);
  PT->makeVector(bc[1],point[nowID],point[backID]);
  PT->makeVector(bc[2],point[nowID],point[nextID]);
  
  
  len = PT->Point3dSize(bc[0]);if(len==0.0)len=1.0;
  ancos1 = PT->InnerProduct(bc[0],bc[1])/(len*PT->Point3dSize(bc[1]));
  ancos2 = PT->InnerProduct(bc[0],bc[2])/(len*PT->Point3dSize(bc[2]));
 
  return ((ancos1/sqrt(1.0-ancos1*ancos1) +ancos2/sqrt(1.0-ancos2*ancos2))/(len*len));
  
}
void Polyhedron::setMVCC(){
  PolarList *nowp;
  int i;
  int nextID=0;
  int backID=0;
  double ddweight=0.0;
  for(i=0;i<numberV;i++){ 
    if(boundary[i]!=1){
      nowp =PHead[i];
      ddweight=0.0;
      while(next(nowp)!=PTail[i]){
	nowp = next(nowp);
	if(nowp->next!=PTail[i]){
	  nextID = nowp->next->ID;
	}else{
	  nextID = PHead[i]->next->ID;
	}
	if(nowp->back!=PHead[i]){
	  backID = nowp->back->ID;
	}else{
	  backID = PTail[i]->back->ID;
	}
	nowp->lambda = getMVCTerm(i,backID,nowp->ID,nextID);
	ddweight += nowp->lambda;
      }
      
      if(ddweight!=0.0){
	nowp =PHead[i];
	while(next(nowp)!=PTail[i]){
	  nowp = next(nowp);
	  nowp->lambda /= ddweight;
	  
	}
      }
    }
  }
  
}
void Polyhedron::setFloaterC(){
  int i;
  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);
  IDList *now=NULL;
  PolarList *nowp=NULL;
  PolarList *nowp2=NULL;
  int checkintersect=0;
  double angle = 0.0;
  double theta = 0.0;
  double interval=0.0;
  double dweight=0.0;
  for(i=0;i<numberV;i++){
    if(boundary[i]!=1){
      dweight=0.0;
      nowp = PHead[i];
      while(next(nowp)!=PTail[i]){
	nowp = next(nowp);
	nowp2 = PHead[i];
	interval = nowp->theta+PI;
	if(interval>=(2.0*PI))interval -= (2.0*PI);
	checkintersect=0;
	while(next(next(nowp2))!=PTail[i]){
	  nowp2 = next(nowp2);
	  if((interval>=nowp2->theta)&&(interval<=next(nowp2)->theta)){
	    v1->x = nowp2->r*cos(nowp2->theta);
	    v1->y = nowp2->r*sin(nowp2->theta);
	    v2->x = next(nowp2)->r*cos(next(nowp2)->theta);
	    v2->y = next(nowp2)->r*sin(next(nowp2)->theta);
	    v3->x = nowp->r*cos(nowp->theta);
	    v3->y = nowp->r*sin(nowp->theta);
	    PT->setC(coef,v0,v1,v2,v3);
	    nowp->lambda += -coef[0]*coef[3];
	    nowp2->lambda += -coef[0]*coef[1];
	    nowp2->next->lambda += -coef[0]*coef[2];
	    dweight += -coef[0]*coef[0];
	    checkintersect=1;
	    break;
	  }
	}
	if(checkintersect==0){
	  nowp2 = next(nowp2);
	  v1->x = nowp2->r*cos(nowp2->theta);
	  v1->y = nowp2->r*sin(nowp2->theta);
	  v2->x = next(PHead[i])->r;
	  v2->y = 0.0;
	  v3->x = nowp->r*cos(nowp->theta);
	  v3->y = nowp->r*sin(nowp->theta);
	  PT->setC(coef,v0,v1,v2,v3);
	  nowp->lambda += -coef[0]*coef[3];
	  nowp2->lambda += -coef[0]*coef[1];
	  PHead[i]->next->lambda += -coef[0]*coef[2];
	  dweight += -coef[0]*coef[0];
	}
      }
      if(dweight!=0.0)
      nowp = PHead[i];
      while(next(nowp)!=PTail[i]){
	nowp = next(nowp);
	nowp->lambda /= -dweight;
      }
    }
  }
  delete coef;
  delete v0;
  delete v1;
  delete v2;
  delete v3;
  
}

double Polyhedron::getMVCTerm(int i,int backID,int nowID,int nextID){
   double angle1,angle2;
  PT->makeVector(bc[0],point[i],point[backID]);
  PT->makeVector(bc[1],point[i],point[nowID]);
  PT->makeVector(bc[2],point[i],point[nextID]);
  double len1 = PT->Point3dSize(bc[0]);
  double len2 = PT->Point3dSize(bc[1]);
  double len3 = PT->Point3dSize(bc[2]);
  angle1 = PT->InnerProduct(bc[0],bc[1])/(len1*len2);
  angle2 = PT->InnerProduct(bc[1],bc[2])/(len2*len3);
  double w1,w2;
  
  
  if((1.0+angle1)==0.0){
    w1 = 1.0;
  }else{
    w1 = sqrt(((1.0-angle1)/(1.0+angle1)));
  }
  if((1.0+angle2)==0.0){
    w2 = 1.0;
  }else{
    w2 = sqrt(((1.0-angle2)/(1.0+angle2)));
  }
  
    
  return ((w1+w2)/len2);
  
 }
 
void Polyhedron::SortIndexP(){
  int i;
  IDList *now=NULL;
  PolarList *nowp=NULL;
  
  PolarList *DummyHead = NULL;
  PolarList *DummyTail = NULL;
  
  for(i=0;i<numberV;i++){
    if(boundary[i]!=1){
      DummyHead = new PolarList();
      DummyTail = new PolarList();
      DummyHead->next = DummyTail;
      DummyTail->back = DummyHead;
      
      now = IHead[i];
      while(next(now)!=ITail[i]){
	now = next(now);
	nowp = PHead[i];
	while(next(nowp)!=PTail[i]){
	  nowp = next(nowp);
	  if(now->ID==nowp->ID){
	    
	    IDtool->AppendPolarI(nowp->ID,DummyTail,nowp->r,nowp->theta,nowp->lambda,nowp->cotw);
	    
	    
	    
	    break;
	  }
	}
      }
      
      IDtool->CleanNeighborPolar(PHead[i],PTail[i]);
      
      PHead[i] = new PolarList();
      PTail[i] = new PolarList();
      PHead[i]->next = PTail[i];
      PTail[i]->back = PHead[i];
      nowp = DummyHead;
      while(next(nowp)!=DummyTail){
	nowp = next(nowp);
	
	IDtool->AppendPolarI(nowp->ID,PTail[i],nowp->r,nowp->theta,nowp->lambda,nowp->cotw);
	
	
      }
      
      
      IDtool->CleanNeighborPolar(DummyHead,DummyTail);
      
    }
  }

}
void Polyhedron::SortIndexPNaturalB(){
  int i;
  IDList *now=NULL;
  PolarList *nowp=NULL;
  
  PolarList *DummyHead = NULL;
  PolarList *DummyTail = NULL;
  
  for(i=0;i<numberV;i++){
    if(PHead[i]->next!=PTail[i]){
      DummyHead = new PolarList();
      DummyTail = new PolarList();
      DummyHead->next = DummyTail;
      DummyTail->back = DummyHead;
      
      now = IHead[i];
      while(next(now)!=ITail[i]){
	now = next(now);
	nowp = PHead[i];
	while(next(nowp)!=PTail[i]){
	  nowp = next(nowp);
	  if(now->ID==nowp->ID){
	    IDtool->AppendPolarI(nowp->ID,DummyTail,nowp->r,nowp->theta,nowp->lambda,0.0);
	    break;
	  }
	}
      }
      
      IDtool->CleanNeighborPolar(PHead[i],PTail[i]);
      
      PHead[i] = new PolarList();
      PTail[i] = new PolarList();
      PHead[i]->next = PTail[i];
      PTail[i]->back = PHead[i];
      nowp = DummyHead;
      while(next(nowp)!=DummyTail){
	nowp = next(nowp);
	IDtool->AppendPolarI(nowp->ID,PTail[i],nowp->r,nowp->theta,nowp->lambda,0.0);
      }
      
      IDtool->CleanNeighborPolar(DummyHead,DummyTail);
      
    }
  }
  
}

  
void Polyhedron::setSigmaZero(){
  int i,j;
  IDList *now=NULL;
  double varphi,ddv,dsize1,sumarea;
  double dddhval=0.0;
  double localsum=0.0;
  for(i=0;i<numberF;i++){
    
    dsize1 = PT->getParametricA(pV[Face[i][0]],
				pV[Face[i][1]],
				pV[Face[i][2]],
				pU[Face[i][0]],
				pU[Face[i][1]],
				pU[Face[i][2]]);
    PT->setParametricDs(bc[0],point[Face[i][0]],
			point[Face[i][1]],point[Face[i][2]],
			pV[Face[i][0]],pV[Face[i][1]],
			pV[Face[i][2]],dsize1);
    PT->setParametricDt(bc[1],point[Face[i][0]],
			point[Face[i][1]],point[Face[i][2]],
			pU[Face[i][0]],pU[Face[i][1]],
			pU[Face[i][2]],dsize1);
    
    E[i] = PT->InnerProduct(bc[0],bc[0]);
    G[i] = PT->InnerProduct(bc[1],bc[1]);
  }
  for(i=0;i<numberV;i++){
         
    sigma[i]=0.0;
    now = FHead[i];
    varphi=0.0;
    localsum=0.0;
    
  
   
    
    
    
    while(next(now)!=FTail[i]){
      now = next(now);
      varphi += (areaMap3D[now->ID]*(0.5*(E[now->ID]+G[now->ID])));
      localsum += (areaMap3D[now->ID]);
      
      
    }
    
    sigma[i] = sqrt((varphi/localsum));
    
    sigma[i] = pow(sigma[i],gammaP);
    
  
  }
  
  
  
} 
  void Polyhedron::setLaplaceC(){
  PolarList *nowp;
  int i;
  
  for(i=0;i<numberV;i++){ 
    if(boundary[i]!=1){
      nowp =PHead[i];
      
      while(next(nowp)!=PTail[i]){
	nowp = next(nowp);
	nowp->lambda = 1.0/((double)(neighborI[i]));
      }
    }
  }
    
}

double Polyhedron::getCurrentE(){
  int i,j;
  IDList *now=NULL;
  double varphi,ddv,dsize1,sumarea;
  double dddhval=0.0;
  double dsum=0.0;
  
  
  
  double localsum=0.0;
  double pV1,pV2,pV3,pU1,pU2,pU3;
  double dE,dG;
  dE = 0.0;
  dG = 0.0;
  if(boundarysigma==0){
  for(i=0;i<numberF;i++){
    if(boundary[Face[i][0]]!=1&&boundary[Face[i][1]]!=1&&boundary[Face[i][2]]!=1){
    pV1 = pV[Face[i][0]];
    pV2 = pV[Face[i][1]];
    pV3 = pV[Face[i][2]];
    pU1 = pU[Face[i][0]];
    pU2 = pU[Face[i][1]];
    pU3 = pU[Face[i][2]];
   
    
    dsize1 = PT->getParametricA(pV1,pV2,pV3,pU1,pU2,pU3);
    
   
      
    PT->setParametricDs(bc[0],point[Face[i][0]],
			point[Face[i][1]],point[Face[i][2]],
			pV1,pV2,pV3,dsize1);
    PT->setParametricDt(bc[1],point[Face[i][0]],
			point[Face[i][1]],point[Face[i][2]],
			pU1,pU2,pU3,dsize1);
    dE = PT->InnerProduct(bc[0],bc[0]);
    
    dG = PT->InnerProduct(bc[1],bc[1]);
    dsum +=  areaMap3D[i]*0.5*(dE+dG);
    }
  }
  }else{
    for(i=0;i<numberF;i++){
      
    pV1 = pV[Face[i][0]];
    pV2 = pV[Face[i][1]];
    pV3 = pV[Face[i][2]];
    pU1 = pU[Face[i][0]];
    pU2 = pU[Face[i][1]];
    pU3 = pU[Face[i][2]];
   
    
    dsize1 = PT->getParametricA(pV1,pV2,pV3,pU1,pU2,pU3);
    
   
      
    PT->setParametricDs(bc[0],point[Face[i][0]],
			point[Face[i][1]],point[Face[i][2]],
			pV1,pV2,pV3,dsize1);
    PT->setParametricDt(bc[1],point[Face[i][0]],
			point[Face[i][1]],point[Face[i][2]],
			pU1,pU2,pU3,dsize1);
    dE = PT->InnerProduct(bc[0],bc[0]);
    
    dG = PT->InnerProduct(bc[1],bc[1]);
    dsum +=  areaMap3D[i]*0.5*(dE+dG);
   
  }
  
  
  }
  dsum =constsumarea3D*sqrt(dsum/sumarea3D);
  return dsum;

}

};