[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
load "msh3"
load "PETSc"                        // PETSc plugin
macro dimension()3// EOM            // '2', '3', '3S' or '3L'
include "macro_ddm.idp"             // additional DDM functions

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

	ShPltMoved = movemesh3(ShPlt,transfo=[x+gu1,y+gu2,z+gu3]);

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.

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

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

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.

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.