mesh_smoother_vsmoother.c

一个用来实现偏微分方程中网格的计算库

    {

      std::cout<<"Parent 1: "<<(it->second)[0]<<std::endl;
      std::cout<<"Parent 2: "<<(it->second)[1]<<std::endl;
      std::cout<<"Hanging Node: "<<it->first<<std::endl<<std::endl;


      //First Parent
      edges[2*i]=(it->second)[1];

      //Second Parent
      edges[2*i+1]=(it->second)[0];

      //Hanging Node
      hnodes[i]=it->first;

      i++;
    }
  }
  std::cout<<"Finished readgr"<<std::endl;

return 0;
}

/**
 * Read Metrics
 */
//int VariationalMeshSmoother::readmetr(char *name, LPLPLPDOUBLE H, int ncells, int n, FILE *sout)
int VariationalMeshSmoother::readmetr(char *name, LPLPLPDOUBLE H, int ncells, int n, FILE *)
{
int i, j, k;
double d;
FILE *stream;

   stream=fopen(name,"r");
   for(i=0;i<ncells;i++)
	  for(j=0;j<n;j++){
	  for(k=0;k<n;k++){fscanf(stream,"%le ",&d); H[i][j][k]=d;}
		  fscanf(stream,"\n");
	  }

   fclose(stream);

return 0;
}

// Stolen from ErrorVector!
float VariationalMeshSmoother::adapt_minimum() const
{
  std::vector<float>& adapt_data = *_adapt_data;

  const unsigned int n = adapt_data.size();
  float min = 1.e30;

  for (unsigned int i=0; i<n; i++)
  {
      // Only positive (or zero) values in the error vector
    libmesh_assert(adapt_data[i] >= 0.);
    min = std::min (min, adapt_data[i]);
  }

  // ErrorVectors are for positive values
  libmesh_assert (min >= 0.);

  return min;
}

void VariationalMeshSmoother::adjust_adapt_data()
{
  //For convenience
  const UnstructuredMesh & aoe_mesh=*_area_of_interest;
  std::vector<float>& adapt_data = *_adapt_data;

  float min=adapt_minimum();

  MeshBase::const_element_iterator       el     = _mesh.elements_begin();
  const MeshBase::const_element_iterator end_el = _mesh.elements_end();

  MeshBase::const_element_iterator       aoe_el     = aoe_mesh.elements_begin();
  const MeshBase::const_element_iterator aoe_end_el = aoe_mesh.elements_end();

  //Counter to keep track of which spot in adapt_data we're looking at
  int i=0;

  for(; el != end_el; el++)
  {
    //Only do this for active elements
    if(adapt_data[i])
    {
      Point centroid=(*el)->centroid();
      bool in_aoe=false;

      //See if the elements centroid lies in the aoe mesh
      //for(aoe_el=aoe_mesh.elements_begin(); aoe_el != aoe_end_el; ++aoe_el)
      //{
	if((*aoe_el)->contains_point(centroid))
	  in_aoe=true;
      //}

      //If the element is not in the area of interest... then set
      //it's adaptivity value to the minimum.
      if(!in_aoe)
	adapt_data[i]=min;
    }
    i++;
  }
}

/**
 * Read Adaptivity
 */
// int VariationalMeshSmoother::read_adp(LPDOUBLE afun, char adap[], FILE *sout)
int VariationalMeshSmoother::read_adp(LPDOUBLE afun, char [], FILE *)
{
  int i, m, j;

  std::vector<float>& adapt_data = *_adapt_data;

  if(_area_of_interest)
    adjust_adapt_data();

  m=adapt_data.size();

  j=0;
  for(i=0;i<m;i++)
  {
    if(adapt_data[i]!=0)
    {
      afun[j]=adapt_data[i];
      j++;
    }
  }

  return 0;
}

