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.