103 ShapeFunctionPressure, DisplacementDim>::
104 assembleWithJacobian(
double const t,
double const dt,
105 std::vector<double>
const& local_x,
106 std::vector<double>
const& local_x_prev,
107 std::vector<double>& local_rhs_data,
108 std::vector<double>& local_Jac_data)
110 assert(local_x.size() == pressure_size + displacement_size);
112 auto p = Eigen::Map<
typename ShapeMatricesTypePressure::template VectorType<
113 pressure_size>
const>(local_x.data() + pressure_index, pressure_size);
116 Eigen::Map<
typename ShapeMatricesTypeDisplacement::template VectorType<
117 displacement_size>
const>(local_x.data() + displacement_index,
121 Eigen::Map<
typename ShapeMatricesTypePressure::template VectorType<
122 pressure_size>
const>(local_x_prev.data() + pressure_index,
125 Eigen::Map<
typename ShapeMatricesTypeDisplacement::template VectorType<
126 displacement_size>
const>(local_x_prev.data() + displacement_index,
130 typename ShapeMatricesTypeDisplacement::template MatrixType<
131 displacement_size + pressure_size,
132 displacement_size + pressure_size>>(
133 local_Jac_data, displacement_size + pressure_size,
134 displacement_size + pressure_size);
137 typename ShapeMatricesTypeDisplacement::template VectorType<
138 displacement_size + pressure_size>>(
139 local_rhs_data, displacement_size + pressure_size);
142 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
146 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
150 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
153 typename ShapeMatricesTypeDisplacement::template MatrixType<
154 displacement_size, pressure_size>
155 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
156 displacement_size, pressure_size>::Zero(displacement_size,
159 typename ShapeMatricesTypeDisplacement::template MatrixType<
160 pressure_size, displacement_size>
161 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
162 pressure_size, displacement_size>::Zero(pressure_size,
165 typename ShapeMatricesTypeDisplacement::template MatrixType<
166 pressure_size, displacement_size>
167 Kpu_k = ShapeMatricesTypeDisplacement::template MatrixType<
168 pressure_size, displacement_size>::Zero(pressure_size,
171 auto const& solid_material =
173 _process_data.solid_materials, _process_data.material_ids,
179 unsigned const n_integration_points =
180 _integration_method.getNumberOfPoints();
182 auto const& b = _process_data.specific_body_force;
183 auto const& medium = _process_data.media_map.getMedium(_element.getID());
184 auto const& solid = medium->phase(
"Solid");
185 auto const& fluid = fluidPhase(*medium);
191 medium->property(MPL::PropertyType::reference_temperature)
192 .template value<double>(vars, x_position, t, dt);
195 auto const& identity2 = Invariants::identity2;
197 for (
unsigned ip = 0; ip < n_integration_points; ip++)
199 auto const& w = _ip_data[ip].integration_weight;
201 auto const& N_u = _ip_data[ip].N_u;
202 auto const& dNdx_u = _ip_data[ip].dNdx_u;
204 auto const& N_p = _ip_data[ip].N_p;
205 auto const& dNdx_p = _ip_data[ip].dNdx_p;
213 ShapeFunctionDisplacement::NPOINTS,
215 dNdx_u, N_u, x_coord, _is_axially_symmetric);
217 auto& eps = _ip_data[ip].eps;
218 eps.noalias() = B * u;
219 auto const& sigma_eff = _ip_data[ip].sigma_eff;
221 double const p_int_pt = N_p.dot(p);
223 phase_pressure = p_int_pt;
225 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
226 t, x_position, dt, T_ref);
227 auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
229 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
230 .template value<double>(vars, x_position, t, dt);
233 solid.property(MPL::PropertyType::density)
234 .template value<double>(vars, x_position, t, dt);
235 auto const porosity =
236 medium->property(MPL::PropertyType::porosity)
237 .template value<double>(vars, x_position, t, dt);
239 auto const [rho_fr, mu] =
244 fluid.property(MPL::PropertyType::density)
245 .template dValue<double>(vars, _process_data.phase_variable,
252 auto const sigma_total =
253 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
262 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
268 medium->property(MPL::PropertyType::permeability)
269 .value(vars, x_position, t, dt));
272 (*medium)[MPL::PropertyType::permeability].dValue(
273 vars, MPL::Variable::mechanical_strain, x_position, t, dt));
275 auto const K_over_mu = K / mu;
277 auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
284 .template block<displacement_size, displacement_size>(
285 displacement_index, displacement_index)
286 .noalias() += B.transpose() * C * B * w;
288 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
289 local_rhs.template segment<displacement_size>(displacement_index)
291 (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
296 Kup.noalias() += B.transpose() * alpha * identity2 * N_p * w;
301 laplace_p.noalias() +=
302 rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
304 storage_p.noalias() +=
305 rho_fr * N_p.transpose() * N_p * w *
306 ((alpha - porosity) * (1.0 - alpha) / K_S + porosity * beta_p);
310 add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
312 (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
314 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
315 dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
321 rho_fr * alpha * N_p.transpose() * identity2.transpose() * B * w;
326 dNdx_p * p - rho_fr * b) *
327 dkde * B * rho_fr / mu * w;
331 .template block<displacement_size, pressure_size>(displacement_index,
335 if (_process_data.mass_lumping)
337 storage_p = storage_p.colwise().sum().eval().asDiagonal();
339 if constexpr (pressure_size == displacement_size)
341 Kpu = Kpu.colwise().sum().eval().asDiagonal();
342 Kpu_k = Kpu_k.colwise().sum().eval().asDiagonal();
348 .template block<pressure_size, pressure_size>(pressure_index,
350 .noalias() += laplace_p + storage_p / dt + add_p_derivative;
354 .template block<pressure_size, displacement_size>(pressure_index,
356 .noalias() += Kpu / dt + Kpu_k;
359 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
360 laplace_p * p + storage_p * (p - p_prev) / dt + Kpu * (u - u_prev) / dt;
363 local_rhs.template segment<displacement_size>(displacement_index)
364 .noalias() += Kup * p;
370 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
371 getIntPtDarcyVelocity(
373 std::vector<GlobalVector*>
const& x,
374 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
375 std::vector<double>& cache)
const
377 int const hydraulic_process_id = _process_data.hydraulic_process_id;
380 assert(!indices.empty());
381 auto const local_x = x[hydraulic_process_id]->get(indices);
383 unsigned const n_integration_points =
384 _integration_method.getNumberOfPoints();
387 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
388 cache, DisplacementDim, n_integration_points);
390 auto p = Eigen::Map<
typename ShapeMatricesTypePressure::template VectorType<
391 pressure_size>
const>(local_x.data() + pressure_index, pressure_size);
396 auto const& medium = _process_data.media_map.getMedium(_element.getID());
397 auto const& fluid = fluidPhase(*medium);
404 double const dt = std::numeric_limits<double>::quiet_NaN();
406 medium->property(MPL::PropertyType::reference_temperature)
407 .template value<double>(vars, x_position, t, dt);
409 auto const& identity2 = Invariants::identity2;
411 for (
unsigned ip = 0; ip < n_integration_points; ip++)
413 double const p_int_pt = _ip_data[ip].N_p.dot(p);
415 phase_pressure = p_int_pt;
417 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
418 .template value<double>(vars, x_position, t, dt);
422 auto const sigma_total =
423 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
430 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
436 medium->property(MPL::PropertyType::permeability)
437 .value(vars, x_position, t, dt));
439 auto const [rho_fr, mu] =
443 auto const K_over_mu = K / mu;
445 auto const& b = _process_data.specific_body_force;
448 auto const& dNdx_p = _ip_data[ip].dNdx_p;
449 cache_matrix.col(ip).noalias() =
450 -K_over_mu * dNdx_p * p + K_over_mu * rho_fr * b;
459 ShapeFunctionPressure, DisplacementDim>::
460 assembleWithJacobianForPressureEquations(
461 const double t,
double const dt, Eigen::VectorXd
const& local_x,
462 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
463 std::vector<double>& local_Jac_data)
467 template VectorType<pressure_size>>(
468 local_b_data, pressure_size);
473 auto const p = local_x.template segment<pressure_size>(pressure_index);
476 local_x_prev.template segment<pressure_size>(pressure_index);
479 typename ShapeMatricesTypeDisplacement::template MatrixType<
480 pressure_size, pressure_size>>(local_Jac_data, pressure_size,
484 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
488 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
492 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
495 auto const& solid_material =
497 _process_data.solid_materials, _process_data.material_ids,
503 auto const& medium = _process_data.media_map.getMedium(_element.getID());
504 auto const& fluid = fluidPhase(*medium);
510 medium->property(MPL::PropertyType::reference_temperature)
511 .template value<double>(vars, x_position, t, dt);
514 auto const& identity2 = Invariants::identity2;
516 auto const staggered_scheme =
517 std::get<Staggered>(_process_data.coupling_scheme);
518 auto const fixed_stress_stabilization_parameter =
519 staggered_scheme.fixed_stress_stabilization_parameter;
520 auto const fixed_stress_over_time_step =
521 staggered_scheme.fixed_stress_over_time_step;
523 int const n_integration_points = _integration_method.getNumberOfPoints();
524 for (
int ip = 0; ip < n_integration_points; ip++)
526 auto const& w = _ip_data[ip].integration_weight;
528 auto const& N_p = _ip_data[ip].N_p;
529 auto const& dNdx_p = _ip_data[ip].dNdx_p;
531 double const p_int_pt = N_p.dot(p);
533 phase_pressure = p_int_pt;
535 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
536 t, x_position, dt, T_ref);
537 auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
540 medium->property(MPL::PropertyType::biot_coefficient)
541 .template value<double>(vars, x_position, t, dt);
545 auto const sigma_total =
546 (_ip_data[ip].sigma_eff - alpha_b * p_int_pt * identity2).eval();
553 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
559 medium->property(MPL::PropertyType::permeability)
560 .value(vars, x_position, t, dt));
561 auto const porosity =
562 medium->property(MPL::PropertyType::porosity)
563 .template value<double>(vars, x_position, t, dt);
565 auto const [rho_fr, mu] =
570 fluid.property(MPL::PropertyType::density)
571 .template dValue<double>(vars, _process_data.phase_variable,
575 auto const K_over_mu = K / mu;
578 rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
582 fixed_stress_stabilization_parameter * alpha_b * alpha_b / K_S;
584 storage.noalias() += rho_fr * N_p.transpose() * N_p * w *
585 ((alpha_b - porosity) * (1.0 - alpha_b) / K_S +
586 porosity * beta_p + beta_FS);
588 auto const& b = _process_data.specific_body_force;
591 local_rhs.noalias() +=
592 dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
597 add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
599 (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
601 if (!fixed_stress_over_time_step)
603 auto const& eps = _ip_data[ip].eps;
604 auto const& eps_prev = _ip_data[ip].eps_prev;
605 const double eps_v_dot =
606 (Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
609 double const strain_rate_b =
610 alpha_b * eps_v_dot -
611 beta_FS * _ip_data[ip].strain_rate_variable;
613 local_rhs.noalias() -= strain_rate_b * rho_fr * N_p * w;
618 local_rhs.noalias() -=
619 alpha_b * _ip_data[ip].strain_rate_variable * rho_fr * N_p * w;
622 local_Jac.noalias() = laplace + storage / dt + add_p_derivative;
624 local_rhs.noalias() -= laplace * p + storage * (p - p_prev) / dt;
630 ShapeFunctionPressure, DisplacementDim>::
631 assembleWithJacobianForDeformationEquations(
632 const double t,
double const dt, Eigen::VectorXd
const& local_x,
633 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
635 auto const p = local_x.template segment<pressure_size>(pressure_index);
637 local_x.template segment<displacement_size>(displacement_index);
640 typename ShapeMatricesTypeDisplacement::template MatrixType<
641 displacement_size, displacement_size>>(
642 local_Jac_data, displacement_size, displacement_size);
646 template VectorType<displacement_size>>(
647 local_b_data, displacement_size);
652 auto const& medium = _process_data.media_map.getMedium(_element.getID());
653 auto const& solid = medium->phase(
"Solid");
654 auto const& fluid = fluidPhase(*medium);
660 medium->property(MPL::PropertyType::reference_temperature)
661 .template value<double>(vars, x_position, t, dt);
664 int const n_integration_points = _integration_method.getNumberOfPoints();
665 for (
int ip = 0; ip < n_integration_points; ip++)
667 auto const& w = _ip_data[ip].integration_weight;
669 auto const& N_u = _ip_data[ip].N_u;
670 auto const& dNdx_u = _ip_data[ip].dNdx_u;
672 auto const& N_p = _ip_data[ip].N_p;
680 ShapeFunctionDisplacement::NPOINTS,
682 dNdx_u, N_u, x_coord, _is_axially_symmetric);
684 auto& eps = _ip_data[ip].eps;
685 auto const& sigma_eff = _ip_data[ip].sigma_eff;
687 phase_pressure = N_p.dot(p);
689 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
690 .template value<double>(vars, x_position, t, dt);
692 solid.property(MPL::PropertyType::density)
693 .template value<double>(vars, x_position, t, dt);
694 auto const porosity =
695 medium->property(MPL::PropertyType::porosity)
696 .template value<double>(vars, x_position, t, dt);
699 t, dt, x_position, fluid, vars);
701 auto const& b = _process_data.specific_body_force;
704 DisplacementDim)>::identity2;
706 eps.noalias() = B * u;
711 auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
714 local_Jac.noalias() += B.transpose() * C * B * w;
719 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
720 local_rhs.noalias() -=
721 (B.transpose() * (sigma_eff - alpha * identity2 * p_at_xi) -
722 N_u_op(N_u).transpose() * rho * b) *
812 ShapeFunctionPressure, DisplacementDim>::
813 postNonLinearSolverConcrete(Eigen::VectorXd
const& local_x,
814 Eigen::VectorXd
const& local_x_prev,
815 double const t,
double const dt,
816 int const process_id)
822 int const n_integration_points = _integration_method.getNumberOfPoints();
824 auto const staggered_scheme_ptr =
825 std::get_if<Staggered>(&_process_data.coupling_scheme);
827 if (staggered_scheme_ptr &&
828 process_id == _process_data.hydraulic_process_id)
830 if (!staggered_scheme_ptr->fixed_stress_over_time_step)
833 local_x.template segment<pressure_size>(pressure_index);
835 local_x_prev.template segment<pressure_size>(pressure_index);
837 for (
int ip = 0; ip < n_integration_points; ip++)
839 auto& ip_data = _ip_data[ip];
841 auto const& N_p = ip_data.N_p;
843 ip_data.strain_rate_variable = N_p.dot(p - p_prev) / dt;
848 if (!staggered_scheme_ptr ||
849 process_id == _process_data.mechanics_related_process_id)
855 _process_data.media_map.getMedium(_element.getID());
858 medium->property(MPL::PropertyType::reference_temperature)
863 local_x.template segment<displacement_size>(displacement_index);
868 for (
int ip = 0; ip < n_integration_points; ip++)
870 auto const& N_u = _ip_data[ip].N_u;
871 auto const& dNdx_u = _ip_data[ip].dNdx_u;
878 DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
880 _is_axially_symmetric);
882 auto& eps = _ip_data[ip].eps;
883 eps.noalias() = B * u;
887 _ip_data[ip].updateConstitutiveRelation(vars, t, x_position, dt, u,
897 DisplacementDim>::postTimestepConcrete(Eigen::VectorXd
const& local_x,
898 Eigen::VectorXd
const& local_x_prev,
899 double const t,
double const dt,
900 int const process_id)
902 auto const staggered_scheme_ptr =
903 std::get_if<Staggered>(&_process_data.coupling_scheme);
905 if (staggered_scheme_ptr &&
906 process_id == _process_data.hydraulic_process_id)
908 if (staggered_scheme_ptr->fixed_stress_over_time_step)
910 auto const fixed_stress_stabilization_parameter =
911 staggered_scheme_ptr->fixed_stress_stabilization_parameter;
914 local_x.template segment<pressure_size>(pressure_index);
916 local_x_prev.template segment<pressure_size>(pressure_index);
921 auto const& solid_material =
923 _process_data.solid_materials, _process_data.material_ids,
927 _process_data.media_map.getMedium(_element.getID());
931 medium->property(MPL::PropertyType::reference_temperature)
932 .template value<double>(vars, x_position, t, dt);
935 int const n_integration_points =
936 _integration_method.getNumberOfPoints();
937 for (
int ip = 0; ip < n_integration_points; ip++)
939 auto& ip_data = _ip_data[ip];
941 auto const& N_p = ip_data.N_p;
943 auto const& eps = ip_data.eps;
944 auto const& eps_prev = ip_data.eps_prev;
945 const double eps_v_dot =
946 (Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
948 auto const C_el = ip_data.computeElasticTangentStiffness(
949 t, x_position, dt, T_ref);
951 solid_material.getBulkModulus(t, x_position, &C_el);
954 medium->property(MPL::PropertyType::biot_coefficient)
955 .template value<double>(vars, x_position, t, dt);
957 ip_data.strain_rate_variable =
958 eps_v_dot - fixed_stress_stabilization_parameter * alpha_b *
959 N_p.dot(p - p_prev) / dt / K_S;
964 unsigned const n_integration_points =
965 _integration_method.getNumberOfPoints();
967 for (
unsigned ip = 0; ip < n_integration_points; ip++)
969 _ip_data[ip].pushBackState();
1079 ShapeFunctionPressure, DisplacementDim>::
1080 computeSecondaryVariableConcrete(
double const t,
double const dt,
1081 Eigen::VectorXd
const& local_x,
1082 Eigen::VectorXd
const& )
1084 auto const p = local_x.template segment<pressure_size>(pressure_index);
1087 ShapeFunctionPressure,
typename ShapeFunctionDisplacement::MeshElement,
1088 DisplacementDim>(_element, _is_axially_symmetric, p,
1089 *_process_data.pressure_interpolated);
1091 int const elem_id = _element.getID();
1094 unsigned const n_integration_points =
1095 _integration_method.getNumberOfPoints();
1097 auto const& medium = _process_data.media_map.getMedium(elem_id);
1104 Eigen::Matrix<double, 3, 3>::Zero());
1106 auto const& identity2 = Invariants::identity2;
1108 for (
unsigned ip = 0; ip < n_integration_points; ip++)
1110 auto const& eps = _ip_data[ip].eps;
1111 sigma_eff_sum += _ip_data[ip].sigma_eff;
1113 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
1114 .template value<double>(vars, x_position, t, dt);
1115 double const p_int_pt = _ip_data[ip].N_p.dot(p);
1117 phase_pressure = p_int_pt;
1122 auto const sigma_total =
1123 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
1131 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
1137 medium->property(MPL::PropertyType::permeability)
1138 .value(vars, x_position, t, dt));
1141 Eigen::Map<Eigen::VectorXd>(
1142 &(*_process_data.permeability)[elem_id * KelvinVectorSize],
1143 KelvinVectorSize) = k_sum / n_integration_points;
1145 Eigen::Matrix<double, 3, 3, 0, 3, 3>
const sigma_avg =
1147 n_integration_points;
1149 Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(sigma_avg);
1151 Eigen::Map<Eigen::Vector3d>(
1152 &(*_process_data.principal_stress_values)[elem_id * 3], 3) =
1155 auto eigen_vectors = e_s.eigenvectors();
1157 for (
auto i = 0; i < 3; i++)
1159 Eigen::Map<Eigen::Vector3d>(
1160 &(*_process_data.principal_stress_vector[i])[elem_id * 3], 3) =
1161 eigen_vectors.col(i);