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;
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)
122 auto const local_matrix_size = local_x.size();
123 assert(local_matrix_size == temperature_size + displacement_size);
125 auto T = Eigen::Map<
typename ShapeMatricesType::template VectorType<
126 temperature_size>
const>(local_x.data() + temperature_index,
129 auto u = Eigen::Map<
typename ShapeMatricesType::template VectorType<
130 displacement_size>
const>(local_x.data() + displacement_index,
132 bool const is_u_non_zero = u.norm() > 0.0;
134 auto T_prev = Eigen::Map<
typename ShapeMatricesType::template VectorType<
135 temperature_size>
const>(local_x_prev.data() + temperature_index,
139 local_Jac_data, local_matrix_size, local_matrix_size);
144 typename ShapeMatricesType::template MatrixType<displacement_size,
147 KuT.setZero(displacement_size, temperature_size);
150 KTT.setZero(temperature_size, temperature_size);
153 DTT.setZero(temperature_size, temperature_size);
155 unsigned const n_integration_points =
156 _integration_method.getNumberOfPoints();
163 auto const& medium = _process_data.media_map.getMedium(_element.getID());
164 auto const& solid_phase = medium->phase(
"Solid");
166 for (
unsigned ip = 0; ip < n_integration_points; ip++)
168 auto const& w = _ip_data[ip].integration_weight;
169 auto const& N = _ip_data[ip].N;
170 auto const& dNdx = _ip_data[ip].dNdx;
173 std::nullopt, _element.getID(),
181 ShapeFunction::NPOINTS,
183 dNdx, N, x_coord, _is_axially_symmetric);
185 auto& sigma = _ip_data[ip].sigma;
186 auto const& sigma_prev = _ip_data[ip].sigma_prev;
188 auto& eps = _ip_data[ip].eps;
189 auto const& eps_prev = _ip_data[ip].eps_prev;
191 auto& eps_m = _ip_data[ip].eps_m;
192 auto const& eps_m_prev = _ip_data[ip].eps_m_prev;
194 auto& state = _ip_data[ip].material_state_variables;
196 const double T_ip = N.dot(T);
197 double const T_prev_ip = N.dot(T_prev);
200 auto const solid_linear_thermal_expansivity_vector =
205 .value(variables, x_position, t, dt));
209 solid_linear_thermal_expansivity_vector * (T_ip - T_prev_ip);
219 eps.noalias() = B * u;
222 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
236 auto&& solution = _ip_data[ip].solid_material.integrateStress(
237 variables_prev, variables, t, x_position, dt, *state);
241 OGS_FATAL(
"Computation of local constitutive relation failed.");
245 std::tie(sigma, state, C) = std::move(*solution);
248 .template block<displacement_size, displacement_size>(
249 displacement_index, displacement_index)
250 .noalias() += B.transpose() * C * B * w;
252 typename ShapeMatricesType::template MatrixType<DisplacementDim,
254 N_u = ShapeMatricesType::template MatrixType<
255 DisplacementDim, displacement_size>::Zero(DisplacementDim,
258 for (
int i = 0; i < DisplacementDim; ++i)
260 N_u.template block<1, displacement_size / DisplacementDim>(
261 i, i * displacement_size / DisplacementDim)
266 solid_phase.property(MPL::PropertyType::density)
267 .template value<double>(variables, x_position, t, dt);
269 auto const& b = _process_data.specific_body_force;
270 local_rhs.template block<displacement_size, 1>(displacement_index, 0)
272 (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
277 auto const alpha_T_tensor =
279 solid_linear_thermal_expansivity_vector);
280 KuT.noalias() += B.transpose() * (C * alpha_T_tensor) * N * w;
282 if (_ip_data[ip].solid_material.getConstitutiveModel() ==
285 auto const s = Invariants::deviatoric_projection * sigma;
286 double const norm_s = Invariants::FrobeniusNorm(s);
287 const double creep_coefficient =
288 _ip_data[ip].solid_material.getTemperatureRelatedCoefficient(
289 t, dt, x_position, T_ip, norm_s);
290 KuT.noalias() += creep_coefficient * B.transpose() * s * N * w;
300 .value(variables, x_position, t, dt);
305 KTT.noalias() += dNdx.transpose() * thermal_conductivity * dNdx * w;
311 .template value<double>(variables, x_position, t, dt);
312 DTT.noalias() += N.transpose() * rho_s * c * N * w;
317 .template block<temperature_size, temperature_size>(temperature_index,
319 .noalias() += KTT + DTT / dt;
323 .template block<displacement_size, temperature_size>(displacement_index,
327 local_rhs.template block<temperature_size, 1>(temperature_index, 0)
328 .noalias() -= KTT * T + DTT * (T - T_prev) / dt;
356 const double t,
double const dt, Eigen::VectorXd
const& local_x,
357 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
358 std::vector<double>& local_Jac_data)
361 local_x.template segment<temperature_size>(temperature_index);
363 auto const local_T_prev =
364 local_x_prev.template segment<temperature_size>(temperature_index);
367 local_x.template segment<displacement_size>(displacement_index);
368 bool const is_u_non_zero = local_u.norm() > 0.0;
371 typename ShapeMatricesType::template MatrixType<displacement_size,
373 local_Jac_data, displacement_size, displacement_size);
376 typename ShapeMatricesType::template VectorType<displacement_size>>(
377 local_b_data, displacement_size);
379 unsigned const n_integration_points =
380 _integration_method.getNumberOfPoints();
384 auto const& medium = _process_data.media_map.getMedium(_element.getID());
385 auto const& solid_phase = medium->phase(
"Solid");
387 for (
unsigned ip = 0; ip < n_integration_points; 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;
394 std::nullopt, _element.getID(),
399 auto const x_coord = x_position.getCoordinates().value()[0];
403 ShapeFunction::NPOINTS,
405 dNdx, N, x_coord, _is_axially_symmetric);
407 auto& sigma = _ip_data[ip].sigma;
408 auto const& sigma_prev = _ip_data[ip].sigma_prev;
410 auto& eps = _ip_data[ip].eps;
411 auto const& eps_prev = _ip_data[ip].eps_prev;
413 auto& eps_m = _ip_data[ip].eps_m;
414 auto const& eps_m_prev = _ip_data[ip].eps_m_prev;
416 auto& state = _ip_data[ip].material_state_variables;
418 const double T_ip = N.dot(local_T);
420 double const dT_ip = T_ip - N.dot(local_T_prev);
430 eps.noalias() = B * local_u;
434 auto const solid_linear_thermal_expansivity_vector =
439 .value(variables, x_position, t, dt));
442 dthermal_strain = solid_linear_thermal_expansivity_vector * dT_ip;
444 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
457 auto&& solution = _ip_data[ip].solid_material.integrateStress(
458 variables_prev, variables, t, x_position, dt, *state);
462 OGS_FATAL(
"Computation of local constitutive relation failed.");
466 std::tie(sigma, state, C) = std::move(*solution);
468 local_Jac.noalias() += B.transpose() * C * B * w;
470 typename ShapeMatricesType::template MatrixType<DisplacementDim,
472 N_u = ShapeMatricesType::template MatrixType<
473 DisplacementDim, displacement_size>::Zero(DisplacementDim,
476 for (
int i = 0; i < DisplacementDim; ++i)
478 N_u.template block<1, displacement_size / DisplacementDim>(
479 i, i * displacement_size / DisplacementDim)
484 solid_phase.property(MPL::PropertyType::density)
485 .template value<double>(variables, x_position, t, dt);
487 auto const& b = _process_data.specific_body_force;
488 local_rhs.noalias() -=
489 (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
496 const double t,
double const dt, Eigen::VectorXd
const& local_x,
497 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
498 std::vector<double>& local_Jac_data)
501 local_x.template segment<temperature_size>(temperature_index);
503 auto const local_T_prev =
504 local_x_prev.template segment<temperature_size>(temperature_index);
507 typename ShapeMatricesType::template MatrixType<temperature_size,
509 local_Jac_data, temperature_size, temperature_size);
512 typename ShapeMatricesType::template VectorType<temperature_size>>(
513 local_b_data, temperature_size);
516 mass.setZero(temperature_size, temperature_size);
519 laplace.setZero(temperature_size, temperature_size);
521 unsigned const n_integration_points =
522 _integration_method.getNumberOfPoints();
524 auto const& medium = _process_data.media_map.getMedium(_element.getID());
525 auto const& solid_phase = medium->phase(
"Solid");
528 for (
unsigned ip = 0; ip < n_integration_points; ip++)
530 auto const& w = _ip_data[ip].integration_weight;
531 auto const& N = _ip_data[ip].N;
532 auto const& dNdx = _ip_data[ip].dNdx;
535 std::nullopt, _element.getID(),
540 const double T_ip = N.dot(local_T);
544 solid_phase.property(MPL::PropertyType::density)
545 .template value<double>(variables, x_position, t, dt);
547 solid_phase.property(MPL::PropertyType::specific_heat_capacity)
548 .template value<double>(variables, x_position, t, dt);
550 mass.noalias() += N.transpose() * rho_s * c_p * N * w;
556 .value(variables, x_position, t, dt);
561 laplace.noalias() += dNdx.transpose() * thermal_conductivity * dNdx * w;
563 local_Jac.noalias() += laplace + mass / dt;
565 local_rhs.noalias() -=
566 laplace * local_T + mass * (local_T - local_T_prev) / dt;