New example with cross-point treatment (PML and HABC)

examples/helmholtz/crossPoints/CMakeLists.txt

+cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
+
+project(demo CXX)
+
+include(CheckCXXCompilerFlag)
+check_cxx_compiler_flag("-std=c++11" HAVE_CXX11)
+if(HAVE_CXX11)
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+else()
+  message(ERROR "The compiler ${CMAKE_CXX_COMPILER} has no C++11 support. Please use a different C++ compiler.")
+endif()
+
+check_cxx_compiler_flag("-Wall" HAVE_WALL)
+if(HAVE_WALL)
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall")
+endif()
+
+macro(find_all_libraries VARNAME LISTNAME PATH SUFFIX)
+  set(${VARNAME})
+  list(LENGTH ${LISTNAME} NUM_LIST)
+  foreach(LIB ${${LISTNAME}})
+    if("${PATH}" STREQUAL "")
+      find_library(FOUND_LIB ${LIB} PATH_SUFFIXES ${SUFFIX})
+    else()
+      find_library(FOUND_LIB ${LIB} PATHS ${PATH} NO_DEFAULT_PATH)
+    endif()
+    if(FOUND_LIB)
+      list(APPEND ${VARNAME} ${FOUND_LIB})
+    endif()
+    unset(FOUND_LIB CACHE)
+  endforeach()
+  list(LENGTH ${VARNAME} NUM_FOUND_LIBRARIES)
+  if(NUM_FOUND_LIBRARIES LESS NUM_LIST)
+    set(${VARNAME})
+  endif()
+endmacro()
+
+
+macro(set_config_option STRING VARNAME)
+  set(${VARNAME} TRUE)
+  list(APPEND CONFIG_OPTIONS ${STRING})
+  message(STATUS "Found " ${STRING})
+endmacro()
+
+macro(opt OPTION HELP VALUE)
+  option(ENABLE_${OPTION} ${HELP} ${VALUE})
+  set(OPT_TEXI "${OPT_TEXI}\n@item ENABLE_${OPTION}\n${HELP} (default: ${VALUE})")
+endmacro()
+
+opt(MPI "Enable MPI" OFF)
+
+####
+#   Mandatory libs
+####
+
+# Gmsh (Mandatory)
+
+find_library(GMSH_LIB gmsh)
+if(NOT GMSH_LIB)
+  message(FATAL_ERROR "Could not find libgmsh")
+endif(NOT GMSH_LIB)
+
+find_path(GMSH_INC gmsh.h)
+if(NOT GMSH_INC)
+  message(FATAL_ERROR "Could not find gmsh.h")
+endif(NOT GMSH_INC)
+
+list(APPEND EXTRA_INCS ${GMSH_INC})
+list(APPEND EXTRA_LIBS ${GMSH_LIB})
+
+# GmshFem
+
+find_library(GMSHFEM_LIB gmshfem)
+if(NOT GMSHFEM_LIB)
+  message(FATAL_ERROR "Could not find libgmshfem")
+endif(NOT GMSHFEM_LIB)
+
+find_path(GMSHFEM_INC gmshfem/GmshFem.h)
+if(NOT GMSHFEM_INC)
+  message(FATAL_ERROR "Could not find 'gmshfem/GmshFem.h'")
+endif(NOT GMSHFEM_INC)
+
+list(APPEND EXTRA_LIBS ${GMSHFEM_LIB})
+list(APPEND EXTRA_INCS ${GMSHFEM_INC}/gmshfem/contrib)
+list(APPEND EXTRA_INCS ${GMSHFEM_INC})
+
+# GmshDdm
+
+find_library(GMSHDDM_LIB gmshddm)
+if(NOT GMSHDDM_LIB)
+  message(FATAL_ERROR "Could not find libgmshddm")
+endif(NOT GMSHDDM_LIB)
+
+find_path(GMSHDDM_INC gmshddm/GmshDdm.h)
+if(NOT GMSHDDM_INC)
+  message(FATAL_ERROR "Could not find 'gmshddm/GmshDdm.h'")
+endif(NOT GMSHDDM_INC)
+
+list(APPEND EXTRA_LIBS ${GMSHDDM_LIB})
+list(APPEND EXTRA_INCS ${GMSHDDM_INC})
+
+# Eigen
+
+find_package(Eigen3 REQUIRED)
+list(APPEND EXTRA_INCS "${EIGEN3_INCLUDE_DIR}")
+
+# OpenMP
+
+find_package(OpenMP)
+if(OPENMP_FOUND)
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+endif(OPENMP_FOUND)
+
+if(APPLE AND EXISTS "/usr/local/opt/libomp")
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Xpreprocessor -fopenmp -I/usr/local/opt/libomp/include")
+  list(APPEND EXTRA_LIBS "-L/usr/local/opt/libomp/lib -lomp")
+endif()
+
+# MPI
+
+if(ENABLE_MPI)
+  find_package(MPI)
+  if(MPI_FOUND)
+    list(APPEND EXTRA_INCS ${MPI_CXX_INCLUDE_PATH})
+    list(APPEND EXTRA_LIBS ${MPI_CXX_LIBRARIES})
+    set(CMAKE_C_COMPILER ${MPI_C_COMPILER})
+    set(CMAKE_CXX_COMPILER ${MPI_CXX_COMPILER})
+    set(CMAKE_Fortran_COMPILER ${MPI_Fortran_COMPILER})
+    set_config_option("MPI" HAVE_MPI)
+  endif()
+endif()
+
+# PETSc
+
+if(PETSC_DIR)
+set(ENV_PETSC_DIR ${PETSC_DIR})
+else()
+  set(ENV_PETSC_DIR $ENV{PETSC_DIR})
+endif()
+
+if(PETSC_ARCH)
+  set(ENV_PETSC_ARCH ${PETSC_ARCH})
+else()
+  set(ENV_PETSC_ARCH $ENV{PETSC_ARCH})
+endif()
+
+set(PETSC_POSSIBLE_CONF_FILES
+    ${ENV_PETSC_DIR}/${ENV_PETSC_ARCH}/conf/petscvariables
+    ${ENV_PETSC_DIR}/${ENV_PETSC_ARCH}/lib/petsc-conf/petscvariables
+    ${ENV_PETSC_DIR}/${ENV_PETSC_ARCH}/lib/petsc/conf/petscvariables)
+
+foreach(FILE ${PETSC_POSSIBLE_CONF_FILES})
+  if(EXISTS ${FILE})
+    # old-style PETSc installations (using PETSC_DIR and PETSC_ARCH)
+    message(STATUS "Using PETSc dir: ${ENV_PETSC_DIR}")
+    message(STATUS "Using PETSc arch: ${ENV_PETSC_ARCH}")
+    # find includes by parsing the petscvariables file
+    file(STRINGS ${FILE} PETSC_VARIABLES NEWLINE_CONSUME)
+  endif()
+endforeach()
+
+if(PETSC_VARIABLES)
+  # try to find PETSC_CC_INCLUDES for PETSc >= 3.4
+  string(REGEX MATCH "PETSC_CC_INCLUDES = [^\n\r]*" PETSC_PACKAGES_INCLUDES ${PETSC_VARIABLES})
+  if(PETSC_PACKAGES_INCLUDES)
+    string(REPLACE "PETSC_CC_INCLUDES = " "" PETSC_PACKAGES_INCLUDES ${PETSC_PACKAGES_INCLUDES})
+  else()
+    # try to find PETSC_PACKAGES_INCLUDES in older versions
+    list(APPEND EXTRA_INCS ${ENV_PETSC_DIR}/include)
+    list(APPEND EXTRA_INCS ${ENV_PETSC_DIR}/${ENV_PETSC_ARCH}/include)
+    string(REGEX MATCH "PACKAGES_INCLUDES = [^\n\r]*" PETSC_PACKAGES_INCLUDES ${PETSC_VARIABLES})
+    string(REPLACE "PACKAGES_INCLUDES = " "" PETSC_PACKAGES_INCLUDES ${PETSC_PACKAGES_INCLUDES})
+  endif()
+
+  if(PETSC_PACKAGES_INCLUDES)
+    if(PETSC_PACKAGES_INCLUDES)
+      string(REPLACE "-I" "" PETSC_PACKAGES_INCLUDES ${PETSC_PACKAGES_INCLUDES})
+      string(REPLACE " " ";" PETSC_PACKAGES_INCLUDES ${PETSC_PACKAGES_INCLUDES})
+      foreach(VAR ${PETSC_PACKAGES_INCLUDES})
+        list(APPEND EXTRA_INCS ${VAR})
+      endforeach()
+    endif()
+  endif()
+
+  # find libraries (<= 3.0)
+  set(PETSC_LIBS_REQUIRED petscksp petscdm petscmat petscvec petsc)
+  find_all_libraries(PETSC_LIBS PETSC_LIBS_REQUIRED ${ENV_PETSC_DIR}/${ENV_PETSC_ARCH}/lib "")
+  # petsc 3.1 creates only one library (libpetsc)
+  if(NOT PETSC_LIBS)
+    find_library(PETSC_LIBS petsc PATHS ${ENV_PETSC_DIR}/${ENV_PETSC_ARCH}/lib NO_DEFAULT_PATH)
+  endif()
+
+  list(APPEND EXTRA_LIBS ${PETSC_LIBS})
+
+  string(REGEX MATCH "PACKAGES_LIBS = [^\n\r]*" PLIBS ${PETSC_VARIABLES})
+  if(PLIBS)
+    string(REPLACE "PACKAGES_LIBS = " "" PLIBS ${PLIBS})
+    string(STRIP ${PLIBS} PLIBS)
+    list(APPEND EXTRA_LIBS "${PLIBS}")
+  endif()
+
+  string(REGEX MATCH "PETSC_EXTERNAL_LIB_BASIC = [^\n\r]*" PLIBS_BASIC ${PETSC_VARIABLES})
+  if(PLIBS_BASIC)
+    string(REPLACE "PETSC_EXTERNAL_LIB_BASIC = " "" PLIBS_BASIC ${PLIBS_BASIC})
+    string(STRIP ${PLIBS_BASIC} PLIBS_BASIC)
+    list(APPEND EXTRA_LIBS "${PLIBS_BASIC}")
+  endif()
+
+  string(REGEX MATCH "PCC_LINKER_LIBS = [^\n\r]*" LLIBS ${PETSC_VARIABLES})
+  if(LLIBS)
+    string(REPLACE "PCC_LINKER_LIBS = " "" LLIBS ${LLIBS})
+    string(STRIP ${LLIBS} LLIBS)
+    list(APPEND EXTRA_LIBS "${LLIBS}")
+  endif()
+
+  set_config_option("PETSc" HAVE_PETSC)
+  list(APPEND EXTRA_LIBS ${PETSC_LIBS})
+else()
+  find_library(PETSC_LIBS petsc)
+  find_path(PETSC_INC "petsc.h" PATH_SUFFIXES include/petsc)
+  if(PETSC_LIBS AND PETSC_INC)
+    set_config_option("PETSc" HAVE_PETSC)
+    list(APPEND EXTRA_LIBS ${PETSC_LIBS})
+    list(APPEND EXTRA_INCS ${PETSC_INC})
+  endif()
+endif()
+
+include_directories(${EXTRA_INCS})
+add_executable(demo main.cpp mesh.cpp ddm2D.cpp ddm3D.cpp Subproblem2D.cpp Subproblem3D.cpp SubproblemDomains.cpp SubproblemParameters.cpp monoDomain.cpp ${EXTRA_INCS})
+target_link_libraries(demo ${EXTRA_LIBS})
+

