first.edp

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

verbosity=2;
mesh3 Th("dodecaedre01");
border cc(t=0,2*pi){x=cos(t);y=sin(t);label=1;}
mesh Th2=buildmesh(cc(50));
fespace Vh2(Th2,P2);
int nbtets=Th.nt;
cout << " Th mes " << Th.mesure << " border mes " << Th.bordermesure << endl;
cout << " nb of Tets = " << nbtets << endl;
if(1) {
  nbtets=2;
  for (int i=0;i<nbtets;i++)
    for (int j=0; j <4; j++)
      cout << i << " " << j << " Th[i][j] = "
	   << Th[i][j] << "  x = "<< Th[i][j].x  << " , y= "<< Th[i][j].y 
	   << ",  label=" << Th[i][j].label << endl;
	    
//   Th(i)   return   the vextex i of Th
//   Th[k]   return   the tet k of Th.

  // get vertices information : 
  int nbvertices=Th.nv;
  nbvertices=2;
  cout << " nb of vertices = " << nbvertices << endl;
  for (int i=0;i<nbvertices;i++)
	cout << "Th(" <<i  << ") : "   // << endl;	
	     << Th(i).x << " " << Th(i).y  << " " << Th(i).z << " " << Th(i).label // version 2.19 
	  << endl;
savemesh(Th,"dd.meshb");
 }
fespace Vh(Th,P23d);
Vh xx=x;
if(xx[].n == Th.nv)
  for(int i=0;i<Th.nv;++i)
    assert(abs(Th(i).x-xx[][i])<1e-6);

func ue =   2*x*x + 3*y*y + 4*z*z + 5*x*y+6*x*z+1;
func f= -18. ;
Vh u=f,b,d,uhe=ue,bc;
cout << " Vh.ndof =  " <<  Vh.ndof << endl;
cout << "  Vh.ndofK " << Vh.ndofK << endl;
cout << Th[5].region << endl;
// cout << Th(0,0,0).region << endl;  a faire ...
cout << Th[5][3] << endl;  // ok.. 


for(int i=0;i<Vh.ndofK;++i )
  cout << Vh(11,i) << " ";
 cout << endl;

cout << ue(0.1,0.2,0.3)<< "  == " << f(0.1,0.2,0.3) << endl; ;
macro Grad3(u) [dx(u),dy(u),dz(u)]  // EOM

varf vbc(u,v) =  on(0,u=1);
varf vlap(u,v) = int3d(Th)(Grad3(v)' *Grad3(u)) + int3d(Th)(f*v) + on(0,u=ue);
varf vBord(u,v,solver=CG) = int2d(Th)(u*v) ;
verbosity=10; 
matrix A=vlap(Vh,Vh);
matrix B=vBord(Vh,Vh);
verbosity=1; 
b[]=vlap(0,Vh);
//bc[]=vbc(0,Vh);
//cout << bc[] <<endl;
{
ofstream fa("A.txt");
ofstream fb("B.txt");
fa << A ;
fb << b[] ;
}


cout << b[]. min << " " << b[].max << endl;
u[]=A^-1*b[];
cout << u[]. min << " " << u[].max << endl;
real err= int3d(Th)( square(u-ue) );
d= ue-u;
cout <<  " err = " << err <<  " " << d[].linfty << endl;
cout << " u (0,0,0) "<< u(0.,0.,0.) << endl;
cout << " dx(u) (0,0,0) "<< dx(u)(0.,0.,0.) << endl;
cout << " dy u (0,0,0) "<< dy(u)(0.,0.,0.) << endl;
cout << " dz u (0,0,0) "<< dz(u)(0.,0.,0.) << endl;
Vh2 u2=u(x,y,0.);
plot(u2,wait=1);
plot(u2,wait=1);
	{ ofstream file("dd.bb"); 
	file << "3 1 1 "<< u[].n << " 2 \n";
	int j;
	for (j=0;j<u[].n ; j++)  
	  file << d[][j] << endl; 
    }