first draft

tutorials/poisson_robin_no_overlap/CMakeLists.txt

+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})

tutorials/poisson_robin_no_overlap/README.md

+# Poisson Equation with overlap and Dirichlet 
+
+We solve a homogeneous Poisson (i.e. Laplace) equation with domain decomposition.
+The domain is a 2 by one rectangle, with boundary conditions $`u(0, y) = 0, u(2, y) = 2`$ and homogeneous von Neumann conditions on the other edges. The analytical solution is $`u(x,y) = x`$.
+
+This tutorial illustrates overlapping DDM. The first domain solves the equation for $`x`$ ranging from 0 to 1.2, and the other domain from 0.8 to 2. They exchange data over the center of the domain (0.8 to 1.2), which provide the boundary conditions on the inner boundaries.
\ No newline at end of file

tutorials/poisson_robin_no_overlap/domain.geo

+// Gmsh project created on Fri Oct  7 15:45:08 2022
+//+
+Point(1) = {0, 0, 0, 0.05};
+//+
+Point(2) = {1, 0, 0, 0.05};
+//+
+Point(3) = {2, 0, 0, 0.05};
+//+
+Point(4) = {2, 1, 0, 0.05};
+//+
+Point(5) = {1, 1, 0, 0.05};
+//+
+Point(6) = {0, 1, 0, 0.05};
+//+
+Line(1) = {1, 2};
+//+
+Line(2) = {2, 3};
+//+
+Line(3) = {4, 3};
+//+
+Line(4) = {4, 5};
+//+
+Line(5) = {5, 6};
+//+
+Line(6) = {6, 1};
+//+
+Line(7) = {2, 5};
+//+
+Physical Curve("left", 8) = {6};
+//+
+Physical Curve("interface", 9) = {7};
+//+
+Physical Curve("right", 10) = {3};
+//+
+Curve Loop(1) = {6, 1, 7, 5};
+//+
+Plane Surface(1) = {1};
+//+
+Curve Loop(2) = {4, -7, 2, -3};
+//+
+Plane Surface(2) = {2};
+//+
+Physical Surface("omega0", 11) = {1};
+//+
+Physical Surface("omega1", 12) = {2};

tutorials/poisson_robin_no_overlap/main.cpp

