1 /*******************************************************
 2  *  This file is part of the FreeFEM tutorial made by  *
 3  *      Pierre Jolivet <pierre@joliv.et>               *
 4  *                                                     *
 5  *  See https://joliv.et/FreeFem-tutorial for more     *
 6  *                                                     *
 7  *   Description: system of linear elasticity          *
 8  *******************************************************/
  
10 load "PETSc"
11 macro dimension()3// EOM
12 include "macro_ddm.idp"
 
14 func Pk = [P2, P2, P2];
 
16 macro def(i)[i, i#B, i#C]//
17 macro init(i)[i, i, i]//
18 real Sqrt = sqrt(2.0);
19 macro epsilon(u)[dx(u), dy(u#B), dz(u#C), (dz(u#B) + dy(u#C)) / Sqrt, (dz(u) + dx(u#C)) / Sqrt, (dy(u) + dx(u#B)) / Sqrt]//
20 macro div(u)(dx(u) + dy(u#B) + dz(u#C))//
 
22 int[int] LL = [2,3, 2,1, 2,2];
23 real n = getARGV("-n", 5);
24 mesh3 Th = cube(10 * n, n, n, [10 * x, y, z], label = LL);
25 Mat A;
26 MatCreate(Th, A, Pk);
 
28 real f = -9000.0;
29 real Young = 1.0e8;
30 real poisson = 0.4;
31 real tmp = 1.0 + poisson;
32 real mu = Young  / (2.0 * tmp);
33 real lambda = Young * poisson / (tmp * (1.0 - 2.0 * poisson));
34 varf vPb(def(u), def(v)) = intN(Th)(lambda * div(u) * div(v) + 2.0 * mu * (epsilon(u)' * epsilon(v))) + intN(Th)(f * vC) + on(1, u = 0, uB = 0, uC = 0);
35 fespace Vh(Th, Pk);
36 matrix Loc = vPb(Vh, Vh);
37 A = Loc;
38 real[int] rhs = vPb(0, Vh);
 
40 Vh<real> def(sol);
41 set(A, sparams = "-ksp_converged_reason -ksp_max_it 50 -pc_gamg_threshold 0.1 -pc_type gamg");
42 sol[] = A^-1 * rhs;
 
44 Vh<real> def(Rb)[6];
45 [Rb[0], RbB[0], RbC[0]] = [1, 0, 0];
46 [Rb[1], RbB[1], RbC[1]] = [0, 1, 0];
47 [Rb[2], RbB[2], RbC[2]] = [0, 0, 1];
48 [Rb[3], RbB[3], RbC[3]] = [y, -x, 0];
49 [Rb[4], RbB[4], RbC[4]] = [-z, 0, x];
50 [Rb[5], RbB[5], RbC[5]] = [0, z, -y];
51 set(A, sparams = "-pc_type none");
52 set(A, sparams = "-pc_type gamg", bs = 3, nearnullspace = Rb);
53 sol[] = 0;
54 sol[] = A^-1 * rhs;
55 real alpha = 1.0;
56 meshN ThMoved = movemesh3(Th, transfo = [x + alpha * sol, y + alpha * solB, z + alpha * solC]);
57 plotDmesh(ThMoved, cmm = "Moved mesh");
58 savevtk("moved.vtu", ThMoved, bin = 1);