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

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

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

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

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

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

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

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

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

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

630{
631 return ProcessLib::getIntegrationPointScalarData(_ip_data,
632 &IpData::width_ip, cache);
633}
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 431 of file HMPhaseFieldFEM-impl.h.

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

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: