assert(version>=1.24); //  for big bug is non symetric matrix see HISTORY, and sign in int1d 
include "beam.edp" // 蓋の変形
//変位ベクトルは $(\overline{u},\overline{v})$
// 水槽の上は蓋の歪によって直線ではない
border e(t=0,10) { x=t; y=-10; label= 1; };      //  bottom
border f(t=0,10) { x=10; y=-10+t ; label= 1; };   //  right
border g(t=0,10) { x=0; y=-t ;label= 2;};       //  left
border h(t=0,10) { x=t; y=vv(t,0)*( t>=0.001 )*(t <= 9.999); label=3;};   //  top of cavity deforme 


mesh sh = buildmesh(h(-20)+f(10)+e(10)+g(10));
plot(sh,wait=1);

fespace Xh(sh,P2),Mh(sh,P1);
fespace V1h(sh,P2);

Xh u1,u2,v1,v2;
Mh p,q,ppp;
/* ここで, 発散定理を使い変分形式を定義する
\begin{eqnarray*}
\int_{\Omega}\textrm{div}u q &=& \int_{\partial\Omega}u\cdot n q+b_x(u_1,q)+b_y(u_2,q)\\
b_x(u_1,q) &=& - \int_{\Omega}\partial_1u_1q,
b_y(u_2,q) = - \int_{\Omega}\partial_2u_2q,
\end{eqnarray*} 
*/
varf bx(u1,q) = int2d(sh)( -(dx(u1)*q)); 

varf by(u1,q) = int2d(sh)( -(dy(u1)*q));
//$Lap(u1,u2)=\int_{\Omega}\nabla u_1\cdot \nabla u_2$

varf Lap(u1,u2)= int2d(sh)(  dx(u1)*dx(u2) + dy(u1)*dy(u2) )
                   +  on(2,u1=1) +  on(1,3,u1=0)   ;

varf Mass(p,q)=int2d(sh)(p*q);

Xh bc1; 
Xh brhs;
                   
/*作用素 $\Delta_D$ を作る
$$
\Delta_D(u)=-\Delta u \quad \textrm{with }
u=1\quad\textrm{on }
$$
*/
matrix A= Lap(Xh,Xh,solver=CG); /*
\begin{eqnarray*}
A_{ij}&=&\int_{\Omega}\nabla \varphi_j\nabla \varphi_i (i\ne j または i=j で対応する要素 q^i\in \Omega\\
A_{ii}&=&E\qquad \textrm{if }q^i\in \Gamma_2
\end{eqnarray*} 
*/
matrix M= Mass(Mh,Mh,solver=CG); //  $M_{ij}=\int_{\Omega}\varphi_j\varphi_i$
matrix Bx= bx(Xh,Mh);//  $B_{x,ij}=\int_{\Omega}\partial_1\varphi_j\psi_i$
matrix By= by(Xh,Mh);//  $B_{y,ij}=\int_{\Omega}\partial_2\varphi_j\psi_i$

Xh bcx,bcy;
/* 圧力 pp=$p^m$ が求められているとき
\begin{eqnarray*}
(-\Delta + E\gamma_{\Gamma_2})u_1^{m+1} + &=& \partial_1 p^m + E\gamma_{\Gamma_2}u_{1ex} + f_1\\
(-\Delta + E\gamma_{\Gamma_2})u_2^{m+1} &=& \partial_2 p^m + E\gamma_{\Gamma_2}u_{2ex} + f_2
\end{eqnarray*}
を解く。この弱形式は
\begin{eqnarray*}
\int_{\Omega}\nabla \vec{u}^{m+1}\cdot \nabla\vec{v}+ E\gamma_{\Gamma_2}1&=& 
\int_{\Omega}p^m\textrm{div}\vec{v}+\int_{\Omega}\vec{f}\cdot \vec{v}
+ E(\gamma_{\Gamma_2}\vec{u}_{ex})\\
&=&-\int_{\Omega}\nabla p^m\cdot\vec{v}+\int_{\Omega_2}\vec{f}\cdot \vec{v}
+ E(\gamma_{\Gamma_2}\vec{u}_{ex})
\end{eqnarray*}
*/
func real[int] divup(real[int] & pp)
{ //  
  int verb=verbosity;
   verbosity=1;
   brhs[]  = Bx'*pp;  // $b_j = p_i^{m}B_{x,ij}$
   brhs[] += bc1[] .*bcx[]; // $b_j = b_j + (E\gamma_{\Gamma}u_{1ex})_j $
   u1[] = A^-1*brhs[]; // $u_1^{m+1}=\Delta_D^{-1}\partial_x p^{m}$
   brhs[]  = By'*pp; //$b_j = p_i^{m}B_{y,ij}$
   brhs[] += bc1[] .*bcy[];
   u2[] = A^-1*brhs[]; //$u_2^{m+1}=\Delta_D^{-1}\partial_y p^{m}$
   ppp[] = M^-1*pp; //$p^{m}|_{q^{i}}$
   ppp[] = 1.e-6*ppp[]; //$p^{m+1}=10^{-6}p^{m}$
   ppp[] =   Bx*u1[];//$-\textrm{div}u_1^{m+1}|_{q^{i}}$
   ppp[] +=  By*u2[];
   verbosity=verb;
   return ppp[] ; // return $-\textrm{div}\vec{u}^{m+1}$
};

