31template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
33TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
39 bool const is_axially_symmetric,
41 : _process_data(process_data),
42 _integration_method(integration_method),
44 _is_axially_symmetric(is_axially_symmetric)
46 unsigned const n_integration_points =
49 _ip_data.reserve(n_integration_points);
52 auto const shape_matrices_u =
55 DisplacementDim>(e, is_axially_symmetric,
58 auto const shape_matrices_p =
63 auto const& solid_material =
68 for (
unsigned ip = 0; ip < n_integration_points; ip++)
70 _ip_data.emplace_back(solid_material);
72 auto const& sm_u = shape_matrices_u[ip];
73 ip_data.integration_weight =
75 sm_u.integralMeasure * sm_u.detJ;
78 ip_data.dNdx_u = sm_u.dNdx;
80 ip_data.N_p = shape_matrices_p[ip].N;
81 ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
87template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
90 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
91 updateConstitutiveVariables(Eigen::VectorXd
const& local_x,
92 Eigen::VectorXd
const& local_x_prev,
93 double const t,
double const dt)
95 [[maybe_unused]]
auto const matrix_size =
96 gas_pressure_size + capillary_pressure_size + temperature_size +
99 assert(local_x.size() == matrix_size);
101 auto const gas_pressure =
102 local_x.template segment<gas_pressure_size>(gas_pressure_index);
103 auto const capillary_pressure =
104 local_x.template segment<capillary_pressure_size>(
105 capillary_pressure_index);
107 auto const temperature =
108 local_x.template segment<temperature_size>(temperature_index);
109 auto const temperature_prev =
110 local_x_prev.template segment<temperature_size>(temperature_index);
112 auto const displacement =
113 local_x.template segment<displacement_size>(displacement_index);
115 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
116 auto const& gas_phase = medium.phase(
"Gas");
117 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
118 auto const& solid_phase = medium.phase(
"Solid");
120 unsigned const n_integration_points =
121 _integration_method.getNumberOfPoints();
123 std::vector<ConstitutiveVariables<DisplacementDim>>
124 ip_constitutive_variables(n_integration_points);
128 for (
unsigned ip = 0; ip < n_integration_points; ip++)
130 auto& ip_data = _ip_data[ip];
131 auto& ip_cv = ip_constitutive_variables[ip];
133 auto const& Np = ip_data.N_p;
135 auto const& Nu = ip_data.N_u;
136 auto const& gradNu = ip_data.dNdx_u;
137 auto const& gradNp = ip_data.dNdx_p;
139 std::nullopt, _element.getID(), ip,
149 double const T = NT.dot(temperature);
150 double const T_prev = NT.dot(temperature_prev);
151 double const pGR = Np.dot(gas_pressure);
152 double const pCap = Np.dot(capillary_pressure);
153 double const pLR = pGR - pCap;
167 ip_data.computeElasticTangentStiffness(t, pos, dt, T_prev, T);
168 auto const K_S = ip_data.solid_material.getBulkModulus(t, pos, &C_el);
170 ip_data.alpha_B = medium.property(MPL::PropertyType::biot_coefficient)
171 .template value<double>(vars, pos, t, dt);
175 ShapeFunctionDisplacement::NPOINTS,
177 gradNu, Nu, x_coord, _is_axially_symmetric);
179 auto& eps = ip_data.eps;
180 eps.noalias() = Bu * displacement;
186 medium.property(MPL::PropertyType::saturation)
187 .template value<double>(
188 vars, pos, t, std::numeric_limits<double>::quiet_NaN());
193 auto const chi = [&medium, pos, t, dt](
double const s_L)
197 return medium.property(MPL::PropertyType::bishops_effective_stress)
198 .template value<double>(vs, pos, t, dt);
200 ip_cv.chi_s_L = chi(ip_data.s_L);
202 medium.property(MPL::PropertyType::bishops_effective_stress)
203 .template dValue<double>(vars, MPL::Variable::liquid_saturation,
210 auto const sigma_total =
211 (_ip_data[ip].sigma_eff - ip_data.alpha_B *
212 (pGR - ip_cv.chi_s_L * pCap) *
213 Invariants::identity2)
221 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
227 ip_data.k_S = MPL::formEigenTensor<DisplacementDim>(
228 medium.property(MPL::PropertyType::permeability)
229 .value(vars, pos, t, dt));
234 MPL::PropertyType::relative_permeability_nonwetting_phase)
235 .template value<double>(vars, pos, t, dt);
237 auto const dk_rel_G_ds_L =
238 medium[MPL::PropertyType::relative_permeability_nonwetting_phase]
239 .template dValue<double>(vars, MPL::Variable::liquid_saturation,
243 medium.property(MPL::PropertyType::relative_permeability)
244 .template value<double>(vars, pos, t, dt);
246 auto const dk_rel_L_ds_L =
247 medium[MPL::PropertyType::relative_permeability]
248 .template dValue<double>(vars, MPL::Variable::liquid_saturation,
252 ip_data.beta_p_SR = (1. - ip_data.alpha_B) / K_S;
255 if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
257 auto& sigma_sw = ip_data.sigma_sw;
258 auto const& sigma_sw_prev = ip_data.sigma_sw_prev;
260 sigma_sw = sigma_sw_prev;
262 auto const sigma_sw_dot =
263 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
265 solid_phase[MPL::PropertyType::swelling_stress_rate]
266 .value(vars, vars_prev, pos, t, dt)));
267 sigma_sw += sigma_sw_dot * dt;
276 .value(vars, pos, t, dt)));
279 ip_data.beta_T_SR = Invariants::trace(ip_data.alpha_T_SR);
282 dthermal_strain = ip_data.alpha_T_SR * (T - T_prev);
284 auto& eps_prev = ip_data.eps_prev;
285 auto& eps_m = ip_data.eps_m;
286 auto& eps_m_prev = ip_data.eps_m_prev;
288 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
290 if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
293 C_el.inverse() * (ip_data.sigma_sw - ip_data.sigma_sw_prev);
300 auto const rho_ref_SR =
301 solid_phase.property(MPL::PropertyType::density)
302 .template value<double>(
303 vars, pos, t, std::numeric_limits<double>::quiet_NaN());
305 double const T0 = _process_data.reference_temperature(t, pos)[0];
306 double const delta_T(T - T0);
307 ip_data.thermal_volume_strain = ip_data.beta_T_SR * delta_T;
310 auto const phi_0 = medium.property(MPL::PropertyType::porosity)
311 .template value<double>(vars, pos, t, dt);
313 auto const phi_S_0 = 1. - phi_0;
315#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
316 auto const& m = Invariants::identity2;
317 double const div_u = m.transpose() * eps;
319 const double phi_S = phi_S_0 * (1. + ip_data.thermal_volume_strain -
320 ip_data.alpha_B * div_u);
322 const double phi_S = phi_S_0;
326 ip_data.phi = 1. - phi_S;
330#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
331 auto const rhoSR = rho_ref_SR * (1. - ip_data.thermal_volume_strain +
332 (ip_data.alpha_B - 1.) * div_u);
334 auto const rhoSR = rho_ref_SR;
337 ip_cv.C = ip_data.updateConstitutiveRelation(vars, t, pos, dt, T_prev);
340 auto& ptm = *_process_data.phase_transition_model_;
341 ptmv = ptm.updateConstitutiveVariables(ptmv, &medium, vars, pos, t, dt);
342 auto const& c = ptmv;
344 auto const phi_L = ip_data.s_L * ip_data.phi;
345 auto const phi_G = (1. - ip_data.s_L) * ip_data.phi;
348 ip_data.lambda = MaterialPropertyLib::formEigenTensor<DisplacementDim>(
352 .value(vars, pos, t, dt));
355 solid_phase.property(MPL::PropertyType::specific_heat_capacity)
356 .template value<double>(vars, pos, t, dt);
357 ip_data.h_S = cpS * T;
358 auto const u_S = ip_data.h_S;
360 ip_data.rho_u_eff = phi_G * c.rhoGR * c.uG + phi_L * c.rhoLR * c.uL +
363 ip_data.rho_G_h_G = phi_G * c.rhoGR * c.hG;
364 ip_data.rho_L_h_L = phi_L * c.rhoLR * c.hL;
365 ip_data.rho_S_h_S = phi_S * rhoSR * ip_data.h_S;
367 ip_data.muGR = c.muGR;
368 ip_data.muLR = c.muLR;
370 ip_data.rhoGR = c.rhoGR;
371 ip_data.rhoLR = c.rhoLR;
372 ip_data.rhoSR = rhoSR;
374 ip_data.rhoCGR = c.rhoCGR;
375 ip_data.rhoCLR = c.rhoCLR;
376 ip_data.rhoWGR = c.rhoWGR;
377 ip_data.rhoWLR = c.rhoWLR;
379 ip_data.dxmWG_dpGR = c.dxmWG_dpGR;
380 ip_data.dxmWG_dpCap = c.dxmWG_dpCap;
381 ip_data.dxmWG_dT = c.dxmWG_dT;
383 ip_data.dxmWL_dpGR = c.dxmWL_dpGR;
384 ip_data.dxmWL_dpCap = c.dxmWL_dpCap;
385 ip_data.dxmWL_dT = c.dxmWL_dT;
387 ip_data.dxmWL_dpLR = c.dxmWL_dpLR;
390 ip_data.xnCG = 1. - c.xnWG;
391 ip_data.xmCG = 1. - c.xmWG;
392 ip_data.xmWG = c.xmWG;
393 ip_data.xmWL = c.xmWL;
394 auto const xmCL = 1. - c.xmWL;
396 ip_data.diffusion_coefficient_vapour = c.diffusion_coefficient_vapour;
397 ip_data.diffusion_coefficient_solute = c.diffusion_coefficient_solute;
400 ip_data.h_CG = c.hCG;
401 ip_data.h_WG = c.hWG;
403 ip_data.pWGR = c.pWGR;
406 ip_data.dxmWG_dpCap * gradpCap +
407 ip_data.dxmWG_dT * gradT;
411 ip_data.dxmWL_dpCap * gradpCap +
412 ip_data.dxmWL_dT * gradT;
419 ip_data.d_CG = ip_data.xmCG == 0.
422 : -phi_G / ip_data.xmCG *
423 ip_data.diffusion_coefficient_vapour *
426 ip_data.d_WG = ip_data.xmWG == 0.
429 : -phi_G / ip_data.xmWG *
430 ip_data.diffusion_coefficient_vapour *
433 ip_data.d_CL = xmCL == 0. ? 0. * gradxmCL
436 ip_data.diffusion_coefficient_solute *
439 ip_data.d_WL = ip_data.xmWL == 0.
442 : -phi_L / ip_data.xmWL *
443 ip_data.diffusion_coefficient_solute *
449 auto const drho_LR_dT =
450 liquid_phase.property(MPL::PropertyType::density)
451 .template dValue<double>(vars, MPL::Variable::temperature, pos,
453 auto const drho_SR_dT =
454 solid_phase.property(MPL::PropertyType::density)
455 .template dValue<double>(vars, MPL::Variable::temperature,
457#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
458 * (1. - ip_data.thermal_volume_strain +
459 (ip_data.alpha_B - 1.) * div_u) -
460 rho_ref_SR * ip_data.beta_T_SR
465 auto const dphi_0_dT =
466 medium[MPL::PropertyType::porosity].template dValue<double>(
467 vars, MPL::Variable::temperature, pos, t, dt);
469 auto const dphi_S_0_dT = -dphi_0_dT;
470 const double dphi_S_dT = dphi_S_0_dT
471#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
472 * (1. + ip_data.thermal_volume_strain -
473 ip_data.alpha_B * div_u) +
474 phi_S_0 * ip_data.beta_T_SR
478 ip_cv.drho_u_eff_dT =
479 phi_G * c.drho_GR_dT * c.uG + phi_G * c.rhoGR * c.du_G_dT +
480 phi_L * drho_LR_dT * c.uL + phi_L * c.rhoLR * c.du_L_dT +
481 phi_S * drho_SR_dT * u_S + phi_S * rhoSR * cpS +
482 dphi_S_dT * rhoSR * u_S;
485 medium[MPL::PropertyType::saturation].template dValue<double>(
486 vars, MPL::Variable::capillary_pressure, pos, t, dt);
490 double const ds_G_dp_cap = -ip_cv.ds_L_dp_cap;
493 double const dphi_G_dp_cap = -ip_cv.ds_L_dp_cap * ip_data.phi;
495 double const dphi_L_dp_cap = ip_cv.ds_L_dp_cap * ip_data.phi;
497 auto const lambdaGR =
498 gas_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
499 ? MPL::formEigenTensor<DisplacementDim>(
501 .property(MPL::PropertyType::thermal_conductivity)
502 .value(vars, pos, t, dt))
503 : MPL::formEigenTensor<DisplacementDim>(0.);
505 auto const dlambda_GR_dT =
506 gas_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
507 ? MPL::formEigenTensor<DisplacementDim>(
508 gas_phase[MPL::PropertyType::thermal_conductivity].dValue(
509 vars, MPL::Variable::temperature, pos, t, dt))
510 : MPL::formEigenTensor<DisplacementDim>(0.);
512 auto const lambdaLR =
513 liquid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
514 ? MPL::formEigenTensor<DisplacementDim>(
516 .property(MPL::PropertyType::thermal_conductivity)
517 .value(vars, pos, t, dt))
518 : MPL::formEigenTensor<DisplacementDim>(0.);
520 auto const dlambda_LR_dT =
521 liquid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
522 ? MPL::formEigenTensor<DisplacementDim>(
523 liquid_phase[MPL::PropertyType::thermal_conductivity]
524 .dValue(vars, MPL::Variable::temperature, pos, t, dt))
525 : MPL::formEigenTensor<DisplacementDim>(0.);
527 auto const lambdaSR =
528 solid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
529 ? MPL::formEigenTensor<DisplacementDim>(
531 .property(MPL::PropertyType::thermal_conductivity)
532 .value(vars, pos, t, dt))
533 : MPL::formEigenTensor<DisplacementDim>(0.);
535 auto const dlambda_SR_dT =
536 solid_phase.hasProperty(MPL::PropertyType::thermal_conductivity)
537 ? MPL::formEigenTensor<DisplacementDim>(
538 solid_phase[MPL::PropertyType::thermal_conductivity]
539 .dValue(vars, MPL::Variable::temperature, pos, t, dt))
540 : MPL::formEigenTensor<DisplacementDim>(0.);
542 ip_cv.dlambda_dp_cap =
543 dphi_G_dp_cap * lambdaGR + dphi_L_dp_cap * lambdaLR;
545 ip_cv.dlambda_dT = phi_G * dlambda_GR_dT + phi_L * dlambda_LR_dT +
546 phi_S * dlambda_SR_dT + dphi_S_dT * lambdaSR;
552 double const drho_LR_dp_GR = c.drho_LR_dp_LR;
553 double const drho_LR_dp_cap = -c.drho_LR_dp_LR;
556 ip_cv.drho_h_eff_dp_GR =
557 phi_G * c.drho_GR_dp_GR *
559 phi_L * drho_LR_dp_GR *
561 ip_cv.drho_h_eff_dp_cap = dphi_G_dp_cap * c.rhoGR * c.hG +
563 dphi_L_dp_cap * c.rhoLR * c.hL +
564 phi_L * drho_LR_dp_cap * c.hL;
567 constexpr double dphi_G_dT = 0;
568 constexpr double dphi_L_dT = 0;
569 ip_cv.drho_h_eff_dT =
570 dphi_G_dT * c.rhoGR * c.hG + phi_G * c.drho_GR_dT * c.hG +
571 phi_G * c.rhoGR * c.dh_G_dT + dphi_L_dT * c.rhoLR * c.hL +
572 phi_L * drho_LR_dT * c.hL + phi_L * c.rhoLR * c.dh_L_dT +
573 dphi_S_dT * rhoSR * ip_data.h_S + phi_S * drho_SR_dT * ip_data.h_S +
576 ip_cv.drho_u_eff_dp_GR =
578 phi_G * c.drho_GR_dp_GR * c.uG + phi_G * c.rhoGR * c.du_G_dp_GR +
580 phi_L * drho_LR_dp_GR * c.uL + phi_L * c.rhoLR * c.du_L_dp_GR;
582 ip_cv.drho_u_eff_dp_cap = dphi_G_dp_cap * c.rhoGR * c.uG +
584 dphi_L_dp_cap * c.rhoLR * c.uL +
585 phi_L * drho_LR_dp_cap * c.uL +
586 phi_L * c.rhoLR * c.du_L_dp_cap;
588 auto const& b = _process_data.specific_body_force;
590 ip_data.k_S * ip_data.k_rel_G / ip_data.muGR;
592 ip_data.k_S * ip_data.k_rel_L / ip_data.muLR;
600 ip_cv.dk_over_mu_G_dp_cap =
601 ip_data.k_S * dk_rel_G_ds_L * ip_cv.ds_L_dp_cap / ip_data.muGR;
602 ip_cv.dk_over_mu_L_dp_cap =
603 ip_data.k_S * dk_rel_L_ds_L * ip_cv.ds_L_dp_cap / ip_data.muLR;
605 ip_data.w_GS = k_over_mu_G * c.rhoGR * b - k_over_mu_G * gradpGR;
606 ip_data.w_LS = k_over_mu_L * gradpCap + k_over_mu_L * c.rhoLR * b -
607 k_over_mu_L * gradpGR;
609 ip_cv.drho_GR_h_w_eff_dp_GR_Npart =
610 c.drho_GR_dp_GR * c.hG * ip_data.w_GS +
611 c.rhoGR * c.hG * k_over_mu_G * c.drho_GR_dp_GR * b;
612 ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart =
613 -c.rhoGR * c.hG * k_over_mu_G - c.rhoLR * c.hL * k_over_mu_L;
615 ip_cv.drho_LR_h_w_eff_dp_cap_Npart =
616 -drho_LR_dp_cap * c.hL * ip_data.w_LS -
617 c.rhoLR * c.hL * k_over_mu_L * drho_LR_dp_cap * b;
618 ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart =
620 -c.rhoLR * c.hL * k_over_mu_L;
622 ip_cv.drho_GR_h_w_eff_dT = c.drho_GR_dT * c.hG * ip_data.w_GS +
623 c.rhoGR * c.dh_G_dT * ip_data.w_GS +
624 drho_LR_dT * c.hL * ip_data.w_LS +
625 c.rhoLR * c.dh_L_dT * ip_data.w_LS;
631 double const s_L = ip_data.s_L;
632 double const s_G = 1. - ip_data.s_L;
633 double const rho_C_FR = s_G * ip_data.rhoCGR + s_L * ip_data.rhoCLR;
634 double const rho_W_FR = s_G * ip_data.rhoWGR + s_L * ip_data.rhoWLR;
636 constexpr double drho_C_GR_dp_cap = 0;
639 ip_cv.dfC_3a_dp_GR = 0.;
640 ip_cv.dfC_3a_dp_cap = 0.;
641 ip_cv.dfC_3a_dT = 0.;
645 double const rho_C_GR_dot =
646 (ip_data.rhoCGR - ip_data.rhoCGR_prev) / dt;
647 double const rho_C_LR_dot =
648 (ip_data.rhoCLR - ip_data.rhoCLR_prev) / dt;
650 s_G * c.drho_C_GR_dp_GR /
652 s_L * c.drho_C_LR_dp_GR /
654 ip_cv.dfC_3a_dp_cap =
655 ds_G_dp_cap * rho_C_GR_dot + s_G * drho_C_GR_dp_cap / dt +
656 ip_cv.ds_L_dp_cap * rho_C_LR_dot - s_L * c.drho_C_LR_dp_LR / dt;
658 s_G * c.drho_C_GR_dT / dt + s_L * c.drho_C_LR_dT / dt;
661 double const drho_C_FR_dp_GR =
662 s_G * c.drho_C_GR_dp_GR +
663 s_L * c.drho_C_LR_dp_GR;
664 ip_cv.dfC_4_MCpG_dp_GR = drho_C_FR_dp_GR *
665 (ip_data.alpha_B - ip_data.phi) *
668 double const drho_C_FR_dT = s_G * c.drho_C_GR_dT + s_L * c.drho_C_LR_dT;
669 ip_cv.dfC_4_MCpG_dT =
670 drho_C_FR_dT * (ip_data.alpha_B - ip_data.phi) * ip_data.beta_p_SR
671#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
672 - rho_C_FR * ip_data.dphi_dT * ip_data.beta_p_SR
677 drho_C_FR_dT * (ip_data.alpha_B - ip_data.phi) * ip_data.beta_T_SR
678#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
679 + rho_C_FR * (ip_data.alpha_B - ip_data.dphi_dT) * ip_data.beta_T_SR
683 ip_cv.dfC_4_MCu_dT = drho_C_FR_dT * ip_data.alpha_B;
685 ip_cv.dfC_2a_dp_GR = -ip_data.phi * c.drho_C_GR_dp_GR -
686 drho_C_FR_dp_GR * pCap *
687 (ip_data.alpha_B - ip_data.phi) *
690 double const drho_C_FR_dp_cap =
691 ds_G_dp_cap * ip_data.rhoCGR + s_G * drho_C_GR_dp_cap +
692 ip_cv.ds_L_dp_cap * ip_data.rhoCLR - s_L * c.drho_C_LR_dp_LR;
694 ip_cv.dfC_2a_dp_cap =
695 ip_data.phi * (-c.drho_C_LR_dp_LR - drho_C_GR_dp_cap) -
696 drho_C_FR_dp_cap * pCap * (ip_data.alpha_B - ip_data.phi) *
698 rho_C_FR * (ip_data.alpha_B - ip_data.phi) * ip_data.beta_p_SR;
701#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
702 ip_data.dphi_dT * (ip_data.rhoCLR - ip_data.rhoCGR) +
704 ip_data.phi * (c.drho_C_LR_dT - c.drho_C_GR_dT) -
705 drho_C_FR_dT * pCap * (ip_data.alpha_B - ip_data.phi) *
707#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
708 + rho_C_FR * pCap * ip_data.dphi_dT * ip_data.beta_p_SR
712 ip_cv.dadvection_C_dp_GR = c.drho_C_GR_dp_GR * k_over_mu_G
715 + c.drho_C_LR_dp_GR * k_over_mu_L;
717 ip_cv.dadvection_C_dp_cap =
719 ip_data.rhoCGR * ip_cv.dk_over_mu_G_dp_cap +
720 (-c.drho_C_LR_dp_LR) * k_over_mu_L +
721 ip_data.rhoCLR * ip_cv.dk_over_mu_L_dp_cap;
723 ip_cv.dfC_4_LCpG_dT =
724 c.drho_C_GR_dT * k_over_mu_G + c.drho_C_LR_dT * k_over_mu_L
728 double const drho_W_FR_dp_GR =
729 s_G * c.drho_W_GR_dp_GR +
730 s_L * c.drho_W_LR_dp_GR;
731 double const drho_W_FR_dp_cap =
732 ds_G_dp_cap * ip_data.rhoWGR + s_G * c.drho_W_GR_dp_cap +
733 ip_cv.ds_L_dp_cap * ip_data.rhoWLR - s_L * c.drho_W_LR_dp_LR;
734 double const drho_W_FR_dT = s_G * c.drho_W_GR_dT + s_L * c.drho_W_LR_dT;
737 ip_data.phi * (c.drho_W_LR_dp_GR - c.drho_W_GR_dp_GR);
738 ip_cv.dfW_2b_dp_GR = drho_W_FR_dp_GR * pCap *
739 (ip_data.alpha_B - ip_data.phi) *
741 ip_cv.dfW_2a_dp_cap =
742 ip_data.phi * (-c.drho_W_LR_dp_LR - c.drho_W_GR_dp_cap);
743 ip_cv.dfW_2b_dp_cap =
744 drho_W_FR_dp_cap * pCap * (ip_data.alpha_B - ip_data.phi) *
746 rho_W_FR * (ip_data.alpha_B - ip_data.phi) * ip_data.beta_p_SR;
749#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
750 ip_data.dphi_dT * (ip_data.rhoWLR - ip_data.rhoWGR) +
752 ip_data.phi * (c.drho_W_LR_dT - c.drho_W_GR_dT);
754 drho_W_FR_dT * pCap * (ip_data.alpha_B - ip_data.phi) *
756#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
757 - rho_W_FR * pCap * ip_data.dphi_dT * ip_data.beta_p_SR
763 ip_cv.dfW_3a_dp_GR = 0.;
764 ip_cv.dfW_3a_dp_cap = 0.;
765 ip_cv.dfW_3a_dT = 0.;
769 double const rho_W_GR_dot =
770 (ip_data.rhoWGR - ip_data.rhoWGR_prev) / dt;
771 double const rho_W_LR_dot =
772 (ip_data.rhoWLR - ip_data.rhoWLR_prev) / dt;
775 s_G * c.drho_W_GR_dp_GR /
777 s_L * c.drho_W_LR_dp_GR /
779 ip_cv.dfW_3a_dp_cap =
780 ds_G_dp_cap * rho_W_GR_dot + s_G * c.drho_W_GR_dp_cap / dt +
781 ip_cv.ds_L_dp_cap * rho_W_LR_dot - s_L * c.drho_W_LR_dp_LR / dt;
783 s_G * c.drho_W_GR_dT / dt + s_L * c.drho_W_LR_dT / dt;
786 ip_cv.dfW_4_LWpG_a_dp_GR = c.drho_W_GR_dp_GR * k_over_mu_G
788 + c.drho_W_LR_dp_GR * k_over_mu_L
791 ip_cv.dfW_4_LWpG_a_dp_cap = c.drho_W_GR_dp_cap * k_over_mu_G +
792 ip_data.rhoWGR * ip_cv.dk_over_mu_G_dp_cap +
793 -c.drho_W_LR_dp_LR * k_over_mu_L +
794 ip_data.rhoWLR * ip_cv.dk_over_mu_L_dp_cap;
796 ip_cv.dfW_4_LWpG_a_dT =
797 c.drho_W_GR_dT * k_over_mu_G
799 + c.drho_W_LR_dT * k_over_mu_L
804 ip_cv.dfW_4_LWpG_d_dp_GR =
805 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
806 ip_cv.dfW_4_LWpG_d_dp_cap =
807 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
808 ip_cv.dfW_4_LWpG_d_dT =
809 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
811 ip_cv.dfW_4_LWpC_a_dp_GR = c.drho_W_LR_dp_GR * k_over_mu_L
814 ip_cv.dfW_4_LWpC_a_dp_cap = -c.drho_W_LR_dp_LR * k_over_mu_L +
815 ip_data.rhoWLR * ip_cv.dk_over_mu_L_dp_cap;
816 ip_cv.dfW_4_LWpC_a_dT = c.drho_W_LR_dT * k_over_mu_L
821 ip_cv.dfW_4_LWpC_d_dp_GR =
822 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
823 ip_cv.dfW_4_LWpC_d_dp_cap =
824 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
825 ip_cv.dfW_4_LWpC_d_dT =
826 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
828 ip_cv.dfC_4_LCpC_a_dp_GR = c.drho_C_LR_dp_GR * k_over_mu_L
831 ip_cv.dfC_4_LCpC_a_dp_cap = -c.drho_C_LR_dp_LR * k_over_mu_L +
832 ip_data.rhoCLR * ip_cv.dk_over_mu_L_dp_cap;
833 ip_cv.dfC_4_LCpC_a_dT = c.drho_W_LR_dT * k_over_mu_L
838 ip_cv.dfC_4_LCpC_d_dp_GR =
839 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
840 ip_cv.dfC_4_LCpC_d_dp_cap =
841 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
842 ip_cv.dfC_4_LCpC_d_dT =
843 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Zero();
846 return ip_constitutive_variables;
849template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
852 ShapeFunctionDisplacement, ShapeFunctionPressure,
854 double const* values,
855 int const integration_order)
857 if (integration_order !=
858 static_cast<int>(_integration_method.getIntegrationOrder()))
861 "Setting integration point initial conditions; The integration "
862 "order of the local assembler for element {:d} is different "
863 "from the integration order in the initial condition.",
867 if (name ==
"sigma_ip")
869 if (_process_data.initial_stress !=
nullptr)
872 "Setting initial conditions for stress from integration "
873 "point data and from a parameter '{:s}' is not possible "
875 _process_data.initial_stress->name);
877 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
878 values, _ip_data, &IpData::sigma_eff);
881 if (name ==
"saturation_ip")
886 if (name ==
"swelling_stress_ip")
888 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
889 values, _ip_data, &IpData::sigma_sw);
891 if (name ==
"epsilon_ip")
893 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
894 values, _ip_data, &IpData::eps);
896 if (name.starts_with(
"material_state_variable_") && name.ends_with(
"_ip"))
898 std::string
const variable_name = name.substr(24, name.size() - 24 - 3);
899 DBUG(
"Setting material state variable '{:s}'", variable_name);
903 auto const& internal_variables =
904 _ip_data[0].solid_material.getInternalVariables();
906 std::find_if(begin(internal_variables), end(internal_variables),
907 [&variable_name](
auto const& iv)
908 {
return iv.name == variable_name; });
909 iv != end(internal_variables))
912 values, _ip_data, &IpData::material_state_variables,
919template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
923 setInitialConditionsConcrete(std::vector<double>
const& local_x_data,
928 [[maybe_unused]]
auto const matrix_size =
929 gas_pressure_size + capillary_pressure_size + temperature_size +
932 assert(local_x_data.size() == matrix_size);
934 auto const local_x = Eigen::Map<Eigen::VectorXd const>(local_x_data.data(),
935 local_x_data.size());
936 auto const capillary_pressure =
937 local_x.template segment<capillary_pressure_size>(
938 capillary_pressure_index);
940 auto const temperature =
941 local_x.template segment<temperature_size>(temperature_index);
943 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
944 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
945 auto const& solid_phase = medium.phase(
"Solid");
947 unsigned const n_integration_points =
948 _integration_method.getNumberOfPoints();
950 for (
unsigned ip = 0; ip < n_integration_points; ip++)
954 auto& ip_data = _ip_data[ip];
955 auto const& Np = ip_data.N_p;
958 std::nullopt, _element.getID(), ip,
962 _element, ip_data.N_u))};
964 double const pCap = Np.dot(capillary_pressure);
967 double const T = NT.dot(temperature);
970 auto& eps = ip_data.eps;
976 medium.property(MPL::PropertyType::saturation)
977 .template value<double>(
978 vars, pos, t, std::numeric_limits<double>::quiet_NaN());
983 ip_data.computeElasticTangentStiffness(t, pos, dt, T, T);
984 auto& sigma_sw = ip_data.sigma_sw;
985 ip_data.eps_m_prev.noalias() =
986 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
987 ? eps + C_el.inverse() * sigma_sw
991 updateConstitutiveVariables(local_x, local_x, t, 0);
993 for (
unsigned ip = 0; ip < n_integration_points; ip++)
995 auto& ip_data = _ip_data[ip];
996 ip_data.pushBackState();
1000template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
1001 int DisplacementDim>
1003 ShapeFunctionDisplacement, ShapeFunctionPressure,
1004 DisplacementDim>::assemble(
double const t,
double const dt,
1005 std::vector<double>
const& local_x,
1006 std::vector<double>
const& local_x_prev,
1007 std::vector<double>& local_M_data,
1008 std::vector<double>& local_K_data,
1009 std::vector<double>& local_rhs_data)
1011 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
1012 temperature_size + displacement_size;
1013 assert(local_x.size() == matrix_size);
1015 auto const capillary_pressure =
1016 Eigen::Map<VectorType<capillary_pressure_size>
const>(
1017 local_x.data() + capillary_pressure_index, capillary_pressure_size);
1019 auto const capillary_pressure_prev =
1020 Eigen::Map<VectorType<capillary_pressure_size>
const>(
1021 local_x_prev.data() + capillary_pressure_index,
1022 capillary_pressure_size);
1026 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1027 local_M_data, matrix_size, matrix_size);
1031 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1032 local_K_data, matrix_size, matrix_size);
1035 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
1036 local_rhs_data, matrix_size);
1042 auto MCpG = local_M.template block<C_size, gas_pressure_size>(
1043 C_index, gas_pressure_index);
1044 auto MCpC = local_M.template block<C_size, capillary_pressure_size>(
1045 C_index, capillary_pressure_index);
1046 auto MCT = local_M.template block<C_size, temperature_size>(
1047 C_index, temperature_index);
1048 auto MCu = local_M.template block<C_size, displacement_size>(
1049 C_index, displacement_index);
1052 auto LCpG = local_K.template block<C_size, gas_pressure_size>(
1053 C_index, gas_pressure_index);
1054 auto LCpC = local_K.template block<C_size, capillary_pressure_size>(
1055 C_index, capillary_pressure_index);
1056 auto LCT = local_K.template block<C_size, temperature_size>(
1057 C_index, temperature_index);
1060 auto MWpG = local_M.template block<W_size, gas_pressure_size>(
1061 W_index, gas_pressure_index);
1062 auto MWpC = local_M.template block<W_size, capillary_pressure_size>(
1063 W_index, capillary_pressure_index);
1064 auto MWT = local_M.template block<W_size, temperature_size>(
1065 W_index, temperature_index);
1066 auto MWu = local_M.template block<W_size, displacement_size>(
1067 W_index, displacement_index);
1070 auto LWpG = local_K.template block<W_size, gas_pressure_size>(
1071 W_index, gas_pressure_index);
1072 auto LWpC = local_K.template block<W_size, capillary_pressure_size>(
1073 W_index, capillary_pressure_index);
1074 auto LWT = local_K.template block<W_size, temperature_size>(
1075 W_index, temperature_index);
1078 auto MTu = local_M.template block<temperature_size, displacement_size>(
1079 temperature_index, displacement_index);
1082 auto KTT = local_K.template block<temperature_size, temperature_size>(
1083 temperature_index, temperature_index);
1086 auto KUpG = local_K.template block<displacement_size, gas_pressure_size>(
1087 displacement_index, gas_pressure_index);
1089 local_K.template block<displacement_size, capillary_pressure_size>(
1090 displacement_index, capillary_pressure_index);
1093 auto fC = local_f.template segment<C_size>(C_index);
1095 auto fW = local_f.template segment<W_size>(W_index);
1097 auto fT = local_f.template segment<temperature_size>(temperature_index);
1099 auto fU = local_f.template segment<displacement_size>(displacement_index);
1101 unsigned const n_integration_points =
1102 _integration_method.getNumberOfPoints();
1104 auto const ip_constitutive_variables = updateConstitutiveVariables(
1105 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1106 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1107 local_x_prev.size()),
1110 for (
unsigned int_point = 0; int_point < n_integration_points; int_point++)
1112 auto& ip = _ip_data[int_point];
1113 auto& ip_cv = ip_constitutive_variables[int_point];
1115 auto const& Np = ip.N_p;
1116 auto const& NT = Np;
1117 auto const& Nu = ip.N_u;
1119 std::nullopt, _element.getID(), int_point,
1125 auto const& NpT = Np.transpose().eval();
1126 auto const& NTT = NT.transpose().eval();
1128 auto const& gradNp = ip.dNdx_p;
1129 auto const& gradNT = gradNp;
1130 auto const& gradNu = ip.dNdx_u;
1132 auto const& gradNpT = gradNp.transpose().eval();
1133 auto const& gradNTT = gradNT.transpose().eval();
1135 auto const& w = ip.integration_weight;
1137 auto const& m = Invariants::identity2;
1139 auto const mT = m.transpose().eval();
1141 auto const x_coord =
1148 ShapeFunctionDisplacement::NPOINTS,
1150 gradNu, Nu, x_coord, _is_axially_symmetric);
1152 auto const BuT = Bu.transpose().eval();
1154 double const pCap = Np.dot(capillary_pressure);
1155 double const pCap_prev = Np.dot(capillary_pressure_prev);
1157 auto& beta_T_SR = ip.beta_T_SR;
1160 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
1162 const double sD_G = ip.diffusion_coefficient_vapour;
1163 const double sD_L = ip.diffusion_coefficient_solute;
1173 auto const s_G = 1. - s_L;
1174 auto const s_L_dot = (s_L - ip.s_L_prev) / dt;
1176 auto& alpha_B = ip.alpha_B;
1177 auto& beta_p_SR = ip.beta_p_SR;
1179 auto const& b = _process_data.specific_body_force;
1185 auto const phi_G = s_G * phi;
1186 auto const phi_L = s_L * phi;
1187 auto const phi_S = 1. - phi;
1190 auto& rho_SR = ip.rhoSR;
1192 auto const rho = phi_G * ip.rhoGR + phi_L * ip.rhoLR + phi_S * rho_SR;
1195 const double rho_C_FR = s_G * ip.rhoCGR + s_L * ip.rhoCLR;
1196 const double rho_W_FR = s_G * ip.rhoWGR + s_L * ip.rhoWLR;
1202 auto const rho_C_GR_dot = (ip.rhoCGR - ip.rhoCGR_prev) / dt;
1203 auto const rho_C_LR_dot = (ip.rhoCLR - ip.rhoCLR_prev) / dt;
1204 auto const rho_W_GR_dot = (ip.rhoWGR - ip.rhoWGR_prev) / dt;
1205 auto const rho_W_LR_dot = (ip.rhoWLR - ip.rhoWLR_prev) / dt;
1207 auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
1209 auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / dt;
1218 MCpG.noalias() += NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1220 NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1222 if (_process_data.apply_mass_lumping)
1224 if (pCap - pCap_prev != 0.)
1228 (phi * (ip.rhoCLR - ip.rhoCGR) -
1229 rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1230 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1234 MCT.noalias() -= NpT * rho_C_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1235 MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
1237 using DisplacementDimMatrix =
1238 Eigen::Matrix<double, DisplacementDim, DisplacementDim>;
1240 DisplacementDimMatrix
const advection_C_G = ip.rhoCGR * k_over_mu_G;
1241 DisplacementDimMatrix
const advection_C_L = ip.rhoCLR * k_over_mu_L;
1243 DisplacementDimMatrix
const diffusion_CGpGR =
1244 -phi_G * ip.rhoGR * D_C_G * ip.dxmWG_dpGR;
1245 DisplacementDimMatrix
const diffusion_CLpGR =
1246 -phi_L * ip.rhoLR * D_C_L * ip.dxmWL_dpGR;
1248 DisplacementDimMatrix
const diffusion_CGpCap =
1249 -phi_G * ip.rhoGR * D_C_G * ip.dxmWG_dpCap;
1250 DisplacementDimMatrix
const diffusion_CLpCap =
1251 -phi_L * ip.rhoLR * D_C_L * ip.dxmWL_dpCap;
1253 DisplacementDimMatrix
const diffusion_CGT =
1254 -phi_G * ip.rhoGR * D_C_G * ip.dxmWG_dT;
1255 DisplacementDimMatrix
const diffusion_CLT =
1256 -phi_L * ip.rhoLR * D_C_L * ip.dxmWL_dT;
1258 DisplacementDimMatrix
const advection_C = advection_C_G + advection_C_L;
1259 DisplacementDimMatrix
const diffusion_C_pGR =
1260 diffusion_CGpGR + diffusion_CLpGR;
1261 DisplacementDimMatrix
const diffusion_C_pCap =
1262 diffusion_CGpCap + diffusion_CLpCap;
1264 DisplacementDimMatrix
const diffusion_C_T =
1265 diffusion_CGT + diffusion_CLT;
1268 gradNpT * (advection_C + diffusion_C_pGR) * gradNp * w;
1271 gradNpT * (diffusion_C_pCap - advection_C_L) * gradNp * w;
1273 LCT.noalias() += gradNpT * (diffusion_C_T)*gradNp * w;
1275 fC.noalias() += gradNpT *
1276 (advection_C_G * ip.rhoGR + advection_C_L * ip.rhoLR) *
1279 if (!_process_data.apply_mass_lumping)
1281 fC.noalias() -= NpT *
1282 (phi * (ip.rhoCLR - ip.rhoCGR) -
1283 rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1288 NpT * phi * (s_G * rho_C_GR_dot + s_L * rho_C_LR_dot) * w;
1294 MWpG.noalias() += NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1296 NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1298 if (_process_data.apply_mass_lumping)
1300 if (pCap - pCap_prev != 0.)
1304 (phi * (ip.rhoWLR - ip.rhoWGR) -
1305 rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1306 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1310 MWT.noalias() -= NpT * rho_W_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1312 MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
1314 DisplacementDimMatrix
const advection_W_G = ip.rhoWGR * k_over_mu_G;
1315 DisplacementDimMatrix
const advection_W_L = ip.rhoWLR * k_over_mu_L;
1317 DisplacementDimMatrix
const diffusion_WGpGR =
1318 phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dpGR;
1319 DisplacementDimMatrix
const diffusion_WLpGR =
1320 phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dpGR;
1322 DisplacementDimMatrix
const diffusion_WGpCap =
1323 phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dpCap;
1324 DisplacementDimMatrix
const diffusion_WLpCap =
1325 phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dpCap;
1327 DisplacementDimMatrix
const diffusion_WGT =
1328 phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dT;
1329 DisplacementDimMatrix
const diffusion_WLT =
1330 phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dT;
1332 DisplacementDimMatrix
const advection_W = advection_W_G + advection_W_L;
1333 DisplacementDimMatrix
const diffusion_W_pGR =
1334 diffusion_WGpGR + diffusion_WLpGR;
1335 DisplacementDimMatrix
const diffusion_W_pCap =
1336 diffusion_WGpCap + diffusion_WLpCap;
1338 DisplacementDimMatrix
const diffusion_W_T =
1339 diffusion_WGT + diffusion_WLT;
1342 gradNpT * (advection_W + diffusion_W_pGR) * gradNp * w;
1345 gradNpT * (diffusion_W_pCap - advection_W_L) * gradNp * w;
1347 LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
1349 fW.noalias() += gradNpT *
1350 (advection_W_G * ip.rhoGR + advection_W_L * ip.rhoLR) *
1353 if (!_process_data.apply_mass_lumping)
1355 fW.noalias() -= NpT *
1356 (phi * (ip.rhoWLR - ip.rhoWGR) -
1357 rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1362 NpT * phi * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
1368 MTu.noalias() += NTT * rho_h_eff * mT * Bu * w;
1370 KTT.noalias() += gradNTT * ip.lambda * gradNT * w;
1372 fT.noalias() -= NTT * rho_u_eff_dot * w;
1375 gradNTT * (ip.rhoGR * h_G * ip.w_GS + ip.rhoLR * h_L * ip.w_LS) * w;
1379 (ip.rhoCGR * ip.h_CG * ip.d_CG + ip.rhoWGR * ip.h_WG * ip.d_WG) * w;
1383 (ip.rhoGR * ip.w_GS.transpose() + ip.rhoLR * ip.w_LS.transpose()) *
1390 KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
1392 KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_s_L * m * Np) * w;
1395 (BuT * ip.sigma_eff - N_u_op(Nu).transpose() * rho * b) * w;
1397 if (_process_data.apply_mass_lumping)
1399 MCpG = MCpG.colwise().sum().eval().asDiagonal();
1400 MCpC = MCpC.colwise().sum().eval().asDiagonal();
1401 MWpG = MWpG.colwise().sum().eval().asDiagonal();
1402 MWpC = MWpC.colwise().sum().eval().asDiagonal();
1409template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
1410 int DisplacementDim>
1413 assembleWithJacobian(
double const t,
double const dt,
1414 std::vector<double>
const& local_x,
1415 std::vector<double>
const& local_x_prev,
1416 std::vector<double>& ,
1417 std::vector<double>& ,
1418 std::vector<double>& local_rhs_data,
1419 std::vector<double>& local_Jac_data)
1421 auto const matrix_size = gas_pressure_size + capillary_pressure_size +
1422 temperature_size + displacement_size;
1423 assert(local_x.size() == matrix_size);
1425 auto const temperature = Eigen::Map<VectorType<temperature_size>
const>(
1426 local_x.data() + temperature_index, temperature_size);
1428 auto const gas_pressure = Eigen::Map<VectorType<gas_pressure_size>
const>(
1429 local_x.data() + gas_pressure_index, gas_pressure_size);
1431 auto const capillary_pressure =
1432 Eigen::Map<VectorType<capillary_pressure_size>
const>(
1433 local_x.data() + capillary_pressure_index, capillary_pressure_size);
1435 auto const displacement = Eigen::Map<VectorType<displacement_size>
const>(
1436 local_x.data() + displacement_index, displacement_size);
1438 auto const gas_pressure_prev =
1439 Eigen::Map<VectorType<gas_pressure_size>
const>(
1440 local_x_prev.data() + gas_pressure_index, gas_pressure_size);
1442 auto const capillary_pressure_prev =
1443 Eigen::Map<VectorType<capillary_pressure_size>
const>(
1444 local_x_prev.data() + capillary_pressure_index,
1445 capillary_pressure_size);
1447 auto const temperature_prev =
1448 Eigen::Map<VectorType<temperature_size>
const>(
1449 local_x_prev.data() + temperature_index, temperature_size);
1451 auto const displacement_prev =
1452 Eigen::Map<VectorType<displacement_size>
const>(
1453 local_x_prev.data() + displacement_index, displacement_size);
1456 MathLib::createZeroedMatrix<MatrixType<matrix_size, matrix_size>>(
1457 local_Jac_data, matrix_size, matrix_size);
1459 auto local_f = MathLib::createZeroedVector<VectorType<matrix_size>>(
1460 local_rhs_data, matrix_size);
1471 C_size, capillary_pressure_size);
1481 C_size, capillary_pressure_size);
1490 W_size, capillary_pressure_size);
1501 W_size, capillary_pressure_size);
1508 temperature_size, displacement_size);
1518 displacement_size, gas_pressure_size);
1521 displacement_size, capillary_pressure_size);
1524 auto fC = local_f.template segment<C_size>(C_index);
1526 auto fW = local_f.template segment<W_size>(W_index);
1528 auto fT = local_f.template segment<temperature_size>(temperature_index);
1530 auto fU = local_f.template segment<displacement_size>(displacement_index);
1532 unsigned const n_integration_points =
1533 _integration_method.getNumberOfPoints();
1535 auto const ip_constitutive_variables = updateConstitutiveVariables(
1536 Eigen::Map<Eigen::VectorXd const>(local_x.data(), local_x.size()),
1537 Eigen::Map<Eigen::VectorXd const>(local_x_prev.data(),
1538 local_x_prev.size()),
1541 for (
unsigned int_point = 0; int_point < n_integration_points; int_point++)
1543 auto& ip = _ip_data[int_point];
1544 auto& ip_cv = ip_constitutive_variables[int_point];
1546 auto const& Np = ip.N_p;
1547 auto const& NT = Np;
1548 auto const& Nu = ip.N_u;
1550 std::nullopt, _element.getID(), int_point,
1556 auto const& NpT = Np.transpose().eval();
1557 auto const& NTT = NT.transpose().eval();
1559 auto const& gradNp = ip.dNdx_p;
1560 auto const& gradNT = gradNp;
1561 auto const& gradNu = ip.dNdx_u;
1563 auto const& gradNpT = gradNp.transpose().eval();
1564 auto const& gradNTT = gradNT.transpose().eval();
1566 auto const& w = ip.integration_weight;
1568 auto const& m = Invariants::identity2;
1569 auto const mT = m.transpose().eval();
1571 auto const x_coord =
1578 ShapeFunctionDisplacement::NPOINTS,
1580 gradNu, Nu, x_coord, _is_axially_symmetric);
1582 auto const BuT = Bu.transpose().eval();
1584 double const div_u_dot =
1585 Invariants::trace(Bu * (displacement - displacement_prev) / dt);
1587 double const pGR = Np.dot(gas_pressure);
1588 double const pCap = Np.dot(capillary_pressure);
1589 double const T = NT.dot(temperature);
1595 double const pGR_prev = Np.dot(gas_pressure_prev);
1596 double const pCap_prev = Np.dot(capillary_pressure_prev);
1597 double const T_prev = NT.dot(temperature_prev);
1598 auto& beta_T_SR = ip.beta_T_SR;
1601 Eigen::Matrix<double, DisplacementDim, DisplacementDim>::Identity();
1603 const double sD_G = ip.diffusion_coefficient_vapour;
1604 const double sD_L = ip.diffusion_coefficient_solute;
1614 auto const s_G = 1. - s_L;
1615 auto const s_L_dot = (s_L - ip.s_L_prev) / dt;
1617 auto& alpha_B = ip.alpha_B;
1618 auto& beta_p_SR = ip.beta_p_SR;
1620 auto const& b = _process_data.specific_body_force;
1626 auto const phi_G = s_G * phi;
1627 auto const phi_L = s_L * phi;
1628 auto const phi_S = 1. - phi;
1631 auto& rho_SR = ip.rhoSR;
1633 auto const rho = phi_G * ip.rhoGR + phi_L * ip.rhoLR + phi_S * rho_SR;
1636 const double rho_C_FR = s_G * ip.rhoCGR + s_L * ip.rhoCLR;
1637 const double rho_W_FR = s_G * ip.rhoWGR + s_L * ip.rhoWLR;
1643 auto const rho_C_GR_dot = (ip.rhoCGR - ip.rhoCGR_prev) / dt;
1644 auto const rho_C_LR_dot = (ip.rhoCLR - ip.rhoCLR_prev) / dt;
1645 auto const rho_W_GR_dot = (ip.rhoWGR - ip.rhoWGR_prev) / dt;
1646 auto const rho_W_LR_dot = (ip.rhoWLR - ip.rhoWLR_prev) / dt;
1648 auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
1650 auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / dt;
1659 MCpG.noalias() += NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1661 NpT * rho_C_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1663 if (_process_data.apply_mass_lumping)
1665 if (pCap - pCap_prev != 0.)
1669 (phi * (ip.rhoCLR - ip.rhoCGR) -
1670 rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1671 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1675 MCT.noalias() -= NpT * rho_C_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1678 .template block<C_size, temperature_size>(C_index,
1680 .noalias() += NpT * ip_cv.dfC_4_MCT_dT * (T - T_prev) / dt * NT * w;
1682 MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w;
1685 .template block<C_size, temperature_size>(C_index,
1687 .noalias() += NpT * ip_cv.dfC_4_MCu_dT * div_u_dot * NT * w;
1692 -phi_G * ip.rhoGR * D_C_G * ip.dxmWG_dpGR;
1694 -phi_L * ip.rhoLR * D_C_L * ip.dxmWL_dpLR;
1696 -phi_G * ip.rhoGR * D_C_G * ip.dxmWG_dT;
1698 -phi_L * ip.rhoLR * D_C_L * ip.dxmWL_dT;
1702 diffusion_C_G_p + diffusion_C_L_p;
1704 diffusion_C_G_T + diffusion_C_L_T;
1706 LCpG.noalias() += gradNpT * (advection_C + diffusion_C_p) * gradNp * w;
1709 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1711 (ip_cv.dadvection_C_dp_GR
1717 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
1719 (ip_cv.dadvection_C_dp_cap
1726 .template block<C_size, temperature_size>(C_index,
1728 .noalias() += gradNpT * ip_cv.dfC_4_LCpG_dT * gradpGR * NT * w;
1731 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
1732 NpT * ip_cv.dfC_4_MCpG_dp_GR * (pGR - pGR_prev) / dt * Np * w;
1736 .template block<C_size, temperature_size>(C_index,
1739 NpT * ip_cv.dfC_4_MCpG_dT * (pGR - pGR_prev) / dt * NT * w;
1742 gradNpT * (advection_C_L + diffusion_C_L_p) * gradNp * w;
1770 LCT.noalias() += gradNpT * diffusion_C_T * gradNp * w;
1773 fC.noalias() += gradNpT *
1774 (advection_C_G * ip.rhoGR + advection_C_L * ip.rhoLR) *
1777 if (!_process_data.apply_mass_lumping)
1780 auto const a = phi * (ip.rhoCLR - ip.rhoCGR) -
1781 rho_C_FR * pCap * (alpha_B - phi) * beta_p_SR;
1782 fC.noalias() -= NpT * a * s_L_dot * w;
1784 local_Jac.template block<C_size, C_size>(C_index, C_index)
1787 (ip_cv.dfC_2a_dp_GR * s_L_dot ) *
1790 local_Jac.template block<C_size, W_size>(C_index, W_index)
1793 (ip_cv.dfC_2a_dp_cap * s_L_dot + a * ip_cv.ds_L_dp_cap / dt) *
1797 .template block<C_size, temperature_size>(C_index,
1799 .noalias() += NpT * ip_cv.dfC_2a_dT * s_L_dot * NT * w;
1803 double const a = s_G * rho_C_GR_dot + s_L * rho_C_LR_dot;
1804 fC.noalias() -= NpT * phi * a * w;
1806 local_Jac.template block<C_size, C_size>(C_index, C_index)
1807 .noalias() += NpT * phi * ip_cv.dfC_3a_dp_GR * Np * w;
1809 local_Jac.template block<C_size, W_size>(C_index, W_index)
1810 .noalias() += NpT * phi * ip_cv.dfC_3a_dp_cap * Np * w;
1813 .template block<C_size, temperature_size>(C_index,
1817#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
1820 phi * ip_cv.dfC_3a_dT) *
1827 MWpG.noalias() += NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * Np * w;
1829 NpT * rho_W_FR * (alpha_B - phi) * beta_p_SR * s_L * Np * w;
1831 if (_process_data.apply_mass_lumping)
1833 if (pCap - pCap_prev != 0.)
1837 (phi * (ip.rhoWLR - ip.rhoWGR) -
1838 rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR) *
1839 s_L_dot * dt / (pCap - pCap_prev) * Np * w;
1843 MWT.noalias() -= NpT * rho_W_FR * (alpha_B - phi) * beta_T_SR * Np * w;
1845 MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w;
1850 phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dpGR;
1852 phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dpLR;
1854 phi_G * ip.rhoGR * D_W_G * ip.dxmWG_dT;
1856 phi_L * ip.rhoLR * D_W_L * ip.dxmWL_dT;
1860 diffusion_W_G_p + diffusion_W_L_p;
1862 diffusion_W_G_T + diffusion_W_L_T;
1864 LWpG.noalias() += gradNpT * (advection_W + diffusion_W_p) * gradNp * w;
1867 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1868 gradNpT * (ip_cv.dfW_4_LWpG_a_dp_GR + ip_cv.dfW_4_LWpG_d_dp_GR) *
1871 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1872 gradNpT * (ip_cv.dfW_4_LWpG_a_dp_cap + ip_cv.dfW_4_LWpG_d_dp_cap) *
1876 .template block<W_size, temperature_size>(W_index,
1878 .noalias() += gradNpT *
1879 (ip_cv.dfW_4_LWpG_a_dT + ip_cv.dfW_4_LWpG_d_dT) *
1883 gradNpT * (advection_W_L + diffusion_W_L_p) * gradNp * w;
1886 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() -=
1887 gradNpT * (ip_cv.dfW_4_LWpC_a_dp_GR + ip_cv.dfW_4_LWpC_d_dp_GR) *
1890 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() -=
1891 gradNpT * (ip_cv.dfW_4_LWpC_a_dp_cap + ip_cv.dfW_4_LWpC_d_dp_cap) *
1895 .template block<W_size, temperature_size>(W_index,
1897 .noalias() -= gradNpT *
1898 (ip_cv.dfW_4_LWpC_a_dT + ip_cv.dfW_4_LWpC_d_dT) *
1901 LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w;
1904 fW.noalias() += gradNpT *
1905 (advection_W_G * ip.rhoGR + advection_W_L * ip.rhoLR) *
1909 if (!_process_data.apply_mass_lumping)
1911 double const f = phi * (ip.rhoWLR - ip.rhoWGR);
1912 double const g = rho_W_FR * pCap * (alpha_B - phi) * beta_p_SR;
1914 fW.noalias() -= NpT * (f - g) * s_L_dot * w;
1916 local_Jac.template block<W_size, C_size>(W_index, C_index)
1917 .noalias() += NpT * (ip_cv.dfW_2a_dp_GR - ip_cv.dfW_2b_dp_GR) *
1923 local_Jac.template block<W_size, W_size>(W_index, W_index)
1926 ((ip_cv.dfW_2a_dp_cap - ip_cv.dfW_2b_dp_cap) * s_L_dot +
1927 (f - g) * ip_cv.ds_L_dp_cap / dt) *
1931 .template block<W_size, temperature_size>(W_index,
1934 NpT * (ip_cv.dfW_2a_dT - ip_cv.dfW_2b_dT) * s_L_dot * Np * w;
1939 NpT * phi * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w;
1941 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
1942 NpT * phi * ip_cv.dfW_3a_dp_GR * Np * w;
1944 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
1945 NpT * phi * ip_cv.dfW_3a_dp_cap * Np * w;
1948 .template block<W_size, temperature_size>(W_index,
1953#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
1954 ip.dphi_dT * (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) +
1956 phi * ip_cv.dfW_3a_dT) *
1963 MTu.noalias() += NTT * rho_h_eff * mT * Bu * w;
1968 .template block<temperature_size, C_size>(temperature_index,
1970 .noalias() += NTT * ip_cv.drho_h_eff_dp_GR * div_u_dot * NT * w;
1975 .template block<temperature_size, W_size>(temperature_index,
1977 .noalias() -= NTT * ip_cv.drho_h_eff_dp_cap * div_u_dot * NT * w;
1982 .template block<temperature_size, temperature_size>(
1983 temperature_index, temperature_index)
1984 .noalias() += NTT * ip_cv.drho_h_eff_dT * div_u_dot * NT * w;
1986 KTT.noalias() += gradNTT * ip.lambda * gradNT * w;
2004 .template block<temperature_size, W_size>(temperature_index,
2006 .noalias() += gradNTT * ip_cv.dlambda_dp_cap * gradT * Np * w;
2010 .template block<temperature_size, temperature_size>(
2011 temperature_index, temperature_index)
2012 .noalias() += gradNTT * ip_cv.dlambda_dT * gradT * NT * w;
2015 fT.noalias() -= NTT * rho_u_eff_dot * w;
2019 .template block<temperature_size, C_size>(temperature_index,
2021 .noalias() += NpT / dt * ip_cv.drho_u_eff_dp_GR * Np * w;
2025 .template block<temperature_size, W_size>(temperature_index,
2027 .noalias() += NpT / dt * ip_cv.drho_u_eff_dp_cap * Np * w;
2032 .template block<temperature_size, temperature_size>(
2033 temperature_index, temperature_index)
2034 .noalias() += NTT * ip_cv.drho_u_eff_dT / dt * NT * w;
2038 gradNTT * (ip.rhoGR * h_G * ip.w_GS + ip.rhoLR * h_L * ip.w_LS) * w;
2042 .template block<temperature_size, C_size>(temperature_index,
2046 gradNTT * ip_cv.drho_GR_h_w_eff_dp_GR_Npart * Np * w +
2048 gradNTT * ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart * gradNp * w;
2052 .template block<temperature_size, W_size>(temperature_index,
2056 gradNTT * (-ip_cv.drho_LR_h_w_eff_dp_cap_Npart) * Np * w +
2058 gradNTT * (-ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart) * gradNp * w;
2062 .template block<temperature_size, temperature_size>(
2063 temperature_index, temperature_index)
2064 .noalias() -= gradNTT * ip_cv.drho_GR_h_w_eff_dT * NT * w;
2069 (ip.rhoGR * ip.w_GS.transpose() + ip.rhoLR * ip.w_LS.transpose()) *
2074 (ip.rhoCGR * ip.h_CG * ip.d_CG + ip.rhoWGR * ip.h_WG * ip.d_WG) * w;
2080 KUpG.noalias() -= (BuT * alpha_B * m * Np) * w;
2085 KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_s_L * m * Np) * w;
2090 .template block<displacement_size, W_size>(displacement_index,
2092 .noalias() += BuT * alpha_B * ip_cv.dchi_ds_L * ip_cv.ds_L_dp_cap *
2096 .template block<displacement_size, displacement_size>(
2097 displacement_index, displacement_index)
2098 .noalias() += BuT * ip_cv.C * Bu * w;
2102 (BuT * ip.sigma_eff - N_u_op(Nu).transpose() * rho * b) * w;
2106 .template block<displacement_size, temperature_size>(
2107 displacement_index, temperature_index)
2108 .noalias() -= BuT * (ip_cv.C * ip.alpha_T_SR) * NT * w;
2118 if (_process_data.apply_mass_lumping)
2120 MCpG = MCpG.colwise().sum().eval().asDiagonal();
2121 MCpC = MCpC.colwise().sum().eval().asDiagonal();
2122 MWpG = MWpG.colwise().sum().eval().asDiagonal();
2123 MWpC = MWpC.colwise().sum().eval().asDiagonal();
2129 fC.noalias() -= LCpG * gas_pressure + LCpC * capillary_pressure +
2131 MCpG * (gas_pressure - gas_pressure_prev) / dt +
2132 MCpC * (capillary_pressure - capillary_pressure_prev) / dt +
2133 MCT * (temperature - temperature_prev) / dt +
2134 MCu * (displacement - displacement_prev) / dt;
2136 local_Jac.template block<C_size, C_size>(C_index, C_index).noalias() +=
2138 local_Jac.template block<C_size, W_size>(C_index, W_index).noalias() +=
2141 .template block<C_size, temperature_size>(C_index, temperature_index)
2142 .noalias() += LCT + MCT / dt;
2144 .template block<C_size, displacement_size>(C_index, displacement_index)
2145 .noalias() += MCu / dt;
2149 fW.noalias() -= LWpG * gas_pressure + LWpC * capillary_pressure +
2151 MWpG * (gas_pressure - gas_pressure_prev) / dt +
2152 MWpC * (capillary_pressure - capillary_pressure_prev) / dt +
2153 MWT * (temperature - temperature_prev) / dt +
2154 MWu * (displacement - displacement_prev) / dt;
2156 local_Jac.template block<W_size, W_size>(W_index, W_index).noalias() +=
2158 local_Jac.template block<W_size, C_size>(W_index, C_index).noalias() +=
2161 .template block<W_size, temperature_size>(W_index, temperature_index)
2162 .noalias() += LWT + MWT / dt;
2164 .template block<W_size, displacement_size>(W_index, displacement_index)
2165 .noalias() += MWu / dt;
2170 KTT * temperature + MTu * (displacement - displacement_prev) / dt;
2173 .template block<temperature_size, temperature_size>(temperature_index,
2177 .template block<temperature_size, displacement_size>(temperature_index,
2179 .noalias() += MTu / dt;
2183 fU.noalias() -= KUpG * gas_pressure + KUpC * capillary_pressure;
2186 .template block<displacement_size, C_size>(displacement_index, C_index)
2189 .template block<displacement_size, W_size>(displacement_index, W_index)
2193template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
2194 int DisplacementDim>
2196 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
2197 getIntPtDarcyVelocityGas(
2199 std::vector<GlobalVector*>
const& ,
2200 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
2201 std::vector<double>& cache)
const
2203 unsigned const n_integration_points =
2204 _integration_method.getNumberOfPoints();
2208 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2209 cache, DisplacementDim, n_integration_points);
2211 for (
unsigned ip = 0; ip < n_integration_points; ip++)
2213 cache_matrix.col(ip).noalias() = _ip_data[ip].w_GS;
2219template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
2220 int DisplacementDim>
2222 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
2223 getIntPtDarcyVelocityLiquid(
2225 std::vector<GlobalVector*>
const& ,
2226 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
2227 std::vector<double>& cache)
const
2229 unsigned const n_integration_points =
2230 _integration_method.getNumberOfPoints();
2234 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2235 cache, DisplacementDim, n_integration_points);
2237 for (
unsigned ip = 0; ip < n_integration_points; ip++)
2239 cache_matrix.col(ip).noalias() = _ip_data[ip].w_LS;
2245template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
2246 int DisplacementDim>
2248 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
2249 getIntPtDiffusionVelocityVapourGas(
2251 std::vector<GlobalVector*>
const& ,
2252 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
2253 std::vector<double>& cache)
const
2255 unsigned const n_integration_points =
2256 _integration_method.getNumberOfPoints();
2260 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2261 cache, DisplacementDim, n_integration_points);
2263 for (
unsigned ip = 0; ip < n_integration_points; ip++)
2265 cache_matrix.col(ip).noalias() = _ip_data[ip].d_WG;
2271template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
2272 int DisplacementDim>
2274 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
2275 getIntPtDiffusionVelocityGasGas(
2277 std::vector<GlobalVector*>
const& ,
2278 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
2279 std::vector<double>& cache)
const
2281 unsigned const n_integration_points =
2282 _integration_method.getNumberOfPoints();
2286 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2287 cache, DisplacementDim, n_integration_points);
2289 for (
unsigned ip = 0; ip < n_integration_points; ip++)
2291 cache_matrix.col(ip).noalias() = _ip_data[ip].d_CG;
2297template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
2298 int DisplacementDim>
2300 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
2301 getIntPtDiffusionVelocitySoluteLiquid(
2303 std::vector<GlobalVector*>
const& ,
2304 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
2305 std::vector<double>& cache)
const
2307 unsigned const n_integration_points =
2308 _integration_method.getNumberOfPoints();
2312 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2313 cache, DisplacementDim, n_integration_points);
2315 for (
unsigned ip = 0; ip < n_integration_points; ip++)
2317 cache_matrix.col(ip).noalias() = _ip_data[ip].d_CL;
2323template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
2324 int DisplacementDim>
2326 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
2327 getIntPtDiffusionVelocityLiquidLiquid(
2329 std::vector<GlobalVector*>
const& ,
2330 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
2331 std::vector<double>& cache)
const
2333 unsigned const n_integration_points =
2334 _integration_method.getNumberOfPoints();
2338 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
2339 cache, DisplacementDim, n_integration_points);
2341 for (
unsigned ip = 0; ip < n_integration_points; ip++)
2343 cache_matrix.col(ip).noalias() = _ip_data[ip].d_WL;
2349template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
2350 int DisplacementDim>
2353 computeSecondaryVariableConcrete(
double const t,
double const dt,
2354 Eigen::VectorXd
const& local_x,
2355 Eigen::VectorXd
const& local_x_prev)
2357 auto const gas_pressure =
2358 local_x.template segment<gas_pressure_size>(gas_pressure_index);
2359 auto const capillary_pressure =
2360 local_x.template segment<capillary_pressure_size>(
2361 capillary_pressure_index);
2362 auto const liquid_pressure = (gas_pressure - capillary_pressure).eval();
2365 ShapeFunctionPressure,
typename ShapeFunctionDisplacement::MeshElement,
2366 DisplacementDim>(_element, _is_axially_symmetric, gas_pressure,
2367 *_process_data.gas_pressure_interpolated);
2370 ShapeFunctionPressure,
typename ShapeFunctionDisplacement::MeshElement,
2371 DisplacementDim>(_element, _is_axially_symmetric, capillary_pressure,
2372 *_process_data.capillary_pressure_interpolated);
2375 ShapeFunctionPressure,
typename ShapeFunctionDisplacement::MeshElement,
2376 DisplacementDim>(_element, _is_axially_symmetric, liquid_pressure,
2377 *_process_data.liquid_pressure_interpolated);
2379 auto const temperature =
2380 local_x.template segment<temperature_size>(temperature_index);
2383 ShapeFunctionPressure,
typename ShapeFunctionDisplacement::MeshElement,
2384 DisplacementDim>(_element, _is_axially_symmetric, temperature,
2385 *_process_data.temperature_interpolated);
2387 unsigned const n_integration_points =
2388 _integration_method.getNumberOfPoints();
2390 double saturation_avg = 0;
2392 updateConstitutiveVariables(local_x, local_x_prev, t, dt);
2394 for (
unsigned ip = 0; ip < n_integration_points; ip++)
2396 saturation_avg += _ip_data[ip].s_L;
2398 saturation_avg /= n_integration_points;
2399 (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > mechanical_strain
double equivalent_plastic_strain
double capillary_pressure
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > total_stress
double liquid_phase_pressure
virtual std::size_t getID() const final
Returns the ID of the element.
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
unsigned getNumberOfPoints() const
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
typename ShapeMatricesTypePressure::template MatrixType< M, N > MatrixType
Eigen::Matrix< double, KelvinVectorSize, 1 > SymmetricTensor
ShapeMatrixPolicyType< ShapeFunctionPressure, DisplacementDim > ShapeMatricesTypePressure
SecondaryData< typename ShapeMatricesTypeDisplacement::ShapeMatrices::ShapeType > _secondary_data
TH2MProcessData< DisplacementDim > & _process_data
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
typename ShapeMatricesTypePressure::GlobalDimMatrixType GlobalDimMatrixType
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
typename ShapeMatricesTypePressure::GlobalDimVectorType GlobalDimVectorType
NumLib::GenericIntegrationMethod const & _integration_method
auto & selectSolidConstitutiveRelation(SolidMaterialsMap const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
template Eigen::Matrix< double, 3, 3 > formEigenTensor< 3 >(MaterialPropertyLib::PropertyDataType const &values)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
KelvinVectorType< DisplacementDim > tensorToKelvin(Eigen::Matrix< double, 3, 3 > const &m)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
std::size_t setIntegrationPointDataMaterialStateVariables(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType member, std::function< std::span< double >(MaterialStateVariables &)> get_values_span)
void setIPDataInitialConditions(std::vector< std::unique_ptr< MeshLib::IntegrationPointWriter > > const &_integration_point_writer, MeshLib::Properties const &mesh_properties, LocalAssemblersVector &local_assemblers, bool const remove_name_suffix=false)
std::size_t setIntegrationPointScalarData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N_u