OGS
TH2MFEM-impl.h
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
2// SPDX-License-Identifier: BSD-3-Clause
3
4#pragma once
5
6#include <Eigen/LU>
7
8#include "ConstitutiveRelations/ConstitutiveData.h"
9#include "MaterialLib/MPL/Medium.h"
10#include "MaterialLib/MPL/Phase.h"
11#include "MaterialLib/MPL/Property.h"
12#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
13#include "MaterialLib/PhysicalConstant.h"
14#include "MathLib/EigenBlockMatrixView.h"
15#include "MathLib/KelvinVector.h"
16#include "NumLib/Fem/Interpolation.h"
17#include "ProcessLib/Reflection/ReflectionSetIPData.h"
18#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
19#include "TH2MProcessData.h"
20
21namespace ProcessLib
22{
23namespace TH2M
24{
25namespace MPL = MaterialPropertyLib;
26
27template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
28 int DisplacementDim>
29TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
30 DisplacementDim>::
31 TH2MLocalAssembler(
32 MeshLib::Element const& e,
33 std::size_t const /*local_matrix_size*/,
34 NumLib::GenericIntegrationMethod const& integration_method,
35 bool const is_axially_symmetric,
36 TH2MProcessData<DisplacementDim>& process_data)
37 : LocalAssemblerInterface<DisplacementDim>(
38 e, integration_method, is_axially_symmetric, process_data)
39{
40 unsigned const n_integration_points =
41 this->integration_method_.getNumberOfPoints();
42
43 _ip_data.resize(n_integration_points);
44 _secondary_data.N_u.resize(n_integration_points);
45
46 auto const shape_matrices_u =
47 NumLib::initShapeMatrices<ShapeFunctionDisplacement,
48 ShapeMatricesTypeDisplacement,
49 DisplacementDim>(e, is_axially_symmetric,
50 this->integration_method_);
51
52 auto const shape_matrices_p =
53 NumLib::initShapeMatrices<ShapeFunctionPressure,
54 ShapeMatricesTypePressure, DisplacementDim>(
55 e, is_axially_symmetric, this->integration_method_);
56
57 for (unsigned ip = 0; ip < n_integration_points; ip++)
58 {
59 auto& ip_data = _ip_data[ip];
60 auto const& sm_u = shape_matrices_u[ip];
61 ip_data.integration_weight =
62 this->integration_method_.getWeightedPoint(ip).getWeight() *
63 sm_u.integralMeasure * sm_u.detJ;
64
65 ip_data.N_u = sm_u.N;
66 ip_data.dNdx_u = sm_u.dNdx;
67
68 ip_data.N_p = shape_matrices_p[ip].N;
69 ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
70
71 _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
72 }
73}
74
75template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
76 int DisplacementDim>
77std::tuple<
78 std::vector<ConstitutiveRelations::ConstitutiveData<DisplacementDim>>,
79 std::vector<ConstitutiveRelations::ConstitutiveTempData<DisplacementDim>>>
80TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
81 DisplacementDim>::
82 updateConstitutiveVariables(
83 Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_x_prev,
84 double const t, double const dt,
85 ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const&
86 models)
87{
88 [[maybe_unused]] auto const matrix_size =
89 gas_pressure_size + capillary_pressure_size + temperature_size +
90 displacement_size;
91
92 assert(local_x.size() == matrix_size);
93
94 auto const gas_pressure =
95 local_x.template segment<gas_pressure_size>(gas_pressure_index);
96 auto const gas_pressure_prev =
97 local_x_prev.template segment<gas_pressure_size>(gas_pressure_index);
98 auto const capillary_pressure =
99 local_x.template segment<capillary_pressure_size>(
100 capillary_pressure_index);
101 auto const capillary_pressure_prev =
102 local_x_prev.template segment<capillary_pressure_size>(
103 capillary_pressure_index);
104
105 auto const temperature =
106 local_x.template segment<temperature_size>(temperature_index);
107 auto const temperature_prev =
108 local_x_prev.template segment<temperature_size>(temperature_index);
109
110 auto const displacement =
111 local_x.template segment<displacement_size>(displacement_index);
112 auto const displacement_prev =
113 local_x_prev.template segment<displacement_size>(displacement_index);
114
115 auto const& medium =
116 *this->process_data_.media_map.getMedium(this->element_.getID());
117 ConstitutiveRelations::MediaData media_data{medium};
118
119 unsigned const n_integration_points =
120 this->integration_method_.getNumberOfPoints();
121
122 std::vector<ConstitutiveRelations::ConstitutiveData<DisplacementDim>>
123 ip_constitutive_data(n_integration_points);
124 std::vector<ConstitutiveRelations::ConstitutiveTempData<DisplacementDim>>
125 ip_constitutive_variables(n_integration_points);
126
127 for (unsigned ip = 0; ip < n_integration_points; ip++)
128 {
129 auto& ip_data = _ip_data[ip];
130 auto& ip_cv = ip_constitutive_variables[ip];
131 auto& ip_cd = ip_constitutive_data[ip];
132 auto& ip_out = this->output_data_[ip];
133 auto& current_state = this->current_states_[ip];
134 auto& prev_state = this->prev_states_[ip];
135
136 auto const& Np = ip_data.N_p;
137 auto const& NT = Np;
138 auto const& Nu = ip_data.N_u;
139 auto const& gradNu = ip_data.dNdx_u;
140 auto const& gradNp = ip_data.dNdx_p;
141 ParameterLib::SpatialPosition const pos{
142 std::nullopt, this->element_.getID(),
143 MathLib::Point3d(
144 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
145 ShapeMatricesTypeDisplacement>(
146 this->element_, Nu))};
147
148 auto const x_coord =
149 pos.getCoordinates().value()[0]; // r for axisymmetry
150
151 double const T = NT.dot(temperature);
152 double const T_prev = NT.dot(temperature_prev);
153 double const pG = Np.dot(gas_pressure);
154 double const pG_prev = Np.dot(gas_pressure_prev);
155 double const pCap = Np.dot(capillary_pressure);
156 double const pCap_prev = Np.dot(capillary_pressure_prev);
157 ConstitutiveRelations::TemperatureData const T_data{T, T_prev};
158 ConstitutiveRelations::GasPressureData const pGR_data{pG, pG_prev};
159 ConstitutiveRelations::CapillaryPressureData const pCap_data{pCap,
160 pCap_prev};
161 ConstitutiveRelations::ReferenceTemperatureData const T0{
162 this->process_data_.reference_temperature(t, pos)[0]};
163 ConstitutiveRelations::GasPressureGradientData<DisplacementDim> const
164 grad_p_GR{gradNp * gas_pressure};
165 ConstitutiveRelations::CapillaryPressureGradientData<
166 DisplacementDim> const grad_p_cap{gradNp * capillary_pressure};
167 ConstitutiveRelations::TemperatureGradientData<DisplacementDim> const
168 grad_T{gradNp * temperature};
169
170 // medium properties
171 models.elastic_tangent_stiffness_model.eval({pos, t, dt}, T_data,
172 ip_cv.C_el_data);
173
174 models.biot_model.eval({pos, t, dt}, media_data, ip_cv.biot_data);
175
176 auto const Bu =
177 LinearBMatrix::computeBMatrix<DisplacementDim,
178 ShapeFunctionDisplacement::NPOINTS,
179 typename BMatricesType::BMatrixType>(
180 gradNu, Nu, x_coord, this->is_axially_symmetric_);
181
182 ip_out.eps_data.eps.noalias() = Bu * displacement;
183 models.S_L_model.eval({pos, t, dt}, media_data, pCap_data,
184 current_state.S_L_data);
185
186 models.chi_S_L_model.eval({pos, t, dt}, media_data,
187 current_state.S_L_data,
188 current_state.chi_S_L);
189
190 models.chi_S_L_prev_model.eval({pos, t, dt}, media_data,
191 prev_state.S_L_data, prev_state.chi_S_L);
192
193 // solid phase compressibility
194 models.beta_p_SR_model.eval({pos, t, dt}, ip_cv.biot_data,
195 ip_cv.C_el_data, ip_cv.beta_p_SR);
196
197 // If there is swelling stress rate, compute swelling stress.
198 models.swelling_model.eval(
199 {pos, t, dt}, media_data, ip_cv.C_el_data, current_state.S_L_data,
200 prev_state.S_L_data, prev_state.swelling_data,
201 current_state.swelling_data, ip_cv.swelling_data);
202
203 // solid phase linear thermal expansion coefficient
204 models.s_therm_exp_model.eval({pos, t, dt}, media_data, T_data, T0,
205 ip_cv.s_therm_exp_data);
206
207 models.mechanical_strain_model.eval(
208 T_data, ip_cv.s_therm_exp_data, ip_out.eps_data,
209 Bu * displacement_prev, prev_state.mechanical_strain_data,
210 ip_cv.swelling_data, current_state.mechanical_strain_data);
211
212 models.s_mech_model.eval(
213 {pos, t, dt}, T_data, current_state.mechanical_strain_data,
214 prev_state.mechanical_strain_data, prev_state.eff_stress_data,
215 current_state.eff_stress_data, this->material_states_[ip],
216 ip_cd.s_mech_data, ip_cv.equivalent_plastic_strain_data);
217
218 models.total_stress_model.eval(current_state.eff_stress_data,
219 ip_cv.biot_data, current_state.chi_S_L,
220 pGR_data, pCap_data,
221 ip_out.total_stress_data);
222
223 models.pure_liquid_density_model.eval({pos, t, dt}, media_data,
224 pGR_data, pCap_data, T_data,
225 current_state.rho_W_LR);
226
227 models.phase_transition_model.eval(
228 {pos, t, dt}, media_data, pGR_data, pCap_data, T_data,
229 current_state.rho_W_LR, ip_out.fluid_enthalpy_data,
230 ip_out.mass_mole_fractions_data, ip_out.fluid_density_data,
231 ip_out.vapour_pressure_data, current_state.constituent_density_data,
232 ip_cv.phase_transition_data);
233
234 models.viscosity_model.eval({pos, t, dt}, media_data, T_data,
235 ip_out.mass_mole_fractions_data,
236 ip_cv.viscosity_data);
237
238 models.porosity_model.eval(
239 {pos, t, dt}, media_data, current_state.S_L_data,
240 prev_state.S_L_data, pCap_data, pGR_data, current_state.chi_S_L,
241 prev_state.chi_S_L, ip_cv.beta_p_SR, ip_out.eps_data,
242 Bu * displacement_prev, prev_state.porosity_data,
243 current_state.porosity_data);
244
245 if (medium.hasProperty(MPL::PropertyType::transport_porosity))
246 {
247 models.transport_porosity_model.eval(
248 {pos, t, dt}, media_data, current_state.S_L_data,
249 prev_state.S_L_data, pCap_data, pGR_data, current_state.chi_S_L,
250 prev_state.chi_S_L, ip_cv.beta_p_SR,
251 current_state.mechanical_strain_data,
252 prev_state.mechanical_strain_data,
253 prev_state.transport_porosity_data, current_state.porosity_data,
254 current_state.transport_porosity_data);
255 }
256 else
257 {
258 current_state.transport_porosity_data.phi =
259 current_state.porosity_data.phi;
260 }
261
262 models.permeability_model.eval(
263 {pos, t, dt}, media_data, current_state.S_L_data, pGR_data,
264 pCap_data, T_data, current_state.transport_porosity_data,
265 ip_out.total_stress_data, current_state.mechanical_strain_data,
266 ip_out.eps_data, ip_cv.equivalent_plastic_strain_data,
267 ip_out.permeability_data);
268
269 models.solid_density_model.eval(
270 {pos, t, dt}, media_data, T_data, current_state.eff_stress_data,
271 pCap_data, pGR_data, current_state.chi_S_L,
272 current_state.porosity_data, ip_out.solid_density_data);
273
274 models.solid_heat_capacity_model.eval({pos, t, dt}, media_data, T_data,
275 ip_cv.solid_heat_capacity_data);
276
277 models.thermal_conductivity_model.eval(
278 {pos, t, dt}, media_data, T_data, current_state.porosity_data,
279 current_state.S_L_data, ip_cv.thermal_conductivity_data);
280
281 models.advection_model.eval(current_state.constituent_density_data,
282 ip_out.permeability_data,
283 current_state.rho_W_LR,
284 ip_cv.viscosity_data,
285 ip_cv.advection_data);
286
287 models.gravity_model.eval(
288 ip_out.fluid_density_data,
289 current_state.porosity_data,
290 current_state.S_L_data,
291 ip_out.solid_density_data,
292 ConstitutiveRelations::SpecificBodyForce<DisplacementDim>(
293 this->process_data_.specific_body_force),
294 ip_cv.volumetric_body_force);
295
296 models.diffusion_velocity_model.eval(grad_p_cap,
297 grad_p_GR,
298 ip_out.mass_mole_fractions_data,
299 ip_cv.phase_transition_data,
300 current_state.porosity_data,
301 current_state.S_L_data,
302 grad_T,
303 ip_out.diffusion_velocity_data);
304
305 models.solid_enthalpy_model.eval(ip_cv.solid_heat_capacity_data, T_data,
306 ip_out.solid_enthalpy_data);
307
308 models.internal_energy_model.eval(ip_out.fluid_density_data,
309 ip_cv.phase_transition_data,
310 current_state.porosity_data,
311 current_state.S_L_data,
312 ip_out.solid_density_data,
313 ip_out.solid_enthalpy_data,
314 current_state.internal_energy_data);
315
316 models.effective_volumetric_enthalpy_model.eval(
317 ip_out.fluid_density_data,
318 ip_out.fluid_enthalpy_data,
319 current_state.porosity_data,
320 current_state.S_L_data,
321 ip_out.solid_density_data,
322 ip_out.solid_enthalpy_data,
323 ip_cv.effective_volumetric_enthalpy_data);
324
325 models.fC_1_model.eval(ip_cv.advection_data, ip_out.fluid_density_data,
326 ip_cv.fC_1);
327
328 if (!this->process_data_.apply_mass_lumping)
329 {
330 models.fC_2a_model.eval(ip_cv.biot_data,
331 pCap_data,
332 current_state.constituent_density_data,
333 current_state.porosity_data,
334 current_state.S_L_data,
335 ip_cv.beta_p_SR,
336 ip_cv.fC_2a);
337 }
338 models.fC_3a_model.eval(dt,
339 current_state.constituent_density_data,
340 prev_state.constituent_density_data,
341 current_state.S_L_data,
342 ip_cv.fC_3a);
343
344 models.fC_4_LCpG_model.eval(ip_cv.advection_data,
345 ip_out.fluid_density_data,
346 ip_cv.phase_transition_data,
347 current_state.porosity_data,
348 current_state.S_L_data,
349 ip_cv.fC_4_LCpG);
350
351 models.fC_4_LCpC_model.eval(ip_cv.advection_data,
352 ip_out.fluid_density_data,
353 ip_cv.phase_transition_data,
354 current_state.porosity_data,
355 current_state.S_L_data,
356 ip_cv.fC_4_LCpC);
357
358 models.fC_4_LCT_model.eval(ip_out.fluid_density_data,
359 ip_cv.phase_transition_data,
360 current_state.porosity_data,
361 current_state.S_L_data,
362 ip_cv.fC_4_LCT);
363
364 models.fC_4_MCpG_model.eval(ip_cv.biot_data,
365 current_state.constituent_density_data,
366 current_state.porosity_data,
367 current_state.S_L_data,
368 ip_cv.beta_p_SR,
369 ip_cv.fC_4_MCpG);
370
371 models.fC_4_MCpC_model.eval(ip_cv.biot_data,
372 pCap_data,
373 current_state.constituent_density_data,
374 current_state.porosity_data,
375 prev_state.S_L_data,
376 current_state.S_L_data,
377 ip_cv.beta_p_SR,
378 ip_cv.fC_4_MCpC);
379
380 models.fC_4_MCT_model.eval(ip_cv.biot_data,
381 current_state.constituent_density_data,
382 current_state.porosity_data,
383 current_state.S_L_data,
384 ip_cv.s_therm_exp_data,
385 ip_cv.fC_4_MCT);
386
387 models.fC_4_MCu_model.eval(ip_cv.biot_data,
388 current_state.constituent_density_data,
389 current_state.S_L_data,
390 ip_cv.fC_4_MCu);
391
392 models.fW_1_model.eval(ip_cv.advection_data, ip_out.fluid_density_data,
393 ip_cv.fW_1);
394
395 if (!this->process_data_.apply_mass_lumping)
396 {
397 models.fW_2_model.eval(ip_cv.biot_data,
398 pCap_data,
399 current_state.constituent_density_data,
400 current_state.porosity_data,
401 current_state.rho_W_LR,
402 current_state.S_L_data,
403 ip_cv.beta_p_SR,
404 ip_cv.fW_2);
405 }
406 models.fW_3a_model.eval(dt,
407 current_state.constituent_density_data,
408 prev_state.constituent_density_data,
409 prev_state.rho_W_LR,
410 current_state.rho_W_LR,
411 current_state.S_L_data,
412 ip_cv.fW_3a);
413
414 models.fW_4_LWpG_model.eval(ip_cv.advection_data,
415 ip_out.fluid_density_data,
416 ip_cv.phase_transition_data,
417 current_state.porosity_data,
418 current_state.S_L_data,
419 ip_cv.fW_4_LWpG);
420
421 models.fW_4_LWpC_model.eval(ip_cv.advection_data,
422 ip_out.fluid_density_data,
423 ip_cv.phase_transition_data,
424 current_state.porosity_data,
425 current_state.S_L_data,
426 ip_cv.fW_4_LWpC);
427
428 models.fW_4_LWT_model.eval(ip_out.fluid_density_data,
429 ip_cv.phase_transition_data,
430 current_state.porosity_data,
431 current_state.S_L_data,
432 ip_cv.fW_4_LWT);
433
434 models.fW_4_MWpG_model.eval(ip_cv.biot_data,
435 current_state.constituent_density_data,
436 current_state.porosity_data,
437 current_state.rho_W_LR,
438 current_state.S_L_data,
439 ip_cv.beta_p_SR,
440 ip_cv.fW_4_MWpG);
441
442 models.fW_4_MWpC_model.eval(ip_cv.biot_data,
443 pCap_data,
444 current_state.constituent_density_data,
445 current_state.porosity_data,
446 prev_state.S_L_data,
447 current_state.rho_W_LR,
448 current_state.S_L_data,
449 ip_cv.beta_p_SR,
450 ip_cv.fW_4_MWpC);
451
452 models.fW_4_MWT_model.eval(ip_cv.biot_data,
453 current_state.constituent_density_data,
454 current_state.porosity_data,
455 current_state.rho_W_LR,
456 current_state.S_L_data,
457 ip_cv.s_therm_exp_data,
458 ip_cv.fW_4_MWT);
459
460 models.fW_4_MWu_model.eval(ip_cv.biot_data,
461 current_state.constituent_density_data,
462 current_state.rho_W_LR,
463 current_state.S_L_data,
464 ip_cv.fW_4_MWu);
465
466 models.fT_1_model.eval(dt,
467 current_state.internal_energy_data,
468 prev_state.internal_energy_data,
469 ip_cv.fT_1);
470
471 // ---------------------------------------------------------------------
472 // Derivatives for Jacobian
473 // ---------------------------------------------------------------------
474
475 models.darcy_velocity_model.eval(
476 grad_p_cap,
477 ip_out.fluid_density_data,
478 grad_p_GR,
479 ip_out.permeability_data,
480 ConstitutiveRelations::SpecificBodyForce<DisplacementDim>(
481 this->process_data_.specific_body_force),
482 ip_cv.viscosity_data,
483 ip_out.darcy_velocity_data);
484
485 models.fT_2_model.eval(ip_out.darcy_velocity_data,
486 ip_out.fluid_density_data,
487 ip_out.fluid_enthalpy_data,
488 ip_cv.fT_2);
489
490 models.fT_3_model.eval(
491 current_state.constituent_density_data,
492 ip_out.darcy_velocity_data,
493 ip_out.diffusion_velocity_data,
494 ip_out.fluid_density_data,
495 ip_cv.phase_transition_data,
496 ConstitutiveRelations::SpecificBodyForce<DisplacementDim>(
497 this->process_data_.specific_body_force),
498 ip_cv.fT_3);
499
500 models.fu_2_KupC_model.eval(ip_cv.biot_data, current_state.chi_S_L,
501 ip_cv.fu_2_KupC);
502 }
503
504 return {ip_constitutive_data, ip_constitutive_variables};
505}
506
507template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
508 int DisplacementDim>
509std::vector<ConstitutiveRelations::DerivativesData<DisplacementDim>>
510TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
511 DisplacementDim>::
512 updateConstitutiveVariablesDerivatives(
513 Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_x_prev,
514 double const t, double const dt,
515 std::vector<
516 ConstitutiveRelations::ConstitutiveData<DisplacementDim>> const&
517 ip_constitutive_data,
518 std::vector<
519 ConstitutiveRelations::ConstitutiveTempData<DisplacementDim>> const&
520 ip_constitutive_variables,
521 ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const&
522 models)
523{
524 [[maybe_unused]] auto const matrix_size =
525 gas_pressure_size + capillary_pressure_size + temperature_size +
526 displacement_size;
527
528 assert(local_x.size() == matrix_size);
529
530 auto const gas_pressure =
531 local_x.template segment<gas_pressure_size>(gas_pressure_index);
532 auto const gas_pressure_prev =
533 local_x_prev.template segment<gas_pressure_size>(gas_pressure_index);
534 auto const temperature =
535 local_x.template segment<temperature_size>(temperature_index);
536 auto const temperature_prev =
537 local_x_prev.template segment<temperature_size>(temperature_index);
538 auto const displacement_prev =
539 local_x_prev.template segment<displacement_size>(displacement_index);
540
541 auto const capillary_pressure =
542 local_x.template segment<capillary_pressure_size>(
543 capillary_pressure_index);
544 auto const capillary_pressure_prev =
545 local_x_prev.template segment<capillary_pressure_size>(
546 capillary_pressure_index);
547
548 auto const& medium =
549 *this->process_data_.media_map.getMedium(this->element_.getID());
550 ConstitutiveRelations::MediaData media_data{medium};
551
552 unsigned const n_integration_points =
553 this->integration_method_.getNumberOfPoints();
554
555 std::vector<ConstitutiveRelations::DerivativesData<DisplacementDim>>
556 ip_d_data(n_integration_points);
557
558 for (unsigned ip = 0; ip < n_integration_points; ip++)
559 {
560 auto const& ip_data = _ip_data[ip];
561 auto& ip_dd = ip_d_data[ip];
562 auto const& ip_cd = ip_constitutive_data[ip];
563 auto const& ip_cv = ip_constitutive_variables[ip];
564 auto const& ip_out = this->output_data_[ip];
565 auto const& current_state = this->current_states_[ip];
566 auto const& prev_state = this->prev_states_[ip];
567
568 auto const& Nu = ip_data.N_u;
569 auto const& Np = ip_data.N_p;
570 auto const& NT = Np;
571 auto const& gradNu = ip_data.dNdx_u;
572
573 ParameterLib::SpatialPosition const pos{
574 std::nullopt, this->element_.getID(),
575 MathLib::Point3d(
576 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
577 ShapeMatricesTypeDisplacement>(
578 this->element_, Nu))};
579
580 auto const x_coord = pos.getCoordinates().value()[0];
581
582 double const T = NT.dot(temperature);
583 double const T_prev = NT.dot(temperature_prev);
584 double const pG = Np.dot(gas_pressure);
585 double const pG_prev = Np.dot(gas_pressure_prev);
586 double const pCap = Np.dot(capillary_pressure);
587 double const pCap_prev = Np.dot(capillary_pressure_prev);
588 ConstitutiveRelations::TemperatureData const T_data{T, T_prev};
589 ConstitutiveRelations::GasPressureData const pGR_data{pG, pG_prev};
590 ConstitutiveRelations::CapillaryPressureData const pCap_data{pCap,
591 pCap_prev};
592
593 auto const Bu =
594 LinearBMatrix::computeBMatrix<DisplacementDim,
595 ShapeFunctionDisplacement::NPOINTS,
596 typename BMatricesType::BMatrixType>(
597 gradNu, Nu, x_coord, this->is_axially_symmetric_);
598
599 models.S_L_model.dEval({pos, t, dt}, media_data, pCap_data,
600 ip_dd.dS_L_dp_cap);
601
602 models.advection_model.dEval(current_state.constituent_density_data,
603 ip_out.permeability_data,
604 ip_cv.viscosity_data,
605 ip_dd.dS_L_dp_cap,
606 ip_cv.phase_transition_data,
607 ip_dd.advection_d_data);
608
609 models.porosity_model.dEval(
610 {pos, t, dt}, media_data, current_state.S_L_data,
611 prev_state.S_L_data, pCap_data, pGR_data, current_state.chi_S_L,
612 prev_state.chi_S_L, ip_cv.beta_p_SR, ip_out.eps_data,
613 Bu * displacement_prev, prev_state.porosity_data,
614 ip_dd.porosity_d_data);
615
616 models.thermal_conductivity_model.dEval(
617 {pos, t, dt}, media_data, T_data, current_state.porosity_data,
618 ip_dd.porosity_d_data, current_state.S_L_data,
619 ip_dd.thermal_conductivity_d_data);
620
621 models.solid_density_model.dEval(
622 {pos, t, dt}, media_data, T_data, current_state.eff_stress_data,
623 pCap_data, pGR_data, current_state.chi_S_L,
624 current_state.porosity_data, ip_dd.solid_density_d_data);
625
626 models.internal_energy_model.dEval(
627 ip_out.fluid_density_data,
628 ip_cv.phase_transition_data,
629 current_state.porosity_data,
630 ip_dd.porosity_d_data,
631 current_state.S_L_data,
632 ip_out.solid_density_data,
633 ip_dd.solid_density_d_data,
634 ip_out.solid_enthalpy_data,
635 ip_cv.solid_heat_capacity_data,
636 ip_dd.effective_volumetric_internal_energy_d_data);
637
638 models.effective_volumetric_enthalpy_model.dEval(
639 ip_out.fluid_density_data,
640 ip_out.fluid_enthalpy_data,
641 ip_cv.phase_transition_data,
642 current_state.porosity_data,
643 ip_dd.porosity_d_data,
644 current_state.S_L_data,
645 ip_out.solid_density_data,
646 ip_dd.solid_density_d_data,
647 ip_out.solid_enthalpy_data,
648 ip_cv.solid_heat_capacity_data,
649 ip_dd.effective_volumetric_enthalpy_d_data);
650 if (!this->process_data_.apply_mass_lumping)
651 {
652 models.fC_2a_model.dEval(ip_cv.biot_data,
653 pCap_data,
654 current_state.constituent_density_data,
655 ip_cv.phase_transition_data,
656 current_state.porosity_data,
657 ip_dd.porosity_d_data,
658 current_state.S_L_data,
659 ip_dd.dS_L_dp_cap,
660 ip_cv.beta_p_SR,
661 ip_dd.dfC_2a);
662 }
663 models.fC_3a_model.dEval(dt,
664 current_state.constituent_density_data,
665 prev_state.constituent_density_data,
666 ip_cv.phase_transition_data,
667 current_state.S_L_data,
668 ip_dd.dS_L_dp_cap,
669 ip_dd.dfC_3a);
670
671 models.fC_4_LCpG_model.dEval(ip_out.permeability_data,
672 ip_cv.viscosity_data,
673 ip_cv.phase_transition_data,
674 ip_dd.advection_d_data,
675 ip_dd.dfC_4_LCpG);
676
677 models.fC_4_LCpC_model.dEval(current_state.constituent_density_data,
678 ip_out.permeability_data,
679 ip_cv.phase_transition_data,
680 ip_dd.dS_L_dp_cap,
681 ip_cv.viscosity_data,
682 ip_dd.dfC_4_LCpC);
683
684 models.fC_4_MCpG_model.dEval(ip_cv.biot_data,
685 current_state.constituent_density_data,
686 ip_cv.phase_transition_data,
687 current_state.porosity_data,
688 ip_dd.porosity_d_data,
689 current_state.S_L_data,
690 ip_cv.beta_p_SR,
691 ip_dd.dfC_4_MCpG);
692
693 models.fC_4_MCT_model.dEval(ip_cv.biot_data,
694 current_state.constituent_density_data,
695 ip_cv.phase_transition_data,
696 current_state.porosity_data,
697 ip_dd.porosity_d_data,
698 current_state.S_L_data,
699 ip_cv.s_therm_exp_data,
700 ip_dd.dfC_4_MCT);
701
702 models.fC_4_MCu_model.dEval(ip_cv.biot_data,
703 ip_cv.phase_transition_data,
704 current_state.S_L_data,
705 ip_dd.dfC_4_MCu);
706
707 if (!this->process_data_.apply_mass_lumping)
708 {
709 models.fW_2_model.dEval(ip_cv.biot_data,
710 pCap_data,
711 current_state.constituent_density_data,
712 ip_cv.phase_transition_data,
713 current_state.porosity_data,
714 ip_dd.porosity_d_data,
715 current_state.rho_W_LR,
716 current_state.S_L_data,
717 ip_dd.dS_L_dp_cap,
718 ip_cv.beta_p_SR,
719 ip_dd.dfW_2);
720 }
721
722 models.fW_3a_model.dEval(dt,
723 current_state.constituent_density_data,
724 ip_cv.phase_transition_data,
725 prev_state.constituent_density_data,
726 prev_state.rho_W_LR,
727 current_state.rho_W_LR,
728 current_state.S_L_data,
729 ip_dd.dS_L_dp_cap,
730 ip_dd.dfW_3a);
731
732 models.fW_4_LWpG_model.dEval(current_state.constituent_density_data,
733 ip_out.permeability_data,
734 ip_cv.phase_transition_data,
735 current_state.rho_W_LR,
736 ip_dd.dS_L_dp_cap,
737 ip_cv.viscosity_data,
738 ip_dd.dfW_4_LWpG);
739
740 models.fW_4_LWpC_model.dEval(ip_cv.advection_data,
741 ip_out.fluid_density_data,
742 ip_out.permeability_data,
743 ip_cv.phase_transition_data,
744 current_state.porosity_data,
745 current_state.rho_W_LR,
746 current_state.S_L_data,
747 ip_dd.dS_L_dp_cap,
748 ip_cv.viscosity_data,
749 ip_dd.dfW_4_LWpC);
750
751 models.fT_1_model.dEval(
752 dt, ip_dd.effective_volumetric_internal_energy_d_data, ip_dd.dfT_1);
753
754 models.fT_2_model.dEval(
755 ip_out.darcy_velocity_data,
756 ip_out.fluid_density_data,
757 ip_out.fluid_enthalpy_data,
758 ip_out.permeability_data,
759 ip_cv.phase_transition_data,
760 ConstitutiveRelations::SpecificBodyForce<DisplacementDim>(
761 this->process_data_.specific_body_force),
762 ip_cv.viscosity_data,
763 ip_dd.dfT_2);
764
765 models.fu_1_KuT_model.dEval(ip_cd.s_mech_data, ip_cv.s_therm_exp_data,
766 ip_dd.dfu_1_KuT);
767
768 models.fu_2_KupC_model.dEval(ip_cv.biot_data,
769 current_state.chi_S_L,
770 pCap_data,
771 ip_dd.dS_L_dp_cap,
772 ip_dd.dfu_2_KupC);
773 }
774
775 return ip_d_data;
776}
777
778template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
779 int DisplacementDim>
780std::size_t TH2MLocalAssembler<
781 ShapeFunctionDisplacement, ShapeFunctionPressure,
782 DisplacementDim>::setIPDataInitialConditions(std::string_view name,
783 double const* values,
784 int const integration_order)
785{
786 if (integration_order !=
787 static_cast<int>(this->integration_method_.getIntegrationOrder()))
788 {
789 OGS_FATAL(
790 "Setting integration point initial conditions; The integration "
791 "order of the local assembler for element {:d} is different "
792 "from the integration order in the initial condition.",
793 this->element_.getID());
794 }
795
796 if (name == "sigma" && this->process_data_.initial_stress.value)
797 {
798 OGS_FATAL(
799 "Setting initial conditions for stress from integration "
800 "point data and from a parameter '{:s}' is not possible "
801 "simultaneously.",
802 this->process_data_.initial_stress.value->name);
803 }
804
805 if (name.starts_with("material_state_variable_"))
806 {
807 name.remove_prefix(24);
808 DBUG("Setting material state variable '{:s}'", name);
809
810 auto const& internal_variables =
811 this->solid_material_.getInternalVariables();
812 if (auto const iv = std::find_if(
813 begin(internal_variables), end(internal_variables),
814 [&name](auto const& iv) { return iv.name == name; });
815 iv != end(internal_variables))
816 {
817 DBUG("Setting material state variable '{:s}'", name);
818 return ProcessLib::setIntegrationPointDataMaterialStateVariables(
819 values, this->material_states_,
820 &ConstitutiveRelations::MaterialStateData<
821 DisplacementDim>::material_state_variables,
822 iv->reference);
823 }
824
825 WARN(
826 "Could not find variable {:s} in solid material model's "
827 "internal variables.",
828 name);
829 return 0;
830 }
831
832 // TODO this logic could be pulled out of the local assembler into the
833 // process. That might lead to a slightly better performance due to less
834 // string comparisons.
835 return ProcessLib::Reflection::reflectSetIPData<DisplacementDim>(
836 name, values, this->current_states_);
837}
838
839template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
840 int DisplacementDim>
841void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
842 DisplacementDim>::
843 setInitialConditionsConcrete(Eigen::VectorXd const local_x,
844 double const t,
845 int const /*process_id*/)
846{
847 [[maybe_unused]] auto const matrix_size =
848 gas_pressure_size + capillary_pressure_size + temperature_size +
849 displacement_size;
850
851 assert(local_x.size() == matrix_size);
852
853 auto const capillary_pressure =
854 local_x.template segment<capillary_pressure_size>(
855 capillary_pressure_index);
856
857 auto const p_GR =
858 local_x.template segment<gas_pressure_size>(gas_pressure_index);
859
860 auto const temperature =
861 local_x.template segment<temperature_size>(temperature_index);
862
863 auto const displacement =
864 local_x.template segment<displacement_size>(displacement_index);
865
866 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
867 auto const& medium =
868 *this->process_data_.media_map.getMedium(this->element_.getID());
869 auto const& solid_phase =
870 medium.phase(MaterialPropertyLib::PhaseName::Solid);
871
872 ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const models{
873 this->solid_material_, *this->process_data_.phase_transition_model_};
874
875 unsigned const n_integration_points =
876 this->integration_method_.getNumberOfPoints();
877
878 for (unsigned ip = 0; ip < n_integration_points; ip++)
879 {
880 MPL::VariableArray vars;
881
882 auto& ip_data = _ip_data[ip];
883 auto& ip_out = this->output_data_[ip];
884 auto& prev_state = this->prev_states_[ip];
885 auto const& Np = ip_data.N_p;
886 auto const& NT = Np;
887 auto const& Nu = ip_data.N_u;
888 auto const& gradNu = ip_data.dNdx_u;
889
890 ParameterLib::SpatialPosition const pos{
891 std::nullopt, this->element_.getID(),
892 MathLib::Point3d(
893 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
894 ShapeMatricesTypeDisplacement>(
895 this->element_, ip_data.N_u))};
896
897 auto const x_coord = pos.getCoordinates().value()[0];
898
899 double const pCap = Np.dot(capillary_pressure);
900 vars.capillary_pressure = pCap;
901
902 double const T = NT.dot(temperature);
903 ConstitutiveRelations::TemperatureData const T_data{
904 T, T}; // T_prev = T in initialization.
905 vars.temperature = T;
906
907 auto const Bu =
908 LinearBMatrix::computeBMatrix<DisplacementDim,
909 ShapeFunctionDisplacement::NPOINTS,
910 typename BMatricesType::BMatrixType>(
911 gradNu, Nu, x_coord, this->is_axially_symmetric_);
912
913 auto& eps = ip_out.eps_data.eps;
914 eps.noalias() = Bu * displacement;
915
916 // Set volumetric strain rate for the general case without swelling.
917 vars.volumetric_strain = Invariants::trace(eps);
918
919 double const S_L =
920 medium.property(MPL::PropertyType::saturation)
921 .template value<double>(
922 vars, pos, t, std::numeric_limits<double>::quiet_NaN());
923 this->prev_states_[ip].S_L_data->S_L = S_L;
924
925 // TODO (naumov) Double computation of C_el might be avoided if
926 // updateConstitutiveVariables is called before. But it might interfere
927 // with eps_m initialization.
928 ConstitutiveRelations::ElasticTangentStiffnessData<DisplacementDim>
929 C_el_data;
930 models.elastic_tangent_stiffness_model.eval({pos, t, dt}, T_data,
931 C_el_data);
932 auto const& C_el = C_el_data.stiffness_tensor;
933
934 // Set eps_m_prev from potentially non-zero eps and sigma_sw from
935 // restart.
936 auto const& sigma_sw = this->current_states_[ip].swelling_data.sigma_sw;
937 prev_state.mechanical_strain_data->eps_m.noalias() =
938 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
939 ? eps + C_el.inverse() * sigma_sw
940 : eps;
941
942 if (this->process_data_.initial_stress.isTotalStress())
943 {
944 auto const alpha_b =
945 medium.property(MPL::PropertyType::biot_coefficient)
946 .template value<double>(vars, pos, t, 0.0 /*dt*/);
947
948 vars.liquid_saturation = S_L;
949 double const bishop =
950 medium.property(MPL::PropertyType::bishops_effective_stress)
951 .template value<double>(vars, pos, t, 0.0 /*dt*/);
952
953 this->current_states_[ip].eff_stress_data.sigma_eff.noalias() +=
954 alpha_b * Np.dot(p_GR - bishop * capillary_pressure) *
955 Invariants::identity2;
956 this->prev_states_[ip].eff_stress_data =
957 this->current_states_[ip].eff_stress_data;
958 }
959 }
960
961 // local_x_prev equal to local_x s.t. the local_x_dot is zero.
962 updateConstitutiveVariables(local_x, local_x, t, 0, models);
963
964 for (unsigned ip = 0; ip < n_integration_points; ip++)
965 {
966 this->material_states_[ip].pushBackState();
967 this->prev_states_[ip] = this->current_states_[ip];
968 }
969}
970
971template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
972 int DisplacementDim>
973void TH2MLocalAssembler<
974 ShapeFunctionDisplacement, ShapeFunctionPressure,
975 DisplacementDim>::assemble(double const t, double const dt,
976 std::vector<double> const& local_x,
977 std::vector<double> const& local_x_prev,
978 std::vector<double>& local_M_data,
979 std::vector<double>& local_K_data,
980 std::vector<double>& local_rhs_data)
981{
982 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
983 temperature_size + displacement_size;
984 assert(local_x.size() == matrix_size);
985
986 auto const capillary_pressure =
987 Eigen::Map<VectorType<capillary_pressure_size> const>(
988 local_x.data() + capillary_pressure_index, capillary_pressure_size);
989
990 auto const capillary_pressure_prev =
991 Eigen::Map<VectorType<capillary_pressure_size> const>(
992 local_x_prev.data() + capillary_pressure_index,
993 capillary_pressure_size);
994
995 // pointer to local_M_data vector
996 auto local_M =
997 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
998 local_M_data, matrix_size, matrix_size);
999
1000 // pointer to local_K_data vector
1001 auto local_K =
1002 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1003 local_K_data, matrix_size, matrix_size);
1004
1005 // pointer to local_rhs_data vector
1006 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
1007 local_rhs_data, matrix_size);
1008
1009 // component-formulation
1010 // W - liquid phase main component
1011 // C - gas phase main component
1012 // pointer-matrices to the mass matrix - C component equation
1013 auto MCpG = local_M.template block<C_size, gas_pressure_size>(
1014 C_index, gas_pressure_index);
1015 auto MCpC = local_M.template block<C_size, capillary_pressure_size>(
1016 C_index, capillary_pressure_index);
1017 auto MCT = local_M.template block<C_size, temperature_size>(
1018 C_index, temperature_index);
1019 auto MCu = local_M.template block<C_size, displacement_size>(
1020 C_index, displacement_index);
1021
1022 // pointer-matrices to the stiffness matrix - C component equation
1023 auto LCpG = local_K.template block<C_size, gas_pressure_size>(
1024 C_index, gas_pressure_index);
1025 auto LCpC = local_K.template block<C_size, capillary_pressure_size>(
1026 C_index, capillary_pressure_index);
1027 auto LCT = local_K.template block<C_size, temperature_size>(
1028 C_index, temperature_index);
1029
1030 // pointer-matrices to the mass matrix - W component equation
1031 auto MWpG = local_M.template block<W_size, gas_pressure_size>(
1032 W_index, gas_pressure_index);
1033 auto MWpC = local_M.template block<W_size, capillary_pressure_size>(
1034 W_index, capillary_pressure_index);
1035 auto MWT = local_M.template block<W_size, temperature_size>(
1036 W_index, temperature_index);
1037 auto MWu = local_M.template block<W_size, displacement_size>(
1038 W_index, displacement_index);
1039
1040 // pointer-matrices to the stiffness matrix - W component equation
1041 auto LWpG = local_K.template block<W_size, gas_pressure_size>(
1042 W_index, gas_pressure_index);
1043 auto LWpC = local_K.template block<W_size, capillary_pressure_size>(
1044 W_index, capillary_pressure_index);
1045 auto LWT = local_K.template block<W_size, temperature_size>(
1046 W_index, temperature_index);
1047
1048 // pointer-matrices to the mass matrix - temperature equation
1049 auto MTu = local_M.template block<temperature_size, displacement_size>(
1050 temperature_index, displacement_index);
1051
1052 // pointer-matrices to the stiffness matrix - temperature equation
1053 auto KTT = local_K.template block<temperature_size, temperature_size>(
1054 temperature_index, temperature_index);
1055
1056 // pointer-matrices to the stiffness matrix - displacement equation
1057 auto KUpG = local_K.template block<displacement_size, gas_pressure_size>(
1058 displacement_index, gas_pressure_index);
1059 auto KUpC =
1060 local_K.template block<displacement_size, capillary_pressure_size>(
1061 displacement_index, capillary_pressure_index);
1062
1063 auto KUU = local_K.template block<displacement_size, displacement_size>(
1064 displacement_index, displacement_index);
1065
1066 // pointer-vectors to the right hand side terms - C-component equation
1067 auto fC = local_f.template segment<C_size>(C_index);
1068 // pointer-vectors to the right hand side terms - W-component equation
1069 auto fW = local_f.template segment<W_size>(W_index);
1070 // pointer-vectors to the right hand side terms - temperature equation
1071 auto fT = local_f.template segment<temperature_size>(temperature_index);
1072 // pointer-vectors to the right hand side terms - displacement equation
1073 auto fU = local_f.template segment<displacement_size>(displacement_index);
1074
1075 unsigned const n_integration_points =
1076 this->integration_method_.getNumberOfPoints();
1077
1078 ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const models{
1079 this->solid_material_, *this->process_data_.phase_transition_model_};
1080
1081 auto const [ip_constitutive_data, ip_constitutive_variables] =
1082 updateConstitutiveVariables(
1083 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1084 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1085 local_x_prev.size()),
1086 t, dt, models);
1087
1088 for (unsigned int_point = 0; int_point < n_integration_points; int_point++)
1089 {
1090 auto& ip = _ip_data[int_point];
1091 auto& ip_cv = ip_constitutive_variables[int_point];
1092 auto& ip_cd = ip_constitutive_data[int_point];
1093
1094 auto& current_state = this->current_states_[int_point];
1095 auto const& prev_state = this->prev_states_[int_point];
1096
1097 auto const& Np = ip.N_p;
1098 auto const& Nu = ip.N_u;
1099 ParameterLib::SpatialPosition const pos{
1100 std::nullopt, this->element_.getID(),
1101 MathLib::Point3d(
1102 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
1103 ShapeMatricesTypeDisplacement>(
1104 this->element_, Nu))};
1105
1106 auto const& NpT = Np.transpose().eval();
1107 auto const& NTT = NpT;
1108
1109 auto const& gradNp = ip.dNdx_p;
1110 auto const& gradNT = gradNp;
1111 auto const& gradNu = ip.dNdx_u;
1112
1113 auto const& gradNpT = gradNp.transpose().eval();
1114 auto const& gradNTT = gradNpT;
1115
1116 auto const& w = ip.integration_weight;
1117
1118 auto const x_coord =
1119 pos.getCoordinates().value()[0]; // r for axisymetric
1120 auto const Bu =
1121 LinearBMatrix::computeBMatrix<DisplacementDim,
1122 ShapeFunctionDisplacement::NPOINTS,
1123 typename BMatricesType::BMatrixType>(
1124 gradNu, Nu, x_coord, this->is_axially_symmetric_);
1125
1126 auto const NTN = (Np.transpose() * Np).eval();
1127 auto const BTI2N = (Bu.transpose() * Invariants::identity2 * Np).eval();
1128
1129 double const pCap = Np.dot(capillary_pressure);
1130 double const pCap_prev = Np.dot(capillary_pressure_prev);
1131
1132 auto const s_L = current_state.S_L_data.S_L;
1133 auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
1134
1135 auto const& b = this->process_data_.specific_body_force;
1136
1137 // ---------------------------------------------------------------------
1138 // C-component equation
1139 // ---------------------------------------------------------------------
1140
1141 MCpG.noalias() += NTN * (ip_cv.fC_4_MCpG.m * w);
1142 MCpC.noalias() += NTN * (ip_cv.fC_4_MCpC.m * w);
1143
1144 if (this->process_data_.apply_mass_lumping)
1145 {
1146 if (pCap - pCap_prev != 0.) // avoid division by Zero
1147 {
1148 MCpC.noalias() +=
1149 NTN * (ip_cv.fC_4_MCpC.ml / (pCap - pCap_prev) * w);
1150 }
1151 }
1152
1153 MCT.noalias() += NTN * (ip_cv.fC_4_MCT.m * w);
1154 MCu.noalias() += BTI2N.transpose() * (ip_cv.fC_4_MCu.m * w);
1155
1156 LCpG.noalias() += gradNpT * ip_cv.fC_4_LCpG.L * gradNp * w;
1157
1158 LCpC.noalias() += gradNpT * ip_cv.fC_4_LCpC.L * gradNp * w;
1159
1160 LCT.noalias() += gradNpT * ip_cv.fC_4_LCT.L * gradNp * w;
1161
1162 fC.noalias() += gradNpT * ip_cv.fC_1.A * b * w;
1163
1164 if (!this->process_data_.apply_mass_lumping)
1165 {
1166 fC.noalias() -= NpT * (ip_cv.fC_2a.a * s_L_dot * w);
1167 }
1168 // fC_III
1169 fC.noalias() -=
1170 NpT * (current_state.porosity_data.phi * ip_cv.fC_3a.a * w);
1171
1172 // ---------------------------------------------------------------------
1173 // W-component equation
1174 // ---------------------------------------------------------------------
1175
1176 MWpG.noalias() += NTN * (ip_cv.fW_4_MWpG.m * w);
1177 MWpC.noalias() += NTN * (ip_cv.fW_4_MWpC.m * w);
1178
1179 if (this->process_data_.apply_mass_lumping)
1180 {
1181 if (pCap - pCap_prev != 0.) // avoid division by Zero
1182 {
1183 MWpC.noalias() +=
1184 NTN * (ip_cv.fW_4_MWpC.ml / (pCap - pCap_prev) * w);
1185 }
1186 }
1187
1188 MWT.noalias() += NTN * (ip_cv.fW_4_MWT.m * w);
1189
1190 MWu.noalias() += BTI2N.transpose() * (ip_cv.fW_4_MWu.m * w);
1191
1192 LWpG.noalias() += gradNpT * ip_cv.fW_4_LWpG.L * gradNp * w;
1193
1194 LWpC.noalias() += gradNpT * ip_cv.fW_4_LWpC.L * gradNp * w;
1195
1196 LWT.noalias() += gradNpT * ip_cv.fW_4_LWT.L * gradNp * w;
1197
1198 fW.noalias() += gradNpT * ip_cv.fW_1.A * b * w;
1199
1200 if (!this->process_data_.apply_mass_lumping)
1201 {
1202 fW.noalias() -= NpT * (ip_cv.fW_2.a * s_L_dot * w);
1203 }
1204
1205 fW.noalias() -=
1206 NpT * (current_state.porosity_data.phi * ip_cv.fW_3a.a * w);
1207
1208 // ---------------------------------------------------------------------
1209 // - temperature equation
1210 // ---------------------------------------------------------------------
1211
1212 MTu.noalias() +=
1213 BTI2N.transpose() *
1214 (ip_cv.effective_volumetric_enthalpy_data.rho_h_eff * w);
1215
1216 KTT.noalias() +=
1217 gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w;
1218
1219 fT.noalias() -= NTT * (ip_cv.fT_1.m * w);
1220
1221 fT.noalias() += gradNTT * ip_cv.fT_2.A * w;
1222
1223 fT.noalias() += gradNTT * ip_cv.fT_3.gradN * w;
1224
1225 fT.noalias() += NTT * (ip_cv.fT_3.N * w);
1226
1227 // ---------------------------------------------------------------------
1228 // - displacement equation
1229 // ---------------------------------------------------------------------
1230
1231 KUpG.noalias() -= BTI2N * (ip_cv.biot_data() * w);
1232
1233 KUpC.noalias() += BTI2N * (ip_cv.fu_2_KupC.m * w);
1234
1235 auto const BuTC =
1236 (Bu.transpose() * ip_cd.s_mech_data.stiffness_tensor).eval();
1237 KUU.noalias() += BuTC * Bu * w;
1238
1239 fU.noalias() -=
1240 (Bu.transpose() * current_state.eff_stress_data.sigma_eff -
1241 N_u_op(Nu).transpose() * ip_cv.volumetric_body_force()) *
1242 w;
1243
1244 if (this->process_data_.apply_mass_lumping)
1245 {
1246 MCpG = MCpG.colwise().sum().eval().asDiagonal();
1247 MCpC = MCpC.colwise().sum().eval().asDiagonal();
1248 MWpG = MWpG.colwise().sum().eval().asDiagonal();
1249 MWpC = MWpC.colwise().sum().eval().asDiagonal();
1250 }
1251 } // int_point-loop
1252}
1253
1254// Assembles the local Jacobian matrix. So far, the linearisation of HT part is
1255// not considered as that in HT process.
1256template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1257 int DisplacementDim>
1258void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
1259 DisplacementDim>::
1260 assembleWithJacobian(double const t, double const dt,
1261 std::vector<double> const& local_x,
1262 std::vector<double> const& local_x_prev,
1263 std::vector<double>& local_rhs_data,
1264 std::vector<double>& local_Jac_data)
1265{
1266 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
1267 temperature_size + displacement_size;
1268 assert(local_x.size() == matrix_size);
1269
1270 auto const temperature = Eigen::Map<VectorType<temperature_size> const>(
1271 local_x.data() + temperature_index, temperature_size);
1272
1273 auto const gas_pressure = Eigen::Map<VectorType<gas_pressure_size> const>(
1274 local_x.data() + gas_pressure_index, gas_pressure_size);
1275
1276 auto const capillary_pressure =
1277 Eigen::Map<VectorType<capillary_pressure_size> const>(
1278 local_x.data() + capillary_pressure_index, capillary_pressure_size);
1279
1280 auto const displacement = Eigen::Map<VectorType<displacement_size> const>(
1281 local_x.data() + displacement_index, displacement_size);
1282
1283 auto const gas_pressure_prev =
1284 Eigen::Map<VectorType<gas_pressure_size> const>(
1285 local_x_prev.data() + gas_pressure_index, gas_pressure_size);
1286
1287 auto const capillary_pressure_prev =
1288 Eigen::Map<VectorType<capillary_pressure_size> const>(
1289 local_x_prev.data() + capillary_pressure_index,
1290 capillary_pressure_size);
1291
1292 auto const temperature_prev =
1293 Eigen::Map<VectorType<temperature_size> const>(
1294 local_x_prev.data() + temperature_index, temperature_size);
1295
1296 auto const displacement_prev =
1297 Eigen::Map<VectorType<displacement_size> const>(
1298 local_x_prev.data() + displacement_index, displacement_size);
1299
1300 auto local_Jac =
1301 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1302 local_Jac_data, matrix_size, matrix_size);
1303
1304 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
1305 local_rhs_data, matrix_size);
1306
1307 // component-formulation
1308 // W - liquid phase main component
1309 // C - gas phase main component
1310
1311 // C component equation matrices
1312 MatrixType<C_size, gas_pressure_size> MCpG =
1313 MatrixType<C_size, gas_pressure_size>::Zero(C_size, gas_pressure_size);
1314 MatrixType<C_size, capillary_pressure_size> MCpC =
1315 MatrixType<C_size, capillary_pressure_size>::Zero(
1316 C_size, capillary_pressure_size);
1317 MatrixType<C_size, temperature_size> MCT =
1318 MatrixType<C_size, temperature_size>::Zero(C_size, temperature_size);
1319 MatrixType<C_size, displacement_size> MCu =
1320 MatrixType<C_size, displacement_size>::Zero(C_size, displacement_size);
1321
1322 MatrixType<C_size, gas_pressure_size> LCpG =
1323 MatrixType<C_size, gas_pressure_size>::Zero(C_size, gas_pressure_size);
1324 MatrixType<C_size, capillary_pressure_size> LCpC =
1325 MatrixType<C_size, capillary_pressure_size>::Zero(
1326 C_size, capillary_pressure_size);
1327 MatrixType<C_size, temperature_size> LCT =
1328 MatrixType<C_size, temperature_size>::Zero(C_size, temperature_size);
1329
1330 // mass matrix - W component equation
1331 MatrixType<W_size, gas_pressure_size> MWpG =
1332 MatrixType<W_size, gas_pressure_size>::Zero(W_size, gas_pressure_size);
1333 MatrixType<W_size, capillary_pressure_size> MWpC =
1334 MatrixType<W_size, capillary_pressure_size>::Zero(
1335 W_size, capillary_pressure_size);
1336 MatrixType<W_size, temperature_size> MWT =
1337 MatrixType<W_size, temperature_size>::Zero(W_size, temperature_size);
1338 MatrixType<W_size, displacement_size> MWu =
1339 MatrixType<W_size, displacement_size>::Zero(W_size, displacement_size);
1340
1341 // stiffness matrix - W component equation
1342 MatrixType<W_size, gas_pressure_size> LWpG =
1343 MatrixType<W_size, gas_pressure_size>::Zero(W_size, gas_pressure_size);
1344 MatrixType<W_size, capillary_pressure_size> LWpC =
1345 MatrixType<W_size, capillary_pressure_size>::Zero(
1346 W_size, capillary_pressure_size);
1347 MatrixType<W_size, temperature_size> LWT =
1348 MatrixType<W_size, temperature_size>::Zero(W_size, temperature_size);
1349
1350 // mass matrix - temperature equation
1351 MatrixType<temperature_size, displacement_size> MTu =
1352 MatrixType<temperature_size, displacement_size>::Zero(
1353 temperature_size, displacement_size);
1354
1355 // stiffness matrix - temperature equation
1356 MatrixType<temperature_size, temperature_size> KTT =
1357 MatrixType<temperature_size, temperature_size>::Zero(temperature_size,
1358 temperature_size);
1359
1360 // stiffness matrices - displacement equation coupling into pressures
1361 MatrixType<displacement_size, gas_pressure_size> KUpG =
1362 MatrixType<displacement_size, gas_pressure_size>::Zero(
1363 displacement_size, gas_pressure_size);
1364 MatrixType<displacement_size, capillary_pressure_size> KUpC =
1365 MatrixType<displacement_size, capillary_pressure_size>::Zero(
1366 displacement_size, capillary_pressure_size);
1367
1368 // pointer-vectors to the right hand side terms - C-component equation
1369 auto fC = local_f.template segment<C_size>(C_index);
1370 // pointer-vectors to the right hand side terms - W-component equation
1371 auto fW = local_f.template segment<W_size>(W_index);
1372 // pointer-vectors to the right hand side terms - temperature equation
1373 auto fT = local_f.template segment<temperature_size>(temperature_index);
1374 // pointer-vectors to the right hand side terms - displacement equation
1375 auto fU = local_f.template segment<displacement_size>(displacement_index);
1376
1377 unsigned const n_integration_points =
1378 this->integration_method_.getNumberOfPoints();
1379
1380 ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const models{
1381 this->solid_material_, *this->process_data_.phase_transition_model_};
1382
1383 auto const [ip_constitutive_data, ip_constitutive_variables] =
1384 updateConstitutiveVariables(
1385 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1386 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1387 local_x_prev.size()),
1388 t, dt, models);
1389
1390 auto const ip_d_data = updateConstitutiveVariablesDerivatives(
1391 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1392 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1393 local_x_prev.size()),
1394 t, dt, ip_constitutive_data, ip_constitutive_variables, models);
1395
1396 for (unsigned int_point = 0; int_point < n_integration_points; int_point++)
1397 {
1398 auto& ip = _ip_data[int_point];
1399 auto& ip_cd = ip_constitutive_data[int_point];
1400 auto& ip_dd = ip_d_data[int_point];
1401 auto& ip_cv = ip_constitutive_variables[int_point];
1402 auto& current_state = this->current_states_[int_point];
1403 auto& prev_state = this->prev_states_[int_point];
1404
1405 auto const& Np = ip.N_p;
1406 auto const& NT = Np;
1407 auto const& Nu = ip.N_u;
1408 ParameterLib::SpatialPosition const pos{
1409 std::nullopt, this->element_.getID(),
1410 MathLib::Point3d(
1411 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
1412 ShapeMatricesTypeDisplacement>(
1413 this->element_, Nu))};
1414
1415 auto const& NpT = Np.transpose().eval();
1416 auto const& NTT = NpT;
1417
1418 auto const& gradNp = ip.dNdx_p;
1419 auto const& gradNT = gradNp;
1420 auto const& gradNu = ip.dNdx_u;
1421
1422 auto const& gradNpT = gradNp.transpose().eval();
1423 auto const& gradNTT = gradNpT;
1424
1425 auto const& w = ip.integration_weight;
1426
1427 auto const x_coord =
1428 pos.getCoordinates().value()[0]; // r for axisymetric
1429 auto const Bu =
1430 LinearBMatrix::computeBMatrix<DisplacementDim,
1431 ShapeFunctionDisplacement::NPOINTS,
1432 typename BMatricesType::BMatrixType>(
1433 gradNu, Nu, x_coord, this->is_axially_symmetric_);
1434
1435 auto const NTN = (Np.transpose() * Np).eval();
1436 auto const BTI2N = (Bu.transpose() * Invariants::identity2 * Np).eval();
1437
1438 double const div_u_dot =
1439 Invariants::trace(Bu * (displacement - displacement_prev) / dt);
1440
1441 double const pGR = Np.dot(gas_pressure);
1442 double const pCap = Np.dot(capillary_pressure);
1443 double const T = NT.dot(temperature);
1444
1445 GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
1446 GlobalDimVectorType const gradpCap = gradNp * capillary_pressure;
1447 GlobalDimVectorType const gradT = gradNT * temperature;
1448
1449 double const pGR_prev = Np.dot(gas_pressure_prev);
1450 double const pCap_prev = Np.dot(capillary_pressure_prev);
1451 double const T_prev = NT.dot(temperature_prev);
1452
1453 auto const& s_L = current_state.S_L_data.S_L;
1454 auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
1455
1456 auto const& b = this->process_data_.specific_body_force;
1457
1458 // ---------------------------------------------------------------------
1459 // C-component equation
1460 // ---------------------------------------------------------------------
1461
1462 MCpG.noalias() += NTN * (ip_cv.fC_4_MCpG.m * w);
1463 MCpC.noalias() += NTN * (ip_cv.fC_4_MCpC.m * w);
1464
1465 if (this->process_data_.apply_mass_lumping)
1466 {
1467 if (pCap - pCap_prev != 0.) // avoid division by Zero
1468 {
1469 MCpC.noalias() +=
1470 NTN * (ip_cv.fC_4_MCpC.ml / (pCap - pCap_prev) * w);
1471 }
1472 }
1473
1474 MCT.noalias() += NTN * (ip_cv.fC_4_MCT.m * w);
1475 // d (fC_4_MCT * T_dot)/d T
1476 local_Jac
1477 .template block<C_size, temperature_size>(C_index,
1478 temperature_index)
1479 .noalias() += NTN * (ip_dd.dfC_4_MCT.dT * (T - T_prev) / dt * w);
1480
1481 MCu.noalias() += BTI2N.transpose() * (ip_cv.fC_4_MCu.m * w);
1482 // d (fC_4_MCu * u_dot)/d T
1483 local_Jac
1484 .template block<C_size, temperature_size>(C_index,
1485 temperature_index)
1486 .noalias() += NTN * (ip_dd.dfC_4_MCu.dT * div_u_dot * w);
1487
1488 LCpG.noalias() += gradNpT * ip_cv.fC_4_LCpG.L * gradNp * w;
1489
1490 // d (fC_4_LCpG * grad p_GR)/d p_GR
1491 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1492 gradNpT * ip_dd.dfC_4_LCpG.dp_GR * gradpGR * Np * w;
1493
1494 // d (fC_4_LCpG * grad p_GR)/d p_cap
1495 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1496 gradNpT * ip_dd.dfC_4_LCpG.dp_cap * gradpGR * Np * w;
1497
1498 // d (fC_4_LCpG * grad p_GR)/d T
1499 local_Jac
1500 .template block<C_size, temperature_size>(C_index,
1501 temperature_index)
1502 .noalias() += gradNpT * ip_dd.dfC_4_LCpG.dT * gradpGR * NT * w;
1503
1504 // d (fC_4_MCpG * p_GR_dot)/d p_GR
1505 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1506 NTN * (ip_dd.dfC_4_MCpG.dp_GR * (pGR - pGR_prev) / dt * w);
1507
1508 // d (fC_4_MCpG * p_GR_dot)/d T
1509 local_Jac
1510 .template block<C_size, temperature_size>(C_index,
1511 temperature_index)
1512 .noalias() +=
1513 NTN * (ip_dd.dfC_4_MCpG.dT * (pGR - pGR_prev) / dt * w);
1514
1515 LCpC.noalias() -= gradNpT * ip_cv.fC_4_LCpC.L * gradNp * w;
1516
1517 /* TODO (naumov) This part is not tested by any of the current ctests.
1518 // d (fC_4_LCpC * grad p_cap)/d p_GR
1519 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1520 gradNpT * ip_dd.dfC_4_LCpC.dp_GR * gradpCap * Np * w;
1521 // d (fC_4_LCpC * grad p_cap)/d p_cap
1522 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1523 gradNpT * ip_dd.dfC_4_LCpC.dp_cap * gradpCap * Np * w;
1524
1525 local_Jac
1526 .template block<C_size, temperature_size>(C_index,
1527 temperature_index)
1528 .noalias() += gradNpT * ip_dd.dfC_4_LCpC.dT * gradpCap * Np * w;
1529 */
1530
1531 LCT.noalias() += gradNpT * ip_cv.fC_4_LCT.L * gradNp * w;
1532
1533 // fC_1
1534 fC.noalias() += gradNpT * ip_cv.fC_1.A * b * w;
1535
1536 if (!this->process_data_.apply_mass_lumping)
1537 {
1538 // fC_2 = \int a * s_L_dot
1539 fC.noalias() -= NpT * (ip_cv.fC_2a.a * s_L_dot * w);
1540
1541 local_Jac.template block<C_size, C_size>(C_index, C_index)
1542 .noalias() +=
1543 NTN * ((ip_dd.dfC_2a.dp_GR * s_L_dot
1544 /*- ip_cv.fC_2a.a * (ds_L_dp_GR = 0) / dt*/) *
1545 w);
1546
1547 local_Jac.template block<C_size, W_size>(C_index, W_index)
1548 .noalias() +=
1549 NTN * ((ip_dd.dfC_2a.dp_cap * s_L_dot +
1550 ip_cv.fC_2a.a * ip_dd.dS_L_dp_cap() / dt) *
1551 w);
1552
1553 local_Jac
1554 .template block<C_size, temperature_size>(C_index,
1555 temperature_index)
1556 .noalias() += NTN * (ip_dd.dfC_2a.dT * s_L_dot * w);
1557 }
1558 {
1559 // fC_3 = \int phi * a
1560 fC.noalias() -=
1561 NpT * (current_state.porosity_data.phi * ip_cv.fC_3a.a * w);
1562
1563 local_Jac.template block<C_size, C_size>(C_index, C_index)
1564 .noalias() += NTN * (current_state.porosity_data.phi *
1565 ip_dd.dfC_3a.dp_GR * w);
1566
1567 local_Jac.template block<C_size, W_size>(C_index, W_index)
1568 .noalias() += NTN * (current_state.porosity_data.phi *
1569 ip_dd.dfC_3a.dp_cap * w);
1570
1571 local_Jac
1572 .template block<C_size, temperature_size>(C_index,
1573 temperature_index)
1574 .noalias() +=
1575 NTN * ((ip_dd.porosity_d_data.dphi_dT * ip_cv.fC_3a.a +
1576 current_state.porosity_data.phi * ip_dd.dfC_3a.dT) *
1577 w);
1578 }
1579 // ---------------------------------------------------------------------
1580 // W-component equation
1581 // ---------------------------------------------------------------------
1582
1583 MWpG.noalias() += NTN * (ip_cv.fW_4_MWpG.m * w);
1584 MWpC.noalias() += NTN * (ip_cv.fW_4_MWpC.m * w);
1585
1586 if (this->process_data_.apply_mass_lumping)
1587 {
1588 if (pCap - pCap_prev != 0.) // avoid division by Zero
1589 {
1590 MWpC.noalias() +=
1591 NTN * (ip_cv.fW_4_MWpC.ml / (pCap - pCap_prev) * w);
1592 }
1593 }
1594
1595 MWT.noalias() += NTN * (ip_cv.fW_4_MWT.m * w);
1596
1597 MWu.noalias() += BTI2N.transpose() * (ip_cv.fW_4_MWu.m * w);
1598
1599 LWpG.noalias() += gradNpT * ip_cv.fW_4_LWpG.L * gradNp * w;
1600
1601 // fW_4 LWpG' parts; LWpG = \int grad (a + d) grad
1602 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1603 gradNpT * ip_dd.dfW_4_LWpG.dp_GR * gradpGR * Np * w;
1604
1605 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1606 gradNpT * ip_dd.dfW_4_LWpG.dp_cap * gradpGR * Np * w;
1607
1608 local_Jac
1609 .template block<W_size, temperature_size>(W_index,
1610 temperature_index)
1611 .noalias() += gradNpT * ip_dd.dfW_4_LWpG.dT * gradpGR * NT * w;
1612
1613 LWpC.noalias() += gradNpT * ip_cv.fW_4_LWpC.L * gradNp * w;
1614
1615 // fW_4 LWp_cap' parts; LWpC = \int grad (a + d) grad
1616 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() -=
1617 gradNpT * ip_dd.dfW_4_LWpC.dp_GR * gradpCap * Np * w;
1618
1619 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() -=
1620 gradNpT * ip_dd.dfW_4_LWpC.dp_cap * gradpCap * Np * w;
1621
1622 local_Jac
1623 .template block<W_size, temperature_size>(W_index,
1624 temperature_index)
1625 .noalias() -= gradNpT * ip_dd.dfW_4_LWpC.dT * gradpCap * NT * w;
1626
1627 LWT.noalias() += gradNpT * ip_cv.fW_4_LWT.L * gradNp * w;
1628
1629 // fW_1
1630 fW.noalias() += gradNpT * ip_cv.fW_1.A * b * w;
1631
1632 // fW_2 = \int a * s_L_dot
1633 if (!this->process_data_.apply_mass_lumping)
1634 {
1635 fW.noalias() -= NpT * (ip_cv.fW_2.a * s_L_dot * w);
1636
1637 local_Jac.template block<W_size, C_size>(W_index, C_index)
1638 .noalias() += NTN * (ip_dd.dfW_2.dp_GR * s_L_dot * w);
1639
1640 // sign negated because of dp_cap = -dp_LR
1641 // TODO (naumov) Had to change the sign to get equal Jacobian WW
1642 // blocks in A2 Test. Where is the error?
1643 local_Jac.template block<W_size, W_size>(W_index, W_index)
1644 .noalias() += NTN * ((ip_dd.dfW_2.dp_cap * s_L_dot +
1645 ip_cv.fW_2.a * ip_dd.dS_L_dp_cap() / dt) *
1646 w);
1647
1648 local_Jac
1649 .template block<W_size, temperature_size>(W_index,
1650 temperature_index)
1651 .noalias() += NTN * (ip_dd.dfW_2.dT * s_L_dot * w);
1652 }
1653
1654 // fW_3 = \int phi * a
1655 fW.noalias() -=
1656 NpT * (current_state.porosity_data.phi * ip_cv.fW_3a.a * w);
1657
1658 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1659 NTN * (current_state.porosity_data.phi * ip_dd.dfW_3a.dp_GR * w);
1660
1661 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1662 NTN * (current_state.porosity_data.phi * ip_dd.dfW_3a.dp_cap * w);
1663
1664 local_Jac
1665 .template block<W_size, temperature_size>(W_index,
1666 temperature_index)
1667 .noalias() +=
1668 NTN * ((ip_dd.porosity_d_data.dphi_dT * ip_cv.fW_3a.a +
1669 current_state.porosity_data.phi * ip_dd.dfW_3a.dT) *
1670 w);
1671
1672 // ---------------------------------------------------------------------
1673 // - temperature equation
1674 // ---------------------------------------------------------------------
1675
1676 MTu.noalias() +=
1677 BTI2N.transpose() *
1678 (ip_cv.effective_volumetric_enthalpy_data.rho_h_eff * w);
1679
1680 // dfT_4/dp_GR
1681 // d (MTu * u_dot)/dp_GR
1682 local_Jac
1683 .template block<temperature_size, C_size>(temperature_index,
1684 C_index)
1685 .noalias() +=
1686 NTN * (ip_dd.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_GR *
1687 div_u_dot * w);
1688
1689 // dfT_4/dp_cap
1690 // d (MTu * u_dot)/dp_cap
1691 local_Jac
1692 .template block<temperature_size, W_size>(temperature_index,
1693 W_index)
1694 .noalias() -=
1695 NTN *
1696 (ip_dd.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_cap *
1697 div_u_dot * w);
1698
1699 // dfT_4/dT
1700 // d (MTu * u_dot)/dT
1701 local_Jac
1702 .template block<temperature_size, temperature_size>(
1703 temperature_index, temperature_index)
1704 .noalias() +=
1705 NTN * (ip_dd.effective_volumetric_enthalpy_d_data.drho_h_eff_dT *
1706 div_u_dot * w);
1707
1708 KTT.noalias() +=
1709 gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w;
1710
1711 // d KTT/dp_GR * T
1712 // TODO (naumov) always zero if lambda_xR have no derivatives wrt. p_GR.
1713 // dlambda_dp_GR =
1714 // (dphi_G_dp_GR = 0) * lambdaGR + phi_G * dlambda_GR_dp_GR +
1715 // (dphi_L_dp_GR = 0) * lambdaLR + phi_L * dlambda_LR_dp_GR +
1716 // (dphi_S_dp_GR = 0) * lambdaSR + phi_S * dlambda_SR_dp_GR +
1717 // = 0
1718 //
1719 // Since dlambda/dp_GR is 0 the derivative is omitted:
1720 // local_Jac
1721 // .template block<temperature_size, C_size>(temperature_index,
1722 // C_index)
1723 // .noalias() += gradNTT * dlambda_dp_GR * gradT * Np * w;
1724
1725 // d KTT/dp_cap * T
1726 local_Jac
1727 .template block<temperature_size, W_size>(temperature_index,
1728 W_index)
1729 .noalias() += gradNTT *
1730 ip_dd.thermal_conductivity_d_data.dlambda_dp_cap *
1731 gradT * Np * w;
1732
1733 // d KTT/dT * T
1734 local_Jac
1735 .template block<temperature_size, temperature_size>(
1736 temperature_index, temperature_index)
1737 .noalias() += gradNTT *
1738 ip_dd.thermal_conductivity_d_data.dlambda_dT * gradT *
1739 NT * w;
1740
1741 // fT_1
1742 fT.noalias() -= NTT * (ip_cv.fT_1.m * w);
1743
1744 // dfT_1/dp_GR
1745 local_Jac
1746 .template block<temperature_size, C_size>(temperature_index,
1747 C_index)
1748 .noalias() += NTN * (ip_dd.dfT_1.dp_GR * w);
1749
1750 // dfT_1/dp_cap
1751 local_Jac
1752 .template block<temperature_size, W_size>(temperature_index,
1753 W_index)
1754 .noalias() += NTN * (ip_dd.dfT_1.dp_cap * w);
1755
1756 // dfT_1/dT
1757 // MTT
1758 local_Jac
1759 .template block<temperature_size, temperature_size>(
1760 temperature_index, temperature_index)
1761 .noalias() += NTN * (ip_dd.dfT_1.dT * w);
1762
1763 // fT_2
1764 fT.noalias() += gradNTT * ip_cv.fT_2.A * w;
1765
1766 // dfT_2/dp_GR
1767 local_Jac
1768 .template block<temperature_size, C_size>(temperature_index,
1769 C_index)
1770 .noalias() -=
1771 // dfT_2/dp_GR first part
1772 gradNTT * ip_dd.dfT_2.dp_GR_Npart * Np * w +
1773 // dfT_2/dp_GR second part
1774 gradNTT * ip_dd.dfT_2.dp_GR_gradNpart * gradNp * w;
1775
1776 // dfT_2/dp_cap
1777 local_Jac
1778 .template block<temperature_size, W_size>(temperature_index,
1779 W_index)
1780 .noalias() -=
1781 // first part of dfT_2/dp_cap
1782 gradNTT * (-ip_dd.dfT_2.dp_cap_Npart) * Np * w +
1783 // second part of dfT_2/dp_cap
1784 gradNTT * (-ip_dd.dfT_2.dp_cap_gradNpart) * gradNp * w;
1785
1786 // dfT_2/dT
1787 local_Jac
1788 .template block<temperature_size, temperature_size>(
1789 temperature_index, temperature_index)
1790 .noalias() -= gradNTT * ip_dd.dfT_2.dT * NT * w;
1791
1792 // fT_3
1793 fT.noalias() += NTT * (ip_cv.fT_3.N * w);
1794
1795 fT.noalias() += gradNTT * ip_cv.fT_3.gradN * w;
1796
1797 // ---------------------------------------------------------------------
1798 // - displacement equation
1799 // ---------------------------------------------------------------------
1800
1801 KUpG.noalias() -= BTI2N * (ip_cv.biot_data() * w);
1802
1803 // dfU_2/dp_GR = dKUpG/dp_GR * p_GR + KUpG. The former is zero, the
1804 // latter is handled below.
1805
1806 KUpC.noalias() += BTI2N * (ip_cv.fu_2_KupC.m * w);
1807
1808 // dfU_2/dp_cap = dKUpC/dp_cap * p_cap + KUpC. The former is handled
1809 // here, the latter below.
1810 local_Jac
1811 .template block<displacement_size, W_size>(displacement_index,
1812 W_index)
1813 .noalias() += BTI2N * (ip_dd.dfu_2_KupC.dp_cap * w);
1814
1815 auto const BuTC =
1816 (Bu.transpose() * ip_cd.s_mech_data.stiffness_tensor).eval();
1817 local_Jac
1818 .template block<displacement_size, displacement_size>(
1819 displacement_index, displacement_index)
1820 .noalias() += BuTC * Bu * w;
1821
1822 // fU_1
1823 fU.noalias() -=
1824 (Bu.transpose() * current_state.eff_stress_data.sigma_eff -
1825 N_u_op(Nu).transpose() * ip_cv.volumetric_body_force()) *
1826 w;
1827
1828 // KuT
1829 local_Jac
1830 .template block<displacement_size, temperature_size>(
1831 displacement_index, temperature_index)
1832 .noalias() -= Bu.transpose() * ip_dd.dfu_1_KuT.dT * NT * w;
1833
1834 /* TODO (naumov) Test with gravity needed to check this Jacobian part.
1835 local_Jac
1836 .template block<displacement_size, temperature_size>(
1837 displacement_index, temperature_index)
1838 .noalias() += N_u_op(Nu).transpose() * ip_cv.drho_dT * b *
1839 N_u_op(Nu).transpose() * w;
1840 */
1841
1842 if (this->process_data_.apply_mass_lumping)
1843 {
1844 MCpG = MCpG.colwise().sum().eval().asDiagonal();
1845 MCpC = MCpC.colwise().sum().eval().asDiagonal();
1846 MWpG = MWpG.colwise().sum().eval().asDiagonal();
1847 MWpC = MWpC.colwise().sum().eval().asDiagonal();
1848 }
1849 } // int_point-loop
1850
1851 // --- Gas ---
1852 // fC_4
1853 fC.noalias() -= LCpG * gas_pressure + LCpC * capillary_pressure +
1854 LCT * temperature +
1855 MCpG * (gas_pressure - gas_pressure_prev) / dt +
1856 MCpC * (capillary_pressure - capillary_pressure_prev) / dt +
1857 MCT * (temperature - temperature_prev) / dt +
1858 MCu * (displacement - displacement_prev) / dt;
1859
1860 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1861 LCpG + MCpG / dt;
1862 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1863 LCpC + MCpC / dt;
1864 local_Jac
1865 .template block<C_size, temperature_size>(C_index, temperature_index)
1866 .noalias() += LCT + MCT / dt;
1867 local_Jac
1868 .template block<C_size, displacement_size>(C_index, displacement_index)
1869 .noalias() += MCu / dt;
1870
1871 // --- Capillary pressure ---
1872 // fW_4
1873 fW.noalias() -= LWpG * gas_pressure + LWpC * capillary_pressure +
1874 LWT * temperature +
1875 MWpG * (gas_pressure - gas_pressure_prev) / dt +
1876 MWpC * (capillary_pressure - capillary_pressure_prev) / dt +
1877 MWT * (temperature - temperature_prev) / dt +
1878 MWu * (displacement - displacement_prev) / dt;
1879
1880 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1881 LWpC + MWpC / dt;
1882 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1883 LWpG + MWpG / dt;
1884 local_Jac
1885 .template block<W_size, temperature_size>(W_index, temperature_index)
1886 .noalias() += LWT + MWT / dt;
1887 local_Jac
1888 .template block<W_size, displacement_size>(W_index, displacement_index)
1889 .noalias() += MWu / dt;
1890
1891 // --- Temperature ---
1892 // fT_4
1893 fT.noalias() -=
1894 KTT * temperature + MTu * (displacement - displacement_prev) / dt;
1895
1896 local_Jac
1897 .template block<temperature_size, temperature_size>(temperature_index,
1898 temperature_index)
1899 .noalias() += KTT;
1900 local_Jac
1901 .template block<temperature_size, displacement_size>(temperature_index,
1902 displacement_index)
1903 .noalias() += MTu / dt;
1904
1905 // --- Displacement ---
1906 // fU_2
1907 fU.noalias() -= KUpG * gas_pressure + KUpC * capillary_pressure;
1908
1909 local_Jac
1910 .template block<displacement_size, C_size>(displacement_index, C_index)
1911 .noalias() += KUpG;
1912 local_Jac
1913 .template block<displacement_size, W_size>(displacement_index, W_index)
1914 .noalias() += KUpC;
1915}
1916
1917template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
1918 int DisplacementDim>
1919void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
1920 DisplacementDim>::
1921 computeSecondaryVariableConcrete(double const t, double const dt,
1922 Eigen::VectorXd const& local_x,
1923 Eigen::VectorXd const& local_x_prev)
1924{
1925 auto const gas_pressure =
1926 local_x.template segment<gas_pressure_size>(gas_pressure_index);
1927 auto const capillary_pressure =
1928 local_x.template segment<capillary_pressure_size>(
1929 capillary_pressure_index);
1930 auto const liquid_pressure = (gas_pressure - capillary_pressure).eval();
1931
1932 NumLib::interpolateToHigherOrderNodes<
1933 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1934 DisplacementDim>(this->element_, this->is_axially_symmetric_,
1935 gas_pressure,
1936 *this->process_data_.gas_pressure_interpolated);
1937
1938 NumLib::interpolateToHigherOrderNodes<
1939 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1940 DisplacementDim>(this->element_, this->is_axially_symmetric_,
1941 capillary_pressure,
1942 *this->process_data_.capillary_pressure_interpolated);
1943
1944 NumLib::interpolateToHigherOrderNodes<
1945 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1946 DisplacementDim>(this->element_, this->is_axially_symmetric_,
1947 liquid_pressure,
1948 *this->process_data_.liquid_pressure_interpolated);
1949
1950 auto const temperature =
1951 local_x.template segment<temperature_size>(temperature_index);
1952
1953 NumLib::interpolateToHigherOrderNodes<
1954 ShapeFunctionPressure, typename ShapeFunctionDisplacement::MeshElement,
1955 DisplacementDim>(this->element_, this->is_axially_symmetric_,
1956 temperature,
1957 *this->process_data_.temperature_interpolated);
1958
1959 ConstitutiveRelations::ConstitutiveModels<DisplacementDim> const models{
1960 this->solid_material_, *this->process_data_.phase_transition_model_};
1961
1962 updateConstitutiveVariables(local_x, local_x_prev, t, dt, models);
1963}
1964
1965} // namespace TH2M
1966} // namespace ProcessLib
OGS_FATAL
#define OGS_FATAL(...)
Definition Error.h:19
DBUG
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:22
WARN
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:34
ParameterLib::SpatialPosition::getCoordinates
std::optional< MathLib::Point3d > const getCoordinates() const
Definition SpatialPosition.h:56
ProcessLib::BMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim >::BMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition BMatrixPolicy.h:48
ProcessLib::TH2M::TH2MLocalAssembler::assemble
void assemble(double const, double const, std::vector< double > const &, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
Definition TH2MFEM-impl.h:975
ProcessLib::TH2M::TH2MLocalAssembler::N_u_op
static constexpr auto & N_u_op
Definition TH2MFEM.h:70
ProcessLib::TH2M::TH2MLocalAssembler::capillary_pressure_index
static const int capillary_pressure_index
Definition TH2MFEM.h:255
ProcessLib::TH2M::TH2MLocalAssembler::ShapeMatricesTypeDisplacement
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
Definition TH2MFEM.h:46
ProcessLib::TH2M::TH2MLocalAssembler::updateConstitutiveVariablesDerivatives
std::vector< ConstitutiveRelations::DerivativesData< DisplacementDim > > updateConstitutiveVariablesDerivatives(Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, std::vector< ConstitutiveRelations::ConstitutiveData< DisplacementDim > > const &ip_constitutive_data, std::vector< ConstitutiveRelations::ConstitutiveTempData< DisplacementDim > > const &ip_constitutive_variables, ConstitutiveRelations::ConstitutiveModels< DisplacementDim > const &models)
Definition TH2MFEM-impl.h:512
ProcessLib::TH2M::TH2MLocalAssembler::_secondary_data
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
Definition TH2MFEM.h:249
ProcessLib::TH2M::TH2MLocalAssembler::capillary_pressure_size
static const int capillary_pressure_size
Definition TH2MFEM.h:256
ProcessLib::TH2M::TH2MLocalAssembler::MatrixType
typename ShapeMatricesTypePressure::template MatrixType< M, N > MatrixType
Definition TH2MFEM.h:57
ProcessLib::TH2M::TH2MLocalAssembler::gas_pressure_index
static const int gas_pressure_index
Definition TH2MFEM.h:253
ProcessLib::TH2M::TH2MLocalAssembler::ShapeMatricesTypePressure
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
Definition TH2MFEM.h:49
ProcessLib::TH2M::TH2MLocalAssembler::updateConstitutiveVariables
std::tuple< std::vector< ConstitutiveRelations::ConstitutiveData< DisplacementDim > >, std::vector< ConstitutiveRelations::ConstitutiveTempData< DisplacementDim > > > updateConstitutiveVariables(Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, ConstitutiveRelations::ConstitutiveModels< DisplacementDim > const &models)
Definition TH2MFEM-impl.h:82
ProcessLib::TH2M::TH2MLocalAssembler::GlobalDimVectorType
typename ShapeMatricesTypePressure::GlobalDimVectorType GlobalDimVectorType
Definition TH2MFEM.h:63
ProcessLib::TH2M::TH2MLocalAssembler::displacement_index
static const int displacement_index
Definition TH2MFEM.h:259
ProcessLib::TH2M::TH2MLocalAssembler::_ip_data
std::vector< IpData > _ip_data
Definition TH2MFEM.h:245
MathLib::createZeroedVector
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Definition EigenMapTools.h:146
MathLib::createZeroedMatrix
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition EigenMapTools.h:25
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:100
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:48
NumLib::interpolateCoordinates
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition InitShapeMatrices.h:93
ProcessLib::ConstitutiveRelations::SpecificBodyForce
BaseLib::StrongType< GlobalDimVector< DisplacementDim >, struct SpecificBodyForceTag > SpecificBodyForce
Definition ConstitutiveRelations/Base.h:32
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:37
ProcessLib::Reflection::reflectSetIPData
std::size_t reflectSetIPData(std::string_view const name, double const *values, std::vector< IPData > &ip_data_vector)
Definition ReflectionSetIPData.h:163
ProcessLib::TH2M::ConstitutiveRelations::GasPressureGradientData
BaseLib::StrongType< GlobalDimVector< DisplacementDim >, struct GasPressureGradientTag > GasPressureGradientData
Definition TH2M/ConstitutiveRelations/Base.h:102
ProcessLib::TH2M::ConstitutiveRelations::TemperatureGradientData
BaseLib::StrongType< GlobalDimVector< DisplacementDim >, struct TemperatureGradientTag > TemperatureGradientData
Definition TH2M/ConstitutiveRelations/Base.h:110
ProcessLib::TH2M::ConstitutiveRelations::CapillaryPressureGradientData
BaseLib::StrongType< GlobalDimVector< DisplacementDim >, struct CapillaryPressureGradientTag > CapillaryPressureGradientData
Definition TH2M/ConstitutiveRelations/Base.h:106
ProcessLib::TH2M::ConstitutiveRelations::ReferenceTemperatureData
BaseLib::StrongType< double, struct ReferenceTemperatureTag > ReferenceTemperatureData
Definition TH2M/ConstitutiveRelations/Base.h:98
ProcessLib::setIntegrationPointDataMaterialStateVariables
std::size_t setIntegrationPointDataMaterialStateVariables(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType member, std::function< std::span< double >(MaterialStateVariables &)> get_values_span)
Definition SetOrGetIntegrationPointData.h:262
ProcessLib::setIPDataInitialConditions
void setIPDataInitialConditions(std::vector< std::unique_ptr< MeshLib::IntegrationPointWriter > > const &_integration_point_writer, MeshLib::Properties const &mesh_properties, LocalAssemblersVector &local_assemblers)
Definition SetIPDataInitialConditions.h:29
MathLib::KelvinVector::Invariants::identity2
static Eigen::Matrix< double, KelvinVectorSize, 1 > const identity2
Kelvin mapping of 2nd order identity tensor.
Definition KelvinVector.h:101
MathLib::KelvinVector::Invariants::trace
static double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.
Definition KelvinVector-impl.h:55
ProcessLib::TH2M::ConstitutiveRelations::BiotModel::eval
void eval(SpaceTimeData const &x_t, MediaData const &media_data, BiotData &out) const
Definition TH2M/ConstitutiveRelations/Biot.cpp:10
ProcessLib::TH2M::ConstitutiveRelations::BishopsModel::eval
void eval(SpaceTimeData const &x_t, MediaData const &media_data, SaturationData const &S_L_data, BishopsData &out) const
Definition Bishops.cpp:28
ProcessLib::TH2M::ConstitutiveRelations::BishopsPrevModel::eval
void eval(SpaceTimeData const &x_t, MediaData const &media_data, PrevState< SaturationData > const &S_L_data, PrevState< BishopsData > &out) const
Definition Bishops.cpp:33
ProcessLib::TH2M::ConstitutiveRelations::ConstitutiveModels::beta_p_SR_model
SolidCompressibilityModel< DisplacementDim, SolidConstitutiveRelation< DisplacementDim > > beta_p_SR_model
Definition TH2M/ConstitutiveRelations/ConstitutiveModels.h:58
ProcessLib::TH2M::ConstitutiveRelations::EffectiveVolumetricEnthalpyModel::dEval
void dEval(FluidDensityData const &fluid_density_data, FluidEnthalpyData const &fluid_enthalpy_data, PhaseTransitionData const &phase_transition_data, PorosityData const &porosity_data, PorosityDerivativeData const &porosity_d_data, SaturationData const &S_L_data, SolidDensityData const &solid_density_data, SolidDensityDerivativeData const &solid_density_d_data, SolidEnthalpyData const &solid_enthalpy_data, SolidHeatCapacityData const &solid_heat_capacity_data, EffectiveVolumetricEnthalpyDerivatives &effective_volumetric_enthalpy_d_data) const
Definition Enthalpy.cpp:29
ProcessLib::TH2M::ConstitutiveRelations::EffectiveVolumetricEnthalpyModel::eval
void eval(FluidDensityData const &fluid_density_data, FluidEnthalpyData const &fluid_enthalpy_data, PorosityData const &porosity_data, SaturationData const &S_L_data, SolidDensityData const &solid_density_data, SolidEnthalpyData const &solid_enthalpy_data, EffectiveVolumetricEnthalpy &effective_volumetric_enthalpy_data) const
Definition Enthalpy.cpp:10
ProcessLib::TH2M::ConstitutiveRelations::FC2aModel::dEval
void dEval(BiotData const &biot_data, CapillaryPressureData const pCap, ConstituentDensityData const &constituent_density_data, PhaseTransitionData const &phase_transition_data, PorosityData const &porosity_data, PorosityDerivativeData const &porosity_d_data, SaturationData const &S_L_data, SaturationDataDeriv const &dS_L_dp_cap, SolidCompressibilityData const &beta_p_SR, FC2aDerivativeData &dfC_2a) const
Definition CEquation.cpp:41
ProcessLib::TH2M::ConstitutiveRelations::FC2aModel::eval
void eval(BiotData const biot_data, CapillaryPressureData const pCap, ConstituentDensityData const &constituent_density_data, PorosityData const &porosity_data, SaturationData const &S_L_data, SolidCompressibilityData const beta_p_SR, FC2aData &fC_2a) const
Definition CEquation.cpp:23
ProcessLib::TH2M::ConstitutiveRelations::FC3aModel::eval
void eval(double const dt, ConstituentDensityData const &constituent_density_data, PrevState< ConstituentDensityData > const &constituent_density_data_prev, SaturationData const &S_L_data, FC3aData &fC_3a) const
Definition CEquation.cpp:96
ProcessLib::TH2M::ConstitutiveRelations::FC3aModel::dEval
void dEval(double const dt, ConstituentDensityData const &constituent_density_data, PrevState< ConstituentDensityData > const &constituent_density_data_prev, PhaseTransitionData const &phase_transition_data, SaturationData const &S_L_data, SaturationDataDeriv const &dS_L_dp_cap, FC3aDerivativeData &dfC_3a) const
Definition CEquation.cpp:120
ProcessLib::TH2M::ConstitutiveRelations::FC4MCpCModel::eval
void eval(BiotData const &biot_data, CapillaryPressureData const pCap, ConstituentDensityData const &constituent_density_data, PorosityData const &porosity_data, PrevState< SaturationData > const &S_L_data_prev, SaturationData const &S_L_data, SolidCompressibilityData const &beta_p_SR, FC4MCpCData &fC_4_MCpC) const
Definition CEquation.cpp:361
ProcessLib::TH2M::ConstitutiveRelations::FC4MCpGModel::eval
void eval(BiotData const &biot_data, ConstituentDensityData const &constituent_density_data, PorosityData const &porosity_data, SaturationData const &S_L_data, SolidCompressibilityData const &beta_p_SR, FC4MCpGData &fC_4_MCpG) const
Definition CEquation.cpp:316
ProcessLib::TH2M::ConstitutiveRelations::FC4MCpGModel::dEval
void dEval(BiotData const &biot_data, ConstituentDensityData const &constituent_density_data, PhaseTransitionData const &phase_transition_data, PorosityData const &porosity_data, PorosityDerivativeData const &porosity_d_data, SaturationData const &S_L_data, SolidCompressibilityData const &beta_p_SR, FC4MCpGDerivativeData &dfC_4_MCpG) const
Definition CEquation.cpp:331
ProcessLib::TH2M::ConstitutiveRelations::FC4MCuModel::eval
void eval(BiotData const &biot_data, ConstituentDensityData const &constituent_density_data, SaturationData const &S_L_data, FC4MCuData &fC_4_MCu) const
Definition CEquation.cpp:436
ProcessLib::TH2M::ConstitutiveRelations::FC4MCuModel::dEval
void dEval(BiotData const &biot_data, PhaseTransitionData const &phase_transition_data, SaturationData const &S_L_data, FC4MCuDerivativeData &dfC_4_MCu) const
Definition CEquation.cpp:449
ProcessLib::TH2M::ConstitutiveRelations::FT1Model::dEval
void dEval(double const dt, EffectiveVolumetricInternalEnergyDerivatives const &effective_volumetric_internal_energy_d_data, FT1DerivativeData &dfT_1) const
Definition TEquation.cpp:27
ProcessLib::TH2M::ConstitutiveRelations::FT1Model::eval
void eval(double const dt, InternalEnergyData const &internal_energy_data, PrevState< InternalEnergyData > const &internal_energy_data_prev, FT1Data &fT_1) const
Definition TEquation.cpp:10
ProcessLib::TH2M::ConstitutiveRelations::FU2KUpCModel::dEval
void dEval(BiotData const &biot_data, BishopsData const &chi_S_L, CapillaryPressureData const &p_cap, SaturationDataDeriv const &dS_L_dp_cap, FU2KUpCDerivativeData &dfu_2_KupC) const
Definition UEquation.cpp:30
ProcessLib::TH2M::ConstitutiveRelations::FU2KUpCModel::eval
void eval(BiotData const &biot_data, BishopsData const &chi_S_L, FU2KUpCData &fu_2_KupC) const
Definition UEquation.cpp:23
ProcessLib::TH2M::ConstitutiveRelations::FW2Model::dEval
void dEval(BiotData const &biot_data, CapillaryPressureData const pCap, ConstituentDensityData const &constituent_density_data, PhaseTransitionData const &phase_transition_data, PorosityData const &porosity_data, PorosityDerivativeData const &porosity_d_data, PureLiquidDensityData const &rho_W_LR, SaturationData const &S_L_data, SaturationDataDeriv const &dS_L_dp_cap, SolidCompressibilityData const &beta_p_SR, FW2DerivativeData &dfW_2) const
Definition WEquation.cpp:42
ProcessLib::TH2M::ConstitutiveRelations::FW2Model::eval
void eval(BiotData const biot_data, CapillaryPressureData const pCap, ConstituentDensityData const &constituent_density_data, PorosityData const &porosity_data, PureLiquidDensityData const &rho_W_LR, SaturationData const &S_L_data, SolidCompressibilityData const beta_p_SR, FW2Data &fW_2) const
Definition WEquation.cpp:23
ProcessLib::TH2M::ConstitutiveRelations::FW3aModel::dEval
void dEval(double const dt, ConstituentDensityData const &constituent_density_data, PhaseTransitionData const &phase_transition_data, PrevState< ConstituentDensityData > const &constituent_density_data_prev, PrevState< PureLiquidDensityData > const &rho_W_LR_prev, PureLiquidDensityData const &rho_W_LR, SaturationData const &S_L_data, SaturationDataDeriv const &dS_L_dp_cap, FW3aDerivativeData &dfW_3a) const
Definition WEquation.cpp:126
ProcessLib::TH2M::ConstitutiveRelations::FW3aModel::eval
void eval(double const dt, ConstituentDensityData const &constituent_density_data, PrevState< ConstituentDensityData > const &constituent_density_data_prev, PrevState< PureLiquidDensityData > const &rho_W_LR_prev, PureLiquidDensityData const &rho_W_LR, SaturationData const &S_L_data, FW3aData &fW_3a) const
Definition WEquation.cpp:102
ProcessLib::TH2M::ConstitutiveRelations::FW4MWpCModel::eval
void eval(BiotData const &biot_data, CapillaryPressureData const pCap, ConstituentDensityData const &constituent_density_data, PorosityData const &porosity_data, PrevState< SaturationData > const &S_L_data_prev, PureLiquidDensityData const &rho_W_LR, SaturationData const &S_L_data, SolidCompressibilityData const &beta_p_SR, FW4MWpCData &fW_4_MWpC) const
Definition WEquation.cpp:382
ProcessLib::TH2M::ConstitutiveRelations::FW4MWpGModel::eval
void eval(BiotData const &biot_data, ConstituentDensityData const &constituent_density_data, PorosityData const &porosity_data, PureLiquidDensityData const &rho_W_LR, SaturationData const &S_L_data, SolidCompressibilityData const &beta_p_SR, FW4MWpGData &fW_4_MWpG) const
Definition WEquation.cpp:365
ProcessLib::TH2M::ConstitutiveRelations::FW4MWuModel::eval
void eval(BiotData const &biot_data, ConstituentDensityData const &constituent_density_data, PureLiquidDensityData const &rho_W_LR, SaturationData const &S_L_data, FW4MWuData &fW_4_MWu) const
Definition WEquation.cpp:429
ProcessLib::TH2M::ConstitutiveRelations::InternalEnergyModel::eval
void eval(FluidDensityData const &fluid_density_data, PhaseTransitionData const &phase_transition_data, PorosityData const &porosity_data, SaturationData const &S_L_data, SolidDensityData const &solid_density_data, SolidEnthalpyData const &solid_enthalpy_data, InternalEnergyData &internal_energy_data) const
Definition InternalEnergy.cpp:10
ProcessLib::TH2M::ConstitutiveRelations::InternalEnergyModel::dEval
void dEval(FluidDensityData const &fluid_density_data, PhaseTransitionData const &phase_transition_data, PorosityData const &porosity_data, PorosityDerivativeData const &porosity_d_data, SaturationData const &S_L_data, SolidDensityData const &solid_density_data, SolidDensityDerivativeData const &solid_density_d_data, SolidEnthalpyData const &solid_enthalpy_data, SolidHeatCapacityData const &solid_heat_capacity_data, EffectiveVolumetricInternalEnergyDerivatives &effective_volumetric_internal_energy_d_data) const
Definition InternalEnergy.cpp:30
ProcessLib::TH2M::ConstitutiveRelations::PhaseTransitionModel::eval
virtual void eval(SpaceTimeData const &x_t, MediaData const &media_data, GasPressureData const &p_GR, CapillaryPressureData const &p_cap, TemperatureData const &T_data, PureLiquidDensityData const &rho_W_LR, FluidEnthalpyData &fluid_enthalpy_data, MassMoleFractionsData &mass_mole_fractions_data, FluidDensityData &fluid_density_data, VapourPartialPressureData &vapour_pressure_data, ConstituentDensityData &constituent_density_data, PhaseTransitionData &cv) const =0
ProcessLib::TH2M::ConstitutiveRelations::PureLiquidDensityModel::eval
void eval(SpaceTimeData const &x_t, MediaData const &media_data, GasPressureData const &p_GR, CapillaryPressureData const &p_cap, TemperatureData const &T_data, PureLiquidDensityData &out) const
Definition PureLiquidDensity.h:18
ProcessLib::TH2M::ConstitutiveRelations::SaturationModel::dEval
void dEval(SpaceTimeData const &x_t, MediaData const &media_data, CapillaryPressureData const &p_cap, SaturationDataDeriv &dS_L_data) const
Definition Saturation.cpp:23
ProcessLib::TH2M::ConstitutiveRelations::SaturationModel::eval
void eval(SpaceTimeData const &x_t, MediaData const &media_data, CapillaryPressureData const &p_cap, SaturationData &S_L_data) const
Definition Saturation.cpp:10
ProcessLib::TH2M::ConstitutiveRelations::SolidEnthalpyModel::eval
void eval(SolidHeatCapacityData const &solid_heat_capacity_data, TemperatureData const &T_data, SolidEnthalpyData &solid_enthalpy_data) const
Definition Enthalpy.cpp:86
ProcessLib::TH2M::ConstitutiveRelations::SolidHeatCapacityModel::eval
void eval(SpaceTimeData const &x_t, MediaData const &media_data, TemperatureData const &T_data, SolidHeatCapacityData &solid_heat_capacity) const
Definition SolidHeatCapacity.cpp:11
ProcessLib::TH2M::ConstitutiveRelations::ViscosityModel::eval
void eval(SpaceTimeData const &x_t, MediaData const &media_data, TemperatureData const &T_data, MassMoleFractionsData const &mass_mole_fractions_data, ViscosityData &viscosity_data) const
Definition Viscosity.cpp:10
ProcessLib::TH2M::LocalAssemblerInterface::integration_method_
NumLib::GenericIntegrationMethod const & integration_method_
Definition TH2M/LocalAssemblerInterface.h:120
ProcessLib::TH2M::LocalAssemblerInterface::LocalAssemblerInterface
LocalAssemblerInterface(MeshLib::Element const &e, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, TH2MProcessData< DisplacementDim > &process_data)
Definition TH2M/LocalAssemblerInterface.h:25
ProcessLib::TH2M::LocalAssemblerInterface::process_data_
TH2MProcessData< DisplacementDim > & process_data_
Definition TH2M/LocalAssemblerInterface.h:104
ProcessLib::TH2M::LocalAssemblerInterface::solid_material_
ConstitutiveRelations::SolidConstitutiveRelation< DisplacementDim > const & solid_material_
Definition TH2M/LocalAssemblerInterface.h:124
ProcessLib::TH2M::LocalAssemblerInterface::current_states_
std::vector< typename ConstitutiveRelations::StatefulData< DisplacementDim > > current_states_
Definition TH2M/LocalAssemblerInterface.h:107
ProcessLib::TH2M::LocalAssemblerInterface::prev_states_
std::vector< typename ConstitutiveRelations::StatefulDataPrev< DisplacementDim > > prev_states_
Definition TH2M/LocalAssemblerInterface.h:111
ProcessLib::TH2M::LocalAssemblerInterface::output_data_
std::vector< ConstitutiveRelations::OutputData< DisplacementDim > > output_data_
Definition TH2M/LocalAssemblerInterface.h:118
ProcessLib::TH2M::LocalAssemblerInterface::material_states_
std::vector< ConstitutiveRelations::MaterialStateData< DisplacementDim > > material_states_
Definition TH2M/LocalAssemblerInterface.h:116