-- FreeFem++ v4.12 (Sat 06 May 2023 07:56:14 PM CEST - git v4.12-72-gb0015903) 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) << " " <