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;
177 ShapeFunction::NPOINTS,
179 dNdx, N, x_coord, _is_axially_symmetric);
181 auto& sigma = _ip_data[ip].sigma;
182 auto const& sigma_prev = _ip_data[ip].sigma_prev;
184 auto& eps = _ip_data[ip].eps;
185 auto const& eps_prev = _ip_data[ip].eps_prev;
187 auto& eps_m = _ip_data[ip].eps_m;
188 auto const& eps_m_prev = _ip_data[ip].eps_m_prev;
190 auto& state = _ip_data[ip].material_state_variables;
192 const double T_ip = N.dot(T);
193 double const T_prev_ip = N.dot(T_prev);
196 auto const solid_linear_thermal_expansivity_vector =
201 .value(variables, x_position, t, dt));
205 solid_linear_thermal_expansivity_vector * (T_ip - T_prev_ip);
215 eps.noalias() = B * u;
218 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
232 auto&& solution = _ip_data[ip].solid_material.integrateStress(
233 variables_prev, variables, t, x_position, dt, *state);
237 OGS_FATAL(
"Computation of local constitutive relation failed.");
241 std::tie(sigma, state, C) = std::move(*solution);
244 .template block<displacement_size, displacement_size>(
245 displacement_index, displacement_index)
246 .noalias() += B.transpose() * C * B * w;
248 typename ShapeMatricesType::template MatrixType<DisplacementDim,
250 N_u = ShapeMatricesType::template MatrixType<
251 DisplacementDim, displacement_size>::Zero(DisplacementDim,
254 for (
int i = 0; i < DisplacementDim; ++i)
256 N_u.template block<1, displacement_size / DisplacementDim>(
257 i, i * displacement_size / DisplacementDim)
262 solid_phase.property(MPL::PropertyType::density)
263 .template value<double>(variables, x_position, t, dt);
265 auto const& b = _process_data.specific_body_force;
266 local_rhs.template block<displacement_size, 1>(displacement_index, 0)
268 (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
273 auto const alpha_T_tensor =
275 solid_linear_thermal_expansivity_vector);
276 KuT.noalias() += B.transpose() * (C * alpha_T_tensor) * N * w;
278 if (_ip_data[ip].solid_material.getConstitutiveModel() ==
281 auto const s = Invariants::deviatoric_projection * sigma;
282 double const norm_s = Invariants::FrobeniusNorm(s);
283 const double creep_coefficient =
284 _ip_data[ip].solid_material.getTemperatureRelatedCoefficient(
285 t, dt, x_position, T_ip, norm_s);
286 KuT.noalias() += creep_coefficient * B.transpose() * s * N * w;
296 .value(variables, x_position, t, dt);
301 KTT.noalias() += dNdx.transpose() * thermal_conductivity * dNdx * w;
307 .template value<double>(variables, x_position, t, dt);
308 DTT.noalias() += N.transpose() * rho_s * c * N * w;
313 .template block<temperature_size, temperature_size>(temperature_index,
315 .noalias() += KTT + DTT / dt;
319 .template block<displacement_size, temperature_size>(displacement_index,
323 local_rhs.template block<temperature_size, 1>(temperature_index, 0)
324 .noalias() -= KTT * T + DTT * (T - T_prev) / dt;
352 const double t,
double const dt, Eigen::VectorXd
const& local_x,
353 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
354 std::vector<double>& local_Jac_data)
357 local_x.template segment<temperature_size>(temperature_index);
359 auto const local_T_prev =
360 local_x_prev.template segment<temperature_size>(temperature_index);
363 local_x.template segment<displacement_size>(displacement_index);
364 bool const is_u_non_zero = local_u.norm() > 0.0;
367 typename ShapeMatricesType::template MatrixType<displacement_size,
369 local_Jac_data, displacement_size, displacement_size);
372 typename ShapeMatricesType::template VectorType<displacement_size>>(
373 local_b_data, displacement_size);
375 unsigned const n_integration_points =
376 _integration_method.getNumberOfPoints();
382 auto const& medium = _process_data.media_map.getMedium(_element.getID());
383 auto const& solid_phase = medium->phase(
"Solid");
385 for (
unsigned ip = 0; ip < n_integration_points; ip++)
387 auto const& w = _ip_data[ip].integration_weight;
388 auto const& N = _ip_data[ip].N;
389 auto const& dNdx = _ip_data[ip].dNdx;
396 ShapeFunction::NPOINTS,
398 dNdx, N, x_coord, _is_axially_symmetric);
400 auto& sigma = _ip_data[ip].sigma;
401 auto const& sigma_prev = _ip_data[ip].sigma_prev;
403 auto& eps = _ip_data[ip].eps;
404 auto const& eps_prev = _ip_data[ip].eps_prev;
406 auto& eps_m = _ip_data[ip].eps_m;
407 auto const& eps_m_prev = _ip_data[ip].eps_m_prev;
409 auto& state = _ip_data[ip].material_state_variables;
411 const double T_ip = N.dot(local_T);
413 double const dT_ip = T_ip - N.dot(local_T_prev);
423 eps.noalias() = B * local_u;
427 auto const solid_linear_thermal_expansivity_vector =
432 .value(variables, x_position, t, dt));
435 dthermal_strain = solid_linear_thermal_expansivity_vector * dT_ip;
437 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
450 auto&& solution = _ip_data[ip].solid_material.integrateStress(
451 variables_prev, variables, t, x_position, dt, *state);
455 OGS_FATAL(
"Computation of local constitutive relation failed.");
459 std::tie(sigma, state, C) = std::move(*solution);
461 local_Jac.noalias() += B.transpose() * C * B * w;
463 typename ShapeMatricesType::template MatrixType<DisplacementDim,
465 N_u = ShapeMatricesType::template MatrixType<
466 DisplacementDim, displacement_size>::Zero(DisplacementDim,
469 for (
int i = 0; i < DisplacementDim; ++i)
471 N_u.template block<1, displacement_size / DisplacementDim>(
472 i, i * displacement_size / DisplacementDim)
477 solid_phase.property(MPL::PropertyType::density)
478 .template value<double>(variables, x_position, t, dt);
480 auto const& b = _process_data.specific_body_force;
481 local_rhs.noalias() -=
482 (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
489 const double t,
double const dt, Eigen::VectorXd
const& local_x,
490 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
491 std::vector<double>& local_Jac_data)
494 local_x.template segment<temperature_size>(temperature_index);
496 auto const local_T_prev =
497 local_x_prev.template segment<temperature_size>(temperature_index);
500 typename ShapeMatricesType::template MatrixType<temperature_size,
502 local_Jac_data, temperature_size, temperature_size);
505 typename ShapeMatricesType::template VectorType<temperature_size>>(
506 local_b_data, temperature_size);
509 mass.setZero(temperature_size, temperature_size);
512 laplace.setZero(temperature_size, temperature_size);
514 unsigned const n_integration_points =
515 _integration_method.getNumberOfPoints();
519 auto const& medium = _process_data.media_map.getMedium(_element.getID());
520 auto const& solid_phase = medium->phase(
"Solid");
523 for (
unsigned ip = 0; ip < n_integration_points; ip++)
525 auto const& w = _ip_data[ip].integration_weight;
526 auto const& N = _ip_data[ip].N;
527 auto const& dNdx = _ip_data[ip].dNdx;
529 const double T_ip = N.dot(local_T);
533 solid_phase.property(MPL::PropertyType::density)
534 .template value<double>(variables, x_position, t, dt);
536 solid_phase.property(MPL::PropertyType::specific_heat_capacity)
537 .template value<double>(variables, x_position, t, dt);
539 mass.noalias() += N.transpose() * rho_s * c_p * N * w;
545 .value(variables, x_position, t, dt);
550 laplace.noalias() += dNdx.transpose() * thermal_conductivity * dNdx * w;
552 local_Jac.noalias() += laplace + mass / dt;
554 local_rhs.noalias() -=
555 laplace * local_T + mass * (local_T - local_T_prev) / dt;