/usr/share/doc/freefem++/examples/examples++-3d/Poisson-cube-ballon.edp is in freefem++-doc 3.19.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 | 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 = movemesh23(Thx,transfo=[x0,x,y],orientation=-1,label=refx);
mesh3 Thx1 = movemesh23(Thx,transfo=[x1,x,y],orientation=1,label=refX);
mesh3 Thy0 = movemesh23(Thy,transfo=[x,y0,y],orientation=+1,label=refy);
mesh3 Thy1 = movemesh23(Thy,transfo=[x,y1,y],orientation=-1,label=refY);
mesh3 Thz0 = movemesh23(Thz,transfo=[x,y,z0],orientation=-1,label=refz);
mesh3 Thz1 = movemesh23(Thz,transfo=[x,y,z1],orientation=+1,label=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 = movemesh23(Th,transfo=[f1*R,f2*R,f3*R],orientation=-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,wait=1);
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,wait=1);
|