main.cpp

#include "AnalyticalFunction.h"

#include <gmshddm/Formulation.h>
#include <gmshddm/GmshDdm.h>
#include <pybind11/complex.h>
#include <pybind11/embed.h>
#include <pybind11/stl.h>
using namespace gmshfem;
using namespace gmshfem::field;
using namespace gmshfem::function;
using namespace gmshfem::analytics;
using namespace gmshfem::post;
using namespace gmshfem::equation;
//***************************************
// Create surface mesh
//***************************************
class BemppMesh
{
 private:
  int _tagSurface;
  int _dim;
  std::unordered_map< unsigned int, unsigned int >
    _nodeConnectBemppMeshToGmshMesh;
  std::vector< unsigned int > _elementsNode0;
  std::vector< unsigned int > _elementsNode1;
  std::vector< unsigned int > _elementsNode2;
  std::vector< double > _nodesCoordX;
  std::vector< double > _nodesCoordY;
  std::vector< double > _nodesCoordZ;
  unsigned int _renameNode(unsigned int nodeNum,
                           std::unordered_map< unsigned int, unsigned int >
                             &nodeConnectGmshMeshToBemppMesh)
  {
    auto it = nodeConnectGmshMeshToBemppMesh.find(nodeNum);
    if(it != nodeConnectGmshMeshToBemppMesh.end()) {
      return it->second;
    }
    else {
      unsigned int size = nodeConnectGmshMeshToBemppMesh.size();
      nodeConnectGmshMeshToBemppMesh.insert(
        std::pair< unsigned int, unsigned int >(nodeNum, size));
      _nodeConnectBemppMeshToGmshMesh.insert(
        std::pair< unsigned int, unsigned int >(size, nodeNum));
      return size;
    }
  }
  void _createSurfaceMesh()
  {
    gmshfem::common::Timer time;
    time.tick();
    std::vector< int > vEntities;
    gmsh::model::getEntitiesForPhysicalGroup(_dim, _tagSurface, vEntities);
    std::unordered_map< unsigned int, unsigned int >
      nodeConnectGmshMeshToBemppMesh;
    for(unsigned int ent = 0; ent < vEntities.size(); ent++) {
      std::vector< int > elementTypes;
      std::vector< std::vector< std::size_t > > elementTags;
      std::vector< std::vector< std::size_t > > nodeTags;
      gmsh::model::mesh::getElements(elementTypes, elementTags, nodeTags, _dim,
                                     vEntities[ent]);
      std::vector< std::size_t > nodeElem = nodeTags[0];

      int numElem = elementTags[0].size();
      _elementsNode0.reserve(_elementsNode0.size() + numElem);
      _elementsNode1.reserve(_elementsNode1.size() + numElem);
      _elementsNode2.reserve(_elementsNode2.size() + numElem);
      unsigned int numNode0;
      unsigned int numNode1;
      unsigned int numNode2;
      for(int i = 0; i < numElem; i++) {
        numNode0 = _renameNode(nodeElem[3 * i], nodeConnectGmshMeshToBemppMesh);
        numNode1 =
          _renameNode(nodeElem[3 * i + 1], nodeConnectGmshMeshToBemppMesh);
        numNode2 =
          _renameNode(nodeElem[3 * i + 2], nodeConnectGmshMeshToBemppMesh);
        _elementsNode0.push_back(numNode0);
        _elementsNode1.push_back(numNode1);
        _elementsNode2.push_back(numNode2);
      }
    }
    _nodesCoordX.reserve(nodeConnectGmshMeshToBemppMesh.size());
    _nodesCoordY.reserve(nodeConnectGmshMeshToBemppMesh.size());
    _nodesCoordZ.reserve(nodeConnectGmshMeshToBemppMesh.size());
    for(unsigned int i = 0; i < nodeConnectGmshMeshToBemppMesh.size(); i++) {
      std::vector< double > coord;
      std::vector< double > parametricCoord;
      int doo;
      int dii;
      gmsh::model::mesh::getNode(_nodeConnectBemppMeshToGmshMesh[i], coord,
                                 parametricCoord, doo, dii);
      _nodesCoordX.push_back(coord[0]);
      _nodesCoordY.push_back(coord[1]);
      _nodesCoordZ.push_back(coord[2]);
    }
    time.tock();
    msg::info << "bem mesh created in " << time << msg::endl;
  }

