Merge branch 'overlap_jacobi' into 'master'

New tutorial case: overlapping DDM (Dirichlet) for Laplace equation

See merge request !17

tutorials/poisson_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_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_overlap/domain.geo

+// Gmsh project created on Thu Oct  6 12:29:07 2022
+//SetFactory("OpenCASCADE");
+//+
+Point(1) = {0, 0, 0, 0.05};
+//+
+Point(2) = {0.8, 0, 0, 0.05};
+//+
+Point(3) = {1.2, 0, 0, 0.05};
+//+
+Point(4) = {2.0, 0, 0, 0.05};
+//+
+Point(5) = {2, 1, 0, 0.05};
+//+
+Point(6) = {1.2, 1, 0, 0.05};
+//+
+Point(7) = {0.8, 1, 0, 0.05};
+//+
+Point(8) = {0.0, 1, 0, 0.05};
+//+
+Line(1) = {1, 2};
+//+
+Line(2) = {2, 7};
+//+
+Line(3) = {7, 8};
+//+
+Line(4) = {8, 1};
+//+
+Line(5) = {3, 2};
+//+
+Line(6) = {6, 3};
+//+
+Line(7) = {6, 7};
+//+
+Line(8) = {3, 4};
+//+
+Line(9) = {4, 5};
+//+
+Line(10) = {5, 6};
+//+
+Curve Loop(1) = {4, 1, 2, 3};
+//+
+Plane Surface(1) = {1};
+//+
+Curve Loop(2) = {7, -2, -5, -6};
+//+
+Plane Surface(2) = {2};
+//+
+Curve Loop(3) = {10, 6, 8, 9};
+//+
+Plane Surface(3) = {3};
+//+
+Physical Surface("left", 11) = {1};
+//+
+Physical Surface("center", 12) = {2};
+//+
+Physical Surface("right", 13) = {3};
+//+
+Physical Curve("gammaLeft", 14) = {4};
+//+
+Physical Curve("gammaRight", 15) = {9};
+//+
+Physical Curve("interface1", 16) = {2};
+//+
+Physical Curve("interface2", 17) = {6};

tutorials/poisson_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 left("left"); // 0 to 1.2
+  Domain center("center"); // 0.8 to 1.2
+  Domain right("right"); // 0.8 to 2
+
+  Domain gammaLeft("gammaLeft"); // 0
+  Domain interface1("interface1"); // 0.8
+  Domain interface2("interface2"); // 1.2
+  Domain gammaRight("gammaRight"); // 2.0
+
+  Subdomain omega(nDom);
+  Interface sigma(nDom);
+  omega(0) = gammaLeft | left | interface1 | center;
+  omega(1) = gammaRight | center | interface2 | right;
+
+  sigma(0, 1) = interface1 | interface2;
+  sigma(1, 0) = interface2 | interface1;
+
+  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(gammaLeft, 0);
+  u(1).addConstraint(gammaRight, 2);
+  auto trueSolution = x<std::complex< double >>();
+
+
+  InterfaceField< std::complex< double >, Form::Form0 > g("g", sigma, FunctionSpaceTypeForm0::HierarchicalH1, order);
+  formulation.addInterfaceField(g);
+
+  // First domain
+  // Laplacian(u) = 0 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), gauss);
+  formulation(0).integral((-g(1, 0)), tf(u(0)), interface2, gauss);
+  formulation(0).integral(dof(u(0)), tf(u(0)), interface2, 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(-(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);
+
+  // Second domain
+  formulation(1).integral(grad(dof(u(1))), grad(tf(u(1))), omega(1), gauss);
+  formulation(1).integral((-g(0, 1)), tf(u(1)), interface1, gauss);
+  formulation(1).integral(dof(u(1)), tf(u(1)), interface1, gauss);
+
+  // Second interface
+  formulation(1, 0).integral(-dof(g(1, 0)), tf(g(1, 0)), sigma(1, 0), gauss);
+  formulation(1, 0).integral(-(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.pre();
+  int iterMax = 1000;
+  gmshDdm.userDefinedParameter(iterMax, "iter");
+  formulation.solve("jacobi", 1e-10, 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));
+    }
+  }
+
+
+  return 0;
+}