Loading [MathJax]/extensions/tex2jax.js
OGS
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim > Class Template Reference

Detailed Description

template<typename ShapeFunction, int DisplacementDim>
class ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >

Definition at line 113 of file HMPhaseFieldFEM.h.

#include <HMPhaseFieldFEM.h>

Inheritance diagram for ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >:
[legend]
Collaboration diagram for ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >:
[legend]

Public Types

using ShapeMatricesType
using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices
using BMatricesType = BMatrixPolicyType<ShapeFunction, DisplacementDim>
using NodalForceVectorType = typename BMatricesType::NodalForceVectorType
using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType
using DeformationVector
using DeformationMatrix
using PhaseFieldVector
using PhaseFieldMatrix
using PressureVector
using PressureMatrix
using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType
using NodalVectorType = typename ShapeMatricesType::NodalVectorType
using IpData
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>

Public Member Functions

 HMPhaseFieldLocalAssembler (HMPhaseFieldLocalAssembler const &)=delete
 HMPhaseFieldLocalAssembler (HMPhaseFieldLocalAssembler &&)=delete
 HMPhaseFieldLocalAssembler (MeshLib::Element const &e, std::size_t const, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, HMPhaseFieldProcessData< DisplacementDim > &process_data)
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
void initializeConcrete () override
void postTimestepConcrete (Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, int const) override
void postNonLinearSolverConcrete (Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, int const process_id) override
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
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
double heaviside (double const v)
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix (const unsigned integration_point) const override
 Provides the shape matrix at the given integration point.
std::vector< double > const & getIntPtWidth (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
virtual ~LocalAssemblerInterface ()=default
virtual void setInitialConditions (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, int const process_id)
virtual void initialize (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
virtual void preAssemble (double const, double const, std::vector< double > const &)
virtual void assemble (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
virtual void assembleForStaggeredScheme (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_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
virtual void assembleWithJacobian (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
virtual void computeSecondaryVariable (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id)
virtual void preTimestep (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
virtual void postTimestep (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
void postNonLinearSolver (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< double > const &) const
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< std::vector< double > > const &) const
 Fits to staggered scheme.
virtual std::optional< VectorSegment > getVectorDeformationSegment () const
Public Member Functions inherited from NumLib::ExtrapolatableElement
virtual ~ExtrapolatableElement ()=default

Static Public Attributes

static int const KelvinVectorSize
static constexpr auto & N_u_op

Private Member Functions

std::vector< double > const & getIntPtSigma (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
std::vector< double > const & getIntPtEpsilon (const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
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 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 assembleWithJacobianForDeformationEquations (double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)

Private Attributes

HMPhaseFieldProcessData< DisplacementDim > & _process_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
NumLib::GenericIntegrationMethod const & _integration_method
MeshLib::Element const & _element
SecondaryData< typename ShapeMatrices::ShapeType > _secondary_data
bool const _is_axially_symmetric

Static Private Attributes

static constexpr int phasefield_index = 0
static constexpr int phasefield_size = ShapeFunction::NPOINTS
static constexpr int pressure_index = phasefield_index + phasefield_size
static constexpr int pressure_size = ShapeFunction::NPOINTS
static constexpr int displacement_index = pressure_index + pressure_size
static constexpr int displacement_size

Member Typedef Documentation

◆ BMatricesType

template<typename ShapeFunction, int DisplacementDim>
using ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::BMatricesType = BMatrixPolicyType<ShapeFunction, DisplacementDim>

Definition at line 130 of file HMPhaseFieldFEM.h.

◆ DeformationMatrix

template<typename ShapeFunction, int DisplacementDim>
using ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::DeformationMatrix
Initial value:
typename ShapeMatricesType::template MatrixType<displacement_size,
displacement_size>

Definition at line 138 of file HMPhaseFieldFEM.h.

◆ DeformationVector

template<typename ShapeFunction, int DisplacementDim>
using ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::DeformationVector
Initial value:
typename ShapeMatricesType::template VectorType<displacement_size>

Definition at line 136 of file HMPhaseFieldFEM.h.

◆ GlobalDimVectorType

template<typename ShapeFunction, int DisplacementDim>
using ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType

Definition at line 134 of file HMPhaseFieldFEM.h.

◆ Invariants

template<typename ShapeFunction, int DisplacementDim>
using ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>

Definition at line 163 of file HMPhaseFieldFEM.h.

◆ IpData

template<typename ShapeFunction, int DisplacementDim>
using ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::IpData
Initial value:
IntegrationPointData<BMatricesType, ShapeMatricesType, DisplacementDim>

Definition at line 154 of file HMPhaseFieldFEM.h.

◆ NodalForceVectorType

template<typename ShapeFunction, int DisplacementDim>
using ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::NodalForceVectorType = typename BMatricesType::NodalForceVectorType

Definition at line 132 of file HMPhaseFieldFEM.h.

◆ NodalMatrixType

template<typename ShapeFunction, int DisplacementDim>
using ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::NodalMatrixType = typename ShapeMatricesType::NodalMatrixType

Definition at line 152 of file HMPhaseFieldFEM.h.

◆ NodalVectorType

template<typename ShapeFunction, int DisplacementDim>
using ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::NodalVectorType = typename ShapeMatricesType::NodalVectorType

Definition at line 153 of file HMPhaseFieldFEM.h.

◆ PhaseFieldMatrix

template<typename ShapeFunction, int DisplacementDim>
using ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::PhaseFieldMatrix
Initial value:
typename ShapeMatricesType::template MatrixType<phasefield_size,
phasefield_size>

Definition at line 143 of file HMPhaseFieldFEM.h.

◆ PhaseFieldVector

template<typename ShapeFunction, int DisplacementDim>
using ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::PhaseFieldVector
Initial value:
typename ShapeMatricesType::template VectorType<phasefield_size>

Definition at line 141 of file HMPhaseFieldFEM.h.

◆ PressureMatrix

template<typename ShapeFunction, int DisplacementDim>
using ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::PressureMatrix
Initial value:
typename ShapeMatricesType::template MatrixType<pressure_size,
pressure_size>

Definition at line 148 of file HMPhaseFieldFEM.h.

◆ PressureVector

template<typename ShapeFunction, int DisplacementDim>
using ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::PressureVector
Initial value:
typename ShapeMatricesType::template VectorType<pressure_size>

Definition at line 146 of file HMPhaseFieldFEM.h.

◆ ShapeMatrices

template<typename ShapeFunction, int DisplacementDim>
using ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::ShapeMatrices = typename ShapeMatricesType::ShapeMatrices

Definition at line 128 of file HMPhaseFieldFEM.h.

◆ ShapeMatricesType

template<typename ShapeFunction, int DisplacementDim>
using ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::ShapeMatricesType
Initial value:
ShapeMatrixPolicyType<ShapeFunction, DisplacementDim>
ShapeMatrixPolicyType
EigenFixedShapeMatrixPolicy< ShapeFunction, GlobalDim > ShapeMatrixPolicyType
Definition ShapeMatrixPolicy.h:141

Definition at line 125 of file HMPhaseFieldFEM.h.

Constructor & Destructor Documentation

◆ HMPhaseFieldLocalAssembler() [1/3]

template<typename ShapeFunction, int DisplacementDim>
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::HMPhaseFieldLocalAssembler ( HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim > const & )
delete

References HMPhaseFieldLocalAssembler().

Referenced by HMPhaseFieldLocalAssembler(), and HMPhaseFieldLocalAssembler().

◆ HMPhaseFieldLocalAssembler() [2/3]

template<typename ShapeFunction, int DisplacementDim>
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::HMPhaseFieldLocalAssembler ( HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim > && )
delete

References HMPhaseFieldLocalAssembler().

◆ HMPhaseFieldLocalAssembler() [3/3]

template<typename ShapeFunction, int DisplacementDim>
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::HMPhaseFieldLocalAssembler ( MeshLib::Element const & e,
std::size_t const ,
NumLib::GenericIntegrationMethod const & integration_method,
bool const is_axially_symmetric,
HMPhaseFieldProcessData< DisplacementDim > & process_data )
inline

Definition at line 168 of file HMPhaseFieldFEM.h.

174 : _process_data(process_data),
175 _integration_method(integration_method),
176 _element(e),
177 _is_axially_symmetric(is_axially_symmetric)
178 {
179 unsigned const n_integration_points =
180 _integration_method.getNumberOfPoints();
181
182 _ip_data.reserve(n_integration_points);
183 _secondary_data.N.resize(n_integration_points);
184
185 auto& solid_material =
186 MaterialLib::Solids::selectSolidConstitutiveRelation(
187 _process_data.solid_materials,
188 _process_data.material_ids,
189 e.getID());
190
191 auto const shape_matrices =
192 NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
193 DisplacementDim>(e, is_axially_symmetric,
194 _integration_method);
195
196 ParameterLib::SpatialPosition x_position;
197 x_position.setElementID(_element.getID());
198
199 for (unsigned ip = 0; ip < n_integration_points; ip++)
200 {
201 _ip_data.emplace_back(solid_material);
202 auto& ip_data = _ip_data[ip];
203 ip_data.integration_weight =
204 _integration_method.getWeightedPoint(ip).getWeight() *
205 shape_matrices[ip].integralMeasure * shape_matrices[ip].detJ;
206
207 static const int kelvin_vector_size =
208 MathLib::KelvinVector::kelvin_vector_dimensions(
209 DisplacementDim);
210 ip_data.eps_tensile.setZero(kelvin_vector_size);
211 ip_data.eps.setZero(kelvin_vector_size);
212 ip_data.eps_prev.resize(kelvin_vector_size);
213 ip_data.D.setZero(kelvin_vector_size, kelvin_vector_size);
214 ip_data.C_tensile.setZero(kelvin_vector_size, kelvin_vector_size);
215 ip_data.C_compressive.setZero(kelvin_vector_size,
216 kelvin_vector_size);
217
218 ip_data.sigma_tensile.setZero(kelvin_vector_size);
219 ip_data.sigma_compressive.setZero(kelvin_vector_size);
220 ip_data.sigma.setZero(kelvin_vector_size);
221 ip_data.strain_energy_tensile = 0.0;
222 ip_data.elastic_energy = 0.0;
223 ip_data.width_ip = 0.0;
224
225 ip_data.N = shape_matrices[ip].N;
226 ip_data.dNdx = shape_matrices[ip].dNdx;
227
228 _secondary_data.N[ip] = shape_matrices[ip].N;
229 }
230 }
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler::_ip_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
Definition HMPhaseFieldFEM.h:347
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler::_process_data
HMPhaseFieldProcessData< DisplacementDim > & _process_data
Definition HMPhaseFieldFEM.h:345
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler::ShapeMatricesType
ShapeMatrixPolicyType< ShapeFunction, DisplacementDim > ShapeMatricesType
Definition HMPhaseFieldFEM.h:125
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler::_integration_method
NumLib::GenericIntegrationMethod const & _integration_method
Definition HMPhaseFieldFEM.h:349
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler::_secondary_data
SecondaryData< typename ShapeMatrices::ShapeType > _secondary_data
Definition HMPhaseFieldFEM.h:351
MaterialLib::Solids::selectSolidConstitutiveRelation
auto & selectSolidConstitutiveRelation(SolidMaterialsMap const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
Definition SelectSolidConstitutiveRelation.h:36
MathLib::KelvinVector::kelvin_vector_dimensions
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition KelvinVector.h:24

References _element, _integration_method, _ip_data, _is_axially_symmetric, _process_data, _secondary_data, MeshLib::Element::getID(), NumLib::initShapeMatrices(), MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialLib::Solids::selectSolidConstitutiveRelation(), and ParameterLib::SpatialPosition::setElementID().

Member Function Documentation

◆ approximateFractureWidth()

template<typename ShapeFunction, int DisplacementDim>
void ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::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 )
overridevirtual

Implements ProcessLib::HMPhaseField::HMPhaseFieldLocalAssemblerInterface.

Definition at line 455 of file HMPhaseFieldFEM-impl.h.

461{
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); });
468
469 auto local_coupled_xs = getCoupledLocalSolutions(x, indices_of_processes);
470 assert(local_coupled_xs.size() ==
471 phasefield_size + displacement_size + pressure_size);
472
473 auto const d = Eigen::Map<PhaseFieldVector const>(
474 &local_coupled_xs[phasefield_index], phasefield_size);
475
476 ParameterLib::SpatialPosition x_position;
477 x_position.setElementID(_element.getID());
478
479 double const ele_d = std::clamp(d.sum() / d.size(), 0.0, 1.0);
480 (*_process_data.ele_d)[_element.getID()] = ele_d;
481
482 if ((*_process_data.ele_d)[_element.getID()] <
483 _process_data.irreversible_threshold)
484 {
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;
489
490 int const n_integration_points =
491 _integration_method.getNumberOfPoints();
492 double width = 0.0;
493 for (int ip = 0; ip < n_integration_points; ip++)
494 {
495 auto eps_tensor =
496 MathLib::KelvinVector::kelvinVectorToTensor(_ip_data[ip].eps);
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);
503
504 // Fracture aperture estimation
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;
508 width += width_ip;
509
510 // Fracture direction estimation
511 auto& normal_ip = _ip_data[ip].normal_ip;
512 if (std::abs(max_principal_strain) > k)
513 {
514 for (int i = 0; i < DisplacementDim; i++)
515 {
516 normal_ip[i] = max_eigen_vector[i];
517 }
518 }
519
520 // Fracture enhanced porosity
521 auto& fracture_enhanced_porosity =
522 _ip_data[ip].fracture_enhanced_porosity;
523 fracture_enhanced_porosity = width_ip / he;
524 }
525
526 // Update aperture for the fractured element
527 (*_process_data.width)[_element.getID()] = width / n_integration_points;
528 }
529}
MathLib::KelvinVector::kelvinVectorToTensor
Eigen::Matrix< double, 3, 3 > kelvinVectorToTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Definition KelvinVector.cpp:64
ProcessLib::getCoupledLocalSolutions
std::vector< double > getCoupledLocalSolutions(std::vector< GlobalVector * > const &global_solutions, std::vector< std::vector< GlobalIndexType > > const &indices)
Definition CoupledSolutionsForStaggeredScheme.cpp:21

References _element, _integration_method, _ip_data, _process_data, displacement_size, ProcessLib::getCoupledLocalSolutions(), MathLib::KelvinVector::kelvinVectorToTensor(), phasefield_index, phasefield_size, pressure_size, and ParameterLib::SpatialPosition::setElementID().

◆ assembleWithJacobianForDeformationEquations()

template<typename ShapeFunction, int DisplacementDim>
void ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::assembleWithJacobianForDeformationEquations ( double const t,
double const dt,
Eigen::VectorXd const & local_x,
std::vector< double > & local_b_data,
std::vector< double > & local_Jac_data )
private

Definition at line 61 of file HMPhaseFieldFEM-impl.h.

65{
66 auto const d = local_x.template segment<phasefield_size>(phasefield_index);
67 auto const u =
68 local_x.template segment<displacement_size>(displacement_index);
69 auto const p = local_x.template segment<pressure_size>(pressure_index);
70
71 auto local_Jac = MathLib::createZeroedMatrix<DeformationMatrix>(
72 local_Jac_data, displacement_size, displacement_size);
73
74 auto local_rhs = MathLib::createZeroedVector<DeformationVector>(
75 local_b_data, displacement_size);
76
77 ParameterLib::SpatialPosition x_position;
78 x_position.setElementID(_element.getID());
79
80 auto const& medium = _process_data.media_map.getMedium(_element.getID());
81 auto const& solid = medium->phase("Solid");
82 auto const& fluid = fluidPhase(*medium);
83 MPL::VariableArray vars;
84
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];
87
88 auto const& identity2 = Invariants::identity2;
89
90 int const n_integration_points = _integration_method.getNumberOfPoints();
91 for (int ip = 0; ip < n_integration_points; ip++)
92 {
93 auto const& w = _ip_data[ip].integration_weight;
94 auto const& N = _ip_data[ip].N;
95 auto const& dNdx = _ip_data[ip].dNdx;
96 double const d_ip = N.dot(d);
97 double const p_ip = N.dot(p);
98
99 auto const x_coord =
100 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
101 _element, N);
102
103 auto const& B =
104 LinearBMatrix::computeBMatrix<DisplacementDim,
105 ShapeFunction::NPOINTS,
106 typename BMatricesType::BMatrixType>(
107 dNdx, N, x_coord, _is_axially_symmetric);
108
109 auto& eps = _ip_data[ip].eps;
110 eps.noalias() = B * u;
111
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,
116 _process_data.energy_split_model);
117
118 auto const& sigma = _ip_data[ip].sigma;
119 auto const& D = _ip_data[ip].D;
120
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;
125
126 // Update the effective bulk modulus
127 auto const& P_sph = Invariants::spherical_projection;
128 auto const D_sph = P_sph * D * identity2;
129 double const bulk_modulus_eff = Invariants::trace(D_sph) / 9.;
130
131 auto const& solid_material =
132 MaterialLib::Solids::selectSolidConstitutiveRelation(
133 _process_data.solid_materials, _process_data.material_ids,
134 _element.getID());
135 auto const bulk_modulus = solid_material.getBulkModulus(t, x_position);
136 auto const alpha_0 =
137 solid.property(MPL::PropertyType::biot_coefficient)
138 .template value<double>(vars, x_position, t, dt);
139 auto const porosity_0 =
140 medium->property(MPL::PropertyType::porosity)
141 .template value<double>(vars, x_position, t, dt);
142 double const bulk_modulus_degradation = bulk_modulus_eff / bulk_modulus;
143
144 // Update Biot's coefficient
145 biot_coefficient = 1. - bulk_modulus_degradation * (1. - alpha_0);
146
147 // The reference porosity
148 auto const porosity_reference = porosity_0 + fracture_enhanced_porosity;
149
150 // Update Biot's modulus
151 biot_modulus_inv = bulk_modulus_degradation * (alpha_0 - porosity_0) *
152 (1. - alpha_0) / bulk_modulus;
153
154 auto const rho_sr =
155 solid.property(MPL::PropertyType::density)
156 .template value<double>(vars, x_position, t, dt);
157 auto const rho_fr =
158 fluid.property(MPL::PropertyType::density)
159 .template value<double>(vars, x_position, t, dt);
160
161 auto const rho =
162 rho_sr * (1. - porosity_reference) + porosity_reference * rho_fr;
163 auto const& b = _process_data.specific_body_force;
164
165 local_rhs.noalias() -=
166 (B.transpose() * (sigma - biot_coefficient * p_ip * identity2) -
167 N_u_op(N).transpose() * rho * b) *
168 w;
169 local_Jac.noalias() += B.transpose() * D * B * w;
170 }
171}
ProcessLib::BMatrixPolicyType::BMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition BMatrixPolicy.h:55
MathLib::createZeroedVector
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Definition EigenMapTools.h:153
MathLib::createZeroedMatrix
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition EigenMapTools.h:32
NumLib::interpolateXCoordinate
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
Definition InitShapeMatrices.h:91
MathLib::KelvinVector::Invariants::spherical_projection
static Eigen::Matrix< double, KelvinVectorSize, KelvinVectorSize > const spherical_projection
Definition KelvinVector.h:105
MathLib::KelvinVector::Invariants::identity2
static Eigen::Matrix< double, KelvinVectorSize, 1 > const identity2
Kelvin mapping of 2nd order identity tensor.
Definition KelvinVector.h:107
MathLib::KelvinVector::Invariants::trace
static double trace(Eigen::Matrix< double, KelvinVectorSize, 1 > const &v)
Trace of the corresponding tensor.
Definition KelvinVector-impl.h:62

References _element, _integration_method, _ip_data, _is_axially_symmetric, _process_data, MaterialPropertyLib::biot_coefficient, ProcessLib::LinearBMatrix::computeBMatrix(), MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, displacement_index, displacement_size, MathLib::KelvinVector::Invariants< KelvinVectorSize >::identity2, NumLib::interpolateXCoordinate(), N_u_op, phasefield_index, MaterialPropertyLib::porosity, pressure_index, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), MathLib::KelvinVector::Invariants< KelvinVectorSize >::spherical_projection, and MathLib::KelvinVector::Invariants< KelvinVectorSize >::trace().

