OGS
ThermoMechanicsFEM.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <memory>
14#include <vector>
15
29
30namespace ProcessLib
31{
32namespace ThermoMechanics
33{
34namespace MPL = MaterialPropertyLib;
35
36template <typename BMatricesType, typename ShapeMatricesType,
37 int DisplacementDim>
39{
48
49 // total stress
50 typename BMatricesType::KelvinVectorType sigma, sigma_prev;
51
52 // total strain
53 typename BMatricesType::KelvinVectorType eps, eps_prev;
54
55 // mechanical strain
56 typename BMatricesType::KelvinVectorType eps_m, eps_m_prev;
57
59 std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
60 DisplacementDim>::MaterialStateVariables>
62
64 typename ShapeMatricesType::NodalRowVectorType N;
65 typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx;
66
68 {
69 eps_prev = eps;
72 material_state_variables->pushBackState();
73 }
74
76};
77
80template <typename ShapeMatrixType>
82{
83 std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
84};
85
126template <typename ShapeFunction, int DisplacementDim>
128 : public ThermoMechanicsLocalAssemblerInterface<DisplacementDim>
129{
130public:
133
134 // Types for displacement.
135 // (Higher order elements = ShapeFunction).
139
140 static int const KelvinVectorSize =
143
145 using RhsVector = typename ShapeMatricesType::template VectorType<
146 ShapeFunction::NPOINTS + ShapeFunction::NPOINTS * DisplacementDim>;
147 using JacobianMatrix = typename ShapeMatricesType::template MatrixType<
148 ShapeFunction::NPOINTS + ShapeFunction::NPOINTS * DisplacementDim,
149 ShapeFunction::NPOINTS + ShapeFunction::NPOINTS * DisplacementDim>;
150 using IpData =
152
154 delete;
156
158 MeshLib::Element const& e,
159 std::size_t const /*local_matrix_size*/,
160 NumLib::GenericIntegrationMethod const& integration_method,
161 bool const is_axially_symmetric,
163
165 std::size_t setIPDataInitialConditions(
166 std::string_view const name,
167 double const* values,
168 int const integration_order) override;
169
170 void assemble(double const /*t*/, double const /*dt*/,
171 std::vector<double> const& /*local_x*/,
172 std::vector<double> const& /*local_x_prev*/,
173 std::vector<double>& /*local_M_data*/,
174 std::vector<double>& /*local_K_data*/,
175 std::vector<double>& /*local_rhs_data*/) override
176 {
177 OGS_FATAL(
178 "ThermoMechanicsLocalAssembler: assembly without jacobian is not "
179 "implemented.");
180 }
181
183 double const t, double const dt, Eigen::VectorXd const& local_x,
184 Eigen::VectorXd const& local_x_prev, int const process_id,
185 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
186 std::vector<double>& local_b_data,
187 std::vector<double>& local_Jac_data) override;
188
189 void assembleWithJacobian(double const t, double const dt,
190 std::vector<double> const& local_x,
191 std::vector<double> const& local_x_prev,
192 std::vector<double>& local_M_data,
193 std::vector<double>& local_K_data,
194 std::vector<double>& local_rhs_data,
195 std::vector<double>& local_Jac_data) override;
196
197 void initializeConcrete() override
198 {
199 unsigned const n_integration_points =
201 for (unsigned ip = 0; ip < n_integration_points; ip++)
202 {
203 auto& ip_data = _ip_data[ip];
204
205 ParameterLib::SpatialPosition const x_position{
206 std::nullopt, _element.getID(), ip,
208 NumLib::interpolateCoordinates<ShapeFunction,
210 _element, ip_data.N))};
211
213 if (_process_data.initial_stress != nullptr)
214 {
215 ip_data.sigma =
217 DisplacementDim>((*_process_data.initial_stress)(
218 std::numeric_limits<
219 double>::quiet_NaN() /* time independent */,
220 x_position));
221 }
222
223 double const t = 0; // TODO (naumov) pass t from top
224 ip_data.solid_material.initializeInternalStateVariables(
225 t, x_position, *ip_data.material_state_variables);
226
227 ip_data.pushBackState();
228 }
229 }
230
231 void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
232 Eigen::VectorXd const& /*local_x_prev*/,
233 double const /*t*/, double const /*dt*/,
234 int const /*process_id*/) override
235 {
236 unsigned const n_integration_points =
238
239 for (unsigned ip = 0; ip < n_integration_points; ip++)
240 {
241 _ip_data[ip].pushBackState();
242 }
243 }
244
245 Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
246 const unsigned integration_point) const override
247 {
248 auto const& N = _secondary_data.N[integration_point];
249
250 // assumes N is stored contiguously in memory
251 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
252 }
253
254private:
278 const double t, double const dt, Eigen::VectorXd const& local_x,
279 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
280 std::vector<double>& local_Jac_data);
281
303 const double t, double const dt, Eigen::VectorXd const& local_x,
304 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
305 std::vector<double>& local_Jac_data);
306
307 std::size_t setSigma(double const* values);
308
309 // TODO (naumov) This method is same as getIntPtSigma but for arguments and
310 // the ordering of the cache_mat.
311 // There should be only one.
312 std::vector<double> getSigma() const override;
313
314 std::vector<double> const& getIntPtSigma(
315 const double t,
316 std::vector<GlobalVector*> const& x,
317 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
318 std::vector<double>& cache) const override;
319
320 std::size_t setEpsilon(double const* values);
321
322 std::vector<double> getEpsilon() const override;
323
324 std::vector<double> const& getIntPtEpsilon(
325 const double t,
326 std::vector<GlobalVector*> const& x,
327 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
328 std::vector<double>& cache) const override;
329
330 std::vector<double> const& getIntPtEpsilonMechanical(
331 const double t,
332 std::vector<GlobalVector*> const& x,
333 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
334 std::vector<double>& cache) const override;
335
336 std::size_t setEpsilonMechanical(double const* values);
337
338 std::vector<double> getEpsilonMechanical() const override;
339
340 unsigned getNumberOfIntegrationPoints() const override;
341
342 int getMaterialID() const override;
343
345 DisplacementDim>::MaterialStateVariables const&
346 getMaterialStateVariablesAt(unsigned integration_point) const override;
347
348private:
350
351 std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
352
357
358 static const int temperature_index = 0;
359 static const int temperature_size = ShapeFunction::NPOINTS;
360 static const int displacement_index = ShapeFunction::NPOINTS;
361 static const int displacement_size =
362 ShapeFunction::NPOINTS * DisplacementDim;
363};
364
365} // namespace ThermoMechanics
366} // namespace ProcessLib
367
#define OGS_FATAL(...)
Definition Error.h:26
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
VectorType< _number_of_dof > NodalForceVectorType
Rhs residual.
std::vector< double > const & getIntPtEpsilonMechanical(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
typename BMatricesType::NodalForceVectorType NodalForceVectorType
SecondaryData< typename ShapeMatrices::ShapeType > _secondary_data
std::size_t setIPDataInitialConditions(std::string_view const name, double const *values, int const integration_order) override
Returns number of read integration points.
typename ShapeMatricesType::template MatrixType< ShapeFunction::NPOINTS+ShapeFunction::NPOINTS *DisplacementDim, ShapeFunction::NPOINTS+ShapeFunction::NPOINTS *DisplacementDim > JacobianMatrix
std::vector< double > const & getIntPtSigma(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
void postTimestepConcrete(Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, int const) override
ThermoMechanicsLocalAssembler(ThermoMechanicsLocalAssembler const &)=delete
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt(unsigned integration_point) const override
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_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
typename ShapeMatricesType::template VectorType< ShapeFunction::NPOINTS+ShapeFunction::NPOINTS *DisplacementDim > RhsVector
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
ThermoMechanicsProcessData< DisplacementDim > & _process_data
NumLib::GenericIntegrationMethod const & _integration_method
void assembleWithJacobianForHeatConductionEquations(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(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 assemble(double const, double const, std::vector< double > const &, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &) override
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
ThermoMechanicsLocalAssembler(ThermoMechanicsLocalAssembler &&)=delete
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_M_data, std::vector< double > &local_K_data, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
std::vector< double > const & getIntPtEpsilon(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
NumLib::ShapeMatrices< NodalRowVectorType, DimNodalMatrixType, DimMatrixType, GlobalDimNodalMatrixType > ShapeMatrices
MatrixType< GlobalDim, GlobalDim > GlobalDimMatrixType
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material
ShapeMatricesType::NodalRowVectorType N
std::unique_ptr< typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables > material_state_variables
ShapeMatricesType::GlobalDimNodalMatrixType dNdx
IntegrationPointData(MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N