============================================== FreeFEM 4.14: examples/3d/test-suite.log ============================================== # TOTAL: 41 # PASS: 32 # SKIP: 0 # XFAIL: 0 # FAIL: 9 # XPASS: 0 # ERROR: 0 .. contents:: :depth: 2 FAIL: beam-3d.edp ================= -- FreeFem++ v4.14 (Thu Dec 21 17:55:20 EST 2023 - git v4.14) file : ./beam-3d.edp Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : load "medit" 2 : include "cube.idp"load "msh3" 2 : load "medit" (already loaded: medit) 3 : // ! basic functions to build regular mesh of a cube 4 : /* 5 : mesh3 Cube(NN,BB,L); 6 : -- build the surface mesh of a 3d box 7 : where: for exqmple: 8 : int[int] NN=[nx,ny,nz]; // the number of seg in the 3 direction 9 : real [int,int] BB=[[xmin,xmax],[ymin,ymax],[zmin,zmax]]; // bounding bax 10 : int [int,int] L=[[1,2],[3,4],[5,6]]; // the label of the 6 face left,right, front, back, down, right 11 : */ 12 : func mesh3 Cube(int[int] & NN,real[int,int] &BB ,int[int,int] & L) 13 : { 14 : // first build the 6 faces of the hex. 15 : real x0=BB(0,0),x1=BB(0,1); 16 : real y0=BB(1,0), *** Warning The identifier y0 hide a Global identifier y1=BB(1,1); *** Warning The identifier y1 hide a Global identifier 17 : real z0=BB(2,0),z1=BB(2,1); 18 : 19 : int nx=NN[0],ny=NN[1],nz=NN[2]; 20 : mesh Thx = square(nx,ny,[x0+(x1-x0)*x,y0+(y1-y0)*y]); 21 : 22 : int[int] rup=[0,L(2,1)], rdown=[0,L(2,0)], 23 : rmid=[1,L(1,0), 2,L(0,1), 3, L(1,1), 4, L(0,0) ]; 24 : mesh3 Th=buildlayers(Thx,nz, zbound=[z0,z1], 25 : labelmid=rmid, labelup = rup, 26 : labeldown = rdown); 27 : 28 : return Th; 29 : } 30 : func mesh3 Cube(int Nx,int Ny,int Nz) 31 : { 32 : int[int] NN=[Nx,Ny,Nz]; 33 : real [int,int] BB=[[0,1],[0,1],[0,1]]; 34 : int[int,int] LL=[[1,2],[3,4],[5,6]]; 35 : return Cube(NN,BB,LL); 36 : } 37 : 38 : 39 : 40 : 3 : int[int] Nxyz=[20,5,5]; 4 : real [int,int] Bxyz=[[0.,5.],[0.,1.],[0.,1.]]; 5 : int [int,int] Lxyz=[[1,2],[2,2],[2,2]]; 6 : mesh3 Th=Cube(Nxyz,Bxyz,Lxyz); 7 : 8 : real E = 21.5e4; 9 : real sigma = 0.29; 10 : real mu = E/(2*(1+sigma)); 11 : real lambda = E*sigma/((1+sigma)*(1-2*sigma)); 12 : real gravity = -0.05; 13 : 14 : fespace Vh(Th,[P1,P1,P1]); 15 : Vh [u1,u2,u3], [v1,v2,v3]; 16 : cout << "lambda,mu,gravity ="<4) cout << bb << endl; 24 : real RR = bb[3]; 25 : // cone using buildlayers with a triangle 26 : int MaxLayersT=(int(2*pi*RR/h)/4)*4; 27 : func zminT = 0; 28 : func zmaxT = 2*pi; 29 : func fx= y*cos(z);// / max( abs(cos(z) ), abs(sin(z))); 30 : func fy= y*sin(z);// / max( abs(cos(z) ), abs(sin(z))); 31 : func fz= x; 32 : mesh3 Th=buildlayers(Th2,coef= min(max(.01,y/RR*1.3),1.) , MaxLayersT,zbound=[zminT,zmaxT],transfo=[fx,fy,fz],facemerge=1); 33 : return Th; 34 : } 3 : load "medit" 4 : // cone using buildlayers with a triangle 5 : real RR=1,HH=1; 6 : border Taxe(t=0,HH){x=t;y=0;label=0;}; 7 : border Hypo(t=1,0){x=HH*t;y=RR*t;label=1;}; 8 : border Vert(t=0,RR){x=HH;y=t;label=2;}; 9 : 10 : int nn=10; 11 : real h= 1./nn; 12 : mesh Th2=buildmesh( Taxe(HH*nn)+ Hypo(sqrt(HH*HH+RR*RR)*nn) + Vert(RR*nn) ) ; 13 : 14 : mesh3 Th3T=BuildAxiOx(Th2,h); 15 : medit("cone",Th3T,wait=1); 16 : // FFCS: testing 3d plots 17 : plot(Th3T,cmm="cone"); 18 : sizestack + 1024 =2176 ( 1152 ) -- mesh: Nb of Triangles = 112, Nb of Vertices 74 FAIL: convect-3d.edp ==================== -- FreeFem++ v4.14 (Thu Dec 21 17:55:20 EST 2023 - git v4.14) file : ./convect-3d.edp Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : load "msh3" 2 : 3 : int nn=8; 4 : 5 : mesh Th2=square(nn,nn,[x*2.-1.,y*2.-1.]); 6 : fespace Vh2(Th2,P1); 7 : int[int] rup=[0,2], rdown=[0,1], rmid=[1,1,2,1,3,1,4,1]; 8 : real zmin=-1,zmax=1.; 9 : 10 : mesh3 Th=buildlayers(Th2,nn, 11 : zbound=[zmin,zmax], 12 : // region=r1, 13 : labelmid=rmid, 14 : reffaceup = rup, 15 : reffacelow = rdown); 16 : func real hill(real r2){return exp(-10.*(r2));}; 17 : 18 : fespace Vh(Th,P13d); 19 : 20 : macro Grad(u) [dx(u),dy(u),dz(u)] ) // EOM 21 : macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) ) //EOM 22 : 23 : Vh v,vo; 24 : Vh2 v2; 25 : real x0=0.3,y0=0.3, *** Warning The identifier y0 hide a Global identifier z0=0; 26 : vo=hill(square(x-x0)+square(y-y0)+square(z-z0)); 27 : 28 : real t=0; 29 : v2=vo(x,y,0); 30 : plot(v2,cmm=" cut y = 0.5, time ="+t,wait=1); 31 : real dt=0.1; 32 : func u1=1.; 33 : func u2=2.; 34 : func u3=3.; 35 : verbosity = 1; 36 : v=convect([u1,u2,u3],-dt,vo); 37 : verbosity = 1; 38 : v2=v(x,y,0); 39 : t += dt; 40 : plot(v2,cmm=" cut y = 0.5, time ="+t,wait=1); 41 : // verification ... 42 : int err=0; 43 : macro Verif(w,val) 44 # { 45 # real so= int3d(Th)(vo); 46 # real soi= int3d(Th)(vo*w); 47 # real sv= int3d(Th)(v); 48 # real svi= int3d(Th)(v*w); 49 # 50 # cout << Stringification(w) << " old " << soi/so << " new " << svi/sv << " delta " << (svi/sv - soi/so )/dt << " ~ " << val << endl; 51 # err += (abs((svi/sv - soi/so )/dt - val)> 0.2); 52 # } 53 # ) //EOF 54 : 55 : Verif(x,1) 44 @ 45 @ 46 @ 47 @ 48 @ 49 @ 50 @ 51 @ 52 @ 53 @ 44 @ { 45 @ real so= int3d(Th)(vo); 46 @ real soi= int3d(Th)(vo*x); 47 @ real sv= int3d(Th)(v); 48 @ real svi= int3d(Th)(v*x); 49 @ 50 @ cout << Stringification((xx)) << " old " << soi/so << " new " << svi/sv << " delta " << (svi/sv - soi/so )/dt << " ~ " << 1 << endl; 51 @ err += (abs((svi/sv - soi/so )/dt - 1)> 0.2); 52 @ } 53 @ 56 : Verif(y,2) 44 @ 45 @ 46 @ 47 @ 48 @ 49 @ 50 @ 51 @ 52 @ 53 @ 44 @ { 45 @ real so= int3d(Th)(vo); 46 @ real soi= int3d(Th)(vo*y); 47 @ real sv= int3d(Th)(v); 48 @ real svi= int3d(Th)(v*y); 49 @ 50 @ cout << Stringification((yy)) << " old " << soi/so << " new " << svi/sv << " delta " << (svi/sv - soi/so )/dt << " ~ " << 2 << endl; 51 @ err += (abs((svi/sv - soi/so )/dt - 2)> 0.2); 52 @ } 53 @ 57 : Verif(z,3) 44 @ 45 @ 46 @ 47 @ 48 @ 49 @ 50 @ 51 @ 52 @ 53 @ 44 @ { 45 @ real so= int3d(Th)(vo); 46 @ real soi= int3d(Th)(vo*z); 47 @ real sv= int3d(Th)(v); 48 @ real svi= int3d(Th)(v*z); 49 @ 50 @ cout << Stringification((zz)) << " old " << soi/so << " new " << svi/sv << " delta " << (svi/sv - soi/so )/dt << " ~ " << 3 << endl; 51 @ err += (abs((svi/sv - soi/so )/dt - 3)> 0.2); 52 @ } 53 @ 58 : assert(err==0); 59 : sizestack + 1024 =2216 ( 1192 ) -- Square mesh : nb vertices =81 , nb triangles = 128 , nb boundary edges 32 -- FESpace: Nb of Nodes 729 Nb of DoF 729 FAIL: intlevelset3d.edp ======================= -- FreeFem++ v4.14 (Thu Dec 21 17:55:20 EST 2023 - git v4.14) file : ./intlevelset3d.edp Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : load "medit" 2 : include "cube.idp"load "msh3" 2 : load "medit" (already loaded: medit) 3 : // ! basic functions to build regular mesh of a cube 4 : /* 5 : mesh3 Cube(NN,BB,L); 6 : -- build the surface mesh of a 3d box 7 : where: for exqmple: 8 : int[int] NN=[nx,ny,nz]; // the number of seg in the 3 direction 9 : real [int,int] BB=[[xmin,xmax],[ymin,ymax],[zmin,zmax]]; // bounding bax 10 : int [int,int] L=[[1,2],[3,4],[5,6]]; // the label of the 6 face left,right, front, back, down, right 11 : */ 12 : func mesh3 Cube(int[int] & NN,real[int,int] &BB ,int[int,int] & L) 13 : { 14 : // first build the 6 faces of the hex. 15 : real x0=BB(0,0),x1=BB(0,1); 16 : real y0=BB(1,0), *** Warning The identifier y0 hide a Global identifier y1=BB(1,1); *** Warning The identifier y1 hide a Global identifier 17 : real z0=BB(2,0),z1=BB(2,1); 18 : 19 : int nx=NN[0],ny=NN[1],nz=NN[2]; 20 : mesh Thx = square(nx,ny,[x0+(x1-x0)*x,y0+(y1-y0)*y]); 21 : 22 : int[int] rup=[0,L(2,1)], rdown=[0,L(2,0)], 23 : rmid=[1,L(1,0), 2,L(0,1), 3, L(1,1), 4, L(0,0) ]; 24 : mesh3 Th=buildlayers(Thx,nz, zbound=[z0,z1], 25 : labelmid=rmid, labelup = rup, 26 : labeldown = rdown); 27 : 28 : return Th; 29 : } 30 : func mesh3 Cube(int Nx,int Ny,int Nz) 31 : { 32 : int[int] NN=[Nx,Ny,Nz]; 33 : real [int,int] BB=[[0,1],[0,1],[0,1]]; 34 : int[int,int] LL=[[1,2],[3,4],[5,6]]; 35 : return Cube(NN,BB,LL); 36 : } 37 : 38 : 39 : 40 : 3 : real surfS1 = 4*pi; 4 : real volS1 =surfS1/3.; 5 : int nn= 16; 6 : int[int] Nxyz=[nn,nn,nn]; 7 : real [int,int] Bxyz=[[-2.,2.],[-2.,2.],[-2.,2.]]; 8 : int [int,int] Lxyz=[[1,1],[1,1],[1,1]]; 9 : mesh3 Th=Cube(Nxyz,Bxyz,Lxyz); 10 : 11 : int err=0; 12 : real eps = 0.5; 13 : func r = sqrt(x*x +y*y+z*z); 14 : 15 : real lc ; 16 : verbosity=3; 17 : lc = int2d(Th,levelset=r-1.)(1.) ; 18 : cout << " area of the level set = " << lc << " = surfS1 " << surfS1 ; 19 : cout << ", Ok = " << (abs(lc-surfS1) < eps) << endl; 20 : if( abs(lc-surfS1) > eps) err++; 21 : fespace Vh(Th,P1); 22 : // test linear and bilinear ... 23 : varf vl(u,v) = int2d(Th,levelset=r-1.)(v) + int2d(Th,levelset=r-1.)(u*v); 24 : real[int] vv=vl(0,Vh); 25 : 26 : cout << " area of the level set (varf linear ) = " << (lc=vv.sum) << "= surfS1 " << surfS1 ; 27 : cout << ", Ok = " << (abs(lc-surfS1) < eps) << endl; 28 : if( abs(lc-surfS1) > eps) err++; 29 : real[int] one(Vh.ndof); 30 : one=1.; 31 : matrix VV=vl(Vh,Vh); // matrix with levelset 32 : vv = VV*one; 33 : cout << " area of the level set (varf bilinear same) = " << (lc=vv.sum) << "= surfS1 " << surfS1; 34 : cout << ", Ok = " << (abs(lc-surfS1) < eps) << endl;; 35 : if( abs(lc-surfS1) > eps) err++; 36 : 37 : // just for test a idea approximation of int of negative part of levelset 38 : // to we just change the measure of the element not the quadrature point 39 : { // test new stuff for level set ... 40 : macro grad(u) [dx(u),dy(u),dz(u)] ) // 41 : Vh u,v; 42 : solve Pxx(u,v) = int3d(Th) ( grad(u) [dx(u),dy(u),dz(u)] '*grad(v) [dx(v),dy(v),dz(v)] *1e-8 ) + int3d(Th, levelset= 1-r) ( grad(u) [dx(u),dy(u),dz(u)] '*grad(v) [dx(v),dy(v),dz(v)] ) + on(1,u=0) + int3d(Th, levelset= 1-r) ( 1*v); 43 : plot(u,wait=1); 44 : varf vxx(u,v) = int3d(Th, levelset= 1-r) ( u*v ) + int3d(Th, levelset= 1-r) ( 1*v); 45 : matrix XX=vxx(Vh,Vh); 46 : real[int] xx=vxx(0,Vh); 47 : real vol1= int3d(Th, levelset= 1-r)(1.); 48 : cout << " vol1 = " << vol1 << " ~= " << Th.measure - volS1 << endl; 49 : err += (abs(vol1-(Th.measure - volS1)) > eps); 50 : cout << " xx.sum = " << xx.sum << " == " << vol1 < 1e-8); 52 : 53 : real[int] yy(Vh.ndof); yy=1; 54 : xx= XX*yy; 55 : cout << " XX.sum = " << xx.sum << " == " << vol1 << endl; 56 : err += (abs(vol1-xx.sum) > 1e-8); 57 : 58 : } 59 : 60 : 61 : 62 : if(0) 63 : {// test on diff mesh3 not wet implemented (FH frev 2014) 64 : mesh3 Th1=Cube(Nxyz,Bxyz,Lxyz); 65 : mesh3 Th2=Cube(Nxyz,Bxyz,Lxyz); 66 : fespace Vh1(Th1,P1); 67 : fespace Vh2(Th2,P1); 68 : 69 : varf vl(u,v) = int2d(Th,levelset=r-1.)(v) + int2d(Th,levelset=r-1.)(u*v); 70 : real[int] vv=vl(0,Vh2); 71 : 72 : cout << " area of the level set (varf linear diff ) = " << (lc=vv.sum) << "= surfS1 " << surfS1 ; 73 : cout << ", Ok = " << (abs(lc-surfS1) < eps) << endl; 74 : if( abs(lc-surfS1) > eps) err++; 75 : real[int] one(Vh1.ndof); 76 : one=1.; 77 : // sorry not implemented to day ... FH 78 : //verbosity=10000; 79 : matrix VV=vl(Vh1,Vh2); // no build of matrix with levelset 80 : vv = VV*one; 81 : cout << " area of the level set (varf bilinear diff ) = " << (lc=vv.sum) << "= surfS1 " << surfS1; 82 : cout << ", Ok = " << (abs(lc-surfS1) < eps) << endl;; 83 : if( abs(lc-surfS1) > eps) err++; 84 : 85 : } 86 : cout << " Nb err " << err << endl; 87 : assert(err==0); 88 : 89 : sizestack + 1024 =4376 ( 3352 ) -- Square mesh : nb vertices =289 , nb triangles = 512 , nb boundary edges 64 FAIL: LaplaceRT-3d.edp ====================== -- FreeFem++ v4.14 (Thu Dec 21 17:55:20 EST 2023 - git v4.14) file : ./LaplaceRT-3d.edp Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : /* 2 : $ - \Delta p = f $ on $\Omega$, 3 : $ dp / dn = (g1d,g2d). n $ on $\Gamma_{1}$ 4 : $ p = gd $ on $\Gamma_{2}$ 5 : with de Mixte finite element formulation 6 : 7 : Find $p\in L^2(\Omega) $ and $u\in H(div) $ such than 8 : $$ u - Grad p = 0 $$ 9 : $$ - div u = f $$ 10 : $$ u. n = (g1d,g2d). n \mbox{ on } \Gamma_{2}$$ 11 : $$ p = gd \mbox{ on }\Gamma_{1}$$ 12 : the variationnel form is: 13 : 14 : $\forall v\in H(div)$; $v.n = 0$ on $\Gamma_{2} $: 15 : 16 : $ \int_\Omega u v + p div v -\int_{\Gamma_{1}} gd* v.n = 0 $ 17 : $\forall q\in L^2$: $ +\int_\Omega q div u = -\int_Omega f q $ 18 : and $ u.n = (g1n,g2n).n$ on $\Gamma_2$ 19 : 20 : */ 21 : include "cube.idp"load "msh3" 2 : load "medit" 3 : // ! basic functions to build regular mesh of a cube 4 : /* 5 : mesh3 Cube(NN,BB,L); 6 : -- build the surface mesh of a 3d box 7 : where: for exqmple: 8 : int[int] NN=[nx,ny,nz]; // the number of seg in the 3 direction 9 : real [int,int] BB=[[xmin,xmax],[ymin,ymax],[zmin,zmax]]; // bounding bax 10 : int [int,int] L=[[1,2],[3,4],[5,6]]; // the label of the 6 face left,right, front, back, down, right 11 : */ 12 : func mesh3 Cube(int[int] & NN,real[int,int] &BB ,int[int,int] & L) 13 : { 14 : // first build the 6 faces of the hex. 15 : real x0=BB(0,0),x1=BB(0,1); 16 : real y0=BB(1,0), *** Warning The identifier y0 hide a Global identifier y1=BB(1,1); *** Warning The identifier y1 hide a Global identifier 17 : real z0=BB(2,0),z1=BB(2,1); 18 : 19 : int nx=NN[0],ny=NN[1],nz=NN[2]; 20 : mesh Thx = square(nx,ny,[x0+(x1-x0)*x,y0+(y1-y0)*y]); 21 : 22 : int[int] rup=[0,L(2,1)], rdown=[0,L(2,0)], 23 : rmid=[1,L(1,0), 2,L(0,1), 3, L(1,1), 4, L(0,0) ]; 24 : mesh3 Th=buildlayers(Thx,nz, zbound=[z0,z1], 25 : labelmid=rmid, labelup = rup, 26 : labeldown = rdown); 27 : 28 : return Th; 29 : } 30 : func mesh3 Cube(int Nx,int Ny,int Nz) 31 : { 32 : int[int] NN=[Nx,Ny,Nz]; 33 : real [int,int] BB=[[0,1],[0,1],[0,1]]; 34 : int[int,int] LL=[[1,2],[3,4],[5,6]]; 35 : return Cube(NN,BB,LL); 36 : } 37 : 38 : 39 : 40 : 22 : int[int] Nxyz=[10,10,10]; 23 : real [int,int] Bxyz=[[0,1],[0,1],[0,1]]; 24 : int [int,int] Lxyz=[[1,1],[1,1],[2,1]]; 25 : mesh3 Th=Cube(Nxyz,Bxyz,Lxyz); 26 : fespace Vh(Th,P1); 27 : fespace Rh(Th,RT03d); 28 : fespace Nh(Th,Edge03d);// Nedelec Finite element. 29 : fespace Ph(Th,P0); 30 : 31 : func gd = 1.; 32 : 33 : func g1n = 2.; 34 : func g2n = 3.; 35 : func g3n = 4.; 36 : 37 : func f = 1.; 38 : 39 : Rh [u1,u2,u3],[v1,v2,v3]; 40 : Nh [e1,e2,e3]; 41 : [u1,u2,u3]=[1+100*x,2+100*y,3+100*z]; 42 : 43 : // a + b ^ x = 44 : /* 45 : b1 x a1 + b2*z - b3*y 46 : b2 ^ y = a2 - b1*z + b3*x 47 : b3 z a3 + b1*y - b2*x 48 : */ 49 : real b1=30,b2=10,b3=20; 50 : func ex1=100+b2*z-b3*y; 51 : 52 : func ex1x=0.; 53 : func ex1y=-b3+0; 54 : func ex1z=b2+0; 55 : 56 : func ex2=200.- b1*z + b3*x ; 57 : func ex2x= b3 +0; 58 : func ex2y= 0. ; 59 : func ex2z= -b1 +0; 60 : func ex3=300.+b1*y - b2*x ; 61 : func ex3x= -b2 +0; 62 : func ex3y= b1 +0; 63 : func ex3z= 0. ; 64 : [e1,e2,e3]=[ex1,ex2,ex3]; 65 : 66 : int k=Th(.1,.2,.3).nuTriangle ; 67 : cout << " u = " << u1(.1,.2,.3) << " " << u2(.1,.2,.3) << " " << u3(.1,.2,.3) << endl; 68 : cout << " dx u = " << dx(u1)(.1,.2,.3) << " " << dy(u2)(.1,.2,.3) << " " << dz(u3)(.1,.2,.3) << endl; 69 : 70 : cout << " e = " << e1(.1,.2,.3) << " " << e2(.1,.2,.3) << " " << e3(.1,.2,.3) << endl; 71 : cout << " ex = " << ex1(.1,.2,.3) << " " << ex2(.1,.2,.3) << " " << ex3(.1,.2,.3) << endl; 72 : 73 : 74 : cout << " dx,dy,dz e1x= " << ex1x(.1,.2,.3) << " " << ex1y(.1,.2,.3) << " " << ex1z(.1,.2,.3) << endl; 75 : cout << " dx,dy,dz e2x= " << ex2x(.1,.2,.3) << " " << ex2y(.1,.2,.3) << " " << ex2z(.1,.2,.3) << endl; 76 : cout << " dx,dy,dz e3x= " << ex3x(.1,.2,.3) << " " << ex3y(.1,.2,.3) << " " << ex3z(.1,.2,.3) << endl; 77 : 78 : cout << " dx,dy,dz e1 = " << dx(e1)(.1,.2,.3) << " " << dy(e1)(.1,.2,.3) << " " << dz(e1)(.1,.2,.3) << endl; 79 : cout << " dx,dy,dz e2 = " << dx(e2)(.1,.2,.3) << " " << dy(e2)(.1,.2,.3) << " " << dz(e2)(.1,.2,.3) << endl; 80 : cout << " dx,dy,dz e3 = " << dx(e3)(.1,.2,.3) << " " << dy(e3)(.1,.2,.3) << " " << dz(e3)(.1,.2,.3) << endl; 81 : 82 : 83 : cout << " k = " << k << endl; 84 : cout << Rh(k,0) << " " < " << int(Th[k].adj((ee=e))) << " " << ee 37 : << " adj: " << ( Th[k].adj((ee=e)) != Th[k]) << endl; 38 : } 39 : // note : if k == int(Th[k].adj(ee=e)) not adjacent element 40 : 41 : 42 : int nbboundaryelement = Th.nbe; 43 : Th.be; 44 : for (int k=0;k 5 1 adj: 1 0 1 <=> -2 1 adj: 1 0 2 <=> 1 2 adj: 1 0 3 <=> -4 3 adj: 1 1 0 <=> -1 0 adj: 1 1 1 <=> 2 0 adj: 1 1 2 <=> 0 2 adj: 1 1 3 <=> -4 3 adj: 1 2 0 <=> 1 1 adj: 1 2 1 <=> -2 1 adj: 1 2 2 <=> 3 2 adj: 1 2 3 <=> -4 3 adj: 1 3 0 <=> -1 0 adj: 1 3 1 <=> 4 0 adj: 1 3 2 <=> 2 2 adj: 1 3 3 <=> -4 3 adj: 1 4 0 <=> 3 1 adj: 1 4 1 <=> -2 1 adj: 1 4 2 <=> 5 2 adj: 1 4 3 <=> -4 3 adj: 1 5 0 <=> -1 0 adj: 1 5 1 <=> 0 0 adj: 1 5 2 <=> 4 2 adj: 1 5 3 <=> -4 3 adj: 1 0 : 4 6 7 , label 6 tet 0 1 N -0 -0 1 1 : 4 0 6 , label 4 tet 0 3 N -1 0 -0 2 : 7 5 4 , label 6 tet 1 0 N -0 0 1 3 : 0 4 5 , label 1 tet 1 3 N -0 -1 -0 4 : 1 5 7 , label 2 tet 2 1 N 1 -0 -0 5 : 1 0 5 , label 1 tet 2 3 N -0 -1 0 6 : 7 3 1 , label 2 tet 3 0 N 1 -0 0 7 : 0 1 3 , label 5 tet 3 3 N -0 -0 -1 8 : 2 3 7 , label 3 tet 4 1 N -0 1 -0 9 : 2 0 3 , label 5 tet 4 3 N 0 -0 -1 10 : 7 6 2 , label 3 tet 5 0 N 0 1 -0 11 : 0 2 6 , label 4 tet 5 3 N -1 -0 -0 boundingbox xmin: 0 xmax: 0 ymin: 0 ymax: 0 solid angle = 12.5664 == 4*pi == 12.5664 -- Cube nv=64 nt=162 nbe=108 kind= 6 P 0 0.111111 0.222222 0 3 A = 0.5 0.5 0.5 FAIL: schwarz-nm-3d.edp ======================= -- FreeFem++ v4.14 (Thu Dec 21 17:55:20 EST 2023 - git v4.14) file : ./schwarz-nm-3d.edp Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : bool withmetis=1; 2 : bool RAS=0; 3 : int sizeoverlaps=2; // size off overlap 4 : int nnx=2,nny=2,nnz=2; 5 : 6 : func bool AddLayers(mesh3 & Th,real[int] &ssd,int n,real[int] &unssd) 7 : { 8 : // build a continuous function uussd (P1) : 9 : // ssd in the caracteristics function on the input sub domain. 10 : // such that : 11 : // unssd = 1 when ssd =1; 12 : // add n layer of element (size of the overlap) 13 : // and unssd = 0 ouside of this layer ... 14 : // --------------------------------- 15 : fespace Vh(Th,P1); 16 : fespace Ph(Th,P0); 17 : Ph s; 18 : assert(ssd.n==Ph.ndof); 19 : assert(unssd.n==Vh.ndof); 20 : unssd=0; 21 : s[]= ssd; 22 : // plot(s,wait=1,fill=1); 23 : Vh u; 24 : varf vM(u,v)=int3d(Th,qforder=1)(u*v/volume); 25 : matrix M=vM(Ph,Vh); 26 : 27 : for(int i=0;i.1; 32 : // plot(u,wait=1); 33 : unssd+= u[]; 34 : s[]= M'*u[];//'; 35 : s = s >0.1; 36 : } 37 : unssd /= (n); 38 : u[]=unssd; 39 : ssd=s[]; 40 : return true; 41 : } 42 : 43 : int withplot=3; 44 : include "cube.idp"load "msh3" 2 : load "medit" 3 : // ! basic functions to build regular mesh of a cube 4 : /* 5 : mesh3 Cube(NN,BB,L); 6 : -- build the surface mesh of a 3d box 7 : where: for exqmple: 8 : int[int] NN=[nx,ny,nz]; // the number of seg in the 3 direction 9 : real [int,int] BB=[[xmin,xmax],[ymin,ymax],[zmin,zmax]]; // bounding bax 10 : int [int,int] L=[[1,2],[3,4],[5,6]]; // the label of the 6 face left,right, front, back, down, right 11 : */ 12 : func mesh3 Cube(int[int] & NN,real[int,int] &BB ,int[int,int] & L) 13 : { 14 : // first build the 6 faces of the hex. 15 : real x0=BB(0,0),x1=BB(0,1); 16 : real y0=BB(1,0), *** Warning The identifier y0 hide a Global identifier y1=BB(1,1); *** Warning The identifier y1 hide a Global identifier 17 : real z0=BB(2,0),z1=BB(2,1); 18 : 19 : int nx=NN[0],ny=NN[1],nz=NN[2]; 20 : mesh Thx = square(nx,ny,[x0+(x1-x0)*x,y0+(y1-y0)*y]); 21 : 22 : int[int] rup=[0,L(2,1)], rdown=[0,L(2,0)], 23 : rmid=[1,L(1,0), 2,L(0,1), 3, L(1,1), 4, L(0,0) ]; 24 : mesh3 Th=buildlayers(Thx,nz, zbound=[z0,z1], 25 : labelmid=rmid, labelup = rup, 26 : labeldown = rdown); 27 : 28 : return Th; 29 : } 30 : func mesh3 Cube(int Nx,int Ny,int Nz) 31 : { 32 : int[int] NN=[Nx,Ny,Nz]; 33 : real [int,int] BB=[[0,1],[0,1],[0,1]]; 34 : int[int,int] LL=[[1,2],[3,4],[5,6]]; 35 : return Cube(NN,BB,LL); 36 : } 37 : 38 : 39 : 40 : 45 : int[int] NN=[25,25,25]; // the number of step in each direction 46 : real [int,int] BB=[[0,1],[0,1],[0,1]]; // bounding box 47 : int [int,int] L=[[1,1],[1,1],[1,1]]; // the label of the 6 face left,right, 48 : // front, back, down, right 49 : mesh3 Th=Cube(NN,BB,L); 50 : int npart= nnx*nny*nnz; 51 : fespace Ph(Th,P0); 52 : fespace Vh(Th,P1); 53 : 54 : Ph part; 55 : Vh sun=0,unssd=0; 56 : Ph xx=x,yy=y,zz=z,nupp; 57 : //part = int(xx*nnx)*nny + int(yy*nny); 58 : part = int(xx*nnx)*nny*nnz + int(yy*nny)*nnz + int(zz*nnz); 59 : //plot(part,wait=1); 60 : if(withmetis) 61 : { 62 : load "metis" load: init metis (v 5 ) ; 63 : int[int] nupart(Th.nt); 64 : metisdual(nupart,Th,npart); 65 : for(int i=0;i1) 69 : plot(part,fill=1,cmm="dual",wait=1); 70 : mesh3[int] aTh(npart); 71 : mesh3 Thi=Th; 72 : fespace Vhi(Thi,P1); 73 : Vhi[int] au(npart),pun(npart); 74 : matrix[int] Rih(npart); 75 : matrix[int] Dih(npart); 76 : matrix[int] aA(npart); 77 : Vhi[int] auntgv(npart); 78 : Vhi[int] rhsi(npart); 79 : 80 : for(int i=0;i0,label=10,split=1); 85 : Rih[i]=interpolate(Vhi,Vh,inside=1); // Vh -> Vhi 86 : if(RAS) 87 : { 88 : suppi= abs(part-i)<0.1; 89 : varf vSuppi(u,v)=int3d(Th,qforder=1)(suppi*v/volume); 90 : unssd[]= vSuppi(0,Vh); 91 : unssd = unssd>0.; 92 : if(withplot>19) 93 : plot(unssd,wait=1); 94 : } 95 : pun[i][]=Rih[i]*unssd[]; 96 : sun[] += Rih[i]'*pun[i][];//'; 97 : if(withplot>9) 98 : plot(part,aTh[i],fill=1,wait=1); 99 : } 100 : real[int] viso=[0,0.1,0.2,0.3]; 101 : plot(sun,wait=1,dim=3,fill=1,viso=viso); 102 : for(int i=0;i8) 107 : plot(pun[i],wait=1); 108 : } 109 : 110 : // verif partition of unite 111 : 112 : macro Grad(u) [dx(u),dy(u),dz(u)] ) //EOM 113 : sun=0; 114 : 115 : for(int i=0;i5) 139 : plot(sun,fill=1,wait=1); 140 : cout << sun[].max << " " << sun[].min<< endl; 141 : // verification of the partition of the unite. 142 : assert( 1.-1e-9 <= sun[].min && 1.+1e-9 >= sun[].max); 143 : 144 : int nitermax=1000; 145 : { 146 : Vh un=0; 147 : for(int iter=0;iter2) plot (cmm="Triangle",Triangle,wait=1); 147 : 148 : func f = 1; 149 : 150 : 151 : meshS Triangle3 = movemesh23(Triangle,transfo=[x,0,y]);//trianglesup 152 : if(wplot>2) plot (cmm="Triangle3",Triangle3,wait=1); 153 : 154 : meshS TriangleS = change(fregion=1,movemeshS(Triangle3,transfo=[x,sinico*y+cosico*z,-cosico*y+sinico*z]));//rotation de -(pi - diedre) par rapport à l'axe des x pour former une face de la pyramide pentagonale 155 : if(wplot>2) plot (cmm="TriangleS",TriangleS,wait=1); 156 : //medit("face pyramide pentagonale",TriangleS); 157 : 158 : meshS TriangleI = change(fregion=2,movemeshS(TriangleS,transfo=[x,-cosdiedre*y+sindiedre*z,-sindiedre*y-cosdiedre*z],orientation=-1));//triangle inf rotation de l'angle diedre par rapport au triangle sup 159 : if(wplot>2) plot (cmm="TriangleI",TriangleI,wait=1); 160 : 161 : meshS Triangles = TriangleI+TriangleS; 162 : 163 : if(wplot>1) plot (cmm="Triangles",Triangles,wait=1); 164 : //medit("face pyramide pentagonale + face antiprisme d'ordre 5",Triangles); 165 : 166 : meshS T1 = change(movemeshS(Triangles,transfo=[x-0.5,y-sinpi3*cosico,z],orientation=1),fregion = region);//translation pour que la figure soit sur le bord du pentagone 167 : meshS T2 = change(movemeshS(T1,transfo=[cos2pi5*x-sin2pi5*y,sin2pi5*x+cos2pi5*y,z],orientation=1),fregion = region+2);; 168 : meshS T3 = change(movemeshS(T2,transfo=[cos2pi5*x-sin2pi5*y,sin2pi5*x+cos2pi5*y,z]),fregion = region+2);; 169 : meshS T4 = change(movemeshS(T3,transfo=[cos2pi5*x-sin2pi5*y,sin2pi5*x+cos2pi5*y,z]),fregion = region+2);; 170 : meshS T5 = change(movemeshS(T4,transfo=[cos2pi5*x-sin2pi5*y,sin2pi5*x+cos2pi5*y,z]),fregion = region+2);; 171 : 172 : meshS Tdemi= T1+T2+T3+T4+T5;//moitié de l'icosaedre 173 : meshS Tdemi0 = movemeshS(Tdemi,transfo=[x,y,z+0.5*h]);//moitié supérieure 174 : meshS Tdemi1 = movemeshS(Tdemi0,transfo=[x,y,-z]);//moitié inférieure 175 : meshS Tdemi1rot = change(movemeshS(Tdemi1,transfo=[cospi5*x-sinpi5*y,sinpi5*x+cospi5*y,z],orientation=-1),fregion = region+10);//rotation de la moitié inférieure pour les emboiter 176 : meshS Ticosaedre = Tdemi0+Tdemi1rot; 177 : //Ticosaedre=trunc(Ticosaedre,1,split=5); 178 : if(wplot) plot(Ticosaedre,wait=1); 179 : //cout << regions(Ticosaedre) << endl; 180 : return Ticosaedre; 181 : } 182 : 183 : func meshS Icosahedron (int orientation) 184 : { 185 : return Icosahedron(orientation,0); 186 : } 187 : func meshS Sphere20 (real R,int N,int orientation,int wplot) *** Warning The identifier N hide a Global identifier 188 : {// Isocaedre regulier !!!! Thank G. Vergez .. 189 : meshS Ticosaedre = Icosahedron(orientation,wplot); 190 : Ticosaedre=trunc(Ticosaedre,1,split=N); 191 : if(wplot) plot(cmm="Icosaedre",Ticosaedre,wait=1); 192 : 193 : //================================================================================= 194 : //Construction de la sphere 3D 195 : //================================================================================= 196 : func metric =dist(x,y,z)/R; 197 : meshS Th = movemeshS(Ticosaedre,transfo=[x/metric,y/metric,z/metric]); 198 : if(wplot) plot (cmm="Th",Th,wait=1); 199 : return Th; 200 : } 201 : 202 : func meshS Sphere20 (real R,int N,int orientation) *** Warning The identifier N hide a Global identifier 203 : { 204 : return Sphere20(R,N,orientation,0); 205 : } 206 : 207 : /* test: 208 : load "tetgen" 209 : { 210 : real hs = 0.1; // mesh size on sphere 211 : int[int] N=[20,20,20]; 212 : real [int,int] B=[[-1,1],[-1,1],[-1,1]]; 213 : int [int,int] L=[[1,2],[3,4],[5,6]]; 214 : 215 : //////////////////////////////// 216 : meshS ThH = 217 : (N,B,L,1); 218 : meshS ThS =Sphere(0.5,hs,7,1); // "gluing" surface meshs to tolat boundary meshes 219 : cout << " xxxx" << ThH.nv << " " << ThS.nv << endl; 220 : 221 : meshS ThHS=ThH+ThS; 222 : savemesh(ThHS,"Hex-Sphere.mesh"); 223 : exec("ffmedit Hex-Sphere.mesh;rm Hex-Sphere.mesh"); 224 : 225 : real voltet=(hs^3)/6.; 226 : cout << " voltet = " << voltet << endl; 227 : real[int] domaine = [0,0,0,1,voltet,0,0,0.7,2,voltet]; 228 : 229 : mesh3 Th = tetg(ThHS,switch="pqaAAYYQ",nbofregions=2,regionlist=domaine); 230 : medit("Cube-With-Ball",Th); 231 : } 232 : 233 : */ 7 : 8 : // build mesh of a box (311) wit 2 holes (300,310) 9 : 10 : real hs = 0.8; 11 : int[int] N=[4/hs,8/hs,11.5/hs]; *** Warning The identifier N hide a Global identifier 12 : real [int,int] B=[[-2,2],[-2,6],[-10,1.5]]; 13 : int [int,int] L=[[311,311],[311,311],[311,311]]; 14 : meshS ThH = SurfaceHex(N,B,L,1); 15 : meshS ThSg = Sphere(1,hs,300,-1); // "gluing" surface meshs to tolat boundary meshes 16 : meshS ThSd = Sphere(1,hs,310,-1); 17 : 18 : plot (ThH,ThSg,ThSd); 19 : ThSd=movemesh(ThSd,[x,4+y,z]); 20 : 21 : meshS ThHS=ThH+ThSg+ThSd; 22 : medit("ThHS", ThHS); 23 : 24 : real voltet=(hs^3)/6.; 25 : cout << " voltet = " << voltet << endl; 26 : real[int] domaine = [0,0,-4,1,voltet]; 27 : real [int] holes=[0,0,0,0,4,0]; 28 : mesh3 Th = tetg(ThHS,switch="pqaAAYYQ",regionlist=domaine,holelist=holes); 29 : medit("Box-With-two-Ball",Th); 30 : // End build mesh 31 : 32 : // int[int] opt=[9,0,64,0,0,3]; // options of mmg3d see freeem++ doc 33 : 34 : real[int] vit=[0,0,-0.3]; 35 : func zero = 0.; 36 : func dep = vit[2]; 37 : 38 : fespace Vh(Th,P1); 39 : macro Grad(u) [dx(u),dy(u),dz(u)] ) // 40 : 41 : Vh uh,vh; // to compute the displacemnt field 42 : problem Lap(uh,vh,solver=CG) = int3d(Th)(Grad(uh) [dx(uh),dy(uh),dz(uh)] '*Grad(vh) [dx(vh),dy(vh),dz(vh)] ) //') for emacs 43 : + on(310,300,uh=dep) +on(311,uh=0.); 44 : 45 : for(int it=0; it<29; it++){ 46 : cout<<" ITERATION "<