diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index afe758b7ebb500884aa67674dcffebbbb75dd260..9d066191a12a2b8d1d095835d88f3424068fad0d 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -59,6 +59,19 @@ void Process::initialize()
     DBUG("Initialize boundary conditions.");
     _boundary_conditions.addBCsForProcessVariables(
         _process_variables, *_local_to_global_index_map, _integration_order);
+
+    for (int variable_id = 0;
+         variable_id < static_cast<int>(_process_variables.size());
+         ++variable_id)
+    {
+        ProcessVariable& pv = _process_variables[variable_id];
+        auto sts =
+            pv.createSourceTerms(*_local_to_global_index_map, variable_id,
+                                 _integration_order);
+
+        std::move(sts.begin(), sts.end(),
+                  std::back_inserter(_source_terms));
+    }
 }
 
 void Process::setInitialConditions(double const t, GlobalVector& x)
@@ -128,6 +141,11 @@ void Process::assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
     assembleConcreteProcess(t, x, M, K, b);
 
     _boundary_conditions.applyNaturalBC(t, x, K, b);
+
+    for (auto const& st : _source_terms)
+    {
+        st->integrateNodalSourceTerm(t, b);
+    }
 }
 
 void Process::assembleWithJacobian(const double t, GlobalVector const& x,
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 88ef33bf5ee60e203fba6c30ad7822211c3ec475..b9fba6a81f6c91bb7e5d1fcbafa1d0274f0fb353 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -14,6 +14,7 @@
 #include "NumLib/ODESolver/TimeDiscretization.h"
 #include "NumLib/NamedFunctionCaller.h"
 #include "ProcessLib/BoundaryCondition/BoundaryConditionCollection.h"
+#include "ProcessLib/SourceTerms/NodalSourceTerm.h"
 #include "ProcessLib/Parameter/Parameter.h"
 
 #include "ExtrapolatorData.h"
@@ -218,6 +219,7 @@ private:
     std::vector<std::reference_wrapper<ProcessVariable>> _process_variables;
 
     BoundaryConditionCollection _boundary_conditions;
+    std::vector<std::unique_ptr<NodalSourceTerm>> _source_terms;
 
     ExtrapolatorData _extrapolator_data;
 };