diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
index 2206073dfd50d206e82f728665ff21efe07b726f..4a5294b938bdf5b3a178149238042f4752919517 100644
--- a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
@@ -230,9 +230,6 @@ void ThermoHydroMechanicsLocalAssembler<
             _process_data.solid_materials, _process_data.material_ids,
             _element.getID());
 
-    ParameterLib::SpatialPosition x_position;
-    x_position.setElementID(_element.getID());
-
     auto const& medium = _process_data.media_map->getMedium(_element.getID());
     auto const& liquid_phase = medium->phase("AqueousLiquid");
     auto const& solid_phase = medium->phase("Solid");
@@ -255,7 +252,6 @@ void ThermoHydroMechanicsLocalAssembler<
         _integration_method.getNumberOfPoints();
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
-        x_position.setIntegrationPoint(ip);
         auto const& w = _ip_data[ip].integration_weight;
 
         auto const& N_u_op = _ip_data[ip].N_u_op;
@@ -273,6 +269,12 @@ void ThermoHydroMechanicsLocalAssembler<
         auto const T_int_pt = N_T.dot(T);
         double const dT_int_pt = N_T.dot(T_dot) * dt;
 
+        ParameterLib::SpatialPosition const x_position{
+            std::nullopt, _element.getID(), ip,
+            MathLib::Point3d(
+                NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
+                                               ShapeMatricesTypeDisplacement>(
+                    _element, N_u))};
         auto const x_coord =
             NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
                                            ShapeMatricesTypeDisplacement>(
@@ -648,9 +650,6 @@ std::vector<double> const& ThermoHydroMechanicsLocalAssembler<
         temperature_size> const>(local_x.data() + temperature_index,
                                  temperature_size);
 
-    ParameterLib::SpatialPosition x_position;
-    x_position.setElementID(_element.getID());
-
     auto const& medium = _process_data.media_map->getMedium(_element.getID());
     auto const& liquid_phase = medium->phase("AqueousLiquid");
     auto const& solid_phase = medium->phase("Solid");
@@ -660,10 +659,14 @@ std::vector<double> const& ThermoHydroMechanicsLocalAssembler<
 
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
-        x_position.setIntegrationPoint(ip);
-
         auto const& N_p = _ip_data[ip].N_p;
 
+        ParameterLib::SpatialPosition const x_position{
+            std::nullopt, _element.getID(), ip,
+            MathLib::Point3d(
+                NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
+                                               ShapeMatricesTypeDisplacement>(
+                    _element, _ip_data[ip].N_u))};
         vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
             N_p.dot(T);  // N_p = N_T
         double const p_int_pt = N_p.dot(p);
@@ -761,8 +764,6 @@ void ThermoHydroMechanicsLocalAssembler<
             temperature_size> const>(local_xdot.data() + temperature_index,
                                      temperature_size);
 
-    ParameterLib::SpatialPosition x_position;
-    x_position.setElementID(_element.getID());
     auto const& medium = _process_data.media_map->getMedium(_element.getID());
     auto const& solid_phase = medium->phase("Solid");
     MaterialPropertyLib::VariableArray vars;
@@ -770,11 +771,16 @@ void ThermoHydroMechanicsLocalAssembler<
     int const n_integration_points = _integration_method.getNumberOfPoints();
     for (int ip = 0; ip < n_integration_points; ip++)
     {
-        x_position.setIntegrationPoint(ip);
         auto const& N_u = _ip_data[ip].N_u;
         auto const& N_T = _ip_data[ip].N_p;
         auto const& dNdx_u = _ip_data[ip].dNdx_u;
 
+        ParameterLib::SpatialPosition const x_position{
+            std::nullopt, _element.getID(), ip,
+            MathLib::Point3d(
+                NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
+                                               ShapeMatricesTypeDisplacement>(
+                    _element, N_u))};
         auto const x_coord =
             NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
                                            ShapeMatricesTypeDisplacement>(
diff --git a/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h
index 66377dbd4c1ec07cb8c3c661a002540974e290b1..73f934e13560120f8b63ed4f177fd22f706b0866 100644
--- a/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h
+++ b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h
@@ -54,12 +54,9 @@ ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::
 
     auto const& medium = *_process_data.media_map->getMedium(_element.getID());
 
-    ParameterLib::SpatialPosition x_position;
-    x_position.setElementID(_element.getID());
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
         auto const& sm = shape_matrices[ip];
-        x_position.setIntegrationPoint(ip);
         _ip_data.emplace_back();
         auto& ip_data = _ip_data[ip];
         _ip_data[ip].integration_weight =
@@ -69,6 +66,11 @@ ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::
         ip_data.N = sm.N;
         ip_data.dNdx = sm.dNdx;
 
+        ParameterLib::SpatialPosition const x_position{
+            std::nullopt, _element.getID(), ip,
+            MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
+                                                            ShapeMatricesType>(
+                _element, sm.N))};
         // Initial porosity. Could be read from integration point data or mesh.
         ip_data.porosity = medium[MPL::porosity].template initialValue<double>(
             x_position,
@@ -121,17 +123,18 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::
     auto const& medium = *_process_data.media_map->getMedium(_element.getID());
     MPL::VariableArray variables;
 
-    ParameterLib::SpatialPosition x_position;
-    x_position.setElementID(_element.getID());
-
     unsigned const n_integration_points =
         _integration_method.getNumberOfPoints();
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
-        x_position.setIntegrationPoint(ip);
-
         auto const& N = _ip_data[ip].N;
 
+        ParameterLib::SpatialPosition const x_position{
+            std::nullopt, _element.getID(), ip,
+            MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
+                                                            ShapeMatricesType>(
+                _element, N))};
+
         double p_cap_ip;
         NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
 
@@ -225,19 +228,21 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::
     MPL::VariableArray variables;
     MPL::VariableArray variables_prev;
 
-    ParameterLib::SpatialPosition x_position;
-    x_position.setElementID(_element.getID());
-
     unsigned const n_integration_points =
         _integration_method.getNumberOfPoints();
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
-        x_position.setIntegrationPoint(ip);
         auto const& w = _ip_data[ip].integration_weight;
 
         auto const& N = _ip_data[ip].N;
         auto const& dNdx = _ip_data[ip].dNdx;
 
+        ParameterLib::SpatialPosition const x_position{
+            std::nullopt, _element.getID(), ip,
+            MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
+                                                            ShapeMatricesType>(
+                _element, N))};
+
         double T_ip;
         NumLib::shapeFunctionInterpolate(T, N, T_ip);
         double T_dot_ip;
@@ -733,19 +738,21 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::assemble(
     MPL::VariableArray variables;
     MPL::VariableArray variables_prev;
 
-    ParameterLib::SpatialPosition x_position;
-    x_position.setElementID(_element.getID());
-
     unsigned const n_integration_points =
         _integration_method.getNumberOfPoints();
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
-        x_position.setIntegrationPoint(ip);
         auto const& w = _ip_data[ip].integration_weight;
 
         auto const& N = _ip_data[ip].N;
         auto const& dNdx = _ip_data[ip].dNdx;
 
+        ParameterLib::SpatialPosition const x_position{
+            std::nullopt, _element.getID(), ip,
+            MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
+                                                            ShapeMatricesType>(
+                _element, N))};
+
         double T_ip;
         NumLib::shapeFunctionInterpolate(T, N, T_ip);
         double T_dot_ip;
@@ -1231,9 +1238,6 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::
     MPL::VariableArray variables;
     MPL::VariableArray variables_prev;
 
-    ParameterLib::SpatialPosition x_position;
-    x_position.setElementID(_element.getID());
-
     unsigned const n_integration_points =
         _integration_method.getNumberOfPoints();
 
@@ -1242,9 +1246,14 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::
 
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
-        x_position.setIntegrationPoint(ip);
         auto const& N = _ip_data[ip].N;
 
+        ParameterLib::SpatialPosition const x_position{
+            std::nullopt, _element.getID(), ip,
+            MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
+                                                            ShapeMatricesType>(
+                _element, N))};
+
         double T_ip;
         NumLib::shapeFunctionInterpolate(T, N, T_ip);
         double T_dot_ip;
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
index fa06f3b21868dac80f6f39d5c381a9c0e83e5f23..329301297ce713504be95a96f132e4ff23e74a24 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
@@ -61,11 +61,8 @@ ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunction,
                                   DisplacementDim>(e, is_axially_symmetric,
                                                    integration_method);
 