 public:
  BemppMesh(const int &tagSurface) :
    _tagSurface(tagSurface), _dim(2)
  {
    _createSurfaceMesh();
  }
  ~BemppMesh() {}

  const std::vector< unsigned int > &getElementsNode0() const
  {
    return _elementsNode0;
  }
  const std::vector< unsigned int > &getElementsNode1() const
  {
    return _elementsNode1;
  }
  const std::vector< unsigned int > &getElementsNode2() const
  {
    return _elementsNode2;
  }
  const std::unordered_map< unsigned int, unsigned int > &
  getnodeConnectBemppMeshToGmshMesh() const
  {
    return _nodeConnectBemppMeshToGmshMesh;
  }
  const std::vector< double > &getnodesCoordX() const
  {
    return _nodesCoordX;
  }
  const std::vector< double > &getnodesCoordY() const
  {
    return _nodesCoordY;
  }
  const std::vector< double > &getnodesCoordZ() const
  {
    return _nodesCoordZ;
  }
};
//***************************************
// BemFormulation
//***************************************
#pragma GCC visibility push(hidden)
class BemFormulation
  : public gmshfem::problem::Formulation< std::complex< double > >
{
 private:
  bool _firstIt;
  pybind11::module _bem;
  gmshfem::field::FieldInterface< std::complex< double > > *_E;
  gmshfem::field::FieldInterface< std::complex< double > > *_H;
  BemppMesh _surfaceMesh;
  std::vector< int > _typeKeys;
  std::vector< unsigned long long > _entityKeys;
  std::complex< double > _Zi;
  std::complex< double > _Ze;
  double _ki;
  double _ke;
  Field< std::complex< double >, Form::Form2 > *_interfaceFieldG;
  std::vector< double > _ChangeOfBasisVectorBemToGmsh;
  double *_bembemL2errorM;
  double *_bembemL2errorJ;
  double _opT;

 public:
  BemFormulation(const std::string &name,
                 gmshfem::field::FieldInterface< std::complex< double > > *E,
                 gmshfem::field::FieldInterface< std::complex< double > > *H,
                 const BemppMesh &surfaceMesh, const std::complex< double > &Zi,
                 const std::complex< double > &Ze, const double &ki,
                 const double &ke, Field< std::complex< double >, Form::Form2 > *g,
                 double *bembemL2errorM, double *bembemL2errorJ,
                 const double &opT) :
    gmshfem::problem::Formulation< std::complex< double > >(name),
    _firstIt(true), _E(E), _H(H), _surfaceMesh(surfaceMesh), _Zi(Zi),
    _Ze(Ze), _ki(ki), _ke(ke), _interfaceFieldG(g),
    _ChangeOfBasisVectorBemToGmsh(), _bembemL2errorM(bembemL2errorM),
    _bembemL2errorJ(bembemL2errorJ), _opT(opT)
  {
    // load our bem.py module
    _bem = pybind11::module::import("bem");
  }

  gmshfem::common::Timer
  pre(const gmshfem::problem::DofsSort::Algorithm algo =
        gmshfem::common::Options::instance()->dofsSortAlgorithm)
  {
    msg::info << "Bem pre-process" << _name << "..." << msg::endl;
    gmshfem::common::Timer time;
    time.tick();

    if(_firstIt) {
      _bem.attr("pre")(
        _ke, _surfaceMesh.getElementsNode0(), _surfaceMesh.getElementsNode1(),
        _surfaceMesh.getElementsNode2(), _surfaceMesh.getnodesCoordX(),
        _surfaceMesh.getnodesCoordY(), _surfaceMesh.getnodesCoordZ(), _opT,
        _Ze, _surfaceMesh.getnodeConnectBemppMeshToGmshMesh());

      pybind11::object edgeNumFollowingBemDofOrderPython =
        _bem.attr("getEdgeNumInBemDofOrder")();
      _entityKeys = pybind11::cast< std::vector< unsigned long long > >(
        edgeNumFollowingBemDofOrderPython);
      _typeKeys.resize(_entityKeys.size(), 1);
      pybind11::object ChangeOfBasisVectorPython =
        _bem.attr("getChangeOfBasisVector")();
      _ChangeOfBasisVectorBemToGmsh =
        pybind11::cast< std::vector< double > >(ChangeOfBasisVectorPython);
      _bem.attr("clearGlobalVariable")();

      for(unsigned int i = 0; i < _entityKeys.size(); ++i) {
        dofs::UnknownDof *dof = new dofs::UnknownDof(
          _typeKeys[i] + GMSHFEM_DOF_FIELD_OFFSET * _E->tag(),
          _entityKeys[i]);
        dof->numDof(i + 1);
        _E->setDof(dof);
      }
      for(unsigned int i = 0; i < _entityKeys.size(); ++i) {
        dofs::UnknownDof *dof = new dofs::UnknownDof(
          _typeKeys[i] + GMSHFEM_DOF_FIELD_OFFSET * _H->tag(),
          _entityKeys[i]);
        dof->numDof(i + 1);
        _H->setDof(dof);
      }
      _firstIt = false;
      time.tock();
      msg::info << _entityKeys.size() << " unknown dofs created in " << time
                << "s:" << msg::endl;
    }
    else {
      time.tock();
      msg::info << "Bem pre-process" << _name << "...Done" << msg::endl;
    }

    return time;
  }

  gmshfem::common::Timer solve(const bool reusePreconditioner = false)
  {
    msg::info << "bem solve" << msg::endl;
    gmshfem::common::Timer time;
    time.tick();
    unsigned int nDof = _entityKeys.size();
    std::vector< std::complex< double > > g_gmshfem(nDof);
    std::vector< std::complex< double > > g_bem;
    g_bem.resize(nDof);
    _interfaceFieldG->getValues(_typeKeys, _entityKeys, g_gmshfem, 0, nDof);
    // low order basis Hdiv in gmshfem (ei low order hcurl basis): n^ei
    // low order basis Hdiv in bempp (ei low order hcurl basis): ei^n
    for(unsigned int i = 0; i < nDof; i++) {
      g_bem[i] = -(g_gmshfem[i]) * _ChangeOfBasisVectorBemToGmsh[i];
    }
    bool physicalSource;
    bool artificialSource;
    getAttribute("ddm::physicalCommutator", physicalSource);
    getAttribute("ddm::artificialCommutator", artificialSource);
    pybind11::object pythonVect =
      _bem.attr("solveIE")(g_bem, physicalSource, artificialSource, _opT);
    std::vector< std::complex< double > > EH_bem =
      pybind11::cast< std::vector< std::complex< double > > >(pythonVect);
    std::vector< std::complex< double > > E_gmshfem(nDof);
    std::vector< std::complex< double > > H_gmshfem(nDof);
    for(unsigned int i = 0; i < nDof; i++) {
      E_gmshfem[i] = EH_bem[i] / _ChangeOfBasisVectorBemToGmsh[i];
    }
    for(unsigned int i = 0; i < nDof; i++) {
      H_gmshfem[i] = EH_bem[i + nDof] / _ChangeOfBasisVectorBemToGmsh[i];
    }
    _E->assignValues(E_gmshfem);
    _H->assignValues(H_gmshfem);
    if(artificialSource && physicalSource) {
      pybind11::object errorPython =
        _bem.attr("getL2errorJM")(g_bem, _Zi, _ki, _ke, _Ze, _opT);
      std::vector< double > error(2);
      error = pybind11::cast< std::vector< double > >(errorPython);
      *_bembemL2errorM = error[0];
      *_bembemL2errorJ = error[1];
    }
    time.tock();
    msg::info << "bem solve in " << time << "s" << msg::endl;
    return time;
  }
};
#pragma GCC visibility pop

