Dear FreeFem++ developers

I’d like to analyze the elasticity problem in which Young’s modulus varies in the field.

Dirichlet conditions were set for a 3-dimensional cube with a predefined displacement of 1.0 and 0.0 vertically outward on one face and the other faces.

The Young’s modulus of the elastic body is then multiplied by a coefficient that varies with the field. For example, here we consider a coefficient such that a hole is made in the center of the cube (i.e., Young’s modulus becomes very small) using a Gaussian function.

The analysis code is as follows:

load "medit" // for data save load "msh3" load "PETSc" // PETSc plugin macro dimension()3// EOM // '2', '3', '3S' or '3L' include "macro_ddm.idp" // additional DDM functions verbosity=3; // Solid Mechanics real E; E = 1.0; // Young's modulus real nu; nu = 0.3; // Poisson's ratio real lambda; lambda = nu*1./(1+nu)/(1-2.*nu); // Lame coefficient real mu; mu = 1./2./(1+nu); // Lame coefficient // Mesh int[int] labs(6); // labels on boundaries of squre: (bottom, right, top, left) int[int] MeshNum(3); // number of mesh: (x cordinate, y cordinate) real[int] Pos0(3),Pos1(3); // Positions of diagonal points of squre: (x cordinate, y cordinate) real meshdiv; meshdiv = 16; mesh3 Shc1; // Mesh of the squre 1 labs = [1,2,3,4,5,6]; MeshNum = [meshdiv,meshdiv,meshdiv]; Pos0 = [ 0.0 , 0.0 , 0.0]; Pos1 = [ 1.0 , 1.0 , 1.0]; Shc1 = cube(MeshNum(0),MeshNum(1),MeshNum(2),[Pos0(0)+(Pos1(0)-Pos0(0))*x, Pos0(1)+(Pos1(1)-Pos0(1))*y, Pos0(2)+(Pos1(2)-Pos0(2))*z],label=labs,flags=1,region=1); mesh3 Sh=Shc1; mesh3 ShPlt=Sh; mesh3 ShPltMoved; int[int] n2o; macro ShN2O()n2o // EOM // Fespaces fespace Wh(Sh, [P2, P2, P2]), WhPlt(ShPlt, [P2, P2, P2]); fespace WhCoefficient(Sh, P0); // Creat matrix buildDmesh(Sh) Mat A; { macro def(i)[i, i#B, i#C]// EOM // vector field definition macro init(i)[i, i, i]// EOM // vector field initialization createMat(Sh, A, [P2, P2, P2]) } // Solid Mechanics macro defSTR (i)[i, i#B, i#C] // EOM // vector field definition macro epsilon(u) [dx(u),dy(u#B),dz(u#C),(dz(u#B)+dy(u#C)), (dz(u)+dx(u#C)), (dy(u)+dx(u#B))] // strain tensor macro D [ [2.*mu+lambda, lambda, lambda, 0, 0, 0 ], [lambda, 2.*mu+lambda, lambda, 0, 0, 0 ], [lambda, lambda, 2.*mu+lambda, 0, 0, 0 ], [0, 0, 0, mu, 0, 0 ], [0, 0, 0, 0, mu, 0 ], [0, 0, 0, 0, 0, mu ] ] //elastic tensor // define u, v and Coefficient Wh defSTR(u), defSTR(v); WhCoefficient Gauss, Character, Coefficient; // Coefficient real sigma; sigma = 0.05; real width; width = 0.3; real minimun; minimun = 1e-7; Gauss = 1.0-( 1/sqrt(2*pi*sigma^2))*exp(-((x-0.5)^2+(y-0.5)^2+(z-0.5)^2)/(2*sigma^2)); Character = (0.5+Gauss/width*(15./16-Gauss^2/width^2*(5./8-3./16*Gauss^2/width^2)))*(Gauss>=-width)*(Gauss<=width)+1.*(Gauss>width); Coefficient = (1.-minimun)*Character+minimun; // Variational formulation varf vPb(defSTR(u), defSTR(v)) = int3d(Sh)((D*epsilon(u))'*epsilon(v)*E*Coefficient) + on(1, uB = 0.0) + on(2, u = 1.0) + on(3, uB = 0.0) + on(4, u = 0.0) + on(5, uC = 0.0) + on(6, uC = 0.0); // ' // Components of the matrix A A = vPb(Wh, Wh); // Right hand side vector real[int] rhs = vPb(0, Wh); // Solver setting set(A, sparams = "-pc_type gamg -ksp_type cg -ksp_rtol 1e-6 -ksp_max_it 1000 -pc_gamg_threshold 0.01"); // Solve u[] = A^-1 * rhs; // Save the solution to the global solution macro guFirs [gu1,gu2,gu3] // Glabal displacement vector WhPlt [pltx, plty, pltz], [reducex, reducey, reducez], [gu1, gu2, gu3]; real[int] sol; ChangeNumbering(A, u[], sol); ChangeNumbering(A, u[], sol, inverse = true); int[int] rest=restrict(Wh, WhPlt, n2o); for[i, v : rest] pltx[][v] = u[][i]; mpiAllReduce(pltx[], gu1[], mpiCommWorld, mpiSUM); // Deformation visualization if(mpirank==0){ ShPltMoved = movemesh3(ShPlt,transfo=[x+gu1,y+gu2,z+gu3]); plot(ShPltMoved,cmm="Deformation"); }

**Problem**

**The results of the analysis are very strange as follows:**

I assumed that this phenomenon could be avoided by selecting and configuring the solver.

In particular, I would like to know if there is a preferred method for this kind of problem.

Thank you for your cooperation.