diff --git a/NumLib/ODESolver/ODESystem.h b/NumLib/ODESolver/ODESystem.h
index b362403f9b547db0410269267dec5cbb2d550565..7dd004d058d15f88f1d8bf4a4a521201a8eea9c9 100644
--- a/NumLib/ODESolver/ODESystem.h
+++ b/NumLib/ODESolver/ODESystem.h
@@ -44,6 +44,9 @@ public:
     static const ODESystemTag ODETag =
         ODESystemTag::FirstOrderImplicitQuasilinear;
 
+    //! Calls process' pre-assembly with the provided state (\c t, \c x).
+    virtual void preAssemble(const double t, GlobalVector const& x) = 0;
+
     //! Assemble \c M, \c K and \c b at the provided state (\c t, \c x).
     virtual void assemble(
         const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
@@ -73,6 +76,10 @@ class ODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
                        NonlinearSolverTag::Picard>
 {
 public:
+
+    //! Calls process' pre-assembly with the provided state (\c t, \c x).
+    virtual void preAssemble(const double t, GlobalVector const& x) = 0;
+
     /*! Assemble \c M, \c K, \c b and the Jacobian
      * \f$ \mathtt{Jac} := \partial r/\partial x_N \f$
      * at the provided state (\c t, \c x).
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index a1b6f29685cc0820daf6011fb53fd7615a9f08ae..064c3d40c5625a97fe8125066cc2869074d22bae 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -84,6 +84,7 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
     _b->setZero();
     _Jac->setZero();
 
+    _ode.preAssemble(t, x_curr);
     _ode.assembleWithJacobian(t, x_curr, xdot, dxdot_dx, dx_dx, *_M, *_K, *_b,
                               *_Jac);
 
@@ -186,6 +187,7 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
     _K->setZero();
     _b->setZero();
 
+    _ode.preAssemble(t, x_curr);
     _ode.assemble(t, x_curr, *_M, *_K, *_b);
 
     LinAlg::finalizeAssembly(*_M);
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index bfb983f87d408f22fe43dbf92bc63c6ac1e02be3..e3fc5cb7983663113db12b0db010403dcd158563 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -35,6 +35,9 @@ class LocalAssemblerInterface
 public:
     virtual ~LocalAssemblerInterface() = default;
 
+    virtual void preAssemble(double const /*t*/,
+                             std::vector<double> const& /*local_x*/){};
+
     virtual void assemble(double const t, std::vector<double> const& local_x,
                           std::vector<double>& local_M_data,
                           std::vector<double>& local_K_data,
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 808bb909c59a4d959589d5df4d18e06e1c9e2859..119e9063dd9e36bdfd5d030ab33c6b01751870be 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -132,6 +132,11 @@ MathLib::MatrixSpecifications Process::getMatrixSpecifications() const
             &l.getGhostIndices(), &_sparsity_pattern};
 }
 
+void Process::preAssemble(const double t, GlobalVector const& x)
+{
+    preAssembleConcreteProcess(t, x);
+}
+
 void Process::assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
                        GlobalMatrix& K, GlobalVector& b)
 {
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index d067d0ff747ba5b10b5807aee426df99c332e0e5..ef2ef4d1b096652af79a874469fe49e9c07e1513 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -77,6 +77,8 @@ public:
         _coupled_solutions = coupled_solutions;
     }
 
+    void preAssemble(const double t, GlobalVector const& x) override final;
+
     virtual void setCoupledSolutionsForStaggeredSchemeToLocalAssemblers() {}
     void assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
                   GlobalMatrix& K, GlobalVector& b) final;
@@ -142,6 +144,11 @@ private:
         MeshLib::Mesh const& mesh,
         unsigned const integration_order) = 0;
 
+    virtual void preAssembleConcreteProcess(const double /*t*/,
+                                            GlobalVector const& /*x*/)
+    {
+    }
+
     virtual void assembleConcreteProcess(const double t, GlobalVector const& x,
                                          GlobalMatrix& M, GlobalMatrix& K,
                                          GlobalVector& b) = 0;
diff --git a/ProcessLib/VectorMatrixAssembler.cpp b/ProcessLib/VectorMatrixAssembler.cpp
index 398c7a5d45bc84d7ae2b429ee0ccc66a579cfe2f..690c9bcc23ce11b2615169a55e970cbd5cf65bbf 100644
--- a/ProcessLib/VectorMatrixAssembler.cpp
+++ b/ProcessLib/VectorMatrixAssembler.cpp
@@ -55,6 +55,17 @@ VectorMatrixAssembler::VectorMatrixAssembler(
 {
 }
 
+void VectorMatrixAssembler::preAssemble(
+    const std::size_t mesh_item_id, LocalAssemblerInterface& local_assembler,
+    const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
+    const GlobalVector& x)
+{
+    auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
+    auto const local_x = x.get(indices);
+
+    local_assembler.preAssemble(t, local_x);
+}
+
 void VectorMatrixAssembler::assemble(
     const std::size_t mesh_item_id, LocalAssemblerInterface& local_assembler,
     const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
diff --git a/ProcessLib/VectorMatrixAssembler.h b/ProcessLib/VectorMatrixAssembler.h
index 5f305cc482b328e08c8f9231efbcb1dcb2c8e04a..a1203361b4d6084aec87bfd4fca0aa54b22aaa6d 100644
--- a/ProcessLib/VectorMatrixAssembler.h
+++ b/ProcessLib/VectorMatrixAssembler.h
@@ -33,6 +33,11 @@ public:
     explicit VectorMatrixAssembler(
         std::unique_ptr<AbstractJacobianAssembler>&& jacobian_assembler);
 
+    void preAssemble(const std::size_t mesh_item_id,
+                     LocalAssemblerInterface& local_assembler,
+                     const NumLib::LocalToGlobalIndexMap& dof_table,
+                     const double t, const GlobalVector& x);
+
     //! Assembles\c M, \c K, and \c b.
     //! \remark Jacobian is not assembled here, see assembleWithJacobian().
     void assemble(std::size_t const mesh_item_id,
diff --git a/Tests/NumLib/ODEs.h b/Tests/NumLib/ODEs.h
index f5326f712fff0c120106d90fa2ed8cd2a455683c..fbe3f27156363517a9fefb610f3f415184be8300 100644
--- a/Tests/NumLib/ODEs.h
+++ b/Tests/NumLib/ODEs.h
@@ -26,6 +26,8 @@ class ODE1 final : public NumLib::ODESystem<
                        NumLib::NonlinearSolverTag::Newton>
 {
 public:
+    void preAssemble(const double t, GlobalVector const& x) override {}
+
     void assemble(const double /*t*/, GlobalVector const& /*x*/,
                   GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) override
     {
@@ -102,6 +104,8 @@ class ODE2 final : public NumLib::ODESystem<
                        NumLib::NonlinearSolverTag::Newton>
 {
 public:
+    void preAssemble(const double t, GlobalVector const& x) override {}
+
     void assemble(const double /*t*/, GlobalVector const& x, GlobalMatrix& M,
                   GlobalMatrix& K, GlobalVector& b) override
     {
@@ -187,6 +191,8 @@ class ODE3 final : public NumLib::ODESystem<
                        NumLib::NonlinearSolverTag::Newton>
 {
 public:
+    void preAssemble(const double t, GlobalVector const& x) override {}
+
     void assemble(const double /*t*/, GlobalVector const& x_curr, GlobalMatrix& M,
                   GlobalMatrix& K, GlobalVector& b) override
     {