24 namespace ThermoMechanics
26 template <
typename ShapeFunction,
typename IntegrationMethod,
28 ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
30 ThermoMechanicsLocalAssembler(
33 bool const is_axially_symmetric,
34 unsigned const integration_order,
36 : _process_data(process_data),
37 _integration_method(integration_order),
39 _is_axially_symmetric(is_axially_symmetric)
41 unsigned const n_integration_points =
44 _ip_data.reserve(n_integration_points);
47 auto const shape_matrices =
49 DisplacementDim>(e, is_axially_symmetric,
55 for (
unsigned ip = 0; ip < n_integration_points; ip++)
57 _ip_data.emplace_back(solid_material);
59 ip_data.integration_weight =
61 shape_matrices[ip].integralMeasure * shape_matrices[ip].detJ;
63 static const int kelvin_vector_size =
66 ip_data.sigma.setZero(kelvin_vector_size);
67 ip_data.eps.setZero(kelvin_vector_size);
70 ip_data.sigma_prev.resize(kelvin_vector_size);
71 ip_data.eps_prev.resize(kelvin_vector_size);
73 ip_data.eps_m.setZero(kelvin_vector_size);
74 ip_data.eps_m_prev.setZero(kelvin_vector_size);
77 ip_data.N = shape_matrices[ip].N;
78 ip_data.dNdx = shape_matrices[ip].dNdx;
84 template <
typename ShapeFunction,
typename IntegrationMethod,
87 ShapeFunction, IntegrationMethod,
88 DisplacementDim>::setIPDataInitialConditions(std::string
const&
name,
90 int const integration_order)
92 if (integration_order !=
93 static_cast<int>(_integration_method.getIntegrationOrder()))
96 "Setting integration point initial conditions; The integration "
97 "order of the local assembler for element {:d} is different from "
98 "the integration order in the initial condition.",
102 if (
name ==
"sigma_ip")
104 return setSigma(values);
106 if (
name ==
"epsilon_ip")
108 return setEpsilon(values);
110 if (
name ==
"epsilon_m_ip")
112 return setEpsilonMechanical(values);
118 template <
typename ShapeFunction,
typename IntegrationMethod,
122 assembleWithJacobian(
double const t,
double const dt,
123 std::vector<double>
const& local_x,
124 std::vector<double>
const& local_xdot,
125 const double ,
const double ,
126 std::vector<double>& ,
127 std::vector<double>& ,
128 std::vector<double>& local_rhs_data,
129 std::vector<double>& local_Jac_data)
131 auto const local_matrix_size = local_x.size();
132 assert(local_matrix_size == temperature_size + displacement_size);
134 auto T = Eigen::Map<
typename ShapeMatricesType::template VectorType<
135 temperature_size>
const>(local_x.data() + temperature_index,
138 auto u = Eigen::Map<
typename ShapeMatricesType::template VectorType<
139 displacement_size>
const>(local_x.data() + displacement_index,
141 bool const is_u_non_zero = u.norm() > 0.0;
143 auto T_dot = Eigen::Map<
typename ShapeMatricesType::template VectorType<
144 temperature_size>
const>(local_xdot.data() + temperature_index,
147 auto local_Jac = MathLib::createZeroedMatrix<JacobianMatrix>(
148 local_Jac_data, local_matrix_size, local_matrix_size);
150 auto local_rhs = MathLib::createZeroedVector<RhsVector>(local_rhs_data,
153 typename ShapeMatricesType::template MatrixType<displacement_size,
156 KuT.setZero(displacement_size, temperature_size);
159 KTT.setZero(temperature_size, temperature_size);
162 DTT.setZero(temperature_size, temperature_size);
164 unsigned const n_integration_points =
165 _integration_method.getNumberOfPoints();
172 auto const& medium = _process_data.media_map->getMedium(_element.getID());
173 auto const& solid_phase = medium->phase(
"Solid");
175 for (
unsigned ip = 0; ip < n_integration_points; ip++)
178 auto const& w = _ip_data[ip].integration_weight;
179 auto const& N = _ip_data[ip].N;
180 auto const& dNdx = _ip_data[ip].dNdx;
183 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
187 ShapeFunction::NPOINTS,
189 dNdx, N, x_coord, _is_axially_symmetric);
191 auto& sigma = _ip_data[ip].sigma;
192 auto const& sigma_prev = _ip_data[ip].sigma_prev;
194 auto& eps = _ip_data[ip].eps;
195 auto const& eps_prev = _ip_data[ip].eps_prev;
197 auto& eps_m = _ip_data[ip].eps_m;
198 auto const& eps_m_prev = _ip_data[ip].eps_m_prev;
200 auto& state = _ip_data[ip].material_state_variables;
202 const double T_ip = N.dot(T);
203 double const dT = N.dot(T_dot) * dt;
206 auto const solid_linear_thermal_expansivity_vector =
207 MPL::formKelvinVectorFromThermalExpansivity<DisplacementDim>(
211 .value(variables, x_position, t, dt));
215 solid_linear_thermal_expansivity_vector.eval() * dT;
225 eps.noalias() = B * u;
228 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
230 variables_prev[
static_cast<int>(MPL::Variable::stress)]
233 variables_prev[
static_cast<int>(MPL::Variable::mechanical_strain)]
237 .emplace<double>(T_ip);
238 variables[
static_cast<int>(MPL::Variable::mechanical_strain)]
244 auto&& solution = _ip_data[ip].solid_material.integrateStress(
245 variables_prev, variables, t, x_position, dt, *state);
249 OGS_FATAL(
"Computation of local constitutive relation failed.");
253 std::tie(sigma, state, C) = std::move(*solution);
256 .template block<displacement_size, displacement_size>(
257 displacement_index, displacement_index)
258 .noalias() += B.transpose() * C * B * w;
260 typename ShapeMatricesType::template MatrixType<DisplacementDim,
262 N_u = ShapeMatricesType::template MatrixType<
263 DisplacementDim, displacement_size>::Zero(DisplacementDim,
266 for (
int i = 0; i < DisplacementDim; ++i)
268 N_u.template block<1, displacement_size / DisplacementDim>(
269 i, i * displacement_size / DisplacementDim)
275 .template value<double>(variables, x_position, t, dt);
277 auto const& b = _process_data.specific_body_force;
278 local_rhs.template block<displacement_size, 1>(displacement_index, 0)
280 (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
285 auto const alpha_T_tensor =
287 solid_linear_thermal_expansivity_vector.eval());
288 KuT.noalias() += B.transpose() * (C * alpha_T_tensor.eval()) * N * w;
290 if (_ip_data[ip].solid_material.getConstitutiveModel() ==
293 auto const s = Invariants::deviatoric_projection * sigma;
294 double const norm_s = Invariants::FrobeniusNorm(s);
295 const double creep_coefficient =
296 _ip_data[ip].solid_material.getTemperatureRelatedCoefficient(
297 t, dt, x_position, T_ip, norm_s);
298 KuT.noalias() += creep_coefficient * B.transpose() * s * N * w;
308 .value(variables, x_position, t, dt);
311 MaterialPropertyLib::formEigenTensor<DisplacementDim>(lambda);
319 .template value<double>(variables, x_position, t, dt);
320 DTT.noalias() += N.transpose() * rho_s *
c * N * w;
325 .template block<temperature_size, temperature_size>(temperature_index,
327 .noalias() += KTT + DTT / dt;
331 .template block<displacement_size, temperature_size>(displacement_index,
335 local_rhs.template block<temperature_size, 1>(temperature_index, 0)
336 .noalias() -= KTT * T + DTT * T_dot;
339 template <
typename ShapeFunction,
typename IntegrationMethod,
343 assembleWithJacobianForStaggeredScheme(
344 const double t,
double const dt, Eigen::VectorXd
const& local_x,
345 Eigen::VectorXd
const& local_xdot,
const double ,
346 const double ,
int const process_id,
347 std::vector<double>& ,
348 std::vector<double>& ,
349 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
352 if (process_id == _process_data.heat_conduction_process_id)
354 assembleWithJacobianForHeatConductionEquations(
355 t, dt, local_x, local_xdot, local_b_data, local_Jac_data);
360 assembleWithJacobianForDeformationEquations(t, dt, local_x, local_xdot,
361 local_b_data, local_Jac_data);
364 template <
typename ShapeFunction,
typename IntegrationMethod,
368 assembleWithJacobianForDeformationEquations(
369 const double t,
double const dt, Eigen::VectorXd
const& local_x,
370 Eigen::VectorXd
const& local_xdot, std::vector<double>& local_b_data,
371 std::vector<double>& local_Jac_data)
374 local_x.template segment<temperature_size>(temperature_index);
376 auto const local_Tdot =
377 local_xdot.template segment<temperature_size>(temperature_index);
380 local_x.template segment<displacement_size>(displacement_index);
381 bool const is_u_non_zero = local_u.norm() > 0.0;
384 typename ShapeMatricesType::template MatrixType<displacement_size,
386 local_Jac_data, displacement_size, displacement_size);
389 typename ShapeMatricesType::template VectorType<displacement_size>>(
390 local_b_data, displacement_size);
392 unsigned const n_integration_points =
393 _integration_method.getNumberOfPoints();
399 auto const& medium = _process_data.media_map->getMedium(_element.getID());
400 auto const& solid_phase = medium->phase(
"Solid");
402 for (
unsigned ip = 0; ip < n_integration_points; ip++)
405 auto const& w = _ip_data[ip].integration_weight;
406 auto const& N = _ip_data[ip].N;
407 auto const& dNdx = _ip_data[ip].dNdx;
410 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
414 ShapeFunction::NPOINTS,
416 dNdx, N, x_coord, _is_axially_symmetric);
418 auto& sigma = _ip_data[ip].sigma;
419 auto const& sigma_prev = _ip_data[ip].sigma_prev;
421 auto& eps = _ip_data[ip].eps;
422 auto const& eps_prev = _ip_data[ip].eps_prev;
424 auto& eps_m = _ip_data[ip].eps_m;
425 auto const& eps_m_prev = _ip_data[ip].eps_m_prev;
427 auto& state = _ip_data[ip].material_state_variables;
429 const double T_ip = N.dot(local_T);
430 double const dT_ip = N.dot(local_Tdot) * dt;
442 eps.noalias() = B * local_u;
446 auto const solid_linear_thermal_expansivity_vector =
447 MPL::formKelvinVectorFromThermalExpansivity<DisplacementDim>(
451 .value(variables, x_position, t, dt));
455 solid_linear_thermal_expansivity_vector.eval() * dT_ip;
457 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
459 variables_prev[
static_cast<int>(MPL::Variable::stress)]
462 variables_prev[
static_cast<int>(MPL::Variable::mechanical_strain)]
466 .emplace<double>(T_ip);
467 variables[
static_cast<int>(MPL::Variable::mechanical_strain)]
471 auto&& solution = _ip_data[ip].solid_material.integrateStress(
472 variables_prev, variables, t, x_position, dt, *state);
476 OGS_FATAL(
"Computation of local constitutive relation failed.");
480 std::tie(sigma, state, C) = std::move(*solution);
482 local_Jac.noalias() += B.transpose() * C * B * w;
484 typename ShapeMatricesType::template MatrixType<DisplacementDim,
486 N_u = ShapeMatricesType::template MatrixType<
487 DisplacementDim, displacement_size>::Zero(DisplacementDim,
490 for (
int i = 0; i < DisplacementDim; ++i)
492 N_u.template block<1, displacement_size / DisplacementDim>(
493 i, i * displacement_size / DisplacementDim)
499 .template value<double>(variables, x_position, t, dt);
501 auto const& b = _process_data.specific_body_force;
502 local_rhs.noalias() -=
503 (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
507 template <
typename ShapeFunction,
typename IntegrationMethod,
511 assembleWithJacobianForHeatConductionEquations(
512 const double t,
double const dt, Eigen::VectorXd
const& local_x,
513 Eigen::VectorXd
const& local_xdot, std::vector<double>& local_b_data,
514 std::vector<double>& local_Jac_data)
517 local_x.template segment<temperature_size>(temperature_index);
519 auto const local_dT =
520 local_xdot.template segment<temperature_size>(temperature_index) * dt;
523 typename ShapeMatricesType::template MatrixType<temperature_size,
525 local_Jac_data, temperature_size, temperature_size);
528 typename ShapeMatricesType::template VectorType<temperature_size>>(
529 local_b_data, temperature_size);
532 mass.setZero(temperature_size, temperature_size);
535 laplace.setZero(temperature_size, temperature_size);
537 unsigned const n_integration_points =
538 _integration_method.getNumberOfPoints();
542 auto const& medium = _process_data.media_map->getMedium(_element.getID());
543 auto const& solid_phase = medium->phase(
"Solid");
546 for (
unsigned ip = 0; ip < n_integration_points; ip++)
549 auto const& w = _ip_data[ip].integration_weight;
550 auto const& N = _ip_data[ip].N;
551 auto const& dNdx = _ip_data[ip].dNdx;
553 const double T_ip = N.dot(local_T);
559 .template value<double>(variables, x_position, t, dt);
562 .template value<double>(variables, x_position, t, dt);
564 mass.noalias() += N.transpose() * rho_s * c_p * N * w;
570 .value(variables, x_position, t, dt);
573 MaterialPropertyLib::formEigenTensor<DisplacementDim>(lambda);
577 local_Jac.noalias() += laplace + mass / dt;
579 local_rhs.noalias() -= laplace * local_T + mass * local_dT / dt;
582 template <
typename ShapeFunction,
typename IntegrationMethod,
586 DisplacementDim>::setSigma(
double const* values)
588 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
589 values, _ip_data, &IpData::sigma);
592 template <
typename ShapeFunction,
typename IntegrationMethod,
595 ShapeFunction, IntegrationMethod, DisplacementDim>::getSigma()
const
597 constexpr
int kelvin_vector_size =
600 return transposeInPlace<kelvin_vector_size>(
601 [
this](std::vector<double>& values)
602 {
return getIntPtSigma(0, {}, {}, values); });
605 template <
typename ShapeFunction,
typename IntegrationMethod,
608 ShapeFunction, IntegrationMethod, DisplacementDim>::
611 std::vector<GlobalVector*>
const& ,
612 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
613 std::vector<double>& cache)
const
615 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
616 _ip_data, &IpData::sigma, cache);
619 template <
typename ShapeFunction,
typename IntegrationMethod,
623 DisplacementDim>::setEpsilon(
double const* values)
625 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
626 values, _ip_data, &IpData::eps);
629 template <
typename ShapeFunction,
typename IntegrationMethod,
632 ShapeFunction, IntegrationMethod, DisplacementDim>::getEpsilon()
const
634 constexpr
int kelvin_vector_size =
637 return transposeInPlace<kelvin_vector_size>(
638 [
this](std::vector<double>& values)
639 {
return getIntPtEpsilon(0, {}, {}, values); });
642 template <
typename ShapeFunction,
typename IntegrationMethod,
645 ShapeFunction, IntegrationMethod, DisplacementDim>::
648 std::vector<GlobalVector*>
const& ,
649 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
650 std::vector<double>& cache)
const
652 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
653 _ip_data, &IpData::eps, cache);
656 template <
typename ShapeFunction,
typename IntegrationMethod,
659 ShapeFunction, IntegrationMethod, DisplacementDim>::
660 getIntPtEpsilonMechanical(
662 std::vector<GlobalVector*>
const& ,
663 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
664 std::vector<double>& cache)
const
666 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
667 _ip_data, &IpData::eps_m, cache);
670 template <
typename ShapeFunction,
typename IntegrationMethod,
673 ShapeFunction, IntegrationMethod,
674 DisplacementDim>::setEpsilonMechanical(
double const* values)
676 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
677 values, _ip_data, &IpData::eps_m);
679 template <
typename ShapeFunction,
typename IntegrationMethod,
683 DisplacementDim>::getEpsilonMechanical()
const
685 constexpr
int kelvin_vector_size =
688 return transposeInPlace<kelvin_vector_size>(
689 [
this](std::vector<double>& values)
690 {
return getIntPtEpsilonMechanical(0, {}, {}, values); });
693 template <
typename ShapeFunction,
typename IntegrationMethod,
696 ShapeFunction, IntegrationMethod,
697 DisplacementDim>::getNumberOfIntegrationPoints()
const
699 return _integration_method.getNumberOfPoints();
702 template <
typename ShapeFunction,
typename IntegrationMethod,
705 DisplacementDim>::MaterialStateVariables
const&
708 getMaterialStateVariablesAt(
unsigned integration_point)
const
710 return *_ip_data[integration_point].material_state_variables;
virtual std::size_t getID() const final
Returns the ID of the element.
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Local assembler of ThermoMechanics process.
MeshLib::Element const & _element
SecondaryData< typename ShapeMatrices::ShapeType > _secondary_data
ShapeMatrixPolicyType< ShapeFunction, DisplacementDim > ShapeMatricesType
IntegrationMethod _integration_method
ThermoMechanicsProcessData< DisplacementDim > & _process_data
std::vector< IpData, Eigen::aligned_allocator< IpData > > _ip_data
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
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)
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
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
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
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