Skip to content
Snippets Groups Projects
Select Git revision
  • c478e9b29a8a29f19538690022f70a5741c3956f
  • master default protected
  • steplayer
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • alphashapes
  • relaying
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 3115-issue-fix
  • 3023-Fillet2D-Update
  • convert_fdivs
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • 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
41 results

GFace.h

Blame
  • Subproblem3D.cpp 57.43 KiB
    #include "Subproblem3D.h"
    
    using gmshfem::equation::dof;
    using gmshfem::equation::tf;
    using gmshfem::function::operator-;
    using gmshfem::function::operator+;
    using gmshfem::function::operator*;
    
    namespace D3 {
    
    
      // **********************************
      // Boundary
      // **********************************
    
      Boundary::Boundary(const unsigned int boundary) : _boundary(boundary)
      {
      }
    
      Boundary::~Boundary()
      {
      }
    
      std::string Boundary::orientation() const
      {
        switch (_boundary) {
          case 0:
            return "E";
            break;
          case 1:
            return "N";
            break;
          case 2:
            return "W";
            break;
          case 3:
            return "S";
            break;
          case 4:
            return "D";
            break;
          case 5:
            return "U";
            break;
          default:
            break;
        }
        return "null";
      }
      
      // **********************************
      // Edge
      // **********************************
    
      Edge::Edge(const unsigned int edge, Boundary *first, Boundary *second) : _edge(edge), _bnd {first, second}
      {
      }
    
      Edge::~Edge()
      {
      }
    
      std::string Edge::orientation() const
      {
        switch (_edge) {
          case 0:
            return "ED";
            break;
          case 1:
            return "ND";
            break;
          case 2:
            return "WD";
            break;
          case 3:
            return "SD";
            break;
            
          case 4:
            return "EN";
            break;
          case 5:
            return "NW";
            break;
          case 6:
            return "WS";
            break;
          case 7:
            return "SE";
            break;
            
          case 8:
            return "EU";
            break;
          case 9:
            return "NU";
            break;
          case 10:
            return "WU";
            break;
          case 11:
            return "SU";
            break;
          default:
            break;
        }
        return "null";
      }
    
      // **********************************
      // Corner
      // **********************************
    
      Corner::Corner(const unsigned int corner, Edge *first, Edge *second, Edge *third) : _corner(corner), _bnd {first, second, third}
      {
      }
    
      Corner::~Corner()
      {
      }
    
      std::string Corner::orientation() const
      {
        switch (_corner) {
          case 0:
            return "END";
            break;
          case 1:
            return "NWD";
            break;
          case 2:
            return "WSD";
            break;
          case 3:
            return "SED";
            break;
            
          case 4:
            return "ENU";
            break;
          case 5:
            return "NWU";
            break;
          case 6:
            return "WSU";
            break;
          case 7:
            return "SEU";
            break;
          default:
            break;
        }
        return "null";
      }
    
      // **********************************
      // Subproblem
      // **********************************
    
      void Subproblem::_parseParameters(std::string &method, std::vector< std::string > &parameters, const std::string &str) const
      {
        method.clear();
        parameters.clear();
        std::string *current = &method;
        for(unsigned int i = 0; i < str.size(); ++i) {
          if(str[i] == '_') {
            parameters.push_back(std::string());
            current = &parameters.back();
          }
          else {
            (*current) += str[i];
          }
        }
      }
    
      Subproblem::Subproblem(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters, const std::string &E, const std::string &N, const std::string &W, const std::string &S, const std::string &D, const std::string &U) : _domains(domains), _parameters(parameters), _formulation(formulation), _fieldName(fieldName), _boundary(6, nullptr), _edge(12, nullptr), _corner(8, nullptr)
      {
        std::string type[6] = {E, N, W, S, D, U};
        std::string method[6];
        std::vector< std::string > bndParameters[6];
        for(unsigned int i = 0; i < 6; ++i) {
          _parseParameters(method[i], bndParameters[i], type[i]);
        }
        
        for(unsigned int b = 0; b < 6; ++b) {
          if(method[b] == "sommerfeld") { // str: sommerfeld
            _boundary[b] = new Sommerfeld(b);
          }
          else if(method[b] == "habc") { // str: habc_[N]_[theta]
            const unsigned int N = std::stoi(bndParameters[b][0]);
            const double theta = std::stof(bndParameters[b][1]);
            _boundary[b] = new HABC(b, N, theta);
          }
          else if(method[b] == "pmlContinuous") { // str: pmlContinuous_[N_pml]_[Type]
            const double N_pml  = std::stof(bndParameters[b][0]);
            const std::string type  = bndParameters[b][1];
            _boundary[b] = new PmlContinuous(b, N_pml, type);
          }
          // marche pas, il faudrait un vrai H(div)
          else if(method[b] == "pmlDiscontinuous") { // str: pmlDiscontinuous_[N_pml]_[Type]
            const double N_pml = std::stof(bndParameters[b][0]);
            const std::string type  = bndParameters[b][1];
            _boundary[b] = new PmlDiscontinuous(b, N_pml, type);
          }
          else {
            gmshfem::msg::error << "Unknown method '" << method[b] << "'" << gmshfem::msg::endl;
          }
        }
        
        for(unsigned int e = 0; e < 12; ++e) {
          if(e < 4) {
            if(method[e%4] == "habc" && method[4] == "habc") {
              _edge[e] = new HABC_HABC(e, static_cast< HABC * >(_boundary[e%4]), static_cast< HABC * >(_boundary[4]));
            }
            else if(method[e%4] == "pmlContinuous" && method[4] == "pmlContinuous") {
              _edge[e] = new PmlContinuous_PmlContinuous(e, static_cast< PmlContinuous * >(_boundary[e%4]), static_cast< PmlContinuous * >(_boundary[4]));
            }
          }
          else if(e >= 4 && e < 8) {
            if(method[e%4] == "habc" && method[(e+1)%4] == "habc") {
              _edge[e] = new HABC_HABC(e, static_cast< HABC * >(_boundary[e%4]), static_cast< HABC * >(_boundary[(e+1)%4]));
            }
            else if(method[e%4] == "pmlContinuous" && method[(e+1)%4] == "pmlContinuous") {
              _edge[e] = new PmlContinuous_PmlContinuous(e, static_cast< PmlContinuous * >(_boundary[e%4]), static_cast< PmlContinuous * >(_boundary[(e+1)%4]));
            }
          }
          else {
            if(method[e%4] == "habc" && method[5] == "habc") {
              _edge[e] = new HABC_HABC(e, static_cast< HABC * >(_boundary[e%4]), static_cast< HABC * >(_boundary[5]));
            }
            else if(method[e%4] == "pmlContinuous" && method[5] == "pmlContinuous") {
              _edge[e] = new PmlContinuous_PmlContinuous(e, static_cast< PmlContinuous * >(_boundary[e%4]), static_cast< PmlContinuous * >(_boundary[5]));
            }
          }
        }
        
        for(unsigned int c = 0; c < 8; ++c) {
          if(c < 4) {
            if(method[c%4] == "habc" && method[(c+1)%4] == "habc" && method[4] == "habc") {
              _corner[c] = new HABC_HABC_HABC(c, static_cast< HABC_HABC * >(_edge[c%4]), static_cast< HABC_HABC * >(_edge[(c+1)%4]), static_cast< HABC_HABC * >(_edge[c+4]));
            }
            else if(method[c%4] == "pmlContinuous" && method[(c+1)%4] == "pmlContinuous" && method[4] == "pmlContinuous") {
              _corner[c] = new PmlContinuous_PmlContinuous_PmlContinuous(c, static_cast< PmlContinuous_PmlContinuous * >(_edge[c%4]), static_cast< PmlContinuous_PmlContinuous * >(_edge[(c+1)%4]), static_cast< PmlContinuous_PmlContinuous * >(_edge[c+4]));
            }
          }
          else {
            if(method[c%4] == "habc" && method[(c+1)%4] == "habc" && method[5] == "habc") {
              _corner[c] = new HABC_HABC_HABC(c, static_cast< HABC_HABC * >(_edge[c%4+8]), static_cast< HABC_HABC * >(_edge[(c+1)%4+8]), static_cast< HABC_HABC * >(_edge[c]));
            }
            else if(method[c%4] == "pmlContinuous" && method[(c+1)%4] == "pmlContinuous" && method[5] == "pmlContinuous") {
              _corner[c] = new PmlContinuous_PmlContinuous_PmlContinuous(c, static_cast< PmlContinuous_PmlContinuous * >(_edge[c%4+8]), static_cast< PmlContinuous_PmlContinuous * >(_edge[(c+1)%4+8]), static_cast< PmlContinuous_PmlContinuous * >(_edge[c]));
            }
          }
        }
      }
    
      Subproblem::~Subproblem()
      {
        for(unsigned int c = 0; c < 8; ++c) {
          if(_corner[c] != nullptr) {
            delete _corner[c];
          }
        }
        for(unsigned int e = 0; e < 12; ++e) {
          if(_edge[e] != nullptr) {
            delete _edge[e];
          }
        }
        for(unsigned int b = 0; b < 6; ++b) {
          if(_boundary[b] != nullptr) {
            delete _boundary[b];
          }
        }
      }
    
      void Subproblem::writeFormulation()
      {
        for(unsigned int b = 0; b < 6; ++b) {
          _boundary[b]->writeFormulation(_formulation, _fieldName, _domains, _parameters);
        }
        
        for(unsigned int e = 0; e < 12; ++e) {
          if(_edge[e] != nullptr) {
            _edge[e]->writeFormulation(_formulation, _domains, _parameters);
          }
        }
        
        for(unsigned int c = 0; c < 8; ++c) {
          if(_corner[c] != nullptr) {
            _corner[c]->writeFormulation(_formulation, _domains, _parameters);
          }
        }
      }
    
      Boundary *Subproblem::getBoundary(const unsigned int b) const
      {
        return _boundary[b];
      }
      
      Edge *Subproblem::getEdge(const unsigned int e) const
      {
        return _edge[e];
      }
    
      Corner *Subproblem::getCorner(const unsigned int c) const
      {
        return _corner[c];
      }
    
      // **********************************
      // Sommerfeld
      // **********************************
    
      Sommerfeld::Sommerfeld(const unsigned int boundary) : Boundary(boundary)
      {
      }
    
      Sommerfeld::~Sommerfeld()
      {
        delete _v;
      }
    
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* Sommerfeld::getV() const
      {
        return _v;
      }
    
      void Sommerfeld::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters)
      {
        gmshfem::domain::Domain sigma = domains.getSigma(_boundary);
        
        const unsigned int neumannOrder = parameters.getNeumannOrder();
        const std::string gauss = parameters.getGauss();
        const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
        
        _v = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("v_" + orientation(), sigma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
        
        gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *u = static_cast< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >(formulation.getField(fieldName));
          
        const std::complex< double > im(0., 1.);
        formulation.integral(-dof(*_v), tf(*u), sigma, gauss);
    
        formulation.integral(dof(*_v), tf(*_v), sigma, gauss);
        formulation.integral(im * kappa * dof(*u), tf(*_v), sigma, gauss);
      }
      
      // **********************************
      // HABC
      // **********************************
    
      HABC::HABC(const unsigned int boundary, const unsigned int N, const double theta) : Boundary(boundary), _N(N), _theta(theta)
      {
      }
    
      HABC::~HABC()
      {
        delete _v;
        for(unsigned int m = 0; m < _N; ++m) {
          delete _uHABC[m];
        }
      }
    
      unsigned int HABC::getN() const
      {
        return _N;
      }
    
      double HABC::getTheta() const
      {
        return _theta;
      }
    
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * HABC::getUHABC(const unsigned int m) const
      {
        if(m > _N) throw gmshfem::common::Exception("Try to get uHABC field " + std::to_string(m));
        return _uHABC[m];
      }
    
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* HABC::getV() const
      {
        return _v;
      }
    
      void HABC::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters)
      {
        gmshfem::domain::Domain sigma = domains.getSigma(_boundary);
        
        const unsigned int neumannOrder = parameters.getNeumannOrder();
        const unsigned int fieldOrder = parameters.getFieldOrder();
        const std::string gauss = parameters.getGauss();
        const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
        
        _uHABC.resize(_N);
        for(unsigned int m = 0; m < _N; ++m) {
          _uHABC[m] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("uHABC_" + orientation() + "_" + std::to_string(m), sigma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
        }
        
        _v = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("v_" + orientation(), sigma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
        
        gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *u = static_cast< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >(formulation.getField(fieldName));
        
        const double pi = 3.141592653589793238462643383279;
        const std::complex< double > expPTheta(std::cos(_theta), std::sin(_theta));
        const std::complex< double > expPTheta2(std::cos(_theta/2.), std::sin(_theta/2.));
        const double M = 2 * _N + 1;
        double cPade[_N];
        for(unsigned int m = 0; m < _N; ++m) {
          cPade[m] = std::tan((m+1) * pi / M) * std::tan((m+1) * pi / M);
        }
    
        const std::complex< double > im(0., 1.);
        
        formulation.integral(-dof(*_v), tf(*u), sigma, gauss);
          
        formulation.integral( dof(*_v), tf(*_v), sigma, gauss);
        formulation.integral( im * kappa * expPTheta2 * dof(*u), tf(*_v), sigma, gauss);
        
        for(unsigned int m = 0; m < _N; ++m) {
          formulation.integral( im * kappa * expPTheta2 * 2./M * cPade[m] * dof(*u), tf(*_v), sigma, gauss);
          formulation.integral( im * kappa * expPTheta2 * 2./M * cPade[m] * dof(*_uHABC[m]), tf(*_v), sigma, gauss);
        }
        
        // auxiliary
        for(unsigned int m = 0; m < _N; ++m) {
          formulation.integral(grad(dof(*_uHABC[m])), grad(tf(*_uHABC[m])), sigma, gauss);
          formulation.integral(- kappa * kappa * (expPTheta * cPade[m] + 1.) * dof(*_uHABC[m]), tf(*_uHABC[m]), sigma, gauss);
          formulation.integral(- kappa * kappa * expPTheta * (cPade[m] + 1.) * dof(*u), tf(*_uHABC[m]), sigma, gauss);
        }
      }
      
      // **********************************
      // PML
      // **********************************
    
      Pml::Pml(const unsigned int boundary, const double pmlSize, const std::string type) : Boundary(boundary), _size(pmlSize), _type(type)
      {
      }
    
      Pml::~Pml()
      {
      }
    
      double Pml::getSize() const
      {
        return _size;
      }
    
      gmshfem::function::ScalarFunction< std::complex< double > > Pml::pmlCoefficients(gmshfem::function::TensorFunction< std::complex< double > > &D, gmshfem::function::ScalarFunction< std::complex< double > > &E,  gmshfem::domain::Domain &bnd, const gmshfem::function::ScalarFunction< std::complex< double > > &kappa) const
      {
        gmshfem::function::ScalarFunction< std::complex< double > > kappaMod;
        const std::complex< double > im(0., 1.);
    
        gmshfem::function::ScalarFunction< std::complex< double > > distSigma;
        gmshfem::function::ScalarFunction< std::complex< double > > sigma;
        if(orientation() == "E" || orientation() == "W") {
          distSigma = abs(gmshfem::function::x< std::complex< double > >() - gmshfem::function::x< std::complex< double > >(bnd));
          if(bnd.maxDim() == 0 || bnd.maxDim() == 1) {
            kappaMod = changeOfCoordinates(kappa, gmshfem::function::x< double >(bnd), gmshfem::function::y< double >(bnd), gmshfem::function::z< double >(bnd));
          }
          else {
            kappaMod = changeOfCoordinates(kappa, gmshfem::function::x< double >(bnd), gmshfem::function::y< double >(), gmshfem::function::z< double >());
          }
        }
        else if(orientation() == "N" || orientation() == "S") {
          distSigma = abs(gmshfem::function::y< std::complex< double > >() - gmshfem::function::y< std::complex< double > >(bnd));
          if(bnd.maxDim() == 0 || bnd.maxDim() == 1) {
            kappaMod = changeOfCoordinates(kappa, gmshfem::function::x< double >(bnd), gmshfem::function::y< double >(bnd), gmshfem::function::z< double >(bnd));
          }
          else {
            kappaMod = changeOfCoordinates(kappa, gmshfem::function::x< double >(), gmshfem::function::y< double >(bnd), gmshfem::function::z< double >());
          }
        }
        else if(orientation() == "D" || orientation() == "U") {
          distSigma = abs(gmshfem::function::z< std::complex< double > >() - gmshfem::function::z< std::complex< double > >(bnd));
          if(bnd.maxDim() == 0 || bnd.maxDim() == 1) {
            kappaMod = changeOfCoordinates(kappa, gmshfem::function::x< double >(bnd), gmshfem::function::y< double >(bnd), gmshfem::function::z< double >(bnd));
          }
          else {
            kappaMod = changeOfCoordinates(kappa, gmshfem::function::x< double >(), gmshfem::function::y< double >(), gmshfem::function::z< double >(bnd));
          }
        }
        
        if(_type == "hs") {
          sigma = 1. / (_size - distSigma) - 1./ _size;
        }
        else if(_type == "h") {
          sigma = 1. / (_size - distSigma);
        }
        else if(_type == "q") {
          const double sigmaStar = 664.55;
          sigma = sigmaStar * distSigma * distSigma / (_size * _size);
        }
        
        gmshfem::function::ScalarFunction< std::complex< double > > kMod = 1. + im * sigma / kappaMod;
        if(orientation() == "E" || orientation() == "W") {
          D = gmshfem::function::tensorDiag< std::complex< double > >(1./kMod, kMod, kMod);
          E = kMod;
        }
        else if(orientation() == "N" || orientation() == "S") {
          D = gmshfem::function::tensorDiag< std::complex< double > >(kMod, 1./kMod, kMod);
          E = kMod;
        }
        else if(orientation() == "D" || orientation() == "U") {
          D = gmshfem::function::tensorDiag< std::complex< double > >(kMod, kMod, 1./kMod);
          E = kMod;
        }
        
        return kappaMod;
      }
      
      // **********************************
      // PML continuous
      // **********************************
    
      PmlContinuous::PmlContinuous(const unsigned int boundary, const double pmlSize, const std::string type) : Pml(boundary, pmlSize, type)
      {
      }
    
      PmlContinuous::~PmlContinuous()
      {
        delete _v;
        delete _uPml;
      }
    
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * PmlContinuous::getUPml() const
      {
        return _uPml;
      }
    
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * PmlContinuous::getV() const
      {
        return _v;
      }
    
      void PmlContinuous::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters)
      {
        gmshfem::domain::Domain sigma = domains.getSigma(_boundary);
        gmshfem::domain::Domain pml = domains.getPml(_boundary);
        
        const unsigned int neumannOrder = parameters.getNeumannOrder();
        const unsigned int fieldOrder = parameters.getFieldOrder();
        const std::string gauss = parameters.getGauss();
        const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
        
        _uPml = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("uPml_" + orientation(), pml, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
        _v = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("v_" + orientation(), sigma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
        
        gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *u = static_cast< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >(formulation.getField(fieldName));
        
        gmshfem::function::TensorFunction< std::complex< double > > D;
        gmshfem::function::ScalarFunction< std::complex< double > > E;
        gmshfem::function::ScalarFunction< std::complex< double > > kappaMod = pmlCoefficients(D, E, sigma, kappa);
            
        formulation.integral(- D * grad(dof(*_uPml)), grad(tf(*_uPml)), pml, gauss);
        formulation.integral(kappaMod * kappaMod * E * dof(*_uPml), tf(*_uPml), pml, gauss);
          
        formulation.integral(dof(*_v), tf(*_uPml), sigma, gauss);
        formulation.integral(- dof(*_v), tf(*u), sigma, gauss);
          
        formulation.integral(dof(*_uPml), tf(*_v), sigma, gauss);
        formulation.integral(- dof(*u), tf(*_v), sigma, gauss);
      }
      
      // **********************************
      // PML discontinuous
      // **********************************
    
      PmlDiscontinuous::PmlDiscontinuous(const unsigned int boundary, const double pmlSize, const std::string type) : Pml(boundary, pmlSize, type)
      {
      }
    
      PmlDiscontinuous::~PmlDiscontinuous()
      {
        delete _v;
        delete _uPml;
      }
    
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * PmlDiscontinuous::getUPml() const
      {
        return _uPml;
      }
    
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form2 > * PmlDiscontinuous::getV() const
      {
        return _v;
      }
    
      void PmlDiscontinuous::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters)
      {
        // Marche pas, il faudrait du vrai H(div) ...
    //    gmshfem::domain::Domain sigma = domains.getSigma(_boundary);
    //    gmshfem::domain::Domain pml = domains.getPml(_boundary);
    //    
    //    const unsigned int neumannOrder = parameters.getNeumannOrder();
    //    const unsigned int fieldOrder = parameters.getFieldOrder();
    //    const std::string gauss = parameters.getGauss();
    //    const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
    //    
    //    _uPml = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("uPml_" + orientation(), pml, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
    //   
    //    _v = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form2 >("v_" + orientation(), sigma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
    //    
    //    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *u = static_cast< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >(formulation.getField(fieldName));
    //
    //    gmshfem::function::TensorFunction< std::complex< double > > D;
    //    gmshfem::function::ScalarFunction< std::complex< double > > E;
    //    gmshfem::function::ScalarFunction< std::complex< double > > kappaMod = pmlCoefficients(D, E, sigma, kappa);
    //
    //    formulation.integral(- D * grad(dof(*_uPml)), grad(tf(*_uPml)), pml, gauss);
    //    formulation.integral(kappaMod * kappaMod * E * dof(*_uPml), tf(*_uPml), pml, gauss);
    //
    //    formulation.integral(dof(*_v), tf(*_uPml), sigma, gauss);
    //    formulation.integral(- dof(*_v), tf(*u), sigma, gauss);
    //
    //    formulation.integral(dof(*_uPml), tf(*_v), sigma, gauss);
    //    formulation.integral(- dof(*u), tf(*_v), sigma, gauss);
      }
      
      // **********************************
      // PML without multiplier
      // **********************************
    
      PmlWithoutMultiplier::PmlWithoutMultiplier(const unsigned int boundary, const double pmlSize, const std::string type) : Pml(boundary, pmlSize, type)
      {
      }
    
      PmlWithoutMultiplier::~PmlWithoutMultiplier()
      {
        delete _v;
        delete _uPml;
      }
    
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * PmlWithoutMultiplier::getUPml() const
      {
        return _uPml;
      }
    
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * PmlWithoutMultiplier::getV() const
      {
        return _v;
      }
    
      void PmlWithoutMultiplier::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters)
      {
        gmshfem::domain::Domain sigma = domains.getSigma(_boundary);
        gmshfem::domain::Domain pml = domains.getPml(_boundary);
        
        const unsigned int neumannOrder = parameters.getNeumannOrder();
        const unsigned int fieldOrder = parameters.getFieldOrder();
        const std::string gauss = parameters.getGauss();
        const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
        
        _uPml = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("uPml_" + orientation(), pml, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
        _v = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("v_" + orientation(), sigma, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
        
        gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *u = static_cast< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >(formulation.getField(fieldName));
        
        gmshfem::function::TensorFunction< std::complex< double > > D;
        gmshfem::function::ScalarFunction< std::complex< double > > E;
        gmshfem::function::ScalarFunction< std::complex< double > > kappaMod = pmlCoefficients(D, E, sigma, kappa);
            
        formulation.integral(- D * grad(dof(*_uPml)), grad(tf(*_uPml)), pml, gauss);
        formulation.integral(kappaMod * kappaMod * E * dof(*_uPml), tf(*_uPml), pml, gauss);
          
        formulation.integral(dof(*_v), tf(*_uPml), sigma, gauss);
        formulation.integral(- dof(*_v), tf(*u), sigma, gauss);
          
        formulation.integral(dof(*_uPml), tf(*_v), sigma, gauss);
        formulation.integral(- dof(*u), tf(*_v), sigma, gauss);
      }
      
    
    
    
      // **********************************
      // HABC_HABC
      // **********************************
    
      HABC_HABC::HABC_HABC(const unsigned int edge, HABC *first, HABC *second) : Edge(edge, first, second)
      {
      }
    
      HABC_HABC::~HABC_HABC()
      {
        const unsigned int N0 = static_cast< HABC * >(_bnd[0])->getN();
        const unsigned int N1 = static_cast< HABC * >(_bnd[1])->getN();
        for(unsigned int m = 0; m < N0; ++m) {
          for(unsigned int n = 0; n < N1; ++n) {
            delete _uHABC_E[m][n];
          }
        }
        for(unsigned int m = 0; m < N0; ++m) {
          delete _vE[0][m];
        }
        for(unsigned int n = 0; n < N1; ++n) {
          delete _vE[1][n];
        }
      }
      
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *HABC_HABC::getUHABC(const unsigned int m, const unsigned int n) const
      {
        return _uHABC_E[m][n];
      }
      
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *HABC_HABC::getV(const unsigned int m, const unsigned int i) const
      {
        return _vE[i][m];
      }
      
      HABC *HABC_HABC::firstBoundary() const
      {
        return static_cast< HABC * >(_bnd[0]);
      }
    
      HABC *HABC_HABC::secondBoundary() const
      {
        return static_cast< HABC * >(_bnd[1]);
      }
    
      void HABC_HABC::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters)
      {
        gmshfem::domain::Domain edge = domains.getEdge(_edge);
        
        const unsigned int neumannOrder = parameters.getNeumannOrder();
        const unsigned int fieldOrder = parameters.getFieldOrder();
        const std::string gauss = parameters.getGauss();
        const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
        
        const unsigned int N0 = static_cast< HABC * >(_bnd[0])->getN();
        const unsigned int N1 = static_cast< HABC * >(_bnd[1])->getN();
        
        const double theta0 = static_cast< HABC * >(_bnd[0])->getTheta();
        const double theta1 = static_cast< HABC * >(_bnd[1])->getTheta();
      
        _uHABC_E.resize(N0);
        for(unsigned int m = 0; m < N0; ++m) {
          _uHABC_E[m].resize(N1);
          for(unsigned int n = 0; n < N1; ++n) {
            _uHABC_E[m][n] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("uHABC_Edge_" + orientation() + "_" + std::to_string(m) + "_" + std::to_string(n), edge, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
          }
        }
        
        _vE[0].resize(N0);
        for(unsigned int m = 0; m < N0; ++m) {
          _vE[0][m] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vE_" + orientation() + "_" + std::to_string(m) + "_first", edge, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
        }
        _vE[1].resize(N1);
        for(unsigned int n = 0; n < N1; ++n) {
          _vE[1][n] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vE_" + orientation() + "_" + std::to_string(n) + "_second", edge, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
        }
        
        const double pi = 3.141592653589793238462643383279;
        
        const std::complex< double > expPTheta_0(std::cos(theta0), std::sin(theta0));
        const std::complex< double > expPTheta2_0(std::cos(theta0/2.), std::sin(theta0/2.));
        
        const std::complex< double > expPTheta_1(std::cos(theta1), std::sin(theta1));
        const std::complex< double > expPTheta2_1(std::cos(theta1/2.), std::sin(theta1/2.));
        
        const double M0 = 2 * N0 + 1;
        double cPade0[N0];
        for(unsigned int m = 0; m < N0; ++m) {
          cPade0[m] = std::tan((m+1) * pi / M0) * std::tan((m+1) * pi / M0);
        }
        
        const double M1 = 2 * N1 + 1;
        double cPade1[N1];
        for(unsigned int n = 0; n < N1; ++n) {
          cPade1[n] = std::tan((n+1) * pi / M1) * std::tan((n+1) * pi / M1);
        }
    
        const std::complex< double > im(0., 1.);
        
        std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * > uHABC[2];
        uHABC[0].resize(N0);
        uHABC[1].resize(N1);
        
        for(unsigned int m = 0; m < N0; ++m) {
          uHABC[0][m] = static_cast< HABC * >(_bnd[0])->getUHABC(m);
          formulation.integral(-dof(*_vE[0][m]), tf(*uHABC[0][m]), edge, gauss);
    
          formulation.integral( dof(*_vE[0][m]), tf(*_vE[0][m]), edge, gauss);
          formulation.integral(-im * kappa * expPTheta2_1 * dof(*uHABC[0][m]), tf(*_vE[0][m]), edge, gauss);
          
          for(unsigned int n = 0; n < N1; ++n) {
            formulation.integral(-im * kappa * expPTheta2_1 * 2./M1 * cPade1[n] * dof(*uHABC[0][m]), tf(*_vE[0][m]), edge, gauss);
            formulation.integral(-im * kappa * expPTheta2_1 * 2./M1 * cPade1[n] * dof(*_uHABC_E[m][n]), tf(*_vE[0][m]), edge, gauss);
          }
        }
        
        for(unsigned int n = 0; n < N1; ++n) {
          uHABC[1][n] = static_cast< HABC * >(_bnd[1])->getUHABC(n);
          formulation.integral(-dof(*_vE[1][n]), tf(*uHABC[1][n]), edge, gauss);
    
          formulation.integral( dof(*_vE[1][n]), tf(*_vE[1][n]), edge, gauss);
          formulation.integral(-im * kappa * expPTheta2_0 * dof(*uHABC[1][n]), tf(*_vE[1][n]), edge, gauss);
          
          for(unsigned int m = 0; m < N0; ++m) {
            formulation.integral(-im * kappa * expPTheta2_0 * 2./M0 * cPade0[m] * dof(*uHABC[1][n]), tf(*_vE[1][n]), edge, gauss);
            formulation.integral(-im * kappa * expPTheta2_0 * 2./M0 * cPade0[m] * dof(*_uHABC_E[m][n]), tf(*_vE[1][n]), edge, gauss);
          }
        }
    
        // auxiliary
        for(unsigned int m = 0; m < N0; ++m) {
          for(unsigned int n = 0; n < N1; ++n) {
            formulation.integral(grad(dof(*_uHABC_E[m][n])), grad(tf(*_uHABC_E[m][n])), edge, gauss);
            formulation.integral(- kappa * kappa * (expPTheta_0 * cPade0[m] + expPTheta_1 * cPade1[n] + 1.) * dof(*_uHABC_E[m][n]), tf(*_uHABC_E[m][n]), edge, gauss);
            formulation.integral(- kappa * kappa * expPTheta_1 * (cPade1[n] + 1.) * dof(*uHABC[0][m]), tf(*_uHABC_E[m][n]), edge, gauss);
            formulation.integral(- kappa * kappa * expPTheta_0 * (cPade0[m] + 1.) * dof(*uHABC[1][n]), tf(*_uHABC_E[m][n]), edge, gauss);
          }
        }
      }
      
      // **********************************
      // PML continuous _ PML continuous
      // **********************************
    
      PmlContinuous_PmlContinuous::PmlContinuous_PmlContinuous(const unsigned int edge, PmlContinuous *first, PmlContinuous *second) : Edge(edge, first, second)
      {
      }
    
      PmlContinuous_PmlContinuous::~PmlContinuous_PmlContinuous()
      {
        delete _uPmlEdge;
        delete _vC[0];
        delete _vC[1];
        delete _vEdge;
      }
    
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* PmlContinuous_PmlContinuous::getUPml() const
      {
        return _uPmlEdge;
      }
    
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* PmlContinuous_PmlContinuous::getV(const unsigned int i) const
      {
        return _vC[i];
      }
      
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* PmlContinuous_PmlContinuous::getVEdge() const
      {
        return _vEdge;
      }
    
      PmlContinuous *PmlContinuous_PmlContinuous::firstBoundary() const
      {
        return static_cast< PmlContinuous * >(_bnd[0]);
      }
    
      PmlContinuous *PmlContinuous_PmlContinuous::secondBoundary() const
      {
        return static_cast< PmlContinuous * >(_bnd[1]);
      }
    
      void PmlContinuous_PmlContinuous::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters)
      {
        gmshfem::domain::Domain edge = domains.getEdge(_edge);
        gmshfem::domain::Domain pmlEdge = domains.getPmlEdge(_edge);
        std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > pmlBnd = domains.getPmlBnd(_edge);
    
        const unsigned int neumannOrder = parameters.getNeumannOrder();
        const unsigned int fieldOrder = parameters.getFieldOrder();
        const std::string gauss = parameters.getGauss();
        const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
    
        _uPmlEdge = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("uPmlEdge_" + orientation(), pmlEdge, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
    
        _vC[0] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_first", pmlBnd.first, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
        _vC[1] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_second", pmlBnd.second, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
    
        gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *uPml[2];
        uPml[0] = static_cast< PmlContinuous * >(_bnd[0])->getUPml();
        uPml[1] = static_cast< PmlContinuous * >(_bnd[1])->getUPml();
    
        gmshfem::function::TensorFunction< std::complex< double > > D[2];
        gmshfem::function::ScalarFunction< std::complex< double > > E[2];
        gmshfem::function::ScalarFunction< std::complex< double > > kappaMod = static_cast< PmlContinuous * >(_bnd[0])->pmlCoefficients(D[0], E[0], edge, kappa);
        static_cast< PmlContinuous * >(_bnd[1])->pmlCoefficients(D[1], E[1], edge, kappa);
    
        formulation.integral(- D[0] * D[1] * grad(dof(*_uPmlEdge)), grad(tf(*_uPmlEdge)), pmlEdge, gauss);
        formulation.integral(kappaMod * kappaMod * E[0] * E[1] * dof(*_uPmlEdge), tf(*_uPmlEdge), pmlEdge, gauss);
    
        formulation.integral(dof(*_vC[0]), tf(*_uPmlEdge), pmlBnd.first, gauss);
        formulation.integral(- dof(*_vC[0]), tf(*uPml[0]), pmlBnd.first, gauss);
    
        formulation.integral(dof(*_uPmlEdge), tf(*_vC[0]), pmlBnd.first, gauss);
        formulation.integral(- dof(*uPml[0]), tf(*_vC[0]), pmlBnd.first, gauss);
    
        formulation.integral(dof(*_vC[1]), tf(*_uPmlEdge), pmlBnd.second, gauss);
        formulation.integral(- dof(*_vC[1]), tf(*uPml[1]), pmlBnd.second, gauss);
    
        formulation.integral(dof(*_uPmlEdge), tf(*_vC[1]), pmlBnd.second, gauss);
        formulation.integral(- dof(*uPml[1]), tf(*_vC[1]), pmlBnd.second, gauss);
        
        // edge equation
        _vEdge = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vEdge_" + orientation(), edge, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
        
        formulation.integral(dof(*_vC[0]), tf(*_vEdge), edge, gauss);
        formulation.integral(dof(*_vEdge), tf(*_vC[0]), edge, gauss);
    
        formulation.integral(-dof(*_vC[1]), tf(*_vEdge), edge, gauss);
        formulation.integral(-dof(*_vEdge), tf(*_vC[1]), edge, gauss);
    
        formulation.integral(dof(*static_cast< PmlContinuous * >(_bnd[0])->getV()), tf(*_vEdge), edge, gauss);
        formulation.integral(dof(*_vEdge), tf(*static_cast< PmlContinuous * >(_bnd[0])->getV()), edge, gauss);
    
        formulation.integral(-dof(*static_cast< PmlContinuous * >(_bnd[1])->getV()), tf(*_vEdge), edge, gauss);
        formulation.integral(-dof(*_vEdge), tf(*static_cast< PmlContinuous * >(_bnd[1])->getV()), edge, gauss);
      }
      
    
    
    
      // **********************************
      // HABC_HABC_HABC
      // **********************************
    
      HABC_HABC_HABC::HABC_HABC_HABC(const unsigned int corner, HABC_HABC *first, HABC_HABC *second, HABC_HABC *third) : Corner(corner, first, second, third)
      {
      }
    
      HABC_HABC_HABC::~HABC_HABC_HABC()
      {
        std::vector< HABC * > boundaries(3, nullptr);
        std::map< std::string, HABC * > boundaryMap;
        for(unsigned int i = 0; i < 3; ++i) {
          boundaryMap.insert(std::pair< std::string, HABC * >(static_cast< HABC_HABC * >(_bnd[i])->firstBoundary()->orientation(), static_cast< HABC_HABC * >(_bnd[i])->firstBoundary()));
          boundaryMap.insert(std::pair< std::string, HABC * >(static_cast< HABC_HABC * >(_bnd[i])->secondBoundary()->orientation(), static_cast< HABC_HABC * >(_bnd[i])->secondBoundary()));
        }
        
        for(auto it = boundaryMap.begin(); it != boundaryMap.end(); ++it) {
          if(it->first == "D" || it->first == "U") {
            boundaries[2] = it->second;
          }
          else {
            if(boundaries[0] == nullptr) {
              boundaries[0] = it->second;
            }
            else {
              boundaries[1] = it->second;
            }
          }
        }
        
        if((boundaries[0]->orientation() == "N" && boundaries[1]->orientation() == "E") ||
           (boundaries[0]->orientation() == "W" && boundaries[1]->orientation() == "N") ||
           (boundaries[0]->orientation() == "S" && boundaries[1]->orientation() == "W") ||
           (boundaries[0]->orientation() == "E" && boundaries[1]->orientation() == "S")) {
          std::swap(boundaries[0], boundaries[1]);
        }
        
        const unsigned int N0 = static_cast< HABC * >(boundaries[0])->getN();
        const unsigned int N1 = static_cast< HABC * >(boundaries[1])->getN();
        const unsigned int N2 = static_cast< HABC * >(boundaries[2])->getN();
        
        for(unsigned int m = 0; m < N0; ++m) {
          for(unsigned int n = 0; n < N1; ++n) {
            for(unsigned int o = 0; o < N2; ++o) {
              delete _uHABC_C[m][n][o];
            }
          }
        }
        
        for(unsigned int m = 0; m < N0; ++m) {
          for(unsigned int n = 0; n < N2; ++n) {
            delete _vC[0][m][n];
          }
        }
        for(unsigned int n = 0; n < N1; ++n) {
          for(unsigned int o = 0; o < N2; ++o) {
            delete _vC[1][n][o];
          }
        }
        for(unsigned int m = 0; m < N0; ++m) {
          for(unsigned int n = 0; n < N1; ++n) {
            delete _vC[2][m][n];
          }
        }
      }
      
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *HABC_HABC_HABC::getUHABC(const unsigned int m, const unsigned int n, const unsigned int o) const
      {
        return _uHABC_C[m][n][o];
      }
      
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *HABC_HABC_HABC::getV(const unsigned int m, const unsigned int n, const unsigned int i) const
      {
        return _vC[i][m][n];
      }
      
      HABC_HABC *HABC_HABC_HABC::firstEdge() const
      {
        return static_cast< HABC_HABC * >(_bnd[0]);
      }
    
      HABC_HABC *HABC_HABC_HABC::secondEdge() const
      {
        return static_cast< HABC_HABC * >(_bnd[1]);
      }
      
      HABC_HABC *HABC_HABC_HABC::thirdEdge() const
      {
        return static_cast< HABC_HABC * >(_bnd[2]);
      }
    
      void HABC_HABC_HABC::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters)
      {
        gmshfem::domain::Domain corner = domains.getCorner(_corner);
        
        const unsigned int neumannOrder = parameters.getNeumannOrder();
        const unsigned int fieldOrder = parameters.getFieldOrder();
        const std::string gauss = parameters.getGauss();
        const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
        
        std::vector< HABC * > boundaries(3, nullptr);
        std::map< std::string, HABC * > boundaryMap;
        for(unsigned int i = 0; i < 3; ++i) {
          boundaryMap.insert(std::pair< std::string, HABC * >(static_cast< HABC_HABC * >(_bnd[i])->firstBoundary()->orientation(), static_cast< HABC_HABC * >(_bnd[i])->firstBoundary()));
          boundaryMap.insert(std::pair< std::string, HABC * >(static_cast< HABC_HABC * >(_bnd[i])->secondBoundary()->orientation(), static_cast< HABC_HABC * >(_bnd[i])->secondBoundary()));
        }
        
        for(auto it = boundaryMap.begin(); it != boundaryMap.end(); ++it) {
          if(it->first == "D" || it->first == "U") {
            boundaries[2] = it->second;
          }
          else {
            if(boundaries[0] == nullptr) {
              boundaries[0] = it->second;
            }
            else {
              boundaries[1] = it->second;
            }
          }
        }
        
        if((boundaries[0]->orientation() == "N" && boundaries[1]->orientation() == "E") ||
           (boundaries[0]->orientation() == "W" && boundaries[1]->orientation() == "N") ||
           (boundaries[0]->orientation() == "S" && boundaries[1]->orientation() == "W") ||
           (boundaries[0]->orientation() == "E" && boundaries[1]->orientation() == "S")) {
          std::swap(boundaries[0], boundaries[1]);
        }
        
        const unsigned int N0 = static_cast< HABC * >(boundaries[0])->getN();
        const unsigned int N1 = static_cast< HABC * >(boundaries[1])->getN();
        const unsigned int N2 = static_cast< HABC * >(boundaries[2])->getN();
        
        const double theta0 = static_cast< HABC * >(boundaries[0])->getTheta();
        const double theta1 = static_cast< HABC * >(boundaries[1])->getTheta();
        const double theta2 = static_cast< HABC * >(boundaries[2])->getTheta();
        
        _uHABC_C.resize(N0);
        for(unsigned int m = 0; m < N0; ++m) {
          _uHABC_C[m].resize(N1);
          for(unsigned int n = 0; n < N1; ++n) {
            _uHABC_C[m][n].resize(N2);
            for(unsigned int o = 0; o < N2; ++o) {
              _uHABC_C[m][n][o] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("uHABC_Corner_" + orientation() + "_" + std::to_string(m) + "_" + std::to_string(n) + "_" + std::to_string(o), corner, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
            }
          }
        }
        
        _vC[0].resize(N0);
        for(unsigned int m = 0; m < N0; ++m) {
          _vC[0][m].resize(N2);
          for(unsigned int o = 0; o < N2; ++o) {
            _vC[0][m][o] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_" + std::to_string(m) + "_" + std::to_string(o) + "_first", corner, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
          }
        }
        _vC[1].resize(N1);
        for(unsigned int n = 0; n < N1; ++n) {
          _vC[1][n].resize(N2);
          for(unsigned int o = 0; o < N2; ++o) {
            _vC[1][n][o] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_" + std::to_string(n) + "_" + std::to_string(o) + "_second", corner, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
          }
        }
        _vC[2].resize(N0);
        for(unsigned int m = 0; m < N0; ++m) {
          _vC[2][m].resize(N1);
          for(unsigned int n = 0; n < N1; ++n) {
            _vC[2][m][n] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_" + std::to_string(m) + "_" + std::to_string(n) + "_third", corner, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
          }
        }
        
        const double pi = 3.141592653589793238462643383279;
        
        const std::complex< double > expPTheta_0(std::cos(theta0), std::sin(theta0));
        const std::complex< double > expPTheta2_0(std::cos(theta0/2.), std::sin(theta0/2.));
        
        const std::complex< double > expPTheta_1(std::cos(theta1), std::sin(theta1));
        const std::complex< double > expPTheta2_1(std::cos(theta1/2.), std::sin(theta1/2.));
        
        const std::complex< double > expPTheta_2(std::cos(theta2), std::sin(theta2));
        const std::complex< double > expPTheta2_2(std::cos(theta2/2.), std::sin(theta2/2.));
    
        const double M0 = 2 * N0 + 1;
        double cPade0[N0];
        for(unsigned int m = 0; m < N0; ++m) {
          cPade0[m] = std::tan((m+1) * pi / M0) * std::tan((m+1) * pi / M0);
        }
        
        const double M1 = 2 * N1 + 1;
        double cPade1[N1];
        for(unsigned int n = 0; n < N1; ++n) {
          cPade1[n] = std::tan((n+1) * pi / M1) * std::tan((n+1) * pi / M1);
        }
        
        const double M2 = 2 * N2 + 1;
        double cPade2[N2];
        for(unsigned int o = 0; o < N2; ++o) {
          cPade2[o] = std::tan((o+1) * pi / M2) * std::tan((o+1) * pi / M2);
        }
    
        const std::complex< double > im(0., 1.);
        
        std::vector< std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * > > uHABC[3];
        uHABC[0].resize(N0, std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >(N2, nullptr));
        uHABC[1].resize(N1, std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >(N2, nullptr));
        uHABC[2].resize(N0, std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >(N1, nullptr));
        
        for(unsigned int m = 0; m < N0; ++m) {
          for(unsigned int o = 0; o < N2; ++o) {
            uHABC[0][m][o] = static_cast< HABC_HABC * >(_bnd[0])->getUHABC(m,o);
            formulation.integral(-dof(*_vC[0][m][o]), tf(*uHABC[0][m][o]), corner, gauss);
    
            formulation.integral( dof(*_vC[0][m][o]), tf(*_vC[0][m][o]), corner, gauss);
            formulation.integral(-im * kappa * expPTheta2_1 * dof(*uHABC[0][m][o]), tf(*_vC[0][m][o]), corner, gauss);
            
            for(unsigned int n = 0; n < N1; ++n) {
              formulation.integral(-im * kappa * expPTheta2_1 * 2./M1 * cPade1[n] * dof(*uHABC[0][m][o]), tf(*_vC[0][m][o]), corner, gauss);
              formulation.integral(-im * kappa * expPTheta2_1 * 2./M1 * cPade1[n] * dof(*_uHABC_C[m][n][o]), tf(*_vC[0][m][o]), corner, gauss);
            }
          }
        }
        
        for(unsigned int n = 0; n < N1; ++n) {
          for(unsigned int o = 0; o < N2; ++o) {
            uHABC[1][n][o] = static_cast< HABC_HABC * >(_bnd[1])->getUHABC(n,o);
            formulation.integral(-dof(*_vC[1][n][o]), tf(*uHABC[1][n][o]), corner, gauss);
    
            formulation.integral( dof(*_vC[1][n][o]), tf(*_vC[1][n][o]), corner, gauss);
            formulation.integral(-im * kappa * expPTheta2_0 * dof(*uHABC[1][n][o]), tf(*_vC[1][n][o]), corner, gauss);
            
            for(unsigned int m = 0; m < N0; ++m) {
              formulation.integral(-im * kappa * expPTheta2_0 * 2./M0 * cPade0[m] * dof(*uHABC[1][n][o]), tf(*_vC[1][n][o]), corner, gauss);
              formulation.integral(-im * kappa * expPTheta2_0 * 2./M0 * cPade0[m] * dof(*_uHABC_C[m][n][o]), tf(*_vC[1][n][o]), corner, gauss);
            }
          }
        }
        
        for(unsigned int m = 0; m < N0; ++m) {
          for(unsigned int n = 0; n < N1; ++n) {
            uHABC[2][m][n] = static_cast< HABC_HABC * >(_bnd[2])->getUHABC(m,n);
            formulation.integral(-dof(*_vC[2][m][n]), tf(*uHABC[2][m][n]), corner, gauss);
    
            formulation.integral( dof(*_vC[2][m][n]), tf(*_vC[2][m][n]), corner, gauss);
            formulation.integral(-im * kappa * expPTheta2_2 * dof(*uHABC[2][m][n]), tf(*_vC[2][m][n]), corner, gauss);
            
            for(unsigned int o = 0; o < N2; ++o) {
              formulation.integral(-im * kappa * expPTheta2_2 * 2./M2 * cPade2[o] * dof(*uHABC[2][m][n]), tf(*_vC[2][m][n]), corner, gauss);
              formulation.integral(-im * kappa * expPTheta2_2 * 2./M2 * cPade2[o] * dof(*_uHABC_C[m][n][o]), tf(*_vC[2][m][n]), corner, gauss);
            }
          }
        }
    
        // auxiliary
        for(unsigned int m = 0; m < N0; ++m) {
          for(unsigned int n = 0; n < N1; ++n) {
            for(unsigned int o = 0; o < N2; ++o) {
              formulation.integral( (expPTheta_0 * cPade0[m] + expPTheta_1 * cPade1[n] + expPTheta_2 * cPade2[o] + 1.) * dof(*_uHABC_C[m][n][o]), tf(*_uHABC_C[m][n][o]), corner, gauss);
              formulation.integral( expPTheta_1 * (cPade1[n] + 1.) * dof(*uHABC[0][m][o]), tf(*_uHABC_C[m][n][o]), corner, gauss);
              formulation.integral( expPTheta_0 * (cPade0[m] + 1.) * dof(*uHABC[1][n][o]), tf(*_uHABC_C[m][n][o]), corner, gauss);
              formulation.integral( expPTheta_2 * (cPade2[o] + 1.) * dof(*uHABC[2][m][n]), tf(*_uHABC_C[m][n][o]), corner, gauss);
            }
          }
        }
      }
      
      // **********************************
      // PML continuous _ PML continuous _ PML continuous
      // **********************************
    
      PmlContinuous_PmlContinuous_PmlContinuous::PmlContinuous_PmlContinuous_PmlContinuous(const unsigned int corner, PmlContinuous_PmlContinuous *first, PmlContinuous_PmlContinuous *second, PmlContinuous_PmlContinuous *third) : Corner(corner, first, second, third)
      {
      }
    
      PmlContinuous_PmlContinuous_PmlContinuous::~PmlContinuous_PmlContinuous_PmlContinuous()
      {
        delete _uPmlCorner;
        delete _vC[0];
        delete _vC[1];
        delete _vC[2];
        delete _vEdge[0];
        delete _vEdge[1];
        delete _vEdge[2];
      }
    
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* PmlContinuous_PmlContinuous_PmlContinuous::getUPml() const
      {
        return _uPmlCorner;
      }
    
      gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* PmlContinuous_PmlContinuous_PmlContinuous::getV(const unsigned int i) const
      {
        return _vC[i];
      }
    
      PmlContinuous_PmlContinuous *PmlContinuous_PmlContinuous_PmlContinuous::firstEdge() const
      {
        return static_cast< PmlContinuous_PmlContinuous * >(_bnd[0]);
      }
    
      PmlContinuous_PmlContinuous *PmlContinuous_PmlContinuous_PmlContinuous::secondEdge() const
      {
        return static_cast< PmlContinuous_PmlContinuous * >(_bnd[1]);
      }
      
      PmlContinuous_PmlContinuous *PmlContinuous_PmlContinuous_PmlContinuous::thirdEdge() const
      {
        return static_cast< PmlContinuous_PmlContinuous * >(_bnd[2]);
      }
    
      void PmlContinuous_PmlContinuous_PmlContinuous::writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters)
      {
        gmshfem::domain::Domain corner = domains.getCorner(_corner);
        gmshfem::domain::Domain pmlCorner = domains.getPmlCorner(_corner);
        std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > pmlBnd = domains.getPmlEdgeBnd(_corner);
        std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > cornerEdge = domains.getCornerEdge(_corner);
    
        const unsigned int neumannOrder = parameters.getNeumannOrder();
        const unsigned int fieldOrder = parameters.getFieldOrder();
        const std::string gauss = parameters.getGauss();
        const gmshfem::function::ScalarFunction< std::complex< double > > kappa = parameters.getKappa();
    
        _uPmlCorner = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("uPmlCorner_" + orientation(), pmlCorner, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, fieldOrder);
    
        _vC[0] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_first", std::get<0>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
        _vC[1] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_second", std::get<1>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
        _vC[2] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_third", std::get<2>(pmlBnd), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
    
    //    gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<0>(pmlBnd), "vC_" + orientation() + "_first");
    //    gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<1>(pmlBnd), "vC_" + orientation() + "_second");
    //    gmshfem::post::save(gmshfem::function::ScalarFunction< double >(1.), std::get<2>(pmlBnd), "vC_" + orientation() + "_third");
    
        gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *uPml[3];
        uPml[0] = static_cast< PmlContinuous_PmlContinuous * >(_bnd[0])->getUPml();
        uPml[1] = static_cast< PmlContinuous_PmlContinuous * >(_bnd[1])->getUPml();
        uPml[2] = static_cast< PmlContinuous_PmlContinuous * >(_bnd[2])->getUPml();
    
        std::vector< PmlContinuous * > boundaries(6);
        for(unsigned int i = 0; i < 3; ++i) {
          boundaries[2 * i] = static_cast< PmlContinuous_PmlContinuous * >(_bnd[i])->firstBoundary();
          boundaries[2 * i + 1] = static_cast< PmlContinuous_PmlContinuous * >(_bnd[i])->secondBoundary();
        }
        std::sort(boundaries.begin(), boundaries.end());
        gmshfem::function::TensorFunction< std::complex< double > > D[3];
        gmshfem::function::ScalarFunction< std::complex< double > > E[3];
        gmshfem::function::ScalarFunction< std::complex< double > > kappaMod = boundaries[0]->pmlCoefficients(D[0], E[0], corner, kappa);
        boundaries[2]->pmlCoefficients(D[1], E[1], corner, kappa);
        boundaries[4]->pmlCoefficients(D[2], E[2], corner, kappa);
    
        formulation.integral(- D[0] * D[1] * D[2] * grad(dof(*_uPmlCorner)), grad(tf(*_uPmlCorner)), pmlCorner, gauss);
        formulation.integral(kappaMod * kappaMod * E[0] * E[1] * E[2] * dof(*_uPmlCorner), tf(*_uPmlCorner), pmlCorner, gauss);
    
        formulation.integral(dof(*_vC[0]), tf(*_uPmlCorner), std::get<0>(pmlBnd), gauss);
        formulation.integral(- dof(*_vC[0]), tf(*uPml[0]), std::get<0>(pmlBnd), gauss);
    
        formulation.integral(dof(*_uPmlCorner), tf(*_vC[0]), std::get<0>(pmlBnd), gauss);
        formulation.integral(- dof(*uPml[0]), tf(*_vC[0]), std::get<0>(pmlBnd), gauss);
    
        formulation.integral(dof(*_vC[1]), tf(*_uPmlCorner), std::get<1>(pmlBnd), gauss);
        formulation.integral(- dof(*_vC[1]), tf(*uPml[1]), std::get<1>(pmlBnd), gauss);
    
        formulation.integral(dof(*_uPmlCorner), tf(*_vC[1]), std::get<1>(pmlBnd), gauss);
        formulation.integral(- dof(*uPml[1]), tf(*_vC[1]), std::get<1>(pmlBnd), gauss);
        
        formulation.integral(dof(*_vC[2]), tf(*_uPmlCorner), std::get<2>(pmlBnd), gauss);
        formulation.integral(- dof(*_vC[2]), tf(*uPml[2]), std::get<2>(pmlBnd), gauss);
    
        formulation.integral(dof(*_uPmlCorner), tf(*_vC[2]), std::get<2>(pmlBnd), gauss);
        formulation.integral(- dof(*uPml[2]), tf(*_vC[2]), std::get<2>(pmlBnd), gauss);
        
        // corner equation
        _vEdge[0] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vCornerEdge_" + orientation() + "_first", std::get<0>(cornerEdge), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
    
        formulation.integral(dof(*_vC[0]), tf(*_vEdge[0]), std::get<0>(cornerEdge), gauss);
        formulation.integral(dof(*_vEdge[0]), tf(*_vC[0]), std::get<0>(cornerEdge), gauss);
    
        formulation.integral(-dof(*_vC[2]), tf(*_vEdge[0]), std::get<0>(cornerEdge), gauss);
        formulation.integral(-dof(*_vEdge[0]), tf(*_vC[2]), std::get<0>(cornerEdge), gauss);
    
        formulation.integral(dof(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[0])->getV(0)), tf(*_vEdge[0]), std::get<0>(cornerEdge), gauss);
        formulation.integral(dof(*_vEdge[0]), tf(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[0])->getV(0)), std::get<0>(cornerEdge), gauss);
    
        formulation.integral(-dof(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[2])->getV(0)), tf(*_vEdge[0]), std::get<0>(cornerEdge), gauss);
        formulation.integral(-dof(*_vEdge[0]), tf(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[2])->getV(0)), std::get<0>(cornerEdge), gauss);
    
        _vEdge[1] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vCornerEdge_" + orientation() + "_second", std::get<1>(cornerEdge), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
    
        formulation.integral(dof(*_vC[1]), tf(*_vEdge[1]), std::get<1>(cornerEdge), gauss);
        formulation.integral(dof(*_vEdge[1]), tf(*_vC[1]), std::get<1>(cornerEdge), gauss);
    
        formulation.integral(-dof(*_vC[2]), tf(*_vEdge[1]), std::get<1>(cornerEdge), gauss);
        formulation.integral(-dof(*_vEdge[1]), tf(*_vC[2]), std::get<1>(cornerEdge), gauss);
    
        formulation.integral(dof(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[1])->getV(0)), tf(*_vEdge[1]), std::get<1>(cornerEdge), gauss);
        formulation.integral(dof(*_vEdge[1]), tf(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[1])->getV(0)), std::get<1>(cornerEdge), gauss);
    
        formulation.integral(-dof(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[2])->getV(1)), tf(*_vEdge[1]), std::get<1>(cornerEdge), gauss);
        formulation.integral(-dof(*_vEdge[1]), tf(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[2])->getV(1)), std::get<1>(cornerEdge), gauss);
    
        _vEdge[2] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vCornerEdge_" + orientation() + "_third", std::get<2>(cornerEdge), gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
    
        formulation.integral(dof(*_vC[0]), tf(*_vEdge[2]), std::get<2>(cornerEdge), gauss);
        formulation.integral(dof(*_vEdge[2]), tf(*_vC[0]), std::get<2>(cornerEdge), gauss);
    
        formulation.integral(-dof(*_vC[1]), tf(*_vEdge[2]), std::get<2>(cornerEdge), gauss);
        formulation.integral(-dof(*_vEdge[2]), tf(*_vC[1]), std::get<2>(cornerEdge), gauss);
    
        formulation.integral(dof(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[0])->getV(1)), tf(*_vEdge[2]), std::get<2>(cornerEdge), gauss);
        formulation.integral(dof(*_vEdge[2]), tf(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[0])->getV(1)), std::get<2>(cornerEdge), gauss);
    
        formulation.integral(-dof(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[1])->getV(1)), tf(*_vEdge[2]), std::get<2>(cornerEdge), gauss);
        formulation.integral(-dof(*_vEdge[2]), tf(*static_cast< PmlContinuous_PmlContinuous * >(_bnd[1])->getV(1)), std::get<2>(cornerEdge), gauss);
        
    
        _vCorner = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vCorner_" + orientation(), corner, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
    
        formulation.integral(dof(*_vEdge[0]), tf(*_vCorner), corner, gauss);
        formulation.integral(dof(*_vCorner), tf(*_vEdge[0]), corner, gauss);
        
        formulation.integral(dof(*_vEdge[1]), tf(*_vCorner), corner, gauss);
        formulation.integral(dof(*_vCorner), tf(*_vEdge[1]), corner, gauss);
        
        formulation.integral(dof(*_vEdge[2]), tf(*_vCorner), corner, gauss);
        formulation.integral(dof(*_vCorner), tf(*_vEdge[2]), corner, gauss);
      }
    
    
    }