diff --git a/ProcessLib/HeatConduction/Tests.cmake b/ProcessLib/HeatConduction/Tests.cmake
index 42874e37d524b10ee94cc6556f6faffb51555c13..636b6ded332688087a6a662f8889d862f98524c3 100644
--- a/ProcessLib/HeatConduction/Tests.cmake
+++ b/ProcessLib/HeatConduction/Tests.cmake
@@ -55,3 +55,28 @@ AddTest(
     standard_solution_bhe2d_pcs_0_ts_840_t_72576000.000000.vtu bhe2d_pcs_0_ts_840_t_72576000.000000.vtu temperature temperature 1e-12 0.0
     REQUIREMENTS NOT OGS_USE_MPI
 )
+
+# test the source term on a subdomain
+AddTest(
+    NAME 1D_HeatConduction_dirichlet_SourceTerm
+    PATH Parabolic/T/1D_dirichlet_source-term
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS line_1_line_1e2_source_term.prj
+    TESTER vtkdiff
+    DIFF_DATA
+    line_1_line_1e2_pcs_0_ts_500_t_39062500.000000_reference.vtu line_1_line_1e2_pcs_0_ts_500_t_39062500.000000.vtu temperature temperature 1e-11 0.0
+    REQUIREMENTS NOT OGS_USE_MPI
+)
+
+# test the source term on a subdomain with the PETSc embedded executable file
+AddTest(
+    NAME 1D_HeatConduction_dirichlet_SourceTerm_PETSc
+    PATH Parabolic/T/1D_dirichlet_source-term
+    EXECUTABLE_ARGS line_1_line_1e2_source_term.prj
+    WRAPPER mpirun
+    WRAPPER_ARGS -np 1
+    TESTER vtkdiff
+    REQUIREMENTS OGS_USE_MPI
+    DIFF_DATA
+    line_1_line_1e2_pcs_0_ts_500_t_39062500.000000_reference.vtu line_1_line_1e2_pcs_0_ts_500_t_39062500_000000_0.vtu temperature temperature 1e-10 0.0
+)
diff --git a/ProcessLib/ProcessVariable.cpp b/ProcessLib/ProcessVariable.cpp
index 475bf1e570fee7a26e26c2689a4228ddbaef0c0f..28b79663ed25facabf688164d75ee0bff11e5808 100644
--- a/ProcessLib/ProcessVariable.cpp
+++ b/ProcessLib/ProcessVariable.cpp
@@ -222,7 +222,7 @@ ProcessVariable::createSourceTerms(
 
     for (auto& config : _source_term_configs)
         source_terms.emplace_back(createSourceTerm(
-            config, dof_table, _mesh, variable_id, integration_order,
+            config, dof_table, config.mesh, variable_id, integration_order,
             _shapefunction_order, parameters));
 
     return source_terms;
diff --git a/ProcessLib/SourceTerms/CreateNodalSourceTerm.cpp b/ProcessLib/SourceTerms/CreateNodalSourceTerm.cpp
index 04c2227954dc8594f1ce424b00c8dc5c4d088530..1d0bb77294d23ca07765b043d010d56d599ff325 100644
--- a/ProcessLib/SourceTerms/CreateNodalSourceTerm.cpp
+++ b/ProcessLib/SourceTerms/CreateNodalSourceTerm.cpp
@@ -20,8 +20,8 @@ namespace ProcessLib
 {
 std::unique_ptr<SourceTerm> createNodalSourceTerm(
     BaseLib::ConfigTree const& config, MeshLib::Mesh const& st_mesh,
-    const NumLib::LocalToGlobalIndexMap& dof_table,
-    std::size_t const bulk_mesh_id, const int variable_id,
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap> dof_table,
+    std::size_t const source_term_mesh_id, const int variable_id,
     const int component_id,
     std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters)
 {
@@ -35,7 +35,8 @@ std::unique_ptr<SourceTerm> createNodalSourceTerm(
 
     auto& param = findParameter<double>(param_name, parameters, 1);
 
-    return std::make_unique<NodalSourceTerm>(dof_table, bulk_mesh_id, st_mesh,
+    return std::make_unique<NodalSourceTerm>(std::move(dof_table),
+                                             source_term_mesh_id, st_mesh,
                                              variable_id, component_id, param);
 }
 
diff --git a/ProcessLib/SourceTerms/CreateNodalSourceTerm.h b/ProcessLib/SourceTerms/CreateNodalSourceTerm.h
index 9f0c97169452ce473e89ff5c00f7e183afc10b63..8922dc94e2d0e0ed36d7e5682629c60db0fe1ace 100644
--- a/ProcessLib/SourceTerms/CreateNodalSourceTerm.h
+++ b/ProcessLib/SourceTerms/CreateNodalSourceTerm.h
@@ -34,8 +34,8 @@ namespace ProcessLib
 {
 std::unique_ptr<SourceTerm> createNodalSourceTerm(
     BaseLib::ConfigTree const& config, MeshLib::Mesh const& st_mesh,
-    const NumLib::LocalToGlobalIndexMap& dof_table, std::size_t mesh_id,
-    const int variable_id, const int component_id,
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap> dof_table,
+    std::size_t mesh_id, const int variable_id, const int component_id,
     std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters);
 
 }   // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/CreateSourceTerm.cpp b/ProcessLib/SourceTerms/CreateSourceTerm.cpp
index b33cd640db6d86103e7bd2772411260ac2aca2d4..b6f4c2709546987fe388a24af0696699f5d074af 100644
--- a/ProcessLib/SourceTerms/CreateSourceTerm.cpp
+++ b/ProcessLib/SourceTerms/CreateSourceTerm.cpp
@@ -21,36 +21,83 @@ namespace ProcessLib
 {
 std::unique_ptr<SourceTerm> createSourceTerm(
     const SourceTermConfig& config,
-    const NumLib::LocalToGlobalIndexMap& dof_table, const MeshLib::Mesh& mesh,
-    const int variable_id, const unsigned integration_order,
-    const unsigned shapefunction_order,
+    const NumLib::LocalToGlobalIndexMap& dof_table_bulk,
+    const MeshLib::Mesh& source_term_mesh, const int variable_id,
+    const unsigned integration_order, const unsigned shapefunction_order,
     std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters)
 {
     //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__type}
     auto const type = config.config.peekConfigParameter<std::string>("type");
 
+    // check basic data consistency
+    if (variable_id >=
+            static_cast<int>(dof_table_bulk.getNumberOfVariables()) ||
+        *config.component_id >=
+            dof_table_bulk.getNumberOfVariableComponents(variable_id))
+    {
+        OGS_FATAL(
+            "Variable id or component id too high. Actual values: (%d, "
+            "%d), maximum values: (%d, %d).",
+            variable_id, *config.component_id,
+            dof_table_bulk.getNumberOfVariables(),
+            dof_table_bulk.getNumberOfVariableComponents(variable_id));
+    }
+
+    if (!source_term_mesh.getProperties()
+             .template existsPropertyVector<std::size_t>("bulk_node_ids"))
+    {
+        OGS_FATAL(
+            "The required bulk node ids map does not exist in the "
+            "source term mesh '%s'.",
+            source_term_mesh.getName().c_str());
+    }
+    std::vector<MeshLib::Node*> const& source_term_nodes =
+        source_term_mesh.getNodes();
+    DBUG(
+        "Found %d nodes for source term at mesh '%s' for the variable %d and "
+        "component %d",
+        source_term_nodes.size(), source_term_mesh.getName().c_str(),
+        variable_id, *config.component_id);
+
+    MeshLib::MeshSubset source_term_mesh_subset(source_term_mesh,
+                                                source_term_nodes);
+
+
     if (type == "Nodal")
     {
+        std::unique_ptr<NumLib::LocalToGlobalIndexMap> dof_table_source_term(
+            dof_table_bulk.deriveBoundaryConstrainedMap(
+                variable_id, {*config.component_id},
+                std::move(source_term_mesh_subset)));
         return ProcessLib::createNodalSourceTerm(
-            config.config, config.mesh, dof_table, mesh.getID(), variable_id,
-            *config.component_id, parameters);
+            config.config, config.mesh, std::move(dof_table_source_term),
+            source_term_mesh.getID(), variable_id, *config.component_id,
+            parameters);
     }
 
     if (type == "Volumetric")
     {
+        std::unique_ptr<NumLib::LocalToGlobalIndexMap> dof_table_source_term(
+            dof_table_bulk.deriveBoundaryConstrainedMap(
+                variable_id, {*config.component_id},
+                std::move(source_term_mesh_subset)));
         return ProcessLib::createVolumetricSourceTerm(
-            config.config, config.mesh, dof_table, parameters,
-            integration_order, shapefunction_order, variable_id,
-            *config.component_id);
+            config.config, config.mesh, std::move(dof_table_source_term),
+            parameters, integration_order, shapefunction_order);
     }
 
     if (type == "Python")
     {
 #ifdef OGS_USE_PYTHON
+        std::unique_ptr<NumLib::LocalToGlobalIndexMap> dof_table_source_term(
+            dof_table_bulk.deriveBoundaryConstrainedMap(
+                std::move(source_term_mesh_subset)));
+
         return ProcessLib::createPythonSourceTerm(
-            config.config, config.mesh, dof_table, mesh.getID(), variable_id,
-            *config.component_id, integration_order, shapefunction_order,
-            mesh.getDimension());
+            config.config, config.mesh, std::move(dof_table_source_term),
+            source_term_mesh.getID(), variable_id, *config.component_id,
+            integration_order, shapefunction_order,
+            source_term_mesh.getDimension());
 #else
         OGS_FATAL("OpenGeoSys has not been built with Python support.");
 #endif
diff --git a/ProcessLib/SourceTerms/CreateVolumetricSourceTerm.cpp b/ProcessLib/SourceTerms/CreateVolumetricSourceTerm.cpp
index 93a1f5d0adf090026a4bc787f4b9edc402e1700d..d75500f9e2054acb99bc7a4e28124b77865dbf76 100644
--- a/ProcessLib/SourceTerms/CreateVolumetricSourceTerm.cpp
+++ b/ProcessLib/SourceTerms/CreateVolumetricSourceTerm.cpp
@@ -21,10 +21,9 @@ namespace ProcessLib
 std::unique_ptr<SourceTerm> createVolumetricSourceTerm(
     BaseLib::ConfigTree const& config,
     MeshLib::Mesh const& source_term_mesh,
-    NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap> source_term_dof_table,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
-    unsigned const integration_order, unsigned const shapefunction_order,
-    int const variable_id, int const component_id)
+    unsigned const integration_order, unsigned const shapefunction_order)
 {
     //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__type}
     config.checkConfigParameter("type", "Volumetric");
@@ -42,9 +41,8 @@ std::unique_ptr<SourceTerm> createVolumetricSourceTerm(
          volumetric_source_term.name.c_str());
 
     return std::make_unique<VolumetricSourceTerm>(
-        source_term_mesh, source_term_dof_table, integration_order,
-        shapefunction_order, variable_id, component_id,
-        volumetric_source_term);
+        source_term_mesh, std::move(source_term_dof_table), integration_order,
+        shapefunction_order, volumetric_source_term);
 }
 
 }  // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/CreateVolumetricSourceTerm.h b/ProcessLib/SourceTerms/CreateVolumetricSourceTerm.h
index 2f1fcb62c54f9b02116d790d44a1021d71b7bd9f..e300528c09baca54e22834d1d0dcb155e4ceb908 100644
--- a/ProcessLib/SourceTerms/CreateVolumetricSourceTerm.h
+++ b/ProcessLib/SourceTerms/CreateVolumetricSourceTerm.h
@@ -35,9 +35,8 @@ class SourceTerm;
 std::unique_ptr<SourceTerm> createVolumetricSourceTerm(
     BaseLib::ConfigTree const& config,
     MeshLib::Mesh const& source_term_mesh,
-    NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap> source_term_dof_table,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
-    unsigned const integration_order, unsigned const shapefunction_order,
-    int const variable_id, int const component_id);
+    unsigned const integration_order, unsigned const shapefunction_order);
 
 }   // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/NodalSourceTerm.cpp b/ProcessLib/SourceTerms/NodalSourceTerm.cpp
index 723b342672cfd1c04021b3591759a802c5d5608e..5737be3aacdebb9ad08a4b03a5b7ba1a362d127f 100644
--- a/ProcessLib/SourceTerms/NodalSourceTerm.cpp
+++ b/ProcessLib/SourceTerms/NodalSourceTerm.cpp
@@ -14,27 +14,20 @@
 namespace ProcessLib
 {
 NodalSourceTerm::NodalSourceTerm(
-    const NumLib::LocalToGlobalIndexMap& source_term_dof_table,
-    std::size_t const bulk_mesh_id,
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap> source_term_dof_table,
+    std::size_t const source_term_mesh_id,
     MeshLib::Mesh const& st_mesh,
     const int variable_id,
     const int component_id,
     Parameter<double> const& parameter)
-    : SourceTerm(source_term_dof_table),
-      _bulk_mesh_id(bulk_mesh_id),
+    : SourceTerm(std::move(source_term_dof_table)),
+      _source_term_mesh_id(source_term_mesh_id),
       _st_mesh(st_mesh),
       _variable_id(variable_id),
       _component_id(component_id),
       _parameter(parameter)
 {
     DBUG("Create NodalSourceTerm.");
-    if (!_st_mesh.getProperties().template existsPropertyVector<std::size_t>(
-            "bulk_node_ids"))
-    {
-        OGS_FATAL(
-            "Required mesh property \"bulk_node_ids\" does not exists on the "
-            "source term mesh '%s'.", _st_mesh.getName().c_str());
-    }
 }
 
 void NodalSourceTerm::integrate(const double t, GlobalVector const& /*x*/,
@@ -42,15 +35,12 @@ void NodalSourceTerm::integrate(const double t, GlobalVector const& /*x*/,
 {
     DBUG("Assemble NodalSourceTerm.");
 
-    auto const& bulk_node_ids_map =
-        *_st_mesh.getProperties().template getPropertyVector<std::size_t>(
-            "bulk_node_ids");
     for (MeshLib::Node const* const node : _st_mesh.getNodes())
     {
         auto const node_id = node->getID();
-        MeshLib::Location const l{_bulk_mesh_id, MeshLib::MeshItemType::Node,
-                                  bulk_node_ids_map[node_id]};
-        auto const index = _source_term_dof_table.getGlobalIndex(
+        MeshLib::Location const l{_source_term_mesh_id,
+                                  MeshLib::MeshItemType::Node, node_id};
+        auto const index = _source_term_dof_table->getGlobalIndex(
             l, _variable_id, _component_id);
 
         SpatialPosition pos;
diff --git a/ProcessLib/SourceTerms/NodalSourceTerm.h b/ProcessLib/SourceTerms/NodalSourceTerm.h
index 4df7f58b07fdb9da6cdbb82bbcec93bdd973aa92..9b04ef6527ed14d25ff99e5a11b6c827f9ae8c2b 100644
--- a/ProcessLib/SourceTerms/NodalSourceTerm.h
+++ b/ProcessLib/SourceTerms/NodalSourceTerm.h
@@ -17,8 +17,8 @@ class NodalSourceTerm final : public SourceTerm
 {
 public:
     explicit NodalSourceTerm(
-        const NumLib::LocalToGlobalIndexMap& source_term_dof_table,
-        std::size_t const bulk_mesh_id, MeshLib::Mesh const& st_mesh,
+        std::unique_ptr<NumLib::LocalToGlobalIndexMap> dof_table,
+        std::size_t const source_term_mesh_id, MeshLib::Mesh const& st_mesh,
         const int variable_id, const int component_id,
         Parameter<double> const& parameter);
 
@@ -26,7 +26,7 @@ public:
                    GlobalMatrix* jac) const override;
 
 private:
-    std::size_t const _bulk_mesh_id;
+    std::size_t const _source_term_mesh_id;
     MeshLib::Mesh const& _st_mesh;
     int const _variable_id;
     int const _component_id;
diff --git a/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.cpp b/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.cpp
index 3707a3801e30d63d68c7193e29df0cbfbbd77cb4..e8cb1635fbe7300b2bd9054b2ef1983be1a3a855 100644
--- a/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.cpp
+++ b/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.cpp
@@ -21,7 +21,7 @@ namespace ProcessLib
 {
 std::unique_ptr<SourceTerm> createPythonSourceTerm(
     BaseLib::ConfigTree const& config, MeshLib::Mesh const& source_term_mesh,
-    NumLib::LocalToGlobalIndexMap const& dof_table,
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap> dof_table,
     std::size_t const bulk_mesh_id, int const variable_id,
     int const component_id, unsigned const integration_order,
     unsigned const shapefunction_order, unsigned const global_dim)
@@ -50,16 +50,6 @@ std::unique_ptr<SourceTerm> createPythonSourceTerm(
                             .cast<ProcessLib::SourceTerms::Python::
                                       PythonSourceTermPythonSideInterface*>();
 
-    if (variable_id >= static_cast<int>(dof_table.getNumberOfVariables()) ||
-        component_id >= dof_table.getNumberOfVariableComponents(variable_id))
-    {
-        OGS_FATAL(
-            "Variable id or component id too high. Actual values: (%d, %d), "
-            "maximum values: (%d, %d).",
-            variable_id, component_id, dof_table.getNumberOfVariables(),
-            dof_table.getNumberOfVariableComponents(variable_id));
-    }
-
     // In case of partitioned mesh the source_term could be empty, i.e. there is
     // no source_term condition.
 #ifdef USE_PETSC
@@ -75,12 +65,12 @@ std::unique_ptr<SourceTerm> createPythonSourceTerm(
     }
 #endif  // USE_PETSC
 
+    auto const global_component_id =
+        dof_table->getGlobalComponent(variable_id, component_id);
     return std::make_unique<ProcessLib::SourceTerms::Python::PythonSourceTerm>(
-        dof_table,
+        std::move(dof_table),
         ProcessLib::SourceTerms::Python::PythonSourceTermData{
-            source_term, dof_table, bulk_mesh_id,
-            dof_table.getGlobalComponent(variable_id, component_id),
-            source_term_mesh},
+            source_term, bulk_mesh_id, global_component_id, source_term_mesh},
         integration_order, shapefunction_order, global_dim, flush_stdout);
 }
 
diff --git a/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.h b/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.h
index b7b2d6a4b7b328ba235a38ef7f6e12714b428f1a..910f54bc3bf2e9f3b64f48a4fabdf64377bd32c3 100644
--- a/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.h
+++ b/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.h
@@ -30,7 +30,7 @@ class SourceTerm;
 
 std::unique_ptr<SourceTerm> createPythonSourceTerm(
     BaseLib::ConfigTree const& config, MeshLib::Mesh const& source_term_mesh,
-    NumLib::LocalToGlobalIndexMap const& dof_table,
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap> dof_table,
     std::size_t const bulk_mesh_id, int const variable_id,
     int const component_id, unsigned const integration_order,
     unsigned const shapefunction_order, unsigned const global_dim);
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTerm.cpp b/ProcessLib/SourceTerms/Python/PythonSourceTerm.cpp
index 733dd183140c9761507a12640b2220274651e9e0..aff419c4b87fa7219dddb4e2900cdcda0bb007b9 100644
--- a/ProcessLib/SourceTerms/Python/PythonSourceTerm.cpp
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTerm.cpp
@@ -58,28 +58,17 @@ namespace SourceTerms
 namespace Python
 {
 PythonSourceTerm::PythonSourceTerm(
-    NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap> source_term_dof_table,
     PythonSourceTermData&& source_term_data, unsigned const integration_order,
     unsigned const shapefunction_order, unsigned const global_dim,
     bool const flush_stdout)
-    : SourceTerm(source_term_dof_table),
+    : SourceTerm(std::move(source_term_dof_table)),
       _source_term_data(std::move(source_term_data)),
       _flush_stdout(flush_stdout)
 {
-    std::vector<MeshLib::Node*> const& source_term_nodes =
-        _source_term_data.source_term_mesh.getNodes();
-    MeshLib::MeshSubset source_term_mesh_subset(
-        _source_term_data.source_term_mesh, source_term_nodes);
-
-    // Create local DOF table from the source term mesh subset for the given
-    // variable and component id.
-    _dof_table_source_term =
-        _source_term_data.dof_table_bulk.deriveBoundaryConstrainedMap(
-            std::move(source_term_mesh_subset));
-
     createLocalAssemblers<PythonSourceTermLocalAssembler>(
         global_dim, _source_term_data.source_term_mesh.getElements(),
-        *_dof_table_source_term, shapefunction_order, _local_assemblers,
+        *_source_term_dof_table, shapefunction_order, _local_assemblers,
         _source_term_data.source_term_mesh.isAxiallySymmetric(),
         integration_order, _source_term_data);
 }
@@ -91,7 +80,7 @@ void PythonSourceTerm::integrate(const double t, const GlobalVector& x,
 
     GlobalExecutor::executeMemberOnDereferenced(
         &PythonSourceTermLocalAssemblerInterface::assemble, _local_assemblers,
-        *_dof_table_source_term, t, x, b, Jac);
+        *_source_term_dof_table, t, x, b, Jac);
 }
 
 }  // namespace Python
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTerm.h b/ProcessLib/SourceTerms/Python/PythonSourceTerm.h
index 9ec77f3952f533ed64171d7a508f735d6256654c..cac510ece8d4e2a8a29c6f3d69089cbbbc4fdfdf 100644
--- a/ProcessLib/SourceTerms/Python/PythonSourceTerm.h
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTerm.h
@@ -30,9 +30,6 @@ struct PythonSourceTermData final
     //! Python object computing source term values.
     PythonSourceTermPythonSideInterface* source_term_object;
 
-    //! DOF table of the entire domain.
-    NumLib::LocalToGlobalIndexMap const& dof_table_bulk;
-
     //! Mesh ID of the entire domain.
     std::size_t const bulk_mesh_id;
 
@@ -49,7 +46,7 @@ class PythonSourceTerm final : public ProcessLib::SourceTerm
 {
 public:
     explicit PythonSourceTerm(
-        NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+        std::unique_ptr<NumLib::LocalToGlobalIndexMap> source_term_dof_table,
         PythonSourceTermData&& source_term_data,
         unsigned const integration_order, unsigned const shapefunction_order,
         unsigned const global_dim, bool const flush_stdout);
@@ -61,9 +58,6 @@ private:
     //! Auxiliary data.
     PythonSourceTermData _source_term_data;
 
-    //! Local dof table for the boundary mesh.
-    std::unique_ptr<NumLib::LocalToGlobalIndexMap> _dof_table_source_term;
-
     //! Local assemblers for all elements of the source term mesh.
     std::vector<std::unique_ptr<PythonSourceTermLocalAssemblerInterface>>
         _local_assemblers;
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h b/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h
index 9a8039ad6ce99ca501a78ba53069e4eb54b05da3..ad569ebeaccce5c6a804162c73a758047d22efa8 100644
--- a/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h
@@ -99,10 +99,10 @@ public:
 
         unsigned const num_integration_points =
             _integration_method.getNumberOfPoints();
-        auto const num_var = _data.dof_table_bulk.getNumberOfVariables();
+        auto const num_var = dof_table_source_term.getNumberOfVariables();
         auto const num_nodes = ShapeFunction::NPOINTS;
         auto const num_comp_total =
-            _data.dof_table_bulk.getNumberOfComponents();
+            dof_table_source_term.getNumberOfComponents();
 
         auto const& bulk_node_ids_map =
             *_data.source_term_mesh.getProperties()
@@ -115,7 +115,7 @@ public:
         for (int var = 0; var < num_var; ++var)
         {
             auto const num_comp =
-                _data.dof_table_bulk.getNumberOfVariableComponents(var);
+                dof_table_source_term.getNumberOfVariableComponents(var);
             for (int comp = 0; comp < num_comp; ++comp)
             {
                 auto const global_component =
@@ -132,7 +132,7 @@ public:
                                           MeshLib::MeshItemType::Node,
                                           bulk_node_id};
                     auto const dof_idx =
-                        _data.dof_table_bulk.getGlobalIndex(loc, var, comp);
+                        dof_table_source_term.getGlobalIndex(loc, var, comp);
                     if (dof_idx == NumLib::MeshComponentMap::nop)
                     {
                         // TODO extend Python BC to mixed FEM ansatz functions
diff --git a/ProcessLib/SourceTerms/SourceTerm.h b/ProcessLib/SourceTerms/SourceTerm.h
index de87c32c8431ef95ef70a44a6345fe5f86c1e9da..9de42d55336e027f1030278dc3d43a8c8a8505f9 100644
--- a/ProcessLib/SourceTerms/SourceTerm.h
+++ b/ProcessLib/SourceTerms/SourceTerm.h
@@ -9,6 +9,8 @@
 
 #pragma once
 
+#include <memory>
+
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
 #include "ProcessLib/Parameter/Parameter.h"
 
@@ -18,8 +20,8 @@ class SourceTerm
 {
 public:
     explicit SourceTerm(
-        const NumLib::LocalToGlobalIndexMap& source_term_dof_table)
-        : _source_term_dof_table(source_term_dof_table)
+        std::unique_ptr<NumLib::LocalToGlobalIndexMap> source_term_dof_table)
+        : _source_term_dof_table{std::move(source_term_dof_table)}
     {
     }
 
@@ -29,7 +31,7 @@ public:
     virtual ~SourceTerm() = default;
 
 protected:
-    NumLib::LocalToGlobalIndexMap const& _source_term_dof_table;
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap> const _source_term_dof_table;
 };
 
 }  // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/VolumetricSourceTerm.cpp b/ProcessLib/SourceTerms/VolumetricSourceTerm.cpp
index c05f41bdd3d342efbc027c66729ccd0d776ecda9..785ca74e9576920462064896c734522a152ad3b9 100644
--- a/ProcessLib/SourceTerms/VolumetricSourceTerm.cpp
+++ b/ProcessLib/SourceTerms/VolumetricSourceTerm.cpp
@@ -15,34 +15,15 @@ namespace ProcessLib
 {
 VolumetricSourceTerm::VolumetricSourceTerm(
     MeshLib::Mesh const& source_term_mesh,
-    NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap> source_term_dof_table,
     unsigned const integration_order, unsigned const shapefunction_order,
-    int const variable_id, int const component_id,
     Parameter<double> const& volumetric_source_term)
-    : SourceTerm(source_term_dof_table),
+    : SourceTerm(std::move(source_term_dof_table)),
       _volumetric_source_term(volumetric_source_term)
 {
-    // check basic data consistency
-    if (variable_id >=
-        static_cast<int>(source_term_dof_table.getNumberOfVariables()))
-    {
-        OGS_FATAL(
-            "Variable id too high. Actual value: %d, maximum value: %d.",
-            variable_id,
-            source_term_dof_table.getNumberOfVariables());
-    }
-    if (component_id >=
-        source_term_dof_table.getNumberOfVariableComponents(variable_id))
-    {
-        OGS_FATAL(
-            "Component id too high. Actual value: %d, maximum value: %d.",
-            component_id,
-            source_term_dof_table.getNumberOfVariableComponents(variable_id));
-    }
-
     ProcessLib::createLocalAssemblers<VolumetricSourceTermLocalAssembler>(
         source_term_mesh.getDimension(), source_term_mesh.getElements(),
-        source_term_dof_table, shapefunction_order, _local_assemblers,
+        *_source_term_dof_table, shapefunction_order, _local_assemblers,
         source_term_mesh.isAxiallySymmetric(), integration_order,
         _volumetric_source_term);
 }
@@ -56,7 +37,7 @@ void VolumetricSourceTerm::integrate(const double t, GlobalVector const& /*x*/,
     // Call global assembler for each local assembly item.
     GlobalExecutor::executeMemberOnDereferenced(
         &VolumetricSourceTermLocalAssemblerInterface::integrate,
-        _local_assemblers, _source_term_dof_table, t, b);
+        _local_assemblers, *_source_term_dof_table, t, b);
 }
 
 }   // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/VolumetricSourceTerm.h b/ProcessLib/SourceTerms/VolumetricSourceTerm.h
index f3d785f310199d87d2cc086e72acaa07790089a7..8225a2dbda9dfb51c0c9eb98592b60f0d6d35569 100644
--- a/ProcessLib/SourceTerms/VolumetricSourceTerm.h
+++ b/ProcessLib/SourceTerms/VolumetricSourceTerm.h
@@ -22,9 +22,8 @@ class VolumetricSourceTerm final : public SourceTerm
 public:
     VolumetricSourceTerm(
         MeshLib::Mesh const& source_term_mesh,
-        NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+        std::unique_ptr<NumLib::LocalToGlobalIndexMap> source_term_dof_table,
         unsigned const integration_order, unsigned const shapefunction_order,
-        int const variable_id, int const component_id,
         Parameter<double> const& volumetric_source_term);
 
     void integrate(const double t, GlobalVector const& x, GlobalVector& b,
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e2_volumetricsourceterm.prj b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e2_volumetricsourceterm.prj
index e673cb7de144b71f7a93fca9bc6c3db77275d0cb..8bdf93d1883aee1f03a3ea278c88f8f8409e8ecc 100644
--- a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e2_volumetricsourceterm.prj
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e2_volumetricsourceterm.prj
@@ -2,6 +2,7 @@
 <OpenGeoSysProject>
     <meshes>
         <mesh>square_1x1_quad_1e2.vtu</mesh>
+        <mesh>square_1x1_quad_1e2_volumetric_source_term.vtu</mesh>
         <mesh>square_1x1_quad_1e2_geometry_right.vtu</mesh>
         <mesh>square_1x1_quad_1e2_geometry_left.vtu</mesh>
     </meshes>
@@ -118,7 +119,7 @@
             </boundary_conditions>
             <source_terms>
                 <source_term>
-                    <mesh>square_1x1_quad_1e2</mesh>
+                    <mesh>square_1x1_quad_1e2_volumetric_source_term</mesh>
                     <type>Volumetric</type>
                     <parameter>volumetric_source_term_parameter</parameter>
                 </source_term>
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e3_volumetricsourceterm.prj b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e3_volumetricsourceterm.prj
index 23892db304c4d425f47e5de5e62d20f9460653ef..5a694b93262f6a83d7a748158224f58481c64064 100644
--- a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e3_volumetricsourceterm.prj
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e3_volumetricsourceterm.prj
@@ -2,6 +2,7 @@
 <OpenGeoSysProject>
     <meshes>
         <mesh>square_1x1_quad_1e3.vtu</mesh>
+        <mesh>square_1x1_quad_1e3_volumetric_source_term.vtu</mesh>
         <mesh>square_1x1_quad_1e3_geometry_right.vtu</mesh>
         <mesh>square_1x1_quad_1e3_geometry_left.vtu</mesh>
     </meshes>
@@ -118,7 +119,7 @@
             </boundary_conditions>
             <source_terms>
                 <source_term>
-                    <mesh>square_1x1_quad_1e3</mesh>
+                    <mesh>square_1x1_quad_1e3_volumetric_source_term</mesh>
                     <type>Volumetric</type>
                     <parameter>volumetric_source_term_parameter</parameter>
                 </source_term>
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e3_volumetricsourcetermdataarray.prj b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e3_volumetricsourcetermdataarray.prj
index 5d3737ffe45db317997ae21aebc28c4495a5a0dd..dacd8025981c1997e9e4a225c9bc159f754d72df 100644
--- a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e3_volumetricsourcetermdataarray.prj
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e3_volumetricsourcetermdataarray.prj
@@ -2,6 +2,7 @@
 <OpenGeoSysProject>
     <meshes>
         <mesh>square_1x1_quad_1e3_volumetricsourcetermdataarray.vtu</mesh>
+        <mesh>square_1x1_quad_1e3_volumetricsourcetermdataarray_volumetric_source_term.vtu</mesh>
         <mesh>square_1x1_quad_1e3_volumetricsourcetermdataarray_geometry_ll.vtu</mesh>
         <mesh>square_1x1_quad_1e3_volumetricsourcetermdataarray_geometry_lr.vtu</mesh>
         <mesh>square_1x1_quad_1e3_volumetricsourcetermdataarray_geometry_ul.vtu</mesh>
@@ -125,7 +126,7 @@
             </boundary_conditions>
             <source_terms>
                 <source_term>
-                    <mesh>square_1x1_quad_1e3_volumetricsourcetermdataarray</mesh>
+                    <mesh>square_1x1_quad_1e3_volumetricsourcetermdataarray_volumetric_source_term</mesh>
                     <type>Volumetric</type>
                     <parameter>volumetric_source_term_parameter</parameter>
                 </source_term>
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e6_with_nodal_sources.prj b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e6_with_nodal_sources.prj
index b33d408f85cce5d4622079f44d3d4204081e6765..4bc7475277d837ba7bbc96fce3321a3c944ad13c 100644
--- a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e6_with_nodal_sources.prj
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1e6_with_nodal_sources.prj
@@ -1,7 +1,11 @@
 <?xml version="1.0" encoding="ISO-8859-1"?>
 <OpenGeoSysProject>
-    <mesh>square_1x1_quad_1e6.vtu</mesh>
-    <geometry>square_1x1_midpoint.gml</geometry>
+    <meshes>
+        <mesh>square_1x1_quad_1e6.vtu</mesh>
+        <mesh>square_1x1_geometry_middle_point.vtu</mesh>
+        <mesh>square_1x1_geometry_left.vtu</mesh>
+        <mesh>square_1x1_geometry_right.vtu</mesh>
+    </meshes>
     <processes>
         <process>
             <name>GW23</name>
@@ -79,22 +83,19 @@
             <initial_condition>p0</initial_condition>
             <boundary_conditions>
                 <boundary_condition>
-                    <geometrical_set>square_1x1_geometry</geometrical_set>
-                    <geometry>left</geometry>
+                    <mesh>square_1x1_geometry_left</mesh>
                     <type>Dirichlet</type>
                     <parameter>p_Dirichlet_left</parameter>
                 </boundary_condition>
                 <boundary_condition>
-                    <geometrical_set>square_1x1_geometry</geometrical_set>
-                    <geometry>right</geometry>
+                    <mesh>square_1x1_geometry_right</mesh>
                     <type>Dirichlet</type>
                     <parameter>p_Dirichlet_right</parameter>
                 </boundary_condition>
             </boundary_conditions>
             <source_terms>
                 <source_term>
-                    <geometrical_set>square_1x1_geometry</geometrical_set>
-                    <geometry>middle_point</geometry>
+                    <mesh>square_1x1_geometry_middle_point</mesh>
                     <type>Nodal</type>
                     <parameter>pressure_source_term</parameter>
                 </source_term>
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_geometry_left.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_geometry_left.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..a13635eec703115fcd1cdd8f7e205c09fa3df626
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_geometry_left.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:738bf064eb885af4ab9d1e8a59489be3b5c2f0e87338b6975ad86548a656976e
+size 87912
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_geometry_middle_point.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_geometry_middle_point.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..243f3dc36caf914516c651435afaf1edb5ff9c9f
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_geometry_middle_point.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:dfdb0edd85c31b5007e94245779ec7cdce9f3db606363e5651b829b7df5918b0
+size 1297
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_geometry_right.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_geometry_right.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..a89ebce793b0d69b7daa8d421b58cfc5580dc1f8
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_geometry_right.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7dc4ea2f8b696874ab4a1437a79479f767cddf736d1e09dd19be207e636b7834
+size 87928
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_quad_1e2_volumetric_source_term.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_quad_1e2_volumetric_source_term.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..c740d4c3bcee30a5710cb11d3117f50112503f11
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_quad_1e2_volumetric_source_term.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:893bcc7f201daccd32ceea681a4d53c90062dc8164ced81d61450cb6c980f1d1
+size 16440
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_quad_1e3_volumetric_source_term.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_quad_1e3_volumetric_source_term.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..9cdee33c196515abeafc61df3adb4a59f9ca6688
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_quad_1e3_volumetric_source_term.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5cb72aebd16f1ed3e2194639df981592801f2af6a6a5b946ae294c35ed51d6ba
+size 143703
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_quad_1e3_volumetricsourcetermdataarray_volumetric_source_term.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_quad_1e3_volumetricsourcetermdataarray_volumetric_source_term.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..b8a8cc2616665243bfb392e6c0106d4378ec2023
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow/square_1x1_quad_1e3_volumetricsourcetermdataarray_volumetric_source_term.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:3b06d86b4b608361e25a5858542e1ad9aec75f34b1d242abf30dba2e28d647e0
+size 178400
diff --git a/Tests/Data/Parabolic/T/1D_dirichlet_source-term/SourceTermDomain.vtu b/Tests/Data/Parabolic/T/1D_dirichlet_source-term/SourceTermDomain.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..5c390d1e2931494609ffef75895040616c0687f9
--- /dev/null
+++ b/Tests/Data/Parabolic/T/1D_dirichlet_source-term/SourceTermDomain.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f2990eedb2204d4bdc58002902b551e4bafc3ee186bf2ed4c05dfb4b1cf4907f
+size 1407
diff --git a/Tests/Data/Parabolic/T/1D_dirichlet_source-term/line_1_geometry_left.vtu b/Tests/Data/Parabolic/T/1D_dirichlet_source-term/line_1_geometry_left.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..2acfaf9fb2956770259901a765efd8d9d4938718
--- /dev/null
+++ b/Tests/Data/Parabolic/T/1D_dirichlet_source-term/line_1_geometry_left.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:612e6cc7d6f49da22c69cf3b4987541f2c40c77e3abe14db8af4b20f8e966742
+size 1253
diff --git a/Tests/Data/Parabolic/T/1D_dirichlet_source-term/line_1_geometry_right.vtu b/Tests/Data/Parabolic/T/1D_dirichlet_source-term/line_1_geometry_right.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..59276ed37cede25d847e42cdd54ca2ee8ba7fc76
--- /dev/null
+++ b/Tests/Data/Parabolic/T/1D_dirichlet_source-term/line_1_geometry_right.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fc509d0b26cc4a0af1e1f92bf7941db2b9cc7e31819ed6349bd5ed1796d83c88
+size 1259
diff --git a/Tests/Data/Parabolic/T/1D_dirichlet_source-term/line_1_line_1e2_pcs_0_ts_500_t_39062500.000000_reference.vtu b/Tests/Data/Parabolic/T/1D_dirichlet_source-term/line_1_line_1e2_pcs_0_ts_500_t_39062500.000000_reference.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..7245032601d3e3efac54bf3b3443fdf51d5a4003
--- /dev/null
+++ b/Tests/Data/Parabolic/T/1D_dirichlet_source-term/line_1_line_1e2_pcs_0_ts_500_t_39062500.000000_reference.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ceb89b82bc5b0c6ede77a01da677c90a5f43aae6c1eb49ea61e72e4f64b44264
+size 3816
diff --git a/Tests/Data/Parabolic/T/1D_dirichlet_source-term/line_1_line_1e2_source_term.prj b/Tests/Data/Parabolic/T/1D_dirichlet_source-term/line_1_line_1e2_source_term.prj
new file mode 100644
index 0000000000000000000000000000000000000000..d92919aae4bc534a8bf2424b575a03b541887efa
--- /dev/null
+++ b/Tests/Data/Parabolic/T/1D_dirichlet_source-term/line_1_line_1e2_source_term.prj
@@ -0,0 +1,150 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <meshes>
+        <mesh>line_1_lines_1e2.vtu</mesh>
+        <mesh>SourceTermDomain.vtu</mesh>
+        <mesh>line_1_geometry_left.vtu</mesh>
+        <mesh>line_1_geometry_right.vtu</mesh>
+    </meshes>
+    <processes>
+        <process>
+            <name>HeatConduction</name>
+            <type>HEAT_CONDUCTION</type>
+            <integration_order>2</integration_order>
+            <thermal_conductivity>K</thermal_conductivity>
+            <heat_capacity>Cp</heat_capacity>
+            <density>rho</density>
+            <process_variables>
+                <process_variable>temperature</process_variable>
+            </process_variables>
+            <secondary_variables>
+                <secondary_variable type="static" internal_name="heat_flux_x" output_name="heat_flux_x"/>
+            </secondary_variables>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="HeatConduction">
+                <nonlinear_solver>basic_picard</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <abstol>1.e-6</abstol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <output>
+                    <variables>
+                        <variable> temperature </variable>
+                        <variable> heat_flux_x </variable>
+                    </variables>
+                </output>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial> 0.0 </t_initial>
+                    <t_end> 39062500 </t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>500</repeat>
+                            <delta_t>78125</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>line_1_line_1e2</prefix>
+            <timesteps>
+                <pair>
+                    <repeat> 1 </repeat>
+                    <each_steps> 500 </each_steps>
+                </pair>
+            </timesteps>
+        </output>
+    </time_loop>
+    <parameters>
+        <parameter>
+            <name>K</name>
+            <type>Constant</type>
+            <value>5</value>
+        </parameter>
+        <parameter>
+            <name>Cp</name>
+            <type>Constant</type>
+            <value>1000</value>
+        </parameter>
+        <parameter>
+            <name>rho</name>
+            <type>Constant</type>
+            <value>2500</value>
+        </parameter>
+        <parameter>
+            <name>T0</name>
+            <type>Constant</type>
+            <value>273.15</value>
+        </parameter>
+        <parameter>
+            <name>T1</name>
+            <type>Constant</type>
+            <value>274.15</value>
+        </parameter>
+        <parameter>
+            <name>st_param</name>
+            <type>Constant</type>
+            <value>5</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>temperature</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>T0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>line_1_geometry_left</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>T1</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <mesh>line_1_geometry_right</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>T1</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+            <source_terms>
+                <source_term>
+                    <mesh>SourceTermDomain</mesh>
+                    <type>Volumetric</type>
+                    <parameter>st_param</parameter>
+                </source_term>
+            </source_terms>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_picard</name>
+            <type>Picard</type>
+            <max_iter>10</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+            <eigen>
+                <solver_type>CG</solver_type>
+                <precon_type>DIAGONAL</precon_type>
+                <max_iteration_step>10000</max_iteration_step>
+                <error_tolerance>1e-16</error_tolerance>
+            </eigen>
+            <petsc>
+                <prefix>gw</prefix>
+                <parameters>-gw_ksp_type cg -gw_pc_type bjacobi -gw_ksp_rtol 1e-16 -gw_ksp_max_it 10000</parameters>
+            </petsc>
+        </linear_solver>
+    </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Parabolic/T/1D_dirichlet_source-term/line_1_lines_1e2.vtu b/Tests/Data/Parabolic/T/1D_dirichlet_source-term/line_1_lines_1e2.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..fa08f9e37f628254798978e37c6c01b458a29cb3
--- /dev/null
+++ b/Tests/Data/Parabolic/T/1D_dirichlet_source-term/line_1_lines_1e2.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:44ce4d6db0975e9fb108401c3b80bed86473ff06fa7a3d92f895c8a567392430
+size 7173