88 updateConstitutiveVariables(
89 Eigen::VectorXd
const& local_x, Eigen::VectorXd
const& local_x_prev,
90 double const t,
double const dt,
94 [[maybe_unused]]
auto const matrix_size =
95 gas_pressure_size + capillary_pressure_size + temperature_size +
98 assert(local_x.size() == matrix_size);
100 auto const gas_pressure =
101 local_x.template segment<gas_pressure_size>(gas_pressure_index);
102 auto const gas_pressure_prev =
103 local_x_prev.template segment<gas_pressure_size>(gas_pressure_index);
104 auto const capillary_pressure =
105 local_x.template segment<capillary_pressure_size>(
106 capillary_pressure_index);
107 auto const capillary_pressure_prev =
108 local_x_prev.template segment<capillary_pressure_size>(
109 capillary_pressure_index);
111 auto const temperature =
112 local_x.template segment<temperature_size>(temperature_index);
113 auto const temperature_prev =
114 local_x_prev.template segment<temperature_size>(temperature_index);
116 auto const displacement =
117 local_x.template segment<displacement_size>(displacement_index);
118 auto const displacement_prev =
119 local_x_prev.template segment<displacement_size>(displacement_index);
122 *this->process_data_.media_map.getMedium(this->element_.getID());
125 unsigned const n_integration_points =
126 this->integration_method_.getNumberOfPoints();
128 std::vector<ConstitutiveRelations::ConstitutiveData<DisplacementDim>>
129 ip_constitutive_data(n_integration_points);
130 std::vector<ConstitutiveRelations::ConstitutiveTempData<DisplacementDim>>
131 ip_constitutive_variables(n_integration_points);
133 for (
unsigned ip = 0; ip < n_integration_points; ip++)
135 auto& ip_data = _ip_data[ip];
136 auto& ip_cv = ip_constitutive_variables[ip];
137 auto& ip_cd = ip_constitutive_data[ip];
138 auto& ip_out = this->output_data_[ip];
139 auto& current_state = this->current_states_[ip];
140 auto& prev_state = this->prev_states_[ip];
142 auto const& Np = ip_data.N_p;
144 auto const& Nu = ip_data.N_u;
145 auto const& gradNu = ip_data.dNdx_u;
146 auto const& gradNp = ip_data.dNdx_p;
148 std::nullopt, this->element_.getID(),
152 this->element_, Nu))};
158 double const T = NT.dot(temperature);
159 double const T_prev = NT.dot(temperature_prev);
160 double const pG = Np.dot(gas_pressure);
161 double const pG_prev = Np.dot(gas_pressure_prev);
162 double const pCap = Np.dot(capillary_pressure);
163 double const pCap_prev = Np.dot(capillary_pressure_prev);
169 this->process_data_.reference_temperature(t, pos)[0]};
171 grad_p_GR{gradNp * gas_pressure};
173 DisplacementDim>
const grad_p_cap{gradNp * capillary_pressure};
175 grad_T{gradNp * temperature};
181 models.
biot_model.
eval({pos, t, dt}, media_data, ip_cv.biot_data);
185 ShapeFunctionDisplacement::NPOINTS,
187 gradNu, Nu, x_coord, this->is_axially_symmetric_);
189 ip_out.eps_data.eps.noalias() = Bu * displacement;
191 current_state.S_L_data);
194 current_state.S_L_data,
195 current_state.chi_S_L);
198 prev_state.S_L_data, prev_state.chi_S_L);
202 ip_cv.C_el_data, ip_cv.beta_p_SR);
206 {pos, t, dt}, media_data, ip_cv.C_el_data, current_state.S_L_data,
207 prev_state.S_L_data, prev_state.swelling_data,
208 current_state.swelling_data, ip_cv.swelling_data);
212 ip_cv.s_therm_exp_data);
215 T_data, ip_cv.s_therm_exp_data, ip_out.eps_data,
216 Bu * displacement_prev, prev_state.mechanical_strain_data,
217 ip_cv.swelling_data, current_state.mechanical_strain_data);
220 {pos, t, dt}, T_data, current_state.mechanical_strain_data,
221 prev_state.mechanical_strain_data, prev_state.eff_stress_data,
222 current_state.eff_stress_data, this->material_states_[ip],
223 ip_cd.s_mech_data, ip_cv.equivalent_plastic_strain_data);
226 ip_cv.biot_data, current_state.chi_S_L,
228 ip_cv.total_stress_data);
231 pGR_data, pCap_data, T_data,
232 current_state.rho_W_LR);
235 {pos, t, dt}, media_data, pGR_data, pCap_data, T_data,
236 current_state.rho_W_LR, ip_out.fluid_enthalpy_data,
237 ip_out.mass_mole_fractions_data, ip_out.fluid_density_data,
238 ip_out.vapour_pressure_data, current_state.constituent_density_data,
239 ip_cv.phase_transition_data);
242 ip_out.mass_mole_fractions_data,
243 ip_cv.viscosity_data);
246 {pos, t, dt}, media_data, current_state.S_L_data,
247 prev_state.S_L_data, pCap_data, pGR_data, current_state.chi_S_L,
248 prev_state.chi_S_L, ip_cv.beta_p_SR, ip_out.eps_data,
249 Bu * displacement_prev, prev_state.porosity_data,
250 current_state.porosity_data);
252 if (medium.hasProperty(MPL::PropertyType::transport_porosity))
255 {pos, t, dt}, media_data, current_state.S_L_data,
256 prev_state.S_L_data, pCap_data, pGR_data, current_state.chi_S_L,
257 prev_state.chi_S_L, ip_cv.beta_p_SR,
258 current_state.mechanical_strain_data,
259 prev_state.mechanical_strain_data,
260 prev_state.transport_porosity_data, current_state.porosity_data,
261 current_state.transport_porosity_data);
265 current_state.transport_porosity_data.phi =
266 current_state.porosity_data.phi;
270 {pos, t, dt}, media_data, current_state.S_L_data, pGR_data,
271 pCap_data, T_data, current_state.transport_porosity_data,
272 ip_cv.total_stress_data, current_state.mechanical_strain_data,
273 ip_out.eps_data, ip_cv.equivalent_plastic_strain_data,
274 ip_out.permeability_data);
277 {pos, t, dt}, media_data, T_data, current_state.eff_stress_data,
278 pCap_data, pGR_data, current_state.chi_S_L,
279 current_state.porosity_data, ip_out.solid_density_data);
282 ip_cv.solid_heat_capacity_data);
285 {pos, t, dt}, media_data, T_data, current_state.porosity_data,
286 current_state.S_L_data, ip_cv.thermal_conductivity_data);
289 ip_out.permeability_data,
290 current_state.rho_W_LR,
291 ip_cv.viscosity_data,
292 ip_cv.advection_data);
295 ip_out.fluid_density_data,
296 current_state.porosity_data,
297 current_state.S_L_data,
298 ip_out.solid_density_data,
300 this->process_data_.specific_body_force},
301 ip_cv.volumetric_body_force);
305 ip_out.mass_mole_fractions_data,
306 ip_cv.phase_transition_data,
307 current_state.porosity_data,
308 current_state.S_L_data,
310 ip_out.diffusion_velocity_data);
313 ip_out.solid_enthalpy_data);
316 ip_cv.phase_transition_data,
317 current_state.porosity_data,
318 current_state.S_L_data,
319 ip_out.solid_density_data,
320 ip_out.solid_enthalpy_data,
321 current_state.internal_energy_data);
324 ip_out.fluid_density_data,
325 ip_out.fluid_enthalpy_data,
326 current_state.porosity_data,
327 current_state.S_L_data,
328 ip_out.solid_density_data,
329 ip_out.solid_enthalpy_data,
330 ip_cv.effective_volumetric_enthalpy_data);
332 models.
fC_1_model.eval(ip_cv.advection_data, ip_out.fluid_density_data,
335 if (!this->process_data_.apply_mass_lumping)
339 current_state.constituent_density_data,
340 current_state.porosity_data,
341 current_state.S_L_data,
346 current_state.constituent_density_data,
347 prev_state.constituent_density_data,
348 current_state.S_L_data,
352 ip_out.fluid_density_data,
353 ip_cv.phase_transition_data,
354 current_state.porosity_data,
355 current_state.S_L_data,
359 ip_out.fluid_density_data,
360 ip_cv.phase_transition_data,
361 current_state.porosity_data,
362 current_state.S_L_data,
366 ip_cv.phase_transition_data,
367 current_state.porosity_data,
368 current_state.S_L_data,
372 current_state.constituent_density_data,
373 current_state.porosity_data,
374 current_state.S_L_data,
380 current_state.constituent_density_data,
381 current_state.porosity_data,
383 current_state.S_L_data,
388 current_state.constituent_density_data,
389 current_state.porosity_data,
390 current_state.S_L_data,
391 ip_cv.s_therm_exp_data,
395 current_state.constituent_density_data,
396 current_state.S_L_data,
399 models.
fW_1_model.eval(ip_cv.advection_data, ip_out.fluid_density_data,
402 if (!this->process_data_.apply_mass_lumping)
406 current_state.constituent_density_data,
407 current_state.porosity_data,
408 current_state.rho_W_LR,
409 current_state.S_L_data,
414 current_state.constituent_density_data,
415 prev_state.constituent_density_data,
417 current_state.rho_W_LR,
418 current_state.S_L_data,
422 ip_out.fluid_density_data,
423 ip_cv.phase_transition_data,
424 current_state.porosity_data,
425 current_state.S_L_data,
429 ip_out.fluid_density_data,
430 ip_cv.phase_transition_data,
431 current_state.porosity_data,
432 current_state.S_L_data,
436 ip_cv.phase_transition_data,
437 current_state.porosity_data,
438 current_state.S_L_data,
442 current_state.constituent_density_data,
443 current_state.porosity_data,
444 current_state.rho_W_LR,
445 current_state.S_L_data,
451 current_state.constituent_density_data,
452 current_state.porosity_data,
454 current_state.rho_W_LR,
455 current_state.S_L_data,
460 current_state.constituent_density_data,
461 current_state.porosity_data,
462 current_state.rho_W_LR,
463 current_state.S_L_data,
464 ip_cv.s_therm_exp_data,
468 current_state.constituent_density_data,
469 current_state.rho_W_LR,
470 current_state.S_L_data,
474 current_state.internal_energy_data,
475 prev_state.internal_energy_data,
484 ip_out.fluid_density_data,
486 ip_out.permeability_data,
488 this->process_data_.specific_body_force},
489 ip_cv.viscosity_data,
490 ip_out.darcy_velocity_data);
492 models.
fT_2_model.eval(ip_out.darcy_velocity_data,
493 ip_out.fluid_density_data,
494 ip_out.fluid_enthalpy_data,
498 current_state.constituent_density_data,
499 ip_out.darcy_velocity_data,
500 ip_out.diffusion_velocity_data,
501 ip_out.fluid_density_data,
502 ip_cv.phase_transition_data,
504 this->process_data_.specific_body_force},
511 return {ip_constitutive_data, ip_constitutive_variables};
519 updateConstitutiveVariablesDerivatives(
520 Eigen::VectorXd
const& local_x, Eigen::VectorXd
const& local_x_prev,
521 double const t,
double const dt,
524 ip_constitutive_data,
527 ip_constitutive_variables,
531 [[maybe_unused]]
auto const matrix_size =
532 gas_pressure_size + capillary_pressure_size + temperature_size +
535 assert(local_x.size() == matrix_size);
537 auto const gas_pressure =
538 local_x.template segment<gas_pressure_size>(gas_pressure_index);
539 auto const gas_pressure_prev =
540 local_x_prev.template segment<gas_pressure_size>(gas_pressure_index);
541 auto const temperature =
542 local_x.template segment<temperature_size>(temperature_index);
543 auto const temperature_prev =
544 local_x_prev.template segment<temperature_size>(temperature_index);
545 auto const displacement_prev =
546 local_x_prev.template segment<displacement_size>(displacement_index);
548 auto const capillary_pressure =
549 local_x.template segment<capillary_pressure_size>(
550 capillary_pressure_index);
551 auto const capillary_pressure_prev =
552 local_x_prev.template segment<capillary_pressure_size>(
553 capillary_pressure_index);
556 *this->process_data_.media_map.getMedium(this->element_.getID());
559 unsigned const n_integration_points =
560 this->integration_method_.getNumberOfPoints();
562 std::vector<ConstitutiveRelations::DerivativesData<DisplacementDim>>
563 ip_d_data(n_integration_points);
565 for (
unsigned ip = 0; ip < n_integration_points; ip++)
567 auto const& ip_data = _ip_data[ip];
568 auto& ip_dd = ip_d_data[ip];
569 auto const& ip_cd = ip_constitutive_data[ip];
570 auto const& ip_cv = ip_constitutive_variables[ip];
571 auto const& ip_out = this->output_data_[ip];
572 auto const& current_state = this->current_states_[ip];
573 auto const& prev_state = this->prev_states_[ip];
575 auto const& Nu = ip_data.N_u;
576 auto const& Np = ip_data.N_p;
578 auto const& gradNu = ip_data.dNdx_u;
585 std::nullopt, this->element_.getID(),
589 this->element_, Nu))};
591 double const T = NT.dot(temperature);
592 double const T_prev = NT.dot(temperature_prev);
593 double const pG = Np.dot(gas_pressure);
594 double const pG_prev = Np.dot(gas_pressure_prev);
595 double const pCap = Np.dot(capillary_pressure);
596 double const pCap_prev = Np.dot(capillary_pressure_prev);
604 ShapeFunctionDisplacement::NPOINTS,
606 gradNu, Nu, x_coord, this->is_axially_symmetric_);
612 ip_out.permeability_data,
613 ip_cv.viscosity_data,
615 ip_cv.phase_transition_data,
616 ip_dd.advection_d_data);
619 {pos, t, dt}, media_data, current_state.S_L_data,
620 prev_state.S_L_data, pCap_data, pGR_data, current_state.chi_S_L,
621 prev_state.chi_S_L, ip_cv.beta_p_SR, ip_out.eps_data,
622 Bu * displacement_prev, prev_state.porosity_data,
623 ip_dd.porosity_d_data);
626 {pos, t, dt}, media_data, T_data, current_state.porosity_data,
627 ip_dd.porosity_d_data, current_state.S_L_data,
628 ip_dd.thermal_conductivity_d_data);
631 {pos, t, dt}, media_data, T_data, current_state.eff_stress_data,
632 pCap_data, pGR_data, current_state.chi_S_L,
633 current_state.porosity_data, ip_dd.solid_density_d_data);
636 ip_out.fluid_density_data,
637 ip_cv.phase_transition_data,
638 current_state.porosity_data,
639 ip_dd.porosity_d_data,
640 current_state.S_L_data,
641 ip_out.solid_density_data,
642 ip_dd.solid_density_d_data,
643 ip_out.solid_enthalpy_data,
644 ip_cv.solid_heat_capacity_data,
645 ip_dd.effective_volumetric_internal_energy_d_data);
648 ip_out.fluid_density_data,
649 ip_out.fluid_enthalpy_data,
650 ip_cv.phase_transition_data,
651 current_state.porosity_data,
652 ip_dd.porosity_d_data,
653 current_state.S_L_data,
654 ip_out.solid_density_data,
655 ip_dd.solid_density_d_data,
656 ip_out.solid_enthalpy_data,
657 ip_cv.solid_heat_capacity_data,
658 ip_dd.effective_volumetric_enthalpy_d_data);
659 if (!this->process_data_.apply_mass_lumping)
663 current_state.constituent_density_data,
664 ip_cv.phase_transition_data,
665 current_state.porosity_data,
666 ip_dd.porosity_d_data,
667 current_state.S_L_data,
673 current_state.constituent_density_data,
674 prev_state.constituent_density_data,
675 ip_cv.phase_transition_data,
676 current_state.S_L_data,
681 ip_cv.viscosity_data,
682 ip_cv.phase_transition_data,
683 ip_dd.advection_d_data,
687 ip_out.permeability_data,
688 ip_cv.phase_transition_data,
690 ip_cv.viscosity_data,
694 current_state.constituent_density_data,
695 ip_cv.phase_transition_data,
696 current_state.porosity_data,
697 ip_dd.porosity_d_data,
698 current_state.S_L_data,
703 current_state.constituent_density_data,
704 ip_cv.phase_transition_data,
705 current_state.porosity_data,
706 ip_dd.porosity_d_data,
707 current_state.S_L_data,
708 ip_cv.s_therm_exp_data,
712 ip_cv.phase_transition_data,
713 current_state.S_L_data,
716 if (!this->process_data_.apply_mass_lumping)
720 current_state.constituent_density_data,
721 ip_cv.phase_transition_data,
722 current_state.porosity_data,
723 ip_dd.porosity_d_data,
724 current_state.rho_W_LR,
725 current_state.S_L_data,
732 current_state.constituent_density_data,
733 ip_cv.phase_transition_data,
734 prev_state.constituent_density_data,
736 current_state.rho_W_LR,
737 current_state.S_L_data,
742 ip_out.permeability_data,
743 ip_cv.phase_transition_data,
744 current_state.rho_W_LR,
746 ip_cv.viscosity_data,
750 ip_out.fluid_density_data,
751 ip_out.permeability_data,
752 ip_cv.phase_transition_data,
753 current_state.porosity_data,
754 current_state.rho_W_LR,
755 current_state.S_L_data,
757 ip_cv.viscosity_data,
761 dt, ip_dd.effective_volumetric_internal_energy_d_data, ip_dd.dfT_1);
764 ip_out.darcy_velocity_data,
765 ip_out.fluid_density_data,
766 ip_out.fluid_enthalpy_data,
767 ip_out.permeability_data,
768 ip_cv.phase_transition_data,
770 this->process_data_.specific_body_force},
771 ip_cv.viscosity_data,
774 models.
fu_1_KuT_model.dEval(ip_cd.s_mech_data, ip_cv.s_therm_exp_data,
778 current_state.chi_S_L,
984 DisplacementDim>::assemble(
double const t,
double const dt,
985 std::vector<double>
const& local_x,
986 std::vector<double>
const& local_x_prev,
987 std::vector<double>& local_M_data,
988 std::vector<double>& local_K_data,
989 std::vector<double>& local_rhs_data)
991 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
992 temperature_size + displacement_size;
993 assert(local_x.size() == matrix_size);
995 auto const capillary_pressure =
996 Eigen::Map<VectorType<capillary_pressure_size>
const>(
997 local_x.data() + capillary_pressure_index, capillary_pressure_size);
999 auto const capillary_pressure_prev =
1000 Eigen::Map<VectorType<capillary_pressure_size>
const>(
1001 local_x_prev.data() + capillary_pressure_index,
1002 capillary_pressure_size);
1007 local_M_data, matrix_size, matrix_size);
1012 local_K_data, matrix_size, matrix_size);
1016 local_rhs_data, matrix_size);
1022 auto MCpG = local_M.template block<C_size, gas_pressure_size>(
1023 C_index, gas_pressure_index);
1024 auto MCpC = local_M.template block<C_size, capillary_pressure_size>(
1025 C_index, capillary_pressure_index);
1026 auto MCT = local_M.template block<C_size, temperature_size>(
1027 C_index, temperature_index);
1028 auto MCu = local_M.template block<C_size, displacement_size>(
1029 C_index, displacement_index);
1032 auto LCpG = local_K.template block<C_size, gas_pressure_size>(
1033 C_index, gas_pressure_index);
1034 auto LCpC = local_K.template block<C_size, capillary_pressure_size>(
1035 C_index, capillary_pressure_index);
1036 auto LCT = local_K.template block<C_size, temperature_size>(
1037 C_index, temperature_index);
1040 auto MWpG = local_M.template block<W_size, gas_pressure_size>(
1041 W_index, gas_pressure_index);
1042 auto MWpC = local_M.template block<W_size, capillary_pressure_size>(
1043 W_index, capillary_pressure_index);
1044 auto MWT = local_M.template block<W_size, temperature_size>(
1045 W_index, temperature_index);
1046 auto MWu = local_M.template block<W_size, displacement_size>(
1047 W_index, displacement_index);
1050 auto LWpG = local_K.template block<W_size, gas_pressure_size>(
1051 W_index, gas_pressure_index);
1052 auto LWpC = local_K.template block<W_size, capillary_pressure_size>(
1053 W_index, capillary_pressure_index);
1054 auto LWT = local_K.template block<W_size, temperature_size>(
1055 W_index, temperature_index);
1058 auto MTu = local_M.template block<temperature_size, displacement_size>(
1059 temperature_index, displacement_index);
1062 auto KTT = local_K.template block<temperature_size, temperature_size>(
1063 temperature_index, temperature_index);
1066 auto KUpG = local_K.template block<displacement_size, gas_pressure_size>(
1067 displacement_index, gas_pressure_index);
1069 local_K.template block<displacement_size, capillary_pressure_size>(
1070 displacement_index, capillary_pressure_index);
1072 auto KUU = local_K.template block<displacement_size, displacement_size>(
1073 displacement_index, displacement_index);
1076 auto fC = local_f.template segment<C_size>(C_index);
1078 auto fW = local_f.template segment<W_size>(W_index);
1080 auto fT = local_f.template segment<temperature_size>(temperature_index);
1082 auto fU = local_f.template segment<displacement_size>(displacement_index);
1084 unsigned const n_integration_points =
1085 this->integration_method_.getNumberOfPoints();
1088 this->solid_material_, *this->process_data_.phase_transition_model_};
1090 auto const [ip_constitutive_data, ip_constitutive_variables] =
1091 updateConstitutiveVariables(
1092 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1093 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1094 local_x_prev.size()),
1097 for (
unsigned int_point = 0; int_point < n_integration_points; int_point++)
1099 auto& ip = _ip_data[int_point];
1100 auto& ip_cv = ip_constitutive_variables[int_point];
1101 auto& ip_cd = ip_constitutive_data[int_point];
1103 auto& current_state = this->current_states_[int_point];
1104 auto const& prev_state = this->prev_states_[int_point];
1106 auto const& Np = ip.N_p;
1107 auto const& NT = Np;
1108 auto const& Nu = ip.N_u;
1110 std::nullopt, this->element_.getID(),
1114 this->element_, Nu))};
1116 auto const& NpT = Np.transpose().eval();
1117 auto const& NTT = NT.transpose().eval();
1119 auto const& gradNp = ip.dNdx_p;
1120 auto const& gradNT = gradNp;
1121 auto const& gradNu = ip.dNdx_u;
1123 auto const& gradNpT = gradNp.transpose().eval();
1124 auto const& gradNTT = gradNT.transpose().eval();
1126 auto const& w = ip.integration_weight;
1128 auto const x_coord =
1131 this->element_, Nu);
1135 ShapeFunctionDisplacement::NPOINTS,
1137 gradNu, Nu, x_coord, this->is_axially_symmetric_);
1139 auto const NTN = (Np.transpose() * Np).eval();
1140 auto const BTI2N = (Bu.transpose() * Invariants::identity2 * Np).eval();
1142 double const pCap = Np.dot(capillary_pressure);
1143 double const pCap_prev = Np.dot(capillary_pressure_prev);
1145 auto const s_L = current_state.S_L_data.S_L;
1146 auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
1148 auto const& b = this->process_data_.specific_body_force;
1154 MCpG.noalias() += NTN * (ip_cv.fC_4_MCpG.m * w);
1155 MCpC.noalias() += NTN * (ip_cv.fC_4_MCpC.m * w);
1157 if (this->process_data_.apply_mass_lumping)
1159 if (pCap - pCap_prev != 0.)
1162 NTN * (ip_cv.fC_4_MCpC.ml / (pCap - pCap_prev) * w);
1166 MCT.noalias() += NTN * (ip_cv.fC_4_MCT.m * w);
1167 MCu.noalias() += BTI2N.transpose() * (ip_cv.fC_4_MCu.m * w);
1169 LCpG.noalias() += gradNpT * ip_cv.fC_4_LCpG.L * gradNp * w;
1171 LCpC.noalias() += gradNpT * ip_cv.fC_4_LCpC.L * gradNp * w;
1173 LCT.noalias() += gradNpT * ip_cv.fC_4_LCT.L * gradNp * w;
1175 fC.noalias() += gradNpT * ip_cv.fC_1.A * b * w;
1177 if (!this->process_data_.apply_mass_lumping)
1179 fC.noalias() -= NpT * (ip_cv.fC_2a.a * s_L_dot * w);
1183 NpT * (current_state.porosity_data.phi * ip_cv.fC_3a.a * w);
1189 MWpG.noalias() += NTN * (ip_cv.fW_4_MWpG.m * w);
1190 MWpC.noalias() += NTN * (ip_cv.fW_4_MWpC.m * w);
1192 if (this->process_data_.apply_mass_lumping)
1194 if (pCap - pCap_prev != 0.)
1197 NTN * (ip_cv.fW_4_MWpC.ml / (pCap - pCap_prev) * w);
1201 MWT.noalias() += NTN * (ip_cv.fW_4_MWT.m * w);
1203 MWu.noalias() += BTI2N.transpose() * (ip_cv.fW_4_MWu.m * w);
1205 LWpG.noalias() += gradNpT * ip_cv.fW_4_LWpG.L * gradNp * w;
1207 LWpC.noalias() += gradNpT * ip_cv.fW_4_LWpC.L * gradNp * w;
1209 LWT.noalias() += gradNpT * ip_cv.fW_4_LWT.L * gradNp * w;
1211 fW.noalias() += gradNpT * ip_cv.fW_1.A * b * w;
1213 if (!this->process_data_.apply_mass_lumping)
1215 fW.noalias() -= NpT * (ip_cv.fW_2.a * s_L_dot * w);
1219 NpT * (current_state.porosity_data.phi * ip_cv.fW_3a.a * w);
1227 (ip_cv.effective_volumetric_enthalpy_data.rho_h_eff * w);
1230 gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w;
1232 fT.noalias() -= NTT * (ip_cv.fT_1.m * w);
1234 fT.noalias() += gradNTT * ip_cv.fT_2.A * w;
1236 fT.noalias() += gradNTT * ip_cv.fT_3.gradN * w;
1238 fT.noalias() += NTT * (ip_cv.fT_3.N * w);
1244 KUpG.noalias() -= BTI2N * (ip_cv.biot_data() * w);
1246 KUpC.noalias() += BTI2N * (ip_cv.fu_2_KupC.m * w);
1249 Bu.transpose() * ip_cd.s_mech_data.stiffness_tensor * Bu * w;
1252 (Bu.transpose() * current_state.eff_stress_data.sigma_eff -
1253 N_u_op(Nu).transpose() * ip_cv.volumetric_body_force()) *
1256 if (this->process_data_.apply_mass_lumping)
1258 MCpG = MCpG.colwise().sum().eval().asDiagonal();
1259 MCpC = MCpC.colwise().sum().eval().asDiagonal();
1260 MWpG = MWpG.colwise().sum().eval().asDiagonal();
1261 MWpC = MWpC.colwise().sum().eval().asDiagonal();
1272 assembleWithJacobian(
double const t,
double const dt,
1273 std::vector<double>
const& local_x,
1274 std::vector<double>
const& local_x_prev,
1275 std::vector<double>& local_rhs_data,
1276 std::vector<double>& local_Jac_data)
1278 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
1279 temperature_size + displacement_size;
1280 assert(local_x.size() == matrix_size);
1282 auto const temperature = Eigen::Map<VectorType<temperature_size>
const>(
1283 local_x.data() + temperature_index, temperature_size);
1285 auto const gas_pressure = Eigen::Map<VectorType<gas_pressure_size>
const>(
1286 local_x.data() + gas_pressure_index, gas_pressure_size);
1288 auto const capillary_pressure =
1289 Eigen::Map<VectorType<capillary_pressure_size>
const>(
1290 local_x.data() + capillary_pressure_index, capillary_pressure_size);
1292 auto const displacement = Eigen::Map<VectorType<displacement_size>
const>(
1293 local_x.data() + displacement_index, displacement_size);
1295 auto const gas_pressure_prev =
1296 Eigen::Map<VectorType<gas_pressure_size>
const>(
1297 local_x_prev.data() + gas_pressure_index, gas_pressure_size);
1299 auto const capillary_pressure_prev =
1300 Eigen::Map<VectorType<capillary_pressure_size>
const>(
1301 local_x_prev.data() + capillary_pressure_index,
1302 capillary_pressure_size);
1304 auto const temperature_prev =
1305 Eigen::Map<VectorType<temperature_size>
const>(
1306 local_x_prev.data() + temperature_index, temperature_size);
1308 auto const displacement_prev =
1309 Eigen::Map<VectorType<displacement_size>
const>(
1310 local_x_prev.data() + displacement_index, displacement_size);
1314 local_Jac_data, matrix_size, matrix_size);
1317 local_rhs_data, matrix_size);
1328 C_size, capillary_pressure_size);
1338 C_size, capillary_pressure_size);
1347 W_size, capillary_pressure_size);
1358 W_size, capillary_pressure_size);
1365 temperature_size, displacement_size);
1375 displacement_size, gas_pressure_size);
1378 displacement_size, capillary_pressure_size);
1381 auto fC = local_f.template segment<C_size>(C_index);
1383 auto fW = local_f.template segment<W_size>(W_index);
1385 auto fT = local_f.template segment<temperature_size>(temperature_index);
1387 auto fU = local_f.template segment<displacement_size>(displacement_index);
1389 unsigned const n_integration_points =
1390 this->integration_method_.getNumberOfPoints();
1393 this->solid_material_, *this->process_data_.phase_transition_model_};
1395 auto const [ip_constitutive_data, ip_constitutive_variables] =
1396 updateConstitutiveVariables(
1397 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1398 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1399 local_x_prev.size()),
1402 auto const ip_d_data = updateConstitutiveVariablesDerivatives(
1403 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1404 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1405 local_x_prev.size()),
1406 t, dt, ip_constitutive_data, ip_constitutive_variables, models);
1408 for (
unsigned int_point = 0; int_point < n_integration_points; int_point++)
1410 auto& ip = _ip_data[int_point];
1411 auto& ip_cd = ip_constitutive_data[int_point];
1412 auto& ip_dd = ip_d_data[int_point];
1413 auto& ip_cv = ip_constitutive_variables[int_point];
1414 auto& current_state = this->current_states_[int_point];
1415 auto& prev_state = this->prev_states_[int_point];
1417 auto const& Np = ip.N_p;
1418 auto const& NT = Np;
1419 auto const& Nu = ip.N_u;
1421 std::nullopt, this->element_.getID(),
1425 this->element_, Nu))};
1427 auto const& NpT = Np.transpose().eval();
1428 auto const& NTT = NT.transpose().eval();
1430 auto const& gradNp = ip.dNdx_p;
1431 auto const& gradNT = gradNp;
1432 auto const& gradNu = ip.dNdx_u;
1434 auto const& gradNpT = gradNp.transpose().eval();
1435 auto const& gradNTT = gradNT.transpose().eval();
1437 auto const& w = ip.integration_weight;
1439 auto const x_coord =
1442 this->element_, Nu);
1446 ShapeFunctionDisplacement::NPOINTS,
1448 gradNu, Nu, x_coord, this->is_axially_symmetric_);
1450 auto const NTN = (Np.transpose() * Np).eval();
1451 auto const BTI2N = (Bu.transpose() * Invariants::identity2 * Np).eval();
1453 double const div_u_dot =
1454 Invariants::trace(Bu * (displacement - displacement_prev) / dt);
1456 double const pGR = Np.dot(gas_pressure);
1457 double const pCap = Np.dot(capillary_pressure);
1458 double const T = NT.dot(temperature);
1464 double const pGR_prev = Np.dot(gas_pressure_prev);
1465 double const pCap_prev = Np.dot(capillary_pressure_prev);
1466 double const T_prev = NT.dot(temperature_prev);
1468 auto const& s_L = current_state.S_L_data.S_L;
1469 auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
1471 auto const& b = this->process_data_.specific_body_force;
1477 MCpG.noalias() += NTN * (ip_cv.fC_4_MCpG.m * w);
1478 MCpC.noalias() += NTN * (ip_cv.fC_4_MCpC.m * w);
1480 if (this->process_data_.apply_mass_lumping)
1482 if (pCap - pCap_prev != 0.)
1485 NTN * (ip_cv.fC_4_MCpC.ml / (pCap - pCap_prev) * w);
1489 MCT.noalias() += NTN * (ip_cv.fC_4_MCT.m * w);
1492 .template block<C_size, temperature_size>(C_index,
1494 .noalias() += NTN * (ip_dd.dfC_4_MCT.dT * (T - T_prev) / dt * w);
1496 MCu.noalias() += BTI2N.transpose() * (ip_cv.fC_4_MCu.m * w);
1499 .template block<C_size, temperature_size>(C_index,
1501 .noalias() += NTN * (ip_dd.dfC_4_MCu.dT * div_u_dot * w);
1503 LCpG.noalias() += gradNpT * ip_cv.fC_4_LCpG.L * gradNp * w;
1506 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1507 gradNpT * ip_dd.dfC_4_LCpG.dp_GR * gradpGR * Np * w;
1510 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1511 gradNpT * ip_dd.dfC_4_LCpG.dp_cap * gradpGR * Np * w;
1515 .template block<C_size, temperature_size>(C_index,
1517 .noalias() += gradNpT * ip_dd.dfC_4_LCpG.dT * gradpGR * NT * w;
1520 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1521 NTN * (ip_dd.dfC_4_MCpG.dp_GR * (pGR - pGR_prev) / dt * w);
1525 .template block<C_size, temperature_size>(C_index,
1528 NTN * (ip_dd.dfC_4_MCpG.dT * (pGR - pGR_prev) / dt * w);
1530 LCpC.noalias() -= gradNpT * ip_cv.fC_4_LCpC.L * gradNp * w;
1546 LCT.noalias() += gradNpT * ip_cv.fC_4_LCT.L * gradNp * w;
1549 fC.noalias() += gradNpT * ip_cv.fC_1.A * b * w;
1551 if (!this->process_data_.apply_mass_lumping)
1554 fC.noalias() -= NpT * (ip_cv.fC_2a.a * s_L_dot * w);
1556 local_Jac.template block<C_size, C_size>(C_index, C_index)
1558 NTN * ((ip_dd.dfC_2a.dp_GR * s_L_dot
1562 local_Jac.template block<C_size, W_size>(C_index, W_index)
1564 NTN * ((ip_dd.dfC_2a.dp_cap * s_L_dot +
1565 ip_cv.fC_2a.a * ip_dd.dS_L_dp_cap() / dt) *
1569 .template block<C_size, temperature_size>(C_index,
1571 .noalias() += NTN * (ip_dd.dfC_2a.dT * s_L_dot * w);
1576 NpT * (current_state.porosity_data.phi * ip_cv.fC_3a.a * w);
1578 local_Jac.template block<C_size, C_size>(C_index, C_index)
1579 .noalias() += NTN * (current_state.porosity_data.phi *
1580 ip_dd.dfC_3a.dp_GR * w);
1582 local_Jac.template block<C_size, W_size>(C_index, W_index)
1583 .noalias() += NTN * (current_state.porosity_data.phi *
1584 ip_dd.dfC_3a.dp_cap * w);
1587 .template block<C_size, temperature_size>(C_index,
1590 NTN * ((ip_dd.porosity_d_data.dphi_dT * ip_cv.fC_3a.a +
1591 current_state.porosity_data.phi * ip_dd.dfC_3a.dT) *
1598 MWpG.noalias() += NTN * (ip_cv.fW_4_MWpG.m * w);
1599 MWpC.noalias() += NTN * (ip_cv.fW_4_MWpC.m * w);
1601 if (this->process_data_.apply_mass_lumping)
1603 if (pCap - pCap_prev != 0.)
1606 NTN * (ip_cv.fW_4_MWpC.ml / (pCap - pCap_prev) * w);
1610 MWT.noalias() += NTN * (ip_cv.fW_4_MWT.m * w);
1612 MWu.noalias() += BTI2N.transpose() * (ip_cv.fW_4_MWu.m * w);
1614 LWpG.noalias() += gradNpT * ip_cv.fW_4_LWpG.L * gradNp * w;
1617 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1618 gradNpT * ip_dd.dfW_4_LWpG.dp_GR * gradpGR * Np * w;
1620 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1621 gradNpT * ip_dd.dfW_4_LWpG.dp_cap * gradpGR * Np * w;
1624 .template block<W_size, temperature_size>(W_index,
1626 .noalias() += gradNpT * ip_dd.dfW_4_LWpG.dT * gradpGR * NT * w;
1628 LWpC.noalias() += gradNpT * ip_cv.fW_4_LWpC.L * gradNp * w;
1631 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() -=
1632 gradNpT * ip_dd.dfW_4_LWpC.dp_GR * gradpCap * Np * w;
1634 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() -=
1635 gradNpT * ip_dd.dfW_4_LWpC.dp_cap * gradpCap * Np * w;
1638 .template block<W_size, temperature_size>(W_index,
1640 .noalias() -= gradNpT * ip_dd.dfW_4_LWpC.dT * gradpCap * NT * w;
1642 LWT.noalias() += gradNpT * ip_cv.fW_4_LWT.L * gradNp * w;
1645 fW.noalias() += gradNpT * ip_cv.fW_1.A * b * w;
1648 if (!this->process_data_.apply_mass_lumping)
1650 fW.noalias() -= NpT * (ip_cv.fW_2.a * s_L_dot * w);
1652 local_Jac.template block<W_size, C_size>(W_index, C_index)
1653 .noalias() += NTN * (ip_dd.dfW_2.dp_GR * s_L_dot * w);
1658 local_Jac.template block<W_size, W_size>(W_index, W_index)
1659 .noalias() += NTN * ((ip_dd.dfW_2.dp_cap * s_L_dot +
1660 ip_cv.fW_2.a * ip_dd.dS_L_dp_cap() / dt) *
1664 .template block<W_size, temperature_size>(W_index,
1666 .noalias() += NTN * (ip_dd.dfW_2.dT * s_L_dot * w);
1671 NpT * (current_state.porosity_data.phi * ip_cv.fW_3a.a * w);
1673 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1674 NTN * (current_state.porosity_data.phi * ip_dd.dfW_3a.dp_GR * w);
1676 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1677 NTN * (current_state.porosity_data.phi * ip_dd.dfW_3a.dp_cap * w);
1680 .template block<W_size, temperature_size>(W_index,
1683 NTN * ((ip_dd.porosity_d_data.dphi_dT * ip_cv.fW_3a.a +
1684 current_state.porosity_data.phi * ip_dd.dfW_3a.dT) *
1693 (ip_cv.effective_volumetric_enthalpy_data.rho_h_eff * w);
1698 .template block<temperature_size, C_size>(temperature_index,
1701 NTN * (ip_dd.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_GR *
1707 .template block<temperature_size, W_size>(temperature_index,
1711 (ip_dd.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_cap *
1717 .template block<temperature_size, temperature_size>(
1718 temperature_index, temperature_index)
1720 NTN * (ip_dd.effective_volumetric_enthalpy_d_data.drho_h_eff_dT *
1724 gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w;
1742 .template block<temperature_size, W_size>(temperature_index,
1744 .noalias() += gradNTT *
1745 ip_dd.thermal_conductivity_d_data.dlambda_dp_cap *
1750 .template block<temperature_size, temperature_size>(
1751 temperature_index, temperature_index)
1752 .noalias() += gradNTT *
1753 ip_dd.thermal_conductivity_d_data.dlambda_dT * gradT *
1757 fT.noalias() -= NTT * (ip_cv.fT_1.m * w);
1761 .template block<temperature_size, C_size>(temperature_index,
1763 .noalias() += NTN * (ip_dd.dfT_1.dp_GR * w);
1767 .template block<temperature_size, W_size>(temperature_index,
1769 .noalias() += NTN * (ip_dd.dfT_1.dp_cap * w);
1774 .template block<temperature_size, temperature_size>(
1775 temperature_index, temperature_index)
1776 .noalias() += NTN * (ip_dd.dfT_1.dT * w);
1779 fT.noalias() += gradNTT * ip_cv.fT_2.A * w;
1783 .template block<temperature_size, C_size>(temperature_index,
1787 gradNTT * ip_dd.dfT_2.dp_GR_Npart * Np * w +
1789 gradNTT * ip_dd.dfT_2.dp_GR_gradNpart * gradNp * w;
1793 .template block<temperature_size, W_size>(temperature_index,
1797 gradNTT * (-ip_dd.dfT_2.dp_cap_Npart) * Np * w +
1799 gradNTT * (-ip_dd.dfT_2.dp_cap_gradNpart) * gradNp * w;
1803 .template block<temperature_size, temperature_size>(
1804 temperature_index, temperature_index)
1805 .noalias() -= gradNTT * ip_dd.dfT_2.dT * NT * w;
1808 fT.noalias() += NTT * (ip_cv.fT_3.N * w);
1810 fT.noalias() += gradNTT * ip_cv.fT_3.gradN * w;
1816 KUpG.noalias() -= BTI2N * (ip_cv.biot_data() * w);
1821 KUpC.noalias() += BTI2N * (ip_cv.fu_2_KupC.m * w);
1826 .template block<displacement_size, W_size>(displacement_index,
1828 .noalias() += BTI2N * (ip_dd.dfu_2_KupC.dp_cap * w);
1831 .template block<displacement_size, displacement_size>(
1832 displacement_index, displacement_index)
1834 Bu.transpose() * ip_cd.s_mech_data.stiffness_tensor * Bu * w;
1838 (Bu.transpose() * current_state.eff_stress_data.sigma_eff -
1839 N_u_op(Nu).transpose() * ip_cv.volumetric_body_force()) *
1844 .template block<displacement_size, temperature_size>(
1845 displacement_index, temperature_index)
1846 .noalias() -= Bu.transpose() * ip_dd.dfu_1_KuT.dT * NT * w;
1856 if (this->process_data_.apply_mass_lumping)
1858 MCpG = MCpG.colwise().sum().eval().asDiagonal();
1859 MCpC = MCpC.colwise().sum().eval().asDiagonal();
1860 MWpG = MWpG.colwise().sum().eval().asDiagonal();
1861 MWpC = MWpC.colwise().sum().eval().asDiagonal();
1867 fC.noalias() -= LCpG * gas_pressure + LCpC * capillary_pressure +
1869 MCpG * (gas_pressure - gas_pressure_prev) / dt +
1870 MCpC * (capillary_pressure - capillary_pressure_prev) / dt +
1871 MCT * (temperature - temperature_prev) / dt +
1872 MCu * (displacement - displacement_prev) / dt;
1874 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1876 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1879 .template block<C_size, temperature_size>(C_index, temperature_index)
1880 .noalias() += LCT + MCT / dt;
1882 .template block<C_size, displacement_size>(C_index, displacement_index)
1883 .noalias() += MCu / dt;
1887 fW.noalias() -= LWpG * gas_pressure + LWpC * capillary_pressure +
1889 MWpG * (gas_pressure - gas_pressure_prev) / dt +
1890 MWpC * (capillary_pressure - capillary_pressure_prev) / dt +
1891 MWT * (temperature - temperature_prev) / dt +
1892 MWu * (displacement - displacement_prev) / dt;
1894 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1896 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1899 .template block<W_size, temperature_size>(W_index, temperature_index)
1900 .noalias() += LWT + MWT / dt;
1902 .template block<W_size, displacement_size>(W_index, displacement_index)
1903 .noalias() += MWu / dt;
1908 KTT * temperature + MTu * (displacement - displacement_prev) / dt;
1911 .template block<temperature_size, temperature_size>(temperature_index,
1915 .template block<temperature_size, displacement_size>(temperature_index,
1917 .noalias() += MTu / dt;
1921 fU.noalias() -= KUpG * gas_pressure + KUpC * capillary_pressure;
1924 .template block<displacement_size, C_size>(displacement_index, C_index)
1927 .template block<displacement_size, W_size>(displacement_index, W_index)