// run with MPI: ff-mpirun -np 4 script.edp // NBPROC 4 include "cube.idp" load "PETSc" // PETSc plugin load "gmsh" macro dimension()3// EOM // 2D or 3D macro vectorialfe()P1// EOM include "macro_ddm.idp" // additional DDM functions macro def(i)[i, i#B, i#C]// EOM // vector field definition macro init(i)[i, i, i]// EOM // vector field initialization real Sqrt = sqrt(2.0); 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]// EOM macro div(u)(dx(u) + dy(u#B) + dz(u#C))// EOM func Pk = [vectorialfe, vectorialfe, vectorialfe]; // finite element space // // int[int] LL = [0,1, 2,3, 4,5]; // mesh3 Th = cube(10 * getARGV("-global", 5), getARGV("-global", 5), getARGV("-global", 5), [10 * x, y, z], label = LL); mesh3 Th = gmshload3("wing_unst_15W.msh"); plot(Th,wait=1); // int Y = 1;real k = 1;real k1 = 1;real convergeN = 1; // int[int] Nxyz=[20,5,5*k];real [int,int] Bxyz=[[0.,15*k1],[0.,5*k1],[0.,5*k1]];int [int,int] Lxyz=[[10,5],[4,3],[2,1]]; mesh3 Th10=Cube(Nxyz,Bxyz,Lxyz); // Nxyz=[20,1,2*k]; Bxyz=[[0.,15*k1],[5*k1,6*k1],[0.,2*k1]]; Lxyz=[[11,5],[4,3],[2,1]]; mesh3 Th11=Cube(Nxyz,Bxyz,Lxyz); // Nxyz=[20,1,1*k]; Bxyz=[[0.,15*k1],[5*k1,6*k1],[2*k1,3*k1]]; Lxyz=[[12,5],[4,3],[2,1]]; mesh3 Th12=Cube(Nxyz,Bxyz,Lxyz); // Nxyz=[20,1,2*k]; Bxyz=[[0.,15],[5*k1,6*k1],[3*k1,5*k1]]; Lxyz=[[13,5],[4,3],[2,1]]; mesh3 Th13=Cube(Nxyz,Bxyz,Lxyz); // Nxyz=[20,5,5*k]; Bxyz=[[0.,15*k1],[6*k1,11*k1],[0.,5*k1]]; Lxyz=[[14,5],[4,3],[2,1]]; mesh3 Th14=Cube(Nxyz,Bxyz,Lxyz); // mesh3 Th = Th10+Th11+Th12+Th13+Th14; Mat A; buildMat(Th, getARGV("-split", 1), A, Pk, mpiCommWorld, 3) real f = -900.0; real strain = 100.0; real Young = 1.0e8; real poisson = 0.45; real tmp = 1.0 + poisson; real mu = Young / (2.0 * tmp); real lambda = Young * poisson / (tmp * (1.0 - 2.0 * poisson)); varf vPb(def(u), def(v)) = int3d(Th)(lambda * div(u) * div(v) + 2.0 * mu * (epsilon(u)' * epsilon(v))) + int2d(Th,3,4)(f * vC) + on(2, u = 0.0, uB = 0.0, uC = 0.0); fespace Wh(Th, Pk); // local finite element space matrix Loc = vPb(Wh, Wh); real[int] rhs = vPb(0, Wh); set(A, sparams = "-ksp_max_it 100"); Wh def(u); // local solution A = Loc; u[] = A^-1 * rhs; real[int] err = A * u[]; // global matrix-vector product exchange(A, rhs, scaled = true); err -= rhs; macro def1(u)u// EOM plotMPI(Th, u, vectorialfe, def1, real, cmm = "Global solution") u[] = err; plotMPI(Th, u, vectorialfe, def1, real, cmm = "Global residual") Wh def(Rb)[6]; [Rb[0], RbB[0], RbC[0]] = [1, 0, 0]; [Rb[1], RbB[1], RbC[1]] = [0, 1, 0]; [Rb[2], RbB[2], RbC[2]] = [0, 0, 1]; [Rb[3], RbB[3], RbC[3]] = [y, x, 0]; [Rb[4], RbB[4], RbC[4]] = [z, 0, x]; [Rb[5], RbB[5], RbC[5]] = [0, z, y]; set(A, sparams = "-pc_type gamg -ksp_type gmres -ksp_max_it 200 -pc_gamg_threshold 0.01", nearnullspace = Rb); u[] = 0; u[] = A^-1 * rhs; plotMPI(Th, u, vectorialfe, def1, real, cmm = "Global solution") real alpha = 1.0; mesh3 ThMoved = movemesh3(Th, transfo = [x + alpha * u, y + alpha * uB, z + alpha * uC]); u[] = mpirank; plotMPI(ThMoved, u, vectorialfe, def1, real, cmm = "Global moved solution")