//初期条件 圧力,流速がゼロ
p=0;q=0;u1=0;v1=0;

for(int i=0;i<2;i++) {    
  bcx=0;  //境界での流速 (0,(-y)*(10.+y)/25.)
  bcy= (-y)*(10.+y)/25.;

  cout << " loop " << i << endl; 
  bc1[] = Lap(0,Xh); //境界条件作用素
  q=0;   
  verbosity=50;
  // $p^{m}\to -\textrm{div}\vec{u}^{m+1}, -\textrm{div}u^{\infty}=0$ を解く
  LinearCG(divup,p[],eps=1.e-3,nbiter=50);
 divup(p[]); // $\vec{u}^{\infty}$ 解を得る
 verbosity=1;
 // 流速ベクトル $(u_1,u_2)$ と圧力 $p$ の表示
 plot([u1,u2],p,wait=1,value=true,coef=0.1,cmm="[u1,u2],p");

 real coef=0.2;

V1h sigma11,sigma22,sigma12; //応力 $\sigma_{ij}$
Vh [uu1,vv1]; //$(\overline{u}_1,\overline{v}_1)$
  [uu1,vv1]=[uu,vv];

/*蓋の変形後形状に流速による力が加わる\\
$\sigma_{11}(x+\overline{u},y+\overline{v})=2\partial_x u_1(x,y)-p(x,y)$,
$\sigma_{22}(x+\overline{u},y+\overline{v})=2\partial_y u_2(x,y)-p(x,y)$\\
$\sigma_{12}(x+\overline{u},y+\overline{v})=\textrm{div}\vec{u}(x,y)$ */
  sigma11([x+uu,y+vv]) = (2*dx(u1)-p);
  sigma22([x+uu,y+vv]) = (2*dy(u2)-p);
  sigma12([x+uu,y+vv]) = (dx(u1)+dy(u2));

//流体の影響を考慮していたの変形を解く
solve  bbst([uu,vv],[w,s],init=i)  = 
    int2d(th)(
    2*mu*(dx(uu)*dx(w)+dy(vv)*dy(s)+ ((dx(vv)+dy(uu))*(dx(s)+dy(w)))/2 )
               + lambda*(dx(uu)+dy(vv))*(dx(w)+dy(s))
             )
  + int2d(th) (-gravity*s)  // 内部力として重力が作用
  + int1d(th,bottombeam)( coef*(sigma11*N.x*w + sigma22*N.y*s + sigma12*(N.y*w+N.x*s) )  )
  + on(1,uu=0,vv=0) //左側を固定
  ;
 plot([uu,vv],wait=1);
 real  err = sqrt(int2d(th)( (uu-uu1)^2 + (vv-vv1)^2 )); 
 cout << " ----- Iteration = " << i <<  " Erreur L2 = " << err << "----------\n";
 th1 = movemesh(th, [x+uu, y+vv]);
 plot(th1,wait=1);
 sh = buildmesh(h(-20)+f(10)+e(10)+g(10));
 plot(sh);
 p=p;q=q;u1=u1;u2=u2;brhs=brhs;ppp=ppp;
 A= Lap(Xh,Xh,solver=CG); 
 M= Mass(Mh,Mh,solver=CG); 
 Bx= bx(Xh,Mh);
 By= by(Xh,Mh);
 bc1=0;  // for resize
 bc1[] = Lap(0,Xh);
}