13 #include <Eigen/Dense>
37 namespace ComponentTransport
39 template <
typename NodalRowVectorType,
typename GlobalDimNodalMatrixType>
43 GlobalDimNodalMatrixType
const& dNdx_,
44 double const& integration_weight_)
51 NodalRowVectorType
const N;
52 GlobalDimNodalMatrixType
const dNdx;
57 double porosity = std::numeric_limits<double>::quiet_NaN();
79 std::size_t
const mesh_item_id,
80 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
81 std::vector<GlobalVector*>
const& x,
double const t)
83 std::vector<double> local_x_vec;
85 auto const n_processes = x.size();
86 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
90 assert(!indices.empty());
91 auto const local_solution = x[process_id]->get(indices);
92 local_x_vec.insert(std::end(local_x_vec),
93 std::begin(local_solution),
94 std::end(local_solution));
102 std::size_t
const mesh_item_id,
103 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
104 std::vector<GlobalVector*>
const& x,
double const t,
double const dt)
106 std::vector<double> local_x_vec;
108 auto const n_processes = x.size();
109 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
113 assert(!indices.empty());
114 auto const local_solution = x[process_id]->get(indices);
115 local_x_vec.insert(std::end(local_x_vec),
116 std::begin(local_solution),
117 std::end(local_solution));
125 std::size_t
const mesh_item_id,
126 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
127 std::vector<GlobalVector*>
const& x,
double const t,
double const dt,
130 std::vector<double> local_x_vec;
132 auto const n_processes = x.size();
133 for (std::size_t pcs_id = 0; pcs_id < n_processes; ++pcs_id)
137 assert(!indices.empty());
138 auto const local_solution = x[pcs_id]->get(indices);
139 local_x_vec.insert(std::end(local_x_vec),
140 std::begin(local_solution),
141 std::end(local_solution));
147 auto const num_r_c = indices.size();
149 std::vector<double> local_M_data;
150 local_M_data.reserve(num_r_c * num_r_c);
151 std::vector<double> local_K_data;
152 local_K_data.reserve(num_r_c * num_r_c);
153 std::vector<double> local_b_data;
154 local_b_data.reserve(num_r_c);
157 local_K_data, local_b_data,
160 auto const r_c_indices =
162 if (!local_M_data.empty())
166 M.
add(r_c_indices, local_M);
168 if (!local_K_data.empty())
172 K.
add(r_c_indices, local_K);
174 if (!local_b_data.empty())
176 b.
add(indices, local_b_data);
181 double const t,
double const dt) = 0;
185 std::vector<GlobalVector*>
const& x,
186 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
187 std::vector<double>& cache)
const = 0;
191 std::vector<GlobalVector*>
const& int_pt_x,
192 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
193 std::vector<double>& cache)
const = 0;
197 Eigen::VectorXd
const& ,
double const ) = 0;
204 double const t,
double const dt, Eigen::VectorXd
const& local_x,
205 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
206 std::vector<double>& local_b_data,
int const transport_process_id) = 0;
212 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
222 ShapeFunction::NPOINTS;
228 typename ShapeMatricesType::template MatrixType<
pressure_size,
231 typename ShapeMatricesType::template VectorType<pressure_size>;
234 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
248 std::size_t
const local_matrix_size,
249 bool is_axially_symmetric,
250 unsigned const integration_order,
252 std::vector<std::reference_wrapper<ProcessVariable>>
const&
253 transport_process_variables)
259 (void)local_matrix_size;
261 unsigned const n_integration_points =
263 _ip_data.reserve(n_integration_points);
268 auto const shape_matrices =
270 GlobalDim>(element, is_axially_symmetric,
274 for (
unsigned ip = 0; ip < n_integration_points; ip++)
277 shape_matrices[ip].N, shape_matrices[ip].dNdx,
279 shape_matrices[ip].integralMeasure *
280 shape_matrices[ip].detJ);
286 .template initialValue<double>(
287 pos, std::numeric_limits<double>::quiet_NaN() );
297 auto& chemical_system_index_map =
300 unsigned const n_integration_points =
302 for (
unsigned ip = 0; ip < n_integration_points; ip++)
305 chemical_system_index_map.empty()
307 : chemical_system_index_map.back() + 1;
308 chemical_system_index_map.push_back(
314 double const t)
override
324 unsigned const n_integration_points =
326 for (
unsigned ip = 0; ip < n_integration_points; ip++)
329 auto const& N = ip_data.N;
330 auto const& chemical_system_id = ip_data.chemical_system_id;
333 std::vector<double> C_int_pt(n_component);
334 for (
unsigned component_id = 0; component_id < n_component;
337 auto const concentration_index =
341 local_x.template segment<concentration_size>(
342 concentration_index);
345 C_int_pt[component_id]);
355 double const t,
double dt)
override
368 unsigned const n_integration_points =
370 for (
unsigned ip = 0; ip < n_integration_points; ip++)
373 auto const& N = ip_data.N;
374 auto& porosity = ip_data.porosity;
375 auto const& porosity_prev = ip_data.porosity_prev;
376 auto const& chemical_system_id = ip_data.chemical_system_id;
379 std::vector<double> C_int_pt(n_component);
380 for (
unsigned component_id = 0; component_id < n_component;
383 auto const concentration_index =
387 local_x.template segment<concentration_size>(
388 concentration_index);
391 C_int_pt[component_id]);
395 vars_prev[
static_cast<int>(
404 .template value<double>(vars, vars_prev, pos, t,
407 vars[
static_cast<int>(
412 C_int_pt, chemical_system_id, medium, vars, pos, t, dt);
417 double const dt)
override
431 ip_data.porosity = ip_data.porosity_prev;
436 ip_data.porosity, t, dt);
439 ip_data.chemical_system_id, medium, ip_data.porosity);
444 std::vector<double>
const& local_x,
445 std::vector<double>
const& ,
446 std::vector<double>& local_M_data,
447 std::vector<double>& local_K_data,
448 std::vector<double>& local_b_data)
override
450 auto const local_matrix_size = local_x.size();
455 assert(local_matrix_size == ShapeFunction::NPOINTS * num_nodal_dof);
457 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
458 local_M_data, local_matrix_size, local_matrix_size);
459 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
460 local_K_data, local_matrix_size, local_matrix_size);
461 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
462 local_b_data, local_matrix_size);
465 auto Kpp = local_K.template block<pressure_size, pressure_size>(
467 auto Mpp = local_M.template block<pressure_size, pressure_size>(
469 auto Bp = local_b.template segment<pressure_size>(
pressure_index);
471 auto local_p = Eigen::Map<const NodalVectorType>(
474 auto const number_of_components = num_nodal_dof - 1;
475 for (
int component_id = 0; component_id < number_of_components;
487 auto concentration_index =
491 local_K.template block<concentration_size, concentration_size>(
492 concentration_index, concentration_index);
494 local_M.template block<concentration_size, concentration_size>(
495 concentration_index, concentration_index);
497 local_M.template block<concentration_size, pressure_size>(
500 local_M.template block<pressure_size, concentration_size>(
503 auto local_C = Eigen::Map<const NodalVectorType>(
507 MCC, MCp, MpC, Kpp, Mpp, Bp);
512 int const component_id,
double const t,
double const dt,
513 Eigen::Ref<const NodalVectorType>
const& C_nodal_values,
514 Eigen::Ref<const NodalVectorType>
const& p_nodal_values,
515 Eigen::Ref<LocalBlockMatrixType> KCC,
516 Eigen::Ref<LocalBlockMatrixType> MCC,
517 Eigen::Ref<LocalBlockMatrixType> MCp,
518 Eigen::Ref<LocalBlockMatrixType> MpC,
519 Eigen::Ref<LocalBlockMatrixType> Kpp,
520 Eigen::Ref<LocalBlockMatrixType> Mpp,
521 Eigen::Ref<LocalSegmentVectorType> Bp)
523 unsigned const n_integration_points =
534 GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
540 auto const& phase = medium.phase(
"AqueousLiquid");
545 auto const& component = phase.component(
548 for (
unsigned ip(0); ip < n_integration_points; ++ip)
553 auto const& N = ip_data.N;
554 auto const& dNdx = ip_data.dNdx;
555 auto const& w = ip_data.integration_weight;
556 auto& porosity = ip_data.porosity;
558 double C_int_pt = 0.0;
559 double p_int_pt = 0.0;
564 vars[
static_cast<int>(
566 vars[
static_cast<int>(
571 .template value<double>(vars, pos, t, dt);
577 .template value<double>(vars, pos, t, dt);
579 auto const& solute_dispersivity_transverse = medium.template value<
583 auto const& solute_dispersivity_longitudinal =
584 medium.template value<double>(
593 .template value<double>(vars, pos, t, dt);
597 .template value<double>(vars, pos, t, dt);
599 auto const& pore_diffusion_coefficient =
600 MaterialPropertyLib::formEigenTensor<GlobalDim>(
602 .value(vars, pos, t, dt));
604 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
610 .template value<double>(vars, pos, t, dt);
616 (dNdx * p_nodal_values - density * b))
619 const double drho_dp =
621 .template dValue<double>(
625 const double drho_dC =
627 .template dValue<double>(
631 double const velocity_magnitude = velocity.norm();
633 velocity_magnitude != 0.0
635 pore_diffusion_coefficient +
636 solute_dispersivity_transverse *
637 velocity_magnitude * I +
638 (solute_dispersivity_longitudinal -
639 solute_dispersivity_transverse) /
640 velocity_magnitude * velocity *
641 velocity.transpose())
643 pore_diffusion_coefficient +
644 solute_dispersivity_transverse *
645 velocity_magnitude * I);
648 auto const N_t_N = (N.transpose() * N).eval();
651 MCp.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dp * w);
652 MCC.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dC * w);
653 KCC.noalias() -= dNdx.transpose() * mass_density_flow * N * w;
658 N.transpose() * mass_density_flow.transpose() * dNdx * w;
660 MCC.noalias() += N_t_N * (R_times_phi * density * w);
661 KCC.noalias() += dNdx.transpose() * hydrodynamic_dispersion * dNdx *
663 N_t_N * (
decay_rate * R_times_phi * density * w);
665 MpC.noalias() += N_t_N * (porosity * drho_dC * w);
668 if (component_id == 0)
670 Mpp.noalias() += N_t_N * (porosity * drho_dp * w);
672 dNdx.transpose() * K_over_mu * dNdx * (density * w);
676 Bp.noalias() += dNdx.transpose() * K_over_mu * b *
677 (density * density * w);
684 Eigen::VectorXd
const& local_x,
685 Eigen::VectorXd
const& local_xdot,
686 int const process_id,
687 std::vector<double>& local_M_data,
688 std::vector<double>& local_K_data,
689 std::vector<double>& local_b_data)
override
694 local_K_data, local_b_data);
700 local_M_data, local_K_data,
701 local_b_data, process_id);
707 Eigen::VectorXd
const& local_x,
708 Eigen::VectorXd
const& local_xdot,
709 std::vector<double>& local_M_data,
710 std::vector<double>& local_K_data,
711 std::vector<double>& local_b_data)
715 auto const local_C = local_x.template segment<concentration_size>(
717 auto const local_Cdot =
720 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
722 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
724 auto local_b = MathLib::createZeroedVector<LocalSegmentVectorType>(
727 unsigned const n_integration_points =
737 auto const& phase = medium.phase(
"AqueousLiquid");
742 for (
unsigned ip(0); ip < n_integration_points; ++ip)
747 auto const& N = ip_data.N;
748 auto const& dNdx = ip_data.dNdx;
749 auto const& w = ip_data.integration_weight;
750 auto& porosity = ip_data.porosity;
751 auto const& porosity_prev = ip_data.porosity_prev;
753 double C_int_pt = 0.0;
754 double p_int_pt = 0.0;
759 vars[
static_cast<int>(
761 vars[
static_cast<int>(
766 vars_prev[
static_cast<int>(
773 .template value<double>(vars, vars_prev, pos, t,
776 vars[
static_cast<int>(
785 .template value<double>(vars, pos, t, dt);
787 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
793 .template value<double>(vars, pos, t, dt);
797 const double drho_dp =
799 .template dValue<double>(
802 const double drho_dC =
804 .template dValue<double>(
809 local_M.noalias() += w * N.transpose() * porosity * drho_dp * N;
811 w * dNdx.transpose() * density * K_over_mu * dNdx;
816 w * density * density * dNdx.transpose() * K_over_mu * b;
821 double dot_C_int_pt = 0.0;
825 w * N.transpose() * porosity * drho_dC * dot_C_int_pt;
831 double const t,
double const dt, Eigen::VectorXd
const& local_x,
832 Eigen::VectorXd
const& local_xdot, std::vector<double>& local_M_data,
833 std::vector<double>& local_K_data,
834 std::vector<double>& ,
int const transport_process_id)
838 auto const local_C = local_x.template segment<concentration_size>(
841 auto const local_pdot =
851 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
853 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
856 unsigned const n_integration_points =
868 GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
872 auto const& phase = medium.phase(
"AqueousLiquid");
875 auto const component_id = transport_process_id - 1;
876 auto const& component = phase.component(
879 for (
unsigned ip(0); ip < n_integration_points; ++ip)
884 auto const& N = ip_data.N;
885 auto const& dNdx = ip_data.dNdx;
886 auto const& w = ip_data.integration_weight;
887 auto& porosity = ip_data.porosity;
888 auto const& porosity_prev = ip_data.porosity_prev;
890 double C_int_pt = 0.0;
891 double p_int_pt = 0.0;
896 vars[
static_cast<int>(
898 vars[
static_cast<int>(
903 vars[
static_cast<int>(
910 vars_prev[
static_cast<int>(
917 .template value<double>(vars, vars_prev, pos, t,
920 vars[
static_cast<int>(
926 .template value<double>(vars, pos, t, dt);
928 auto const& solute_dispersivity_transverse = medium.template value<
931 auto const& solute_dispersivity_longitudinal =
932 medium.template value<double>(
939 .template value<double>(vars, pos, t, dt);
942 .template value<double>(vars, pos, t, dt);
944 auto const& pore_diffusion_coefficient =
945 MaterialPropertyLib::formEigenTensor<GlobalDim>(
947 .value(vars, pos, t, dt));
949 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
954 .template value<double>(vars, pos, t, dt);
960 (dNdx * local_p - density * b))
963 double const velocity_magnitude = velocity.norm();
965 velocity_magnitude != 0.0
967 pore_diffusion_coefficient +
968 solute_dispersivity_transverse *
969 velocity_magnitude * I +
970 (solute_dispersivity_longitudinal -
971 solute_dispersivity_transverse) /
972 velocity_magnitude * velocity *
973 velocity.transpose())
975 pore_diffusion_coefficient +
976 solute_dispersivity_transverse *
977 velocity_magnitude * I);
980 auto const N_t_N = (N.transpose() * N).eval();
984 const double drho_dC =
986 .template dValue<double>(
990 N_t_N * (R_times_phi * C_int_pt * drho_dC * w);
993 local_M.noalias() += N_t_N * (R_times_phi * density * w);
999 double dot_p_int_pt = 0.0;
1002 const double drho_dp =
1004 .template dValue<double>(
1008 local_K.noalias() +=
1009 N_t_N * ((R_times_phi * drho_dp * dot_p_int_pt) * w) -
1010 dNdx.transpose() * velocity * N * (density * w);
1014 local_K.noalias() +=
1015 N.transpose() * velocity.transpose() * dNdx * (density * w);
1017 local_K.noalias() +=
1018 dNdx.transpose() * hydrodynamic_dispersion * dNdx *
1020 N_t_N * (
decay_rate * R_times_phi * density * w);
1025 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1026 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
1027 std::vector<double>& local_b_data,
1028 int const transport_process_id)
override
1030 auto const local_C = local_x.template segment<concentration_size>(
1034 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1036 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1038 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
1041 unsigned const n_integration_points =
1050 auto const& medium =
1052 auto const component_id = transport_process_id - 1;
1053 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1058 auto const& N = ip_data.N;
1059 auto const w = ip_data.integration_weight;
1060 auto& porosity = ip_data.porosity;
1061 auto const& porosity_prev = ip_data.porosity_prev;
1062 auto const chemical_system_id = ip_data.chemical_system_id;
1064 double C_int_pt = 0.0;
1067 vars[
static_cast<int>(
1070 auto const porosity_dot = (porosity - porosity_prev) / dt;
1074 vars_prev[
static_cast<int>(
1081 .template value<double>(vars, vars_prev, pos, t,
1085 local_M.noalias() += w * N.transpose() * porosity * N;
1087 local_K.noalias() += w * N.transpose() * porosity_dot * N;
1089 auto const C_post_int_pt =
1091 component_id, chemical_system_id);
1093 local_b.noalias() +=
1094 w * N.transpose() * porosity * (C_post_int_pt - C_int_pt) / dt;
1100 std::vector<GlobalVector*>
const& x,
1101 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
1102 std::vector<double>& cache)
const override
1104 assert(x.size() == dof_table.size());
1106 auto const n_processes = x.size();
1107 std::vector<std::vector<double>> local_x;
1108 local_x.reserve(n_processes);
1110 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1112 auto const indices =
1114 assert(!indices.empty());
1115 local_x.push_back(x[process_id]->get(indices));
1119 if (n_processes == 1)
1121 auto const local_p = Eigen::Map<const NodalVectorType>(
1123 auto const local_C = Eigen::Map<const NodalVectorType>(
1130 constexpr
int pressure_process_id = 0;
1131 constexpr
int concentration_process_id = 1;
1132 auto const local_p = Eigen::Map<const NodalVectorType>(
1134 auto const local_C = Eigen::Map<const NodalVectorType>(
1142 Eigen::Ref<const NodalVectorType>
const& p_nodal_values,
1143 Eigen::Ref<const NodalVectorType>
const& C_nodal_values,
1144 std::vector<double>& cache)
const
1146 auto const n_integration_points =
1151 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1152 cache, GlobalDim, n_integration_points);
1159 auto const& medium =
1161 auto const& phase = medium.phase(
"AqueousLiquid");
1163 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
1165 auto const& ip_data =
_ip_data[ip];
1166 auto const& N = ip_data.N;
1167 auto const& dNdx = ip_data.dNdx;
1168 auto const& porosity = ip_data.porosity;
1170 pos.setIntegrationPoint(ip);
1172 double C_int_pt = 0.0;
1173 double p_int_pt = 0.0;
1178 vars[
static_cast<int>(
1180 vars[
static_cast<int>(
1187 double const dt = std::numeric_limits<double>::quiet_NaN();
1188 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1192 .template value<double>(vars, pos, t, dt);
1195 cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
1200 .template value<double>(vars, pos, t, dt);
1203 cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
1211 const unsigned integration_point)
const override
1213 auto const& N =
_ip_data[integration_point].N;
1216 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
1221 std::vector<double>
const& local_x)
const override
1223 auto const local_p = Eigen::Map<const NodalVectorType>(
1225 auto const local_C = Eigen::Map<const NodalVectorType>(
1231 auto const shape_matrices =
1235 std::array{pnt_local_coords})[0];
1242 auto const& medium =
1244 auto const& phase = medium.phase(
"AqueousLiquid");
1250 .emplace<double>(c_int_pt);
1255 .emplace<double>(p_int_pt);
1259 double const dt = std::numeric_limits<double>::quiet_NaN();
1260 auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1265 .template value<double>(vars, pos, t, dt);
1270 .template value<double>(vars, pos, t, dt);
1274 q += K_over_mu * rho_w * b;
1276 Eigen::Vector3d flux(0.0, 0.0, 0.0);
1277 flux.head<GlobalDim>() = rho_w *
q;
1282 std::vector<double>
const& local_x)
override
1284 unsigned const n_integration_points =
1287 std::vector<double> interpolated_values(n_integration_points);
1288 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1291 interpolated_values[ip]);
1293 return interpolated_values;
1298 std::vector<GlobalVector*>
const& int_pt_x,
1299 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
1300 std::vector<double>& cache)
const override
1303 assert(int_pt_x.size() == 1);
1307 unsigned const n_integration_points =
1309 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1311 auto const& chemical_system_id =
_ip_data[ip].chemical_system_id;
1312 auto const c_int_pt = int_pt_x[0]->get(chemical_system_id);
1313 cache.push_back(c_int_pt);
1322 Eigen::VectorXd
const& local_x,
1323 Eigen::VectorXd
const& )
override
1325 auto const local_p =
1327 auto const local_C = local_x.template segment<concentration_size>(
1330 std::vector<double> ele_velocity;
1333 auto const n_integration_points =
1335 auto const ele_velocity_mat =
1339 Eigen::Map<LocalVectorType>(
1342 ele_velocity_mat.rowwise().sum() / n_integration_points;
1348 auto const& medium =
1356 ip_data.porosity = ip_data.porosity_prev;
1360 medium, ip_data.porosity);
1364 std::vector<GlobalIndexType> chemical_system_indices;
1365 chemical_system_indices.reserve(n_integration_points);
1367 std::back_inserter(chemical_system_indices),
1368 [](
auto const& ip_data)
1369 { return ip_data.chemical_system_id; });
1372 ele_id, chemical_system_indices);
1377 double const ,
double const )
override
1379 unsigned const n_integration_points =
1382 for (
unsigned ip = 0; ip < n_integration_points; ip++)
1393 std::vector<std::reference_wrapper<ProcessVariable>>
const
1398 Eigen::aligned_allocator<
GlobalMatrix::IndexType GlobalIndexType
virtual void updatePorosityPostReaction(GlobalIndexType const &, MaterialPropertyLib::Medium const &, double &)
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 double getConcentration(int const, GlobalIndexType const) const
virtual void initializeChemicalSystemConcrete(std::vector< double > const &, GlobalIndexType const &, MaterialPropertyLib::Medium const &, ParameterLib::SpatialPosition const &, double const)
int add(IndexType row, IndexType col, double val)
Global vector based on Eigen vector.
void add(IndexType rowId, double v)
add entry
virtual std::size_t getID() const final
Returns the ID of the element.
MathLib::RowColumnIndices< GlobalIndexType > RowColumnIndices
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
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 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 & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
void setStaggeredCoupledSolutions(std::size_t const, CoupledSolutionsForStaggeredScheme *const coupling_term)
CoupledSolutionsForStaggeredScheme * _coupled_solutions
ComponentTransportLocalAssemblerInterface()=default
virtual std::vector< double > const & getInterpolatedLocalSolution(const double, std::vector< GlobalVector * > const &int_pt_x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const =0
std::vector< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType >, Eigen::aligned_allocator< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType > > > _ip_data
typename ShapeMatricesType::template MatrixType< pressure_size, pressure_size > LocalBlockMatrixType
static const int first_concentration_index
void assembleHydraulicEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
std::vector< double > interpolateNodalValuesToIntegrationPoints(std::vector< double > const &local_x) override
void assembleBlockMatrices(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::GlobalDimVectorType GlobalDimVectorType
void setChemicalSystemConcrete(Eigen::VectorXd const &local_x, double const t, double dt) 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 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
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
Eigen::Matrix< double, Eigen::Dynamic, 1 > LocalVectorType
void assembleComponentTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &, int const transport_process_id)
std::vector< double > const & getInterpolatedLocalSolution(const double, std::vector< GlobalVector * > const &int_pt_x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &cache) const override
Eigen::Vector3d getFlux(MathLib::Point3d const &pnt_local_coords, double const t, std::vector< double > const &local_x) const override
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > LocalMatrixType
typename ShapeMatricesType::NodalRowVectorType NodalRowVectorType
static const int concentration_size
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
typename ShapeMatricesType::template VectorType< pressure_size > LocalSegmentVectorType
IntegrationMethod const _integration_method
void initializeChemicalSystemConcrete(Eigen::VectorXd const &local_x, double const t) override
static const int pressure_index
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
LocalAssemblerData(MeshLib::Element const &element, std::size_t const local_matrix_size, bool is_axially_symmetric, unsigned const integration_order, ComponentTransportProcessData const &process_data, std::vector< std::reference_wrapper< ProcessVariable >> const &transport_process_variables)
typename ShapeMatricesType::NodalVectorType NodalVectorType
ComponentTransportProcessData const & _process_data
std::vector< std::reference_wrapper< ProcessVariable > > const _transport_process_variables
typename ShapeMatricesType::GlobalDimNodalMatrixType GlobalDimNodalMatrixType
static const int pressure_size
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
void postSpeciationCalculation(std::size_t const ele_id, double const t, double const dt) override
MeshLib::Element const & _element
void postTimestepConcrete(Eigen::VectorXd const &, double const, double const) override
void computeSecondaryVariableConcrete(double const t, double const, Eigen::VectorXd const &local_x, Eigen::VectorXd const &) override
void setChemicalSystemID(std::size_t const) override
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 assembleForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
@ 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)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > computeShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, PointContainer const &points)
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
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)
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.
Eigen::VectorXd const specific_body_force
ChemistryLib::ChemicalSolverInterface *const chemical_solver_interface
const int hydraulic_process_id
MeshLib::PropertyVector< double > * mesh_prop_velocity
bool const non_advective_form
ParameterLib::Parameter< double > const *const temperature
bool const chemically_induced_porosity_change
std::unique_ptr< MaterialPropertyLib::MaterialSpatialDistributionMap > media_map
GlobalDimNodalMatrixType const dNdx
GlobalIndexType chemical_system_id
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
NodalRowVectorType const N
IntegrationPointData(NodalRowVectorType const &N_, GlobalDimNodalMatrixType const &dNdx_, double const &integration_weight_)
double const integration_weight