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