diff --git a/tutorials/poisson_robin_no_overlap/CMakeLists.txt b/tutorials/poisson_robin_no_overlap/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3ed63f13623e9d7ed45010b2a7ffd418461a2b0c
--- /dev/null
+++ b/tutorials/poisson_robin_no_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_robin_no_overlap/README.md b/tutorials/poisson_robin_no_overlap/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..cc49cc5f69f2ef2c94642fb22c58c033bb24a62c
--- /dev/null
+++ b/tutorials/poisson_robin_no_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_robin_no_overlap/domain.geo b/tutorials/poisson_robin_no_overlap/domain.geo
new file mode 100644
index 0000000000000000000000000000000000000000..805d938ef162b3b014b9b1b633cd9fd9f44ac643
--- /dev/null
+++ b/tutorials/poisson_robin_no_overlap/domain.geo
@@ -0,0 +1,45 @@
+// 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};
diff --git a/tutorials/poisson_robin_no_overlap/main.cpp b/tutorials/poisson_robin_no_overlap/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0e3690747079d48695b0650b1897694db3d5a180
--- /dev/null
+++ b/tutorials/poisson_robin_no_overlap/main.cpp
@@ -0,0 +1,128 @@
+#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;
+}