Referenced by assembleWithJacobianForStaggeredScheme().

◆ assembleWithJacobianForStaggeredScheme()

template<typename ShapeFunction, int DisplacementDim>
void ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::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 )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 31 of file HMPhaseFieldFEM-impl.h.

38{
39 // For the equations with phase field.
40 if (process_id == _process_data._phasefield_process_id)
41 {
42 assembleWithJacobianPhaseFieldEquations(t, dt, local_x, local_b_data,
43 local_Jac_data);
44 return;
45 }
46
47 // For the equations for hydro
48 if (process_id == _process_data._hydro_process_id)
49 {
50 assembleWithJacobianHydroEquations(t, dt, local_x, local_x_prev,
51 local_b_data, local_Jac_data);
52 return;
53 }
54
55 // For the equations with deformation
56 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_b_data,
57 local_Jac_data);
58}
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler::assembleWithJacobianPhaseFieldEquations
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)
Definition HMPhaseFieldFEM-impl.h:334
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler::assembleWithJacobianHydroEquations
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)
Definition HMPhaseFieldFEM-impl.h:175
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler::assembleWithJacobianForDeformationEquations
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)
Definition HMPhaseFieldFEM-impl.h:62

References _process_data, assembleWithJacobianForDeformationEquations(), assembleWithJacobianHydroEquations(), and assembleWithJacobianPhaseFieldEquations().

