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: primal substructuring method         *
 8  *******************************************************/
  
10 load "hpddm_substructuring"
11 macro dimension()2// EOM
12 include "macro_ddm_substructuring.idp"
 
14 macro def(i)[i, i#B]// EOM
15 macro init(i)[i, i]// EOM
16 real Sqrt = sqrt(2.0);
17 macro epsilon(u)[dx(u), dy(u#B), (dy(u) + dx(u#B)) / Sqrt]// EOM
18 macro div(u)(dx(u) + dy(u#B))// EOM
19 macro BC(u, val)u = val, u#B = val// EOM
20 func Pk = [P2, P2];
 
22 meshN Th;
23 fespace Wh(Th, Pk);
24 int[int][int] intersection;
25 int[int] interfaceNb;
26 {
27     int[int] l = [2, 1, 2, 2];
28     Th = square(2 * getARGV("-global", 10), getARGV("-global", 10), [2 * x, y], label = l);    // global mesh
29     buildSubstructuring(Th, interfaceNb, 10, 1, 2, 1, intersection, Pk, BC, mpiCommWorld, false);
30 }
 
32 real f = -90000.0;
33 real strain = 100.0;
34 real Young = 1.0e8;
35 real poisson = 0.45;
36 real tmp = 1.0 + poisson;
37 real mu = Young  / (2.0 * tmp);
38 real lambda = Young * poisson / (tmp * (1.0 - 2.0 * poisson));
39 varf vPb(def(u), def(v)) = intN(Th)(lambda * div(u) * div(v) + 2.0 * mu * (epsilon(u)' * epsilon(v))) + intN(Th)(f * vB) + on(1, u = 0.0, uB = 0.0);
40 matrix<real> Loc = vPb(Wh, Wh, sym = 1);
41 real[int] rhs = vPb(0, Wh);
 
43 bdd A(Loc, intersection);
44 Wh[int] def(Rb)(0);
45 real[int] float(Wh.ndof);
46 varf floatingPb(def(u), def(v)) = on(1, BC(u, 1.0));
47 float = floatingPb(0, Wh);
48 if(float.max < 0.9) {
49     Rb.resize(3);
50     [Rb[0], RbB[0]] = [1, 0];
51     [Rb[1], RbB[1]] = [0, 1];
52     [Rb[2], RbB[2]] = [y, -x];
53 }
54 renumber(A, Loc, interfaceNb, R = Rb, effort = rhs);
55 AttachCoarseOperator(mpiCommWorld, A, R = Rb);
56 Wh def(u);
57 u[] = A^-1 * rhs;
58 OriginalNumbering(A, u[], interfaceNb);
59 plotMPI(Th, def(u), Pk, def, real, cmm = "Solution");