double VariationalMeshSmoother::jac3(double x1,double y1,double z1,double x2,double y2,
	    double z2,double x3,double y3,double z3)
{
  return( x1*(y2*z3 - y3*z2) + y1*(z2*x3 - z3*x2) + z1*(x2*y3 - x3*y2) );
}

double VariationalMeshSmoother::jac2(double x1,double y1,double x2,double y2)
{
  return( x1*y2-x2*y1 );
}

/**
 * BasisA determines matrix H^(-T)Q on one Jacobian matrix
 */
int VariationalMeshSmoother::basisA(int n, LPLPDOUBLE Q, int nvert, LPDOUBLE K, LPLPDOUBLE H, int me)
{
 LPLPDOUBLE U;
 int i,j,k;
 U=alloc_d_n1_n2(n,nvert);
 for(i=0;i<n;i++) U[i]=alloc_d_n1(nvert);

 if(n==2){
    if(nvert==4){//quad
	U[0][0]=(-(1-K[1]));
	U[0][1]=( (1-K[1]));
	U[0][2]=(-    K[1]);
	U[0][3]=(     K[1]);
	U[1][0]=(-(1-K[0]));
	U[1][1]=(-    K[0]);
	U[1][2]=( (1-K[0]));
	U[1][3]=(     K[0]);
    }
    if(nvert==3){//tri
	/*---for right target triangle---
	U[0][0]=-1.0; U[1][0]=-1.0;
	U[0][1]=1.0;  U[1][1]=0.0;
	U[0][2]=0.0;  U[1][2]=1.0;
	-----for regular triangle---*/
	U[0][0]=-1.0; U[1][0]=-1.0/sqrt(3.0);
	U[0][1]= 1.0; U[1][1]=-1.0/sqrt(3.0);
	U[0][2]=   0; U[1][2]= 2.0/sqrt(3.0);
    }
    if(nvert==6){/*--curved triangle*/
	U[0][0]=-3*(1-K[0]-K[1]*0.5)+(K[0]-0.5*K[1])+K[1];
	U[0][1]=-(1-K[0]-K[1]*0.5)+(K[0]-0.5*K[1])*3-K[1];
	U[0][2]=0;
	U[0][3]=K[1]*4;
	U[0][4]=-4*K[1];
	U[0][5]=4*(1-K[0]-K[1]*0.5)-4*(K[0]-0.5*K[1]);
	U[1][0]=(1-K[2]-K[3]*0.5)*(-1.5)+(K[2]-0.5*K[3])*(0.5)+K[3]*(0.5);
	U[1][1]=(1-K[2]-K[3]*0.5)*(0.5)+(K[2]-0.5*K[3])*(-1.5)+K[3]*(0.5);
	U[1][2]=(1-K[2]-K[3]*0.5)*(-1)+(K[2]-0.5*K[3])*(-1)+K[3]*(3);
	U[1][3]=(K[2]-0.5*K[3])*(4.0)+K[3]*(-2.0);
	U[1][4]=(1-K[2]-K[3]*0.5)*(4.0)+K[3]*(-2.0);
	U[1][5]=(1-K[2]-K[3]*0.5)*(-2.0)+(K[2]-0.5*K[3])*(-2.0);
    }
 }
 if(n==3){
     if(nvert==8){//hex
       U[0][0]=(-(1-K[0])*(1-K[1]));
       U[0][1]=( (1-K[0])*(1-K[1]));
       U[0][2]=(-    K[0]*(1-K[1]));
       U[0][3]=(     K[0]*(1-K[1]));
       U[0][4]=(-(1-K[0])*    K[1]);
       U[0][5]=( (1-K[0])*    K[1]);
       U[0][6]=(-    K[0]*    K[1]);
       U[0][7]=(     K[0]*    K[1]);
       U[1][0]=(-(1-K[2])*(1-K[3]));
       U[1][1]=(-    K[2]*(1-K[3]));
       U[1][2]=( (1-K[2])*(1-K[3]));
       U[1][3]=(     K[2]*(1-K[3]));
       U[1][4]=(-(1-K[2])*    K[3]);
       U[1][5]=(-    K[2]*    K[3]);
       U[1][6]=( (1-K[2])*    K[3]);
       U[1][7]=(     K[2]*    K[3]);
       U[2][0]=(-(1-K[4])*(1-K[5]));
       U[2][1]=(-    K[4]*(1-K[5]));
       U[2][2]=(-(1-K[4])*    K[5]);
       U[2][3]=(-    K[4]*    K[5]);
       U[2][4]=( (1-K[4])*(1-K[5]));
       U[2][5]=(     K[4]*(1-K[5]));
       U[2][6]=( (1-K[4])*    K[5]);
       U[2][7]=(     K[4]*    K[5]);
     }
     if(nvert==4){//linear tetr
      /*---for regular reference tetrahedron--------*/
      U[0][0]=-1; U[1][0]=-1.0/sqrt(3.0); U[2][0]=-1.0/sqrt(6.0);
      U[0][1]=1;  U[1][1]=-1.0/sqrt(3.0); U[2][1]=-1.0/sqrt(6.0);
      U[0][2]=0;  U[1][2]=2.0/sqrt(3.0);  U[2][2]=-1.0/sqrt(6.0);
      U[0][3]=0;  U[1][3]=0;            U[2][3]=3.0/sqrt(6.0);
      /*---for right corner reference tetrahedron----
      U[0][0]=-1; U[1][0]=-1; U[2][0]=-1;
      U[0][1]=1;  U[1][1]=0;  U[2][1]=0;
      U[0][2]=0;  U[1][2]=1;  U[2][2]=0;
      U[0][3]=0;  U[1][3]=0;  U[2][3]=1;*/
     }
     if(nvert==6){//prism
	/* for regular triangle in the prism base */
	U[0][0]=-1+K[0]; U[1][0]=(-1+K[1])/2.0; U[2][0]=-1+K[2]+K[3]/2.0;
	U[0][1]=1-K[0]; U[1][1]=(-1+K[1])/2.0; U[2][1]=-K[2]+K[3]/2.0;
	U[0][2]=0; U[1][2]=(1-K[1]); U[2][2]=-K[3];
	U[0][3]=-K[0]; U[1][3]=(-K[1])/2.0; U[2][3]=1-K[2]-K[3]/2.0;
	U[0][4]=K[0]; U[1][4]=(-K[1])/2.0; U[2][4]=K[2]-K[3]/2.0;
	U[0][5]=0; U[1][5]=K[1]; U[2][5]=K[3];
     }
     if(nvert==10){//quad tetr
	 U[0][0]=(1-K[0]-K[1]/2-K[2]/3)*(-3)+(K[0]-K[1]/2-K[2]/3)+(K[1]-K[2]/3)+K[2];
	 U[0][1]=(1-K[0]-K[1]/2-K[2]/3)*(-1)+(K[0]-K[1]/2-K[2]/3)*3-(K[1]-K[2]/3)-K[2];
	 U[0][2]=0; U[0][3]=0;
	 U[0][4]=4*(K[1]-K[2]/3); U[0][5]=-4*(K[1]-K[2]/3);
	 U[0][6]=4*(1-K[0]-K[1]/2-K[2]/3)-4*(K[0]-K[1]/2-K[2]/3);
	 U[0][7]=4*K[2]; U[0][8]=0; U[0][9]=-4*K[2];
	 U[1][0]=(1-K[3]-K[4]/2-K[5]/3)*(-3/2)+(K[3]-K[4]/2-K[5]/3)/2+(K[4]-K[5]/3)/2+K[5]/2;
	 U[1][1]=(1-K[3]-K[4]/2-K[5]/3)/2+(K[3]-K[4]/2-K[5]/3)*(-3/2)+(K[4]-K[5]/3)/2+K[5]/2;
	 U[1][2]=-(1-K[3]-K[4]/2-K[5]/3)-(K[3]-K[4]/2-K[5]/3)+3*(K[4]-K[5]/3)-K[5];
	 U[1][3]=0; U[1][4]=4*(K[3]-K[4]/2-K[5]/3)-2*(K[4]-K[5]/3);
	 U[1][5]=4*(1-K[3]-K[4]/2-K[5]/3)-2*(K[4]-K[5]/3);
	 U[1][6]=-2*(1-K[3]-K[4]/2-K[5]/3)-2*(K[3]-K[4]/2-K[5]/3);
	 U[1][7]=-2*K[5]; U[1][8]=4*K[5]; U[1][9]=-2*K[5];
	 U[2][0]=-(1-K[6]-K[7]/2-K[8]/3)+(K[6]-K[7]/2-K[8]/3)/3+(K[7]-K[8]/3)/3+K[8]/3;
	 U[2][1]=(1-K[6]-K[7]/2-K[8]/3)/3-(K[6]-K[7]/2-K[8]/3)+(K[7]-K[8]/3)/3+K[8]/3;
	 U[2][2]=(1-K[6]-K[7]/2-K[8]/3)/3+(K[6]-K[7]/2-K[8]/3)/3-(K[7]-K[8]/3)+K[8]/3;
	 U[2][3]=-(1-K[6]-K[7]/2-K[8]/3)-(K[6]-K[7]/2-K[8]/3)-(K[7]-K[8]/3)+3*K[8];
	 U[2][4]=-4*(K[6]-K[7]/2-K[8]/3)/3-4*(K[7]-K[8]/3)/3;
	 U[2][5]=-4*(1-K[6]-K[7]/2-K[8]/3)/3-4*(K[7]-K[8]/3)/3;
	 U[2][6]=-4*(1-K[6]-K[7]/2-K[8]/3)/3-4*(K[6]-K[7]/2-K[8]/3)/3;
	 U[2][7]=4*(K[6]-K[7]/2-K[8]/3)-4*K[8]/3;
	 U[2][8]=4*(K[7]-K[8]/3)-4*K[8]/3; U[2][9]=4*(1-K[6]-K[7]/2-K[8]/3)-4*K[8]/3;
     }
 }

 if(me==1){
    for(i=0;i<n;i++)
	for(j=0;j<nvert;j++) Q[i][j]=U[i][j];

 } else {
    for(i=0;i<n;i++){
       for(j=0;j<nvert;j++){ Q[i][j]=0;
	  for(k=0;k<n;k++) Q[i][j]+=H[k][i]*U[k][j];
       }
    }
 }

 for(i=0;i<n;i++) free(U[i]);
 free(U);
return 0;
}