◆ assembleWithJacobianHydroEquations()

template<typename ShapeFunction, int DisplacementDim>
void ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::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 )
private

Definition at line 174 of file HMPhaseFieldFEM-impl.h.

180{
181 auto const d = local_x.template segment<phasefield_size>(phasefield_index);
182
183 auto const p = local_x.template segment<pressure_size>(pressure_index);
184 auto const p_prev =
185 local_x_prev.template segment<pressure_size>(pressure_index);
186
187 auto local_Jac = MathLib::createZeroedMatrix<PressureMatrix>(
188 local_Jac_data, pressure_size, pressure_size);
189
190 auto local_rhs = MathLib::createZeroedVector<PressureVector>(local_b_data,
191 pressure_size);
192
193 typename ShapeMatricesType::NodalMatrixType mass =
194 ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
195
196 typename ShapeMatricesType::NodalMatrixType laplace =
197 ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
198
199 typename ShapeMatricesType::NodalMatrixType stablizing =
200 ShapeMatricesType::NodalMatrixType::Zero(pressure_size, pressure_size);
201
202 ParameterLib::SpatialPosition x_position;
203 x_position.setElementID(_element.getID());
204
205 auto const& medium = _process_data.media_map.getMedium(_element.getID());
206 auto const& fluid = fluidPhase(*medium);
207 MPL::VariableArray vars;
208
209 auto const& P_sph = Invariants::spherical_projection;
210 auto const& P_dev = Invariants::deviatoric_projection;
211 auto const& identity2 = Invariants::identity2;
212 auto const& ones2 = Invariants::ones2;
213
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];
216 double const width = (*_process_data.width)[_element.getID()];
217 double const fracture_threshold = _process_data.fracture_threshold;
218 double const fracture_permeability_parameter =
219 _process_data.fracture_permeability_parameter;
220 double const fixed_stress_stabilization_parameter =
221 _process_data.fixed_stress_stabilization_parameter;
222 double const spatial_stabilization_parameter =
223 _process_data.spatial_stabilization_parameter;
224 auto const he =
225 ls / _process_data.diffused_range_parameter; // element size
226
227 int const n_integration_points = _integration_method.getNumberOfPoints();
228 for (int ip = 0; ip < n_integration_points; ip++)
229 {
230 auto const& w = _ip_data[ip].integration_weight;
231 auto const& N = _ip_data[ip].N;
232 auto const& dNdx = _ip_data[ip].dNdx;
233 double const d_ip = N.dot(d);
234
235 auto const rho_fr =
236 fluid.property(MPL::PropertyType::density)
237 .template value<double>(vars, x_position, t, dt);
238 double const cf = _process_data.fluid_compressibility;
239 auto const Km = medium->property(MPL::PropertyType::permeability)
240 .template value<double>(vars, x_position, t, dt);
241 auto const K = MPL::formEigenTensor<DisplacementDim>(Km);
242 auto const mu = fluid.property(MPL::PropertyType::viscosity)
243 .template value<double>(vars, x_position, t, dt);
244
245 auto const vol_strain = Invariants::trace(_ip_data[ip].eps);
246 auto const vol_strain_prev = Invariants::trace(_ip_data[ip].eps_prev);
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;
251
252 // The reference porosity
253 auto const porosity_0 =
254 medium->property(MPL::PropertyType::porosity)
255 .template value<double>(vars, x_position, t, dt);
256 auto const porosity_reference = porosity_0 + fracture_enhanced_porosity;
257
258 double const dv_dt = (vol_strain - vol_strain_prev) / dt;
259
260 // The degraded shear modulus
261 auto const& D = _ip_data[ip].D;
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 =
265 Invariants::FrobeniusNorm(D_dev) / 2.;
266 auto const degraded_bulk_modulus = Invariants::trace(D_sph) / 9.;
267
268 double const residual_bulk_modulus = [&]
269 {
270 if ((*_process_data.ele_d)[_element.getID()] <
271 _process_data.fracture_threshold)
272 {
273 // The residual bulk modulus in the fractured element
274 double const degradation_threshold =
275 _process_data.degradation_derivative->degradation(
276 fracture_threshold, k, ls);
277 auto const& D_threshold =
278 degradation_threshold * _ip_data[ip].C_tensile +
279 _ip_data[ip].C_compressive;
280 auto const D_sph_threshold = P_sph * D_threshold * identity2;
281
282 return Invariants::trace(D_sph_threshold) / 9.;
283 }
284 return degraded_bulk_modulus;
285 }();
286
287 double const modulus_rm = fixed_stress_stabilization_parameter *
288 biot_coefficient * biot_coefficient /
289 residual_bulk_modulus;
290
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;
296
297 mass.noalias() +=
298 (biot_modulus_inv + cf * porosity_reference + modulus_rm) *
299 N.transpose() * N * w;
300
301 auto const K_over_mu = K / mu;
302 laplace.noalias() += dNdx.transpose() * K_over_mu * dNdx * w;
303 auto const& b = _process_data.specific_body_force;
304
305 // bodyforce-driven Darcy flow
306 local_rhs.noalias() += dNdx.transpose() * rho_fr * K_over_mu * b * w;
307
308 local_rhs.noalias() -= (biot_coefficient * dv_dt -
309 modulus_rm * _ip_data[ip].coupling_pressure) *
310 N.transpose() * w;
311
312 if ((*_process_data.ele_d)[_element.getID()] <
313 _process_data.irreversible_threshold)
314 {
315 // Fracture-enhanced permeability
316 auto const& normal_ip = _ip_data[ip].normal_ip;
317 auto const Kf =
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;
324 }
325 }
326 local_Jac.noalias() = laplace + mass / dt + stablizing / dt;
327
328 local_rhs.noalias() -=
329 laplace * p + mass * (p - p_prev) / dt + stablizing * (p - p_prev) / dt;
330}
MaterialPropertyLib::formEigenTensor
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Definition FormEigenTensor.cpp:126
EigenFixedShapeMatrixPolicy::NodalMatrixType
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
Definition ShapeMatrixPolicy.h:74
MathLib::KelvinVector::Invariants::FrobeniusNorm
static double FrobeniusNorm(Eigen::Matrix< double, KelvinVectorSize, 1 > const &deviatoric_v)
Get the norm of the deviatoric stress.
Definition KelvinVector-impl.h:27
MathLib::KelvinVector::Invariants::deviatoric_projection
static Eigen::Matrix< double, KelvinVectorSize, KelvinVectorSize > const deviatoric_projection
Definition KelvinVector.h:101
MathLib::KelvinVector::Invariants::ones2
static Eigen::Matrix< double, KelvinVectorSize, 1 > const ones2
Definition KelvinVector.h:111