-    ParameterLib::SpatialPosition x_position;
-    x_position.setElementID(e.getID());
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
-        x_position.setIntegrationPoint(ip);
         auto& ip_data = ip_data_[ip];
         auto const& sm_u = shape_matrices_u[ip];
         ip_data_[ip].integration_weight =
@@ -117,20 +114,21 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         this->process_data_.media_map->getMedium(this->element_.getID());
     MPL::VariableArray variables;
 
-    ParameterLib::SpatialPosition x_position;
-    x_position.setElementID(this->element_.getID());
-
     auto const& solid_phase = medium->phase("Solid");
 
     unsigned const n_integration_points =
         this->integration_method_.getNumberOfPoints();
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
-        x_position.setIntegrationPoint(ip);
-
         // N is used for both T and p variables.
         auto const& N = ip_data_[ip].N_p;
 
+        ParameterLib::SpatialPosition const x_position{
+            std::nullopt, this->element_.getID(), ip,
+            MathLib::Point3d(
+                NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
+                                               ShapeMatricesTypeDisplacement>(
+                    this->element_, ip_data_[ip].N_u))};
         double p_cap_ip;
         NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
 
@@ -177,9 +175,6 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
     auto& medium =
         *this->process_data_.media_map->getMedium(this->element_.getID());
 
