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;
188 std::vector<GlobalVector*>
const& x,
189 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
190 std::vector<double>& cache)
const = 0;
194 Eigen::VectorXd
const& ,
double const ) = 0;
201 double const t,
double const dt, Eigen::VectorXd
const& local_x,
202 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
203 std::vector<double>& local_b_data,
int const transport_process_id) = 0;
206template <
typename ShapeFunction,
int GlobalDim>
218 ShapeFunction::NPOINTS;
224 typename ShapeMatricesType::template MatrixType<
pressure_size,
227 typename ShapeMatricesType::template VectorType<pressure_size>;
230 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
244 std::size_t
const local_matrix_size,
246 bool is_axially_symmetric,
248 std::vector<std::reference_wrapper<ProcessVariable>>
const&
249 transport_process_variables)
251 : ShapeFunction::NPOINTS),
253 ? ShapeFunction::NPOINTS
254 : 2 * ShapeFunction::NPOINTS),
260 (void)local_matrix_size;
262 unsigned const n_integration_points =
264 _ip_data.reserve(n_integration_points);
271 auto const shape_matrices =
273 GlobalDim>(element, is_axially_symmetric,
277 for (
unsigned ip = 0; ip < n_integration_points; ip++)
280 shape_matrices[ip].dNdx,
282 shape_matrices[ip].integralMeasure *
283 shape_matrices[ip].detJ * aperture_size);
289 .template initialValue<double>(
290 pos, std::numeric_limits<double>::quiet_NaN() );
300 auto& chemical_system_index_map =
303 unsigned const n_integration_points =
305 for (
unsigned ip = 0; ip < n_integration_points; ip++)
308 chemical_system_index_map.empty()
310 : chemical_system_index_map.back() + 1;
311 chemical_system_index_map.push_back(
317 double const t)
override
331 unsigned const n_integration_points =
334 for (
unsigned ip = 0; ip < n_integration_points; ip++)
337 auto const& N = Ns[ip];
338 auto const& chemical_system_id = ip_data.chemical_system_id;
341 std::vector<double> C_int_pt(n_component);
342 for (
unsigned component_id = 0; component_id < n_component;
345 auto const concentration_index =
349 local_x.template segment<concentration_size>(
350 concentration_index);
353 C_int_pt[component_id]);
363 double const t,
double dt)
override
380 unsigned const n_integration_points =
383 for (
unsigned ip = 0; ip < n_integration_points; ip++)
386 auto const& N = Ns[ip];
387 auto& porosity = ip_data.porosity;
388 auto const& porosity_prev = ip_data.porosity_prev;
389 auto const& chemical_system_id = ip_data.chemical_system_id;
392 std::vector<double> C_int_pt(n_component);
393 for (
unsigned component_id = 0; component_id < n_component;
396 auto const concentration_index =
400 local_x.template segment<concentration_size>(
401 concentration_index);
404 C_int_pt[component_id]);
416 .template value<double>(vars, vars_prev, pos, t,
423 C_int_pt, chemical_system_id, medium, vars, pos, t, dt);
428 double const dt)
override
442 ip_data.porosity = ip_data.porosity_prev;
447 ip_data.porosity, t, dt);
450 ip_data.chemical_system_id, medium, ip_data.porosity);
455 std::vector<double>
const& local_x,
456 std::vector<double>
const& ,
457 std::vector<double>& local_M_data,
458 std::vector<double>& local_K_data,
459 std::vector<double>& local_b_data)
override
461 auto const local_matrix_size = local_x.size();
466 assert(local_matrix_size == ShapeFunction::NPOINTS * num_nodal_dof);
468 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
469 local_M_data, local_matrix_size, local_matrix_size);
470 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
471 local_K_data, local_matrix_size, local_matrix_size);
472 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
473 local_b_data, local_matrix_size);
476 auto Kpp = local_K.template block<pressure_size, pressure_size>(
478 auto Mpp = local_M.template block<pressure_size, pressure_size>(
480 auto Bp = local_b.template segment<pressure_size>(
pressure_index);
482 auto local_p = Eigen::Map<const NodalVectorType>(
489 auto const number_of_components = num_nodal_dof - 1;
490 for (
int component_id = 0; component_id < number_of_components;
502 auto concentration_index =
506 local_K.template block<concentration_size, concentration_size>(
507 concentration_index, concentration_index);
509 local_M.template block<concentration_size, concentration_size>(
510 concentration_index, concentration_index);
512 local_M.template block<concentration_size, pressure_size>(
515 local_M.template block<pressure_size, concentration_size>(
518 auto local_C = Eigen::Map<const NodalVectorType>(
522 MCC, MCp, MpC, Kpp, Mpp, Bp);
526 auto const stoichiometric_matrix =
530 assert(stoichiometric_matrix);
532 for (Eigen::SparseMatrix<double>::InnerIterator it(
533 *stoichiometric_matrix, component_id);
537 auto const stoichiometric_coefficient = it.value();
538 auto const coupled_component_id = it.row();
539 auto const kinetic_prefactor =
543 auto const concentration_index =
545 auto const coupled_concentration_index =
550 concentration_index, coupled_concentration_index);
554 stoichiometric_coefficient,
564 Eigen::Ref<const NodalVectorType>
const& C_nodal_values,
565 Eigen::Ref<const NodalVectorType>
const& p_nodal_values,
566 Eigen::Ref<LocalBlockMatrixType> KCC,
567 Eigen::Ref<LocalBlockMatrixType> MCC,
568 Eigen::Ref<LocalBlockMatrixType> MCp,
569 Eigen::Ref<LocalBlockMatrixType> MpC,
570 Eigen::Ref<LocalBlockMatrixType> Kpp,
571 Eigen::Ref<LocalBlockMatrixType> Mpp,
572 Eigen::Ref<LocalSegmentVectorType> Bp)
574 unsigned const n_integration_points =
586 auto const& phase = medium.
phase(
"AqueousLiquid");
597 std::vector<GlobalDimVectorType> ip_flux_vector;
598 double average_velocity_norm = 0.0;
601 ip_flux_vector.reserve(n_integration_points);
608 for (
unsigned ip(0); ip < n_integration_points; ++ip)
613 auto const& dNdx = ip_data.dNdx;
614 auto const& N = Ns[ip];
615 auto const& w = ip_data.integration_weight;
616 auto& porosity = ip_data.porosity;
618 double C_int_pt = 0.0;
619 double p_int_pt = 0.0;
629 .template value<double>(vars, pos, t, dt);
632 auto const& retardation_factor =
634 .template value<double>(vars, pos, t, dt);
636 auto const& solute_dispersivity_transverse = medium.template value<
640 auto const& solute_dispersivity_longitudinal =
641 medium.template value<double>(
643 longitudinal_dispersivity);
650 .template value<double>(vars, pos, t, dt);
652 auto const decay_rate =
654 .template value<double>(vars, pos, t, dt);
656 auto const& pore_diffusion_coefficient =
657 MaterialPropertyLib::formEigenTensor<GlobalDim>(
659 .value(vars, pos, t, dt));
661 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
667 .template value<double>(vars, pos, t, dt);
673 (dNdx * p_nodal_values - density * b))
676 const double drho_dp =
678 .template dValue<double>(
683 const double drho_dC =
685 .template dValue<double>(
692 pore_diffusion_coefficient, velocity, porosity,
693 solute_dispersivity_transverse,
694 solute_dispersivity_longitudinal);
696 const double R_times_phi(retardation_factor * porosity);
698 auto const N_t_N = (N.transpose() * N).eval();
702 MCp.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dp * w);
703 MCC.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dC * w);
704 KCC.noalias() -= dNdx.transpose() * mass_density_flow * N * w;
708 ip_flux_vector.emplace_back(mass_density_flow);
709 average_velocity_norm += velocity.norm();
711 MCC.noalias() += N_t_N * (R_times_phi * density * w);
712 KCC.noalias() += N_t_N * (decay_rate * R_times_phi * density * w);
713 KCC_Laplacian.noalias() +=
714 dNdx.transpose() * hydrodynamic_dispersion * dNdx * density * w;
716 MpC.noalias() += N_t_N * (porosity * drho_dC * w);
719 if (component_id == 0)
721 Mpp.noalias() += N_t_N * (porosity * drho_dp * w);
723 dNdx.transpose() * K_over_mu * dNdx * (density * w);
727 Bp.noalias() += dNdx.transpose() * K_over_mu * b *
728 (density * density * w);
736 typename ShapeFunction::MeshElement>(
741 average_velocity_norm /
742 static_cast<double>(n_integration_points),
746 KCC.noalias() += KCC_Laplacian;
750 Eigen::Ref<LocalBlockMatrixType> KCmCn,
751 double const stoichiometric_coefficient,
752 double const kinetic_prefactor)
754 unsigned const n_integration_points =
764 auto const& phase = medium.
phase(
"AqueousLiquid");
772 for (
unsigned ip(0); ip < n_integration_points; ++ip)
775 auto const& w = ip_data.integration_weight;
776 auto const& N = Ns[ip];
777 auto& porosity = ip_data.porosity;
779 auto const retardation_factor =
781 .template value<double>(vars, pos, t, dt);
784 .template value<double>(vars, pos, t, dt);
788 .template value<double>(vars, pos, t, dt);
790 KCmCn.noalias() -= w * N.transpose() * stoichiometric_coefficient *
791 kinetic_prefactor * retardation_factor *
792 porosity * density * N;
797 Eigen::VectorXd
const& local_x,
798 Eigen::VectorXd
const& local_x_prev,
799 int const process_id,
800 std::vector<double>& local_M_data,
801 std::vector<double>& local_K_data,
802 std::vector<double>& local_b_data)
override
807 local_M_data, local_K_data, local_b_data);
812 local_M_data, local_K_data,
819 local_M_data, local_K_data,
820 local_b_data, process_id);
826 Eigen::VectorXd
const& local_x,
827 Eigen::VectorXd
const& local_x_prev,
828 std::vector<double>& local_M_data,
829 std::vector<double>& local_K_data,
830 std::vector<double>& local_b_data)
834 auto const local_C = local_x.template segment<concentration_size>(
836 auto const local_C_prev =
841 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
843 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
845 auto local_b = MathLib::createZeroedVector<LocalSegmentVectorType>(
848 unsigned const n_integration_points =
860 auto const& phase = medium.
phase(
"AqueousLiquid");
869 for (
unsigned ip(0); ip < n_integration_points; ++ip)
874 auto const& dNdx = ip_data.dNdx;
875 auto const& w = ip_data.integration_weight;
876 auto const& N = Ns[ip];
877 auto& porosity = ip_data.porosity;
878 auto const& porosity_prev = ip_data.porosity_prev;
880 double const C_int_pt = N.dot(local_C);
881 double const p_int_pt = N.dot(local_p);
882 double const T_int_pt = N.dot(local_T);
896 .template value<double>(vars, vars_prev, pos, t,
907 .template value<double>(vars, pos, t, dt);
909 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
915 .template value<double>(vars, pos, t, dt);
919 const double drho_dp =
921 .template dValue<double>(
925 const double drho_dC =
927 .template dValue<double>(
932 local_M.noalias() += w * N.transpose() * porosity * drho_dp * N;
934 w * dNdx.transpose() * density * K_over_mu * dNdx;
939 w * density * density * dNdx.transpose() * K_over_mu * b;
944 double const C_dot = (C_int_pt - N.dot(local_C_prev)) / dt;
947 w * N.transpose() * porosity * drho_dC * C_dot;
953 Eigen::VectorXd
const& local_x,
954 Eigen::VectorXd
const& ,
955 std::vector<double>& local_M_data,
956 std::vector<double>& local_K_data,
957 std::vector<double>& )
959 assert(local_x.size() ==
967 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
969 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
978 auto const& liquid_phase = medium.
phase(
"AqueousLiquid");
986 unsigned const n_integration_points =
989 std::vector<GlobalDimVectorType> ip_flux_vector;
990 double average_velocity_norm = 0.0;
991 ip_flux_vector.reserve(n_integration_points);
997 for (
unsigned ip(0); ip < n_integration_points; ip++)
1001 auto const& ip_data = this->
_ip_data[ip];
1002 auto const& dNdx = ip_data.dNdx;
1003 auto const& w = ip_data.integration_weight;
1004 auto const& N = Ns[ip];
1006 double p_at_xi = 0.;
1008 double T_at_xi = 0.;
1016 auto const porosity =
1018 .template value<double>(vars, pos, t, dt);
1022 auto const fluid_density =
1025 .template value<double>(vars, pos, t, dt);
1027 auto const specific_heat_capacity_fluid =
1030 .template value<double>(vars, pos, t, dt);
1033 local_M.noalias() += w *
1035 vars, porosity, fluid_density,
1036 specific_heat_capacity_fluid, pos, t, dt) *
1040 auto const viscosity =
1043 .template value<double>(vars, pos, t, dt);
1045 auto const intrinsic_permeability =
1046 MaterialPropertyLib::formEigenTensor<GlobalDim>(
1050 .value(vars, pos, t, dt));
1053 intrinsic_permeability / viscosity;
1055 process_data.has_gravity
1057 (dNdx * local_p - fluid_density * b))
1062 vars, fluid_density, specific_heat_capacity_fluid, velocity,
1065 local_K.noalias() +=
1066 w * dNdx.transpose() * thermal_conductivity_dispersivity * dNdx;
1068 ip_flux_vector.emplace_back(velocity * fluid_density *
1069 specific_heat_capacity_fluid);
1070 average_velocity_norm += velocity.norm();
1073 NumLib::assembleAdvectionMatrix<typename ShapeFunction::MeshElement>(
1074 process_data.stabilizer, this->_ip_data,
1076 average_velocity_norm /
static_cast<double>(n_integration_points),
1081 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1082 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_M_data,
1083 std::vector<double>& local_K_data,
1084 std::vector<double>& ,
int const transport_process_id)
1086 assert(
static_cast<int>(local_x.size()) ==
1092 auto const local_p =
1097 auto const local_C = local_x.template segment<concentration_size>(
1101 auto const local_p_prev =
1104 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1106 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1112 unsigned const n_integration_points =
1115 std::vector<GlobalDimVectorType> ip_flux_vector;
1116 double average_velocity_norm = 0.0;
1119 ip_flux_vector.reserve(n_integration_points);
1132 auto const& medium =
1134 auto const& phase = medium.
phase(
"AqueousLiquid");
1135 auto const component_id =
1137 auto const& component = phase.
component(
1144 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1149 auto const& dNdx = ip_data.dNdx;
1150 auto const& w = ip_data.integration_weight;
1151 auto const& N = Ns[ip];
1152 auto& porosity = ip_data.porosity;
1153 auto const& porosity_prev = ip_data.porosity_prev;
1155 double const C_int_pt = N.dot(local_C);
1156 double const p_int_pt = N.dot(local_p);
1157 double const T_int_pt = N.dot(local_T);
1170 vars_prev.
porosity = porosity_prev;
1176 .template value<double>(vars, vars_prev, pos, t,
1182 auto const& retardation_factor =
1184 .template value<double>(vars, pos, t, dt);
1186 auto const& solute_dispersivity_transverse = medium.template value<
1189 auto const& solute_dispersivity_longitudinal =
1190 medium.template value<double>(
1192 longitudinal_dispersivity);
1195 auto const density =
1197 .template value<double>(vars, pos, t, dt);
1198 auto const decay_rate =
1200 .template value<double>(vars, pos, t, dt);
1202 auto const& pore_diffusion_coefficient =
1203 MaterialPropertyLib::formEigenTensor<GlobalDim>(
1205 .value(vars, pos, t, dt));
1207 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1212 .template value<double>(vars, pos, t, dt);
1218 (dNdx * local_p - density * b))
1224 pore_diffusion_coefficient, velocity, porosity,
1225 solute_dispersivity_transverse,
1226 solute_dispersivity_longitudinal);
1228 double const R_times_phi = retardation_factor * porosity;
1229 auto const N_t_N = (N.transpose() * N).eval();
1233 const double drho_dC =
1235 .template dValue<double>(
1238 local_M.noalias() +=
1239 N_t_N * (R_times_phi * C_int_pt * drho_dC * w);
1242 local_M.noalias() += N_t_N * (R_times_phi * density * w);
1247 double const p_dot = (p_int_pt - N.dot(local_p_prev)) / dt;
1249 const double drho_dp =
1251 .template dValue<double>(vars,
1253 liquid_phase_pressure,
1256 local_K.noalias() +=
1257 N_t_N * ((R_times_phi * drho_dp * p_dot) * w) -
1258 dNdx.transpose() * velocity * N * (density * w);
1262 ip_flux_vector.emplace_back(velocity * density);
1263 average_velocity_norm += velocity.norm();
1265 local_K.noalias() +=
1266 N_t_N * (decay_rate * R_times_phi * density * w);
1268 KCC_Laplacian.noalias() += dNdx.transpose() *
1269 hydrodynamic_dispersion * dNdx *
1276 typename ShapeFunction::MeshElement>(
1279 average_velocity_norm /
1280 static_cast<double>(n_integration_points),
1283 local_K.noalias() += KCC_Laplacian;
1287 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1288 Eigen::VectorXd
const& local_x_prev,
int const process_id,
1289 std::vector<double>& ,
1290 std::vector<double>& ,
1291 std::vector<double>& local_b_data,
1292 std::vector<double>& local_Jac_data)
override
1297 local_b_data, local_Jac_data);
1301 int const component_id = process_id - 1;
1303 t, dt, local_x, local_x_prev, local_b_data, local_Jac_data,
1309 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1310 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
1311 std::vector<double>& local_Jac_data)
1313 auto const p = local_x.template segment<pressure_size>(
pressure_index);
1314 auto const c = local_x.template segment<concentration_size>(
1321 auto local_Jac = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1323 auto local_rhs = MathLib::createZeroedVector<LocalSegmentVectorType>(
1326 unsigned const n_integration_points =
1335 auto const& medium =
1337 auto const& phase = medium.
phase(
"AqueousLiquid");
1346 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1351 auto const& dNdx = ip_data.dNdx;
1352 auto const& w = ip_data.integration_weight;
1353 auto const& N = Ns[ip];
1354 auto& phi = ip_data.porosity;
1355 auto const& phi_prev = ip_data.porosity_prev;
1357 double const p_ip = N.dot(p);
1358 double const c_ip = N.dot(c);
1360 double const cdot_ip = (c_ip - N.dot(c_prev)) / dt;
1372 .template value<double>(vars, vars_prev, pos, t,
1379 .template value<double>(vars, pos, t, dt);
1381 auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1386 .template value<double>(vars, pos, t, dt);
1388 auto const drho_dp =
1390 .template dValue<double>(
1394 auto const drho_dc =
1396 .template dValue<double>(
1401 local_Jac.noalias() += w * N.transpose() * phi * drho_dp / dt * N +
1402 w * dNdx.transpose() * rho * k / mu * dNdx;
1404 local_rhs.noalias() -=
1405 w * N.transpose() * phi *
1406 (drho_dp * N * p_prev + drho_dc * cdot_ip) +
1407 w * rho * dNdx.transpose() * k / mu * dNdx * p;
1411 local_rhs.noalias() +=
1412 w * rho * dNdx.transpose() * k / mu * rho * b;
1418 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1419 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
1420 std::vector<double>& local_Jac_data,
int const component_id)
1422 auto const concentration_index =
1425 auto const p = local_x.template segment<pressure_size>(
pressure_index);
1427 local_x.template segment<concentration_size>(concentration_index);
1437 auto local_Jac = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1439 auto local_rhs = MathLib::createZeroedVector<LocalSegmentVectorType>(
1445 unsigned const n_integration_points =
1448 std::vector<GlobalDimVectorType> ip_flux_vector;
1449 double average_velocity_norm = 0.0;
1450 ip_flux_vector.reserve(n_integration_points);
1462 auto const& medium =
1464 auto const& phase = medium.
phase(
"AqueousLiquid");
1465 auto const& component = phase.
component(
1472 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1477 auto const& dNdx = ip_data.dNdx;
1478 auto const& w = ip_data.integration_weight;
1479 auto const& N = Ns[ip];
1480 auto& phi = ip_data.porosity;
1481 auto const& phi_prev = ip_data.porosity_prev;
1483 double const p_ip = N.dot(p);
1484 double const c_ip = N.dot(c);
1501 .template value<double>(vars, vars_prev, pos, t,
1509 .template value<double>(vars, pos, t, dt);
1511 auto const alpha_T = medium.template value<double>(
1513 auto const alpha_L = medium.template value<double>(
1517 .template value<double>(vars, pos, t, dt);
1521 .template value<double>(vars, pos, t, dt);
1523 auto const Dp = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1525 .value(vars, pos, t, dt));
1527 auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1531 .template value<double>(vars, pos, t, dt);
1543 local_Jac.noalias() +=
1544 w * rho * N.transpose() * phi * R * (alpha + 1 / dt) * N;
1546 KCC_Laplacian.noalias() += w * rho * dNdx.transpose() * D * dNdx;
1548 auto const cdot = (c - c_prev) / dt;
1549 local_rhs.noalias() -=
1550 w * rho * N.transpose() * phi * R * N * (cdot + alpha * c);
1552 ip_flux_vector.emplace_back(q * rho);
1553 average_velocity_norm += q.norm();
1556 NumLib::assembleAdvectionMatrix<typename ShapeFunction::MeshElement>(
1559 average_velocity_norm /
static_cast<double>(n_integration_points),
1562 local_rhs.noalias() -= KCC_Laplacian * c;
1564 local_Jac.noalias() += KCC_Laplacian;
1568 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1569 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
1570 std::vector<double>& local_b_data,
1571 int const transport_process_id)
override
1573 auto const local_C = local_x.template segment<concentration_size>(
1577 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1579 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1581 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
1584 unsigned const n_integration_points =
1593 auto const& medium =
1595 auto const component_id = transport_process_id - 1;
1601 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1606 auto const w = ip_data.integration_weight;
1607 auto const& N = Ns[ip];
1608 auto& porosity = ip_data.porosity;
1609 auto const& porosity_prev = ip_data.porosity_prev;
1610 auto const chemical_system_id = ip_data.chemical_system_id;
1612 double C_int_pt = 0.0;
1617 auto const porosity_dot = (porosity - porosity_prev) / dt;
1621 vars_prev.
porosity = porosity_prev;
1627 .template value<double>(vars, vars_prev, pos, t,
1631 local_M.noalias() += w * N.transpose() * porosity * N;
1633 local_K.noalias() += w * N.transpose() * porosity_dot * N;
1635 if (chemical_system_id == -1)
1640 auto const C_post_int_pt =
1642 component_id, chemical_system_id);
1644 local_b.noalias() +=
1645 w * N.transpose() * porosity * (C_post_int_pt - C_int_pt) / dt;
1651 std::vector<GlobalVector*>
const& x,
1652 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
1653 std::vector<double>& cache)
const override
1655 assert(x.size() == dof_table.size());
1657 auto const n_processes = x.size();
1658 std::vector<std::vector<double>> local_x;
1659 local_x.reserve(n_processes);
1661 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1663 auto const indices =
1665 assert(!indices.empty());
1666 local_x.push_back(x[process_id]->get(indices));
1670 if (n_processes == 1)
1672 auto const local_p = Eigen::Map<const NodalVectorType>(
1674 auto const local_C = Eigen::Map<const NodalVectorType>(
1681 constexpr int pressure_process_id = 0;
1682 constexpr int concentration_process_id = 1;
1683 auto const local_p = Eigen::Map<const NodalVectorType>(
1685 auto const local_C = Eigen::Map<const NodalVectorType>(
1693 Eigen::Ref<const NodalVectorType>
const& p_nodal_values,
1694 Eigen::Ref<const NodalVectorType>
const& C_nodal_values,
1695 std::vector<double>& cache)
const
1697 auto const n_integration_points =
1702 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1703 cache, GlobalDim, n_integration_points);
1714 auto const& medium =
1716 auto const& phase = medium.
phase(
"AqueousLiquid");
1722 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
1724 auto const& ip_data =
_ip_data[ip];
1725 auto const& dNdx = ip_data.dNdx;
1726 auto const& N = Ns[ip];
1727 auto const& porosity = ip_data.porosity;
1729 pos.setIntegrationPoint(ip);
1731 double C_int_pt = 0.0;
1732 double p_int_pt = 0.0;
1743 double const dt = std::numeric_limits<double>::quiet_NaN();
1744 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1748 .template value<double>(vars, pos, t, dt);
1751 cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
1756 .template value<double>(vars, pos, t, dt);
1758 cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
1766 const unsigned integration_point)
const override
1769 typename ShapeFunction::MeshElement>()[integration_point];
1772 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
1777 std::vector<double>
const& local_x)
const override
1779 auto const local_p = Eigen::Map<const NodalVectorType>(
1781 auto const local_C = Eigen::Map<const NodalVectorType>(
1787 auto const shape_matrices =
1791 std::array{pnt_local_coords})[0];
1801 auto const& medium =
1803 auto const& phase = medium.
phase(
"AqueousLiquid");
1816 double const dt = std::numeric_limits<double>::quiet_NaN();
1817 auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1822 .template value<double>(vars, pos, t, dt);
1827 .template value<double>(vars, pos, t, dt);
1830 q += K_over_mu * rho_w * b;
1832 Eigen::Vector3d flux(0.0, 0.0, 0.0);
1833 flux.head<GlobalDim>() = rho_w * q;
1840 Eigen::VectorXd
const& local_x,
1841 Eigen::VectorXd
const& )
override
1843 auto const local_p =
1845 auto const local_C = local_x.template segment<concentration_size>(
1848 std::vector<double> ele_velocity;
1851 auto const n_integration_points =
1853 auto const ele_velocity_mat =
1857 Eigen::Map<LocalVectorType>(
1860 ele_velocity_mat.rowwise().sum() / n_integration_points;
1864 std::size_t
const ele_id)
override
1866 auto const n_integration_points =
1875 ip_data.porosity = ip_data.porosity_prev;
1879 medium, ip_data.porosity);
1884 [](
double const s,
auto const& ip)
1885 {
return s + ip.porosity; }) /
1886 n_integration_points;
1889 std::vector<GlobalIndexType> chemical_system_indices;
1890 chemical_system_indices.reserve(n_integration_points);
1892 std::back_inserter(chemical_system_indices),
1893 [](
auto const& ip_data)
1894 { return ip_data.chemical_system_id; });
1897 ele_id, chemical_system_indices);
1902 std::vector<GlobalVector*>
const& x,
1903 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
1904 std::vector<double>& cache)
const override
1906 std::vector<double> local_x_vec;
1908 auto const n_processes = x.size();
1909 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1911 auto const indices =
1913 assert(!indices.empty());
1914 auto const local_solution = x[process_id]->get(indices);
1915 local_x_vec.insert(std::end(local_x_vec),
1916 std::begin(local_solution),
1917 std::end(local_solution));
1921 auto const p = local_x.template segment<pressure_size>(
pressure_index);
1922 auto const c = local_x.template segment<concentration_size>(
1925 auto const n_integration_points =
1930 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1931 cache, GlobalDim, n_integration_points);
1942 auto const& medium =
1944 auto const& phase = medium.
phase(
"AqueousLiquid");
1946 int const component_id = 0;
1947 auto const& component = phase.
component(
1954 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
1956 auto const& ip_data =
_ip_data[ip];
1957 auto const& dNdx = ip_data.dNdx;
1958 auto const& N = Ns[ip];
1959 auto const& phi = ip_data.porosity;
1961 pos.setIntegrationPoint(ip);
1963 double const p_ip = N.dot(p);
1964 double const c_ip = N.dot(c);
1970 double const dt = std::numeric_limits<double>::quiet_NaN();
1972 auto const& k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1976 .template value<double>(vars, pos, t, dt);
1978 .template value<double>(vars, pos, t, dt);
1986 auto const alpha_T = medium.template value<double>(
1988 auto const alpha_L = medium.template value<double>(
1990 auto const Dp = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1992 .value(vars, pos, t, dt));
1999 cache_mat.col(ip).noalias() = q * c_ip - phi * D * dNdx * c;
2006 Eigen::VectorXd
const& ,
2007 double const ,
double const ,
2008 int const )
override
2010 unsigned const n_integration_points =
2013 for (
unsigned ip = 0; ip < n_integration_points; ip++)
2024 std::vector<std::reference_wrapper<ProcessVariable>>
const
2027 std::vector<IntegrationPointData<GlobalDimNodalMatrixType>>
_ip_data;
2031 const double fluid_density,
const double specific_heat_capacity_fluid,
2035 auto const& medium =
2037 auto const& solid_phase = medium.
phase(
"Solid");
2039 auto const specific_heat_capacity_solid =
2043 .template value<double>(vars, pos, t, dt);
2045 auto const solid_density =
2047 .template value<double>(vars, pos, t, dt);
2049 return solid_density * specific_heat_capacity_solid * (1 - porosity) +
2050 fluid_density * specific_heat_capacity_fluid * porosity;
2055 const double fluid_density,
const double specific_heat_capacity_fluid,
2060 auto const& medium =
2063 auto thermal_conductivity =
2064 MaterialPropertyLib::formEigenTensor<GlobalDim>(
2068 .value(vars, pos, t, dt));
2070 auto const thermal_dispersivity_transversal =
2073 thermal_transversal_dispersivity)
2074 .
template value<double>();
2076 auto const thermal_dispersivity_longitudinal =
2079 thermal_longitudinal_dispersivity)
2080 .
template value<double>();
2085 return thermal_conductivity +
2086 fluid_density * specific_heat_capacity_fluid *
2089 GlobalDimMatrixType::Zero(GlobalDim, GlobalDim),
2090 velocity, 0 , thermal_dispersivity_transversal,
2091 thermal_dispersivity_longitudinal);
2095 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 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) 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
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
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 & getIntPtMolarFlux(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< double > &cache) const 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 > &, std::vector< double > &, 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