+#include "gmsh.h"
+
+#include <gmshddm/Formulation.h>
+#include <gmshddm/GmshDdm.h>
+#include <gmshddm/MPIInterface.h>
+#include <gmshfem/Message.h>
+
+
+using namespace gmshddm;
+using namespace gmshddm::common;
+using namespace gmshddm::domain;
+using namespace gmshddm::problem;
+using namespace gmshddm::field;
+using namespace gmshddm::mpi;
+
+using namespace gmshfem;
+using namespace gmshfem::domain;
+using namespace gmshfem::equation;
+using namespace gmshfem::problem;
+using namespace gmshfem::field;
+using namespace gmshfem::function;
+using namespace gmshfem::post;
+
+
+int main(int argc, char **argv)
+{
+  GmshDdm gmshDdm(argc, argv);
+
+  unsigned int nDom = 2;
+  int order = 1;
+  double lc = 0.05;
+  std::string gauss = "Gauss10";
+
+  gmshDdm.userDefinedParameter(nDom, "nDom");
+  gmshDdm.userDefinedParameter(order, "order");
+  gmshDdm.userDefinedParameter(lc, "lc");
+  gmshDdm.userDefinedParameter(gauss, "gauss");
+
+  gmsh::open("../domain.geo"); // Assume we run from the "build" subfolder
+  gmsh::model::mesh::generate();
+
+  // Domains are partitions of the x-axis.
+  Domain omega0("omega0");
+  Domain omega1("omega1");
+
+  // Curves
+  Domain left("left"); // 0 to 1.2
+  Domain interface("interface"); // 0.8 to 1.2
+  Domain right("right"); // 0.8 to 2
+
+  
+
+  Subdomain omega(nDom);
+  Interface sigma(nDom);
+  omega(0) = omega0;// | interface;
+  omega(1) = omega1;// | interface;
+
+  sigma(0, 1) = interface;
+  sigma(1, 0) = interface;
+
+  std::vector< std::vector< unsigned int > > topology(nDom);
+  topology[0] = {1};
+  topology[1] = {0};
+
+  // Create DDM formulation
+  // And BCs on the subdomain problems
+  gmshddm::problem::Formulation< std::complex< double > > formulation("Laplace", topology);
+  SubdomainField< std::complex< double >, Form::Form0 > u("u", omega,  FunctionSpaceTypeForm0::HierarchicalH1, order);
+
+  // Toy problem with solution u(x,y) = x
+  u(0).addConstraint(left, 0);
+  u(1).addConstraint(right, 0);
+  auto trueSolution = 0.5 * x<std::complex< double >>() * (2.0 - x<std::complex< double >>());
+
+
+  InterfaceField< std::complex< double >, Form::Form0 > g("g", sigma, FunctionSpaceTypeForm0::HierarchicalH1, order);
+  formulation.addInterfaceField(g);
+
+
+  // First domain
+  // Laplacian(u) = -1 on the domain, and u matches g(1,0) on the interface
+  formulation(0).integral(grad(dof(u(0))), grad(tf(u(0))), omega(0)|interface, gauss);
+  formulation(0).integral(-1.0, tf(u(0)), omega(0)|interface, gauss);
+
+  formulation(0).integral((formulation.artificialSource(-g(1, 0))), tf(u(0)), interface, gauss);
+  //formulation(0).integral(dof(u(0)), tf(u(0)), interface, gauss);
+
+  // First interface (What 0 sends to 1)
+  // Compute g(0,1) such that g(0,1) + g(1,0) = 2 u(0)
+  formulation(0, 1).integral(-dof(g(0, 1)), tf(g(0, 1)), sigma(0, 1), gauss);
+  formulation(0, 1).integral(-formulation.artificialSource(g(1, 0)), tf(g(0, 1)), sigma(0, 1), gauss);
+  formulation(0, 1).integral((2 * u(0)), tf(g(0, 1)), sigma(0, 1), gauss);
+  //formulation(0, 1).integral((2* grad(u(0))) * normal<std::complex<double>>(), tf(g(0, 1)), sigma(0, 1), gauss);
+
+
+  // Second domain
+  formulation(1).integral(grad(dof(u(1))), grad(tf(u(1))), omega(1)|interface, gauss);
+  formulation(1).integral(-1.0, tf(u(1)), omega(1)|interface, gauss);
+
+  formulation(1).integral(formulation.artificialSource(-g(0, 1)), tf(u(1)), interface, gauss);
+  //formulation(1).integral(dof(u(1)), tf(u(1)), interface, gauss);
+
+  // Second interface
+  formulation(1, 0).integral(-dof(g(1, 0)), tf(g(1, 0)), sigma(1, 0), gauss);
+  formulation(1, 0).integral(-formulation.artificialSource(g(0, 1)), tf(g(1, 0)), sigma(1, 0), gauss);
+  formulation(1, 0).integral((2*u(1)), tf(g(1, 0)), sigma(1, 0), gauss);
+  //formulation(1, 0).integral((2* grad(u(1))) * normal<std::complex<double>>(), tf(g(0, 1)), sigma(0, 1), gauss);
+
+
+  formulation.pre();
+  int iterMax = 1000;
+  gmshDdm.userDefinedParameter(iterMax, "iter");
+  formulation.solve("gmres", 1e-6, iterMax);
+
+
+  for(unsigned int i = 0; i < nDom; ++i) {
+    if(isItMySubdomain(i)) {
+      save(u(i), omega(i), "u_" + std::to_string(i));
+      save(u(i) - trueSolution, omega(i), "err_" + std::to_string(i));
+      save(g(i, 1-i), sigma(i,1-i), "g");
+    }
+  }
+
+  save(trueSolution, omega(0) | omega(1) | interface, "true_sol");
+
+
+  return 0;
+}