44 ShapeFunctionPressure, DisplacementDim>::
45 HydroMechanicsLocalAssemblerFracture(
48 std::vector<unsigned>
const& dofIndex_to_localIndex,
50 bool const is_axially_symmetric,
53 e, is_axially_symmetric, integration_method,
54 ShapeFunctionDisplacement::NPOINTS * DisplacementDim +
55 ShapeFunctionPressure::NPOINTS,
56 dofIndex_to_localIndex),
57 _process_data(process_data)
61 unsigned const n_integration_points =
64 _ip_data.reserve(n_integration_points);
67 auto const shape_matrices_u =
70 DisplacementDim>(e, is_axially_symmetric,
73 auto const shape_matrices_p =
76 e, is_axially_symmetric, integration_method);
79 auto frac_id =
_process_data.map_materialID_to_fractureID[mat_id];
85 aperture0_node_values =
89 for (
unsigned ip = 0; ip < n_integration_points; ip++)
92 auto const& sm_u = shape_matrices_u[ip];
93 auto const& sm_p = shape_matrices_p[ip];
102 ip_data.integration_weight =
103 sm_u.detJ * sm_u.integralMeasure *
108 ShapeFunctionDisplacement::NPOINTS * DisplacementDim);
110 DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
113 ip_data.N_p = sm_p.N;
114 ip_data.dNdx_p = sm_p.dNdx;
119 ip_data.w.setZero(DisplacementDim);
120 ip_data.sigma_eff.setZero(DisplacementDim);
123 ip_data.w_prev.resize(DisplacementDim);
124 ip_data.sigma_eff_prev.resize(DisplacementDim);
126 ip_data.C.resize(DisplacementDim, DisplacementDim);
128 ip_data.aperture0 = aperture0_node_values.dot(sm_u.N);
129 ip_data.aperture = ip_data.aperture0;
131 auto const initial_effective_stress =
132 _process_data.initial_fracture_effective_stress(0, x_position);
133 for (
int i = 0; i < DisplacementDim; i++)
135 ip_data.sigma_eff[i] = initial_effective_stress[i];
136 ip_data.sigma_eff_prev[i] = initial_effective_stress[i];
144 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
145 assembleWithJacobianConcrete(
double const t,
double const dt,
146 Eigen::VectorXd
const& local_x,
147 Eigen::VectorXd
const& local_x_prev,
148 Eigen::VectorXd& local_b,
149 Eigen::MatrixXd& local_J)
151 auto const p = local_x.segment(pressure_index, pressure_size);
152 auto const p_prev = local_x_prev.segment(pressure_index, pressure_size);
153 auto const g = local_x.segment(displacement_index, displacement_size);
155 local_x_prev.segment(displacement_index, displacement_size);
157 auto rhs_p = local_b.segment(pressure_index, pressure_size);
158 auto rhs_g = local_b.segment(displacement_index, displacement_size);
159 auto J_pp = local_J.block(pressure_index, pressure_index, pressure_size,
161 auto J_pg = local_J.block(pressure_index, displacement_index, pressure_size,
163 auto J_gp = local_J.block(displacement_index, pressure_index,
164 displacement_size, pressure_size);
165 auto J_gg = local_J.block(displacement_index, displacement_index,
166 displacement_size, displacement_size);
168 assembleBlockMatricesWithJacobian(t, dt, p, p_prev, g, g_prev, rhs_p, rhs_g,
169 J_pp, J_pg, J_gg, J_gp);
175 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
176 assembleBlockMatricesWithJacobian(
177 double const t,
double const dt,
178 Eigen::Ref<const Eigen::VectorXd>
const& p,
179 Eigen::Ref<const Eigen::VectorXd>
const& p_prev,
180 Eigen::Ref<const Eigen::VectorXd>
const& g,
181 Eigen::Ref<const Eigen::VectorXd>
const& g_prev,
182 Eigen::Ref<Eigen::VectorXd> rhs_p, Eigen::Ref<Eigen::VectorXd> rhs_g,
183 Eigen::Ref<Eigen::MatrixXd> J_pp, Eigen::Ref<Eigen::MatrixXd> J_pg,
184 Eigen::Ref<Eigen::MatrixXd> J_gg, Eigen::Ref<Eigen::MatrixXd> J_gp)
186 auto const& R = _fracture_property->R;
190 auto constexpr index_normal = DisplacementDim - 1;
193 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
197 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
200 typename ShapeMatricesTypeDisplacement::template MatrixType<
201 displacement_size, pressure_size>
202 Kgp = ShapeMatricesTypeDisplacement::template MatrixType<
203 displacement_size, pressure_size>::Zero(displacement_size,
207 Eigen::MatrixXd
const global2local_rotation =
208 R.template topLeftCorner<ShapeFunctionPressure::DIM, DisplacementDim>();
211 global2local_rotation *
212 _process_data.specific_body_force;
218 auto const& medium = _process_data.media_map.getMedium(_element.getID());
219 auto const& liquid_phase = medium->phase(
"AqueousLiquid");
221 medium->property(MPL::PropertyType::reference_temperature)
222 .template value<double>(variables, x_position, t, dt);
225 unsigned const n_integration_points = _ip_data.size();
226 for (
unsigned ip = 0; ip < n_integration_points; ip++)
228 auto& ip_data = _ip_data[ip];
229 auto const& ip_w = ip_data.integration_weight;
230 auto const& N_p = ip_data.N_p;
231 auto const& dNdx_p = ip_data.dNdx_p;
232 auto const& H_g = ip_data.H_u;
233 auto const& identity2 =
237 std::nullopt, _element.getID(),
243 auto& mat = ip_data.fracture_material;
244 auto& effective_stress = ip_data.sigma_eff;
245 auto const& effective_stress_prev = ip_data.sigma_eff_prev;
247 auto const& w_prev = ip_data.w_prev;
249 auto& state = *ip_data.material_state_variables;
250 auto& b_m = ip_data.aperture;
253 liquid_phase.property(MPL::PropertyType::density)
254 .template value<double>(variables, x_position, t, dt);
258 medium->property(MPL::PropertyType::biot_coefficient)
259 .template value<double>(variables, x_position, t, dt);
262 medium->property(MPL::PropertyType::storage)
263 .template value<double>(variables, x_position, t, dt);
266 liquid_phase.property(MPL::PropertyType::viscosity)
267 .template value<double>(variables, x_position, t, dt);
270 w.noalias() = R * H_g * g;
273 b_m = ip_data.aperture0 + w[index_normal];
277 "Element {:d}, gp {:d}: Fracture aperture is {:g}, but it must "
279 "non-negative. Setting it to zero.",
280 _element.getID(), ip, b_m);
284 auto const initial_effective_stress =
285 _process_data.initial_fracture_effective_stress(0, x_position);
287 Eigen::Map<typename HMatricesType::ForceVectorType const>
const stress0(
288 initial_effective_stress.data(), initial_effective_stress.size());
291 mat.computeConstitutiveRelation(
292 t, x_position, ip_data.aperture0, stress0, w_prev, w,
293 effective_stress_prev, effective_stress, C, state);
299 H_g.transpose() * R.transpose() * effective_stress * ip_w;
300 J_gg.noalias() += H_g.transpose() * R.transpose() * C * R * H_g * ip_w;
306 H_g.transpose() * R.transpose() * alpha * identity2 * N_p * ip_w;
314 auto const permeability =
315 medium->property(MPL::PropertyType::permeability)
316 .value(variables, x_position, t, dt);
318 auto& k = ip_data.permeability;
319 k = std::get<double>(permeability);
320 double const k_over_mu = k / mu;
321 storage_p.noalias() += N_p.transpose() * b_m * S * N_p * ip_w;
322 laplace_p.noalias() +=
323 dNdx_p.transpose() * b_m * k_over_mu * dNdx_p * ip_w;
325 dNdx_p.transpose() * b_m * k_over_mu * rho_fr * gravity_vec * ip_w;
331 Eigen::Matrix<double, 1, displacement_size>
const mT_R_Hg =
332 identity2.transpose() * R * H_g;
334 ip_data.darcy_velocity = -k_over_mu * grad_head;
336 N_p.transpose() * S * N_p * (p - p_prev) / dt * mT_R_Hg * ip_w;
339 double const dk_db_over_mu =
340 medium->property(MPL::PropertyType::permeability)
341 .template dValue<double>(variables,
342 MPL::Variable::fracture_aperture,
346 dNdx_p.transpose() * k_over_mu * grad_head * mT_R_Hg * ip_w;
347 J_pg.noalias() += dNdx_p.transpose() * b_m * dk_db_over_mu * grad_head *
352 J_gp.noalias() -= Kgp;
355 J_pp.noalias() += laplace_p + storage_p / dt;
358 J_pg.noalias() += Kgp.transpose() / dt;
361 rhs_p.noalias() -= laplace_p * p + storage_p * (p - p_prev) / dt +
362 Kgp.transpose() * (g - g_prev) / dt;
365 rhs_g.noalias() -= -Kgp * p;
372 DisplacementDim>::postTimestepConcreteWithVector(
const double t,
374 Eigen::VectorXd
const&
377 auto const nodal_g = local_x.segment(displacement_index, displacement_size);
379 auto const& R = _fracture_property->R;
382 auto constexpr index_normal = DisplacementDim - 1;
385 auto const e_id = _element.getID();
388 unsigned const n_integration_points = _ip_data.size();
389 for (
unsigned ip = 0; ip < n_integration_points; ip++)
391 auto& ip_data = _ip_data[ip];
392 auto const& H_g = ip_data.H_u;
393 auto& mat = ip_data.fracture_material;
394 auto& effective_stress = ip_data.sigma_eff;
395 auto const& effective_stress_prev = ip_data.sigma_eff_prev;
397 auto const& w_prev = ip_data.w_prev;
399 auto& state = *ip_data.material_state_variables;
400 auto& b_m = ip_data.aperture;
407 _element, ip_data.N_p))};
410 w.noalias() = R * H_g * nodal_g;
413 b_m = ip_data.aperture0 + w[index_normal];
417 "Element {:d}, gp {:d}: Fracture aperture is {:g}, but it is "
418 "expected to be non-negative. Setting it to zero now.",
419 _element.getID(), ip, b_m);
423 auto const initial_effective_stress =
424 _process_data.initial_fracture_effective_stress(0, x_position);
426 Eigen::Map<typename HMatricesType::ForceVectorType const>
const stress0(
427 initial_effective_stress.data(), initial_effective_stress.size());
430 mat.computeConstitutiveRelation(
431 t, x_position, ip_data.aperture0, stress0, w_prev, w,
432 effective_stress_prev, effective_stress, C, state);
438 HMatricesType::ForceVectorType::Zero(DisplacementDim);
440 HMatricesType::ForceVectorType::Zero(DisplacementDim);
443 double ele_Fs = -std::numeric_limits<double>::max();
444 for (
auto const& ip : _ip_data)
446 ele_b += ip.aperture;
447 ele_k += ip.permeability;
449 ele_sigma_eff += ip.sigma_eff;
450 ele_velocity += ip.darcy_velocity;
452 ele_Fs, ip.material_state_variables->getShearYieldFunctionValue());
454 ele_b /=
static_cast<double>(n_integration_points);
455 ele_k /=
static_cast<double>(n_integration_points);
456 ele_w /=
static_cast<double>(n_integration_points);
457 ele_sigma_eff /=
static_cast<double>(n_integration_points);
458 ele_velocity /=
static_cast<double>(n_integration_points);
459 auto const element_id = _element.getID();
460 (*_process_data.mesh_prop_b)[element_id] = ele_b;
461 (*_process_data.mesh_prop_k_f)[element_id] = ele_k;
463 Eigen::Map<GlobalDimVectorType>(
464 &(*_process_data.element_fracture_stresses)[e_id * DisplacementDim]) =
467 Eigen::Map<GlobalDimVectorType>(
468 &(*_process_data.element_fracture_velocities)[e_id * DisplacementDim]) =
471 Eigen::Map<GlobalDimVectorType>(
472 &(*_process_data.element_local_jumps)[e_id * DisplacementDim]) = ele_w;
474 (*_process_data.mesh_prop_fracture_shear_failure)[element_id] = ele_Fs;