OGS
SmallDeformationFEM.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <limits>
14#include <memory>
15#include <optional>
16#include <vector>
17
22#include "NumLib/DOF/LocalDOF.h"
33
34namespace ProcessLib
35{
36namespace SmallDeformation
37{
38namespace MPL = MaterialPropertyLib;
39
40template <typename BMatricesType, typename ShapeMatricesType,
41 int DisplacementDim>
43{
45 typename ShapeMatricesType::NodalRowVectorType N_u;
46 typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx_u;
47
49};
50
53template <typename ShapeMatrixType>
55{
56 std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
57};
58
59template <typename ShapeFunction, int DisplacementDim>
61 : public SmallDeformationLocalAssemblerInterface<DisplacementDim>
62{
63public:
70
77
81 using IpData =
83
84 static constexpr auto& N_u_op = MathLib::eigenBlockMatrixView<
85 DisplacementDim, typename ShapeMatricesType::NodalRowVectorType>;
86
88 delete;
90
92 MeshLib::Element const& e,
93 std::size_t const /*local_matrix_size*/,
94 NumLib::GenericIntegrationMethod const& integration_method,
95 bool const is_axially_symmetric,
98 e, integration_method, is_axially_symmetric, process_data)
99 {
100 unsigned const n_integration_points =
102
103 _ip_data.resize(n_integration_points);
104 _secondary_data.N.resize(n_integration_points);
105
106 auto const shape_matrices =
108 DisplacementDim>(
109 e, is_axially_symmetric, this->integration_method_);
110
111 for (unsigned ip = 0; ip < n_integration_points; ip++)
112 {
113 auto& ip_data = _ip_data[ip];
114 auto const& sm = shape_matrices[ip];
115 _ip_data[ip].integration_weight =
117 sm.integralMeasure * sm.detJ;
118
119 ip_data.N_u = sm.N;
120 ip_data.dNdx_u = sm.dNdx;
121
122 _secondary_data.N[ip] = shape_matrices[ip].N;
123 }
124 }
125
126 void initializeConcrete() override
127 {
128 unsigned const n_integration_points =
130 for (unsigned ip = 0; ip < n_integration_points; ip++)
131 {
132 auto const& ip_data = _ip_data[ip];
133
134 ParameterLib::SpatialPosition const x_position{
135 std::nullopt, this->element_.getID(), ip,
137 NumLib::interpolateCoordinates<ShapeFunction,
139 this->element_, ip_data.N_u))};
140
142 if (this->process_data_.initial_stress != nullptr)
143 {
144 this->current_states_[ip].stress_data.sigma.noalias() =
146 DisplacementDim>((*this->process_data_.initial_stress)(
147 std::numeric_limits<
148 double>::quiet_NaN() /* time independent */,
149 x_position));
150 }
151
152 double const t = 0; // TODO (naumov) pass t from top
153 auto& material_state = this->material_states_[ip];
154 this->solid_material_.initializeInternalStateVariables(
155 t, x_position, *material_state.material_state_variables);
156
157 material_state.pushBackState();
158 }
159
160 for (unsigned ip = 0; ip < n_integration_points; ip++)
161 {
162 this->prev_states_[ip] = this->current_states_[ip];
163 }
164 }
165
168 BMatrixType const& B, Eigen::Ref<Eigen::VectorXd const> const& u,
169 Eigen::Ref<Eigen::VectorXd const> const& u_prev,
170 ParameterLib::SpatialPosition const& x_position, double const t,
171 double const dt,
173 CS,
174 MaterialPropertyLib::Medium const& medium,
176 current_state,
178 prev_state,
181 output_data) const
182 {
183 double const T_ref =
184 this->process_data_.reference_temperature
185 ? (*this->process_data_.reference_temperature)(t, x_position)[0]
186 : std::numeric_limits<double>::quiet_NaN();
187
189 models{this->process_data_, this->solid_material_};
191 tmp;
193
194 CS.eval(models, t, dt, x_position, //
195 medium, //
196 T_ref, B * u, B * u_prev, //
197 current_state, prev_state, material_state, tmp, output_data,
198 CD);
199
200 return CD;
201 }
202
203 void assemble(double const /*t*/, double const /*dt*/,
204 std::vector<double> const& /*local_x*/,
205 std::vector<double> const& /*local_x_prev*/,
206 std::vector<double>& /*local_M_data*/,
207 std::vector<double>& /*local_K_data*/,
208 std::vector<double>& /*local_b_data*/) override
209 {
210 OGS_FATAL(
211 "SmallDeformationLocalAssembler: assembly without jacobian is not "
212 "implemented.");
213 }
214
215 std::optional<BBarMatrixType> getDilatationalBBarMatrix() const
216 {
217 if (!(this->process_data_.use_b_bar))
218 {
219 return std::nullopt;
220 }
221
223 DisplacementDim, ShapeFunction::NPOINTS, ShapeFunction,
227 }
228
229 void assembleWithJacobian(double const t, double const dt,
230 std::vector<double> const& local_x,
231 std::vector<double> const& local_x_prev,
232 std::vector<double>& /*local_M_data*/,
233 std::vector<double>& /*local_K_data*/,
234 std::vector<double>& local_b_data,
235 std::vector<double>& local_Jac_data) override
236 {
237 auto const local_matrix_size = local_x.size();
238
239 auto local_Jac = MathLib::createZeroedMatrix<StiffnessMatrixType>(
240 local_Jac_data, local_matrix_size, local_matrix_size);
241
242 auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
243 local_b_data, local_matrix_size);
244
245 auto [u] = localDOF(local_x);
246 auto [u_prev] = localDOF(local_x_prev);
247
248 unsigned const n_integration_points =
250
252 x_position.setElementID(this->element_.getID());
253
256 auto const& medium =
257 *this->process_data_.media_map.getMedium(this->element_.getID());
258
259 auto const B_dil_bar = getDilatationalBBarMatrix();
260
261 for (unsigned ip = 0; ip < n_integration_points; ip++)
262 {
263 x_position.setIntegrationPoint(ip);
264 auto const& w = _ip_data[ip].integration_weight;
265 auto const& N = _ip_data[ip].N_u;
266 auto const& dNdx = _ip_data[ip].dNdx_u;
267
268 auto const x_coord =
269 NumLib::interpolateXCoordinate<ShapeFunction,
271 this->element_, N);
273 DisplacementDim, ShapeFunction::NPOINTS, BBarMatrixType,
275 dNdx, N, B_dil_bar, x_coord, this->is_axially_symmetric_);
276
277 auto const CD = updateConstitutiveRelations(
278 B, u, u_prev, x_position, t, dt, constitutive_setting, medium,
279 this->current_states_[ip], this->prev_states_[ip],
280 this->material_states_[ip], this->output_data_[ip]);
281
282 auto const& sigma = this->current_states_[ip].stress_data.sigma;
283 auto const& b = *CD.volumetric_body_force;
284 auto const& C = CD.s_mech_data.stiffness_tensor;
285
286 local_b.noalias() -=
287 (B.transpose() * sigma - N_u_op(N).transpose() * b) * w;
288 local_Jac.noalias() += B.transpose() * C * B * w;
289 }
290 }
291
292 void postTimestepConcrete(Eigen::VectorXd const& local_x,
293 Eigen::VectorXd const& local_x_prev,
294 double const t, double const dt,
295 int const /*process_id*/) override
296 {
297 unsigned const n_integration_points =
299
301 x_position.setElementID(this->element_.getID());
302
305
306 auto const& medium =
307 *this->process_data_.media_map.getMedium(this->element_.getID());
308
309 auto const B_dil_bar = getDilatationalBBarMatrix();
310
311 for (unsigned ip = 0; ip < n_integration_points; ip++)
312 {
313 x_position.setIntegrationPoint(ip);
314 auto const& N = _ip_data[ip].N_u;
315 auto const& dNdx = _ip_data[ip].dNdx_u;
316
317 auto const x_coord =
318 NumLib::interpolateXCoordinate<ShapeFunction,
320 this->element_, N);
322 DisplacementDim, ShapeFunction::NPOINTS, BBarMatrixType,
324 dNdx, N, B_dil_bar, x_coord, this->is_axially_symmetric_);
325
327 B, local_x, local_x_prev, x_position, t, dt,
328 constitutive_setting, medium, this->current_states_[ip],
329 this->prev_states_[ip], this->material_states_[ip],
330 this->output_data_[ip]);
331
332 this->material_states_[ip].pushBackState();
333 }
334
335 for (unsigned ip = 0; ip < n_integration_points; ip++)
336 {
337 this->prev_states_[ip] = this->current_states_[ip];
338 }
339 }
340
341 std::vector<double> const& getMaterialForces(
342 std::vector<double> const& local_x,
343 std::vector<double>& nodal_values) override
344 {
346 DisplacementDim, ShapeFunction, ShapeMatricesType,
349 GradientMatrixType>(local_x, nodal_values,
351 this->current_states_, this->output_data_,
352 this->element_, this->is_axially_symmetric_);
353 }
354
355 Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
356 const unsigned integration_point) const override
357 {
358 auto const& N = _secondary_data.N[integration_point];
359
360 // assumes N is stored contiguously in memory
361 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
362 }
363
364private:
365 static constexpr auto localDOF(std::vector<double> const& x)
366 {
367 return NumLib::localDOF<
369 }
370
371private:
372 std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
373
375};
376
377} // namespace SmallDeformation
378} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
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)
void setIntegrationPoint(unsigned integration_point)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
MatrixType< 3, ShapeFunction::NPOINTS > BBarMatrixType
MatrixType< _number_of_dof, _number_of_dof > StiffnessMatrixType
VectorType< _number_of_dof > NodalForceVectorType
Rhs residual.
VectorTypeFixedSize< MathLib::VectorizedTensor::size(DisplacementDim)> GradientVectorType
MatrixType< MathLib::VectorizedTensor::size(DisplacementDim), _number_of_dof > GradientMatrixType
typename BMatricesType::NodalForceVectorType NodalForceVectorType
typename ShapeMatricesType::NodalVectorType NodalVectorType
ConstitutiveRelations::ConstitutiveData< DisplacementDim > updateConstitutiveRelations(BMatrixType const &B, Eigen::Ref< Eigen::VectorXd const > const &u, Eigen::Ref< Eigen::VectorXd const > const &u_prev, ParameterLib::SpatialPosition const &x_position, double const t, double const dt, typename ConstitutiveRelations::ConstitutiveSetting< DisplacementDim > &CS, MaterialPropertyLib::Medium const &medium, typename ConstitutiveRelations::StatefulData< DisplacementDim > &current_state, typename ConstitutiveRelations::StatefulDataPrev< DisplacementDim > const &prev_state, MaterialStateData< DisplacementDim > &material_state, typename ConstitutiveRelations::OutputData< DisplacementDim > &output_data) const
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
std::optional< BBarMatrixType > getDilatationalBBarMatrix() const
SecondaryData< typename ShapeMatrices::ShapeType > _secondary_data
typename BMatricesType::NodalForceVectorType NodalDisplacementVectorType
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
void postTimestepConcrete(Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, double const t, double const dt, int const) override
typename ShapeMatricesType::NodalMatrixType NodalMatrixType
void assembleWithJacobian(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &, std::vector< double > &, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
typename GMatricesType::GradientMatrixType GradientMatrixType
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
ShapeMatrixPolicyType< ShapeFunction, DisplacementDim > ShapeMatricesType
std::vector< double > const & getMaterialForces(std::vector< double > const &local_x, std::vector< double > &nodal_values) override
typename BMatricesType::StiffnessMatrixType StiffnessMatrixType
SmallDeformationLocalAssembler(SmallDeformationLocalAssembler &&)=delete
static constexpr auto localDOF(std::vector< double > const &x)
typename GMatricesType::GradientVectorType GradientVectorType
void assemble(double const, double const, std::vector< double > const &, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
SmallDeformationLocalAssembler(MeshLib::Element const &e, std::size_t const, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, SmallDeformationProcessData< DisplacementDim > &process_data)
SmallDeformationLocalAssembler(SmallDeformationLocalAssembler const &)=delete
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
constexpr Eigen::CwiseNullaryOp< EigenBlockMatrixViewFunctor< D, M >, typename EigenBlockMatrixViewFunctor< D, M >::Matrix > eigenBlockMatrixView(const Eigen::MatrixBase< M > &matrix)
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
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)
auto localDOF(ElementDOFVector const &x)
Definition LocalDOF.h:64
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
BBarMatrixType computeDilatationalBbar(std::vector< IpData, Eigen::aligned_allocator< IpData > > const &ip_data, MeshLib::Element const &element, NumLib::GenericIntegrationMethod const &integration_method, const bool is_axially_symmetric)
BMatrixType computeBMatrixPossiblyWithBbar(DNDX_Type const &dNdx, N_Type const &N, std::optional< BBarMatrixType > const &B_dil_bar, const double radius, const bool is_axially_symmetric)
Fills a B matrix, or a B bar matrix if required.
std::vector< double > const & getMaterialForces(std::vector< double > const &local_x, std::vector< double > &nodal_values, IntegrationMethod const &_integration_method, IPData const &_ip_data, StressData const &stress_data, OutputData const &output_data, MeshLib::Element const &element, bool const is_axially_symmetric)
NumLib::ShapeMatrices< NodalRowVectorType, DimNodalMatrixType, DimMatrixType, GlobalDimNodalMatrixType > ShapeMatrices
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
VectorType< ShapeFunction::NPOINTS > NodalVectorType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
Data that is needed for the equation system assembly.
void eval(ConstitutiveModels< DisplacementDim > &models, double const t, double const dt, ParameterLib::SpatialPosition const &x_position, MaterialPropertyLib::Medium const &medium, double const T_ref, KelvinVector< DisplacementDim > const &eps, KelvinVector< DisplacementDim > const &eps_prev, StatefulData< DisplacementDim > &state, StatefulDataPrev< DisplacementDim > const &prev_state, MaterialStateData< DisplacementDim > &mat_state, ConstitutiveTempData< DisplacementDim > &tmp, OutputData< DisplacementDim > &out, ConstitutiveData< DisplacementDim > &cd) const
Evaluate the constitutive setting.
Data that is needed for output purposes, but not directly for the assembly.
Data whose state must be tracked by the process.
ShapeMatricesType::GlobalDimNodalMatrixType dNdx_u
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
ConstitutiveRelations::ConstitutiveSetting< DisplacementDim > constitutive_setting
std::vector< typename ConstitutiveRelations::StatefulDataPrev< DisplacementDim > > prev_states_
std::vector< typename ConstitutiveRelations::StatefulData< DisplacementDim > > current_states_
std::vector< typename ConstitutiveRelations::OutputData< DisplacementDim > > output_data_
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material_
std::vector< MaterialStateData< DisplacementDim > > material_states_