bamg.cpp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

	     << " file" << endl;
	exit(3);
      }
  // some verification
  NoMeshReconstruction= fmeshr !=0;
  if (!fmeshback) fmeshback=fmeshr;
  fileout = fmeshout || fam || fnopo || fftq || fam || famdba ||famfmt || wbb || wBB  ;
  if (!fileout && !foM )
    {
      cerr << " No Output file a given " << endl;
      MeshError(1);
    }
  if ( maxsubdiv > boundmaxsubdiv || maxsubdiv <= 1.0)
    {
      cerr << " -maxsubdiv " << maxsubdiv << " is not in ]1,"<< boundmaxsubdiv << "]" << endl;
      exit(3);
    }
  if (iso) 
    anisomax=1;
  if (!(fmetrix||fMbb))
    NbSmooth=0; // no metric -> no smoothing 
  if( ! rbb) // to set the rbb filename by default 
    rbb=fMbb;
  if( ! rBB) // to set the rbb filename by default 
    rBB=fMBB;

  if (fMbb && rbb) 
    rbbeqMbb = !strcmp(rbb,fMbb);
  if (fMBB && rBB) 
    rBBeqMBB = !strcmp(rBB,fMBB);

#ifdef DRAWING 
  if(initgraphique)
    {
      initgraphique();
      initgraph=1;
      withrgraphique=0; 
    }
#endif
  
     if (verbosity) 
       cout << bamgversion <<" : ";
     if (verbosity)
       if (fgeom && fileout)
	 cout << " Construction of a mesh from the geometry : " << fgeom << endl
	      << "     with less than " << nbvx << " vertices " << endl;
       else if(fmeshback && fileout)
	 if (NoMeshReconstruction)
	   cout <<  " Modification of a adpated mesh " 
		<< fmeshback 
		<< "    with less than " << nbvx << " vertices" << endl;
	 else
	   cout <<  " Construction of a adpated mesh from the background mesh: " 
		<< fmeshback 
		<< "    with less than " << nbvx << " vertices" <<  endl;
       else if (fmeshback && foM)
	 cout << " Construction of the metric file on the background mesh: " 
	      << fmeshback << endl;
       

     if (fgeom && fileout)
       {
	 if (verbosity) 
	   cout << " Construction of a mesh from the geometry : " << fgeom << endl
		<< "     with less than " << nbvx << " vertices " << endl;

	 time0 = CPUtime();
	 Geometry Gh(fgeom);
	 hmin = Max(hmin,Gh.MinimalHmin());
	 hmax = Min(hmax,Gh.MaximalHmax());
	 if(fmetrix) 
	   Gh.ReadMetric(fmetrix,hmin,hmax,coef);
	 else
	   {
	     Int4 iv;
	     for(iv=0;iv<Gh.nbv;iv++)
	       {
		 MetricAnIso M=Gh[iv];
		 MatVVP2x2 Vp(M/coef);
		 Vp.Maxh(hmin);
		 Vp.Minh(hmax);
		 Gh.vertices[iv].m = Vp;
	       }
	   }
	 
	 time1 = CPUtime();
	 if(verbosity>3)
	   cout << " Cpu for read the geometry " << time1-time0 << "s" << endl;
	 Triangles Th(nbvx,Gh);
	 if(costheta<=1)
	   Th.MakeQuadrangles(costheta);
	 if (allquad)
	      Th.SplitElement(allquad==2);
	   
	 if(SplitEdgeWith2Boundary)
	   Th.SplitInternalEdgeWithBorderVertices();
	 Th.ReNumberingTheTriangleBySubDomain();
	 time2 = CPUtime();
	 if(verbosity>3) Th.ShowHistogram();
	  if(NbSmooth>0) Th.SmoothingVertex(NbSmooth,omega);
	  if(verbosity>3 && NbSmooth>0) Th.ShowHistogram();
#ifdef DRAWING
	  if(initgraph)
	    {
	      Th.InitDraw();
	      Th.Draw();
	      Th.inquire();
	    }
#endif
	 time3 = CPUtime();
	 if(verbosity> 1)
	   cout << " Cpu for meshing         :" << setw(8) << time2-time1 
		<< "s, for  Smoothing " <<  time3-time2 
		<< "s  Nb Vertices/s = " <<  (Th.nbv) /(time2-time1)
		<< setw(0) <<endl ;
	 if (fmeshout) Th.Write(fmeshout  ,Triangles::BDMesh);
	 if (famfmt)   Th.Write(famfmt    ,Triangles::am_fmtMesh);
	 if (fam)      Th.Write(fam       ,Triangles::amMesh);
	 if (famdba)   Th.Write(famdba    ,Triangles::amdbaMesh);
	 if (fftq)     Th.Write(fftq      ,Triangles::ftqMesh);
	 if (fmsh)     Th.Write(fmsh      ,Triangles::mshMesh);
	 if (fnopo)    Th.Write(fnopo     ,Triangles::NOPOMesh);
#ifdef DRAWING 
	  if(initgraph)
	    {
	      
	      reffecran();
	      Th.InitDraw();
	      Th.Draw();
	      Th.inquire();
	    }
#endif

	 time3 = CPUtime();
	 if(verbosity>0) 
	       {
		 cout << " Cpu for meshing with io :" << setw(8) <<time3-time0 
		      << "s  Nb Vertices/s = " <<  (Th.nbv) /(time3-time0)	 << setw(0)<< endl;
		 cout << " Nb vertices = " << Th.nbv;
	       if (Th.nbt-Th.NbOutT-Th.NbOfQuad*2) 
		 cout  << " Nb Triangles = " << (Th.nbt-Th.NbOutT-Th.NbOfQuad*2);
	       if (Th.NbOfQuad ) 
		 cout  << " Nb Quadrilaterals = " << Th.NbOfQuad  ;
	       cout   << endl;
	       }	   

       }
     else if (fmeshback && (fmetrix||fMbb||fMBB||NoMeshReconstruction) && 
	      (fileout || foM  || ( rbb && wbb)  ||( rBB && wBB) ) ) 
       {
	 // to be sure the value was not stupide
	 
	 time0 = CPUtime();
	 Triangles BTh(fmeshback,cutoffradian); // read the background mesh 
	 hmin = Max(hmin,BTh.MinimalHmin());
	 hmax = Min(hmax,BTh.MaximalHmax());

	 BTh.MakeQuadTree();
	 time1 = CPUtime();
	 if (fmetrix) 
	   BTh.ReadMetric(fmetrix,hmin,hmax,coef);
	 else
	   { // init with hmax 
	     Metric Mhmax(hmax);
	     for (Int4 iv=0;iv<BTh.nbv;iv++)
		 BTh[iv].m = Mhmax;
	   }
	 if (fMbb)
	   {
	     solMbb = ReadbbFile(fMbb,nbsolbb,lsolbb,2,2);
	     if (lsolbb != BTh.nbv) 
	       {
	         cerr << "fatal error  nbsol " << nbsolbb << " " << lsolbb<< " =! " << BTh.nbv << endl;
		 cerr << "  size of sol incorrect " << endl;
		 MeshError(99);
	       }
	     assert(lsolbb==BTh.nbv);
	     BTh.IntersectConsMetric(solMbb,nbsolbb,0,hmin,hmax,sqrt(err)*coef,1e30,AbsError?0.0:CutOff,
				     nbjacoby,Rescaling,power,ChoiseHessien);
	    if(!rbbeqMbb)
	      delete [] solMbb,solMbb =0;

	   }
	 if (fMBB)
	   {
	     solMBB = ReadBBFile(fMBB,nbsolBB,lsolBB,typesolsBB,2,2);
	     if (lsolBB != BTh.nbv) 
	       {
	         cerr << "fatal error  nbsol " << nbsolBB << " " << lsolBB << " =! " << BTh.nbv << endl;
		 cerr << "  size of sol incorrect " << endl;
		 MeshError(99);
	       }
	     assert(lsolBB==BTh.nbv);
	     BTh.IntersectConsMetric(solMBB,nbsolBB,0,hmin,hmax,sqrt(err)*coef,1e30,AbsError?0.0:CutOff,
				     nbjacoby,Rescaling,ChoiseHessien);
	    if(!rBBeqMBB)
	      delete [] solMBB,solMBB =0;

	   }
	 
	 BTh.IntersectGeomMetric(errg,iso);
	 if(raison) 
	   BTh.SmoothMetric(raison);
	 BTh.MaxSubDivision(maxsubdiv);
	 //
	 if (iso) 
	   anisomax=1;
	 BTh.BoundAnisotropy(anisomax,hminaniso);
	 time2 = CPUtime();
	 if(foM)
	   {
	     if(verbosity >2)
	       cout << " -- write Metric  file " << foM <<endl;
		 
	     ofstream f(foM);
	     if(f) BTh.WriteMetric(f,iso);
	   }

	 if ( fileout) 
	   {
	     if ( NoMeshReconstruction)
               if (( fmeshback == fmeshr) || (!strcmp(fmeshback,fmeshr)))
		 Thr=&BTh,Thb=0; // back and r mesh are the same 
	       else
                 Thr= new Triangles(fmeshr,cutoffradian),Thb=&BTh; // read the new 

		
	     
	     Triangles & Th( *(NoMeshReconstruction 
			     ?  new Triangles(*Thr,&Thr->Gh,Thb,nbvx) // copy the mesh + free space to modification  
			     :  new Triangles(nbvx,BTh,KeepBackVertices)     // construct a new mesh 
			     ));
             if (Thr != &BTh) delete Thr;
             
	     if(costheta<=1.0)
	       Th.MakeQuadrangles(costheta);
	     if (allquad)
	       Th.SplitElement(allquad==2);
	     if(SplitEdgeWith2Boundary)
	       Th.SplitInternalEdgeWithBorderVertices();
	     Th.ReNumberingTheTriangleBySubDomain();
	     time3 = CPUtime();
	     
	     if(verbosity>0) 
	       cout << " Cpu time for meshing          : " << time3-time2 
		    << "s  Nb Vertices/s = " <<  (Th.nbv) /(time3-time2)
		    << endl;
	     if(verbosity>3) Th.ShowHistogram();
	     if(NbSmooth>0) Th.SmoothingVertex(NbSmooth,omega);
	     if(verbosity>3 && NbSmooth>0) Th.ShowHistogram();
#ifdef DRAWING
	  if(initgraph)
	    {
	      Th.InitDraw();
	      Th.Draw();
	      Th.inquire();
	    }
#endif
	 Th.BTh.ReMakeTriangleContainingTheVertex();

	 if (fmeshout) Th.Write(fmeshout  ,Triangles::BDMesh);
	 if (famfmt)   Th.Write(famfmt    ,Triangles::am_fmtMesh);
	 if (fam)      Th.Write(fam       ,Triangles::amMesh);
	 if (famdba)   Th.Write(famdba    ,Triangles::amdbaMesh);
	 if (fftq)     Th.Write(fftq      ,Triangles::ftqMesh);
	 if (fmsh)     Th.Write(fmsh      ,Triangles::mshMesh);
	 if (fnopo)    Th.Write(fnopo     ,Triangles::NOPOMesh);
	 
	     if ( ( rbb && wbb)  ||( rBB && wBB))  // the code for interpolation 
	       {
		 if(verbosity >1)
		   {
		     if (rbb ) 
		       cout << " -- interpolation P1  files : " << rbb << " -> " << wbb <<endl;
		     if (rBB ) 
		       cout << " -- interpolation P1  files : " << rBB<< " -> " << wBB <<endl;
		   }
		 double time00  = CPUtime();	   
		 const int dim=2;
		 // optimisation read only si rbb != fMbb
		
		 double *solbb=0;
		 double *solBB=0;

		 if (rbb) 
		   solbb =  rbbeqMbb? solMbb : ReadbbFile(rbb,nbsolbb,lsolbb,2,2);
		 if (rBB) 
		   solBB =  rBBeqMBB? solMBB : ReadBBFile(rBB,nbsolBB,lsolBB,typesolsBB,2,2);

		 // cout << " " << rbbeqMbb << " " <<  sol << " " << nbsol << " " << lsol << endl; 
		 if (!solBB && !solbb ) 
		   {
		     cerr << " Fatal Error "  << "solBB =  " <<  solBB << " solbb= " << solbb << endl;
		     exit(2);}
		 
		 ofstream *fbb = wbb ? new ofstream(wbb) :0;
		 ofstream *fBB = wBB ? new ofstream(wBB) :0;
		 Int4   nbfieldBB = 0, nbfieldbb = nbsolbb;
		 if (fbb)
		   *fbb  << dim << " " << nbsolbb << " " << Th.nbv <<" " << 2 << endl; 
		 if (fBB)
		   {
		     int i;
		     *fBB  << dim << " " << nbsolBB ;
		     for ( i=0;i<nbsolBB;i++)
		       *fBB << " " << (typesolsBB ?typesolsBB[i]+1:1) ;
		     *fBB << " " << Th.nbv <<" " << 2 << endl; 
		     assert(dim==2);
		     for ( i=0;i<nbsolBB;i++)
		     nbfieldBB += typesolsBB ? typesolsBB[i]+1 : 1;
		     // this code is good only if dim == 2 
		   }
		 cout << "nb of field BB " << nbfieldBB << endl;
		 //		 if(fBB) fBB->precision(15);
		 //if(fbb) fbb->precision(15);
		 for(i=0;i<Th.nbv;i++)
		   {
		     Int4 i0,i1,i2;
		     double a[3];
		     Icoor2 dete[3];
		     I2 I = Th.BTh.toI2(Th.vertices[i].r);
		     Triangle & tb = *Th.BTh.FindTriangleContening(I,dete);
		     
		     if (tb.det>0) { // internal point 
		       a[0]= (Real8) dete[0]/ tb.det;
		       a[1]= (Real8) dete[1] / tb.det;
		       a[2] = (Real8) dete[2] / tb.det;
		       i0=Th.BTh.Number(tb[0]);
		       i1=Th.BTh.Number(tb[1]);
		       i2=Th.BTh.Number(tb[2]);}
		     else {
		       double aa,bb;
		       
		       TriangleAdjacent  ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj();
		       
		       int k = ta;
		       Triangle & tc = *(Triangle *)ta;
		       i0=Th.BTh.Number(tc[0]);
		       i1=Th.BTh.Number(tc[1]);
		       i2=Th.BTh.Number(tc[2]);
		       a[VerticesOfTriangularEdge[k][1]] =aa;
		       a[VerticesOfTriangularEdge[k][0]] = bb;
		       a[OppositeVertex[k]] = 1- aa -bb;}
		     
		     Int4  ibb0 = nbfieldbb*i0;
		     Int4  ibb1 = nbfieldbb*i1;
		     Int4  ibb2 = nbfieldbb*i2;

		     Int4  iBB0 = nbfieldBB*i0;
		     Int4  iBB1 = nbfieldBB*i1;
		     Int4  iBB2 = nbfieldBB*i2;
		     Int4 j=0;

		     for ( j=0;j<nbfieldbb;j++) 
		       *fbb  << " " << ( a[0] * solbb[ibb0++] + a[1] * solbb[ibb1++] + a[2]* solbb[ibb2++]) ;

		     for (j=0;j<nbfieldBB;j++) 
		       *fBB  << " " << ( a[0] * solBB[iBB0++] + a[1] * solBB[iBB1++] + a[2]* solBB[iBB2++]) ;
		     
		     if (fbb) *fbb  << endl;
		     if (fBB) *fBB  << endl;
		   }
		 if (fbb)
		   delete fbb  ; // close
		 if (fBB)
		   delete fBB  ; // close
		 if (solbb) 
		   delete [] solbb;
		 if (solBB) 
		   delete [] solBB;
		 double  time04 = CPUtime();
		 if(verbosity>0) 
		   cout << " Cpu time for interpolation " << time04-time00 << " s" <<endl;
		 
	       }
	     time4 = CPUtime();
	     if(verbosity>0) 
	       {
	       cout << " Cpu time for meshing with io  : " << time4-time0
		    << "s  Nb Vertices/s = " << Th.nbv/(time4-time0)
		    << endl
		    << " Nb vertices = " << Th.nbv;
	       if (Th.nbt-Th.NbOutT-Th.NbOfQuad*2) 
		 cout  << " Nb Triangles = " << (Th.nbt-Th.NbOutT-Th.NbOfQuad*2);
	       if (Th.NbOfQuad ) 
		   cout  << " Nb Quadrilaterals = " << Th.NbOfQuad  ;
		     cout   << endl;
	       }
	     
#ifdef DRAWING
	     if(initgraph)
	       {
		 reffecran();
		 Th.InitDraw();
		 Th.Draw();
		 Th.inquire();
	       }
#endif
	     
	     // cout << "delete &Th " ;
	     delete &Th;
	     //cout << " end " << endl;
	   }
       }// if (fgeom && fmeshout)
     // clean the 
     for(i=1;i<datargc;i++)
       delete [] datargv[i] ;

#ifdef DRAWING 
     if(initgraph)
       closegraphique();
     initgraph=0;
#endif
     cout << flush;
     return(0);
     }