// file Lame-Evol.edp
//
// SOLUTION OF THE EVOLUTION LAME SYSTEM OF ELASTICITY:
//
// 

//
// 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 a=0.;                                              // Damping coefficient
real T=10.;                                             // Final time


func g1 = 0.;                                           // Applied surface forces, Comp. 1
func g2 = 0.;

func u0 = 0.;                                           // Inital data
func v0 = 0.;
func u1 = 0.;
func v1 = 0.;

real f0 = -5., omega = 10.*pi, k0 = 1.1;                 // Data for external forces

// --------------------------------------------------------------------------------------------
// ROUTINE: fext (THE RIGHT HAND SIDE - EXTERNAL FORCES)
//---------------------------------------------------------------------------------------------
// IN:   time
// OUT:  ftime (THE RIGHT HAND SIDE AT TIME time)
//---------------------------------------------------------------------------------------------
func real fext(real time)
{
// real ftime = f0*sin(omega*time);
// real ftime = exp(k0*time)*f0*sin(omega*time);
real ftime = exp(-k0*time)*f0*sin(omega*time);

return ftime;
}
// --------------------------------------------------------------------------------------------
// END OF ROUTINE
//---------------------------------------------------------------------------------------------

// DATA FOR TIME APPROXIMATION
int Mtime=200;

//
// AUXILIAR VARIABLES
//
//real hx=x2-x1, hy=y2-y1;
//real c2=c*c;
real sqrt2=sqrt(2.), usqrt2=1./sqrt2, mu2=2.*mu, scale=100.;
real dt=T/Mtime, dt2=dt*dt, atime0=a/dt, atime1=1./dt2, alpha=atime0+atime1,
    atime2=2.*atime1+atime0, tt=0.;
int Mtime1=Mtime+1, kiter, jiter=0;

//
// THE MESH:
//

mesh Th=square(20,10,[20*x,2*y-1]);    // The domain: (0.,20) X (-1.,1.)
// mesh Th=square(100,50,[20*x,2*y-1]);    // The domain: (0.,20) X (-1.,1.)

//
// THE FINITE ELEMENT SPACE:
//

fespace Vh(Th,P1);
// fespace Vh(Th,P2);
Vh u, v, uu, vv, uold, vold, uvold, vvold, fext1, fext2;
Vh[int] U(Mtime1), V(Mtime1);

//
// SOME MACRO DEFINITIONS:
//

macro epsilon(u1,u2) [dx(u1),dy(u2),usqrt2*(dy(u1)+dx(u2))] // EOM
		
macro div(u,v) (dx(u)+dy(v)) // EOM

// -------------------------------------
// THE LINEAR STATIONARY PROBLEMS
// -------------------------------------
problem lamet([u,v],[uu,vv],init=jiter) = int2d(Th)(alpha*(u*uu+v*vv)
        + lambda*div(u,v)*div(uu,vv)
        + mu2*(epsilon(u,v)'*epsilon(uu,vv)))
        - int2d(Th)(fext1*uu+fext2*vv) - int1d(Th,1)(g1*uu+g2*vv)
        + on(4,u=0.,v=0.);

//
// SOLVING THE PROBLEM
//

// STARTING

kiter=0;
U[0]=u0; V[0]=v0;
plot(U[0],cmm="Comp.1 - ITERATE "+kiter,value=true,fill=true,wait=true);
plot(V[0],cmm="Comp.2 - ITERATE "+kiter,value=true,fill=true,wait=true);

kiter=1;
U[1]=u0+dt*u1; V[1]=v0+dt*v1;
plot(U[1],cmm="Comp.1 - ITERATE "+kiter,value=true,fill=true,wait=true);
plot(V[1],cmm="Comp.2 - ITERATE "+kiter,value=true,fill=true,wait=true);

uvold=U[0]; vvold=V[0];
uold=U[1]; vold=V[1];

// -------------------------------------
// TIME ITERATES
// -------------------------------------

for (int k=2;k<=Mtime;k++)
{
	kiter=k;
	tt=k*dt;
//	z=k*dt;
//	fext=force;
	fext1=atime2*uold-atime1*uvold;
	fext2=fext(tt)+atime2*vold-atime1*vvold;
	lamet;
	
/*
	plot(u,cmm="Comp.1 - ITERATE "+kiter,value=true,fill=true,wait=true,
	    dim=3,nbiso=40);
	plot(v,cmm="Comp.2 - ITERATE "+kiter,value=true,fill=true,wait=true,
	    dim=3,nbiso=40);
*/

	mesh Thd = movemesh(Th,[x+scale*u, y+scale*v]);
//    plot(Thd,wait=1,ps="lamedeform.eps",cmm="THE DEFORMED MESH");
    plot(Th,Thd,wait=1,cmm="THE DEFORMED MESH");

	U[kiter]=u; V[kiter]=v;

/*
	if (kiter % nadapt == 0)
	{
		nmesh=nmesh+1; jiter=-1;
		Th=adaptmesh(Th,[u,v],err=errmesh);
		plot(Th,cmm=" MESH NO. "+nmesh,wait=true);
		u=u;
	};
*/
	uvold=uold; vvold=vold;
	uold=u; vold=v;
	jiter=1;
}

/*
// HERE HERE HERE
//
// DISPLAY, RECORD AND PRINT
//

plot(u,value=true,fill=true,cmm="FINAL ITERATE",wait=1);

// WRITING


{ // save solution
    ofstream file("sol.txt");
    file << Mtime << endl;
	for (int kwrite=0;kwrite<=Mtime;kwrite++)
	{
		time=kwrite*dt;
		file << kwrite << time << endl;
		file << U[kwrite] << endl;
	}
} // close writing file (end of block)


{ // read
    ifstream file("sol.txt");
    file >> Mtime >> endl;
	for (int kwrite=0;kwrite<=Mtime;kwrite++)
	{
		file >> kwrite >> time >> endl;
		file >> U[kwrite] >> endl;
	}
/} // close reading file (end of block)





//
// FIGURES:
//

plot([u,v],wait=1,coef=scale,cmm="THE DISPLACEMENTS");
//plot([u,v],wait=1,ps="lamevect.eps",coef=scale,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");
mesh Thd = movemesh(Th,[x+scale*u, y+scale*v]);
//plot(Thd,wait=1,ps="lamedeform.eps",cmm="THE DEFORMED MESH");
plot(Thd,wait=1,cmm="THE DEFORMED MESH");

real dxmin  = u[].min; 
real dymin  = v[].min; 

cout << " min displacements:    x = "<< dxmin<< " y=" << dymin << endl; 
cout << " displacement at (20,0)  = "<< u(20,0) << " " << v(20,0) << endl;
*/