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 $ 18and $ u.n = (g1n,g2n).n$ on $\Gamma_2$ 19 20*/ 21include "cube.idp" 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]]; 25mesh3 Th=Cube(Nxyz,Bxyz,Lxyz); 26fespace Vh(Th,P1); 27fespace Rh(Th,RT03d); 28fespace Nh(Th,Edge03d);// Nedelec Finite element. 29fespace Ph(Th,P0); 30 31func gd = 1.; 32 33func g1n = 2.; 34func g2n = 3.; 35func g3n = 4.; 36 37func f = 1.; 38 39Rh [u1,u2,u3],[v1,v2,v3]; 40Nh [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*/ 49real b1=30,b2=10,b3=20; 50func ex1=100+b2*z-b3*y; 51 52func ex1x=0.; 53func ex1y=-b3+0; 54func ex1z=b2+0; 55 56func ex2=200.- b1*z + b3*x ; 57func ex2x= b3 +0; 58func ex2y= 0. ; 59func ex2z= -b1 +0; 60func ex3=300.+b1*y - b2*x ; 61func ex3x= -b2 +0; 62func ex3y= b1 +0; 63func ex3z= 0. ; 64[e1,e2,e3]=[ex1,ex2,ex3]; 65 66int k=Th(.1,.2,.3).nuTriangle ; 67cout << " u = " << u1(.1,.2,.3) << " " << u2(.1,.2,.3) << " " << u3(.1,.2,.3) << endl; 68cout << " dx u = " << dx(u1)(.1,.2,.3) << " " << dy(u2)(.1,.2,.3) << " " << dz(u3)(.1,.2,.3) << endl; 69 70cout << " e = " << e1(.1,.2,.3) << " " << e2(.1,.2,.3) << " " << e3(.1,.2,.3) << endl; 71cout << " ex = " << ex1(.1,.2,.3) << " " << ex2(.1,.2,.3) << " " << ex3(.1,.2,.3) << endl; 72 73 74cout << " dx,dy,dz e1x= " << ex1x(.1,.2,.3) << " " << ex1y(.1,.2,.3) << " " << ex1z(.1,.2,.3) << endl; 75cout << " dx,dy,dz e2x= " << ex2x(.1,.2,.3) << " " << ex2y(.1,.2,.3) << " " << ex2z(.1,.2,.3) << endl; 76cout << " dx,dy,dz e3x= " << ex3x(.1,.2,.3) << " " << ex3y(.1,.2,.3) << " " << ex3z(.1,.2,.3) << endl; 77 78cout << " dx,dy,dz e1 = " << dx(e1)(.1,.2,.3) << " " << dy(e1)(.1,.2,.3) << " " << dz(e1)(.1,.2,.3) << endl; 79cout << " dx,dy,dz e2 = " << dx(e2)(.1,.2,.3) << " " << dy(e2)(.1,.2,.3) << " " << dz(e2)(.1,.2,.3) << endl; 80cout << " dx,dy,dz e3 = " << dx(e3)(.1,.2,.3) << " " << dy(e3)(.1,.2,.3) << " " << dz(e3)(.1,.2,.3) << endl; 81 82 83cout << " k = " << k << endl; 84cout << Rh(k,0) << " " <<Rh(k,1) << " " <<Rh(k,2) << " " <<Rh(k,3) << endl; 85cout << " df = " << u1[][Rh(k,0)] << " " << u1[][Rh(k,1)] <<" " << u1[][Rh(k,2)] << " " << u1[][Rh(k,2)] << endl; 86// cout << u1[] << endl; 87 88Vh P,Q; 89Ph p,q; 90macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) // 91macro Grad(u) [dx(u),dy(u),dz(u)] // 92 problem laplace(P,Q,solver=CG) = 93 int3d(Th) ( Grad(P)'*Grad(Q)) //') for emacs 94 - int3d(Th)(f*Q) 95 + on(1,P=gd) 96 - int2d(Th,2) ( (g1n*N.x+g2n*N.y+g3n*N.z)*Q); 97 98fespace RPh(Th,[RT03d,P0]); 99varf von1([u1,u2,u3,p],[v1,v2,v3,q]) = 100 int3d(Th)( p*q*1e-15+ u1*v1 + u2*v2 + u3*v3 + p*div(v1,v2,v3) + div(u1,u2,u3)*q ) 101 - int3d(Th) ( f*q) 102 + int2d(Th,1)( gd*(v1*N.x +v2*N.y + v3*N.z) ) // int on gamma 103 + on(2,u1=g1n,u2=g2n,u3=g3n); 104 105RPh [vv1,vv2,vv3,qq]; 106// some verification Boundary Condition 107// and interpolation ... 108real[int] ron=von1(0,RPh); 109vv3[]=von1(0,RPh); 110cout << " vv: = " << vv1(.1,.2,.001) << " " << vv2(.1,.2,.001) << " " << vv3(.1,.2,.001) << endl; 111[vv1,vv2,vv3,qq]=[g1n,g2n,g3n,100]; 112[v1,v2,v3]=[g1n,g2n,g3n]; 113 114cout << " vv: = " << vv1(.1,.2,.001) << " " << vv2(.1,.2,.001) << " " << vv3(.1,.2,.001) << " " << qq(.1,.2,.001) << endl; 115cout << " v : = " << v1(.1,.2,.001) << " " << v2(.1,.2,.001) << " " << v3(.1,.2,.001) << endl; 116 117// end of verification of Boundary Condition ... 118 119problem laplaceMixte([u1,u2,u3,p],[v1,v2,v3,q],eps=1.0e-10,tgv=1e30,dimKrylov=1000) = 120 int3d(Th)( p*q*1e-15+ u1*v1 + u2*v2 + u3*v3 + p*div(v1,v2,v3) + div(u1,u2,u3)*q ) 121 + int3d(Th) ( f*q) 122 - int2d(Th,1)( gd*(v1*N.x +v2*N.y + v3*N.z) ) // int on gamma 123 + on(2,u1=g1n,u2=g2n,u3=g3n); 124 125laplace; 126 127// FFCS: add 3D view parameters 128real[int] CameraPositionValue = [0.0165449,3.23891,-0.991528]; 129real[int] CameraFocalPointValue = [0.5,0.5,0.5]; 130real[int] CameraViewUpValue = [0.671735,0.442219,0.594318]; 131real[int] CutPlaneOriginValue = [0.5,0.5,1.01]; 132real[int] CutPlaneNormalValue = [0.689523,0.722423,0.0516115]; 133plot(P,fill=0,boundary=0,ShowMeshes=1,CutPlane=1, 134 CameraPosition=CameraPositionValue, 135 CameraFocalPoint=CameraFocalPointValue, 136 CameraViewUp=CameraViewUpValue, 137 CutPlaneOrigin=CutPlaneOriginValue, 138 CutPlaneNormal = CutPlaneNormalValue); 139 140laplaceMixte; 141 142real errL2=sqrt(int3d(Th)(square(P-p))) ; 143cout << " int 2 x,yz "<<int2d(Th,2)(x) << " " << int2d(Th,2)(y) << " " << int2d(Th,2)(z) << endl; 144cout << " int 2 gn "<<int2d(Th,2)(g1n) << " " << int2d(Th,2)(g2n) << " " << int2d(Th,2)(g3n) << endl; 145cout << " int 2 U "<<int2d(Th,2)(u1) << " " << int2d(Th,2)(u2) << " " << int2d(Th,2)(u3) << endl; 146cout << " int 2 V "<<int2d(Th,2)(vv1) << " " << int2d(Th,2)(vv2) << " " << int2d(Th,2)(vv3) << endl; 147cout << " int 2 DP "<<int2d(Th,2)(dx(P)) << " " << int2d(Th,2)(dy(P)) << " " << int2d(Th,2)(dz(P)) << endl; 148 149cout << " diff: u Gamma_2 " << sqrt(int2d(Th,2) ( square((g1n*N.x+g2n*N.y+g3n*N.z) - (u1*N.x +u2*N.y + u3*N.z) ) ) ) <<endl; 150cout << " diff: P Gamma_2 " << sqrt(int2d(Th,2) ( square((g1n*N.x+g2n*N.y+g3n*N.z) - (dx(P)*N.x +dy(P)*N.y + dz(P)*N.z) ) ) ) <<endl; 151cout << " diff err L2 :" << errL2 << endl; 152cout << " P L2 :" <<sqrt(int3d(Th)(square(P))) << endl; 153cout << " p L2 :" <<sqrt(int3d(Th)(square(p))) << endl; 154assert(errL2<0.05);