Skip to content
Snippets Groups Projects
Select Git revision
  • 19bbbe2de922ac084f07c6e400a59190112dde5e
  • master default protected
  • rgpu
  • oras_vs_osm
  • refactor_coupled
  • lumi-stable
  • fix-compile-without-mpi
  • clean_multirhs
  • oras_comp
  • hpddm_integration
  • blockProduct
  • multiSrcs
  • splitPrePro
  • reuseGCR
  • helmholtz_2d_ddm
  • fix-template-instanciantion-clang-macos
  • customSchwarz
  • hp-convergence-test
  • fix_krylov
  • solverCorrection
  • boris-martin-master-patch-52103
  • gmshddm_1_0_0
22 results

monoDomain.cpp

Blame
  • monoDomain.cpp 31.94 KiB
    #include "mesh.h"
    #include "SubproblemDomains.h"
    #include "Subproblem2D.h"
    #include "Subproblem3D.h"
    
    #include <gmshfem/GmshFem.h>
    #include <gmshfem/FieldInterface.h>
    #include <gmshfem/Formulation.h>
    #include <gmshfem/AnalyticalFunction.h>
    #include <gmshfem/Post.h>
    #include <gmshfem/Function.h>
    #include <gmshfem/io.h>
    
    #include <gmshddm/GmshDdm.h>
    #include <gmshddm/Subdomain.h>
    #include <gmshddm/Interface.h>
    #include <gmshddm/SubdomainField.h>
    #include <gmshddm/InterfaceField.h>
    #include <gmshddm/Formulation.h>
    #include <gmshddm/MPIInterface.h>
    
    #include <algorithm>
    #include <fstream>
    
    using gmshfem::equation::dof;
    using gmshfem::equation::tf;
    using gmshfem::function::operator-;
    using gmshfem::function::operator*;
    using gmshfem::function::abs;
    
    namespace D2 {
    
    
      void monoDomain()
      {
        gmshddm::common::GmshDdm *gmshDdm = gmshddm::common::GmshDdm::currentInstance();
        
        // ************************
        // P H Y S I C S
        // ************************
        double pi = 3.141592653589793238462643383279;
        double k = 4. * pi;
        gmshDdm->userDefinedParameter(k, "k");
        double R = 0.5;
        gmshDdm->userDefinedParameter(R, "R");
        
        // ************************
        // M E S H
        // ************************
        double lc = 0.03333333333333; // Marmousi = 20; other = 0.03333333333333
        gmshDdm->userDefinedParameter(lc, "lc");
        int meshOrder = 1;
        gmshDdm->userDefinedParameter(meshOrder, "meshOrder");
        
        // ************************
        // M O D E L I N G
        // ************************
        std::string boundary = "sommerfeld"; // sommerfeld, pml, habc
        gmshDdm->userDefinedParameter(boundary, "boundary");
        std::string pmlMethod = "continuous"; // continuous, discontinuous, withoutMultiplier
        gmshDdm->userDefinedParameter(pmlMethod, "pmlMethod");
        std::string pmlType = "hs"; // hs, h, q
        gmshDdm->userDefinedParameter(pmlType, "pmlType");
        
        double N = 6;
        gmshDdm->userDefinedParameter(N, "N");
        double pmlSize = N * lc;
        if(pmlSize == 0.) {
          pmlSize = 0.3;
        }
        
        std::string gauss = "Gauss10";
        gmshDdm->userDefinedParameter(gauss, "gauss");
        int fieldOrder = 1;
        gmshDdm->userDefinedParameter(fieldOrder, "fieldOrder");
        int neumannOrder = 1;
        gmshDdm->userDefinedParameter(neumannOrder, "neumannOrder");
        double stab = 0.11;
        gmshDdm->userDefinedParameter(stab, "stab");
        double thetaPade = 0.;
        gmshDdm->userDefinedParameter(thetaPade, "thetaPade");
        thetaPade *= pi;
        
        // ************************
        // P O S T
        // ************************
        bool computeError = false;
        gmshDdm->userDefinedParameter(computeError, "error");
        bool computeInterfaceError = false;
        gmshDdm->userDefinedParameter(computeInterfaceError, "interfaceError");
        bool getMatrices = false;
        gmshDdm->userDefinedParameter(getMatrices, "getMatrices");
        bool convergence = false;
        gmshDdm->userDefinedParameter(convergence, "convergence");
        std::string fileName = "none";
        gmshDdm->userDefinedParameter(fileName, "file");
        bool saveNeumannTraces = false;
        gmshDdm->userDefinedParameter(saveNeumannTraces, "saveNeumannTraces");
        bool spy = false;
        gmshDdm->userDefinedParameter(spy, "spy");
        
        // mesh
        monoDomain(lc, pmlSize, R, meshOrder);
        // source
        gmshfem::analytics::AnalyticalFunction< gmshfem::analytics::helmholtz2D::ScatteringByASoftCylinder< std::complex< double > > > fAnalytic(k, R, 1., 1., 2 * k);
    
        gmshfem::msg::info << "Running 'monoDomain2D'" << gmshfem::msg::endl;
        gmshfem::msg::info << "Parameters:" << gmshfem::msg::endl;
        gmshfem::msg::info << " * physics:" << gmshfem::msg::endl;
        gmshfem::msg::info << "   - k: " << k << " (" << k/pi << "*pi" << ")" << gmshfem::msg::endl;
        gmshfem::msg::info << "   - R: " << R << gmshfem::msg::endl;
        gmshfem::msg::info << " * mesh:" << gmshfem::msg::endl;
        gmshfem::msg::info << "   - lc: " << lc << " (eta_h = " << (2*pi/k)/lc << ")" << gmshfem::msg::endl;
        gmshfem::msg::info << "   - meshOrder: " << meshOrder << gmshfem::msg::endl;
        gmshfem::msg::info << " * modeling:" << gmshfem::msg::endl;
        gmshfem::msg::info << "   - boundary: " << boundary << gmshfem::msg::endl;
        gmshfem::msg::info << "   - N: " << N << gmshfem::msg::endl;
        if(boundary == "pml") {
          gmshfem::msg::info << "   - pmlSize: " << pmlSize << gmshfem::msg::endl;
          gmshfem::msg::info << "   - pmlMethod: " << pmlMethod << gmshfem::msg::endl;
          gmshfem::msg::info << "   - pmlType: " << pmlType << gmshfem::msg::endl;
        }
        else if(boundary == "habc") {
          gmshfem::msg::info << "   - thetaPade: " << thetaPade << " (" << thetaPade/pi << "*pi" << ")" << gmshfem::msg::endl;
        }
        gmshfem::msg::info << "   - gauss: " << gauss << gmshfem::msg::endl;
        gmshfem::msg::info << "   - fieldOrder: " << fieldOrder << gmshfem::msg::endl;
        gmshfem::msg::info << "   - neumannOrder: " << neumannOrder << gmshfem::msg::endl;
        gmshfem::msg::info << "   - stab: " << stab << " * lc (= " << stab * lc << ")" << gmshfem::msg::endl;
        
          
        gmshfem::domain::Domain omega("omega");
        gmshfem::domain::Domain gamma("gamma");
        gmshfem::domain::Domain sigmas;
        std::vector< gmshfem::domain::Domain > pml(4);
        std::vector< gmshfem::domain::Domain > pmlCorner(4);
        std::vector< gmshfem::domain::Domain > sigma(4);
        std::vector< std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > > pmlBnd(4);
        std::vector< gmshfem::domain::Domain > corner(4);
        std::vector< std::string > dir {"E", "N", "W", "S"};
        for(unsigned int i = 0; i < 4; ++i) {
          pml[i] = gmshfem::domain::Domain("pml" + dir[i]);
          pmlCorner[i] = gmshfem::domain::Domain("pmlCorner"+ dir[i] + dir[(i+1)%4]);
          sigma[i] = gmshfem::domain::Domain("sigma" + dir[i]);
          pmlBnd[i].first = gmshfem::domain::Domain("pmlBnd" + dir[i] + "_second");
          pmlBnd[i].second = gmshfem::domain::Domain("pmlBnd" + dir[(i+1)%4] + "_first");
          corner[i] = gmshfem::domain::Domain("corner" + dir[i] + dir[(i+1)%4]);
            
          sigmas |= sigma[i];
        }
          
        gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > u("u", omega | sigmas | gamma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
        u.addConstraint(gamma, fAnalytic);
          
        gmshfem::problem::Formulation< std::complex< double > > formulation("HelmholtzMonoDomain");
    
        formulation.integral(-grad(dof(u)), grad(tf(u)), omega, gauss);
        formulation.integral(k * k * dof(u), tf(u), omega, gauss);
        
        SubproblemDomains domains;
        domains.setOmega(omega);
        domains.setSigma(sigma);
        domains.setPml(pml);
        domains.setPmlBnd(pmlBnd);
        domains.setPmlCorner(pmlCorner);
        domains.setCorner(corner);
        
        SubproblemParameters parameters;
        parameters.setGauss(gauss);
        parameters.setKappa(k);
        parameters.setNeumannOrder(neumannOrder);
        parameters.setFieldOrder(fieldOrder);
        parameters.setStab(stab * lc);
        
        std::string bnd = boundary;
        if(boundary == "habc") {
          bnd += "_" + std::to_string(N) + "_" + std::to_string(thetaPade);
        }
        else if(boundary == "pml") {
          if(pmlMethod == "continuous") {
            bnd += "Continuous_" + std::to_string(pmlSize) + "_" + pmlType;
          }
          else if(pmlMethod == "discontinuous") {
            bnd += "Discontinuous_" + std::to_string(pmlSize) + "_" + pmlType;
          }
          else if(pmlMethod == "withoutMultiplier") {
            bnd += "WithoutMultiplier_" + std::to_string(pmlSize) + "_" + pmlType;
          }
        }
        Subproblem subproblem(formulation, u.name(), domains, parameters, bnd, bnd, bnd, bnd);
        subproblem.writeFormulation();
          
        if(spy) {
          gmshfem::common::Options::instance()->dofsSortAlgorithm = gmshfem::problem::DofsSort::Algorithm::None;
        }
        formulation.pre();
        gmshfem::common::Timer assembleTime = formulation.assemble();
        gmshfem::common::Timer solveTime = formulation.solve();
              
        // POST
        
        double errorL2 = 0.;
        if(computeError) {
          const double num = std::real(gmshfem::post::integrate(pow(abs(u - fAnalytic), 2), omega, gauss));
          const double den = std::real(gmshfem::post::integrate(pow(abs(fAnalytic), 2), omega, gauss));
          errorL2 = std::sqrt(num/den);
          
          gmshfem::msg::info << "Error L2 = " << errorL2 << gmshfem::msg::endl;
          
          gmshfem::post::save(u - fAnalytic, omega, "e");
        }
    
        if(computeInterfaceError) {
          if(boundary == "pml") {
            double numInterface = 0.;
            double denInterface = 0.;
            
            for(unsigned int b = 0; b < 4; ++b) {
              const Boundary *boundary = subproblem.getBoundary(b);
              if(auto *bnd = dynamic_cast< const PmlContinuous * >(boundary)) {
                numInterface += std::real(gmshfem::post::integrate(pow(abs(u - *bnd->getUPml()), 2), sigma[b], gauss));
                denInterface += std::real(gmshfem::post::integrate(pow(abs(u), 2), sigma[b], gauss));
              }
              else if(auto *bnd = dynamic_cast< const PmlDiscontinuous * >(boundary)) {
                numInterface += std::real(gmshfem::post::integrate(pow(abs(u - *bnd->getUPml()), 2), sigma[b], gauss));
                denInterface += std::real(gmshfem::post::integrate(pow(abs(u), 2), sigma[b], gauss));
              }
              else if(auto *bnd = dynamic_cast< const PmlWithoutMultiplier * >(boundary)) {
                numInterface += std::real(gmshfem::post::integrate(pow(abs(u - *bnd->getUPml()), 2), sigma[b], gauss));
                denInterface += std::real(gmshfem::post::integrate(pow(abs(u), 2), sigma[b], gauss));
              }
            }
          
            for(unsigned int c = 0; c < 4; ++c) {
              const Corner *corner = subproblem.getCorner(c);
              if(auto *cr = dynamic_cast< const PmlContinuous_PmlContinuous * >(corner)) {
                numInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->firstBoundary()->getUPml() - *cr->getUPml()), 2), pmlBnd[c].first, gauss));
                denInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->firstBoundary()->getUPml()), 2), pmlBnd[c].first, gauss));
    
                numInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->secondBoundary()->getUPml() - *cr->getUPml()), 2), pmlBnd[c].second, gauss));
                denInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->secondBoundary()->getUPml()), 2), pmlBnd[c].second, gauss));
              }
              else if(auto *cr = dynamic_cast< const PmlDiscontinuous_PmlDiscontinuous * >(corner)) {
                numInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->firstBoundary()->getUPml() - *cr->getUPml()), 2), pmlBnd[c].first, gauss));
                denInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->firstBoundary()->getUPml()), 2), pmlBnd[c].first, gauss));
    
                numInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->secondBoundary()->getUPml() - *cr->getUPml()), 2), pmlBnd[c].second, gauss));
                denInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->secondBoundary()->getUPml()), 2), pmlBnd[c].second, gauss));
              }
              else if(auto *cr = dynamic_cast< const PmlWithoutMultiplier_PmlWithoutMultiplier * >(corner)) {
                numInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->firstBoundary()->getUPml() - *cr->getUPml()), 2), pmlBnd[c].first, gauss));
                denInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->firstBoundary()->getUPml()), 2), pmlBnd[c].first, gauss));
    
                numInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->secondBoundary()->getUPml() - *cr->getUPml()), 2), pmlBnd[c].second, gauss));
                denInterface += std::real(gmshfem::post::integrate(pow(abs(*cr->secondBoundary()->getUPml()), 2), pmlBnd[c].second, gauss));
              }
            }
            
            gmshfem::msg::info << "Interface L2 error = " << std::sqrt(numInterface/denInterface) << gmshfem::msg::endl;
            
            if(fileName != "none") {
              gmshfem::common::CSVio file(fileName, ';', gmshfem::common::OpeningMode::Append);
              file << (2*pi/k)/lc << lc << 1./lc << std::sqrt(numInterface/denInterface) << gmshfem::csv::endl;
              file.close();
            }
          }
          else {
            gmshfem::msg::error << "Unable to define a boundary error with 'boundary = " << boundary << "'" << gmshfem::msg::endl;
          }
        }
    
        if(getMatrices) {
    //      gmshfem::algebra::MatrixCRS< std::complex< double > > U;
    //      gmshfem::algebra::MatrixCRS< std::complex< double > > L;
    //      gmshfem::algebra::MatrixCRS< std::complex< double > > M;
    //      gmshfem::algebra::MatrixCRS< std::complex< double > > C;
          gmshfem::algebra::MatrixCRS< std::complex< double > > A;
    //
    //      std::vector< const gmshfem::field::FieldInterface< std::complex< double > > * > fieldsU;
    //      fieldsU.push_back(&u);
    //      for(unsigned int i = 0; i < 4; ++i) {
    //        fieldsU.push_back(subproblem.getU(i));
    //      }
    //      for(unsigned int i = 0; i < 4; ++i) {
    //        fieldsU.push_back(subproblem.getUC(i));
    //      }
    //
    //      std::vector< const gmshfem::field::FieldInterface< std::complex< double > > * > fieldsL;
    //      for(unsigned int i = 0; i < 4; ++i) {
    //        fieldsL.push_back(subproblem.getV(i));
    //      }
    //      if(boundary == "pml") {
    //        for(unsigned int i = 0; i < 4; ++i) {
    //          fieldsL.push_back(subproblem.getVC(i, 0));
    //          fieldsL.push_back(subproblem.getVC(i, 1));
    //        }
    //      }
    
      //    std::vector< const gmshfem::field::FieldInterface< std::complex< double > > * > fieldsC;
      //    if(boundary == "pml") {
      //      if(form == "form0") {
      //        for(unsigned int i = 0; i < 4; ++i) {
      //          fieldsC.push_back(static_cast< NeumannTrace< gmshfem::field::Form::Form0 > * >(neumannTrace)->get("vC_" + orientation[(i-1)%4] + "_" + orientation[i]));
      //        }
      //      }
      //    }
      //    else if (boundary == "habc") {
      //      for(unsigned int i = 0; i < 4; ++i) {
      //        for(unsigned int n = 0; n < N; ++n) {
      //          fieldsC.push_back(static_cast< NeumannTrace< gmshfem::field::Form::Form0 > * >(neumannTrace)->get("vC_" + orientation[i] + "_" + std::to_string(n) + "_first"));
      //          fieldsC.push_back(static_cast< NeumannTrace< gmshfem::field::Form::Form0 > * >(neumannTrace)->get("vC_" + orientation[(i-1)%4] + "_" + std::to_string(n) + "_second"));
      //        }
      //      }
      //    }
    
    //      formulation.getLHSBlock(U, fieldsU, fieldsU);
    //      formulation.getLHSBlock(L, fieldsL, fieldsU);
    //      formulation.getLHSBlock(C, fieldsC, fieldsL);
    //      formulation.getLHSBlock(M, fieldsL, fieldsL);
          formulation.getLHS(A);
    //
    //      U.save("U");
    //      L.save("L");
    //      C.save("C");
    //      M.save("M");
          A.save("A");
        }
        
        gmshfem::post::save(u, omega, "u");
            
        if(convergence) {
          gmshfem::common::CSVio file(fileName, ';',  gmshfem::common::OpeningMode::Append);
    //      file << lc << (2*pi/k)/lc << errorL2 << assembleTime << solveTime << formulation.getTotalNumberOfDof() << gmshfem::csv::endl;
          file << k << (2*pi/k)/lc << errorL2 << assembleTime << solveTime << gmshfem::csv::endl;
        }
        
        if(saveNeumannTraces) {
          for(unsigned int b = 0; b < 4; ++b) {
            const Boundary *boundary = subproblem.getBoundary(b);
            if(auto *bnd = dynamic_cast< const Sommerfeld * >(boundary)) {
              gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
            }
            else if(auto *bnd = dynamic_cast< const HABC * >(boundary)) {
              gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
            }
            else if(auto *bnd = dynamic_cast< const PmlContinuous * >(boundary)) {
              gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
            }
            else if(auto *bnd = dynamic_cast< const PmlDiscontinuous * >(boundary)) {
              auto n = gmshfem::function::normal< std::complex< double > >();
              gmshfem::post::save(*(bnd->getV()) * n, bnd->getV()->domain(), bnd->getV()->name());
            }
            else if(auto *bnd = dynamic_cast< const PmlWithoutMultiplier * >(boundary)) {
              gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
            }
          }
          
          for(unsigned int c = 0; c < 4; ++c) {
            const Corner *corner = subproblem.getCorner(c);
            if(auto *cr = dynamic_cast< const PmlContinuous_PmlContinuous * >(corner)) {
              gmshfem::post::save(*(cr->getV(0)), cr->getV(0)->domain(), cr->getV(0)->name());
              gmshfem::post::save(*(cr->getV(1)), cr->getV(1)->domain(), cr->getV(1)->name());
            }
            else if(auto *cr = dynamic_cast< const PmlDiscontinuous_PmlDiscontinuous * >(corner)) {
              auto n = gmshfem::function::normal< std::complex< double > >();
              gmshfem::post::save(*(cr->getV(0)) * n, cr->getV(0)->domain(), cr->getV(0)->name());
              gmshfem::post::save(*(cr->getV(1)) * n, cr->getV(1)->domain(), cr->getV(1)->name());
            }
            else if(auto *cr = dynamic_cast< const PmlWithoutMultiplier_PmlWithoutMultiplier * >(corner)) {
              gmshfem::post::save(*(cr->getV(0)), cr->getV(0)->domain(), cr->getV(0)->name());
              gmshfem::post::save(*(cr->getV(1)), cr->getV(1)->domain(), cr->getV(1)->name());
            }
          }
        }
        
        if(spy) {
          gmshfem::algebra::MatrixCRSFast< std::complex< double > > A;
          formulation.getLHS(A);
          A.saveSpyPlot("spy", 8);
        }
        
      }
    
      
    }
    
    
    namespace D3 {
    
    
      void monoDomain()
      {
        gmshddm::common::GmshDdm *gmshDdm = gmshddm::common::GmshDdm::currentInstance();
        
        // ************************
        // P H Y S I C S
        // ************************
        double pi = 3.141592653589793238462643383279;
        double k = 2. * pi;
        gmshDdm->userDefinedParameter(k, "k");
        double R = 0.5;
        gmshDdm->userDefinedParameter(R, "R");
        
        // ************************
        // M E S H
        // ************************
        double lc = 0.06666666666666; // Marmousi = 20; other = 0.06666666666666
        gmshDdm->userDefinedParameter(lc, "lc");
        int meshOrder = 1;
        gmshDdm->userDefinedParameter(meshOrder, "meshOrder");
        
        // ************************
        // M O D E L I N G
        // ************************
        std::string boundary = "sommerfeld"; // sommerfeld, pml, habc
        gmshDdm->userDefinedParameter(boundary, "boundary");
        std::string pmlMethod = "continuous"; // continuous, discontinuous
        gmshDdm->userDefinedParameter(pmlMethod, "pmlMethod");
        std::string pmlType = "hs"; // hs, h, q
        gmshDdm->userDefinedParameter(pmlType, "pmlType");
        
        double N = 6;
        gmshDdm->userDefinedParameter(N, "N");
        double pmlSize = N * lc;
        if(pmlSize == 0.) {
          pmlSize = 0.3;
        }
        
        std::string gauss = "Gauss10";
        gmshDdm->userDefinedParameter(gauss, "gauss");
        int fieldOrder = 1;
        gmshDdm->userDefinedParameter(fieldOrder, "fieldOrder");
        int neumannOrder = 1;
        gmshDdm->userDefinedParameter(neumannOrder, "neumannOrder");
        double stab = 0.11;
        gmshDdm->userDefinedParameter(stab, "stab");
        double thetaPade = 0.;
        gmshDdm->userDefinedParameter(thetaPade, "thetaPade");
        thetaPade *= pi;
        
        // ************************
        // P O S T
        // ************************
        bool onPlane = false;
        gmshDdm->userDefinedParameter(onPlane, "onPlane");
        bool computeError = false;
        gmshDdm->userDefinedParameter(computeError, "error");
        bool computeInterfaceError = false;
        gmshDdm->userDefinedParameter(computeInterfaceError, "interfaceError");
        bool getMatrices = false;
        gmshDdm->userDefinedParameter(getMatrices, "getMatrices");
        bool saveNeumannTraces = false;
        gmshDdm->userDefinedParameter(saveNeumannTraces, "saveNeumannTraces");
        
        // mesh
        monoDomain(lc, pmlSize, R, meshOrder);
        // source
        gmshfem::analytics::AnalyticalFunction< gmshfem::analytics::helmholtz3D::ScatteringByASoftSphere< std::complex< double > > > fAnalytic(k, R, 1., 1., 1., 2 * k);
        
        gmshfem::msg::info << "Running 'monoDomain3D'" << gmshfem::msg::endl;
        gmshfem::msg::info << "Parameters:" << gmshfem::msg::endl;
        gmshfem::msg::info << " * physics:" << gmshfem::msg::endl;
        gmshfem::msg::info << "   - k: " << k << " (" << k/pi << "*pi" << ")" << gmshfem::msg::endl;
        gmshfem::msg::info << "   - R: " << R << gmshfem::msg::endl;
        gmshfem::msg::info << " * mesh:" << gmshfem::msg::endl;
        gmshfem::msg::info << "   - lc: " << lc << " (eta_h = " << (2*pi/k)/lc << ")" << gmshfem::msg::endl;
        gmshfem::msg::info << "   - meshOrder: " << meshOrder << gmshfem::msg::endl;
        gmshfem::msg::info << " * modeling:" << gmshfem::msg::endl;
        gmshfem::msg::info << "   - boundary: " << boundary << gmshfem::msg::endl;
        gmshfem::msg::info << "   - N: " << N << gmshfem::msg::endl;
        if(boundary == "pml") {
          gmshfem::msg::info << "   - pmlSize: " << pmlSize << gmshfem::msg::endl;
          gmshfem::msg::info << "   - pmlMethod: " << pmlMethod << gmshfem::msg::endl;
          gmshfem::msg::info << "   - pmlType: " << pmlType << gmshfem::msg::endl;
        }
        else if(boundary == "habc") {
          gmshfem::msg::info << "   - thetaPade: " << thetaPade << " (" << thetaPade/pi << "*pi" << ")" << gmshfem::msg::endl;
        }
        gmshfem::msg::info << "   - gauss: " << gauss << gmshfem::msg::endl;
        gmshfem::msg::info << "   - fieldOrder: " << fieldOrder << gmshfem::msg::endl;
        gmshfem::msg::info << "   - neumannOrder: " << neumannOrder << gmshfem::msg::endl;
        gmshfem::msg::info << "   - stab: " << stab << " * lc (= " << stab * lc << ")" << gmshfem::msg::endl;
        
        gmshfem::domain::Domain omega("omega");
        gmshfem::domain::Domain gamma("gamma");
        gmshfem::domain::Domain sigmas;
        std::vector< gmshfem::domain::Domain > pml(6);
        std::vector< gmshfem::domain::Domain > pmlEdge(12);
        std::vector< gmshfem::domain::Domain > pmlCorner(8);
        std::vector< gmshfem::domain::Domain > sigma(6);
        std::vector< std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > > pmlBnd(12);
        std::vector< std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > > pmlEdgeBnd(8);
        std::vector< gmshfem::domain::Domain > edge(12);
        std::vector< gmshfem::domain::Domain > corner(8);
        std::vector< std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > > cornerEdge(8);
        std::vector< std::string > dir {"E", "N", "W", "S", "D", "U"};
        // face
        for(unsigned int i = 0; i < 6; ++i) {
          pml[i] = gmshfem::domain::Domain("pml" + dir[i]);
          sigma[i] = gmshfem::domain::Domain("sigma" + dir[i]);
          sigmas |= sigma[i];
        }
        // edge
        for(unsigned int i = 0; i < 12; ++i) {
          std::vector< std::string > count {"first", "second", "third", "fourth"};
          if(i < 4) {
            pmlEdge[i] = gmshfem::domain::Domain("pml" + dir[i%4] + "D");
            pmlBnd[i].first = gmshfem::domain::Domain("pmlBnd" + dir[i%4] + "_first");
            pmlBnd[i].second = gmshfem::domain::Domain("pmlBndD_" + count[i%4]);
            edge[i] = gmshfem::domain::Domain("edge" + dir[i%4] + "D");
          }
          else if(i >= 4 && i < 8) {
            pmlEdge[i] = gmshfem::domain::Domain("pml" + dir[i%4] + dir[(i+1)%4]);
            pmlBnd[i].first = gmshfem::domain::Domain("pmlBnd" + dir[i%4] + "_second");
            pmlBnd[i].second = gmshfem::domain::Domain("pmlBnd" + dir[(i+1)%4] + "_fourth");
            edge[i] = gmshfem::domain::Domain("edge" + dir[i%4] + dir[(i+1)%4]);
          }
          else {
            pmlEdge[i] = gmshfem::domain::Domain("pml" + dir[i%4] + "U");
            pmlBnd[i].first = gmshfem::domain::Domain("pmlBnd" + dir[i%4] + "_third");
            pmlBnd[i].second = gmshfem::domain::Domain("pmlBndU_" + count[i%4]);
            edge[i] = gmshfem::domain::Domain("edge" + dir[i%4] + "U");
          }
        }
        // corner
        for(unsigned int i = 0; i < 8; ++i) {
          if(i < 4) {
            pmlCorner[i] = gmshfem::domain::Domain("pml" + dir[i%4] + dir[(i+1)%4] + "D");
            std::get<0>(pmlEdgeBnd[i]) = gmshfem::domain::Domain("pmlBnd" + dir[i%4] + "D_second");
            std::get<1>(pmlEdgeBnd[i]) = gmshfem::domain::Domain("pmlBnd" + dir[(i+1)%4] + "D_first");
            std::get<2>(pmlEdgeBnd[i]) = gmshfem::domain::Domain("pmlBnd" + dir[i%4] + dir[(i+1)%4] + "_first");
            corner[i] = gmshfem::domain::Domain("corner" + dir[i%4] + dir[(i+1)%4] + "D");
            std::get<0>(cornerEdge[i]) = gmshfem::domain::Domain("cornerEdge" + dir[i%4] + dir[(i+1)%4] + "D_first");
            std::get<1>(cornerEdge[i]) = gmshfem::domain::Domain("cornerEdge" + dir[i%4] + dir[(i+1)%4] + "D_second");
            std::get<2>(cornerEdge[i]) = gmshfem::domain::Domain("cornerEdge" + dir[i%4] + dir[(i+1)%4] + "D_third");
          }
          else {
            pmlCorner[i] = gmshfem::domain::Domain("pml" + dir[i%4] + dir[(i+1)%4] + "U");
            std::get<0>(pmlEdgeBnd[i]) = gmshfem::domain::Domain("pmlBnd" + dir[i%4] + "U_second");
            std::get<1>(pmlEdgeBnd[i]) = gmshfem::domain::Domain("pmlBnd" + dir[(i+1)%4] + "U_first");
            std::get<2>(pmlEdgeBnd[i]) = gmshfem::domain::Domain("pmlBnd" + dir[i%4] + dir[(i+1)%4] + "_second");
            corner[i] = gmshfem::domain::Domain("corner" + dir[i%4] + dir[(i+1)%4] + "U");
            std::get<0>(cornerEdge[i]) = gmshfem::domain::Domain("cornerEdge" + dir[i%4] + dir[(i+1)%4] + "U_first");
            std::get<1>(cornerEdge[i]) = gmshfem::domain::Domain("cornerEdge" + dir[i%4] + dir[(i+1)%4] + "U_second");
            std::get<2>(cornerEdge[i]) = gmshfem::domain::Domain("cornerEdge" + dir[i%4] + dir[(i+1)%4] + "U_third");
          }
        }
          
        gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > u("u", omega | sigmas | gamma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
        u.addConstraint(gamma, fAnalytic);
          
        gmshfem::problem::Formulation< std::complex< double > > formulation("HelmholtzMonoDomain");
    
        formulation.integral(-grad(dof(u)), grad(tf(u)), omega, gauss);
        formulation.integral(k * k * dof(u), tf(u), omega, gauss);
        
        SubproblemDomains domains;
        domains.setOmega(omega);
        domains.setSigma(sigma);
        domains.setPml(pml);
        domains.setPmlBnd(pmlBnd);
        domains.setPmlEdge(pmlEdge);
        domains.setEdge(edge);
        domains.setPmlEdgeBnd(pmlEdgeBnd);
        domains.setPmlCorner(pmlCorner);
        domains.setCorner(corner);
        domains.setCornerEdge(cornerEdge);
    
        SubproblemParameters parameters;
        parameters.setGauss(gauss);
        parameters.setKappa(k);
        parameters.setNeumannOrder(neumannOrder);
        parameters.setFieldOrder(fieldOrder);
        parameters.setStab(stab * lc);
        
        std::string bnd = boundary;
        if(boundary == "habc") {
          bnd += "_" + std::to_string(N) + "_" + std::to_string(thetaPade);
        }
        else if(boundary == "pml") {
          if(pmlMethod == "continuous") {
            bnd += "Continuous_" + std::to_string(pmlSize) + "_" + pmlType;
          }
          else if(pmlMethod == "discontinuous") {
            bnd += "Discontinuous_" + std::to_string(pmlSize) + "_" + pmlType;
          }
        }
        Subproblem subproblem(formulation, u.name(), domains, parameters, bnd, bnd, bnd, bnd, bnd, bnd);
        subproblem.writeFormulation();
        
        formulation.pre();
        formulation.assemble();
        formulation.solve();
    
        // POST
        const double eps = 0.01;
        gmshfem::post::Plane planeXY("planeXY", -1.+eps, -1.+eps, 0., 2.-2.*eps, 0., 0., 0., 2.-2.*eps, 0., 100., 100.);
        gmshfem::post::Plane planeYZ("planeYZ", -1.+eps, -1.+eps, -1.+eps, 2.-2.*eps, 2.-2.*eps, 0., 0., 0., 2.-2.*eps, 100., 100.);
        gmshfem::post::Plane planeZX("planeZX", -1.+eps, 1.-eps, -1.+eps, 2.-2.*eps, -2.+2.*eps, 0., 0., 0., 2.-2.*eps, 100., 100.);
        gmshfem::post::Sphere scatter("scatter", 0., 0., 0., R+eps, 100. * R);
        
        if(computeError) {
          const double num = std::real(gmshfem::post::integrate(pow(abs(u - fAnalytic), 2), omega, gauss));
          const double den = std::real(gmshfem::post::integrate(pow(abs(fAnalytic), 2), omega, gauss));
    
          gmshfem::msg::info << "Error L2 = " << std::sqrt(num/den) << gmshfem::msg::endl;
          
          if(onPlane) {
            gmshfem::post::save(gmshfem::function::heaviside(gmshfem::function::r3d< std::complex< double > >() - R) * (u - fAnalytic), planeXY, "e_XY", "pos");
            gmshfem::post::save(gmshfem::function::heaviside(gmshfem::function::r3d< std::complex< double > >() - R) * (u - fAnalytic), planeYZ, "e_YZ", "pos");
            gmshfem::post::save(gmshfem::function::heaviside(gmshfem::function::r3d< std::complex< double > >() - R) * (u - fAnalytic), planeZX, "e_ZX", "pos");
            gmshfem::post::save(u - fAnalytic, scatter, "e_SCAT", "pos");
          }
          else {
            gmshfem::post::save(gmshfem::function::norm(u - fAnalytic), omega, "e");
          }
        }
    
        if(onPlane) {
          gmshfem::post::save(u, planeXY, "u_XY", "pos");
          gmshfem::post::save(u, planeYZ, "u_YZ", "pos");
          gmshfem::post::save(u, planeZX, "u_ZX", "pos");
          gmshfem::post::save(u, scatter, "u_SCAT", "pos");
        }
        else {
          gmshfem::post::save(u, omega, "u");
        }
        
        if(getMatrices) {
          gmshfem::algebra::MatrixCRS< std::complex< double > > A;
          formulation.getLHS(A);
          A.save("A");
        }
        
        if(saveNeumannTraces) {
          for(unsigned int b = 0; b < 6; ++b) {
            const Boundary *boundary = subproblem.getBoundary(b);
            if(auto *bnd = dynamic_cast< const Sommerfeld * >(boundary)) {
              gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
            }
            else if(auto *bnd = dynamic_cast< const HABC * >(boundary)) {
              gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
            }
            else if(auto *bnd = dynamic_cast< const PmlContinuous * >(boundary)) {
              gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
            }
            else if(auto *bnd = dynamic_cast< const PmlDiscontinuous * >(boundary)) {
              auto n = gmshfem::function::normal< std::complex< double > >();
              gmshfem::post::save(*(bnd->getV()) * n, bnd->getV()->domain(), bnd->getV()->name());
            }
            else if(auto *bnd = dynamic_cast< const PmlWithoutMultiplier * >(boundary)) {
              gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
            }
          }
          
          for(unsigned int e = 0; e < 12; ++e) {
            const Edge *edge = subproblem.getEdge(e);
            if(auto *eg = dynamic_cast< const PmlContinuous_PmlContinuous * >(edge)) {
              gmshfem::post::save(*(eg->getV(0)), eg->getV(0)->domain(), eg->getV(0)->name());
              gmshfem::post::save(*(eg->getV(1)), eg->getV(1)->domain(), eg->getV(1)->name());
            }
    //        else if(auto *eg = dynamic_cast< const PmlDiscontinuous_PmlDiscontinuous * >(edge)) {
    //          auto n = gmshfem::function::normal< std::complex< double > >();
    //          gmshfem::post::save(*(eg->getV(0)) * n, eg->getV(0)->domain(), eg->getV(0)->name());
    //          gmshfem::post::save(*(eg->getV(1)) * n, eg->getV(1)->domain(), eg->getV(1)->name());
    //        }
    //        else if(auto *eg = dynamic_cast< const PmlWithoutMultiplier_PmlWithoutMultiplier * >(boundary)) {
    //          gmshfem::post::save(*(bnd->getV()), bnd->getV()->domain(), bnd->getV()->name());
    //        }
          }
          
          for(unsigned int c = 0; c < 8; ++c) {
            const Corner *corner = subproblem.getCorner(c);
            if(auto *cr = dynamic_cast< const PmlContinuous_PmlContinuous_PmlContinuous * >(corner)) {
              gmshfem::post::save(*(cr->getV(0)), cr->getV(0)->domain(), cr->getV(0)->name());
              gmshfem::post::save(*(cr->getV(1)), cr->getV(1)->domain(), cr->getV(1)->name());
              gmshfem::post::save(*(cr->getV(2)), cr->getV(2)->domain(), cr->getV(2)->name());
            }
    //        else if(auto *cr = dynamic_cast< const PmlDiscontinuous_PmlDiscontinuous_PmlDiscontinuous * >(corner)) {
    //          auto n = gmshfem::function::normal< std::complex< double > >();
    //          gmshfem::post::save(*(cr->getV(0)) * n, cr->getV(0)->domain(), cr->getV(0)->name());
    //          gmshfem::post::save(*(cr->getV(1)) * n, cr->getV(1)->domain(), cr->getV(1)->name());
    //          gmshfem::post::save(*(cr->getV(2)) * n, cr->getV(2)->domain(), cr->getV(2)->name());
    //        }
    //        else if(auto *cr = dynamic_cast< const PmlWithoutMultiplier_PmlWithoutMultiplier_PmlWithoutMultiplier * >(corner)) {
    //          gmshfem::post::save(*(cr->getV(0)), cr->getV(0)->domain(), cr->getV(0)->name());
    //          gmshfem::post::save(*(cr->getV(1)), cr->getV(1)->domain(), cr->getV(1)->name());
    //        }
          }
        }
      }
    
    
    }