62 double const t,
double const dt, std::vector<double>
const& local_x,
63 std::vector<double>
const& ,
64 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
65 std::vector<double>& local_b_data)
67 auto const local_matrix_size = local_x.size();
70 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
73 local_M_data, local_matrix_size, local_matrix_size);
75 local_K_data, local_matrix_size, local_matrix_size);
77 local_b_data, local_matrix_size);
79 unsigned const n_integration_points =
80 _integration_method.getNumberOfPoints();
82 auto p_nodal_values = Eigen::Map<const NodalVectorType>(
83 &local_x[pressure_index], pressure_size);
85 auto const& b = _process_data.specific_body_force;
90 GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
93 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
94 auto const& phase = medium.phase(
"AqueousLiquid");
95 auto const& component =
96 phase.component(_transport_process_variable.getName());
98 auto KCC = local_K.template block<concentration_size, concentration_size>(
99 concentration_index, concentration_index);
100 auto MCC = local_M.template block<concentration_size, concentration_size>(
101 concentration_index, concentration_index);
102 auto Kpp = local_K.template block<pressure_size, pressure_size>(
103 pressure_index, pressure_index);
104 auto Mpp = local_M.template block<pressure_size, pressure_size>(
105 pressure_index, pressure_index);
106 auto Bp = local_b.template block<pressure_size, 1>(pressure_index, 0);
108 for (std::size_t ip(0); ip < n_integration_points; ++ip)
110 auto const& ip_data = _ip_data[ip];
111 auto const& N = ip_data.N;
112 auto const& dNdx = ip_data.dNdx;
113 auto const& w = ip_data.integration_weight;
116 std::nullopt, _element.getID(),
121 double C_int_pt = 0.0;
122 double p_int_pt = 0.0;
128 .template value<double>(vars, pos, t, dt);
130 double const dSw_dpc =
132 .template dValue<double>(
144 auto const specific_storage =
146 .template value<double>(vars, pos, t, dt);
149 auto const porosity =
151 .template value<double>(vars, pos, t, dt);
153 auto const retardation_factor =
155 .template value<double>(vars, pos, t, dt);
157 auto const solute_dispersivity_transverse =
159 .template value<double>(vars, pos, t, dt);
160 auto const solute_dispersivity_longitudinal =
162 .template value<double>(vars, pos, t, dt);
166 .template value<double>(vars, pos, t, dt);
169 auto const decay_rate =
171 .template value<double>(vars, pos, t, dt);
172 auto const pore_diffusion_coefficient =
175 .value(vars, pos, t, dt));
183 .template value<double>(vars, pos, t, dt);
186 .template value<double>(vars, pos, t, dt);
187 auto const K_times_k_rel_over_mu = K * (k_rel / mu);
190 _process_data.has_gravity
192 (dNdx * p_nodal_values - density * b))
196 double const velocity_magnitude = velocity.norm();
198 velocity_magnitude != 0.0
200 porosity * pore_diffusion_coefficient +
201 solute_dispersivity_transverse * velocity_magnitude * I +
202 (solute_dispersivity_longitudinal -
203 solute_dispersivity_transverse) /
204 velocity_magnitude * velocity * velocity.transpose())
206 solute_dispersivity_transverse *
207 velocity_magnitude * I);
211 (dNdx.transpose() * hydrodynamic_dispersion * dNdx +
212 N.transpose() * velocity.transpose() * dNdx +
213 N.transpose() * decay_rate * porosity * retardation_factor * N) *
215 MCC.noalias() += w * N.transpose() * porosity * retardation_factor * N;
216 Kpp.noalias() += w * dNdx.transpose() * K_times_k_rel_over_mu * dNdx;
218 double const drhow_dp(0.0);
219 Mpp.noalias() += (specific_storage * Sw + porosity * Sw * drhow_dp -
220 porosity * dSw_dpc) *
221 ip_data.mass_operator;
223 if (_process_data.has_gravity)
225 Bp += w * density * dNdx.transpose() * K_times_k_rel_over_mu * b;
236 std::vector<GlobalVector*>
const& x,
237 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
238 std::vector<double>& cache)
const
240 unsigned const n_integration_points =
241 _integration_method.getNumberOfPoints();
243 constexpr int process_id = 0;
246 assert(!indices.empty());
247 auto const local_x = x[process_id]->get(indices);
251 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
252 cache, GlobalDim, n_integration_points);
257 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
258 auto const& phase = medium.phase(
"AqueousLiquid");
260 auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
261 &local_x[ShapeFunction::NPOINTS], ShapeFunction::NPOINTS);
263 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
265 auto const& ip_data = _ip_data[ip];
266 auto const& N = ip_data.N;
267 auto const& dNdx = ip_data.dNdx;
270 std::nullopt, _element.getID(),
275 auto const dt = std::numeric_limits<double>::quiet_NaN();
280 .template value<double>(vars, pos, t, dt);
282 double C_int_pt = 0.0;
283 double p_int_pt = 0.0;
287 vars.capillary_pressure = -p_int_pt;
289 .template value<double>(vars, pos, t, dt);
291 vars.liquid_saturation = Sw;
294 .template value<double>(vars, pos, t, dt);
296 cache_mat.col(ip).noalias() = -dNdx * p_nodal_values;
297 if (_process_data.has_gravity)
299 vars.concentration = C_int_pt;
300 vars.liquid_phase_pressure = p_int_pt;
302 .template value<double>(vars, pos, t, dt);
303 auto const b = _process_data.specific_body_force;
305 cache_mat.col(ip).noalias() += rho_w * b;
307 cache_mat.col(ip).noalias() = k_rel / mu * (K * cache_mat.col(ip));
327 std::vector<GlobalVector*>
const& x,
328 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
329 std::vector<double>& cache)
const
333 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
335 unsigned const n_integration_points =
336 _integration_method.getNumberOfPoints();
338 constexpr int process_id = 0;
341 assert(!indices.empty());
342 auto const local_x = x[process_id]->get(indices);
346 Eigen::Matrix<double, 1, Eigen::Dynamic, Eigen::RowMajor>>(
347 cache, n_integration_points);
349 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
351 auto const& ip_data = _ip_data[ip];
352 auto const& N = ip_data.N;
355 std::nullopt, _element.getID(),
360 double C_int_pt = 0.0;
361 double p_int_pt = 0.0;
366 auto const dt = std::numeric_limits<double>::quiet_NaN();
368 .template value<double>(vars, pos, t, dt);