34namespace SmallDeformation
36template <
typename ShapeFunction,
int DisplacementDim>
42 bool const is_axially_symmetric,
44 : _process_data(process_data),
45 _integration_method(integration_method),
47 _is_axially_symmetric(is_axially_symmetric)
49 unsigned const n_integration_points =
52 _ip_data.reserve(n_integration_points);
55 auto const shape_matrices =
57 DisplacementDim>(e, is_axially_symmetric,
63 for (
unsigned ip = 0; ip < n_integration_points; ip++)
65 _ip_data.emplace_back(solid_material);
67 auto const& sm = shape_matrices[ip];
69 ip_data.dNdx = sm.dNdx;
70 ip_data.integration_weight =
72 sm.integralMeasure * sm.detJ;
75 static const int kelvin_vector_size =
77 ip_data._sigma.setZero(kelvin_vector_size);
78 ip_data._eps.setZero(kelvin_vector_size);
81 ip_data._sigma_prev.resize(kelvin_vector_size);
82 ip_data._eps_prev.resize(kelvin_vector_size);
84 ip_data._C.resize(kelvin_vector_size, kelvin_vector_size);
90template <
typename ShapeFunction,
int DisplacementDim>
93 std::vector<double>
const& local_x,
94 std::vector<double>
const& ,
95 std::vector<double>& ,
96 std::vector<double>& ,
97 std::vector<double>& local_b_data,
98 std::vector<double>& local_Jac_data)
100 assert(_element.getDimension() == DisplacementDim);
102 auto const local_matrix_size = local_x.size();
104 auto local_Jac = MathLib::createZeroedMatrix<StiffnessMatrixType>(
105 local_Jac_data, local_matrix_size, local_matrix_size);
107 auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
108 local_b_data, local_matrix_size);
110 unsigned const n_integration_points =
111 _integration_method.getNumberOfPoints();
118 for (
unsigned ip = 0; ip < n_integration_points; ip++)
121 auto const& w = _ip_data[ip].integration_weight;
123 auto const& N = _ip_data[ip].N;
124 auto const& dNdx = _ip_data[ip].dNdx;
126 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
130 ShapeFunction::NPOINTS,
132 dNdx, N, x_coord, _is_axially_symmetric);
134 auto const& eps_prev = _ip_data[ip]._eps_prev;
135 auto const& sigma_prev = _ip_data[ip]._sigma_prev;
137 auto& eps = _ip_data[ip]._eps;
138 auto& sigma = _ip_data[ip]._sigma;
139 auto& state = _ip_data[ip]._material_state_variables;
142 B * Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
143 local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
154 variables_prev.
temperature = _process_data._reference_temperature;
156 auto&& solution = _ip_data[ip]._solid_material.integrateStress(
157 variables_prev, variables, t, x_position, dt, *state);
161 OGS_FATAL(
"Computation of local constitutive relation failed.");
165 std::tie(sigma, state, C) = std::move(*solution);
167 local_b.noalias() -= B.transpose() * sigma * w;
168 local_Jac.noalias() += B.transpose() * C * B * w;
172template <
typename ShapeFunction,
int DisplacementDim>
175 double const , Eigen::VectorXd
const& )
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);
182 unsigned const n_integration_points =
183 _integration_method.getNumberOfPoints();
184 for (
unsigned ip = 0; ip < n_integration_points; ip++)
186 auto const& ip_data = _ip_data[ip];
188 ele_stress += ip_data._sigma;
189 ele_strain += ip_data._eps;
191 ele_stress /= n_integration_points;
192 ele_strain /= n_integration_points;
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)
200 (*_process_data._mesh_prop_stress_yz)[_element.getID()] = ele_stress[4];
201 (*_process_data._mesh_prop_stress_xz)[_element.getID()] = ele_stress[5];
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)
210 (*_process_data._mesh_prop_strain_yz)[_element.getID()] = ele_strain[4];
211 (*_process_data._mesh_prop_strain_xz)[_element.getID()] = ele_strain[5];
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > mechanical_strain
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > stress
virtual std::size_t getID() const final
Returns the ID of the element.
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
unsigned getNumberOfPoints() const
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
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
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
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.
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N