text: /**
   * \file
   * \copyright
   * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
   *            Distributed under a Modified BSD License.
   *              See accompanying file LICENSE.txt or
   *              http://www.opengeosys.org/project/license
   *
   */
  
  #pragma once
  
  #include "ThermoHydroMechanicsFEM.h"
  
  #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
  
  #include "MaterialLib/MPL/Utils/GetLiquidThermalExpansivity.h"
  #include "MaterialLib/MPL/Medium.h"
  #include "MaterialLib/MPL/Property.h"
  #include "MaterialLib/MPL/Utils/FormEffectiveThermalConductivity.h"
  #include "MaterialLib/MPL/Utils/FormEigenTensor.h"
  
  #include "MathLib/KelvinVector.h"
  #include "NumLib/Function/Interpolation.h"
  #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
  
  namespace ProcessLib
  {
  namespace ThermoHydroMechanics
  {
  template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
            typename IntegrationMethod, int DisplacementDim>
  ThermoHydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
                                     ShapeFunctionPressure, IntegrationMethod,
                                     DisplacementDim>::
      ThermoHydroMechanicsLocalAssembler(
          MeshLib::Element const& e,
          std::size_t const /*local_matrix_size*/,
          bool const is_axially_symmetric,
          unsigned const integration_order,
          ThermoHydroMechanicsProcessData<DisplacementDim>& process_data)
      : _process_data(process_data),
        _integration_method(integration_order),
        _element(e),
        _is_axially_symmetric(is_axially_symmetric)
  {
      unsigned const n_integration_points =
          _integration_method.getNumberOfPoints();
  
      _ip_data.reserve(n_integration_points);
      _secondary_data.N_u.resize(n_integration_points);
  
      auto const shape_matrices_u =
          NumLib::initShapeMatrices<ShapeFunctionDisplacement,
                                    ShapeMatricesTypeDisplacement,
                                    DisplacementDim>(e, is_axially_symmetric,
                                                     _integration_method);
  
      auto const shape_matrices_p =
          NumLib::initShapeMatrices<ShapeFunctionPressure,
                                    ShapeMatricesTypePressure, DisplacementDim>(
              e, is_axially_symmetric, _integration_method);
  
      auto const& solid_material =
          MaterialLib::Solids::selectSolidConstitutiveRelation(
              _process_data.solid_materials, _process_data.material_ids,
              e.getID());
  
      for (unsigned ip = 0; ip < n_integration_points; ip++)
      {
          _ip_data.emplace_back(solid_material);
          auto& ip_data = _ip_data[ip];
          auto const& sm_u = shape_matrices_u[ip];
          ip_data.integration_weight =
              _integration_method.getWeightedPoint(ip).getWeight() *
              sm_u.integralMeasure * sm_u.detJ;
  
          ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
              DisplacementDim, displacement_size>::Zero(DisplacementDim,
                                                        displacement_size);
          for (int i = 0; i < DisplacementDim; ++i)
              ip_data.N_u_op
                  .template block<1, displacement_size / DisplacementDim>(
                      i, i * displacement_size / DisplacementDim)
                  .noalias() = sm_u.N;
  
          ip_data.N_u = sm_u.N;
          ip_data.dNdx_u = sm_u.dNdx;
  
          ip_data.N_p = shape_matrices_p[ip].N;
          ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
  
          _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
      }
  }
  
  // Assembles the local Jacobian matrix. So far, the linearisation of HT part is
  // not considered as that in HT process.
  template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
            typename IntegrationMethod, int DisplacementDim>
  void ThermoHydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
                                          ShapeFunctionPressure,
                                          IntegrationMethod, DisplacementDim>::
      assembleWithJacobian(double const t, double const dt,
                           std::vector<double> const& local_x,
                           std::vector<double> const& local_xdot,
                           const double /*dxdot_dx*/, const double /*dx_dx*/,
                           std::vector<double>& /*local_M_data*/,
                           std::vector<double>& /*local_K_data*/,
                           std::vector<double>& local_rhs_data,
                           std::vector<double>& local_Jac_data)
  {
      assert(local_x.size() ==
             pressure_size + displacement_size + temperature_size);
  
      auto T = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
          temperature_size> const>(local_x.data() + temperature_index,
                                   temperature_size);
  
      auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
          pressure_size> const>(local_x.data() + pressure_index, pressure_size);
  
      auto u =
          Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
              displacement_size> const>(local_x.data() + displacement_index,
                                        displacement_size);
  
      auto T_dot =
          Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
              temperature_size> const>(local_xdot.data() + temperature_index,
                                       temperature_size);
  
      auto p_dot =
          Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
              pressure_size> const>(local_xdot.data() + pressure_index,
                                    pressure_size);
      auto u_dot =
          Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
              displacement_size> const>(local_xdot.data() + displacement_index,
                                        displacement_size);
  
      auto local_Jac = MathLib::createZeroedMatrix<
          typename ShapeMatricesTypeDisplacement::template MatrixType<
              temperature_size + displacement_size + pressure_size,
              temperature_size + displacement_size + pressure_size>>(
          local_Jac_data, displacement_size + pressure_size + temperature_size,
          displacement_size + pressure_size + temperature_size);
  
      auto local_rhs = MathLib::createZeroedVector<
          typename ShapeMatricesTypeDisplacement::template VectorType<
              displacement_size + pressure_size + temperature_size>>(
          local_rhs_data, displacement_size + pressure_size + temperature_size);
  
      typename ShapeMatricesTypePressure::NodalMatrixType MTT;
      MTT.setZero(temperature_size, temperature_size);
  
      typename ShapeMatricesTypePressure::NodalMatrixType KTT;
      KTT.setZero(temperature_size, temperature_size);
  
      typename ShapeMatricesTypePressure::NodalMatrixType KTp;
      KTp.setZero(temperature_size, pressure_size);
  
      typename ShapeMatricesTypePressure::NodalMatrixType laplace_p;
      laplace_p.setZero(pressure_size, pressure_size);
  
      typename ShapeMatricesTypePressure::NodalMatrixType storage_p;
      storage_p.setZero(pressure_size, pressure_size);
  
      typename ShapeMatricesTypePressure::NodalMatrixType storage_T;
      storage_T.setZero(pressure_size, temperature_size);
  
      typename ShapeMatricesTypeDisplacement::template MatrixType<
          displacement_size, pressure_size>
          Kup;
      Kup.setZero(displacement_size, pressure_size);
      typename ShapeMatricesTypeDisplacement::template MatrixType<
          displacement_size, temperature_size>
          KuT;
      KuT.setZero(displacement_size, 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");
      MaterialPropertyLib::VariableArray vars;
  
      auto const& identity2 = MathLib::KelvinVector::Invariants<
          MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value>::
          identity2;
  
      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_u_op = _ip_data[ip].N_u_op;
  
          auto const& N_u = _ip_data[ip].N_u;
          auto const& dNdx_u = _ip_data[ip].dNdx_u;
  
          auto const& N_p = _ip_data[ip].N_p;
          auto const& dNdx_p = _ip_data[ip].dNdx_p;
  
          // same shape function for pressure and temperature since they have the
          // same order
          auto const& N_T = N_p;
          auto const& dNdx_T = dNdx_p;
          auto const T_int_pt = N_T.dot(T);
  
          auto const x_coord =
              NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
                                             ShapeMatricesTypeDisplacement>(
                  _element, N_u);
          auto const B =
              LinearBMatrix::computeBMatrix<DisplacementDim,
                                            ShapeFunctionDisplacement::NPOINTS,
                                            typename BMatricesType::BMatrixType>(
                  dNdx_u, N_u, x_coord, _is_axially_symmetric);
  
          auto& eps = _ip_data[ip].eps;
          auto const& sigma_eff = _ip_data[ip].sigma_eff;
  
          vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
              T_int_pt;
          double const p_int_pt = N_p.dot(p);
          vars[static_cast<int>(MaterialPropertyLib::Variable::phase_pressure)] =
              p_int_pt;
  
          auto const solid_density =
              solid_phase.property(MaterialPropertyLib::PropertyType::density)
                  .template value<double>(vars, x_position, t, dt);
          auto const specific_storage =
              solid_phase.property(MaterialPropertyLib::PropertyType::storage)
                  .template value<double>(vars, x_position, t, dt);
          // Consider anisotropic thermal expansion.
          // Read in 3x3 tensor. 2D case also requires expansion coeff. for z-
          // component.
          Eigen::Matrix<double, 3,
                        3> const solid_linear_thermal_expansion_coefficient =
              MaterialPropertyLib::formEigenTensor<3>(
                  solid_phase
                      .property(
                          MaterialPropertyLib::PropertyType::thermal_expansivity)
                      .value(vars, x_position, t, dt));
  
          auto const porosity =
              solid_phase.property(MaterialPropertyLib::PropertyType::porosity)
                  .template value<double>(vars, x_position, t, dt);
  
          auto const alpha =
              solid_phase
                  .property(MaterialPropertyLib::PropertyType::biot_coefficient)
                  .template value<double>(vars, x_position, t, dt);
          // For stress dependent permeability.
          vars[static_cast<int>(MaterialPropertyLib::Variable::stress)]
              .emplace<SymmetricTensor>(
                  MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
                      (_ip_data[ip].sigma_eff - alpha * identity2 * p_int_pt)
                          .eval()));
  
          auto const intrinsic_permeability =
              MaterialPropertyLib::formEigenTensor<DisplacementDim>(
                  solid_phase
                      .property(MaterialPropertyLib::PropertyType::permeability)
                      .value(vars, x_position, t, dt));
  
          auto const fluid_density =
              liquid_phase.property(MaterialPropertyLib::PropertyType::density)
                  .template value<double>(vars, x_position, t, dt);
  
          double const fluid_volumetric_thermal_expansion_coefficient =
              MaterialPropertyLib::getLiquidThermalExpansivity(
                  liquid_phase, vars, fluid_density, x_position, t, dt);
  
          // Use the viscosity model to compute the viscosity
          auto const viscosity =
              liquid_phase.property(MaterialPropertyLib::PropertyType::viscosity)
                  .template value<double>(vars, x_position, t, dt);
          GlobalDimMatrixType K_over_mu = intrinsic_permeability / viscosity;
  
          double const T0 = _process_data.reference_temperature(t, x_position)[0];
  
          auto const& b = _process_data.specific_body_force;
  
          // TODO (Wenqing) : Change dT to time step wise increment
          double const delta_T(T_int_pt - T0);
          MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
              thermal_strain =
                  MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
                      solid_linear_thermal_expansion_coefficient) *
                  delta_T;
  
          double const rho_s = solid_density * (1 -
                  solid_linear_thermal_expansion_coefficient.trace() * delta_T);
  
          auto velocity = (-K_over_mu * dNdx_p * p).eval();
          velocity += K_over_mu * fluid_density * b;
  
          //
          // displacement equation, displacement part
          //
          eps.noalias() = B * u;
          auto C = _ip_data[ip].updateConstitutiveRelationThermal(
              t, x_position, dt, u,
              _process_data.reference_temperature(t, x_position)[0],
              thermal_strain);
  
          local_Jac
              .template block<displacement_size, displacement_size>(
                  displacement_index, displacement_index)
              .noalias() += B.transpose() * C * B * w;
  
          auto const rho = rho_s * (1 - porosity) + porosity * fluid_density;
          local_rhs.template segment<displacement_size>(displacement_index)
              .noalias() -=
              (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w;
  
          //
          // displacement equation, pressure part (K_up)
          //
  
          Kup.noalias() += B.transpose() * alpha * identity2 * N_p * w;
  
          //
          // pressure equation, pressure part (K_pp and M_pp).
          //
          laplace_p.noalias() += dNdx_p.transpose() * K_over_mu * dNdx_p * w;
  
          storage_p.noalias() += N_p.transpose() * specific_storage * N_p * w;
          //
          //  RHS, pressure part
          //
          local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
              dNdx_p.transpose() * fluid_density * K_over_mu * b * w;
          //
          // pressure equation, temperature part (M_pT)
          //
          auto const beta =
              porosity * fluid_volumetric_thermal_expansion_coefficient +
              (1 - porosity) * solid_linear_thermal_expansion_coefficient.trace();
          storage_T.noalias() += N_T.transpose() * beta * N_T * w;
  
          //
          // pressure equation, displacement part.
          //
          // Reusing Kup.transpose().
  
          //
          // temperature equation, temperature part.
          //
          const double c_f =
              liquid_phase
                  .property(
                      MaterialPropertyLib::PropertyType::specific_heat_capacity)
                  .template value<double>(vars, x_position, t, dt);
          auto const fluid_thermal_conductivity =
              liquid_phase
                  .property(
                      MaterialPropertyLib::PropertyType::thermal_conductivity)
                  .template value<double>(vars, x_position, t, dt);
          GlobalDimMatrixType effective_thermal_condictivity =
              MaterialPropertyLib::formEffectiveThermalConductivity<
                  DisplacementDim>(
                  solid_phase
                      .property(
                          MaterialPropertyLib::PropertyType::thermal_conductivity)
                      .value(vars, x_position, t, dt),
                  fluid_thermal_conductivity, porosity);
  
          KTT.noalias() +=
              (dNdx_T.transpose() * effective_thermal_condictivity * dNdx_T +
               dNdx_T.transpose() * velocity * N_p * fluid_density * c_f) *
              w;
  
          auto const effective_volumetric_heat_capacity =
              porosity * fluid_density * c_f +
              (1.0 - porosity) * solid_density *
                  solid_phase
                      .property(MaterialPropertyLib::PropertyType::
                                    specific_heat_capacity)
                      .template value<double>(vars, x_position, t, dt);
  
          MTT.noalias() +=
              N_T.transpose() * effective_volumetric_heat_capacity * N_T * w;
  
          //
          // temperature equation, pressure part
          //
          KTp.noalias() += fluid_density * c_f * N_T.transpose() *
                           (dNdx_T * T).transpose() * K_over_mu * dNdx_p * w;
      }
      // temperature equation, temperature part
      local_Jac
          .template block<temperature_size, temperature_size>(temperature_index,
                                                              temperature_index)
          .noalias() += KTT + MTT / dt;
  
      // temperature equation, pressure part
      local_Jac
          .template block<temperature_size, pressure_size>(temperature_index,
                                                           pressure_index)
          .noalias() -= KTp;
      // displacement equation, temperature part
      local_Jac
          .template block<displacement_size, temperature_size>(displacement_index,
                                                               temperature_index)
          .noalias() -= KuT;
  
      // displacement equation, pressure part
      local_Jac
          .template block<displacement_size, pressure_size>(displacement_index,
                                                            pressure_index)
          .noalias() -= Kup;
  
      // pressure equation, temperature part.
      local_Jac
          .template block<pressure_size, temperature_size>(pressure_index,
                                                           temperature_index)
          .noalias() -= storage_T / dt;
  
      // pressure equation, pressure part.
      local_Jac
          .template block<pressure_size, pressure_size>(pressure_index,
                                                        pressure_index)
          .noalias() += laplace_p + storage_p / dt;
  
      // pressure equation, displacement part.
      local_Jac
          .template block<pressure_size, displacement_size>(pressure_index,
                                                            displacement_index)
          .noalias() += Kup.transpose() / dt;
  
      // pressure equation (f_p)
      local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
          laplace_p * p + storage_p * p_dot - storage_T * T_dot +
          Kup.transpose() * u_dot;
  
      // displacement equation (f_u)
      local_rhs.template segment<displacement_size>(displacement_index)
          .noalias() += Kup * p;
  
      // temperature equation (f_T)
      local_rhs.template segment<temperature_size>(temperature_index).noalias() -=
          KTT * T + MTT * T_dot;
  }
  
  template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
            typename IntegrationMethod, int DisplacementDim>
  std::vector<double> const& ThermoHydroMechanicsLocalAssembler<
      ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
      DisplacementDim>::
      getIntPtDarcyVelocity(
          const double t,
          std::vector<GlobalVector*> const& x,
          std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
          std::vector<double>& cache) const
  {
      unsigned const n_integration_points =
          _integration_method.getNumberOfPoints();
  
      constexpr int process_id = 0;  // monolithic scheme;
      auto const indices =
          NumLib::getIndices(_element.getID(), *dof_table[process_id]);
      assert(!indices.empty());
      auto const local_x = x[process_id]->get(indices);
  
      cache.clear();
      auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
          double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
          cache, DisplacementDim, n_integration_points);
  
      auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
          pressure_size> const>(local_x.data() + pressure_index, pressure_size);
      auto T = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
          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");
      MaterialPropertyLib::VariableArray vars;
  
      auto const& identity2 = MathLib::KelvinVector::Invariants<
          MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value>::
          identity2;
  
      for (unsigned ip = 0; ip < n_integration_points; ip++)
      {
          x_position.setIntegrationPoint(ip);
  
          auto const& N_p = _ip_data[ip].N_p;
  
          vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
              N_p.dot(T);  // N_p = N_T
          double const p_int_pt = N_p.dot(p);
          vars[static_cast<int>(MaterialPropertyLib::Variable::phase_pressure)] =
              p_int_pt;
  
          // TODO (naumov) Temporary value not used by current material models.
          // Need extension of secondary variables interface.
          double const dt = std::numeric_limits<double>::quiet_NaN();
          auto const viscosity =
              liquid_phase.property(MaterialPropertyLib::PropertyType::viscosity)
                  .template value<double>(vars, x_position, t, dt);
  
          auto const alpha =
              solid_phase
                  .property(MaterialPropertyLib::PropertyType::biot_coefficient)
                  .template value<double>(vars, x_position, t, dt);
          // For stress dependent permeability.
          vars[static_cast<int>(MaterialPropertyLib::Variable::stress)]
              .emplace<SymmetricTensor>(
                  MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
                      (_ip_data[ip].sigma_eff - alpha * identity2 * p_int_pt)
                          .eval()));
  
          GlobalDimMatrixType K_over_mu =
              MaterialPropertyLib::formEigenTensor<DisplacementDim>(
                  solid_phase
                      .property(MaterialPropertyLib::PropertyType::permeability)
                      .value(vars, x_position, t, dt)) /
              viscosity;
  
          auto const fluid_density =
              liquid_phase.property(MaterialPropertyLib::PropertyType::density)
                  .template value<double>(vars, x_position, t, dt);
          auto const& b = _process_data.specific_body_force;
  
          // Compute the velocity
          auto const& dNdx_p = _ip_data[ip].dNdx_p;
          cache_matrix.col(ip).noalias() =
              -K_over_mu * dNdx_p * p + K_over_mu * fluid_density * b;
      }
  
      return cache;
  }
  
  template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
            typename IntegrationMethod, int DisplacementDim>
  void ThermoHydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
                                          ShapeFunctionPressure,
                                          IntegrationMethod, DisplacementDim>::
      postNonLinearSolverConcrete(std::vector<double> const& local_x,
                                  double const t, double const dt,
                                  bool const use_monolithic_scheme)
  {
      const int displacement_offset =
          use_monolithic_scheme ? displacement_index : 0;
  
      auto u =
          Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
              displacement_size> const>(local_x.data() + displacement_offset,
                                        displacement_size);
  
      auto T = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
          temperature_size> const>(local_x.data() + temperature_index,
                                   temperature_size);
      auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
          pressure_size> const>(local_x.data() + pressure_index, pressure_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;
  
      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;
  
          auto const x_coord =
              NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
                                             ShapeMatricesTypeDisplacement>(
                  _element, N_u);
          auto const B =
              LinearBMatrix::computeBMatrix<DisplacementDim,
                                            ShapeFunctionDisplacement::NPOINTS,
                                            typename BMatricesType::BMatrixType>(
                  dNdx_u, N_u, x_coord, _is_axially_symmetric);
  
          double const T0 = _process_data.reference_temperature(t, x_position)[0];
  
          double const T_int_pt = N_T.dot(T);
          vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
              T_int_pt;
          vars[static_cast<int>(MaterialPropertyLib::Variable::phase_pressure)] =
              N_T.dot(p);  // N_T = N_p
  
          // Read in 3x3 tensor. 2D case also requires expansion coeff. for z-
          // component.
          Eigen::Matrix<double, 3,
                        3> const solid_linear_thermal_expansion_coefficient =
              MaterialPropertyLib::formEigenTensor<3>(
                  solid_phase
                      .property(
                          MaterialPropertyLib::PropertyType::thermal_expansivity)
                      .value(vars, x_position, t, dt));
  
          double const delta_T(T_int_pt - T0);
          MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
              thermal_strain =
                  MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
                      solid_linear_thermal_expansion_coefficient) *
                  delta_T;
  
          auto& eps = _ip_data[ip].eps;
          eps.noalias() = B * u;
  
          _ip_data[ip].updateConstitutiveRelationThermal(
              t, x_position, dt, u,
              _process_data.reference_temperature(t, x_position)[0],
              thermal_strain);
      }
  }
  
  template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
            typename IntegrationMethod, int DisplacementDim>
  void ThermoHydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
                                          ShapeFunctionPressure,
                                          IntegrationMethod, DisplacementDim>::
      computeSecondaryVariableConcrete(double const /*t*/, double const /*dt*/,
                                       std::vector<double> const& local_x,
                                       std::vector<double> const& /*local_x_dot*/)
  {
      auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
          pressure_size> const>(local_x.data() + pressure_index, pressure_size);
  
      NumLib::interpolateToHigherOrderNodes<
          ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
          DisplacementDim>(_element, _is_axially_symmetric, p,
                           *_process_data.pressure_interpolated);
  
      auto T = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
          pressure_size> const>(local_x.data() + temperature_index,
                                temperature_size);
  
      NumLib::interpolateToHigherOrderNodes<
          ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
          DisplacementDim>(_element, _is_axially_symmetric, T,
                           *_process_data.temperature_interpolated);
  }
  
  }  // namespace ThermoHydroMechanics
  }  // namespace ProcessLib
  