mesh_smoother_vsmoother.c

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

  }

  free(Gamma);

  /*-----BN correction - 2D only!---*/
  epsilon=0.000000001;
  if(NBN>0) for(counter=0;counter<iter[2];counter++)
  {
    //update adaptation function before each iteration
    if((adp!=0)&&(gr==0))
      adp_renew(n, N, R, ncells, cells, afun, adp, sout);
    Jk=minJ_BC(N, R, mask, ncells, cells, mcells, eps, w, me, H, vol, msglev, &Vmin, &emax, &qmin, adp, afun, NBN, sout);
    if(msglev>=1)
      fprintf(sout, "NBC niter=%d, qmin*G/vol=%e, Vmin=%e, emax=%e  \n",counter,qmin,Vmin,emax);

    //Outrageous Enm1 to make sure we hit this atleast once
    Enm1=99999;

    //Now that we've moved the boundary nodes (or not) we need to resmoooth
    for(j=0;(j<iter[1]);j++)
    {
      if(fabs(emax-Enm1)<1e-2)
	break;

      //Save off the error from the previous smoothing step
      Enm1=emax;

      //update adaptation function before each iteration
      if((adp!=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, adp, afun, sout);

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

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

  /*----------free memory-----------------*/
  //free(maskf);
  //if(adp!=0) free(afun);


  return;
}

/**
 * Determines the values of maxE_theta
 */
// double VariationalMeshSmoother::maxE(int n, int N, LPLPDOUBLE R, int ncells, LPLPINT cells, LPINT mcells,
//             int me, LPLPLPDOUBLE H, double v, double epsilon, double w, LPDOUBLE Gamma,
//			double *qmin, FILE *sout)
double VariationalMeshSmoother::maxE(int n, int, LPLPDOUBLE R, int ncells, LPLPINT cells, LPINT mcells,
	    int me, LPLPLPDOUBLE H, double v, double epsilon, double w, LPDOUBLE Gamma,
			double *qmin, FILE *)
{
  LPLPDOUBLE Q;
  double K[9], a1[3], a2[3], a3[3];
  double  gemax, det, tr, E=0.0, sigma=0.0, chi, vmin;
  int  ii, i, j, k, l, m, nn, kk, ll;

  Q = alloc_d_n1_n2(3, 10);
  for(i=0;i<n;i++) Q[i]=alloc_d_n1(3*n+n%2);

  gemax=-1e32; vmin=1e32;

for(ii=0; ii<ncells; ii++) if(mcells[ii]>=0){
    if(n==2){
		if(cells[ii][3]==-1){//tri
	  basisA(2,Q,3,K,H[ii],me);
	  for(k=0;k<2;k++){a1[k]=0; a2[k]=0;
	     for(l=0;l<3;l++){
		 a1[k]=a1[k]+Q[k][l]*R[cells[ii][l]][0];
		 a2[k]=a2[k]+Q[k][l]*R[cells[ii][l]][1];
	     }
	  }
	  det=jac2(a1[0],a1[1],a2[0],a2[1]);
	  tr=0.5*(a1[0]*a1[0]+a2[0]*a2[0]+a1[1]*a1[1]+a2[1]*a2[1]);
	  chi=0.5*(det+sqrt(det*det+epsilon*epsilon));
	  E=(1-w)*tr/chi+0.5*w*(v+det*det/v)/chi;
	  if(E>gemax) gemax=E;
		  if(vmin>det) vmin=det;
		} else if(cells[ii][4]==-1){ E=0; //quad
	  for(i=0; i<2; i++){ K[0]=i;
	    for(j=0; j<2; j++){ K[1]=j;
			   basisA(2,Q,4,K,H[ii],me);
	       for(k=0;k<2;k++){a1[k]=0; a2[k]=0;
		  for(l=0;l<4;l++){
		      a1[k]=a1[k]+Q[k][l]*R[cells[ii][l]][0];
		      a2[k]=a2[k]+Q[k][l]*R[cells[ii][l]][1];
		  }
	       }
	       det=jac2(a1[0],a1[1],a2[0],a2[1]);
	       tr=0.5*(a1[0]*a1[0]+a2[0]*a2[0]+a1[1]*a1[1]+a2[1]*a2[1]);
	       chi=0.5*(det+sqrt(det*det+epsilon*epsilon));
	       E+=0.25*((1-w)*tr/chi+0.5*w*(v+det*det/v)/chi);
			   if(vmin>det) vmin=det;
			}
	  }
	  if(E>gemax) gemax=E;
		} else { E=0; //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;
			      basisA(2,Q,6,K,H[ii],me);
		  for(k=0;k<2;k++){a1[k]=0; a2[k]=0;
		     for(l=0;l<6;l++){
			a1[k]=a1[k]+Q[k][l]*R[cells[ii][l]][0];
			a2[k]=a2[k]+Q[k][l]*R[cells[ii][l]][1];
		     }
		  }
		  det=jac2(a1[0],a1[1],a2[0],a2[1]);
		  if(i==j) sigma=1.0/12; else sigma=1.0/24;
		  tr=0.5*(a1[0]*a1[0]+a2[0]*a2[0]+a1[1]*a1[1]+a2[1]*a2[1]);
		  chi=0.5*(det+sqrt(det*det+epsilon*epsilon));
		  E+=sigma*((1-w)*tr/chi+0.5*w*(v+det*det/v)/chi);
				  if(vmin>det) vmin=det;
			   }
			}
	    if(E>gemax) gemax=E;
		}
	}
    if(n==3){
       if(cells[ii][4]==-1){//tetr
	   basisA(3,Q,4,K,H[ii],me);
	   for(k=0;k<3;k++){a1[k]=0; a2[k]=0; a3[k]=0;
	      for(l=0;l<4;l++){
		  a1[k]=a1[k]+Q[k][l]*R[cells[ii][l]][0];
		  a2[k]=a2[k]+Q[k][l]*R[cells[ii][l]][1];
		  a3[k]=a3[k]+Q[k][l]*R[cells[ii][l]][2];
	      }
	   }
	   det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
	   tr=0; for(k=0;k<3;k++) tr+=(a1[k]*a1[k]+a2[k]*a2[k]+a3[k]*a3[k])/3.0;
	   chi=0.5*(det+sqrt(det*det+epsilon*epsilon));
	   E=(1-w)*pow(tr,1.5)/chi+0.5*w*(v+det*det/v)/chi;
	   if(E>gemax) gemax=E;
		   if(vmin>det) vmin=det;
       } else if(cells[ii][6]==-1){//prism
	  E=0;
	  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);
		      basisA(3,Q,6,K,H[ii],me);
		      for(kk=0;kk<3;kk++){a1[kk]=0; a2[kk]=0; a3[kk]=0;
			  for(ll=0;ll<6;ll++){
			      a1[kk]=a1[kk]+Q[kk][ll]*R[cells[ii][ll]][0];
			      a2[kk]=a2[kk]+Q[kk][ll]*R[cells[ii][ll]][1];
			      a3[kk]=a3[kk]+Q[kk][ll]*R[cells[ii][ll]][2];
			  }
		      }
		      det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
		      tr=0; for(kk=0;kk<3;kk++) tr+=(a1[kk]*a1[kk]+a2[kk]*a2[kk]+a3[kk]*a3[kk])/3.0;
		      chi=0.5*(det+sqrt(det*det+epsilon*epsilon));
		      E+=((1-w)*pow(tr,1.5)/chi+0.5*w*(v+det*det/v)/chi)/12.0;
					  if(vmin>det) vmin=det;
				  }
			  }
		  }
	  if(E>gemax) gemax=E;
       } else if(cells[ii][8]==-1){ E=0; //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;
				basisA(3,Q,8,K,H[ii],me);
				for(kk=0;kk<3;kk++){a1[kk]=0; a2[kk]=0; a3[kk]=0;
				   for(ll=0;ll<8;ll++){
				      a1[kk]=a1[kk]+Q[kk][ll]*R[cells[ii][ll]][0];
				      a2[kk]=a2[kk]+Q[kk][ll]*R[cells[ii][ll]][1];
				      a3[kk]=a3[kk]+Q[kk][ll]*R[cells[ii][ll]][2];
				   }
				}
				det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
								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);
				tr=0; for(kk=0;kk<3;kk++) tr+=(a1[kk]*a1[kk]+a2[kk]*a2[kk]+a3[kk]*a3[kk])/3.0;
				chi=0.5*(det+sqrt(det*det+epsilon*epsilon));
				E+=((1-w)*pow(tr,1.5)/chi+0.5*w*(v+det*det/v)/chi)*sigma;
								if(vmin>det) vmin=det;
							 }
			 }
		      }
		  }
	      }
	  }
	  if(E>gemax) gemax=E;
       } else {//quad tetr
		   E=0;
	  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; }
		    basisA(3,Q,10,K,H[ii],me);
		    for(kk=0;kk<3;kk++){a1[kk]=0; a2[kk]=0; a3[kk]=0;
		       for(ll=0;ll<10;ll++){
			  a1[kk]=a1[kk]+Q[kk][ll]*R[cells[ii][ll]][0];
			  a2[kk]=a2[kk]+Q[kk][ll]*R[cells[ii][ll]][1];
			  a3[kk]=a3[kk]+Q[kk][ll]*R[cells[ii][ll]][2];
		       }
		    }
		    det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
		    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;
		    tr=0; for(kk=0;kk<3;kk++) tr+=(a1[kk]*a1[kk]+a2[kk]*a2[kk]+a3[kk]*a3[kk])/3.0;
		    chi=0.5*(det+sqrt(det*det+epsilon*epsilon));
		    E+=((1-w)*pow(tr,1.5)/chi+0.5*w*(v+det*det/v)/chi)*sigma;
					if(vmin>det) vmin=det;
			  }
			}
		  }
	  if(E>gemax) gemax=E;
	   }
	}
	Gamma[ii]=E;
}

(*qmin)=vmin;

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

return gemax;
}

/**
  * Compute min Jacobian determinant (minq), min cell volume (Vmin), and average cell volume (vol).
  */
// double VariationalMeshSmoother::minq(int n, int N, LPLPDOUBLE R, int ncells, LPLPINT cells, LPINT mcells,
//                                      int me, LPLPLPDOUBLE H, double *vol, double *Vmin, FILE *sout)
double VariationalMeshSmoother::minq(int n, int, LPLPDOUBLE R, int ncells, LPLPINT cells, LPINT mcells,
				     int me, LPLPLPDOUBLE H, double *vol, double *Vmin, FILE *)
{
  LPLPDOUBLE Q;
  double K[9], a1[3], a2[3], a3[3];
  double  v, vmin, gqmin, det, vcell, sigma=0.0;
  int  ii, i, j, k, l, m, nn, kk, ll;

  Q = alloc_d_n1_n2(3, 10);
  for(i=0;i<n;i++) Q[i]=alloc_d_n1(3*n+n%2);

  v=0; vmin=1e32; gqmin=1e32;

for(ii=0; ii<ncells; ii++) if(mcells[ii]>=0){
    if(n==2){//2D
	if(cells[ii][3]==-1){//tri
	  basisA(2,Q,3,K,H[ii],me);
	  for(k=0;k<2;k++){a1[k]=0; a2[k]=0;
	     for(l=0;l<3;l++){
		 a1[k]=a1[k]+Q[k][l]*R[cells[ii][l]][0];
		 a2[k]=a2[k]+Q[k][l]*R[cells[ii][l]][1];
	     }
	  }
	  det=jac2(a1[0],a1[1],a2[0],a2[1]);
	  if(gqmin>det) gqmin=det;
	  if(vmin>det) vmin=det;
	  v+=det;
	} else if(cells[ii][4]==-1){ vcell=0; //quad
	  for(i=0; i<2; i++){ K[0]=i;
	    for(j=0; j<2; j++){ K[1]=j;
	       basisA(2,Q,4,K,H[ii],me);
	       for(k=0;k<2;k++){a1[k]=0; a2[k]=0;
		  for(l=0;l<4;l++){
		      a1[k]=a1[k]+Q[k][l]*R[cells[ii][l]][0];
		      a2[k]=a2[k]+Q[k][l]*R[cells[ii][l]][1];
		  }
	       }
	       det=jac2(a1[0],a1[1],a2[0],a2[1]);
	       if(gqmin>det) gqmin=det;
			   v+=0.25*det;
			   vcell+=0.25*det;
	    }
	  }
	  if(vmin>vcell) vmin=vcell;
		} else { vcell=0; //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;
			      basisA(2,Q,6,K,H[ii],me);
		  for(k=0;k<2;k++){a1[k]=0; a2[k]=0;
		     for(l=0;l<6;l++){
			a1[k]=a1[k]+Q[k][l]*R[cells[ii][l]][0];
			a2[k]=a2[k]+Q[k][l]*R[cells[ii][l]][1];
		     }
		  }
		  det=jac2(a1[0],a1[1],a2[0],a2[1]);
		  if(gqmin>det) gqmin=det;
				  if(i==j) sigma=1.0/12; else sigma=1.0/24;
				  v+=sigma*det;
				  vcell+=sigma*det;
			   }
			}
			if(vmin>vcell) vmin=vcell;
		}
	}
    if(n==3){//3D
       if(cells[ii][4]==-1){//tetr
	   basisA(3,Q,4,K,H[ii],me);
	   for(k=0;k<3;k++){a1[k]=0; a2[k]=0; a3[k]=0;
	      for(l=0;l<4;l++){
		  a1[k]=a1[k]+Q[k][l]*R[cells[ii][l]][0];
		  a2[k]=a2[k]+Q[k][l]*R[cells[ii][l]][1];
		  a3[k]=a3[k]+Q[k][l]*R[cells[ii][l]][2];
	      }
	   }
	   det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
	   if(gqmin>det) gqmin=det;
		   if(vmin>det) vmin=det;
	   v+=det;
       } else if(cells[ii][6]==-1){//prism
	  vcell=0;
	  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);
		      basisA(3,Q,6,K,H[ii],me);
		      for(kk=0;kk<3;kk++){a1[kk]=0; a2[kk]=0; a3[kk]=0;
			  for(ll=0;ll<6;ll++){
			      a1[kk]=a1[kk]+Q[kk][ll]*R[cells[ii][ll]][0];
			      a2[kk]=a2[kk]+Q[kk][ll]*R[cells[ii][ll]][1];
			      a3[kk]=a3[kk]+Q[kk][ll]*R[cells[ii][ll]][2];
			  }
		      }
		      det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
		      if(gqmin>det) gqmin=det; sigma=1.0/12.0;
		      v+=sigma*det; vcell+=sigma*det;
				  }
			  }
		  }
	  if(vmin>vcell) vmin=vcell;
	   } else if(cells[ii][8]==-1){ vcell=0; //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;
				basisA(3,Q,8,K,H[ii],me);
				for(kk=0;kk<3;kk++){a1[kk]=0; a2[kk]=0; a3[kk]=0;
				   for(ll=0;ll<8;ll++){
				      a1[kk]=a1[kk]+Q[kk][ll]*R[cells[ii][ll]][0];
				      a2[kk]=a2[kk]+Q[kk][ll]*R[cells[ii][ll]][1];
				      a3[kk]=a3[kk]+Q[kk][ll]*R[cells[ii][ll]][2];
				   }
				}
				det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
				if(gqmin>det) gqmin=det;
				    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);
				v+=sigma*det;
				vcell+=sigma*det;
			    }
			 }
		      }
		  }
	      }
	  }
	  if(vmin>vcell) vmin=vcell;
       } else {//quad tetr
	  vcell=0;
	  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; }
		    basisA(3,Q,10,K,H[ii],me);
		    for(kk=0;kk<3;kk++){a1[kk]=0; a2[kk]=0; a3[kk]=0;
		       for(ll=0;ll<10;ll++){
			  a1[kk]=a1[kk]+Q[kk][ll]*R[cells[ii][ll]][0];
			  a2[kk]=a2[kk]+Q[kk][ll]*R[cells[ii][ll]][1];
			  a3[kk]=a3[kk]+Q[kk][ll]*R[cells[ii][ll]][2];
		       }
		    }
		    det=jac3(a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],a3[0],a3[1],a3[2]);
		    if(gqmin>det) gqmin=det;
		    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;
		    v+=sigma*det;   vcell+=sigma*det;
			  }
			  }
		  }
		  if(vmin>vcell) vmin=vcell;
       }
    }
}


(*vol)=v/(double)ncells;
(*Vmin)=vmin;