OGS
SmallDeformationLocalAssemblerMatrix-impl.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <Eigen/Core>
14#include <valarray>
15#include <vector>
16
26#include "SecondaryData.h"
29
30namespace ProcessLib
31{
32namespace LIE
33{
34namespace SmallDeformation
35{
36template <typename ShapeFunction, int DisplacementDim>
39 MeshLib::Element const& e,
40 std::size_t const /*local_matrix_size*/,
41 NumLib::GenericIntegrationMethod const& integration_method,
42 bool const is_axially_symmetric,
44 : _process_data(process_data),
45 _integration_method(integration_method),
46 _element(e),
47 _is_axially_symmetric(is_axially_symmetric)
48{
49 unsigned const n_integration_points =
51
52 _ip_data.reserve(n_integration_points);
53 _secondary_data.N.resize(n_integration_points);
54
55 auto const shape_matrices =
57 DisplacementDim>(e, is_axially_symmetric,
59
61 _process_data.solid_materials, _process_data.material_ids, e.getID());
62
63 for (unsigned ip = 0; ip < n_integration_points; ip++)
64 {
65 _ip_data.emplace_back(solid_material);
66 auto& ip_data = _ip_data[ip];
67 auto const& sm = shape_matrices[ip];
68 ip_data.N = sm.N;
69 ip_data.dNdx = sm.dNdx;
70 ip_data.integration_weight =
72 sm.integralMeasure * sm.detJ;
73
74 // Initialize current time step values
75 static const int kelvin_vector_size =
77 ip_data._sigma.setZero(kelvin_vector_size);
78 ip_data._eps.setZero(kelvin_vector_size);
79
80 // Previous time step values are not initialized and are set later.
81 ip_data._sigma_prev.resize(kelvin_vector_size);
82 ip_data._eps_prev.resize(kelvin_vector_size);
83
84 ip_data._C.resize(kelvin_vector_size, kelvin_vector_size);
85
86 _secondary_data.N[ip] = sm.N;
87 }
88}
89
90template <typename ShapeFunction, int DisplacementDim>
92 assembleWithJacobian(double const t, double const dt,
93 std::vector<double> const& local_x,
94 std::vector<double> const& /*local_x_prev*/,
95 std::vector<double>& /*local_M_data*/,
96 std::vector<double>& /*local_K_data*/,
97 std::vector<double>& local_b_data,
98 std::vector<double>& local_Jac_data)
99{
100 assert(_element.getDimension() == DisplacementDim);
101
102 auto const local_matrix_size = local_x.size();
103
104 auto local_Jac = MathLib::createZeroedMatrix<StiffnessMatrixType>(
105 local_Jac_data, local_matrix_size, local_matrix_size);
106
107 auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
108 local_b_data, local_matrix_size);
109
110 unsigned const n_integration_points =
111 _integration_method.getNumberOfPoints();
112
113 MPL::VariableArray variables;
114 MPL::VariableArray variables_prev;
116 x_position.setElementID(_element.getID());
117
118 for (unsigned ip = 0; ip < n_integration_points; ip++)
119 {
120 x_position.setIntegrationPoint(ip);
121 auto const& w = _ip_data[ip].integration_weight;
122
123 auto const& N = _ip_data[ip].N;
124 auto const& dNdx = _ip_data[ip].dNdx;
125 auto const x_coord =
126 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
127 _element, N);
128 auto const B =
129 LinearBMatrix::computeBMatrix<DisplacementDim,
130 ShapeFunction::NPOINTS,
132 dNdx, N, x_coord, _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
178 const int n = DisplacementDim == 2 ? 4 : 6;
179 Eigen::VectorXd ele_stress = Eigen::VectorXd::Zero(n);
180 Eigen::VectorXd ele_strain = Eigen::VectorXd::Zero(n);
181
182 unsigned const n_integration_points =
183 _integration_method.getNumberOfPoints();
184 for (unsigned ip = 0; ip < n_integration_points; ip++)
185 {
186 auto const& ip_data = _ip_data[ip];
187
188 ele_stress += ip_data._sigma;
189 ele_strain += ip_data._eps;
190 }
191 ele_stress /= n_integration_points;
192 ele_strain /= n_integration_points;
193
194 (*_process_data._mesh_prop_stress_xx)[_element.getID()] = ele_stress[0];
195 (*_process_data._mesh_prop_stress_yy)[_element.getID()] = ele_stress[1];
196 (*_process_data._mesh_prop_stress_zz)[_element.getID()] = ele_stress[2];
197 (*_process_data._mesh_prop_stress_xy)[_element.getID()] = ele_stress[3];
198 if (DisplacementDim == 3)
199 {
200 (*_process_data._mesh_prop_stress_yz)[_element.getID()] = ele_stress[4];
201 (*_process_data._mesh_prop_stress_xz)[_element.getID()] = ele_stress[5];
202 }
203
204 (*_process_data._mesh_prop_strain_xx)[_element.getID()] = ele_strain[0];
205 (*_process_data._mesh_prop_strain_yy)[_element.getID()] = ele_strain[1];
206 (*_process_data._mesh_prop_strain_zz)[_element.getID()] = ele_strain[2];
207 (*_process_data._mesh_prop_strain_xy)[_element.getID()] = ele_strain[3];
208 if (DisplacementDim == 3)
209 {
210 (*_process_data._mesh_prop_strain_yz)[_element.getID()] = ele_strain[4];
211 (*_process_data._mesh_prop_strain_xz)[_element.getID()] = ele_strain[5];
212 }
213}
214
215} // namespace SmallDeformation
216} // namespace LIE
217} // namespace ProcessLib
#define OGS_FATAL(...)
Definition: Error.h:26
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > mechanical_strain
Definition: VariableType.h:182
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > stress
Definition: VariableType.h:192
double getWeight() const
Definition: WeightedPoint.h:80
virtual std::size_t getID() const final
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
Definition: BMatrixPolicy.h:55
void computeSecondaryVariableConcreteWithVector(double const t, Eigen::VectorXd const &local_x) override
void assembleWithJacobian(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &, std::vector< double > &, std::vector< double > &, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
std::vector< IntegrationPointDataMatrix< ShapeMatricesType, BMatricesType, DisplacementDim >, Eigen::aligned_allocator< IntegrationPointDataMatrix< ShapeMatricesType, BMatricesType, DisplacementDim > > > _ip_data
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, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:24
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
Definition: KelvinVector.h:57
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 computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
Definition: LinearBMatrix.h:42
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
Definition: SecondaryData.h:27