-
Notifications
You must be signed in to change notification settings - Fork 188
/
elasticity-3d-simple.edp
88 lines (56 loc) · 2.42 KB
/
elasticity-3d-simple.edp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
//ff-mpirun -np 8 elasticity-3d-simple.edp -wg -ffddm_schwarz_method ras -ffddm_schwarz_coarse_correction BNN -ffddm_geneo_nu 20 -global 10
// If you have openmpi you may need to add the option --oversubscribe to allow more processes than the number of cores available on your computer
// for the make check:
// NBPROC 4
// PARAM -raspart -ffddm_schwarz_method ras -ffddm_schwarz_coarse_correction BNN -ffddm_geneo_nu 20 -global 5
macro dimension 3// EOM // 2D or 3D
include "ffddm.idp"
load "msh3"
macro def(i)[i, i#B, i#C]// EOM // vector field definition
macro init(i)[i, i, i]// EOM // vector field initialization
real Sqrt = sqrt(2.0);
macro epsilon(u)[dx(u), dy(u#B), dz(u#C), (dz(u#B) + dy(u#C)) / Sqrt, (dz(u) + dx(u#C)) / Sqrt, (dy(u) + dx(u#B)) / Sqrt]// EOM
macro div(u)(dx(u) + dy(u#B) + dz(u#C))// EOM
func Pk = [P1,P1,P1]; // finite element space
int[int] LL = [2,3, 2,1, 2,2];
mesh3 ThGlobal = cube(6 * getARGV("-global", 5), getARGV("-global", 5), getARGV("-global", 5), [6 * x, y, z], label = LL);
real f = -9000.0;
real strain = 100.0;
real Young = 2.0e11;
real poisson = 0.35;
real tmp = 1.0 + poisson;
real mu = Young / (2.0 * tmp);
real lambda = Young * poisson / (tmp * (1.0 - 2.0 * poisson));
int dirichlet = 1;
macro Varf(varfName, meshName, PhName)
varf varfName(def(u), def(v)) = int3d(meshName)(lambda * div(u) * div(v) + 2.0 * mu * (epsilon(u)' * epsilon(v))) + int3d(meshName)(f * vC) + on(dirichlet, u = 0.0, uB = 0.0, uC = 0.0); // EOM
vtgv = -2;
vtgvelim = -2;
vsym = 1;
ffddmbuild(E,ThGlobal,real,def,init,Pk,mpiCommWorld)
//macro Ewithhpddm()1//
macro Edefmplot(u)u//
ffddmsetup(E,E,Varf,null)
real[int] rhs(1);
ffddmbuildrhs(E,Varf,rhs)
real[int] x0(rhs.n);
x0 = 0;
EVhi def(u), def(err);
//set(EhpddmOP,sparams="-hpddm_reuse_preconditioner 1");
if (mpirank == 0) cout << "RAS :" << endl;
u[] = EfGMRES(x0, rhs, 1.e-6, 200, "right");
Ewritesummary
if (mpirank == 0) cout << endl << "RAS + GENEO :" << endl;
ffddmgeneosetup(E,Varf)
u[] = EfGMRES(x0, rhs, 1.e-6, 200, "right");
Ewritesummary
err[] = EA(u[]);
err[] -= rhs;
ffddmplot(E,u, "Global solution");
ffddmplot(E,err, "Global residual");
EVhglob def(uglob);
EfromVhi(u[],EVhglob,uglob[])
real alpha = 20000.0;
EThglob = movemesh3(EThglob, transfo = [x + alpha * uglob, y + alpha * uglobB, z + alpha * uglobC]);
u[] = mpirank;
ffddmplot(E,u, "Deformed mesh");