examples/helmholtz/crossPoints/README.md

+# Helmholtz cross-points
+> Anthony Royer
+
+Dependency(ies):
+* GmshDDM : 1.0.0
+* GmshFEM : 1.0.0
+* ...
+
+## Project description
+
+In this work, we present a non-overlapping substructured DDM with PML transmission conditions or HABC 
+transmission conditionz forvcheckerboard (Cartesian) decompositions that takes cross-points into account. 
+In such decompositions, each subdomain is surrounded by PMLs associated to edges and corners. 
+
+## Installation
+
+```
+mkdir build && cd build
+cmake ..
+make
+sudo make install
+```
+
+## Usage
+
+Simply run:
+```
+  ./example [PARAM]
+```
+
+## References
+
+A. Royer, C. Geuzaine, E. Béchet and A. Modave, "A non-overlapping domain decomposition method with perfectly matched layer transmission conditions for the Helmholtz equation", in review, 2021

examples/helmholtz/crossPoints/Subproblem2D.h

+#ifndef H_SUBPROBLEM2D
+#define H_SUBPROBLEM2D
+
+#include "SubproblemDomains.h"
+#include "SubproblemParameters.h"
+
+#include <gmshfem/Formulation.h>
+
+namespace D2 {
+
+
+  class Boundary
+  {
+   protected:
+    const unsigned int _boundary;
+    
+   public:
+    Boundary(const unsigned int boundary);
+    virtual ~Boundary();
+    
+    virtual void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters) = 0;
+    
+    std::string orientation() const;
+  };
+
+  class Corner
+  {
+   protected:
+    const unsigned int _corner;
+    Boundary *_bnd[2];
+    
+   public:
+    Corner(const unsigned int corner, Boundary *first, Boundary *second);
+    virtual ~Corner();
+    
+    virtual void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) = 0;
+    
+    std::string orientation() const;
+  };
+
+  class Subproblem
+  {
+   private:
+    const SubproblemDomains _domains;
+    const SubproblemParameters _parameters;
+    gmshfem::problem::Formulation< std::complex< double > > &_formulation;
+    const std::string _fieldName;
+    std::vector< Boundary * > _boundary;
+    std::vector< Corner * > _corner;
+    
+    void _parseParameters(std::string &method, std::vector< std::string > &parameters, const std::string &str) const;
+    
+   public:
+    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);
+    ~Subproblem();
+    
+    void writeFormulation();
+    
+    Boundary *getBoundary(const unsigned int b) const;
+    Corner *getCorner(const unsigned int c) const;
+  };
+
+
+
+
+  class Sommerfeld : public Boundary
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _v;
+    
+   public:
+    Sommerfeld(const unsigned int boundary);
+    ~Sommerfeld();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getV() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+
+  class HABC : public Boundary
+  {
+   protected:
+    const unsigned int _N;
+    const double _theta;
+    std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * > _uHABC;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _v;
+    
+   public:
+    HABC(const unsigned int boundary, const unsigned int N, const double theta);
+    ~HABC();
+    
+    unsigned int getN() const;
+    double getTheta() const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getUHABC(const unsigned int m) const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getV() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+  
+  class Pml : public Boundary
+  {
+   protected:
+    const double _size;
+    const std::string _type;
+    
+   public:
+    Pml(const unsigned int boundary, const double pmlSize, const std::string type);
+    virtual ~Pml();
+    
+    virtual double getSize() const;
+    virtual gmshfem::function::ScalarFunction< std::complex< double > >  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;
+    
+    virtual void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters) = 0;
+  };
+
+  class PmlContinuous : public Pml // with Lagrange multiplier
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _uPml;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _v;
+    
+   public:
+    PmlContinuous(const unsigned int boundary, const double pmlSize, const std::string type);
+    ~PmlContinuous();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getUPml() const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getV() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+
+  class PmlDiscontinuous : public Pml // with Lagrange multiplier
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _uPml;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form2 > * _v;
+    
+   public:
+    PmlDiscontinuous(const unsigned int boundary, const double pmlSize, const std::string type);
+    ~PmlDiscontinuous();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getUPml() const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form2 >* getV() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+  
+  class PmlWithoutMultiplier : public Pml // without Lagrange multiplier
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _uPml;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _v;
+    
+   public:
+    PmlWithoutMultiplier(const unsigned int boundary, const double pmlSize, const std::string type);
+    ~PmlWithoutMultiplier();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getUPml() const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getV() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+
+
+
+
+  class HABC_HABC : public Corner
+  {
+   protected:
+    std::vector< std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * > > _uHABC_C;
+    std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * > _vC[2];
+    
+   public:
+    HABC_HABC(const unsigned int corner, HABC *first, HABC *second);
+    ~HABC_HABC();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *getV(const unsigned int m, const unsigned int i) const;
+    
+    HABC *firstBoundary() const;
+    HABC *secondBoundary() const;
+      
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+
+  class PmlContinuous_PmlContinuous : public Corner
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _uPmlCorner;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _vC[2];
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _vCorner;
+    
+   public:
+    PmlContinuous_PmlContinuous(const unsigned int corner, PmlContinuous *first, PmlContinuous *second);
+    ~PmlContinuous_PmlContinuous();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getUPml() const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getV(const unsigned int i) const;
+    
+    PmlContinuous *firstBoundary() const;
+    PmlContinuous *secondBoundary() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+
+  class PmlDiscontinuous_PmlDiscontinuous : public Corner
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _uPmlCorner;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form2 > * _vC[2];
+    
+   public:
+    PmlDiscontinuous_PmlDiscontinuous(const unsigned int corner, PmlDiscontinuous *first, PmlDiscontinuous *second);
+    ~PmlDiscontinuous_PmlDiscontinuous();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getUPml() const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form2 >* getV(const unsigned int i) const;
+    
+    PmlDiscontinuous *firstBoundary() const;
+    PmlDiscontinuous *secondBoundary() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+  
+  class PmlWithoutMultiplier_PmlWithoutMultiplier : public Corner
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _uPmlCorner;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _vC[2];
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _vCorner;
+    
+   public:
+    PmlWithoutMultiplier_PmlWithoutMultiplier(const unsigned int corner, PmlWithoutMultiplier *first, PmlWithoutMultiplier *second);
+    ~PmlWithoutMultiplier_PmlWithoutMultiplier();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getUPml() const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getV(const unsigned int i) const;
+    
+    PmlWithoutMultiplier *firstBoundary() const;
+    PmlWithoutMultiplier *secondBoundary() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+
+  class PmlContinuous_Sommerfeld : public Corner
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _vC;
+    
+   public:
+    PmlContinuous_Sommerfeld(const unsigned int corner, PmlContinuous *first, Sommerfeld *second);
+    ~PmlContinuous_Sommerfeld();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getV() const;
+    
+    PmlContinuous *firstBoundary() const;
+    Sommerfeld *secondBoundary() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+
+  class Sommerfeld_PmlContinuous : public Corner
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _vC;
+    
+   public:
+    Sommerfeld_PmlContinuous(const unsigned int corner, Sommerfeld *first, PmlContinuous *second);
+    ~Sommerfeld_PmlContinuous();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getV() const;
+    
+    Sommerfeld *firstBoundary() const;
+    PmlContinuous *secondBoundary() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+  
+  class PmlDiscontinuous_Sommerfeld : public Corner
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form2 > * _vC;
+    
+   public:
+    PmlDiscontinuous_Sommerfeld(const unsigned int corner, PmlDiscontinuous *first, Sommerfeld *second);
+    ~PmlDiscontinuous_Sommerfeld();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form2 >* getV() const;
+    
+    PmlDiscontinuous *firstBoundary() const;
+    Sommerfeld *secondBoundary() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+
+  class Sommerfeld_PmlDiscontinuous : public Corner
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form2 > * _vC;
+    
+   public:
+    Sommerfeld_PmlDiscontinuous(const unsigned int corner, Sommerfeld *first, PmlDiscontinuous *second);
+    ~Sommerfeld_PmlDiscontinuous();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form2 >* getV() const;
+    
+    Sommerfeld *firstBoundary() const;
+    PmlDiscontinuous *secondBoundary() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+  
+  class HABC_Sommerfeld : public Corner
+  {
+   protected:
+    std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * > _vC;
+    
+   public:
+    HABC_Sommerfeld(const unsigned int corner, HABC *first, Sommerfeld *second);
+    ~HABC_Sommerfeld();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *getV(const unsigned int m) const;
+    
+    HABC *firstBoundary() const;
+    Sommerfeld *secondBoundary() const;
+      
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+  
+  class Sommerfeld_HABC : public Corner
+  {
+   protected:
+    std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * > _vC;
+    
+   public:
+    Sommerfeld_HABC(const unsigned int corner, Sommerfeld *first, HABC *second);
+    ~Sommerfeld_HABC();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *getV(const unsigned int m) const;
+    
+    Sommerfeld *firstBoundary() const;
+    HABC *secondBoundary() const;
+      
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+
+
+}
+
+#endif // H_SUBPROBLEM2D

examples/helmholtz/crossPoints/Subproblem3D.h

+#ifndef H_SUBPROBLEM3D
+#define H_SUBPROBLEM3D
+
+#include "SubproblemDomains.h"
+#include "SubproblemParameters.h"
+
+#include <gmshfem/Formulation.h>
+
+namespace D3 {
+
+
+  class Boundary
+  {
+   protected:
+    const unsigned int _boundary;
+    
+   public:
+    Boundary(const unsigned int boundary);
+    virtual ~Boundary();
+    
+    virtual void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters) = 0;
+    
+    std::string orientation() const;
+  };
+  
+  class Edge
+  {
+   protected:
+    const unsigned int _edge;
+    Boundary *_bnd[2];
+    
+   public:
+    Edge(const unsigned int edge, Boundary *first, Boundary *second);
+    virtual ~Edge();
+    
+    virtual void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) = 0;
+    
+    std::string orientation() const;
+  };
+
+  class Corner
+  {
+   protected:
+    const unsigned int _corner;
+    Edge *_bnd[3];
+    
+   public:
+    Corner(const unsigned int corner, Edge *first, Edge *second, Edge *third);
+    virtual ~Corner();
+    
+    virtual void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) = 0;
+    
+    std::string orientation() const;
+  };
+
+  class Subproblem
+  {
+   private:
+    const SubproblemDomains _domains;
+    const SubproblemParameters _parameters;
+    gmshfem::problem::Formulation< std::complex< double > > &_formulation;
+    const std::string _fieldName;
+    std::vector< Boundary * > _boundary;
+    std::vector< Edge * > _edge;
+    std::vector< Corner * > _corner;
+    
+    void _parseParameters(std::string &method, std::vector< std::string > &parameters, const std::string &str) const;
+    
+   public:
+    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);
+    ~Subproblem();
+    
+    void writeFormulation();
+    
+    Boundary *getBoundary(const unsigned int b) const;
+    Edge *getEdge(const unsigned int e) const;
+    Corner *getCorner(const unsigned int c) const;
+  };
+
+
+
+
+  class Sommerfeld : public Boundary
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _v;
+    
+   public:
+    Sommerfeld(const unsigned int boundary);
+    ~Sommerfeld();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getV() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+  
+  class HABC : public Boundary
+  {
+   protected:
+    const unsigned int _N;
+    const double _theta;
+    std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * > _uHABC;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _v;
+    
+   public:
+    HABC(const unsigned int boundary, const unsigned int N, const double theta);
+    ~HABC();
+    
+    unsigned int getN() const;
+    double getTheta() const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getUHABC(const unsigned int m) const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getV() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+  
+  class Pml : public Boundary
+  {
+   protected:
+    const double _size;
+    const std::string _type;
+    
+   public:
+    Pml(const unsigned int boundary, const double pmlSize, const std::string type);
+    virtual ~Pml();
+    
+    double getSize() const;
+    gmshfem::function::ScalarFunction< std::complex< double > >  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;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters) = 0;
+  };
+  
+  class PmlContinuous : public Pml
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _uPml;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _v;
+    
+   public:
+    PmlContinuous(const unsigned int boundary, const double pmlSize, const std::string type);
+    ~PmlContinuous();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getUPml() const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getV() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+  
+  // Marche pas, il faudrait du vrai H(div)
+  class PmlDiscontinuous : public Pml
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _uPml;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form2 > * _v;
+    
+   public:
+    PmlDiscontinuous(const unsigned int boundary, const double pmlSize, const std::string type);
+    ~PmlDiscontinuous();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getUPml() const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form2 >* getV() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+  
+  class PmlWithoutMultiplier : public Pml
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _uPml;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _v;
+    
+   public:
+    PmlWithoutMultiplier(const unsigned int boundary, const double pmlSize, const std::string type);
+    ~PmlWithoutMultiplier();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getUPml() const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getV() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const std::string &fieldName, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+
+
+
+
+  class HABC_HABC : public Edge
+  {
+   protected:
+    std::vector< std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * > > _uHABC_E;
+    std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * > _vE[2];
+    
+   public:
+    HABC_HABC(const unsigned int corner, HABC *first, HABC *second);
+    ~HABC_HABC();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *getUHABC(const unsigned int m, const unsigned int n) const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *getV(const unsigned int m, const unsigned int i) const;
+    
+    HABC *firstBoundary() const;
+    HABC *secondBoundary() const;
+      
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+  
+  class PmlContinuous_PmlContinuous : public Edge
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _uPmlEdge;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _vC[2];
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _vEdge;
+    
+   public:
+    PmlContinuous_PmlContinuous(const unsigned int edge, PmlContinuous *first, PmlContinuous *second);
+    ~PmlContinuous_PmlContinuous();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getUPml() const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getV(const unsigned int i) const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getVEdge() const;
+    
+    PmlContinuous *firstBoundary() const;
+    PmlContinuous *secondBoundary() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+  
+  
+  
+  
+  class HABC_HABC_HABC : public Corner
+  {
+   protected:
+    std::vector< std::vector< std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * > > > _uHABC_C;
+    std::vector< std::vector< gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * > > _vC[3];
+    
+   public:
+    HABC_HABC_HABC(const unsigned int corner, HABC_HABC *first, HABC_HABC *second, HABC_HABC *third);
+    ~HABC_HABC_HABC();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *getUHABC(const unsigned int m, const unsigned int n, const unsigned int o) const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > *getV(const unsigned int m, const unsigned int n, const unsigned int i) const;
+    
+    HABC_HABC *firstEdge() const;
+    HABC_HABC *secondEdge() const;
+    HABC_HABC *thirdEdge() const;
+      
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+  
+  class PmlContinuous_PmlContinuous_PmlContinuous : public Corner
+  {
+   protected:
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _uPmlCorner;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _vC[3];
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _vEdge[3];
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 > * _vCorner;
+    
+   public:
+    PmlContinuous_PmlContinuous_PmlContinuous(const unsigned int corner, PmlContinuous_PmlContinuous *first, PmlContinuous_PmlContinuous *second, PmlContinuous_PmlContinuous *third);
+    ~PmlContinuous_PmlContinuous_PmlContinuous();
+    
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getUPml() const;
+    gmshfem::field::Field< std::complex< double >, gmshfem::field::Form::Form0 >* getV(const unsigned int i) const;
+    
+    PmlContinuous_PmlContinuous *firstEdge() const;
+    PmlContinuous_PmlContinuous *secondEdge() const;
+    PmlContinuous_PmlContinuous *thirdEdge() const;
+    
+    void writeFormulation(gmshfem::problem::Formulation< std::complex< double > > &formulation, const SubproblemDomains &domains, const SubproblemParameters &parameters) override;
+  };
+
+}
+
+#endif // H_SUBPROBLEM3D

examples/helmholtz/crossPoints/SubproblemDomains.cpp

+#include "SubproblemDomains.h"
+
+namespace D2 {
+
+
+  SubproblemDomains::SubproblemDomains() : _omega(), _sigma(), _pml(), _pmlBnd(), _pmlCorner(), _corner()
+  {
+
+  }
+
+  SubproblemDomains::~SubproblemDomains()
+  {
+  }
+    
+  void SubproblemDomains::setOmega(const gmshfem::domain::Domain &omega)
+  {
+    _omega = omega;
+  }
+
+  void SubproblemDomains::setSigma(const std::vector< gmshfem::domain::Domain > &sigma)
+  {
+    _sigma = sigma;
+  }
+
+  gmshfem::domain::Domain SubproblemDomains::getSigma(const unsigned int b) const
+  {
+    return _sigma[b];
+  }
+
+  void SubproblemDomains::setPml(const std::vector< gmshfem::domain::Domain > &pml)
+  {
+    _pml = pml;
+  }
+
+  gmshfem::domain::Domain SubproblemDomains::getPml(const unsigned int b) const
+  {
+    return _pml[b];
+  }
+
+  void SubproblemDomains::setPmlBnd(const std::vector< std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > > &pmlBnd)
+  {
+    _pmlBnd = pmlBnd;
+  }
+
+  std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > SubproblemDomains::getPmlBnd(const unsigned int c) const
+  {
+    return _pmlBnd[c];
+  }
+
+  void SubproblemDomains::setPmlCorner(const std::vector< gmshfem::domain::Domain > &pmlCorner)
+  {
+    _pmlCorner = pmlCorner;
+  }
+
+  gmshfem::domain::Domain SubproblemDomains::getPmlCorner(const unsigned int c) const
+  {
+    return _pmlCorner[c];
+  }
+
+  void SubproblemDomains::setCorner(const std::vector< gmshfem::domain::Domain > &corner)
+  {
+    _corner = corner;
+  }
+
+  gmshfem::domain::Domain SubproblemDomains::getCorner(const unsigned int c) const
+  {
+    return _corner[c];
+  }
+
+
+}
+
+
+namespace D3 {
+
+
+  SubproblemDomains::SubproblemDomains() : _omega(), _sigma(), _pml(), _pmlBnd(), _pmlEdge(), _edge(), _pmlEdgeBnd(), _pmlCorner(), _corner()
+  {
+
+  }
+
+  SubproblemDomains::~SubproblemDomains()
+  {
+  }
+    
+  void SubproblemDomains::setOmega(const gmshfem::domain::Domain &omega)
+  {
+    _omega = omega;
+  }
+
+  void SubproblemDomains::setSigma(const std::vector< gmshfem::domain::Domain > &sigma)
+  {
+    _sigma = sigma;
+  }
+
+  gmshfem::domain::Domain SubproblemDomains::getSigma(const unsigned int b) const
+  {
+    return _sigma[b];
+  }
+
+  void SubproblemDomains::setPml(const std::vector< gmshfem::domain::Domain > &pml)
+  {
+    _pml = pml;
+  }
+
+  gmshfem::domain::Domain SubproblemDomains::getPml(const unsigned int b) const
+  {
+    return _pml[b];
+  }
+
+  void SubproblemDomains::setPmlBnd(const std::vector< std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > > &pmlBnd)
+  {
+    _pmlBnd = pmlBnd;
+  }
+
+  std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > SubproblemDomains::getPmlBnd(const unsigned int e) const
+  {
+    return _pmlBnd[e];
+  }
+
+  void SubproblemDomains::setPmlEdge(const std::vector< gmshfem::domain::Domain > &pmlEdge)
+  {
+    _pmlEdge = pmlEdge;
+  }
+
+  gmshfem::domain::Domain SubproblemDomains::getPmlEdge(const unsigned int e) const
+  {
+    return _pmlEdge[e];
+  }
+
+  void SubproblemDomains::setEdge(const std::vector< gmshfem::domain::Domain > &edge)
+  {
+    _edge = edge;
+  }
+
+  gmshfem::domain::Domain SubproblemDomains::getEdge(const unsigned int e) const
+  {
+    return _edge[e];
+  }
+  
+  void SubproblemDomains::setPmlEdgeBnd(const std::vector< std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > > &pmlEdgeBnd)
+  {
+    _pmlEdgeBnd = pmlEdgeBnd;
+  }
+
+  std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > SubproblemDomains::getPmlEdgeBnd(const unsigned int c) const
+  {
+    return _pmlEdgeBnd[c];
+  }
+
+  void SubproblemDomains::setPmlCorner(const std::vector< gmshfem::domain::Domain > &pmlCorner)
+  {
+    _pmlCorner = pmlCorner;
+  }
+
+  gmshfem::domain::Domain SubproblemDomains::getPmlCorner(const unsigned int c) const
+  {
+    return _pmlCorner[c];
+  }
+
+  void SubproblemDomains::setCorner(const std::vector< gmshfem::domain::Domain > &corner)
+  {
+    _corner = corner;
+  }
+
+  gmshfem::domain::Domain SubproblemDomains::getCorner(const unsigned int c) const
+  {
+    return _corner[c];
+  }
+  
+  void SubproblemDomains::setCornerEdge(const std::vector< std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > > &cornerEdge)
+  {
+    _cornerEdge = cornerEdge;
+  }
+  
+  std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > SubproblemDomains::getCornerEdge(const unsigned int c) const
+  {
+    return _cornerEdge[c];
+  }
+  
+  void SubproblemDomains::saveDebug() const
+  {
+    const std::string ori[6] = {"E", "N", "W", "S", "D", "U"};
+    
+    _omega.saveDebug("omega");
+    for(unsigned int b = 0; b < 6; ++b) {
+      _sigma[b].saveDebug("sigma" + ori[b]);
+      _pml[b].saveDebug("pml" + ori[b]);
+    }
+    for(unsigned int e = 0; e < 12; ++e) {
+      if(e < 4) {
+        _pmlBnd[e].first.saveDebug("pmlBnd" + ori[e%4] + "D_first");
+        _pmlBnd[e].second.saveDebug("pmlBnd" + ori[e%4] + "D_second");
+        _pmlEdge[e].saveDebug("pmlEdge" + ori[e%4] + "D");
+        _edge[e].saveDebug("edge" + ori[e%4] + "D");
+      }
+      else if(e >= 4 && e < 8) {
+        _pmlBnd[e].first.saveDebug("pmlBnd" + ori[e%4] + ori[(e+1)%4] + "_first");
+        _pmlBnd[e].second.saveDebug("pmlBnd" + ori[e%4] + ori[(e+1)%4] + "_second");
+        _pmlEdge[e].saveDebug("pmlEdge" + ori[e%4] + ori[(e+1)%4]);
+        _edge[e].saveDebug("edge" + ori[e%4] + ori[(e+1)%4]);
+      }
+      else {
+        _pmlBnd[e].first.saveDebug("pmlBnd" + ori[e%4] + "U_first");
+        _pmlBnd[e].second.saveDebug("pmlBnd" + ori[e%4] + "U_second");
+        _pmlEdge[e].saveDebug("pmlEdge" + ori[e%4] + "U");
+        _edge[e].saveDebug("edge" + ori[e%4] + "U");
+      }
+    }
+    for(unsigned int c = 0; c < 8; ++c) {
+      if(c < 4) {
+        std::get<0>(_pmlEdgeBnd[c]).saveDebug("pmlEdgeBnd" + ori[c%4] + ori[(c+1)%4] + "D_first");
+        std::get<1>(_pmlEdgeBnd[c]).saveDebug("pmlEdgeBnd" + ori[c%4] + ori[(c+1)%4] + "D_second");
+        std::get<2>(_pmlEdgeBnd[c]).saveDebug("pmlEdgeBnd" + ori[c%4] + ori[(c+1)%4] + "D_third");
+        _pmlCorner[c].saveDebug("pmlCorner" + ori[c%4] + ori[(c+1)%4] + "D");
+        _corner[c].saveDebug("corner" + ori[c%4] + ori[(c+1)%4] + "D");
+        std::get<0>(_cornerEdge[c]).saveDebug("cornerEdge" + ori[c%4] + ori[(c+1)%4] + "D_first");
+        std::get<1>(_cornerEdge[c]).saveDebug("cornerEdge" + ori[c%4] + ori[(c+1)%4] + "D_second");
+        std::get<2>(_cornerEdge[c]).saveDebug("cornerEdge" + ori[c%4] + ori[(c+1)%4] + "D_third");
+      }
+      else {
+        std::get<0>(_pmlEdgeBnd[c]).saveDebug("pmlEdgeBnd" + ori[c%4] + ori[(c+1)%4] + "U_first");
+        std::get<1>(_pmlEdgeBnd[c]).saveDebug("pmlEdgeBnd" + ori[c%4] + ori[(c+1)%4] + "U_second");
+        std::get<2>(_pmlEdgeBnd[c]).saveDebug("pmlEdgeBnd" + ori[c%4] + ori[(c+1)%4] + "U_third");
+        _pmlCorner[c].saveDebug("pmlCorner" + ori[c%4] + ori[(c+1)%4] + "U");
+        _corner[c].saveDebug("corner" + ori[c%4] + ori[(c+1)%4] + "U");
+        std::get<0>(_cornerEdge[c]).saveDebug("cornerEdge" + ori[c%4] + ori[(c+1)%4] + "U_first");
+        std::get<1>(_cornerEdge[c]).saveDebug("cornerEdge" + ori[c%4] + ori[(c+1)%4] + "U_second");
+        std::get<2>(_cornerEdge[c]).saveDebug("cornerEdge" + ori[c%4] + ori[(c+1)%4] + "U_third");
+      }
+    }
+  }
+
+
+}

examples/helmholtz/crossPoints/SubproblemDomains.h

+#ifndef H_SUBPROBLEMDOMAINS
+#define H_SUBPROBLEMDOMAINS
+
+#include <gmshfem/Domain.h>
+
+#include <tuple>
+
+namespace D2 {
+
+
+  class SubproblemDomains
+  {
+   private:
+    gmshfem::domain::Domain _omega;
+    // boundary
+    std::vector< gmshfem::domain::Domain > _sigma;
+    std::vector< gmshfem::domain::Domain > _pml;
+    // corner
+    std::vector< std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > > _pmlBnd;
+    std::vector< gmshfem::domain::Domain > _pmlCorner;
+    std::vector< gmshfem::domain::Domain > _corner;
+    
+   public:
+    SubproblemDomains();
+    ~SubproblemDomains();
+    
+    void setOmega(const gmshfem::domain::Domain &omega);
+    
+    void setSigma(const std::vector< gmshfem::domain::Domain > &sigma);
+    gmshfem::domain::Domain getSigma(const unsigned int b) const;
+    
+    void setPml(const std::vector< gmshfem::domain::Domain > &pml);
+    gmshfem::domain::Domain getPml(const unsigned int b) const;
+    
+    void setPmlBnd(const std::vector< std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > > &pmlBnd);
+    std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > getPmlBnd(const unsigned int c) const;
+    
+    void setPmlCorner(const std::vector< gmshfem::domain::Domain > &pmlCorner);
+    gmshfem::domain::Domain getPmlCorner(const unsigned int c) const;
+    
+    void setCorner(const std::vector< gmshfem::domain::Domain > &corner);
+    gmshfem::domain::Domain getCorner(const unsigned int c) const;
+  };
+
+
+}
+
+namespace D3 {
+
+
+  class SubproblemDomains
+  {
+   private:
+    gmshfem::domain::Domain _omega;
+    // boundary
+    std::vector< gmshfem::domain::Domain > _sigma;
+    std::vector< gmshfem::domain::Domain > _pml;
+    // edge
+    std::vector< std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > > _pmlBnd;
+    std::vector< gmshfem::domain::Domain > _pmlEdge;
+    std::vector< gmshfem::domain::Domain > _edge;
+    // corner
+    std::vector< std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > > _pmlEdgeBnd;
+    std::vector< gmshfem::domain::Domain > _pmlCorner;
+    std::vector< gmshfem::domain::Domain > _corner;
+    std::vector< std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > > _cornerEdge;
+    
+   public:
+    SubproblemDomains();
+    ~SubproblemDomains();
+    
+    void setOmega(const gmshfem::domain::Domain &omega);
+    
+    void setSigma(const std::vector< gmshfem::domain::Domain > &sigma);
+    gmshfem::domain::Domain getSigma(const unsigned int b) const;
+    
+    void setPml(const std::vector< gmshfem::domain::Domain > &pml);
+    gmshfem::domain::Domain getPml(const unsigned int b) const;
+    
+    void setPmlBnd(const std::vector< std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > > &pmlBnd);
+    std::pair< gmshfem::domain::Domain, gmshfem::domain::Domain > getPmlBnd(const unsigned int e) const;
+    
+    void setPmlEdge(const std::vector< gmshfem::domain::Domain > &pmlEdge);
+    gmshfem::domain::Domain getPmlEdge(const unsigned int e) const;
+  
+    void setEdge(const std::vector< gmshfem::domain::Domain > &edge);
+    gmshfem::domain::Domain getEdge(const unsigned int e) const;
+    
+    void setPmlEdgeBnd(const std::vector< std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > > &pmlEdgeBnd);
+    std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > getPmlEdgeBnd(const unsigned int c) const;
+    
+    void setPmlCorner(const std::vector< gmshfem::domain::Domain > &pmlCorner);
+    gmshfem::domain::Domain getPmlCorner(const unsigned int c) const;
+    
+    void setCorner(const std::vector< gmshfem::domain::Domain > &corner);
+    gmshfem::domain::Domain getCorner(const unsigned int c) const;
+    
+    void setCornerEdge(const std::vector< std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > > &cornerEdge);
+    std::tuple< gmshfem::domain::Domain, gmshfem::domain::Domain, gmshfem::domain::Domain > getCornerEdge(const unsigned int c) const;
+    
+    void saveDebug() const;
+  };
+
+
+}
+
+#endif // H_SUBPROBLEMDOMAINS

examples/helmholtz/crossPoints/SubproblemParameters.cpp

+#include "SubproblemParameters.h"
+
+SubproblemParameters::SubproblemParameters()
+{
+
+}
+
+SubproblemParameters::~SubproblemParameters()
+{
+}
+
+void SubproblemParameters::setGauss(const std::string &gauss)
+{
+  _gauss = gauss;
+}
+
+std::string SubproblemParameters::getGauss() const
+{
+  return _gauss;
+}
+
+void SubproblemParameters::setKappa(const gmshfem::function::ScalarFunction< std::complex< double > > &kappa)
+{
+  _kappa = kappa;
+}
+
+gmshfem::function::ScalarFunction< std::complex< double > > SubproblemParameters::getKappa() const
+{
+  return _kappa;
+}
+
+void SubproblemParameters::setNeumannOrder(const unsigned int neumannOrder)
+{
+  _neumannOrder = neumannOrder;
+}
+
+unsigned int SubproblemParameters::getNeumannOrder() const
+{
+  return _neumannOrder;
+}
+
+void SubproblemParameters::setFieldOrder(const unsigned int fieldOrder)
+{
+  _fieldOrder = fieldOrder;
+}
+
+unsigned int SubproblemParameters::getFieldOrder() const
+{
+  return _fieldOrder;
+}
+
+void SubproblemParameters::setStab(const double stab)
+{
+  _stab = stab;
+}
+
+double SubproblemParameters::getStab() const
+{
+  return _stab;
+}

examples/helmholtz/crossPoints/SubproblemParameters.h

+#ifndef H_SUBPROBLEMPARAMETERS
+#define H_SUBPROBLEMPARAMETERS
+
+#include <gmshfem/Function.h>
+
+class SubproblemParameters
+{
+ private:
+  std::string _gauss;
+  gmshfem::function::ScalarFunction< std::complex< double > > _kappa;
+  unsigned int _neumannOrder;
+  unsigned int _fieldOrder;
+  double _stab;
+  
+ public:
+  SubproblemParameters();
+  ~SubproblemParameters();
+  
+  void setGauss(const std::string &gauss);
+  std::string getGauss() const;
+  
+  void setKappa(const gmshfem::function::ScalarFunction< std::complex< double > > &kappa);
+  gmshfem::function::ScalarFunction< std::complex< double > > getKappa() const;
+  
+  void setNeumannOrder(const unsigned int neumannOrder);
+  unsigned int getNeumannOrder() const;
+  
+  void setFieldOrder(const unsigned int fieldOrder);
+  unsigned int getFieldOrder() const;
+  
+  void setStab(const double stab);
+  double getStab() const;
+};
+
+#endif // H_SUBPROBLEMPARAMETERS

examples/helmholtz/crossPoints/ddm2D.h

+#ifndef H_DDM2D
+#define H_DDM2D
+
+namespace D2 {
+
+  void ddm();
+
+}
+
+
+#endif // H_DDM2D

examples/helmholtz/crossPoints/ddm3D.h

+#ifndef H_DDM3D
+#define H_DDM3D
+
+namespace D3 {
+
+  void ddm();
+
+}
+
+
+#endif // H_DDM3D

examples/helmholtz/crossPoints/main.cpp

+#include "monoDomain.h"
+#include "ddm2D.h"
+#include "ddm3D.h"
+
+#include <gmshddm/GmshDdm.h>
+
+#include <gmshfem/Message.h>
+
+int main(int argc, char **argv)
+{
+  gmshddm::common::GmshDdm gmshDdm(argc, argv);
+  
+  std::string problem = "monoDomain2D"; // monoDomain, ddm
+  gmshDdm.userDefinedParameter(problem, "problem");
+  
+  if(problem == "monoDomain2D") {
+    D2::monoDomain();
+  }
+  else if(problem == "monoDomain3D") {
+    D3::monoDomain();
+  }
+  else if(problem == "ddm2D") {
+    D2::ddm();
+  }
+  else if(problem == "ddm3D") {
+    D3::ddm();
+  }
+  else {
+    gmshfem::msg::error << "Unknown 'problem = " << problem << "'" << gmshfem::msg::endl;
+  }
+  
+  return 0;
+}

examples/helmholtz/crossPoints/mesh.h

+#ifndef H_MESH
+#define H_MESH
+
+#include <gmshfem/Function.h>
+
+namespace D2 {
+
+  void subdomain(const int index, const int jndex, const double xSize, const double ySize, const double circleSize, const double lc, const double pmlSizeE, const double pmlSizeN, const double pmlSizeW, const double pmlSizeS, bool transfinite);
+  void monoDomain(const double lc, const double pmlSize, const double R, const int order);
+  void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const double xSize, const double ySize, const double circleSize, const double lc, const bool pml, const double pmlSize, const bool pmlExt, const double pmlSizeExt, const int order, const unsigned int circlePosX = 0, const unsigned int circlePosY = 0, bool transfinite = false);
+
+}
+
+namespace D3 {
+
+  void subdomain(const int index, const int jndex, const int kndex, const double xSize, const double ySize, const double zSize, const double sphereSize, const double lc, const double pmlSizeE, const double pmlSizeN, const double pmlSizeW, const double pmlSizeS, const double pmlSizeD, const double pmlSizeU, bool transfinite);
+  void monoDomain(const double lc, const double pmlSize, const double R, const int order);
+  void checkerboard(const unsigned int nDomX, const unsigned int nDomY, const unsigned int nDomZ, const double xSize, const double ySize, const double zSize, const double sphereSize, const double lc, const bool pml, const double pmlSize, const bool pmlExt, const double pmlSizeExt, const int order, const unsigned int spherePosX = 0, const unsigned int spherePosY = 0, const unsigned int spherePosZ = 0, bool transfinite = false);
+  
+}
+
+#endif // H_MESH

examples/helmholtz/crossPoints/monoDomain.h

+#ifndef H_MONODOMAIN
+#define H_MONODOMAIN
+
+namespace D2 {
+
+  void monoDomain();
+
+}
+
+namespace D3 {
+
+  void monoDomain();
+
+}
+
+
+#endif // H_MONODOMAIN