working!

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(1) = {0, 0, 0, 0.02};
 //+
-Point(2) = {1, 0, 0, 0.05};
+Point(2) = {1, 0, 0, 0.02};
 //+
-Point(3) = {2, 0, 0, 0.05};
+Point(3) = {2, 0, 0, 0.02};
 //+
-Point(4) = {2, 1, 0, 0.05};
+Point(4) = {2, 1, 0, 0.02};
 //+
-Point(5) = {1, 1, 0, 0.05};
+Point(5) = {1, 1, 0, 0.02};
 //+
-Point(6) = {0, 1, 0, 0.05};
+Point(6) = {0, 1, 0, 0.02};
 //+
 Line(1) = {1, 2};
 //+

tutorials/poisson_robin_no_overlap/main.cpp

@@ -52,8 +52,8 @@ int main(int argc, char **argv)
 
   Subdomain omega(nDom);
   Interface sigma(nDom);
-  omega(0) = omega0;// | interface;
-  omega(1) = omega1;// | interface;
+  omega(0) = omega0;
+  omega(1) = omega1;
 
   sigma(0, 1) = interface;
   sigma(1, 0) = interface;
@@ -76,35 +76,26 @@ int main(int argc, char **argv)
   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);
+  double alpha = 1;
 
   // 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);
+  for (int i = 0; i <= 1; ++i) {
+    std::complex<double> im(0.0, 1.0);
+    auto other = 1-i;
+    // Domain
+    formulation(i).integral(grad(dof(u(i))), grad(tf(u(i))), omega(i), gauss);
+    formulation(i).integral(formulation.physicalSource(-1.0), tf(u(i)), omega(i), gauss);
+    formulation(i).integral(formulation.artificialSource(-g(other, i)), tf(u(i)), interface, gauss);
+    formulation(i).integral(alpha * (dof(u(i))), tf(u(i)), interface, gauss);
+
+    // Interface
+    formulation(i, other).integral(dof(g(i, other)), tf(g(i, other)), sigma(i, other), gauss);
+    formulation(i, other).integral(formulation.artificialSource(g(other, i)), tf(g(i, other)), sigma(i, other), gauss);
+    // It is critical that u(i) be not an artificial source (why ?)
+    formulation(i, other).integral(((-2.0 * alpha) * u(i)), tf(g(i, other)), sigma(i, other), gauss);
+  }
 
 
   formulation.pre();
@@ -117,7 +108,6 @@ int main(int argc, char **argv)
     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");
     }
   }