39namespace ComponentTransport
41template <
typename GlobalDimNodalMatrixType>
45 double const& integration_weight_)
51 GlobalDimNodalMatrixType
const dNdx;
58 double porosity = std::numeric_limits<double>::quiet_NaN();
73 std::size_t
const mesh_item_id,
74 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
75 std::vector<GlobalVector*>
const& x,
double const t)
77 std::vector<double> local_x_vec;
79 auto const n_processes = x.size();
80 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
84 assert(!indices.empty());
85 auto const local_solution = x[process_id]->get(indices);
86 local_x_vec.insert(std::end(local_x_vec),
87 std::begin(local_solution),
88 std::end(local_solution));
96 std::size_t
const mesh_item_id,
97 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
98 std::vector<GlobalVector*>
const& x,
double const t,
double const dt)
100 std::vector<double> local_x_vec;
102 auto const n_processes = x.size();
103 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
107 assert(!indices.empty());
108 auto const local_solution = x[process_id]->get(indices);
109 local_x_vec.insert(std::end(local_x_vec),
110 std::begin(local_solution),
111 std::end(local_solution));
119 std::size_t
const mesh_item_id,
120 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
121 std::vector<GlobalVector*>
const& x,
double const t,
double const dt,
124 std::vector<double> local_x_vec;
126 auto const n_processes = x.size();
127 for (std::size_t pcs_id = 0; pcs_id < n_processes; ++pcs_id)
131 assert(!indices.empty());
132 auto const local_solution = x[pcs_id]->get(indices);
133 local_x_vec.insert(std::end(local_x_vec),
134 std::begin(local_solution),
135 std::end(local_solution));
141 auto const num_r_c = indices.size();
143 std::vector<double> local_M_data;
144 local_M_data.reserve(num_r_c * num_r_c);
145 std::vector<double> local_K_data;
146 local_K_data.reserve(num_r_c * num_r_c);
147 std::vector<double> local_b_data;
148 local_b_data.reserve(num_r_c);
151 local_K_data, local_b_data,
154 auto const r_c_indices =
156 if (!local_M_data.empty())
160 M.
add(r_c_indices, local_M);
162 if (!local_K_data.empty())
166 K.
add(r_c_indices, local_K);
168 if (!local_b_data.empty())
170 b.
add(indices, local_b_data);
175 double const t,
double const dt) = 0;
178 std::size_t
const ele_id) = 0;
182 std::vector<GlobalVector*>
const& x,
183 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
184 std::vector<double>& cache)
const = 0;
187 const double t, std::vector<GlobalVector*>
const& x,
188 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
189 std::vector<double>& cache,
int const component_id)
const = 0;
193 Eigen::VectorXd
const& ,
double const ) = 0;
200 double const t,
double const dt, Eigen::VectorXd
const& local_x,
201 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
202 std::vector<double>& local_b_data,
int const transport_process_id) = 0;
205template <
typename ShapeFunction,
int GlobalDim>
217 ShapeFunction::NPOINTS;
223 typename ShapeMatricesType::template MatrixType<
pressure_size,
226 typename ShapeMatricesType::template VectorType<pressure_size>;
229 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
243 std::size_t
const local_matrix_size,
245 bool is_axially_symmetric,
247 std::vector<std::reference_wrapper<ProcessVariable>>
const&
248 transport_process_variables)
250 : ShapeFunction::NPOINTS),
252 ? ShapeFunction::NPOINTS
253 : 2 * ShapeFunction::NPOINTS),
259 (void)local_matrix_size;
261 unsigned const n_integration_points =
263 _ip_data.reserve(n_integration_points);
270 auto const shape_matrices =
272 GlobalDim>(element, is_axially_symmetric,
276 for (
unsigned ip = 0; ip < n_integration_points; ip++)
279 shape_matrices[ip].dNdx,
281 shape_matrices[ip].integralMeasure *
282 shape_matrices[ip].detJ * aperture_size);
288 .template initialValue<double>(
289 pos, std::numeric_limits<double>::quiet_NaN() );
299 auto& chemical_system_index_map =
302 unsigned const n_integration_points =
304 for (
unsigned ip = 0; ip < n_integration_points; ip++)
307 chemical_system_index_map.empty()
309 : chemical_system_index_map.back() + 1;
310 chemical_system_index_map.push_back(
316 double const t)
override
330 unsigned const n_integration_points =
333 for (
unsigned ip = 0; ip < n_integration_points; ip++)
336 auto const& N = Ns[ip];
337 auto const& chemical_system_id = ip_data.chemical_system_id;
340 std::vector<double> C_int_pt(n_component);
341 for (
unsigned component_id = 0; component_id < n_component;
344 auto const concentration_index =
348 local_x.template segment<concentration_size>(
349 concentration_index);
352 C_int_pt[component_id]);
362 double const t,
double dt)
override
379 unsigned const n_integration_points =
382 for (
unsigned ip = 0; ip < n_integration_points; ip++)
385 auto const& N = Ns[ip];
386 auto& porosity = ip_data.porosity;
387 auto const& porosity_prev = ip_data.porosity_prev;
388 auto const& chemical_system_id = ip_data.chemical_system_id;
391 std::vector<double> C_int_pt(n_component);
392 for (
unsigned component_id = 0; component_id < n_component;
395 auto const concentration_index =
399 local_x.template segment<concentration_size>(
400 concentration_index);
403 C_int_pt[component_id]);
415 .template value<double>(vars, vars_prev, pos, t,
422 C_int_pt, chemical_system_id, medium, vars, pos, t, dt);
427 double const dt)
override
441 ip_data.porosity = ip_data.porosity_prev;
446 ip_data.porosity, t, dt);
449 ip_data.chemical_system_id, medium, ip_data.porosity);
454 std::vector<double>
const& local_x,
455 std::vector<double>
const& ,
456 std::vector<double>& local_M_data,
457 std::vector<double>& local_K_data,
458 std::vector<double>& local_b_data)
override
460 auto const local_matrix_size = local_x.size();
465 assert(local_matrix_size == ShapeFunction::NPOINTS * num_nodal_dof);
467 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
468 local_M_data, local_matrix_size, local_matrix_size);
469 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
470 local_K_data, local_matrix_size, local_matrix_size);
471 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
472 local_b_data, local_matrix_size);
475 auto Kpp = local_K.template block<pressure_size, pressure_size>(
477 auto Mpp = local_M.template block<pressure_size, pressure_size>(
479 auto Bp = local_b.template segment<pressure_size>(
pressure_index);
481 auto local_p = Eigen::Map<const NodalVectorType>(
488 auto const number_of_components = num_nodal_dof - 1;
489 for (
int component_id = 0; component_id < number_of_components;
501 auto concentration_index =
505 local_K.template block<concentration_size, concentration_size>(
506 concentration_index, concentration_index);
508 local_M.template block<concentration_size, concentration_size>(
509 concentration_index, concentration_index);
511 local_M.template block<concentration_size, pressure_size>(
514 local_M.template block<pressure_size, concentration_size>(
517 auto local_C = Eigen::Map<const NodalVectorType>(
521 MCC, MCp, MpC, Kpp, Mpp, Bp);
525 auto const stoichiometric_matrix =
529 assert(stoichiometric_matrix);
531 for (Eigen::SparseMatrix<double>::InnerIterator it(
532 *stoichiometric_matrix, component_id);
536 auto const stoichiometric_coefficient = it.value();
537 auto const coupled_component_id = it.row();
538 auto const kinetic_prefactor =
542 auto const concentration_index =
544 auto const coupled_concentration_index =
549 concentration_index, coupled_concentration_index);
553 stoichiometric_coefficient,
563 Eigen::Ref<const NodalVectorType>
const& C_nodal_values,
564 Eigen::Ref<const NodalVectorType>
const& p_nodal_values,
565 Eigen::Ref<LocalBlockMatrixType> KCC,
566 Eigen::Ref<LocalBlockMatrixType> MCC,
567 Eigen::Ref<LocalBlockMatrixType> MCp,
568 Eigen::Ref<LocalBlockMatrixType> MpC,
569 Eigen::Ref<LocalBlockMatrixType> Kpp,
570 Eigen::Ref<LocalBlockMatrixType> Mpp,
571 Eigen::Ref<LocalSegmentVectorType> Bp)
573 unsigned const n_integration_points =
585 auto const& phase = medium.
phase(
"AqueousLiquid");
596 std::vector<GlobalDimVectorType> ip_flux_vector;
597 double average_velocity_norm = 0.0;
600 ip_flux_vector.reserve(n_integration_points);
607 for (
unsigned ip(0); ip < n_integration_points; ++ip)
612 auto const& dNdx = ip_data.dNdx;
613 auto const& N = Ns[ip];
614 auto const& w = ip_data.integration_weight;
615 auto& porosity = ip_data.porosity;
617 double C_int_pt = 0.0;
618 double p_int_pt = 0.0;
628 .template value<double>(vars, pos, t, dt);
631 auto const& retardation_factor =
633 .template value<double>(vars, pos, t, dt);
635 auto const& solute_dispersivity_transverse = medium.template value<
639 auto const& solute_dispersivity_longitudinal =
640 medium.template value<double>(
642 longitudinal_dispersivity);
649 .template value<double>(vars, pos, t, dt);
651 auto const decay_rate =
653 .template value<double>(vars, pos, t, dt);
655 auto const& pore_diffusion_coefficient =
656 MaterialPropertyLib::formEigenTensor<GlobalDim>(
658 .value(vars, pos, t, dt));
660 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
666 .template value<double>(vars, pos, t, dt);
672 (dNdx * p_nodal_values - density * b))
675 const double drho_dp =
677 .template dValue<double>(
682 const double drho_dC =
684 .template dValue<double>(
691 pore_diffusion_coefficient, velocity, porosity,
692 solute_dispersivity_transverse,
693 solute_dispersivity_longitudinal);
695 const double R_times_phi(retardation_factor * porosity);
697 auto const N_t_N = (N.transpose() * N).eval();
701 MCp.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dp * w);
702 MCC.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dC * w);
703 KCC.noalias() -= dNdx.transpose() * mass_density_flow * N * w;
707 ip_flux_vector.emplace_back(mass_density_flow);
708 average_velocity_norm += velocity.norm();
710 MCC.noalias() += N_t_N * (R_times_phi * density * w);
711 KCC.noalias() += N_t_N * (decay_rate * R_times_phi * density * w);
712 KCC_Laplacian.noalias() +=
713 dNdx.transpose() * hydrodynamic_dispersion * dNdx * density * w;
715 MpC.noalias() += N_t_N * (porosity * drho_dC * w);
718 if (component_id == 0)
720 Mpp.noalias() += N_t_N * (porosity * drho_dp * w);
722 dNdx.transpose() * K_over_mu * dNdx * (density * w);
726 Bp.noalias() += dNdx.transpose() * K_over_mu * b *
727 (density * density * w);
735 typename ShapeFunction::MeshElement>(
740 average_velocity_norm /
741 static_cast<double>(n_integration_points),
745 KCC.noalias() += KCC_Laplacian;
749 Eigen::Ref<LocalBlockMatrixType> KCmCn,
750 double const stoichiometric_coefficient,
751 double const kinetic_prefactor)
753 unsigned const n_integration_points =
763 auto const& phase = medium.
phase(
"AqueousLiquid");
771 for (
unsigned ip(0); ip < n_integration_points; ++ip)
774 auto const& w = ip_data.integration_weight;
775 auto const& N = Ns[ip];
776 auto& porosity = ip_data.porosity;
778 auto const retardation_factor =
780 .template value<double>(vars, pos, t, dt);
783 .template value<double>(vars, pos, t, dt);
787 .template value<double>(vars, pos, t, dt);
789 KCmCn.noalias() -= w * N.transpose() * stoichiometric_coefficient *
790 kinetic_prefactor * retardation_factor *
791 porosity * density * N;
796 Eigen::VectorXd
const& local_x,
797 Eigen::VectorXd
const& local_x_prev,
798 int const process_id,
799 std::vector<double>& local_M_data,
800 std::vector<double>& local_K_data,
801 std::vector<double>& local_b_data)
override
806 local_M_data, local_K_data, local_b_data);
811 local_M_data, local_K_data,
818 local_M_data, local_K_data,
819 local_b_data, process_id);
825 Eigen::VectorXd
const& local_x,
826 Eigen::VectorXd
const& local_x_prev,
827 std::vector<double>& local_M_data,
828 std::vector<double>& local_K_data,
829 std::vector<double>& local_b_data)
833 auto const local_C = local_x.template segment<concentration_size>(
835 auto const local_C_prev =
840 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
842 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
844 auto local_b = MathLib::createZeroedVector<LocalSegmentVectorType>(
847 unsigned const n_integration_points =
859 auto const& phase = medium.
phase(
"AqueousLiquid");
868 for (
unsigned ip(0); ip < n_integration_points; ++ip)
873 auto const& dNdx = ip_data.dNdx;
874 auto const& w = ip_data.integration_weight;
875 auto const& N = Ns[ip];
876 auto& porosity = ip_data.porosity;
877 auto const& porosity_prev = ip_data.porosity_prev;
879 double const C_int_pt = N.dot(local_C);
880 double const p_int_pt = N.dot(local_p);
881 double const T_int_pt = N.dot(local_T);
895 .template value<double>(vars, vars_prev, pos, t,
906 .template value<double>(vars, pos, t, dt);
908 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
914 .template value<double>(vars, pos, t, dt);
918 const double drho_dp =
920 .template dValue<double>(
924 const double drho_dC =
926 .template dValue<double>(
931 local_M.noalias() += w * N.transpose() * porosity * drho_dp * N;
933 w * dNdx.transpose() * density * K_over_mu * dNdx;
938 w * density * density * dNdx.transpose() * K_over_mu * b;
943 double const C_dot = (C_int_pt - N.dot(local_C_prev)) / dt;
946 w * N.transpose() * porosity * drho_dC * C_dot;
952 Eigen::VectorXd
const& local_x,
953 Eigen::VectorXd
const& ,
954 std::vector<double>& local_M_data,
955 std::vector<double>& local_K_data,
956 std::vector<double>& )
958 assert(local_x.size() ==
966 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
968 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
977 auto const& liquid_phase = medium.
phase(
"AqueousLiquid");
985 unsigned const n_integration_points =
988 std::vector<GlobalDimVectorType> ip_flux_vector;
989 double average_velocity_norm = 0.0;
990 ip_flux_vector.reserve(n_integration_points);
996 for (
unsigned ip(0); ip < n_integration_points; ip++)
1000 auto const& ip_data = this->
_ip_data[ip];
1001 auto const& dNdx = ip_data.dNdx;
1002 auto const& w = ip_data.integration_weight;
1003 auto const& N = Ns[ip];
1005 double p_at_xi = 0.;
1007 double T_at_xi = 0.;
1015 auto const porosity =
1017 .template value<double>(vars, pos, t, dt);
1021 auto const fluid_density =
1024 .template value<double>(vars, pos, t, dt);
1026 auto const specific_heat_capacity_fluid =
1029 .template value<double>(vars, pos, t, dt);
1032 local_M.noalias() += w *
1034 vars, porosity, fluid_density,
1035 specific_heat_capacity_fluid, pos, t, dt) *
1039 auto const viscosity =
1042 .template value<double>(vars, pos, t, dt);
1044 auto const intrinsic_permeability =
1045 MaterialPropertyLib::formEigenTensor<GlobalDim>(
1049 .value(vars, pos, t, dt));
1052 intrinsic_permeability / viscosity;
1054 process_data.has_gravity
1056 (dNdx * local_p - fluid_density * b))
1061 vars, fluid_density, specific_heat_capacity_fluid, velocity,
1064 local_K.noalias() +=
1065 w * dNdx.transpose() * thermal_conductivity_dispersivity * dNdx;
1067 ip_flux_vector.emplace_back(velocity * fluid_density *
1068 specific_heat_capacity_fluid);
1069 average_velocity_norm += velocity.norm();
1072 NumLib::assembleAdvectionMatrix<typename ShapeFunction::MeshElement>(
1073 process_data.stabilizer, this->_ip_data,
1075 average_velocity_norm /
static_cast<double>(n_integration_points),
1080 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1081 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_M_data,
1082 std::vector<double>& local_K_data,
1083 std::vector<double>& ,
int const transport_process_id)
1085 assert(
static_cast<int>(local_x.size()) ==
1091 auto const local_p =
1096 auto const local_C = local_x.template segment<concentration_size>(
1100 auto const local_p_prev =
1103 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1105 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1111 unsigned const n_integration_points =
1114 std::vector<GlobalDimVectorType> ip_flux_vector;
1115 double average_velocity_norm = 0.0;
1118 ip_flux_vector.reserve(n_integration_points);
1131 auto const& medium =
1133 auto const& phase = medium.
phase(
"AqueousLiquid");
1134 auto const component_id =
1136 auto const& component = phase.
component(
1143 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1148 auto const& dNdx = ip_data.dNdx;
1149 auto const& w = ip_data.integration_weight;
1150 auto const& N = Ns[ip];
1151 auto& porosity = ip_data.porosity;
1152 auto const& porosity_prev = ip_data.porosity_prev;
1154 double const C_int_pt = N.dot(local_C);
1155 double const p_int_pt = N.dot(local_p);
1156 double const T_int_pt = N.dot(local_T);
1169 vars_prev.
porosity = porosity_prev;
1175 .template value<double>(vars, vars_prev, pos, t,
1181 auto const& retardation_factor =
1183 .template value<double>(vars, pos, t, dt);
1185 auto const& solute_dispersivity_transverse = medium.template value<
1188 auto const& solute_dispersivity_longitudinal =
1189 medium.template value<double>(
1191 longitudinal_dispersivity);
1194 auto const density =
1196 .template value<double>(vars, pos, t, dt);
1197 auto const decay_rate =
1199 .template value<double>(vars, pos, t, dt);
1201 auto const& pore_diffusion_coefficient =
1202 MaterialPropertyLib::formEigenTensor<GlobalDim>(
1204 .value(vars, pos, t, dt));
1206 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1211 .template value<double>(vars, pos, t, dt);
1217 (dNdx * local_p - density * b))
1223 pore_diffusion_coefficient, velocity, porosity,
1224 solute_dispersivity_transverse,
1225 solute_dispersivity_longitudinal);
1227 double const R_times_phi = retardation_factor * porosity;
1228 auto const N_t_N = (N.transpose() * N).eval();
1232 const double drho_dC =
1234 .template dValue<double>(
1237 local_M.noalias() +=
1238 N_t_N * (R_times_phi * C_int_pt * drho_dC * w);
1241 local_M.noalias() += N_t_N * (R_times_phi * density * w);
1246 double const p_dot = (p_int_pt - N.dot(local_p_prev)) / dt;
1248 const double drho_dp =
1250 .template dValue<double>(vars,
1252 liquid_phase_pressure,
1255 local_K.noalias() +=
1256 N_t_N * ((R_times_phi * drho_dp * p_dot) * w) -
1257 dNdx.transpose() * velocity * N * (density * w);
1261 ip_flux_vector.emplace_back(velocity * density);
1262 average_velocity_norm += velocity.norm();
1264 local_K.noalias() +=
1265 N_t_N * (decay_rate * R_times_phi * density * w);
1267 KCC_Laplacian.noalias() += dNdx.transpose() *
1268 hydrodynamic_dispersion * dNdx *
1275 typename ShapeFunction::MeshElement>(
1278 average_velocity_norm /
1279 static_cast<double>(n_integration_points),
1282 local_K.noalias() += KCC_Laplacian;
1286 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1287 Eigen::VectorXd
const& local_x_prev,
int const process_id,
1288 std::vector<double>& local_b_data,
1289 std::vector<double>& local_Jac_data)
override
1294 local_b_data, local_Jac_data);
1298 int const component_id = process_id - 1;
1300 t, dt, local_x, local_x_prev, local_b_data, local_Jac_data,
1306 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1307 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
1308 std::vector<double>& local_Jac_data)
1310 auto const p = local_x.template segment<pressure_size>(
pressure_index);
1311 auto const c = local_x.template segment<concentration_size>(
1318 auto local_Jac = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1320 auto local_rhs = MathLib::createZeroedVector<LocalSegmentVectorType>(
1323 unsigned const n_integration_points =
1332 auto const& medium =
1334 auto const& phase = medium.
phase(
"AqueousLiquid");
1343 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1348 auto const& dNdx = ip_data.dNdx;
1349 auto const& w = ip_data.integration_weight;
1350 auto const& N = Ns[ip];
1351 auto& phi = ip_data.porosity;
1352 auto const& phi_prev = ip_data.porosity_prev;
1354 double const p_ip = N.dot(p);
1355 double const c_ip = N.dot(c);
1357 double const cdot_ip = (c_ip - N.dot(c_prev)) / dt;
1369 .template value<double>(vars, vars_prev, pos, t,
1376 .template value<double>(vars, pos, t, dt);
1378 auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1383 .template value<double>(vars, pos, t, dt);
1385 auto const drho_dp =
1387 .template dValue<double>(
1391 auto const drho_dc =
1393 .template dValue<double>(
1398 local_Jac.noalias() += w * N.transpose() * phi * drho_dp / dt * N +
1399 w * dNdx.transpose() * rho * k / mu * dNdx;
1401 local_rhs.noalias() -=
1402 w * N.transpose() * phi *
1403 (drho_dp * N * p_prev + drho_dc * cdot_ip) +
1404 w * rho * dNdx.transpose() * k / mu * dNdx * p;
1408 local_rhs.noalias() +=
1409 w * rho * dNdx.transpose() * k / mu * rho * b;
1415 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1416 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
1417 std::vector<double>& local_Jac_data,
int const component_id)
1419 auto const concentration_index =
1422 auto const p = local_x.template segment<pressure_size>(
pressure_index);
1424 local_x.template segment<concentration_size>(concentration_index);
1434 auto local_Jac = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1436 auto local_rhs = MathLib::createZeroedVector<LocalSegmentVectorType>(
1442 unsigned const n_integration_points =
1445 std::vector<GlobalDimVectorType> ip_flux_vector;
1446 double average_velocity_norm = 0.0;
1447 ip_flux_vector.reserve(n_integration_points);
1459 auto const& medium =
1461 auto const& phase = medium.
phase(
"AqueousLiquid");
1462 auto const& component = phase.
component(
1469 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1474 auto const& dNdx = ip_data.dNdx;
1475 auto const& w = ip_data.integration_weight;
1476 auto const& N = Ns[ip];
1477 auto& phi = ip_data.porosity;
1478 auto const& phi_prev = ip_data.porosity_prev;
1480 double const p_ip = N.dot(p);
1481 double const c_ip = N.dot(c);
1498 .template value<double>(vars, vars_prev, pos, t,
1506 .template value<double>(vars, pos, t, dt);
1508 auto const alpha_T = medium.template value<double>(
1510 auto const alpha_L = medium.template value<double>(
1514 .template value<double>(vars, pos, t, dt);
1518 .template value<double>(vars, pos, t, dt);
1520 auto const Dp = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1522 .value(vars, pos, t, dt));
1524 auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1528 .template value<double>(vars, pos, t, dt);
1540 local_Jac.noalias() +=
1541 w * rho * N.transpose() * phi * R * (alpha + 1 / dt) * N;
1543 KCC_Laplacian.noalias() += w * rho * dNdx.transpose() * D * dNdx;
1545 auto const cdot = (c - c_prev) / dt;
1546 local_rhs.noalias() -=
1547 w * rho * N.transpose() * phi * R * N * (cdot + alpha * c);
1549 ip_flux_vector.emplace_back(q * rho);
1550 average_velocity_norm += q.norm();
1553 NumLib::assembleAdvectionMatrix<typename ShapeFunction::MeshElement>(
1556 average_velocity_norm /
static_cast<double>(n_integration_points),
1559 local_rhs.noalias() -= KCC_Laplacian * c;
1561 local_Jac.noalias() += KCC_Laplacian;
1565 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1566 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
1567 std::vector<double>& local_b_data,
1568 int const transport_process_id)
override
1570 auto const local_C = local_x.template segment<concentration_size>(
1574 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1576 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1578 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
1581 unsigned const n_integration_points =
1590 auto const& medium =
1592 auto const component_id = transport_process_id - 1;
1598 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1603 auto const w = ip_data.integration_weight;
1604 auto const& N = Ns[ip];
1605 auto& porosity = ip_data.porosity;
1606 auto const& porosity_prev = ip_data.porosity_prev;
1607 auto const chemical_system_id = ip_data.chemical_system_id;
1609 double C_int_pt = 0.0;
1614 auto const porosity_dot = (porosity - porosity_prev) / dt;
1618 vars_prev.
porosity = porosity_prev;
1624 .template value<double>(vars, vars_prev, pos, t,
1628 local_M.noalias() += w * N.transpose() * porosity * N;
1630 local_K.noalias() += w * N.transpose() * porosity_dot * N;
1632 if (chemical_system_id == -1)
1637 auto const C_post_int_pt =
1639 component_id, chemical_system_id);
1641 local_b.noalias() +=
1642 w * N.transpose() * porosity * (C_post_int_pt - C_int_pt) / dt;
1648 std::vector<GlobalVector*>
const& x,
1649 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
1650 std::vector<double>& cache)
const override
1652 assert(x.size() == dof_table.size());
1654 auto const n_processes = x.size();
1655 std::vector<std::vector<double>> local_x;
1656 local_x.reserve(n_processes);
1658 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1660 auto const indices =
1662 assert(!indices.empty());
1663 local_x.push_back(x[process_id]->get(indices));
1667 if (n_processes == 1)
1669 auto const local_p = Eigen::Map<const NodalVectorType>(
1671 auto const local_C = Eigen::Map<const NodalVectorType>(
1678 constexpr int pressure_process_id = 0;
1679 constexpr int concentration_process_id = 1;
1680 auto const local_p = Eigen::Map<const NodalVectorType>(
1682 auto const local_C = Eigen::Map<const NodalVectorType>(
1690 Eigen::Ref<const NodalVectorType>
const& p_nodal_values,
1691 Eigen::Ref<const NodalVectorType>
const& C_nodal_values,
1692 std::vector<double>& cache)
const
1694 auto const n_integration_points =
1699 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1700 cache, GlobalDim, n_integration_points);
1711 auto const& medium =
1713 auto const& phase = medium.
phase(
"AqueousLiquid");
1719 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
1721 auto const& ip_data =
_ip_data[ip];
1722 auto const& dNdx = ip_data.dNdx;
1723 auto const& N = Ns[ip];
1724 auto const& porosity = ip_data.porosity;
1726 pos.setIntegrationPoint(ip);
1728 double C_int_pt = 0.0;
1729 double p_int_pt = 0.0;
1740 double const dt = std::numeric_limits<double>::quiet_NaN();
1741 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1745 .template value<double>(vars, pos, t, dt);
1748 cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
1753 .template value<double>(vars, pos, t, dt);
1755 cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
1763 const unsigned integration_point)
const override
1766 typename ShapeFunction::MeshElement>()[integration_point];
1769 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
1774 std::vector<double>
const& local_x)
const override
1776 auto const local_p = Eigen::Map<const NodalVectorType>(
1778 auto const local_C = Eigen::Map<const NodalVectorType>(
1784 auto const shape_matrices =
1788 std::array{pnt_local_coords})[0];
1798 auto const& medium =
1800 auto const& phase = medium.
phase(
"AqueousLiquid");
1813 double const dt = std::numeric_limits<double>::quiet_NaN();
1814 auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1819 .template value<double>(vars, pos, t, dt);
1824 .template value<double>(vars, pos, t, dt);
1827 q += K_over_mu * rho_w * b;
1829 Eigen::Vector3d flux(0.0, 0.0, 0.0);
1830 flux.head<GlobalDim>() = rho_w * q;
1837 Eigen::VectorXd
const& local_x,
1838 Eigen::VectorXd
const& )
override
1840 auto const local_p =
1842 auto const local_C = local_x.template segment<concentration_size>(
1845 std::vector<double> ele_velocity;
1848 auto const n_integration_points =
1850 auto const ele_velocity_mat =
1854 Eigen::Map<LocalVectorType>(
1857 ele_velocity_mat.rowwise().sum() / n_integration_points;
1861 std::size_t
const ele_id)
override
1863 auto const n_integration_points =
1872 ip_data.porosity = ip_data.porosity_prev;
1876 medium, ip_data.porosity);
1881 [](
double const s,
auto const& ip)
1882 {
return s + ip.porosity; }) /
1883 n_integration_points;
1886 std::vector<GlobalIndexType> chemical_system_indices;
1887 chemical_system_indices.reserve(n_integration_points);
1889 std::back_inserter(chemical_system_indices),
1890 [](
auto const& ip_data)
1891 { return ip_data.chemical_system_id; });
1894 ele_id, chemical_system_indices);
1898 const double t, std::vector<GlobalVector*>
const& x,
1899 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
1900 std::vector<double>& cache,
int const component_id)
const override
1902 std::vector<double> local_x_vec;
1904 auto const n_processes = x.size();
1905 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1907 auto const indices =
1909 assert(!indices.empty());
1910 auto const local_solution = x[process_id]->get(indices);
1911 local_x_vec.insert(std::end(local_x_vec),
1912 std::begin(local_solution),
1913 std::end(local_solution));
1917 auto const p = local_x.template segment<pressure_size>(
pressure_index);
1918 auto const c = local_x.template segment<concentration_size>(
1921 auto const n_integration_points =
1926 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1927 cache, GlobalDim, n_integration_points);
1938 auto const& medium =
1940 auto const& phase = medium.
phase(
"AqueousLiquid");
1942 auto const& component = phase.
component(
1949 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
1951 auto const& ip_data =
_ip_data[ip];
1952 auto const& dNdx = ip_data.dNdx;
1953 auto const& N = Ns[ip];
1954 auto const& phi = ip_data.porosity;
1956 pos.setIntegrationPoint(ip);
1958 double const p_ip = N.dot(p);
1959 double const c_ip = N.dot(c);
1965 double const dt = std::numeric_limits<double>::quiet_NaN();
1967 auto const& k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1971 .template value<double>(vars, pos, t, dt);
1973 .template value<double>(vars, pos, t, dt);
1981 auto const alpha_T = medium.template value<double>(
1983 auto const alpha_L = medium.template value<double>(
1985 auto const Dp = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1987 .value(vars, pos, t, dt));
1994 cache_mat.col(ip).noalias() = q * c_ip - D * dNdx * c;
2001 Eigen::VectorXd
const& ,
2002 double const ,
double const ,
2003 int const )
override
2005 unsigned const n_integration_points =
2008 for (
unsigned ip = 0; ip < n_integration_points; ip++)
2019 std::vector<std::reference_wrapper<ProcessVariable>>
const
2022 std::vector<IntegrationPointData<GlobalDimNodalMatrixType>>
_ip_data;
2026 const double fluid_density,
const double specific_heat_capacity_fluid,
2030 auto const& medium =
2032 auto const& solid_phase = medium.
phase(
"Solid");
2034 auto const specific_heat_capacity_solid =
2038 .template value<double>(vars, pos, t, dt);
2040 auto const solid_density =
2042 .template value<double>(vars, pos, t, dt);
2044 return solid_density * specific_heat_capacity_solid * (1 - porosity) +
2045 fluid_density * specific_heat_capacity_fluid * porosity;
2050 const double fluid_density,
const double specific_heat_capacity_fluid,
2055 auto const& medium =
2058 auto thermal_conductivity =
2059 MaterialPropertyLib::formEigenTensor<GlobalDim>(
2063 .value(vars, pos, t, dt));
2065 auto const thermal_dispersivity_transversal =
2068 thermal_transversal_dispersivity)
2069 .
template value<double>();
2071 auto const thermal_dispersivity_longitudinal =
2074 thermal_longitudinal_dispersivity)
2075 .
template value<double>();
2080 return thermal_conductivity +
2081 fluid_density * specific_heat_capacity_fluid *
2084 GlobalDimMatrixType::Zero(GlobalDim, GlobalDim),
2085 velocity, 0 , thermal_dispersivity_transversal,
2086 thermal_dispersivity_longitudinal);
2090 Eigen::VectorXd
const& local_x)
GlobalMatrix::IndexType GlobalIndexType
virtual void updatePorosityPostReaction(GlobalIndexType const &, MaterialPropertyLib::Medium const &, double &)
virtual double getKineticPrefactor(std::size_t reaction_id) const
virtual void computeSecondaryVariable(std::size_t const, std::vector< GlobalIndexType > const &)
virtual void setChemicalSystemConcrete(std::vector< double > const &, GlobalIndexType const &, MaterialPropertyLib::Medium const *, MaterialPropertyLib::VariableArray const &, ParameterLib::SpatialPosition const &, double const, double const)
virtual void updateVolumeFractionPostReaction(GlobalIndexType const &, MaterialPropertyLib::Medium const &, ParameterLib::SpatialPosition const &, double const, double const, double const)
std::vector< GlobalIndexType > chemical_system_index_map
virtual Eigen::SparseMatrix< double > const * getStoichiometricMatrix() const
virtual double getConcentration(int const, GlobalIndexType const) const
virtual void initializeChemicalSystemConcrete(std::vector< double > const &, GlobalIndexType const &, MaterialPropertyLib::Medium const &, ParameterLib::SpatialPosition const &, double const)
Medium * getMedium(std::size_t element_id)
Phase const & phase(std::size_t index) const
Property const & property(PropertyType const &p) const
Component const & component(std::size_t const &index) const
double liquid_phase_pressure
int add(IndexType row, IndexType col, double val)
Global vector based on Eigen vector.
void add(IndexType rowId, double v)
add entry
std::size_t getID() const
Returns the ID of the element.
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
unsigned getNumberOfPoints() const
MathLib::RowColumnIndices< GlobalIndexType > RowColumnIndices
auto const & NsHigherOrder() const
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
virtual std::vector< double > const & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
virtual void setChemicalSystemConcrete(Eigen::VectorXd const &, double const, double const)=0
virtual void postSpeciationCalculation(std::size_t const ele_id, double const t, double const dt)=0
virtual void computeReactionRelatedSecondaryVariable(std::size_t const ele_id)=0
virtual void setChemicalSystemID(std::size_t const)=0
void initializeChemicalSystem(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t)
void setChemicalSystem(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double const dt)
virtual void assembleReactionEquationConcrete(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, int const transport_process_id)=0
void assembleReactionEquation(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double const dt, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, int const process_id)
virtual void initializeChemicalSystemConcrete(Eigen::VectorXd const &, double const)=0
virtual std::vector< double > const & getIntPtMolarFlux(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache, int const component_id) const =0
ComponentTransportLocalAssemblerInterface()=default
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
void assembleForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
MeshLib::Element const & _element
void assembleWithJacobianHydraulicEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
const int first_concentration_index
void assembleHeatTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &)
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > LocalMatrixType
void postTimestepConcrete(Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, int const) override
std::vector< double > const & getIntPtMolarFlux(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< double > &cache, int const component_id) const override
void postSpeciationCalculation(std::size_t const ele_id, double const t, double const dt) override
NumLib::GenericIntegrationMethod const & _integration_method
void initializeChemicalSystemConcrete(Eigen::VectorXd const &local_x, double const t) override
void computeReactionRelatedSecondaryVariable(std::size_t const ele_id) override
std::vector< double > const & calculateIntPtDarcyVelocity(const double t, Eigen::Ref< const NodalVectorType > const &p_nodal_values, Eigen::Ref< const NodalVectorType > const &C_nodal_values, std::vector< double > &cache) const
void assembleBlockMatrices(GlobalDimVectorType const &b, int const component_id, double const t, double const dt, Eigen::Ref< const NodalVectorType > const &C_nodal_values, Eigen::Ref< const NodalVectorType > const &p_nodal_values, Eigen::Ref< LocalBlockMatrixType > KCC, Eigen::Ref< LocalBlockMatrixType > MCC, Eigen::Ref< LocalBlockMatrixType > MCp, Eigen::Ref< LocalBlockMatrixType > MpC, Eigen::Ref< LocalBlockMatrixType > Kpp, Eigen::Ref< LocalBlockMatrixType > Mpp, Eigen::Ref< LocalSegmentVectorType > Bp)
typename ShapeMatricesType::template MatrixType< pressure_size, pressure_size > LocalBlockMatrixType
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
void assembleWithJacobianForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
static const int temperature_size
std::vector< IntegrationPointData< GlobalDimNodalMatrixType > > _ip_data
ComponentTransportProcessData const & _process_data
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
void computeSecondaryVariableConcrete(double const t, double const, Eigen::VectorXd const &local_x, Eigen::VectorXd const &) override
std::vector< std::reference_wrapper< ProcessVariable > > const _transport_process_variables
void assembleHydraulicEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
typename ShapeMatricesType::GlobalDimNodalMatrixType GlobalDimNodalMatrixType
typename ShapeMatricesType::NodalRowVectorType NodalRowVectorType
std::vector< double > const & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const override
const int temperature_index
static const int concentration_size
void assembleComponentTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &, int const transport_process_id)
NodalVectorType getLocalTemperature(double const t, Eigen::VectorXd const &local_x)
Eigen::Vector3d getFlux(MathLib::Point3d const &pnt_local_coords, double const t, std::vector< double > const &local_x) const override
GlobalDimMatrixType getThermalConductivityDispersivity(MaterialPropertyLib::VariableArray const &vars, const double fluid_density, const double specific_heat_capacity_fluid, const GlobalDimVectorType &velocity, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
void assemble(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
void assembleWithJacobianComponentTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data, int const component_id)
Eigen::Matrix< double, Eigen::Dynamic, 1 > LocalVectorType
static const int pressure_size
double getHeatEnergyCoefficient(MaterialPropertyLib::VariableArray const &vars, const double porosity, const double fluid_density, const double specific_heat_capacity_fluid, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
void assembleReactionEquationConcrete(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, int const transport_process_id) override
LocalAssemblerData(MeshLib::Element const &element, std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const &integration_method, bool is_axially_symmetric, ComponentTransportProcessData const &process_data, std::vector< std::reference_wrapper< ProcessVariable > > const &transport_process_variables)
typename ShapeMatricesType::template VectorType< pressure_size > LocalSegmentVectorType
void assembleKCmCn(int const component_id, double const t, double const dt, Eigen::Ref< LocalBlockMatrixType > KCmCn, double const stoichiometric_coefficient, double const kinetic_prefactor)
static const int pressure_index
void setChemicalSystemConcrete(Eigen::VectorXd const &local_x, double const t, double dt) override
void setChemicalSystemID(std::size_t const) override
typename ShapeMatricesType::NodalVectorType NodalVectorType
@ longitudinal_dispersivity
used to compute the hydrodynamic dispersion tensor.
@ transversal_dispersivity
used to compute the hydrodynamic dispersion tensor.
@ retardation_factor
specify retardation factor used in component transport process.
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
void assembleAdvectionMatrix(IPData const &ip_data_vector, NumLib::ShapeMatrixCache const &shape_matrix_cache, std::vector< FluxVectorType > const &ip_flux_vector, Eigen::MatrixBase< Derived > &laplacian_matrix)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > computeShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, PointContainer const &points)
Eigen::MatrixXd computeHydrodynamicDispersion(NumericalStabilization const &stabilizer, std::size_t const element_id, Eigen::MatrixXd const &pore_diffusion_coefficient, Eigen::VectorXd const &velocity, double const porosity, double const solute_dispersivity_transverse, double const solute_dispersivity_longitudinal)
NumLib::ShapeMatrices< NodalRowVectorType, DimNodalMatrixType, DimMatrixType, GlobalDimNodalMatrixType > ShapeMatrices
MatrixType< GlobalDim, ShapeFunction::NPOINTS > GlobalDimNodalMatrixType
MatrixType< GlobalDim, GlobalDim > GlobalDimMatrixType
VectorType< GlobalDim > GlobalDimVectorType
VectorType< ShapeFunction::NPOINTS > NodalVectorType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
virtual Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > getNodalValuesOnElement(MeshLib::Element const &element, double const t) const
Returns a matrix of values for all nodes of the given element.
NumLib::NumericalStabilization stabilizer
std::vector< Eigen::VectorXd > const projected_specific_body_force_vectors
Projected specific body force vector: R * R^T * b.
MaterialPropertyLib::MaterialSpatialDistributionMap media_map
ChemistryLib::ChemicalSolverInterface *const chemical_solver_interface
static const int hydraulic_process_id
ParameterLib::Parameter< double > const & aperture_size
MeshLib::PropertyVector< double > * mesh_prop_velocity
bool const non_advective_form
NumLib::ShapeMatrixCache shape_matrix_cache
caches for each mesh element type the shape matrix
ParameterLib::Parameter< double > const *const temperature
const int thermal_process_id
bool const chemically_induced_porosity_change
MeshLib::PropertyVector< double > * mesh_prop_porosity
GlobalIndexType chemical_system_id
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
GlobalDimNodalMatrixType const dNdx
IntegrationPointData(GlobalDimNodalMatrixType const &dNdx_, double const &integration_weight_)
double const integration_weight