OGS
SmallDeformationLocalAssemblerMatrix-impl.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <Eigen/Core>
14#include <vector>
15
27#include "SecondaryData.h"
30
31namespace ProcessLib
32{
33namespace LIE
34{
35namespace SmallDeformation
36{
37template <typename ShapeFunction, int DisplacementDim>
40 MeshLib::Element const& e,
41 std::size_t const /*local_matrix_size*/,
42 NumLib::GenericIntegrationMethod const& integration_method,
43 bool const is_axially_symmetric,
45 : _process_data(process_data),
46 _integration_method(integration_method),
47 _element(e),
48 _is_axially_symmetric(is_axially_symmetric)
49{
50 unsigned const n_integration_points =
52
53 _ip_data.reserve(n_integration_points);
54 _secondary_data.N.resize(n_integration_points);
55
56 auto const shape_matrices =
58 DisplacementDim>(e, is_axially_symmetric,
60
62 _process_data.solid_materials, _process_data.material_ids, e.getID());
63
64 for (unsigned ip = 0; ip < n_integration_points; ip++)
65 {
66 _ip_data.emplace_back(solid_material);
67 auto& ip_data = _ip_data[ip];
68 auto const& sm = shape_matrices[ip];
69 ip_data.N_u = sm.N;
70 ip_data.dNdx_u = sm.dNdx;
71 ip_data.integration_weight =
73 sm.integralMeasure * sm.detJ;
74
75 // Initialize current time step values
76 static const int kelvin_vector_size =
78 ip_data._sigma.setZero(kelvin_vector_size);
79 ip_data._eps.setZero(kelvin_vector_size);
80
81 // Previous time step values are not initialized and are set later.
82 ip_data._sigma_prev.resize(kelvin_vector_size);
83 ip_data._eps_prev.resize(kelvin_vector_size);
84
85 ip_data._C.resize(kelvin_vector_size, kelvin_vector_size);
86
87 _secondary_data.N[ip] = sm.N;
88 }
89}
90
91template <typename ShapeFunction, int DisplacementDim>
93 assembleWithJacobian(double const t, double const dt,
94 std::vector<double> const& local_x,
95 std::vector<double> const& /*local_x_prev*/,
96 std::vector<double>& local_b_data,
97 std::vector<double>& local_Jac_data)
98{
99 assert(_element.getDimension() == DisplacementDim);
100
101 auto const local_matrix_size = local_x.size();
102
104 local_Jac_data, local_matrix_size, local_matrix_size);
105
107 local_b_data, local_matrix_size);
108
109 unsigned const n_integration_points =
110 _integration_method.getNumberOfPoints();
111
112 MPL::VariableArray variables;
113 MPL::VariableArray variables_prev;
115 x_position.setElementID(_element.getID());
116
117 auto const B_dil_bar = getDilatationalBBarMatrix();
118
119 for (unsigned ip = 0; ip < n_integration_points; ip++)
120 {
121 auto const& w = _ip_data[ip].integration_weight;
122
123 auto const& N = _ip_data[ip].N_u;
124 auto const& dNdx = _ip_data[ip].dNdx_u;
125 auto const x_coord =
127 _element, N);
128
130 DisplacementDim, ShapeFunction::NPOINTS, BBarMatrixType,
131 typename BMatricesType::BMatrixType>(dNdx, N, B_dil_bar, x_coord,
132 this->_is_axially_symmetric);
133
134 auto const& eps_prev = _ip_data[ip]._eps_prev;
135 auto const& sigma_prev = _ip_data[ip]._sigma_prev;
136
137 auto& eps = _ip_data[ip]._eps;
138 auto& sigma = _ip_data[ip]._sigma;
139 auto& state = _ip_data[ip]._material_state_variables;
140
141 eps.noalias() =
142 B * Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
143 local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
144
145 variables.mechanical_strain
147 eps);
148 variables_prev.stress
150 sigma_prev);
151 variables_prev.mechanical_strain
153 eps_prev);
154 variables_prev.temperature = _process_data.reference_temperature;
155
156 auto&& solution = _ip_data[ip]._solid_material.integrateStress(
157 variables_prev, variables, t, x_position, dt, *state);
158
159 if (!solution)
160 {
161 OGS_FATAL("Computation of local constitutive relation failed.");
162 }
163
165 std::tie(sigma, state, C) = std::move(*solution);
166
167 local_b.noalias() -= B.transpose() * sigma * w;
168 local_Jac.noalias() += B.transpose() * C * B * w;
169 }
170}
171
172template <typename ShapeFunction, int DisplacementDim>
175 double const /*t*/, Eigen::VectorXd const& /*local_x*/)
176{
177 // Compute average value per element
179 KV sigma_avg = KV::Zero();
180 auto const e_id = _element.getID();
181
182 unsigned const n_integration_points =
183 _integration_method.getNumberOfPoints();
184 for (unsigned ip = 0; ip < n_integration_points; ip++)
185 {
186 sigma_avg += _ip_data[ip]._sigma;
187 }
188 sigma_avg /= n_integration_points;
189
190 Eigen::Map<KV>(
191 &(*_process_data.element_stresses)[e_id * KV::RowsAtCompileTime]) =
193}
194
195template <typename ShapeFunction, int DisplacementDim>
196std::vector<double> const&
199 const double /*t*/,
200 std::vector<GlobalVector*> const& /*x*/,
201 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
202 std::vector<double>& cache) const
203{
205 _ip_data, &IntegrationPointDataType::_sigma, cache);
206}
207template <typename ShapeFunction, int DisplacementDim>
208std::vector<double> const&
211 const double /*t*/,
212 std::vector<GlobalVector*> const& /*x*/,
213 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
214 std::vector<double>& cache) const
215{
217 _ip_data, &IntegrationPointDataType::_eps, cache);
218}
219
220} // namespace SmallDeformation
221} // namespace LIE
222} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
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)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
void computeSecondaryVariableConcreteWithVector(double const t, Eigen::VectorXd const &local_x) override
std::vector< IntegrationPointDataType, Eigen::aligned_allocator< IntegrationPointDataType > > _ip_data
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 assembleWithJacobian(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
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
SmallDeformationLocalAssemblerMatrix(SmallDeformationLocalAssemblerMatrix const &)=delete
auto & selectSolidConstitutiveRelation(SolidMaterialsMap const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
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)
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 & getIntegrationPointKelvinVectorData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N