diff --git a/CMakeLists.txt b/CMakeLists.txt
index 4961ed0bc28eb8e4aa6b533af4f2a00e87cac3c4..263bfb4013fcab1c4f9c46c3a332541137f19638 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -225,6 +225,8 @@ include_directories(src/algebra src/analytics src/common src/common/io src/equat
add_subdirectory(examples)
add_subdirectory(tutorials)
+add_subdirectory(unitTests/correctness)
+
#----------------------------------
#
diff --git a/examples/distributed/HelmholtzCheckerboard/CMakeLists.txt b/examples/distributed/HelmholtzCheckerboard/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3e8d28b40a80596978df92599f594c10f009d6e8
--- /dev/null
+++ b/examples/distributed/HelmholtzCheckerboard/CMakeLists.txt
@@ -0,0 +1,18 @@
+cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
+
+project(tuto CXX)
+
+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})
+endif()
+
+
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/../../example.cmake)
+
+add_executable(tuto main.cpp ${EXTRA_INCS})
+target_link_libraries(tuto ${EXTRA_LIBS})
diff --git a/examples/distributed/HelmholtzCheckerboard/main.cpp b/examples/distributed/HelmholtzCheckerboard/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..666068f379a2877131014ba765849ecaf5a3ddf1
--- /dev/null
+++ b/examples/distributed/HelmholtzCheckerboard/main.cpp
@@ -0,0 +1,595 @@
+#include <PetscInterface.h>
+#include <algorithm>
+#include <chrono>
+#include <gmsh.h>
+#include <gmshfem/CSVio.h>
+#include <gmshfem/DistributedField.h>
+#include <gmshfem/Formulation.h>
+#include <gmshfem/GmshFem.h>
+#include <gmshfem/MatrixFactory.h>
+#include <gmshfem/Solver.h>
+#include <gmshfem/SolverMonitor.h>
+#include <gmshfem/VectorFactory.h>
+#include <thread>
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+using namespace gmshfem::problem;
+using namespace gmshfem::domain;
+using namespace gmshfem::field;
+using namespace gmshfem::function;
+using namespace gmshfem::post;
+using namespace gmshfem::equation;
+using namespace gmshfem::msg;
+
+
+class MonitorL2Error : public gmshfem::system::AbstractMonitor
+{
+ private:
+ Domain dom;
+ DistributedField< std::complex< double >, Form::Form0 > &field;
+ DistributedField< std::complex< double >, Form::Form0 > &UnknownField;
+ double normSol;
+ Formulation< std::complex< double > > &formulation;
+
+ public:
+ MonitorL2Error(Domain dom, DistributedField< std::complex< double >, Form::Form0 > &field, DistributedField< std::complex< double >, Form::Form0 > &UnknownField, Formulation< std::complex< double > > &form) :
+ dom(dom), field(field), UnknownField(UnknownField), formulation(form)
+ {
+ // COMPUTE SOL NORM, REDUCE, IN CALLBACK: set values into field
+ normSol = gmshfem::common::GmshFem::reductionSum(real(integrate((field)*conj(field), dom, "Gauss6")));
+ if(gmshfem::common::GmshFem::getMPIRank() == 0)
+ msg::info << "Initalize L2 monitor with " << normSol << msg::endl;
+ }
+
+ virtual PetscErrorCode Monitor(KSP ksp, PetscInt n, PetscReal rnorm)
+ {
+ auto MPI_Rank = gmshfem::common::GmshFem::getMPIRank();
+
+ Vec x;
+ VecType t;
+ PetscInt size;
+ KSPBuildSolution(ksp, nullptr, &x);
+ VecGetType(x, &t);
+ std::vector< std::complex< double > > values;
+
+ values.resize(formulation.getNumberOfUnknownDof());
+ formulation.getDistributedContext()->readScatteredData(values, x);
+ formulation.setSolutionIntoFields(values);
+
+ double normErr = gmshfem::common::GmshFem::reductionSum(real(integrate((field - UnknownField) * conj(field - UnknownField), dom, "Gauss6")));
+
+
+ if(MPI_Rank == 0)
+ msg::info << "Monitor called " << n << " and rnorm " << rnorm << ". Rel L² error is: " << sqrt(normErr / normSol) << msg::endl;
+
+ return PETSC_SUCCESS;
+ }
+};
+
+void genMesh(double L, double H, unsigned nx, unsigned ny, double lc, unsigned overlapSize, std::vector< std::pair< double, double > > sources)
+{
+ namespace geo = gmsh::model::geo;
+ std::vector< int > sourcePoints, sourceDoms;
+ CSVio out_source_domains("source_subdomains");
+
+
+ gmsh::model::add("genMesh");
+ double dx = L / nx;
+ double dy = H / ny;
+ double delta = lc * overlapSize;
+ if(delta > dx || delta > dy)
+ throw Exception("Overlap too big!");
+
+
+ // points[i][j] with i the x axis and j the y axis. Each subdomain has 4 points in a direction, but offset of 3
+ // because of sharing
+ std::vector< std::vector< int > > points(1 + 3 * nx, std::vector< int >(1 + 3 * ny));
+
+ // All points
+ for(unsigned i = 0; i < nx; ++i) {
+ for(unsigned j = 0; j < ny; ++j) {
+ double x = i * dx;
+ double y = j * dy;
+ points[3 * i][3 * j] = gmsh::model::geo::addPoint(x, y, 0., lc);
+ points[3 * i + 1][3 * j] = gmsh::model::geo::addPoint(x + delta, y, 0., lc);
+ points[3 * i + 2][3 * j] = gmsh::model::geo::addPoint(x + dx - delta, y, 0., lc);
+
+ points[3 * i][3 * j + 1] = gmsh::model::geo::addPoint(x, y + delta, 0., lc);
+ points[3 * i + 1][3 * j + 1] = gmsh::model::geo::addPoint(x + delta, y + delta, 0., lc);
+ points[3 * i + 2][3 * j + 1] = gmsh::model::geo::addPoint(x + dx - delta, y + delta, 0., lc);
+
+ points[3 * i][3 * j + 2] = gmsh::model::geo::addPoint(x, y + dy - delta, 0., lc);
+ points[3 * i + 1][3 * j + 2] = gmsh::model::geo::addPoint(x + delta, y + dy - delta, 0., lc);
+ points[3 * i + 2][3 * j + 2] = gmsh::model::geo::addPoint(x + dx - delta, y + dy - delta, 0., lc);
+
+ if(j == (ny - 1)) {
+ points[3 * i][3 * j + 3] = gmsh::model::geo::addPoint(x, y + dy, 0., lc);
+ points[3 * i + 1][3 * j + 3] = gmsh::model::geo::addPoint(x + delta, y + dy, 0., lc);
+ points[3 * i + 2][3 * j + 3] = gmsh::model::geo::addPoint(x + dx - delta, y + dy, 0., lc);
+ }
+
+ if(i == (nx - 1)) {
+ points[3 * i + 3][3 * j] = gmsh::model::geo::addPoint(x + dx, y, 0., lc);
+ points[3 * i + 3][3 * j + 1] = gmsh::model::geo::addPoint(x + dx, y + delta, 0., lc);
+ points[3 * i + 3][3 * j + 2] = gmsh::model::geo::addPoint(x + dx, y + dy - delta, 0., lc);
+ if(j == (ny - 1)) {
+ points[3 * i + 3][3 * j + 3] = gmsh::model::geo::addPoint(x + dx, y + dy, 0., lc);
+ }
+ }
+ }
+ }
+
+ // There are 24 lines on each subdomain. Some are duplicated!
+ struct SubdomainLines {
+ // All ids in order [bottom, right, up, left]
+ // Lines are always pointing to the right or upward
+ int inner[4]; // Inne square
+ int outer[4]; // Interior parts (corner excluded) of outer square
+ int cornerLines[16]; // 4 * cornerId + edgeId,
+ };
+
+ std::vector< std::vector< SubdomainLines > > lines(nx, std::vector< SubdomainLines >(ny));
+
+
+ // All lines
+ for(unsigned i = 0; i < nx; ++i) {
+ for(unsigned j = 0; j < ny; ++j) {
+ auto &sub = lines[i][j];
+ sub.inner[0] = geo::addLine(points[3 * i + 1][3 * j + 1], points[3 * i + 2][3 * j + 1]);
+ sub.inner[1] = geo::addLine(points[3 * i + 2][3 * j + 1], points[3 * i + 2][3 * j + 2]);
+ sub.inner[2] = geo::addLine(points[3 * i + 1][3 * j + 2], points[3 * i + 2][3 * j + 2]);
+ sub.inner[3] = geo::addLine(points[3 * i + 1][3 * j + 1], points[3 * i + 1][3 * j + 2]);
+
+ // Always create bottom boundary, give it to previous subdomain
+ sub.outer[0] = geo::addLine(points[3 * i + 1][3 * j], points[3 * i + 2][3 * j]);
+ if(j > 0)
+ lines[i][j - 1].outer[2] = sub.outer[0];
+
+ // If I'm at the end, I create the last ones, otherwise I'll read them.
+ if(i == (nx - 1))
+ sub.outer[1] = geo::addLine(points[3 * i + 3][3 * j + 1], points[3 * i + 3][3 * j + 2]);
+ if(j == (ny - 1))
+ sub.outer[2] = geo::addLine(points[3 * i + 1][3 * j + 3], points[3 * i + 2][3 * j + 3]);
+
+ sub.outer[3] = geo::addLine(points[3 * i][3 * j + 1], points[3 * i][3 * j + 2]);
+ if(i > 0)
+ lines[i - 1][j].outer[1] = sub.outer[3];
+
+
+ // Always define the inner corners
+
+ // Bottom left: right and up parts
+ sub.cornerLines[1] = geo::addLine(points[3 * i + 1][3 * j], points[3 * i + 1][3 * j + 1]);
+ sub.cornerLines[2] = geo::addLine(points[3 * i][3 * j + 1], points[3 * i + 1][3 * j + 1]);
+
+ // Bottom right: up and left parts
+ sub.cornerLines[4 + 2] = geo::addLine(points[3 * i + 2][3 * j + 1], points[3 * i + 3][3 * j + 1]);
+ sub.cornerLines[4 + 3] = geo::addLine(points[3 * i + 2][3 * j], points[3 * i + 2][3 * j + 1]);
+
+ // top right: down and left parts
+ sub.cornerLines[8 + 0] = geo::addLine(points[3 * i + 2][3 * j + 2], points[3 * i + 3][3 * j + 2]);
+ sub.cornerLines[8 + 3] = geo::addLine(points[3 * i + 2][3 * j + 2], points[3 * i + 2][3 * j + 3]);
+
+ // top left: down and right parts
+ sub.cornerLines[12 + 0] = geo::addLine(points[3 * i][3 * j + 2], points[3 * i + 1][3 * j + 2]);
+ sub.cornerLines[12 + 1] = geo::addLine(points[3 * i + 1][3 * j + 2], points[3 * i + 1][3 * j + 3]);
+
+ // Shared parts of corners
+ // Left: always create, share if i > 0
+ sub.cornerLines[3] = geo::addLine(points[3 * i][3 * j + 0], points[3 * i][3 * j + 1]);
+ sub.cornerLines[12 + 3] = geo::addLine(points[3 * i][3 * j + 2], points[3 * i][3 * j + 3]);
+ if(i > 0) {
+ lines[i - 1][j].cornerLines[4 + 1] = sub.cornerLines[3];
+ lines[i - 1][j].cornerLines[8 + 1] = sub.cornerLines[12 + 3];
+ }
+
+ // Bottom: always create
+ sub.cornerLines[0] = geo::addLine(points[3 * i][3 * j + 0], points[3 * i + 1][3 * j + 0]);
+ sub.cornerLines[4] = geo::addLine(points[3 * i + 2][3 * j + 0], points[3 * i + 3][3 * j + 0]);
+ if(j > 0) {
+ lines[i][j - 1].cornerLines[12 + 2] = sub.cornerLines[0];
+ lines[i][j - 1].cornerLines[8 + 2] = sub.cornerLines[4];
+ }
+
+ // IF i'm the last col, add the right
+ if(i == (nx - 1)) {
+ sub.cornerLines[4 + 1] = geo::addLine(points[3 * i + 3][3 * j + 0], points[3 * i + 3][3 * j + 1]);
+ sub.cornerLines[8 + 1] = geo::addLine(points[3 * i + 3][3 * j + 2], points[3 * i + 3][3 * j + 3]);
+ }
+
+ if(j == (ny - 1)) {
+ sub.cornerLines[12 + 2] = geo::addLine(points[3 * i + 0][3 * j + 3], points[3 * i + 1][3 * j + 3]);
+ sub.cornerLines[8 + 2] = geo::addLine(points[3 * i + 2][3 * j + 3], points[3 * i + 3][3 * j + 3]);
+ }
+ }
+ }
+
+ // Make the surfaces
+ struct SubdomaineFaces {
+ int inner;
+ int corner[4];
+ int edges[4];
+ };
+
+ std::vector< std::vector< SubdomaineFaces > > surfaces(nx, std::vector< SubdomaineFaces >(ny));
+
+ for(unsigned i = 0; i < nx; ++i) {
+ for(unsigned j = 0; j < ny; ++j) {
+ auto &theseLines = lines[i][j];
+ int loop = geo::addCurveLoop({theseLines.inner[0], theseLines.inner[1], -theseLines.inner[2], -theseLines.inner[3]});
+ surfaces[i][j].inner = geo::addPlaneSurface({loop});
+
+ // "Edges"
+ loop = geo::addCurveLoop({theseLines.outer[0], theseLines.cornerLines[7], -theseLines.inner[0], -theseLines.cornerLines[1]});
+ surfaces[i][j].edges[0] = geo::addPlaneSurface({loop});
+
+ loop = geo::addCurveLoop({-theseLines.inner[1], theseLines.cornerLines[6], theseLines.outer[1], -theseLines.cornerLines[8]});
+ surfaces[i][j].edges[1] = geo::addPlaneSurface({loop});
+
+ loop = geo::addCurveLoop({theseLines.inner[2], theseLines.cornerLines[11], -theseLines.outer[2], -theseLines.cornerLines[13]});
+ surfaces[i][j].edges[2] = geo::addPlaneSurface({loop});
+
+ loop = geo::addCurveLoop({theseLines.inner[3], -theseLines.cornerLines[12], -theseLines.outer[3], theseLines.cornerLines[2]});
+ surfaces[i][j].edges[3] = geo::addPlaneSurface({loop});
+
+ // "Corners"
+ for(int offset = 0; offset < 16; offset += 4) {
+ loop = geo::addCurveLoop({theseLines.cornerLines[offset + 0], theseLines.cornerLines[offset + 1], -theseLines.cornerLines[offset + 2], -theseLines.cornerLines[offset + 3]});
+ surfaces[i][j].corner[offset / 4] = geo::addPlaneSurface({loop});
+ }
+ }
+ }
+
+ // Add source
+ for(auto &pos : sources) {
+ double x = pos.first;
+ double y = pos.second;
+ int iSrc = static_cast< int >(x / dx);
+ int jSrc = static_cast< int >(y / dy);
+ int srcDom = iSrc + jSrc * nx;
+ sourceDoms.push_back(srcDom);
+ sourcePoints.push_back(geo::addPoint(pos.first, pos.second, 0, lc));
+
+ out_source_domains << srcDom << csv::endl;
+
+ gmsh::model::geo::synchronize();
+
+
+ double xloc = x - iSrc * dx;
+ double yloc = y - jSrc * dy;
+
+ auto point = sourcePoints.back();
+ if(xloc > delta && xloc < (dx - delta)) {
+ // Inner
+ if(yloc > delta && yloc < (dy - delta)) {
+ gmsh::model::mesh::embed(0, {point}, 2, surfaces[iSrc][jSrc].inner);
+ }
+ else if(yloc < delta) {
+ // S edge
+ gmsh::model::mesh::embed(0, {point}, 2, surfaces[iSrc][jSrc].edges[0]);
+ }
+ else if(yloc > (dy - delta)) {
+ // N edge
+ gmsh::model::mesh::embed(0, {point}, 2, surfaces[iSrc][jSrc].edges[2]);
+ }
+ else {
+ throw common::Exception("Bad source position");
+ }
+ }
+ else if(xloc < delta) {
+ // W edge
+ if(yloc > delta && yloc < (dy - delta)) {
+ gmsh::model::mesh::embed(0, {point}, 2, surfaces[iSrc][jSrc].edges[3]);
+ }
+ else if(yloc < delta) {
+ // SW corner
+ gmsh::model::mesh::embed(0, {point}, 2, surfaces[iSrc][jSrc].corner[0]);
+ }
+ else if(yloc > (dy - delta)) {
+ // NW corner
+ gmsh::model::mesh::embed(0, {point}, 2, surfaces[iSrc][jSrc].corner[3]);
+ }
+ else {
+ throw common::Exception("Bad source position");
+ }
+ }
+ else if(xloc > (dx - delta)) {
+ // E edge
+ if(yloc > delta && yloc < (dy - delta)) {
+ gmsh::model::mesh::embed(0, {point}, 2, surfaces[iSrc][jSrc].edges[1]);
+ }
+ else if(yloc < delta) {
+ // SE corner
+ gmsh::model::mesh::embed(0, {point}, 2, surfaces[iSrc][jSrc].corner[1]);
+ }
+ else if(yloc > (dy - delta)) {
+ // NE corner
+ gmsh::model::mesh::embed(0, {point}, 2, surfaces[iSrc][jSrc].corner[2]);
+ }
+ else {
+ throw common::Exception("Bad source position");
+ }
+ }
+ }
+
+ //int src = geo::addPoint(xsrc, ysrc, 0, lc);
+ //int iSrc = nx / 2;
+ //int jSrc = ny / 2;
+
+ //gmsh::model::geo::synchronize();
+ //gmsh::model::mesh::embed(0, {src}, 2, surfaces[0][0].inner);
+
+ gmsh::model::mesh::generate();
+ //gmsh::write("full_mesh.msh");
+
+
+ for(unsigned i = 0; i < nx; ++i) {
+ for(unsigned j = 0; j < ny; ++j) {
+ //gmsh::model::removePhysicalGroups();
+
+ std::vector< int > inner, overlap, boundary;
+ // Fill interior, overlap and boundary
+ inner.push_back(surfaces[i][j].inner);
+ for(int k = 0; k < 4; ++k) {
+ inner.push_back(surfaces[i][j].corner[k]);
+ inner.push_back(surfaces[i][j].edges[k]);
+ }
+
+ if(i > 0) {
+ overlap.push_back(surfaces[i - 1][j].edges[1]);
+ overlap.push_back(surfaces[i - 1][j].corner[1]);
+ overlap.push_back(surfaces[i - 1][j].corner[2]);
+
+ // Left bnd
+ boundary.push_back(lines[i - 1][j].inner[1]);
+ boundary.push_back(lines[i - 1][j].cornerLines[4 + 3]);
+ boundary.push_back(lines[i - 1][j].cornerLines[8 + 3]);
+ }
+ else {
+ boundary.push_back(lines[i][j].outer[3]);
+ boundary.push_back(lines[i][j].cornerLines[3]);
+ boundary.push_back(lines[i][j].cornerLines[12 + 3]);
+ }
+
+ if(j > 0) {
+ overlap.push_back(surfaces[i][j - 1].edges[2]);
+ overlap.push_back(surfaces[i][j - 1].corner[2]);
+ overlap.push_back(surfaces[i][j - 1].corner[3]);
+
+ // Bottom bnd
+ boundary.push_back(lines[i][j - 1].inner[2]);
+ boundary.push_back(lines[i][j - 1].cornerLines[8 + 0]);
+ boundary.push_back(lines[i][j - 1].cornerLines[12 + 0]);
+ }
+ else {
+ boundary.push_back(lines[i][j].outer[0]);
+ boundary.push_back(lines[i][j].cornerLines[0]);
+ boundary.push_back(lines[i][j].cornerLines[4]);
+ }
+
+ if(i != (nx - 1)) {
+ overlap.push_back(surfaces[i + 1][j].edges[3]);
+ overlap.push_back(surfaces[i + 1][j].corner[3]);
+ overlap.push_back(surfaces[i + 1][j].corner[0]);
+
+ // Right bnd
+ boundary.push_back(lines[i + 1][j].inner[3]);
+ boundary.push_back(lines[i + 1][j].cornerLines[4 + 1]);
+ boundary.push_back(lines[i + 1][j].cornerLines[8 + 1]);
+ }
+ else {
+ boundary.push_back(lines[i][j].outer[1]);
+ boundary.push_back(lines[i][j].cornerLines[4 + 1]);
+ boundary.push_back(lines[i][j].cornerLines[8 + 1]);
+ }
+
+ if(j != (ny - 1)) {
+ overlap.push_back(surfaces[i][j + 1].edges[0]);
+ overlap.push_back(surfaces[i][j + 1].corner[1]);
+ overlap.push_back(surfaces[i][j + 1].corner[0]);
+
+ // top bnd
+ boundary.push_back(lines[i][j + 1].inner[0]);
+ boundary.push_back(lines[i][j + 1].cornerLines[2]);
+ boundary.push_back(lines[i][j + 1].cornerLines[4 + 2]);
+ }
+ else {
+ // top bnd
+ boundary.push_back(lines[i][j].outer[2]);
+ boundary.push_back(lines[i][j].cornerLines[8 + 2]);
+ boundary.push_back(lines[i][j].cornerLines[12 + 2]);
+ }
+ gmsh::model::geo::synchronize();
+
+ /*
+ Corner of overlaps
+ */
+ if(i != 0 && j != 0) {
+ // SW corner
+ overlap.push_back(surfaces[i - 1][j - 1].corner[2]);
+ }
+ if(i != 0 && j != (ny - 1)) {
+ // NW corner
+ overlap.push_back(surfaces[i - 1][j + 1].corner[1]);
+ }
+ if(i != (nx - 1) && j != 0) {
+ // SE corner
+ overlap.push_back(surfaces[i + 1][j - 1].corner[3]);
+ }
+ if(i != (nx - 1) && j != (ny - 1)) {
+ // NE corner
+ overlap.push_back(surfaces[i + 1][j + 1].corner[0]);
+ }
+
+ gmsh::model::addPhysicalGroup(2, inner, 1 + 100 * (i + j * nx));
+ gmsh::model::addPhysicalGroup(2, overlap, 2 + 100 * (i + j * nx));
+ gmsh::model::addPhysicalGroup(1, boundary, 3 + 100 * (i + j * nx));
+
+ /*gmsh::model::addPhysicalGroup(2, {3}, 2, "overlap");
+ gmsh::model::addPhysicalGroup(1, {l110}, 1, "gammaLat");
+ gmsh::model::addPhysicalGroup(1, {l89, l910}, 2, "gammaTop");
+ gmsh::model::addPhysicalGroup(1, {l12, l23}, 3, "gammaBottom");
+ gmsh::model::addPhysicalGroup(0, {srcLoc}, 1, "src");*/
+ for(int srcIdx = 0; srcIdx < sourceDoms.size(); ++srcIdx) {
+ if(i + nx * j == sourceDoms[srcIdx])
+ gmsh::model::addPhysicalGroup(0, {sourcePoints[srcIdx]}, srcIdx + 1);
+ }
+
+ gmsh::model::geo::synchronize();
+
+ //gmsh::write("mesh_" + std::to_string(i) + "_" + std::to_string(j) + ".msh");
+ }
+ }
+
+ gmsh::write("full_mesh.msh");
+}
+/*
+void solveFullProblem(unsigned nCells, unsigned order, double k)
+{
+ gmsh::open("full_mesh.msh");
+ Domain omega("omega");
+ Domain gamma("gamma");
+ Domain src("src");
+
+ Formulation< std::complex< double > > formulation("formulation_full");
+ Field< std::complex< double >, Form::Form0 > v("v", omega, FunctionSpaceTypeForm0::HierarchicalH1, order);
+
+ formulation.integral(grad(dof(v)), grad(tf(v)), v.domain(), "Gauss6");
+ formulation.integral(-k * k * dof(v), tf(v), v.domain(), "Gauss6");
+ formulation.integral(std::complex< double >(0, 1) * k * dof(v), tf(v), gamma, "Gauss6");
+ formulation.integral(-1, tf(v), src, "Gauss6");
+
+ formulation.pre();
+ formulation.assemble();
+ formulation.solve(false);
+ save(v, omega, "v_full");
+}*/
+
+int main(int argc, char **argv)
+{
+ using namespace std::chrono_literals;
+
+ GmshFem gmshFem(argc, argv);
+ unsigned int MPI_Rank = gmshFem.getMPIRank();
+ unsigned MPI_Size = gmshFem.getMPISize();
+
+ int numSolsToSave = 1;
+ gmshFem.userDefinedParameter(numSolsToSave, "numSolsToSave");
+
+ int order = 1;
+ gmshFem.userDefinedParameter(order, "order");
+ unsigned nx{2}, ny{1};
+
+ //if(MPI_Size != 2)
+ // throw Exception("Must run with 2 ranks! (Got " + std::to_string(MPI_Size) + ")");
+ gmshFem.userDefinedParameter(nx, "nx");
+ gmshFem.userDefinedParameter(ny, "ny");
+ int overlapSize = 1;
+ gmshFem.userDefinedParameter(overlapSize, "overlapSize");
+
+ double wavelength = 0.05;
+ gmshFem.userDefinedParameter(wavelength, "wavelength");
+
+ double k = 2 * 3.1415 / wavelength;
+
+ int iDom = MPI_Rank % nx;
+ int jDom = MPI_Rank / nx;
+
+
+ /****
+ * Gmsh part
+ ****/
+ std::vector<double> x_src = {0.17, 0.22, 0.27, 0.32, 0.37, 0.42, 0.47, 0.52, 0.57, 0.62, 0.67, 0.72, 0.77, 0.82, 0.87};
+ std::vector< std::pair< double, double > > sources;
+ for (double x: x_src)
+ {
+ sources.push_back({x, 0.8});
+ sources.push_back({x, 0.7});
+ sources.push_back({x, 0.1});
+ }
+
+ if(MPI_Rank == 0) {
+ genMesh(1, 1, nx, ny, wavelength / 10, overlapSize, sources);
+ }
+ gmshFem.BarrierMPI();
+
+ std::vector< int > srcDomains, sourcesOwned;
+ CSVio out_source_domains("source_subdomains", ';', OpeningMode::Reading);
+ int tmp;
+ while((out_source_domains >> tmp).isEOF() == false) {
+ srcDomains.push_back(tmp);
+ if(tmp == MPI_Rank)
+ sourcesOwned.push_back(srcDomains.size() - 1);
+ if(MPI_Rank == 0)
+ msg::debug << "Source " << srcDomains.size() - 1 << " goes in " << tmp << msg::endl;
+ }
+
+
+
+ //gmsh::open("mesh_" + std::to_string(iDom) + "_" + std::to_string(jDom) + ".msh");
+ gmsh::open("full_mesh.msh");
+ gmsh::model::mesh::createEdges();
+ gmsh::model::mesh::createFaces();
+
+
+ Formulation< std::complex< double > > formulation("formulation_" + std::to_string(MPI_Rank));
+ Domain omega(2, 1 + 100 * (iDom + jDom * nx));
+ Domain overlap(2, 2 + 100 * (iDom + jDom * nx));
+ Domain boundary(1, 3 + 100 * (iDom + jDom * nx));
+
+ DistributedField< std::complex< double >, Form::Form0 > v("v", omega, overlap, FunctionSpaceTypeForm0::HierarchicalH1, order);
+
+ formulation.integral(grad(dof(v)), grad(tf(v)), v.domain(), "Gauss6");
+ std::complex< double > k2Complex(k * k, -k * k * 0.0);
+ std::complex< double > im(0., 1.0);
+ formulation.integral(-k2Complex * dof(v), tf(v), v.domain(), "Gauss6");
+ formulation.integral(im * k * dof(v), tf(v), boundary, "Gauss6");
+ for(int srcIdx = 0; srcIdx < srcDomains.size(); ++srcIdx) {
+ formulation.setCurrentRHS(srcIdx);
+
+ if(std::find(sourcesOwned.begin(), sourcesOwned.end(), srcIdx) != sourcesOwned.end()) {
+ Domain thisSrc(0, srcIdx + 1);
+ formulation.integral(-1, tf(v), thisSrc, "Gauss6");
+ }
+ }
+
+ formulation.pre();
+ formulation.preDistributed();
+ formulation.assemble();
+ formulation.solveAllDistributed(system::DISTRIBUTED_SOLVE_TYPE::ORAS, false);
+
+ for (unsigned kSol = 0; kSol < formulation.numRHS(); ++kSol) {
+ formulation.loadSolution(kSol);
+ if (kSol < numSolsToSave)
+ save(v, omega, "v_rhs_" + std::to_string(kSol) + "_rank_" + std::to_string(MPI_Rank));
+ }
+
+ // Error analysis on first RHS
+ formulation.loadSolution(0);
+ decltype(v) v_backup = v;
+ std::unique_ptr< system::AbstractMonitor > monitor = std::make_unique< MonitorL2Error >(omega, v_backup, v, formulation);
+ formulation.setMonitor(std::move(monitor));
+ formulation.solveDistributed(system::DISTRIBUTED_SOLVE_TYPE::ORAS, true);
+
+
+ save(v, omega, "v_" + std::to_string(MPI_Rank));
+ save(v - v_backup, omega, "delta_" + std::to_string(MPI_Rank));
+ auto vnorm = real(integrate((v)*conj(v), omega, "Gauss6"));
+ auto error1 = real(integrate((v - v_backup) * conj(v - v_backup), omega, "Gauss6"));
+
+ msg::print << "L2 absolute error in domain " << MPI_Rank << " : " << error1 << " and sol has norm " << vnorm << " so rel error is " << error1 / vnorm << msg::endl;
+ double v_tot = gmshFem.reductionSum(vnorm);
+ double err_tot = gmshFem.reductionSum(error1);
+ if(MPI_Rank == 0)
+ msg::print << "L2 absolute error in domain (total)): " << err_tot << " and sol has norm " << v_tot << " so rel error is " << err_tot / v_tot << msg::endl;
+
+
+ //save(v - solution, omega, "err_" + std::to_string(MPI_Rank));
+ /*
+ if(error1 > 1e-2) {
+ throw common::Exception("Error above tolerance : " + std::to_string(error1));
+ }*/
+
+ return 0;
+}
diff --git a/examples/distributed/Poisson2DRobin/CMakeLists.txt b/examples/distributed/Poisson2DRobin/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..af8f060b31fe272fe3efb2f6f7cb71f56422b74a
--- /dev/null
+++ b/examples/distributed/Poisson2DRobin/CMakeLists.txt
@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
+
+project(tuto CXX)
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/../../example.cmake)
+
+add_executable(tuto main.cpp ${EXTRA_INCS})
+target_link_libraries(tuto ${EXTRA_LIBS})
diff --git a/examples/distributed/Poisson2DRobin/main.cpp b/examples/distributed/Poisson2DRobin/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c088e954ff79d878c26628c28c1d9699c1c216c0
--- /dev/null
+++ b/examples/distributed/Poisson2DRobin/main.cpp
@@ -0,0 +1,437 @@
+#include <algorithm>
+#include <chrono>
+#include <gmsh.h>
+#include <gmshfem/DistributedField.h>
+#include <gmshfem/Formulation.h>
+#include <gmshfem/GmshFem.h>
+#include <gmshfem/MatrixFactory.h>
+#include <gmshfem/VectorFactory.h>
+#include <thread>
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+using namespace gmshfem::problem;
+using namespace gmshfem::domain;
+using namespace gmshfem::field;
+using namespace gmshfem::function;
+using namespace gmshfem::post;
+using namespace gmshfem::equation;
+using namespace gmshfem::msg;
+
+
+void genMesh(double L, double H, unsigned nx, unsigned ny, double lc, unsigned overlapSize)
+{
+ namespace geo = gmsh::model::geo;
+
+
+ gmsh::model::add("genMesh");
+ double dx = L / nx;
+ double dy = H / ny;
+ double delta = lc * overlapSize;
+ if(delta > dx || delta > dy)
+ throw Exception("Overlap too big!");
+
+
+ // points[i][j] with i the x axis and j the y axis. Each subdomain has 4 points in a direction, but offset of 3
+ // because of sharing
+ std::vector< std::vector< int > > points(1 + 3 * nx, std::vector< int >(1 + 3 * ny));
+
+ // All points
+ for(unsigned i = 0; i < nx; ++i) {
+ for(unsigned j = 0; j < ny; ++j) {
+ double x = i * dx;
+ double y = j * dy;
+ points[3 * i][3 * j] = gmsh::model::geo::addPoint(x, y, 0., lc);
+ points[3 * i + 1][3 * j] = gmsh::model::geo::addPoint(x + delta, y, 0., lc);
+ points[3 * i + 2][3 * j] = gmsh::model::geo::addPoint(x + dx - delta, y, 0., lc);
+
+ points[3 * i][3 * j + 1] = gmsh::model::geo::addPoint(x, y + delta, 0., lc);
+ points[3 * i + 1][3 * j + 1] = gmsh::model::geo::addPoint(x + delta, y + delta, 0., lc);
+ points[3 * i + 2][3 * j + 1] = gmsh::model::geo::addPoint(x + dx - delta, y + delta, 0., lc);
+
+ points[3 * i][3 * j + 2] = gmsh::model::geo::addPoint(x, y + dy - delta, 0., lc);
+ points[3 * i + 1][3 * j + 2] = gmsh::model::geo::addPoint(x + delta, y + dy - delta, 0., lc);
+ points[3 * i + 2][3 * j + 2] = gmsh::model::geo::addPoint(x + dx - delta, y + dy - delta, 0., lc);
+
+ if(j == (ny - 1)) {
+ points[3 * i][3 * j + 3] = gmsh::model::geo::addPoint(x, y + dy, 0., lc);
+ points[3 * i + 1][3 * j + 3] = gmsh::model::geo::addPoint(x + delta, y + dy, 0., lc);
+ points[3 * i + 2][3 * j + 3] = gmsh::model::geo::addPoint(x + dx - delta, y + dy, 0., lc);
+ }
+
+ if(i == (nx - 1)) {
+ points[3 * i + 3][3 * j] = gmsh::model::geo::addPoint(x + dx, y, 0., lc);
+ points[3 * i + 3][3 * j + 1] = gmsh::model::geo::addPoint(x + dx, y + delta, 0., lc);
+ points[3 * i + 3][3 * j + 2] = gmsh::model::geo::addPoint(x + dx, y + dy - delta, 0., lc);
+ if(j == (ny - 1)) {
+ points[3 * i + 3][3 * j + 3] = gmsh::model::geo::addPoint(x + dx, y + dy, 0., lc);
+ }
+ }
+ }
+ }
+
+ // There are 24 lines on each subdomain. Some are duplicated!
+ struct SubdomainLines {
+ // All ids in order [bottom, right, up, left]
+ // Lines are always pointing to the right or upward
+ int inner[4]; // Inne square
+ int outer[4]; // Interior parts (corner excluded) of outer square
+ int cornerLines[16]; // 4 * cornerId + edgeId,
+ };
+
+ std::vector< std::vector< SubdomainLines > > lines(nx, std::vector< SubdomainLines >(ny));
+
+
+ // All lines
+ for(unsigned i = 0; i < nx; ++i) {
+ for(unsigned j = 0; j < ny; ++j) {
+ auto &sub = lines[i][j];
+ sub.inner[0] = geo::addLine(points[3 * i + 1][3 * j + 1], points[3 * i + 2][3 * j + 1]);
+ sub.inner[1] = geo::addLine(points[3 * i + 2][3 * j + 1], points[3 * i + 2][3 * j + 2]);
+ sub.inner[2] = geo::addLine(points[3 * i + 1][3 * j + 2], points[3 * i + 2][3 * j + 2]);
+ sub.inner[3] = geo::addLine(points[3 * i + 1][3 * j + 1], points[3 * i + 1][3 * j + 2]);
+
+ // Always create bottom boundary, give it to previous subdomain
+ sub.outer[0] = geo::addLine(points[3 * i + 1][3 * j], points[3 * i + 2][3 * j]);
+ if(j > 0)
+ lines[i][j - 1].outer[2] = sub.outer[0];
+
+ // If I'm at the end, I create the last ones, otherwise I'll read them.
+ if(i == (nx - 1))
+ sub.outer[1] = geo::addLine(points[3 * i + 3][3 * j + 1], points[3 * i + 3][3 * j + 2]);
+ if(j == (ny - 1))
+ sub.outer[2] = geo::addLine(points[3 * i + 1][3 * j + 3], points[3 * i + 2][3 * j + 3]);
+
+ sub.outer[3] = geo::addLine(points[3 * i][3 * j + 1], points[3 * i][3 * j + 2]);
+ if(i > 0)
+ lines[i - 1][j].outer[1] = sub.outer[3];
+
+
+ // Always define the inner corners
+
+ // Bottom left: right and up parts
+ sub.cornerLines[1] = geo::addLine(points[3 * i + 1][3 * j], points[3 * i + 1][3 * j + 1]);
+ sub.cornerLines[2] = geo::addLine(points[3 * i][3 * j + 1], points[3 * i + 1][3 * j + 1]);
+
+ // Bottom right: up and left parts
+ sub.cornerLines[4 + 2] = geo::addLine(points[3 * i + 2][3 * j + 1], points[3 * i + 3][3 * j + 1]);
+ sub.cornerLines[4 + 3] = geo::addLine(points[3 * i + 2][3 * j], points[3 * i + 2][3 * j + 1]);
+
+ // top right: down and left parts
+ sub.cornerLines[8 + 0] = geo::addLine(points[3 * i + 2][3 * j + 2], points[3 * i + 3][3 * j + 2]);
+ sub.cornerLines[8 + 3] = geo::addLine(points[3 * i + 2][3 * j + 2], points[3 * i + 2][3 * j + 3]);
+
+ // top left: down and right parts
+ sub.cornerLines[12 + 0] = geo::addLine(points[3 * i][3 * j + 2], points[3 * i + 1][3 * j + 2]);
+ sub.cornerLines[12 + 1] = geo::addLine(points[3 * i + 1][3 * j + 2], points[3 * i + 1][3 * j + 3]);
+
+ // Shared parts of corners
+ // Left: always create, share if i > 0
+ sub.cornerLines[3] = geo::addLine(points[3 * i][3 * j + 0], points[3 * i][3 * j + 1]);
+ sub.cornerLines[12 + 3] = geo::addLine(points[3 * i][3 * j + 2], points[3 * i][3 * j + 3]);
+ if(i > 0) {
+ lines[i - 1][j].cornerLines[4 + 1] = sub.cornerLines[3];
+ lines[i - 1][j].cornerLines[8 + 1] = sub.cornerLines[12 + 3];
+ }
+
+ // Bottom: always create
+ sub.cornerLines[0] = geo::addLine(points[3 * i][3 * j + 0], points[3 * i + 1][3 * j + 0]);
+ sub.cornerLines[4] = geo::addLine(points[3 * i + 2][3 * j + 0], points[3 * i + 3][3 * j + 0]);
+ if(j > 0) {
+ lines[i][j - 1].cornerLines[12 + 2] = sub.cornerLines[0];
+ lines[i][j - 1].cornerLines[8 + 2] = sub.cornerLines[4];
+ }
+
+ // IF i'm the last col, add the right
+ if(i == (nx - 1)) {
+ sub.cornerLines[4 + 1] = geo::addLine(points[3 * i + 3][3 * j + 0], points[3 * i + 3][3 * j + 1]);
+ sub.cornerLines[8 + 1] = geo::addLine(points[3 * i + 3][3 * j + 2], points[3 * i + 3][3 * j + 3]);
+ }
+
+ if(j == (ny - 1)) {
+ sub.cornerLines[12 + 2] = geo::addLine(points[3 * i + 0][3 * j + 3], points[3 * i + 1][3 * j + 3]);
+ sub.cornerLines[8 + 2] = geo::addLine(points[3 * i + 2][3 * j + 3], points[3 * i + 3][3 * j + 3]);
+ }
+ }
+ }
+
+ // Make the surfaces
+ struct SubdomaineFaces {
+ int inner;
+ int corner[4];
+ int edges[4];
+ };
+
+ std::vector< std::vector< SubdomaineFaces > > surfaces(nx, std::vector< SubdomaineFaces >(ny));
+
+ for(unsigned i = 0; i < nx; ++i) {
+ for(unsigned j = 0; j < ny; ++j) {
+ auto &theseLines = lines[i][j];
+ int loop = geo::addCurveLoop({theseLines.inner[0], theseLines.inner[1], -theseLines.inner[2], -theseLines.inner[3]});
+ surfaces[i][j].inner = geo::addPlaneSurface({loop});
+
+ // "Edges"
+ loop = geo::addCurveLoop({theseLines.outer[0], theseLines.cornerLines[7], -theseLines.inner[0], -theseLines.cornerLines[1]});
+ surfaces[i][j].edges[0] = geo::addPlaneSurface({loop});
+
+ loop = geo::addCurveLoop({-theseLines.inner[1], theseLines.cornerLines[6], theseLines.outer[1], -theseLines.cornerLines[8]});
+ surfaces[i][j].edges[1] = geo::addPlaneSurface({loop});
+
+ loop = geo::addCurveLoop({theseLines.inner[2], theseLines.cornerLines[11], -theseLines.outer[2], -theseLines.cornerLines[13]});
+ surfaces[i][j].edges[2] = geo::addPlaneSurface({loop});
+
+ loop = geo::addCurveLoop({theseLines.inner[3], -theseLines.cornerLines[12], -theseLines.outer[3], theseLines.cornerLines[2]});
+ surfaces[i][j].edges[3] = geo::addPlaneSurface({loop});
+
+ // "Corners"
+ for(int offset = 0; offset < 16; offset += 4) {
+ loop = geo::addCurveLoop({theseLines.cornerLines[offset + 0], theseLines.cornerLines[offset + 1], -theseLines.cornerLines[offset + 2], -theseLines.cornerLines[offset + 3]});
+ surfaces[i][j].corner[offset / 4] = geo::addPlaneSurface({loop});
+ }
+ }
+ }
+
+
+ gmsh::model::geo::synchronize();
+
+ gmsh::model::mesh::generate();
+ //gmsh::write("full_mesh.msh");
+
+
+ for(unsigned i = 0; i < nx; ++i) {
+ for(unsigned j = 0; j < ny; ++j) {
+ //gmsh::model::removePhysicalGroups();
+
+ std::vector< int > inner, overlap, boundary, coupling;
+ // Fill interior, overlap and boundary
+ inner.push_back(surfaces[i][j].inner);
+ for(int k = 0; k < 4; ++k) {
+ inner.push_back(surfaces[i][j].corner[k]);
+ inner.push_back(surfaces[i][j].edges[k]);
+ }
+
+ if(i > 0) {
+ overlap.push_back(surfaces[i - 1][j].edges[1]);
+ overlap.push_back(surfaces[i - 1][j].corner[1]);
+ overlap.push_back(surfaces[i - 1][j].corner[2]);
+
+ // Left bnd
+ coupling.push_back(lines[i - 1][j].inner[1]);
+ coupling.push_back(lines[i - 1][j].cornerLines[4 + 3]);
+ coupling.push_back(lines[i - 1][j].cornerLines[8 + 3]);
+
+ if(j != 0)
+ coupling.push_back(lines[i - 1][j - 1].cornerLines[8 + 3]);
+ if(j != (ny - 1))
+ coupling.push_back(lines[i - 1][j + 1].cornerLines[4 + 3]);
+ }
+ else {
+ boundary.push_back(lines[i][j].outer[3]);
+ boundary.push_back(lines[i][j].cornerLines[3]);
+ boundary.push_back(lines[i][j].cornerLines[12 + 3]);
+
+ if(j != 0)
+ boundary.push_back(lines[i][j - 1].cornerLines[15]);
+ if(j != (ny - 1))
+ boundary.push_back(lines[i][j + 1].cornerLines[3]);
+ }
+
+ if(j > 0) {
+ overlap.push_back(surfaces[i][j - 1].edges[2]);
+ overlap.push_back(surfaces[i][j - 1].corner[2]);
+ overlap.push_back(surfaces[i][j - 1].corner[3]);
+
+ // Bottom bnd
+ coupling.push_back(lines[i][j - 1].inner[2]);
+ coupling.push_back(lines[i][j - 1].cornerLines[8 + 0]);
+ coupling.push_back(lines[i][j - 1].cornerLines[12 + 0]);
+
+ if(i != 0)
+ coupling.push_back(lines[i - 1][j - 1].cornerLines[8]);
+ if(i != (nx - 1))
+ coupling.push_back(lines[i + 1][j - 1].cornerLines[12]);
+ }
+ else {
+ boundary.push_back(lines[i][j].outer[0]);
+ boundary.push_back(lines[i][j].cornerLines[0]);
+ boundary.push_back(lines[i][j].cornerLines[4]);
+
+ if(i != 0)
+ boundary.push_back(lines[i - 1][j].cornerLines[4]);
+ if(i != (nx - 1))
+ boundary.push_back(lines[i + 1][j].cornerLines[0]);
+ }
+
+ if(i != (nx - 1)) {
+ overlap.push_back(surfaces[i + 1][j].edges[3]);
+ overlap.push_back(surfaces[i + 1][j].corner[3]);
+ overlap.push_back(surfaces[i + 1][j].corner[0]);
+
+ // Right bnd
+ coupling.push_back(lines[i + 1][j].inner[3]);
+ coupling.push_back(lines[i + 1][j].cornerLines[1]);
+ coupling.push_back(lines[i + 1][j].cornerLines[12 + 1]);
+
+ if(j != 0)
+ coupling.push_back(lines[i + 1][j - 1].cornerLines[12 + 1]);
+ if(j != (ny - 1))
+ coupling.push_back(lines[i + 1][j + 1].cornerLines[1]);
+ }
+ else {
+ boundary.push_back(lines[i][j].outer[1]);
+ boundary.push_back(lines[i][j].cornerLines[4 + 1]);
+ boundary.push_back(lines[i][j].cornerLines[8 + 1]);
+
+ if(j != 0)
+ boundary.push_back(lines[i][j - 1].cornerLines[8 + 1]);
+ if(j != (ny - 1))
+ boundary.push_back(lines[i][j + 1].cornerLines[4 + 1]);
+ }
+
+ if(j != (ny - 1)) {
+ overlap.push_back(surfaces[i][j + 1].edges[0]);
+ overlap.push_back(surfaces[i][j + 1].corner[1]);
+ overlap.push_back(surfaces[i][j + 1].corner[0]);
+
+ // top bnd
+ coupling.push_back(lines[i][j + 1].inner[0]);
+ coupling.push_back(lines[i][j + 1].cornerLines[2]);
+ coupling.push_back(lines[i][j + 1].cornerLines[4 + 2]);
+
+ if(i != 0)
+ coupling.push_back(lines[i - 1][j + 1].cornerLines[4 + 2]);
+ if(i != (nx - 1))
+ coupling.push_back(lines[i + 1][j + 1].cornerLines[2]);
+ }
+ else {
+ // top bnd
+ boundary.push_back(lines[i][j].outer[2]);
+ boundary.push_back(lines[i][j].cornerLines[8 + 2]);
+ boundary.push_back(lines[i][j].cornerLines[12 + 2]);
+
+ if(i != 0)
+ boundary.push_back(lines[i - 1][j].cornerLines[10]);
+ if(i != (nx - 1))
+ boundary.push_back(lines[i + 1][j].cornerLines[14]);
+ }
+
+ /*
+ Corner of overlaps
+ */
+ if(i != 0 && j != 0) {
+ // SW corner
+ overlap.push_back(surfaces[i - 1][j - 1].corner[2]);
+ }
+ if(i != 0 && j != (ny - 1)) {
+ // NW corner
+ overlap.push_back(surfaces[i - 1][j + 1].corner[1]);
+ }
+ if(i != (nx - 1) && j != 0) {
+ // SE corner
+ overlap.push_back(surfaces[i + 1][j - 1].corner[3]);
+ }
+ if(i != (nx - 1) && j != (ny - 1)) {
+ // NE corner
+ overlap.push_back(surfaces[i + 1][j + 1].corner[0]);
+ }
+
+ gmsh::model::geo::synchronize();
+
+ gmsh::model::addPhysicalGroup(2, inner, 1 + 100 * (i + j * nx));
+ gmsh::model::addPhysicalGroup(2, overlap, 2 + 100 * (i + j * nx));
+ gmsh::model::addPhysicalGroup(1, coupling, 3 + 100 * (i + j * nx));
+ gmsh::model::addPhysicalGroup(1, boundary, 4 + 100 * (i + j * nx));
+
+ gmsh::model::geo::synchronize();
+
+ //gmsh::write("mesh_" + std::to_string(i) + "_" + std::to_string(j) + ".msh");
+ }
+ }
+
+ gmsh::write("full_mesh.msh");
+}
+
+int main(int argc, char **argv)
+{
+ using namespace std::chrono_literals;
+
+ GmshFem gmshFem(argc, argv);
+ unsigned int MPI_Rank = gmshFem.getMPIRank();
+ unsigned MPI_Size = gmshFem.getMPISize();
+
+
+ int order = 1;
+ gmshFem.userDefinedParameter(order, "order");
+ unsigned nx{2}, ny{1};
+ gmshFem.userDefinedParameter(nx, "nx");
+ gmshFem.userDefinedParameter(ny, "ny");
+ if(MPI_Size != nx * ny)
+ throw Exception("Must run with same number of subdomains and ranks! (Got " + std::to_string(MPI_Size) + ")");
+
+ int overlapSize = 1;
+ gmshFem.userDefinedParameter(overlapSize, "overlapSize");
+
+ double h = 0.01;
+ gmshFem.userDefinedParameter(h, "h");
+
+ int iDom = MPI_Rank % nx;
+ int jDom = MPI_Rank / nx;
+
+
+ /****
+ * Gmsh part
+ ****/
+ if(MPI_Rank == 0) {
+ genMesh(1, 1, nx, ny, h, overlapSize);
+ }
+ gmshFem.BarrierMPI();
+
+ //gmsh::open("mesh_" + std::to_string(iDom) + "_" + std::to_string(jDom) + ".msh");
+ gmsh::open("full_mesh.msh");
+ gmsh::model::mesh::createEdges();
+ gmsh::model::mesh::createFaces();
+
+ auto analytical = 1 + x< std::complex< double > >() + 2 * x< std::complex< double > >() * y< std::complex< double > >() + pow(y< std::complex< double > >(), 3);
+
+ Formulation< std::complex< double > > formulation("formulation_" + std::to_string(MPI_Rank));
+ Domain omega(2, 1 + 100 * (iDom + jDom * nx));
+ Domain overlap(2, 2 + 100 * (iDom + jDom * nx));
+ Domain coupling(1, 3 + 100 * (iDom + jDom * nx));
+ Domain boundary(1, 4 + 100 * (iDom + jDom * nx));
+
+ DistributedField< std::complex< double >, Form::Form0 > v("v", omega, overlap, FunctionSpaceTypeForm0::HierarchicalH1, order);
+
+ // "Convection" coefficient
+ double alpha = 1;
+ gmshFem.userDefinedParameter(alpha, "alpha");
+
+ formulation.integral(grad(dof(v)), grad(tf(v)), v.domain(), "Gauss6");
+ formulation.integral(alpha * dof(v), tf(v), coupling, "Gauss6");
+ formulation.integral(6. * y< std::complex< double > >(), (tf(v)), v.domain(), "Gauss6");
+
+ v.addConstraint(boundary, analytical);
+
+ formulation.pre();
+ formulation.preDistributed();
+ formulation.assemble();
+ formulation.solveDistributed(system::DISTRIBUTED_SOLVE_TYPE::ORAS, false);
+
+
+ auto errorLoc = real(integrate((v - analytical) * conj(v - analytical), omega, "Gauss6"));
+ auto solNormLoc = real(integrate((analytical) * conj(analytical), omega, "Gauss6"));
+ double l2_rel_error = gmshFem.reductionSum(errorLoc) / gmshFem.reductionSum(solNormLoc);
+
+ if(MPI_Rank == 0)
+ msg::print << "L2 relative error in the whole domain: " << l2_rel_error << msg::endl;
+ save(v, omega, "v_" + std::to_string(MPI_Rank));
+ save(analytical, omega, "sol_" + std::to_string(MPI_Rank));
+
+ double tol = 1e-2;
+ if(order > 1)
+ tol = 1e-4;
+ if(l2_rel_error > 1e-2) {
+ throw common::Exception("Error above tolerance : " + std::to_string(tol));
+ }
+
+ return 0;
+}
diff --git a/src/field/DistributedField.cpp b/src/field/DistributedField.cpp
index dc3068525834c9bf1e56a454087bcccdf12e317c..94d4c029e02210ce1852e92de5aa223815f451a2 100644
--- a/src/field/DistributedField.cpp
+++ b/src/field/DistributedField.cpp
@@ -121,7 +121,7 @@ namespace gmshfem::field
std::vector< ShareDofInfo > local, global;
local.reserve(_sharedDofs.size());
-
+
unsigned rank = gmshfem::common::GmshFem::getMPIRank();
unsigned commSize = gmshfem::common::GmshFem::getMPISize();
@@ -236,7 +236,6 @@ namespace gmshfem::field
_toSend.clear();
_toRead.clear();
-
msg::debug << "On rank " << rank << " indexToDof has size: " << _indexToDof.size() << msg::endl;
MPI_Barrier(MPI_COMM_WORLD);
@@ -354,7 +353,6 @@ namespace gmshfem::field
void DistributedField< T_Scalar, T_Form >::_computeToSend()
{
// Send all owned interfaces and inner dofs that are read by someone else
-
std::vector< DofGlobalIndex > local, global;
local.reserve(_sharedDofs.size());
diff --git a/src/problem/Formulation.cpp b/src/problem/Formulation.cpp
index 14b7dea98066e6fb99ee57f353bf88fe3418cbe4..0e141d2d1876edfef26e4018af6d8c3863e62f1a 100644
--- a/src/problem/Formulation.cpp
+++ b/src/problem/Formulation.cpp
@@ -219,9 +219,28 @@ namespace gmshfem::problem
void Formulation< T_Scalar >::setCurrentRHS(unsigned idx)
{
_currentRHS = idx;
+ auto oldNumRHS = _numRHS;
_numRHS = std::max(idx + 1, _numRHS);
}
+ template< class T_Scalar >
+ void Formulation< T_Scalar >::forceRHSCount(unsigned idx)
+ {
+ _numRHS = idx;
+ _multiB.clear();
+ _multiB.resize(numRHS());
+ for(auto &bi : _multiB) {
+ bi.init(_dofs.nbrUnknownDofs());
+ }
+ }
+
+ template< class T_Scalar >
+ void Formulation< T_Scalar >::setMonitor(std::unique_ptr<system::AbstractMonitor> &&ptr)
+ {
+ _solver->setMonitor(std::move(ptr));
+ }
+
+
// Bilinear terms
template< class T_Scalar >
unsigned int Formulation< T_Scalar >::integral(const equation::EquationInterface< T_Scalar, Degree::Degree0, field::Form::Form0 > &equationLhs, const equation::EquationInterface< T_Scalar, Degree::Degree0, field::Form::Form0 > &equationRhs, const domain::GeometricObject &domain, const std::string &integrationType, const term::ProductType productType)
@@ -1272,12 +1291,14 @@ namespace gmshfem::problem
msg::info << "Assembling " << _name << "..." << msg::endl;
if(numRHS() != _multiB.size()) {
- _multiB.clear();
_multiB.resize(numRHS());
}
for(auto &bi : _multiB) {
- bi.init(_dofs.nbrUnknownDofs());
+ // If previous RHS have the good size, don't touch them
+ if(bi.size() != _dofs.nbrUnknownDofs())
+ bi.init(_dofs.nbrUnknownDofs());
}
+
s_setLinkedIndices(_dofs.nbrLinkedDofs(), _unknownFields, _A, _multiB);
@@ -1607,6 +1628,85 @@ namespace gmshfem::problem
return time;
}
+ template< class T_Scalar >
+ common::Timer Formulation< T_Scalar >::solveDistributed(system::DISTRIBUTED_SOLVE_TYPE type, const bool reusePreconditioner)
+ {
+ if (gmshfem::common::GmshFem::getMPIRank() == 0)
+ msg::info << "Solving distributed system " << _name << "..." << msg::endl;
+
+ common::Timer time;
+ time.tick();
+
+ std::vector< T_Scalar > values, valuesReordered(_dofs.nbrUnknownDofs(), 0);
+
+ try {
+ if (!_distributedContext)
+ throw common::Exception("Missing distributed context");
+ _solver->solveDistributed(values, type, *_distributedContext, reusePreconditioner);
+ }
+ catch(const std::exception &exc) {
+ msg::error << "Unable to solve the system" << msg::endl;
+ msg::error << exc.what() << msg::endl;
+ msg::info << "Solving aborted" << msg::endl;
+ time.tock();
+ return time;
+ }
+
+ // Copy all dofs
+ if(values.size() != _dofs.nbrUnknownDofs()) {
+ msg::error << "An error occurred while solving the system. Expected a vector of size " << _dofs.nbrUnknownDofs() << " but got " << values.size() << "." << msg::endl;
+ msg::info << "Solving aborted" << msg::endl;
+ time.tock();
+ return time;
+ }
+
+ setSolutionIntoFields((values));
+
+ time.tock();
+ if(gmshfem::common::GmshFem::getMPIRank() == 0)
+ msg::info << "Done solving in " << time << "s" << msg::endl;
+ return time;
+ }
+
+ template< class T_Scalar >
+ common::Timer Formulation< T_Scalar >::solveAllDistributed(system::DISTRIBUTED_SOLVE_TYPE solveType, const bool reusePreconditioner)
+ {
+ if(gmshfem::common::GmshFem::getMPIRank() == 0)
+
+ msg::info << "Solving distributed system" << _name << " with " << _numRHS << " right-hand sides..." << msg::endl;
+
+ common::Timer time;
+ time.tick();
+
+ _solutions.resize(_numRHS);
+
+ try {
+ if (!_distributedContext)
+ throw common::Exception("Missing distributed context");
+ _solver->solveAllDistributed(_solutions, solveType, *_distributedContext, reusePreconditioner);
+ }
+ catch(const std::exception &exc) {
+ msg::error << "Unable to solve the system" << msg::endl;
+ msg::error << exc.what() << msg::endl;
+ msg::info << "Solving aborted" << msg::endl;
+ time.tock();
+ return time;
+ }
+
+ if(_solutions.back().size() != _dofs.nbrUnknownDofs()) {
+ msg::error << "An error occurred while solving the system. Solution size is " << _solutions.back().size() << " but expected " << _dofs.nbrUnknownDofs() << msg::endl;
+ msg::info << "Solving aborted" << msg::endl;
+ time.tock();
+ return time;
+ }
+
+
+ time.tock();
+ if(gmshfem::common::GmshFem::getMPIRank() == 0)
+ msg::info << "Done solving in " << time << "s" << msg::endl;
+ return time;
+ }
+
template< class T_Scalar >
void Formulation< T_Scalar >::setSolutionIntoFields(const std::vector< T_Scalar > &solution)
{
@@ -1629,6 +1729,7 @@ namespace gmshfem::problem
}
}
+
template< class T_Scalar >
const system::DistributedContext< T_Scalar >* Formulation< T_Scalar >::getDistributedContext() const
{
@@ -2049,6 +2150,12 @@ namespace gmshfem::problem
s_getRHSBlock(vector, tfTags, &_multiB.at(0));
}
+ template< class T_Scalar >
+ const std::vector< std::vector< T_Scalar > > &Formulation< T_Scalar >::getSolutionsVector() const
+ {
+ return _solutions;
+ }
+
template< class T_Scalar >
void Formulation< T_Scalar >::setRHS(const algebra::Vector< T_Scalar > &vector)
{
diff --git a/src/problem/Formulation.h b/src/problem/Formulation.h
index 736e81552f4d55f9ee1bf0e0bb364211d963e79e..4e5de7b34cb501dc531ece00235613a529b4cfc9 100644
--- a/src/problem/Formulation.h
+++ b/src/problem/Formulation.h
@@ -108,6 +108,10 @@ namespace gmshfem::problem
// Multi RHS interface
unsigned numRHS() const;
void setCurrentRHS(unsigned idx);
+ void forceRHSCount(unsigned idx);
+
+ void setMonitor(std::unique_ptr<system::AbstractMonitor>&& ptr);
+ void disableSolverShellOptions();
// Bilinear term
unsigned int integral(const equation::EquationInterface< T_Scalar, Degree::Degree0, field::Form::Form0 > &equationLhs, const equation::EquationInterface< T_Scalar, Degree::Degree0, field::Form::Form0 > &equationRhs, const domain::GeometricObject &domain, const std::string &integrationType, const term::ProductType productType = term::ProductType::Hermitian);
diff --git a/src/system/MatrixFactory.cpp b/src/system/MatrixFactory.cpp
index 626c5009f664a746f7598d581f229d94703de6ee..a4ae130d988b1604f95557a5715c7eeec4864cc5 100644
--- a/src/system/MatrixFactory.cpp
+++ b/src/system/MatrixFactory.cpp
@@ -8,13 +8,13 @@
#include "CSVio.h"
#include "Dof.h"
#include "Exception.h"
+#include "GmshFem.h"
#include "KahanSum.h"
#include "Message.h"
#include "PPMio.h"
#include "PetscInterface.h"
#include "gmshfemDefines.h"
#include "instantiate.h"
-
namespace gmshfem::system
{
@@ -569,6 +569,111 @@ namespace gmshfem::system
}
return matPetsc;
}
+
+#ifdef HAVE_MPI
+
+ template< class T_Scalar >
+ Mat MatrixFactory< T_Scalar >::getPetscDistributed(const DistributedContext<T_Scalar>& distributedContext) const
+ {
+ const auto& localDofs = distributedContext.localIDofOwned();
+ const auto &localToGlobal = distributedContext.localToGlobal();
+ unsigned nLocalRows = localDofs.size();
+
+ // PETSc can determine the number of columns because the matrix is square
+ std::vector< PetscInt > ai, aj;
+ ai.reserve(nLocalRows + 1);
+ ai.push_back(0);
+ std::vector< PetscScalar > values;
+ unsigned long long k = 0;
+ for(unsigned i = 0; i < nLocalRows; ++i) {
+ // Copy the row
+ unsigned long long rowToCopy = localDofs[i] - 1; //numDofs start at 1
+ ai.push_back(k + _ai[rowToCopy + 1] - _ai[rowToCopy]);
+ for(unsigned l = 0; l < _ai[rowToCopy + 1] - _ai[rowToCopy]; ++l) {
+ // All entries from the row get copied, and columns are adapted
+ values.push_back((*_module)[_ai[rowToCopy] + l]);
+ aj.push_back(localToGlobal.at(_aj[_ai[rowToCopy] + l]) - 1);
+ ++k;
+ }
+ }
+
+ Mat matPetsc;
+ msg::debug << "Distributed matrix. Local rows: " << nLocalRows << " ai.size(), aj.size(), values.size() : " << ai.size() << ',' << aj.size() << ',' << values.size() << msg::endl;
+ MatCreateMPIAIJWithArrays(PETSC_COMM_WORLD, nLocalRows, nLocalRows, PETSC_DETERMINE, PETSC_DETERMINE, ai.data(), aj.data(), values.data(), &matPetsc);
+
+ return matPetsc;
+ }
+ template< class T_Scalar >
+ void MatrixFactory< T_Scalar >::modifyPetscORAS(Mat mat, const std::vector< unsigned long long > &localIdOfOwned, const std::vector< unsigned long long > &localIdOfNonOwned, const std::vector< unsigned long long > &localToGlobalMap, const IS row, const IS col) const
+ {
+ unsigned nLocalRows = localIdOfOwned.size();
+
+ auto rank = common::GmshFem::getMPIRank();
+
+ const PetscInt *rowsIdx;
+ PetscInt nRows, nrowsGlob;
+ ISGetIndices(row, &rowsIdx);
+ ISGetLocalSize(row, &nRows);
+ ISGetSize(row, &nrowsGlob);
+
+ PetscInt n, m;
+ MatGetSize(mat, &n, &m);
+
+ common::GmshFem::BarrierMPI();
+
+ MatAssemblyBegin(mat, MAT_FLUSH_ASSEMBLY);
+ MatAssemblyEnd(mat, MAT_FLUSH_ASSEMBLY);
+
+ common::GmshFem::BarrierMPI();
+
+ /**
+ * TODO optimisation:
+ * - Sort ISSort pour qu'il soit ordonnée -> std::lower_bound etc
+ * - MatSetValues
+ */
+
+ for(unsigned k = 0; k < nRows; ++k) {
+ PetscInt globalDof = rowsIdx[k] + 1; // 1-Indexing in FEM code
+ auto it = std::find_if(localIdOfNonOwned.cbegin(), localIdOfNonOwned.cend(), [globalDof, &localToGlobalMap](unsigned long long loc) { return localToGlobalMap.at(loc - 1) == globalDof; });
+ if(it != localIdOfNonOwned.cend()) {
+ auto globalIdx = it - localIdOfNonOwned.begin(); // i if dof is localIdOfNonOwned[i]
+ //msg::print << rank << ": Modifying global DOF " << globalDof << " which is the " << globalIdx << "th global dof ( k is " << k << ") and local ID is " << *it << msg::endl;
+
+ // "k" is a local line to be updated, *it is the local ID of the line (in 1-indexing)
+ auto rowToCopy = *it - 1;
+ for(unsigned l = 0; l < _ai[rowToCopy + 1] - _ai[rowToCopy]; ++l) {
+ // All entries from the row get copied, and columns are adapted
+ //values.push_back((*_module)[_ai[rowToCopy] + l]);
+ auto val = (*_module)[_ai[rowToCopy] + l];
+ auto globCol = localToGlobalMap.at(_aj[_ai[rowToCopy] + l]) - 1;
+ auto colPtr = std::find(rowsIdx, rowsIdx + nRows, globCol);
+ if(colPtr == rowsIdx + nRows)
+ throw common::Exception("Missing column ?");
+ auto j = colPtr - rowsIdx;
+ //aj.push_back(localToGlobal.at(_aj[_ai[rowToCopy] + l]) - 1);
+ //msg::print << '[' << rank << "] changes entry " << k <<','<< j << msg::endl;
+ PetscScalar replaced;
+ //MatGetValue(mat, k, j, &replaced);
+ MatSetValue(mat, k, j, val, INSERT_VALUES);
+
+ //msg::print << "Replacing in " << k << ',' << j << " value " << replaced << " by " << val << msg::endl;
+ }
+ }
+ else {
+ //msg::print << rank << ": Global dof untouched " << globalDof << msg::endl;
+ }
+
+ }
+
+ // Reassemble the matrix
+ MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
+ common::GmshFem::BarrierMPI();
+
+ ISRestoreIndices(row, &rowsIdx);
+ }
+#endif
+
#endif
template< class T_Scalar >
diff --git a/src/system/MatrixFactory.h b/src/system/MatrixFactory.h
index ba74876a3ec38e109964b89b547bbacd20b76e39..3ab9ca90f40a85ff19290ad337133af7e6484e92 100644
--- a/src/system/MatrixFactory.h
+++ b/src/system/MatrixFactory.h
@@ -15,10 +15,12 @@
#include "OmpInterface.h"
#include "gmshfemDefines.h"
#include "scalar.h"
+#include "DistributedContext.h"
#include <vector>
typedef struct _p_Mat *Mat;
+typedef struct _p_IS *IS;
namespace gmshfem::system
{
@@ -72,6 +74,10 @@ namespace gmshfem::system
bool havePattern() const;
#ifdef HAVE_PETSC
Mat getPetsc() const;
+ #ifdef HAVE_MPI
+ Mat getPetscDistributed(const DistributedContext<T_Scalar>& distributedContext) const;
+ void modifyPetscORAS(Mat mat, const std::vector< unsigned long long > &localIdOfOwned, const std::vector< unsigned long long > &localIdOfNonOwned, const std::vector< unsigned long long >& localToGlobalMap, const IS row, const IS col) const;
+ #endif
#endif
common::Memory memory() const;
diff --git a/src/system/Solver.cpp b/src/system/Solver.cpp
index ae9dbe1f3e78fd13aa07588820c7b7e160a7b9ac..665b2f25e9d98917b35389ec5a3500af828df44a 100644
--- a/src/system/Solver.cpp
+++ b/src/system/Solver.cpp
@@ -25,6 +25,49 @@
namespace gmshfem::system
{
+ // Code to use ASM preconditioner
+ template< class T_Scalar >
+ struct ModifySubmatrixContext {
+ std::vector< PetscInt > *ghostsIndices;
+ size_t nLocalRows;
+ const std::vector< unsigned long long > *localIdOfNonOwned;
+ const std::vector< unsigned long long > *localIdOfOwned;
+ const std::vector< unsigned long long > *localToGlobalMap;
+ MatrixFactory< T_Scalar > *A;
+ };
+
+ template< class T_Scalar >
+ static PetscErrorCode modifySubMatrix(PC pc, PetscInt nsub, const IS row[], const IS col[], Mat submat[], void *ctxPtr)
+ {
+ ModifySubmatrixContext< T_Scalar > *ctx = static_cast< ModifySubmatrixContext< T_Scalar > * >(ctxPtr);
+ if(ctx == nullptr)
+ return PETSC_ERR_ARG_WRONG;
+
+ auto &ghosts = *ctx->ghostsIndices;
+ PetscInt m, n;
+ MatGetSize(submat[0], &m, &n);
+ MatType type;
+ MatGetType(submat[0], &type);
+
+ PetscInt nRows;
+ ISGetSize(row[0], &nRows);
+
+ //ISRestoreIndices(row[0], &rowsIdx);
+
+
+ ctx->A->modifyPetscORAS(submat[0], *ctx->localIdOfOwned, *ctx->localIdOfNonOwned, *ctx->localToGlobalMap, row[0], col[0]);
+
+ return PETSC_SUCCESS;
+ };
+
+ template< class T_Scalar >
+ void Solver< T_Scalar >::setMonitor(std::unique_ptr< system::AbstractMonitor > &&ptr)
+ {
+ _customMonitor = std::move(ptr);
+ }
+
+
+
template< class T_Scalar >
Solver< T_Scalar >::Solver(MatrixFactory< T_Scalar > *const A, std::vector< VectorFactory< T_Scalar > > &b) :
@@ -127,6 +170,11 @@ namespace gmshfem::system
// use direct sparse solver by default
KSPSetType(_ksp, "preonly");
PCSetType(pc, (_isCholesky< T_Scalar >(_A->getOptions().symmetric(), _A->getOptions().hermitian()) ? "cholesky" : "lu"));
+
+ PetscViewerAndFormat *vf;
+ //PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD), PETSC_VIEWER_DEFAULT, &vf);
+ //KSPMonitorSet(_ksp, (PetscErrorCode(*)(KSP, PetscInt, PetscReal, void *))KSPMonitorResidual, vf, (PetscErrorCode(*)(void **))PetscViewerAndFormatDestroy);
+
#if defined(PETSC_HAVE_MUMPS)
PCFactorSetMatSolverType(pc, "mumps");
#elif defined(PETSC_HAVE_MKL_PARDISO)
@@ -181,6 +229,210 @@ namespace gmshfem::system
}
}
+ template< class T_Scalar >
+ void Solver< T_Scalar >::solveDistributed(std::vector< T_Scalar > &values, system::DISTRIBUTED_SOLVE_TYPE solveType, const DistributedContext<T_Scalar>& distributedContext, const bool reusePreconditioner)
+ {
+ // Unpack distributed Context
+ const std::vector< unsigned long long > &localIdOfOwned = distributedContext.localIDofOwned();
+ const std::vector< unsigned long long > &localIdOfNonOwned = distributedContext.localIDofNonOwned();
+ const std::vector< unsigned long long > &localToGlobalMap = distributedContext.localToGlobal();
+
+ //*******************
+ // Solve :
+ // A x = b for rhs index 0
+ //*******************
+ try {
+#ifdef HAVE_PETSC
+ if(_A->getModule()->name() != "A" && _A->getModule()->name() != "AFrequency") {
+ throw common::Exception("A separately assembled system can not be solved");
+ }
+
+ PC pc;
+ KSPSetReusePreconditioner(_kspDistributed, reusePreconditioner ? PETSC_TRUE : PETSC_FALSE);
+
+ Mat A, B;
+ KSPGetOperators(_kspDistributed, &A, &B);
+ PetscInt m, n;
+ MatGetSize(A, &m, &n);
+ if((scalar::IsComplex< T_Scalar >::value == true && scalar::IsComplex< PetscScalar >::value == false ? 2 * _A->size() : _A->size()) != unsigned(m)) {
+ KSPSetOperators(_kspDistributed, NULL, NULL);
+ }
+
+ if((_A->size() != static_cast< unsigned long long >(m) || _A->size() != static_cast< unsigned long long >(n)) || !reusePreconditioner) {
+ Mat mat = _A->getPetscDistributed(distributedContext);
+ KSPSetOperators(_kspDistributed, mat, mat);
+ MatDestroy(&mat);
+ }
+ KSPGetPC(_kspDistributed, &pc);
+
+ std::string directSolver = "petsc";
+#if defined(PETSC_HAVE_MUMPS)
+ msg::debug << "Using MUMPS." << msg::endl;
+ directSolver = "mumps";
+#elif defined(PETSC_HAVE_MKL_PARDISO)
+ directSolver = "mkl_pardiso";
+#elif defined(PETSC_HAVE_UMFPACK) || defined(PETSC_HAVE_SUITESPARSE)
+ directSolver = "umfpack";
+#endif
+
+
+ bool useDDM = false;
+ switch(solveType)
+ {
+ case DISTRIBUTED_SOLVE_TYPE::DIRECT:
+ KSPSetType(_kspDistributed, "preonly");
+ PCSetType(pc, "lu");
+ PCFactorSetMatSolverType(pc, directSolver.c_str());
+ break;
+ case DISTRIBUTED_SOLVE_TYPE::RAS:
+ KSPSetType(_kspDistributed, "gmres");
+ PCSetType(pc, "asm");
+ useDDM = true;
+ break;
+ case DISTRIBUTED_SOLVE_TYPE::ORAS:
+ KSPSetType(_kspDistributed, "gmres");
+ PCSetType(pc, "asm");
+ useDDM = true;
+ break;
+
+ }
+ //KSPSetType(_kspDistributed, "gmres");
+ //PCSetType(pc, "asm");
+ //IS * loc, *glob;
+ //PetscInt nsub;
+
+ //KSPSetType(_kspDistributed, "preonly");
+ //PCSetType(pc, (_isCholesky< T_Scalar >(_A->getOptions().symmetric(), _A->getOptions().hermitian()) ? "cholesky" : "lu"));
+
+ PetscViewerAndFormat *vf;
+ PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD), PETSC_VIEWER_DEFAULT, &vf);
+ // Montior res
+ KSPMonitorSet(_kspDistributed, (PetscErrorCode(*)(KSP, PetscInt, PetscReal, void *))KSPMonitorResidual, vf, (PetscErrorCode(*)(void **))PetscViewerAndFormatDestroy);
+ if (_customMonitor) {
+ KSPMonitorSet(_kspDistributed, [](KSP ksp, PetscInt it, PetscReal rn, void * ctx)
+ {system::AbstractMonitor* monitor= static_cast<system::AbstractMonitor*> (ctx); return monitor->Monitor(ksp, it, rn);}, _customMonitor.get(), nullptr);
+ }
+
+ KSPAppendOptionsPrefix(_kspDistributed, "fem_");
+ KSPSetFromOptions(_kspDistributed);
+
+ Vec sol;
+ Vec b = _b.at(0).getPetscDistributed(distributedContext);
+ VecScale(b, -1.0);
+
+ // Create ghosted vector
+ std::vector< PetscInt > ghostsIndices;
+ for(auto sharedLocalIdx : localIdOfNonOwned) {
+ ghostsIndices.push_back(localToGlobalMap[sharedLocalIdx - 1] - 1);
+ }
+ VecCreateGhost(PETSC_COMM_WORLD, localIdOfOwned.size(), PETSC_DETERMINE, ghostsIndices.size(), ghostsIndices.data(), &sol);
+ msg::debug << "Local vector with ghosts : " << localIdOfOwned.size() << " owned DOFs and " << ghostsIndices.size() << " ghosts." << msg::endl;
+ PetscInt low, high;
+ VecGetOwnershipRange(sol, &low, &high);
+ PetscInt size;
+ VecGetLocalSize(sol, &size);
+ std::vector<PetscInt> isdata(localIdOfOwned.size() + ghostsIndices.size());
+ for (unsigned k = 0; k < localIdOfOwned.size(); ++k){
+ isdata[k] = low + k;
+ }
+ for (unsigned k = 0; k < ghostsIndices.size(); ++k){
+ isdata[k + localIdOfOwned.size()] = ghostsIndices[k];
+ }
+
+
+ IS mySubdomain;
+ if(useDDM && !reusePreconditioner) {
+ ISCreateGeneral(PETSC_COMM_WORLD, isdata.size(), isdata.data(), PETSC_COPY_VALUES, &mySubdomain);
+ PCASMSetLocalSubdomains(pc, 1, &mySubdomain, nullptr);
+ }
+
+
+ KSPType ksptype;
+ KSPGetType(_kspDistributed, &ksptype);
+ PCType pctype;
+ PCGetType(pc, &pctype);
+ MatSolverType stype;
+ PCFactorGetMatSolverType(pc, &stype);
+ if(common::GmshFem::getMPIRank() == 0) {
+ msg::info << "System size: " << size << " - "
+ << (ksptype ? ksptype : "none") << " "
+ << (pctype ? pctype : "none") << " "
+ << (stype ? stype : "none") << msg::endl;
+ }
+
+
+ ModifySubmatrixContext<T_Scalar> ctx;
+ ctx.ghostsIndices = &ghostsIndices;
+ ctx.nLocalRows = localIdOfOwned.size();
+ ctx.localIdOfNonOwned = &localIdOfNonOwned;
+ ctx.localIdOfOwned = &localIdOfOwned;
+ ctx.localToGlobalMap = &localToGlobalMap;
+ ctx.A = _A;
+
+
+
+ if (useDDM && solveType == DISTRIBUTED_SOLVE_TYPE::ORAS)
+ PCSetModifySubMatrices(pc, modifySubMatrix<T_Scalar>, &ctx);
+
+
+ KSPSetUp(_kspDistributed);
+ // SUB KSP
+ PetscInt nloc, firstloc;
+ KSP *subksp;
+ if(useDDM) {
+ PCASMGetSubKSP(pc, &nloc, &firstloc, &subksp);
+ msg::debug << "nloc: " << nloc << ", firstloc: " << firstloc << msg::endl;
+ PC subpc;
+ //subksp[0];
+ KSPGetType(subksp[0], &ksptype);
+ KSPGetPC(subksp[0], &subpc);
+ PCSetType(subpc, PCLU);
+ PCGetType(subpc, &pctype);
+ PCFactorSetMatSolverType(subpc, "mumps");
+ PCFactorGetMatSolverType(subpc, &stype);
+ if (gmshfem::common::GmshFem::getMPIRank() == 0)
+ msg::info << "Solving subdomains using : " << ksptype << ", " << (pctype ? pctype : "none") << ", " << (stype ? stype : "none") << msg::endl;
+ }
+
+
+ KSPSolve(_kspDistributed, b, sol);
+ // VISUALIZE SUB
+ /*PCASMGetLocalSubdomains(pc, &nsub, &loc, &glob);
+ msg::print << "Num subdomains:" << nsub << msg::endl;
+ int isSize;
+ ISGetLocalSize(loc[0], &isSize);
+ msg::print << "Local size:" << isSize << msg::endl;
+ */
+
+ // Visualize solution
+ VecGhostUpdateBegin(sol, INSERT_VALUES, SCATTER_FORWARD);
+ VecGhostUpdateEnd(sol, INSERT_VALUES, SCATTER_FORWARD);
+ Vec x_local;
+ VecGhostGetLocalForm(sol, &x_local);
+
+
+ VecGetLocalSize(x_local, &size);
+
+ values.resize((scalar::IsComplex< T_Scalar >::value == true && scalar::IsComplex< PetscScalar >::value == false ? size / 2 : size));
+ distributedContext.readScatteredData(values, sol);
+ //PetscScalar *array;
+ //VecGetArray(x_local, &array);
+ //PetscInterface< T_Scalar, PetscScalar > interface;
+ //interface(values, array, size);
+ //VecRestoreArray(x_local, &array);
+
+ VecDestroy(&sol);
+ VecDestroy(&b);
+#else
+ msg::error << "Unable to solve the linear system without PETSc" << msg::endl;
+#endif
+ }
+ catch(const std::exception &exc) {
+ values.clear();
+ throw;
+ }
+ }
+
template< class T_Scalar >
void Solver< T_Scalar >::solveAll(std::vector< std::vector< T_Scalar > > &values, const bool reusePreconditioner)
{
@@ -303,6 +555,219 @@ namespace gmshfem::system
#endif
}
+ template< class T_Scalar >
+ void Solver< T_Scalar >::solveAllDistributed(std::vector< std::vector< T_Scalar > > &values, DISTRIBUTED_SOLVE_TYPE solveType, const DistributedContext<T_Scalar>& distributedContext, const bool reusePreconditioner)
+ {
+ //*******************
+// Solve with MPI:
+// A X = B for a thin rectangular matrix B (multi RHS)
+//*******************
+
+ // Unpack distributed Context
+ const std::vector< unsigned long long > &localIdOfOwned = distributedContext.localIDofOwned();
+ const std::vector< unsigned long long > &localIdOfNonOwned = distributedContext.localIDofNonOwned();
+ const std::vector< unsigned long long > &localToGlobalMap = distributedContext.localToGlobal();
+
+#if !defined(PETSC_HAVE_HPDDM) || !defined(HAVE_PETSC)
+ throw gmshfem::common::Exception("Solving of multiple RHS requires PETSc with HPDDM");
+#else
+
+
+ try {
+ if(_A->getModule()->name() != "A" && _A->getModule()->name() != "AFrequency") {
+ throw common::Exception("A separately assembled system can not be solved");
+ }
+
+ PC pc;
+ KSPSetReusePreconditioner(_kspDistributed, reusePreconditioner ? PETSC_TRUE : PETSC_FALSE);
+
+ Mat A, B;
+ KSPGetOperators(_kspDistributed, &A, &B);
+ PetscInt m, n, mloc, nloc;
+ MatGetSize(A, &m, &n);
+ if((scalar::IsComplex< T_Scalar >::value == true && scalar::IsComplex< PetscScalar >::value == false ? 2 * _A->size() : _A->size()) != unsigned(m)) {
+ KSPSetOperators(_kspDistributed, NULL, NULL);
+ }
+
+ if((_A->size() != static_cast< unsigned long long >(m) || _A->size() != static_cast< unsigned long long >(n)) || !reusePreconditioner) {
+
+ Mat mat = _A->getPetscDistributed(distributedContext);
+ KSPSetOperators(_kspDistributed, mat, mat);
+ MatGetSize(mat, &m, &n);
+ MatGetLocalSize(mat, &mloc, &nloc);
+ MatDestroy(&mat);
+ }
+ KSPGetPC(_kspDistributed, &pc);
+
+ std::string directSolver = "petsc";
+#if defined(PETSC_HAVE_MUMPS)
+ msg::debug << "Using MUMPS." << msg::endl;
+ directSolver = "mumps";
+#elif defined(PETSC_HAVE_MKL_PARDISO)
+ directSolver = "mkl_pardiso";
+#elif defined(PETSC_HAVE_UMFPACK) || defined(PETSC_HAVE_SUITESPARSE)
+ directSolver = "umfpack";
+#endif
+
+
+ KSPSetType(_kspDistributed, KSPHPDDM);
+
+
+ bool useDDM = false;
+ switch(solveType)
+ {
+ case DISTRIBUTED_SOLVE_TYPE::DIRECT:
+ KSPHPDDMSetType(_kspDistributed, KSP_HPDDM_TYPE_PREONLY);
+ PCSetType(pc, "lu");
+ PCFactorSetMatSolverType(pc, directSolver.c_str());
+ break;
+ case DISTRIBUTED_SOLVE_TYPE::RAS:
+ KSPHPDDMSetType(_kspDistributed, KSP_HPDDM_TYPE_BGMRES);
+ PCSetType(pc, "asm");
+ useDDM = true;
+ break;
+ case DISTRIBUTED_SOLVE_TYPE::ORAS:
+ KSPHPDDMSetType(_kspDistributed, KSP_HPDDM_TYPE_BGMRES);
+ PCSetType(pc, "asm");
+ useDDM = true;
+ break;
+
+ }
+
+
+ KSPAppendOptionsPrefix(_kspDistributed, "fem_");
+ KSPSetFromOptions(_kspDistributed);
+ PetscViewerAndFormat *vf;
+ PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD), PETSC_VIEWER_DEFAULT, &vf);
+ KSPMonitorSet(_kspDistributed, (PetscErrorCode(*)(KSP, PetscInt, PetscReal, void *))KSPMonitorResidual, vf, (PetscErrorCode(*)(void **))PetscViewerAndFormatDestroy);
+
+
+ // Define block matrices
+ unsigned nRHS = _b.size();
+ Mat BlockRHS, BlockY;
+ BlockRHS = VectorFactory< T_Scalar >::generateDistributedBlockRHS(_b, distributedContext);
+ MatDuplicate(BlockRHS, MAT_DO_NOT_COPY_VALUES, &BlockY);
+ MatScale(BlockRHS, -1);
+
+
+ // If DDM is used, define subdomains
+ std::vector< PetscInt > ghostsIndices;
+ for(auto sharedLocalIdx : localIdOfNonOwned) {
+ ghostsIndices.push_back(localToGlobalMap[sharedLocalIdx - 1] - 1);
+ }
+
+ PetscInt low, high;
+ MatGetOwnershipRange(BlockRHS, &low, &high);
+ PetscInt size;
+ //VecGetLocalSize(sol, &size);
+ std::vector<PetscInt> isdata(localIdOfOwned.size() + ghostsIndices.size());
+ for (unsigned k = 0; k < localIdOfOwned.size(); ++k){
+ isdata[k] = low + k;
+ }
+ for (unsigned k = 0; k < ghostsIndices.size(); ++k){
+ isdata[k + localIdOfOwned.size()] = ghostsIndices[k];
+ }
+
+ IS mySubdomain;
+ if(useDDM && !reusePreconditioner) {
+ ISCreateGeneral(PETSC_COMM_WORLD, isdata.size(), isdata.data(), PETSC_COPY_VALUES, &mySubdomain);
+ PCASMSetLocalSubdomains(pc, 1, &mySubdomain, nullptr);
+ }
+
+ // Setup DDM and subsolvers
+ ModifySubmatrixContext<T_Scalar> ctx;
+ ctx.ghostsIndices = &ghostsIndices;
+ ctx.nLocalRows = localIdOfOwned.size();
+ ctx.localIdOfNonOwned = &localIdOfNonOwned;
+ ctx.localIdOfOwned = &localIdOfOwned;
+ ctx.localToGlobalMap = &localToGlobalMap;
+ ctx.A = _A;
+
+
+ if (useDDM && solveType == DISTRIBUTED_SOLVE_TYPE::ORAS)
+ PCSetModifySubMatrices(pc, modifySubMatrix<T_Scalar>, &ctx);
+
+ KSPSetUp(_kspDistributed);
+ // SUB KSP
+ PetscInt nLocalSubdomains, firstloc;
+ KSP *subksp;
+ if(useDDM) {
+ PCASMGetSubKSP(pc, &nLocalSubdomains, &firstloc, &subksp);
+ msg::debug << "nloc: " << nLocalSubdomains << ", firstloc: " << firstloc << msg::endl;
+ PC subpc;
+ KSPType ksptype;
+ PCType pctype;
+ MatSolverType stype;
+ //subksp[0];
+ KSPGetType(subksp[0], &ksptype);
+ KSPGetPC(subksp[0], &subpc);
+ PCSetType(subpc, PCLU);
+ PCGetType(subpc, &pctype);
+ PCFactorSetMatSolverType(subpc, "mumps");
+ PCFactorGetMatSolverType(subpc, &stype);
+ if(gmshfem::common::GmshFem::getMPIRank() == 0)
+ msg::info << "Solving subdomains using : " << ksptype << ", " << (pctype ? pctype : "none") << ", " << (stype ? stype : "none") << msg::endl;
+ }
+
+
+ KSPType ksptype;
+ KSPGetType(_kspDistributed, &ksptype);
+ PCType pctype;
+ PCGetType(pc, &pctype);
+ MatSolverType stype;
+ PCFactorGetMatSolverType(pc, &stype);
+ if(common::GmshFem::getMPIRank() == 0) {
+ msg::info << "System size: " << m << " - "
+ << (ksptype ? ksptype : "none") << " "
+ << (pctype ? pctype : "none") << " "
+ << (stype ? stype : "none") << msg::endl;
+ }
+
+ KSPMatSolve(_kspDistributed, BlockRHS, BlockY);
+
+#warning "Handle FEM complex + real PETSc ?"
+ //values.resize((scalar::IsComplex< T_Scalar >::value == true && scalar::IsComplex< PetscScalar >::value == false ? size / 2 : size));
+
+
+ for (unsigned c = 0; c < nRHS; ++c) {
+ values[c].resize(mloc + ghostsIndices.size());
+ Vec col;
+ MatDenseGetColumnVecRead(BlockY, c, &col);
+ distributedContext.readScatteredData(values[c], col);
+ MatDenseRestoreColumnVecRead(BlockY, c, &col);
+}
+
+ // Create a single Vec that contains the block solution, column major format.
+ /*PetscScalar *array;
+ MatDenseGetArray(BlockY, &array);
+ PetscInterface< T_Scalar, PetscScalar > interface;
+ //MatView(BlockY, PETSC_VIEWER_STDOUT_WORLD);
+
+ for(unsigned iRHS = 0; iRHS < nRHS; ++iRHS) {
+ // Get the pointer to the beginning of the iRHS-th column.
+ PetscScalar *col_start = array + iRHS * mloc;
+
+
+ // TODO: have value read ghosts
+ values[iRHS].resize(mloc);
+ interface(values[iRHS], col_start, mloc);
+ }
+
+ MatDenseRestoreArray(BlockY, &array);*/
+
+ MatDestroy(&BlockRHS);
+ MatDestroy(&BlockY);
+
+ // Don't set fields to the solved system; loadSolution() should be used.
+ }
+
+ catch(const std::exception &exc) {
+ values.clear();
+ throw;
+ }
+#endif
+ }
+
template< class T_Scalar >
void Solver< T_Scalar >::eigensolve(std::vector< scalar::ComplexPrecision< T_Scalar > > &eigenvalues, std::vector< std::vector< scalar::ComplexPrecision< T_Scalar > > > &eigenvectors, const bool computeEigenvectors, const unsigned int numberOfEigenvalues, const scalar::ComplexPrecision< T_Scalar > target)
{
diff --git a/src/system/VectorFactory.cpp b/src/system/VectorFactory.cpp
index ca568cf56e2addc839d417bec2d623c16c4f24a6..239e0dd0408403321c696f631335348cecf52e5f 100644
--- a/src/system/VectorFactory.cpp
+++ b/src/system/VectorFactory.cpp
@@ -47,6 +47,28 @@ namespace gmshfem::system
_vec.resize(_size, 0.);
}
+ template< class T_Scalar >
+ void VectorFactory< T_Scalar >::save(const std::string &name) const
+ {
+ common::CSVio file(name);
+ if(!file.isOpen()) {
+ throw common::Exception("Cannot open file " + name + ".csv");
+ }
+
+ file << csv::precision(scalar::PrecisionDigits< T_Scalar >::value) << csv::scientific;
+ if(scalar::IsComplex< T_Scalar >::value) {
+ for(auto i = 0ULL; i < _vec.size(); ++i) {
+ file << std::real(_vec[i]) << std::imag(_vec[i]) << csv::endl;
+ }
+ }
+ else {
+ for(auto i = 0ULL; i < _vec.size(); ++i) {
+ file << _vec[i] << csv::endl;
+ }
+ }
+ file.close();
+ }
+
template< class T_Scalar >
void VectorFactory< T_Scalar >::removeVector()
{
@@ -268,6 +290,60 @@ namespace gmshfem::system
interface(_size, &_vec[0], &vecPetsc, _shouldBeDestroyedWithPetsc);
return vecPetsc;
}
+
+#ifdef HAVE_MPI
+ template< class T_Scalar >
+ Vec VectorFactory< T_Scalar >::getPetscDistributed(const DistributedContext<T_Scalar>& distributedContext) const
+ {
+ PetscScalar* data = static_cast<PetscScalar*>(generateRawLocalData(distributedContext));
+ unsigned nLocDofs = distributedContext.localIDofOwned().size();
+
+ Vec vecPetsc;
+ VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, nLocDofs, PETSC_DETERMINE, data, &vecPetsc);
+ _shouldBeDestroyedWithPetsc.push_back(data);
+
+ return vecPetsc;
+ }
+ template< class T_Scalar >
+ void *VectorFactory< T_Scalar >::generateRawLocalData(const DistributedContext<T_Scalar>& distributedContext) const
+ {
+ const auto& localDofs = distributedContext.localIDofOwned();
+ unsigned nLocalRows = localDofs.size();
+ // In the distributed vector, copy the owned DOFs
+ //std::vector<PetscScalar> data(nLocalRows);
+ PetscScalar* data = (PetscScalar *)std::malloc(nLocalRows * sizeof(PetscScalar));
+ for (unsigned i = 0; i < nLocalRows; ++i) {
+ //msg::print << "vec: " << i << ", " << localDofs[i] << ", " << localToGlobal[i] << msg::endl;
+ data[i] = _vec.at(localDofs[i]-1);
+ }
+
+ return data; // Must be deleted by caller with std::free
+ }
+ template< class T_Scalar >
+ Mat VectorFactory< T_Scalar >::generateDistributedBlockRHS(const std::vector< VectorFactory< T_Scalar > > &allRHS, const DistributedContext<T_Scalar>& distributedContext)
+ {
+ unsigned nLocDofs = distributedContext.localIDofOwned().size();
+ unsigned nRHS = allRHS.size();
+
+ Mat blockRHS;
+ MatCreateDense(PETSC_COMM_WORLD, nLocDofs, PETSC_DECIDE, PETSC_DECIDE, nRHS, nullptr, &blockRHS);
+ PetscScalar *array;
+ MatDenseGetArray(blockRHS, &array);
+ PetscInterface< T_Scalar, PetscScalar > interface;
+
+ for(unsigned iRHS = 0; iRHS < nRHS; ++iRHS) {
+ // Get the pointer to the beginning of the iRHS-th column.
+ PetscScalar *col_start = array + iRHS * nLocDofs;
+ PetscScalar *colData = static_cast< PetscScalar * >(allRHS[iRHS].generateRawLocalData(distributedContext));
+ std::copy(colData, colData + nLocDofs, col_start);
+ std::free(colData);
+ }
+
+ MatDenseRestoreArray(blockRHS, &array);
+
+ return blockRHS;
+ }
+#endif
#endif
template< class T_Scalar >
diff --git a/src/system/VectorFactory.h b/src/system/VectorFactory.h
index 215b567eb37767c2705ee0c6a66be90462233506..e27d5013b4a8440ba78d75180b7c174b0e5fce70 100644
--- a/src/system/VectorFactory.h
+++ b/src/system/VectorFactory.h
@@ -40,6 +40,7 @@ namespace gmshfem::system
~VectorFactory();
void init(const unsigned long long size);
+ void save(const std::string& name) const;
void removeVector();
void removePetscData();
@@ -56,6 +57,13 @@ namespace gmshfem::system
#ifdef HAVE_PETSC
Vec getPetsc() const;
+ #ifdef HAVE_MPI
+ Vec getPetscDistributed(const DistributedContext<T_Scalar>& distributedContext) const;
+ // Return PetscScalar* but not declared yet
+ void* generateRawLocalData(const DistributedContext<T_Scalar>& distributedContext) const;
+
+ static Mat generateDistributedBlockRHS(const std::vector<VectorFactory<T_Scalar>>& allRHS, const DistributedContext<T_Scalar>& distributedContext);
+ #endif
#endif
common::Memory memory() const;
unsigned long long size() const;
diff --git a/unitTests/correctness/CMakeLists.txt b/unitTests/correctness/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..aad68ce18c83b60d1b323b51852be434960ba31c
--- /dev/null
+++ b/unitTests/correctness/CMakeLists.txt
@@ -0,0 +1,39 @@
+
+# List of all test cases for iteration
+set(TEST_CASES
+ distributedPoisson1DDirichlet
+ distributedPoisson1DLagrangeMultiplier
+ distributedPoisson2DTwoRanksDirect
+ distributedPoisson2DORAS
+)
+
+# Define test cases with associated ranks
+set(distributedPoisson1DDirichlet_RANKS 4)
+set(distributedPoisson1DLagrangeMultiplier_RANKS 2)
+set(distributedPoisson2DTwoRanksDirect_RANKS 2)
+set(distributedPoisson2DORAS_RANKS 4)
+
+foreach(TUTO IN LISTS TEST_CASES)
+ # Retrieve the number of ranks for this test case
+ set(NRANKS_VAR "${TUTO}_RANKS")
+ set(NRANKS "${${NRANKS_VAR}}")
+
+ # Define test commands
+ set(TEST_DIR "${CMAKE_CURRENT_SOURCE_DIR}/${TUTO}/build")
+
+ # Ensure the directory is created
+ add_test(NAME "correctness:${TUTO}:mkdir" COMMAND "${CMAKE_COMMAND}" -E make_directory "${TEST_DIR}")
+
+ # Configure the test
+ add_test(NAME "correctness:${TUTO}:cmake" COMMAND "${CMAKE_COMMAND}" -DCMAKE_PREFIX_PATH=${CMAKE_INSTALL_PREFIX} .. WORKING_DIRECTORY "${TEST_DIR}")
+
+ # Build the test
+ add_test(NAME "correctness:${TUTO}:build" COMMAND "${CMAKE_MAKE_PROGRAM}" WORKING_DIRECTORY "${TEST_DIR}")
+
+ # Determine if the test is MPI based on the number of ranks
+ if(NOT "${NRANKS}" STREQUAL "1")
+ add_test(NAME "correctness:${TUTO}:run" COMMAND mpirun -n ${NRANKS} ./tuto -maxThreads 1 WORKING_DIRECTORY "${TEST_DIR}")
+ else()
+ add_test(NAME "correctness:${TUTO}:run" COMMAND ./tuto -maxThreads 1 WORKING_DIRECTORY "${TEST_DIR}")
+ endif()
+endforeach()
diff --git a/unitTests/correctness/distributedPoisson1DDirichlet/CMakeLists.txt b/unitTests/correctness/distributedPoisson1DDirichlet/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3ed63f13623e9d7ed45010b2a7ffd418461a2b0c
--- /dev/null
+++ b/unitTests/correctness/distributedPoisson1DDirichlet/CMakeLists.txt
@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
+
+project(tuto CXX)
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/../tutorial.cmake)
+
+add_executable(tuto main.cpp ${EXTRA_INCS})
+target_link_libraries(tuto ${EXTRA_LIBS})
diff --git a/unitTests/correctness/distributedPoisson1DDirichlet/main.cpp b/unitTests/correctness/distributedPoisson1DDirichlet/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..66d8e5c8e95c3097bdc1837a18b866781f7da54d
--- /dev/null
+++ b/unitTests/correctness/distributedPoisson1DDirichlet/main.cpp
@@ -0,0 +1,153 @@
+#include <algorithm>
+#include <gmsh.h>
+#include <gmshfem/Formulation.h>
+#include <gmshfem/GmshFem.h>
+#include <gmshfem/DistributedField.h>
+#include <gmshfem/MatrixFactory.h>
+#include <gmshfem/VectorFactory.h>
+#include <chrono>
+#include <thread>
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+using namespace gmshfem::problem;
+using namespace gmshfem::domain;
+using namespace gmshfem::field;
+using namespace gmshfem::function;
+using namespace gmshfem::post;
+using namespace gmshfem::equation;
+using namespace gmshfem::msg;
+
+void genMesh(double lc, unsigned nDom, unsigned nOverlap)
+{
+ namespace geo = gmsh::model::geo;
+
+ gmsh::model::add("genMesh");
+ gmsh::option::setNumber("General.ExpertMode", 1);
+ double L = 2;
+ double domSize = L / nDom;
+ double delta = nOverlap * lc;
+ if (delta > domSize)
+ {
+ throw common::Exception("Subdomains too small!");
+ }
+
+ std::vector<int> points;
+ points.reserve(nDom * 3 + 1);
+ points.push_back(geo::addPoint(0, 0, 0, lc));
+ for (unsigned iDom = 0; iDom < nDom; ++iDom)
+ {
+ double x = iDom * domSize;
+ points.push_back(geo::addPoint(x + delta, 0, 0, lc));
+ points.push_back(geo::addPoint(x + domSize - delta, 0, 0, lc));
+ points.push_back(geo::addPoint(x + domSize, 0, 0, lc));
+ }
+
+ std::vector<int> lines;
+ lines.reserve(3 * nDom);
+ for (unsigned iDom = 0; iDom < nDom; ++iDom)
+ {
+ lines.push_back(geo::addLine(points[iDom * 3], points[iDom * 3 + 1]));
+ lines.push_back(geo::addLine(points[iDom * 3 + 1], points[iDom * 3 + 2]));
+ lines.push_back(geo::addLine(points[iDom * 3 + 2], points[iDom * 3 + 3]));
+ }
+
+ gmsh::model::geo::synchronize();
+ gmsh::model::mesh::generate();
+
+ for (unsigned iDom = 0; iDom < nDom; ++iDom)
+ {
+ gmsh::model::addPhysicalGroup(1, {lines[iDom*3], lines[iDom*3+1], lines[iDom*3+2]}, 100 * iDom + 1, "omega");
+ std::vector<int> overlap;
+ if (iDom > 0)
+ overlap.push_back(lines[3*iDom-1]);
+ if (iDom != (nDom - 1))
+ overlap.push_back(lines[3*iDom+3]);
+ gmsh::model::addPhysicalGroup(1, {overlap}, 100 * iDom + 2, "overlap");
+ if (iDom == 0)
+ gmsh::model::addPhysicalGroup(0, {points[0]}, iDom * 100 + 3, "boundary");
+ if (iDom == (nDom - 1))
+ gmsh::model::addPhysicalGroup(0, {points.back()}, iDom * 100 + 4, "boundary");
+ gmsh::model::geo::synchronize();
+ }
+
+ // Everyone needs the full mesh (so far)
+ gmsh::model::addPhysicalGroup(1, {lines}, 1234, "omega_all");
+ gmsh::write("mesh_full.msh");
+
+ gmsh::model::remove();
+}
+
+int main(int argc, char **argv)
+{
+ using namespace std::chrono_literals;
+
+ GmshFem gmshFem(argc, argv);
+ unsigned int MPI_Rank = gmshFem.getMPIRank();
+ unsigned MPI_Size = gmshFem.getMPISize();
+ //if (MPI_Size != 2)
+ // throw Exception("Must run with 2 ranks! (Got " + std::to_string(MPI_Size) + ")");
+
+ int order = 1;
+ gmshFem.userDefinedParameter(order, "order");
+
+ double lc = 2.0 / MPI_Size / 3; // 3 elements per domain
+ gmshFem.userDefinedParameter(lc, "lc");
+
+ /****
+ * Gmsh part
+ ****/
+ if (MPI_Rank == 0)
+ genMesh(lc, MPI_Size, 1);
+
+ std::this_thread::sleep_for(100ms);
+ bool volumeSrc = false;
+ gmshFem.userDefinedParameter(volumeSrc, "volumeSrc");
+
+ // -X²/2 + X has laplacian -1, is 0 on 0 and 2 on 2
+ // function::ScalarFunction<std::complex<double>> analytical = -x<std::complex<double>>()*x<std::complex<double>>()/2 + 2*x<std::complex<double>>();
+ auto analytical = real(-x<std::complex<double>>() * x<std::complex<double>>() / 2 + 2 * x<std::complex<double>>());
+ if (!volumeSrc)
+ {
+ analytical = 1 * x<std::complex<double>>();
+ }
+ gmsh::open("mesh_full.msh");
+ gmsh::model::mesh::createEdges();
+
+ Formulation<std::complex<double>> formulation("formulation_" + std::to_string(MPI_Rank));
+ Domain omega(1, MPI_Rank * 100 + 1);
+ Domain overlap(1, MPI_Rank * 100 + 2);
+ Domain boundary;
+
+ DistributedField<std::complex<double>, Form::Form0> v("v", omega, overlap, FunctionSpaceTypeForm0::HierarchicalH1, order);
+ if (MPI_Rank == 0 || MPI_Rank == (MPI_Size - 1))
+ {
+ int localIdx = MPI_Rank == 0 ? 3 : 4;
+ boundary = Domain(0, MPI_Rank * 100 + localIdx);
+ v.addConstraint(boundary, x<std::complex<double>>());
+ }
+
+ formulation.integral(grad(dof(v)), grad(tf(v)), v.domain(), "Gauss6");
+ if (volumeSrc)
+ formulation.integral(-1, tf(v), v.domain(), "Gauss6");
+
+ formulation.pre();
+ formulation.preDistributed();
+ formulation.assemble();
+ gmshFem.BarrierMPI();
+
+ formulation.solveDistributed(system::DISTRIBUTED_SOLVE_TYPE::DIRECT, false);
+ save(v, omega, "v_" + std::to_string(MPI_Rank));
+
+ auto error1 = real(integrate(pow(real(v - analytical), 2), omega, "Gauss6")) / real(integrate(pow(real(analytical), 2), omega, "Gauss6"));
+ msg::print << "L2 relative error in domain " << MPI_Rank << " : " << error1 << msg::endl;
+ save(analytical, omega, "v_truth_" + std::to_string(MPI_Rank));
+ save(analytical - v, omega, "v_err_" + std::to_string(MPI_Rank));
+
+ if (abs(error1) > 1e-6)
+ {
+ throw common::Exception("Error above tolerance. Log10 is : " + std::to_string(log10(error1)));
+ }
+
+ return 0;
+}
diff --git a/unitTests/correctness/distributedPoisson1DLagrangeMultiplier/CMakeLists.txt b/unitTests/correctness/distributedPoisson1DLagrangeMultiplier/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3ed63f13623e9d7ed45010b2a7ffd418461a2b0c
--- /dev/null
+++ b/unitTests/correctness/distributedPoisson1DLagrangeMultiplier/CMakeLists.txt
@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
+
+project(tuto CXX)
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/../tutorial.cmake)
+
+add_executable(tuto main.cpp ${EXTRA_INCS})
+target_link_libraries(tuto ${EXTRA_LIBS})
diff --git a/unitTests/correctness/distributedPoisson1DLagrangeMultiplier/main.cpp b/unitTests/correctness/distributedPoisson1DLagrangeMultiplier/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..25fc98d1318578601d3d16786942e0ca9693840d
--- /dev/null
+++ b/unitTests/correctness/distributedPoisson1DLagrangeMultiplier/main.cpp
@@ -0,0 +1,151 @@
+#include <algorithm>
+#include <gmsh.h>
+#include <gmshfem/Formulation.h>
+#include <gmshfem/GmshFem.h>
+#include <gmshfem/DistributedField.h>
+#include <gmshfem/MatrixFactory.h>
+#include <gmshfem/VectorFactory.h>
+#include <chrono>
+#include <thread>
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+using namespace gmshfem::problem;
+using namespace gmshfem::domain;
+using namespace gmshfem::field;
+using namespace gmshfem::function;
+using namespace gmshfem::post;
+using namespace gmshfem::equation;
+using namespace gmshfem::msg;
+
+void genMesh(double lc, unsigned nDom, unsigned nOverlap)
+{
+ namespace geo = gmsh::model::geo;
+
+ gmsh::model::add("genMesh");
+ gmsh::option::setNumber("General.ExpertMode", 1);
+ double L = 2;
+ double domSize = L / nDom;
+ double delta = nOverlap * lc;
+ if (delta > domSize)
+ {
+ throw common::Exception("Subdomains too small!");
+ }
+
+ std::vector<int> points;
+ points.reserve(nDom * 3 + 1);
+ points.push_back(geo::addPoint(0, 0, 0, lc));
+ for (unsigned iDom = 0; iDom < nDom; ++iDom)
+ {
+ double x = iDom * domSize;
+ points.push_back(geo::addPoint(x + delta, 0, 0, lc));
+ points.push_back(geo::addPoint(x + domSize - delta, 0, 0, lc));
+ points.push_back(geo::addPoint(x + domSize, 0, 0, lc));
+ }
+
+ std::vector<int> lines;
+ lines.reserve(3 * nDom);
+ for (unsigned iDom = 0; iDom < nDom; ++iDom)
+ {
+ lines.push_back(geo::addLine(points[iDom * 3], points[iDom * 3 + 1]));
+ lines.push_back(geo::addLine(points[iDom * 3 + 1], points[iDom * 3 + 2]));
+ lines.push_back(geo::addLine(points[iDom * 3 + 2], points[iDom * 3 + 3]));
+ }
+
+ gmsh::model::geo::synchronize();
+ gmsh::model::mesh::generate();
+
+ for (unsigned iDom = 0; iDom < nDom; ++iDom)
+ {
+ gmsh::model::addPhysicalGroup(1, {lines[iDom*3], lines[iDom*3+1], lines[iDom*3+2]}, 100 * iDom + 1, "omega");
+ std::vector<int> overlap;
+ if (iDom > 0)
+ overlap.push_back(lines[3*iDom-1]);
+ if (iDom != (nDom - 1))
+ overlap.push_back(lines[3*iDom+3]);
+ gmsh::model::addPhysicalGroup(1, {overlap}, 100 * iDom + 2, "overlap");
+ if (iDom == 0)
+ gmsh::model::addPhysicalGroup(0, {points[0]}, iDom * 100 + 3, "boundary");
+ if (iDom == (nDom - 1))
+ gmsh::model::addPhysicalGroup(0, {points.back()}, iDom * 100 + 4, "boundary");
+ gmsh::model::geo::synchronize();
+ }
+
+ // Everyone needs the full mesh (so far)
+ gmsh::model::addPhysicalGroup(1, {lines}, 1234, "omega_all");
+ gmsh::write("mesh_full.msh");
+
+ gmsh::model::remove();
+}
+
+int main(int argc, char **argv)
+{
+ using namespace std::chrono_literals;
+
+ GmshFem gmshFem(argc, argv);
+ unsigned int MPI_Rank = gmshFem.getMPIRank();
+ unsigned MPI_Size = gmshFem.getMPISize();
+ if (MPI_Size == 1)
+ throw Exception("Must run with multiple ranks! (Got " + std::to_string(MPI_Size) + ")");
+
+ int order = 1;
+ gmshFem.userDefinedParameter(order, "order");
+
+ double lc = 2.0 / MPI_Size / 3; // 3 elements per domain
+ gmshFem.userDefinedParameter(lc, "lc");
+
+ /****
+ * Gmsh part
+ ****/
+ if (MPI_Rank == 0)
+ genMesh(lc, MPI_Size, 1);
+
+ std::this_thread::sleep_for(100ms);
+
+ auto analytical = 1 * x<std::complex<double>>();
+ gmsh::open("mesh_full.msh");
+ gmsh::model::mesh::createEdges();
+
+ Formulation<std::complex<double>> formulation("formulation_" + std::to_string(MPI_Rank));
+ Domain omega(1, MPI_Rank * 100 + 1);
+ Domain overlap(1, MPI_Rank * 100 + 2);
+ Domain boundary;
+
+ DistributedField<std::complex<double>, Form::Form0> v("v", omega, overlap, FunctionSpaceTypeForm0::HierarchicalH1, order);
+ std::unique_ptr<Field<std::complex<double>, Form::Form0>> lambda_ptr; // Optional Lagrange multiplier
+ if (MPI_Rank == 0 || MPI_Rank == (MPI_Size - 1))
+ {
+ int localIdx = MPI_Rank == 0 ? 3 : 4;
+ boundary = Domain(0, MPI_Rank * 100 + localIdx);
+ lambda_ptr = std::make_unique<Field<std::complex<double>, Form::Form0>>("lambda", boundary, FunctionSpaceTypeForm0::Lagrange);
+ // Add inhomogenuous Neumann :
+ formulation.integral(dof(*lambda_ptr), tf(v), boundary, "Gauss0");
+ // Apply constraint
+ formulation.integral(dof(v), tf(*lambda_ptr), boundary, "Gauss0");
+ formulation.integral(-analytical, tf(*lambda_ptr), boundary, "Gauss0");
+
+ }
+
+ formulation.integral(grad(dof(v)), grad(tf(v)), v.domain(), "Gauss6");
+
+
+ formulation.pre();
+ formulation.preDistributed();
+ formulation.assemble();
+ gmshFem.BarrierMPI();
+
+ formulation.solveDistributed(system::DISTRIBUTED_SOLVE_TYPE::DIRECT, false);
+ save(v, omega, "v_" + std::to_string(MPI_Rank));
+
+ auto error1 = real(integrate(pow(real(v - analytical), 2), omega, "Gauss6")) / real(integrate(pow(real(analytical), 2), omega, "Gauss6"));
+ msg::print << "L2 relative error in domain " << MPI_Rank << " : " << error1 << msg::endl;
+ save(analytical, omega, "v_truth_" + std::to_string(MPI_Rank));
+ save(analytical - v, omega, "v_err_" + std::to_string(MPI_Rank));
+
+ if (abs(error1) > 1e-6)
+ {
+ throw common::Exception("Error above tolerance. Log10 is : " + std::to_string(log10(error1)));
+ }
+
+ return 0;
+}
diff --git a/unitTests/correctness/distributedPoisson2DORAS/CMakeLists.txt b/unitTests/correctness/distributedPoisson2DORAS/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3ed63f13623e9d7ed45010b2a7ffd418461a2b0c
--- /dev/null
+++ b/unitTests/correctness/distributedPoisson2DORAS/CMakeLists.txt
@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
+
+project(tuto CXX)
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/../tutorial.cmake)
+
+add_executable(tuto main.cpp ${EXTRA_INCS})
+target_link_libraries(tuto ${EXTRA_LIBS})
diff --git a/unitTests/correctness/distributedPoisson2DORAS/main.cpp b/unitTests/correctness/distributedPoisson2DORAS/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5923c183c0182d5ffca9bdc8ba47c28e7221ee7d
--- /dev/null
+++ b/unitTests/correctness/distributedPoisson2DORAS/main.cpp
@@ -0,0 +1,166 @@
+#include <algorithm>
+#include <gmsh.h>
+#include <gmshfem/Formulation.h>
+#include <gmshfem/GmshFem.h>
+#include <gmshfem/DistributedField.h>
+#include <gmshfem/MatrixFactory.h>
+#include <gmshfem/VectorFactory.h>
+#include <chrono>
+#include <thread>
+#include <gmshfem/CSVio.h>
+using namespace gmshfem;
+using namespace gmshfem::common;
+using namespace gmshfem::problem;
+using namespace gmshfem::domain;
+using namespace gmshfem::field;
+using namespace gmshfem::function;
+using namespace gmshfem::post;
+using namespace gmshfem::equation;
+using namespace gmshfem::msg;
+
+void genMesh()
+{
+ gmsh::model::add("genMesh");
+
+ namespace geo = gmsh::model::geo;
+
+ std::vector<int> points, hlines, vlines, surfaces;
+ std::vector<double> xs{0, 0.4, 0.5, 0.6, 1.0};
+ std::vector<double> ys = xs;
+ for (double y : ys)
+ {
+ for (double x : xs)
+ points.push_back(geo::addPoint(x, y, 0));
+ }
+
+ for (unsigned j = 0; j < 5; ++j)
+ {
+ for (unsigned i = 0; i < 4; ++i)
+ {
+ hlines.push_back(geo::addLine(points[j * 5 + i], points[j * 5 + i + 1]));
+ gmsh::model::geo::mesh::setTransfiniteCurve(hlines.back(), 2);
+
+ }
+ }
+
+ for (unsigned j = 0; j < 4; ++j)
+ {
+ for (unsigned i = 0; i < 5; ++i)
+ {
+ vlines.push_back(geo::addLine(points[j * 5 + i], points[(j+1) * 5 + i]));
+ gmsh::model::geo::mesh::setTransfiniteCurve(vlines.back(), 2);
+
+ }
+ }
+
+
+ for (unsigned j = 0; j < 4; ++j)
+ {
+ for (unsigned i = 0; i < 4; ++i)
+ {
+ int curve = geo::addCurveLoop({hlines[4 * j + i], vlines[5 * j + i + 1], -hlines[4 * (j + 1) + i], -vlines[5 * j + i]});
+ surfaces.push_back(geo::addPlaneSurface({curve}));
+ geo::mesh::setTransfiniteSurface(surfaces.back());
+ //geo::mesh::setRecombine(2, surfaces.back());
+ }
+ }
+
+
+ gmsh::model::geo::synchronize();
+
+ bool makeFourDomains = true;
+ using gmsh::model::addPhysicalGroup;
+ if (makeFourDomains) {
+ // Each subdomain has an "omega" (1), an "overlap" (2), a part of the physical boundary (3) (overlap inc.), and a coupling boundary (4)
+ // ID: above + 100*domain index
+ for (int i = 0; i < 2; ++i) {
+ for (int j = 0; j < 2; ++j) {
+ int domIdx = 2*i+j;
+ int offset = 8*i + 2*j;
+ addPhysicalGroup(2, {surfaces[offset+ 0], surfaces[offset+1], surfaces[offset+4], surfaces[offset+4+1]}, 1 + 100*domIdx);
+ }
+ }
+
+ // Define overlaps
+ addPhysicalGroup(2, {surfaces[2], surfaces[6], surfaces[10], surfaces[9], surfaces[8]}, 2 + 100 * 0);
+ addPhysicalGroup(2, {surfaces[1], surfaces[5], surfaces[9], surfaces[10], surfaces[11]}, 2 + 100 * 1);
+ addPhysicalGroup(2, {surfaces[4], surfaces[5], surfaces[6], surfaces[10], surfaces[14]}, 2 + 100 * 2);
+ addPhysicalGroup(2, {surfaces[5], surfaces[6], surfaces[7], surfaces[9], surfaces[13]}, 2 + 100 * 3);
+
+ // Define physical boundaries (with overlap)
+ addPhysicalGroup(1, {hlines[0], hlines[1], hlines[2], vlines[0], vlines[5], vlines[10]}, 3 + 100 * 0);
+ addPhysicalGroup(1, {hlines[1], hlines[2], hlines[3], vlines[4], vlines[9], vlines[14]}, 3 + 100 * 1);
+ addPhysicalGroup(1, {hlines[16], hlines[17], hlines[18], vlines[5], vlines[10], vlines[15]}, 3 + 100 * 2);
+ addPhysicalGroup(1, {hlines[17], hlines[18], hlines[19], vlines[9], vlines[14], vlines[19]}, 3 + 100 * 3);
+
+ // Define coupling boundaries
+ addPhysicalGroup(1, {hlines[12], hlines[13], hlines[14], vlines[3], vlines[8], vlines[13]}, 4 + 100 * 0);
+ addPhysicalGroup(1, {hlines[13], hlines[14], hlines[15], vlines[1], vlines[6], vlines[11]}, 4 + 100 * 1);
+ addPhysicalGroup(1, {hlines[4], hlines[5], hlines[6], vlines[8], vlines[13], vlines[18]}, 4 + 100 * 2);
+ addPhysicalGroup(1, {hlines[5], hlines[6], hlines[7], vlines[6], vlines[11], vlines[16]}, 4 + 100 * 3);
+ }
+
+ //gmsh::model::addPhysicalGroup(1, hlines);
+ //gmsh::model::addPhysicalGroup(2, surfaces);
+ gmsh::model::mesh::generate();
+
+ gmsh::write("mesh.msh");
+ gmsh::model::remove();
+}
+
+int main(int argc, char **argv)
+{
+ using namespace std::chrono_literals;
+
+ GmshFem gmshFem(argc, argv);
+ unsigned int MPI_Rank = gmshFem.getMPIRank();
+ unsigned MPI_Size = gmshFem.getMPISize();
+ if (MPI_Size != 4)
+ throw Exception("Must run with 4 ranks! (Got " + std::to_string(MPI_Size) + ")");
+
+ /****
+ * Gmsh part
+ ****/
+ if (MPI_Rank == 0)
+ genMesh();
+ gmshFem.BarrierMPI();
+
+ // x + 2y + 1, harmonic fct
+ auto analytical = real(x<std::complex<double>>() + 2 * y<std::complex<double>>() + 1);
+ gmsh::open("mesh.msh");
+ gmsh::model::mesh::createEdges();
+ gmsh::model::mesh::createFaces();
+
+ Formulation<std::complex<double>> formulation("formulation_" + std::to_string(MPI_Rank));
+ //*
+ int shift = MPI_Rank * 100;
+ Domain omega(2, 1 + shift);
+ Domain overlap(2, 2 + shift);
+ Domain boundary(1, 3 + shift);
+ Domain coupling(1, 4 + shift);
+
+ DistributedField<std::complex<double>, Form::Form0> v("v", omega, overlap, FunctionSpaceTypeForm0::HierarchicalH1, 1);
+
+ formulation.integral(grad(dof(v)), grad(tf(v)), v.domain(), "Gauss6");
+ formulation.integral(dof(v), tf(v), coupling, "Gauss6");
+
+
+ v.addConstraint(boundary, analytical);
+
+ formulation.pre();
+ formulation.preDistributed();
+ formulation.assemble();
+ formulation.solveDistributed(system::DISTRIBUTED_SOLVE_TYPE::ORAS, false);
+
+ auto error1 = real(integrate(pow(real(v - analytical), 2), omega, "Gauss6")) / real(integrate(pow(real(analytical), 2), omega, "Gauss6"));
+ msg::print << "L2 absolute error in domain " << MPI_Rank << " : " << error1 << msg::endl;
+ save(v, omega, "v_" + std::to_string(MPI_Rank));
+ save(analytical, omega, "sol_" + std::to_string(MPI_Rank));
+
+ if (error1 > 1e-4)
+ {
+ throw common::Exception("Error above tolerance : " + std::to_string(error1));
+ }
+
+ return 0;
+}