32 ShapeFunctionPressure, GlobalDim>::
33 HydroMechanicsLocalAssemblerMatrix(
35 std::size_t
const n_variables,
37 std::vector<unsigned>
const& dofIndex_to_localIndex,
39 bool const is_axially_symmetric,
42 e, is_axially_symmetric,
43 (n_variables - 1) * ShapeFunctionDisplacement::NPOINTS * GlobalDim +
44 ShapeFunctionPressure::NPOINTS,
45 dofIndex_to_localIndex),
46 _process_data(process_data)
48 unsigned const n_integration_points =
51 _ip_data.reserve(n_integration_points);
54 auto const shape_matrices_u =
57 e, is_axially_symmetric, integration_method);
59 auto const shape_matrices_p =
62 e, is_axially_symmetric, integration_method);
69 for (
unsigned ip = 0; ip < n_integration_points; ip++)
73 _ip_data.emplace_back(solid_material);
75 auto const& sm_u = shape_matrices_u[ip];
76 auto const& sm_p = shape_matrices_p[ip];
78 ip_data.integration_weight =
79 sm_u.detJ * sm_u.integralMeasure *
81 ip_data.darcy_velocity.setZero();
84 ip_data.dNdx_u = sm_u.dNdx;
86 for (
int i = 0; i < GlobalDim; ++i)
91 .noalias() = ip_data.N_u;
95 ip_data.dNdx_p = sm_p.dNdx;
106 auto const initial_effective_stress =
110 ip_data.sigma_eff[i] = initial_effective_stress[i];
111 ip_data.sigma_eff_prev[i] = initial_effective_stress[i];
119 ShapeFunctionPressure, GlobalDim>::
120 assembleWithJacobianConcrete(
double const t,
double const dt,
121 Eigen::VectorXd
const& local_x,
122 Eigen::VectorXd
const& local_x_prev,
123 Eigen::VectorXd& local_rhs,
124 Eigen::MatrixXd& local_Jac)
126 auto p =
const_cast<Eigen::VectorXd&
>(local_x).segment(pressure_index,
128 auto p_prev =
const_cast<Eigen::VectorXd&
>(local_x_prev)
129 .segment(pressure_index, pressure_size);
131 if (_process_data.deactivate_matrix_in_flow)
133 setPressureOfInactiveNodes(t, p);
136 auto u = local_x.segment(displacement_index, displacement_size);
137 auto u_prev = local_x_prev.segment(displacement_index, displacement_size);
139 auto rhs_p = local_rhs.template segment<pressure_size>(pressure_index);
141 local_rhs.template segment<displacement_size>(displacement_index);
143 auto J_pp = local_Jac.template block<pressure_size, pressure_size>(
144 pressure_index, pressure_index);
145 auto J_pu = local_Jac.template block<pressure_size, displacement_size>(
146 pressure_index, displacement_index);
147 auto J_uu = local_Jac.template block<displacement_size, displacement_size>(
148 displacement_index, displacement_index);
149 auto J_up = local_Jac.template block<displacement_size, pressure_size>(
150 displacement_index, pressure_index);
152 assembleBlockMatricesWithJacobian(t, dt, p, p_prev, u, u_prev, rhs_p, rhs_u,
153 J_pp, J_pu, J_uu, J_up);
159 ShapeFunctionPressure, GlobalDim>::
160 assembleBlockMatricesWithJacobian(
161 double const t,
double const dt,
162 Eigen::Ref<const Eigen::VectorXd>
const& p,
163 Eigen::Ref<const Eigen::VectorXd>
const& p_prev,
164 Eigen::Ref<const Eigen::VectorXd>
const& u,
165 Eigen::Ref<const Eigen::VectorXd>
const& u_prev,
166 Eigen::Ref<Eigen::VectorXd> rhs_p, Eigen::Ref<Eigen::VectorXd> rhs_u,
167 Eigen::Ref<Eigen::MatrixXd> J_pp, Eigen::Ref<Eigen::MatrixXd> J_pu,
168 Eigen::Ref<Eigen::MatrixXd> J_uu, Eigen::Ref<Eigen::MatrixXd> J_up)
170 assert(this->_element.getDimension() == GlobalDim);
173 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
177 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
180 typename ShapeMatricesTypeDisplacement::template MatrixType<
181 displacement_size, pressure_size>
182 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
183 displacement_size, pressure_size>::Zero(displacement_size,
186 auto const& gravity_vec = _process_data.specific_body_force;
193 unsigned const n_integration_points = _ip_data.size();
194 for (
unsigned ip = 0; ip < n_integration_points; ip++)
198 auto& ip_data = _ip_data[ip];
199 auto const& ip_w = ip_data.integration_weight;
200 auto const& N_u = ip_data.N_u;
201 auto const& dNdx_u = ip_data.dNdx_u;
202 auto const& N_p = ip_data.N_p;
203 auto const& dNdx_p = ip_data.dNdx_p;
204 auto const& H_u = ip_data.H_u;
212 ShapeFunctionDisplacement::NPOINTS,
214 dNdx_u, N_u, x_coord, _is_axially_symmetric);
216 auto const& eps_prev = ip_data.eps_prev;
217 auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
218 auto& sigma_eff = ip_data.sigma_eff;
220 auto& eps = ip_data.eps;
221 auto& state = ip_data.material_state_variables;
223 auto const alpha = _process_data.biot_coefficient(t, x_position)[0];
224 auto const rho_sr = _process_data.solid_density(t, x_position)[0];
225 auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
226 auto const porosity = _process_data.porosity(t, x_position)[0];
228 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
229 auto const& identity2 =
232 eps.noalias() = B * u;
243 variables_prev.
temperature = _process_data.reference_temperature;
245 auto&& solution = _ip_data[ip].solid_material.integrateStress(
246 variables_prev, variables, t, x_position, dt, *state);
250 OGS_FATAL(
"Computation of local constitutive relation failed.");
254 std::tie(sigma_eff, state, C) = std::move(*solution);
256 J_uu.noalias() += B.transpose() * C * B * ip_w;
258 rhs_u.noalias() -= B.transpose() * sigma_eff * ip_w;
259 rhs_u.noalias() -= -H_u.transpose() * rho * gravity_vec * ip_w;
265 if (!_process_data.deactivate_matrix_in_flow)
268 Kup.noalias() += B.transpose() * alpha * identity2 * N_p * ip_w;
270 double const k_over_mu =
271 _process_data.intrinsic_permeability(t, x_position)[0] /
272 _process_data.fluid_viscosity(t, x_position)[0];
273 double const S = _process_data.specific_storage(t, x_position)[0];
275 auto q = ip_data.darcy_velocity.head(GlobalDim);
276 q.noalias() = -k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
278 laplace_p.noalias() +=
279 dNdx_p.transpose() * k_over_mu * dNdx_p * ip_w;
280 storage_p.noalias() += N_p.transpose() * S * N_p * ip_w;
283 dNdx_p.transpose() * rho_fr * k_over_mu * gravity_vec * ip_w;
288 J_up.noalias() -= Kup;
291 J_pp.noalias() += laplace_p + storage_p / dt;
294 J_pu.noalias() += Kup.transpose() / dt;
297 rhs_p.noalias() -= laplace_p * p + storage_p * (p - p_prev) / dt +
298 Kup.transpose() * (u - u_prev) / dt;
301 rhs_u.noalias() -= -Kup * p;
325 ShapeFunctionPressure, GlobalDim>::
326 postTimestepConcreteWithBlockVectors(
327 double const t,
double const dt,
328 Eigen::Ref<const Eigen::VectorXd>
const& p,
329 Eigen::Ref<const Eigen::VectorXd>
const& u)
335 auto const e_id = _element.getID();
339 KV sigma_avg = KV::Zero();
341 velocity_avg.setZero();
343 unsigned const n_integration_points = _ip_data.size();
344 for (
unsigned ip = 0; ip < n_integration_points; ip++)
348 auto& ip_data = _ip_data[ip];
350 auto const& eps_prev = ip_data.eps_prev;
351 auto const& sigma_eff_prev = ip_data.sigma_eff_prev;
353 auto& eps = ip_data.eps;
354 auto& sigma_eff = ip_data.sigma_eff;
355 auto& state = ip_data.material_state_variables;
357 auto const& N_u = ip_data.N_u;
358 auto const& dNdx_u = ip_data.dNdx_u;
366 ShapeFunctionDisplacement::NPOINTS,
368 dNdx_u, N_u, x_coord, _is_axially_symmetric);
370 eps.noalias() = B * u;
381 variables_prev.
temperature = _process_data.reference_temperature;
383 auto&& solution = _ip_data[ip].solid_material.integrateStress(
384 variables_prev, variables, t, x_position, dt, *state);
388 OGS_FATAL(
"Computation of local constitutive relation failed.");
392 std::tie(sigma_eff, state, C) = std::move(*solution);
394 sigma_avg += ip_data.sigma_eff;
396 if (!_process_data.deactivate_matrix_in_flow)
399 double const k_over_mu =
400 _process_data.intrinsic_permeability(t, x_position)[0] /
401 _process_data.fluid_viscosity(t, x_position)[0];
402 auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
403 auto const& gravity_vec = _process_data.specific_body_force;
404 auto const& dNdx_p = ip_data.dNdx_p;
406 ip_data.darcy_velocity.head(GlobalDim).noalias() =
407 -k_over_mu * (dNdx_p * p + rho_fr * gravity_vec);
408 velocity_avg += ip_data.darcy_velocity.head(GlobalDim);
412 sigma_avg /= n_integration_points;
413 velocity_avg /= n_integration_points;
416 &(*_process_data.element_stresses)[e_id * KV::RowsAtCompileTime]) =
419 Eigen::Map<GlobalDimVector>(
420 &(*_process_data.element_velocities)[e_id * GlobalDim]) = velocity_avg;
423 ShapeFunctionPressure,
typename ShapeFunctionDisplacement::MeshElement,
424 GlobalDim>(_element, _is_axially_symmetric, p,
425 *_process_data.mesh_prop_nodal_p);
479 ShapeFunctionDisplacement, ShapeFunctionPressure, GlobalDim>::
480 getIntPtDarcyVelocity(
482 std::vector<GlobalVector*>
const& ,
483 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
484 std::vector<double>& cache)
const
486 unsigned const n_integration_points = _ip_data.size();
490 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
491 cache, GlobalDim, n_integration_points);
493 for (
unsigned ip = 0; ip < n_integration_points; ip++)
495 cache_matrix.col(ip).noalias() = _ip_data[ip].darcy_velocity;