39 DisplacementDim, ConstitutiveTraits>::
40 ThermoRichardsMechanicsLocalAssembler(
44 bool const is_axially_symmetric,
48 e, integration_method, is_axially_symmetric, process_data)
50 unsigned const n_integration_points =
53 ip_data_.resize(n_integration_points);
55 auto const shape_matrices_u =
58 DisplacementDim>(e, is_axially_symmetric,
61 auto const shape_matrices =
63 DisplacementDim>(e, is_axially_symmetric,
66 for (
unsigned ip = 0; ip < n_integration_points; ip++)
69 auto const& sm_u = shape_matrices_u[ip];
72 sm_u.integralMeasure * sm_u.detJ;
75 ip_data.dNdx_u = sm_u.dNdx;
78 ip_data.N_p = shape_matrices[ip].N;
79 ip_data.dNdx_p = shape_matrices[ip].dNdx;
87 ConstitutiveTraits>::setInitialConditionsConcrete(Eigen::VectorXd
const
92 assert(local_x.size() ==
93 temperature_size + pressure_size + displacement_size);
95 auto const p_L = local_x.template segment<pressure_size>(pressure_index);
97 local_x.template segment<temperature_size>(temperature_index);
99 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
101 *this->process_data_.media_map.getMedium(this->element_.getID());
105 typename ConstitutiveTraits::ConstitutiveSetting
const constitutive_setting;
106 auto models = ConstitutiveTraits::createConstitutiveModels(
107 this->process_data_, this->solid_material_);
109 unsigned const n_integration_points =
110 this->integration_method_.getNumberOfPoints();
111 for (
unsigned ip = 0; ip < n_integration_points; ip++)
114 auto const& N = ip_data_[ip].N_p;
117 std::nullopt, this->element_.getID(), ip,
121 this->element_, ip_data_[ip].N_u))};
138 medium.property(MPL::PropertyType::saturation)
139 .template value<double>(variables, x_position, t, dt);
140 std::get<PrevState<SaturationData>>(this->prev_states_[ip])->S_L = S_L;
142 constitutive_setting.init(models, t, dt, x_position, media_data,
143 {T_ip, 0, {}}, this->current_states_[ip],
144 this->prev_states_[ip]);
146 if (this->process_data_.initial_stress.value)
149 convertInitialStressType(ip, t, x_position, medium, variables,
159 ConstitutiveTraits>::
160 convertInitialStressType(
unsigned const ip,
165 double const p_at_ip)
167 bool constexpr is_strain_temperature_constitutive =
170 ConstitutiveTraits>::value;
171 if (is_strain_temperature_constitutive &&
172 this->process_data_.initial_stress.type ==
178 if (!is_strain_temperature_constitutive &&
184 double const alpha_b =
185 medium.
property(MPL::PropertyType::biot_coefficient)
186 .template value<double>(variables, x_position, t, 0.0 );
188 double const bishop =
189 medium.
property(MPL::PropertyType::bishops_effective_stress)
190 .template value<double>(variables, x_position, t, 0.0 );
192 ConstitutiveTraits::ConstitutiveSetting::convertInitialStressType(
193 this->current_states_[ip], this->prev_states_[ip],
194 bishop * alpha_b * p_at_ip * Invariants::identity2);
201 ConstitutiveTraits>::
202 assembleWithJacobian(
double const t,
double const dt,
203 std::vector<double>
const& local_x,
204 std::vector<double>
const& local_x_prev,
205 std::vector<double>& local_rhs_data,
206 std::vector<double>& local_Jac_data)
209 *this->process_data_.media_map.getMedium(this->element_.getID());
216 typename ConstitutiveTraits::ConstitutiveSetting constitutive_setting;
218 for (
unsigned ip = 0; ip < this->integration_method_.getNumberOfPoints();
222 std::nullopt, this->element_.getID(), ip,
226 this->element_, ip_data_[ip].N_u))};
228 assembleWithJacobianSingleIP(
230 local_x, local_x_prev,
231 ip_data_[ip], constitutive_setting,
234 this->current_states_[ip], this->prev_states_[ip],
235 this->material_states_[ip], this->output_data_[ip]);
236 loc_mat += loc_mat_current_ip;
239 massLumping(loc_mat);
241 addToLocalMatrixData(dt, local_x, local_x_prev, loc_mat, local_rhs_data,
269 ConstitutiveTraits>::
270 addToLocalMatrixData(
272 std::vector<double>
const& local_x,
273 std::vector<double>
const& local_x_prev,
275 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
276 ConstitutiveTraits>::LocalMatrices
const& loc_mat,
277 std::vector<double>& local_rhs_data,
278 std::vector<double>& local_Jac_data)
const
280 constexpr auto local_matrix_dim =
281 displacement_size + pressure_size + temperature_size;
282 assert(local_x.size() == local_matrix_dim);
285 typename ShapeMatricesTypeDisplacement::template MatrixType<
286 local_matrix_dim, local_matrix_dim>>(
287 local_Jac_data, local_matrix_dim, local_matrix_dim);
291 template VectorType<local_matrix_dim>>(
292 local_rhs_data, local_matrix_dim);
294 local_Jac.noalias() = loc_mat.Jac;
295 local_rhs.noalias() = -loc_mat.res;
300 block_TT(local_Jac).noalias() += loc_mat.M_TT / dt + loc_mat.K_TT;
301 block_Tp(local_Jac).noalias() +=
302 loc_mat.K_Tp + loc_mat.dK_TT_dp + loc_mat.M_Tp / dt;
304 block_pT(local_Jac).noalias() += loc_mat.M_pT / dt + loc_mat.K_pT;
305 block_pp(local_Jac).noalias() +=
306 loc_mat.K_pp + loc_mat.storage_p_a_p / dt + loc_mat.storage_p_a_S_Jpp;
307 block_pu(local_Jac).noalias() = loc_mat.M_pu / dt;
312 auto const [T, p_L, u] = localDOF(local_x);
313 auto const [T_prev, p_L_prev, u_prev] = localDOF(local_x_prev);
315 block_T(local_rhs).noalias() -= loc_mat.M_TT * (T - T_prev) / dt +
316 loc_mat.K_TT * T + loc_mat.K_Tp * p_L +
317 loc_mat.M_Tp * (p_L - p_L_prev) / dt;
318 block_p(local_rhs).noalias() -=
319 loc_mat.K_pp * p_L + loc_mat.K_pT * T +
320 (loc_mat.storage_p_a_p + loc_mat.storage_p_a_S) * (p_L - p_L_prev) /
322 loc_mat.M_pu * (u - u_prev) / dt + loc_mat.M_pT * (T - T_prev) / dt;
329 ConstitutiveTraits>::
330 assembleWithJacobianSingleIP(
331 double const t,
double const dt,
333 std::vector<double>
const& local_x,
334 std::vector<double>
const& local_x_prev,
336 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
337 ConstitutiveTraits>::IpData
const& ip_data,
338 typename ConstitutiveTraits::ConstitutiveSetting& CS,
341 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
342 ConstitutiveTraits>::LocalMatrices& out,
343 typename ConstitutiveTraits::StatefulData& current_state,
344 typename ConstitutiveTraits::StatefulDataPrev
const& prev_state,
346 typename ConstitutiveTraits::OutputData& output_data)
const
348 auto const& N_u = ip_data.N_u;
349 auto const& dNdx_u = ip_data.dNdx_u;
352 auto const& N = ip_data.N_p;
353 auto const& dNdx = ip_data.dNdx_p;
357 ShapeFunctionDisplacement::NPOINTS,
360 this->is_axially_symmetric_);
362 typename ConstitutiveTraits::ConstitutiveData CD;
364 auto const [T, p_L, u] = localDOF(local_x);
365 auto const [T_prev, p_L_prev, u_prev] = localDOF(local_x_prev);
368 auto models = ConstitutiveTraits::createConstitutiveModels(
369 this->process_data_, this->solid_material_);
370 typename ConstitutiveTraits::ConstitutiveTempData tmp;
372 double const T_ip = N * T;
373 double const T_prev_ip = N * T_prev;
376 double const p_cap_ip = -N * p_L;
377 double const p_cap_prev_ip = -N * p_L_prev;
382 CS.eval(models, t, dt, x_position,
384 {T_ip, T_prev_ip, grad_T_ip},
385 {p_cap_ip, p_cap_prev_ip, grad_p_cap_ip},
386 eps, current_state, prev_state, mat_state, tmp, output_data,
391 NodalMatrix
const NTN = N.transpose() * N;
392 NodalMatrix
const dNTdN = dNdx.transpose() * dNdx;
397 DisplacementDim)>::identity2;
398 typename ShapeMatricesTypeDisplacement::template MatrixType<
399 displacement_size, pressure_size>
400 BTI2N = B.transpose() * identity2 * N;
422 block_p(out.res).noalias() =
423 dNdx.transpose() * std::get<EqPData<DisplacementDim>>(CD).rhs_p_dNT_V;
424 block_u(out.res).noalias() =
426 ProcessLib::Graph::get<TotalStressData<DisplacementDim>>(
429 static_cast<int>(this->process_data_.apply_body_force_for_deformation) *
430 N_u_op(N_u).transpose() *
431 std::get<GravityData<DisplacementDim>>(CD).volumetric_body_force;
434 out.storage_p_a_p.noalias() =
435 std::get<EqPData<DisplacementDim>>(CD).storage_p_a_p_X_NTN * NTN;
436 out.storage_p_a_S.noalias() =
437 std::get<TRMStorageData>(CD).storage_p_a_S_X_NTN * NTN;
438 out.storage_p_a_S_Jpp.noalias() =
439 std::get<TRMStorageData>(CD).storage_p_a_S_Jpp_X_NTN * NTN;
443 std::get<EqTData<DisplacementDim>>(CD).M_TT_X_NTN * NTN;
445 std::get<TRMVaporDiffusionData<DisplacementDim>>(CD).M_Tp_X_NTN * NTN;
448 std::get<EqPData<DisplacementDim>>(CD).M_pT_X_NTN * NTN;
450 std::get<EqPData<DisplacementDim>>(CD).M_pu_X_BTI2N * BTI2N.transpose();
455 std::get<TRMHeatStorageAndFluxData<DisplacementDim>>(CD)
459 (std::get<EqTData<DisplacementDim>>(CD).K_TT_NT_V_dN.transpose() *
464 out.dK_TT_dp.noalias() =
466 (std::get<TRMHeatStorageAndFluxData<DisplacementDim>>(CD)
467 .K_Tp_NT_V_dN.transpose() *
473 std::get<ThermoOsmosisData<DisplacementDim>>(CD).K_Tp_Laplace *
479 dNdx.transpose() * std::get<EqPData<DisplacementDim>>(CD).K_pp_Laplace *
485 std::get<ThermoOsmosisData<DisplacementDim>>(CD).K_pT_Laplace * dNdx;
488 block_pT(out.Jac).noalias() =
489 std::get<TRMVaporDiffusionData<DisplacementDim>>(CD).J_pT_X_dNTdN *
491 block_pp(out.Jac).noalias() =
492 std::get<TRMStorageData>(CD).J_pp_X_NTN * NTN +
493 std::get<EqPData<DisplacementDim>>(CD).J_pp_X_BTI2NT_u_dot_N *
494 BTI2N.transpose() * (u - u_prev) / dt *
497 std::get<EqPData<DisplacementDim>>(CD).J_pp_dNT_V_N * N;
499 block_uT(out.Jac).noalias() =
501 std::get<SolidMechanicsDataStateless<DisplacementDim>>(CD).J_uT_BT_K_N *
503 block_up(out.Jac).noalias() =
505 std::get<SolidMechanicsDataStateless<DisplacementDim>>(CD)
508 N_u_op(N_u).transpose() *
509 std::get<GravityData<DisplacementDim>>(CD).J_up_HT_V_N * N;
510 block_uu(out.Jac).noalias() =
512 std::get<SolidMechanicsDataStateless<DisplacementDim>>(CD)
516 out *= ip_data.integration_weight;
523 ConstitutiveTraits>::
524 computeSecondaryVariableConcrete(
double const t,
double const dt,
525 Eigen::VectorXd
const& local_x,
526 Eigen::VectorXd
const& local_x_prev)
528 auto const T = block_T(local_x);
529 auto const p_L = block_p(local_x);
530 auto const u = block_u(local_x);
532 auto const T_prev = block_T(local_x_prev);
533 auto const p_L_prev = block_p(local_x_prev);
535 auto const e_id = this->element_.getID();
536 auto const& process_data = this->process_data_;
537 auto& medium = *process_data.media_map.getMedium(e_id);
539 unsigned const n_integration_points =
540 this->integration_method_.getNumberOfPoints();
542 typename ConstitutiveTraits::ConstitutiveSetting constitutive_setting;
544 auto models = ConstitutiveTraits::createConstitutiveModels(
545 process_data, this->solid_material_);
546 typename ConstitutiveTraits::ConstitutiveTempData tmp;
547 typename ConstitutiveTraits::ConstitutiveData CD;
549 for (
unsigned ip = 0; ip < n_integration_points; ip++)
551 auto& current_state = this->current_states_[ip];
552 auto& output_data = this->output_data_[ip];
554 auto const& ip_data = ip_data_[ip];
557 auto const& N = ip_data.N_p;
558 auto const& N_u = ip_data.N_u;
559 auto const& dNdx_u = ip_data.dNdx_u;
560 auto const& dNdx = ip_data.dNdx_p;
563 std::nullopt, this->element_.getID(), ip,
567 this->element_, N_u))};
571 this->element_, N_u);
574 ShapeFunctionDisplacement::NPOINTS,
576 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
578 double const T_ip = N * T;
579 double const T_prev_ip = N * T_prev;
582 double const p_cap_ip = -N * p_L;
583 double const p_cap_prev_ip = -N * p_L_prev;
588 constitutive_setting.eval(models,
591 {T_ip, T_prev_ip, grad_T_ip},
592 {p_cap_ip, p_cap_prev_ip, grad_p_cap_ip},
593 eps, current_state, this->prev_states_[ip],
594 this->material_states_[ip], tmp, output_data,
599 ShapeFunction,
typename ShapeFunctionDisplacement::MeshElement,
600 DisplacementDim>(this->element_, this->is_axially_symmetric_, p_L,
601 *process_data.pressure_interpolated);
603 ShapeFunction,
typename ShapeFunctionDisplacement::MeshElement,
604 DisplacementDim>(this->element_, this->is_axially_symmetric_, T,
605 *process_data.temperature_interpolated);