References _element, _integration_method, _ip_data, _process_data, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, MathLib::KelvinVector::Invariants< KelvinVectorSize >::deviatoric_projection, MaterialPropertyLib::formEigenTensor(), MathLib::KelvinVector::Invariants< KelvinVectorSize >::FrobeniusNorm(), MathLib::KelvinVector::Invariants< KelvinVectorSize >::identity2, MathLib::KelvinVector::Invariants< KelvinVectorSize >::ones2, MaterialPropertyLib::permeability, phasefield_index, MaterialPropertyLib::porosity, pressure_index, pressure_size, ParameterLib::SpatialPosition::setElementID(), MathLib::KelvinVector::Invariants< KelvinVectorSize >::spherical_projection, MathLib::KelvinVector::Invariants< KelvinVectorSize >::trace(), and MaterialPropertyLib::viscosity.

Referenced by assembleWithJacobianForStaggeredScheme().

◆ assembleWithJacobianPhaseFieldEquations()

template<typename ShapeFunction, int DisplacementDim>
void ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::assembleWithJacobianPhaseFieldEquations ( double const t,
double const dt,
Eigen::VectorXd const & local_x,
std::vector< double > & local_b_data,
std::vector< double > & local_Jac_data )
private

Definition at line 333 of file HMPhaseFieldFEM-impl.h.

