39 ShapeFunctionPressure, GlobalDim>::
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 * GlobalDim +
51 ShapeFunctionPressure::NPOINTS,
52 dofIndex_to_localIndex),
53 _process_data(process_data)
55 unsigned const n_integration_points =
58 _ip_data.reserve(n_integration_points);
61 auto const shape_matrices_u =
64 e, is_axially_symmetric, integration_method);
66 auto const shape_matrices_p =
69 e, is_axially_symmetric, integration_method);
76 for (
unsigned ip = 0; ip < n_integration_points; ip++)
78 _ip_data.emplace_back(solid_material);
80 auto const& sm_u = shape_matrices_u[ip];
81 auto const& sm_p = shape_matrices_p[ip];
83 ip_data.integration_weight =
84 sm_u.detJ * sm_u.integralMeasure *
86 ip_data.darcy_velocity.setZero();
89 ip_data.dNdx_u = sm_u.dNdx;
91 for (
int i = 0; i < GlobalDim; ++i)
96 .noalias() = ip_data.N_u;
100 ip_data.dNdx_p = sm_p.dNdx;
111 auto const initial_effective_stress =
115 ip_data.sigma_eff[i] = initial_effective_stress[i];
116 ip_data.sigma_eff_prev[i] = initial_effective_stress[i];
124 ShapeFunctionPressure, GlobalDim>::
125 assembleWithJacobianConcrete(
double const t,
double const dt,
126 Eigen::VectorXd
const& local_x,
127 Eigen::VectorXd
const& local_x_prev,
128 Eigen::VectorXd& local_rhs,
129 Eigen::MatrixXd& local_Jac)
131 auto p =
const_cast<Eigen::VectorXd&
>(local_x).segment(pressure_index,
133 auto const p_prev = local_x_prev.segment(pressure_index, pressure_size);
135 if (_process_data.deactivate_matrix_in_flow)
137 setPressureOfInactiveNodes(t, p);
140 auto u = local_x.segment(displacement_index, displacement_size);
141 auto u_prev = local_x_prev.segment(displacement_index, displacement_size);
143 auto rhs_p = local_rhs.template segment<pressure_size>(pressure_index);
145 local_rhs.template segment<displacement_size>(displacement_index);
147 auto J_pp = local_Jac.template block<pressure_size, pressure_size>(
148 pressure_index, pressure_index);
149 auto J_pu = local_Jac.template block<pressure_size, displacement_size>(
150 pressure_index, displacement_index);
151 auto J_uu = local_Jac.template block<displacement_size, displacement_size>(
152 displacement_index, displacement_index);
153 auto J_up = local_Jac.template block<displacement_size, pressure_size>(
154 displacement_index, pressure_index);
156 assembleBlockMatricesWithJacobian(t, dt, p, p_prev, u, u_prev, rhs_p, rhs_u,
157 J_pp, J_pu, J_uu, J_up);
163 ShapeFunctionPressure, GlobalDim>::
164 assembleBlockMatricesWithJacobian(
165 double const t,
double const dt,
166 Eigen::Ref<const Eigen::VectorXd>
const& p,
167 Eigen::Ref<const Eigen::VectorXd>
const& p_prev,
168 Eigen::Ref<const Eigen::VectorXd>
const& u,
169 Eigen::Ref<const Eigen::VectorXd>
const& u_prev,
170 Eigen::Ref<Eigen::VectorXd> rhs_p, Eigen::Ref<Eigen::VectorXd> rhs_u,
171 Eigen::Ref<Eigen::MatrixXd> J_pp, Eigen::Ref<Eigen::MatrixXd> J_pu,
172 Eigen::Ref<Eigen::MatrixXd> J_uu, Eigen::Ref<Eigen::MatrixXd> J_up)
174 assert(this->_element.getDimension() == GlobalDim);
177 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
181 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
184 typename ShapeMatricesTypeDisplacement::template MatrixType<
185 displacement_size, pressure_size>
186 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
187 displacement_size, pressure_size>::Zero(displacement_size,
190 auto const& gravity_vec = _process_data.specific_body_force;
197 auto const& medium = _process_data.media_map.getMedium(_element.getID());
198 auto const& liquid_phase = medium->phase(
"AqueousLiquid");
199 auto const& solid_phase = medium->phase(
"Solid");
202 medium->property(MPL::PropertyType::reference_temperature)
203 .template value<double>(variables, x_position, t, dt);
207 bool const has_storage_property =
208 medium->hasProperty(MPL::PropertyType::storage);
210 auto const B_dil_bar = getDilatationalBBarMatrix();
212 unsigned const n_integration_points = _ip_data.size();
213 for (
unsigned ip = 0; ip < n_integration_points; ip++)
215 auto& ip_data = _ip_data[ip];
216 auto const& ip_w = ip_data.integration_weight;
217 auto const& N_u = ip_data.N_u;
218 auto const& dNdx_u = ip_data.dNdx_u;
219 auto const& N_p = ip_data.N_p;
220 auto const& dNdx_p = ip_data.dNdx_p;
221 auto const& H_u = ip_data.H_u;
233 dNdx_u, N_u, B_dil_bar, x_coord, this->_is_axially_symmetric)
236 auto const& eps_prev = ip_data.eps_prev;
237 auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
238 auto& sigma_eff = ip_data.sigma_eff;
240 auto& eps = ip_data.eps;
241 auto& state = ip_data.material_state_variables;
244 solid_phase.property(MPL::PropertyType::density)
245 .template value<double>(variables, x_position, t, dt);
247 liquid_phase.property(MPL::PropertyType::density)
248 .template value<double>(variables, x_position, t, dt);
251 medium->property(MPL::PropertyType::biot_coefficient)
252 .template value<double>(variables, x_position, t, dt);
253 auto const porosity =
254 medium->property(MPL::PropertyType::porosity)
255 .template value<double>(variables, x_position, t, dt);
257 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
258 auto const& identity2 =
261 eps.noalias() = B * u;
274 auto&& solution = _ip_data[ip].solid_material.integrateStress(
275 variables_prev, variables, t, x_position, dt, *state);
279 OGS_FATAL(
"Computation of local constitutive relation failed.");
283 std::tie(sigma_eff, state, C) = std::move(*solution);
285 J_uu.noalias() += B.transpose() * C * B * ip_w;
287 rhs_u.noalias() -= B.transpose() * sigma_eff * ip_w;
288 rhs_u.noalias() -= -H_u.transpose() * rho * gravity_vec * ip_w;
294 if (!_process_data.deactivate_matrix_in_flow)
297 Kup.noalias() += B.transpose() * alpha * identity2 * N_p * ip_w;
301 liquid_phase.property(MPL::PropertyType::viscosity)
302 .template value<double>(variables, x_position, t, dt);
304 auto const k_over_mu =
306 medium->property(MPL::PropertyType::permeability)
307 .value(variables, x_position, t, dt)) /
313 ? medium->property(MPL::PropertyType::storage)
314 .template value<double>(variables, x_position, t, dt)
316 (liquid_phase.property(MPL::PropertyType::density)
317 .
template dValue<double>(
319 MPL::Variable::liquid_phase_pressure,
322 (alpha - porosity) * (1.0 - alpha) /
323 _ip_data[ip].solid_material.getBulkModulus(
326 auto q = ip_data.darcy_velocity.head(GlobalDim);
327 q.noalias() = -k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
329 laplace_p.noalias() +=
330 dNdx_p.transpose() * k_over_mu * dNdx_p * ip_w;
331 storage_p.noalias() += N_p.transpose() * S * N_p * ip_w;
334 dNdx_p.transpose() * rho_fr * k_over_mu * gravity_vec * ip_w;
339 J_up.noalias() -= Kup;
342 J_pp.noalias() += laplace_p + storage_p / dt;
345 J_pu.noalias() += Kup.transpose() / dt;
348 rhs_p.noalias() -= laplace_p * p + storage_p * (p - p_prev) / dt +
349 Kup.transpose() * (u - u_prev) / dt;
352 rhs_u.noalias() -= -Kup * p;
376 ShapeFunctionPressure, GlobalDim>::
377 postTimestepConcreteWithBlockVectors(
378 double const t,
double const dt,
379 Eigen::Ref<const Eigen::VectorXd>
const& p,
380 Eigen::Ref<const Eigen::VectorXd>
const& u)
386 auto const e_id = _element.getID();
390 KV sigma_avg = KV::Zero();
392 velocity_avg.setZero();
394 unsigned const n_integration_points = _ip_data.size();
396 auto const& medium = _process_data.media_map.getMedium(_element.getID());
397 auto const& liquid_phase = medium->phase(
"AqueousLiquid");
399 medium->property(MPL::PropertyType::reference_temperature)
400 .template value<double>(variables, x_position, t, dt);
404 auto const B_dil_bar = getDilatationalBBarMatrix();
406 for (
unsigned ip = 0; ip < n_integration_points; ip++)
408 auto& ip_data = _ip_data[ip];
410 auto const& eps_prev = ip_data.eps_prev;
411 auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
413 auto& eps = ip_data.eps;
414 auto& sigma_eff = ip_data.sigma_eff;
415 auto& state = ip_data.material_state_variables;
417 auto const& N_u = ip_data.N_u;
418 auto const& N_p = ip_data.N_p;
419 auto const& dNdx_u = ip_data.dNdx_u;
431 dNdx_u, N_u, B_dil_bar, x_coord, this->_is_axially_symmetric)
434 eps.noalias() = B * u;
446 auto&& solution = _ip_data[ip].solid_material.integrateStress(
447 variables_prev, variables, t, x_position, dt, *state);
451 OGS_FATAL(
"Computation of local constitutive relation failed.");
455 std::tie(sigma_eff, state, C) = std::move(*solution);
457 sigma_avg += ip_data.sigma_eff;
459 if (!_process_data.deactivate_matrix_in_flow)
463 liquid_phase.property(MPL::PropertyType::density)
464 .template value<double>(variables, x_position, t, dt);
467 liquid_phase.property(MPL::PropertyType::viscosity)
468 .template value<double>(variables, x_position, t, dt);
471 medium->property(MPL::PropertyType::permeability)
472 .value(variables, x_position, t, dt))
477 auto const& gravity_vec = _process_data.specific_body_force;
478 auto const& dNdx_p = ip_data.dNdx_p;
480 ip_data.darcy_velocity.head(GlobalDim).noalias() =
481 -k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
482 velocity_avg.noalias() += ip_data.darcy_velocity.head(GlobalDim);
486 sigma_avg /= n_integration_points;
487 velocity_avg /= n_integration_points;
490 &(*_process_data.element_stresses)[e_id * KV::RowsAtCompileTime]) =
493 Eigen::Map<GlobalDimVector>(
494 &(*_process_data.element_velocities)[e_id * GlobalDim]) = velocity_avg;
497 ShapeFunctionPressure,
typename ShapeFunctionDisplacement::MeshElement,
498 GlobalDim>(_element, _is_axially_symmetric, p,
499 *_process_data.mesh_prop_nodal_p);