Skip to content
Snippets Groups Projects
Select Git revision
  • a3221785cab590b387209068b407cbd5e61faf80
  • master default protected
  • dof-renumbering
  • test-dof-hash
  • gdemesy-master-patch-30528
  • eval-space-time
  • oscillating_multiharm
  • MH_movement
  • axisqu
  • write_vtu_and_ensight_formats
  • movingband
  • CP_1972_add_vtu_file_writing
  • mortar
  • fast_freq_sweep_Resolution
  • applyresolvent_again
  • marteaua-master-patch-54323
  • patch-1
  • binde-master-patch-08072
  • binde-master-patch-52461
  • BCGSL
  • resolvent
  • getdp_3_5_0
  • getdp_3_4_0
  • getdp_3_3_0
  • getdp_3_2_0
  • getdp_3_1_0
  • getdp_3_0_4
  • getdp_3_0_3
  • getdp_3_0_2
  • getdp_3_0_1
  • getdp_3_0_0
  • onelab_mobile_2.1.0
  • getdp_2_11_3 protected
  • getdp_2_11_2 protected
  • getdp_2_11_1 protected
  • getdp_2_11_0 protected
  • getdp_2_10_0 protected
  • getdp_2_9_2 protected
  • getdp_2_9_1 protected
  • getdp_2_9_0 protected
  • getdp_2_8_0 protected
41 results

F_PeWe.cpp

