/usr/share/doc/freefem++/examples/examples++-mpi/VG.edp is in freefem++-doc 3.19.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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 | // other
load "medit"
load "MUMPS_FreeFem"
include "cube.idp"
real ttgv=1e10;
string ssparams="nprow=1, npcol="+mpisize;
int[int] Nxyz=[40,8,8];
real [int,int] Bxyz=[[0.,5.],[0.,1.],[0.,1.]];
int [int,int] Lxyz=[[1,1],[2,2],[2,2]];
mesh3 Th=Cube(Nxyz,Bxyz,Lxyz);
real E = 21.5e4;
real sigma = 0.29;
real mu = E/(2*(1+sigma));
real lambda = E*sigma/((1+sigma)*(1-2*sigma));
real gravity = -0.05;
fespace Vh(Th,[P2,P2,P2]);
//fespace Vh(Th,[P1,P1,P1]);
Vh [u1,u2,u3], [v1,v2,v3];
cout << "lambda,mu,gravity ="<<lambda<< " " << mu << " " << gravity << endl;
real sqrt2=sqrt(2.);
macro epsilon(u1,u2,u3) [dx(u1),dy(u2),dz(u3),(dz(u2)+dy(u3))/sqrt2,(dz(u1)+dx(u3))/sqrt2,(dy(u1)+dx(u2))/sqrt2] // EOM
macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM
real time=clock();
real tMatrix = clock();
varf vLame([u1,u2,u3],[v1,v2,v3]) =
int3d(Th)(
lambda*div(u1,u2,u3)*div(v1,v2,v3)
+ 2.*mu*( epsilon(u1,u2,u3)'*epsilon(v1,v2,v3) ) //') for emacs
)
+ on(1,u1=0,u2=0,u3=0);
matrix MLame=vLame(Vh,Vh,tgv=ttgv);
tMatrix = clock()-tMatrix;
real tFact = clock();
set(MLame,solver=sparsesolver,tgv=ttgv,sparams=ssparams);
tFact = clock()-tFact;
real tsdc = clock();
varf vsdc([u1,u2,u3],[v1,v2,v3])=int3d(Th) (gravity*v3)+ on(1,u1=0,u2=0,u3=0);
real[int] sdc= vsdc(0,Vh);
tsdc = clock()-tsdc;
real tsolve=clock();
u1[] = MLame^-1*sdc;
tsolve= clock()-tsolve;
cout << "===============================================" << endl;
cout << "==== CPU time =====" << endl;
cout << "===============================================" << endl;
cout << " ALL solving steps :::: " << clock()-time << endl;
cout << " Matrix :::: " << tMatrix << endl;
cout << " Fact :::: " << tFact << endl;
cout << " Second member :::: " << tsdc << endl;
cout << " Solve :::: " << tsolve << endl;
cout << "===============================================" << endl;
real dmax= u1[].max;
cout << " max deplacement = " << dmax << endl;
real coef= 0.1/dmax;
int[int] ref2=[1,0,2,0];
mesh3 Thm=movemesh3(Th,transfo=[x+u1*coef,y+u2*coef,z+u3*coef],label=ref2);
savemesh(Thm,"beam-deformed-mumps.mesh");
|