diff --git a/src/problem/Formulation.cpp b/src/problem/Formulation.cpp
index 0f326a37fac8a6ed0909016ba775240d1c253879..59bbff9a78cb56baef74f6769867c8aeb728b154 100644
--- a/src/problem/Formulation.cpp
+++ b/src/problem/Formulation.cpp
@@ -344,7 +344,7 @@ namespace gmshddm
     }
 
     template< class T_Scalar >
-    gmshfem::common::Timer Formulation< T_Scalar >::pre()
+    gmshfem::common::Timer Formulation< T_Scalar >::pre(bool mustAssemble)
     {
       const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
       const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
@@ -798,7 +798,51 @@ namespace gmshddm
 
       _infoNumberOfDofs();
 
-      s_printResources("DDM pre-processing before volume assembly");
+
+      if(mustAssemble) {
+        s_printResources("DDM pre-processing before volume assembly");
+
+        assemble();
+      }
+
+      time.tock();
+      if(MPI_Rank == 0) {
+        gmshfem::msg::info << "Done pre-processing in " << time << "s" << gmshfem::msg::endl;
+      }
+
+      return time;
+    }
+
+
+    template< class T_Scalar >
+    gmshfem::common::Timer Formulation< T_Scalar >::assemble()
+    {
+      gmshfem::common::Timer time;
+      time.tick();
+
+      const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
+      const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
+
+      _physicalCommutator = true;
+      _artificialCommutator = false;
+      togglePhysicalAndArtificialSourceTerms();
+
+      // Clear volume
+      for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
+        if(mpi::isItMySubdomain(idom)) {
+          // Clear formulation
+          _volume[idom]->setSystemToZero();
+        }
+      }
+
+      for(auto idom = 0ULL; idom < _surface.size(); ++idom) {
+        if(mpi::isItMySubdomain(idom)) {
+          for(auto it = _surface[idom].begin(); it != _surface[idom].end(); ++it) {
+            gmshfem::msg::info << " - Clearing " << it->second->name() << gmshfem::msg::endl;
+            it->second->setSystemToZero();
+          }
+        }
+      }
 
       // Assemble locally on each processor
       _assembleAllVolume();
@@ -820,12 +864,13 @@ namespace gmshddm
       _extractRHS();
 
       MPI_Barrier(PETSC_COMM_WORLD);
-
+      
       s_printResources("DDM pre-processing end");
 
+      
       time.tock();
       if(MPI_Rank == 0) {
-        gmshfem::msg::info << "Done pre-processing in " << time << "s" << gmshfem::msg::endl;
+        gmshfem::msg::info << "Done assembling in " << time << "s" << gmshfem::msg::endl;
       }
 
       return time;
diff --git a/src/problem/Formulation.h b/src/problem/Formulation.h
index 3ce7aaf4b329ce1168333360a149278e4d59b1cd..ec02567d8ffdc6d29fc62e4cba5f01404ba20f27 100644
--- a/src/problem/Formulation.h
+++ b/src/problem/Formulation.h
@@ -90,7 +90,8 @@ namespace gmshddm
       gmshfem::problem::Formulation< T_Scalar > &operator()(unsigned int i);
       gmshfem::problem::Formulation< T_Scalar > &operator()(unsigned int i, unsigned int j);
 
-      gmshfem::common::Timer pre();
+      gmshfem::common::Timer pre(bool mustAssemble = true);
+      gmshfem::common::Timer assemble();
 
       gmshfem::common::Timer solve(const std::string &solver, const double tolerance = 1e-6, const int iterMax = 1000, const bool sameMatrixWithArtificialAndPhysicalSources = false, const bool skipFinalSolutionComputation = false);
       gmshfem::common::Timer solve(AbstractIterativeSolver<T_Scalar> &solver, const double tolerance = 1e-6, const int iterMax = 1000, const bool sameMatrixWithArtificialAndPhysicalSources = false, const bool skipFinalSolutionComputation = false);