/**
 * Specify adaptive function
 */
// void VariationalMeshSmoother::adp_renew(int n, int N, LPLPDOUBLE R, int ncells, LPLPINT cells, LPDOUBLE afun, int adp, FILE *sout)
void VariationalMeshSmoother::adp_renew(int n, int N, LPLPDOUBLE R, int ncells, LPLPINT cells, LPDOUBLE afun, int adp, FILE *)
{
  int i, j, nvert;
  double x, y, z;

  //evaluates given adaptive function on the provided mesh
  if(adp<0){
    for(i=0;i<ncells;i++)
    {
      x=y=z=0;
      nvert=0;
      while(cells[i][nvert]>=0)
      {
	x+=R[cells[i][nvert]][0];
	y+=R[cells[i][nvert]][1];
	if(n==3)
	  z+=R[cells[i][nvert]][2];
	nvert++;
      }
      afun[i]=5*(x+y+z); //adaptive function, cell based
    }
  }
  if(adp>0)
  {
    for(i=0;i<N;i++)
    {
      z=0; for(j=0;j<n;j++) z+=R[i][j];
      afun[i]=5*sin(R[i][0]); //adaptive function, node based
    }
  }

  return;
}

/**
 * Preprocess mesh data and control smoothing/untangling iterations
 */
void VariationalMeshSmoother::full_smooth(int n, int N, LPLPDOUBLE R, LPINT mask, int ncells, LPLPINT cells, LPINT mcells,
		 int nedges, LPINT edges, LPINT hnodes, double w, LPINT iter, int me,
		     LPLPLPDOUBLE H, int adp, char adap[], int gr, FILE *sout)
{
  LPDOUBLE afun=NULL, Gamma;
  LPINT maskf;
  double  Jk, epsilon, eps, qmin, vol, emax, Vmin, Enm1;
  int i, ii, j, counter, NBN, ladp, msglev=1;

  //Adaptive function is on cells
  if(adp<0)
    afun=alloc_d_n1(ncells);

  //Adaptive function is on nodes
  if(adp>0)
    afun=alloc_d_n1(N);

  maskf=alloc_i_n1(N);
  Gamma=alloc_d_n1(ncells);

  if(msglev>=1)
    fprintf(sout,"N=%d ncells=%d nedges=%d \n",N,ncells,nedges);


  //Boundary node counting
  NBN=0;
  for(i=0;i<N;i++)
    if(mask[i]==2 || mask[i]==1)
      NBN++;

  if(NBN>0)
  {
    if(msglev>=1) fprintf(sout,"# of Boundary Nodes=%d  \n",NBN);

    NBN=0;
    for(i=0;i<N;i++)
      if(mask[i]==2)
	NBN++;
    if(msglev>=1) fprintf(sout,"# of moving Boundary Nodes=%d  \n",NBN);
  }

  for(i=0;i<N;i++)
  {
    if((mask[i]==1)||(mask[i]==2))
      maskf[i]=1;
    else
      maskf[i]=0;
  }

  /*-------determination of min jacobian-------*/
  qmin=minq(n, N, R, ncells, cells, mcells, me, H, &vol, &Vmin, sout);
  if(me>1) vol=1.0;
  if(msglev>=1) fprintf(sout,"vol=%e  qmin=%e min volume = %e\n",vol,qmin,Vmin);

  epsilon=0.000000001;
  //compute max distortion measure over all cells
  if(qmin<0){eps=sqrt(epsilon*epsilon+0.004*qmin*qmin*vol*vol);} else eps=epsilon;
  emax=maxE(n, N, R, ncells, cells, mcells, me, H, vol, eps, w, Gamma, &qmin, sout);
  if(msglev>=1) fprintf(sout," emax=%e \n",emax);

  /*-------unfolding/smoothing-----------*/

  //iteration tolerance
  ii=0; counter=0;
  Enm1=1.0;

  //read adaptive function from file
  if(adp*gr!=0) read_adp(afun, adap, sout);

  while(((qmin<=0)||(counter<iter[0])||(fabs(emax-Enm1)>1e-3))&&(ii<iter[1])&&(counter<iter[1]))
  {
    libmesh_assert(counter<iter[1]);

    Enm1=emax;

    if((ii>=0)&&(ii%2==0))
    {
      if(qmin<0)
	eps=sqrt(epsilon*epsilon+0.004*qmin*qmin*vol*vol);
      else
	eps=epsilon;
    }

    if((qmin<=0)||(counter<ii)) ladp=0; else ladp=adp;

    //update adaptation function before each iteration
    if((ladp!=0)&&(gr==0))
      adp_renew(n, N, R, ncells, cells, afun, adp, sout);

    Jk=minJ(n, N, R, maskf, ncells, cells, mcells, eps, w, me, H, vol, nedges, edges, hnodes,
	    msglev, &Vmin, &emax, &qmin, ladp, afun, sout);

    if(qmin>0)
      counter++;
    else
      ii++;

    if(msglev>=1)
    {
      fprintf(sout, "niter=%d, qmin*G/vol=%e, Vmin=%e, emax=%e  Jk=%e \n",counter,qmin,Vmin,emax, Jk);
      fprintf(sout," emax=%e, Enm1=%e \n",emax,Enm1);
    }