Finite Element

 1 load "msh3"
 2 load "medit"
 3 
 4 // Parameters
 5 searchMethod=1; // More safe seach algo
 6 real a = 1, d = 0.5, h = 0.5;
 7 int nnb = 7, nni = 10;
 8 int nz = 3;
 9 func zmin = 0;
10 func zmax = h;
11 
12 // Mesh 2D
13 border b1(t=0.5, -0.5){x=a*t; y=-a/2; label=1;}
14 border b2(t=0.5, -0.5){x=a/2; y=a*t; label=2;}
15 border b3(t=0.5, -0.5){x=a*t; y=a/2; label=3;}
16 border b4(t=0.5, -0.5){x=-a/2; y=a*t; label=4;}
17 border i1(t=0, 2.*pi){x=d/2*cos(t); y=-d/2*sin(t); label=7;}
18 mesh Th = buildmesh(b1(-nnb) + b3(nnb) + b2(-nnb) + b4(nnb) + i1(nni));
19 
20 { // Cleaning the memory correctly
21     int[int] old2new(0:Th.nv-1);
22     fespace Vh2(Th, P1);
23     Vh2 sorder = x + y;
24     sort(sorder[], old2new);
25     int[int] new2old = old2new^-1; // Inverse permutation
26     Th = change(Th, renumv=new2old);
27     sorder[] = 0:Th.nv-1;
28 }
29 {
30     fespace Vh2(Th, P1);
31     Vh2 nu;
32     nu[] = 0:Th.nv-1;
33     plot(nu, cmm="nu=", wait=true);
34 }
35 
36 // Mesh 3D
37 int[int] rup = [0, 5], rlow = [0, 6], rmid = [1, 1, 2, 2, 3, 3, 4, 4, 7, 7], rtet = [0, 41];
38 mesh3 Th3 = buildlayers(Th, nz, zbound=[zmin, zmax],
39     reftet=rtet, reffacemid=rmid, reffaceup=rup, reffacelow=rlow);
40 for(int i = 1; i <= 6; ++i)
41     cout << " int " << i << " : " << int2d(Th3,i)(1.) << " " << int2d(Th3,i)(1./area) << endl;
42 
43 plot(Th3, wait=true);
44 medit("Th3", Th3);
45 
46 fespace Vh(Th3, P2, periodic=[[1, x, z], [3, x, z], [2, y, z], [4, y, z], [5, x, y], [6, x, y]]);

 1 // Parameters
 2 func f = 1 + x - y;
 3 
 4 // Mesh
 5 mesh Th = square(10, 10);
 6 
 7 // Fespace
 8 fespace Vh(Th, P1);
 9 int n = Vh.ndof;
10 int n1 = n+1;
11 Vh uh, vh;
12 
13 // Problem
14 varf va (uh, vh)
15     = int2d(Th)(
16           dx(uh)*dx(vh)
17         + dy(uh)*dy(vh)
18     )
19     ;
20 
21 varf vL (uh, vh) = int2d(Th)(f*vh);
22 varf vb (uh, vh) = int2d(Th)(1.*vh);
23 
24 matrix A = va(Vh, Vh);
25 real[int] b = vL(0, Vh);
26 real[int] B = vb(0, Vh);
27 
28 // Block matrix
29 matrix AA = [ [ A, B ], [ B', 0 ] ];
30 set(AA, solver=sparsesolver);
31 
32 real[int] bb(n+1), xx(n+1), b1(1), l(1);
33 b1 = 0;
34 // Builds the right hand side block
35 bb = [b, b1];
36 
37 // Solve
38 xx = AA^-1 * bb;
39 
40 // Set values
41 [uh[],l] = xx;
42 
43 // Display
44 cout << " l = " << l(0) << " , b(u, 1) =" << B'*uh[] << endl;
45 
46 // Plot
47 plot(uh);