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