# \[PETSc\] 3D elasticity problem in which Young's modulus varies in the field **URL:** \ **Category:** General Discussion\ **Created:** [May 16, 2023, 12:58pm UTC](https://community.freefem.org/t/petsc-3d-elasticity-problem-in-which-youngs-modulus-varies-in-the-field/2461 "2023-05-16T12:58:37Z")\ **Posts on this page:** 7\ **Page:** 1
**Author:** ![Keita](https://avatars.discourse-cdn.com/v4/letter/k/977dab/32.png) [@Keita](https://community.freefem.org/u/Keita)\ **Post date:** [May 16, 2023, 12:58pm UTC](https://community.freefem.org/t/petsc-3d-elasticity-problem-in-which-youngs-modulus-varies-in-the-field/2461/1 "2023-05-16T12:58:37Z")
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:** ![230516_AnalysisResult](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/b/b3b88077836c853510a15f198be8782d6bd5be7a.png) 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. ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [May 16, 2023, 1:25pm UTC](https://community.freefem.org/t/petsc-3d-elasticity-problem-in-which-youngs-modulus-varies-in-the-field/2461/2 "2023-05-16T13:25:42Z")
1. Use LU until you figure things out 2. Do not use the default `tgv` value of `1E30` ---
**Author:** ![Keita](https://avatars.discourse-cdn.com/v4/letter/k/977dab/32.png) [@Keita](https://community.freefem.org/u/Keita)\ **Post date:** [May 16, 2023, 2:01pm UTC](https://community.freefem.org/t/petsc-3d-elasticity-problem-in-which-youngs-modulus-varies-in-the-field/2461/3 "2023-05-16T14:01:10Z")
Dear Prof. Jolivet Thank you very much for your reply! First, I will try to understand the basic theoretical part… I plan to gain knowledge mainly from the following: [KSP: Linear System Solvers — PETSc 3.19.1 documentation](https://petsc.org/release/manual/ksp/) Is it enough? If you have any suggestions on the best way to get to grips with it in the shortest possible time, please let me know. Keywords are fine … I hope you will consider it. ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [May 16, 2023, 6:17pm UTC](https://community.freefem.org/t/petsc-3d-elasticity-problem-in-which-youngs-modulus-varies-in-the-field/2461/4 "2023-05-16T18:17:34Z")
Is it enough for what? ---
**Author:** ![Keita](https://avatars.discourse-cdn.com/v4/letter/k/977dab/32.png) [@Keita](https://community.freefem.org/u/Keita)\ **Post date:** [May 17, 2023, 2:02am UTC](https://community.freefem.org/t/petsc-3d-elasticity-problem-in-which-youngs-modulus-varies-in-the-field/2461/5 "2023-05-17T02:02:58Z")
To figure things out. I’m sorry for asking such a stupid question… ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [May 17, 2023, 4:35am UTC](https://community.freefem.org/t/petsc-3d-elasticity-problem-in-which-youngs-modulus-varies-in-the-field/2461/6 "2023-05-17T04:35:31Z")
No, I told you the way to figure things out is to start with a simple LU method. ---
**Author:** ![Keita](https://avatars.discourse-cdn.com/v4/letter/k/977dab/32.png) [@Keita](https://community.freefem.org/u/Keita)\ **Post date:** [May 17, 2023, 6:04am UTC](https://community.freefem.org/t/petsc-3d-elasticity-problem-in-which-youngs-modulus-varies-in-the-field/2461/7 "2023-05-17T06:04:15Z")
OK, I understand. Thank you very much for your guidance.