-    ParameterLib::SpatialPosition x_position;
-    x_position.setElementID(this->element_.getID());
-
     LocalMatrices loc_mat;
     loc_mat.setZero();
     LocalMatrices loc_mat_current_ip;
@@ -191,7 +186,12 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
     for (unsigned ip = 0; ip < this->integration_method_.getNumberOfPoints();
          ++ip)
     {
-        x_position.setIntegrationPoint(ip);
+        ParameterLib::SpatialPosition const x_position{
+            std::nullopt, this->element_.getID(), ip,
+            MathLib::Point3d(
+                NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
+                                               ShapeMatricesTypeDisplacement>(
+                    this->element_, ip_data_[ip].N_u))};
 
         assembleWithJacobianSingleIP(
             t, dt, x_position,    //
@@ -463,9 +463,6 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
     auto const& process_data = this->process_data_;
     auto& medium = *process_data.media_map->getMedium(e_id);
 
-    ParameterLib::SpatialPosition x_position;
-    x_position.setElementID(e_id);
-
     unsigned const n_integration_points =
         this->integration_method_.getNumberOfPoints();
 
@@ -490,8 +487,6 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         auto& current_state = this->current_states_[ip];
         auto& output_data = this->output_data_[ip];
 
-        x_position.setIntegrationPoint(ip);
-
         auto const& ip_data = ip_data_[ip];
 
         // N is used for both p and T variables
@@ -500,6 +495,12 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         auto const& dNdx_u = ip_data.dNdx_u;
         auto const& dNdx = ip_data.dNdx_p;
 
+        ParameterLib::SpatialPosition const x_position{
+            std::nullopt, this->element_.getID(), ip,
+            MathLib::Point3d(
+                NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
+                                               ShapeMatricesTypeDisplacement>(
+                    this->element_, N_u))};
         auto const x_coord =
             NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
                                            ShapeMatricesTypeDisplacement>(