Revert demos/navier

demos/navier/main.cpp

@@ -47,7 +47,7 @@ int main(int argc, char **argv)
   gmshFem.userDefinedParameter(pointsPerSWavelength, "pointsPerSWavelength");
   double lc = 2 * pi / (pointsPerSWavelength * kS);        //.............[m]        characteristic length
 
-  const double rExt = 5.;                                  //.............[m]        external radius of the annulus geometry
+  const double rExt = 2.;                                  //.............[m]        external radius of the annulus geometry
   const double rInt = 1.;                                  //.............[m]        internal radius of the annulus geometry
       
   int ABCOrder = 2;
@@ -203,6 +203,149 @@ int main(int argc, char **argv)
   formulation.galerkin(-w * w * rho * dof(ux), tf(ux), surface, "Gauss2");
   formulation.galerkin(-w * w * rho * dof(uy), tf(uy), surface, "Gauss2");
 
+  //---------------------------------
+  //- Absorbing Boundary Conditions -
+  //---------------------------------
+  if(ABCOrder == 0) {
+    msg::info << "Use zero order Lysmer Kuhlemeyer ABC." << msg::endl;
+      
+    // On ABC boundary int -i w cp rho In u.u' -i w cs rho It u.u'
+    formulation.galerkin(-im * w * cP * rho * nX * nX * dof(ux), tf(ux), gammaInf, "Gauss2");
+    formulation.galerkin(-im * w * cP * rho * nX * nY * dof(uy), tf(ux), gammaInf, "Gauss2");
+    formulation.galerkin(-im * w * cP * rho * nX * nY * dof(ux), tf(uy), gammaInf, "Gauss2");
+    formulation.galerkin(-im * w * cP * rho * nY * nY * dof(uy), tf(uy), gammaInf, "Gauss2");
+
+    formulation.galerkin(-im * w * cS * rho * nY * nY * dof(ux), tf(ux), gammaInf, "Gauss2");
+    formulation.galerkin( im * w * cS * rho * nX * nY * dof(uy), tf(ux), gammaInf, "Gauss2");
+    formulation.galerkin(-im * w * cS * rho * nX * nX * dof(uy), tf(uy), gammaInf, "Gauss2");
+    formulation.galerkin( im * w * cS * rho * nX * nY * dof(ux), tf(uy), gammaInf, "Gauss2");
+  }
+  else if(ABCOrder == 2) {
+    msg::info << "Use high order ABC." << msg::endl;
+    //------------------------------------------------------------------------------------------------------------------------------
+    //- See p 1713 of "A high-order ABC for 2D time-harmonic elastodynamic scattering problems", V. Mattesi, M. Darbas C. Geuzaine -
+    //------------------------------------------------------------------------------------------------------------------------------
+
+
+    //-------------------------------------------------------------------
+    //- STEP 1 : Find v in H^1/2(Gamma) such that v =  Lambda_{1,eps} u -
+    //-------------------------------------------------------------------
+                             
+    // On gammaInf int v.v' - i rho w^2 [  1/kPEps (n vO,v') + 1/kSEps (t v1,v') ] (=0)
+    formulation.galerkin(dof(vx), tf(vx), gammaInf, "Gauss2");
+    formulation.galerkin(dof(vy), tf(vy), gammaInf, "Gauss2");
+
+    formulation.galerkin(-im * rho * pow(w,2) * (1.0/kPEps) * nX * dof(v0), tf(vx), gammaInf, "Gauss2");
+    formulation.galerkin(-im * rho * pow(w,2) * (1.0/kPEps) * nY * dof(v0), tf(vy), gammaInf, "Gauss2");
+
+    formulation.galerkin(-im * rho * pow(w,2) * (1.0/kSEps) * (-nY) * dof(v1), tf(vx), gammaInf, "Gauss2");
+    formulation.galerkin(-im * rho * pow(w,2) * (1.0/kSEps) * ( nX) * dof(v1), tf(vy), gammaInf, "Gauss2");
+
+    // (v0,v0') - sum R_i (h_i,v0') =0
+    formulation.galerkin(dof(v0), tf(v0), gammaInf, "Gauss2");
+
+    for(unsigned int i = 0; i < padeOrder ; ++i) {
+      formulation.galerkin(-padeR[i] * dof(h[i]), tf(v0), gammaInf, "Gauss2");
+
+      // Sl(hl,hl') - (grad_G/kPEps hl,grad_G/kPEps hl') -(n.u,hl') (=0) forall l=1,...,L
+      formulation.galerkin(padeS[i] * dof(h[i]), tf(h[i]), gammaInf, "Gauss2");
+
+      formulation.galerkin(-pow(1.0/kPEps,2) * grad(dof(h[i])), grad(tf(h[i])), gammaInf, "Gauss2");
+      
+      formulation.galerkin(-nX * dof(ux), tf(h[i]), gammaInf, "Gauss2");
+      formulation.galerkin(-nY * dof(uy), tf(h[i]), gammaInf, "Gauss2");
+    }
+    // (v1,v1') - sum Rl (il,v1') =0
+    formulation.galerkin(dof(v1), tf(v1), gammaInf, "Gauss2");
+    for(unsigned int i = 0; i < padeOrder ; ++i) {
+      formulation.galerkin(-padeR[i] * dof(ix[i]), tf(v1), gammaInf, "Gauss2");
+                           
+      // Sl( il,il') -(grad_G/kSEps il,grad_G/kSEps il') - (t. u,il') (=0) forall l=1,...L
+      formulation.galerkin( padeS[i] * dof(ix[i]), tf(ix[i]), gammaInf, "Gauss2");
+
+      formulation.galerkin(-pow(1.0/kSEps,2) * grad(dof(ix[i])), grad(tf(ix[i])), gammaInf, "Gauss2");
+
+      formulation.galerkin(-(-nY) * dof(ux), tf(ix[i]), gammaInf, "Gauss2");
+      formulation.galerkin(-( nX) * dof(uy), tf(ix[i]), gammaInf, "Gauss2");
+    }
+    
+    //-------------------------------------------------------------------
+    //- Step 2 : Find q in H^-1/2(Gamma) st : (I+ Lambda_{2,eps}) q = v -
+    //-------------------------------------------------------------------
+
+    // (q,q') -i[1/ks(grad_G q0,q') - 1/kp (q1 n,q')] - (v,q') (=0)
+    formulation.galerkin( dof(qx), tf(qx), gammaInf, "Gauss2");
+    formulation.galerkin( dof(qy), tf(qy), gammaInf, "Gauss2");
+
+    formulation.galerkin(-im * (1.0/kSEps) * vector< std::complex< double > >(1., 0., 0.) * grad(dof(q0)), tf(qx), gammaInf, "Gauss2");
+    formulation.galerkin(-im * (1.0/kSEps) * vector< std::complex< double > >(0., 1., 0.) * grad(dof(q0)), tf(qy), gammaInf, "Gauss2");
+
+    formulation.galerkin( im * (1.0/kPEps) * nX*(-nY) * vector< std::complex< double > >(1., 0., 0.) * grad(dof(q1)), tf(qx), gammaInf, "Gauss2");
+    formulation.galerkin( im * (1.0/kPEps) * nX * nX  * vector< std::complex< double > >(0., 1., 0.) * grad(dof(q1)), tf(qx), gammaInf, "Gauss2");
+    formulation.galerkin( im * (1.0/kPEps) * nY*(-nY) * vector< std::complex< double > >(1., 0., 0.) * grad(dof(q1)), tf(qy), gammaInf, "Gauss2");
+    formulation.galerkin( im * (1.0/kPEps) * nY * nX  * vector< std::complex< double > >(0., 1., 0.) * grad(dof(q1)), tf(qy), gammaInf, "Gauss2");
+
+    formulation.galerkin(-dof(vx), tf(qx), gammaInf, "Gauss2");
+    formulation.galerkin(-dof(vy), tf(qy), gammaInf, "Gauss2");
+
+    //(q0,q0')=sum Rl (jl,q0') 
+    formulation.galerkin( dof(q0), tf(q0), gammaInf, "Gauss2");
+    for(unsigned int i = 0; i < padeOrder ; ++i) {
+      formulation.galerkin(-padeR[i] * dof(j[i]), tf(q0), gammaInf, "Gauss2");
+
+      //Sl(jl,jl')-(grad/kSEps jl,grad/kSEps jl') - (q.n,jl') =0
+      formulation.galerkin(padeS[i] * dof(j[i]), tf(j[i]), gammaInf, "Gauss2");
+
+        formulation.galerkin(-pow(1.0/kSEps,2) * grad(dof(j[i])), grad(tf(j[i])), gammaInf, "Gauss2");
+
+        formulation.galerkin(-nX * dof(qx), tf(j[i]), gammaInf, "Gauss2");
+        formulation.galerkin(-nY * dof(qy), tf(j[i]), gammaInf, "Gauss2");
+    }
+    // (q1,q1')= sum (Rl kl , q1')
+    formulation.galerkin( dof(q1), tf(q1), gammaInf, "Gauss2"); 
+
+    for(unsigned int i = 0; i < padeOrder ; ++i) {
+      formulation.galerkin(-padeR[i] * dof(k[i]), tf(q1), gammaInf, "Gauss2");
+
+      // (Sl kl,kl') -(grad_G/kPEps kl,grad_G/kPEps kl') - (t.q, kl') (=0)
+      formulation.galerkin(padeS[i] * dof(k[i]), tf(k[i]), gammaInf, "Gauss2");
+
+      formulation.galerkin(-pow(1.0/kPEps,2) * grad(dof(k[i])), grad(tf(k[i])), gammaInf, "Gauss2");
+
+      formulation.galerkin(-(-nY) * dof(qx), tf(k[i]), gammaInf, "Gauss2");
+      formulation.galerkin(-( nX) * dof(qy), tf(k[i]), gammaInf, "Gauss2");
+    }
+
+    //----------------------------------
+    //- Step 3: -tu.u' with t=q+2mu Mu -
+    //----------------------------------
+
+    // - (q,u') - 2 mu (Mu,u')
+    formulation.galerkin(-dof(qx), tf(ux), gammaInf, "Gauss2");
+    formulation.galerkin(-dof(qy), tf(uy), gammaInf, "Gauss2");
+
+    formulation.galerkin(-2.0 * mu * vector< std::complex< double > >(1., 0., 0.) * grad(dof(un)), tf(ux), gammaInf, "Gauss2");
+    formulation.galerkin(-2.0 * mu * vector< std::complex< double > >(0., 1., 0.) * grad(dof(un)), tf(uy), gammaInf, "Gauss2");
+
+    formulation.galerkin( 2.0 * mu * nX*(-nY) * vector< std::complex< double > >(1., 0., 0.) * grad(dof(ut)), tf(ux), gammaInf, "Gauss2");
+    formulation.galerkin( 2.0 * mu * nX * nX  * vector< std::complex< double > >(0., 1., 0.) * grad(dof(ut)), tf(ux), gammaInf, "Gauss2");
+    formulation.galerkin( 2.0 * mu * nY*(-nY) * vector< std::complex< double > >(1., 0., 0.) * grad(dof(ut)), tf(uy), gammaInf, "Gauss2");
+    formulation.galerkin( 2.0 * mu * nY * nX  * vector< std::complex< double > >(0., 1., 0.) * grad(dof(ut)), tf(uy), gammaInf, "Gauss2");
+
+    // (un,un')=(u.n,un')
+    formulation.galerkin( dof(un), tf(un), gammaInf, "Gauss2");
+
+    formulation.galerkin(- nX * dof(ux), tf(un), gammaInf, "Gauss2");
+    formulation.galerkin(- nY * dof(uy), tf(un), gammaInf, "Gauss2");
+
+    // (ut,ut')=(u.t,ut')
+    formulation.galerkin( dof(ut), tf(ut), gammaInf, "Gauss2");
+
+    formulation.galerkin(-(-nY) * dof(ux), tf(ut), gammaInf, "Gauss2");
+    formulation.galerkin(-( nX) * dof(uy), tf(ut), gammaInf, "Gauss2");
+
+  }
+
   // Prepro
   formulation.pre();
   
@@ -214,5 +357,17 @@ int main(int argc, char **argv)
   save(ux);
   save(uy);
 
+  save(xComp(solution), surface, "ux_exact");
+  save(yComp(solution), surface, "uy_exact");
+
+  std::complex< double > num = integrate(pow(abs(xComp(solution) - ux), 2) + pow(abs(yComp(solution) - uy), 2), surface, "Gauss2");
+  std::complex< double > den = integrate(pow(abs(xComp(solution) ), 2) + pow(abs(yComp(solution) ), 2), surface, "Gauss2");
+  
+  std::complex<double> L2error = sqrt(num / den);
+  msg::info << "L_2 error = " << L2error << msg::endl;
+
+  CSVio file("convergence.txt", ' ', OpeningMode::Append);
+  file << 1. / pointsPerSWavelength << std::real(sqrt(num / den)) << csv::endl;
+  
   return 0;
 }