//***************************************
// Geometry
//***************************************
void meshSphere(double const &lc, const int &tagSurface, const int &tagVolume,
                const double &r)
{
  gmsh::model::add("sphere");
  gmsh::model::occ::addPoint(0.5, 0, 0, 1);
  gmsh::model::occ::addSphere(0, 0, 0, r, 1);
  gmsh::model::occ::synchronize();
  gmsh::model::addPhysicalGroup(3, {1}, tagVolume);
  gmsh::model::setPhysicalName(3, 1, "volume");
  gmsh::model::addPhysicalGroup(2, {1}, tagSurface);
  gmsh::model::setPhysicalName(2, tagSurface, "surface");
  gmsh::option::setNumber("Mesh.CharacteristicLengthMin", lc);
  gmsh::option::setNumber("Mesh.CharacteristicLengthMax", lc);
  gmsh::model::mesh::generate(3);
}

int main(int argc, char **argv)
{

  gmshddm::common::GmshDdm gmshDdm(argc, argv);
  /****
   *  start python interpreter
   ****/
  pybind11::scoped_interpreter guard{};
  /****
   * input Data
   ****/
  double pi = 3.14159265359;
  double epsrext = 1;
  double murext = 1;
  double epsrint = 4;
  double murint = 1;
  double mu0 = 4 * pi * 1e-7;
  double eps0 = 8.854187817e-12;
  double muext = murext * mu0;
  double epsext = eps0 * epsrext;
  double muint = murint * mu0;
  double epsint = eps0 * epsrint;
  double frequency =
    4 / (2. * pi * std::sqrt(epsint * muint));
  double w = 2 * pi * frequency;
  std::complex< double > im(0., 1.);
  std::complex< double > k1e = im * w * epsext;
  std::complex< double > k2e = im * w * muext;
  double ke = std::sqrt(real(-k1e * k2e));
  std::complex< double > k1i = im * w * epsint;
  std::complex< double > k2i = im * w * muint;
  double ki = std::sqrt(real(-k1i * k2i));
  double Zi = std::sqrt(muint / epsint);
  double Ze = std::sqrt(muext / epsext);
  int tagSurface = 1;
  int tagVolume = 3;
  double r = 1;
  int FEorderAlpha = 2;
  int FEorderBeta = 3;
  double GMRESTol = 4;
  double pointsByWl = 10;
  int pade = 1;
  double angleDivisor = 2.;
  int padeOrder = 2;
  gmshDdm.userDefinedParameter(FEorderAlpha, "FEorderAlpha");
  gmshDdm.userDefinedParameter(FEorderBeta, "FEorderBeta");
  gmshDdm.userDefinedParameter(pointsByWl, "pointsByWl");
  gmshDdm.userDefinedParameter(pade, "pade");
  double lc = 2 * pi / (pointsByWl * std::max(ki, ke));
  std::string gauss =
    "Gauss" + std::to_string(2 * (std::max(FEorderAlpha, FEorderBeta) + 1));
  gmshDdm.userDefinedParameter(gauss, "gauss");
  std::vector< unsigned int > degree = {0, unsigned(FEorderBeta)};
  double mL2Error = 0;
  double jL2Error = 0;
  double opT = Ze;
  double Tp = 1 / opT;
  double Tm = 1 / Ze;
  /****
   * Generate mesh
   ****/
  meshSphere(lc, tagSurface, tagVolume, r);
  /****
   * create surface mesh
   ****/
  BemppMesh surfaceMesh(tagSurface);
  /****
   * ddm part
   ******/
  msg::info << "********************************" << msg::endl;
  msg::info << " Wave properties:" << msg::endl;
  msg::info << " - kint = " << ki << " [1/m]." << msg::endl;
  msg::info << " - kout = " << ke << " [1/m]." << msg::endl;
  msg::info << " Problem properties:" << msg::endl;
  msg::info << " - Ze : " << Ze << msg::endl;
  msg::info << " - Zi : " << Zi << msg::endl;
  msg::info << " - FEorderAlpha = " << FEorderAlpha << "" << msg::endl;
  msg::info << " - FEorderBeta = " << FEorderBeta << "" << msg::endl;
  msg::info << " Mesh properties:" << msg::endl;
  msg::info << " - Approximate number of points by wavelength = " << pointsByWl
            << "" << msg::endl;
  msg::info << " GMRES tolerance:" << msg::endl;
  msg::info << " - " << GMRESTol << msg::endl;
  msg::info << " Use Pade:" << msg::endl;
  msg::info << " - " << pade << msg::endl;
  msg::info << "********************************" << msg::endl;
  // create domains
  gmshddm::domain::Subdomain omega(2);
  gmshddm::domain::Interface sigma(2);
  omega(0) = gmshfem::domain::Domain(3, tagVolume);
  omega(1) = gmshfem::domain::Domain(2, tagSurface);
  sigma(0, 1) = gmshfem::domain::Domain(2, tagSurface);
  sigma(1, 0) = gmshfem::domain::Domain(2, tagSurface);
  // Create fields
  gmshddm::field::SubdomainField< std::complex< double >,
                                  gmshfem::field::Form::Form1 >
    E("E", omega | sigma,
      gmshfem::field::FunctionSpaceTypeForm1::HierarchicalHCurl, FEorderAlpha);
  Field< std::complex< double >, Form::Form1 > Hbem(
    "Hbem", omega(1), FunctionSpaceTypeForm1::HierarchicalHCurl, 0);
  gmshddm::field::InterfaceField< std::complex< double >,
                                  gmshfem::field::Form::Form2 >
    g("g", sigma, gmshfem::field::FunctionSpaceTypeForm2::P_HierarchicalHCurl,
      degree);

  std::vector< gmshfem::problem::Formulation< std::complex< double > > * > vol(2);
  vol[0] = new gmshfem::problem::Formulation< std::complex< double > >("int");
  vol[1] = new BemFormulation("ext", &E(1), &Hbem, surfaceMesh, Zi, Ze, ki, ke,
                              &g(0, 1), &mL2Error, &jL2Error, opT);
  // Define interface formulations
  std::vector< std::unordered_map<
    unsigned int, gmshfem::problem::Formulation< std::complex< double > > * > >
    sur(2);
  sur[0][1] =
    new gmshfem::problem::Formulation< std::complex< double > >("interface_01");
  sur[1][0] =
    new gmshfem::problem::Formulation< std::complex< double > >("interface_10");
  // Define ddm formulation
  gmshddm::problem::Formulation< std::complex< double > > formulation(
    "weak coupling", vol, sur);
  VectorFunction< std::complex< double > > n = normal< std::complex< double > >();
  // Tell to the formulation that g is field that have to be exchanged
  // between subdomains
  formulation.addInterfaceField(g);
  // interior problem formulation
  formulation(0).integral(curl(dof(E(0))), curl(tf(E(0))), omega(0), gauss);
  formulation(0).integral(-ki * ki * dof(E(0)), tf(E(0)), omega(0), gauss);
  // Boundary conditions
  formulation(0).integral(-im * ki * Zi * formulation.artificialSource(g(1, 0)),
                          tf(E(0)), sigma(0, 1), gauss);
  formulation(0, 1).integral(dof(g(0, 1)), tf(g(0, 1)), sigma(0, 1), gauss);
  formulation(0, 1).integral(-Tp * n % (E(0) % n), tf(g(0, 1)), sigma(0, 1),
                             gauss);
  formulation(0, 1).integral(-formulation.artificialSource(g(1, 0)),
                             tf(g(0, 1)), sigma(0, 1), gauss);
  formulation(1, 0).integral(dof(g(1, 0)), tf(g(1, 0)), sigma(1, 0), gauss);
  formulation(1, 0).integral(Tp * E(1), tf(g(1, 0)), sigma(1, 0), gauss);
  formulation(1, 0).integral(-formulation.artificialSource(g(0, 1)), tf(g(1, 0)), sigma(1, 0), gauss);
  std::vector< FieldInterface< std::complex< double > > * > fieldBucket;
  if(pade == 0) {
    formulation(0).integral(-im * ki * Zi * Tm * dof(E(0)), tf(E(0)),
                            sigma(0, 1), gauss);
    formulation(1, 0).integral(Tm * E(1), tf(g(1, 0)), sigma(1, 0), gauss);
    formulation(0, 1).integral(-Tm * n % (E(0) % n), tf(g(0, 1)), sigma(0, 1),
                               gauss);
  }
  else {
    // all what we need for pade approximations
    const double angle = pi / angleDivisor;
    const double M = 2. * padeOrder + 1.;
    std::complex< double > exp1 =
      std::complex< double >(std::cos(-angle / 2), std::sin(-angle / 2));
    std::complex< double > expC =
      std::complex< double >(std::cos(angle / 2.), std::sin(angle / 2.));
    std::complex< double > exp3 =
      std::complex< double >(std::cos(-angle), std::sin(-angle));
    std::complex< double > C0 = expC;
    std::vector< std::complex< double > > Aj(padeOrder, 0.);
    std::vector< std::complex< double > > Bj(padeOrder, 0.);
    std::complex< double > kepsilonsquare = std::complex< double >(
      ke, 0.4 * pow(ke, 1. / 3.) * pow((1 / r), 2. / 3.));
    kepsilonsquare *= kepsilonsquare;
    for(int i = 0; i < padeOrder; i++) {
      std::complex< double > aj = std::sin((i + 1) * pi / M);
      aj = aj * aj * 2. / M;
      std::complex< double > bj = std::cos((i + 1) * pi / M);
      bj *= bj;
      std::complex< double > den = 1. + bj * (exp3 - 1.);
      Aj[i] = exp1 * aj / (den * den);
      Bj[i] = exp3 * bj / den;
      C0 += expC * aj * (exp3 - 1.) / den;
    }
    // Declare Auxiliary  fields
    Field< std::complex< double >, Form::Form1 > *E_auxi0 =
      new Field< std::complex< double >, Form::Form1 >(
        "E_auxi0", sigma(0, 1),
        gmshfem::field::FunctionSpaceTypeForm1::HierarchicalHCurl,
        FEorderBeta);
    Field< std::complex< double >, Form::Form1 > *E_auxi1 =
      new Field< std::complex< double >, Form::Form1 >(
        "E_auxi1", sigma(1, 0),
        gmshfem::field::FunctionSpaceTypeForm1::HierarchicalHCurl,
        0);
    fieldBucket.push_back(E_auxi0);
    fieldBucket.push_back(E_auxi1);
    std::vector< Field< std::complex< double >, Form::Form0 > * > zeta01;
    std::vector< Field< std::complex< double >, Form::Form1 > * > phi01;
    std::vector< Field< std::complex< double >, Form::Form1 > * > phi10;
    std::vector< Field< std::complex< double >, Form::Form0 > * > zeta10;
    for(int i = 0; i < padeOrder; ++i) {
      phi01.push_back(new Field< std::complex< double >, Form::Form1 >(
        "phi01_" + std::to_string(i), sigma(0, 1),
        gmshfem::field::FunctionSpaceTypeForm1::HierarchicalHCurl,
        0));
      zeta01.push_back(new Field< std::complex< double >, Form::Form0 >(
        "zeta01_" + std::to_string(i), sigma(0, 1),
        gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, 1));
      phi10.push_back(new Field< std::complex< double >, Form::Form1 >(
        "phi10_" + std::to_string(i), sigma(1, 0),
        gmshfem::field::FunctionSpaceTypeForm1::HierarchicalHCurl,
        0));
      zeta10.push_back(new Field< std::complex< double >, Form::Form0 >(
        "zeta10_" + std::to_string(i), sigma(1, 0),
        gmshfem::field::FunctionSpaceTypeForm0::HierarchicalH1, 1));
      fieldBucket.push_back(zeta10.back());
      fieldBucket.push_back(phi10.back());
      fieldBucket.push_back(zeta01.back());
      fieldBucket.push_back(phi01.back());
    }

    // Boundary conditions
    formulation(0).integral(im * ki * Zi / Ze * dof(*E_auxi0), tf(E(0)),
                            sigma(0, 1), gauss);
    // Auxiliary equations E_auxi0
    formulation(0).integral(C0 * dof(*E_auxi0), tf(*E_auxi0), sigma(0, 1),
                            gauss);
    formulation(0).integral(dof(E(0)), tf(*E_auxi0), sigma(0, 1), gauss);
    formulation(0).integral(-1. / kepsilonsquare * curl(dof(E(0))),
                            curl(tf(*E_auxi0)), sigma(0, 1), gauss);

    for(int i = 0; i < padeOrder; ++i) {
      formulation(0).integral(Aj[i] * grad(dof(*zeta01[i])), tf(*E_auxi0),
                              sigma(0, 1), gauss);
      formulation(0).integral(-Aj[i] / kepsilonsquare * curl(dof(*phi01[i])),
                              curl(tf(*E_auxi0)), sigma(0, 1), gauss);
      // Auxiliary equations phi01
      formulation(0).integral(dof(*phi01[i]), tf(*phi01[i]), sigma(0, 1),
                              gauss);
      formulation(0).integral(Bj[i] * grad(dof(*zeta01[i])), tf(*phi01[i]),
                              sigma(0, 1), gauss);
      formulation(0).integral(-Bj[i] / kepsilonsquare * curl(dof(*phi01[i])),
                              curl(tf(*phi01[i])), sigma(0, 1), gauss);

      formulation(0).integral(-dof(*E_auxi0), tf(*phi01[i]), sigma(0, 1),
                              gauss);
      // Auxiliary equations zeta01
      formulation(0).integral(dof(*zeta01[i]), tf(*zeta01[i]), sigma(0, 1),
                              gauss);
      formulation(0).integral(1. / kepsilonsquare * dof(*phi01[i]),
                              grad(tf(*zeta01[i])), sigma(0, 1), gauss);
    }
    formulation(0, 1).integral(1. / Ze * (*E_auxi0), tf(g(0, 1)), sigma(0, 1),
                               gauss);
    formulation(1, 0).integral(-1. / Ze * dof(*E_auxi1), tf(g(1, 0)),
                               sigma(1, 0), gauss);
    // Auxiliary equations E_auxi1
    formulation(1, 0).integral(C0 * dof(*E_auxi1), tf(*E_auxi1), sigma(1, 0),
                               gauss);
    formulation(1, 0).integral(-(E(1) % n) % n, tf(*E_auxi1), sigma(1, 0),
                               gauss);
    formulation(1, 0).integral(-1. / kepsilonsquare * curl(E(1)),
                               curl(tf(*E_auxi1)), sigma(1, 0), gauss);
    for(int i = 0; i < padeOrder; ++i) {

      formulation(1, 0).integral(Aj[i] * grad(dof(*zeta10[i])), tf(*E_auxi1),
                                 sigma(1, 0), gauss);

      formulation(1, 0).integral(-Aj[i] / kepsilonsquare * curl(dof(*phi10[i])),
                                 curl(tf(*E_auxi1)), sigma(1, 0), gauss);
      // Auxiliary equations phi10
      formulation(1, 0).integral(dof(*phi10[i]), tf(*phi10[i]), sigma(1, 0),
                                 gauss);
      formulation(1, 0).integral(Bj[i] * grad(dof(*zeta10[i])), tf(*phi10[i]),
                                 sigma(1, 0), gauss);

      formulation(1, 0).integral(-Bj[i] / kepsilonsquare * curl(dof(*phi10[i])),
                                 curl(tf(*phi10[i])), sigma(1, 0), gauss);
      formulation(1, 0).integral(-dof(*E_auxi1), tf(*phi10[i]), sigma(1, 0),
                                 gauss);
      // Auxiliary equations zeta10
      formulation(1, 0).integral(dof(*zeta10[i]), tf(*zeta10[i]), sigma(1, 0),
                                 gauss);
      formulation(1, 0).integral(1. / kepsilonsquare * dof(*phi10[i]),
                                 grad(tf(*zeta10[i])), sigma(1, 0), gauss);
    }
  }

  // Solve the DDM formulation
  formulation.pre();
  formulation.solve("gmres", pow(10, -GMRESTol), 100, true);

  for(unsigned int i = 0; i < fieldBucket.size(); ++i) {
    delete fieldBucket[i];
  }
  msg::info << "***************************************************************************************" <<msg::endl;
  msg::info << "Differences in L2-norm between the strong BEM-BEM coupling and the weak FEM-BEM coupling"<<msg::endl;
  msg::info << "***************************************************************************************" <<msg::endl;
  msg::info <<"e_bembem_M: "  << mL2Error << msg::endl;
  msg::info << "e_bembem_J: " << jL2Error << msg::endl;

  msg::info << "***************************************************************************************" <<msg::endl;
  msg::info << "Errors in L2-norm between the analytial solution and the weak FEM-BEM coupling"<<msg::endl;
  msg::info << "***************************************************************************************" <<msg::endl;
  AnalyticalFunction< maxwell3D::MMagneticSurfaceCurrent< std::complex< double > > >
    Mexact(ki, ke, k2i, k2e, r);
  msg::info << "e_analy_M: "
            << real(sqrt(
                 integrate(pow(abs(xComp(Mexact) - xComp(E(1) % n)), 2) +
                             pow(abs(yComp(Mexact) - yComp(E(1) % n)), 2) +
                             pow(abs(zComp(Mexact) - zComp(E(1) % n)), 2),
                           omega(1), gauss) /
                 integrate(pow(abs(xComp(Mexact)), 2) +
                             pow(abs(yComp(Mexact)), 2) +
                             pow(abs(zComp(Mexact)), 2),
                           omega(1), gauss)))
            << msg::endl;

  AnalyticalFunction<
    maxwell3D::JElectricSurfaceCurrent< std::complex< double > > >
    Jexact(ki, ke, k2i, k2e, r);
  msg::info << "e_analy_J: "<<real(
    sqrt(integrate(pow(abs(xComp(Jexact) - xComp(Hbem % n)), 2) +
                     pow(abs(yComp(Jexact) - yComp(Hbem % n)), 2) +
                     pow(abs(zComp(Jexact) - zComp(Hbem % n)), 2),
                   omega(1), gauss) /
         integrate(pow(abs(xComp(Jexact)), 2) +
                     pow(abs(yComp(Jexact)), 2) +
                     pow(abs(zComp(Jexact)), 2),
                   omega(1), gauss)))<< msg::endl;;

  return 0;
}