40 DisplacementDim, ConstitutiveTraits>::
41 ThermoRichardsMechanicsLocalAssembler(
45 bool const is_axially_symmetric,
49 e, integration_method, is_axially_symmetric, process_data)
51 unsigned const n_integration_points =
54 ip_data_.resize(n_integration_points);
56 auto const shape_matrices_u =
59 DisplacementDim>(e, is_axially_symmetric,
62 auto const shape_matrices =
64 DisplacementDim>(e, is_axially_symmetric,
67 for (
unsigned ip = 0; ip < n_integration_points; ip++)
70 auto const& sm_u = shape_matrices_u[ip];
73 sm_u.integralMeasure * sm_u.detJ;
76 ip_data.dNdx_u = sm_u.dNdx;
79 ip_data.N_p = shape_matrices[ip].N;
80 ip_data.dNdx_p = shape_matrices[ip].dNdx;
88 ConstitutiveTraits>::setInitialConditionsConcrete(Eigen::VectorXd
const
93 assert(local_x.size() ==
94 temperature_size + pressure_size + displacement_size);
96 auto const p_L = local_x.template segment<pressure_size>(pressure_index);
98 local_x.template segment<temperature_size>(temperature_index);
100 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
102 *this->process_data_.media_map.getMedium(this->element_.getID());
106 typename ConstitutiveTraits::ConstitutiveSetting
const constitutive_setting;
107 typename ConstitutiveTraits::ConstitutiveModels models(
108 this->process_data_, this->solid_material_);
110 unsigned const n_integration_points =
111 this->integration_method_.getNumberOfPoints();
112 for (
unsigned ip = 0; ip < n_integration_points; ip++)
115 auto const& N = ip_data_[ip].N_p;
118 std::nullopt, this->element_.getID(), ip,
122 this->element_, ip_data_[ip].N_u))};
139 medium.property(MPL::PropertyType::saturation)
140 .template value<double>(variables, x_position, t, dt);
141 std::get<PrevState<SaturationData>>(this->prev_states_[ip])->S_L = S_L;
143 constitutive_setting.init(models, t, dt, x_position, media_data,
144 {T_ip, 0, {}}, this->current_states_[ip],
145 this->prev_states_[ip]);
147 if (this->process_data_.initial_stress.value)
150 convertInitialStressType(ip, t, x_position, medium, variables,
160 ConstitutiveTraits>::
161 convertInitialStressType(
unsigned const ip,
166 double const p_at_ip)
168 bool constexpr is_strain_temperature_constitutive =
171 ConstitutiveTraits>::value;
172 if (is_strain_temperature_constitutive &&
173 this->process_data_.initial_stress.type ==
179 if (!is_strain_temperature_constitutive &&
185 double const alpha_b =
186 medium.
property(MPL::PropertyType::biot_coefficient)
187 .template value<double>(variables, x_position, t, 0.0 );
189 double const bishop =
190 medium.
property(MPL::PropertyType::bishops_effective_stress)
191 .template value<double>(variables, x_position, t, 0.0 );
193 ConstitutiveTraits::ConstitutiveSetting::convertInitialStressType(
194 this->current_states_[ip], this->prev_states_[ip],
195 bishop * alpha_b * p_at_ip * Invariants::identity2);
202 ConstitutiveTraits>::
203 assembleWithJacobian(
double const t,
double const dt,
204 std::vector<double>
const& local_x,
205 std::vector<double>
const& local_x_prev,
206 std::vector<double>& ,
207 std::vector<double>& ,
208 std::vector<double>& local_rhs_data,
209 std::vector<double>& local_Jac_data)
212 *this->process_data_.media_map.getMedium(this->element_.getID());
219 typename ConstitutiveTraits::ConstitutiveSetting constitutive_setting;
221 for (
unsigned ip = 0; ip < this->integration_method_.getNumberOfPoints();
225 std::nullopt, this->element_.getID(), ip,
229 this->element_, ip_data_[ip].N_u))};
231 assembleWithJacobianSingleIP(
233 local_x, local_x_prev,
234 ip_data_[ip], constitutive_setting,
237 this->current_states_[ip], this->prev_states_[ip],
238 this->material_states_[ip], this->output_data_[ip]);
239 loc_mat += loc_mat_current_ip;
242 massLumping(loc_mat);
244 addToLocalMatrixData(dt, local_x, local_x_prev, loc_mat, local_rhs_data,
272 ConstitutiveTraits>::
273 addToLocalMatrixData(
275 std::vector<double>
const& local_x,
276 std::vector<double>
const& local_x_prev,
278 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
279 ConstitutiveTraits>::LocalMatrices
const& loc_mat,
280 std::vector<double>& local_rhs_data,
281 std::vector<double>& local_Jac_data)
const
283 constexpr auto local_matrix_dim =
284 displacement_size + pressure_size + temperature_size;
285 assert(local_x.size() == local_matrix_dim);
288 typename ShapeMatricesTypeDisplacement::template MatrixType<
289 local_matrix_dim, local_matrix_dim>>(
290 local_Jac_data, local_matrix_dim, local_matrix_dim);
294 template VectorType<local_matrix_dim>>(
295 local_rhs_data, local_matrix_dim);
297 local_Jac.noalias() = loc_mat.Jac;
298 local_rhs.noalias() = -loc_mat.res;
303 block_TT(local_Jac).noalias() += loc_mat.M_TT / dt + loc_mat.K_TT;
304 block_Tp(local_Jac).noalias() +=
305 loc_mat.K_Tp + loc_mat.dK_TT_dp + loc_mat.M_Tp / dt;
307 block_pT(local_Jac).noalias() += loc_mat.M_pT / dt + loc_mat.K_pT;
308 block_pp(local_Jac).noalias() +=
309 loc_mat.K_pp + loc_mat.storage_p_a_p / dt + loc_mat.storage_p_a_S_Jpp;
310 block_pu(local_Jac).noalias() = loc_mat.M_pu / dt;
315 auto const [T, p_L, u] = localDOF(local_x);
316 auto const [T_prev, p_L_prev, u_prev] = localDOF(local_x_prev);
318 block_T(local_rhs).noalias() -= loc_mat.M_TT * (T - T_prev) / dt +
319 loc_mat.K_TT * T + loc_mat.K_Tp * p_L +
320 loc_mat.M_Tp * (p_L - p_L_prev) / dt;
321 block_p(local_rhs).noalias() -=
322 loc_mat.K_pp * p_L + loc_mat.K_pT * T +
323 (loc_mat.storage_p_a_p + loc_mat.storage_p_a_S) * (p_L - p_L_prev) /
325 loc_mat.M_pu * (u - u_prev) / dt + loc_mat.M_pT * (T - T_prev) / dt;
332 ConstitutiveTraits>::
333 assembleWithJacobianSingleIP(
334 double const t,
double const dt,
336 std::vector<double>
const& local_x,
337 std::vector<double>
const& local_x_prev,
339 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
340 ConstitutiveTraits>::IpData
const& ip_data,
341 typename ConstitutiveTraits::ConstitutiveSetting& CS,
344 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
345 ConstitutiveTraits>::LocalMatrices& out,
346 typename ConstitutiveTraits::StatefulData& current_state,
347 typename ConstitutiveTraits::StatefulDataPrev
const& prev_state,
349 typename ConstitutiveTraits::OutputData& output_data)
const
351 auto const& N_u = ip_data.N_u;
352 auto const& dNdx_u = ip_data.dNdx_u;
355 auto const& N = ip_data.N_p;
356 auto const& dNdx = ip_data.dNdx_p;
360 ShapeFunctionDisplacement::NPOINTS,
363 this->is_axially_symmetric_);
365 auto const [T, p_L, u] = localDOF(local_x);
366 auto const [T_prev, p_L_prev, u_prev] = localDOF(local_x_prev);
370 typename ConstitutiveTraits::ConstitutiveModels models(
371 this->process_data_, this->solid_material_);
372 typename ConstitutiveTraits::ConstitutiveTempData tmp;
373 typename ConstitutiveTraits::ConstitutiveData CD;
376 double const T_ip = N * T;
377 double const T_prev_ip = N * T_prev;
379 double const p_cap_ip = -N * p_L;
380 double const p_cap_prev_ip = -N * p_L_prev;
385 CS.eval(models, t, dt, x_position,
387 {T_ip, T_prev_ip, grad_T_ip},
388 {p_cap_ip, p_cap_prev_ip, grad_p_cap_ip},
389 eps, current_state, prev_state, mat_state, tmp, output_data,
394 NodalMatrix
const NTN = N.transpose() * N;
395 NodalMatrix
const dNTdN = dNdx.transpose() * dNdx;
400 DisplacementDim)>::identity2;
401 typename ShapeMatricesTypeDisplacement::template MatrixType<
402 displacement_size, pressure_size>
403 BTI2N = B.transpose() * identity2 * N;
425 block_p(out.res).noalias() =
426 dNdx.transpose() * std::get<EqPData<DisplacementDim>>(CD).rhs_p_dNT_V;
427 block_u(out.res).noalias() =
429 ProcessLib::Graph::get<TotalStressData<DisplacementDim>>(
432 static_cast<int>(this->process_data_.apply_body_force_for_deformation) *
433 N_u_op(N_u).transpose() *
434 std::get<GravityData<DisplacementDim>>(CD).volumetric_body_force;
437 out.storage_p_a_p.noalias() =
438 std::get<EqPData<DisplacementDim>>(CD).storage_p_a_p_X_NTN * NTN;
439 out.storage_p_a_S.noalias() =
440 std::get<TRMStorageData>(CD).storage_p_a_S_X_NTN * NTN;
441 out.storage_p_a_S_Jpp.noalias() =
442 std::get<TRMStorageData>(CD).storage_p_a_S_Jpp_X_NTN * NTN;
446 std::get<EqTData<DisplacementDim>>(CD).M_TT_X_NTN * NTN;
448 std::get<TRMVaporDiffusionData<DisplacementDim>>(CD).M_Tp_X_NTN * NTN;
451 std::get<EqPData<DisplacementDim>>(CD).M_pT_X_NTN * NTN;
453 std::get<EqPData<DisplacementDim>>(CD).M_pu_X_BTI2N * BTI2N.transpose();
458 std::get<TRMHeatStorageAndFluxData<DisplacementDim>>(CD)
462 (std::get<EqTData<DisplacementDim>>(CD).K_TT_NT_V_dN.transpose() *
467 out.dK_TT_dp.noalias() =
469 (std::get<TRMHeatStorageAndFluxData<DisplacementDim>>(CD)
470 .K_Tp_NT_V_dN.transpose() *
476 std::get<ThermoOsmosisData<DisplacementDim>>(CD).K_Tp_Laplace *
482 dNdx.transpose() * std::get<EqPData<DisplacementDim>>(CD).K_pp_Laplace *
488 std::get<ThermoOsmosisData<DisplacementDim>>(CD).K_pT_Laplace * dNdx;
491 block_pT(out.Jac).noalias() =
492 std::get<TRMVaporDiffusionData<DisplacementDim>>(CD).J_pT_X_dNTdN *
494 block_pp(out.Jac).noalias() =
495 std::get<TRMStorageData>(CD).J_pp_X_NTN * NTN +
496 std::get<EqPData<DisplacementDim>>(CD).J_pp_X_BTI2NT_u_dot_N *
497 BTI2N.transpose() * (u - u_prev) / dt *
500 std::get<EqPData<DisplacementDim>>(CD).J_pp_dNT_V_N * N;
502 block_uT(out.Jac).noalias() =
504 std::get<SolidMechanicsDataStateless<DisplacementDim>>(CD).J_uT_BT_K_N *
506 block_up(out.Jac).noalias() =
508 std::get<SolidMechanicsDataStateless<DisplacementDim>>(CD)
511 N_u_op(N_u).transpose() *
512 std::get<GravityData<DisplacementDim>>(CD).J_up_HT_V_N * N;
513 block_uu(out.Jac).noalias() =
515 std::get<SolidMechanicsDataStateless<DisplacementDim>>(CD)
519 out *= ip_data.integration_weight;
526 ConstitutiveTraits>::
527 computeSecondaryVariableConcrete(
double const t,
double const dt,
528 Eigen::VectorXd
const& local_x,
529 Eigen::VectorXd
const& local_x_prev)
531 auto const T = block_T(local_x);
532 auto const p_L = block_p(local_x);
533 auto const u = block_u(local_x);
535 auto const T_prev = block_T(local_x_prev);
536 auto const p_L_prev = block_p(local_x_prev);
538 auto const e_id = this->element_.getID();
539 auto const& process_data = this->process_data_;
540 auto& medium = *process_data.media_map.getMedium(e_id);
542 unsigned const n_integration_points =
543 this->integration_method_.getNumberOfPoints();
545 double saturation_avg = 0;
546 double porosity_avg = 0;
547 double liquid_density_avg = 0;
548 double viscosity_avg = 0;
551 KV sigma_avg = KV::Zero();
553 typename ConstitutiveTraits::ConstitutiveSetting constitutive_setting;
555 typename ConstitutiveTraits::ConstitutiveModels models(
556 process_data, this->solid_material_);
557 typename ConstitutiveTraits::ConstitutiveTempData tmp;
558 typename ConstitutiveTraits::ConstitutiveData CD;
560 for (
unsigned ip = 0; ip < n_integration_points; ip++)
562 auto& current_state = this->current_states_[ip];
563 auto& output_data = this->output_data_[ip];
565 auto const& ip_data = ip_data_[ip];
568 auto const& N = ip_data.N_p;
569 auto const& N_u = ip_data.N_u;
570 auto const& dNdx_u = ip_data.dNdx_u;
571 auto const& dNdx = ip_data.dNdx_p;
574 std::nullopt, this->element_.getID(), ip,
578 this->element_, N_u))};
582 this->element_, N_u);
585 ShapeFunctionDisplacement::NPOINTS,
587 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
589 double const T_ip = N * T;
590 double const T_prev_ip = N * T_prev;
593 double const p_cap_ip = -N * p_L;
594 double const p_cap_prev_ip = -N * p_L_prev;
599 constitutive_setting.eval(models,
602 {T_ip, T_prev_ip, grad_T_ip},
603 {p_cap_ip, p_cap_prev_ip, grad_p_cap_ip},
604 eps, current_state, this->prev_states_[ip],
605 this->material_states_[ip], tmp, output_data,
608 saturation_avg += std::get<SaturationData>(current_state).S_L;
609 porosity_avg += std::get<PorosityData>(current_state).phi;
611 liquid_density_avg += std::get<LiquidDensityData>(output_data).rho_LR;
612 viscosity_avg += std::get<LiquidViscosityData>(output_data).viscosity;
613 sigma_avg += ConstitutiveTraits::ConstitutiveSetting::statefulStress(
616 saturation_avg /= n_integration_points;
617 porosity_avg /= n_integration_points;
618 viscosity_avg /= n_integration_points;
619 liquid_density_avg /= n_integration_points;
620 sigma_avg /= n_integration_points;
622 (*process_data.element_saturation)[e_id] = saturation_avg;
623 (*process_data.element_porosity)[e_id] = porosity_avg;
624 (*process_data.element_liquid_density)[e_id] = liquid_density_avg;
625 (*process_data.element_viscosity)[e_id] = viscosity_avg;
628 &(*process_data.element_stresses)[e_id * KV::RowsAtCompileTime]) =
632 ShapeFunction,
typename ShapeFunctionDisplacement::MeshElement,
633 DisplacementDim>(this->element_, this->is_axially_symmetric_, p_L,
634 *process_data.pressure_interpolated);
636 ShapeFunction,
typename ShapeFunctionDisplacement::MeshElement,
637 DisplacementDim>(this->element_, this->is_axially_symmetric_, T,
638 *process_data.temperature_interpolated);