# [PETSc] 3D elasticity problem in which Young's modulus varies in the field

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
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.

1. Use LU until you figure things out
2. Do not use the default `tgv` value of `1E30`

Dear Prof. Jolivet

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

Is it enough?
If you have any suggestions on the best way to get to grips with it in the shortest possible time,
Keywords are fine …

I hope you will consider it.

Is it enough for what?

To figure things out.

I’m sorry for asking such a stupid question…

No, I told you the way to figure things out is to start with a simple LU method.

OK, I understand.
Thank you very much for your guidance.