Loading [MathJax]/extensions/MathMenu.js
OGS
HMPhaseFieldFEM.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <memory>
14#include <vector>
15
36
37namespace ProcessLib
38{
39namespace HMPhaseField
40{
41namespace MPL = MaterialPropertyLib;
42template <typename BMatricesType, typename ShapeMatrixType, int DisplacementDim>
44{
53
54 typename ShapeMatrixType::NodalRowVectorType N;
55 typename ShapeMatrixType::GlobalDimNodalMatrixType dNdx;
56
57 typename BMatricesType::KelvinVectorType eps, eps_prev, eps_tensile;
58
59 typename BMatricesType::KelvinVectorType sigma_tensile, sigma_compressive,
62 double width_ip;
63 Eigen::Vector<double, DisplacementDim> normal_ip;
64
66 std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
67 DisplacementDim>::MaterialStateVariables>
69
70 typename BMatricesType::KelvinMatrixType D, C_tensile, C_compressive;
71
73 double coupling_pressure = 0.0;
77
85
86 template <typename DisplacementVectorType>
88 double const t,
90 double const /*dt*/,
91 DisplacementVectorType const& /*u*/,
92 double const degradation,
94 energy_split_model)
95 {
97 decltype(sigma), decltype(D), DisplacementDim>(
100 eps, energy_split_model, t, x, solid_material);
101 }
102};
103
106template <typename ShapeMatrixType>
108{
109 std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
110};
111
112template <typename ShapeFunction, int DisplacementDim>
114{
115private:
116 static constexpr int phasefield_index = 0;
117 static constexpr int phasefield_size = ShapeFunction::NPOINTS;
119 static constexpr int pressure_size = ShapeFunction::NPOINTS;
121 static constexpr int displacement_size =
122 ShapeFunction::NPOINTS * DisplacementDim;
123
124public:
127
129
131
133
135
137 typename ShapeMatricesType::template VectorType<displacement_size>;
139 typename ShapeMatricesType::template MatrixType<displacement_size,
142 typename ShapeMatricesType::template VectorType<phasefield_size>;
144 typename ShapeMatricesType::template MatrixType<phasefield_size,
147 typename ShapeMatricesType::template VectorType<pressure_size>;
149 typename ShapeMatricesType::template MatrixType<pressure_size,
151
154 using IpData =
156
157 static int const KelvinVectorSize =
159
160 static constexpr auto& N_u_op = MathLib::eigenBlockMatrixView<
161 DisplacementDim, typename ShapeMatricesType::NodalRowVectorType>;
162
164
167
169 MeshLib::Element const& e,
170 std::size_t const /*local_matrix_size*/,
171 NumLib::GenericIntegrationMethod const& integration_method,
172 bool const is_axially_symmetric,
174 : _process_data(process_data),
175 _integration_method(integration_method),
176 _element(e),
177 _is_axially_symmetric(is_axially_symmetric)
178 {
179 unsigned const n_integration_points =
181
182 _ip_data.reserve(n_integration_points);
183 _secondary_data.N.resize(n_integration_points);
184
185 auto& solid_material =
187 _process_data.solid_materials,
188 _process_data.material_ids,
189 e.getID());
190
191 auto const shape_matrices =
193 DisplacementDim>(e, is_axially_symmetric,
195
197 x_position.setElementID(_element.getID());
198
199 for (unsigned ip = 0; ip < n_integration_points; ip++)
200 {
201 _ip_data.emplace_back(solid_material);
202 auto& ip_data = _ip_data[ip];
203 ip_data.integration_weight =
205 shape_matrices[ip].integralMeasure * shape_matrices[ip].detJ;
206
207 static const int kelvin_vector_size =
209 DisplacementDim);
210 ip_data.eps_tensile.setZero(kelvin_vector_size);
211 ip_data.eps.setZero(kelvin_vector_size);
212 ip_data.eps_prev.resize(kelvin_vector_size);
213 ip_data.D.setZero(kelvin_vector_size, kelvin_vector_size);
214 ip_data.C_tensile.setZero(kelvin_vector_size, kelvin_vector_size);
215 ip_data.C_compressive.setZero(kelvin_vector_size,
216 kelvin_vector_size);
217
218 ip_data.sigma_tensile.setZero(kelvin_vector_size);
219 ip_data.sigma_compressive.setZero(kelvin_vector_size);
220 ip_data.sigma.setZero(kelvin_vector_size);
221 ip_data.strain_energy_tensile = 0.0;
222 ip_data.elastic_energy = 0.0;
223 ip_data.width_ip = 0.0;
224
225 ip_data.N = shape_matrices[ip].N;
226 ip_data.dNdx = shape_matrices[ip].dNdx;
227
228 _secondary_data.N[ip] = shape_matrices[ip].N;
229 }
230 }
231
233 double const t, double const dt, Eigen::VectorXd const& local_x,
234 Eigen::VectorXd const& local_x_prev, int const process_id,
235 std::vector<double>& local_b_data,
236 std::vector<double>& local_Jac_data) override;
237
238 void initializeConcrete() override
239 {
240 unsigned const n_integration_points =
242
243 for (unsigned ip = 0; ip < n_integration_points; ip++)
244 {
245 _ip_data[ip].pushBackState();
246
247 // Specify the direction of preexisting fracture if it exists. The
248 // default value is zero.
249 auto& normal_ip = _ip_data[ip].normal_ip;
250 auto const fracture_normal =
251 _process_data.specific_fracture_direction;
252 normal_ip = fracture_normal;
253 }
254 // Specify the aperture of preexisting fractures if they exist. The
255 // default value is zero.
257 x_position.setElementID(_element.getID());
258 auto const width_init = _process_data.width_init(0, x_position)[0];
259 (*_process_data.width)[_element.getID()] = width_init;
260 }
261
262 void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
263 Eigen::VectorXd const& /*local_x_prev*/,
264 double const /*t*/, double const /*dt*/,
265 int const /*process_id*/) override
266 {
267 unsigned const n_integration_points =
269
270 for (unsigned ip = 0; ip < n_integration_points; ip++)
271 {
272 _ip_data[ip].pushBackState();
273 }
274 }
275
276 void postNonLinearSolverConcrete(Eigen::VectorXd const& local_x,
277 Eigen::VectorXd const& local_x_prev,
278 double const t, double const dt,
279 int const process_id) override;
280
282 std::size_t mesh_item_id,
283 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
284 std::vector<GlobalVector*> const& x, double const t,
285 double const dt) override;
286
287 void computeEnergy(
288 std::size_t mesh_item_id,
289 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
290 std::vector<GlobalVector*> const& x, double const t,
291 double& elastic_energy, double& surface_energy,
292 double& pressure_work) override;
293
294 inline double heaviside(double const v) { return (v < 0) ? 0.0 : 1.0; }
295
296 Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
297 const unsigned integration_point) const override
298 {
299 auto const& N = _secondary_data.N[integration_point];
300
301 // assumes N is stored contiguously in memory
302 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
303 }
304
305 std::vector<double> const& getIntPtWidth(
306 const double /*t*/,
307 std::vector<GlobalVector*> const& /*x*/,
308 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
309 std::vector<double>& cache) const override;
310
311private:
312 std::vector<double> const& getIntPtSigma(
313 const double /*t*/,
314 std::vector<GlobalVector*> const& /*x*/,
315 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
316 std::vector<double>& cache) const override
317 {
319 _ip_data, &IpData::sigma, cache);
320 }
321
322 std::vector<double> const& getIntPtEpsilon(
323 const double /*t*/,
324 std::vector<GlobalVector*> const& /*x*/,
325 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
326 std::vector<double>& cache) const override
327 {
329 _ip_data, &IpData::eps, cache);
330 }
331
333 const double t, double const dt, Eigen::VectorXd const& local_x,
334 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
335 std::vector<double>& local_Jac_data);
336
338 double const t, double const dt, Eigen::VectorXd const& local_x,
339 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data);
340
342 double const t, double const dt, Eigen::VectorXd const& local_x,
343 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data);
344
346
347 std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
348
353};
354
355} // namespace HMPhaseField
356} // namespace ProcessLib
357
358#include "HMPhaseFieldFEM-impl.h"
Definition of the Element class.
constexpr double getWeight() const
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
void setElementID(std::size_t element_id)
VectorType< _number_of_dof > NodalForceVectorType
Rhs residual.
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
typename BMatricesType::NodalForceVectorType NodalForceVectorType
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
void computeEnergy(std::size_t mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double &elastic_energy, double &surface_energy, double &pressure_work) override
HMPhaseFieldLocalAssembler(HMPhaseFieldLocalAssembler &&)=delete
HMPhaseFieldLocalAssembler(HMPhaseFieldLocalAssembler const &)=delete
std::vector< double > const & getIntPtSigma(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
typename ShapeMatricesType::template VectorType< displacement_size > DeformationVector
typename ShapeMatricesType::NodalVectorType NodalVectorType
void assembleWithJacobianForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
HMPhaseFieldProcessData< DisplacementDim > & _process_data
typename ShapeMatricesType::NodalMatrixType NodalMatrixType
void postNonLinearSolverConcrete(Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, int const process_id) override
std::vector< double > const & getIntPtWidth(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
ShapeMatrixPolicyType< ShapeFunction, DisplacementDim > ShapeMatricesType
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)
typename ShapeMatricesType::template MatrixType< phasefield_size, phasefield_size > PhaseFieldMatrix
void assembleWithJacobianHydroEquations(const double t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void assembleWithJacobianForDeformationEquations(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
typename ShapeMatricesType::template MatrixType< pressure_size, pressure_size > PressureMatrix
typename ShapeMatricesType::template MatrixType< displacement_size, displacement_size > DeformationMatrix
NumLib::GenericIntegrationMethod const & _integration_method
typename ShapeMatricesType::template VectorType< pressure_size > PressureVector
typename ShapeMatricesType::template VectorType< phasefield_size > PhaseFieldVector
void postTimestepConcrete(Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, int const) override
HMPhaseFieldLocalAssembler(MeshLib::Element const &e, std::size_t const, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, HMPhaseFieldProcessData< DisplacementDim > &process_data)
SecondaryData< typename ShapeMatrices::ShapeType > _secondary_data
std::vector< double > const & getIntPtEpsilon(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
void calculateStress(T_VECTOR &sigma, T_VECTOR &sigma_tensile, T_VECTOR &sigma_compressive, T_VECTOR &eps_tensile, T_MATRIX &D, T_MATRIX &C_tensile, T_MATRIX &C_compressive, double &strain_energy_tensile, double &elastic_energy, double const degradation, T_VECTOR const &eps, EnergySplitModel const &energy_split_model, double const t, ParameterLib::SpatialPosition const &x, MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
auto & selectSolidConstitutiveRelation(SolidMaterialsMap const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
constexpr Eigen::CwiseNullaryOp< EigenBlockMatrixViewFunctor< D, M >, typename EigenBlockMatrixViewFunctor< D, M >::Matrix > eigenBlockMatrixView(const Eigen::MatrixBase< M > &matrix)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
std::vector< double > const & getIntegrationPointKelvinVectorData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
NumLib::ShapeMatrices< NodalRowVectorType, DimNodalMatrixType, DimMatrixType, GlobalDimNodalMatrixType > ShapeMatrices
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
VectorType< GlobalDim > GlobalDimVectorType
VectorType< ShapeFunction::NPOINTS > NodalVectorType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
BMatricesType::KelvinMatrixType C_compressive
void updateConstitutiveRelation(double const t, ParameterLib::SpatialPosition const &x, double const, DisplacementVectorType const &, double const degradation, MaterialLib::Solids::Phasefield::EnergySplitModel const energy_split_model)
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material
BMatricesType::KelvinMatrixType C_tensile
BMatricesType::KelvinVectorType eps_prev
BMatricesType::KelvinVectorType sigma_compressive
BMatricesType::KelvinVectorType eps_tensile
std::unique_ptr< typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables > material_state_variables
IntegrationPointData(MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
Eigen::Vector< double, DisplacementDim > normal_ip
ShapeMatrixType::NodalRowVectorType N
ShapeMatrixType::GlobalDimNodalMatrixType dNdx
BMatricesType::KelvinVectorType sigma_tensile
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N