mesh_smoother_vsmoother.c

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

		  }
		  if(i==cells[l][1]){
			  if(mask[cells[l][0]]>0) {if(ind==0) {j=cells[l][5]; ind++;} else {k=cells[l][5]; ind++;}}
			  if(mask[cells[l][2]]>0) {if(ind==0) {j=cells[l][3]; ind++;} else {k=cells[l][3]; ind++;}}
		  }
		  if(i==cells[l][2]){
			  if(mask[cells[l][1]]>0) {if(ind==0) {j=cells[l][3]; ind++;} else {k=cells[l][3]; ind++;}}
			  if(mask[cells[l][0]]>0) {if(ind==0) {j=cells[l][4]; ind++;} else {k=cells[l][4]; ind++;}}
		  }
		  if(i==cells[l][3]){j=cells[l][1]; k=cells[l][2];}
	  if(i==cells[l][4]){j=cells[l][2]; k=cells[l][0];}
	  if(i==cells[l][5]){j=cells[l][0]; k=cells[l][1];}
		  break;}
	  }//end boundary connectivity
      //lines
      del1=R[j][0]-R[i][0]; del2=R[i][0]-R[k][0];
      if((fabs(del1)<eps)&&(fabs(del2)<eps)) {constr[I*4]=1; constr[I*4+1]=0;
	  constr[I*4+2]=R[i][0]; constr[I*4+3]=R[i][1];} else {
	  del1=R[j][1]-R[i][1]; del2=R[i][1]-R[k][1];
	  if((fabs(del1)<eps)&&(fabs(del2)<eps)) {constr[I*4]=0; constr[I*4+1]=1;
	      constr[I*4+2]=R[i][0]; constr[I*4+3]=R[i][1];
	      } else {
		  J=(R[i][0]-R[j][0])*(R[k][1]-R[j][1])-(R[k][0]-R[j][0])*(R[i][1]-R[j][1]);
		  if(fabs(J)<eps) {
		      constr[I*4]=1.0/(R[k][0]-R[j][0]); constr[I*4+1]=-1.0/(R[k][1]-R[j][1]);
		      constr[I*4+2]=R[i][0]; constr[I*4+3]=R[i][1];
} else //circle
{
		      x0=((R[k][1]-R[j][1])*(R[i][0]*R[i][0]-R[j][0]*R[j][0]+R[i][1]*R[i][1]-R[j][1]*R[j][1])+
			  (R[j][1]-R[i][1])*(R[k][0]*R[k][0]-R[j][0]*R[j][0]+R[k][1]*R[k][1]-R[j][1]*R[j][1]))/(2*J);
		      y0=((R[j][0]-R[k][0])*(R[i][0]*R[i][0]-R[j][0]*R[j][0]+R[i][1]*R[i][1]-R[j][1]*R[j][1])+
			  (R[i][0]-R[j][0])*(R[k][0]*R[k][0]-R[j][0]*R[j][0]+R[k][1]*R[k][1]-R[j][1]*R[j][1]))/(2*J);
		      constr[I*4]=x0; constr[I*4+1]=y0;
		      constr[I*4+2]=(R[i][0]-x0)*(R[i][0]-x0)+(R[i][1]-y0)*(R[i][1]-y0);
		  }
	      }
      }
  }
/*for(i=0;i<NCN;i++){
    for(j=0;j<4;j++) fprintf(stdout," %e ",constr[4*i+j]);
    fprintf(stdout, "constr %d node %d \n",i,Bind[i]);
}*/



//----initial values----
for(i=0;i<NCN;i++) lam[i]=0;
for(nz=0;nz<5;nz++){
  //---------find H and -grad J-----------------
   nonzero=0; Jpr=0;
   for(i=0; i<2*N; i++){ b[i] = 0; hm[i]=0; }
   for(i=0; i<ncells; i++){
	nvert=0; while(cells[i][nvert]>=0) nvert++;
	for(j=0;j<nvert;j++){ G[i][j]=0;
	if(adp<0) for(k=0;k<abs(adp);k++) G[i][j]=G[i][j]+afun[i*(-adp)+k]; }
	for(index=0;index<2;index++){
	    for(k=0;k<nvert;k++){ F[index][k]=0;
		for(j=0;j<nvert;j++) W[index][k][j]=0;
	    }
	}
	if(mcells[i]>=0){
		     Jpr+=localP(2, W, F, R, cells[i], mask, epsilon, w, nvert, H[i],
			 me, vol, 0, &lVmin, &lqmin, adp, afun, G[i], sout);
	} else {
	    for(index=0;index<2;index++) for(j=0;j<nvert;j++) W[index][j][j]=1;
	}
	for(index=0;index<2;index++){
	    for(l=0; l<nvert; l++){//-diagonal Hessian-
		hm[cells[i][l]+index*N]=hm[cells[i][l]+index*N]+W[index][l][l];
		b[cells[i][l]+index*N]=b[cells[i][l]+index*N]-F[index][l];
	    }
	}
    }
    //------||grad J||_2------
    for(i=0;i<2*N;i++){ nonzero=nonzero+b[i]*b[i];}
    //-----solve for Plam--------
    for(I=0;I<NCN;I++){i=Bind[I];
       if(constr[4*I+3]<0.5/eps) {Bx=constr[4*I]; By=constr[4*I+1];
	   g=(R[i][0]-constr[4*I+2])*constr[4*I]+(R[i][1]-constr[4*I+3])*constr[4*I+1];
       } else {
	  Bx=2*(R[i][0]-constr[4*I]); By=2*(R[i][1]-constr[4*I+1]);
	      hm[i]+=2*lam[I]; hm[i+N]+=2*lam[I];
	      g=(R[i][0]-constr[4*I])*(R[i][0]-constr[4*I])+(R[i][1]-constr[4*I+1])*(R[i][1]-constr[4*I+1])-constr[4*I+2];
       }
       Jpr+=lam[I]*g;
       qq=Bx*b[i]/hm[i]+By*b[i+N]/hm[i+N]-g;
       a=Bx*Bx/hm[i]+By*By/hm[i+N];
       if(a!=0) Plam[I]=qq/a; else fprintf(sout,"error: B^TH-1B is degenerate \n");
       b[i]-=Plam[I]*Bx; b[i+N]-=Plam[I]*By;
       Plam[I]-=lam[I];
   }
   //-----------solve for P------------
   for(i=0;i<N;i++) {P[i][0]=b[i]/hm[i]; P[i][1]=b[i+N]/hm[i+N];}
   //-----correct solution-----
   for(i=0;i<N;i++) for(j=0;j<2;j++) if((fabs(P[i][j])<eps)||(mask[i]==1)) P[i][j]=0;
   //----P is determined--------
   if(msglev>=3){
       for(i=0;i<N;i++){
	   for(j=0;j<2;j++)  if(P[i][j]!=0) fprintf(sout, "P[%d][%d]=%f  ",i,j,P[i][j]);
       }
   }
   //-------local minimization problem, determination of tau----------
   if(msglev>=3) fprintf(sout, "dJ=%e L0=%e \n",sqrt(nonzero), Jpr);
   L=1e32; j=1;
   while((Jpr<=L)&&(j>-30)){
      j=j-1;
      tau=pow(2.0,j);
      for(i=0;i<N;i++){
	  for(k=0;k<2;k++) Rpr[i][k]=R[i][k]+tau*P[i][k];}
      J=0; gVmin=1e32; gemax=-1e32; gqmin=1e32;
      for(i=0; i<ncells; i++) if(mcells[i]>=0){
	  nvert=0; while(cells[i][nvert]>=0) nvert++;
		  lemax=localP(2, W, F, Rpr, cells[i], mask, epsilon, w, nvert, H[i], me, vol, 1, &lVmin,
		       &lqmin, adp, afun, G[i], sout);
	      J+=lemax;
	  if(gVmin>lVmin) gVmin=lVmin;
	  if(gemax<lemax) gemax=lemax;
	  if(gqmin>lqmin) gqmin=lqmin;
	  }
      L=J;
      //----constraints contribution----
      for(I=0;I<NCN;I++){i=Bind[I];
	  if(constr[4*I+3]<0.5/eps) g=(Rpr[i][0]-constr[4*I+2])*constr[4*I]+(Rpr[i][1]-constr[4*I+3])*constr[4*I+1];
	      else g=(Rpr[i][0]-constr[4*I])*(Rpr[i][0]-constr[4*I])+
		   (Rpr[i][1]-constr[4*I+1])*(Rpr[i][1]-constr[4*I+1])-constr[4*I+2];
	  L+=(lam[I]+tau*Plam[I])*g;
	  }
      //----end of constraints----
      if(msglev>=3) fprintf(sout," tau=%f J=%f \n",tau,J);
   }
   if(j==-30) T=0; else {
      Jpr=L; qq=J;
      for(i=0;i<N;i++){
	  for(k=0;k<2;k++) Rpr[i][k]=R[i][k]+tau*0.5*P[i][k];}
      J=0; gVmin0=1e32; gemax0=-1e32; gqmin0=1e32;
      for(i=0; i<ncells; i++) if(mcells[i]>=0){
	  nvert=0; while(cells[i][nvert]>=0) nvert++;
		  lemax=localP(2, W, F, Rpr, cells[i], mask, epsilon, w, nvert, H[i], me, vol, 1, &lVmin,
		       &lqmin, adp, afun, G[i], sout);
		  J+=lemax;
	      if(gVmin0>lVmin) gVmin0=lVmin;
	  if(gemax0<lemax) gemax0=lemax;
	  if(gqmin0>lqmin) gqmin0=lqmin;
	  }
      L=J;
      //----constraints contribution----
      for(I=0;I<NCN;I++){i=Bind[I];
	  if(constr[4*I+3]<0.5/eps) g=(Rpr[i][0]-constr[4*I+2])*constr[4*I]+(Rpr[i][1]-constr[4*I+3])*constr[4*I+1];
	      else g=(Rpr[i][0]-constr[4*I])*(Rpr[i][0]-constr[4*I])+
		     (Rpr[i][1]-constr[4*I+1])*(Rpr[i][1]-constr[4*I+1])-constr[4*I+2];
	  L+=(lam[I]+tau*0.5*Plam[I])*g;
	  }
      //----end of constraints----
   }
   if(Jpr>L) {
      T=0.5*tau;
      (*Vmin)=gVmin0;
      (*emax)=gemax0;
      (*qmin)=gqmin0;
   } else {
      T=tau; L=Jpr; J=qq;
      (*Vmin)=gVmin;
      (*emax)=gemax;
      (*qmin)=gqmin;
   }

   for(j=0;j<N;j++) for(k=0;k<2;k++) R[j][k]+=T*P[j][k];
   for(i=0;i<NCN;i++) lam[i]+=T*Plam[i];

}//end Lagrangian iter

 if(msglev>=2) fprintf(sout, "tau=%e, J=%e, L=%e  \n",T,J,L);

free(lam);
free(b);
for(i=0;i<2;i++){
   for(j=0;j<6;j++) free(W[i][j]);
   free(W[i]);
   free(F[i]);}
free(W);
free(F);
for(i=0;i<N;i++) {free(Rpr[i]);
    free(P[i]);}
free(Rpr);
free(P);
for(i=0;i<ncells;i++) free(G[i]);
free(G);
free(Bind);
free(constr);
free(hm);
free(Plam);

return (sqrt(nonzero));

}

/**
 * composes local matrix W and right side F from all quadrature nodes of one cell
 */
double VariationalMeshSmoother::localP(int n, LPLPLPDOUBLE W, LPLPDOUBLE F, LPLPDOUBLE R, LPINT cell, LPINT mask, double epsilon,
	      double w, int nvert, LPLPDOUBLE H, int me, double vol, int f, double *Vmin,
	      double *qmin, int adp, LPDOUBLE afun, LPDOUBLE Gloc, FILE *sout)
{

  int  ii, jj, kk, i, j, k, l, m, nn;
  double  sigma=0.0, fun, lqmin, gqmin, g;
  double K[9];
  /*f - flag, f=0 for determination of Hessian and gradient of the functional,
    f=1 for determination of the functional value only;
    K - determines approximation rule for local integral over the cell*/


  /*-------adaptivity, determined on the first step for adp>0 (nodal based)------*/
  if(f==0)
  {
    if(adp>0)
      avertex(n, afun, Gloc, R, cell, nvert, adp, sout);
    if(adp==0)
    {
      for(i=0;i<n;i++)
	Gloc[i]=1.0;
    }
  }

fun=0; gqmin=1e32; g=0;//Vmin

//cell integration depending on cell type
if(n==2){//2D
	if(nvert==3){//tri
		sigma=1.0;
	fun+=vertex(n, W, F, R, cell, epsilon, w, nvert, K, H, me, vol, f, &lqmin, adp, Gloc, sigma, sout);
		g+=sigma*lqmin;
	if(gqmin>lqmin) gqmin=lqmin;
	}
    if(nvert==4){//quad
	for(i=0; i<2; i++){ K[0]=i;
	    for(j=0; j<2; j++){ K[1]=j;
			    sigma=0.25;
			    fun+=vertex(n, W, F, R, cell, epsilon, w, nvert, K, H, me,
				vol, f, &lqmin, adp, Gloc, sigma, sout);
		g+=sigma*lqmin;
		if(gqmin>lqmin) gqmin=lqmin;
	    }
	  }
	} else {//quad tri
	    for(i=0; i<3; i++){ K[0]=i*0.5; k=i/2; K[1]=(double)k;
	       for(j=0; j<3; j++){  K[2]=j*0.5; k=j/2; K[3]=(double)k;
			       if(i==j) sigma=1.0/12; else sigma=1.0/24;
			       fun+=vertex(n, W, F, R, cell, epsilon, w, nvert, K, H, me,
				   vol, f, &lqmin, adp, Gloc, sigma, sout);
			   g+=sigma*lqmin;
		   if(gqmin>lqmin) gqmin=lqmin;
			   }
			}
	}
}
if(n==3){//3D
       if(nvert==4){//tetr
		   sigma=1.0;
		   fun+=vertex(n, W, F, R, cell, epsilon, w, nvert, K, H, me,
			   vol, f, &lqmin, adp, Gloc, sigma, sout);
		   g+=sigma*lqmin;
	   if(gqmin>lqmin) gqmin=lqmin;
       }
       if(nvert==6){//prism
	  for(i=0;i<2;i++){ K[0]=i;
	      for(j=0;j<2;j++){ K[1]=j;
		  for(k=0;k<3;k++) {K[2]=(double)k/2.0; K[3]=(double)(k%2);
				      sigma=1.0/12.0;
				      fun+=vertex(n, W, F, R, cell, epsilon, w, nvert, K, H, me,
				      vol, f, &lqmin, adp, Gloc, sigma, sout);
			      g+=sigma*lqmin;
		      if(gqmin>lqmin) gqmin=lqmin;
				  }
			  }
		  }
	   }
       if(nvert==8){//hex
	 for(i=0; i<2; i++){ K[0]=i;
	     for(j=0; j<2; j++){ K[1]=j;
		 for(k=0; k<2; k++){ K[2]=k;
		     for(l=0; l<2; l++){ K[3]=l;
			 for(m=0; m<2; m++){ K[4]=m;
			     for(nn=0; nn<2; nn++){ K[5]=nn;
							    if((i==nn)&&(j==l)&&(k==m)) sigma=(double)1/27;
				if(((i==nn)&&(j==l)&&(k!=m))||((i==nn)&&(j!=l)&&(k==m))||((i!=nn)&&(j==l)&&(k==m))) sigma=(double)1/(27*2);
				if(((i==nn)&&(j!=l)&&(k!=m))||((i!=nn)&&(j!=l)&&(k==m))||((i!=nn)&&(j==l)&&(k!=m))) sigma=(double)1/(27*4);
				if((i!=nn)&&(j!=l)&&(k!=m)) sigma=(double)1/(27*8);
							    fun+=vertex(n, W, F, R, cell, epsilon, w, nvert, K, H, me,
						vol, f, &lqmin, adp, Gloc, sigma, sout);
					g+=sigma*lqmin;
				if(gqmin>lqmin) gqmin=lqmin;
							 }
						 }
					 }
				 }
			 }
		 }
	   } else {//quad tetr
	  for(i=0;i<4;i++){
	    for(j=0;j<4;j++){
		  for(k=0;k<4;k++){
		    switch (i)
		    { case 0 : K[0]=0; K[1]=0; K[2]=0; break;
		      case 1 : K[0]=1; K[1]=0; K[2]=0; break;
		      case 2 : K[0]=1.0/2; K[1]=1; K[2]=0; break;
		      case 3 : K[0]=1.0/2; K[1]=1.0/3; K[2]=1; break; }
		    switch (j)
		    { case 0 : K[3]=0; K[4]=0; K[5]=0; break;
		      case 1 : K[3]=1; K[4]=0; K[5]=0; break;
		      case 2 : K[3]=1.0/2; K[4]=1; K[5]=0; break;
		      case 3 : K[3]=1.0/2; K[4]=1.0/3; K[5]=1; break; }
		    switch (k)
		    { case 0 : K[6]=0; K[7]=0; K[8]=0; break;
		      case 1 : K[6]=1; K[7]=0; K[8]=0; break;
		      case 2 : K[6]=1.0/2; K[7]=1; K[8]=0; break;
		      case 3 : K[6]=1.0/2; K[7]=1.0/3; K[8]=1; break; }
			   if((i==j)&&(j==k)) sigma=1.0/120; else
		       if((i==j)||(j==k)||(i==k)) sigma=1.0/360; else sigma=1.0/720;
			   fun+=vertex(n, W, F, R, cell, epsilon, w, nvert, K, H, me,
			       vol, f, &lqmin, adp, Gloc, sigma, sout);
		       g+=sigma*lqmin;
	       if(gqmin>lqmin) gqmin=lqmin;
			  }
			}
		  }
	   }
}

/*---fixed nodes correction---*/
for(ii=0;ii<nvert;ii++)
{
    if(mask[cell[ii]]==1)
    {
      for(kk=0;kk<n;kk++)
      {
	for(jj=0;jj<nvert;jj++)
	{
	  W[kk][ii][jj]=0;
	  W[kk][jj][ii]=0;
	}

	W[kk][ii][ii]=1;
	F[kk][ii]=0;
      }
    }
}
/*---end of fixed nodes correction---*/

(*Vmin)=g;
(*qmin)=gqmin/vol;

return fun;
}

/**
 * avertex - assembly of adaptivity metric on a cell
 */
// double VariationalMeshSmoother::avertex(int n, LPDOUBLE afun, LPDOUBLE G, LPLPDOUBLE R, LPINT cell, int nvert, int adp, FILE *sout)
double VariationalMeshSmoother::avertex(int n, LPDOUBLE afun, LPDOUBLE G, LPLPDOUBLE R, LPINT cell, int nvert, int adp, FILE *)
{
  LPLPDOUBLE Q;
  LPDOUBLE K;
  double a1[3], a2[3], a3[3], qu[3];
  int i,j;
  double det, g, df0, df1, df2;
  K=alloc_d_n1(8);
  Q = alloc_d_n1_n2(3, nvert);
  for(i=0;i<n;i++) Q[i]=alloc_d_n1(nvert);

  for(i=0;i<8;i++) K[i]=0.5; //cell center

  basisA(n,Q,nvert,K,Q,1);
  for(i=0;i<n;i++){a1[i]=0; a2[i]=0; a3[i]=0; qu[i]=0;
     for(j=0;j<nvert;j++){
	a1[i]+=Q[i][j]*R[cell[j]][0];
	a2[i]+=Q[i][j]*R[cell[j]][1];
		if(n==3) a3[i]+=Q[i][j]*R[cell[j]][2];
		qu[i]+=Q[i][j]*afun[cell[j]];
     }
  }
  if(n==2) det=jac2(a1[0],a1[1],a2[0],a2[1]); else det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
  if(det!=0){
	  if(n==2){
	  df0=jac2(qu[0],qu[1],a2[0],a2[1])/det;
	  df1=jac2(a1[0],a1[1],qu[0],qu[1])/det;
		  g=(1+df0*df0+df1*df1);
		  if(adp==2){//directional adaptation G=diag(g_i)
			  G[0]=1+df0*df0; G[1]=1+df1*df1;
		  } else for(i=0;i<n;i++) G[i]=g; //simple adaptation G=gI
	  } else {
		  df0=(qu[0]*(a2[1]*a3[2]-a2[2]*a3[1])+qu[1]*(a2[2]*a3[0]-a2[0]*a3[2])+qu[2]*(a2[0]*a3[1]-a2[1]*a3[0]))/det;
	  df1=(qu[0]*(a3[1]*a1[2]-a3[2]*a1[1])+qu[1]*(a3[2]*a1[0]-a3[0]*a1[2])+qu[2]*(a3[0]*a1[1]-a3[1]*a1[0]))/det;
	  df2=(qu[0]*(a1[1]*a2[2]-a1[2]*a2[1])+qu[1]*(a1[2]*a2[0]-a1[0]*a2[2])+qu[2]*(a1[0]*a2[1]-a1[1]*a2[0]))/det;
		  g=(1+df0*df0+df1*df1+df2*df2);
	  if(adp==2){//directional adaptation G=diag(g_i)
			  G[0]=1+df0*df0; G[1]=1+df1*df1; G[2]=1+df2*df2;
		  } else for(i=0;i<n;i++) G[i]=g; //simple adaptation G=gI
	  }
  } else {g=1.0; for(i=0;i<n;i++) G[i]=g;}


 for(i=0;i<n;i++) free(Q[i]);
 free(Q);

 return(g);

}

/**
 * Computes local matrics W and local rhs F on one basis
 */
double VariationalMeshSmoother::vertex(int n, LPLPLPDOUBLE W, LPLPDOUBLE F, LPLPDOUBLE R, LPINT cell,
	      double epsilon, double w, int nvert, LPDOUBLE K, LPLPDOUBLE H,
		  int me, double vol, int f, double *qmin, int adp, LPDOUBLE g, double sigma, FILE *)
//	          int me, double vol, int f, double *qmin, int adp, LPDOUBLE g, double sigma, FILE *sout)
{
  LPLPLPDOUBLE P = NULL, d2phi = NULL;
  LPLPDOUBLE  gpr = NULL, dphi = NULL, dfe = NULL;
  LPLPDOUBLE Q;
  double a1[3], a2[3], a3[3], av1[3], av2[3], av3[3];
  int i,j,k,i1;
  double tr=0.0, det=0.0, dchi, chi, fet, phit=0.0, G, con=100.0;
  Q = alloc_d_n1_n2(3, nvert);
  for(i=0;i<n;i++) Q[i]=alloc_d_n1(nvert);
if(f==0){
  gpr = alloc_d_n1_n2(3, nve