OGS
ThermoRichardsMechanicsFEM-impl.h
Go to the documentation of this file.
1
12#pragma once
13
14#include <spdlog/fmt/bundled/format.h>
15
16#include <cassert>
17#include <type_traits>
18
19#include "MaterialLib/MPL/Medium.h"
20#include "MathLib/KelvinVector.h"
21#include "NumLib/Fem/Interpolation.h"
22#include "ProcessLib/Deformation/LinearBMatrix.h"
23#include "ProcessLib/Graph/Get.h"
24#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EqP.h"
25#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EqT.h"
26#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/Gravity.h"
27#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/LiquidDensity.h"
28#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/LiquidViscosity.h"
29#include "ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TotalStressData.h"
30#include "ThermoRichardsMechanicsFEM.h"
31
32namespace ProcessLib
33{
34namespace ThermoRichardsMechanics
35{
36template <typename ShapeFunctionDisplacement, typename ShapeFunction,
37 int DisplacementDim, typename ConstitutiveTraits>
38ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunction,
39 DisplacementDim, ConstitutiveTraits>::
40 ThermoRichardsMechanicsLocalAssembler(
41 MeshLib::Element const& e,
42 std::size_t const /*local_matrix_size*/,
43 NumLib::GenericIntegrationMethod const& integration_method,
44 bool const is_axially_symmetric,
45 ThermoRichardsMechanicsProcessData<DisplacementDim, ConstitutiveTraits>&
46 process_data)
47 : LocalAssemblerInterface<DisplacementDim, ConstitutiveTraits>(
48 e, integration_method, is_axially_symmetric, process_data)
49{
50 unsigned const n_integration_points =
51 this->integration_method_.getNumberOfPoints();
52
53 ip_data_.resize(n_integration_points);
54
55 auto const shape_matrices_u =
56 NumLib::initShapeMatrices<ShapeFunctionDisplacement,
57 ShapeMatricesTypeDisplacement,
58 DisplacementDim>(e, is_axially_symmetric,
59 integration_method);
60
61 auto const shape_matrices =
62 NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
63 DisplacementDim>(e, is_axially_symmetric,
64 integration_method);
65
66 for (unsigned ip = 0; ip < n_integration_points; ip++)
67 {
68 auto& ip_data = ip_data_[ip];
69 auto const& sm_u = shape_matrices_u[ip];
70 ip_data_[ip].integration_weight =
71 integration_method.getWeightedPoint(ip).getWeight() *
72 sm_u.integralMeasure * sm_u.detJ;
73
74 ip_data.N_u = sm_u.N;
75 ip_data.dNdx_u = sm_u.dNdx;
76
77 // ip_data.N_p and ip_data.dNdx_p are used for both p and T variables
78 ip_data.N_p = shape_matrices[ip].N;
79 ip_data.dNdx_p = shape_matrices[ip].dNdx;
80 }
81}
82
83template <typename ShapeFunctionDisplacement, typename ShapeFunction,
84 int DisplacementDim, typename ConstitutiveTraits>
85void ThermoRichardsMechanicsLocalAssembler<
86 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
87 ConstitutiveTraits>::setInitialConditionsConcrete(Eigen::VectorXd const
88 local_x,
89 double const t,
90 int const /*process_id*/)
91{
92 assert(local_x.size() ==
93 temperature_size + pressure_size + displacement_size);
94
95 auto const p_L = local_x.template segment<pressure_size>(pressure_index);
96 auto const T =
97 local_x.template segment<temperature_size>(temperature_index);
98
99 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
100 auto const& medium =
101 *this->process_data_.media_map.getMedium(this->element_.getID());
102
103 MediaData const media_data{medium};
104
105 typename ConstitutiveTraits::ConstitutiveSetting const constitutive_setting;
106 auto models = ConstitutiveTraits::createConstitutiveModels(
107 this->process_data_, this->solid_material_);
108
109 unsigned const n_integration_points =
110 this->integration_method_.getNumberOfPoints();
111 for (unsigned ip = 0; ip < n_integration_points; ip++)
112 {
113 // N is used for both T and p variables.
114 auto const& N = ip_data_[ip].N_p;
115
116 ParameterLib::SpatialPosition const x_position{
117 std::nullopt, this->element_.getID(),
118 MathLib::Point3d(
119 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
120 ShapeMatricesTypeDisplacement>(
121 this->element_, ip_data_[ip].N_u))};
122
123 double T_ip;
124 NumLib::shapeFunctionInterpolate(T, N, T_ip);
125
126 double p_cap_ip;
127 NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
128
129 MPL::VariableArray variables;
130
131 if (medium.hasProperty(MPL::PropertyType::saturation))
132 {
133 variables.capillary_pressure = p_cap_ip;
134 variables.liquid_phase_pressure = -p_cap_ip;
135 // setting pG to 1 atm
136 // TODO : rewrite equations s.t. p_L = pG-p_cap
137 variables.gas_phase_pressure = 1.0e5;
138 variables.temperature = T_ip;
139
140 double const S_L =
141 medium.property(MPL::PropertyType::saturation)
142 .template value<double>(variables, x_position, t, dt);
143 std::get<PrevState<SaturationData>>(this->prev_states_[ip])->S_L =
144 S_L;
145
146 variables.liquid_saturation = S_L;
147 }
148
149 constitutive_setting.init(models, t, dt, x_position, media_data,
150 {T_ip, 0, {}}, this->current_states_[ip],
151 this->prev_states_[ip]);
152
153 if (this->process_data_.initial_stress.value)
154 {
155 if (!medium.hasProperty(MPL::PropertyType::saturation))
156 {
157 OGS_FATAL(
158 "The medium has no saturation property required to compute "
159 "initial stress.");
160 }
161 convertInitialStressType(ip, t, x_position, medium, variables,
162 -p_cap_ip);
163 }
164 }
165}
166
167template <typename ShapeFunctionDisplacement, typename ShapeFunction,
168 int DisplacementDim, typename ConstitutiveTraits>
169void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
170 ShapeFunction, DisplacementDim,
171 ConstitutiveTraits>::
172 convertInitialStressType(unsigned const ip,
173 double const t,
174 ParameterLib::SpatialPosition const x_position,
175 MaterialPropertyLib::Medium const& medium,
176 MPL::VariableArray const& variables,
177 double const p_at_ip)
178{
179 bool constexpr is_strain_temperature_constitutive =
180 std::is_same<ConstitutiveStress_StrainTemperature::ConstitutiveTraits<
181 DisplacementDim>,
182 ConstitutiveTraits>::value;
183 if (is_strain_temperature_constitutive &&
184 this->process_data_.initial_stress.type ==
185 InitialStress::Type::Effective)
186 {
187 return;
188 }
189
190 if (!is_strain_temperature_constitutive &&
191 this->process_data_.initial_stress.type == InitialStress::Type::Total)
192 {
193 return;
194 }
195
196 double const alpha_b =
197 medium.property(MPL::PropertyType::biot_coefficient)
198 .template value<double>(variables, x_position, t, 0.0 /*dt*/);
199
200 double const bishop =
201 medium.property(MPL::PropertyType::bishops_effective_stress)
202 .template value<double>(variables, x_position, t, 0.0 /*dt*/);
203
204 ConstitutiveTraits::ConstitutiveSetting::convertInitialStressType(
205 this->current_states_[ip], this->prev_states_[ip],
206 bishop * alpha_b * p_at_ip * Invariants::identity2);
207}
208
209template <typename ShapeFunctionDisplacement, typename ShapeFunction,
210 int DisplacementDim, typename ConstitutiveTraits>
211void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
212 ShapeFunction, DisplacementDim,
213 ConstitutiveTraits>::
214 assembleWithJacobian(double const t, double const dt,
215 std::vector<double> const& local_x,
216 std::vector<double> const& local_x_prev,
217 std::vector<double>& local_rhs_data,
218 std::vector<double>& local_Jac_data)
219{
220 auto& medium =
221 *this->process_data_.media_map.getMedium(this->element_.getID());
222
223 LocalMatrices loc_mat;
224 loc_mat.setZero();
225 LocalMatrices loc_mat_current_ip;
226 loc_mat_current_ip.setZero(); // only to set the right matrix sizes
227
228 typename ConstitutiveTraits::ConstitutiveSetting constitutive_setting;
229
230 for (unsigned ip = 0; ip < this->integration_method_.getNumberOfPoints();
231 ++ip)
232 {
233 ParameterLib::SpatialPosition const x_position{
234 std::nullopt, this->element_.getID(),
235 MathLib::Point3d(
236 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
237 ShapeMatricesTypeDisplacement>(
238 this->element_, ip_data_[ip].N_u))};
239
240 assembleWithJacobianSingleIP(
241 t, dt, x_position, //
242 local_x, local_x_prev, //
243 ip_data_[ip], constitutive_setting,
244 medium, //
245 loc_mat_current_ip, //
246 this->current_states_[ip], this->prev_states_[ip],
247 this->material_states_[ip], this->output_data_[ip]);
248 loc_mat += loc_mat_current_ip;
249 }
250
251 massLumping(loc_mat);
252
253 addToLocalMatrixData(dt, local_x, local_x_prev, loc_mat, local_rhs_data,
254 local_Jac_data);
255}
256
257template <typename ShapeFunctionDisplacement, typename ShapeFunction,
258 int DisplacementDim, typename ConstitutiveTraits>
259void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
260 ShapeFunction, DisplacementDim,
261 ConstitutiveTraits>::
262 massLumping(typename ThermoRichardsMechanicsLocalAssembler<
263 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
264 ConstitutiveTraits>::LocalMatrices& loc_mat) const
265{
266 if (this->process_data_.apply_mass_lumping)
267 {
268 loc_mat.storage_p_a_p =
269 loc_mat.storage_p_a_p.colwise().sum().eval().asDiagonal();
270 loc_mat.storage_p_a_S =
271 loc_mat.storage_p_a_S.colwise().sum().eval().asDiagonal();
272 loc_mat.storage_p_a_S_Jpp =
273 loc_mat.storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
274 }
275}
276
277template <typename ShapeFunctionDisplacement, typename ShapeFunction,
278 int DisplacementDim, typename ConstitutiveTraits>
279void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
280 ShapeFunction, DisplacementDim,
281 ConstitutiveTraits>::
282 addToLocalMatrixData(
283 double const dt,
284 std::vector<double> const& local_x,
285 std::vector<double> const& local_x_prev,
286 typename ThermoRichardsMechanicsLocalAssembler<
287 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
288 ConstitutiveTraits>::LocalMatrices const& loc_mat,
289 std::vector<double>& local_rhs_data,
290 std::vector<double>& local_Jac_data) const
291{
292 constexpr auto local_matrix_dim =
293 displacement_size + pressure_size + temperature_size;
294 assert(local_x.size() == local_matrix_dim);
295
296 auto local_Jac = MathLib::createZeroedMatrix<
297 typename ShapeMatricesTypeDisplacement::template MatrixType<
298 local_matrix_dim, local_matrix_dim>>(
299 local_Jac_data, local_matrix_dim, local_matrix_dim);
300
301 auto local_rhs =
302 MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
303 template VectorType<local_matrix_dim>>(
304 local_rhs_data, local_matrix_dim);
305
306 local_Jac.noalias() = loc_mat.Jac;
307 local_rhs.noalias() = -loc_mat.res;
308
309 //
310 // -- Jacobian
311 //
312 block_TT(local_Jac).noalias() += loc_mat.M_TT / dt + loc_mat.K_TT;
313 block_Tp(local_Jac).noalias() +=
314 loc_mat.K_Tp + loc_mat.dK_TT_dp + loc_mat.M_Tp / dt;
315
316 block_pT(local_Jac).noalias() += loc_mat.M_pT / dt + loc_mat.K_pT;
317 block_pp(local_Jac).noalias() +=
318 loc_mat.K_pp + loc_mat.storage_p_a_p / dt + loc_mat.storage_p_a_S_Jpp;
319 block_pu(local_Jac).noalias() = loc_mat.M_pu / dt;
320
321 //
322 // -- Residual
323 //
324 auto const [T, p_L, u] = localDOF(local_x);
325 auto const [T_prev, p_L_prev, u_prev] = localDOF(local_x_prev);
326
327 block_T(local_rhs).noalias() -= loc_mat.M_TT * (T - T_prev) / dt +
328 loc_mat.K_TT * T + loc_mat.K_Tp * p_L +
329 loc_mat.M_Tp * (p_L - p_L_prev) / dt;
330 block_p(local_rhs).noalias() -=
331 loc_mat.K_pp * p_L + loc_mat.K_pT * T +
332 (loc_mat.storage_p_a_p + loc_mat.storage_p_a_S) * (p_L - p_L_prev) /
333 dt +
334 loc_mat.M_pu * (u - u_prev) / dt + loc_mat.M_pT * (T - T_prev) / dt;
335}
336
337template <typename ShapeFunctionDisplacement, typename ShapeFunction,
338 int DisplacementDim, typename ConstitutiveTraits>
339void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
340 ShapeFunction, DisplacementDim,
341 ConstitutiveTraits>::
342 assembleWithJacobianSingleIP(
343 double const t, double const dt,
344 ParameterLib::SpatialPosition const& x_position,
345 std::vector<double> const& local_x,
346 std::vector<double> const& local_x_prev,
347 typename ThermoRichardsMechanicsLocalAssembler<
348 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
349 ConstitutiveTraits>::IpData const& ip_data,
350 typename ConstitutiveTraits::ConstitutiveSetting& CS,
351 MaterialPropertyLib::Medium& medium,
352 typename ThermoRichardsMechanicsLocalAssembler<
353 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
354 ConstitutiveTraits>::LocalMatrices& out,
355 typename ConstitutiveTraits::StatefulData& current_state,
356 typename ConstitutiveTraits::StatefulDataPrev const& prev_state,
357 MaterialStateData<DisplacementDim>& mat_state,
358 typename ConstitutiveTraits::OutputData& output_data) const
359{
360 auto const& N_u = ip_data.N_u;
361 auto const& dNdx_u = ip_data.dNdx_u;
362
363 // N and dNdx are used for both p and T variables
364 auto const& N = ip_data.N_p;
365 auto const& dNdx = ip_data.dNdx_p;
366
367 auto const B =
368 LinearBMatrix::computeBMatrix<DisplacementDim,
369 ShapeFunctionDisplacement::NPOINTS,
370 typename BMatricesType::BMatrixType>(
371 dNdx_u, N_u, (*x_position.getCoordinates())[0],
372 this->is_axially_symmetric_);
373
374 typename ConstitutiveTraits::ConstitutiveData CD;
375
376 auto const [T, p_L, u] = localDOF(local_x);
377 auto const [T_prev, p_L_prev, u_prev] = localDOF(local_x_prev);
378
379 {
380 auto models = ConstitutiveTraits::createConstitutiveModels(
381 this->process_data_, this->solid_material_);
382 typename ConstitutiveTraits::ConstitutiveTempData tmp;
383
384 double const T_ip = N * T;
385 double const T_prev_ip = N * T_prev;
386 GlobalDimVectorType const grad_T_ip = dNdx * T;
387
388 double const p_cap_ip = -N * p_L;
389 double const p_cap_prev_ip = -N * p_L_prev;
390 GlobalDimVectorType const grad_p_cap_ip = -dNdx * p_L;
391
392 KelvinVectorType eps = B * u;
393
394 CS.eval(models, t, dt, x_position, //
395 medium, //
396 {T_ip, T_prev_ip, grad_T_ip}, //
397 {p_cap_ip, p_cap_prev_ip, grad_p_cap_ip}, //
398 eps, current_state, prev_state, mat_state, tmp, output_data,
399 CD);
400 }
401
402 using NodalMatrix = typename ShapeMatricesType::NodalMatrixType;
403 NodalMatrix const NTN = N.transpose() * N;
404 NodalMatrix const dNTdN = dNdx.transpose() * dNdx;
405
406 // TODO is identity2.transpose() * B something like divergence?
407 auto const& identity2 = MathLib::KelvinVector::Invariants<
408 MathLib::KelvinVector::kelvin_vector_dimensions(
409 DisplacementDim)>::identity2;
410 typename ShapeMatricesTypeDisplacement::template MatrixType<
411 displacement_size, pressure_size>
412 BTI2N = B.transpose() * identity2 * N;
413
414 /*
415 * Conventions:
416 *
417 * * use positive signs exclusively, any negative signs must be included in
418 * the coefficients coming from the constitutive setting
419 * * the used coefficients from the constitutive setting are named after the
420 * terms they appear in, e.g. K_TT_X_dNTdN means it is a scalar (X) that
421 * will be multiplied by dNdx.transpose() * dNdx. Placefolders for the
422 * coefficients are:
423 * * X -> scalar
424 * * V -> vector
425 * * K -> Kelvin vector
426 * * the Laplace terms have a special name, e.g., K_TT_Laplace
427 * * there shall be only one contribution to each of the LocalMatrices,
428 * assigned with = assignment; this point might be relaxed in the future
429 * * this method will overwrite the data in the passed LocalMatrices& out
430 * argument, not add to it
431 */
432
433 // residual, order T, p, u
434 block_p(out.res).noalias() =
435 dNdx.transpose() * std::get<EqPData<DisplacementDim>>(CD).rhs_p_dNT_V;
436 block_u(out.res).noalias() =
437 B.transpose() *
438 ProcessLib::Graph::get<TotalStressData<DisplacementDim>>(
439 CD, current_state)
440 .sigma_total -
441 static_cast<int>(this->process_data_.apply_body_force_for_deformation) *
442 N_u_op(N_u).transpose() *
443 std::get<GravityData<DisplacementDim>>(CD).volumetric_body_force;
444
445 // Storage matrices
446 out.storage_p_a_p.noalias() =
447 std::get<EqPData<DisplacementDim>>(CD).storage_p_a_p_X_NTN * NTN;
448 out.storage_p_a_S.noalias() =
449 std::get<TRMStorageData>(CD).storage_p_a_S_X_NTN * NTN;
450 out.storage_p_a_S_Jpp.noalias() =
451 std::get<TRMStorageData>(CD).storage_p_a_S_Jpp_X_NTN * NTN;
452
453 // M matrices, order T, p, u
454 out.M_TT.noalias() =
455 std::get<EqTData<DisplacementDim>>(CD).M_TT_X_NTN * NTN;
456 out.M_Tp.noalias() =
457 std::get<TRMVaporDiffusionData<DisplacementDim>>(CD).M_Tp_X_NTN * NTN;
458
459 out.M_pT.noalias() =
460 std::get<EqPData<DisplacementDim>>(CD).M_pT_X_NTN * NTN;
461 out.M_pu.noalias() =
462 std::get<EqPData<DisplacementDim>>(CD).M_pu_X_BTI2N * BTI2N.transpose();
463
464 // K matrices, order T, p, u
465 out.K_TT.noalias() =
466 dNdx.transpose() *
467 std::get<TRMHeatStorageAndFluxData<DisplacementDim>>(CD)
468 .K_TT_Laplace *
469 dNdx +
470 N.transpose() *
471 (std::get<EqTData<DisplacementDim>>(CD).K_TT_NT_V_dN.transpose() *
472 dNdx) +
473 std::get<TRMVaporDiffusionData<DisplacementDim>>(CD).K_TT_X_dNTdN *
474 dNTdN;
475
476 out.dK_TT_dp.noalias() =
477 N.transpose() *
478 (std::get<TRMHeatStorageAndFluxData<DisplacementDim>>(CD)
479 .K_Tp_NT_V_dN.transpose() *
480 dNdx) +
481 std::get<TRMHeatStorageAndFluxData<DisplacementDim>>(CD).K_Tp_X_NTN *
482 NTN;
483 out.K_Tp.noalias() =
484 dNdx.transpose() *
485 std::get<ThermoOsmosisData<DisplacementDim>>(CD).K_Tp_Laplace *
486 dNdx +
487 std::get<TRMVaporDiffusionData<DisplacementDim>>(CD).K_Tp_X_dNTdN *
488 dNTdN;
489
490 out.K_pp.noalias() =
491 dNdx.transpose() * std::get<EqPData<DisplacementDim>>(CD).K_pp_Laplace *
492 dNdx +
493 std::get<TRMVaporDiffusionData<DisplacementDim>>(CD).K_pp_X_dNTdN *
494 dNTdN;
495 out.K_pT.noalias() =
496 dNdx.transpose() *
497 std::get<ThermoOsmosisData<DisplacementDim>>(CD).K_pT_Laplace * dNdx;
498
499 // direct Jacobian contributions, order T, p, u
500 block_pT(out.Jac).noalias() =
501 std::get<TRMVaporDiffusionData<DisplacementDim>>(CD).J_pT_X_dNTdN *
502 dNTdN;
503 block_pp(out.Jac).noalias() =
504 std::get<TRMStorageData>(CD).J_pp_X_NTN * NTN +
505 std::get<EqPData<DisplacementDim>>(CD).J_pp_X_BTI2NT_u_dot_N *
506 BTI2N.transpose() * (u - u_prev) / dt *
507 N // TODO something with volumetric strain rate?
508 + dNdx.transpose() *
509 std::get<EqPData<DisplacementDim>>(CD).J_pp_dNT_V_N * N;
510
511 block_uT(out.Jac).noalias() =
512 B.transpose() *
513 std::get<SolidMechanicsDataStateless<DisplacementDim>>(CD).J_uT_BT_K_N *
514 N;
515 block_up(out.Jac).noalias() =
516 B.transpose() *
517 std::get<SolidMechanicsDataStateless<DisplacementDim>>(CD)
518 .J_up_BT_K_N *
519 N +
520 N_u_op(N_u).transpose() *
521 std::get<GravityData<DisplacementDim>>(CD).J_up_HT_V_N * N;
522 block_uu(out.Jac).noalias() =
523 B.transpose() *
524 std::get<SolidMechanicsDataStateless<DisplacementDim>>(CD)
525 .stiffness_tensor *
526 B;
527
528 out *= ip_data.integration_weight;
529}
530
531template <typename ShapeFunctionDisplacement, typename ShapeFunction,
532 int DisplacementDim, typename ConstitutiveTraits>
533void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
534 ShapeFunction, DisplacementDim,
535 ConstitutiveTraits>::
536 computeSecondaryVariableConcrete(double const t, double const dt,
537 Eigen::VectorXd const& local_x,
538 Eigen::VectorXd const& local_x_prev)
539{
540 auto const T = block_T(local_x);
541 auto const p_L = block_p(local_x);
542 auto const u = block_u(local_x);
543
544 auto const T_prev = block_T(local_x_prev);
545 auto const p_L_prev = block_p(local_x_prev);
546
547 auto const e_id = this->element_.getID();
548 auto const& process_data = this->process_data_;
549 auto& medium = *process_data.media_map.getMedium(e_id);
550
551 unsigned const n_integration_points =
552 this->integration_method_.getNumberOfPoints();
553
554 typename ConstitutiveTraits::ConstitutiveSetting constitutive_setting;
555
556 auto models = ConstitutiveTraits::createConstitutiveModels(
557 process_data, this->solid_material_);
558 typename ConstitutiveTraits::ConstitutiveTempData tmp;
559 typename ConstitutiveTraits::ConstitutiveData CD;
560
561 for (unsigned ip = 0; ip < n_integration_points; ip++)
562 {
563 auto& current_state = this->current_states_[ip];
564 auto& output_data = this->output_data_[ip];
565
566 auto const& ip_data = ip_data_[ip];
567
568 // N is used for both p and T variables
569 auto const& N = ip_data.N_p;
570 auto const& N_u = ip_data.N_u;
571 auto const& dNdx_u = ip_data.dNdx_u;
572 auto const& dNdx = ip_data.dNdx_p;
573
574 ParameterLib::SpatialPosition const x_position{
575 std::nullopt, this->element_.getID(),
576 MathLib::Point3d(
577 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
578 ShapeMatricesTypeDisplacement>(
579 this->element_, N_u))};
580 auto const x_coord =
581 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
582 ShapeMatricesTypeDisplacement>(
583 this->element_, N_u);
584 auto const B =
585 LinearBMatrix::computeBMatrix<DisplacementDim,
586 ShapeFunctionDisplacement::NPOINTS,
587 typename BMatricesType::BMatrixType>(
588 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
589
590 double const T_ip = N * T;
591 double const T_prev_ip = N * T_prev;
592 GlobalDimVectorType const grad_T_ip = dNdx * T;
593
594 double const p_cap_ip = -N * p_L;
595 double const p_cap_prev_ip = -N * p_L_prev;
596 GlobalDimVectorType const grad_p_cap_ip = -dNdx * p_L;
597
598 KelvinVectorType eps = B * u;
599
600 constitutive_setting.eval(models, //
601 t, dt, x_position, //
602 medium, //
603 {T_ip, T_prev_ip, grad_T_ip}, //
604 {p_cap_ip, p_cap_prev_ip, grad_p_cap_ip}, //
605 eps, current_state, this->prev_states_[ip],
606 this->material_states_[ip], tmp, output_data,
607 CD);
608 }
609
610 NumLib::interpolateToHigherOrderNodes<
611 ShapeFunction, typename ShapeFunctionDisplacement::MeshElement,
612 DisplacementDim>(this->element_, this->is_axially_symmetric_, p_L,
613 *process_data.pressure_interpolated);
614 NumLib::interpolateToHigherOrderNodes<
615 ShapeFunction, typename ShapeFunctionDisplacement::MeshElement,
616 DisplacementDim>(this->element_, this->is_axially_symmetric_, T,
617 *process_data.temperature_interpolated);
618}
619} // namespace ThermoRichardsMechanics
620} // namespace ProcessLib
OGS_FATAL
#define OGS_FATAL(...)
Definition Error.h:26
MaterialPropertyLib::Medium::property
Property const & property(PropertyType const &p) const
Definition Medium.cpp:54
MathLib::WeightedPoint::getWeight
constexpr double getWeight() const
Definition WeightedPoint.h:82
NumLib::GenericIntegrationMethod::getWeightedPoint
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
Definition GenericIntegrationMethod.h:36
ParameterLib::SpatialPosition::getCoordinates
std::optional< MathLib::Point3d > const getCoordinates() const
Definition SpatialPosition.h:63
ProcessLib::BMatrixPolicyType::BMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition BMatrixPolicy.h:55
ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler::ShapeMatricesTypeDisplacement
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
Definition ThermoRichardsMechanicsFEM.h:49
ProcessLib::ThermoRichardsMechanics::ThermoRichardsMechanicsLocalAssembler::ShapeMatricesType
ShapeMatrixPolicyType< ShapeFunction, DisplacementDim > ShapeMatricesType
Definition ThermoRichardsMechanicsFEM.h:53
MathLib::KelvinVector::kelvin_vector_dimensions
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition KelvinVector.h:24
MathLib::createZeroedVector
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Definition EigenMapTools.h:153
MathLib::createZeroedMatrix
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition EigenMapTools.h:32
NumLib::detail::shapeFunctionInterpolate
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Definition Interpolation.h:27
NumLib::interpolateXCoordinate
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition InitShapeMatrices.h:91
NumLib::interpolateToHigherOrderNodes
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
Definition Interpolation.h:107
NumLib::initShapeMatrices
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
Definition InitShapeMatrices.h:54
NumLib::interpolateCoordinates
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition InitShapeMatrices.h:102
ProcessLib::Graph::get
auto & get(Tuples &... ts)
Definition Get.h:59
ProcessLib::LinearBMatrix::computeBMatrix
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
Definition LinearBMatrix.h:44
EigenFixedShapeMatrixPolicy::NodalMatrixType
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
Definition ShapeMatrixPolicy.h:74
ProcessLib::ThermoRichardsMechanics::ConstitutiveStress_StrainTemperature::ConstitutiveTraits
Definition Traits.h:22