32 bool const is_axially_symmetric,
34 : _process_data(process_data),
35 _integration_method(integration_method),
37 _is_axially_symmetric(is_axially_symmetric)
39 unsigned const n_integration_points =
42 _ip_data.reserve(n_integration_points);
45 auto const shape_matrices =
47 DisplacementDim>(e, is_axially_symmetric,
53 for (
unsigned ip = 0; ip < n_integration_points; ip++)
55 _ip_data.emplace_back(solid_material);
57 ip_data.integration_weight =
59 shape_matrices[ip].integralMeasure * shape_matrices[ip].detJ;
61 static const int kelvin_vector_size =
64 ip_data.sigma.setZero(kelvin_vector_size);
65 ip_data.eps.setZero(kelvin_vector_size);
68 ip_data.sigma_prev.resize(kelvin_vector_size);
69 ip_data.eps_prev.resize(kelvin_vector_size);
71 ip_data.eps_m.setZero(kelvin_vector_size);
72 ip_data.eps_m_prev.setZero(kelvin_vector_size);
75 ip_data.N = shape_matrices[ip].N;
76 ip_data.dNdx = shape_matrices[ip].dNdx;
151 std::vector<double>
const& local_x,
152 std::vector<double>
const& local_x_prev,
153 std::vector<double>& local_rhs_data,
154 std::vector<double>& local_Jac_data)
156 auto const local_matrix_size = local_x.size();
157 assert(local_matrix_size == temperature_size + displacement_size);
159 auto T = Eigen::Map<
typename ShapeMatricesType::template VectorType<
160 temperature_size>
const>(local_x.data() + temperature_index,
163 auto u = Eigen::Map<
typename ShapeMatricesType::template VectorType<
164 displacement_size>
const>(local_x.data() + displacement_index,
166 bool const is_u_non_zero = u.norm() > 0.0;
168 auto T_prev = Eigen::Map<
typename ShapeMatricesType::template VectorType<
169 temperature_size>
const>(local_x_prev.data() + temperature_index,
173 local_Jac_data, local_matrix_size, local_matrix_size);
178 typename ShapeMatricesType::template MatrixType<displacement_size,
181 KuT.setZero(displacement_size, temperature_size);
184 KTT.setZero(temperature_size, temperature_size);
187 DTT.setZero(temperature_size, temperature_size);
189 unsigned const n_integration_points =
190 _integration_method.getNumberOfPoints();
197 auto const& medium = _process_data.media_map.getMedium(_element.getID());
198 auto const& solid_phase = medium->phase(
"Solid");
200 for (
unsigned ip = 0; ip < n_integration_points; ip++)
202 auto const& w = _ip_data[ip].integration_weight;
203 auto const& N = _ip_data[ip].N;
204 auto const& dNdx = _ip_data[ip].dNdx;
207 std::nullopt, _element.getID(),
215 ShapeFunction::NPOINTS,
217 dNdx, N, x_coord, _is_axially_symmetric);
219 auto& sigma = _ip_data[ip].sigma;
220 auto const& sigma_prev = _ip_data[ip].sigma_prev;
222 auto& eps = _ip_data[ip].eps;
223 auto const& eps_prev = _ip_data[ip].eps_prev;
225 auto& eps_m = _ip_data[ip].eps_m;
226 auto const& eps_m_prev = _ip_data[ip].eps_m_prev;
228 auto& state = _ip_data[ip].material_state_variables;
230 const double T_ip = N.dot(T);
231 double const T_prev_ip = N.dot(T_prev);
234 auto const solid_linear_thermal_expansivity_vector =
239 .value(variables, x_position, t, dt));
243 solid_linear_thermal_expansivity_vector * (T_ip - T_prev_ip);
253 eps.noalias() = B * u;
256 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
270 auto&& solution = _ip_data[ip].solid_material.integrateStress(
271 variables_prev, variables, t, x_position, dt, *state);
275 OGS_FATAL(
"Computation of local constitutive relation failed.");
279 std::tie(sigma, state, C) = std::move(*solution);
282 .template block<displacement_size, displacement_size>(
283 displacement_index, displacement_index)
284 .noalias() += B.transpose() * C * B * w;
286 typename ShapeMatricesType::template MatrixType<DisplacementDim,
288 N_u = ShapeMatricesType::template MatrixType<
289 DisplacementDim, displacement_size>::Zero(DisplacementDim,
292 for (
int i = 0; i < DisplacementDim; ++i)
294 N_u.template block<1, displacement_size / DisplacementDim>(
295 i, i * displacement_size / DisplacementDim)
300 solid_phase.property(MPL::PropertyType::density)
301 .template value<double>(variables, x_position, t, dt);
303 auto const& b = _process_data.specific_body_force;
304 local_rhs.template block<displacement_size, 1>(displacement_index, 0)
306 (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
311 auto const alpha_T_tensor =
313 solid_linear_thermal_expansivity_vector);
314 KuT.noalias() += B.transpose() * (C * alpha_T_tensor) * N * w;
316 if (_ip_data[ip].solid_material.getConstitutiveModel() ==
319 auto const s = Invariants::deviatoric_projection * sigma;
320 double const norm_s = Invariants::FrobeniusNorm(s);
321 const double creep_coefficient =
322 _ip_data[ip].solid_material.getTemperatureRelatedCoefficient(
323 t, dt, x_position, T_ip, norm_s);
324 KuT.noalias() += creep_coefficient * B.transpose() * s * N * w;
334 .value(variables, x_position, t, dt);
339 KTT.noalias() += dNdx.transpose() * thermal_conductivity * dNdx * w;
345 .template value<double>(variables, x_position, t, dt);
346 DTT.noalias() += N.transpose() * rho_s * c * N * w;
351 .template block<temperature_size, temperature_size>(temperature_index,
353 .noalias() += KTT + DTT / dt;
357 .template block<displacement_size, temperature_size>(displacement_index,
361 local_rhs.template block<temperature_size, 1>(temperature_index, 0)
362 .noalias() -= KTT * T + DTT * (T - T_prev) / dt;
390 const double t,
double const dt, Eigen::VectorXd
const& local_x,
391 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
392 std::vector<double>& local_Jac_data)
395 local_x.template segment<temperature_size>(temperature_index);
397 auto const local_T_prev =
398 local_x_prev.template segment<temperature_size>(temperature_index);
401 local_x.template segment<displacement_size>(displacement_index);
402 bool const is_u_non_zero = local_u.norm() > 0.0;
405 typename ShapeMatricesType::template MatrixType<displacement_size,
407 local_Jac_data, displacement_size, displacement_size);
410 typename ShapeMatricesType::template VectorType<displacement_size>>(
411 local_b_data, displacement_size);
413 unsigned const n_integration_points =
414 _integration_method.getNumberOfPoints();
418 auto const& medium = _process_data.media_map.getMedium(_element.getID());
419 auto const& solid_phase = medium->phase(
"Solid");
421 for (
unsigned ip = 0; ip < n_integration_points; ip++)
423 auto const& w = _ip_data[ip].integration_weight;
424 auto const& N = _ip_data[ip].N;
425 auto const& dNdx = _ip_data[ip].dNdx;
428 std::nullopt, _element.getID(),
437 ShapeFunction::NPOINTS,
439 dNdx, N, x_coord, _is_axially_symmetric);
441 auto& sigma = _ip_data[ip].sigma;
442 auto const& sigma_prev = _ip_data[ip].sigma_prev;
444 auto& eps = _ip_data[ip].eps;
445 auto const& eps_prev = _ip_data[ip].eps_prev;
447 auto& eps_m = _ip_data[ip].eps_m;
448 auto const& eps_m_prev = _ip_data[ip].eps_m_prev;
450 auto& state = _ip_data[ip].material_state_variables;
452 const double T_ip = N.dot(local_T);
454 double const dT_ip = T_ip - N.dot(local_T_prev);
464 eps.noalias() = B * local_u;
468 auto const solid_linear_thermal_expansivity_vector =
473 .value(variables, x_position, t, dt));
476 dthermal_strain = solid_linear_thermal_expansivity_vector * dT_ip;
478 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
491 auto&& solution = _ip_data[ip].solid_material.integrateStress(
492 variables_prev, variables, t, x_position, dt, *state);
496 OGS_FATAL(
"Computation of local constitutive relation failed.");
500 std::tie(sigma, state, C) = std::move(*solution);
502 local_Jac.noalias() += B.transpose() * C * B * w;
504 typename ShapeMatricesType::template MatrixType<DisplacementDim,
506 N_u = ShapeMatricesType::template MatrixType<
507 DisplacementDim, displacement_size>::Zero(DisplacementDim,
510 for (
int i = 0; i < DisplacementDim; ++i)
512 N_u.template block<1, displacement_size / DisplacementDim>(
513 i, i * displacement_size / DisplacementDim)
518 solid_phase.property(MPL::PropertyType::density)
519 .template value<double>(variables, x_position, t, dt);
521 auto const& b = _process_data.specific_body_force;
522 local_rhs.noalias() -=
523 (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
530 const double t,
double const dt, Eigen::VectorXd
const& local_x,
531 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
532 std::vector<double>& local_Jac_data)
535 local_x.template segment<temperature_size>(temperature_index);
537 auto const local_T_prev =
538 local_x_prev.template segment<temperature_size>(temperature_index);
541 typename ShapeMatricesType::template MatrixType<temperature_size,
543 local_Jac_data, temperature_size, temperature_size);
546 typename ShapeMatricesType::template VectorType<temperature_size>>(
547 local_b_data, temperature_size);
550 mass.setZero(temperature_size, temperature_size);
553 laplace.setZero(temperature_size, temperature_size);
555 unsigned const n_integration_points =
556 _integration_method.getNumberOfPoints();
558 auto const& medium = _process_data.media_map.getMedium(_element.getID());
559 auto const& solid_phase = medium->phase(
"Solid");
562 for (
unsigned ip = 0; ip < n_integration_points; ip++)
564 auto const& w = _ip_data[ip].integration_weight;
565 auto const& N = _ip_data[ip].N;
566 auto const& dNdx = _ip_data[ip].dNdx;
569 std::nullopt, _element.getID(),
574 const double T_ip = N.dot(local_T);
578 solid_phase.property(MPL::PropertyType::density)
579 .template value<double>(variables, x_position, t, dt);
581 solid_phase.property(MPL::PropertyType::specific_heat_capacity)
582 .template value<double>(variables, x_position, t, dt);
584 mass.noalias() += N.transpose() * rho_s * c_p * N * w;
590 .value(variables, x_position, t, dt);
595 laplace.noalias() += dNdx.transpose() * thermal_conductivity * dNdx * w;
597 local_Jac.noalias() += laplace + mass / dt;
599 local_rhs.noalias() -=
600 laplace * local_T + mass * (local_T - local_T_prev) / dt;