diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index c5766420777908d3492593da3bd4adde5519256e..56c65416ba6c9bf73ec25a80b1d5c50ed0c960a7 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -95,8 +95,9 @@ template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
typename IntegrationMethod, int DisplacementDim>
void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
IntegrationMethod, DisplacementDim>::
- updateConstitutiveVariables(Eigen::VectorXd const& local_x, double const t,
- double const dt)
+ updateConstitutiveVariables(Eigen::VectorXd const& local_x,
+ Eigen::VectorXd const& local_x_dot,
+ double const t, double const dt)
{
[[maybe_unused]] auto const matrix_size =
gas_pressure_size + capillary_pressure_size + temperature_size +
@@ -112,6 +113,8 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
auto const temperature =
local_x.template segment<temperature_size>(temperature_index);
+ auto const temperature_dot =
+ local_x_dot.template segment<temperature_size>(temperature_index);
auto const displacement =
local_x.template segment<displacement_size>(displacement_index);
@@ -142,19 +145,11 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ShapeMatricesTypeDisplacement>(
_element, Nu);
- auto const Bu =
- LinearBMatrix::computeBMatrix<DisplacementDim,
- ShapeFunctionDisplacement::NPOINTS,
- typename BMatricesType::BMatrixType>(
- gradNu, Nu, x_coord, _is_axially_symmetric);
-
double const T = NT.dot(temperature);
double const pGR = Np.dot(gas_pressure);
double const pCap = Np.dot(capillary_pressure);
double const pLR = pGR - pCap;
- double const div_u = mT * Bu * displacement;
-
MPL::VariableArray vars;
vars[static_cast<int>(MPL::Variable::temperature)] = T;
vars[static_cast<int>(MPL::Variable::phase_pressure)] = pGR;
@@ -204,6 +199,28 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
// isotropic solid phase volumetric thermal expansion coefficient
ip_data.beta_T_SR = Invariants::trace(ip_data.alpha_T_SR);
+ double const T_dot = NT.dot(temperature_dot);
+ MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
+ dthermal_strain = ip_data.alpha_T_SR * T_dot * dt;
+
+ auto const Bu =
+ LinearBMatrix::computeBMatrix<DisplacementDim,
+ ShapeFunctionDisplacement::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ gradNu, Nu, x_coord, _is_axially_symmetric);
+
+ auto& eps = ip_data.eps;
+ eps.noalias() = Bu * displacement;
+
+ auto& eps_prev = ip_data.eps_prev;
+ auto& eps_m = ip_data.eps_m;
+ auto& eps_m_prev = ip_data.eps_m_prev;
+
+ eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
+ vars[static_cast<int>(MaterialPropertyLib::Variable::mechanical_strain)]
+ .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
+ eps_m);
+
auto const rho_ref_SR =
solid_phase.property(MPL::PropertyType::density)
.template value<double>(
@@ -222,6 +239,7 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
.template value<double>(vars, pos, t, dt);
// solid phase volume fraction
+ double const div_u = mT * eps;
auto const phi_S_0 = 1. - phi_0;
const double phi_S = phi_S_0 * (1. + ip_data.thermal_volume_strain -
ip_data.alpha_B * div_u);
@@ -232,6 +250,8 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
auto const rhoSR = rho_ref_SR * (1. - ip_data.thermal_volume_strain +
(ip_data.alpha_B - 1.) * div_u);
+ ip_data.updateConstitutiveRelation(vars, t, pos, dt, T - T_dot * dt);
+
// constitutive model object as specified in process creation
auto& ptm = *_process_data.phase_transition_model_;
ptm.computeConstitutiveVariables(&medium, vars, pos, t, dt);
@@ -318,8 +338,8 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
assert(local_x.size() == matrix_size);
updateConstitutiveVariables(
- Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()), t,
- std::numeric_limits<double>::quiet_NaN());
+ Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
+ Eigen::VectorXd::Zero(matrix_size), t, 0);
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
@@ -345,9 +365,6 @@ void TH2MLocalAssembler<
temperature_size + displacement_size;
assert(local_x.size() == matrix_size);
- auto const temperature = Eigen::Map<VectorType<temperature_size> const>(
- local_x.data() + temperature_index, temperature_size);
-
auto const gas_pressure = Eigen::Map<VectorType<gas_pressure_size> const>(
local_x.data() + gas_pressure_index, gas_pressure_size);
@@ -355,17 +372,11 @@ void TH2MLocalAssembler<
Eigen::Map<VectorType<capillary_pressure_size> const>(
local_x.data() + capillary_pressure_index, capillary_pressure_size);
- auto const displacement = Eigen::Map<VectorType<displacement_size> const>(
- local_x.data() + displacement_index, displacement_size);
-
auto const capillary_pressure_dot =
Eigen::Map<VectorType<capillary_pressure_size> const>(
local_x_dot.data() + capillary_pressure_index,
capillary_pressure_size);
- auto const temperature_dot = Eigen::Map<VectorType<temperature_size> const>(
- local_x_dot.data() + temperature_index, temperature_size);
-
// pointer to local_M_data vector
auto local_M =
MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
@@ -450,10 +461,10 @@ void TH2MLocalAssembler<
_integration_method.getNumberOfPoints();
updateConstitutiveVariables(
- Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()), t,
- dt);
-
- MPL::VariableArray vars;
+ Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
+ Eigen::Map<Eigen::VectorXd const>(local_x_dot.data(),
+ local_x_dot.size()),
+ t, dt);
for (unsigned int_point = 0; int_point < n_integration_points; int_point++)
{
@@ -494,9 +505,6 @@ void TH2MLocalAssembler<
auto const BuT = Bu.transpose().eval();
- auto& eps = ip.eps;
- auto const& sigma_eff = ip.sigma_eff;
-
double const pCap = Np.dot(capillary_pressure);
GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
@@ -705,25 +713,7 @@ void TH2MLocalAssembler<
KUpC.noalias() += (BuT * alpha_B * s_L * m * Np) * w;
- eps.noalias() = Bu * displacement;
-
- auto& eps_prev = ip.eps_prev;
- auto& eps_m = ip.eps_m;
- auto& eps_m_prev = ip.eps_m_prev;
-
- double const T = NT.dot(temperature);
- double const T_dot = NT.dot(temperature_dot);
-
- MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
- dthermal_strain = ip.alpha_T_SR * T_dot * dt;
- eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
- vars[static_cast<int>(MaterialPropertyLib::Variable::mechanical_strain)]
- .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
- eps_m);
-
- ip.updateConstitutiveRelation(vars, t, pos, dt, T - T_dot * dt);
-
- fU.noalias() -= (BuT * sigma_eff - Nu_op.transpose() * rho * b) * w;
+ fU.noalias() -= (BuT * ip.sigma_eff - Nu_op.transpose() * rho * b) * w;
if (_process_data.apply_mass_lumping)
{
@@ -966,100 +956,9 @@ template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
typename IntegrationMethod, int DisplacementDim>
void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
IntegrationMethod, DisplacementDim>::
- postNonLinearSolverConcrete(std::vector<double> const& local_x,
- std::vector<double> const& local_xdot,
- double const t, double const dt,
- bool const use_monolithic_scheme,
- int const /*process_id*/)
-{
- const int displacement_offset =
- use_monolithic_scheme ? displacement_index : 0;
-
- auto const u =
- Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
- displacement_size> const>(local_x.data() + displacement_offset,
- displacement_size);
-
- auto const temperature =
- Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
- temperature_size> const>(local_x.data() + temperature_index,
- temperature_size);
- auto const p =
- Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
- gas_pressure_size> const>(local_x.data() + gas_pressure_index,
- gas_pressure_size);
-
- auto const dT =
- Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
- temperature_size> const>(local_xdot.data() + temperature_index,
- temperature_size) *
- dt;
-
- 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 T = N_T.dot(temperature);
- vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] = T;
- vars[static_cast<int>(MaterialPropertyLib::Variable::phase_pressure)] =
- N_T.dot(p); // N_T = N_p
-
- // solid phase linear thermal expansion coefficient
- auto const alpha_T_SR = MathLib::KelvinVector::tensorToKelvin<
- DisplacementDim>(MaterialPropertyLib::formEigenTensor<3>(
- solid_phase
- .property(
- MaterialPropertyLib::PropertyType::thermal_expansivity)
- .value(vars, x_position, t, dt)));
-
- double const dT_int_pt = N_T.dot(dT);
-
- MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
- dthermal_strain = alpha_T_SR * dT_int_pt;
-
- auto& eps = _ip_data[ip].eps;
- eps.noalias() = B * u;
-
- auto& eps_prev = _ip_data[ip].eps_prev;
- auto& eps_m = _ip_data[ip].eps_m;
- auto& eps_m_prev = _ip_data[ip].eps_m_prev;
- eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
- vars[static_cast<int>(MaterialPropertyLib::Variable::mechanical_strain)]
- .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
- eps_m);
-
- _ip_data[ip].updateConstitutiveRelation(vars, t, x_position, dt,
- T - dT_int_pt);
- }
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
- typename IntegrationMethod, int DisplacementDim>
-void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
- IntegrationMethod, DisplacementDim>::
- computeSecondaryVariableConcrete(double const t, double const /*dt*/,
+ computeSecondaryVariableConcrete(double const t, double const dt,
Eigen::VectorXd const& local_x,
- Eigen::VectorXd const& /*local_x_dot*/)
+ Eigen::VectorXd const& local_x_dot)
{
auto const gas_pressure =
local_x.template segment<gas_pressure_size>(gas_pressure_index);
@@ -1096,7 +995,7 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
double saturation_avg = 0;
- updateConstitutiveVariables(local_x, t, 0);
+ updateConstitutiveVariables(local_x, local_x_dot, t, dt);
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h
index 55c0496adc1463d93f175d2ef1c44ed68069d37c..d93b5a2af7be26b2ca4411d1bc1a5188de15c10d 100644
--- a/ProcessLib/TH2M/TH2MFEM.h
+++ b/ProcessLib/TH2M/TH2MFEM.h
@@ -70,8 +70,6 @@ public:
MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
- using SymmetricTensor = Eigen::Matrix<double, KelvinVectorSize, 1>;
-
TH2MLocalAssembler(TH2MLocalAssembler const&) = delete;
TH2MLocalAssembler(TH2MLocalAssembler&&) = delete;
@@ -151,12 +149,6 @@ private:
double const t, double const dt, Eigen::VectorXd const& local_x,
Eigen::VectorXd const& local_x_dot) override;
- void postNonLinearSolverConcrete(std::vector<double> const& local_x,
- std::vector<double> const& local_xdot,
- double const t, double const dt,
- bool const use_monolithic_scheme,
- int const process_id) override;
-
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
@@ -178,6 +170,7 @@ private:
std::vector<double>& cache) const override;
void updateConstitutiveVariables(Eigen::VectorXd const& local_x,
+ Eigen::VectorXd const& local_x_dot,
double const t, double const dt);
std::size_t setSigma(double const* values)
diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp
index f944ec6b0a3df3d4f2cc1956533a9e6feada9270..0b8658495b2a2ef5d094b9f59948a96ee02e6b41 100644
--- a/ProcessLib/TH2M/TH2MProcess.cpp
+++ b/ProcessLib/TH2M/TH2MProcess.cpp
@@ -364,24 +364,6 @@ void TH2MProcess<DisplacementDim>::postTimestepConcreteProcess(
x, t, dt);
}
-template <int DisplacementDim>
-void TH2MProcess<DisplacementDim>::postNonLinearSolverConcreteProcess(
- GlobalVector const& x, GlobalVector const& xdot, const double t,
- double const dt, const int process_id)
-{
- if (!hasMechanicalProcess(process_id))
- {
- return;
- }
-
- DBUG("PostNonLinearSolver TH2MProcess.");
- // Calculate strain, stress or other internal variables of mechanics.
- GlobalExecutor::executeMemberOnDereferenced(
- &LocalAssemblerInterface::postNonLinearSolver, _local_assemblers,
- getDOFTable(process_id), x, xdot, t, dt, _use_monolithic_scheme,
- process_id);
-}
-
template <int DisplacementDim>
void TH2MProcess<DisplacementDim>::computeSecondaryVariableConcrete(
double const t, double const dt, std::vector<GlobalVector*> const& x,
diff --git a/ProcessLib/TH2M/TH2MProcess.h b/ProcessLib/TH2M/TH2MProcess.h
index 1983344d2e5a530ffc3e018bdd44d0b52726abd9..599035de8d6972e6f1366a2d92f1583928e64b7c 100644
--- a/ProcessLib/TH2M/TH2MProcess.h
+++ b/ProcessLib/TH2M/TH2MProcess.h
@@ -95,11 +95,6 @@ private:
const double t, const double dt,
int const /*process_id*/) override;
- void postNonLinearSolverConcreteProcess(GlobalVector const& x,
- GlobalVector const& xdot,
- const double t, double const dt,
- int const process_id) override;
-
NumLib::LocalToGlobalIndexMap const& getDOFTable(
const int process_id) const override;