OGS
LocalAssemblerInterface.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <Eigen/Eigenvalues>
14#include <vector>
15
19#include "MaterialForces.h"
28
29namespace ProcessLib
30{
31namespace SmallDeformation
32{
33template <int DisplacementDim>
38{
40 MeshLib::Element const& e,
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 solid_material_(MaterialLib::Solids::selectSolidConstitutiveRelation(
49 process_data_.solid_materials, process_data_.material_ids,
50 element_.getID()))
51 {
52 unsigned const n_integration_points =
54
55 current_states_.resize(n_integration_points);
56 prev_states_.resize(n_integration_points);
57 output_data_.resize(n_integration_points);
58
59 material_states_.reserve(n_integration_points);
60 for (unsigned ip = 0; ip < n_integration_points; ++ip)
61 {
62 material_states_.emplace_back(
63 solid_material_.createMaterialStateVariables());
64
65 // Set initial stress field to zero. Might be overwritten by
66 // integration point data or initial stress.
67 this->current_states_[ip].stress_data.sigma.noalias() =
69 DisplacementDim>::Zero();
70 }
71 }
73 std::size_t setIPDataInitialConditions(std::string_view name,
74 double const* values,
75 int const integration_order)
76 {
77 if (integration_order !=
78 static_cast<int>(integration_method_.getIntegrationOrder()))
79 {
81 "Setting integration point initial conditions; The integration "
82 "order of the local assembler for element {:d} is different "
83 "from the integration order in the initial condition.",
84 element_.getID());
85 }
86
87 if (name.starts_with("material_state_variable_"))
88 {
89 name.remove_prefix(24);
90
91 auto const& internal_variables =
92 solid_material_.getInternalVariables();
93 if (auto const iv = std::find_if(
94 begin(internal_variables), end(internal_variables),
95 [&name](auto const& iv) { return iv.name == name; });
96 iv != end(internal_variables))
97 {
98 DBUG("Setting material state variable '{:s}'", name);
99 return ProcessLib::
101 values, material_states_,
103 DisplacementDim>::material_state_variables,
104 iv->reference);
105 }
106
107 WARN(
108 "Could not find variable {:s} in solid material model's "
109 "internal variables.",
110 name);
111 return 0;
112 }
113
114 // TODO this logic could be pulled out of the local assembler into the
115 // process. That might lead to a slightly better performance due to less
116 // string comparisons.
117 return ProcessLib::Reflection::reflectSetIPData<DisplacementDim>(
118 name, values, current_states_);
119 }
120
122 double const /*t*/, double const /*dt*/, Eigen::VectorXd const& /*x*/,
123 Eigen::VectorXd const& /*x_prev*/) override
124 {
125 int const elem_id = element_.getID();
127 x_position.setElementID(elem_id);
128 unsigned const n_integration_points =
130
131 auto sigma_sum = MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
132 Eigen::Matrix<double, 3, 3>::Zero());
133
134 for (unsigned ip = 0; ip < n_integration_points; ip++)
135 {
136 x_position.setIntegrationPoint(ip);
137 auto const& sigma = current_states_[ip].stress_data.sigma;
138 sigma_sum += sigma;
139 }
140
141 Eigen::Matrix<double, 3, 3, 0, 3, 3> const sigma_avg =
143 n_integration_points;
144
145 Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(
146 sigma_avg);
147
148 Eigen::Map<Eigen::Vector3d>(
149 &(*process_data_.principal_stress_values)[elem_id * 3], 3) =
150 e_s.eigenvalues();
151
152 auto eigen_vectors = e_s.eigenvectors();
153
154 for (auto i = 0; i < 3; i++)
155 {
156 Eigen::Map<Eigen::Vector3d>(
157 &(*process_data_.principal_stress_vector[i])[elem_id * 3], 3) =
158 eigen_vectors.col(i);
159 }
160 }
161
162 // TODO move to NumLib::ExtrapolatableElement
164 {
166 }
167
168 int getMaterialID() const
169 {
170 return process_data_.material_ids == nullptr
171 ? 0
172 : (*process_data_.material_ids)[element_.getID()];
173 }
174
176 std::function<std::span<double>(
178 MaterialStateVariables&)> const& get_values_span,
179 int const& n_components) const
180 {
184 get_values_span, n_components);
185 }
186
188 DisplacementDim>::MaterialStateVariables const&
189 getMaterialStateVariablesAt(unsigned integration_point) const
190 {
191 return *material_states_[integration_point].material_state_variables;
192 }
193
195 {
197
199 &Self::current_states_, &Self::output_data_);
200 }
201
202protected:
204
205 // Material state is special, because it contains both the current and the
206 // old state.
207 std::vector<MaterialStateData<DisplacementDim>> material_states_;
208 std::vector<typename ConstitutiveRelations::StatefulData<DisplacementDim>>
209 current_states_; // TODO maybe do not store but rather re-evaluate for
210 // state update
211 std::vector<
214 std::vector<typename ConstitutiveRelations::OutputData<DisplacementDim>>
216
223};
224
225} // namespace SmallDeformation
226} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Eigen::Matrix< double, 3, 3 > kelvinVectorToTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
auto reflectWithoutName(Accessors &&... accessors)
std::size_t setIntegrationPointDataMaterialStateVariables(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType member, std::function< std::span< double >(MaterialStateVariables &)> get_values_span)
std::vector< double > getIntegrationPointDataMaterialStateVariables(IntegrationPointDataVector const &ip_data_vector, MemberType member, std::function< std::span< double >(MaterialStateVariables &)> get_values_span, int const n_components)
std::vector< double > getMaterialStateVariableInternalState(std::function< std::span< double >(typename MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables &)> const &get_values_span, int const &n_components) const
ConstitutiveRelations::ConstitutiveSetting< DisplacementDim > constitutive_setting
std::vector< typename ConstitutiveRelations::StatefulDataPrev< DisplacementDim > > prev_states_
std::size_t setIPDataInitialConditions(std::string_view name, double const *values, int const integration_order)
Returns number of read integration points.
std::vector< typename ConstitutiveRelations::StatefulData< DisplacementDim > > current_states_
std::vector< typename ConstitutiveRelations::OutputData< DisplacementDim > > output_data_
SmallDeformationLocalAssemblerInterface(MeshLib::Element const &e, NumLib::GenericIntegrationMethod const &integration_method, bool const is_axially_symmetric, SmallDeformationProcessData< DisplacementDim > &process_data)
void computeSecondaryVariableConcrete(double const, double const, Eigen::VectorXd const &, Eigen::VectorXd const &) override
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt(unsigned integration_point) const
MaterialLib::Solids::MechanicsBase< DisplacementDim > const & solid_material_
std::vector< MaterialStateData< DisplacementDim > > material_states_