14#include <Eigen/Eigenvalues>
30template <
typename ShapeFunction,
int DisplacementDim>
33 Eigen::VectorXd
const& local_x,
34 Eigen::VectorXd
const& local_x_prev,
36 std::vector<double>& local_b_data,
37 std::vector<double>& local_Jac_data)
51 local_b_data, local_Jac_data);
60template <
typename ShapeFunction,
int DisplacementDim>
63 double const t,
double const dt, Eigen::VectorXd
const& local_x,
64 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
69 auto const p = local_x.template segment<pressure_size>(
pressure_index);
81 auto const& solid = medium->phase(
"Solid");
82 auto const& fluid = fluidPhase(*medium);
85 double const k =
_process_data.residual_stiffness(t, x_position)[0];
86 double const ls =
_process_data.crack_length_scale(t, x_position)[0];
91 for (
int ip = 0; ip < n_integration_points; ip++)
93 auto const& w =
_ip_data[ip].integration_weight;
95 auto const& dNdx =
_ip_data[ip].dNdx;
96 double const d_ip = N.dot(d);
97 double const p_ip = N.dot(p);
105 ShapeFunction::NPOINTS,
110 eps.noalias() = B * u;
112 double const degradation =
113 _process_data.degradation_derivative->degradation(d_ip, k, ls);
114 _ip_data[ip].updateConstitutiveRelation(
115 t, x_position, dt, u, degradation,
118 auto const& sigma =
_ip_data[ip].sigma;
121 auto& biot_coefficient =
_ip_data[ip].biot_coefficient;
122 auto& biot_modulus_inv =
_ip_data[ip].biot_modulus_inv;
123 auto const& fracture_enhanced_porosity =
124 _ip_data[ip].fracture_enhanced_porosity;
128 auto const D_sph = P_sph * D * identity2;
131 auto const& solid_material =
135 auto const bulk_modulus = solid_material.getBulkModulus(t, x_position);
138 .template value<double>(vars, x_position, t, dt);
139 auto const porosity_0 =
141 .template value<double>(vars, x_position, t, dt);
142 double const bulk_modulus_degradation = bulk_modulus_eff / bulk_modulus;
145 biot_coefficient = 1. - bulk_modulus_degradation * (1. - alpha_0);
148 auto const porosity_reference = porosity_0 + fracture_enhanced_porosity;
151 biot_modulus_inv = bulk_modulus_degradation * (alpha_0 - porosity_0) *
152 (1. - alpha_0) / bulk_modulus;
156 .template value<double>(vars, x_position, t, dt);
159 .template value<double>(vars, x_position, t, dt);
162 rho_sr * (1. - porosity_reference) + porosity_reference * rho_fr;
165 local_rhs.noalias() -=
166 (B.transpose() * (sigma - biot_coefficient * p_ip * identity2) -
167 N_u_op(N).transpose() * rho * b) *
169 local_Jac.noalias() += B.transpose() * D * B * w;
173template <
typename ShapeFunction,
int DisplacementDim>
176 Eigen::VectorXd
const& local_x,
177 Eigen::VectorXd
const& local_x_prev,
178 std::vector<double>& local_b_data,
179 std::vector<double>& local_Jac_data)
183 auto const p = local_x.template segment<pressure_size>(
pressure_index);
206 auto const& fluid = fluidPhase(*medium);
214 double const k =
_process_data.residual_stiffness(t, x_position)[0];
215 double const ls =
_process_data.crack_length_scale(t, x_position)[0];
217 double const fracture_threshold =
_process_data.fracture_threshold;
218 double const fracture_permeability_parameter =
220 double const fixed_stress_stabilization_parameter =
222 double const spatial_stabilization_parameter =
228 for (
int ip = 0; ip < n_integration_points; ip++)
230 auto const& w =
_ip_data[ip].integration_weight;
232 auto const& dNdx =
_ip_data[ip].dNdx;
233 double const d_ip = N.dot(d);
237 .template value<double>(vars, x_position, t, dt);
240 .template value<double>(vars, x_position, t, dt);
243 .template value<double>(vars, x_position, t, dt);
247 auto const& biot_coefficient =
_ip_data[ip].biot_coefficient;
248 auto const& biot_modulus_inv =
_ip_data[ip].biot_modulus_inv;
249 auto const& fracture_enhanced_porosity =
250 _ip_data[ip].fracture_enhanced_porosity;
253 auto const porosity_0 =
255 .template value<double>(vars, x_position, t, dt);
256 auto const porosity_reference = porosity_0 + fracture_enhanced_porosity;
258 double const dv_dt = (vol_strain - vol_strain_prev) / dt;
262 auto const D_sph = P_sph * D * identity2;
263 auto const D_dev = P_dev * D * (ones2 - identity2) / std::sqrt(2.);
264 auto const degraded_shear_modulus =
268 double const residual_bulk_modulus = [&]
274 double const degradation_threshold =
276 fracture_threshold, k, ls);
277 auto const& D_threshold =
278 degradation_threshold *
_ip_data[ip].C_tensile +
280 auto const D_sph_threshold = P_sph * D_threshold * identity2;
284 return degraded_bulk_modulus;
287 double const modulus_rm = fixed_stress_stabilization_parameter *
288 biot_coefficient * biot_coefficient /
289 residual_bulk_modulus;
291 double const stablization_spatial =
292 spatial_stabilization_parameter * 0.25 * he * he /
293 (degraded_bulk_modulus + 4. / 3. * degraded_shear_modulus);
294 stablizing.noalias() +=
295 dNdx.transpose() * stablization_spatial * dNdx * w;
298 (biot_modulus_inv + cf * porosity_reference + modulus_rm) *
299 N.transpose() * N * w;
301 auto const K_over_mu = K / mu;
302 laplace.noalias() += dNdx.transpose() * K_over_mu * dNdx * w;
306 local_rhs.noalias() += dNdx.transpose() * rho_fr * K_over_mu * b * w;
308 local_rhs.noalias() -= (biot_coefficient * dv_dt -
309 modulus_rm *
_ip_data[ip].coupling_pressure) *
316 auto const& normal_ip =
_ip_data[ip].normal_ip;
318 std::pow(1. - d_ip, fracture_permeability_parameter) * width *
319 width * width / he / 12.0 *
320 (Eigen::Matrix<double, DisplacementDim,
321 DisplacementDim>::Identity() -
322 normal_ip * normal_ip.transpose());
323 laplace.noalias() += dNdx.transpose() * Kf / mu * dNdx * w;
326 local_Jac.noalias() = laplace + mass / dt + stablizing / dt;
328 local_rhs.noalias() -=
329 laplace * p + mass * (p - p_prev) / dt + stablizing * (p - p_prev) / dt;
332template <
typename ShapeFunction,
int DisplacementDim>
335 Eigen::VectorXd
const& local_x,
336 std::vector<double>& local_b_data,
337 std::vector<double>& local_Jac_data)
340 auto const p = local_x.template segment<pressure_size>(
pressure_index);
352 auto const& solid_material =
357 auto const bulk_modulus = solid_material.getBulkModulus(t, x_position);
359 auto const& solid = medium->phase(
"Solid");
362 double const k =
_process_data.residual_stiffness(t, x_position)[0];
363 double const ls =
_process_data.crack_length_scale(t, x_position)[0];
364 double const gc =
_process_data.crack_resistance(t, x_position)[0];
369 for (
int ip = 0; ip < n_integration_points; ip++)
371 auto const& w =
_ip_data[ip].integration_weight;
373 auto const& dNdx =
_ip_data[ip].dNdx;
375 double const d_ip = N.dot(d);
376 double const p_ip = N.dot(p);
377 double const degradation =
378 _process_data.degradation_derivative->degradation(d_ip, k, ls);
379 double const degradation_df1 =
380 _process_data.degradation_derivative->degradationDf1(d_ip, k, ls);
381 double const degradation_df2 =
382 _process_data.degradation_derivative->degradationDf2(d_ip, k, ls);
384 _ip_data[ip].updateConstitutiveRelation(
385 t, x_position, dt, u, degradation,
388 auto& biot_coefficient =
_ip_data[ip].biot_coefficient;
389 auto& biot_modulus_inv =
_ip_data[ip].biot_modulus_inv;
393 .template value<double>(vars, x_position, t, dt);
394 auto const porosity_0 =
396 .template value<double>(vars, x_position, t, dt);
400 auto const D_sph = P_sph * D * identity2;
402 double const bulk_modulus_degradation = bulk_modulus_eff / bulk_modulus;
405 biot_coefficient = 1. - bulk_modulus_degradation * (1. - alpha_0);
408 biot_modulus_inv = bulk_modulus_degradation * (alpha_0 - porosity_0) *
409 (1. - alpha_0) / bulk_modulus;
411 auto const& strain_energy_tensile =
_ip_data[ip].strain_energy_tensile;
412 auto const& C_tensile =
_ip_data[ip].C_tensile;
413 auto const C_tensile_sph = P_sph * C_tensile * identity2;
416 auto const driven_energy =
418 (strain_energy_tensile + p_ip * p_ip / 2. * bulk_modulus_plus /
419 bulk_modulus * (alpha_0 - porosity_0) *
420 (1. - alpha_0) / bulk_modulus) *
423 local_Jac.noalias() += driven_energy * N * degradation_df2;
425 local_rhs.noalias() -= driven_energy * degradation_df1;
427 calculateCrackLocalJacobianAndResidual<
428 decltype(dNdx),
decltype(N),
decltype(w),
decltype(d),
429 decltype(local_Jac),
decltype(local_rhs)>(
430 dNdx, N, w, d, local_Jac, local_rhs, gc, ls,
435template <
typename ShapeFunction,
int DisplacementDim>
438 Eigen::VectorXd
const& local_x_prev,
439 double const ,
double const dt,
443 auto const p = local_x.template segment<pressure_size>(
pressure_index);
447 for (
int ip = 0; ip < n_integration_points; ip++)
450 _ip_data[ip].coupling_pressure = N.dot(p - p_prev) / dt;
454template <
typename ShapeFunction,
int DisplacementDim>
457 std::size_t mesh_item_id,
458 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
459 std::vector<GlobalVector*>
const& x,
double const t,
462 std::vector<std::vector<GlobalIndexType>> indices_of_processes;
463 indices_of_processes.reserve(dof_tables.size());
464 std::transform(dof_tables.begin(), dof_tables.end(),
465 std::back_inserter(indices_of_processes),
466 [&](
auto const dof_table)
467 { return NumLib::getIndices(mesh_item_id, *dof_table); });
470 assert(local_coupled_xs.size() ==
473 auto const d = Eigen::Map<PhaseFieldVector const>(
479 double const ele_d = std::clamp(d.sum() / d.size(), 0.0, 1.0);
485 double const width_init =
_process_data.width_init(t, x_position)[0];
486 double const k =
_process_data.residual_stiffness(t, x_position)[0];
487 double const ls =
_process_data.crack_length_scale(t, x_position)[0];
488 double const he = ls /
_process_data.diffused_range_parameter;
490 int const n_integration_points =
493 for (
int ip = 0; ip < n_integration_points; ip++)
497 Eigen::EigenSolver<
decltype(eps_tensor)> eigen_solver(eps_tensor);
498 Eigen::MatrixXf::Index maxIndex;
499 double const max_principal_strain =
500 eigen_solver.eigenvalues().real().maxCoeff(&maxIndex);
501 auto const max_eigen_vector =
502 eigen_solver.eigenvectors().real().col(maxIndex);
505 auto& width_ip =
_ip_data[ip].width_ip;
506 width_ip = max_principal_strain * he;
507 width_ip = width_ip < width_init ? width_init : width_ip;
511 auto& normal_ip =
_ip_data[ip].normal_ip;
512 if (std::abs(max_principal_strain) > k)
514 for (
int i = 0; i < DisplacementDim; i++)
516 normal_ip[i] = max_eigen_vector[i];
521 auto& fracture_enhanced_porosity =
522 _ip_data[ip].fracture_enhanced_porosity;
523 fracture_enhanced_porosity = width_ip / he;
531template <
typename ShapeFunction,
int DisplacementDim>
533 std::size_t mesh_item_id,
534 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
535 std::vector<GlobalVector*>
const& x,
double const t,
double& elastic_energy,
536 double& surface_energy,
double& pressure_work)
538 std::vector<std::vector<GlobalIndexType>> indices_of_processes;
539 indices_of_processes.reserve(dof_tables.size());
540 std::transform(dof_tables.begin(), dof_tables.end(),
541 std::back_inserter(indices_of_processes),
542 [&](
auto const dof_table)
543 { return NumLib::getIndices(mesh_item_id, *dof_table); });
545 auto const local_coupled_xs =
547 assert(local_coupled_xs.size() ==
550 auto const d = Eigen::Map<PhaseFieldVector const>(
552 auto const u = Eigen::Map<DeformationVector const>(
554 auto const p = Eigen::Map<PressureVector const>(
560 double element_elastic_energy = 0.0;
561 double element_surface_energy = 0.0;
562 double element_pressure_work = 0.0;
564 double const gc =
_process_data.crack_resistance(t, x_position)[0];
565 double const ls =
_process_data.crack_length_scale(t, x_position)[0];
567 for (
int ip = 0; ip < n_integration_points; ip++)
569 auto const& w =
_ip_data[ip].integration_weight;
571 auto const& dNdx =
_ip_data[ip].dNdx;
572 double const d_ip = N.dot(d);
573 double const p_ip = N.dot(p);
575 element_elastic_energy +=
_ip_data[ip].elastic_energy * w;
581 element_surface_energy +=
583 ((1 - d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) * w;
589 element_surface_energy += 0.5 * gc *
590 ((1 - d_ip) * (1 - d_ip) / ls +
591 (dNdx * d).dot((dNdx * d)) * ls) *
597 element_surface_energy +=
598 gc / std::numbers::pi *
599 ((1 - d_ip * d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) *
605 element_pressure_work += p_ip * (
N_u_op(N) * u).dot(dNdx * d) * w;
609 int const n_all_nodes = indices_of_processes[1].size();
610 int const n_regular_nodes = std::count_if(
611 begin(indices_of_processes[1]), end(indices_of_processes[1]),
613 if (n_all_nodes != n_regular_nodes)
615 element_elastic_energy *=
616 static_cast<double>(n_regular_nodes) / n_all_nodes;
617 element_surface_energy *=
618 static_cast<double>(n_regular_nodes) / n_all_nodes;
619 element_pressure_work *=
620 static_cast<double>(n_regular_nodes) / n_all_nodes;
623 elastic_energy += element_elastic_energy;
624 surface_energy += element_surface_energy;
625 pressure_work += element_pressure_work;
628template <
typename ShapeFunction,
int DisplacementDim>
629std::vector<double>
const&
632 std::vector<GlobalVector*>
const& ,
633 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
634 std::vector<double>& cache)
const
GlobalMatrix::IndexType GlobalIndexType
void setElementID(std::size_t element_id)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
void approximateFractureWidth(std::size_t mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double const dt) override
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
void computeEnergy(std::size_t mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double &elastic_energy, double &surface_energy, double &pressure_work) override
MeshLib::Element const & _element
static constexpr int displacement_size
static constexpr int phasefield_index
static constexpr auto & N_u_op
void assembleWithJacobianForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
HMPhaseFieldProcessData< DisplacementDim > & _process_data
static constexpr int pressure_size
void postNonLinearSolverConcrete(Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, int const process_id) override
std::vector< double > const & getIntPtWidth(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
void assembleWithJacobianPhaseFieldEquations(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void assembleWithJacobianHydroEquations(const double t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void assembleWithJacobianForDeformationEquations(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
bool const _is_axially_symmetric
NumLib::GenericIntegrationMethod const & _integration_method
static constexpr int displacement_index
static constexpr int pressure_index
static constexpr int phasefield_size
auto & selectSolidConstitutiveRelation(SolidMaterialsMap const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Eigen::Matrix< double, 3, 3 > kelvinVectorToTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
std::vector< double > const & getIntegrationPointScalarData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
std::vector< double > getCoupledLocalSolutions(std::vector< GlobalVector * > const &global_solutions, std::vector< std::vector< GlobalIndexType > > const &indices)
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
static Eigen::Matrix< double, KelvinVectorSize, KelvinVectorSize > const spherical_projection
static double FrobeniusNorm(Eigen::Matrix< double, KelvinVectorSize, 1 > const &deviatoric_v)
Get the norm of the deviatoric stress.
static Eigen::Matrix< double, KelvinVectorSize, 1 > const identity2
Kelvin mapping of 2nd order identity tensor.
static Eigen::Matrix< double, KelvinVectorSize, KelvinVectorSize > const deviatoric_projection
static Eigen::Matrix< double, KelvinVectorSize, 1 > const ones2
static double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.