diff --git a/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.cpp b/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.cpp
index 17d0eef15c71a42246240c052755b4d0757ee79a..bfccdecf8783e55d39c69dd24b7f8df0e48ea88e 100644
--- a/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.cpp
+++ b/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.cpp
@@ -338,37 +338,18 @@ OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeTet4);
 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeTri3);
 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeTri6);
 
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeHex20, 1);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeHex8, 1);
+OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePoint1, 0);
+
 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeLine2, 1);
 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeLine3, 1);
 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePoint1, 1);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePrism15, 1);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePrism6, 1);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePyra13, 1);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePyra5, 1);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad4, 1);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad8, 1);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad9, 1);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTet10, 1);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTet4, 1);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTri3, 1);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTri6, 1);
 
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeHex20, 2);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeHex8, 2);
 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeLine2, 2);
 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeLine3, 2);
 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePoint1, 2);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePrism15, 2);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePrism6, 2);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePyra13, 2);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePyra5, 2);
 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad4, 2);
 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad8, 2);
 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad9, 2);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTet10, 2);
-OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTet4, 2);
 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTri3, 2);
 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTri6, 2);
 
diff --git a/ProcessLib/BoundaryConditionAndSourceTerm/Python/Utils/CreateLocalAssemblers.h b/ProcessLib/BoundaryConditionAndSourceTerm/Python/Utils/CreateLocalAssemblers.h
index 0960b713c2e507b574dccabbf536c0b082e15a2f..e4a629d8913ed0382d971ad3df74a133bff7fc72 100644
--- a/ProcessLib/BoundaryConditionAndSourceTerm/Python/Utils/CreateLocalAssemblers.h
+++ b/ProcessLib/BoundaryConditionAndSourceTerm/Python/Utils/CreateLocalAssemblers.h
@@ -33,10 +33,6 @@ void createLocalAssemblersPython(
     NumLib::IntegrationOrder const integration_order,
     ExtraCtorArgs&&... extra_ctor_args)
 {
-    static_assert(
-        GlobalDim == 1 || GlobalDim == 2 || GlobalDim == 3,
-        "Meshes with dimension greater than three are not supported.");
-
     using LocAsmFactory =
         LocalAssemblerFactoryPython<LocalAssemblerInterface,
                                     LocalAssemblerImplementation, GlobalDim,
@@ -82,6 +78,12 @@ void createLocalAssemblersPython(
 
     switch (dimension)
     {
+        case 0:
+            detail::createLocalAssemblersPython<0,
+                                                LocalAssemblerImplementation>(
+                dof_table, mesh_elements, local_assemblers, integration_order,
+                std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+            break;
         case 1:
             detail::createLocalAssemblersPython<1,
                                                 LocalAssemblerImplementation>(
diff --git a/ProcessLib/BoundaryConditionAndSourceTerm/Python/Utils/NsAndWeight.h b/ProcessLib/BoundaryConditionAndSourceTerm/Python/Utils/NsAndWeight.h
index d6cb88eecaef2eb0e5a2bea9bbbd1fcd7a0468a5..7cbede4e660317c960ea277703ff795b50b02e51 100644
--- a/ProcessLib/BoundaryConditionAndSourceTerm/Python/Utils/NsAndWeight.h
+++ b/ProcessLib/BoundaryConditionAndSourceTerm/Python/Utils/NsAndWeight.h
@@ -91,13 +91,19 @@ struct NsAndWeight<ShapeFunction, ShapeFunction, ShapeMatrix, ShapeMatrix>
 
     Eigen::Ref<const Eigen::RowVectorXd> N(unsigned order) const
     {
-        if (order >= 2)
+        // For point elements shape functions are the same for all orders, so
+        // this specialization will be selected, which is OK for any shape
+        // function order for point elements.
+        if constexpr (!std::is_same_v<ShapeFunction, NumLib::ShapePoint1>)
         {
-            OGS_FATAL(
-                "Only shape functions of order < 2 are available in the "
-                "current setting. You have requested order {}. Maybe there is "
-                "an error in the OGS project file.",
-                order);
+            if (order >= 2)
+            {
+                OGS_FATAL(
+                    "Only shape functions of order < 2 are available in the "
+                    "current setting. You have requested order {}. Maybe there "
+                    "is an error in the OGS project file.",
+                    order);
+            }
         }
 
         return N_;
diff --git a/ProcessLib/BoundaryConditionAndSourceTerm/Utils/CreateLocalAssemblers.h b/ProcessLib/BoundaryConditionAndSourceTerm/Utils/CreateLocalAssemblers.h
index 22b2b005b1057adba78b2437685a016c50d06a04..1d462cc294b80f5504af12dfc2e8f1a1612bff61 100644
--- a/ProcessLib/BoundaryConditionAndSourceTerm/Utils/CreateLocalAssemblers.h
+++ b/ProcessLib/BoundaryConditionAndSourceTerm/Utils/CreateLocalAssemblers.h
@@ -33,10 +33,6 @@ void createLocalAssemblers(
     NumLib::IntegrationOrder const integration_order,
     ExtraCtorArgs&&... extra_ctor_args)
 {
-    static_assert(
-        GlobalDim == 1 || GlobalDim == 2 || GlobalDim == 3,
-        "Meshes with dimension greater than three are not supported.");
-
     using LocalAssemblerFactory =
         LocalAssemblerFactory<LocalAssemblerInterface,
                               LocalAssemblerImplementation, GlobalDim,
@@ -83,6 +79,12 @@ void createLocalAssemblers(
 
     switch (dimension)
     {
+        case 0:
+            detail::createLocalAssemblers<0, LocalAssemblerImplementation>(
+                dof_table, shapefunction_order, mesh_elements, local_assemblers,
+                integration_order,
+                std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+            break;
         case 1:
             detail::createLocalAssemblers<1, LocalAssemblerImplementation>(
                 dof_table, shapefunction_order, mesh_elements, local_assemblers,
diff --git a/ProcessLib/ThermoHydroMechanics/Tests.cmake b/ProcessLib/ThermoHydroMechanics/Tests.cmake
index e8632084e92a4d0ef8da955e31e264339efb6e1e..8ebcf7e202899fc9784c64307c1778e8ac2a4f8e 100644
--- a/ProcessLib/ThermoHydroMechanics/Tests.cmake
+++ b/ProcessLib/ThermoHydroMechanics/Tests.cmake
@@ -125,6 +125,40 @@ AddTest(
     expected_pointheatsource_linear-mesh_ts_10_t_50000.000000.vtu pointheatsource_linear-mesh_ts_10_t_50000.000000.vtu epsilon epsilon 1e-5 1e-5
     expected_pointheatsource_linear-mesh_ts_10_t_50000.000000.vtu pointheatsource_linear-mesh_ts_10_t_50000.000000.vtu sigma sigma 1e-5 1e-5
 )
+# ThermoHydroMechanics; Small deformation, linear poroelastic, point heat source consolidation, variation with a volumetric source term
+AddTest(
+    NAME ThermoHydroMechanics_point_heat_injection_vol_st
+    PATH ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term
+    RUNTIME 45
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS pointheatsource_quadratic-mesh.xml
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    DIFF_DATA
+    expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu displacement displacement 1e-5 1e-5
+    expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu pressure pressure 1e-5 1e-5
+    expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu temperature temperature 1e-5 1e-5
+    expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu epsilon epsilon 1e-5 1e-5
+    expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu sigma sigma 1e-5 1e-5
+)
+# ThermoHydroMechanics; Small deformation, linear poroelastic, point heat source consolidation, variation with a Python source term
+AddTest(
+    NAME ThermoHydroMechanics_point_heat_injection_py_st
+    PATH ThermoHydroMechanics/Linear/Point_injection/with-python-source-term
+    RUNTIME 45
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS pointheatsource_quadratic-mesh.xml
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    DIFF_DATA
+    expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu displacement displacement 1e-5 1e-5
+    expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu pressure pressure 1e-5 1e-5
+    expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu temperature temperature 1e-5 1e-5
+    expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu epsilon epsilon 1e-5 1e-5
+    expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu sigma sigma 1e-5 1e-5
+)
 # ThermoHydroMechanics; Small deformation, linear elastic, porosity=0, anisotropic thermal expansion
 AddTest(
     NAME ThermoHydroMechanics_cube_ortho-thermal-expansion-phi0
diff --git a/ProcessLib/ThermoRichardsMechanics/Tests.cmake b/ProcessLib/ThermoRichardsMechanics/Tests.cmake
index 7cd0815d23805a3ecfd37ccfcc87d63f9be609ed..5dfac37ec31ebc19a7d69027464bf262559d69e7 100644
--- a/ProcessLib/ThermoRichardsMechanics/Tests.cmake
+++ b/ProcessLib/ThermoRichardsMechanics/Tests.cmake
@@ -205,7 +205,7 @@ AddTest(
     expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu PointHeatSource_ts_10_t_50000.000000.vtu NodalForces NodalForces 3.6e2 0
     # submesh residuum output quarter_002_2nd_r_gt_2
     PointHeatSource_quarter_002_2nd_r_gt_2_ts_10_t_50000.000000.vtu PointHeatSource_quarter_002_2nd_r_gt_2_ts_10_t_50000.000000.vtu displacement displacement 1e-15 0
-    PointHeatSource_quarter_002_2nd_r_gt_2_ts_10_t_50000.000000.vtu PointHeatSource_quarter_002_2nd_r_gt_2_ts_10_t_50000.000000.vtu pressure pressure 8.2e-8 0
+    PointHeatSource_quarter_002_2nd_r_gt_2_ts_10_t_50000.000000.vtu PointHeatSource_quarter_002_2nd_r_gt_2_ts_10_t_50000.000000.vtu pressure pressure 8.2e-8 2.2e-10
     PointHeatSource_quarter_002_2nd_r_gt_2_ts_10_t_50000.000000.vtu PointHeatSource_quarter_002_2nd_r_gt_2_ts_10_t_50000.000000.vtu temperature temperature 1.2e-13 0
     PointHeatSource_quarter_002_2nd_r_gt_2_ts_10_t_50000.000000.vtu PointHeatSource_quarter_002_2nd_r_gt_2_ts_10_t_50000.000000.vtu epsilon epsilon 1e-15 0
     PointHeatSource_quarter_002_2nd_r_gt_2_ts_10_t_50000.000000.vtu PointHeatSource_quarter_002_2nd_r_gt_2_ts_10_t_50000.000000.vtu sigma sigma 7e-8 0
@@ -218,7 +218,7 @@ AddTest(
     PointHeatSource_quarter_002_2nd_r_lt_2_ts_10_t_50000.000000.vtu PointHeatSource_quarter_002_2nd_r_lt_2_ts_10_t_50000.000000.vtu temperature temperature 3e-11 0
     PointHeatSource_quarter_002_2nd_r_lt_2_ts_10_t_50000.000000.vtu PointHeatSource_quarter_002_2nd_r_lt_2_ts_10_t_50000.000000.vtu epsilon epsilon 1e-15 0
     PointHeatSource_quarter_002_2nd_r_lt_2_ts_10_t_50000.000000.vtu PointHeatSource_quarter_002_2nd_r_lt_2_ts_10_t_50000.000000.vtu sigma sigma 1.1e-5 0
-    PointHeatSource_quarter_002_2nd_r_lt_2_ts_10_t_50000.000000.vtu PointHeatSource_quarter_002_2nd_r_lt_2_ts_10_t_50000.000000.vtu HeatFlowRate HeatFlowRate 6.2e-12 0
+    PointHeatSource_quarter_002_2nd_r_lt_2_ts_10_t_50000.000000.vtu PointHeatSource_quarter_002_2nd_r_lt_2_ts_10_t_50000.000000.vtu HeatFlowRate HeatFlowRate 6.3e-12 0
     PointHeatSource_quarter_002_2nd_r_lt_2_ts_10_t_50000.000000.vtu PointHeatSource_quarter_002_2nd_r_lt_2_ts_10_t_50000.000000.vtu MassFlowRate MassFlowRate 1e-15 0
     PointHeatSource_quarter_002_2nd_r_lt_2_ts_10_t_50000.000000.vtu PointHeatSource_quarter_002_2nd_r_lt_2_ts_10_t_50000.000000.vtu NodalForces NodalForces 2.1e-8 0
 )
diff --git a/Tests/Data/TH2M/submesh_residuum_assembly/p_cap.xml b/Tests/Data/TH2M/submesh_residuum_assembly/p_cap.xml
index bb662ee5e18199ecf339fcb7858de52a7f43a087..c28e26f8232408ff7c0be8f8930fbdbfa01c99b9 100644
--- a/Tests/Data/TH2M/submesh_residuum_assembly/p_cap.xml
+++ b/Tests/Data/TH2M/submesh_residuum_assembly/p_cap.xml
@@ -60,7 +60,7 @@
     <add sel="/*/time_loop/submesh_residuum_output">
         <prefix>p_cap_{:meshname}</prefix>
     </add>
-    
+
     <add sel="/*/time_loop/output/variables">
         <variable>HeatFlowRate</variable>
         <variable>GasMassFlowRate</variable>
@@ -88,7 +88,7 @@
                 <absolute_tolerance>1e-15</absolute_tolerance>
                 <relative_tolerance>0</relative_tolerance>
             </vtkdiff>
-            
+
             <vtkdiff>
                 <regex>p_cap.*_ts_[0-9]*_.*\.vtu</regex>
                 <field>gas_pressure_interpolated</field>
@@ -101,11 +101,11 @@
                 <absolute_tolerance>1e-15</absolute_tolerance>
                 <relative_tolerance>0</relative_tolerance>
             </vtkdiff>
-            
+
             <vtkdiff>
                 <regex>p_cap.*_ts_[0-9]*_.*\.vtu</regex>
                 <field>capillary_pressure_interpolated</field>
-                <absolute_tolerance>9e-11</absolute_tolerance>
+                <absolute_tolerance>1.2e-10</absolute_tolerance>
                 <relative_tolerance>0</relative_tolerance>
             </vtkdiff>
             <vtkdiff>
@@ -114,7 +114,7 @@
                 <absolute_tolerance>1e-15</absolute_tolerance>
                 <relative_tolerance>0</relative_tolerance>
             </vtkdiff>
-            
+
             <vtkdiff>
                 <regex>p_cap.*_ts_[0-9]*_.*\.vtu</regex>
                 <field>displacement</field>
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/HeatTransportInStationaryFlow/WithFreezingPhase.prj b/Tests/Data/ThermoHydroMechanics/Linear/HeatTransportInStationaryFlow/WithFreezingPhase.prj
index 5f29ac8a8d1db0b594a8f6b48a2e49e79e631a4c..10c859b5bd268746a4b6bfdc1bda1bd8a41d6e96 100644
--- a/Tests/Data/ThermoHydroMechanics/Linear/HeatTransportInStationaryFlow/WithFreezingPhase.prj
+++ b/Tests/Data/ThermoHydroMechanics/Linear/HeatTransportInStationaryFlow/WithFreezingPhase.prj
@@ -378,7 +378,7 @@
         <vtkdiff>
             <regex>WithFreezingPhase_ts_.*.vtu</regex>
             <field>temperature</field>
-            <absolute_tolerance>1e-12</absolute_tolerance>
+            <absolute_tolerance>1.2e-12</absolute_tolerance>
             <relative_tolerance>0</relative_tolerance>
         </vtkdiff>
         <vtkdiff>
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu
new file mode 120000
index 0000000000000000000000000000000000000000..91d0a2ec11a747ab30147ff199367f23de7442a1
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu
@@ -0,0 +1 @@
+../expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu
\ No newline at end of file
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/pointheatsource_quadratic-mesh.xml b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/pointheatsource_quadratic-mesh.xml
new file mode 100644
index 0000000000000000000000000000000000000000..0d4b1c81a5f3128d3bb28b2054906851399eb569
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/pointheatsource_quadratic-mesh.xml
@@ -0,0 +1,15 @@
+<?xml version='1.0' encoding='ISO-8859-1'?>
+<OpenGeoSysProjectDiff base_file="../pointheatsource_quadratic-mesh.prj">
+    <replace sel="//source_term/type/text()">
+        Python
+    </replace>
+    <add sel="//source_term">
+        <source_term_object>source_term</source_term_object>
+    </add>
+    <add sel="/*">
+        <python_script>source_term.py</python_script>
+    </add>
+
+    <!-- axial symmetry would multiply with an integral measure of 0 at the center point, thereby disabling the heat source -->
+    <remove sel="/*/meshes/mesh[text()='quarter_circle_geometry_center.vtu']/@axially_symmetric" />
+</OpenGeoSysProjectDiff>
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/quarter_002_2nd.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/quarter_002_2nd.vtu
new file mode 120000
index 0000000000000000000000000000000000000000..571461b87faa9e664c0ad79cb6b1688c63919d1a
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/quarter_002_2nd.vtu
@@ -0,0 +1 @@
+../quarter_002_2nd.vtu
\ No newline at end of file
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/quarter_circle_geometry_bottom.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/quarter_circle_geometry_bottom.vtu
new file mode 120000
index 0000000000000000000000000000000000000000..2befae203d94b6ca745205283bbb5017c052047c
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/quarter_circle_geometry_bottom.vtu
@@ -0,0 +1 @@
+../quarter_circle_geometry_bottom.vtu
\ No newline at end of file
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/quarter_circle_geometry_center.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/quarter_circle_geometry_center.vtu
new file mode 120000
index 0000000000000000000000000000000000000000..7e3e61a8e92c7b34bd79ff5d9536b906dccec23b
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/quarter_circle_geometry_center.vtu
@@ -0,0 +1 @@
+../quarter_circle_geometry_center.vtu
\ No newline at end of file
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/quarter_circle_geometry_left.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/quarter_circle_geometry_left.vtu
new file mode 120000
index 0000000000000000000000000000000000000000..1619e5104c1d69dcfb8e92363ff495f4e5ad974b
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/quarter_circle_geometry_left.vtu
@@ -0,0 +1 @@
+../quarter_circle_geometry_left.vtu
\ No newline at end of file
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/quarter_circle_geometry_out.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/quarter_circle_geometry_out.vtu
new file mode 120000
index 0000000000000000000000000000000000000000..09740266319fe37dc42fce3eb4268ca9a554f736
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/quarter_circle_geometry_out.vtu
@@ -0,0 +1 @@
+../quarter_circle_geometry_out.vtu
\ No newline at end of file
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/source_term.py b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/source_term.py
new file mode 100644
index 0000000000000000000000000000000000000000..21296f2c4f17067fe3a89b048d917e21049bb694
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-python-source-term/source_term.py
@@ -0,0 +1,15 @@
+try:
+    import ogs.callbacks as OpenGeoSys
+except ModuleNotFoundError:
+    import OpenGeoSys
+
+
+class SourceTerm(OpenGeoSys.SourceTerm):
+    def getFlux(self, _t, coords, _primary_vars):
+        x, y, z = coords
+        value = 150
+        Jac = [0.0] * 4
+        return (value, Jac)
+
+
+source_term = SourceTerm()
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu
new file mode 120000
index 0000000000000000000000000000000000000000..91d0a2ec11a747ab30147ff199367f23de7442a1
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu
@@ -0,0 +1 @@
+../expected_pointheatsource_quadratic-mesh_ts_10_t_50000.000000.vtu
\ No newline at end of file
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/pointheatsource_quadratic-mesh.xml b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/pointheatsource_quadratic-mesh.xml
new file mode 100644
index 0000000000000000000000000000000000000000..085ce9217c96fcbe5da2863a7f617caf48505563
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/pointheatsource_quadratic-mesh.xml
@@ -0,0 +1,9 @@
+<?xml version='1.0' encoding='ISO-8859-1'?>
+<OpenGeoSysProjectDiff base_file="../pointheatsource_quadratic-mesh.prj">
+    <replace sel="//source_term/type/text()">
+        Volumetric
+    </replace>
+
+    <!-- axial symmetry would multiply with an integral measure of 0 at the center point, thereby disabling the heat source -->
+    <remove sel="/*/meshes/mesh[text()='quarter_circle_geometry_center.vtu']/@axially_symmetric" />
+</OpenGeoSysProjectDiff>
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/quarter_002_2nd.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/quarter_002_2nd.vtu
new file mode 120000
index 0000000000000000000000000000000000000000..571461b87faa9e664c0ad79cb6b1688c63919d1a
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/quarter_002_2nd.vtu
@@ -0,0 +1 @@
+../quarter_002_2nd.vtu
\ No newline at end of file
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/quarter_circle_geometry_bottom.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/quarter_circle_geometry_bottom.vtu
new file mode 120000
index 0000000000000000000000000000000000000000..2befae203d94b6ca745205283bbb5017c052047c
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/quarter_circle_geometry_bottom.vtu
@@ -0,0 +1 @@
+../quarter_circle_geometry_bottom.vtu
\ No newline at end of file
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/quarter_circle_geometry_center.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/quarter_circle_geometry_center.vtu
new file mode 120000
index 0000000000000000000000000000000000000000..7e3e61a8e92c7b34bd79ff5d9536b906dccec23b
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/quarter_circle_geometry_center.vtu
@@ -0,0 +1 @@
+../quarter_circle_geometry_center.vtu
\ No newline at end of file
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/quarter_circle_geometry_left.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/quarter_circle_geometry_left.vtu
new file mode 120000
index 0000000000000000000000000000000000000000..1619e5104c1d69dcfb8e92363ff495f4e5ad974b
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/quarter_circle_geometry_left.vtu
@@ -0,0 +1 @@
+../quarter_circle_geometry_left.vtu
\ No newline at end of file
diff --git a/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/quarter_circle_geometry_out.vtu b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/quarter_circle_geometry_out.vtu
new file mode 120000
index 0000000000000000000000000000000000000000..09740266319fe37dc42fce3eb4268ca9a554f736
--- /dev/null
+++ b/Tests/Data/ThermoHydroMechanics/Linear/Point_injection/with-volumetric-source-term/quarter_circle_geometry_out.vtu
@@ -0,0 +1 @@
+../quarter_circle_geometry_out.vtu
\ No newline at end of file