OGS
ThermoMechanicsFEM-impl.h
Go to the documentation of this file.
1
13#pragma once
14
21
22namespace ProcessLib
23{
24namespace ThermoMechanics
25{
26template <typename ShapeFunction, int DisplacementDim>
29 MeshLib::Element const& e,
30 std::size_t const /*local_matrix_size*/,
31 NumLib::GenericIntegrationMethod const& integration_method,
32 bool const is_axially_symmetric,
34 : _process_data(process_data),
35 _integration_method(integration_method),
36 _element(e),
37 _is_axially_symmetric(is_axially_symmetric)
38{
39 unsigned const n_integration_points =
41
42 _ip_data.reserve(n_integration_points);
43 _secondary_data.N.resize(n_integration_points);
44
45 auto const shape_matrices =
47 DisplacementDim>(e, is_axially_symmetric,
49
51 _process_data.solid_materials, _process_data.material_ids, e.getID());
52
53 for (unsigned ip = 0; ip < n_integration_points; ip++)
54 {
55 _ip_data.emplace_back(solid_material);
56 auto& ip_data = _ip_data[ip];
57 ip_data.integration_weight =
59 shape_matrices[ip].integralMeasure * shape_matrices[ip].detJ;
60
61 static const int kelvin_vector_size =
63 // Initialize current time step values
64 ip_data.sigma.setZero(kelvin_vector_size);
65 ip_data.eps.setZero(kelvin_vector_size);
66
67 // Previous time step values are not initialized and are set later.
68 ip_data.sigma_prev.resize(kelvin_vector_size);
69 ip_data.eps_prev.resize(kelvin_vector_size);
70
71 ip_data.eps_m.setZero(kelvin_vector_size);
72 ip_data.eps_m_prev.setZero(kelvin_vector_size);
74 x_position.setElementID(_element.getID());
75 ip_data.N = shape_matrices[ip].N;
76 ip_data.dNdx = shape_matrices[ip].dNdx;
77
78 _secondary_data.N[ip] = shape_matrices[ip].N;
79 }
80}
81
82template <typename ShapeFunction, int DisplacementDim>
84 setIPDataInitialConditions(std::string_view const name,
85 double const* values,
86 int const integration_order)
87{
88 if (integration_order !=
89 static_cast<int>(_integration_method.getIntegrationOrder()))
90 {
92 "Setting integration point initial conditions; The integration "
93 "order of the local assembler for element {:d} is different from "
94 "the integration order in the initial condition.",
95 _element.getID());
96 }
97
98 if (name == "sigma")
99 {
100 return setSigma(values);
101 }
102 if (name == "epsilon")
103 {
104 return setEpsilon(values);
105 }
106 if (name == "epsilon_m")
107 {
108 return setEpsilonMechanical(values);
109 }
110
111 return 0;
112}
113
114template <typename ShapeFunction, int DisplacementDim>
116 assembleWithJacobian(double const t, double const dt,
117 std::vector<double> const& local_x,
118 std::vector<double> const& local_x_prev,
119 std::vector<double>& local_rhs_data,
120 std::vector<double>& local_Jac_data)
121{
122 auto const local_matrix_size = local_x.size();
123 assert(local_matrix_size == temperature_size + displacement_size);
124
125 auto T = Eigen::Map<typename ShapeMatricesType::template VectorType<
126 temperature_size> const>(local_x.data() + temperature_index,
127 temperature_size);
128
129 auto u = Eigen::Map<typename ShapeMatricesType::template VectorType<
130 displacement_size> const>(local_x.data() + displacement_index,
131 displacement_size);
132 bool const is_u_non_zero = u.norm() > 0.0;
133
134 auto T_prev = Eigen::Map<typename ShapeMatricesType::template VectorType<
135 temperature_size> const>(local_x_prev.data() + temperature_index,
136 temperature_size);
137
138 auto local_Jac = MathLib::createZeroedMatrix<JacobianMatrix>(
139 local_Jac_data, local_matrix_size, local_matrix_size);
140
141 auto local_rhs = MathLib::createZeroedVector<RhsVector>(local_rhs_data,
142 local_matrix_size);
143
144 typename ShapeMatricesType::template MatrixType<displacement_size,
145 temperature_size>
146 KuT;
147 KuT.setZero(displacement_size, temperature_size);
148
150 KTT.setZero(temperature_size, temperature_size);
151
153 DTT.setZero(temperature_size, temperature_size);
154
155 unsigned const n_integration_points =
156 _integration_method.getNumberOfPoints();
157
158 MPL::VariableArray variables;
159 MPL::VariableArray variables_prev;
161 x_position.setElementID(_element.getID());
162
163 auto const& medium = _process_data.media_map.getMedium(_element.getID());
164 auto const& solid_phase = medium->phase("Solid");
165
166 for (unsigned ip = 0; ip < n_integration_points; ip++)
167 {
168 x_position.setIntegrationPoint(ip);
169 auto const& w = _ip_data[ip].integration_weight;
170 auto const& N = _ip_data[ip].N;
171 auto const& dNdx = _ip_data[ip].dNdx;
172
173 auto const x_coord =
174 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
175 _element, N);
176 auto const& B =
177 LinearBMatrix::computeBMatrix<DisplacementDim,
178 ShapeFunction::NPOINTS,
180 dNdx, N, x_coord, _is_axially_symmetric);
181
182 auto& sigma = _ip_data[ip].sigma;
183 auto const& sigma_prev = _ip_data[ip].sigma_prev;
184
185 auto& eps = _ip_data[ip].eps;
186 auto const& eps_prev = _ip_data[ip].eps_prev;
187
188 auto& eps_m = _ip_data[ip].eps_m;
189 auto const& eps_m_prev = _ip_data[ip].eps_m_prev;
190
191 auto& state = _ip_data[ip].material_state_variables;
192
193 const double T_ip = N.dot(T); // T at integration point
194 double const T_prev_ip = N.dot(T_prev);
195
196 // Consider also anisotropic thermal expansion.
197 auto const solid_linear_thermal_expansivity_vector =
198 MPL::formKelvinVector<DisplacementDim>(
199 solid_phase
200 .property(
202 .value(variables, x_position, t, dt));
203
205 dthermal_strain =
206 solid_linear_thermal_expansivity_vector * (T_ip - T_prev_ip);
207
208 //
209 // displacement equation, displacement part
210 //
211 // For the restart computation, the displacement may not be
212 // reloaded but the initial strains are always available. For such case,
213 // the following computation is skipped.
214 if (is_u_non_zero)
215 {
216 eps.noalias() = B * u;
217 }
218
219 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
220
221 variables_prev.stress
223 sigma_prev);
224 variables_prev.mechanical_strain
226 eps_m_prev);
227 variables_prev.temperature = T_ip;
228 variables.mechanical_strain
230 eps_m);
231 variables.temperature = T_ip;
232
233 auto&& solution = _ip_data[ip].solid_material.integrateStress(
234 variables_prev, variables, t, x_position, dt, *state);
235
236 if (!solution)
237 {
238 OGS_FATAL("Computation of local constitutive relation failed.");
239 }
240
242 std::tie(sigma, state, C) = std::move(*solution);
243
244 local_Jac
245 .template block<displacement_size, displacement_size>(
246 displacement_index, displacement_index)
247 .noalias() += B.transpose() * C * B * w;
248
249 typename ShapeMatricesType::template MatrixType<DisplacementDim,
250 displacement_size>
251 N_u = ShapeMatricesType::template MatrixType<
252 DisplacementDim, displacement_size>::Zero(DisplacementDim,
253 displacement_size);
254
255 for (int i = 0; i < DisplacementDim; ++i)
256 {
257 N_u.template block<1, displacement_size / DisplacementDim>(
258 i, i * displacement_size / DisplacementDim)
259 .noalias() = N;
260 }
261
262 auto const rho_s =
263 solid_phase.property(MPL::PropertyType::density)
264 .template value<double>(variables, x_position, t, dt);
265
266 auto const& b = _process_data.specific_body_force;
267 local_rhs.template block<displacement_size, 1>(displacement_index, 0)
268 .noalias() -=
269 (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
270
271 //
272 // displacement equation, temperature part
273 // The computation of KuT can be ignored.
274 auto const alpha_T_tensor =
276 solid_linear_thermal_expansivity_vector);
277 KuT.noalias() += B.transpose() * (C * alpha_T_tensor) * N * w;
278
279 if (_ip_data[ip].solid_material.getConstitutiveModel() ==
281 {
282 auto const s = Invariants::deviatoric_projection * sigma;
283 double const norm_s = Invariants::FrobeniusNorm(s);
284 const double creep_coefficient =
285 _ip_data[ip].solid_material.getTemperatureRelatedCoefficient(
286 t, dt, x_position, T_ip, norm_s);
287 KuT.noalias() += creep_coefficient * B.transpose() * s * N * w;
288 }
289
290 //
291 // temperature equation, temperature part;
292 //
293 auto const lambda =
294 solid_phase
295 .property(
297 .value(variables, x_position, t, dt);
298
299 GlobalDimMatrixType const thermal_conductivity =
300 MaterialPropertyLib::formEigenTensor<DisplacementDim>(lambda);
301
302 KTT.noalias() += dNdx.transpose() * thermal_conductivity * dNdx * w;
303
304 auto const c =
305 solid_phase
306 .property(
308 .template value<double>(variables, x_position, t, dt);
309 DTT.noalias() += N.transpose() * rho_s * c * N * w;
310 }
311
312 // temperature equation, temperature part
313 local_Jac
314 .template block<temperature_size, temperature_size>(temperature_index,
315 temperature_index)
316 .noalias() += KTT + DTT / dt;
317
318 // displacement equation, temperature part
319 local_Jac
320 .template block<displacement_size, temperature_size>(displacement_index,
321 temperature_index)
322 .noalias() -= KuT;
323
324 local_rhs.template block<temperature_size, 1>(temperature_index, 0)
325 .noalias() -= KTT * T + DTT * (T - T_prev) / dt;
326}
327
328template <typename ShapeFunction, int DisplacementDim>
330 assembleWithJacobianForStaggeredScheme(const double t, double const dt,
331 Eigen::VectorXd const& local_x,
332 Eigen::VectorXd const& local_x_prev,
333 int const process_id,
334 std::vector<double>& local_b_data,
335 std::vector<double>& local_Jac_data)
336{
337 // For the equations with pressure
338 if (process_id == _process_data.heat_conduction_process_id)
339 {
340 assembleWithJacobianForHeatConductionEquations(
341 t, dt, local_x, local_x_prev, local_b_data, local_Jac_data);
342 return;
343 }
344
345 // For the equations with deformation
346 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_x_prev,
347 local_b_data, local_Jac_data);
348}
349
350template <typename ShapeFunction, int DisplacementDim>
353 const double t, double const dt, Eigen::VectorXd const& local_x,
354 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
355 std::vector<double>& local_Jac_data)
356{
357 auto const local_T =
358 local_x.template segment<temperature_size>(temperature_index);
359
360 auto const local_T_prev =
361 local_x_prev.template segment<temperature_size>(temperature_index);
362
363 auto const local_u =
364 local_x.template segment<displacement_size>(displacement_index);
365 bool const is_u_non_zero = local_u.norm() > 0.0;
366
367 auto local_Jac = MathLib::createZeroedMatrix<
368 typename ShapeMatricesType::template MatrixType<displacement_size,
369 displacement_size>>(
370 local_Jac_data, displacement_size, displacement_size);
371
372 auto local_rhs = MathLib::createZeroedVector<
373 typename ShapeMatricesType::template VectorType<displacement_size>>(
374 local_b_data, displacement_size);
375
376 unsigned const n_integration_points =
377 _integration_method.getNumberOfPoints();
378
379 MPL::VariableArray variables;
380 MPL::VariableArray variables_prev;
382 x_position.setElementID(_element.getID());
383 auto const& medium = _process_data.media_map.getMedium(_element.getID());
384 auto const& solid_phase = medium->phase("Solid");
385
386 for (unsigned ip = 0; ip < n_integration_points; ip++)
387 {
388 x_position.setIntegrationPoint(ip);
389 auto const& w = _ip_data[ip].integration_weight;
390 auto const& N = _ip_data[ip].N;
391 auto const& dNdx = _ip_data[ip].dNdx;
392
393 auto const x_coord =
394 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
395 _element, N);
396 auto const& B =
397 LinearBMatrix::computeBMatrix<DisplacementDim,
398 ShapeFunction::NPOINTS,
400 dNdx, N, x_coord, _is_axially_symmetric);
401
402 auto& sigma = _ip_data[ip].sigma;
403 auto const& sigma_prev = _ip_data[ip].sigma_prev;
404
405 auto& eps = _ip_data[ip].eps;
406 auto const& eps_prev = _ip_data[ip].eps_prev;
407
408 auto& eps_m = _ip_data[ip].eps_m;
409 auto const& eps_m_prev = _ip_data[ip].eps_m_prev;
410
411 auto& state = _ip_data[ip].material_state_variables;
412
413 const double T_ip = N.dot(local_T); // T at integration point
414 variables.temperature = T_ip;
415 double const dT_ip = T_ip - N.dot(local_T_prev);
416
417 //
418 // displacement equation, displacement part
419 //
420 // For the restart computation, the displacement may not be
421 // reloaded but the initial strains are always available. For such case,
422 // the following computation is skipped.
423 if (is_u_non_zero)
424 {
425 eps.noalias() = B * local_u;
426 }
427
428 // Consider also anisotropic thermal expansion.
429 auto const solid_linear_thermal_expansivity_vector =
430 MPL::formKelvinVector<DisplacementDim>(
431 solid_phase
432 .property(
434 .value(variables, x_position, t, dt));
435
437 dthermal_strain = solid_linear_thermal_expansivity_vector * dT_ip;
438
439 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
440
441 variables_prev.stress
443 sigma_prev);
444 variables_prev.mechanical_strain
446 eps_m_prev);
447 variables_prev.temperature = T_ip;
448 variables.mechanical_strain
450 eps_m);
451
452 auto&& solution = _ip_data[ip].solid_material.integrateStress(
453 variables_prev, variables, t, x_position, dt, *state);
454
455 if (!solution)
456 {
457 OGS_FATAL("Computation of local constitutive relation failed.");
458 }
459
461 std::tie(sigma, state, C) = std::move(*solution);
462
463 local_Jac.noalias() += B.transpose() * C * B * w;
464
465 typename ShapeMatricesType::template MatrixType<DisplacementDim,
466 displacement_size>
467 N_u = ShapeMatricesType::template MatrixType<
468 DisplacementDim, displacement_size>::Zero(DisplacementDim,
469 displacement_size);
470
471 for (int i = 0; i < DisplacementDim; ++i)
472 {
473 N_u.template block<1, displacement_size / DisplacementDim>(
474 i, i * displacement_size / DisplacementDim)
475 .noalias() = N;
476 }
477
478 auto const rho_s =
479 solid_phase.property(MPL::PropertyType::density)
480 .template value<double>(variables, x_position, t, dt);
481
482 auto const& b = _process_data.specific_body_force;
483 local_rhs.noalias() -=
484 (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
485 }
486}
487
488template <typename ShapeFunction, int DisplacementDim>
491 const double t, double const dt, Eigen::VectorXd const& local_x,
492 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
493 std::vector<double>& local_Jac_data)
494{
495 auto const local_T =
496 local_x.template segment<temperature_size>(temperature_index);
497
498 auto const local_T_prev =
499 local_x_prev.template segment<temperature_size>(temperature_index);
500
501 auto local_Jac = MathLib::createZeroedMatrix<
502 typename ShapeMatricesType::template MatrixType<temperature_size,
503 temperature_size>>(
504 local_Jac_data, temperature_size, temperature_size);
505
506 auto local_rhs = MathLib::createZeroedVector<
507 typename ShapeMatricesType::template VectorType<temperature_size>>(
508 local_b_data, temperature_size);
509
511 mass.setZero(temperature_size, temperature_size);
512
513 typename ShapeMatricesType::NodalMatrixType laplace;
514 laplace.setZero(temperature_size, temperature_size);
515
516 unsigned const n_integration_points =
517 _integration_method.getNumberOfPoints();
518
520 x_position.setElementID(_element.getID());
521 auto const& medium = _process_data.media_map.getMedium(_element.getID());
522 auto const& solid_phase = medium->phase("Solid");
523 MPL::VariableArray variables;
524
525 for (unsigned ip = 0; ip < n_integration_points; ip++)
526 {
527 x_position.setIntegrationPoint(ip);
528 auto const& w = _ip_data[ip].integration_weight;
529 auto const& N = _ip_data[ip].N;
530 auto const& dNdx = _ip_data[ip].dNdx;
531
532 const double T_ip = N.dot(local_T); // T at integration point
533 variables.temperature = T_ip;
534
535 auto const rho_s =
536 solid_phase.property(MPL::PropertyType::density)
537 .template value<double>(variables, x_position, t, dt);
538 auto const c_p =
539 solid_phase.property(MPL::PropertyType::specific_heat_capacity)
540 .template value<double>(variables, x_position, t, dt);
541
542 mass.noalias() += N.transpose() * rho_s * c_p * N * w;
543
544 auto const lambda =
545 solid_phase
546 .property(
548 .value(variables, x_position, t, dt);
549
550 GlobalDimMatrixType const thermal_conductivity =
551 MaterialPropertyLib::formEigenTensor<DisplacementDim>(lambda);
552
553 laplace.noalias() += dNdx.transpose() * thermal_conductivity * dNdx * w;
554 }
555 local_Jac.noalias() += laplace + mass / dt;
556
557 local_rhs.noalias() -=
558 laplace * local_T + mass * (local_T - local_T_prev) / dt;
559}
560
561template <typename ShapeFunction, int DisplacementDim>
562std::size_t
564 double const* values)
565{
566 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
567 values, _ip_data, &IpData::sigma);
568}
569
570template <typename ShapeFunction, int DisplacementDim>
571std::vector<double>
573{
574 constexpr int kelvin_vector_size =
576
577 return transposeInPlace<kelvin_vector_size>(
578 [this](std::vector<double>& values)
579 { return getIntPtSigma(0, {}, {}, values); });
580}
581
582template <typename ShapeFunction, int DisplacementDim>
583std::vector<double> const&
585 const double /*t*/,
586 std::vector<GlobalVector*> const& /*x*/,
587 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
588 std::vector<double>& cache) const
589{
590 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
591 _ip_data, &IpData::sigma, cache);
592}
593
594template <typename ShapeFunction, int DisplacementDim>
595std::size_t
597 double const* values)
598{
599 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
600 values, _ip_data, &IpData::eps);
601}
602
603template <typename ShapeFunction, int DisplacementDim>
604std::vector<double> ThermoMechanicsLocalAssembler<
605 ShapeFunction, DisplacementDim>::getEpsilon() const
606{
607 constexpr int kelvin_vector_size =
609
610 return transposeInPlace<kelvin_vector_size>(
611 [this](std::vector<double>& values)
612 { return getIntPtEpsilon(0, {}, {}, values); });
613}
614
615template <typename ShapeFunction, int DisplacementDim>
616std::vector<double> const&
618 const double /*t*/,
619 std::vector<GlobalVector*> const& /*x*/,
620 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
621 std::vector<double>& cache) const
622{
623 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
624 _ip_data, &IpData::eps, cache);
625}
626
627template <typename ShapeFunction, int DisplacementDim>
628std::vector<double> const&
631 const double /*t*/,
632 std::vector<GlobalVector*> const& /*x*/,
633 std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
634 std::vector<double>& cache) const
635{
636 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
637 _ip_data, &IpData::eps_m, cache);
638}
639
640template <typename ShapeFunction, int DisplacementDim>
642 ShapeFunction, DisplacementDim>::setEpsilonMechanical(double const* values)
643{
644 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
645 values, _ip_data, &IpData::eps_m);
646}
647template <typename ShapeFunction, int DisplacementDim>
648std::vector<double> ThermoMechanicsLocalAssembler<
649 ShapeFunction, DisplacementDim>::getEpsilonMechanical() const
650{
651 constexpr int kelvin_vector_size =
653
654 return transposeInPlace<kelvin_vector_size>(
655 [this](std::vector<double>& values)
656 { return getIntPtEpsilonMechanical(0, {}, {}, values); });
657}
658
659template <typename ShapeFunction, int DisplacementDim>
661 ShapeFunction, DisplacementDim>::getNumberOfIntegrationPoints() const
662{
663 return _integration_method.getNumberOfPoints();
664}
665
666template <typename ShapeFunction, int DisplacementDim>
667int ThermoMechanicsLocalAssembler<ShapeFunction,
668 DisplacementDim>::getMaterialID() const
669{
670 return _process_data.material_ids == nullptr
671 ? 0
672 : (*_process_data.material_ids)[_element.getID()];
673}
674
675template <typename ShapeFunction, int DisplacementDim>
677 DisplacementDim>::MaterialStateVariables const&
679 getMaterialStateVariablesAt(unsigned integration_point) const
680{
681 return *_ip_data[integration_point].material_state_variables;
682}
683} // namespace ThermoMechanics
684} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
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)
void setIntegrationPoint(unsigned integration_point)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
std::vector< double > const & getIntPtEpsilonMechanical(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
SecondaryData< typename ShapeMatrices::ShapeType > _secondary_data
std::size_t setIPDataInitialConditions(std::string_view const name, double const *values, int const integration_order) override
Returns number of read integration points.
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
ThermoMechanicsLocalAssembler(ThermoMechanicsLocalAssembler const &)=delete
MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt(unsigned integration_point) const override
ShapeMatrixPolicyType< ShapeFunction, DisplacementDim > ShapeMatricesType
ThermoMechanicsProcessData< DisplacementDim > & _process_data
NumLib::GenericIntegrationMethod const & _integration_method
void assembleWithJacobianForHeatConductionEquations(const double t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void assembleWithJacobianForDeformationEquations(const double t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void assembleWithJacobianForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
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
void assembleWithJacobian(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_rhs_data, std::vector< double > &local_Jac_data) override
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)
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.
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N