44 ShapeFunctionPressure, GlobalDim>::
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 * GlobalDim +
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 e, is_axially_symmetric, integration_method);
72 auto const shape_matrices_p =
75 e, is_axially_symmetric, integration_method);
82 aperture0_node_values = frac_prop.aperture0.getNodalValuesOnElement(
87 for (
unsigned ip = 0; ip < n_integration_points; ip++)
90 auto const& sm_u = shape_matrices_u[ip];
91 auto const& sm_p = shape_matrices_p[ip];
93 ip_data.integration_weight =
94 sm_u.detJ * sm_u.integralMeasure *
97 ip_data.H_u.setZero(GlobalDim,
98 ShapeFunctionDisplacement::NPOINTS * GlobalDim);
100 GlobalDim, ShapeFunctionDisplacement::NPOINTS,
103 ip_data.N_p = sm_p.N;
104 ip_data.dNdx_p = sm_p.dNdx;
109 ip_data.w.setZero(GlobalDim);
110 ip_data.sigma_eff.setZero(GlobalDim);
113 ip_data.w_prev.resize(GlobalDim);
114 ip_data.sigma_eff_prev.resize(GlobalDim);
116 ip_data.C.resize(GlobalDim, GlobalDim);
118 ip_data.aperture0 = aperture0_node_values.dot(sm_u.N);
119 ip_data.aperture = ip_data.aperture0;
121 auto const initial_effective_stress =
122 _process_data.initial_fracture_effective_stress(0, x_position);
123 for (
int i = 0; i < GlobalDim; i++)
125 ip_data.sigma_eff[i] = initial_effective_stress[i];
126 ip_data.sigma_eff_prev[i] = initial_effective_stress[i];
134 ShapeFunctionPressure, GlobalDim>::
135 assembleWithJacobianConcrete(
double const t,
double const dt,
136 Eigen::VectorXd
const& local_x,
137 Eigen::VectorXd
const& local_x_prev,
138 Eigen::VectorXd& local_b,
139 Eigen::MatrixXd& local_J)
141 auto const p = local_x.segment(pressure_index, pressure_size);
142 auto const p_prev = local_x_prev.segment(pressure_index, pressure_size);
143 auto const g = local_x.segment(displacement_index, displacement_size);
145 local_x_prev.segment(displacement_index, displacement_size);
147 auto rhs_p = local_b.segment(pressure_index, pressure_size);
148 auto rhs_g = local_b.segment(displacement_index, displacement_size);
149 auto J_pp = local_J.block(pressure_index, pressure_index, pressure_size,
151 auto J_pg = local_J.block(pressure_index, displacement_index, pressure_size,
153 auto J_gp = local_J.block(displacement_index, pressure_index,
154 displacement_size, pressure_size);
155 auto J_gg = local_J.block(displacement_index, displacement_index,
156 displacement_size, displacement_size);
158 assembleBlockMatricesWithJacobian(t, dt, p, p_prev, g, g_prev, rhs_p, rhs_g,
159 J_pp, J_pg, J_gg, J_gp);
165 ShapeFunctionPressure, GlobalDim>::
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& g,
171 Eigen::Ref<const Eigen::VectorXd>
const& g_prev,
172 Eigen::Ref<Eigen::VectorXd> rhs_p, Eigen::Ref<Eigen::VectorXd> rhs_g,
173 Eigen::Ref<Eigen::MatrixXd> J_pp, Eigen::Ref<Eigen::MatrixXd> J_pg,
174 Eigen::Ref<Eigen::MatrixXd> J_gg, Eigen::Ref<Eigen::MatrixXd> J_gp)
176 auto const& frac_prop = *_process_data.fracture_property;
177 auto const& R = frac_prop.R;
181 auto constexpr index_normal = GlobalDim - 1;
184 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
188 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
191 typename ShapeMatricesTypeDisplacement::template MatrixType<
192 displacement_size, pressure_size>
193 Kgp = ShapeMatricesTypeDisplacement::template MatrixType<
194 displacement_size, pressure_size>::Zero(displacement_size,
198 Eigen::MatrixXd
const global2local_rotation =
199 R.template topLeftCorner<ShapeFunctionPressure::DIM, GlobalDim>();
202 global2local_rotation *
203 _process_data.specific_body_force;
209 auto const& medium = _process_data.media_map.getMedium(_element.getID());
210 auto const& liquid_phase = medium->phase(
"AqueousLiquid");
212 medium->property(MPL::PropertyType::reference_temperature)
213 .template value<double>(variables, x_position, t, dt);
216 unsigned const n_integration_points = _ip_data.size();
217 for (
unsigned ip = 0; ip < n_integration_points; ip++)
219 auto& ip_data = _ip_data[ip];
220 auto const& ip_w = ip_data.integration_weight;
221 auto const& N_p = ip_data.N_p;
222 auto const& dNdx_p = ip_data.dNdx_p;
223 auto const& H_g = ip_data.H_u;
224 auto const& identity2 =
227 auto& mat = ip_data.fracture_material;
228 auto& effective_stress = ip_data.sigma_eff;
229 auto const& effective_stress_prev = ip_data.sigma_eff_prev;
231 auto const& w_prev = ip_data.w_prev;
233 auto& state = *ip_data.material_state_variables;
234 auto& b_m = ip_data.aperture;
237 liquid_phase.property(MPL::PropertyType::density)
238 .template value<double>(variables, x_position, t, dt);
242 medium->property(MPL::PropertyType::biot_coefficient)
243 .template value<double>(variables, x_position, t, dt);
246 medium->property(MPL::PropertyType::storage)
247 .template value<double>(variables, x_position, t, dt);
250 liquid_phase.property(MPL::PropertyType::viscosity)
251 .template value<double>(variables, x_position, t, dt);
254 w.noalias() = R * H_g * g;
257 b_m = ip_data.aperture0 + w[index_normal];
261 "Element {:d}, gp {:d}: Fracture aperture is {:g}, but it must "
263 "non-negative. Setting it to zero.",
264 _element.getID(), ip, b_m);
268 auto const initial_effective_stress =
269 _process_data.initial_fracture_effective_stress(0, x_position);
271 Eigen::Map<typename HMatricesType::ForceVectorType const>
const stress0(
272 initial_effective_stress.data(), initial_effective_stress.size());
275 mat.computeConstitutiveRelation(
276 t, x_position, ip_data.aperture0, stress0, w_prev, w,
277 effective_stress_prev, effective_stress, C, state);
283 H_g.transpose() * R.transpose() * effective_stress * ip_w;
284 J_gg.noalias() += H_g.transpose() * R.transpose() * C * R * H_g * ip_w;
290 H_g.transpose() * R.transpose() * alpha * identity2 * N_p * ip_w;
298 auto const permeability =
299 medium->property(MPL::PropertyType::permeability)
300 .value(variables, x_position, t, dt);
302 auto& k = ip_data.permeability;
303 k = std::get<double>(permeability);
304 double const k_over_mu = k / mu;
305 storage_p.noalias() += N_p.transpose() * b_m * S * N_p * ip_w;
306 laplace_p.noalias() +=
307 dNdx_p.transpose() * b_m * k_over_mu * dNdx_p * ip_w;
309 dNdx_p.transpose() * b_m * k_over_mu * rho_fr * gravity_vec * ip_w;
315 Eigen::Matrix<double, 1, displacement_size>
const mT_R_Hg =
316 identity2.transpose() * R * H_g;
318 ip_data.darcy_velocity = -k_over_mu * grad_head;
320 N_p.transpose() * S * N_p * (p - p_prev) / dt * mT_R_Hg * ip_w;
323 double const dk_db_over_mu =
324 medium->property(MPL::PropertyType::permeability)
325 .template dValue<double>(variables,
326 MPL::Variable::fracture_aperture,
330 dNdx_p.transpose() * k_over_mu * grad_head * mT_R_Hg * ip_w;
331 J_pg.noalias() += dNdx_p.transpose() * b_m * dk_db_over_mu * grad_head *
336 J_gp.noalias() -= Kgp;
339 J_pp.noalias() += laplace_p + storage_p / dt;
342 J_pg.noalias() += Kgp.transpose() / dt;
345 rhs_p.noalias() -= laplace_p * p + storage_p * (p - p_prev) / dt +
346 Kgp.transpose() * (g - g_prev) / dt;
349 rhs_g.noalias() -= -Kgp * p;
356 GlobalDim>::postTimestepConcreteWithVector(
const double t,
358 Eigen::VectorXd
const& local_x)
360 auto const nodal_g = local_x.segment(displacement_index, displacement_size);
362 auto const& frac_prop = *_process_data.fracture_property;
363 auto const& R = frac_prop.R;
366 auto constexpr index_normal = GlobalDim - 1;
369 auto const e_id = _element.getID();
372 unsigned const n_integration_points = _ip_data.size();
373 for (
unsigned ip = 0; ip < n_integration_points; ip++)
375 auto& ip_data = _ip_data[ip];
376 auto const& H_g = ip_data.H_u;
377 auto& mat = ip_data.fracture_material;
378 auto& effective_stress = ip_data.sigma_eff;
379 auto const& effective_stress_prev = ip_data.sigma_eff_prev;
381 auto const& w_prev = ip_data.w_prev;
383 auto& state = *ip_data.material_state_variables;
384 auto& b_m = ip_data.aperture;
387 w.noalias() = R * H_g * nodal_g;
390 b_m = ip_data.aperture0 + w[index_normal];
394 "Element {:d}, gp {:d}: Fracture aperture is {:g}, but it is "
395 "expected to be non-negative. Setting it to zero now.",
396 _element.getID(), ip, b_m);
400 auto const initial_effective_stress =
401 _process_data.initial_fracture_effective_stress(0, x_position);
403 Eigen::Map<typename HMatricesType::ForceVectorType const>
const stress0(
404 initial_effective_stress.data(), initial_effective_stress.size());
407 mat.computeConstitutiveRelation(
408 t, x_position, ip_data.aperture0, stress0, w_prev, w,
409 effective_stress_prev, effective_stress, C, state);
415 HMatricesType::ForceVectorType::Zero(GlobalDim);
417 HMatricesType::ForceVectorType::Zero(GlobalDim);
420 double ele_Fs = -std::numeric_limits<double>::max();
421 for (
auto const& ip : _ip_data)
423 ele_b += ip.aperture;
424 ele_k += ip.permeability;
426 ele_sigma_eff += ip.sigma_eff;
427 ele_velocity += ip.darcy_velocity;
429 ele_Fs, ip.material_state_variables->getShearYieldFunctionValue());
431 ele_b /=
static_cast<double>(n_integration_points);
432 ele_k /=
static_cast<double>(n_integration_points);
433 ele_w /=
static_cast<double>(n_integration_points);
434 ele_sigma_eff /=
static_cast<double>(n_integration_points);
435 ele_velocity /=
static_cast<double>(n_integration_points);
436 auto const element_id = _element.getID();
437 (*_process_data.mesh_prop_b)[element_id] = ele_b;
438 (*_process_data.mesh_prop_k_f)[element_id] = ele_k;
440 Eigen::Map<GlobalDimVectorType>(
441 &(*_process_data.element_fracture_stresses)[e_id * GlobalDim]) =
444 Eigen::Map<GlobalDimVectorType>(
445 &(*_process_data.element_fracture_velocities)[e_id * GlobalDim]) =
448 Eigen::Map<GlobalDimVectorType>(
449 &(*_process_data.element_local_jumps)[e_id * GlobalDim]) = ele_w;
451 (*_process_data.mesh_prop_fracture_shear_failure)[element_id] = ele_Fs;