// file Lame3D.edp
//
// SOLUTION OF THE LAME SYSTEM OF ELASTICITY:
//
//

load "msh3"
		
load "medit"
		
//
// THE DATA:
//

real E = 21.e5, nu = 0.28;                               // Young modulus, Poisson ratio
real mu= E/(2*(1+nu)), lambda = E*nu/((1+nu)*(1-2*nu));  // Lame coefficients
real gamma1=mu, gamma2=mu+lambda;                        // Other coefficients

func f10 = 0.;                                           // External forces
func f20 = 0.;                                           // External forces
func f30 = 0.;                                           // External forces
func g10 = -z;                                           // Applied surface forces
func g20 = 0.;                                           // Applied surface forces
func g30 = 0.;                                           // Applied surface forces

//
// AUXILIAR VARIABLES:
//

real unsq2=1./sqrt(2.), mu2=2.*mu, coef=1000.;

//
// SOME MACRO DEFINITIONS:
//

macro epsilon(u,v,w) [dx(u),dy(v),dz(w),
		unsq2*(dy(u)+dx(v)),unsq2*(dz(u)+dx(w)),unsq2*(dz(v)+dy(w))] // EOM
		
macro div(u,v,w) (dx(u)+dy(v)+dz(w)) // EOM

verbosity = 10;

//
// THE 2D MESH:
//

real x0 = 0., x1 = 1.;                          // $x_0$, $x_1$
real y0 = 0., y1 = 2.;                          // $y_0$, $y_1$

int NInt=100, NSigma=99, NBase=98, NTop=97;     // labels: interior, $\Sigma$, $z = 0$, $z = T$
int Nxtot=3, Nytot=7, Nztot=15;                 // steps in directions $x$, $y$, $z$
real errmesh=.1;                                // initial error for mesh adaptation

border C1(s=x0, x1){x=s; y=y0;     label=NSigma;}
border C2(s=y0, y1){x=x1; y=s;     label=NSigma;}
border C3(s=x1, x0){x=s; y=y1;     label=NSigma;}
border C4(s=y1, y0){x=x0; y=s;     label=NSigma;}
plot(C1(Nxtot)+C2(Nytot)+C3(Nxtot)+C4(Nytot), wait=true);

mesh Th2d=buildmesh(C1(Nxtot)+C2(Nytot)+C3(Nxtot)+C4(Nytot));
plot(Th2d, wait=true);
//savemesh(Th2d,"E-Mesh2D.msh");

//errold=errmesh;
//for (int i=0;i< 4;i++)
//{
//	imesh=i+1;
//	Th2d=adaptmesh(Th2d,f2d0,err=errold);
//	plot(Th2d, cmm="2D MESH NO. "+imesh, wait=true);
//	errold = .5*errold;
//};

//
// THE 3D MESH:
//

real z0 = 0., z1 = 4.;                                 // $_z0$, $z_1$
real px=.5*(x0+x1), py=.5*(y0+y1);
int refint = Th2d(px,py).region;
		
int[int] rup=[refint,NTop,NSigma,NSigma],              // upper face --> NTop
		 rdown=[refint,NBase,NSigma,NSigma],           // lower face --> NBase
		 rmid=[NSigma,NSigma];                         // lateral faces, NSigma--> NSigma

mesh3 Th3d = buildlayers(Th2d, Nztot, zbound=[z0,z1],
		 reffacemid=rmid, reffaceup=rup, reffacelow=rdown);
plot(Th3d,wait=true);
//plot(Th3d,wait=1,ps="mixedTh3d.eps");
//savemesh(Th3d,"E-Mesh3D.mesh");
//medit("3D Mesh", Th3d);

//
// THE FINITE ELEMENT SPACE:
//

fespace Vh(Th3d,P1);
// fespace Vh(Th,P2);
Vh u,v,w,uu,vv,ww,f1=f10,f2=f20,f3=f30,g1=g10,g2=g20,g3=g30;

//
// THE PROBLEM:
//

problem lame3d([u,v,w],[uu,vv,ww]) = int3d(Th3d)(lambda*(div(u,v,w)*div(uu,vv,ww))
        + mu2*(epsilon(u,v,w)'*epsilon(uu,vv,ww)))
        - int3d(Th3d)(f1*uu+f2*vv+f3*ww) - int2d(Th3d,NSigma)(g1*uu+g2*vv+g3*ww)
        + on(NBase, u=0.,v=0.,w=0.);

//
// THE COMPUTATIONS:
//

lame3d;

//
// FIGURES:
//

plot([u,v,w],wait=1,ps="lamevect.eps",coef=coef,cmm="THE DISPLACEMENTS");

plot(u,wait=1,fill=1,value=1,cmm="THE x-DISPLACEMENT");

plot(v,wait=1,fill=1,value=1,cmm="THE y-DISPLACEMENT");

plot(w,wait=1,fill=1,value=1,cmm="THE z-DISPLACEMENT");

mesh3 Th3dd = movemesh3(Th3d,transfo=[x+coef*u, y+coef*v, z+coef*w]);
plot(Th3d,Th3dd, wait=1, cmm="THE DEFORMED MESH");

//mesh3 Thm=movemesh3(Th,transfo=[x+u1*coef,y+u2*coef,z+u3*coef],label=ref2);
//Thm=change(Thm,label=ref2);

real dxmin  = u[].min; 
real dymin  = v[].min; 
real dzmin  = w[].min; 
real dxmax  = u[].max; 
real dymax  = v[].max; 
real dzmax  = w[].max; 

cout << " min displacements: x = " << dxmin <<
		" y = " << dymin << " z = " << dzmin << endl;
cout << " max displacements: x = " << dxmax <<
		" y = " << dymax << " z = " << dzmax << endl;