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 106 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 int getNumberOfVectorElementsForDeformation () 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 123 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 131 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 129 of file HMPhaseFieldFEM.h.

◆ GlobalDimVectorType

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

Definition at line 127 of file HMPhaseFieldFEM.h.

◆ Invariants

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

Definition at line 156 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 147 of file HMPhaseFieldFEM.h.

◆ NodalForceVectorType

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

Definition at line 125 of file HMPhaseFieldFEM.h.

◆ NodalMatrixType

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

Definition at line 145 of file HMPhaseFieldFEM.h.

◆ NodalVectorType

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

Definition at line 146 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 136 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 134 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 141 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 139 of file HMPhaseFieldFEM.h.

◆ ShapeMatrices

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

Definition at line 121 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:134

Definition at line 118 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 161 of file HMPhaseFieldFEM.h.

167 : _process_data(process_data),
168 _integration_method(integration_method),
169 _element(e),
170 _is_axially_symmetric(is_axially_symmetric)
171 {
172 unsigned const n_integration_points =
173 _integration_method.getNumberOfPoints();
174
175 _ip_data.reserve(n_integration_points);
176 _secondary_data.N.resize(n_integration_points);
177
178 auto& solid_material =
179 MaterialLib::Solids::selectSolidConstitutiveRelation(
180 _process_data.solid_materials,
181 _process_data.material_ids,
182 e.getID());
183
184 auto const shape_matrices =
185 NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
186 DisplacementDim>(e, is_axially_symmetric,
187 _integration_method);
188
189 ParameterLib::SpatialPosition x_position;
190 x_position.setElementID(_element.getID());
191
192 for (unsigned ip = 0; ip < n_integration_points; ip++)
193 {
194 _ip_data.emplace_back(solid_material);
195 auto& ip_data = _ip_data[ip];
196 ip_data.integration_weight =
197 _integration_method.getWeightedPoint(ip).getWeight() *
198 shape_matrices[ip].integralMeasure * shape_matrices[ip].detJ;
199
200 static const int kelvin_vector_size =
201 MathLib::KelvinVector::kelvin_vector_dimensions(
202 DisplacementDim);
203 ip_data.eps_tensile.setZero(kelvin_vector_size);
204 ip_data.eps.setZero(kelvin_vector_size);
205 ip_data.eps_prev.resize(kelvin_vector_size);
206 ip_data.D.setZero(kelvin_vector_size, kelvin_vector_size);
207 ip_data.C_tensile.setZero(kelvin_vector_size, kelvin_vector_size);
208 ip_data.C_compressive.setZero(kelvin_vector_size,
209 kelvin_vector_size);
210
211 ip_data.sigma_tensile.setZero(kelvin_vector_size);
212 ip_data.sigma_compressive.setZero(kelvin_vector_size);
213 ip_data.sigma.setZero(kelvin_vector_size);
214 ip_data.strain_energy_tensile = 0.0;
215 ip_data.elastic_energy = 0.0;
216 ip_data.width_ip = 0.0;
217
218 ip_data.N = shape_matrices[ip].N;
219 ip_data.dNdx = shape_matrices[ip].dNdx;
220
221 _secondary_data.N[ip] = shape_matrices[ip].N;
222 }
223 }
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler::_ip_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
Definition HMPhaseFieldFEM.h:340
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler::_process_data
HMPhaseFieldProcessData< DisplacementDim > & _process_data
Definition HMPhaseFieldFEM.h:338
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler::ShapeMatricesType
ShapeMatrixPolicyType< ShapeFunction, DisplacementDim > ShapeMatricesType
Definition HMPhaseFieldFEM.h:118
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler::_integration_method
NumLib::GenericIntegrationMethod const & _integration_method
Definition HMPhaseFieldFEM.h:342
ProcessLib::HMPhaseField::HMPhaseFieldLocalAssembler::_secondary_data
SecondaryData< typename ShapeMatrices::ShapeType > _secondary_data
Definition HMPhaseFieldFEM.h:344
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:30
MathLib::KelvinVector::kelvin_vector_dimensions
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition KelvinVector.h:18

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 451 of file HMPhaseFieldFEM-impl.h.

