39 ShapeFunctionPressure, DisplacementDim>::
40 HydroMechanicsLocalAssemblerMatrix(
42 std::size_t
const n_variables,
44 std::vector<unsigned>
const& dofIndex_to_localIndex,
46 bool const is_axially_symmetric,
49 e, is_axially_symmetric, integration_method,
50 (n_variables - 1) * ShapeFunctionDisplacement::NPOINTS *
52 ShapeFunctionPressure::NPOINTS,
53 dofIndex_to_localIndex),
54 _process_data(process_data)
56 unsigned const n_integration_points =
59 _ip_data.reserve(n_integration_points);
62 auto const shape_matrices_u =
65 DisplacementDim>(e, is_axially_symmetric,
68 auto const shape_matrices_p =
71 e, is_axially_symmetric, integration_method);
78 for (
unsigned ip = 0; ip < n_integration_points; ip++)
80 _ip_data.emplace_back(solid_material);
82 auto const& sm_u = shape_matrices_u[ip];
83 auto const& sm_p = shape_matrices_p[ip];
85 ip_data.integration_weight =
86 sm_u.detJ * sm_u.integralMeasure *
88 ip_data.darcy_velocity.setZero();
91 ip_data.dNdx_u = sm_u.dNdx;
93 for (
int i = 0; i < DisplacementDim; ++i)
98 .noalias() = ip_data.N_u;
101 ip_data.N_p = sm_p.N;
102 ip_data.dNdx_p = sm_p.dNdx;
113 auto const initial_effective_stress =
117 ip_data.sigma_eff[i] = initial_effective_stress[i];
118 ip_data.sigma_eff_prev[i] = initial_effective_stress[i];
126 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
127 assembleWithJacobianConcrete(
double const t,
double const dt,
128 Eigen::VectorXd
const& local_x,
129 Eigen::VectorXd
const& local_x_prev,
130 Eigen::VectorXd& local_rhs,
131 Eigen::MatrixXd& local_Jac)
133 auto p =
const_cast<Eigen::VectorXd&
>(local_x).segment(pressure_index,
135 auto const p_prev = local_x_prev.segment(pressure_index, pressure_size);
137 if (_process_data.deactivate_matrix_in_flow)
139 setPressureOfInactiveNodes(t, p);
142 auto u = local_x.segment(displacement_index, displacement_size);
143 auto u_prev = local_x_prev.segment(displacement_index, displacement_size);
145 auto rhs_p = local_rhs.template segment<pressure_size>(pressure_index);
147 local_rhs.template segment<displacement_size>(displacement_index);
149 auto J_pp = local_Jac.template block<pressure_size, pressure_size>(
150 pressure_index, pressure_index);
151 auto J_pu = local_Jac.template block<pressure_size, displacement_size>(
152 pressure_index, displacement_index);
153 auto J_uu = local_Jac.template block<displacement_size, displacement_size>(
154 displacement_index, displacement_index);
155 auto J_up = local_Jac.template block<displacement_size, pressure_size>(
156 displacement_index, pressure_index);
158 assembleBlockMatricesWithJacobian(t, dt, p, p_prev, u, u_prev, rhs_p, rhs_u,
159 J_pp, J_pu, J_uu, J_up);
165 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
166 assembleBlockMatricesWithJacobian(
167 double const t,
double const dt,
168 Eigen::Ref<const Eigen::VectorXd>
const& p,
169 Eigen::Ref<const Eigen::VectorXd>
const& p_prev,
170 Eigen::Ref<const Eigen::VectorXd>
const& u,
171 Eigen::Ref<const Eigen::VectorXd>
const& u_prev,
172 Eigen::Ref<Eigen::VectorXd> rhs_p, Eigen::Ref<Eigen::VectorXd> rhs_u,
173 Eigen::Ref<Eigen::MatrixXd> J_pp, Eigen::Ref<Eigen::MatrixXd> J_pu,
174 Eigen::Ref<Eigen::MatrixXd> J_uu, Eigen::Ref<Eigen::MatrixXd> J_up)
176 assert(this->_element.getDimension() == DisplacementDim);
179 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
183 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
186 typename ShapeMatricesTypeDisplacement::template MatrixType<
187 displacement_size, pressure_size>
188 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
189 displacement_size, pressure_size>::Zero(displacement_size,
192 auto const& gravity_vec = _process_data.specific_body_force;
199 auto const& medium = _process_data.media_map.getMedium(_element.getID());
200 auto const& liquid_phase = medium->phase(
"AqueousLiquid");
201 auto const& solid_phase = medium->phase(
"Solid");
204 medium->property(MPL::PropertyType::reference_temperature)
205 .template value<double>(variables, x_position, t, dt);
209 bool const has_storage_property =
210 medium->hasProperty(MPL::PropertyType::storage);
212 auto const B_dil_bar = getDilatationalBBarMatrix();
214 unsigned const n_integration_points = _ip_data.size();
215 for (
unsigned ip = 0; ip < n_integration_points; ip++)
217 auto& ip_data = _ip_data[ip];
218 auto const& ip_w = ip_data.integration_weight;
219 auto const& N_u = ip_data.N_u;
220 auto const& dNdx_u = ip_data.dNdx_u;
221 auto const& N_p = ip_data.N_p;
222 auto const& dNdx_p = ip_data.dNdx_p;
223 auto const& H_u = ip_data.H_u;
228 std::nullopt, _element.getID(),
238 DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
240 dNdx_u, N_u, B_dil_bar, x_coord, this->_is_axially_symmetric)
243 auto const& eps_prev = ip_data.eps_prev;
244 auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
245 auto& sigma_eff = ip_data.sigma_eff;
247 auto& eps = ip_data.eps;
248 auto& state = ip_data.material_state_variables;
251 solid_phase.property(MPL::PropertyType::density)
252 .template value<double>(variables, x_position, t, dt);
254 liquid_phase.property(MPL::PropertyType::density)
255 .template value<double>(variables, x_position, t, dt);
258 medium->property(MPL::PropertyType::biot_coefficient)
259 .template value<double>(variables, x_position, t, dt);
260 auto const porosity =
261 medium->property(MPL::PropertyType::porosity)
262 .template value<double>(variables, x_position, t, dt);
264 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
265 auto const& identity2 =
268 eps.noalias() = B * u;
282 auto&& solution = _ip_data[ip].solid_material.integrateStress(
283 variables_prev, variables, t, x_position, dt, *state);
287 OGS_FATAL(
"Computation of local constitutive relation failed.");
291 std::tie(sigma_eff, state, C) = std::move(*solution);
293 J_uu.noalias() += B.transpose() * C * B * ip_w;
295 rhs_u.noalias() -= B.transpose() * sigma_eff * ip_w;
296 rhs_u.noalias() -= -H_u.transpose() * rho * gravity_vec * ip_w;
302 if (!_process_data.deactivate_matrix_in_flow)
305 Kup.noalias() += B.transpose() * alpha * identity2 * N_p * ip_w;
309 liquid_phase.property(MPL::PropertyType::viscosity)
310 .template value<double>(variables, x_position, t, dt);
312 auto const k_over_mu =
314 medium->property(MPL::PropertyType::permeability)
315 .value(variables, x_position, t, dt)) /
321 ? medium->property(MPL::PropertyType::storage)
322 .template value<double>(variables, x_position, t, dt)
324 (liquid_phase.property(MPL::PropertyType::density)
325 .
template dValue<double>(
327 MPL::Variable::liquid_phase_pressure,
330 (alpha - porosity) * (1.0 - alpha) /
331 _ip_data[ip].solid_material.getBulkModulus(
334 auto q = ip_data.darcy_velocity.head(DisplacementDim);
335 q.noalias() = -k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
337 laplace_p.noalias() +=
338 dNdx_p.transpose() * k_over_mu * dNdx_p * ip_w;
339 storage_p.noalias() += N_p.transpose() * S * N_p * ip_w;
342 dNdx_p.transpose() * rho_fr * k_over_mu * gravity_vec * ip_w;
347 J_up.noalias() -= Kup;
350 J_pp.noalias() += laplace_p + storage_p / dt;
353 J_pu.noalias() += Kup.transpose() / dt;
356 rhs_p.noalias() -= laplace_p * p + storage_p * (p - p_prev) / dt +
357 Kup.transpose() * (u - u_prev) / dt;
360 rhs_u.noalias() -= -Kup * p;
384 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
385 postTimestepConcreteWithBlockVectors(
386 double const t,
double const dt,
387 Eigen::Ref<const Eigen::VectorXd>
const& p,
388 Eigen::Ref<const Eigen::VectorXd>
const& u)
394 auto const e_id = _element.getID();
398 KV sigma_avg = KV::Zero();
400 velocity_avg.setZero();
402 unsigned const n_integration_points = _ip_data.size();
404 auto const& medium = _process_data.media_map.getMedium(_element.getID());
405 auto const& liquid_phase = medium->phase(
"AqueousLiquid");
407 medium->property(MPL::PropertyType::reference_temperature)
408 .template value<double>(variables, x_position, t, dt);
412 auto const B_dil_bar = getDilatationalBBarMatrix();
414 for (
unsigned ip = 0; ip < n_integration_points; ip++)
416 auto& ip_data = _ip_data[ip];
418 auto const& eps_prev = ip_data.eps_prev;
419 auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
421 auto& eps = ip_data.eps;
422 auto& sigma_eff = ip_data.sigma_eff;
423 auto& state = ip_data.material_state_variables;
425 auto const& N_u = ip_data.N_u;
426 auto const& N_p = ip_data.N_p;
427 auto const& dNdx_u = ip_data.dNdx_u;
432 std::nullopt, _element.getID(),
442 DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
444 dNdx_u, N_u, B_dil_bar, x_coord, this->_is_axially_symmetric)
447 eps.noalias() = B * u;
460 auto&& solution = _ip_data[ip].solid_material.integrateStress(
461 variables_prev, variables, t, x_position, dt, *state);
465 OGS_FATAL(
"Computation of local constitutive relation failed.");
469 std::tie(sigma_eff, state, C) = std::move(*solution);
471 sigma_avg += ip_data.sigma_eff;
473 if (!_process_data.deactivate_matrix_in_flow)
477 liquid_phase.property(MPL::PropertyType::density)
478 .template value<double>(variables, x_position, t, dt);
481 liquid_phase.property(MPL::PropertyType::viscosity)
482 .template value<double>(variables, x_position, t, dt);
485 medium->property(MPL::PropertyType::permeability)
486 .value(variables, x_position, t, dt))
491 auto const& gravity_vec = _process_data.specific_body_force;
492 auto const& dNdx_p = ip_data.dNdx_p;
494 ip_data.darcy_velocity.head(DisplacementDim).noalias() =
495 -k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
496 velocity_avg.noalias() +=
497 ip_data.darcy_velocity.head(DisplacementDim);
501 sigma_avg /= n_integration_points;
502 velocity_avg /= n_integration_points;
505 &(*_process_data.element_stresses)[e_id * KV::RowsAtCompileTime]) =
508 Eigen::Map<DisplacementDimVector>(
509 &(*_process_data.element_velocities)[e_id * DisplacementDim]) =
513 ShapeFunctionPressure,
typename ShapeFunctionDisplacement::MeshElement,
514 DisplacementDim>(_element, _is_axially_symmetric, p,
515 *_process_data.mesh_prop_nodal_p);