poisson-cube-ballon.edp

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

verbosity=1;
load "msh3"
load "tetgen"
load "medit"

 
// 
mesh3 ThHex;
real volumetet;  // use in tetg.
{
	//  first  build the 6 faces of the hex.
real x0=-1,x1=1;
real y0=-1.1,y1=1.1;
real z0=-1.2,z1=1.2;

int nx=19,ny=20,nz=21;
//  a  volume  of  on tet. 
volumetet= (x1-x0)*(y1-y0)*(z1-z0)/ (nx*ny*ny) /6.;

mesh Thx = square(ny,nz,[y0+(y1-y0)*x,z0+(z1-z0)*y]);
mesh Thy = square(nx,nz,[x0+(x1-x0)*x,z0+(z1-z0)*y]);
mesh Thz = square(nx,ny,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

int[int] refz=[0,5];  //  bas
int[int] refZ=[0,6];   //  haut
int[int] refy=[0,3];  //  devant
int[int] refY=[0,4];   // derriere
int[int] refx=[0,1];  // gauche
int[int] refX=[0,2];   // droite 


mesh3 Thx0 = movemesh2D3Dsurf(Thx,transfo=[x0,x,y],normal=-1,refface=refx);
mesh3 Thx1 = movemesh2D3Dsurf(Thx,transfo=[x1,x,y],normal=1,refface=refX);
mesh3 Thy0 = movemesh2D3Dsurf(Thy,transfo=[x,y0,y],normal=+1,refface=refy);
mesh3 Thy1 = movemesh2D3Dsurf(Thy,transfo=[x,y1,y],normal=-1,refface=refY);
mesh3 Thz0 = movemesh2D3Dsurf(Thz,transfo=[x,y,z0],normal=-1,refface=refz);
mesh3 Thz1 = movemesh2D3Dsurf(Thz,transfo=[x,y,z1],normal=+1,refface=refZ);

//medit("  --- ", Thx0,Thx1,Thy0,Thy1,Thz0,Thz1);
 ThHex = Thx0+Thx1+Thy0+Thy1+Thz0+Thz1;
 
}
mesh3 Thsph; // 

{
mesh  Th=square(10,20,[x*pi-pi/2,2*y*pi]);  //  $]\frac{-pi}{2},frac{-pi}{2}[\times]0,2\pi[ $
//  a paratrization of a sphere 
func f1 =cos(x)*cos(y);
func f2 =cos(x)*sin(y);
func f3 = sin(x);
//  de  partiel derivatrive of the parametrization DF
func f1x=sin(x)*cos(y);   
func f1y=-cos(x)*sin(y);
func f2x=-sin(x)*sin(y);
func f2y=cos(x)*cos(y);
func f3x=cos(x);
func f3y=0;
// $  M = DF^t DF $
func m11=f1x^2+f2x^2+f3x^2;
func m21=f1x*f1y+f2x*f2y+f3x*f3y;
func m22=f1y^2+f2y^2+f3y^2;

func perio=[[4,y],[2,y],[1,x],[3,x]];  // to store the periodic condition 

// the intial mesh
savemesh(Th,"sphere",[f1,f2,f3]);

real R=0.5,hh=0.1/R;// hh  taille du maille sur la shere unite. 
real vv= 1/square(hh);
verbosity=2;
Th=adaptmesh(Th,m11*vv,m21*vv,m22*vv,IsMetric=1,inquire=1,periodic=perio);
plot(Th,wait=1);
Th=adaptmesh(Th,m11*vv,m21*vv,m22*vv,IsMetric=1,periodic=perio);
plot(Th,wait=1);
Th=adaptmesh(Th,m11*vv,m21*vv,m22*vv,IsMetric=1,periodic=perio);
plot(Th,wait=1);
Th=adaptmesh(Th,m11*vv,m21*vv,m22*vv,IsMetric=1,periodic=perio);

Thsph = movemesh2D3Dsurf(Th,transfo=[f1*R,f2*R,f3*R],normal=-1);
}


////////////////////////////////
mesh3 ThS = ThHex+Thsph; // "gluing" surface meshs to tolat boundary meshes
medit("Bounday mesh",ThS,wait=1);

// build a mesh of a axis parallel box with TetGen

real[int] domaine = [0,0,0,1,volumetet,0,0,0.7,2,volumetet];
mesh3 Th = tetg(ThS,switch="pqaAAYYQ",nbofregions=2,regionlist=domaine);    
// Tetrahelize the interior of the cube with tetgen
medit("tetg",Th);
savemesh(Th,"Th-hex-sph.mesh");

fespace Ph(Th,P03d);
fespace Vh(Th,P13d);
Ph reg=region;

cout << "  centre = " << reg(0,0,0) << endl;
cout << " exterieur = " << reg(0,0,0.7) << endl;

macro Grad(u) [dx(u),dy(u),dz(u)] // EOM

Vh uh,vh;
real f=1.;
real gn = 1.;
real cf= 1;
problem P(uh,vh)=
   int3d(Th,1)( Grad(uh)'*Grad(vh)*100) 
  +  int3d(Th,2)( Grad(uh)'*Grad(vh)*2) 
  + int3d(Th) (vh*f)
  + on(-1,uh=-1) + on(1,uh=1) 
  + int2d(Th,2,-2)(vh*gn)
  + int2d(Th,3,-3)(cf*vh*uh)
  ; 
  
  P;

plot(uh,wait=1, nbiso=6);
medit("   uh ",Th, uh);