457{
458 std::vector<std::vector<GlobalIndexType>> indices_of_processes;
459 indices_of_processes.reserve(dof_tables.size());
460 std::transform(dof_tables.begin(), dof_tables.end(),
461 std::back_inserter(indices_of_processes),
462 [&](auto const dof_table)
463 { return NumLib::getIndices(mesh_item_id, *dof_table); });
464
465 auto local_coupled_xs = getCoupledLocalSolutions(x, indices_of_processes);
466 assert(local_coupled_xs.size() ==
467 phasefield_size + displacement_size + pressure_size);
468
469 auto const d = Eigen::Map<PhaseFieldVector const>(
470 &local_coupled_xs[phasefield_index], phasefield_size);
471
472 ParameterLib::SpatialPosition x_position;
473 x_position.setElementID(_element.getID());
474
475 double const ele_d = std::clamp(d.sum() / d.size(), 0.0, 1.0);
476 (*_process_data.ele_d)[_element.getID()] = ele_d;
477
478 if ((*_process_data.ele_d)[_element.getID()] <
479 _process_data.irreversible_threshold)
480 {
481 double const width_init = _process_data.width_init(t, x_position)[0];
482 double const k = _process_data.residual_stiffness(t, x_position)[0];
483 double const ls = _process_data.crack_length_scale(t, x_position)[0];
484 double const he = ls / _process_data.diffused_range_parameter;
485
486 int const n_integration_points =
487 _integration_method.getNumberOfPoints();
488 double width = 0.0;
489 for (int ip = 0; ip < n_integration_points; ip++)
490 {
491 auto eps_tensor =
492 MathLib::KelvinVector::kelvinVectorToTensor(_ip_data[ip].eps);
493 Eigen::EigenSolver<decltype(eps_tensor)> eigen_solver(eps_tensor);
494 Eigen::MatrixXf::Index maxIndex;
495 double const max_principal_strain =
496 eigen_solver.eigenvalues().real().maxCoeff(&maxIndex);
497 auto const max_eigen_vector =
498 eigen_solver.eigenvectors().real().col(maxIndex);
499
500 // Fracture aperture estimation
501 auto& width_ip = _ip_data[ip].width_ip;
502 width_ip = max_principal_strain * he;
503 width_ip = width_ip < width_init ? width_init : width_ip;
504 width += width_ip;
505
506 // Fracture direction estimation
507 auto& normal_ip = _ip_data[ip].normal_ip;
508 if (std::abs(max_principal_strain) > k)
509 {
510 for (int i = 0; i < DisplacementDim; i++)
511 {
512 normal_ip[i] = max_eigen_vector[i];
513 }
514 }
515
516 // Fracture enhanced porosity
517 auto& fracture_enhanced_porosity =
518 _ip_data[ip].fracture_enhanced_porosity;
519 fracture_enhanced_porosity = width_ip / he;
520 }
521
522 // Update aperture for the fractured element
523 (*_process_data.width)[_element.getID()] = width / n_integration_points;
524 }
525}
MathLib::KelvinVector::kelvinVectorToTensor
Eigen::Matrix< double, 3, 3 > kelvinVectorToTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
Definition KelvinVector.cpp:57
ProcessLib::getCoupledLocalSolutions
std::vector< double > getCoupledLocalSolutions(std::vector< GlobalVector * > const &global_solutions, std::vector< std::vector< GlobalIndexType > > const &indices)
Definition CoupledSolutionsForStaggeredScheme.cpp:12

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 55 of file HMPhaseFieldFEM-impl.h.

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

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, ParameterLib::SpatialPosition::getCoordinates(), MathLib::KelvinVector::Invariants< KelvinVectorSize >::identity2, NumLib::interpolateCoordinates(), N_u_op, phasefield_index, MaterialPropertyLib::porosity, pressure_index, MaterialLib::Solids::selectSolidConstitutiveRelation(), MaterialPropertyLib::Solid, 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 25 of file HMPhaseFieldFEM-impl.h.

32{
33 // For the equations with phase field.
34 if (process_id == _process_data._phasefield_process_id)
35 {
36 assembleWithJacobianPhaseFieldEquations(t, dt, local_x, local_b_data,
37 local_Jac_data);
38 return;
39 }
40
41 // For the equations for hydro
42 if (process_id == _process_data._hydro_process_id)
43 {
44 assembleWithJacobianHydroEquations(t, dt, local_x, local_x_prev,
45 local_b_data, local_Jac_data);
46 return;
47 }
48
49 // For the equations with deformation
50 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_b_data,
51 local_Jac_data);
52}
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:330
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:171
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:56

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 170 of file HMPhaseFieldFEM-impl.h.

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

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 329 of file HMPhaseFieldFEM-impl.h.

334{
335 auto const d = local_x.template segment<phasefield_size>(phasefield_index);
336 auto const p = local_x.template segment<pressure_size>(pressure_index);
337 auto const u =
338 local_x.template segment<displacement_size>(displacement_index);
339
340 auto local_Jac = MathLib::createZeroedMatrix<PhaseFieldMatrix>(
341 local_Jac_data, phasefield_size, phasefield_size);
342 auto local_rhs = MathLib::createZeroedVector<PhaseFieldVector>(
343 local_b_data, phasefield_size);
344
345 ParameterLib::SpatialPosition x_position;
346 x_position.setElementID(_element.getID());
347
348 auto const& solid_material =
349 MaterialLib::Solids::selectSolidConstitutiveRelation(
350 _process_data.solid_materials, _process_data.material_ids,
351 _element.getID());
352
353 auto const bulk_modulus = solid_material.getBulkModulus(t, x_position);
354 auto const& medium = _process_data.media_map.getMedium(_element.getID());
355 auto const& solid = medium->phase(MaterialPropertyLib::PhaseName::Solid);
356 MPL::VariableArray vars;
357
358 double const k = _process_data.residual_stiffness(t, x_position)[0];
359 double const ls = _process_data.crack_length_scale(t, x_position)[0];
360 double const gc = _process_data.crack_resistance(t, x_position)[0];
361
362 auto const& identity2 = Invariants::identity2;
363
364 int const n_integration_points = _integration_method.getNumberOfPoints();
365 for (int ip = 0; ip < n_integration_points; ip++)
366 {
367 auto const& w = _ip_data[ip].integration_weight;
368 auto const& N = _ip_data[ip].N;
369 auto const& dNdx = _ip_data[ip].dNdx;
370
371 double const d_ip = N.dot(d);
372 double const p_ip = N.dot(p);
373 double const degradation =
374 _process_data.degradation_derivative->degradation(d_ip, k, ls);
375 double const degradation_df1 =
376 _process_data.degradation_derivative->degradationDf1(d_ip, k, ls);
377 double const degradation_df2 =
378 _process_data.degradation_derivative->degradationDf2(d_ip, k, ls);
379
380 _ip_data[ip].updateConstitutiveRelation(
381 t, x_position, dt, u, degradation,
382 _process_data.energy_split_model);
383
384 auto& biot_coefficient = _ip_data[ip].biot_coefficient;
385 auto& biot_modulus_inv = _ip_data[ip].biot_modulus_inv;
386
387 auto const alpha_0 =
388 solid.property(MPL::PropertyType::biot_coefficient)
389 .template value<double>(vars, x_position, t, dt);
390 auto const porosity_0 =
391 medium->property(MPL::PropertyType::porosity)
392 .template value<double>(vars, x_position, t, dt);
393
394 auto const& D = _ip_data[ip].D;
395 auto const& P_sph = Invariants::spherical_projection;
396 auto const D_sph = P_sph * D * identity2;
397 double const bulk_modulus_eff = Invariants::trace(D_sph) / 9.;
398 double const bulk_modulus_degradation = bulk_modulus_eff / bulk_modulus;
399
400 // Update Biot's coefficient
401 biot_coefficient = 1. - bulk_modulus_degradation * (1. - alpha_0);
402
403 // Update Biot's modulus
404 biot_modulus_inv = bulk_modulus_degradation * (alpha_0 - porosity_0) *
405 (1. - alpha_0) / bulk_modulus;
406
407 auto const& strain_energy_tensile = _ip_data[ip].strain_energy_tensile;
408 auto const& C_tensile = _ip_data[ip].C_tensile;
409 auto const C_tensile_sph = P_sph * C_tensile * identity2;
410 double const bulk_modulus_plus = Invariants::trace(C_tensile_sph) / 9.;
411
412 auto const driven_energy =
413 N.transpose() *
414 (strain_energy_tensile + p_ip * p_ip / 2. * bulk_modulus_plus /
415 bulk_modulus * (alpha_0 - porosity_0) *
416 (1. - alpha_0) / bulk_modulus) *
417 w;
418
419 local_Jac.noalias() += driven_energy * N * degradation_df2;
420
421 local_rhs.noalias() -= driven_energy * degradation_df1;
422
423 calculateCrackLocalJacobianAndResidual<
424 decltype(dNdx), decltype(N), decltype(w), decltype(d),
425 decltype(local_Jac), decltype(local_rhs)>(
426 dNdx, N, w, d, local_Jac, local_rhs, gc, ls,
427 _process_data.phasefield_model);
428 }
429}

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(), MaterialPropertyLib::Solid, 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 528 of file HMPhaseFieldFEM-impl.h.

533{
534 std::vector<std::vector<GlobalIndexType>> indices_of_processes;
535 indices_of_processes.reserve(dof_tables.size());
536 std::transform(dof_tables.begin(), dof_tables.end(),
537 std::back_inserter(indices_of_processes),
538 [&](auto const dof_table)
539 { return NumLib::getIndices(mesh_item_id, *dof_table); });
540
541 auto const local_coupled_xs =
542 getCoupledLocalSolutions(x, indices_of_processes);
543 assert(local_coupled_xs.size() ==
544 phasefield_size + displacement_size + pressure_size);
545
546 auto const d = Eigen::Map<PhaseFieldVector const>(
547 &local_coupled_xs[phasefield_index], phasefield_size);
548 auto const u = Eigen::Map<DeformationVector const>(
549 &local_coupled_xs[displacement_index], displacement_size);
550 auto const p = Eigen::Map<PressureVector const>(
551 &local_coupled_xs[pressure_index], pressure_size);
552
553 ParameterLib::SpatialPosition x_position;
554 x_position.setElementID(_element.getID());
555
556 double element_elastic_energy = 0.0;
557 double element_surface_energy = 0.0;
558 double element_pressure_work = 0.0;
559
560 double const gc = _process_data.crack_resistance(t, x_position)[0];
561 double const ls = _process_data.crack_length_scale(t, x_position)[0];
562 int const n_integration_points = _integration_method.getNumberOfPoints();
563 for (int ip = 0; ip < n_integration_points; ip++)
564 {
565 auto const& w = _ip_data[ip].integration_weight;
566 auto const& N = _ip_data[ip].N;
567 auto const& dNdx = _ip_data[ip].dNdx;
568 double const d_ip = N.dot(d);
569 double const p_ip = N.dot(p);
570
571 element_elastic_energy += _ip_data[ip].elastic_energy * w;
572
573 switch (_process_data.phasefield_model)
574 {
575 case MaterialLib::Solids::Phasefield::PhaseFieldModel::AT1:
576 {
577 element_surface_energy +=
578 gc * 0.375 *
579 ((1 - d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) * w;
580
581 break;
582 }
583 case MaterialLib::Solids::Phasefield::PhaseFieldModel::AT2:
584 {
585 element_surface_energy += 0.5 * gc *
586 ((1 - d_ip) * (1 - d_ip) / ls +
587 (dNdx * d).dot((dNdx * d)) * ls) *
588 w;
589 break;
590 }
591 case MaterialLib::Solids::Phasefield::PhaseFieldModel::COHESIVE:
592 {
593 element_surface_energy +=
594 gc / std::numbers::pi *
595 ((1 - d_ip * d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) *
596 w;
597 break;
598 }
599 }
600
601 element_pressure_work += p_ip * (N_u_op(N) * u).dot(dNdx * d) * w;
602 }
603
604#ifdef USE_PETSC
605 int const n_all_nodes = indices_of_processes[1].size();
606 int const n_regular_nodes = std::count_if(
607 begin(indices_of_processes[1]), end(indices_of_processes[1]),
608 [](GlobalIndexType const& index) { return index >= 0; });
609 if (n_all_nodes != n_regular_nodes)
610 {
611 element_elastic_energy *=
612 static_cast<double>(n_regular_nodes) / n_all_nodes;
613 element_surface_energy *=
614 static_cast<double>(n_regular_nodes) / n_all_nodes;
615 element_pressure_work *=
616 static_cast<double>(n_regular_nodes) / n_all_nodes;
617 }
618#endif // USE_PETSC
619 elastic_energy += element_elastic_energy;
620 surface_energy += element_surface_energy;
621 pressure_work += element_pressure_work;
622}

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 315 of file HMPhaseFieldFEM.h.

320 {
321 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
322 _ip_data, &IpData::eps, cache);
323 }
ProcessLib::getIntegrationPointKelvinVectorData
std::vector< double > const & getIntegrationPointKelvinVectorData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
Definition SetOrGetIntegrationPointData.h:60

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 305 of file HMPhaseFieldFEM.h.

310 {
311 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
312 _ip_data, &IpData::sigma, cache);
313 }

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 626 of file HMPhaseFieldFEM-impl.h.

631{
632 return ProcessLib::getIntegrationPointScalarData(_ip_data,
633 &IpData::width_ip, cache);
634}
ProcessLib::getIntegrationPointScalarData
std::vector< double > const & getIntegrationPointScalarData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
Definition SetOrGetIntegrationPointData.h:168

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 289 of file HMPhaseFieldFEM.h.

291 {
292 auto const& N = _secondary_data.N[integration_point];
293
294 // assumes N is stored contiguously in memory
295 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
296 }

References _secondary_data.

◆ heaviside()

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

Definition at line 287 of file HMPhaseFieldFEM.h.

287{ 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 231 of file HMPhaseFieldFEM.h.

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

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 432 of file HMPhaseFieldFEM-impl.h.

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

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 255 of file HMPhaseFieldFEM.h.

259 {
260 unsigned const n_integration_points =
261 _integration_method.getNumberOfPoints();
262
263 for (unsigned ip = 0; ip < n_integration_points; ip++)
264 {
265 _ip_data[ip].pushBackState();
266 }
267 }

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 343 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 342 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 340 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 345 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 338 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 344 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 113 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 114 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 150 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:43
EigenFixedShapeMatrixPolicy::NodalRowVectorType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
Definition ShapeMatrixPolicy.h:71

Definition at line 153 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 109 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 110 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 111 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 112 of file HMPhaseFieldFEM.h.

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

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