117 std::vector<double>
const& local_x,
118 std::vector<double>
const& local_x_prev,
119 std::vector<double>& ,
120 std::vector<double>& ,
121 std::vector<double>& local_rhs_data,
122 std::vector<double>& local_Jac_data)
124 auto const local_matrix_size = local_x.size();
125 assert(local_matrix_size == temperature_size + displacement_size);
127 auto T = Eigen::Map<
typename ShapeMatricesType::template VectorType<
128 temperature_size>
const>(local_x.data() + temperature_index,
131 auto u = Eigen::Map<
typename ShapeMatricesType::template VectorType<
132 displacement_size>
const>(local_x.data() + displacement_index,
134 bool const is_u_non_zero = u.norm() > 0.0;
136 auto T_prev = Eigen::Map<
typename ShapeMatricesType::template VectorType<
137 temperature_size>
const>(local_x_prev.data() + temperature_index,
140 auto local_Jac = MathLib::createZeroedMatrix<JacobianMatrix>(
141 local_Jac_data, local_matrix_size, local_matrix_size);
143 auto local_rhs = MathLib::createZeroedVector<RhsVector>(local_rhs_data,
146 typename ShapeMatricesType::template MatrixType<displacement_size,
149 KuT.setZero(displacement_size, temperature_size);
152 KTT.setZero(temperature_size, temperature_size);
155 DTT.setZero(temperature_size, temperature_size);
157 unsigned const n_integration_points =
158 _integration_method.getNumberOfPoints();
165 auto const& medium = _process_data.media_map.getMedium(_element.getID());
166 auto const& solid_phase = medium->phase(
"Solid");
168 for (
unsigned ip = 0; ip < n_integration_points; ip++)
171 auto const& w = _ip_data[ip].integration_weight;
172 auto const& N = _ip_data[ip].N;
173 auto const& dNdx = _ip_data[ip].dNdx;
176 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
180 ShapeFunction::NPOINTS,
182 dNdx, N, x_coord, _is_axially_symmetric);
184 auto& sigma = _ip_data[ip].sigma;
185 auto const& sigma_prev = _ip_data[ip].sigma_prev;
187 auto& eps = _ip_data[ip].eps;
188 auto const& eps_prev = _ip_data[ip].eps_prev;
190 auto& eps_m = _ip_data[ip].eps_m;
191 auto const& eps_m_prev = _ip_data[ip].eps_m_prev;
193 auto& state = _ip_data[ip].material_state_variables;
195 const double T_ip = N.dot(T);
196 double const T_prev_ip = N.dot(T_prev);
199 auto const solid_linear_thermal_expansivity_vector =
200 MPL::formKelvinVector<DisplacementDim>(
204 .value(variables, x_position, t, dt));
208 solid_linear_thermal_expansivity_vector * (T_ip - T_prev_ip);
218 eps.noalias() = B * u;
221 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
235 auto&& solution = _ip_data[ip].solid_material.integrateStress(
236 variables_prev, variables, t, x_position, dt, *state);
240 OGS_FATAL(
"Computation of local constitutive relation failed.");
244 std::tie(sigma, state, C) = std::move(*solution);
247 .template block<displacement_size, displacement_size>(
248 displacement_index, displacement_index)
249 .noalias() += B.transpose() * C * B * w;
251 typename ShapeMatricesType::template MatrixType<DisplacementDim,
253 N_u = ShapeMatricesType::template MatrixType<
254 DisplacementDim, displacement_size>::Zero(DisplacementDim,
257 for (
int i = 0; i < DisplacementDim; ++i)
259 N_u.template block<1, displacement_size / DisplacementDim>(
260 i, i * displacement_size / DisplacementDim)
265 solid_phase.property(MPL::PropertyType::density)
266 .template value<double>(variables, x_position, t, dt);
268 auto const& b = _process_data.specific_body_force;
269 local_rhs.template block<displacement_size, 1>(displacement_index, 0)
271 (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
276 auto const alpha_T_tensor =
278 solid_linear_thermal_expansivity_vector);
279 KuT.noalias() += B.transpose() * (C * alpha_T_tensor) * N * w;
281 if (_ip_data[ip].solid_material.getConstitutiveModel() ==
284 auto const s = Invariants::deviatoric_projection * sigma;
285 double const norm_s = Invariants::FrobeniusNorm(s);
286 const double creep_coefficient =
287 _ip_data[ip].solid_material.getTemperatureRelatedCoefficient(
288 t, dt, x_position, T_ip, norm_s);
289 KuT.noalias() += creep_coefficient * B.transpose() * s * N * w;
299 .value(variables, x_position, t, dt);
302 MaterialPropertyLib::formEigenTensor<DisplacementDim>(lambda);
304 KTT.noalias() += dNdx.transpose() * thermal_conductivity * dNdx * w;
310 .template value<double>(variables, x_position, t, dt);
311 DTT.noalias() += N.transpose() * rho_s * c * N * w;
316 .template block<temperature_size, temperature_size>(temperature_index,
318 .noalias() += KTT + DTT / dt;
322 .template block<displacement_size, temperature_size>(displacement_index,
326 local_rhs.template block<temperature_size, 1>(temperature_index, 0)
327 .noalias() -= KTT * T + DTT * (T - T_prev) / dt;
355 const double t,
double const dt, Eigen::VectorXd
const& local_x,
356 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
357 std::vector<double>& local_Jac_data)
360 local_x.template segment<temperature_size>(temperature_index);
362 auto const local_T_prev =
363 local_x_prev.template segment<temperature_size>(temperature_index);
366 local_x.template segment<displacement_size>(displacement_index);
367 bool const is_u_non_zero = local_u.norm() > 0.0;
370 typename ShapeMatricesType::template MatrixType<displacement_size,
372 local_Jac_data, displacement_size, displacement_size);
375 typename ShapeMatricesType::template VectorType<displacement_size>>(
376 local_b_data, displacement_size);
378 unsigned const n_integration_points =
379 _integration_method.getNumberOfPoints();
385 auto const& medium = _process_data.media_map.getMedium(_element.getID());
386 auto const& solid_phase = medium->phase(
"Solid");
388 for (
unsigned ip = 0; ip < n_integration_points; ip++)
391 auto const& w = _ip_data[ip].integration_weight;
392 auto const& N = _ip_data[ip].N;
393 auto const& dNdx = _ip_data[ip].dNdx;
396 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
400 ShapeFunction::NPOINTS,
402 dNdx, N, x_coord, _is_axially_symmetric);
404 auto& sigma = _ip_data[ip].sigma;
405 auto const& sigma_prev = _ip_data[ip].sigma_prev;
407 auto& eps = _ip_data[ip].eps;
408 auto const& eps_prev = _ip_data[ip].eps_prev;
410 auto& eps_m = _ip_data[ip].eps_m;
411 auto const& eps_m_prev = _ip_data[ip].eps_m_prev;
413 auto& state = _ip_data[ip].material_state_variables;
415 const double T_ip = N.dot(local_T);
417 double const dT_ip = T_ip - N.dot(local_T_prev);
427 eps.noalias() = B * local_u;
431 auto const solid_linear_thermal_expansivity_vector =
432 MPL::formKelvinVector<DisplacementDim>(
436 .value(variables, x_position, t, dt));
439 dthermal_strain = solid_linear_thermal_expansivity_vector * dT_ip;
441 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
454 auto&& solution = _ip_data[ip].solid_material.integrateStress(
455 variables_prev, variables, t, x_position, dt, *state);
459 OGS_FATAL(
"Computation of local constitutive relation failed.");
463 std::tie(sigma, state, C) = std::move(*solution);
465 local_Jac.noalias() += B.transpose() * C * B * w;
467 typename ShapeMatricesType::template MatrixType<DisplacementDim,
469 N_u = ShapeMatricesType::template MatrixType<
470 DisplacementDim, displacement_size>::Zero(DisplacementDim,
473 for (
int i = 0; i < DisplacementDim; ++i)
475 N_u.template block<1, displacement_size / DisplacementDim>(
476 i, i * displacement_size / DisplacementDim)
481 solid_phase.property(MPL::PropertyType::density)
482 .template value<double>(variables, x_position, t, dt);
484 auto const& b = _process_data.specific_body_force;
485 local_rhs.noalias() -=
486 (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
493 const double t,
double const dt, Eigen::VectorXd
const& local_x,
494 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
495 std::vector<double>& local_Jac_data)
498 local_x.template segment<temperature_size>(temperature_index);
500 auto const local_T_prev =
501 local_x_prev.template segment<temperature_size>(temperature_index);
504 typename ShapeMatricesType::template MatrixType<temperature_size,
506 local_Jac_data, temperature_size, temperature_size);
509 typename ShapeMatricesType::template VectorType<temperature_size>>(
510 local_b_data, temperature_size);
513 mass.setZero(temperature_size, temperature_size);
516 laplace.setZero(temperature_size, temperature_size);
518 unsigned const n_integration_points =
519 _integration_method.getNumberOfPoints();
523 auto const& medium = _process_data.media_map.getMedium(_element.getID());
524 auto const& solid_phase = medium->phase(
"Solid");
527 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;
534 const double T_ip = N.dot(local_T);
538 solid_phase.property(MPL::PropertyType::density)
539 .template value<double>(variables, x_position, t, dt);
541 solid_phase.property(MPL::PropertyType::specific_heat_capacity)
542 .template value<double>(variables, x_position, t, dt);
544 mass.noalias() += N.transpose() * rho_s * c_p * N * w;
550 .value(variables, x_position, t, dt);
553 MaterialPropertyLib::formEigenTensor<DisplacementDim>(lambda);
555 laplace.noalias() += dNdx.transpose() * thermal_conductivity * dNdx * w;
557 local_Jac.noalias() += laplace + mass / dt;
559 local_rhs.noalias() -=
560 laplace * local_T + mass * (local_T - local_T_prev) / dt;