1//  non linear elasticity model
2//
3//  -------------------------------
4//  with huge utilisation of  macro  new version
5// more simple
6// ---------------------------
7//   optimize version
8// ------------
9//  problem is  find $(uu,vn)$  minimizing  $J$
10//  $ min J(un,vn) = \int 1/2 (F2) -  int Pa * un $
11//   $ dJ(u,u,uu,vv) = int dF2(u,v,uu,vv) df(F2(u,v)) $
12//   where $F2 =  (^t {E}  A {E} ) $,
13//   $E(U) =  1/2 (\nabla U + \nabla U^t + \nabla U^t  \nabla U) $
14//         ($u_1$)
15//  with U=(   )
16//         ($u_2$)
17// so:
18//$$ E_ij = 0.5 ( d_i u_j + d_j u_i ) + \sum_k d_i u_k * d_j*u_k  \leqno(1)$$
19//  the symetric tensor $t_{ij}$ are a vector  [t11,2*t12,t22]
20// date from
21/*
22Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow
23
24Stefan Turek and Jaroslav Hron⋆
25Institute for Applied Mathematics and Numerics, University of Dortmund,
26Vogelpothsweg 87, 44227 Dortmund, Germany
27
28http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.550.1689&rep=rep1&type=pdf
29
30section 4.2 CSM test:
31*/
32
33macro EL(u,v) [dx(u),(dx(v)+dy(u)),dy(v)] // is $[\epsilon_{11},2\epsilon_{12},\epsilon_{22}]$
34
35macro ENL(u,v) [
36(dx(u)*dx(u)+dx(v)*dx(v))*0.5,
37(dx(u)*dy(u)+dx(v)*dy(v))    ,
38(dy(u)*dy(u)+dy(v)*dy(v))*0.5 ] // EOM ENL
39
40macro dENL(u,v,uu,vv) [(dx(u)*dx(uu)+dx(v)*dx(vv)),
41 (dx(u)*dy(uu)+dx(v)*dy(vv)+dx(uu)*dy(u)+dx(vv)*dy(v)),
42 (dy(u)*dy(uu)+dy(v)*dy(vv)) ] //
43
44
45macro E(u,v) (EL(u,v)+ENL(u,v)) // is $[\E{11},\sqrt2\E{12},\E{22}]
46macro dE(u,v,uu,vv) (EL(uu,vv)+dENL(u,v,uu,vv)) //
47macro ddE(u,v,uu,vv,uuu,vvv) dENL(uuu,vvv,uu,vv) //
48macro F2(u,v) (E(u,v)'*A*E(u,v)) //
49macro dF2(u,v,uu,vv)  (E(u,v)'*A*dE(u,v,uu,vv) ) //
50macro ddF2(u,v,uu,vv,uuu,vvv) (
51            (dE(u,v,uu,vv)'*A*dE(u,v,uuu,vvv))
52          + (E(u,v)'*A*ddE(u,v,uu,vv,uuu,vvv))  )// EOM
53
54//  for hyper elasticity problem
55//  -----------------------------
56macro f(u) (u) // end of macro
57macro df(u) (1) // end of macro
58macro ddf(u) (0) // end of macro
59
60// -------------------------------
61//
62//   $  \sigma = 2 \mu E + \lambda tr(E) Id $
63//   $   A(u,v)= \sigma(u):\E(v) $
64//
65//   ( a b )
66//   ( b c )
67//
68//  tr*Id : (a,b,c) -> (a+c,0,a+c)
69// so the associed matrix is:
70//   ( 1 0 1 )
71//   ( 0 0 0 )
72//   ( 1 0 1 )
73// ------------------
74verbosity=0;
75
76// data CMS 1
77real xp=0.35,yp=0.0;//
78real upbench= -0.007187, vpbench= -0.066;
79
80real EE = 1.4e6;
81real rho = 1000;
82real sigma = 0.4;
83real mu = EE/(2*(1+sigma));
84real lambda = EE*sigma/((1+sigma)*(1-2*sigma));
85real gravity = -2;
86int nnn = 10;
87
88mesh Th = square(17*nnn,nnn,[0.35*x,0.02*(y-0.5)]);
89cout<<" Th nt : " << Th.nt<<endl;
90int bottom=1, right=2,upper=3,left=4;
91
92plot(Th);
93
94real a11= 2*mu +  lambda  ;
95real a22= mu ; //  because [0,2*t12,0]' A [0,2*s12,0]  = 2*mu*(t12*s12+t21*s21) = 4*mu*t12*s12
96real a33= 2*mu +   lambda ;
97real a12= 0 ;
98real a13= lambda ;
99real a23= 0 ;
100//  symetric part
101real a21= a12 ;
102real a31= a13 ;
103real a32= a23 ;
104func A = [ [ a11,a12,a13],[ a21,a22,a23],[ a31,a32,a33] ];
105
106
107
108
109
110
111fespace Vh(Th,[P2,P2]);
112
113
114
115
116
117
118Vh [uu,vv], [w,s],[un,vn];
119[un,vn]=[0,0];//  intialisation
120[uu,vv]=[0,0];
121
122
123
124real errb=0;
125// Newton's method
126// ---------------
127for (int i=0;i<10;i++)
128{
129  cout << "Loop " << i << endl;
130  solve NonLin([uu,vv],[w,s])=
131    int2d(Th)(        ddF2(un,vn,uu,vv,w,s)  )
132     - int2d(Th)(  dF2(un,vn,w,s) - rho*gravity*s)
133     + on(left,uu=0,vv=0);
134  ;
135  real res = uu[].linfty ; //  norme  L^2 of [uu,vv]
136  errb = max(abs(un(xp,yp)-upbench) , abs(vn(xp,yp)-vpbench));
137
138  un[] -= uu[];
139  cout << "     Linfty residual = " << res << endl;
140  cout << "     u of A =" << un(xp,yp) << " " << vn(xp,yp)  <<  " bench " << upbench  << " "<< vpbench << " err bench =" << errb << endl;
141  plot(movemesh(Th, [x+un, y+vn]),wait=1,cmm=" i =" +i + " err =" + res);
142
143  if (res<1e-10) break;
144}
145
146cout << " err bench = " << errb << endl;
147//plot([un,vn],wait=1);
148mesh th1 = movemesh(Th, [x+un, y+vn]);
149plot(th1,wait=1,ps="nl-elas.eps");
150assert( errb < 1e-2);