22 namespace HydroMechanics
24 template <
int GlobalDim,
typename RotationMatrix>
28 using M = Eigen::Matrix<double, GlobalDim, GlobalDim>;
30 tensor.diagonal().head(GlobalDim - 1).setConstant(value);
31 return (R.transpose() * tensor * R).eval();
34 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
35 typename IntegrationMethod,
int GlobalDim>
36 HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
37 ShapeFunctionPressure, IntegrationMethod,
39 HydroMechanicsLocalAssemblerFracture(
42 std::vector<unsigned>
const& dofIndex_to_localIndex,
43 bool const is_axially_symmetric,
44 unsigned const integration_order,
47 e, is_axially_symmetric,
48 ShapeFunctionDisplacement::NPOINTS * GlobalDim +
49 ShapeFunctionPressure::NPOINTS,
50 dofIndex_to_localIndex),
51 _process_data(process_data)
55 IntegrationMethod integration_method(integration_order);
56 unsigned const n_integration_points =
57 integration_method.getNumberOfPoints();
59 _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 aperture0_node_values = frac_prop.aperture0.getNodalValuesOnElement(
81 for (
unsigned ip = 0; ip < n_integration_points; ip++)
86 auto const& sm_u = shape_matrices_u[ip];
87 auto const& sm_p = shape_matrices_p[ip];
89 ip_data.integration_weight =
90 sm_u.detJ * sm_u.integralMeasure *
91 integration_method.getWeightedPoint(ip).getWeight();
93 ip_data.H_u.setZero(GlobalDim,
94 ShapeFunctionDisplacement::NPOINTS * GlobalDim);
96 GlobalDim, ShapeFunctionDisplacement::NPOINTS,
100 ip_data.dNdx_p = sm_p.dNdx;
103 ip_data.w.setZero(GlobalDim);
104 ip_data.sigma_eff.setZero(GlobalDim);
107 ip_data.w_prev.resize(GlobalDim);
108 ip_data.sigma_eff_prev.resize(GlobalDim);
110 ip_data.C.resize(GlobalDim, GlobalDim);
112 ip_data.aperture0 = aperture0_node_values.dot(sm_u.N);
113 ip_data.aperture = ip_data.aperture0;
115 ip_data.permeability_state =
116 frac_prop.permeability_model->getNewState();
118 auto const initial_effective_stress =
119 _process_data.initial_fracture_effective_stress(0, x_position);
120 for (
int i = 0; i < GlobalDim; i++)
122 ip_data.sigma_eff[i] = initial_effective_stress[i];
123 ip_data.sigma_eff_prev[i] = initial_effective_stress[i];
128 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
129 typename IntegrationMethod,
int GlobalDim>
131 ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
132 GlobalDim>::assembleWithJacobianConcrete(
double const t,
double const dt,
133 Eigen::VectorXd
const& local_x,
134 Eigen::VectorXd
const& local_xdot,
135 Eigen::VectorXd& local_b,
136 Eigen::MatrixXd& local_J)
138 auto const p = local_x.segment(pressure_index, pressure_size);
139 auto const p_dot = local_xdot.segment(pressure_index, pressure_size);
140 auto const g = local_x.segment(displacement_index, displacement_size);
142 local_xdot.segment(displacement_index, displacement_size);
144 auto rhs_p = local_b.segment(pressure_index, pressure_size);
145 auto rhs_g = local_b.segment(displacement_index, displacement_size);
146 auto J_pp = local_J.block(pressure_index, pressure_index, pressure_size,
148 auto J_pg = local_J.block(pressure_index, displacement_index, pressure_size,
150 auto J_gp = local_J.block(displacement_index, pressure_index,
151 displacement_size, pressure_size);
152 auto J_gg = local_J.block(displacement_index, displacement_index,
153 displacement_size, displacement_size);
155 assembleBlockMatricesWithJacobian(t, dt, p, p_dot, g, g_dot, rhs_p, rhs_g,
156 J_pp, J_pg, J_gg, J_gp);
159 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
160 typename IntegrationMethod,
int GlobalDim>
162 ShapeFunctionPressure,
163 IntegrationMethod, 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_dot,
168 Eigen::Ref<const Eigen::VectorXd>
const& g,
169 Eigen::Ref<const Eigen::VectorXd>
const& g_dot,
170 Eigen::Ref<Eigen::VectorXd> rhs_p, Eigen::Ref<Eigen::VectorXd> rhs_g,
171 Eigen::Ref<Eigen::MatrixXd> J_pp, Eigen::Ref<Eigen::MatrixXd> J_pg,
172 Eigen::Ref<Eigen::MatrixXd> J_gg, Eigen::Ref<Eigen::MatrixXd> J_gp)
174 auto const& frac_prop = *_process_data.fracture_property;
175 auto const& R = frac_prop.R;
179 auto constexpr index_normal = GlobalDim - 1;
182 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
186 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
189 typename ShapeMatricesTypeDisplacement::template MatrixType<
190 displacement_size, pressure_size>
191 Kgp = ShapeMatricesTypeDisplacement::template MatrixType<
192 displacement_size, pressure_size>::Zero(displacement_size,
196 Eigen::MatrixXd
const global2local_rotation =
197 R.template topLeftCorner<ShapeFunctionPressure::DIM, GlobalDim>();
200 global2local_rotation * _process_data.specific_body_force;
205 unsigned const n_integration_points = _ip_data.size();
206 for (
unsigned ip = 0; ip < n_integration_points; ip++)
210 auto& ip_data = _ip_data[ip];
211 auto const& ip_w = ip_data.integration_weight;
212 auto const& N_p = ip_data.N_p;
213 auto const& dNdx_p = ip_data.dNdx_p;
214 auto const& H_g = ip_data.H_u;
215 auto const& identity2 =
218 auto& mat = ip_data.fracture_material;
220 auto const& effective_stress_prev = ip_data.sigma_eff_prev;
222 auto const& w_prev = ip_data.w_prev;
224 auto& state = *ip_data.material_state_variables;
225 auto& b_m = ip_data.aperture;
227 double const S = frac_prop.specific_storage(t, x_position)[0];
228 double const mu = _process_data.fluid_viscosity(t, x_position)[0];
229 auto const alpha = frac_prop.biot_coefficient(t, x_position)[0];
230 auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
233 w.noalias() = R * H_g * g;
236 b_m = ip_data.aperture0 + w[index_normal];
240 "Element {:d}, gp {:d}: Fracture aperture is {:g}, but it must "
242 "non-negative. Setting it to zero.",
243 _element.getID(), ip, b_m);
247 auto const initial_effective_stress =
248 _process_data.initial_fracture_effective_stress(0, x_position);
250 Eigen::Map<typename HMatricesType::ForceVectorType const>
const stress0(
251 initial_effective_stress.data(), initial_effective_stress.size());
254 mat.computeConstitutiveRelation(
255 t, x_position, ip_data.aperture0, stress0, w_prev, w,
258 auto& k = ip_data.permeability;
259 k = frac_prop.permeability_model->permeability(
260 ip_data.permeability_state.get(), ip_data.aperture0, b_m);
263 double const local_dk_db =
264 frac_prop.permeability_model->dpermeability_daperture(
265 ip_data.permeability_state.get(), ip_data.aperture0, b_m);
272 J_gg.noalias() += H_g.transpose() * R.transpose() * C * R * H_g * ip_w;
278 H_g.transpose() * R.transpose() *
alpha * identity2 * N_p * ip_w;
283 storage_p.noalias() += N_p.transpose() * b_m * S * N_p * ip_w;
284 laplace_p.noalias() +=
285 dNdx_p.transpose() * b_m * k / mu * dNdx_p * ip_w;
287 dNdx_p.transpose() * b_m * k / mu * rho_fr * gravity_vec * ip_w;
292 auto const grad_head_over_mu =
293 ((dNdx_p * p + rho_fr * gravity_vec) / mu).eval();
294 Eigen::Matrix<double, 1, displacement_size>
const mT_R_Hg =
295 identity2.transpose() * R * H_g;
297 ip_data.darcy_velocity =
298 -global2local_rotation.transpose() * k * grad_head_over_mu;
299 J_pg.noalias() += N_p.transpose() * S * N_p * p_dot * mT_R_Hg * ip_w;
301 dNdx_p.transpose() * k * grad_head_over_mu * mT_R_Hg * ip_w;
302 J_pg.noalias() += dNdx_p.transpose() * b_m * local_dk_db *
303 grad_head_over_mu * mT_R_Hg * ip_w;
307 J_gp.noalias() -= Kgp;
310 J_pp.noalias() += laplace_p + storage_p / dt;
313 J_pg.noalias() += Kgp.transpose() / dt;
317 laplace_p * p + storage_p * p_dot + Kgp.transpose() * g_dot;
320 rhs_g.noalias() -= -Kgp * p;
323 template <
typename ShapeFunctionDisplacement,
typename ShapeFunctionPressure,
324 typename IntegrationMethod,
int GlobalDim>
326 ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
327 GlobalDim>::postTimestepConcreteWithVector(
const double t,
329 Eigen::VectorXd
const& local_x)
331 auto const nodal_g = local_x.segment(displacement_index, displacement_size);
333 auto const& frac_prop = *_process_data.fracture_property;
334 auto const& R = frac_prop.R;
337 auto constexpr index_normal = GlobalDim - 1;
342 unsigned const n_integration_points = _ip_data.size();
343 for (
unsigned ip = 0; ip < n_integration_points; ip++)
347 auto& ip_data = _ip_data[ip];
348 auto const& H_g = ip_data.H_u;
349 auto& mat = ip_data.fracture_material;
351 auto const& effective_stress_prev = ip_data.sigma_eff_prev;
353 auto const& w_prev = ip_data.w_prev;
355 auto& state = *ip_data.material_state_variables;
356 auto& b_m = ip_data.aperture;
359 w.noalias() = R * H_g * nodal_g;
362 b_m = ip_data.aperture0 + w[index_normal];
366 "Element {:d}, gp {:d}: Fracture aperture is {:g}, but it is "
367 "expected to be non-negative.",
368 _element.getID(), ip, b_m);
371 auto const initial_effective_stress =
372 _process_data.initial_fracture_effective_stress(0, x_position);
374 Eigen::Map<typename HMatricesType::ForceVectorType const>
const stress0(
375 initial_effective_stress.data(), initial_effective_stress.size());
378 mat.computeConstitutiveRelation(
379 t, x_position, ip_data.aperture0, stress0, w_prev, w,
386 HMatricesType::ForceVectorType::Zero(GlobalDim);
388 HMatricesType::ForceVectorType::Zero(GlobalDim);
389 double ele_Fs = -std::numeric_limits<double>::max();
391 for (
auto const& ip : _ip_data)
393 ele_b += ip.aperture;
394 ele_k += ip.permeability;
396 ele_sigma_eff += ip.sigma_eff;
398 ele_Fs, ip.material_state_variables->getShearYieldFunctionValue());
399 ele_velocity += ip.darcy_velocity;
401 ele_b /=
static_cast<double>(n_integration_points);
402 ele_k /=
static_cast<double>(n_integration_points);
403 ele_w /=
static_cast<double>(n_integration_points);
404 ele_sigma_eff /=
static_cast<double>(n_integration_points);
405 ele_velocity /=
static_cast<double>(n_integration_points);
406 auto const element_id = _element.getID();
407 (*_process_data.mesh_prop_b)[element_id] = ele_b;
408 (*_process_data.mesh_prop_k_f)[element_id] = ele_k;
409 (*_process_data.mesh_prop_w_n)[element_id] = ele_w[index_normal];
410 (*_process_data.mesh_prop_w_s)[element_id] = ele_w[0];
411 (*_process_data.mesh_prop_fracture_stress_normal)[element_id] =
412 ele_sigma_eff[index_normal];
413 (*_process_data.mesh_prop_fracture_stress_shear)[element_id] =
415 (*_process_data.mesh_prop_fracture_shear_failure)[element_id] = ele_Fs;
419 (*_process_data.mesh_prop_w_s2)[element_id] = ele_w[1];
420 (*_process_data.mesh_prop_fracture_stress_shear2)[element_id] =
424 for (
unsigned i = 0; i < GlobalDim; i++)
426 (*_process_data.mesh_prop_velocity)[element_id * GlobalDim + i] =
void DBUG(char const *fmt, Args const &... args)
virtual constexpr unsigned getDimension() const =0
Get dimension of the mesh element.
virtual std::size_t getID() const final
Returns the ID of the element.
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
VectorType< DisplacementDim > ForceVectorType
HydroMechanicsProcessData< GlobalDim > & _process_data
Eigen::Matrix< double, GlobalDim, 1 > GlobalDimVector
ShapeMatrixPolicyType< ShapeFunctionPressure, ShapeFunctionPressure::DIM > ShapeMatricesTypePressure
std::vector< IntegrationPointDataType, Eigen::aligned_allocator< IntegrationPointDataType > > _ip_data
typename HMatricesType::HMatrixType HMatrixType
typename ShapeMatricesTypePressure::DimVectorType DimVectorType
ShapeMatrixPolicyType< ShapeFunctionDisplacement, ShapeFunctionDisplacement::DIM > ShapeMatricesTypeDisplacement
Eigen::Matrix< double, 3u, 3u, Eigen::RowMajor > RotationMatrix
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
Eigen::Matrix< double, GlobalDim, GlobalDim > createRotatedTensor(RotationMatrix const &R, double const value)
void computeHMatrix(N_Type const &N, HMatrixType &H)
Fills a H-matrix based on given shape function.
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
VectorType< ShapeFunction::NPOINTS > NodalVectorType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType