Subproblem3D.cpp

#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) {
        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) {
        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) {
        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();
  
    if(N0 != N1) {
      throw gmshfem::common::Exception("Edge with HABC having different number of auxilary fields");
    }
    
    _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);
    }
    
    std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * > uHABC[2];
    for(unsigned int m = 0; m < N0; ++m) {
      uHABC[0].push_back(static_cast< HABC * >(_bnd[0])->getUHABC(m));
    }
    
    for(unsigned int n = 0; n < N1; ++n) {
      uHABC[1].push_back(static_cast< HABC * >(_bnd[1])->getUHABC(n));
    }
    
    const double pi = 3.141592653589793238462643383279;
    const std::complex< double > expPTheta(std::cos(theta0), std::sin(theta0));
    const std::complex< double > expMTheta(std::cos(-theta0), std::sin(-theta0));
    const std::complex< double > expPTheta2(std::cos(theta0/2.), std::sin(theta0/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.);
    
    for(unsigned int m = 0; m < N0; ++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 * dof(*uHABC[0][m]), tf(*_vE[0][m]), edge, gauss);
      
      for(unsigned int n = 0; n < N1; ++n) {
        formulation.integral(-im * kappa * expPTheta2 * 2./M1 * cPade1[n] * dof(*uHABC[0][m]), tf(*_vE[0][m]), edge, gauss);
        formulation.integral(-im * kappa * expPTheta2 * 2./M1 * cPade1[n] * dof(*_uHABC_E[m][n]), tf(*_vE[0][m]), edge, gauss);
      }
    }
    
    for(unsigned int n = 0; n < N1; ++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 * dof(*uHABC[1][n]), tf(*_vE[1][n]), edge, gauss);
      
      for(unsigned int m = 0; m < N0; ++m) {
        formulation.integral(-im * kappa * expPTheta2 * 2./M0 * cPade0[m] * dof(*uHABC[1][n]), tf(*_vE[1][n]), edge, gauss);
        formulation.integral(-im * kappa * expPTheta2 * 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 * (cPade0[m] + cPade1[n] + expMTheta) * dof(*_uHABC_E[m][n]), tf(*_uHABC_E[m][n]), edge, gauss);
        formulation.integral(- kappa * kappa * expPTheta * (cPade0[m] + 1.) * dof(*uHABC[1][n]), tf(*_uHABC_E[m][n]), edge, gauss);
        formulation.integral(- kappa * kappa * expPTheta * (cPade1[n] + 1.) * dof(*uHABC[0][m]), 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(6);
    for(unsigned int i = 0; i < 3; ++i) {
      boundaries[2 * i] = static_cast< HABC_HABC * >(_bnd[i])->firstBoundary();
      boundaries[2 * i + 1] = static_cast< HABC_HABC * >(_bnd[i])->secondBoundary();
    }
    std::sort(boundaries.begin(), boundaries.end());
    
    const unsigned int N0 = static_cast< HABC * >(boundaries[0])->getN();
    const unsigned int N1 = static_cast< HABC * >(boundaries[2])->getN();
    const unsigned int N2 = static_cast< HABC * >(boundaries[4])->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 < N1; ++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 o = 0; o < N2; ++o) {
      for(unsigned int m = 0; m < N0; ++m) {
        delete _vC[2][o][m];
      }
    }
  }
  
  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(6);
    for(unsigned int i = 0; i < 3; ++i) {
      boundaries[2 * i] = static_cast< HABC_HABC * >(_bnd[i])->firstBoundary();
      boundaries[2 * i + 1] = static_cast< HABC_HABC * >(_bnd[i])->secondBoundary();
    }
    std::sort(boundaries.begin(), boundaries.end());
    
    const unsigned int N0 = static_cast< HABC * >(boundaries[0])->getN();
    const unsigned int N1 = static_cast< HABC * >(boundaries[2])->getN();
    const unsigned int N2 = static_cast< HABC * >(boundaries[4])->getN();
    
    const double theta0 = static_cast< HABC * >(boundaries[0])->getTheta();
  //  const double theta1 = static_cast< HABC * >(_bnd[1])->getTheta();
  
    if(N0 != N1 || N1 != N2 || N0 != N2) {
      throw gmshfem::common::Exception("Corner with HABC having different number of auxilary fields");
    }
    
    _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(N1);
      for(unsigned int n = 0; n < N1; ++n) {
        _vC[0][m][n] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_" + std::to_string(m) + "_" + std::to_string(n) + "_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(N2);
    for(unsigned int o = 0; o < N2; ++o) {
      _vC[2][o].resize(N0);
      for(unsigned int m = 0; m < N0; ++m) {
        _vC[2][o][m] = new gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >("vC_" + orientation() + "_" + std::to_string(o) + "_" + std::to_string(m) + "_third", corner, gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, neumannOrder);
      }
    }
    
    std::vector< std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * > > uHABC[3];
    for(unsigned int m = 0; m < N0; ++m) {
      uHABC[0].push_back(std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >());
      for(unsigned int n = 0; n < N1; ++n) {
        uHABC[0][m].push_back(static_cast< HABC_HABC * >(_bnd[0])->getUHABC(m, n));
      }
    }
    for(unsigned int n = 0; n < N1; ++n) {
      uHABC[1].push_back(std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >());
      for(unsigned int o = 0; o < N2; ++o) {
        uHABC[1][n].push_back(static_cast< HABC_HABC * >(_bnd[1])->getUHABC(n, o));
      }
    }
    for(unsigned int o = 0; o < N2; ++o) {
      uHABC[2].push_back(std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * >());
      for(unsigned int m = 0; m < N0; ++m) {
        uHABC[2][o].push_back(static_cast< HABC_HABC * >(_bnd[2])->getUHABC(o, m));
      }
    }
    
    const double pi = 3.141592653589793238462643383279;
    const std::complex< double > expPTheta(std::cos(theta0), std::sin(theta0));
    const std::complex< double > expMTheta(std::cos(-theta0), std::sin(-theta0));
    const std::complex< double > expPTheta2(std::cos(theta0/2.), std::sin(theta0/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.);
    
    for(unsigned int m = 0; m < N0; ++m) {
      for(unsigned int n = 0; n < N1; ++n) {
        formulation.integral(-dof(*_vC[0][m][n]), tf(*uHABC[0][m][n]), corner, gauss);

        formulation.integral( dof(*_vC[0][m][n]), tf(*_vC[0][m][n]), corner, gauss);
        formulation.integral(-im * kappa * expPTheta2 * dof(*uHABC[0][m][n]), tf(*_vC[0][m][n]), corner, gauss);
        
        for(unsigned int o = 0; o < N2; ++o) {
          formulation.integral(-im * kappa * expPTheta2 * 2./M2 * cPade2[o] * dof(*uHABC[0][m][n]), tf(*_vC[0][m][n]), corner, gauss);
          formulation.integral(-im * kappa * expPTheta2 * 2./M2 * cPade2[o] * dof(*_uHABC_C[m][n][o]), tf(*_vC[0][m][n]), corner, gauss);
        }
      }
    }
    
    for(unsigned int n = 0; n < N1; ++n) {
      for(unsigned int o = 0; o < N2; ++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 * 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 * 2./M0 * cPade0[m] * dof(*uHABC[1][n][o]), tf(*_vC[1][n][o]), corner, gauss);
          formulation.integral(-im * kappa * expPTheta2 * 2./M0 * cPade0[m] * dof(*_uHABC_C[m][n][o]), tf(*_vC[1][n][o]), corner, gauss);
        }
      }
    }
    
    for(unsigned int o = 0; o < N2; ++o) {
      for(unsigned int m = 0; m < N0; ++m) {
        formulation.integral(-dof(*_vC[2][o][m]), tf(*uHABC[2][o][m]), corner, gauss);

        formulation.integral( dof(*_vC[2][o][m]), tf(*_vC[2][o][m]), corner, gauss);
        formulation.integral(-im * kappa * expPTheta2 * dof(*uHABC[2][o][m]), tf(*_vC[2][o][m]), corner, gauss);
        
        for(unsigned int n = 0; n < N1; ++n) {
          formulation.integral(-im * kappa * expPTheta2 * 2./M1 * cPade1[n] * dof(*uHABC[2][o][m]), tf(*_vC[2][o][m]), corner, gauss);
          formulation.integral(-im * kappa * expPTheta2 * 2./M1 * cPade1[n] * dof(*_uHABC_C[m][n][o]), tf(*_vC[2][o][m]), 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( (cPade0[m] + cPade1[n] + cPade2[o] + expMTheta) * dof(*_uHABC_C[m][n][o]), tf(*_uHABC_C[m][n][o]), corner, gauss);
          formulation.integral( (cPade0[m] + 1.) * dof(*uHABC[1][n][o]), tf(*_uHABC_C[m][n][o]), corner, gauss);
          formulation.integral( (cPade1[n] + 1.) * dof(*uHABC[2][o][m]), tf(*_uHABC_C[m][n][o]), corner, gauss);
          formulation.integral( (cPade2[o] + 1.) * dof(*uHABC[0][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);
  }


}