diff --git a/examples/maxwell/main.cpp b/examples/maxwell/main.cpp
index 0b8ad67676bba5fe1c8432fe7d02be96932d9b0c..89b65a8750bc59e54154d4b0e7cfd429e72bd531 100644
--- a/examples/maxwell/main.cpp
+++ b/examples/maxwell/main.cpp
@@ -140,6 +140,10 @@ int main(int argc, char **argv)
     im * beta * kz / (kc * kc) *
       cos< std::complex< double > >(kz * z< std::complex< double > >()) *
       sin< std::complex< double > >(ky * y< std::complex< double > >()));
+
+  VectorFunction< std::complex< double > > functionZero =
+    vector< std::complex< double > >(0., 0., 0.);
+
   /****
      * FEM
      ****/
@@ -209,8 +213,8 @@ int main(int argc, char **argv)
     // VectorFunction<std::complex<double>> n = normal<std::complex<double>>();
     // VOLUME TERMS
 
-    E.addConstraint(border, vector< std::complex< double > >(0., 0., 0.));
-    lambda.addConstraint(border, vector< std::complex< double > >(0., 0., 0.));
+    E.addConstraint(border, functionZero);
+    lambda.addConstraint(border, functionZero);
 
     formulation.integral(curl(dof(E)), curl(tf(E)), omega, gauss);
     formulation.integral(-k * k * dof(E), tf(E), omega, gauss);
@@ -294,8 +298,7 @@ int main(int argc, char **argv)
     // subdomains
     formulation.addInterfaceField(g);
     VectorFunction< std::complex< double > > n = normal< std::complex< double > >();
-    VectorFunction< std::complex< double > > functionZero =
-      vector< std::complex< double > >(0., 0., 0.);
+
     for(unsigned int i = 0; i < unsigned(nDom); ++i) {
       E(i).addConstraint(border(i), functionZero);
       lambda(i).addConstraint(border(i), functionZero);
@@ -346,9 +349,14 @@ int main(int argc, char **argv)
                                    sigma(i, j), gauss);*/
       }
     }
+
     formulation.pre();
-    // Solve the DDM formulation
-    formulation.solve("gmres", 1e-4, 100, false); // FIXME: last argument should be autodetected!
+    // Solve the DDM formulation; the last argument tells the solver that the
+    // subdomain matrices will be identical with physical and artificial sources
+    // (no bilinear terms involving the sources and no non-homogeneous Dirichlet
+    // BCs) - this saves 2 factorizations per subdomain
+    formulation.solve("gmres", 1e-4, 100, true);
+
     //save field E
     for(int i = 0; i < nDom; ++i) {
       save(E(i), omega(i), "E" + std::to_string(i));
diff --git a/src/problem/Formulation.cpp b/src/problem/Formulation.cpp
index 5b2a561846a2a0f0805f2b82b7d330af14850ce0..c3b5920540def01c9ebe63c3998ae8b95d56a72f 100644
--- a/src/problem/Formulation.cpp
+++ b/src/problem/Formulation.cpp
@@ -664,9 +664,10 @@ namespace gmshddm
     template< class T_Scalar >
     gmshfem::common::Timer Formulation< T_Scalar >::solve(const std::string &solver, const double tolerance, const int iterMax, const bool sameMatrixWithArtificialAndPhysicalSources)
     {
-      // TODO: sameMatrixWithArtificialAndPhysicalSources could be automatically
-      // detected by examining if bilinear terms of volume formulations use the
-      // physical/artifical communicator
+      // sameMatrixWithArtificialAndPhysicalSources could be automatically set
+      // to true if bilinear terms of volume formulations do not use the
+      // physical/artifical communicator and if there are no non-homogeneous
+      // Dirichlet BCs
 
       gmshfem::common::Timer time;
       time.tick();