Finite Element

 1load "msh3"
 2load "medit"
 3
 4// Parameters
 5searchMethod=1; // More safe seach algo
 6real a = 1, d = 0.5, h = 0.5;
 7int nnb = 7, nni = 10;
 8int nz = 3;
 9func zmin = 0;
10func zmax = h;
11
12// Mesh 2D
13border b1(t=0.5, -0.5){x=a*t; y=-a/2; label=1;}
14border b2(t=0.5, -0.5){x=a/2; y=a*t; label=2;}
15border b3(t=0.5, -0.5){x=a*t; y=a/2; label=3;}
16border b4(t=0.5, -0.5){x=-a/2; y=a*t; label=4;}
17border i1(t=0, 2.*pi){x=d/2*cos(t); y=-d/2*sin(t); label=7;}
18mesh 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
37int[int] rup = [0, 5], rlow = [0, 6], rmid = [1, 1, 2, 2, 3, 3, 4, 4, 7, 7], rtet = [0, 41];
38mesh3 Th3 = buildlayers(Th, nz, zbound=[zmin, zmax],
39    reftet=rtet, reffacemid=rmid, reffaceup=rup, reffacelow=rlow);
40for(int i = 1; i <= 6; ++i)
41    cout << " int " << i << " : " << int2d(Th3,i)(1.) << " " << int2d(Th3,i)(1./area) << endl;
42
43plot(Th3, wait=true);
44medit("Th3", Th3);
45
46fespace Vh(Th3, P2, periodic=[[1, x, z], [3, x, z], [2, y, z], [4, y, z], [5, x, y], [6, x, y]]);

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