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,
138 auto local_Jac = MathLib::createZeroedMatrix<JacobianMatrix>(
139 local_Jac_data, local_matrix_size, local_matrix_size);
141 auto local_rhs = MathLib::createZeroedVector<RhsVector>(local_rhs_data,
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++)
169 auto const& w = _ip_data[ip].integration_weight;
170 auto const& N = _ip_data[ip].N;
171 auto const& dNdx = _ip_data[ip].dNdx;
174 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
178 ShapeFunction::NPOINTS,
180 dNdx, N, x_coord, _is_axially_symmetric);
182 auto& sigma = _ip_data[ip].sigma;
183 auto const& sigma_prev = _ip_data[ip].sigma_prev;
185 auto& eps = _ip_data[ip].eps;
186 auto const& eps_prev = _ip_data[ip].eps_prev;
188 auto& eps_m = _ip_data[ip].eps_m;
189 auto const& eps_m_prev = _ip_data[ip].eps_m_prev;
191 auto& state = _ip_data[ip].material_state_variables;
193 const double T_ip = N.dot(T);
194 double const T_prev_ip = N.dot(T_prev);
197 auto const solid_linear_thermal_expansivity_vector =
198 MPL::formKelvinVector<DisplacementDim>(
202 .value(variables, x_position, t, dt));
206 solid_linear_thermal_expansivity_vector * (T_ip - T_prev_ip);
216 eps.noalias() = B * u;
219 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
233 auto&& solution = _ip_data[ip].solid_material.integrateStress(
234 variables_prev, variables, t, x_position, dt, *state);
238 OGS_FATAL(
"Computation of local constitutive relation failed.");
242 std::tie(sigma, state, C) = std::move(*solution);
245 .template block<displacement_size, displacement_size>(
246 displacement_index, displacement_index)
247 .noalias() += B.transpose() * C * B * w;
249 typename ShapeMatricesType::template MatrixType<DisplacementDim,
251 N_u = ShapeMatricesType::template MatrixType<
252 DisplacementDim, displacement_size>::Zero(DisplacementDim,
255 for (
int i = 0; i < DisplacementDim; ++i)
257 N_u.template block<1, displacement_size / DisplacementDim>(
258 i, i * displacement_size / DisplacementDim)
263 solid_phase.property(MPL::PropertyType::density)
264 .template value<double>(variables, x_position, t, dt);
266 auto const& b = _process_data.specific_body_force;
267 local_rhs.template block<displacement_size, 1>(displacement_index, 0)
269 (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
274 auto const alpha_T_tensor =
276 solid_linear_thermal_expansivity_vector);
277 KuT.noalias() += B.transpose() * (C * alpha_T_tensor) * N * w;
279 if (_ip_data[ip].solid_material.getConstitutiveModel() ==
282 auto const s = Invariants::deviatoric_projection * sigma;
283 double const norm_s = Invariants::FrobeniusNorm(s);
284 const double creep_coefficient =
285 _ip_data[ip].solid_material.getTemperatureRelatedCoefficient(
286 t, dt, x_position, T_ip, norm_s);
287 KuT.noalias() += creep_coefficient * B.transpose() * s * N * w;
297 .value(variables, x_position, t, dt);
300 MaterialPropertyLib::formEigenTensor<DisplacementDim>(lambda);
302 KTT.noalias() += dNdx.transpose() * thermal_conductivity * dNdx * w;
308 .template value<double>(variables, x_position, t, dt);
309 DTT.noalias() += N.transpose() * rho_s * c * N * w;
314 .template block<temperature_size, temperature_size>(temperature_index,
316 .noalias() += KTT + DTT / dt;
320 .template block<displacement_size, temperature_size>(displacement_index,
324 local_rhs.template block<temperature_size, 1>(temperature_index, 0)
325 .noalias() -= KTT * T + DTT * (T - T_prev) / dt;
353 const double t,
double const dt, Eigen::VectorXd
const& local_x,
354 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
355 std::vector<double>& local_Jac_data)
358 local_x.template segment<temperature_size>(temperature_index);
360 auto const local_T_prev =
361 local_x_prev.template segment<temperature_size>(temperature_index);
364 local_x.template segment<displacement_size>(displacement_index);
365 bool const is_u_non_zero = local_u.norm() > 0.0;
368 typename ShapeMatricesType::template MatrixType<displacement_size,
370 local_Jac_data, displacement_size, displacement_size);
373 typename ShapeMatricesType::template VectorType<displacement_size>>(
374 local_b_data, displacement_size);
376 unsigned const n_integration_points =
377 _integration_method.getNumberOfPoints();
383 auto const& medium = _process_data.media_map.getMedium(_element.getID());
384 auto const& solid_phase = medium->phase(
"Solid");
386 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 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
398 ShapeFunction::NPOINTS,
400 dNdx, N, x_coord, _is_axially_symmetric);
402 auto& sigma = _ip_data[ip].sigma;
403 auto const& sigma_prev = _ip_data[ip].sigma_prev;
405 auto& eps = _ip_data[ip].eps;
406 auto const& eps_prev = _ip_data[ip].eps_prev;
408 auto& eps_m = _ip_data[ip].eps_m;
409 auto const& eps_m_prev = _ip_data[ip].eps_m_prev;
411 auto& state = _ip_data[ip].material_state_variables;
413 const double T_ip = N.dot(local_T);
415 double const dT_ip = T_ip - N.dot(local_T_prev);
425 eps.noalias() = B * local_u;
429 auto const solid_linear_thermal_expansivity_vector =
430 MPL::formKelvinVector<DisplacementDim>(
434 .value(variables, x_position, t, dt));
437 dthermal_strain = solid_linear_thermal_expansivity_vector * dT_ip;
439 eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
452 auto&& solution = _ip_data[ip].solid_material.integrateStress(
453 variables_prev, variables, t, x_position, dt, *state);
457 OGS_FATAL(
"Computation of local constitutive relation failed.");
461 std::tie(sigma, state, C) = std::move(*solution);
463 local_Jac.noalias() += B.transpose() * C * B * w;
465 typename ShapeMatricesType::template MatrixType<DisplacementDim,
467 N_u = ShapeMatricesType::template MatrixType<
468 DisplacementDim, displacement_size>::Zero(DisplacementDim,
471 for (
int i = 0; i < DisplacementDim; ++i)
473 N_u.template block<1, displacement_size / DisplacementDim>(
474 i, i * displacement_size / DisplacementDim)
479 solid_phase.property(MPL::PropertyType::density)
480 .template value<double>(variables, x_position, t, dt);
482 auto const& b = _process_data.specific_body_force;
483 local_rhs.noalias() -=
484 (B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
491 const double t,
double const dt, Eigen::VectorXd
const& local_x,
492 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
493 std::vector<double>& local_Jac_data)
496 local_x.template segment<temperature_size>(temperature_index);
498 auto const local_T_prev =
499 local_x_prev.template segment<temperature_size>(temperature_index);
502 typename ShapeMatricesType::template MatrixType<temperature_size,
504 local_Jac_data, temperature_size, temperature_size);
507 typename ShapeMatricesType::template VectorType<temperature_size>>(
508 local_b_data, temperature_size);
511 mass.setZero(temperature_size, temperature_size);
514 laplace.setZero(temperature_size, temperature_size);
516 unsigned const n_integration_points =
517 _integration_method.getNumberOfPoints();
521 auto const& medium = _process_data.media_map.getMedium(_element.getID());
522 auto const& solid_phase = medium->phase(
"Solid");
525 for (
unsigned ip = 0; ip < n_integration_points; ip++)
528 auto const& w = _ip_data[ip].integration_weight;
529 auto const& N = _ip_data[ip].N;
530 auto const& dNdx = _ip_data[ip].dNdx;
532 const double T_ip = N.dot(local_T);
536 solid_phase.property(MPL::PropertyType::density)
537 .template value<double>(variables, x_position, t, dt);
539 solid_phase.property(MPL::PropertyType::specific_heat_capacity)
540 .template value<double>(variables, x_position, t, dt);
542 mass.noalias() += N.transpose() * rho_s * c_p * N * w;
548 .value(variables, x_position, t, dt);
551 MaterialPropertyLib::formEigenTensor<DisplacementDim>(lambda);
553 laplace.noalias() += dNdx.transpose() * thermal_conductivity * dNdx * w;
555 local_Jac.noalias() += laplace + mass / dt;
557 local_rhs.noalias() -=
558 laplace * local_T + mass * (local_T - local_T_prev) / dt;