metric.cpp

来自「FreeFem++可以生成高质量的有限元网格。可以用于流体力学」· C++ 代码 · 共 1,214 行 · 第 1/3 页

    {
      Mmass[iv]=0;
      OnBoundary[iv]=0;
      Mmassxx[iv]=0;
    }

  for (i=0;i<nbt;i++) 
    if(triangles[i].link) // the real triangles 
      {
	const Triangle &t=triangles[i];
	// coor of 3 vertices 
	R2 A=t[0];
	R2 B=t[1];
	R2 C=t[2];


	// number of the 3 vertices
	iA = Number(t[0]);
	iB = Number(t[1]);
	iC = Number(t[2]);
	
	Real8 dett = bamg::Area2(A,B,C);
	detT[i]=dett;
	dett /= 6;

	// construction of on boundary 
	int nbb =0;
	for(int j=0;j<3;j++)
          {
	    Triangle *ta=t.Adj(j);
	    if ( ! ta || !ta->link) // no adj triangle => edge on boundary
	      OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1,
		OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1,
		nbb++;
	  }
	
	workT[i] = nbb;
	Mmass[iA] += dett;
	Mmass[iB] += dett;
	Mmass[iC] += dett;
	
	if((nbb==0)|| !choice)
	  {
	    Mmassxx[iA] += dett;
	    Mmassxx[iB] += dett;
	    Mmassxx[iC] += dett;
	  }
      }
  else
    workT[i]=-1;

//  for (Int4 kcount=0;kcount<n;kcount++,ss++)
    for (Int4 nusol=0;nusol<nbsol;nusol++)
    { //for all Solution  

      Real8 smin=ss[0],smax=ss[0];
      
      Real8 h1=1.e30,h2=1e-30,rx=0;
      Real8 coef = 1./(anisomax*anisomax);
      Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30;  
      int nbfield = typsols? sizeoftype[typsols[nusol]] : 1; 
      if (nbfield == 1) 
       for ( iv=0,k=0; iv<nbv; iv++,k+=n )
				{
				  dxdx[iv]=dxdy[iv]=dydy[iv]=0;
				  smin=Min(smin,ss[k]);
				  smax=Max(smax,ss[k]);
				 }
			  else
			   {
         //  cas vectoriel 
          for ( iv=0,k=0; iv<nbv; iv++,k+=n )
          {	
           double v=0;		     
			     for (int i=0;i<nbfield;i++) 
			         v += ss[k+i]*ss[k+i];
			     v = sqrt(v);
				   smin=Min(smin,v);
				   smax=Max(smax,v);
			    }
			   }
      Real8 sdelta = smax-smin;
      Real8 absmax=Max(Abs(smin),Abs(smax));
      Real8 cnorm = DoNormalisation ? coef2/sdelta : coef2;
      
      if(verbosity>2) 
	     cout << "    Solution " << nusol <<  " Min = " << smin << " Max = " 
	       << smax << " Delta =" << sdelta << " cnorm = " << cnorm <<  " Nb of fields =" << nbfield << endl;

      
      if ( sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)) 
				{
				  if (verbosity>2)
				    cout << "      Solution " << nusol << " is constant. We skip. " 
					 << " Min = " << smin << " Max = " << smax << endl;
				continue;
				}
				
	 double *sf  = ss; 
	 for (Int4 nufield=0;nufield<nbfield;nufield++,ss++) 
	   {
	     for ( iv=0,k=0; iv<nbv; iv++,k+=n )
		       dxdx[iv]=dxdy[iv]=dydy[iv]=0;
       for (i=0;i<nbt;i++) 
	      if(triangles[i].link)
	  {// for real all triangles 
	    // coor of 3 vertices 
	    R2 A=triangles[i][0];
	    R2 B=triangles[i][1];
	    R2 C=triangles[i][2];
	    
	    
	    // warning the normal is internal and the 
	    //   size is the length of the edge
	    R2 nAB = Orthogonal(B-A);
	    R2 nBC = Orthogonal(C-B);
	    R2 nCA = Orthogonal(A-C);
	    // remark :  nAB + nBC + nCA == 0 

	    // number of the 3 vertices
	    iA = Number(triangles[i][0]);
	    iB = Number(triangles[i][1]);
	    iC = Number(triangles[i][2]);
	    
	    // for the test of  boundary edge
	    // the 3 adj triangles 
	    Triangle *tBC = triangles[i].TriangleAdj(OppositeEdge[0]);
	    Triangle *tCA = triangles[i].TriangleAdj(OppositeEdge[1]);
	    Triangle *tAB = triangles[i].TriangleAdj(OppositeEdge[2]);

	    // value of the P1 fonction on 3 vertices 
	    sA = ss[iA*n];
	    sB = ss[iB*n];
	    sC = ss[iC*n];

	    R2 Grads = (nAB * sC + nBC * sA + nCA * sB ) /detT[i] ;
	    if(choice) 
	      {
		int nbb = 0;
		Real8 dd = detT[i];
		Real8 lla,llb,llc,llf;
		Real8  taa[3][3],bb[3];
		// construction of the trans of lin system
		for (int j=0;j<3;j++)
		  {
		    int ie = OppositeEdge[j];
		    TriangleAdjacent ta = triangles[i].Adj(ie);
		    Triangle *tt = ta;
		    if (tt && tt->link)
		      {
			Vertex &v = *ta.OppositeVertex();
			R2 V = v;
			Int4 iV = Number(v);
			Real8 lA  = bamg::Area2(V,B,C)/dd;
			Real8 lB  = bamg::Area2(A,V,C)/dd;
			Real8 lC  = bamg::Area2(A,B,V)/dd;
			taa[0][j] =  lB*lC;
			taa[1][j] =  lC*lA;
			taa[2][j] =  lA*lB;
			//Real8 xx = V.x-V.y;
			//Real8 yy = V.x + V.y;
			//cout << " iv " << ss[iV*n] << " == " << (8*xx*xx+yy*yy)
			//     << " l = " << lA << " " << lB << " " << lC 
			//     << " = " << lA+lB+lC << " " <<  V << " == " << A*lA+B*lB+C*lC << endl;
			
			lla = lA,llb=lB,llc=lC,llf=ss[iV*n] ;

			bb[j]     =  ss[iV*n] - ( sA*lA + sB*lB + sC*lC ) ;
		      }
		    else
		      {
			nbb++;
			taa[0][j]=0;
			taa[1][j]=0;
			taa[2][j]=0;
			taa[j][j]=1;
			bb[j]=0;
		      }
		  }

		// resolution of 3x3 lineaire system transpose
		Real8 det33 =  det3x3(taa[0],taa[1],taa[2]);		
		Real8 cBC   =  det3x3(bb,taa[1],taa[2]);
		Real8 cCA   =  det3x3(taa[0],bb,taa[2]);
		Real8 cAB   =  det3x3(taa[0],taa[1],bb);
		
		assert(det33);
		//	det33=1;
		// verif
		//	cout << " " << (taa[0][0]*cBC +  taa[1][0]*cCA + taa[2][0] * cAB)/det33 << " == " << bb[0] ;
		//	cout << " " << (taa[0][1]*cBC +  taa[1][1]*cCA + taa[2][1] * cAB)/det33 << " == " << bb[1];
		//	cout << " " << (taa[0][2]*cBC +  taa[1][2]*cCA + taa[2][2] * cAB)/det33 << " == " << bb[2] 
		//	     << "  -- " ;
		//cout << lla*sA + llb*sB+llc*sC+ (lla*llb* cAB +  llb*llc* cBC + llc*lla*cCA)/det33 
		//   << " == " << llf <<  endl;
		// computation of the gradient in the element 
		
		// H( li*lj) = grad li grad lj + grad lj grad lj
		// grad li = njk  / detT ; with i j k ={A,B,C)
		Real8 Hxx = cAB * ( nBC.x*nCA.x) +  cBC * ( nCA.x*nAB.x) + cCA * (nAB.x*nBC.x);
		Real8 Hyy = cAB * ( nBC.y*nCA.y) +  cBC * ( nCA.y*nAB.y) + cCA * (nAB.y*nBC.y);
		Real8 Hxy = cAB * ( nBC.y*nCA.x) +  cBC * ( nCA.y*nAB.x) + cCA * (nAB.y*nBC.x) 
		          + cAB * ( nBC.x*nCA.y) +  cBC * ( nCA.x*nAB.y) + cCA * (nAB.x*nBC.y);
		Real8 coef = 1.0/(3*dd*det33);
		Real8 coef2 = 2*coef;
		//	cout << " H = " << Hxx << " " << Hyy << " " <<  Hxy/2 << " coef2 = " << coef2 << endl;
		Hxx *= coef2;
		Hyy *= coef2;
		Hxy *= coef2;
		//cout << i  << " H = " << 3*Hxx/dd << " " << 3*Hyy/dd << " " <<  3*Hxy/(dd*2) << " nbb = " << nbb << endl;
		if(nbb==0)
		  {
		    dxdx[iA] += Hxx;
		    dydy[iA] += Hyy;
		    dxdy[iA] += Hxy;
		    
		    dxdx[iB] += Hxx;
		    dydy[iB] += Hyy;
		    dxdy[iB] += Hxy;
		    
		    dxdx[iC] += Hxx;
		    dydy[iC] += Hyy;
		    dxdy[iC] += Hxy;
		  }
		
	      }
	    else
	      {
		
		// if edge on boundary no contribution  => normal = 0
		if ( ! tBC || ! tBC->link ) nBC = O;
		if ( ! tCA || ! tCA->link ) nCA = O;
		if ( ! tAB || ! tAB->link ) nAB = O;
	    
		// remark we forgot a 1/2 because
		//       $\\int_{edge} w_i = 1/2 $ if $i$ is in edge 
		//                          0  if not
		// if we don't take the  boundary 
		// dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x;
		
		dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x;
		dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x;
		dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x;
		
		// warning optimization (1) the divide by 2 is done on the metrix construction
		dxdy[iA] += (( nCA.y + nAB.y ) *Grads.x + ( nCA.x + nAB.x ) *Grads.y) ;
		dxdy[iB] += (( nAB.y + nBC.y ) *Grads.x + ( nAB.x + nBC.x ) *Grads.y) ;
		dxdy[iC] += (( nBC.y + nCA.y ) *Grads.x + ( nBC.x + nCA.x ) *Grads.y) ; 
		
		dydy[iA] += ( nCA.y + nAB.y ) *Grads.y;
		dydy[iB] += ( nAB.y + nBC.y ) *Grads.y;
		dydy[iC] += ( nBC.y + nCA.y ) *Grads.y;
	      }
	    
	  } // for real all triangles 
     Int4 kk=0;
      for ( iv=0,k=0 ; iv<nbv; iv++,k+=n )
	if(Mmassxx[iv]>0) 
	  {
	    dxdx[iv] /= 2*Mmassxx[iv];
	    // warning optimization (1) on term dxdy[iv]*ci/2 
	    dxdy[iv] /= 4*Mmassxx[iv];
	    dydy[iv] /= 2*Mmassxx[iv];
	    // Compute the matrix with abs(eigen value)
	    Metric M(dxdx[iv], dxdy[iv], dydy[iv]);
	    MatVVP2x2 Vp(M);
	    //cout <<iv <<  "  M  = " <<  M <<  " aniso= " << Vp.Aniso() ;
	    Vp.Abs();
	    M = Vp;
	      dxdx[iv] = M.a11;
	      dxdy[iv] = M.a21;
	      dydy[iv] = M.a22;
	      //  cout << " (abs)  iv M  = " <<  M <<  " aniso= " << Vp.Aniso() <<endl;
	  }
	else kk++;
      
      
      // correction of second derivate
      // by a laplacien

      Real8 *d2[3] = { dxdx, dxdy, dydy};
      Real8 *dd;
      for (int xy = 0;xy<3;xy++)
	{
	  dd = d2[xy];
      // do leat 2 iteration for boundary problem
	  for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++)
	    {
	      for (i=0;i<nbt;i++) 
		if(triangles[i].link) // the real triangles 
		  {
		    // number of the 3 vertices
		    iA = Number(triangles[i][0]);
		    iB = Number(triangles[i][1]);
		    iC = Number(triangles[i][2]);
		    Real8 cc=3;
		    if(ijacobi==0)
		      cc = Max((Real8) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.);
		    workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc;
		  }
	      for (iv=0;iv<nbv;iv++)
		workV[iv]=0;

	      for (i=0;i<nbt;i++) 
		if(triangles[i].link) // the real triangles 
		  {
		    // number of the 3 vertices
		    iA = Number(triangles[i][0]);
		    iB = Number(triangles[i][1]);
		    iC = Number(triangles[i][2]);
		    Real8 cc =  workT[i]*detT[i];
		    workV[iA] += cc;
		    workV[iB] += cc;
		    workV[iC] += cc;
		  }

	      for (iv=0;iv<nbv;iv++)
		if( ijacobi<NbJacobi || OnBoundary[iv])
		  dd[iv] = workV[iv]/(Mmass[iv]*6);
	      

	    }

	  
	}

      // constuction  of the metrix from the Hessian dxdx. dxdy,dydy

      Real8 rCutOff=CutOff*absmax;// relative cut off 

      for ( iv=0,k=0 ; iv<nbv; iv++,k+=n )
	{ // for all vertices 
	  //{
	  //Metric M(dxdx[iv], dxdy[iv], dydy[iv]);
	  // MatVVP2x2 Vp(M);	  
	  //cout << " iv M="<<  M << "  Vp = " << Vp << " aniso  " << Vp.Aniso() << endl;
	  //}
	  MetricIso Miso;
// new code to compute ci ---	  
	  Real8 ci ;
	  if (RelativeMetric)
	    { //   compute the norm of the solution
	       double xx =0,*sfk=sf+k; 
	       for (int ifield=0;ifield<nbfield;ifield++,sfk++)
	          xx += *sfk* *sfk;	       
	       xx=sqrt(xx);
	       ci = coef2/Max(xx,rCutOff);
	    }
	  else ci = cnorm;
	  
 // old 
//	  Real8 ci = RelativeMetric ? coef2/(Max(Abs(ss[k]),rCutOff)) : cnorm ;
 //   modif F Hecht 101099
	  Metric Miv(dxdx[iv]*ci, dxdy[iv]*ci,  dydy[iv]*ci);
	  MatVVP2x2 Vp(Miv);

	  Vp.Abs();
	 if(power!=1.0) 
	      Vp.pow(power);
	  


	  h1=Min(h1,Vp.lmin());
	  h2=Max(h2,Vp.lmax());

	  Vp.Maxh(hmin);
	  Vp.Minh(hmax);

	  rx = Max(rx,Vp.Aniso2());

	  Vp.BoundAniso2(coef);

	  hn1=Min(hn1,Vp.lmin());
	  hn2=Max(hn2,Vp.lmax());
	  rnx = Max(rnx,Vp.Aniso2());

	  Metric MVp(Vp);
	  vertices[iv].m.IntersectWith(MVp);
	}// for all vertices 
      if (verbosity>2)
	{ 
	  cout << "              Field " << nufield << " of solution " << nusol  << endl;
	  cout << "              before bounding :  Hmin = " << sqrt(1/h2) << " Hmax = " 
	       << sqrt(1/h1)  << " factor of anisotropy max  = " << sqrt(rx) << endl;
	  cout << "              after  bounding :  Hmin = " << sqrt(1/hn2) << " Hmax = " 
	       << sqrt(1/hn1)  << " factor of anisotropy max  = " << sqrt(rnx) << endl;
	}
	 } //  end of for all field
    }// end for all solution 

  delete [] detT;
  delete [] Mmass;
  delete [] dxdx;
  delete [] dxdy;
  delete [] dydy;
  delete []  workT;
  delete [] workV;
  delete [] Mmassxx;
  delete []  OnBoundary;
 
}


void Triangles::ReadMetric(const char * fmetrix,const Real8 hmin1=1.0e-30,const Real8 hmax1=1.0e30,const Real8 coef=1)
{