OGS 6.1.0-1721-g6382411ad
SmallDeformationLocalAssemblerMatrix-impl.h
Go to the documentation of this file.
1 
10 #pragma once
11 
13 
14 #include <valarray>
15 #include <vector>
16 
17 #include <Eigen/Eigen>
18 
22 
27 
29 
31 #include "SecondaryData.h"
33 
34 namespace ProcessLib
35 {
36 namespace LIE
37 {
38 namespace SmallDeformation
39 {
40 template <typename ShapeFunction, typename IntegrationMethod,
41  int DisplacementDim>
42 SmallDeformationLocalAssemblerMatrix<ShapeFunction, IntegrationMethod,
43  DisplacementDim>::
44  SmallDeformationLocalAssemblerMatrix(
45  MeshLib::Element const& e,
46  std::size_t const /*local_matrix_size*/,
47  bool const is_axially_symmetric,
48  unsigned const integration_order,
50  : _process_data(process_data),
51  _integration_method(integration_order),
52  _element(e),
53  _is_axially_symmetric(is_axially_symmetric)
54 {
55  unsigned const n_integration_points =
56  _integration_method.getNumberOfPoints();
57 
58  _ip_data.reserve(n_integration_points);
59  _secondary_data.N.resize(n_integration_points);
60 
61  auto const shape_matrices =
62  initShapeMatrices<ShapeFunction, ShapeMatricesType, IntegrationMethod,
63  DisplacementDim>(e, is_axially_symmetric,
65 
67  _process_data.solid_materials, _process_data.material_ids, e.getID());
68 
69  for (unsigned ip = 0; ip < n_integration_points; ip++)
70  {
71  _ip_data.emplace_back(solid_material);
72  auto& ip_data = _ip_data[ip];
73  auto const& sm = shape_matrices[ip];
74  ip_data.N = sm.N;
75  ip_data.dNdx = sm.dNdx;
76  ip_data.integration_weight =
77  _integration_method.getWeightedPoint(ip).getWeight() *
78  sm.integralMeasure * sm.detJ;
79 
80  // Initialize current time step values
81  static const int kelvin_vector_size =
83  DisplacementDim>::value;
84  ip_data._sigma.setZero(kelvin_vector_size);
85  ip_data._eps.setZero(kelvin_vector_size);
86 
87  // Previous time step values are not initialized and are set later.
88  ip_data._sigma_prev.resize(kelvin_vector_size);
89  ip_data._eps_prev.resize(kelvin_vector_size);
90 
91  ip_data._C.resize(kelvin_vector_size, kelvin_vector_size);
92 
93  _secondary_data.N[ip] = sm.N;
94  }
95 }
96 
97 template <typename ShapeFunction, typename IntegrationMethod,
98  int DisplacementDim>
99 void SmallDeformationLocalAssemblerMatrix<ShapeFunction, IntegrationMethod,
100  DisplacementDim>::
101  assembleWithJacobian(double const t, std::vector<double> const& local_x,
102  std::vector<double> const& /*local_xdot*/,
103  const double /*dxdot_dx*/, const double /*dx_dx*/,
104  std::vector<double>& /*local_M_data*/,
105  std::vector<double>& /*local_K_data*/,
106  std::vector<double>& local_b_data,
107  std::vector<double>& local_Jac_data)
108 {
109  assert(_element.getDimension() == DisplacementDim);
110 
111  auto const local_matrix_size = local_x.size();
112 
113  auto local_Jac = MathLib::createZeroedMatrix<StiffnessMatrixType>(
114  local_Jac_data, local_matrix_size, local_matrix_size);
115 
116  auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
117  local_b_data, local_matrix_size);
118 
119  unsigned const n_integration_points =
120  _integration_method.getNumberOfPoints();
121 
122  SpatialPosition x_position;
123  x_position.setElementID(_element.getID());
124 
125  for (unsigned ip = 0; ip < n_integration_points; ip++)
126  {
127  x_position.setIntegrationPoint(ip);
128  auto const& w = _ip_data[ip].integration_weight;
129 
130  auto const& N = _ip_data[ip].N;
131  auto const& dNdx = _ip_data[ip].dNdx;
132  auto const x_coord =
133  interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(_element,
134  N);
135  auto const B =
136  LinearBMatrix::computeBMatrix<DisplacementDim,
137  ShapeFunction::NPOINTS,
138  typename BMatricesType::BMatrixType>(
139  dNdx, N, x_coord, _is_axially_symmetric);
140 
141  auto const& eps_prev = _ip_data[ip]._eps_prev;
142  auto const& sigma_prev = _ip_data[ip]._sigma_prev;
143 
144  auto& eps = _ip_data[ip]._eps;
145  auto& sigma = _ip_data[ip]._sigma;
146  auto& state = _ip_data[ip]._material_state_variables;
147 
148  eps.noalias() =
149  B * Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
150  local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
151 
152  auto&& solution = _ip_data[ip]._solid_material.integrateStress(
153  t, x_position, _process_data.dt, eps_prev, eps, sigma_prev, *state,
154  _process_data._reference_temperature);
155 
156  if (!solution)
157  {
158  OGS_FATAL("Computation of local constitutive relation failed.");
159  }
160 
162  std::tie(sigma, state, C) = std::move(*solution);
163 
164  local_b.noalias() -= B.transpose() * sigma * w;
165  local_Jac.noalias() += B.transpose() * C * B * w;
166  }
167 }
168 
169 template <typename ShapeFunction, typename IntegrationMethod,
170  int DisplacementDim>
171 void SmallDeformationLocalAssemblerMatrix<ShapeFunction, IntegrationMethod,
172  DisplacementDim>::
174  double const /*t*/, Eigen::VectorXd const& /*local_x*/)
175 {
176  // Compute average value per element
177  const int n = DisplacementDim == 2 ? 4 : 6;
178  Eigen::VectorXd ele_stress = Eigen::VectorXd::Zero(n);
179  Eigen::VectorXd ele_strain = Eigen::VectorXd::Zero(n);
180 
181  unsigned const n_integration_points =
182  _integration_method.getNumberOfPoints();
183  for (unsigned ip = 0; ip < n_integration_points; ip++)
184  {
185  auto& ip_data = _ip_data[ip];
186 
187  ele_stress += ip_data._sigma;
188  ele_strain += ip_data._eps;
189  }
190  ele_stress /= n_integration_points;
191  ele_strain /= n_integration_points;
192 
193  (*_process_data._mesh_prop_stress_xx)[_element.getID()] = ele_stress[0];
194  (*_process_data._mesh_prop_stress_yy)[_element.getID()] = ele_stress[1];
195  (*_process_data._mesh_prop_stress_zz)[_element.getID()] = ele_stress[2];
196  (*_process_data._mesh_prop_stress_xy)[_element.getID()] = ele_stress[3];
197  if (DisplacementDim == 3)
198  {
199  (*_process_data._mesh_prop_stress_yz)[_element.getID()] = ele_stress[4];
200  (*_process_data._mesh_prop_stress_xz)[_element.getID()] = ele_stress[5];
201  }
202 
203  (*_process_data._mesh_prop_strain_xx)[_element.getID()] = ele_strain[0];
204  (*_process_data._mesh_prop_strain_yy)[_element.getID()] = ele_strain[1];
205  (*_process_data._mesh_prop_strain_zz)[_element.getID()] = ele_strain[2];
206  (*_process_data._mesh_prop_strain_xy)[_element.getID()] = ele_strain[3];
207  if (DisplacementDim == 3)
208  {
209  (*_process_data._mesh_prop_strain_yz)[_element.getID()] = ele_strain[4];
210  (*_process_data._mesh_prop_strain_xz)[_element.getID()] = ele_strain[5];
211  }
212 }
213 
214 } // namespace SmallDeformation
215 } // namespace LIE
216 } // namespace ProcessLib
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:23
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
Definition: SecondaryData.h:29
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool is_axially_symmetric, IntegrationMethod const &integration_method)
MechanicsBase< DisplacementDim > & selectSolidConstitutiveRelation(std::map< int, std::unique_ptr< MechanicsBase< DisplacementDim >>> const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
void computeSecondaryVariableConcreteWithVector(double const t, Eigen::VectorXd const &local_x) override
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, KelvinVectorDimensions< DisplacementDim >::value, Eigen::RowMajor > KelvinMatrixType
Definition: KelvinVector.h:59
void setElementID(std::size_t element_id)
virtual unsigned getDimension() const =0
Get dimension of the mesh element.
std::vector< IntegrationPointDataMatrix< ShapeMatricesType, BMatricesType, DisplacementDim >, Eigen::aligned_allocator< IntegrationPointDataMatrix< ShapeMatricesType, BMatricesType, DisplacementDim > > > _ip_data
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:41
void assembleWithJacobian(double const t, std::vector< double > const &local_x, std::vector< double > const &, const double, const double, std::vector< double > &, std::vector< double > &, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:90
#define OGS_FATAL(fmt,...)
Definition: Error.h:71
void setIntegrationPoint(unsigned integration_point)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:54