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
121
122
123
124
125
126
127
128
129
130
131
132
133
134
|
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;
// FFCS: with 3D view parameters
real[int] CameraPositionValue = [3.50634,-2.51489,2.60313];
real[int] CameraFocalPointValue = [0.0604689,-0.304636,-0.256484];
real[int] CameraViewUpValue = [0.7198,0.502367,-0.479078];
real[int] CutPlaneOriginValue = [-0.5,-0.55,0.0335184];
real[int] CutPlaneNormalValue = [0,0,1];
plot(uh,wait=1, nbiso=6,
BorderAsMesh = 1,
CameraPosition=CameraPositionValue,
CameraFocalPoint=CameraFocalPointValue,
CameraViewUp=CameraViewUpValue,
CutPlaneOrigin=CutPlaneOriginValue,
CutPlaneNormal = CutPlaneNormalValue);
medit(" uh ",Th, uh,wait=1);
|
