/*
  Exemple of Mixte furmulation for Lame equation (Linear Elasticity)
   
  The Stress and Strian are  H(div div),  
  for a Sym Stress S we have 
      ( div(div S) in H^-1 
  and S.nn are continuous

  The Displacement u are H(Curl) (    u.tau is continuous )
  where n,tau are normal / tangent 

  This formulation is call TD-NSS
*/
if (HaveUMFPACK) //  need UMPACK  computation 
{
load "Element_Mixte"
load "Element_P3"
int nn=10;
mesh Th= square(nn*4,nn,[x*4,y]);

// caoutchoouc 
//real mu     = 0.012; 
//real lambda = 0.4; 
// rho rho=1e-16;
// ALu
real  mu = 3.8;
real lambda= 6.2;
real rho = 5e-4;
//   s11, s12, s22 
real twomu=mu*2;
real mu1= 1./(twomu);
real lambda1 = 1./(2*(2*lambda+twomu)) - 1./(twomu*2);
macro BItr(a,b)  /* S-> a S + b tr(S) Id */  
[ [ a + b , 0 , 0 ,b       ], 
  [     0 , a , 0 ,0         ], 
  [     0 , 0 , a ,0         ], 
  [   b   , 0 , 0 , a+ b  ]    ] //

func A = BItr(twomu,lambda) ;
func A1 = BItr(mu1,lambda1) ;

real gravity = -9.81*rho;



macro Strain(u)  [dx(u#1),(dy(u#1)+dx(u#2))/2,(dy(u#1)+dx(u#2))/2,dy(u#2)] //EOM
macro Strain2(s)  [s#11, s#12, s#12, s#22] //EOM
macro S(s) [s#11,s#12,s#22]//
  macro U(u) [u#1,u#2] //
  macro dn(u) (U(u)'*[N.x,N.y]) //')
  macro dnn(s) ( [N.x,N.y]' * [[ s#11,s#12] , [s#12, s#22] ] *[N.x,N.y]) //')

fespace Wh(Th,[P2,P2]);
fespace Sh(Th,TDNNS1);
fespace Vh(Th,BDM1Ortho);
Sh S(s),S(t);
Vh U(u),U(v);;
Wh U(w),U(ww);

  solve Lame2([s11,s12,s22,u1,u2],[t11,t12,t22,v1,v2])=
  int2d(Th)( 
	       Strain2(s)'*(A1*Strain2(t)) //' 
	    - Strain(u)'*Strain2(t)  //'
	    - Strain(v)'*Strain2(s) //'
     //      + 1e-6*(u1*v2+u2*v2)
	     )
  +  intalledges(Th)( dnn(t)*dn(u) + dnn(s)*dn(v) ) 
  //  - int1d(Th,4) ( 0*dnn(t) ) // here  0 <=> U_d .n * dnn(t) 
  +  int2d(Th) ( gravity*v2) //-int2d(Th,1,2,3) ( dNtau(t). v.tau) 
  +on(1,2,3,s11=0,s12=0,s22=0)
  +on(4,u1=0,u2=0); 



  solve Lame1([w1,w2],[ww1,ww2])=
  int2d(Th)( 
	      Strain(w)'*(A*Strain(ww)) //' 
	     )
 -int2d(Th)(gravity*ww2)
  +on(4,w1=0,w2=0); 
plot([w1,w2], wait=1);
plot([u1,u2], wait=1);

real errl2 = sqrt(int2d(Th)( square(w1-u1) + square(w2-w2)));
cout << " errL2 = " << errl2 << endl;
}



