diff --git a/tutorials/poisson_overlap/CMakeLists.txt b/tutorials/poisson_overlap/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3ed63f13623e9d7ed45010b2a7ffd418461a2b0c
--- /dev/null
+++ b/tutorials/poisson_overlap/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/tutorials/poisson_overlap/README.md b/tutorials/poisson_overlap/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..cc49cc5f69f2ef2c94642fb22c58c033bb24a62c
--- /dev/null
+++ b/tutorials/poisson_overlap/README.md
@@ -0,0 +1,6 @@
+# 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
diff --git a/tutorials/poisson_overlap/domain.geo b/tutorials/poisson_overlap/domain.geo
new file mode 100644
index 0000000000000000000000000000000000000000..283f0836ac1f517266a36dcff7cd5df9525c1fca
--- /dev/null
+++ b/tutorials/poisson_overlap/domain.geo
@@ -0,0 +1,64 @@
+// 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};
diff --git a/tutorials/poisson_overlap/main.cpp b/tutorials/poisson_overlap/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..405d84fe9c4bc7b52b18938bde2b10b9b02c948c
--- /dev/null
+++ b/tutorials/poisson_overlap/main.cpp
@@ -0,0 +1,116 @@
+#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;
+}