338{
339 auto const d = local_x.template segment<phasefield_size>(phasefield_index);
340 auto const p = local_x.template segment<pressure_size>(pressure_index);
341 auto const u =
342 local_x.template segment<displacement_size>(displacement_index);
343
344 auto local_Jac = MathLib::createZeroedMatrix<PhaseFieldMatrix>(
345 local_Jac_data, phasefield_size, phasefield_size);
346 auto local_rhs = MathLib::createZeroedVector<PhaseFieldVector>(
347 local_b_data, phasefield_size);
348
349 ParameterLib::SpatialPosition x_position;
350 x_position.setElementID(_element.getID());
351
352 auto const& solid_material =
353 MaterialLib::Solids::selectSolidConstitutiveRelation(
354 _process_data.solid_materials, _process_data.material_ids,
355 _element.getID());
356
357 auto const bulk_modulus = solid_material.getBulkModulus(t, x_position);
358 auto const& medium = _process_data.media_map.getMedium(_element.getID());
359 auto const& solid = medium->phase("Solid");
360 MPL::VariableArray vars;
361
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];
365
366 auto const& identity2 = Invariants::identity2;
367
368 int const n_integration_points = _integration_method.getNumberOfPoints();
369 for (int ip = 0; ip < n_integration_points; ip++)
370 {
371 auto const& w = _ip_data[ip].integration_weight;
372 auto const& N = _ip_data[ip].N;
373 auto const& dNdx = _ip_data[ip].dNdx;
374
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);
383
384 _ip_data[ip].updateConstitutiveRelation(
385 t, x_position, dt, u, degradation,
386 _process_data.energy_split_model);
387
388 auto& biot_coefficient = _ip_data[ip].biot_coefficient;
389 auto& biot_modulus_inv = _ip_data[ip].biot_modulus_inv;
390
391 auto const alpha_0 =
392 solid.property(MPL::PropertyType::biot_coefficient)
393 .template value<double>(vars, x_position, t, dt);
394 auto const porosity_0 =
395 medium->property(MPL::PropertyType::porosity)
396 .template value<double>(vars, x_position, t, dt);
397
398 auto const& D = _ip_data[ip].D;
399 auto const& P_sph = Invariants::spherical_projection;
400 auto const D_sph = P_sph * D * identity2;
401 double const bulk_modulus_eff = Invariants::trace(D_sph) / 9.;
402 double const bulk_modulus_degradation = bulk_modulus_eff / bulk_modulus;
403
404 // Update Biot's coefficient
405 biot_coefficient = 1. - bulk_modulus_degradation * (1. - alpha_0);
406
407 // Update Biot's modulus
408 biot_modulus_inv = bulk_modulus_degradation * (alpha_0 - porosity_0) *
409 (1. - alpha_0) / bulk_modulus;
410
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;
414 double const bulk_modulus_plus = Invariants::trace(C_tensile_sph) / 9.;
415
416 auto const driven_energy =
417 N.transpose() *
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) *
421 w;
422
423 local_Jac.noalias() += driven_energy * N * degradation_df2;
424
425 local_rhs.noalias() -= driven_energy * degradation_df1;
426
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,
431 _process_data.phasefield_model);
432 }
433}

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::biot_coefficient, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), displacement_index, MathLib::KelvinVector::Invariants< KelvinVectorSize >::identity2, phasefield_index, phasefield_size, MaterialPropertyLib::porosity, pressure_index, MaterialLib::Solids::selectSolidConstitutiveRelation(), ParameterLib::SpatialPosition::setElementID(), MathLib::KelvinVector::Invariants< KelvinVectorSize >::spherical_projection, and MathLib::KelvinVector::Invariants< KelvinVectorSize >::trace().

Referenced by assembleWithJacobianForStaggeredScheme().

◆ computeEnergy()

template<typename ShapeFunction, int DisplacementDim>
void ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::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 )
overridevirtual

Implements ProcessLib::HMPhaseField::HMPhaseFieldLocalAssemblerInterface.

Definition at line 532 of file HMPhaseFieldFEM-impl.h.

537{
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); });
544
545 auto const local_coupled_xs =
546 getCoupledLocalSolutions(x, indices_of_processes);
547 assert(local_coupled_xs.size() ==
548 phasefield_size + displacement_size + pressure_size);
549
550 auto const d = Eigen::Map<PhaseFieldVector const>(
551 &local_coupled_xs[phasefield_index], phasefield_size);
552 auto const u = Eigen::Map<DeformationVector const>(
553 &local_coupled_xs[displacement_index], displacement_size);
554 auto const p = Eigen::Map<PressureVector const>(
555 &local_coupled_xs[pressure_index], pressure_size);
556
557 ParameterLib::SpatialPosition x_position;
558 x_position.setElementID(_element.getID());
559
560 double element_elastic_energy = 0.0;
561 double element_surface_energy = 0.0;
562 double element_pressure_work = 0.0;
563
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];
566 int const n_integration_points = _integration_method.getNumberOfPoints();
567 for (int ip = 0; ip < n_integration_points; ip++)
568 {
569 auto const& w = _ip_data[ip].integration_weight;
570 auto const& N = _ip_data[ip].N;
571 auto const& dNdx = _ip_data[ip].dNdx;
572 double const d_ip = N.dot(d);
573 double const p_ip = N.dot(p);
574
575 element_elastic_energy += _ip_data[ip].elastic_energy * w;
576
577 switch (_process_data.phasefield_model)
578 {
579 case MaterialLib::Solids::Phasefield::PhaseFieldModel::AT1:
580 {
581 element_surface_energy +=
582 gc * 0.375 *
583 ((1 - d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) * w;
584
585 break;
586 }
587 case MaterialLib::Solids::Phasefield::PhaseFieldModel::AT2:
588 {
589 element_surface_energy += 0.5 * gc *
590 ((1 - d_ip) * (1 - d_ip) / ls +
591 (dNdx * d).dot((dNdx * d)) * ls) *
592 w;
593 break;
594 }
595 case MaterialLib::Solids::Phasefield::PhaseFieldModel::COHESIVE:
596 {
597 element_surface_energy +=
598 gc / std::numbers::pi *
599 ((1 - d_ip * d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) *
600 w;
601 break;
602 }
603 }
604
605 element_pressure_work += p_ip * (N_u_op(N) * u).dot(dNdx * d) * w;
606 }
607
608#ifdef USE_PETSC
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]),
612 [](GlobalIndexType const& index) { return index >= 0; });
613 if (n_all_nodes != n_regular_nodes)
614 {
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;
621 }
622#endif // USE_PETSC
623 elastic_energy += element_elastic_energy;
624 surface_energy += element_surface_energy;
625 pressure_work += element_pressure_work;
626}

References _element, _integration_method, _ip_data, _process_data, MaterialLib::Solids::Phasefield::AT1, MaterialLib::Solids::Phasefield::AT2, MaterialLib::Solids::Phasefield::COHESIVE, displacement_index, displacement_size, ProcessLib::getCoupledLocalSolutions(), N_u_op, phasefield_index, phasefield_size, pressure_index, pressure_size, and ParameterLib::SpatialPosition::setElementID().

◆ getIntPtEpsilon()

template<typename ShapeFunction, int DisplacementDim>
std::vector< double > const & ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::getIntPtEpsilon ( const double ,
std::vector< GlobalVector * > const & ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & ,
std::vector< double > & cache ) const
inlineoverrideprivatevirtual

Implements ProcessLib::HMPhaseField::HMPhaseFieldLocalAssemblerInterface.

Definition at line 322 of file HMPhaseFieldFEM.h.

327 {
328 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
329 _ip_data, &IpData::eps, cache);
330 }
ProcessLib::getIntegrationPointKelvinVectorData
std::vector< double > const & getIntegrationPointKelvinVectorData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
Definition SetOrGetIntegrationPointData.h:68

References _ip_data, ProcessLib::HMPhaseField::IntegrationPointData< BMatricesType, ShapeMatricesType, DisplacementDim >::eps, and ProcessLib::getIntegrationPointKelvinVectorData().

◆ getIntPtSigma()

template<typename ShapeFunction, int DisplacementDim>
std::vector< double > const & ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::getIntPtSigma ( const double ,
std::vector< GlobalVector * > const & ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & ,
std::vector< double > & cache ) const
inlineoverrideprivatevirtual

Implements ProcessLib::HMPhaseField::HMPhaseFieldLocalAssemblerInterface.

Definition at line 312 of file HMPhaseFieldFEM.h.

317 {
318 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
319 _ip_data, &IpData::sigma, cache);
320 }

References _ip_data, ProcessLib::getIntegrationPointKelvinVectorData(), and ProcessLib::HMPhaseField::IntegrationPointData< BMatricesType, ShapeMatricesType, DisplacementDim >::sigma.

◆ getIntPtWidth()

template<typename ShapeFunction, int DisplacementDim>
std::vector< double > const & ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::getIntPtWidth ( const double ,
std::vector< GlobalVector * > const & ,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & ,
std::vector< double > & cache ) const
overridevirtual

Implements ProcessLib::HMPhaseField::HMPhaseFieldLocalAssemblerInterface.

Definition at line 630 of file HMPhaseFieldFEM-impl.h.

635{
636 return ProcessLib::getIntegrationPointScalarData(_ip_data,
637 &IpData::width_ip, cache);
638}
ProcessLib::getIntegrationPointScalarData
std::vector< double > const & getIntegrationPointScalarData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
Definition SetOrGetIntegrationPointData.h:176

References _ip_data, ProcessLib::getIntegrationPointScalarData(), and ProcessLib::HMPhaseField::IntegrationPointData< BMatricesType, ShapeMatricesType, DisplacementDim >::width_ip.

◆ getShapeMatrix()

template<typename ShapeFunction, int DisplacementDim>
Eigen::Map< const Eigen::RowVectorXd > ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::getShapeMatrix ( const unsigned integration_point) const
inlineoverridevirtual

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 296 of file HMPhaseFieldFEM.h.

298 {
299 auto const& N = _secondary_data.N[integration_point];
300
301 // assumes N is stored contiguously in memory
302 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
303 }

References _secondary_data.

◆ heaviside()

template<typename ShapeFunction, int DisplacementDim>
double ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::heaviside ( double const v)
inline

Definition at line 294 of file HMPhaseFieldFEM.h.

294{ return (v < 0) ? 0.0 : 1.0; }

◆ initializeConcrete()

template<typename ShapeFunction, int DisplacementDim>
void ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::initializeConcrete ( )
inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 238 of file HMPhaseFieldFEM.h.

239 {
240 unsigned const n_integration_points =
241 _integration_method.getNumberOfPoints();
242
243 for (unsigned ip = 0; ip < n_integration_points; ip++)
244 {
245 _ip_data[ip].pushBackState();
246
247 // Specify the direction of preexisting fracture if it exists. The
248 // default value is zero.
249 auto& normal_ip = _ip_data[ip].normal_ip;
250 auto const fracture_normal =
251 _process_data.specific_fracture_direction;
252 normal_ip = fracture_normal;
253 }
254 // Specify the aperture of preexisting fractures if they exist. The
255 // default value is zero.
256 ParameterLib::SpatialPosition x_position;
257 x_position.setElementID(_element.getID());
258 auto const width_init = _process_data.width_init(0, x_position)[0];
259 (*_process_data.width)[_element.getID()] = width_init;
260 }

References _element, _integration_method, _ip_data, _process_data, and ParameterLib::SpatialPosition::setElementID().

◆ postNonLinearSolverConcrete()

template<typename ShapeFunction, int DisplacementDim>
void ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::postNonLinearSolverConcrete ( Eigen::VectorXd const & local_x,
Eigen::VectorXd const & local_x_prev,
double const t,
double const dt,
int const process_id )
overridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 436 of file HMPhaseFieldFEM-impl.h.

441{
442 int const n_integration_points = _integration_method.getNumberOfPoints();
443 auto const p = local_x.template segment<pressure_size>(pressure_index);
444 auto const p_prev =
445 local_x_prev.template segment<pressure_size>(pressure_index);
446
447 for (int ip = 0; ip < n_integration_points; ip++)
448 {
449 auto const& N = _ip_data[ip].N;
450 _ip_data[ip].coupling_pressure = N.dot(p - p_prev) / dt;
451 }
452}

References _integration_method, _ip_data, and pressure_index.

◆ postTimestepConcrete()

template<typename ShapeFunction, int DisplacementDim>
void ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::postTimestepConcrete ( Eigen::VectorXd const & ,
Eigen::VectorXd const & ,
double const ,
double const ,
int const  )
inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 262 of file HMPhaseFieldFEM.h.

266 {
267 unsigned const n_integration_points =
268 _integration_method.getNumberOfPoints();
269
270 for (unsigned ip = 0; ip < n_integration_points; ip++)
271 {
272 _ip_data[ip].pushBackState();
273 }
274 }

References _integration_method, and _ip_data.

Member Data Documentation

◆ _element

template<typename ShapeFunction, int DisplacementDim>
MeshLib::Element const& ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::_element
private

Definition at line 350 of file HMPhaseFieldFEM.h.

Referenced by HMPhaseFieldLocalAssembler(), approximateFractureWidth(), assembleWithJacobianForDeformationEquations(), assembleWithJacobianHydroEquations(), assembleWithJacobianPhaseFieldEquations(), computeEnergy(), and initializeConcrete().

◆ _integration_method

template<typename ShapeFunction, int DisplacementDim>
NumLib::GenericIntegrationMethod const& ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::_integration_method
private

Definition at line 349 of file HMPhaseFieldFEM.h.

Referenced by HMPhaseFieldLocalAssembler(), approximateFractureWidth(), assembleWithJacobianForDeformationEquations(), assembleWithJacobianHydroEquations(), assembleWithJacobianPhaseFieldEquations(), computeEnergy(), initializeConcrete(), postNonLinearSolverConcrete(), and postTimestepConcrete().

◆ _ip_data

template<typename ShapeFunction, int DisplacementDim>
std::vector<IpData, Eigen::aligned_allocator<IpData> > ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::_ip_data
private

Definition at line 347 of file HMPhaseFieldFEM.h.

Referenced by HMPhaseFieldLocalAssembler(), approximateFractureWidth(), assembleWithJacobianForDeformationEquations(), assembleWithJacobianHydroEquations(), assembleWithJacobianPhaseFieldEquations(), computeEnergy(), getIntPtEpsilon(), getIntPtSigma(), getIntPtWidth(), initializeConcrete(), postNonLinearSolverConcrete(), and postTimestepConcrete().

◆ _is_axially_symmetric

template<typename ShapeFunction, int DisplacementDim>
bool const ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::_is_axially_symmetric
private

Definition at line 352 of file HMPhaseFieldFEM.h.

Referenced by HMPhaseFieldLocalAssembler(), and assembleWithJacobianForDeformationEquations().

◆ _process_data

template<typename ShapeFunction, int DisplacementDim>
HMPhaseFieldProcessData<DisplacementDim>& ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::_process_data
private

Definition at line 345 of file HMPhaseFieldFEM.h.

Referenced by HMPhaseFieldLocalAssembler(), approximateFractureWidth(), assembleWithJacobianForDeformationEquations(), assembleWithJacobianForStaggeredScheme(), assembleWithJacobianHydroEquations(), assembleWithJacobianPhaseFieldEquations(), computeEnergy(), and initializeConcrete().

◆ _secondary_data

template<typename ShapeFunction, int DisplacementDim>
SecondaryData<typename ShapeMatrices::ShapeType> ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::_secondary_data
private

Definition at line 351 of file HMPhaseFieldFEM.h.

Referenced by HMPhaseFieldLocalAssembler(), and getShapeMatrix().

◆ displacement_index

template<typename ShapeFunction, int DisplacementDim>
int ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::displacement_index = pressure_index + pressure_size
staticconstexprprivate

Definition at line 120 of file HMPhaseFieldFEM.h.

Referenced by assembleWithJacobianForDeformationEquations(), assembleWithJacobianPhaseFieldEquations(), and computeEnergy().

◆ displacement_size

template<typename ShapeFunction, int DisplacementDim>
int ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::displacement_size
staticconstexprprivate
Initial value:
=
ShapeFunction::NPOINTS * DisplacementDim

Definition at line 121 of file HMPhaseFieldFEM.h.

Referenced by approximateFractureWidth(), assembleWithJacobianForDeformationEquations(), and computeEnergy().

◆ KelvinVectorSize

template<typename ShapeFunction, int DisplacementDim>
int const ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::KelvinVectorSize
static
Initial value:
=
MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim)

Definition at line 157 of file HMPhaseFieldFEM.h.

◆ N_u_op

template<typename ShapeFunction, int DisplacementDim>
auto& ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::N_u_op
staticconstexpr
Initial value:
= MathLib::eigenBlockMatrixView<
DisplacementDim, typename ShapeMatricesType::NodalRowVectorType>
MathLib::eigenBlockMatrixView
constexpr Eigen::CwiseNullaryOp< EigenBlockMatrixViewFunctor< D, M >, typename EigenBlockMatrixViewFunctor< D, M >::Matrix > eigenBlockMatrixView(const Eigen::MatrixBase< M > &matrix)
Definition EigenBlockMatrixView.h:49
EigenFixedShapeMatrixPolicy::NodalRowVectorType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
Definition ShapeMatrixPolicy.h:78

Definition at line 160 of file HMPhaseFieldFEM.h.

Referenced by assembleWithJacobianForDeformationEquations(), and computeEnergy().

◆ phasefield_index

template<typename ShapeFunction, int DisplacementDim>
int ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::phasefield_index = 0
staticconstexprprivate

Definition at line 116 of file HMPhaseFieldFEM.h.

Referenced by approximateFractureWidth(), assembleWithJacobianForDeformationEquations(), assembleWithJacobianHydroEquations(), assembleWithJacobianPhaseFieldEquations(), and computeEnergy().

◆ phasefield_size

template<typename ShapeFunction, int DisplacementDim>
int ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::phasefield_size = ShapeFunction::NPOINTS
staticconstexprprivate

Definition at line 117 of file HMPhaseFieldFEM.h.

Referenced by approximateFractureWidth(), assembleWithJacobianPhaseFieldEquations(), and computeEnergy().

◆ pressure_index

template<typename ShapeFunction, int DisplacementDim>
int ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::pressure_index = phasefield_index + phasefield_size
staticconstexprprivate

Definition at line 118 of file HMPhaseFieldFEM.h.

Referenced by assembleWithJacobianForDeformationEquations(), assembleWithJacobianHydroEquations(), assembleWithJacobianPhaseFieldEquations(), computeEnergy(), and postNonLinearSolverConcrete().

◆ pressure_size

template<typename ShapeFunction, int DisplacementDim>
int ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler< ShapeFunction, DisplacementDim >::pressure_size = ShapeFunction::NPOINTS
staticconstexprprivate

Definition at line 119 of file HMPhaseFieldFEM.h.

Referenced by approximateFractureWidth(), assembleWithJacobianHydroEquations(), and computeEnergy().

The documentation for this class was generated from the following files: