New bicgl solver

demos/helmholtz/biGCL/CMakeLists.txt

+cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
+
+project(demo CXX)
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/../../demos.cmake)
+
+add_executable(demo main.cpp ${EXTRA_INCS})
+target_link_libraries(demo ${EXTRA_LIBS})

demos/helmholtz/biGCL/main.cpp

+#include <gmshddm/GmshDdm.h>
+#include <gmshddm/Formulation.h>
+
+#include <gmshfem/Message.h>
+
+#include <gmsh.h>
+
+using gmshfem::equation::dof;
+using gmshfem::equation::tf;
+using gmshfem::function::operator-;
+using gmshfem::function::operator*;
+
+void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const double xSize, const double ySize, double pmlSize, const double lc, bool pml, const int order, const unsigned int sourcePosX, const unsigned int sourcePosY)
+{
+  //    ***********
+  // j  *         *
+  // ^  *         *
+  // |  *         *
+  // |  ***********
+  // |
+  // -------> i
+  std::vector< std::vector< int > > omega(nDomX);
+  for(unsigned int i = 0; i < nDomX; ++i) {
+    omega[i].resize(nDomY);
+  }
+  
+  // ******************
+  // borders
+  // ******************
+  
+  //  ****(1)****
+  //  *         *
+  // (2)       (0)
+  //  *         *
+  //  ****(3)****
+  std::vector< std::vector< std::vector< int > > > borders(nDomX);
+  for(unsigned int i = 0; i < nDomX; ++i) {
+    borders[i].resize(nDomY);
+    for(unsigned int j = 0; j < nDomY; ++j) {
+      borders[i][j].resize(4);
+    }
+  }
+  
+  // ******************
+  // pmls
+  // ******************
+
+  //   1         0
+  // 0 ****(1)**** 1
+  //   *         *
+  //  (2)       (0)
+  //   *         *
+  // 1 ****(3)**** 0
+  //   0         1
+  
+  std::vector< std::vector< std::vector< int > > > pmls(nDomX);
+  std::vector< std::vector< std::vector< int > > > pmlsInf(nDomX);
+  std::vector< std::vector< std::vector< std::pair< int, int > > > > pmlsBnd(nDomX);
+  for(unsigned int i = 0; i < nDomX; ++i) {
+    pmls[i].resize(nDomY);
+    pmlsInf[i].resize(nDomY);
+    pmlsBnd[i].resize(nDomY);
+    for(unsigned int j = 0; j < nDomY; ++j) {
+      pmls[i][j].resize(4);
+      pmlsInf[i][j].resize(4);
+      pmlsBnd[i][j].resize(4);
+    }
+  }
+  
+  // ******************
+  // pmlsSquare
+  // ******************
+  
+  // (1)*******(0)
+  //  *         *
+  //  *         *
+  //  *         *
+  // (2)*******(3)
+  
+  std::vector< std::vector< std::vector< int > > > pmlsSquare(nDomX);
+  for(unsigned int i = 0; i < nDomX; ++i) {
+    pmlsSquare[i].resize(nDomY);
+    for(unsigned int j = 0; j < nDomY; ++j) {
+      pmlsSquare[i][j].resize(4);
+    }
+  }
+  
+//  std::vector< std::vector< std::vector< int > > > corners(nDomX);
+//  for(unsigned int i = 0; i < nDomX; ++i) {
+//    corners[i].resize(nDomY);
+//    for(unsigned int j = 0; j < nDomY; ++j) {
+//      corners[i][j].resize(4);
+//    }
+//  }
+  
+  std::vector< int > left(nDomY);
+  int bottom = 0;
+  for(unsigned int i = 0; i < nDomX; ++i) {
+    for(unsigned int j = 0; j < nDomY; ++j) {
+      // (p3)*******(p2)
+      //  *          *
+      //  *          *
+      //  *          *
+      // (p0)*******(p1)
+      int p[4];
+      p[0] = gmsh::model::geo::addPoint(i * xSize,       j * ySize,       0., lc);
+      p[1] = gmsh::model::geo::addPoint((i + 1) * xSize, j * ySize,       0., lc);
+      p[2] = gmsh::model::geo::addPoint((i + 1) * xSize, (j + 1) * ySize, 0., lc);
+      p[3] = gmsh::model::geo::addPoint(i * xSize,       (j + 1) * ySize, 0., lc);
+      
+//      corners[i][j][0] = p[2];
+//      corners[i][j][1] = p[3];
+//      corners[i][j][2] = p[0];
+//      corners[i][j][3] = p[1];
+
+      int l[4];
+      borders[i][j][0] = l[0] = gmsh::model::geo::addLine(p[1], p[2]);
+      borders[i][j][1] = l[1] = gmsh::model::geo::addLine(p[2], p[3]);
+      borders[i][j][2] = l[2] = gmsh::model::geo::addLine(p[3], p[0]);
+      borders[i][j][3] = l[3] = gmsh::model::geo::addLine(p[0], p[1]);
+      
+      if(i == 0) {
+        left[j] = l[0];
+      }
+      else {
+        borders[i][j][2] = -left[j];
+        left[j] = l[0];
+      }
+      
+      if(j == 0) {
+        bottom = l[1];
+      }
+      else {
+        borders[i][j][3] = -bottom;
+        bottom = l[1];
+      }
+            
+      int ll = gmsh::model::geo::addCurveLoop({l[0], l[1], l[2], l[3]});
+      
+      int surface = gmsh::model::geo::addPlaneSurface({ll});
+      omega[i][j] = surface;
+      
+      if(i == sourcePosX && j == sourcePosY) {
+        // source
+        int pS = gmsh::model::geo::addPoint(sourcePosX * i + xSize / 2., sourcePosY * j + ySize / 2., 0., lc);
+        
+        gmsh::model::geo::synchronize();
+        gmsh::model::mesh::embed(0, {pS}, 2, surface);
+        
+        gmsh::model::addPhysicalGroup(0, {pS}, 1);
+        gmsh::model::setPhysicalName(0, 1, "source");
+      }
+      
+      if(pml) {
+        // pmls
+        const double extrudeX[4] = {pmlSize, 0., -pmlSize, 0.};
+        const double extrudeY[4] = {0., pmlSize, 0., -pmlSize};
+        
+        for(unsigned int k = 0; k < 4; ++k) {
+          std::vector< std::pair< int, int > > extrudeEntities;
+          gmsh::model::geo::extrude({std::make_pair(1, -borders[i][j][k])}, extrudeX[k], extrudeY[k], 0., extrudeEntities, {static_cast< int >(/*(k % 2 ? 1. : 2.) */ pmlSize / lc)}, std::vector< double >(), true);
+          pmlsInf[i][j][k] = extrudeEntities[0].second;
+          pmls[i][j][k] = extrudeEntities[1].second;
+          if(-borders[i][j][k] == extrudeEntities[2].second) {
+            pmlsBnd[i][j][k].first = extrudeEntities[3].second;
+            pmlsBnd[i][j][k].second = extrudeEntities[4].second;
+          }
+          else {
+            pmlsBnd[i][j][k].first = extrudeEntities[2].second;
+            pmlsBnd[i][j][k].second = extrudeEntities[3].second;
+          }
+        }
+        
+        // pmlsSquare
+        std::vector< std::pair< int, int > > extrudeEntities;
+        
+        gmsh::model::geo::extrude({std::make_pair(1, -pmlsBnd[i][j][0].second)}, 0., pmlSize, 0., extrudeEntities, {static_cast< int >(pmlSize / lc)}, std::vector< double >(), true);
+        pmlsSquare[i][j][0] = extrudeEntities[1].second;
+        extrudeEntities.clear();
+        
+        gmsh::model::geo::extrude({std::make_pair(1, -pmlsBnd[i][j][2].first)}, 0., pmlSize, 0., extrudeEntities, {static_cast< int >(pmlSize / lc)}, std::vector< double >(), true);
+        pmlsSquare[i][j][1] = extrudeEntities[1].second;
+        extrudeEntities.clear();
+        
+        gmsh::model::geo::extrude({std::make_pair(1, -pmlsBnd[i][j][2].second)}, 0., -pmlSize, 0., extrudeEntities, {static_cast< int >(pmlSize / lc)}, std::vector< double >(), true);
+        pmlsSquare[i][j][2] = extrudeEntities[1].second;
+        extrudeEntities.clear();
+        
+        gmsh::model::geo::extrude({std::make_pair(1, -pmlsBnd[i][j][0].first)}, 0., -pmlSize, 0., extrudeEntities, {static_cast< int >(pmlSize / lc)}, std::vector< double >(), true);
+        pmlsSquare[i][j][3] = extrudeEntities[1].second;
+        extrudeEntities.clear();
+      }
+    }
+  }
+  gmsh::model::geo::removeAllDuplicates();
+  
+  gmsh::model::geo::synchronize();
+    
+  // ******************
+  // physicals
+  // ******************
+
+  for(unsigned int i = 0; i < nDomX; ++i) {
+    for(unsigned int j = 0; j < nDomY; ++j) {
+      gmsh::model::addPhysicalGroup(2, {omega[i][j]}, i * nDomY + j + 1);
+      gmsh::model::setPhysicalName(2, i * nDomY + j + 1, "omega_" + std::to_string(i) + "_" + std::to_string(j));
+      
+      gmsh::model::addPhysicalGroup(1, {borders[i][j][0]},  10000 + i * nDomY + j + 1);
+      gmsh::model::setPhysicalName(1, 10000 + i * nDomY + j + 1, "sigmaE_" + std::to_string(i) + "_" + std::to_string(j));
+      
+      gmsh::model::addPhysicalGroup(1, {borders[i][j][1]},  20000 + i * nDomY + j + 1);
+      gmsh::model::setPhysicalName(1, 20000 + i * nDomY + j + 1, "sigmaN_" + std::to_string(i) + "_" + std::to_string(j));
+      
+      gmsh::model::addPhysicalGroup(1, {borders[i][j][2]},  30000 + i * nDomY + j + 1);
+      gmsh::model::setPhysicalName(1, 30000 + i * nDomY + j + 1, "sigmaW_" + std::to_string(i) + "_" + std::to_string(j));
+      
+      gmsh::model::addPhysicalGroup(1, {borders[i][j][3]},  40000 + i * nDomY + j + 1);
+      gmsh::model::setPhysicalName(1, 40000 + i * nDomY + j + 1, "sigmaS_" + std::to_string(i) + "_" + std::to_string(j));
+      
+      
+//      gmsh::model::addPhysicalGroup(0, {corners[i][j][0]},  10000 + i * nDomY + j + 1);
+//      gmsh::model::setPhysicalName(0, 10000 + i * nDomY + j + 1, "crossPointNE_" + std::to_string(i) + "_" + std::to_string(j));
+//
+//      gmsh::model::addPhysicalGroup(0, {corners[i][j][1]},  20000 + i * nDomY + j + 1);
+//      gmsh::model::setPhysicalName(0, 20000 + i * nDomY + j + 1, "crossPointNW_" + std::to_string(i) + "_" + std::to_string(j));
+//
+//      gmsh::model::addPhysicalGroup(0, {corners[i][j][2]},  30000 + i * nDomY + j + 1);
+//      gmsh::model::setPhysicalName(0, 30000 + i * nDomY + j + 1, "crossPointSW_" + std::to_string(i) + "_" + std::to_string(j));
+//
+//      gmsh::model::addPhysicalGroup(0, {corners[i][j][3]},  40000 + i * nDomY + j + 1);
+//      gmsh::model::setPhysicalName(0, 40000 + i * nDomY + j + 1, "crossPointSE_" + std::to_string(i) + "_" + std::to_string(j));
+
+      if(pml) {
+        // PMLs
+        gmsh::model::addPhysicalGroup(2, {pmls[i][j][0]}, 10000 + i * nDomY + j + 1);
+        gmsh::model::setPhysicalName(2, 10000 + i * nDomY + j + 1, "pmlE_" + std::to_string(i) + "_" + std::to_string(j));
+        
+        gmsh::model::addPhysicalGroup(2, {pmls[i][j][1]}, 20000 + i * nDomY + j + 1);
+        gmsh::model::setPhysicalName(2, 20000 + i * nDomY + j + 1, "pmlN_" + std::to_string(i) + "_" + std::to_string(j));
+        
+        gmsh::model::addPhysicalGroup(2, {pmls[i][j][2]}, 30000 + i * nDomY + j + 1);
+        gmsh::model::setPhysicalName(2, 30000 + i * nDomY + j + 1, "pmlW_" + std::to_string(i) + "_" + std::to_string(j));
+        
+        gmsh::model::addPhysicalGroup(2, {pmls[i][j][3]}, 40000 + i * nDomY + j + 1);
+        gmsh::model::setPhysicalName(2, 40000 + i * nDomY + j + 1, "pmlS_" + std::to_string(i) + "_" + std::to_string(j));
+        
+        // PMLs square
+        gmsh::model::addPhysicalGroup(2, {pmlsSquare[i][j][0]}, 100000 + i * nDomY + j + 1);
+        gmsh::model::setPhysicalName(2, 100000 + i * nDomY + j + 1, "pmlNE_" + std::to_string(i) + "_" + std::to_string(j));
+        
+        gmsh::model::addPhysicalGroup(2, {pmlsSquare[i][j][1]}, 200000 + i * nDomY + j + 1);
+        gmsh::model::setPhysicalName(2, 200000 + i * nDomY + j + 1, "pmlNW_" + std::to_string(i) + "_" + std::to_string(j));
+        
+        gmsh::model::addPhysicalGroup(2, {pmlsSquare[i][j][2]}, 300000 + i * nDomY + j + 1);
+        gmsh::model::setPhysicalName(2, 300000 + i * nDomY + j + 1, "pmlSW_" + std::to_string(i) + "_" + std::to_string(j));
+        
+        gmsh::model::addPhysicalGroup(2, {pmlsSquare[i][j][3]}, 400000 + i * nDomY + j + 1);
+        gmsh::model::setPhysicalName(2, 400000 + i * nDomY + j + 1, "pmlSE_" + std::to_string(i) + "_" + std::to_string(j));
+        
+        // Sigmas in PMLs
+        gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][0].first}, 100000 + i * nDomY + j + 1);
+        gmsh::model::setPhysicalName(1, 100000 + i * nDomY + j + 1, "sigmaE_S_" + std::to_string(i) + "_" + std::to_string(j));
+        gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][0].second}, 150000 + i * nDomY + j + 1);
+        gmsh::model::setPhysicalName(1, 150000 + i * nDomY + j + 1, "sigmaE_N_" + std::to_string(i) + "_" + std::to_string(j));
+        
+        gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][1].first}, 200000 + i * nDomY + j + 1);
+        gmsh::model::setPhysicalName(1, 200000 + i * nDomY + j + 1, "sigmaN_E_" + std::to_string(i) + "_" + std::to_string(j));
+        gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][1].second}, 250000 + i * nDomY + j + 1);
+        gmsh::model::setPhysicalName(1, 250000 + i * nDomY + j + 1, "sigmaN_W_" + std::to_string(i) + "_" + std::to_string(j));
+        
+        gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][2].first}, 300000 + i * nDomY + j + 1);
+        gmsh::model::setPhysicalName(1, 300000 + i * nDomY + j + 1, "sigmaW_N_" + std::to_string(i) + "_" + std::to_string(j));
+        gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][2].second}, 350000 + i * nDomY + j + 1);
+        gmsh::model::setPhysicalName(1, 350000 + i * nDomY + j + 1, "sigmaW_S_" + std::to_string(i) + "_" + std::to_string(j));
+        
+        gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][3].first}, 400000 + i * nDomY + j + 1);
+        gmsh::model::setPhysicalName(1, 400000 + i * nDomY + j + 1, "sigmaS_W_" + std::to_string(i) + "_" + std::to_string(j));
+        gmsh::model::addPhysicalGroup(1, {pmlsBnd[i][j][3].second}, 450000 + i * nDomY + j + 1);
+        gmsh::model::setPhysicalName(1, 450000 + i * nDomY + j + 1, "sigmaS_E_" + std::to_string(i) + "_" + std::to_string(j));
+      }
+    }
+  }
+
+  gmsh::model::geo::synchronize();
+  gmsh::model::mesh::generate();
+  gmsh::model::mesh::setOrder(order);
+}
+
+int main(int argc, char **argv)
+{
+  gmshddm::common::GmshDdm gmshDdm(argc, argv);
+  
+  int nDomX = 10, nDomY = 10;
+  gmshDdm.userDefinedParameter(nDomX, "nDomX");
+  gmshDdm.userDefinedParameter(nDomY, "nDomY");
+  double xSize = 0.1, ySize = 0.1;
+  gmshDdm.userDefinedParameter(xSize, "xSize");
+  gmshDdm.userDefinedParameter(ySize, "ySize");
+  double pmlSize = (xSize + ySize) / 10.;
+  gmshDdm.userDefinedParameter(pmlSize, "pmlSize");
+  bool withPml = false;
+  gmshDdm.userDefinedParameter(withPml, "pml");
+  double lc = 0.1;
+  gmshDdm.userDefinedParameter(lc, "lc");
+  int meshOrder = 1;
+  gmshDdm.userDefinedParameter(meshOrder, "meshOrder");
+  int sourcePosX = 0, sourcePosY = 0;
+  gmshDdm.userDefinedParameter(sourcePosX, "sourcePosX");
+  gmshDdm.userDefinedParameter(sourcePosY, "sourcePosY");
+  
+  // Build geometry and mesh using gmsh API
+  checkerboard(nDomX, nDomY, xSize, ySize, pmlSize, lc, withPml, meshOrder, sourcePosX, sourcePosY);
+  
+  const double pi = 3.14159265359;
+  const std::complex< double > im(0., 1.);
+  double k = 2. * pi;
+  gmshDdm.userDefinedParameter(k, "k");
+  std::string interface = "sommerfeld";
+  gmshDdm.userDefinedParameter(interface, "interface");
+  std::string gauss = "Gauss4";
+  gmshDdm.userDefinedParameter(gauss, "gauss");
+  
+  unsigned int nDom = nDomX * nDomY;
+  // Define domain
+  gmshddm::domain::Subdomain omega(nDom);
+  gmshfem::domain::Domain gamma = gmshfem::domain::Domain(0, 1);
+
+  std::vector< gmshddm::domain::Subdomain > sigma(4, nDom);
+  std::vector< gmshddm::domain::Interface > sigmaInterface(4, nDom);
+
+  std::vector< std::pair< gmshddm::domain::Subdomain, gmshddm::domain::Subdomain > > corner(4, std::make_pair(gmshddm::domain::Subdomain(nDom), gmshddm::domain::Subdomain(nDom)) );
+  // Define topology
+  std::vector< std::vector< unsigned int > > topology(nDom);
+  for(unsigned int i = 0; i < static_cast< unsigned int >(nDomX); ++i) {
+    for(unsigned int j = 0; j < static_cast< unsigned int >(nDomY); ++j) {
+      unsigned int index = i * nDomY + j;
+
+      omega(index) = gmshfem::domain::Domain(2, index + 1);
+
+      for(unsigned int f = 0; f < 4; ++f) {
+        sigma[f](index) = gmshfem::domain::Domain(1, (f+1) * 10000 + i * nDomY + j + 1);
+//        corner[f].first(index) = gmshfem::domain::Domain(0, ((f-1)%4 + 1) * 10000 + i * nDomY + j + 1);
+//        corner[f].second(index) = gmshfem::domain::Domain(0, (f + 1) * 10000 + i * nDomY + j + 1);
+      }
+
+      if(i != static_cast< unsigned int >(nDomX) - 1) { // E
+        sigmaInterface[0](index, (i+1) * nDomY + j) = sigma[0](index);
+        topology[index].push_back((i+1) * nDomY + j);
+      }
+
+      if(j != static_cast< unsigned int >(nDomY) - 1) { // N
+        sigmaInterface[1](index, i * nDomY + (j+1)) = sigma[1](index);
+        topology[index].push_back(i * nDomY + (j+1));
+      }
+
+      if(i != 0) { // W
+        sigmaInterface[2](index, (i-1) * nDomY + j) = sigma[2](index);
+        topology[index].push_back((i-1) * nDomY + j);
+      }
+
+      if(j != 0) { // S
+        sigmaInterface[3](index, i * nDomY + (j-1)) = sigma[3](index);
+        topology[index].push_back(i * nDomY + (j-1));
+      }
+    }
+  }
+  
+  gmshddm::problem::Formulation< std::complex< double > > formulation("HelmholtzDDM", topology);
+
+  // Fields definition
+  gmshddm::field::SubdomainField< std::complex< double >, gmshfem::field::Form::Form0 > u("u", omega | sigma[0] | sigma[1] | sigma[2] | sigma[3], gmshfem::field::FunctionSpaceTypeForm0::Lagrange);
+  u(sourcePosX * nDomY + sourcePosY).domain(u(sourcePosX * nDomY + sourcePosY).domain() | gamma);
+
+  gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > lambda("lambda", gamma, gmshfem::field::FunctionSpaceTypeForm0::Lagrange);
+
+  const std::vector< std::string > orientation {"E", "N", "W", "S"};
+  std::vector< gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 > > g;
+  g.reserve(4);
+  for(unsigned int f = 0; f < 4; ++f) {
+    g.push_back(gmshddm::field::InterfaceField< std::complex< double >, gmshfem::field::Form::Form0 >("g_" + orientation[f], sigmaInterface[f], gmshfem::field::FunctionSpaceTypeForm0::Lagrange));
+    formulation.addInterfaceField(g[f]);
+  }
+
+  std::vector< gmshddm::field::SubdomainField< std::complex< double >, gmshfem::field::Form::Form0 > > v;
+  v.reserve(4);
+  for(unsigned int f = 0; f < 4; ++f) {
+    v.push_back(gmshddm::field::SubdomainField< std::complex< double >, gmshfem::field::Form::Form0 >("v_" + orientation[f], sigma[f], gmshfem::field::FunctionSpaceTypeForm0::Lagrange));
+  }
+
+  for(unsigned int i = 0; i < nDomX * nDomY; ++i) {
+    formulation(i).galerkin(-grad(dof(u(i))), grad(tf(u(i))), omega(i), gauss);
+    formulation(i).galerkin(k * k * dof(u(i)), tf(u(i)), omega(i), gauss);
+
+    if(i == sourcePosX * nDomY + sourcePosY) {
+      formulation(i).galerkin(dof(lambda), tf(u(i)), gamma, gauss);
+      formulation(i).galerkin(dof(u(i)), tf(lambda), gamma, gauss);
+      formulation(i).galerkin(formulation.physicalSource(- 1.), tf(lambda), gamma, gauss);
+    }
+    
+    // boundary
+    if(interface == "sommerfeld") {
+      for(unsigned int f = 0; f < 4; ++f) {
+        formulation(i).galerkin(-dof(v[f](i)), tf(u(i)), sigma[f](i), gauss);
+        formulation(i).galerkin(dof(v[f](i)), tf(v[f](i)), sigma[f](i), gauss);
+        formulation(i).galerkin(im * k * dof(u(i)), tf(v[f](i)), sigma[f](i), gauss);
+      }
+    }
+      
+    // coupling
+    for(unsigned int f = 0; f < 4; ++f) {
+      for(unsigned int jj = 0; jj < topology[i].size(); ++jj) {
+        const unsigned int j = topology[i][jj];
+        if(!sigmaInterface[f](i,j).isEmpty()) {
+          formulation(i).galerkin(formulation.artificialSource( - g[(f+2)%4](j,i) ), tf(u(i)), sigmaInterface[f](i,j), gauss);
+        }
+      }
+    }
+
+    // interface
+    for(unsigned int f = 0; f < 4; ++f) {
+      for(unsigned int jj = 0; jj < topology[i].size(); ++jj) {
+        const unsigned int j = topology[i][jj];
+        if(!sigmaInterface[f](i,j).isEmpty()) {
+          formulation(i,j).galerkin(dof(g[f](i,j)), tf(g[f](i,j)), sigmaInterface[f](i,j), gauss);
+          formulation(i,j).galerkin(formulation.artificialSource( g[(f+2)%4](j,i) ), tf(g[f](i,j)), sigmaInterface[f](i,j), gauss);
+          formulation(i,j).galerkin(2. * v[f](i), tf(g[f](i,j)), sigmaInterface[f](i,j), gauss);
+          
+        }
+      }
+    }
+  }
+  
+  double tol = 1e-4;
+  gmshDdm.userDefinedParameter(tol, "tol");
+  std::string solver = "gmres";
+  gmshDdm.userDefinedParameter(solver, "solver");
+
+  formulation.pre();
+  formulation.solve(solver, tol, 1000);
+
+  for(unsigned int i = 0; i < nDom; ++i) {
+    gmshfem::post::save(u(i), omega(i), "u_" + std::to_string(i), "msh");
+  }
+
+  return 0;
+}

src/field/InterfaceField.cpp

@@ -9,40 +9,71 @@ namespace gmshddm
   namespace field
   {
 
+    //
+    // class InterfaceFieldInterface
+    //
+    
+    template< class T_Scalar >
+    InterfaceFieldInterface< T_Scalar >::InterfaceFieldInterface(const std::string &name) :
+      _name(name)
+    {
+    }
+    
+    template< class T_Scalar >
+    InterfaceFieldInterface< T_Scalar >::~InterfaceFieldInterface()
+    {
+    }
+    
+    template< class T_Scalar >
+    std::string InterfaceFieldInterface< T_Scalar >::name() const
+    {
+      return _name;
+    }
+    
+    
+    template class InterfaceFieldInterface< std::complex< double > >;
+    template class InterfaceFieldInterface< double >;
+    template class InterfaceFieldInterface< std::complex< float > >;
+    template class InterfaceFieldInterface< float >;
+    
+    
+    //
+    // class InterfaceField
+    //
     
     template< class T_Scalar, gmshfem::field::Form T_Form >
     InterfaceField< T_Scalar, T_Form >::InterfaceField() :
-      InterfaceFieldInterface< T_Scalar >(), _name(), _fields()
+      InterfaceFieldInterface< T_Scalar >(), _fields()
     {
     }
 
     template< class T_Scalar, gmshfem::field::Form T_Form >
     InterfaceField< T_Scalar, T_Form >::InterfaceField(const std::string &name, const domain::Interface &domains, const gmshfem::field::FunctionSpaceOfForm< T_Form > &type, const unsigned int degree) :
-      InterfaceFieldInterface< T_Scalar >(), _name(name), _fields(domains.numberOfSubdomains())
+      InterfaceFieldInterface< T_Scalar >(name), _fields(domains.numberOfSubdomains())
     {
       for(unsigned int i = 0; i < domains.numberOfSubdomains(); ++i) {
         for(auto it = domains[i].begin(); it != domains[i].end(); ++it) {
-          _fields[i].insert(std::make_pair(it->first, gmshfem::field::Field< T_Scalar, T_Form >(_name + "_" + std::to_string(i) + "_" + std::to_string(it->first), it->second, type, degree)));
+          _fields[i].insert(std::make_pair(it->first, gmshfem::field::Field< T_Scalar, T_Form >(this->_name + "_" + std::to_string(i) + "_" + std::to_string(it->first), it->second, type, degree)));
         }
       }
     }
 
     template< class T_Scalar, gmshfem::field::Form T_Form >
     InterfaceField< T_Scalar, T_Form >::InterfaceField(const InterfaceField< T_Scalar, T_Form > &other) :
-      InterfaceFieldInterface< T_Scalar >(), _name(other._name), _fields(other._fields)
+      InterfaceFieldInterface< T_Scalar >(other._name), _fields(other._fields)
     {
     }
 
     template< class T_Scalar, gmshfem::field::Form T_Form >
     InterfaceField< T_Scalar, T_Form >::InterfaceField(InterfaceField< T_Scalar, T_Form > &&other) :
-      InterfaceFieldInterface< T_Scalar >(), _name(std::move(other._name)), _fields(std::move(other._fields))
+      InterfaceFieldInterface< T_Scalar >(other._name), _fields(std::move(other._fields))
     {
     }
 
     template< class T_Scalar, gmshfem::field::Form T_Form >
     InterfaceField< T_Scalar, T_Form > &InterfaceField< T_Scalar, T_Form >::operator=(const InterfaceField< T_Scalar, T_Form > &other)
     {
-      _name = other._name;
+      this->_name = other._name;
       _fields = other._fields;
       return *this;
     }
@@ -72,11 +103,11 @@ namespace gmshddm
     const gmshfem::field::Field< T_Scalar, T_Form > &InterfaceField< T_Scalar, T_Form >::operator()(const unsigned int i, const unsigned int j) const
     {
       if(i >= _fields.size()) {
-        throw gmshfem::common::Exception("Try to acces to field(" + std::to_string(i) + ", " + std::to_string(j) + ") of field '" + _name + "'");
+        throw gmshfem::common::Exception("Try to acces to field(" + std::to_string(i) + ", " + std::to_string(j) + ") of field '" + this->_name + "'");
       }
       auto it = _fields[i].find(j);
       if(it == _fields[i].end()) {
-        throw gmshfem::common::Exception("Try to acces to field(" + std::to_string(i) + ", " + std::to_string(j) + ") of field '" + _name + "'");
+        throw gmshfem::common::Exception("Try to acces to field(" + std::to_string(i) + ", " + std::to_string(j) + ") of field '" + this->_name + "'");
       }
 
       return it->second;
@@ -86,11 +117,11 @@ namespace gmshddm
     gmshfem::field::Field< T_Scalar, T_Form > &InterfaceField< T_Scalar, T_Form >::operator()(const unsigned int i, const unsigned int j)
     {
       if(i >= _fields.size()) {
-        throw gmshfem::common::Exception("Try to acces to field(" + std::to_string(i) + ", " + std::to_string(j) + ") of field '" + _name + "'");
+        throw gmshfem::common::Exception("Try to acces to field(" + std::to_string(i) + ", " + std::to_string(j) + ") of field '" + this->_name + "'");
       }
       auto it = _fields[i].find(j);
       if(it == _fields[i].end()) {
-        throw gmshfem::common::Exception("Try to acces to field(" + std::to_string(i) + ", " + std::to_string(j) + ") of field '" + _name + "'");
+        throw gmshfem::common::Exception("Try to acces to field(" + std::to_string(i) + ", " + std::to_string(j) + ") of field '" + this->_name + "'");
       }
 
       return it->second;

src/field/InterfaceField.h

@@ -19,20 +19,23 @@ namespace gmshddm
     template< class T_Scalar >
     class InterfaceFieldInterface
     {
+     protected:
+      std::string _name;
+      
      public:
-      InterfaceFieldInterface() {}
-      virtual ~InterfaceFieldInterface() {}
+      InterfaceFieldInterface(const std::string &name = "");
+      virtual ~InterfaceFieldInterface();
 
       virtual bool isDefined(const unsigned int i, const unsigned int j) const = 0;
       virtual const gmshfem::field::FieldInterface< T_Scalar > &operator()(const unsigned int i, const unsigned int j) const = 0;
       virtual gmshfem::field::FieldInterface< T_Scalar > &operator()(const unsigned int i, const unsigned int j) = 0;
+      
+      std::string name() const;
     };
 
     template< class T_Scalar, gmshfem::field::Form T_Form >
     class InterfaceField : public InterfaceFieldInterface< T_Scalar >
     {
-     private:
-      std::string _name;
       std::vector< std::unordered_map< unsigned int, gmshfem::field::Field< T_Scalar, T_Form > > > _fields;
 
      public:

src/problem/Formulation.cpp

@@ -84,6 +84,7 @@ namespace gmshddm
     void Formulation< T_Scalar >::addInterfaceField(field::InterfaceFieldInterface< T_Scalar > &interfaceField)
     {
       _interfaceFields.push_back(&interfaceField);
+      gmshfem::msg::debug << "New InterfaceField named " << _interfaceFields.back()->name() << " is added to the formulation '" << _name << "'." << gmshfem::msg::endl;
     }
 
     template< class T_Scalar >
@@ -177,7 +178,7 @@ namespace gmshddm
       gmshfem::common::Timer time;
       time.tick();
       
-      gmshfem::msg::info << "Solve " << _name << "..." << gmshfem::msg::endl;
+      gmshfem::msg::info << "Solve " << _name << " using " << solver << "..." << gmshfem::msg::endl;
 
       // ******************************************************
       // Iterative solver
@@ -248,6 +249,21 @@ namespace gmshddm
         KSPGMRESSetRestart(ksp, 10000);
         KSPSolve(ksp, _rhs.getPetsc(), Y);
       }
+      else if(solver == "bicgl") { // bicgl
+        _IA = true;
+        KSP ksp;
+        PC pc;
+        KSPCreate(PETSC_COMM_SELF, &ksp);
+        KSPSetFromOptions(ksp);
+        KSPSetOperators(ksp, A, A);
+        KSPSetType(ksp, "bcgsl");
+        KSPMonitorSet(ksp, kspPrint< T_Scalar >, this, PETSC_NULL);
+        KSPSetTolerances(ksp, tolerance, PETSC_DEFAULT, PETSC_DEFAULT, iterMax);
+        KSPGetPC(ksp, &pc);
+        PCSetType(pc, PCNONE);
+        // KSPBCGSLSetEll(ksp, 2);
+        KSPSolve(ksp, _rhs.getPetsc(), Y);
+      }
 
       // ******************************************************
       // Final solution