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
32
33namespace ProcessLib
34{
35namespace ThermoRichardsMechanics
36{
37template <typename ShapeFunctionDisplacement, typename ShapeFunction,
38 int DisplacementDim, typename ConstitutiveTraits>
39ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunction,
40 DisplacementDim, ConstitutiveTraits>::
41 ThermoRichardsMechanicsLocalAssembler(
42 MeshLib::Element const& e,
43 std::size_t const /*local_matrix_size*/,
44 NumLib::GenericIntegrationMethod const& integration_method,
45 bool const is_axially_symmetric,
47 process_data)
48 : LocalAssemblerInterface<DisplacementDim, ConstitutiveTraits>(
49 e, integration_method, is_axially_symmetric, process_data)
50{
51 unsigned const n_integration_points =
53
54 ip_data_.resize(n_integration_points);
55
56 auto const shape_matrices_u =
57 NumLib::initShapeMatrices<ShapeFunctionDisplacement,
59 DisplacementDim>(e, is_axially_symmetric,
60 integration_method);
61
62 auto const shape_matrices =
64 DisplacementDim>(e, is_axially_symmetric,
65 integration_method);
66
67 for (unsigned ip = 0; ip < n_integration_points; ip++)
68 {
69 auto& ip_data = ip_data_[ip];
70 auto const& sm_u = shape_matrices_u[ip];
71 ip_data_[ip].integration_weight =
72 integration_method.getWeightedPoint(ip).getWeight() *
73 sm_u.integralMeasure * sm_u.detJ;
74
75 ip_data.N_u = sm_u.N;
76 ip_data.dNdx_u = sm_u.dNdx;
77
78 // ip_data.N_p and ip_data.dNdx_p are used for both p and T variables
79 ip_data.N_p = shape_matrices[ip].N;
80 ip_data.dNdx_p = shape_matrices[ip].dNdx;
81 }
82}
83
84template <typename ShapeFunctionDisplacement, typename ShapeFunction,
85 int DisplacementDim, typename ConstitutiveTraits>
87 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
88 ConstitutiveTraits>::setInitialConditionsConcrete(Eigen::VectorXd const
89 local_x,
90 double const t,
91 int const /*process_id*/)
92{
93 assert(local_x.size() ==
94 temperature_size + pressure_size + displacement_size);
95
96 auto const p_L = local_x.template segment<pressure_size>(pressure_index);
97 auto const T =
98 local_x.template segment<temperature_size>(temperature_index);
99
100 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
101 auto const& medium =
102 *this->process_data_.media_map.getMedium(this->element_.getID());
103
104 MediaData const media_data{medium};
105
106 typename ConstitutiveTraits::ConstitutiveSetting const constitutive_setting;
107 typename ConstitutiveTraits::ConstitutiveModels models(
108 this->process_data_, this->solid_material_);
109
110 unsigned const n_integration_points =
111 this->integration_method_.getNumberOfPoints();
112 for (unsigned ip = 0; ip < n_integration_points; ip++)
113 {
114 // N is used for both T and p variables.
115 auto const& N = ip_data_[ip].N_p;
116
117 ParameterLib::SpatialPosition const x_position{
118 std::nullopt, this->element_.getID(), ip,
120 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
122 this->element_, ip_data_[ip].N_u))};
123
124 double T_ip;
126
127 double p_cap_ip;
128 NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
129
130 MPL::VariableArray variables;
131 variables.capillary_pressure = p_cap_ip;
132 variables.liquid_phase_pressure = -p_cap_ip;
133 // setting pG to 1 atm
134 // TODO : rewrite equations s.t. p_L = pG-p_cap
135 variables.gas_phase_pressure = 1.0e5;
136 variables.temperature = T_ip;
137
138 double const S_L =
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;
142
143 constitutive_setting.init(models, t, dt, x_position, media_data,
144 {T_ip, 0, {}}, this->current_states_[ip],
145 this->prev_states_[ip]);
146
147 if (this->process_data_.initial_stress.value)
148 {
149 variables.liquid_saturation = S_L;
150 convertInitialStressType(ip, t, x_position, medium, variables,
151 -p_cap_ip);
152 }
153 }
154}
155
156template <typename ShapeFunctionDisplacement, typename ShapeFunction,
157 int DisplacementDim, typename ConstitutiveTraits>
158void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
159 ShapeFunction, DisplacementDim,
160 ConstitutiveTraits>::
161 convertInitialStressType(unsigned const ip,
162 double const t,
163 ParameterLib::SpatialPosition const x_position,
164 MaterialPropertyLib::Medium const& medium,
165 MPL::VariableArray const& variables,
166 double const p_at_ip)
167{
168 bool constexpr is_strain_temperature_constitutive =
170 DisplacementDim>,
171 ConstitutiveTraits>::value;
172 if (is_strain_temperature_constitutive &&
173 this->process_data_.initial_stress.type ==
175 {
176 return;
177 }
178
179 if (!is_strain_temperature_constitutive &&
180 this->process_data_.initial_stress.type == InitialStress::Type::Total)
181 {
182 return;
183 }
184
185 double const alpha_b =
186 medium.property(MPL::PropertyType::biot_coefficient)
187 .template value<double>(variables, x_position, t, 0.0 /*dt*/);
188
189 double const bishop =
190 medium.property(MPL::PropertyType::bishops_effective_stress)
191 .template value<double>(variables, x_position, t, 0.0 /*dt*/);
192
193 ConstitutiveTraits::ConstitutiveSetting::convertInitialStressType(
194 this->current_states_[ip], this->prev_states_[ip],
195 bishop * alpha_b * p_at_ip * Invariants::identity2);
196}
197
198template <typename ShapeFunctionDisplacement, typename ShapeFunction,
199 int DisplacementDim, typename ConstitutiveTraits>
200void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
201 ShapeFunction, DisplacementDim,
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>& /*local_M_data*/,
207 std::vector<double>& /*local_K_data*/,
208 std::vector<double>& local_rhs_data,
209 std::vector<double>& local_Jac_data)
210{
211 auto& medium =
212 *this->process_data_.media_map.getMedium(this->element_.getID());
213
214 LocalMatrices loc_mat;
215 loc_mat.setZero();
216 LocalMatrices loc_mat_current_ip;
217 loc_mat_current_ip.setZero(); // only to set the right matrix sizes
218
219 typename ConstitutiveTraits::ConstitutiveSetting constitutive_setting;
220
221 for (unsigned ip = 0; ip < this->integration_method_.getNumberOfPoints();
222 ++ip)
223 {
224 ParameterLib::SpatialPosition const x_position{
225 std::nullopt, this->element_.getID(), ip,
227 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
229 this->element_, ip_data_[ip].N_u))};
230
231 assembleWithJacobianSingleIP(
232 t, dt, x_position, //
233 local_x, local_x_prev, //
234 ip_data_[ip], constitutive_setting,
235 medium, //
236 loc_mat_current_ip, //
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;
240 }
241
242 massLumping(loc_mat);
243
244 addToLocalMatrixData(dt, local_x, local_x_prev, loc_mat, local_rhs_data,
245 local_Jac_data);
246}
247
248template <typename ShapeFunctionDisplacement, typename ShapeFunction,
249 int DisplacementDim, typename ConstitutiveTraits>
250void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
251 ShapeFunction, DisplacementDim,
252 ConstitutiveTraits>::
253 massLumping(typename ThermoRichardsMechanicsLocalAssembler<
254 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
255 ConstitutiveTraits>::LocalMatrices& loc_mat) const
256{
257 if (this->process_data_.apply_mass_lumping)
258 {
259 loc_mat.storage_p_a_p =
260 loc_mat.storage_p_a_p.colwise().sum().eval().asDiagonal();
261 loc_mat.storage_p_a_S =
262 loc_mat.storage_p_a_S.colwise().sum().eval().asDiagonal();
263 loc_mat.storage_p_a_S_Jpp =
264 loc_mat.storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
265 }
266}
267
268template <typename ShapeFunctionDisplacement, typename ShapeFunction,
269 int DisplacementDim, typename ConstitutiveTraits>
270void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
271 ShapeFunction, DisplacementDim,
272 ConstitutiveTraits>::
273 addToLocalMatrixData(
274 double const dt,
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
282{
283 constexpr auto local_matrix_dim =
284 displacement_size + pressure_size + temperature_size;
285 assert(local_x.size() == local_matrix_dim);
286
287 auto local_Jac = MathLib::createZeroedMatrix<
288 typename ShapeMatricesTypeDisplacement::template MatrixType<
289 local_matrix_dim, local_matrix_dim>>(
290 local_Jac_data, local_matrix_dim, local_matrix_dim);
291
292 auto local_rhs =
293 MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
294 template VectorType<local_matrix_dim>>(
295 local_rhs_data, local_matrix_dim);
296
297 local_Jac.noalias() = loc_mat.Jac;
298 local_rhs.noalias() = -loc_mat.res;
299
300 //
301 // -- Jacobian
302 //
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;
306
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;
311
312 //
313 // -- Residual
314 //
315 auto const [T, p_L, u] = localDOF(local_x);
316 auto const [T_prev, p_L_prev, u_prev] = localDOF(local_x_prev);
317
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) /
324 dt +
325 loc_mat.M_pu * (u - u_prev) / dt + loc_mat.M_pT * (T - T_prev) / dt;
326}
327
328template <typename ShapeFunctionDisplacement, typename ShapeFunction,
329 int DisplacementDim, typename ConstitutiveTraits>
330void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
331 ShapeFunction, DisplacementDim,
332 ConstitutiveTraits>::
333 assembleWithJacobianSingleIP(
334 double const t, double const dt,
335 ParameterLib::SpatialPosition const& x_position,
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
350{
351 auto const& N_u = ip_data.N_u;
352 auto const& dNdx_u = ip_data.dNdx_u;
353
354 // N and dNdx are used for both p and T variables
355 auto const& N = ip_data.N_p;
356 auto const& dNdx = ip_data.dNdx_p;
357
358 auto const B =
359 LinearBMatrix::computeBMatrix<DisplacementDim,
360 ShapeFunctionDisplacement::NPOINTS,
362 dNdx_u, N_u, (*x_position.getCoordinates())[0],
363 this->is_axially_symmetric_);
364
365 auto const [T, p_L, u] = localDOF(local_x);
366 auto const [T_prev, p_L_prev, u_prev] = localDOF(local_x_prev);
367
368 GlobalDimVectorType const grad_T_ip = dNdx * T;
369
370 typename ConstitutiveTraits::ConstitutiveModels models(
371 this->process_data_, this->solid_material_);
372 typename ConstitutiveTraits::ConstitutiveTempData tmp;
373 typename ConstitutiveTraits::ConstitutiveData CD;
374
375 {
376 double const T_ip = N * T;
377 double const T_prev_ip = N * T_prev;
378
379 double const p_cap_ip = -N * p_L;
380 double const p_cap_prev_ip = -N * p_L_prev;
381 GlobalDimVectorType const grad_p_cap_ip = -dNdx * p_L;
382
383 KelvinVectorType eps = B * u;
384
385 CS.eval(models, t, dt, x_position, //
386 medium, //
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,
390 CD);
391 }
392
393 using NodalMatrix = typename ShapeMatricesType::NodalMatrixType;
394 NodalMatrix const NTN = N.transpose() * N;
395 NodalMatrix const dNTdN = dNdx.transpose() * dNdx;
396
397 // TODO is identity2.transpose() * B something like divergence?
398 auto const& identity2 = MathLib::KelvinVector::Invariants<
400 DisplacementDim)>::identity2;
401 typename ShapeMatricesTypeDisplacement::template MatrixType<
402 displacement_size, pressure_size>
403 BTI2N = B.transpose() * identity2 * N;
404
405 /*
406 * Conventions:
407 *
408 * * use positive signs exclusively, any negative signs must be included in
409 * the coefficients coming from the constitutive setting
410 * * the used coefficients from the constitutive setting are named after the
411 * terms they appear in, e.g. K_TT_X_dNTdN means it is a scalar (X) that
412 * will be multiplied by dNdx.transpose() * dNdx. Placefolders for the
413 * coefficients are:
414 * * X -> scalar
415 * * V -> vector
416 * * K -> Kelvin vector
417 * * the Laplace terms have a special name, e.g., K_TT_Laplace
418 * * there shall be only one contribution to each of the LocalMatrices,
419 * assigned with = assignment; this point might be relaxed in the future
420 * * this method will overwrite the data in the passed LocalMatrices& out
421 * argument, not add to it
422 */
423
424 // residual, order T, p, u
425 block_p(out.res).noalias() =
426 dNdx.transpose() * std::get<EqPData<DisplacementDim>>(CD).rhs_p_dNT_V;
427 block_u(out.res).noalias() =
428 B.transpose() *
429 ProcessLib::Graph::get<TotalStressData<DisplacementDim>>(
430 CD, current_state)
431 .sigma_total -
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;
435
436 // Storage matrices
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;
443
444 // M matrices, order T, p, u
445 out.M_TT.noalias() =
446 std::get<EqTData<DisplacementDim>>(CD).M_TT_X_NTN * NTN;
447 out.M_Tp.noalias() =
448 std::get<TRMVaporDiffusionData<DisplacementDim>>(CD).M_Tp_X_NTN * NTN;
449
450 out.M_pT.noalias() =
451 std::get<EqPData<DisplacementDim>>(CD).M_pT_X_NTN * NTN;
452 out.M_pu.noalias() =
453 std::get<EqPData<DisplacementDim>>(CD).M_pu_X_BTI2N * BTI2N.transpose();
454
455 // K matrices, order T, p, u
456 out.K_TT.noalias() =
457 dNdx.transpose() *
458 std::get<TRMHeatStorageAndFluxData<DisplacementDim>>(CD)
459 .K_TT_Laplace *
460 dNdx +
461 N.transpose() *
462 (std::get<EqTData<DisplacementDim>>(CD).K_TT_NT_V_dN.transpose() *
463 dNdx) +
465 dNTdN;
466
467 out.dK_TT_dp.noalias() =
468 N.transpose() *
469 (std::get<TRMHeatStorageAndFluxData<DisplacementDim>>(CD)
470 .K_Tp_NT_V_dN.transpose() *
471 dNdx) +
473 NTN;
474 out.K_Tp.noalias() =
475 dNdx.transpose() *
476 std::get<ThermoOsmosisData<DisplacementDim>>(CD).K_Tp_Laplace *
477 dNdx +
479 dNTdN;
480
481 out.K_pp.noalias() =
482 dNdx.transpose() * std::get<EqPData<DisplacementDim>>(CD).K_pp_Laplace *
483 dNdx +
485 dNTdN;
486 out.K_pT.noalias() =
487 dNdx.transpose() *
488 std::get<ThermoOsmosisData<DisplacementDim>>(CD).K_pT_Laplace * dNdx;
489
490 // direct Jacobian contributions, order T, p, u
491 block_pT(out.Jac).noalias() =
492 std::get<TRMVaporDiffusionData<DisplacementDim>>(CD).J_pT_X_dNTdN *
493 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 *
498 N // TODO something with volumetric strain rate?
499 + dNdx.transpose() *
500 std::get<EqPData<DisplacementDim>>(CD).J_pp_dNT_V_N * N;
501
502 block_uT(out.Jac).noalias() =
503 B.transpose() *
504 std::get<SolidMechanicsDataStateless<DisplacementDim>>(CD).J_uT_BT_K_N *
505 N;
506 block_up(out.Jac).noalias() =
507 B.transpose() *
508 std::get<SolidMechanicsDataStateless<DisplacementDim>>(CD)
509 .J_up_BT_K_N *
510 N +
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() =
514 B.transpose() *
515 std::get<SolidMechanicsDataStateless<DisplacementDim>>(CD)
516 .stiffness_tensor *
517 B;
518
519 out *= ip_data.integration_weight;
520}
521
522template <typename ShapeFunctionDisplacement, typename ShapeFunction,
523 int DisplacementDim, typename ConstitutiveTraits>
524void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
525 ShapeFunction, DisplacementDim,
526 ConstitutiveTraits>::
527 computeSecondaryVariableConcrete(double const t, double const dt,
528 Eigen::VectorXd const& local_x,
529 Eigen::VectorXd const& local_x_prev)
530{
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);
534
535 auto const T_prev = block_T(local_x_prev);
536 auto const p_L_prev = block_p(local_x_prev);
537
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);
541
542 unsigned const n_integration_points =
543 this->integration_method_.getNumberOfPoints();
544
545 double saturation_avg = 0;
546 double porosity_avg = 0;
547 double liquid_density_avg = 0;
548 double viscosity_avg = 0;
549
551 KV sigma_avg = KV::Zero();
552
553 typename ConstitutiveTraits::ConstitutiveSetting constitutive_setting;
554
555 typename ConstitutiveTraits::ConstitutiveModels models(
556 process_data, this->solid_material_);
557 typename ConstitutiveTraits::ConstitutiveTempData tmp;
558 typename ConstitutiveTraits::ConstitutiveData CD;
559
560 for (unsigned ip = 0; ip < n_integration_points; ip++)
561 {
562 auto& current_state = this->current_states_[ip];
563 auto& output_data = this->output_data_[ip];
564
565 auto const& ip_data = ip_data_[ip];
566
567 // N is used for both p and T variables
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;
572
573 ParameterLib::SpatialPosition const x_position{
574 std::nullopt, this->element_.getID(), ip,
576 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
578 this->element_, N_u))};
579 auto const x_coord =
580 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
582 this->element_, N_u);
583 auto const B =
584 LinearBMatrix::computeBMatrix<DisplacementDim,
585 ShapeFunctionDisplacement::NPOINTS,
587 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
588
589 double const T_ip = N * T;
590 double const T_prev_ip = N * T_prev;
591 GlobalDimVectorType const grad_T_ip = dNdx * T;
592
593 double const p_cap_ip = -N * p_L;
594 double const p_cap_prev_ip = -N * p_L_prev;
595 GlobalDimVectorType const grad_p_cap_ip = -dNdx * p_L;
596
597 KelvinVectorType eps = B * u;
598
599 constitutive_setting.eval(models, //
600 t, dt, x_position, //
601 medium, //
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,
606 CD);
607
608 saturation_avg += std::get<SaturationData>(current_state).S_L;
609 porosity_avg += std::get<PorosityData>(current_state).phi;
610
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(
614 current_state);
615 }
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;
621
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;
626
627 Eigen::Map<KV>(
628 &(*process_data.element_stresses)[e_id * KV::RowsAtCompileTime]) =
630
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);
639}
640} // namespace ThermoRichardsMechanics
641} // namespace ProcessLib
Property const & property(PropertyType const &p) const
Definition Medium.cpp:54
double getWeight() const
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
std::optional< MathLib::Point3d > const & getCoordinates() const
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
ShapeMatrixPolicyType< ShapeFunction, DisplacementDim > ShapeMatricesType
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
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)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
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.
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType