only compute exact solution for sphere

tutorials/helmholtz/scattering3D/main.cpp

@@ -102,8 +102,6 @@ int main(int argc, char **argv)
     sphere(Rin, Rext, lc);
     // partition
     gmsh::model::mesh::partition(nDom);
-
-    // debug
     gmsh::write("sphere.msh");
   }
   else {
@@ -113,14 +111,17 @@ int main(int argc, char **argv)
         gmsh::merge(mesh + "_" + std::to_string(i + 1) + ".msh");
       }
     }
-
     if(gmsh::model::getNumberOfPartitions() != nDom)
       msg::error << "Wrong number of partitions in mesh file " << gmsh::model::getNumberOfPartitions() << " vs. " << nDom << msg::endl;
   }
 
-  Subdomain omega = Subdomain::buildSubdomainOf(gmshfem::domain::Domain("omega"));
-  Subdomain gammaScat = Subdomain::buildSubdomainOf(gmshfem::domain::Domain("gamma_scat"));
-  Subdomain gammaExt = Subdomain::buildSubdomainOf(gmshfem::domain::Domain("gamma_ext"));
+  Domain omegaMono = gmshfem::domain::Domain("omega");
+  Domain gammaExtMono = gmshfem::domain::Domain("gamma_scat");
+  Domain gammaScatMono = gmshfem::domain::Domain("gamma_ext");
+
+  Subdomain omega = Subdomain::buildSubdomainOf(omegaMono);
+  Subdomain gammaScat = Subdomain::buildSubdomainOf(gammaExtMono);
+  Subdomain gammaExt = Subdomain::buildSubdomainOf(gammaScatMono);
 
   std::vector< std::vector< unsigned int > > topology;
   Interface sigma = Interface::buildInterface(topology);
@@ -132,11 +133,24 @@ int main(int argc, char **argv)
   formulation.addInterfaceField(g);
 
   // Define Analytical solution
-  Function< std::complex< double >, Degree::Degree0 > *solution = nullptr;
-  solution = new AnalyticalFunction< helmholtz3D::ScatteringByAHardSphere< std::complex< double > > >(-k, Rin, 0., 0., 0., 2 * k, theta); // -k because we use e^iwt convention
-  solution->activateMemorization();
+  Function< std::complex< double >, Degree::Degree0 > *u_exact = nullptr;
+  if(mesh.empty()) {
+    u_exact = new AnalyticalFunction< helmholtz3D::ScatteringByAHardSphere< std::complex< double > > >(-k, Rin, 0., 0., 0., 2 * k, theta); // -k because we use e^iwt convention
+    u_exact->activateMemorization();
+  }
+
   // Define incident field
-  Function< std::complex< double >, Degree::Degree0 > uinc = AnalyticalFunction< helmholtz3D::dr_ScatteringByAHardSphere< std::complex< double > > >(-k, Rin, 0., 0., 0., 2 * k, theta);
+  Function< std::complex< double >, Degree::Degree0 > uinc, dn_uinc;
+
+  // TODO actually define uinc in whole omega, so that we can e.g. compute the total field
+  if(mesh.empty()) {
+    uinc = *u_exact;
+    dn_uinc = AnalyticalFunction< helmholtz3D::dr_ScatteringByAHardSphere< std::complex< double > > >(-k, Rin, 0., 0., 0., 2 * k, theta);
+  }
+  else {
+    uinc = exp< std::complex< double > >(im * k * z< std::complex< double > >());
+    dn_uinc = normal< std::complex< double > >() * vector< std::complex< double > > (0, 0, im * k * uinc);
+  }
 
   ScalarFunction< std::complex< double > > Si = 0., dS = 0.;
   if(Transmission == "Taylor0") {
@@ -165,7 +179,7 @@ int main(int argc, char **argv)
     }
 
     // Incident field
-    formulation(i).integral( formulation.physicalSource(uinc) , tf(u(i)), gammaScat(i), gauss);
+    formulation(i).integral( formulation.physicalSource(dn_uinc) , tf(u(i)), gammaScat(i), gauss);
 
     // Artificial source
     for(unsigned int jj = 0; jj < topology[i].size(); ++jj) {
@@ -197,19 +211,21 @@ int main(int argc, char **argv)
   formulation.pre();
   formulation.solve("gmres", tolerence, maxIter, true);
 
+  if(u_exact) {
     // Compute analytic L2-error
     double num = 0., den = 0.;
     double local_err; // L2 local error
     for(auto i = 0; i < nDom; ++i) {
       if(gmshddm::mpi::isItMySubdomain(i)) {
-      std::complex< double > denLocal = integrate(pow(abs(*solution), 2), omega(i), gauss);
-      std::complex< double > numLocal = integrate(pow(abs(*solution - u(i)), 2), omega(i), gauss);
+        std::complex< double > denLocal = integrate(pow(abs(*u_exact), 2), omega(i), gauss);
+        std::complex< double > numLocal = integrate(pow(abs(*u_exact - u(i)), 2), omega(i), gauss);
         num += numLocal.real();
         den += denLocal.real();
         local_err = sqrt(num / den);
         msg::info << "Analytic-ddm L2 error (%) = " << 100.*local_err << " on subdomain " << i << msg::endl;
       }
     }
+  }
 
   // save solutions
   gmsh::option::setNumber("Mesh.Binary", 1);
@@ -217,7 +233,7 @@ int main(int argc, char **argv)
   gmsh::option::setNumber("PostProcessing.SaveMesh", 0);
   for(auto i = 0; i < nDom; ++i) {
     if(gmshddm::mpi::isItMySubdomain(i)) {
-      int tag = save(u(i) - uinc, omega(i), "u", "", "", true, 0, 0., i + 1);
+      int tag = save(u(i), omega(i), "u", "", "", true, 0, 0., i + 1);
       gmsh::view::write(tag, "u_" + std::to_string(i + 1) + ".msh");
 
       // also save a cut
@@ -242,13 +258,6 @@ int main(int argc, char **argv)
   bool ComputeMono = false; // monodomain solution
   gmshDdm.userDefinedParameter(ComputeMono, "ComputeMono");
   if (ComputeMono) {
-    Domain omegaMono, gammaExtMono, gammaScatMono;
-    for(auto i = 0; i < nDom; ++i) {
-      omegaMono |= omega(i);
-      gammaScatMono |= gammaScat(i);
-      gammaExtMono |= gammaExt(i);
-    }
-
     // monodomain solution
     gmshfem::field::Field< std::complex< double >, Form::Form0 > uMono("uMono", omegaMono | gammaScatMono | gammaExtMono, FunctionSpaceTypeForm0::HierarchicalH1, FEMorder);
     gmshfem::problem::Formulation< std::complex< double > > monodomain("Helmholtz");
@@ -264,7 +273,7 @@ int main(int argc, char **argv)
     }
 
     // Incident field
-    monodomain.integral(uinc, tf(uMono), gammaScatMono, gauss);
+    monodomain.integral(dn_uinc, tf(uMono), gammaScatMono, gauss);
     // solve
     monodomain.pre();
     monodomain.assemble();