Blame
  • Formulation.cpp 73.28 KiB
    // GmshDDM - Copyright (C) 2019-2022, A. Royer, C. Geuzaine, Université de Liège
    //
    // See the LICENSE.txt file for license information. Please report all
    // issues on https://gitlab.onelab.info/gmsh/ddm/issues
    
    #include "Formulation.h"
    #include "Solvers.h"
    
    #include "GmshDdm.h"
    #include "MPIInterface.h"
    #include "Options.h"
    #include "gmshddmDefines.h"
    #include "resources.h"
    
    #include <gmshfem/CSVio.h>
    #include <gmshfem/Dof.h>
    #include <gmshfem/Message.h>
    #include <gmshfem/Vector.h>
    
    #include <regex>
    
    #ifdef HAVE_MPI
    #include <mpi.h>
    #endif
    
    #ifdef HAVE_PETSC
    #include <petsc.h>
    #endif
    
    namespace gmshddm
    {
    
    
      namespace problem
      {
    
        using gmshddm::common::s_printResources;
    
        template< class T_Scalar >
        Formulation< T_Scalar >::Formulation(const std::string &name, const std::vector< std::vector< unsigned int > > &topology) :
          _name(name), _volume(), _surface(topology.size()), _interfaceFields(), _interfaceMapping(), _rhs(), _numberOfIterations(0), _physicalCommutator(false), _artificialCommutator(false), _IA(false)
        {
          for(auto i = 0ULL; i < topology.size(); ++i) {
            gmshfem::msg::debug << "Create volume formulation (" << i << ")." << gmshfem::msg::endl;
            _volume.push_back(new gmshfem::problem::Formulation< T_Scalar >("volume_" + std::to_string(i)));
            for(auto ii = 0ULL; ii < topology[i].size(); ++ii) {
              const unsigned int j = topology[i][ii];
              gmshfem::msg::debug << "Create surface formulation (" << i << ", " << j << ")." << gmshfem::msg::endl;
              _surface[i].emplace(j, new gmshfem::problem::Formulation< T_Scalar >("surface_" + std::to_string(i) + "_" + std::to_string(j)));
            }
          }
        }
    
        template< class T_Scalar >
        Formulation< T_Scalar >::Formulation(const std::string &name, const std::vector< gmshfem::problem::Formulation< T_Scalar > * > &volume, const std::vector< std::unordered_map< unsigned int, gmshfem::problem::Formulation< T_Scalar > * > > &surface) :
          _name(name), _volume(volume), _surface(surface), _interfaceFields(), _interfaceMapping(), _rhs(), _numberOfIterations(0), _physicalCommutator(false), _artificialCommutator(false), _IA(false)
        {
        }
    
        template< class T_Scalar >
        Formulation< T_Scalar >::~Formulation()
        {
          for(auto i = 0ULL; i < _volume.size(); ++i) {
            delete _volume[i];
          }
          for(auto i = 0ULL; i < _surface.size(); ++i) {
            for(auto it = _surface[i].begin(); it != _surface[i].end(); ++it) {
              delete it->second;
            }
          }
        }
    
        template< class T_Scalar >
        typename Formulation< T_Scalar >::sourceStatus Formulation< T_Scalar >::indexStatusInShots(unsigned id) const {
        for (unsigned i = 0; i < _activeShots.size(); ++i) {
            if (_numberedPhysicalSourceTerms[i].count(id) > 0) {
                if (_activeShots[i]) {
                    return sourceStatus::ACTIVE;
                } else {
                    return sourceStatus::INACTIVE;
                }
            }
        }
        return sourceStatus::NOT_IN_NUMBERED_TERMS;
    }
    
        template< class T_Scalar >
        gmshfem::algebra::Vector< T_Scalar > Formulation< T_Scalar >::_sortDofValues(const gmshfem::algebra::Vector< T_Scalar > &vector, const unsigned int fieldTag)
        {
          auto it = _interfaceMapping.find(fieldTag);
          if(it == _interfaceMapping.end()) {
            gmshfem::msg::debug << "There is any mapping for the field tagged '" << fieldTag << "'" << gmshfem::msg::endl;
            return vector;
          }
          gmshfem::algebra::Vector< T_Scalar > newVector(vector.size());
    #pragma omp parallel for
          for(auto i = 0ULL; i < vector.size(); ++i) {
            newVector[i] = vector[it->second[i]];
          }
          return newVector;
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::_fillG(const T_Scalar *array)
        {
          const unsigned int MPI_Rank = mpi::getMPIRank();
          const unsigned int MPI_Size = mpi::getMPISize();
          if(MPI_Size == 1) {
            unsigned int pos = 0;
            for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
              for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
                for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
                  if((*itField->first).isDefined(idom, it->first)) {
                    const unsigned long long size = (*itField->first)(idom, it->first).numberOfUnknownDofs();
                    gmshfem::algebra::Vector< T_Scalar > g_ij(std::vector< T_Scalar >(&array[pos], &array[pos + size]));
                    (*itField->first)(idom, it->first).setUnknownVector(g_ij, gmshfem::dofs::RawOrder::Hash);
                    pos += size;
                  }
                }
              }
            }
          }
          else {
            std::vector< std::vector< std::vector< gmshfem::algebra::Vector< T_Scalar > * > > > g_ji;
            std::vector< MPI_Request > requests;
            std::vector< MPI_Status > status;
    
            // Send and recv g_ij
            for(auto recv : {true, false}) {
              unsigned int pos = 0;
              unsigned int currentSubdomain = 0;
              unsigned int currentInterface = 0;
              unsigned int currentField = 0;
              for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
                if(mpi::isItMySubdomain(idom)) {
                  if(recv) {
                    g_ji.push_back(std::vector< std::vector< gmshfem::algebra::Vector< T_Scalar > * > >());
                  }
                  for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
                    if(mpi::isItMySubdomain(it->first)) {
                      for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
                        if((*itField->first).isDefined(idom, it->first)) {
                          const unsigned long long size = (*itField->first)(idom, it->first).numberOfUnknownDofs();
                          if(!recv) {
                            gmshfem::algebra::Vector< T_Scalar > g_ij(std::vector< T_Scalar >(&array[pos], &array[pos + size]));
                            (*itField->first)(idom, it->first).setUnknownVector(g_ij, gmshfem::dofs::RawOrder::Hash);
                          }
                          pos += size;
                        }
                      }
                    }
                    else {
                      if(recv) {
                        g_ji[currentSubdomain].push_back(std::vector< gmshfem::algebra::Vector< T_Scalar > * >());
                      }
                      for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
                        if(!(*itField->first).isDefined(idom, it->first)) continue;
                        const unsigned int tagIJ = (*itField->first)(idom, it->first).tag();
                        const unsigned int tagJI = (*itField->second)(it->first, idom).tag();
                        const unsigned long long size = (*itField->first)(idom, it->first).numberOfUnknownDofs();
    
                        if(recv) {
                          g_ji[currentSubdomain][currentInterface].push_back(new gmshfem::algebra::Vector< T_Scalar >(size));
                          if(size) {
                            gmshfem::msg::debug << "MPI_Irecv g rank=" << MPI_Rank << " source=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagJI << " size=" << size << gmshfem::msg::endl;
                            MPI_Request req;
                            MPI_Irecv(&((*g_ji[currentSubdomain][currentInterface][currentField])[0]), size, MPI_C_DOUBLE_COMPLEX, mpi::rankOfSubdomain(it->first), tagJI, MPI_COMM_WORLD, &req);
                            requests.push_back(req);
                          }
                        }
                        else {
                          gmshfem::algebra::Vector< T_Scalar > g_ij(std::vector< T_Scalar >(&array[pos], &array[pos + size]));
                          (*itField->first)(idom, it->first).setUnknownVector(g_ij, gmshfem::dofs::RawOrder::Hash);
                          if(size) {
                            gmshfem::msg::debug << "MPI_Send g rank=" << MPI_Rank << " dest=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagIJ << " size=" << size << gmshfem::msg::endl;
                            MPI_Send(&g_ij[0], size, MPI_C_DOUBLE_COMPLEX, mpi::rankOfSubdomain(it->first), tagIJ, MPI_COMM_WORLD);
                          }
                        }
    
                        pos += size;
                        currentField += 1;
                      }
                      currentField = 0;
                      currentInterface += 1;
                    }
                  }
                  currentInterface = 0;
                  currentSubdomain += 1;
                }
              }
            }
            // wait until we receive the data
            gmshfem::msg::debug << "MPI_Wait g starting on rank=" << MPI_Rank << gmshfem::msg::endl;
            status.resize(requests.size());
            MPI_Waitall(requests.size(), &requests[0], &status[0]);
            for(auto s : status)
              if(s.MPI_ERROR != MPI_SUCCESS)
                gmshfem::msg::warning << "MPI_WaitAll g status " << s.MPI_ERROR  << gmshfem::msg::endl;
            gmshfem::msg::debug << "MPI_Wait g done on rank=" << MPI_Rank << gmshfem::msg::endl;
    
            MPI_Barrier(MPI_COMM_WORLD);
    
            {
              // assign g_ji
              unsigned int currentSubdomain = 0;
              unsigned int currentInterface = 0;
              unsigned int currentField = 0;
              for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
                if(mpi::isItMySubdomain(idom)) {
                  for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
                    if(mpi::isItMySubdomain(it->first)) continue;
                    for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
                      if(!(*itField->first).isDefined(idom, it->first)) continue;
                      const unsigned int tagIJ = (*itField->first)(idom, it->first).tag();
    
                      (*itField->second)(it->first, idom).setUnknownVector(_sortDofValues(*g_ji[currentSubdomain][currentInterface][currentField], tagIJ), gmshfem::dofs::RawOrder::Hash);
                      delete g_ji[currentSubdomain][currentInterface][currentField];
    
                      currentField += 1;
                    }
                    currentField = 0;
                    currentInterface += 1;
                  }
                  currentInterface = 0;
                  currentSubdomain += 1;
                }
              }
            }
    
            MPI_Barrier(MPI_COMM_WORLD);
          }
        }
    
        template< class T_Scalar >
        unsigned int Formulation< T_Scalar >::size() const
        {
          return _volume.size();
        }
    
        template< class T_Scalar >
        unsigned int Formulation< T_Scalar >::size(unsigned int i) const
        {
          return _surface.at(i).size();
        }
    
        template< class T_Scalar >
        const gmshfem::function::ScalarFunction< T_Scalar > Formulation< T_Scalar >::physicalSource(const gmshfem::function::ScalarFunction< T_Scalar > &f)
        {
          return commutator(f, &_physicalCommutator);
        }
    
        template< class T_Scalar >
        const gmshfem::function::VectorFunction< T_Scalar > Formulation< T_Scalar >::physicalSource(const gmshfem::function::VectorFunction< T_Scalar > &f)
        {
          return commutator(f, &_physicalCommutator);
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::physicalSourceTerm(const unsigned int termTag)
        {
          _physicalSourceTerms.insert(termTag);
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::numberedPhysicalSourceTerm(const unsigned int id, const unsigned int termTag)
        {
          _numberedPhysicalSourceTerms.at(id).insert(termTag);
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::compareRHSexport()
        {
          Vec RHS;
          Mat BlockRHS;
          _buildPetscRHS(&RHS);
          _buildPetscAllRHS(&BlockRHS);
    
          {
            PetscViewer viewer;
            PetscViewerASCIIOpen(PETSC_COMM_WORLD, "output_vec.txt", &viewer);
            PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_DENSE);
            VecView(RHS, viewer);
            PetscViewerDestroy(&viewer);
          }
    
          {
            PetscViewer viewer;
            PetscViewerASCIIOpen(PETSC_COMM_WORLD, "output_block.txt", &viewer);
            PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_DENSE);
            MatView(BlockRHS, viewer);
            PetscViewerDestroy(&viewer);
          }
        }
    
        template< class T_Scalar >
        const gmshfem::function::ScalarFunction< T_Scalar > Formulation< T_Scalar >::artificialSource(const gmshfem::function::ScalarFunction< T_Scalar > &f)
        {
          return commutator(f, &_artificialCommutator);
        }
    
        template< class T_Scalar >
        const gmshfem::function::VectorFunction< T_Scalar > Formulation< T_Scalar >::artificialSource(const gmshfem::function::VectorFunction< T_Scalar > &f)
        {
          return commutator(f, &_artificialCommutator);
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::artificialSourceTerm(const unsigned int termTag)
        {
          _artificialSourceTerms.insert(termTag);
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::togglePhysicalAndArtificialSourceTerms()
        {
          for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              for(auto it = _volume[idom]->begin(); it != _volume[idom]->end(); ++it) {
                unsigned termTag = (*it)->tag();
                auto termStatus = indexStatusInShots(termTag);
                if(_physicalSourceTerms.count(termTag) > 0 || termStatus != sourceStatus::NOT_IN_NUMBERED_TERMS) {
                  if(_physicalCommutator && termStatus != sourceStatus::INACTIVE) {
                    (*it)->activate();
                  }
                  else {
                    // Disable if all physical terms are disabled OR if it's a disabled shot
                    (*it)->deactivate();
                  }
                }
                else if(_artificialSourceTerms.count(termTag)) {
                  if(_artificialCommutator) {
                    (*it)->activate();
                  }
                  else {
                    (*it)->deactivate();
                  }
                }
              }
              for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
                for(auto itTerm = it->second->begin(); itTerm != it->second->end(); ++itTerm) {
                  if(_physicalSourceTerms.find((*itTerm)->tag()) != _physicalSourceTerms.end()) {
                    if(_physicalCommutator) {
                      (*itTerm)->activate();
                    }
                    else {
                      (*itTerm)->deactivate();
                    }
                  }
                  else if(_artificialSourceTerms.find((*itTerm)->tag()) != _artificialSourceTerms.end()) {
                    if(_artificialCommutator) {
                      (*itTerm)->activate();
                    }
                    else {
                      (*itTerm)->deactivate();
                    }
                  }
                }
              }
            }
          }
        }
    
        template< typename T_Scalar >
        void Formulation< T_Scalar >::setShotNumber(unsigned N)
        {
          _activeShots.resize(N, false);
          _numberedPhysicalSourceTerms.resize(N);
          _multi_rhs.resize(N);
        }
    
        template< typename T_Scalar >
        void Formulation< T_Scalar >::enableShot(unsigned i)
        {
          _activeShots.at(i) = true;
        }
    
        template< typename T_Scalar >
        void Formulation< T_Scalar >::disableShot(unsigned i)
        {
          _activeShots.at(i) = false;
        }
    
        template< typename T_Scalar >
        void Formulation< T_Scalar >::toggleShot(unsigned i)
        {
          _activeShots.at(i) = !_activeShots.at(i);
        }
    
        template< typename T_Scalar >
        bool Formulation< T_Scalar >::isShotEnabled(unsigned i) const
        {
          return _activeShots.at(i);
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::disableAllShots()
        {
          for (unsigned i = 0; i < _activeShots.size(); ++i)
            _activeShots.at(i) = false;
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::addInterfaceField(field::InterfaceFieldInterface< T_Scalar > &interfaceField)
        {
          _interfaceFields.insert(std::make_pair(&interfaceField, nullptr));
          gmshfem::msg::debug << "New InterfaceField named " << interfaceField.name() << " is added to the formulation '" << _name << "'." << gmshfem::msg::endl;
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::addInterfaceField(field::InterfaceFieldInterface< T_Scalar > &interfaceField, field::InterfaceFieldInterface< T_Scalar > &interfaceFieldMapping)
        {
          _interfaceFields.insert(std::make_pair(&interfaceField, &interfaceFieldMapping));
          gmshfem::msg::debug << "New InterfaceField named " << interfaceField.name() << " is added to the formulation '" << _name << "'." << gmshfem::msg::endl;
        }
    
        template< class T_Scalar >
        field::InterfaceFieldInterface< T_Scalar > *Formulation< T_Scalar >::getInterfaceField(const std::string &name) const
        {
          for(auto it = _interfaceFields.begin(); it != _interfaceFields.end(); ++it) {
            if(it->first->name() == name) {
              return it->first;
            }
          }
          return nullptr;
        }
    
        template< class T_Scalar >
        gmshfem::problem::Formulation< T_Scalar > &Formulation< T_Scalar >::operator()(unsigned int i)
        {
          return *_volume[i];
        }
    
        template< class T_Scalar >
        gmshfem::problem::Formulation< T_Scalar > &Formulation< T_Scalar >::operator()(unsigned int i, unsigned int j)
        {
          return *_surface[i].at(j);
        }
    
        template< class T_Scalar >
        gmshfem::common::Timer Formulation< T_Scalar >::pre(bool mustAssemble)
        {
          const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
          const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
    
          const unsigned int MPI_Rank = mpi::getMPIRank();
          const unsigned int MPI_Size = mpi::getMPISize();
    
          MPI_Barrier(MPI_COMM_WORLD);
    
          gmshfem::common::Timer time;
          time.tick();
    
          if(MPI_Rank == 0) {
            gmshfem::msg::info << "Pre-processing " << _name << "..." << gmshfem::msg::endl;
          }
          if(_interfaceFields.size() == 0) {
            gmshfem::msg::warning << "There is no interface field defined. Perhaps you have forgotten to declare them with 'addInterfaceField' function." << gmshfem::msg::endl;
          }
          s_printResources("DDM pre-processing begin", true);
    
          // Build interface field mapping
          if(MPI_Size != 1) {
            const unsigned int nDom = _volume.size();
            for(auto it = _interfaceFields.begin(); it != _interfaceFields.end(); ++it) {
              if(it->second == nullptr) {
                for(auto i = 0U; i < nDom; ++i) {
                  for(auto j = 0U; j < nDom; ++j) {
                    if(!it->first->isDefined(i, j)) continue;
    
                    for(auto it2 = _interfaceFields.begin(); it2 != _interfaceFields.end(); ++it2) {
                      if(!it2->first->isDefined(j, i)) continue;
    
                      if((*it2->first)(j, i).domain() == (*it->first)(i, j).domain()) {
                        it->second = it2->first;
                        gmshfem::msg::debug << "Interface field '" << it->first->name() << "' is mapped with interface field '" << it2->first->name() << "'" << gmshfem::msg::endl;
                        break;
                      }
                    }
                    if(it->second != nullptr) break;
                  }
                  if(it->second != nullptr) break;
                }
    
                if(it->second == nullptr) {
                  throw gmshfem::common::Exception("Unable to match the interface field '" + it->first->name() + "' with another interface field; give it manally");
                }
              }
            }
          }
    
          // ******************************************************
          // Initialize Kyrlov solver (compute b of (I-A) g^n = b)
          // ******************************************************
          _physicalCommutator = true;
          _artificialCommutator = false;
          togglePhysicalAndArtificialSourceTerms();
    
          // Pre-process the problem for all subdomains to get the dofs from the neighbours
          _preProVolume();
    
          if(MPI_Size == 1) {
            for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
              for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
                gmshfem::common::Options::instance()->verbose = innerVerbosity;
                it->second->pre();
                gmshfem::common::Options::instance()->verbose = outerVerbosity;
              }
            }
          }
          else {
            // Create MPI structure
            MPI_Aint struct_data[2] = {offsetof(gmshfem::dofs::RawDofKey, numType),
                                       offsetof(gmshfem::dofs::RawDofKey, entity)};
            int block_lengths[2] = {1, 1};
            MPI_Datatype types[2] = {MPI_UNSIGNED_LONG_LONG, MPI_UNSIGNED_LONG_LONG};
    
            // Commit the new datatype
            MPI_Datatype MPI_DOF;
            MPI_Type_create_struct(2, block_lengths, struct_data, types, &MPI_DOF);
            MPI_Type_commit(&MPI_DOF);
    
            // Send the RawDofKey vector
    
            // If
            //
            // 1) the mesh on the interface is identical on both subdomains (same
            //    entities, same elements and same nodes, with the same tags, and
            //    stored in the same order)
            //
            // 2) for Dofs associated with entities other than nodes (edges and
            //    faces), the edges and face tags correspond on both subdomains
            //    (e.g. they have been prealably generated globally in Gmsh)
            //
            // then exchanging the RawDofKey vector is actually enough, provided
            // that the Dofs are correctly sorted afterwards (see _interfaceMapping
            // below)
            //
            // Since generating the edge and face tags beforehand is cumbersome
            // (e.g. if the mesh was not generated with Gmsh), the code below also
            // exchanges the edge and face tags, based only on the assumption that
            // the meshes match (same entities, same elements and same nodes, with
            // the same tags, and stored in the same order).
    
            std::vector< MPI_Request > requests;
            std::vector< MPI_Status > status;
            std::vector< std::vector< std::vector< std::vector< gmshfem::dofs::RawDofKey > * > > > dofKeys;
    
            for(auto recv : {true, false}) {
              unsigned int currentSubdomain = 0;
              unsigned int currentInterface = 0;
              unsigned int currentField = 0;
              for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
                if(mpi::isItMySubdomain(idom)) {
                  if(recv) {
                    dofKeys.push_back(std::vector< std::vector< std::vector< gmshfem::dofs::RawDofKey > * > >());
                  }
                  for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
                    if(recv) {
                      gmshfem::common::Options::instance()->verbose = innerVerbosity;
                      it->second->pre();
                      gmshfem::common::Options::instance()->verbose = outerVerbosity;
                    }
                    if(mpi::isItMySubdomain(it->first)) continue;
                    if(recv) {
                      dofKeys[currentSubdomain].push_back(std::vector< std::vector< gmshfem::dofs::RawDofKey > * >());
                    }
                    for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
                      if(!(*itField->first).isDefined(idom, it->first)) continue;
                      if(recv) {
                        dofKeys[currentSubdomain][currentInterface].push_back(new std::vector< gmshfem::dofs::RawDofKey >());
                      }
                      const unsigned int tagIJ = (*itField->first)(idom, it->first).tag();
                      const unsigned int tagJI = (*itField->second)(it->first, idom).tag();
                      std::vector< gmshfem::dofs::RawDof > dofs;
                      // get dofs
                      (*itField->first)(idom, it->first).getUnknownDofVector(dofs);
                      if(recv) {
                        dofKeys[currentSubdomain][currentInterface][currentField]->resize(dofs.size());
                        if(dofs.size()) {
                          gmshfem::msg::debug << "MPI_Irecv dofs rank=" << MPI_Rank << " source=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagJI << " size=" << dofs.size() << gmshfem::msg::endl;
                          MPI_Request req;
                          MPI_Irecv(&(*dofKeys[currentSubdomain][currentInterface][currentField])[0], dofs.size(), MPI_DOF, mpi::rankOfSubdomain(it->first), tagJI, MPI_COMM_WORLD, &req);
                          requests.push_back(req);
                        }
                      }
                      else {
                        std::vector< gmshfem::dofs::RawDofKey > myDofKeys(dofs.size());
    #pragma omp parallel for
                        for(auto i = 0ULL; i < dofs.size(); ++i) {
                          myDofKeys[i].numType = dofs[i].numType;
                          myDofKeys[i].entity = dofs[i].entity;
                          dofs[i].numType -= GMSHFEM_DOF_FIELD_OFFSET * tagIJ;
                          dofs[i].numType += GMSHFEM_DOF_FIELD_OFFSET * tagJI;
                        }
                        // set dofs to the associated field
                        (*itField->second)(it->first, idom).setUnknownDofVector(dofs);
                        if(dofs.size()) {
                          gmshfem::msg::debug << "MPI_Send dofs rank=" << MPI_Rank << " dest=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagIJ << " size=" << dofs.size() << gmshfem::msg::endl;
                          MPI_Send(&myDofKeys[0], dofs.size(), MPI_DOF, mpi::rankOfSubdomain(it->first), tagIJ, MPI_COMM_WORLD);
                        }
                      }
    
                      currentField += 1;
                    }
                    currentField = 0;
                    currentInterface += 1;
                  }
                  currentInterface = 0;
                  currentSubdomain += 1;
                }
              }
            }
            // wait until we receive the data
            gmshfem::msg::debug << "MPI_Wait dofs starting on rank=" << MPI_Rank << gmshfem::msg::endl;
            status.resize(requests.size());
            MPI_Waitall(requests.size(), &requests[0], &status[0]);
            for(auto s : status)
              if(s.MPI_ERROR != MPI_SUCCESS)
                gmshfem::msg::warning << "MPI_WaitAll dofs status " << s.MPI_ERROR  << gmshfem::msg::endl;
            requests.clear();
            status.clear();
            gmshfem::msg::debug << "MPI_Wait dofs done on rank=" << MPI_Rank << gmshfem::msg::endl;
            MPI_Barrier(MPI_COMM_WORLD);
    
            // Generate the EdgeTags and FaceTags vectors
            std::vector< std::vector< std::vector< std::vector< unsigned long long > > > > edgeTags, edgeTagsNeighbor, faceTags, faceTagsNeighbor;
            {
              // This is only useful for Dofs associated to edges, in case a global
              // numbering of edges is not available. We might want to detect when
              // this communication step is actually necessary.
    
              unsigned int currentSubdomain = 0;
              unsigned int currentInterface = 0;
              unsigned int currentField = 0;
              for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
                if(mpi::isItMySubdomain(idom)) {
                  edgeTags.push_back(std::vector< std::vector< std::vector< unsigned long long > > >());
                  edgeTagsNeighbor.push_back(std::vector< std::vector< std::vector< unsigned long long > > >());
                  faceTags.push_back(std::vector< std::vector< std::vector< unsigned long long > > >());
                  faceTagsNeighbor.push_back(std::vector< std::vector< std::vector< unsigned long long > > >());
                  for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
                    if(mpi::isItMySubdomain(it->first)) continue;
                    edgeTags[currentSubdomain].push_back(std::vector< std::vector< unsigned long long > >());
                    edgeTagsNeighbor[currentSubdomain].push_back(std::vector< std::vector< unsigned long long > >());
                    faceTags[currentSubdomain].push_back(std::vector< std::vector< unsigned long long > >());
                    faceTagsNeighbor[currentSubdomain].push_back(std::vector< std::vector< unsigned long long > >());
                    for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
                      if(!(*itField->first).isDefined(idom, it->first)) continue;
                      edgeTags[currentSubdomain][currentInterface].push_back(std::vector< unsigned long long >());
                      edgeTagsNeighbor[currentSubdomain][currentInterface].push_back(std::vector< unsigned long long >());
                      faceTags[currentSubdomain][currentInterface].push_back(std::vector< unsigned long long >());
                      faceTagsNeighbor[currentSubdomain][currentInterface].push_back(std::vector< unsigned long long >());
    
                      // build the interface map
                      const gmshfem::domain::Domain domain = (*itField->first)(idom, it->first).domain();
                      for(auto itEntity = domain.cbegin(); itEntity != domain.cend(); ++itEntity) {
                        // edge case
                        if((*itField->first)(idom, it->first).getFunctionSpace()->getNumberOfKeysByEdge()) {
                          std::vector< int > elementTypes, orientations;
                          std::vector< std::size_t > edgeNodes, localEdgeTags;
                          gmsh::model::mesh::getElementTypes(elementTypes, itEntity->first, itEntity->second);
                          for(auto iElementType = 0ULL; iElementType < elementTypes.size(); ++iElementType) {
                            gmsh::model::mesh::getElementEdgeNodes(elementTypes[iElementType], edgeNodes, itEntity->second, true);
                            gmsh::model::mesh::getEdges(edgeNodes, localEdgeTags, orientations);
                            edgeTags[currentSubdomain][currentInterface][currentField].reserve(edgeTags[currentSubdomain][currentInterface][currentField].size() + localEdgeTags.size());
                            for(auto i = 0ULL; i < localEdgeTags.size(); ++i) {
                              edgeTags[currentSubdomain][currentInterface][currentField].push_back(localEdgeTags[i]);
                            }
                          }
                        }
                        // face (triangle) case
                        if((*itField->first)(idom, it->first).getFunctionSpace()->getNumberOfKeysByTriangularFace()) {
                          std::vector< int > elementTypes, orientations;
                          std::vector< std::size_t > faceNodes, localFaceTags;
                          gmsh::model::mesh::getElementTypes(elementTypes, itEntity->first, itEntity->second);
                          for(auto iElementType = 0ULL; iElementType < elementTypes.size(); ++iElementType) {
                            gmsh::model::mesh::getElementFaceNodes(elementTypes[iElementType], 3, faceNodes, itEntity->second, true);
                            gmsh::model::mesh::getFaces(3, faceNodes, localFaceTags, orientations);
                            faceTags[currentSubdomain][currentInterface][currentField].reserve(faceTags[currentSubdomain][currentInterface][currentField].size() + localFaceTags.size());
                            for(auto i = 0ULL; i < localFaceTags.size(); ++i) {
                              faceTags[currentSubdomain][currentInterface][currentField].push_back(localFaceTags[i]);
                            }
                          }
                        }
                        // face (quadrangle) case
                        if((*itField->first)(idom, it->first).getFunctionSpace()->getNumberOfKeysByQuadrangularFace()) {
                          std::vector< int > elementTypes, orientations;
                          std::vector< std::size_t > faceNodes, localFaceTags;
                          gmsh::model::mesh::getElementTypes(elementTypes, itEntity->first, itEntity->second);
                          for(auto iElementType = 0ULL; iElementType < elementTypes.size(); ++iElementType) {
                            gmsh::model::mesh::getElementFaceNodes(elementTypes[iElementType], 4, faceNodes, itEntity->second, true);
                            gmsh::model::mesh::getFaces(4, faceNodes, localFaceTags, orientations);
                            faceTags[currentSubdomain][currentInterface][currentField].reserve(faceTags[currentSubdomain][currentInterface][currentField].size() + localFaceTags.size());
                            for(auto i = 0ULL; i < localFaceTags.size(); ++i) {
                              faceTags[currentSubdomain][currentInterface][currentField].push_back(localFaceTags[i]);
                            }
                          }
                        }
                      }
                      edgeTagsNeighbor[currentSubdomain][currentInterface][currentField].resize(edgeTags[currentSubdomain][currentInterface][currentField].size());
                      faceTagsNeighbor[currentSubdomain][currentInterface][currentField].resize(faceTags[currentSubdomain][currentInterface][currentField].size());
    
                      currentField += 1;
                    }
                    currentField = 0;
                    currentInterface += 1;
                  }
                  currentInterface = 0;
                  currentSubdomain += 1;
                }
              }
            }
    
            // Send the EdgeTags vectors
            for(auto recv : {true, false}) {
              unsigned int currentSubdomain = 0;
              unsigned int currentInterface = 0;
              unsigned int currentField = 0;
              for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
                if(mpi::isItMySubdomain(idom)) {
                  for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
                    if(mpi::isItMySubdomain(it->first)) continue;
                    for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
                      if(!(*itField->first).isDefined(idom, it->first)) continue;
    
                      auto size = edgeTags[currentSubdomain][currentInterface][currentField].size();
                      if(size) {
                        const unsigned int tagIJ = (*itField->first)(idom, it->first).tag();
                        const unsigned int tagJI = (*itField->second)(it->first, idom).tag();
                        if(recv) {
                          gmshfem::msg::debug << "MPI_Irecv edges rank=" << MPI_Rank << " source=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagJI << " size=" << size << gmshfem::msg::endl;
                          MPI_Request req;
                          MPI_Irecv(&edgeTagsNeighbor[currentSubdomain][currentInterface][currentField][0], size, MPI_UNSIGNED_LONG_LONG, mpi::rankOfSubdomain(it->first), tagJI, MPI_COMM_WORLD, &req);
                          requests.push_back(req);
                        }
                        else {
                          gmshfem::msg::debug << "MPI_Send edges rank=" << MPI_Rank << " dest=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagIJ << " size=" << size << gmshfem::msg::endl;
                          MPI_Send(&edgeTags[currentSubdomain][currentInterface][currentField][0], size, MPI_UNSIGNED_LONG_LONG, mpi::rankOfSubdomain(it->first), tagIJ, MPI_COMM_WORLD);
                        }
                      }
    
                      currentField += 1;
                    }
                    currentField = 0;
                    currentInterface += 1;
                  }
                  currentInterface = 0;
                  currentSubdomain += 1;
                }
              }
            }
            // wait until we receive the data
            gmshfem::msg::debug << "MPI_Wait edges starting on rank=" << MPI_Rank << gmshfem::msg::endl;
            status.resize(requests.size());
            MPI_Waitall(requests.size(), &requests[0], &status[0]);
            for(auto s : status)
              if(s.MPI_ERROR != MPI_SUCCESS)
                gmshfem::msg::warning << "MPI_WaitAll edges status " << s.MPI_ERROR  << gmshfem::msg::endl;
            requests.clear();
            status.clear();
            gmshfem::msg::debug << "MPI_Wait edges done on rank=" << MPI_Rank << gmshfem::msg::endl;
            MPI_Barrier(MPI_COMM_WORLD);
    
            // Send the FaceTags vectors
            for(auto recv : {true, false}) {
              unsigned int currentSubdomain = 0;
              unsigned int currentInterface = 0;
              unsigned int currentField = 0;
              for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
                if(mpi::isItMySubdomain(idom)) {
                  for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
                    if(mpi::isItMySubdomain(it->first)) continue;
                    for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
                      if(!(*itField->first).isDefined(idom, it->first)) continue;
    
                      auto size = faceTags[currentSubdomain][currentInterface][currentField].size();
                      if(size) {
                        const unsigned int tagIJ = (*itField->first)(idom, it->first).tag();
                        const unsigned int tagJI = (*itField->second)(it->first, idom).tag();
                        if(recv) {
                          gmshfem::msg::debug << "MPI_Irecv faces rank=" << MPI_Rank << " source=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagJI << " size=" << size << gmshfem::msg::endl;
                          MPI_Request req;
                          MPI_Irecv(&faceTagsNeighbor[currentSubdomain][currentInterface][currentField][0], size, MPI_UNSIGNED_LONG_LONG, mpi::rankOfSubdomain(it->first), tagJI, MPI_COMM_WORLD, &req);
                          requests.push_back(req);
                        }
                        else {
                          gmshfem::msg::debug << "MPI_Send faces rank=" << MPI_Rank << " dest=" << mpi::rankOfSubdomain(it->first) << " tag=" << tagIJ << " size=" << size << gmshfem::msg::endl;
                          MPI_Send(&faceTags[currentSubdomain][currentInterface][currentField][0], size, MPI_UNSIGNED_LONG_LONG, mpi::rankOfSubdomain(it->first), tagIJ, MPI_COMM_WORLD);
                        }
                      }
    
                      currentField += 1;
                    }
                    currentField = 0;
                    currentInterface += 1;
                  }
                  currentInterface = 0;
                  currentSubdomain += 1;
                }
              }
            }
            // wait until we receive the data
            gmshfem::msg::debug << "MPI_Wait faces starting on rank=" << MPI_Rank << gmshfem::msg::endl;
            status.resize(requests.size());
            MPI_Waitall(requests.size(), &requests[0], &status[0]);
            for(auto s : status)
              if(s.MPI_ERROR != MPI_SUCCESS)
                gmshfem::msg::warning << "MPI_WaitAll faces status " << s.MPI_ERROR  << gmshfem::msg::endl;
            requests.clear();
            status.clear();
            gmshfem::msg::debug << "MPI_Wait faces done on rank=" << MPI_Rank << gmshfem::msg::endl;
            MPI_Barrier(MPI_COMM_WORLD);
    
            // Generate the interface mapping
            {
              unsigned int currentSubdomain = 0;
              unsigned int currentInterface = 0;
              unsigned int currentField = 0;
              for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
                if(mpi::isItMySubdomain(idom)) {
                  for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
                    if(mpi::isItMySubdomain(it->first)) continue;
                    for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
                      if(!(*itField->first).isDefined(idom, it->first)) continue;
                      const unsigned int tagIJ = (*itField->first)(idom, it->first).tag();
                      const unsigned int tagJI = (*itField->second)(it->first, idom).tag();
    
                      const unsigned long long numberOfEdges = edgeTags[currentSubdomain][currentInterface][currentField].size();
                      const unsigned long long numberOfFaces = faceTags[currentSubdomain][currentInterface][currentField].size();
    
                      // build maps
                      std::unordered_map< unsigned long long, unsigned long long > edgeToEdge(numberOfEdges);
                      std::unordered_map< unsigned long long, unsigned long long > faceToFace(numberOfFaces);
    
                      for(auto i = 0U; i < numberOfEdges; ++i) {
                        edgeToEdge.insert(std::make_pair(edgeTags[currentSubdomain][currentInterface][currentField][i], edgeTagsNeighbor[currentSubdomain][currentInterface][currentField][i]));
                      }
                      for(auto i = 0U; i < numberOfFaces; ++i) {
                        faceToFace.insert(std::make_pair(faceTags[currentSubdomain][currentInterface][currentField][i], faceTagsNeighbor[currentSubdomain][currentInterface][currentField][i]));
                      }
    
                      std::vector< gmshfem::dofs::RawDof > myDofs;
                      (*itField->second)(it->first, idom).getUnknownDofVector(myDofs);
    
                      _interfaceMapping.insert(std::make_pair(tagIJ, std::vector< unsigned long long >(myDofs.size())));
    
                      const unsigned int numberOfNodeKey = (*itField->first)(idom, it->first).getFunctionSpace()->getNumberOfKeysByNode();
                      const unsigned int numberOfEdgeKey = (*itField->first)(idom, it->first).getFunctionSpace()->getNumberOfKeysByEdge();
                      const unsigned int numberOfFaceKey = (*itField->first)(idom, it->first).getFunctionSpace()->getNumberOfKeysByTriangularFace() + (*itField->first)(idom, it->first).getFunctionSpace()->getNumberOfKeysByQuadrangularFace();
    #pragma omp parallel for
                      for(auto i = 0ULL; i < myDofs.size(); ++i) {
                        unsigned int firstKey = myDofs[i].numType - GMSHFEM_DOF_FIELD_OFFSET * tagJI;
                        if((*itField->first)(idom, it->first).multiplicity() > 1) { // To deal with compound fields
                          unsigned int multiplicity = 1;
                          while(firstKey >= GMSHFEM_DOF_FIELD_OFFSET || multiplicity < (*itField->first)(idom, it->first).multiplicity()) {
                            firstKey -= GMSHFEM_DOF_FIELD_OFFSET;
                            ++multiplicity;
                          }
                        }
                        if(numberOfNodeKey == 0) { // 1-Form as the firstKey for edges always starts with value '1', even if there is no node keys.
                          firstKey--;
                        }
    
                        if(firstKey < numberOfNodeKey) { // if the dof is assocated to a node
                          // nothing to do (note that we might want to handle meshes with different node tags on both sides)
                        }
                        else if(firstKey >= numberOfNodeKey && firstKey < numberOfNodeKey + numberOfEdgeKey) { // if the dof is assocated to a edge
                          myDofs[i].entity = edgeToEdge[myDofs[i].entity];
                        }
                        else if(firstKey >= numberOfNodeKey + numberOfEdgeKey && firstKey < numberOfNodeKey + numberOfEdgeKey + numberOfFaceKey) { // if the dof is assocated to a face
                          myDofs[i].entity = faceToFace[myDofs[i].entity];
                        }
                        else {
                          gmshfem::msg::error << "An unknown key is found (firstKey = " << firstKey << ")" << gmshfem::msg::endl;
                        }
    
                        auto itFind = std::find_if(dofKeys[currentSubdomain][currentInterface][currentField]->begin(), dofKeys[currentSubdomain][currentInterface][currentField]->end(), [&](const gmshfem::dofs::RawDofKey &other) { return myDofs[i].numType == other.numType && myDofs[i].entity == other.entity; });
                        if(itFind == dofKeys[currentSubdomain][currentInterface][currentField]->end()) {
                          gmshfem::msg::error << "Unable to find the corresponding interface" << gmshfem::msg::endl;
                        }
                        _interfaceMapping[tagIJ][i] = std::distance(dofKeys[currentSubdomain][currentInterface][currentField]->begin(), itFind);
                      }
                      delete dofKeys[currentSubdomain][currentInterface][currentField];
    
                      currentField += 1;
                    }
                    currentField = 0;
                    currentInterface += 1;
                  }
                  currentInterface = 0;
                  currentSubdomain += 1;
                }
              }
            }
    
            // Free the MPI_DOF datatype
            MPI_Type_free(&MPI_DOF);
          }
    
          MPI_Barrier(MPI_COMM_WORLD);
    
          _infoNumberOfDofs();
    
          if(mustAssemble) {
    
            assemble();
            time.tock();
            if(MPI_Rank == 0) {
              gmshfem::msg::info << "Done pre-processing and assembling in " << time << "s" << gmshfem::msg::endl;
            }
    
            return time;
          }
          else {
            time.tock();
            if(MPI_Rank == 0) {
              gmshfem::msg::info << "Done pre-processing (no assembly) in " << time << "s" << gmshfem::msg::endl;
            }
            return time;
          }
        }
    
    
        template< class T_Scalar >
        gmshfem::common::Timer Formulation< T_Scalar >::assemble(int storeID)
        {
          gmshfem::common::Timer time;
          time.tick();
    
          const unsigned int MPI_Rank = mpi::getMPIRank();
    
          _physicalCommutator = true;
          _artificialCommutator = false;
          togglePhysicalAndArtificialSourceTerms();
    
          // Clear volume
          for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              // Clear formulation
              _volume[idom]->setSystemToZero();
            }
          }
    
          for(auto idom = 0ULL; idom < _surface.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
                gmshfem::msg::info << " - Clearing " << it->second->name() << gmshfem::msg::endl;
                it->second->setSystemToZero();
              }
            }
          }
    
          // Assemble locally on each processor
          s_printResources("DDM pre-processing before volume assembly");
          _assembleAllVolume();
    
          MPI_Barrier(MPI_COMM_WORLD);
          s_printResources("DDM pre-processing before volume solve");
    
          // solve locally on each processor
          _solveAllVolume();
    
          MPI_Barrier(MPI_COMM_WORLD);
          s_printResources("DDM pre-processing before surface assembly and solve");
    
          _assembleAndSolveSurface();
    
          MPI_Barrier(MPI_COMM_WORLD);
    
          // get the local interface unknowns for each processor
          _extractRHS();
          if (storeID > -1) {
            if (storeID >= _activeShots.size()) {
              throw gmshfem::common::Exception("Invalid storeID argument to gmshddm::problem::Formulation::assemble()");
            }
            _multi_rhs.at(storeID) = _rhs; // Do a copy. Make moving an option ?
          }
    
          MPI_Barrier(MPI_COMM_WORLD);
    
          s_printResources("DDM pre-processing end");
    
    
          time.tock();
          if(MPI_Rank == 0) {
            gmshfem::msg::info << "Done assembling in " << time << "s" << gmshfem::msg::endl;
          }
    
          return time;
        }
    
        template< class T_Scalar >
        gmshfem::common::Timer Formulation< T_Scalar >::solve(const std::string& solver, const double tolerance, const int iterMax, const bool sameMatrixWithArtificialAndPhysicalSources, const bool skipFinalSolutionComputation)
        {
          std::unique_ptr < AbstractIterativeSolver< T_Scalar >> solverObject;
          // For backwards compatibility.
          if(solver == "jacobi") {
            solverObject = std::make_unique< JacobiIterativeSolver< T_Scalar > >();
          }
          else {
            solverObject = std::make_unique< KSPIterativeSolver< T_Scalar > >(solver);
          }
    
          return solve(*solverObject, tolerance, iterMax, sameMatrixWithArtificialAndPhysicalSources, skipFinalSolutionComputation);
        }
    
        template< class T_Scalar >
        gmshfem::common::Timer Formulation< T_Scalar >::solve(AbstractIterativeSolver<T_Scalar>& solver, const double tolerance, const int iterMax, const bool sameMatrixWithArtificialAndPhysicalSources, const bool skipFinalSolutionComputation)
        {
          // sameMatrixWithArtificialAndPhysicalSources could be automatically set
          // to true if bilinear terms of volume formulations do not use the
          // physical/artifical commutator and if there are no non-homogeneous
          // Dirichlet BCs
    
          const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
          const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
    
          const unsigned int MPI_Rank = mpi::getMPIRank();
    
          MPI_Barrier(MPI_COMM_WORLD);
    
          gmshfem::common::Timer time;
          time.tick();
    
          if(MPI_Rank == 0) {
            gmshfem::msg::info << "Solving " << _name << " using " << solver.solverName() << "..." << gmshfem::msg::endl;
          }
          s_printResources("DDM iterative solver begin");
    
          // ******************************************************
          // Iterative solver
          // ******************************************************
    
          // PETSc
          Mat A;
          Vec X, Y, RHS;
    
          // local size per proc and global size
          PetscInt locSize = _rhs.size();
          PetscInt globSize;
    
          // Create Petsc distributed vector - let Petsc compute the global vector size
          _buildPetscRHS(&RHS);
          VecGetSize(RHS, &globSize);
    
          VecDuplicate(RHS, &X); // current iterate
          VecDuplicate(RHS, &Y); // new iterate
    
          MatCreateShell(MPI_COMM_WORLD, locSize, locSize, globSize, globSize, this, &A);
          switch(solver.productType()) {
              case ProductType::A:
                MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) & MatVectProductA< T_Scalar >);
                break;
              case ProductType::I_A:
                MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) & MatVectProductI_A< T_Scalar >);
                break;
              default:
              throw gmshfem::common::Exception("Unknown product type.");
          }
    
          _physicalCommutator = false;
          _artificialCommutator = true;
          togglePhysicalAndArtificialSourceTerms();
    
    
          // Assemble the volumic systems for each subdomain
          for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              _volume[idom]->setAttribute("ddm::physicalCommutator", _physicalCommutator);
              _volume[idom]->setAttribute("ddm::artificialCommutator", _artificialCommutator);
              gmshfem::common::Options::instance()->verbose = innerVerbosity;
              if(sameMatrixWithArtificialAndPhysicalSources) {
                _deactivateBilinear(idom);
                _volume[idom]->setRHSToZero();
                _volume[idom]->assemble();
              }
              else {
                _volume[idom]->removeSystem();
                _volume[idom]->initSystem();
                _volume[idom]->pre();
                _volume[idom]->assemble();
                _deactivateBilinear(idom);
              }
              gmshfem::common::Options::instance()->verbose = outerVerbosity;
            }
          }
    
    
          // Call the AbstractIterativeSolver object
          solver.solve(A, X, Y, RHS, tolerance, iterMax);
          _numberOfIterations = solver.absoluteResidual().size();
    
    
          // ******************************************************
          // Final solution
          // ******************************************************
          if(skipFinalSolutionComputation) {
            time.tock();
            return time;
          }
    
          _physicalCommutator = true;
          _artificialCommutator = true;
          togglePhysicalAndArtificialSourceTerms();
    
          for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              _volume[idom]->setAttribute("ddm::physicalCommutator", _physicalCommutator);
              _volume[idom]->setAttribute("ddm::artificialCommutator", _artificialCommutator);
              gmshfem::common::Options::instance()->verbose = innerVerbosity;
              if(sameMatrixWithArtificialAndPhysicalSources) {
                _volume[idom]->setRHSToZero();
              }
              else {
                _volume[idom]->removeSystem();
                _volume[idom]->initSystem();
                _volume[idom]->pre();
              }
              gmshfem::common::Options::instance()->verbose = outerVerbosity;
            }
          }
    
          _fromPetscToInterfaceFields(Y);
    
          // SolveVolumePDE:
          for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              gmshfem::common::Options::instance()->verbose = innerVerbosity;
              if(!sameMatrixWithArtificialAndPhysicalSources) {
                _activateBilinear(idom);
              }
              _volume[idom]->assemble();
              _volume[idom]->solve(true);
              if(sameMatrixWithArtificialAndPhysicalSources) {
                _activateBilinear(idom);
              }
              gmshfem::common::Options::instance()->verbose = outerVerbosity;
            }
          }
    
          MatDestroy(&A);
          VecDestroy(&X);
          VecDestroy(&Y);
          VecDestroy(&RHS);
    
          MPI_Barrier(MPI_COMM_WORLD);
          s_printResources("DDM iterative solver end");
    
          time.tock();
          if(MPI_Rank == 0) {
            gmshfem::msg::info << "Done solving in " << time << "s" << gmshfem::msg::endl;
          }
    
          return time;
        }
    
        template< class T_Scalar >
        gmshfem::algebra::MatrixCCS< T_Scalar > Formulation< T_Scalar >::computeMatrix()
        {
          gmshfem::common::Timer time;
          time.tick();
    
          const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
          const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
    
          Mat A;
          Vec X, Y, Y_seq;
          VecScatter     ctx;
    
          // Create a distributed vector that acts as a RHS
          PetscInt locSize = _rhs.size();
          VecCreateMPI(MPI_COMM_WORLD, locSize, PETSC_DETERMINE, &Y);
          VecDuplicate(Y, &X);
          
          // Create a sequential full vector and a scatter object to synchronize the global vector of outputs,
          // Needed to build the matrix
          PetscInt globSize;
          VecGetSize(Y, &globSize);
          VecCreateSeq(MPI_COMM_SELF, globSize, &Y_seq);
          VecScatterCreateToAll(Y, &ctx, &Y_seq);
    
          if (mpi::getMPIRank() == 0)
            gmshfem::msg::info << "Local size:" << locSize << " and global size:" << globSize << "." << gmshfem::msg::endl;
    
    
          // Setup the A matrix
          MatCreateShell(MPI_COMM_WORLD, locSize, locSize, globSize, globSize, this, &A);
          MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) & MatVectProductI_A< T_Scalar >);
    
          _physicalCommutator = false;
          _artificialCommutator = true;
          togglePhysicalAndArtificialSourceTerms();
    
      
          // Compute sub matrices then disable bilinear terms (to assemble only RHS)
          for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
            _volume[idom]->setAttribute("ddm::physicalCommutator", _physicalCommutator);
            _volume[idom]->setAttribute("ddm::artificialCommutator", _artificialCommutator);
            gmshfem::common::Options::instance()->verbose = innerVerbosity;
            _volume[idom]->removeSystem();
            _volume[idom]->initSystem();
            _volume[idom]->pre();
            _volume[idom]->assemble();
            _deactivateBilinear(idom);
            gmshfem::common::Options::instance()->verbose = outerVerbosity;
          }
    
    
          gmshfem::algebra::Vector< T_Scalar > vec(globSize);
          std::vector< PetscInt > aiVector(globSize);
          for(auto i = 0ULL; i < globSize; ++i) {
            aiVector[i] = i;
          }
          std::vector< PetscScalar > aVector(globSize);
    
          std::vector< unsigned long long > aj(globSize + 1, 0);
          std::vector< unsigned long long > ai;
          ai.reserve(globSize);
          std::vector< T_Scalar > a;
          a.reserve(globSize);
    
          std::vector<PetscScalar> zeros(globSize, 0);
          // Init RHS vector
          VecSet(X, 0.0);
          VecSetValue(X, 0, 1.0, INSERT_VALUES);
          VecAssemblyBegin(X);
          VecAssemblyEnd(X);
    
          for(auto i = 0ULL; i < globSize; ++i) {
    
            VecSet(X, 0.0);
            VecSetValue(X, i, 1.0, INSERT_VALUES);
            VecAssemblyBegin(X);
            VecAssemblyEnd(X);
    
            if(mpi::getMPIRank() == 0) {
             
              gmshfem::msg::info << "Column number " << i << " over " << globSize << "." << gmshfem::msg::endl;
            }
            MatMult(A, X, Y);
    
            // Scatter to rank 0
            VecScatterBegin(ctx, Y, Y_seq, INSERT_VALUES, SCATTER_FORWARD);
            VecScatterEnd(ctx, Y, Y_seq, INSERT_VALUES, SCATTER_FORWARD);
    
            VecGetValues(Y_seq, globSize, &aiVector[0], &aVector[0]);
    
    
    
            aj[i + 1] = aj[i];
            for(auto j = 0ULL; j < aVector.size(); ++j) {
              if(aVector[j] != 0) {
                aj[i + 1]++;
                ai.push_back(j);
                a.push_back(T_Scalar(aVector[j]));
              }
              
            }
    
            ai.reserve(aj.size() + globSize);
            a.reserve(a.size() + globSize);
    
            // Update RHS
            VecSetValue(X, i, 0.0, INSERT_VALUES);
            if(i + 1 < globSize)
              VecSetValue(X, i + 1, 1.0, INSERT_VALUES);
            VecAssemblyBegin(X);
            VecAssemblyEnd(X);
          }
    
          for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
            _activateBilinear(idom);
          }
    
          VecScatterDestroy(&ctx);
    
    
          gmshfem::algebra::MatrixCCS< T_Scalar > mat(std::move(ai), std::move(aj), std::move(a));
    
          time.tock();
          if(mpi::getMPIRank() == 0) {
            gmshfem::msg::info << "Computed the matrix in " << time.time() << "s." << gmshfem::msg::endl;
          }
          return mat;
        }
    
        template< class T_Scalar >
        gmshfem::algebra::MatrixCCS< T_Scalar > Formulation< T_Scalar >::computeMatrixBlock()
        {
          gmshfem::common::Timer time;
          time.tick();
    
          const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
          const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
    
          Mat A;
          Vec Y, Y_seq;
          VecScatter     ctx;
    
          // Create a distributed vector that acts as a RHS
          PetscInt locSize = _rhs.size();
          VecCreateMPI(MPI_COMM_WORLD, locSize, PETSC_DETERMINE, &Y);
          VecSetUp(Y);
          
          // Create a sequential full vector and a scatter object to synchronize the global vector of outputs,
          // Needed to build the matrix
          PetscInt globSize;
          VecGetSize(Y, &globSize);
    
          // Create the identity matrix
          Mat ID;
          {
            MatCreate(PETSC_COMM_WORLD, &ID);
            MatSetSizes(ID, locSize, PETSC_DECIDE, globSize, globSize);
            MatSetType(ID, MATMPIDENSE);
            MatSetUp(ID);
            MatZeroEntries(ID);
    
            PetscInt start, end;
            MatGetOwnershipRange(ID, &start, &end);
    
            for(PetscInt i = start; i < end; ++i) {
              PetscScalar val = 1.0;
              MatSetValues(ID, 1, &i, 1, &i, &val, INSERT_VALUES);
            }
    
            MatAssemblyBegin(ID, MAT_FINAL_ASSEMBLY);
            MatAssemblyEnd(ID, MAT_FINAL_ASSEMBLY);
          }
    
          VecCreateSeq(MPI_COMM_SELF, globSize, &Y_seq);
          VecScatterCreateToAll(Y, &ctx, &Y_seq);
    
          if (mpi::getMPIRank() == 0)
            gmshfem::msg::info << "Local size:" << locSize << " and global size:" << globSize << "." << gmshfem::msg::endl;
    
    
          // Setup the A matrix
          MatCreateShell(MPI_COMM_WORLD, locSize, locSize, globSize, globSize, this, &A);
          MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) & MatVectProductI_A< T_Scalar >);
    
          _physicalCommutator = false;
          _artificialCommutator = true;
          togglePhysicalAndArtificialSourceTerms();
    
      
          // Compute sub matrices then disable bilinear terms (to assemble only RHS)
          for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
            _volume[idom]->setAttribute("ddm::physicalCommutator", _physicalCommutator);
            _volume[idom]->setAttribute("ddm::artificialCommutator", _artificialCommutator);
            gmshfem::common::Options::instance()->verbose = innerVerbosity;
            _volume[idom]->removeSystem();
            _volume[idom]->initSystem();
            _volume[idom]->pre();
            _volume[idom]->assemble();
            _deactivateBilinear(idom);
            gmshfem::common::Options::instance()->verbose = outerVerbosity;
          }
    
    
          Mat Yblock;
          MatDuplicate(ID, MAT_DO_NOT_COPY_VALUES, &Yblock);
          MatMatProductI_A< T_Scalar >(A, ID, Yblock);
    
          std::vector< PetscInt > aiVector(globSize);
          std::vector< PetscScalar > aVector(globSize);
    
          std::vector< unsigned long long > aj(globSize + 1, 0);
          std::vector< unsigned long long > ai;
          ai.reserve(globSize);
          std::vector< T_Scalar > a;
          a.reserve(globSize);
    
          for(auto i = 0ULL; i < globSize; ++i) {
            aiVector[i] = i;
          }
    
          for(auto i = 0ULL; i < globSize; ++i) {
    
            
            if(mpi::getMPIRank() == 0) {
    
              gmshfem::msg::info << "Column number " << i << " over " << globSize << "." << gmshfem::msg::endl;
            }
    
            // Scatter to rank 0
            MatGetColumnVector(Yblock, Y, i);
            VecScatterBegin(ctx, Y, Y_seq, INSERT_VALUES, SCATTER_FORWARD);
            VecScatterEnd(ctx, Y, Y_seq, INSERT_VALUES, SCATTER_FORWARD);
            VecGetValues(Y_seq, globSize, &aiVector[0], &aVector[0]);
    
           
            aj[i + 1] = aj[i];
            for(auto j = 0ULL; j < aVector.size(); ++j) {
              if(aVector[j] != 0) {
                aj[i + 1]++;
                ai.push_back(j);
                a.push_back(T_Scalar(aVector[j]));
              }
            }
    
            ai.reserve(aj.size() + globSize);
            a.reserve(a.size() + globSize);
    
          }
          
    
          for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
            _activateBilinear(idom);
          }
    
          VecScatterDestroy(&ctx);
          gmshfem::algebra::MatrixCCS< T_Scalar > mat(std::move(ai), std::move(aj), std::move(a));
    
          time.tock();
          if(mpi::getMPIRank() == 0) {
            gmshfem::msg::info << "Computed the matrix in " << time.time() << "s." << gmshfem::msg::endl;
          }
          return mat;
    
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::setSolutionIntoInterfaceFields(const std::vector< T_Scalar > &solution)
        {
          _fillG(solution.data());
        }
    
        template< class T_Scalar >
        const gmshfem::algebra::Vector< T_Scalar > &Formulation< T_Scalar >::getRHS() const
        {
          return _rhs;
        }
    
        template< class T_Scalar >
        unsigned int Formulation< T_Scalar >::numberOfIterations() const
        {
          return _numberOfIterations;
        }
    
        template< class T_Scalar >
        int MatVectProductA(Mat A, Vec X, Vec Y)
        {
          return MatVectProductImpl<T_Scalar>(A, X, Y, false);
        }
    
        template< class T_Scalar >
        int MatVectProductI_A(Mat A, Vec X, Vec Y)
        {
          return MatVectProductImpl<T_Scalar>(A, X, Y, true);
        }
    
        template< class T_Scalar >
        int MatMatProductI_A(Mat A, Mat X, Mat Y)
        {
          return MatMatProductImpl<T_Scalar>(A, X, Y, true);
        }
    
        template< class T_Scalar >
        int MatMatProductA(Mat A, Mat X, Mat Y)
        {
          return MatMatProductImpl<T_Scalar>(A, X, Y, false);
        }
    
        template< class T_Scalar >
        void Formulation<T_Scalar>::_fromPetscToInterfaceFields(Vec G)
        {
          const PetscScalar *petscArray;
          PetscInt petscArraySize;
          VecGetArrayRead(G, &petscArray);
          VecGetLocalSize(G, &petscArraySize);
          const T_Scalar *array = PetscInterface< T_Scalar, PetscScalar, PetscInt >::arrayInterface(petscArray, petscArraySize);
          _fillG(array);
          PetscInterface< T_Scalar, PetscScalar, PetscInt >::freeArray(array);
          VecRestoreArrayRead(G, &petscArray);
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::_preProVolume()
        {
          const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
          const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
    
          for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              gmshfem::common::Options::instance()->verbose = innerVerbosity;
              _volume[idom]->pre();
              gmshfem::common::Options::instance()->verbose = outerVerbosity;
            }
          }
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::_infoNumberOfDofs() const
        {
          for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              gmshfem::msg::info << " - Subdomain " << idom << ":" << gmshfem::msg::endl;
              gmshfem::msg::info << "  * Subdomain size: " << _volume[idom]->getNumberOfUnknownDof() << " dofs" << gmshfem::msg::endl;
              unsigned long long totalInterfaceSize = 0;
              for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
                totalInterfaceSize += it->second->getNumberOfUnknownDof();
              }
              gmshfem::msg::info << "  * Interface size: " << totalInterfaceSize << " dofs:" << gmshfem::msg::endl;
              for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
                gmshfem::msg::info << "   · Interface 'g(" << idom << ", " << it->first << ")' size: " << it->second->getNumberOfUnknownDof() << " dofs" << gmshfem::msg::endl;
              }
            }
          }
        }
    
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::_assembleAllVolume()
        {
          const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
          const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
    
          for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              gmshfem::msg::info << " - Assembling " << _volume[idom]->name() << gmshfem::msg::endl;
              _volume[idom]->setAttribute("ddm::physicalCommutator", _physicalCommutator);
              _volume[idom]->setAttribute("ddm::artificialCommutator", _artificialCommutator);
              gmshfem::common::Options::instance()->verbose = innerVerbosity;
              _volume[idom]->assemble();
              gmshfem::common::Options::instance()->verbose = outerVerbosity;
              if(gmshddm::common::Options::instance()->memory) {
                gmshfem::msg::info << " - Subdomain " << idom << ": number of non zeros : " << _volume[idom]->getNumberOfNonZeros() << gmshfem::msg::endl;
              }
            }
          }
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::_solveAllVolume()
        {
          const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
          const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
    
          for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              gmshfem::msg::info << " - Solving " << _volume[idom]->name() << gmshfem::msg::endl;
              if(gmshddm::common::Options::instance()->memory) {
                gmshfem::msg::info << " - Subdomain " << idom << ": estimated factorization memory : " << _volume[idom]->getEstimatedFactorizationMemoryUsage() << gmshfem::msg::endl;
              }
              gmshfem::common::Options::instance()->verbose = innerVerbosity;
              _volume[idom]->solve();
              gmshfem::common::Options::instance()->verbose = outerVerbosity;
            }
          }
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::_assembleAndSolveSurface()
        {
          const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
          const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
    
          for(auto idom = 0ULL; idom < _surface.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
                gmshfem::msg::info << " - Assembling and solving " << it->second->name() << gmshfem::msg::endl;
                gmshfem::common::Options::instance()->verbose = innerVerbosity;
                it->second->assemble();
                it->second->solve();
                gmshfem::common::Options::instance()->verbose = outerVerbosity;
              }
            }
          }
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::_extractRHS()
        {
          _rhs.clear();
          for(auto idom = 0ULL; idom < _surface.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
                for(auto itField = _interfaceFields.begin(); itField != _interfaceFields.end(); ++itField) {
                  if((*itField->first).isDefined(idom, it->first)) {
                    gmshfem::algebra::Vector< T_Scalar > g_ij;
                    (*itField->first)(idom, it->first).getUnknownVector(g_ij, gmshfem::dofs::RawOrder::Hash);
                    // concatenate the local unknowns g_ij of the current process and set it in the rhs for the iterative solver
                    _rhs.concatenate(g_ij);
                  }
                }
              }
            }
          }
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::_activateBilinear(const unsigned long long idom)
        {
          for(auto formulation : *_volume[idom]) {
            if(formulation->isBilinear()) {
              formulation->activate();
            }
          }
          for(auto formulation : _surface[idom]) {
            for(auto term : *formulation.second) {
              if(term->isBilinear()) {
                term->activate();
              }
            }
          }
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::_deactivateBilinear(const unsigned long long idom)
        {
          for(auto formulation : *_volume[idom]) {
            if(formulation->isBilinear()) {
              formulation->deactivate();
            }
          }
          for(auto formulation : _surface[idom]) {
            for(auto term : *formulation.second) {
              if(term->isBilinear()) {
                term->deactivate();
              }
            }
          }
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::_buildPetscRHS(Vec *RHS)
        {
    
          auto locSize = _rhs.size();
          VecCreateMPI(MPI_COMM_WORLD, locSize, PETSC_DETERMINE, RHS);
          VecCopy(_rhs.getPetsc(), *RHS);
        }
    
        template< class T_Scalar >
        void Formulation< T_Scalar >::_buildPetscAllRHS(Mat *RHS)
        {
          gmshfem::common::Timer time;
          time.tick();
          if(mpi::getMPIRank() == 0) {
            gmshfem::msg::info << "Starting to build a block of RHS." << gmshfem::msg::endl;
          }
    
          auto locSize = _rhs.size();
          auto nRHS = _activeShots.size();
    
          MatCreateDense(MPI_COMM_WORLD, locSize, PETSC_DETERMINE, PETSC_DETERMINE, nRHS, NULL, RHS);
          MatSetUp(*RHS);
    
          std::vector< PetscInt > rowIndices(locSize);
          PetscInt startRow, endRow;
          MatGetOwnershipRange(*RHS, &startRow, &endRow);
    
          for(PetscInt i = 0; i < locSize; i++) {
            rowIndices[i] = startRow + i;
          }
    
          // For each RHS, copy the vector into the column of the block
          for(PetscInt j = 0; j < nRHS; j++) {
            Vec bj = _multi_rhs[j].getPetsc();
    
            PetscScalar *colValues;
            VecGetArray(bj, &colValues);
            MatSetValues(*RHS, locSize, rowIndices.data(), 1, &j, colValues, INSERT_VALUES);
            VecRestoreArray(bj, &colValues);
          }
    
          MatAssemblyBegin(*RHS, MAT_FINAL_ASSEMBLY);
          MatAssemblyEnd(*RHS, MAT_FINAL_ASSEMBLY);
    
    
          time.tock();
          if(mpi::getMPIRank() == 0) {
            gmshfem::msg::info << "Done assembling block of RHS in " << time << "s" << gmshfem::msg::endl;
          }
        }
    
    
        template< class T_Scalar >
        int MatVectProductImpl(Mat A, Vec X, Vec Y, bool IA)
        {
          PetscFunctionBegin;
    
          const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
          const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
    
          Formulation< T_Scalar > *formulation;
          MatShellGetContext(A, (void **)&formulation);
    
          formulation->_fromPetscToInterfaceFields(X);
    
          // SolveVolumePDE
          for(auto idom = 0ULL; idom < formulation->_volume.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              gmshfem::common::Options::instance()->verbose = innerVerbosity;
              formulation->_volume[idom]->setRHSToZero();
              formulation->_volume[idom]->assemble();
              gmshfem::common::Options::instance()->verbose = outerVerbosity;
            }
          }
    
          for(auto idom = 0ULL; idom < formulation->_volume.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              gmshfem::common::Options::instance()->verbose = innerVerbosity;
              formulation->_volume[idom]->solve(true);
              gmshfem::common::Options::instance()->verbose = outerVerbosity;
            }
          }
    
          // SolveSurfacePDE
          for(auto idom = 0ULL; idom < formulation->_surface.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              for(auto it = formulation->_surface[idom].begin(); it != formulation->_surface[idom].end(); ++it) {
                gmshfem::common::Options::instance()->verbose = innerVerbosity;
                it->second->setRHSToZero();
                it->second->assemble();
                gmshfem::common::Options::instance()->verbose = outerVerbosity;
              }
            }
          }
    
          for(auto idom = 0ULL; idom < formulation->_surface.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              for(auto it = formulation->_surface[idom].begin(); it != formulation->_surface[idom].end(); ++it) {
                gmshfem::common::Options::instance()->verbose = innerVerbosity;
                it->second->solve(true);
                gmshfem::common::Options::instance()->verbose = outerVerbosity;
              }
            }
          }
    
          MPI_Barrier(MPI_COMM_WORLD);
    
          gmshfem::algebra::Vector< T_Scalar > g_local;
          for(auto idom = 0ULL; idom < formulation->_surface.size(); ++idom) {
            if(mpi::isItMySubdomain(idom)) {
              for(auto it = formulation->_surface[idom].begin(); it != formulation->_surface[idom].end(); ++it) {
                for(auto itField = formulation->_interfaceFields.begin(); itField != formulation->_interfaceFields.end(); ++itField) {
                  if((*itField->first).isDefined(idom, it->first)) {
                    gmshfem::algebra::Vector< T_Scalar > g_ij;
                    // get the solved g_ij
                    (*itField->first)(idom, it->first).getUnknownVector(g_ij, gmshfem::dofs::RawOrder::Hash);
                    // concatenate to get one vector per process
                    g_local.concatenate(g_ij);
                  }
                }
              }
            }
          }
    
          VecCopy(g_local.getPetsc(), Y);
    
          if(IA) {
            VecAYPX(Y, -1., X);
          }
    
          PetscFunctionReturn(0);
        }
    
        template< class T_Scalar >
        int MatMatProductImpl(Mat A, Mat X, Mat Y, bool IA)
        {
    
          // EXTREMLY UNOPTIMIZED
          PetscFunctionBegin;
    
          PetscInt n, p, n_local, p_local;
          MatGetSize(X, &n, &p);
          MatGetLocalSize(X, &n_local, &p_local);
          Vec Xcol, Ycol;
          VecCreateMPI(MPI_COMM_WORLD, n_local, n, &Xcol);
          VecDuplicate(Xcol, &Ycol);
    
        for (auto j = 0; j < p; ++j) {
          MatGetColumnVector(X, Xcol, j);
          MatVectProductImpl<T_Scalar>(A, Xcol, Ycol, IA);
    
          // Copy to matrix
          PetscInt low, high, i;
          PetscScalar val;
    
          // Get the ownership range for the current rank
          VecGetOwnershipRange(Ycol, &low, &high);
    
          // Loop over local rows and set the value for the specified column
          for (i = low; i < high; ++i) {
              VecGetValues(Ycol, 1, &i, &val);
              MatSetValue(Y, i, j, val, INSERT_VALUES);
          }
    
        }
    
          PetscFunctionReturn(0);
        }
    
        
    
    
        template class Formulation< std::complex< double > >;
        template class Formulation< std::complex< float > >;
    
    
      } // namespace problem
    
    
    } // namespace gmshddm