Skip to content
Snippets Groups Projects
Select Git revision
  • 177f6917e1590a3c2103786eb3408e01b2139ffb
  • master default protected
  • alphashapes
  • quadMeshingTools
  • cygwin_conv_path
  • macos_arm64
  • add-transfiniteautomatic-to-geo
  • patch_releases_4_10
  • HierarchicalHDiv
  • isuruf-master-patch-63355
  • hyperbolic
  • hexdom
  • hxt_update
  • jf
  • 1618-pythonocc-and-gmsh-api-integration
  • octreeSizeField
  • hexbl
  • alignIrregularVertices
  • getEdges
  • patch_releases_4_8
  • isuruf-master-patch-51992
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
  • gmsh_4_8_4
  • gmsh_4_8_3
  • gmsh_4_8_2
  • gmsh_4_8_1
  • gmsh_4_8_0
  • gmsh_4_7_1
  • gmsh_4_7_0
41 results

delaunay_refinement.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    ddm3D.cpp 102.54 KiB
    #include "ddm3D.h"
    
    #include "mesh.h"
    #include "SubproblemDomains.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>
    #include <vector>
    
    using gmshfem::equation::dof;
    using gmshfem::equation::tf;
    using gmshfem::function::operator-;
    using gmshfem::function::operator*;
    using gmshfem::function::abs;
    
    namespace D3 {
    
    
      void ddm()
      {
        gmshddm::common::GmshDdm *gmshDdm = gmshddm::common::GmshDdm::currentInstance();
    
        // ************************
        // P H Y S I C S
        // ************************
        double pi = 3.14159265359;
        double k = 2. * pi;
        gmshDdm->userDefinedParameter(k, "k");
        double R = 0.5;
        gmshDdm->userDefinedParameter(R, "R");
        std::string benchmark = "scattering"; // scattering
        gmshDdm->userDefinedParameter(benchmark, "benchmark");
    
        // ************************
        // M E S H
        // ************************
        double lc = 0.1; // 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 boundaryExt = "sommerfeld"; // sommerfeld, pml, habc
        gmshDdm->userDefinedParameter(boundaryExt, "boundaryExt");
        std::string pmlMethod = "continuous"; // continuous, discontinuous
        gmshDdm->userDefinedParameter(pmlMethod, "pmlMethod");
        std::string pmlMethodExt = "continuous"; // continuous, discontinuous
        gmshDdm->userDefinedParameter(pmlMethodExt, "pmlMethodExt");
        std::string pmlType = "hs"; // hs, h, q
        gmshDdm->userDefinedParameter(pmlType, "pmlType");
        std::string pmlTypeExt = "hs"; // hs, h, q
        gmshDdm->userDefinedParameter(pmlTypeExt, "pmlTypeExt");
    
    
        unsigned int nDomX = 2, nDomY = 2, nDomZ = 2;
        gmshDdm->userDefinedParameter(nDomX, "nDomX");
        gmshDdm->userDefinedParameter(nDomY, "nDomY");
        gmshDdm->userDefinedParameter(nDomZ, "nDomZ");
        unsigned int iterMax = 1000;
        gmshDdm->userDefinedParameter(iterMax, "iterMax");
        double res = 1e-6;
        gmshDdm->userDefinedParameter(res, "res");
        double sizeX = 2., sizeY = 2., sizeZ = 2.;
        gmshDdm->userDefinedParameter(sizeX, "sizeX");
        gmshDdm->userDefinedParameter(sizeY, "sizeY");
        gmshDdm->userDefinedParameter(sizeZ, "sizeZ");
    
        unsigned int N = 6;
        gmshDdm->userDefinedParameter(N, "N");
        double pmlSize = N * lc;
        if(pmlSize == 0.) {
          pmlSize = 0.3;
        }
    
        unsigned int NExt = N;
        gmshDdm->userDefinedParameter(NExt, "NExt");
        double pmlSizeExt = NExt * lc;
        if(pmlSizeExt == 0.) {
          pmlSizeExt = 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 monoDomainError = false;
        gmshDdm->userDefinedParameter(monoDomainError, "monoDomainError");
        std::string fileName = "none";
        gmshDdm->userDefinedParameter(fileName, "file");
        bool saveEMono = false;
        gmshDdm->userDefinedParameter(saveEMono, "saveEMono");
        bool saveUMono = false;
        gmshDdm->userDefinedParameter(saveUMono, "saveUMono");
        bool saveU = false;
        gmshDdm->userDefinedParameter(saveU, "saveU");
        bool wavenumberPlot = false;
        gmshDdm->userDefinedParameter(wavenumberPlot, "wavenumber");
        bool meshPlot = false;
        gmshDdm->userDefinedParameter(meshPlot, "mesh");
        bool saveMesh = false;
        gmshDdm->userDefinedParameter(saveMesh, "saveMesh");
    
        if(benchmark == "scattering") {
          checkerboard(nDomX, nDomY, nDomZ, sizeX, sizeY, sizeZ, R, lc, (boundary == "pml"), pmlSize, (boundaryExt == "pml"), pmlSizeExt, meshOrder, 0, 0, 0);
        }
        if(saveMesh) {
          gmsh::write("mesh.msh");
        }
    
        gmshfem::msg::info << "Running 'ddm3D'" << 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 << "   - grid: (" << nDomX << ", " << nDomY << ", " << nDomZ << ")" << gmshfem::msg::endl;
        gmshfem::msg::info << "   - iterMax: " << iterMax << gmshfem::msg::endl;
        gmshfem::msg::info << "   - res: " << res << gmshfem::msg::endl;
        gmshfem::msg::info << "   - subdomain size: (" << sizeX << ", " << sizeY << ", " << sizeZ << ")" << gmshfem::msg::endl;
        gmshfem::msg::info << "   - boundary: " << boundary << gmshfem::msg::endl;
        gmshfem::msg::info << "   - boundaryExt: " << boundaryExt << gmshfem::msg::endl;
        if(boundary == "pml") {
          gmshfem::msg::info << "   - N: " << N << gmshfem::msg::endl;
          gmshfem::msg::info << "   - pmlSize: " << pmlSize << gmshfem::msg::endl;
          gmshfem::msg::info << "   - pmlType: " << pmlType << gmshfem::msg::endl;
          gmshfem::msg::info << "   - pmlMethod: " << pmlMethod << gmshfem::msg::endl;
        }
        else if(boundary == "habc") {
          gmshfem::msg::info << "   - N: " << N << gmshfem::msg::endl;
          gmshfem::msg::info << "   - thetaPade: " << thetaPade << " (" << thetaPade/pi << "*pi" << ")" << gmshfem::msg::endl;
        }
        if(boundaryExt == "pml") {
          gmshfem::msg::info << "   - NExt: " << NExt << gmshfem::msg::endl;
          gmshfem::msg::info << "   - pmlSizeExt: " << pmlSizeExt << gmshfem::msg::endl;
          gmshfem::msg::info << "   - pmlTypeExt: " << pmlTypeExt << gmshfem::msg::endl;
          gmshfem::msg::info << "   - pmlMethodExt: " << pmlMethodExt << gmshfem::msg::endl;
        }
        else if(boundaryExt == "habc") {
          gmshfem::msg::info << "   - NExt: " << NExt << 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 << gmshfem::msg::endl;
    
    
        // source
        gmshfem::analytics::AnalyticalFunction< gmshfem::analytics::helmholtz3D::ScatteringByASoftSphere< std::complex< double > > > fAnalytic(k, R, sizeX/2., sizeY/2., sizeZ/2., 2 * k);
    
        const unsigned int nDom = nDomX * nDomY * nDomZ;
        // Define domain
        gmshddm::domain::Subdomain omega(nDom);
        gmshfem::domain::Domain gamma;
        if(benchmark == "scattering") {
          gamma = gmshfem::domain::Domain("gamma");
        }
    
        std::vector< gmshddm::domain::Subdomain > pml(6, nDom);
        std::vector< gmshddm::domain::Subdomain > pmlEdge(12, nDom);
        std::vector< gmshddm::domain::Subdomain > pmlCorner(8, nDom);
        std::vector< gmshddm::domain::Subdomain > sigma(6, nDom);
        std::vector< gmshddm::domain::Interface > sigmaInterface(6, nDom);
        std::vector< std::pair< gmshddm::domain::Subdomain, gmshddm::domain::Subdomain > > pmlBnd(12, std::make_pair(gmshddm::domain::Subdomain(nDom), gmshddm::domain::Subdomain(nDom)) );
        std::vector< std::pair< gmshddm::domain::Interface, gmshddm::domain::Interface > > pmlBndInterface(12, std::make_pair(gmshddm::domain::Interface(nDom), gmshddm::domain::Interface(nDom)) );
        std::vector< std::tuple< gmshddm::domain::Subdomain, gmshddm::domain::Subdomain, gmshddm::domain::Subdomain > > pmlEdgeBnd(8, std::make_tuple(gmshddm::domain::Subdomain(nDom), gmshddm::domain::Subdomain(nDom), gmshddm::domain::Subdomain(nDom)) );
        std::vector< std::tuple< gmshddm::domain::Interface, gmshddm::domain::Interface, gmshddm::domain::Interface > > pmlEdgeBndInterface(8, std::make_tuple(gmshddm::domain::Interface(nDom), gmshddm::domain::Interface(nDom), gmshddm::domain::Interface(nDom)) );
    
        std::vector< gmshddm::domain::Subdomain > edge(12, nDom);
        std::vector< std::pair< gmshddm::domain::Interface, gmshddm::domain::Interface > > edgeInterface(12, std::make_pair(gmshddm::domain::Interface(nDom), gmshddm::domain::Interface(nDom)) );
        std::vector< gmshddm::domain::Subdomain > corner(8, nDom);
        std::vector< std::tuple< gmshddm::domain::Interface, gmshddm::domain::Interface, gmshddm::domain::Interface > > cornerInterface(8, std::make_tuple(gmshddm::domain::Interface(nDom), gmshddm::domain::Interface(nDom), gmshddm::domain::Interface(nDom)) );
        std::vector< std::tuple< gmshddm::domain::Subdomain, gmshddm::domain::Subdomain, gmshddm::domain::Subdomain > > cornerEdge(8, std::make_tuple(gmshddm::domain::Subdomain(nDom), gmshddm::domain::Subdomain(nDom), gmshddm::domain::Subdomain(nDom)));
        // Define topology
        std::vector< std::vector< unsigned int > > topology(nDom);
        std::vector< std::string > dir {"E", "N", "W", "S", "D", "U"};
        for(unsigned int i = 0; i < static_cast< unsigned int >(nDomX); ++i) {
          for(unsigned int j = 0; j < static_cast< unsigned int >(nDomY); ++j) {
            for(unsigned int k = 0; k < static_cast< unsigned int >(nDomZ); ++k) {
              unsigned int index = i * nDomZ * nDomY + j * nDomZ + k;
    
              const std::string subdomainTag = "_" + std::to_string(i) + "_" + std::to_string(j) + "_" + std::to_string(k);
              omega(index) = gmshfem::domain::Domain("omega" + subdomainTag);
    
              for(unsigned int b = 0; b < 6; ++b) {
                if(boundary == "pml" || boundaryExt == "pml") {
                  pml[b](index) = gmshfem::domain::Domain("pml" + dir[b] + subdomainTag);
                }
                sigma[b](index) = gmshfem::domain::Domain("sigma" + dir[b] + subdomainTag);
              }
    
              for(unsigned int e = 0; e < 12; ++e) {
                std::vector< std::string > count {"first", "second", "third", "fourth"};
                if(e < 4) {
                  if(boundary == "pml" || boundaryExt == "pml") {
                    pmlEdge[e](index) = gmshfem::domain::Domain("pml" + dir[e%4] + "D" + subdomainTag);
                    pmlBnd[e].first(index) = gmshfem::domain::Domain("pmlBnd" + dir[e%4] + "_first" + subdomainTag);
                    pmlBnd[e].second(index) = gmshfem::domain::Domain("pmlBndD_" + count[e%4] + subdomainTag);
                  }
                  edge[e](index) = gmshfem::domain::Domain("edge" + dir[e%4] + "D" + subdomainTag);
                }
                else if(e >= 4 && e < 8) {
                  if(boundary == "pml" || boundaryExt == "pml") {
                    pmlEdge[e](index) = gmshfem::domain::Domain("pml" + dir[e%4] + dir[(e+1)%4] + subdomainTag);
                    pmlBnd[e].first(index) = gmshfem::domain::Domain("pmlBnd" + dir[e%4] + "_second" + subdomainTag);
                    pmlBnd[e].second(index) = gmshfem::domain::Domain("pmlBnd" + dir[(e+1)%4] + "_fourth" + subdomainTag);
                  }
                  edge[e](index) = gmshfem::domain::Domain("edge" + dir[e%4] + dir[(e+1)%4] + subdomainTag);
                }
                else {
                  if(boundary == "pml" || boundaryExt == "pml") {
                    pmlEdge[e](index) = gmshfem::domain::Domain("pml" + dir[e%4] + "U" + subdomainTag);
                    pmlBnd[e].first(index) = gmshfem::domain::Domain("pmlBnd" + dir[e%4] + "_third" + subdomainTag);
                    pmlBnd[e].second(index) = gmshfem::domain::Domain("pmlBndU_" + count[e%4] + subdomainTag);
                  }
                  edge[e](index) = gmshfem::domain::Domain("edge" + dir[e%4] + "U" + subdomainTag);
                }
              }
              
              for(unsigned int c = 0; c < 8; ++c) {
                if(c < 4) {
                  if(boundary == "pml" || boundaryExt == "pml") {
                    pmlCorner[c](index) = gmshfem::domain::Domain("pml" + dir[c%4] + dir[(c+1)%4] + "D" + subdomainTag);
                    
                    std::get<0>(pmlEdgeBnd[c])(index) = gmshfem::domain::Domain("pmlBnd" + dir[c%4] + "D_second" + subdomainTag);
                    std::get<1>(pmlEdgeBnd[c])(index) = gmshfem::domain::Domain("pmlBnd" + dir[(c+1)%4] + "D_first" + subdomainTag);
                    std::get<2>(pmlEdgeBnd[c])(index) = gmshfem::domain::Domain("pmlBnd" + dir[c%4] + dir[(c+1)%4] + "_first" + subdomainTag);
                    
                    std::get<0>(cornerEdge[c])(index) = gmshfem::domain::Domain("cornerEdge" + dir[c%4] + dir[(c+1)%4] + "D_first" + subdomainTag);
                    std::get<1>(cornerEdge[c])(index) = gmshfem::domain::Domain("cornerEdge" + dir[c%4] + dir[(c+1)%4] + "D_second" + subdomainTag);
                    std::get<2>(cornerEdge[c])(index) = gmshfem::domain::Domain("cornerEdge" + dir[c%4] + dir[(c+1)%4] + "D_third" + subdomainTag);
                  }
                  corner[c](index) = gmshfem::domain::Domain("corner" + dir[c%4] + dir[(c+1)%4] + "D" + subdomainTag);
                }
                else {
                  if(boundary == "pml" || boundaryExt == "pml") {
                    pmlCorner[c](index) = gmshfem::domain::Domain("pml" + dir[c%4] + dir[(c+1)%4] + "U" + subdomainTag);
                    
                    std::get<0>(pmlEdgeBnd[c])(index) = gmshfem::domain::Domain("pmlBnd" + dir[c%4] + "U_second" + subdomainTag);
                    std::get<1>(pmlEdgeBnd[c])(index) = gmshfem::domain::Domain("pmlBnd" + dir[(c+1)%4] + "U_first" + subdomainTag);
                    std::get<2>(pmlEdgeBnd[c])(index) = gmshfem::domain::Domain("pmlBnd" + dir[c%4] + dir[(c+1)%4] + "_second" + subdomainTag);
                    
                    std::get<0>(cornerEdge[c])(index) = gmshfem::domain::Domain("cornerEdge" + dir[c%4] + dir[(c+1)%4] + "U_first" + subdomainTag);
                    std::get<1>(cornerEdge[c])(index) = gmshfem::domain::Domain("cornerEdge" + dir[c%4] + dir[(c+1)%4] + "U_second" + subdomainTag);
                    std::get<2>(cornerEdge[c])(index) = gmshfem::domain::Domain("cornerEdge" + dir[c%4] + dir[(c+1)%4] + "U_third" + subdomainTag);
                  }
                  corner[c](index) = gmshfem::domain::Domain("corner" + dir[c%4] + dir[(c+1)%4] + "U" + subdomainTag);
                }
              }
        
              if(i != static_cast< unsigned int >(nDomX) - 1) { // E
                if(boundary == "pml" || boundaryExt == "pml") {
                  pmlBndInterface[0].second(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = pmlBnd[0].second(index);
                  pmlBndInterface[4].second(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = pmlBnd[4].second(index);
                  pmlBndInterface[8].second(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = pmlBnd[8].second(index);
                  pmlBndInterface[7].first(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = pmlBnd[7].first(index);
                  
                  std::get<1>(pmlEdgeBndInterface[0])(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = std::get<1>(pmlEdgeBnd[0])(index);
                  std::get<1>(pmlEdgeBndInterface[4])(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = std::get<1>(pmlEdgeBnd[4])(index);
                  std::get<0>(pmlEdgeBndInterface[7])(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = std::get<0>(pmlEdgeBnd[7])(index);
                  std::get<0>(pmlEdgeBndInterface[3])(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = std::get<0>(pmlEdgeBnd[3])(index);
    
                }
                sigmaInterface[0](index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = sigma[0](index);
                
                edgeInterface[0].second(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = edge[0](index);
                edgeInterface[4].second(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = edge[4](index);
                edgeInterface[8].second(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = edge[8](index);
                edgeInterface[7].first(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = edge[7](index);
                
                std::get<1>(cornerInterface[0])(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = corner[0](index);
                std::get<1>(cornerInterface[4])(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = corner[4](index);
                std::get<0>(cornerInterface[7])(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = corner[7](index);
                std::get<0>(cornerInterface[3])(index, (i+1) * nDomY * nDomZ + j * nDomZ + k) = corner[3](index);
    
                topology[index].push_back((i+1) * nDomY * nDomZ + j * nDomZ + k);
              }
    
              if(j != static_cast< unsigned int >(nDomY) - 1) { // N
                if(boundary == "pml" || boundaryExt == "pml") {
                  pmlBndInterface[1].second(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = pmlBnd[1].second(index);
                  pmlBndInterface[5].second(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = pmlBnd[5].second(index);
                  pmlBndInterface[9].second(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = pmlBnd[9].second(index);
                  pmlBndInterface[4].first(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = pmlBnd[4].first(index);
                  
                  std::get<1>(pmlEdgeBndInterface[1])(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = std::get<1>(pmlEdgeBnd[1])(index);
                  std::get<1>(pmlEdgeBndInterface[5])(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = std::get<1>(pmlEdgeBnd[5])(index);
                  std::get<0>(pmlEdgeBndInterface[4])(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = std::get<0>(pmlEdgeBnd[4])(index);
                  std::get<0>(pmlEdgeBndInterface[0])(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = std::get<0>(pmlEdgeBnd[0])(index);
                }
                sigmaInterface[1](index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = sigma[1](index);
                
                edgeInterface[1].second(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = edge[1](index);
                edgeInterface[5].second(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = edge[5](index);
                edgeInterface[9].second(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = edge[9](index);
                edgeInterface[4].first(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = edge[4](index);
                
                std::get<1>(cornerInterface[1])(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = corner[1](index);
                std::get<1>(cornerInterface[5])(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = corner[5](index);
                std::get<0>(cornerInterface[4])(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = corner[4](index);
                std::get<0>(cornerInterface[0])(index, i * nDomY * nDomZ + (j+1) * nDomZ + k) = corner[0](index);
    
                topology[index].push_back(i * nDomY * nDomZ + (j+1) * nDomZ + k);
              }
    
              if(i != 0) { // W
                if(boundary == "pml" || boundaryExt == "pml") {
                  pmlBndInterface[2].second(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = pmlBnd[2].second(index);
                  pmlBndInterface[6].second(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = pmlBnd[6].second(index);
                  pmlBndInterface[10].second(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = pmlBnd[10].second(index);
                  pmlBndInterface[5].first(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = pmlBnd[5].first(index);
                  
                  std::get<1>(pmlEdgeBndInterface[2])(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = std::get<1>(pmlEdgeBnd[2])(index);
                  std::get<1>(pmlEdgeBndInterface[6])(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = std::get<1>(pmlEdgeBnd[6])(index);
                  std::get<0>(pmlEdgeBndInterface[5])(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = std::get<0>(pmlEdgeBnd[5])(index);
                  std::get<0>(pmlEdgeBndInterface[1])(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = std::get<0>(pmlEdgeBnd[1])(index);
                }
                sigmaInterface[2](index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = sigma[2](index);
                
                edgeInterface[2].second(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = edge[2](index);
                edgeInterface[6].second(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = edge[6](index);
                edgeInterface[10].second(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = edge[10](index);
                edgeInterface[5].first(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = edge[5](index);
                
                std::get<1>(cornerInterface[2])(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = corner[2](index);
                std::get<1>(cornerInterface[6])(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = corner[6](index);
                std::get<0>(cornerInterface[5])(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = corner[5](index);
                std::get<0>(cornerInterface[1])(index, (i-1) * nDomY * nDomZ + j * nDomZ + k) = corner[1](index);
    
                topology[index].push_back((i-1) * nDomY * nDomZ + j * nDomZ + k);
              }
    
              if(j != 0) { // S
                if(boundary == "pml" || boundaryExt == "pml") {
                  pmlBndInterface[3].second(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = pmlBnd[3].second(index);
                  pmlBndInterface[7].second(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = pmlBnd[7].second(index);
                  pmlBndInterface[11].second(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = pmlBnd[11].second(index);
                  pmlBndInterface[6].first(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = pmlBnd[6].first(index);
                  
                  std::get<1>(pmlEdgeBndInterface[3])(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = std::get<1>(pmlEdgeBnd[3])(index);
                  std::get<1>(pmlEdgeBndInterface[7])(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = std::get<1>(pmlEdgeBnd[7])(index);
                  std::get<0>(pmlEdgeBndInterface[6])(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = std::get<0>(pmlEdgeBnd[6])(index);
                  std::get<0>(pmlEdgeBndInterface[2])(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = std::get<0>(pmlEdgeBnd[2])(index);
                }
                sigmaInterface[3](index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = sigma[3](index);
                
                edgeInterface[3].second(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = edge[3](index);
                edgeInterface[7].second(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = edge[7](index);
                edgeInterface[11].second(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = edge[11](index);
                edgeInterface[6].first(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = edge[6](index);
                
                std::get<1>(cornerInterface[3])(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = corner[3](index);
                std::get<1>(cornerInterface[7])(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = corner[7](index);
                std::get<0>(cornerInterface[6])(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = corner[6](index);
                std::get<0>(cornerInterface[2])(index, i * nDomY * nDomZ + (j-1) * nDomZ + k) = corner[2](index);
    
                topology[index].push_back(i * nDomY * nDomZ + (j-1) * nDomZ + k);
              }
              
              if(k != 0) { // D
                if(boundary == "pml" || boundaryExt == "pml") {
                  pmlBndInterface[0].first(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = pmlBnd[0].first(index);
                  pmlBndInterface[1].first(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = pmlBnd[1].first(index);
                  pmlBndInterface[2].first(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = pmlBnd[2].first(index);
                  pmlBndInterface[3].first(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = pmlBnd[3].first(index);
                  
                  std::get<2>(pmlEdgeBndInterface[0])(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = std::get<2>(pmlEdgeBnd[0])(index);
                  std::get<2>(pmlEdgeBndInterface[1])(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = std::get<2>(pmlEdgeBnd[1])(index);
                  std::get<2>(pmlEdgeBndInterface[2])(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = std::get<2>(pmlEdgeBnd[2])(index);
                  std::get<2>(pmlEdgeBndInterface[3])(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = std::get<2>(pmlEdgeBnd[3])(index);
                }
                sigmaInterface[4](index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = sigma[4](index);
          
                edgeInterface[0].first(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = edge[0](index);
                edgeInterface[1].first(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = edge[1](index);
                edgeInterface[2].first(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = edge[2](index);
                edgeInterface[3].first(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = edge[3](index);
                
                std::get<2>(cornerInterface[0])(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = corner[0](index);
                std::get<2>(cornerInterface[1])(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = corner[1](index);
                std::get<2>(cornerInterface[2])(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = corner[2](index);
                std::get<2>(cornerInterface[3])(index, i * nDomY * nDomZ + j * nDomZ + (k-1)) = corner[3](index);
    
                topology[index].push_back(i * nDomY * nDomZ + j * nDomZ + (k-1));
              }
              
              if(k != static_cast< unsigned int >(nDomZ) - 1) { // U
                if(boundary == "pml" || boundaryExt == "pml") {
                  pmlBndInterface[8].first(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = pmlBnd[8].first(index);
                  pmlBndInterface[9].first(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = pmlBnd[9].first(index);
                  pmlBndInterface[10].first(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = pmlBnd[10].first(index);
                  pmlBndInterface[11].first(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = pmlBnd[11].first(index);
                  
                  std::get<2>(pmlEdgeBndInterface[4])(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = std::get<2>(pmlEdgeBnd[4])(index);
                  std::get<2>(pmlEdgeBndInterface[5])(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = std::get<2>(pmlEdgeBnd[5])(index);
                  std::get<2>(pmlEdgeBndInterface[6])(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = std::get<2>(pmlEdgeBnd[6])(index);
                  std::get<2>(pmlEdgeBndInterface[7])(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = std::get<2>(pmlEdgeBnd[7])(index);
                }
                sigmaInterface[5](index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = sigma[5](index);
                
                edgeInterface[8].first(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = edge[8](index);
                edgeInterface[9].first(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = edge[9](index);
                edgeInterface[10].first(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = edge[10](index);
                edgeInterface[11].first(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = edge[11](index);
                
                std::get<2>(cornerInterface[4])(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = corner[4](index);
                std::get<2>(cornerInterface[5])(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = corner[5](index);
                std::get<2>(cornerInterface[6])(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = corner[6](index);
                std::get<2>(cornerInterface[7])(index, i * nDomY * nDomZ + j * nDomZ + (k+1)) = corner[7](index);
    
                topology[index].push_back(i * nDomY * nDomZ + j * nDomZ + (k+1));
              }
            }
          }
        }
        
        // Kappa definition
        gmshfem::function::ScalarFunction< std::complex< double > > kappa;
        if(benchmark == "scattering") {
          kappa = k;
        }
        
        // Fields definition
        gmshddm::field::SubdomainField< std::complex< double >, gmshfem::field::Form::Form0 > u("u", omega | sigma[0] | sigma[1] | sigma[2] | sigma[3], gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
        if(benchmark == "scattering") {
          u(0).domain(u(0).domain() | gamma);
        }
        gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > lambda("lambda", gamma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
    
        std::vector< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > gContinuous;
    //    std::vector< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 > * > gDiscontinuous;
        for(unsigned int b = 0; b < 6; ++b) {
          gContinuous.push_back(new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gC_" + dir[b], sigmaInterface[b], gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder));
    //      gDiscontinuous.push_back(new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gD_" + dir[b], sigmaInterface[b], gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder));
        }
        
        std::vector< std::pair< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > > gEdgePmlContinuous;
    //    std::vector< std::pair< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 > * > > gEdgePmlDiscontinuous;
        std::vector< std::vector< std::pair< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > > > gEdgeHABC;
        for(unsigned int e = 0; e < 12; ++e) {
          if(e < 4) {
            gEdgePmlContinuous.push_back(std::make_pair(
              new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgePmlContinuous_" + dir[e%4] + "D_first", pmlBndInterface[e].first, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
              new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgePmlContinuous_" + dir[e%4] + "D_second", pmlBndInterface[e].second, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder)));
    //        gEdgePmlDiscontinuous.push_back(std::make_pair(
    //          new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gEdgePmlDiscontinuous_" + dir[e%4] + "D_first", pmlBndInterface[e].first, gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder),
    //          new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gEdgePmlDiscontinuous_" + dir[e%4] + "D_second", pmlBndInterface[e].second, gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder)));
            gEdgeHABC.push_back(std::vector< std::pair< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > >(N));
            for(unsigned int n = 0; n < N; ++n) {
              gEdgeHABC[e][n] = std::make_pair(
                new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgeHABC_" + dir[e%4] + "D_" + std::to_string(n) + "_first", edgeInterface[e].first, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
                new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgeHABC_" + dir[e%4] + "D_" + std::to_string(n) + "_second", edgeInterface[e].second, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder));
            }
          }
          else if(e >= 4 && e < 8) {
            gEdgePmlContinuous.push_back(std::make_pair(
              new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgePmlContinuous_" + dir[e%4] + dir[(e+1)%4] + "_first", pmlBndInterface[e].first, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
              new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgePmlContinuous_" + dir[e%4] + dir[(e+1)%4] + "_second", pmlBndInterface[e].second, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder)));
    //        gEdgePmlDiscontinuous.push_back(std::make_pair(
    //          new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gEdgePmlDiscontinuous_" + dir[e%4] + dir[(e+1)%4] + "_first", pmlBndInterface[e].first, gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder),
    //          new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gEdgePmlDiscontinuous_" + dir[e%4] + dir[(e+1)%4] + "_second", pmlBndInterface[e].second, gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder)));
            gEdgeHABC.push_back(std::vector< std::pair< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > >(N));
            for(unsigned int n = 0; n < N; ++n) {
              gEdgeHABC[e][n] = std::make_pair(
                new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgeHABC_" + dir[e%4] + dir[(e+1)%4] + "_" + std::to_string(n) + "_first", edgeInterface[e].first, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
                new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgeHABC_" + dir[e%4] + dir[(e+1)%4] + "_" + std::to_string(n) + "_second", edgeInterface[e].second, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder));
            }
          }
          else {
            gEdgePmlContinuous.push_back(std::make_pair(
              new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgePmlContinuous_" + dir[e%4] + "U_first", pmlBndInterface[e].first, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
              new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgePmlContinuous_" + dir[e%4] + "U_second", pmlBndInterface[e].second, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder)));
    //        gEdgePmlDiscontinuous.push_back(std::make_pair(
    //          new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gEdgePmlDiscontinuous_" + dir[e%4] + "U_first", pmlBndInterface[e].first, gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder),
    //          new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gEdgePmlDiscontinuous_" + dir[e%4] + "U_second", pmlBndInterface[e].second, gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder)));
            gEdgeHABC.push_back(std::vector< std::pair< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > >(N));
            for(unsigned int n = 0; n < N; ++n) {
              gEdgeHABC[e][n] = std::make_pair(
                new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgeHABC_" + dir[e%4] + "U_" + std::to_string(n) + "_first", edgeInterface[e].first, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
                new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gEdgeHABC_" + dir[e%4] + "U_" + std::to_string(n) + "_second", edgeInterface[e].second, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder));
            }
          }
        }
    
        std::vector< std::tuple< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > > gCornerPmlContinuous;
    //    std::vector< std::pair< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 > * > > gCornerPmlDiscontinuous;
        std::vector< std::vector< std::vector< std::tuple< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > > > > gCornerHABC;
        for(unsigned int c = 0; c < 8; ++c) {
          if(c < 4) {
            gCornerPmlContinuous.push_back(std::make_tuple(
              new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerPmlContinuous_" + dir[c%4] + dir[(c+1)%4] + "D_first", std::get<0>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
              new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerPmlContinuous_" + dir[c%4] + dir[(c+1)%4] + "D_second", std::get<1>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
              new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerPmlContinuous_" + dir[c%4] + dir[(c+1)%4] + "D_third", std::get<2>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder)));
    //        gCornerPmlDiscontinuous.push_back(std::make_tuple(
    //          new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gCornerPmlDiscontinuous_" + dir[c%4] + dir[(c+1)%4] + "D_first", std::get<0>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder),
    //          new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gCornerPmlDiscontinuous_" + dir[c%4] + dir[(c+1)%4] + "D_second", std::get<1>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder),
    //          new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gCornerPmlDiscontinuous_" + dir[c%4] + dir[(c+1)%4] + "D_third", std::get<2>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder)));
            gCornerHABC.push_back(std::vector< std::vector< std::tuple< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > > >(N));
            for(unsigned int n = 0; n < N; ++n) {
              gCornerHABC[c][n] = std::vector< std::tuple< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > >(N);
              for(unsigned int m = 0; m < N; ++m) {
                gCornerHABC[c][n][m] = std::make_tuple(
                  new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "D_" + std::to_string(n) + "_" + std::to_string(m) + "_first", std::get<0>(cornerInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
                  new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "D_" + std::to_string(n) + "_" + std::to_string(m) + "_second", std::get<1>(cornerInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
                  new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "D_" + std::to_string(n) + "_" + std::to_string(m) + "_third", std::get<2>(cornerInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder));
    //            gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<0>(cornerInterface[c])(0,1), "gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "D_" + std::to_string(n) + "_" + std::to_string(m) + "_first");
    //            gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<1>(cornerInterface[c])(0,1), "gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "D_" + std::to_string(n) + "_" + std::to_string(m) + "_second");
    //            gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<2>(cornerInterface[c])(0,1), "gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "D_" + std::to_string(n) + "_" + std::to_string(m) + "_third");
              }
            }
          }
          else {
            gCornerPmlContinuous.push_back(std::make_tuple(
              new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerPmlContinuous_" + dir[c%4] + dir[(c+1)%4] + "U_first", std::get<0>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
              new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerPmlContinuous_" + dir[c%4] + dir[(c+1)%4] + "U_second", std::get<1>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
              new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerPmlContinuous_" + dir[c%4] + dir[(c+1)%4] + "U_third", std::get<2>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder)));
    //        gCornerPmlDiscontinuous.push_back(std::make_tuple(
    //          new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gCornerPmlDiscontinuous_" + dir[c%4] + dir[(c+1)%4] + "U_first", std::get<0>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder),
    //          new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gCornerPmlDiscontinuous_" + dir[c%4] + dir[(c+1)%4] + "U_second", std::get<1>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder),
    //          new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form2 >("gCornerPmlDiscontinuous_" + dir[c%4] + dir[(c+1)%4] + "U_third", std::get<2>(pmlEdgeBndInterface[c]), gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl, neumannOrder)));
            gCornerHABC.push_back(std::vector< std::vector< std::tuple< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > > >(N));
            for(unsigned int n = 0; n < N; ++n) {
              gCornerHABC[c][n] = std::vector< std::tuple< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > *, gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > * > >(N);
              for(unsigned int m = 0; m < N; ++m) {
                gCornerHABC[c][n][m] = std::make_tuple(
                  new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "U_" + std::to_string(n) + "_" + std::to_string(m) + "_first", std::get<0>(cornerInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
                  new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "U_" + std::to_string(n) + "_" + std::to_string(m) + "_second", std::get<1>(cornerInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder),
                  new gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "U_" + std::to_string(n) + "_" + std::to_string(m) + "_third", std::get<2>(cornerInterface[c]), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder));
    //            gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<0>(cornerInterface[c])(0,1), "gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "U_" + std::to_string(n) + "_" + std::to_string(m) + "_first");
    //            gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<1>(cornerInterface[c])(0,1), "gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "U_" + std::to_string(n) + "_" + std::to_string(m) + "_second");
    //            gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<2>(cornerInterface[c])(0,1), "gCornerHABC_" + dir[c%4] + dir[(c+1)%4] + "U_" + std::to_string(n) + "_" + std::to_string(m) + "_third");
              }
            }
          }
        }
    
        gmshddm::problem::Formulation< std::complex< double > > formulation("HelmholtzDDMScattering", topology);
    
        for(unsigned int b = 0; b < 6; ++b) {
          formulation.addInterfaceField(*gContinuous[b]);
    //      formulation.addInterfaceField(*gDiscontinuous[b]);
        }
        
        for(unsigned int e = 0; e < 12; ++e) {
          formulation.addInterfaceField(*gEdgePmlContinuous[e].first);
          formulation.addInterfaceField(*gEdgePmlContinuous[e].second);
          
    //      formulation.addInterfaceField(*gEdgePmlDiscontinuous[e].first);
    //      formulation.addInterfaceField(*gEdgePmlDiscontinuous[e].second);
          for(unsigned int n = 0; n < N; ++n) {
            formulation.addInterfaceField(*gEdgeHABC[e][n].first);
            formulation.addInterfaceField(*gEdgeHABC[e][n].second);
          }
        }
        
        for(unsigned int c = 0; c < 8; ++c) {
          formulation.addInterfaceField(*std::get<0>(gCornerPmlContinuous[c]));
          formulation.addInterfaceField(*std::get<1>(gCornerPmlContinuous[c]));
          formulation.addInterfaceField(*std::get<2>(gCornerPmlContinuous[c]));
          
    //      formulation.addInterfaceField(*std::get<0>(gCornerPmlDiscontinuous[c]));
    //      formulation.addInterfaceField(*std::get<1>(gCornerPmlDiscontinuous[c]));
    //      formulation.addInterfaceField(*std::get<2>(gCornerPmlDiscontinuous[c]));
          for(unsigned int n = 0; n < N; ++n) {
            for(unsigned int m = 0; m < N; ++m) {
              formulation.addInterfaceField(*std::get<0>(gCornerHABC[c][n][m]));
              formulation.addInterfaceField(*std::get<1>(gCornerHABC[c][n][m]));
              formulation.addInterfaceField(*std::get<2>(gCornerHABC[c][n][m]));
            }
          }
        }
        
        std::vector< std::vector< std::vector< Subproblem * > > > subproblem(nDomX);
        for(unsigned int i = 0; i < nDomX; ++i) {
          subproblem[i].resize(nDomY);
          for(unsigned int j = 0; j < nDomY; ++j) {
            subproblem[i][j].resize(nDomZ);
            for(unsigned int k = 0; k < nDomZ; ++k) {
              unsigned int index = i * nDomY * nDomZ + j * nDomZ + k;
              formulation(index).integral(-grad(dof(u(index))), grad(tf(u(index))), omega(index), gauss);
              formulation(index).integral(kappa * kappa * dof(u(index)), tf(u(index)), omega(index), gauss);
    
              if(benchmark == "scattering") {
                if(index == 0) {
                  formulation(index).integral(dof(lambda), tf(u(index)), gamma, gauss);
                  formulation(index).integral(dof(u(index)), tf(lambda), gamma, gauss);
                  formulation(index).integral(formulation.physicalSource(- fAnalytic), tf(lambda), gamma, gauss);
                }
              }
    
              SubproblemDomains domains;
              domains.setOmega(omega(index));
              domains.setSigma({ sigma[0](index), sigma[1](index), sigma[2](index), sigma[3](index), sigma[4](index), sigma[5](index) });
              domains.setPml({ pml[0](index), pml[1](index), pml[2](index), pml[3](index), pml[4](index), pml[5](index) });
              domains.setPmlBnd({ std::make_pair(pmlBnd[0].first(index), pmlBnd[0].second(index)), std::make_pair(pmlBnd[1].first(index), pmlBnd[1].second(index)), std::make_pair(pmlBnd[2].first(index), pmlBnd[2].second(index)), std::make_pair(pmlBnd[3].first(index), pmlBnd[3].second(index)), std::make_pair(pmlBnd[4].first(index), pmlBnd[4].second(index)), std::make_pair(pmlBnd[5].first(index), pmlBnd[5].second(index)), std::make_pair(pmlBnd[6].first(index), pmlBnd[6].second(index)), std::make_pair(pmlBnd[7].first(index), pmlBnd[7].second(index)), std::make_pair(pmlBnd[8].first(index), pmlBnd[8].second(index)), std::make_pair(pmlBnd[9].first(index), pmlBnd[9].second(index)), std::make_pair(pmlBnd[10].first(index), pmlBnd[10].second(index)), std::make_pair(pmlBnd[11].first(index), pmlBnd[11].second(index)) });
              domains.setPmlEdge({ pmlEdge[0](index), pmlEdge[1](index), pmlEdge[2](index), pmlEdge[3](index), pmlEdge[4](index), pmlEdge[5](index), pmlEdge[6](index), pmlEdge[7](index), pmlEdge[8](index), pmlEdge[9](index), pmlEdge[10](index), pmlEdge[11](index) });
              domains.setEdge({ edge[0](index), edge[1](index), edge[2](index), edge[3](index), edge[4](index), edge[5](index), edge[6](index), edge[7](index), edge[8](index), edge[9](index), edge[10](index), edge[11](index) });
              domains.setPmlEdgeBnd({ std::make_tuple(std::get<0>(pmlEdgeBnd[0])(index), std::get<1>(pmlEdgeBnd[0])(index), std::get<2>(pmlEdgeBnd[0])(index)), std::make_tuple(std::get<0>(pmlEdgeBnd[1])(index), std::get<1>(pmlEdgeBnd[1])(index), std::get<2>(pmlEdgeBnd[1])(index)), std::make_tuple(std::get<0>(pmlEdgeBnd[2])(index), std::get<1>(pmlEdgeBnd[2])(index), std::get<2>(pmlEdgeBnd[2])(index)), std::make_tuple(std::get<0>(pmlEdgeBnd[3])(index), std::get<1>(pmlEdgeBnd[3])(index), std::get<2>(pmlEdgeBnd[3])(index)), std::make_tuple(std::get<0>(pmlEdgeBnd[4])(index), std::get<1>(pmlEdgeBnd[4])(index), std::get<2>(pmlEdgeBnd[4])(index)), std::make_tuple(std::get<0>(pmlEdgeBnd[5])(index), std::get<1>(pmlEdgeBnd[5])(index), std::get<2>(pmlEdgeBnd[5])(index)), std::make_tuple(std::get<0>(pmlEdgeBnd[6])(index), std::get<1>(pmlEdgeBnd[6])(index), std::get<2>(pmlEdgeBnd[6])(index)), std::make_tuple(std::get<0>(pmlEdgeBnd[7])(index), std::get<1>(pmlEdgeBnd[7])(index), std::get<2>(pmlEdgeBnd[7])(index)) });
              domains.setPmlCorner({ pmlCorner[0](index), pmlCorner[1](index), pmlCorner[2](index), pmlCorner[3](index), pmlCorner[4](index), pmlCorner[5](index), pmlCorner[6](index), pmlCorner[7](index) });
              domains.setCorner({ corner[0](index), corner[1](index), corner[2](index), corner[3](index), corner[4](index), corner[5](index), corner[6](index), corner[7](index) });
              domains.setCornerEdge({ std::make_tuple(std::get<0>(cornerEdge[0])(index), std::get<1>(cornerEdge[0])(index), std::get<2>(cornerEdge[0])(index)), std::make_tuple(std::get<0>(cornerEdge[1])(index), std::get<1>(cornerEdge[1])(index), std::get<2>(cornerEdge[1])(index)), std::make_tuple(std::get<0>(cornerEdge[2])(index), std::get<1>(cornerEdge[2])(index), std::get<2>(cornerEdge[2])(index)), std::make_tuple(std::get<0>(cornerEdge[3])(index), std::get<1>(cornerEdge[3])(index), std::get<2>(cornerEdge[3])(index)), std::make_tuple(std::get<0>(cornerEdge[4])(index), std::get<1>(cornerEdge[4])(index), std::get<2>(cornerEdge[4])(index)), std::make_tuple(std::get<0>(cornerEdge[5])(index), std::get<1>(cornerEdge[5])(index), std::get<2>(cornerEdge[5])(index)), std::make_tuple(std::get<0>(cornerEdge[6])(index), std::get<1>(cornerEdge[6])(index), std::get<2>(cornerEdge[6])(index)), std::make_tuple(std::get<0>(cornerEdge[7])(index), std::get<1>(cornerEdge[7])(index), std::get<2>(cornerEdge[7])(index)) });
    
              SubproblemParameters parameters;
              parameters.setGauss(gauss);
              parameters.setKappa(kappa);
              parameters.setNeumannOrder(neumannOrder);
              parameters.setFieldOrder(fieldOrder);
              parameters.setStab(stab);
    
              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;
                }
              }
              std::string bndExt = boundaryExt;
              if(boundaryExt == "habc") {
                bndExt += "_" + std::to_string(NExt) + "_" + std::to_string(thetaPade);
              }
              else if(boundaryExt == "pml") {
                if(pmlMethodExt == "continuous") {
                  bndExt += "Continuous_" + std::to_string(pmlSizeExt) + "_" + pmlTypeExt;
                }
                else if(pmlMethodExt == "discontinuous") {
                  bndExt += "Discontinuous_" + std::to_string(pmlSizeExt) + "_" + pmlTypeExt;
                }
              }
              gmshfem::msg::info << "Subdomain (" << i << ", " << j << ", " << k << ") has boundaries [" << (i == nDomX-1 ? bndExt : bnd) << ", " << (j == nDomY-1 ? bndExt : bnd) << ", " << (i == 0 ? bndExt : bnd) << ", " << (j == 0 ? bndExt : bnd) << ", " << (k == 0 ? bndExt : bnd) << ", " << (k == nDomZ-1 ? bndExt : bnd) << "]" << gmshfem::msg::endl;
              subproblem[i][j][k] = new Subproblem(formulation(index), u(index).name(), domains, parameters, (i == nDomX-1 ? bndExt : bnd), (j == nDomY-1 ? bndExt : bnd), (i == 0 ? bndExt : bnd), (j == 0 ? bndExt : bnd), (k == 0 ? bndExt : bnd), (k == nDomZ-1 ? bndExt : bnd));
              subproblem[i][j][k]->writeFormulation();
    
              // coupling
              const int faceNeighbor[6] = {
                2, 3, 0, 1, 5, 4
              };
              for(unsigned int b = 0; b < 6; ++b) {
                for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
                  const unsigned int jndex = topology[index][jj];
                  
                  Boundary *boundary = subproblem[i][j][k]->getBoundary(b);
                  if(dynamic_cast< Sommerfeld * >(boundary)) {
                    if(!sigmaInterface[b](index, jndex).isEmpty()) {
                      formulation(index).integral(formulation.artificialSource( - (*gContinuous[faceNeighbor[b]])(jndex, index)), tf(u(index)), sigmaInterface[b](index, jndex), gauss);
                    }
                  }
                  else if(dynamic_cast< HABC * >(boundary)) {
                    if(!sigmaInterface[b](index, jndex).isEmpty()) {
                      formulation(index).integral(formulation.artificialSource( - (*gContinuous[faceNeighbor[b]])(jndex, index)), tf(u(index)), sigmaInterface[b](index, jndex), gauss);
                    }
                  }
                  else if(dynamic_cast< PmlContinuous * >(boundary)) {
                    if(!sigmaInterface[b](index, jndex).isEmpty()) {
                      formulation(index).integral(formulation.artificialSource( - (*gContinuous[faceNeighbor[b]])(jndex, index)), tf(u(index)), sigmaInterface[b](index, jndex), gauss);
                    }
                  }
                }
              }
              
              const int edgeNeighbor[12][2] = {
                // Down
                {8, 2},
                {9, 3},
                {10, 0},
                {11, 1},
                // Rot
                {7, 5},
                {4, 6},
                {5, 7},
                {6, 4},
                // Up
                {0, 10},
                {1, 11},
                {2, 8},
                {3, 9}
              };
              for(unsigned int e = 0; e < 12; ++e) {
                for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
                  const unsigned int jndex = topology[index][jj];
    
                  Edge *edge = subproblem[i][j][k]->getEdge(e);
                  if(edge == nullptr) {
                    continue;
                  }
                  if(auto eg = dynamic_cast< HABC_HABC * >(edge)) {
                    if(!edgeInterface[e].first(index, jndex).isEmpty()) {
                      if(e >= 4 && e < 8) {
                        for(unsigned int n = 0; n < N; ++n) {
                          formulation(index).integral(formulation.artificialSource( - (*gEdgeHABC[edgeNeighbor[e][0]][n].second)(jndex, index) ), tf(*eg->firstBoundary()->getUHABC(n)), edgeInterface[e].first(index, jndex), gauss);
                        }
                      }
                      else {
                        for(unsigned int n = 0; n < N; ++n) {
                          formulation(index).integral(formulation.artificialSource( - (*gEdgeHABC[edgeNeighbor[e][0]][n].first)(jndex, index) ), tf(*eg->firstBoundary()->getUHABC(n)), edgeInterface[e].first(index, jndex), gauss);
                        }
                      }
                    }
                    if(!edgeInterface[e].second(index, jndex).isEmpty()) {
                      if(e >= 4 && e < 8) {
                        for(unsigned int n = 0; n < N; ++n) {
                          formulation(index).integral(formulation.artificialSource( - (*gEdgeHABC[edgeNeighbor[e][1]][n].first)(jndex, index) ), tf(*eg->secondBoundary()->getUHABC(n)), edgeInterface[e].second(index, jndex), gauss);
                        }
                      }
                      else {
                        for(unsigned int n = 0; n < N; ++n) {
                          formulation(index).integral(formulation.artificialSource( - (*gEdgeHABC[edgeNeighbor[e][1]][n].second)(jndex, index) ), tf(*eg->secondBoundary()->getUHABC(n)), edgeInterface[e].second(index, jndex), gauss);
                        }
                      }
                    }
                  }
                  else if(auto eg = dynamic_cast< PmlContinuous_PmlContinuous * >(edge)) {
                    if(!pmlBndInterface[e].first(index, jndex).isEmpty()) {
                      if(e >= 4 && e < 8) {
                        formulation(index).integral(formulation.artificialSource( - (*gEdgePmlContinuous[edgeNeighbor[e][0]].second)(jndex, index) ), tf(*eg->firstBoundary()->getUPml()), pmlBndInterface[e].first(index, jndex), gauss);
                      }
                      else {
                        formulation(index).integral(formulation.artificialSource( - (*gEdgePmlContinuous[edgeNeighbor[e][0]].first)(jndex, index) ), tf(*eg->firstBoundary()->getUPml()), pmlBndInterface[e].first(index, jndex), gauss);
                      }
                    }
                    if(!pmlBndInterface[e].second(index, jndex).isEmpty()) {
                      if(e >= 4 && e < 8) {
                        formulation(index).integral(formulation.artificialSource( - (*gEdgePmlContinuous[edgeNeighbor[e][1]].first)(jndex, index) ), tf(*eg->secondBoundary()->getUPml()), pmlBndInterface[e].second(index, jndex), gauss);
                      }
                      else {
                        formulation(index).integral(formulation.artificialSource( - (*gEdgePmlContinuous[edgeNeighbor[e][1]].second)(jndex, index) ), tf(*eg->secondBoundary()->getUPml()), pmlBndInterface[e].second(index, jndex), gauss);
                      }
                    }
                  }
                }
              }
    
              const int cornerNeighbor[8][3] = {
                // Down
                {3, 1, 4},
                {0, 2, 5},
                {1, 3, 6},
                {2, 0, 7},
                // Up
                {7, 5, 0},
                {4, 6, 1},
                {5, 7, 2},
                {6, 4, 3}
              };
              for(unsigned int c = 0; c < 8; ++c) {
                for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
                  const unsigned int jndex = topology[index][jj];
    
                  Corner *corner = subproblem[i][j][k]->getCorner(c);
                  if(corner == nullptr) {
                    continue;
                  }
                  if(auto cr = dynamic_cast< HABC_HABC_HABC * >(corner)) {
                    if(!std::get<0>(cornerInterface[c])(index, jndex).isEmpty()) {
                      for(unsigned int n = 0; n < N; ++n) {
                        for(unsigned int m = 0; m < N; ++m) {
                          formulation(index).integral(formulation.artificialSource( - (*std::get<1>(gCornerHABC[cornerNeighbor[c][0]][n][m]))(jndex, index) ), tf(*cr->firstEdge()->getUHABC(n, m)), std::get<0>(cornerInterface[c])(index, jndex), gauss);
                        }
                      }
                    }
                    if(!std::get<1>(cornerInterface[c])(index, jndex).isEmpty()) {
                      for(unsigned int n = 0; n < N; ++n) {
                        for(unsigned int m = 0; m < N; ++m) {
                          formulation(index).integral(formulation.artificialSource( - (*std::get<0>(gCornerHABC[cornerNeighbor[c][1]][n][m]))(jndex, index) ), tf(*cr->secondEdge()->getUHABC(n, m)), std::get<1>(cornerInterface[c])(index, jndex), gauss);
                        }
                      }
                    }
                    if(!std::get<2>(cornerInterface[c])(index, jndex).isEmpty()) {
                      for(unsigned int n = 0; n < N; ++n) {
                        for(unsigned int m = 0; m < N; ++m) {
                          formulation(index).integral(formulation.artificialSource( - (*std::get<2>(gCornerHABC[cornerNeighbor[c][2]][n][m]))(jndex, index) ), tf(*cr->thirdEdge()->getUHABC(n, m)), std::get<2>(cornerInterface[c])(index, jndex), gauss);
                        }
                      }
                    }
                  }
                  else if(auto cr = dynamic_cast< PmlContinuous_PmlContinuous_PmlContinuous * >(corner)) {
                    if(!std::get<0>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
                      formulation(index).integral(formulation.artificialSource( - (*std::get<1>(gCornerPmlContinuous[cornerNeighbor[c][0]]))(jndex, index) ), tf(*cr->firstEdge()->getUPml()), std::get<0>(pmlEdgeBndInterface[c])(index, jndex), gauss);
                    }
                    if(!std::get<1>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
                      formulation(index).integral(formulation.artificialSource( - (*std::get<0>(gCornerPmlContinuous[cornerNeighbor[c][1]]))(jndex, index) ), tf(*cr->secondEdge()->getUPml()), std::get<1>(pmlEdgeBndInterface[c])(index, jndex), gauss);
                    }
                    if(!std::get<2>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
                      formulation(index).integral(formulation.artificialSource( - (*std::get<2>(gCornerPmlContinuous[cornerNeighbor[c][2]]))(jndex, index) ), tf(*cr->thirdEdge()->getUPml()), std::get<2>(pmlEdgeBndInterface[c])(index, jndex), gauss);
                    }
                  }
                }
              }
    
              // interface
              for(unsigned int b = 0; b < 6; ++b) {
                for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
                  const unsigned int jndex = topology[index][jj];
    
                  Boundary *boundary = subproblem[i][j][k]->getBoundary(b);
                  if(auto bnd = dynamic_cast< Sommerfeld * >(boundary)) {
                    if(!sigmaInterface[b](index, jndex).isEmpty()) {
                      formulation(index, jndex).integral(dof((*gContinuous[b])(index, jndex)), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
                      formulation(index, jndex).integral(formulation.artificialSource( (*gContinuous[faceNeighbor[b]])(jndex, index) ), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
                      formulation(index, jndex).integral(2. * *bnd->getV(), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
                    }
                  }
                  else if(auto bnd = dynamic_cast< HABC * >(boundary)) {
                    if(!sigmaInterface[b](index, jndex).isEmpty()) {
                      formulation(index, jndex).integral(dof((*gContinuous[b])(index, jndex)), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
                      formulation(index, jndex).integral(formulation.artificialSource( (*gContinuous[faceNeighbor[b]])(jndex, index) ), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
                      formulation(index, jndex).integral(2. * *bnd->getV(), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
                    }
                  }
                  else if(auto bnd = dynamic_cast< PmlContinuous * >(boundary)) {
                    if(!sigmaInterface[b](index, jndex).isEmpty()) {
                      formulation(index, jndex).integral(dof((*gContinuous[b])(index, jndex)), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
                      formulation(index, jndex).integral(formulation.artificialSource( (*gContinuous[faceNeighbor[b]])(jndex, index) ), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
                      formulation(index, jndex).integral(2. * *bnd->getV(), tf((*gContinuous[b])(index, jndex)), sigmaInterface[b](index, jndex), gauss);
                    }
                  }
                }
              }
              
              for(unsigned int e = 0; e < 12; ++e) {
                for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
                  const unsigned int jndex = topology[index][jj];
    
                  Edge *edge = subproblem[i][j][k]->getEdge(e);
                  if(edge == nullptr) {
                    continue;
                  }
                  if(auto eg = dynamic_cast< HABC_HABC * >(edge)) {
                    if(!edgeInterface[e].first(index, jndex).isEmpty()) {
                      if(e >= 4 && e < 8) {
                        for(unsigned int n = 0; n < N; ++n) {
                          formulation(index, jndex).integral(dof((*gEdgeHABC[e][n].first)(index, jndex)), tf((*gEdgeHABC[e][n].first)(index, jndex)), edgeInterface[e].first(index, jndex), gauss);
                          formulation(index, jndex).integral(formulation.artificialSource( (*gEdgeHABC[edgeNeighbor[e][0]][n].second)(jndex, index) ), tf((*gEdgeHABC[e][n].first)(index, jndex)), edgeInterface[e].first(index, jndex), gauss);
                          formulation(index, jndex).integral(2. * *eg->getV(n, 0), tf((*gEdgeHABC[e][n].first)(index, jndex)), edgeInterface[e].first(index, jndex), gauss);
                        }
                      }
                      else {
                        for(unsigned int n = 0; n < N; ++n) {
                          formulation(index, jndex).integral(dof((*gEdgeHABC[e][n].first)(index, jndex)), tf((*gEdgeHABC[e][n].first)(index, jndex)), edgeInterface[e].first(index, jndex), gauss);
                          formulation(index, jndex).integral(formulation.artificialSource( (*gEdgeHABC[edgeNeighbor[e][0]][n].first)(jndex, index) ), tf((*gEdgeHABC[e][n].first)(index, jndex)), edgeInterface[e].first(index, jndex), gauss);
                          formulation(index, jndex).integral(2. * *eg->getV(n, 0), tf((*gEdgeHABC[e][n].first)(index, jndex)), edgeInterface[e].first(index, jndex), gauss);
                        }
                      }
                    }
                    if(!edgeInterface[e].second(index, jndex).isEmpty()) {
                      if(e >= 4 && e < 8) {
                        for(unsigned int n = 0; n < N; ++n) {
                          formulation(index, jndex).integral(dof((*gEdgeHABC[e][n].second)(index, jndex)), tf((*gEdgeHABC[e][n].second)(index, jndex)), edgeInterface[e].second(index, jndex), gauss);
                          formulation(index, jndex).integral(formulation.artificialSource( (*gEdgeHABC[edgeNeighbor[e][1]][n].first)(jndex, index) ), tf((*gEdgeHABC[e][n].second)(index, jndex)), edgeInterface[e].second(index, jndex), gauss);
                          formulation(index, jndex).integral(2. * *eg->getV(n, 1), tf((*gEdgeHABC[e][n].second)(index, jndex)), edgeInterface[e].second(index, jndex), gauss);
                        }
                      }
                      else {
                        for(unsigned int n = 0; n < N; ++n) {
                          formulation(index, jndex).integral(dof((*gEdgeHABC[e][n].second)(index, jndex)), tf((*gEdgeHABC[e][n].second)(index, jndex)), edgeInterface[e].second(index, jndex), gauss);
                          formulation(index, jndex).integral(formulation.artificialSource( (*gEdgeHABC[edgeNeighbor[e][1]][n].second)(jndex, index) ), tf((*gEdgeHABC[e][n].second)(index, jndex)), edgeInterface[e].second(index, jndex), gauss);
                          formulation(index, jndex).integral(2. * *eg->getV(n, 1), tf((*gEdgeHABC[e][n].second)(index, jndex)), edgeInterface[e].second(index, jndex), gauss);
                        }
                      }
                    }
                  }
                  else if(auto eg = dynamic_cast< PmlContinuous_PmlContinuous * >(edge)) {
                    if(!pmlBndInterface[e].first(index, jndex).isEmpty()) {
                      if(e >= 4 && e < 8) {
                        formulation(index, jndex).integral(dof((*gEdgePmlContinuous[e].first)(index, jndex)), tf((*gEdgePmlContinuous[e].first)(index, jndex)), pmlBndInterface[e].first(index, jndex), gauss);
                        formulation(index, jndex).integral(formulation.artificialSource( (*gEdgePmlContinuous[edgeNeighbor[e][0]].second)(jndex, index) ), tf((*gEdgePmlContinuous[e].first)(index, jndex)), pmlBndInterface[e].first(index, jndex), gauss);
                        formulation(index, jndex).integral(2. * *eg->getV(0), tf((*gEdgePmlContinuous[e].first)(index, jndex)), pmlBndInterface[e].first(index, jndex), gauss);
                      }
                      else {
                        formulation(index, jndex).integral(dof((*gEdgePmlContinuous[e].first)(index, jndex)), tf((*gEdgePmlContinuous[e].first)(index, jndex)), pmlBndInterface[e].first(index, jndex), gauss);
                        formulation(index, jndex).integral(formulation.artificialSource( (*gEdgePmlContinuous[edgeNeighbor[e][0]].first)(jndex, index) ), tf((*gEdgePmlContinuous[e].first)(index, jndex)), pmlBndInterface[e].first(index, jndex), gauss);
                        formulation(index, jndex).integral(2. * *eg->getV(0), tf((*gEdgePmlContinuous[e].first)(index, jndex)), pmlBndInterface[e].first(index, jndex), gauss);
                      }
                    }
                    if(!pmlBndInterface[e].second(index, jndex).isEmpty()) {
                      if(e >= 4 && e < 8) {
                        formulation(index, jndex).integral(dof((*gEdgePmlContinuous[e].second)(index, jndex)), tf((*gEdgePmlContinuous[e].second)(index, jndex)), pmlBndInterface[e].second(index, jndex), gauss);
                        formulation(index, jndex).integral(formulation.artificialSource( (*gEdgePmlContinuous[edgeNeighbor[e][1]].first)(jndex, index) ), tf((*gEdgePmlContinuous[e].second)(index, jndex)), pmlBndInterface[e].second(index, jndex), gauss);
                        formulation(index, jndex).integral(2. * *eg->getV(1), tf((*gEdgePmlContinuous[e].second)(index, jndex)), pmlBndInterface[e].second(index, jndex), gauss);
                      }
                      else {
                        formulation(index, jndex).integral(dof((*gEdgePmlContinuous[e].second)(index, jndex)), tf((*gEdgePmlContinuous[e].second)(index, jndex)), pmlBndInterface[e].second(index, jndex), gauss);
                        formulation(index, jndex).integral(formulation.artificialSource( (*gEdgePmlContinuous[edgeNeighbor[e][1]].second)(jndex, index) ), tf((*gEdgePmlContinuous[e].second)(index, jndex)), pmlBndInterface[e].second(index, jndex), gauss);
                        formulation(index, jndex).integral(2. * *eg->getV(1), tf((*gEdgePmlContinuous[e].second)(index, jndex)), pmlBndInterface[e].second(index, jndex), gauss);
                      }
                    }
                  }
                }
              }
              
              for(unsigned int c = 0; c < 8; ++c) {
                for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
                  const unsigned int jndex = topology[index][jj];
    
                  Corner *corner = subproblem[i][j][k]->getCorner(c);
                  if(corner == nullptr) {
                    continue;
                  }
                  if(auto cr = dynamic_cast< HABC_HABC_HABC * >(corner)) {
                    if(!std::get<0>(cornerInterface[c])(index, jndex).isEmpty()) {
                      for(unsigned int n = 0; n < N; ++n) {
                        for(unsigned int m = 0; m < N; ++m) {
                          formulation(index, jndex).integral(dof((*std::get<0>(gCornerHABC[c][n][m]))(index, jndex)), tf((*std::get<0>(gCornerHABC[c][n][m]))(index, jndex)), std::get<0>(cornerInterface[c])(index, jndex), gauss);
                          formulation(index, jndex).integral(formulation.artificialSource( (*std::get<1>(gCornerHABC[cornerNeighbor[c][0]][n][m]))(jndex, index) ), tf((*std::get<0>(gCornerHABC[c][n][m]))(index, jndex)), std::get<0>(cornerInterface[c])(index, jndex), gauss);
                          formulation(index, jndex).integral(2. * *cr->getV(n, m, 0), tf((*std::get<0>(gCornerHABC[c][n][m]))(index, jndex)), std::get<0>(cornerInterface[c])(index, jndex), gauss);
                        }
                      }
                    }
                    if(!std::get<1>(cornerInterface[c])(index, jndex).isEmpty()) {
                      for(unsigned int n = 0; n < N; ++n) {
                        for(unsigned int m = 0; m < N; ++m) {
                          formulation(index, jndex).integral(dof((*std::get<1>(gCornerHABC[c][n][m]))(index, jndex)), tf((*std::get<1>(gCornerHABC[c][n][m]))(index, jndex)), std::get<1>(cornerInterface[c])(index, jndex), gauss);
                          formulation(index, jndex).integral(formulation.artificialSource( (*std::get<0>(gCornerHABC[cornerNeighbor[c][1]][n][m]))(jndex, index) ), tf((*std::get<1>(gCornerHABC[c][n][m]))(index, jndex)), std::get<1>(cornerInterface[c])(index, jndex), gauss);
                          formulation(index, jndex).integral(2. * *cr->getV(n, m, 1), tf((*std::get<1>(gCornerHABC[c][n][m]))(index, jndex)), std::get<1>(cornerInterface[c])(index, jndex), gauss);
                        }
                      }
                    }
                    if(!std::get<2>(cornerInterface[c])(index, jndex).isEmpty()) {
                      for(unsigned int n = 0; n < N; ++n) {
                        for(unsigned int m = 0; m < N; ++m) {
                          formulation(index, jndex).integral(dof((*std::get<2>(gCornerHABC[c][n][m]))(index, jndex)), tf((*std::get<2>(gCornerHABC[c][n][m]))(index, jndex)), std::get<2>(cornerInterface[c])(index, jndex), gauss);
                          formulation(index, jndex).integral(formulation.artificialSource( (*std::get<2>(gCornerHABC[cornerNeighbor[c][2]][n][m]))(jndex, index) ), tf((*std::get<2>(gCornerHABC[c][n][m]))(index, jndex)), std::get<2>(cornerInterface[c])(index, jndex), gauss);
                          formulation(index, jndex).integral(2. * *cr->getV(n, m, 2), tf((*std::get<2>(gCornerHABC[c][n][m]))(index, jndex)), std::get<2>(cornerInterface[c])(index, jndex), gauss);
                        }
                      }
                    }
                  }
                  else if(auto cr = dynamic_cast< PmlContinuous_PmlContinuous_PmlContinuous * >(corner)) {
                    if(!std::get<0>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
                      formulation(index, jndex).integral(dof((*std::get<0>(gCornerPmlContinuous[c]))(index, jndex)), tf((*std::get<0>(gCornerPmlContinuous[c]))(index, jndex)), std::get<0>(pmlEdgeBndInterface[c])(index, jndex), gauss);
                      formulation(index, jndex).integral(formulation.artificialSource( (*std::get<1>(gCornerPmlContinuous[cornerNeighbor[c][0]]))(jndex, index) ), tf((*std::get<0>(gCornerPmlContinuous[c]))(index, jndex)), std::get<0>(pmlEdgeBndInterface[c])(index, jndex), gauss);
                      formulation(index, jndex).integral(2. * *cr->getV(0), tf((*std::get<0>(gCornerPmlContinuous[c]))(index, jndex)), std::get<0>(pmlEdgeBndInterface[c])(index, jndex), gauss);
                    }
                    if(!std::get<1>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
                      formulation(index, jndex).integral(dof((*std::get<1>(gCornerPmlContinuous[c]))(index, jndex)), tf((*std::get<1>(gCornerPmlContinuous[c]))(index, jndex)), std::get<1>(pmlEdgeBndInterface[c])(index, jndex), gauss);
                      formulation(index, jndex).integral(formulation.artificialSource( (*std::get<0>(gCornerPmlContinuous[cornerNeighbor[c][1]]))(jndex, index) ), tf((*std::get<1>(gCornerPmlContinuous[c]))(index, jndex)), std::get<1>(pmlEdgeBndInterface[c])(index, jndex), gauss);
                      formulation(index, jndex).integral(2. * *cr->getV(1), tf((*std::get<1>(gCornerPmlContinuous[c]))(index, jndex)), std::get<1>(pmlEdgeBndInterface[c])(index, jndex), gauss);
                    }
                    if(!std::get<2>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
                      formulation(index, jndex).integral(dof((*std::get<2>(gCornerPmlContinuous[c]))(index, jndex)), tf((*std::get<2>(gCornerPmlContinuous[c]))(index, jndex)), std::get<2>(pmlEdgeBndInterface[c])(index, jndex), gauss);
                      formulation(index, jndex).integral(formulation.artificialSource( (*std::get<2>(gCornerPmlContinuous[cornerNeighbor[c][2]]))(jndex, index) ), tf((*std::get<2>(gCornerPmlContinuous[c]))(index, jndex)), std::get<2>(pmlEdgeBndInterface[c])(index, jndex), gauss);
                      formulation(index, jndex).integral(2. * *cr->getV(2), tf((*std::get<2>(gCornerPmlContinuous[c]))(index, jndex)), std::get<2>(pmlEdgeBndInterface[c])(index, jndex), gauss);
                    }
                  }
                }
              }
              
    //          for(unsigned int e = 0; e < 12; ++e) {
    //            for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
    //              const unsigned int jndex = topology[index][jj];
    //
    //              Edge *edge = subproblem[i][j][k]->getEdge(e);
    //              if(edge == nullptr) {
    //                continue;
    //              }
    //              if(auto eg = dynamic_cast< HABC_HABC * >(edge)) {
    //                if(!edgeInterface[e].first(index, jndex).isEmpty()) {
    //                  if(e >= 4 && e < 8) {
    //                    for(unsigned int n = 0; n < N; ++n) {
    //                      gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgeHABC[e][n].first)(index, jndex), (*gEdgeHABC[e][n].first)(index, jndex).domain(), (*gEdgeHABC[e][n].first)(index, jndex).name());
    //                      gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1) + (*gEdgeHABC[edgeNeighbor[e][0]][n].second)(jndex, index), (*gEdgeHABC[edgeNeighbor[e][0]][n].second)(jndex, index).domain(), (*gEdgeHABC[edgeNeighbor[e][0]][n].second)(jndex, index).name() + "_pair");
    //                    }
    //                  }
    //                  else {
    //                    for(unsigned int n = 0; n < N; ++n) {
    //                      gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgeHABC[e][n].first)(index, jndex), (*gEdgeHABC[e][n].first)(index, jndex).domain(), (*gEdgeHABC[e][n].first)(index, jndex).name());
    //                      gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgeHABC[edgeNeighbor[e][0]][n].first)(jndex, index), (*gEdgeHABC[edgeNeighbor[e][0]][n].first)(jndex, index).domain(), (*gEdgeHABC[edgeNeighbor[e][0]][n].first)(jndex, index).name() + "_pair");
    //                    }
    //                  }
    //                }
    //                if(!edgeInterface[e].second(index, jndex).isEmpty()) {
    //                  if(e >= 4 && e < 8) {
    //                    for(unsigned int n = 0; n < N; ++n) {
    //                      gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgeHABC[e][n].second)(index, jndex), (*gEdgeHABC[e][n].second)(index, jndex).domain(), (*gEdgeHABC[e][n].second)(index, jndex).name());
    //                      gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgeHABC[edgeNeighbor[e][1]][n].first)(jndex, index), (*gEdgeHABC[edgeNeighbor[e][1]][n].first)(jndex, index).domain(), (*gEdgeHABC[edgeNeighbor[e][1]][n].first)(jndex, index).name() + "_pair");
    //                    }
    //                  }
    //                  else {
    //                    for(unsigned int n = 0; n < N; ++n) {
    //                      gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgeHABC[e][n].second)(index, jndex), (*gEdgeHABC[e][n].second)(index, jndex).domain(), (*gEdgeHABC[e][n].second)(index, jndex).name());
    //                      gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgeHABC[edgeNeighbor[e][1]][n].second)(jndex, index), (*gEdgeHABC[edgeNeighbor[e][1]][n].second)(jndex, index).domain(), (*gEdgeHABC[edgeNeighbor[e][1]][n].second)(jndex, index).name() + "_pair");
    //                    }
    //                  }
    //                }
    //              }
    //              else if(auto eg = dynamic_cast< PmlContinuous_PmlContinuous * >(edge)) {
    //                if(!pmlBndInterface[e].first(index, jndex).isEmpty()) {
    //                  if(e >= 4 && e < 8) {
    //                    gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgePmlContinuous[e].first)(index, jndex), (*gEdgePmlContinuous[e].first)(index, jndex).domain(), (*gEdgePmlContinuous[e].first)(index, jndex).name());
    //                    gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1) + (*gEdgePmlContinuous[edgeNeighbor[e][0]].second)(jndex, index), (*gEdgePmlContinuous[edgeNeighbor[e][0]].second)(jndex, index).domain(), (*gEdgePmlContinuous[edgeNeighbor[e][0]].second)(jndex, index).name() + "_pair");
    //                  }
    //                  else {
    //                    gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgePmlContinuous[e].first)(index, jndex), (*gEdgePmlContinuous[e].first)(index, jndex).domain(), (*gEdgePmlContinuous[e].first)(index, jndex).name());
    //                    gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1) + (*gEdgePmlContinuous[edgeNeighbor[e][0]].first)(jndex, index), (*gEdgePmlContinuous[edgeNeighbor[e][0]].second)(jndex, index).domain(), (*gEdgePmlContinuous[edgeNeighbor[e][0]].first)(jndex, index).name() + "_pair");
    //                  }
    //                }
    //                if(!pmlBndInterface[e].second(index, jndex).isEmpty()) {
    //                  if(e >= 4 && e < 8) {
    //                    gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgePmlContinuous[e].second)(index, jndex), (*gEdgePmlContinuous[e].second)(index, jndex).domain(), (*gEdgePmlContinuous[e].second)(index, jndex).name());
    //                    gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1) + (*gEdgePmlContinuous[edgeNeighbor[e][1]].first)(jndex, index), (*gEdgePmlContinuous[edgeNeighbor[e][1]].first)(jndex, index).domain(), (*gEdgePmlContinuous[edgeNeighbor[e][1]].first)(jndex, index).name() + "_pair");
    //                  }
    //                  else {
    //                    gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*gEdgePmlContinuous[e].second)(index, jndex), (*gEdgePmlContinuous[e].second)(index, jndex).domain(), (*gEdgePmlContinuous[e].second)(index, jndex).name());
    //                    gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1) + (*gEdgePmlContinuous[edgeNeighbor[e][1]].second)(jndex, index), (*gEdgePmlContinuous[edgeNeighbor[e][1]].second)(jndex, index).domain(), (*gEdgePmlContinuous[edgeNeighbor[e][1]].second)(jndex, index).name() + "_pair");
    //                  }
    //                }
    //              }
    //            }
    //          }
    
    //          for(unsigned int c = 0; c < 8; ++c) {
    //            for(unsigned int jj = 0; jj < topology[index].size(); ++jj) {
    //              const unsigned int jndex = topology[index][jj];
    //
    //              Corner *corner = subproblem[i][j][k]->getCorner(c);
    //              if(corner == nullptr) {
    //                continue;
    //              }
    //              if(auto cr = dynamic_cast< HABC_HABC_HABC * >(corner)) {
    //                if(!std::get<0>(cornerInterface[c])(index, jndex).isEmpty()) {
    //                  for(unsigned int n = 0; n < N; ++n) {
    //                    for(unsigned int m = 0; m < N; ++m) {
    //                      gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<0>(gCornerHABC[c][n][m]))(index, jndex), (*std::get<0>(gCornerHABC[c][n][m]))(index, jndex).domain(), (*std::get<0>(gCornerHABC[c][n][m]))(index, jndex).name());
    //                      gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<1>(gCornerHABC[cornerNeighbor[c][0]][n][m]))(jndex, index), (*std::get<1>(gCornerHABC[cornerNeighbor[c][0]][n][m]))(jndex, index).domain(), (*std::get<1>(gCornerHABC[cornerNeighbor[c][0]][n][m]))(jndex, index).name() + "_pair");
    //                    }
    //                  }
    //                }
    //                if(!std::get<1>(cornerInterface[c])(index, jndex).isEmpty()) {
    //                  for(unsigned int n = 0; n < N; ++n) {
    //                    for(unsigned int m = 0; m < N; ++m) {
    //                      gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<1>(gCornerHABC[c][n][m]))(index, jndex), (*std::get<1>(gCornerHABC[c][n][m]))(index, jndex).domain(), (*std::get<1>(gCornerHABC[c][n][m]))(index, jndex).name());
    //                      gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<0>(gCornerHABC[cornerNeighbor[c][1]][n][m]))(jndex, index), (*std::get<0>(gCornerHABC[cornerNeighbor[c][1]][n][m]))(jndex, index).domain(), (*std::get<0>(gCornerHABC[cornerNeighbor[c][1]][n][m]))(jndex, index).name() + "_pair");
    //                    }
    //                  }
    //                }
    //                if(!std::get<2>(cornerInterface[c])(index, jndex).isEmpty()) {
    //                  for(unsigned int n = 0; n < N; ++n) {
    //                    for(unsigned int m = 0; m < N; ++m) {
    //                      gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<2>(gCornerHABC[c][n][m]))(index, jndex), (*std::get<2>(gCornerHABC[c][n][m]))(index, jndex).domain(), (*std::get<2>(gCornerHABC[c][n][m]))(index, jndex).name());
    //                      gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<2>(gCornerHABC[cornerNeighbor[c][2]][n][m]))(jndex, index), (*std::get<2>(gCornerHABC[cornerNeighbor[c][2]][n][m]))(jndex, index).domain(), (*std::get<2>(gCornerHABC[cornerNeighbor[c][2]][n][m]))(jndex, index).name() + "_pair");
    //                    }
    //                  }
    //                }
    //              }
    //              else if(auto cr = dynamic_cast< PmlContinuous_PmlContinuous_PmlContinuous * >(corner)) {
    //                if(!std::get<0>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
    //                  gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<0>(gCornerPmlContinuous[c]))(index, jndex), (*std::get<0>(gCornerPmlContinuous[c]))(index, jndex).domain(), (*std::get<0>(gCornerPmlContinuous[c]))(index, jndex).name());
    //                  gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<1>(gCornerPmlContinuous[cornerNeighbor[c][0]]))(jndex, index), (*std::get<1>(gCornerPmlContinuous[cornerNeighbor[c][0]]))(jndex, index).domain(), (*std::get<1>(gCornerPmlContinuous[cornerNeighbor[c][0]]))(jndex, index).name() + "_pair");
    //                }
    //                if(!std::get<1>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
    //                  gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<1>(gCornerPmlContinuous[c]))(index, jndex), (*std::get<1>(gCornerPmlContinuous[c]))(index, jndex).domain(), (*std::get<1>(gCornerPmlContinuous[c]))(index, jndex).name());
    //                  gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<0>(gCornerPmlContinuous[cornerNeighbor[c][1]]))(jndex, index), (*std::get<0>(gCornerPmlContinuous[cornerNeighbor[c][1]]))(jndex, index).domain(), (*std::get<0>(gCornerPmlContinuous[cornerNeighbor[c][1]]))(jndex, index).name() + "_pair");
    //                }
    //                if(!std::get<2>(pmlEdgeBndInterface[c])(index, jndex).isEmpty()) {
    //                  gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<2>(gCornerPmlContinuous[c]))(index, jndex), (*std::get<2>(gCornerPmlContinuous[c]))(index, jndex).domain(), (*std::get<2>(gCornerPmlContinuous[c]))(index, jndex).name());
    //                  gmshfem::post::save(gmshfem::function::ScalarFunction< std::complex< double > >(1.) + (*std::get<2>(gCornerPmlContinuous[cornerNeighbor[c][2]]))(jndex, index), (*std::get<2>(gCornerPmlContinuous[cornerNeighbor[c][2]]))(jndex, index).domain(), (*std::get<2>(gCornerPmlContinuous[cornerNeighbor[c][2]]))(jndex, index).name() + "_pair");
    //                }
    //              }
    //            }
    //          }
    
    //          exit(0);
            }
          }
        }
    
        formulation.pre();
        formulation.solve("gmres", res, iterMax);
        
        for(unsigned int i = 0; i < nDomX; ++i) {
          for(unsigned int j = 0; j < nDomY; ++j) {
            for(unsigned int k = 0; k < nDomZ; ++k) {
              unsigned int index = i * nDomY * nDomZ + j * nDomZ + k;
              if(gmshddm::mpi::isItMySubdomain(index)) {
                if(saveU) {
                  gmshfem::post::save(u(index), omega(index), "u_" + std::to_string(index), "msh");
                }
              }
            }
          }
        }
    
        if(monoDomainError) {
          gmshfem::problem::Formulation< std::complex< double > > formulationMono("HelmholtzMonoDomain");
          
          std::vector< gmshfem::domain::Domain > sigmaMono(6);
          std::vector< gmshfem::domain::Domain > pmlMono(6);
          std::vector< gmshfem::domain::Domain > pmlEdgeMono(12);
          std::vector< gmshfem::domain::Domain > pmlCornerMono(8);
          std::vector< std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > > pmlBndMono(12);
          std::vector< std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > > pmlEdgeBndMono(8);
          std::vector< gmshfem::domain::Domain > edgeMono(12);
          std::vector< gmshfem::domain::Domain > cornerMono(8);
          std::vector< std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > > cornerEdgeMono(8);
          
          for(unsigned int i = 0; i < nDomX; ++i) {
            for(unsigned int j = 0; j < nDomY; ++j) {
              for(unsigned int k = 0; k < nDomZ; ++k) {
                unsigned int index = i * nDomY * nDomZ + j * nDomZ + k;
                if(i == 0) {
                  sigmaMono[2] |= sigma[2](index);
                  pmlMono[2] |= pml[2](index);
                  if(j == 0) {
                    pmlEdgeMono[6] |= pmlEdge[6](index);
                    edgeMono[6] |= edge[6](index);
                    pmlBndMono[6].first |= pmlBnd[6].first(index);
                    pmlBndMono[6].second |= pmlBnd[6].second(index);
                    if(k == 0) {
                      pmlCornerMono[2] |= pmlCorner[2](index);
                      cornerMono[2] |= corner[2](index);
                      std::get<0>(pmlEdgeBndMono[2]) |= std::get<0>(pmlEdgeBnd[2])(index);
                      std::get<1>(pmlEdgeBndMono[2]) |= std::get<1>(pmlEdgeBnd[2])(index);
                      std::get<2>(pmlEdgeBndMono[2]) |= std::get<2>(pmlEdgeBnd[2])(index);
                      std::get<0>(cornerEdgeMono[2]) |= std::get<0>(cornerEdge[2])(index);
                      std::get<1>(cornerEdgeMono[2]) |= std::get<1>(cornerEdge[2])(index);
                      std::get<2>(cornerEdgeMono[2]) |= std::get<2>(cornerEdge[2])(index);
                    }
                    if(k == nDomZ-1) {
                      pmlCornerMono[6] |= pmlCorner[6](index);
                      cornerMono[6] |= corner[6](index);
                      std::get<0>(pmlEdgeBndMono[6]) |= std::get<0>(pmlEdgeBnd[6])(index);
                      std::get<1>(pmlEdgeBndMono[6]) |= std::get<1>(pmlEdgeBnd[6])(index);
                      std::get<2>(pmlEdgeBndMono[6]) |= std::get<2>(pmlEdgeBnd[6])(index);
                      std::get<0>(cornerEdgeMono[6]) |= std::get<0>(cornerEdge[6])(index);
                      std::get<1>(cornerEdgeMono[6]) |= std::get<1>(cornerEdge[6])(index);
                      std::get<2>(cornerEdgeMono[6]) |= std::get<2>(cornerEdge[6])(index);
                    }
                  }
                  if(k == 0) {
                    pmlEdgeMono[2] |= pmlEdge[2](index);
                    edgeMono[2] |= edge[2](index);
                    pmlBndMono[2].first |= pmlBnd[2].first(index);
                    pmlBndMono[2].second |= pmlBnd[2].second(index);
                  }
                  if(j == nDomY-1) {
                    pmlEdgeMono[5] |= pmlEdge[5](index);
                    edgeMono[5] |= edge[5](index);
                    pmlBndMono[5].first |= pmlBnd[5].first(index);
                    pmlBndMono[5].second |= pmlBnd[5].second(index);
                    if(k == 0) {
                      pmlCornerMono[1] |= pmlCorner[1](index);
                      cornerMono[1] |= corner[1](index);
                      std::get<0>(pmlEdgeBndMono[1]) |= std::get<0>(pmlEdgeBnd[1])(index);
                      std::get<1>(pmlEdgeBndMono[1]) |= std::get<1>(pmlEdgeBnd[1])(index);
                      std::get<2>(pmlEdgeBndMono[1]) |= std::get<2>(pmlEdgeBnd[1])(index);
                      std::get<0>(cornerEdgeMono[1]) |= std::get<0>(cornerEdge[1])(index);
                      std::get<1>(cornerEdgeMono[1]) |= std::get<1>(cornerEdge[1])(index);
                      std::get<2>(cornerEdgeMono[1]) |= std::get<2>(cornerEdge[1])(index);
                    }
                    if(k == nDomZ-1) {
                      pmlCornerMono[5] |= pmlCorner[5](index);
                      cornerMono[5] |= corner[5](index);
                      std::get<0>(pmlEdgeBndMono[5]) |= std::get<0>(pmlEdgeBnd[5])(index);
                      std::get<1>(pmlEdgeBndMono[5]) |= std::get<1>(pmlEdgeBnd[5])(index);
                      std::get<2>(pmlEdgeBndMono[5]) |= std::get<2>(pmlEdgeBnd[5])(index);
                      std::get<0>(cornerEdgeMono[5]) |= std::get<0>(cornerEdge[5])(index);
                      std::get<1>(cornerEdgeMono[5]) |= std::get<1>(cornerEdge[5])(index);
                      std::get<2>(cornerEdgeMono[5]) |= std::get<2>(cornerEdge[5])(index);
                    }
                  }
                  if(k == nDomZ-1) {
                    pmlEdgeMono[10] |= pmlEdge[10](index);
                    edgeMono[10] |= edge[10](index);
                    pmlBndMono[10].first |= pmlBnd[10].first(index);
                    pmlBndMono[10].second |= pmlBnd[10].second(index);
                  }
                }
                if(j == 0) {
                  sigmaMono[3] |= sigma[3](index);
                  pmlMono[3] |= pml[3](index);
                  if(k == 0) {
                    pmlEdgeMono[3] |= pmlEdge[3](index);
                    edgeMono[3] |= edge[3](index);
                    pmlBndMono[3].first |= pmlBnd[3].first(index);
                    pmlBndMono[3].second |= pmlBnd[3].second(index);
                  }
                  if(k == nDomZ-1) {
                    pmlEdgeMono[11] |= pmlEdge[11](index);
                    edgeMono[11] |= edge[11](index);
                    pmlBndMono[11].first |= pmlBnd[11].first(index);
                    pmlBndMono[11].second |= pmlBnd[11].second(index);
                  }
                }
                if(k == 0) {
                  sigmaMono[4] |= sigma[4](index);
                  pmlMono[4] |= pml[4](index);
                }
                if(i == nDomX-1) {
                  sigmaMono[0] |= sigma[0](index);
                  pmlMono[0] |= pml[0](index);
                  if(j == 0) {
                    pmlEdgeMono[7] |= pmlEdge[7](index);
                    edgeMono[7] |= edge[7](index);
                    pmlBndMono[7].first |= pmlBnd[7].first(index);
                    pmlBndMono[7].second |= pmlBnd[7].second(index);
                    if(k == 0) {
                      pmlCornerMono[3] |= pmlCorner[3](index);
                      cornerMono[3] |= corner[3](index);
                      std::get<0>(pmlEdgeBndMono[3]) |= std::get<0>(pmlEdgeBnd[3])(index);
                      std::get<1>(pmlEdgeBndMono[3]) |= std::get<1>(pmlEdgeBnd[3])(index);
                      std::get<2>(pmlEdgeBndMono[3]) |= std::get<2>(pmlEdgeBnd[3])(index);
                      std::get<0>(cornerEdgeMono[3]) |= std::get<0>(cornerEdge[3])(index);
                      std::get<1>(cornerEdgeMono[3]) |= std::get<1>(cornerEdge[3])(index);
                      std::get<2>(cornerEdgeMono[3]) |= std::get<2>(cornerEdge[3])(index);
                    }
                    if(k == nDomZ-1) {
                      pmlCornerMono[7] |= pmlCorner[7](index);
                      cornerMono[7] |= corner[7](index);
                      std::get<0>(pmlEdgeBndMono[7]) |= std::get<0>(pmlEdgeBnd[7])(index);
                      std::get<1>(pmlEdgeBndMono[7]) |= std::get<1>(pmlEdgeBnd[7])(index);
                      std::get<2>(pmlEdgeBndMono[7]) |= std::get<2>(pmlEdgeBnd[7])(index);
                      std::get<0>(cornerEdgeMono[7]) |= std::get<0>(cornerEdge[7])(index);
                      std::get<1>(cornerEdgeMono[7]) |= std::get<1>(cornerEdge[7])(index);
                      std::get<2>(cornerEdgeMono[7]) |= std::get<2>(cornerEdge[7])(index);
                    }
                  }
                  if(k == 0) {
                    pmlEdgeMono[0] |= pmlEdge[0](index);
                    edgeMono[0] |= edge[0](index);
                    pmlBndMono[0].first |= pmlBnd[0].first(index);
                    pmlBndMono[0].second |= pmlBnd[0].second(index);
                  }
                  if(j == nDomY-1) {
                    pmlEdgeMono[4] |= pmlEdge[4](index);
                    edgeMono[4] |= edge[4](index);
                    pmlBndMono[4].first |= pmlBnd[4].first(index);
                    pmlBndMono[4].second |= pmlBnd[4].second(index);
                    if(k == 0) {
                      pmlCornerMono[0] |= pmlCorner[0](index);
                      cornerMono[0] |= corner[0](index);
                      std::get<0>(pmlEdgeBndMono[0]) |= std::get<0>(pmlEdgeBnd[0])(index);
                      std::get<1>(pmlEdgeBndMono[0]) |= std::get<1>(pmlEdgeBnd[0])(index);
                      std::get<2>(pmlEdgeBndMono[0]) |= std::get<2>(pmlEdgeBnd[0])(index);
                      std::get<0>(cornerEdgeMono[0]) |= std::get<0>(cornerEdge[0])(index);
                      std::get<1>(cornerEdgeMono[0]) |= std::get<1>(cornerEdge[0])(index);
                      std::get<2>(cornerEdgeMono[0]) |= std::get<2>(cornerEdge[0])(index);
                    }
                    if(k == nDomZ-1) {
                      pmlCornerMono[4] |= pmlCorner[4](index);
                      cornerMono[4] |= corner[4](index);
                      std::get<0>(pmlEdgeBndMono[4]) |= std::get<0>(pmlEdgeBnd[4])(index);
                      std::get<1>(pmlEdgeBndMono[4]) |= std::get<1>(pmlEdgeBnd[4])(index);
                      std::get<2>(pmlEdgeBndMono[4]) |= std::get<2>(pmlEdgeBnd[4])(index);
                      std::get<0>(cornerEdgeMono[4]) |= std::get<0>(cornerEdge[4])(index);
                      std::get<1>(cornerEdgeMono[4]) |= std::get<1>(cornerEdge[4])(index);
                      std::get<2>(cornerEdgeMono[4]) |= std::get<2>(cornerEdge[4])(index);
                    }
                  }
                  if(k == nDomZ-1) {
                    pmlEdgeMono[8] |= pmlEdge[8](index);
                    edgeMono[8] |= edge[8](index);
                    pmlBndMono[8].first |= pmlBnd[8].first(index);
                    pmlBndMono[8].second |= pmlBnd[8].second(index);
                  }
                }
                if(j == nDomY-1) {
                  sigmaMono[1] |= sigma[1](index);
                  pmlMono[1] |= pml[1](index);
                  if(k == 0) {
                    pmlEdgeMono[1] |= pmlEdge[1](index);
                    edgeMono[1] |= edge[1](index);
                    pmlBndMono[1].first |= pmlBnd[1].first(index);
                    pmlBndMono[1].second |= pmlBnd[1].second(index);
                  }
                  if(k == nDomZ-1) {
                    pmlEdgeMono[9] |= pmlEdge[9](index);
                    edgeMono[9] |= edge[9](index);
                    pmlBndMono[9].first |= pmlBnd[9].first(index);
                    pmlBndMono[9].second |= pmlBnd[9].second(index);
                  }
                }
                if(k == nDomZ-1) {
                  sigmaMono[5] |= sigma[5](index);
                  pmlMono[5] |= pml[5](index);
                }
              }
            }
          }
    
          gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > uMono("uMono", omega.getUnion(), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
          gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > lambdaMono("lamdbaMono", gamma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
    
          formulationMono.integral(-grad(dof(uMono)), grad(tf(uMono)), omega.getUnion(), gauss);
          formulationMono.integral(kappa * kappa * dof(uMono), tf(uMono), omega.getUnion(), gauss);
    
          formulationMono.integral(dof(lambda), tf(uMono), gamma, gauss);
          formulationMono.integral(dof(uMono), tf(lambda), gamma, gauss);
          if(benchmark == "scattering") {
            formulationMono.integral(- fAnalytic, tf(lambda), gamma, gauss);
          }
    
          SubproblemDomains domains;
          domains.setOmega(omega.getUnion());
          domains.setSigma(sigmaMono);
          domains.setPml(pmlMono);
          domains.setPmlBnd(pmlBndMono);
          domains.setPmlEdge(pmlEdgeMono);
          domains.setEdge(edgeMono);
          domains.setPmlEdgeBnd(pmlEdgeBndMono);
          domains.setPmlCorner(pmlCornerMono);
          domains.setCorner(cornerMono);
          domains.setCornerEdge(cornerEdgeMono);
    
          SubproblemParameters parameters;
          parameters.setGauss(gauss);
          parameters.setKappa(kappa);
          parameters.setNeumannOrder(neumannOrder);
          parameters.setFieldOrder(fieldOrder);
          parameters.setStab(stab);
    
          std::string bndExt = boundaryExt;
          if(boundaryExt == "habc") {
            bndExt += "_" + std::to_string(NExt) + "_" + std::to_string(thetaPade);
          }
          else if(boundaryExt == "pml") {
            if(pmlMethodExt == "continuous") {
              bndExt += "Continuous_" + std::to_string(pmlSizeExt) + "_" + pmlTypeExt;
            }
            else if(pmlMethodExt == "discontinuous") {
              bndExt += "Discontinuous_" + std::to_string(pmlSizeExt) + "_" + pmlTypeExt;
            }
          }
          gmshfem::msg::info << "Monodomain has boundaries [" << bndExt << ", " << bndExt << ", " << bndExt << ", " << bndExt << ", " << bndExt << ", " << bndExt  << "]" << gmshfem::msg::endl;
          Subproblem subproblem(formulationMono, uMono.name(), domains, parameters, bndExt, bndExt, bndExt, bndExt, bndExt, bndExt);
          subproblem.writeFormulation();
    
          formulationMono.pre();
          formulationMono.assemble();
          formulationMono.solve();
    
          double num = 0.;
          double den = 0.;
          for(unsigned int i = 0; i < nDomX * nDomY * nDomZ; ++i) {
            num += std::real(gmshfem::post::integrate(pow(abs(u(i) - uMono), 2), omega(i), gauss));
            den += std::real(gmshfem::post::integrate(pow(abs(uMono), 2), omega(i), gauss));
          }
    
          gmshfem::msg::info << "Error mono L2 = " << std::sqrt(num/den) << gmshfem::msg::endl;
    
          for(unsigned int i = 0; i < nDomX * nDomY * nDomZ; ++i) {
            if(gmshddm::mpi::isItMySubdomain(i)) {
              if(saveEMono) {
                gmshfem::post::save(u(i) - uMono, omega(i), "eMono_" + std::to_string(i));
              }
              if(saveUMono) {
                gmshfem::post::save(uMono, omega(i), "uMono_" + std::to_string(i));
              }
            }
          }
    
          if(fileName != "none") {
            gmshfem::common::CSVio file(fileName, ';', gmshfem::common::OpeningMode::Append);
            file << iterMax << formulation.relativeResidual().back() << std::sqrt(num/den) << gmshfem::csv::endl;
            file.close();
          }
        }
    
        if(wavenumberPlot) {
          gmshfem::common::CSVio file(fileName, ';', gmshfem::common::OpeningMode::Append);
          file << k << formulation.relativeResidual().size()-1 << gmshfem::csv::endl;
          file.close();
        }
    
        if(meshPlot) {
          gmshfem::common::CSVio file(fileName, ';', gmshfem::common::OpeningMode::Append);
          file << (2*pi/k)/lc << lc << 1./lc << formulation.relativeResidual().size()-1 << gmshfem::csv::endl;
          file.close();
        }
      }
    
    
    }