mesh_smoother_vsmoother.c

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

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

return gqmin;
}

/**
 * Executes one step of minimization algorithm:
 * finds minimization direction (P=H^{-1} \grad J) and solves approximately
 * local minimization problem for optimal step in this minimization direction (tau=min J(R+tau P))
 */
double VariationalMeshSmoother::minJ(int n, int N, LPLPDOUBLE R, LPINT mask, int ncells, LPLPINT cells, LPINT mcells,
	    double epsilon, double w, int me, LPLPLPDOUBLE H, double vol, int nedges,
	    LPINT edges, LPINT hnodes, int msglev, double *Vmin, double *emax, double *qmin,
	    int adp, LPDOUBLE afun, FILE *sout)
{

LPLPLPDOUBLE W; //local Hessian matrix;
LPLPDOUBLE  F, A, G; //F - local gradient; G - adaptation metric;
LPDOUBLE b, u, a; // matrix & rhs for solver, u - solution vector;
LPINT ia, ja; // matrix connectivity for solver;
LPLPINT JA; //A, JA - internal form of global matrix;
LPLPDOUBLE Rpr, P; //P - minimization direction;
double  tau=0.0, J, T, Jpr, lVmin, lemax, lqmin, gVmin=0.0, gemax=0.0,
gqmin=0.0, gtmin0=0.0, gtmax0=0.0, gqmin0=0.0;
double eps, nonzero, Tau_hn, g_i;
int index, i, ii, j, jj, k, l, m, nz;
int sch, ind, columns, nvert, ind_i, ind_j, ind_k;
/* Jpr - value of functional;
   nonzero - norm of gradient;
   columns - max number of nonzero entries in every row of global matrix;
*/

columns=n*n*10;
W=alloc_d_n1_n2_n3(n,3*n+n%2,3*n+n%2);
F = alloc_d_n1_n2(n,3*n+n%2);
for(i=0;i<n;i++){
   F[i]=alloc_d_n1(3*n+n%2);
   W[i]=alloc_d_n1_n2(3*n+n%2,3*n+n%2);
   for(j=0;j<3*n+n%2;j++) W[i][j]=alloc_d_n1(3*n+n%2);}
Rpr=alloc_d_n1_n2(N,n);
P=alloc_d_n1_n2(N,n);
for(i=0;i<N;i++) {Rpr[i]=alloc_d_n1(n);
 P[i]=alloc_d_n1(n);}
A = alloc_d_n1_n2(n*N, columns);
b = alloc_d_n1(n*N);
u = alloc_d_n1(n*N);
a = alloc_d_n1(n*N*columns);
ia = alloc_i_n1(n*N+1);
ja = alloc_i_n1(n*N*columns);
JA = alloc_i_n1_n2(n*N, columns);
for(i=0;i<n*N;i++) {
   A[i] = alloc_d_n1(columns);
   JA[i] = alloc_i_n1(columns);}
G=alloc_d_n1_n2(ncells,n);
for(i=0;i<ncells;i++) G[i]=alloc_d_n1(n);

//---------find minimization direction P-----------------
nonzero=0; Jpr=0;
for(i=0; i<n*N; i++){//initialise matrix and rhs
    for(j=0; j<columns; j++){ A[i][j] = 0; JA[i][j] =0;}
    b[i] = 0;
    }
for(i=0; i<ncells; i++)
{
  nvert=0;
  while(cells[i][nvert]>=0)
    nvert++;
  //---determination of local matrices on each cell----

  for(j=0;j<n;j++)
  {
    G[i][j]=0; //adaptation metric G is held constant throughout minJ run
    if(adp<0)
    {
      for(k=0;k<abs(adp);k++)
	G[i][j]=G[i][j]+afun[i*(-adp)+k]; //cell-based adaptivity is computed here
    }
  }
     for(index=0;index<n;index++){//initialise local matrices
	 for(k=0;k<3*n+n%2;k++){ F[index][k]=0;
	     for(j=0;j<3*n+n%2;j++) W[index][k][j]=0;
	 }
     }
     if(mcells[i]>=0){//if cell is not excluded
	 Jpr+=localP(n, 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<n;index++) for(j=0;j<nvert;j++) W[index][j][j]=1;
     }
    //----assembly of an upper triangular part of a global matrix A------
    for(index=0;index<n;index++)
    {
      for(l=0; l<nvert; l++)
      {
	for(m=0; m<nvert; m++){
	  if((W[index][l][m]!=0)&&(cells[i][m]>=cells[i][l]))
	  {
	    sch=0; ind=1;
	    while (ind!=0)
	    {
	      if(A[cells[i][l]+index*N][sch]!=0)
	      {
		  if(JA[cells[i][l]+index*N][sch]==(cells[i][m]+index*N))
		  {
		    A[cells[i][l]+index*N][sch] = A[cells[i][l]+index*N][sch] + W[index][l][m];
		    ind=0;
		  }
		  else
		    sch++;
	      }
	      else
	      {
		A[cells[i][l]+index*N][sch] = W[index][l][m];
		JA[cells[i][l]+index*N][sch]=cells[i][m]+index*N;
		ind = 0;
	      }
	      if(sch>columns-1) fprintf(sout,"error: # of nonzero entries in the %d row of Hessian =%d >= columns=%d \n",cells[i][l],sch,columns);
	    }
	  }
	}
	b[cells[i][l]+index*N]=b[cells[i][l]+index*N]-F[index][l];
      }
    }
//------end of matrix A---------
  }
 //-------HN correction--------
 Tau_hn=pow(vol,1.0/(double)n)*1e-10; //tolerance for HN being the mid-edge node for its parents
 for(i=0;i<nedges;i++){
     ind_i=hnodes[i]; ind_j=edges[2*i]; ind_k=edges[2*i+1];
     for(j=0;j<n;j++){
	 g_i=R[ind_i][j]-0.5*(R[ind_j][j]+R[ind_k][j]);
	 Jpr+=g_i*g_i/(2*Tau_hn);
	 b[ind_i+j*N]-=g_i/Tau_hn;
	 b[ind_j+j*N]+=0.5*g_i/Tau_hn;
	 b[ind_k+j*N]+=0.5*g_i/Tau_hn;
     }
     for(ind=0;ind<columns;ind++){
	 if(JA[ind_i][ind]==ind_i) A[ind_i][ind]+=1.0/Tau_hn;
	 if(JA[ind_i+N][ind]==ind_i+N) A[ind_i+N][ind]+=1.0/Tau_hn;
	 if(n==3) if(JA[ind_i+2*N][ind]==ind_i+2*N) A[ind_i+2*N][ind]+=1.0/Tau_hn;
	 if((JA[ind_i][ind]==ind_j)||(JA[ind_i][ind]==ind_k)) A[ind_i][ind]-=0.5/Tau_hn;
	 if((JA[ind_i+N][ind]==ind_j+N)||(JA[ind_i+N][ind]==ind_k+N)) A[ind_i+N][ind]-=0.5/Tau_hn;
	 if(n==3) if((JA[ind_i+2*N][ind]==ind_j+2*N)||(JA[ind_i+2*N][ind]==ind_k+2*N)) A[ind_i+2*N][ind]-=0.5/Tau_hn;
	 if(JA[ind_j][ind]==ind_i) A[ind_j][ind]-=0.5/Tau_hn;
	 if(JA[ind_k][ind]==ind_i) A[ind_k][ind]-=0.5/Tau_hn;
	 if(JA[ind_j+N][ind]==ind_i+N) A[ind_j+N][ind]-=0.5/Tau_hn;
	 if(JA[ind_k+N][ind]==ind_i+N) A[ind_k+N][ind]-=0.5/Tau_hn;
	 if(n==3) if(JA[ind_j+2*N][ind]==ind_i+2*N) A[ind_j+2*N][ind]-=0.5/Tau_hn;
	 if(n==3) if(JA[ind_k+2*N][ind]==ind_i+2*N) A[ind_k+2*N][ind]-=0.5/Tau_hn;
	 if((JA[ind_j][ind]==ind_j)||(JA[ind_j][ind]==ind_k)) A[ind_j][ind]+=0.25/Tau_hn;
	 if((JA[ind_k][ind]==ind_j)||(JA[ind_k][ind]==ind_k)) A[ind_k][ind]+=0.25/Tau_hn;
	 if((JA[ind_j+N][ind]==ind_j+N)||(JA[ind_j+N][ind]==ind_k+N)) A[ind_j+N][ind]+=0.25/Tau_hn;
	 if((JA[ind_k+N][ind]==ind_j+N)||(JA[ind_k+N][ind]==ind_k+N)) A[ind_k+N][ind]+=0.25/Tau_hn;
	 if(n==3) if((JA[ind_j+2*N][ind]==ind_j+2*N)||(JA[ind_j+2*N][ind]==ind_k+2*N)) A[ind_j+2*N][ind]+=0.25/Tau_hn;
	 if(n==3) if((JA[ind_k+2*N][ind]==ind_j+2*N)||(JA[ind_k+2*N][ind]==ind_k+2*N)) A[ind_k+2*N][ind]+=0.25/Tau_hn;
     }
 }
//------||\grad J||_2------
for(i=0;i<n*N;i++){ nonzero=nonzero+b[i]*b[i];}

//-----sort matrix A-----
  for(i=0;i<n*N;i++)
  {
      for(j=columns-1;j>1;j--)
      {
	  for(k=0;k<j;k++)
	  {
	     if(JA[i][k]>JA[i][k+1])
	     {
	       sch=JA[i][k];
	       JA[i][k]=JA[i][k+1];
	       JA[i][k+1]=sch;
	       tau=A[i][k];
	       A[i][k]=A[i][k+1];
	       A[i][k+1]=tau;
	     }
	 }
      }
  }

eps=sqrt(vol)*1e-9;

//-------solver for P (unconstrained)-------
  ia[0]=0; j=0;
  for(i=0;i<n*N;i++)
  {
      u[i]=0;
      nz=0;
      for(k=0;k<columns;k++){
	 if(A[i][k]!=0)
	 {
	   nz++;
	   ja[j]=JA[i][k]+1;
	   a[j]=A[i][k];
	   j++;
	 }
      }
      ia[i+1]=ia[i]+nz;
  }
  nz=ia[n*N];
  m=n*N;
  if(msglev>=3) sch=1; else sch=0;


  nz=solver(m,ia,ja,a,u,b,eps,100,sch, sout);
//    sol_pcg_pj(m,ia,ja,a,u,b,eps,100,sch);

  for(i=0;i<N;i++){ //ensure fixed nodes are not moved
      for(index=0;index<n;index++)
	  if(mask[i]==1) P[i][index]=0; else P[i][index]=u[i+index*N];
  }
 //----P is determined--------
 if(msglev>=4){
 for(i=0;i<N;i++){
    for(j=0;j<n;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 J0=%e \n",sqrt(nonzero),Jpr);
J=1e32; j=1;
while((Jpr<=J)&&(j>-30)){
   j=j-1;
   tau=pow(2.0,j); //search through finite # of values for tau
   for(i=0;i<N;i++)
   {
       for(k=0;k<n;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(n, 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;

	//------HN correction--------
	for(ii=0;ii<nedges;ii++)
	{
	    ind_i=hnodes[ii];
	    ind_j=edges[2*ii];
	    ind_k=edges[2*ii+1];
	    for(jj=0;jj<n;jj++)
	    {
		g_i=Rpr[ind_i][jj]-0.5*(Rpr[ind_j][jj]+Rpr[ind_k][jj]);
		J+=g_i*g_i/(2*Tau_hn);
	    }
	 }
      }
   }
   if(msglev>=3)
     fprintf(sout, "tau=%f J=%f \n",tau,J);
}

if(j==-30)
  T=0;
else
{
  Jpr=J;
  for(i=0;i<N;i++)
  {
    for(k=0;k<n;k++)
      Rpr[i][k]=R[i][k]+tau*0.5*P[i][k];
  }
  J=0; gtmin0=1e32; gtmax0=-1e32; gqmin0=1e32;
  for(i=0; i<ncells; i++)
  {
    if(mcells[i]>=0)
    {
      nvert=0; while(cells[i][nvert]>=0) nvert++;
      lemax=localP(n, W, F, Rpr, cells[i], mask, epsilon, w, nvert, H[i], me, vol, 1, &lVmin, &lqmin, adp, afun, G[i], sout);
      J+=lemax;
      if(gtmin0>lVmin)
	gtmin0=lVmin;
      if(gtmax0<lemax)
	gtmax0=lemax;
      if(gqmin0>lqmin)
	gqmin0=lqmin;
      //-------HN correction--------
       for(ii=0;ii<nedges;ii++){
	   ind_i=hnodes[ii]; ind_j=edges[2*ii]; ind_k=edges[2*ii+1];
	   for(jj=0;jj<n;jj++){
	       g_i=Rpr[ind_i][jj]-0.5*(Rpr[ind_j][jj]+Rpr[ind_k][jj]);
	       J+=g_i*g_i/(2*Tau_hn);
	   }
       }//HN
   }
  }
}
 if(Jpr>J) {
    T=0.5*tau;
    (*Vmin)=gtmin0;
    (*emax)=gtmax0;
    (*qmin)=gqmin0;
 } else {
    T=tau; J=Jpr;
    (*Vmin)=gVmin;
    (*emax)=gemax;
    (*qmin)=gqmin;
 }

 nonzero=0;
for(j=0;j<N;j++) for(k=0;k<n;k++){
     R[j][k]=R[j][k]+T*P[j][k];
     nonzero+=T*P[j][k]*T*P[j][k];
}

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

free(b);
free(u);
free(ia);
free(ja);
free(a);
for(i=0;i<n;i++){
   for(j=0;j<3*n+n%2;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<n*N;i++) {
    free(A[i]);
    free(JA[i]);}
free(A);
free(JA);
for(i=0;i<ncells;i++) free(G[i]);
free(G);

return (sqrt(nonzero));

}

/**
 * minJ() with sliding Boundary Nodes constraints and no account for HN,
 * using Lagrange multiplier formulation: minimize L=J+\sum lam*g;
 * only works in 2D
 */
double VariationalMeshSmoother::minJ_BC(int N, LPLPDOUBLE R, LPINT mask, int ncells, LPLPINT cells, LPINT mcells,
	       double epsilon, double w, int me, LPLPLPDOUBLE H, double vol, int msglev,
		   double *Vmin, double *emax, double *qmin, int adp, LPDOUBLE afun, int NCN, FILE *sout)
{
/*---new form of matrices, 5 iterations for minL---*/
LPLPLPDOUBLE W;
LPLPDOUBLE  F, G;
LPDOUBLE b, hm, Plam, constr, lam;
LPINT Bind;
LPLPDOUBLE Rpr, P;
double  tau=0.0, J=0.0, T, Jpr, L, lVmin, lqmin, gVmin=0.0, gqmin=0.0, gVmin0=0.0,
gqmin0=0.0, lemax, gemax=0.0, gemax0=0.0;
double a, g, qq=0.0, eps, nonzero, x0, y0, del1, del2, Bx, By;
int index, i, j, k=0, l, nz, I;
int ind, nvert;

//----memory------
Bind=alloc_i_n1(NCN); //array of sliding BN
lam=alloc_d_n1(NCN);
hm=alloc_d_n1(2*N);
Plam=alloc_d_n1(NCN);
constr=alloc_d_n1(4*NCN); //holds constraints = local approximation to the boundary
F = alloc_d_n1_n2(2, 6);
W=alloc_d_n1_n2_n3(2, 6, 6);
for(i=0;i<2;i++){
   F[i]=alloc_d_n1(6);
   W[i]=alloc_d_n1_n2(6,6);
   for(j=0;j<6;j++) W[i][j]=alloc_d_n1(6);}
Rpr=alloc_d_n1_n2(N,2);
P=alloc_d_n1_n2(N,2);
for(i=0;i<N;i++) {Rpr[i]=alloc_d_n1(2);
 P[i]=alloc_d_n1(2);}
b = alloc_d_n1(2*N);
G=alloc_d_n1_n2(ncells,6);
for(i=0;i<ncells;i++) G[i]=alloc_d_n1(6);

//-----assembler of constraints-----
  eps=sqrt(vol)*1e-9;
  for(i=0;i<4*NCN;i++) constr[i]=1.0/eps;
  I=0; for(i=0;i<N;i++) if(mask[i]==2) {Bind[I]=i; I++;}
  for(I=0;I<NCN;I++){ i=Bind[I];
      ind=0;
      //---boundary connectivity----
      for(l=0;l<ncells;l++){
	  nvert=0; while(cells[l][nvert]>=0) nvert++;
		  switch(nvert)
		  {case 3: //tri
		  if(i==cells[l][0]){
			  if(mask[cells[l][1]]>0) {if(ind==0) {j=cells[l][1]; ind++;} else {k=cells[l][1]; ind++;}}
			  if(mask[cells[l][2]]>0) {if(ind==0) {j=cells[l][2]; ind++;} else {k=cells[l][2]; ind++;}}
		  }
		  if(i==cells[l][1]){
			  if(mask[cells[l][0]]>0) {if(ind==0) {j=cells[l][0]; ind++;} else {k=cells[l][0]; ind++;}}
			  if(mask[cells[l][2]]>0) {if(ind==0) {j=cells[l][2]; ind++;} else {k=cells[l][2]; ind++;}}
		  }
		  if(i==cells[l][2]){
			  if(mask[cells[l][1]]>0) {if(ind==0) {j=cells[l][1]; ind++;} else {k=cells[l][1]; ind++;}}
			  if(mask[cells[l][0]]>0) {if(ind==0) {j=cells[l][0]; ind++;} else {k=cells[l][0]; ind++;}}
		  }
		  break;
		  case 4: //quad
		  if((i==cells[l][0])||(i==cells[l][3])){
			  if(mask[cells[l][1]]>0) {if(ind==0) {j=cells[l][1]; ind++;} else {k=cells[l][1]; ind++;}}
			  if(mask[cells[l][2]]>0) {if(ind==0) {j=cells[l][2]; ind++;} else {k=cells[l][2]; ind++;}}
		  }
		  if((i==cells[l][1])||(i==cells[l][2])){
			  if(mask[cells[l][0]]>0) {if(ind==0) {j=cells[l][0]; ind++;} else {k=cells[l][0]; ind++;}}
			  if(mask[cells[l][3]]>0) {if(ind==0) {j=cells[l][3]; ind++;} else {k=cells[l][3]; ind++;}}
		  }
		  break;
		  case 6: //quad tri
		  if(i==cells[l][0]){
			  if(mask[cells[l][1]]>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][4]; ind++;} else {k=cells[l][4]; ind++;}}