20 namespace RichardsComponentTransport
22 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
26 std::size_t
const local_matrix_size,
27 bool is_axially_symmetric,
28 unsigned const integration_order,
31 : _element_id(element.getID()),
32 _process_data(process_data),
33 _integration_method(integration_order),
34 _transport_process_variable(transport_process_variable)
38 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
39 (void)local_matrix_size;
41 unsigned const n_integration_points =
43 _ip_data.reserve(n_integration_points);
45 auto const shape_matrices =
46 NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType, GlobalDim>(
49 for (
unsigned ip = 0; ip < n_integration_points; ip++)
51 auto const& sm = shape_matrices[ip];
52 double const integration_factor = sm.integralMeasure * sm.detJ;
57 sm.N.transpose() * sm.N * integration_factor *
62 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
64 double const t,
double const dt, std::vector<double>
const& local_x,
65 std::vector<double>
const& ,
66 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
67 std::vector<double>& local_b_data)
69 auto const local_matrix_size = local_x.size();
72 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
74 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
75 local_M_data, local_matrix_size, local_matrix_size);
76 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
77 local_K_data, local_matrix_size, local_matrix_size);
78 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
79 local_b_data, local_matrix_size);
81 unsigned const n_integration_points =
82 _integration_method.getNumberOfPoints();
87 auto p_nodal_values = Eigen::Map<const NodalVectorType>(
88 &local_x[pressure_index], pressure_size);
90 auto const& b = _process_data.specific_body_force;
95 GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
98 auto const& medium = *_process_data.media_map->getMedium(_element_id);
99 auto const& phase = medium.phase(
"AqueousLiquid");
100 auto const& component =
101 phase.component(_transport_process_variable.getName());
103 auto KCC = local_K.template block<concentration_size, concentration_size>(
104 concentration_index, concentration_index);
105 auto MCC = local_M.template block<concentration_size, concentration_size>(
106 concentration_index, concentration_index);
107 auto Kpp = local_K.template block<pressure_size, pressure_size>(
108 pressure_index, pressure_index);
109 auto Mpp = local_M.template block<pressure_size, pressure_size>(
110 pressure_index, pressure_index);
111 auto Bp = local_b.template block<pressure_size, 1>(pressure_index, 0);
113 for (std::size_t ip(0); ip < n_integration_points; ++ip)
117 auto const& ip_data = _ip_data[ip];
118 auto const& N = ip_data.N;
119 auto const& dNdx = ip_data.dNdx;
120 auto const& w = ip_data.integration_weight;
122 double C_int_pt = 0.0;
123 double p_int_pt = 0.0;
127 vars[
static_cast<int>(
130 .template value<double>(vars, pos, t, dt);
132 double const dSw_dpc =
134 .template dValue<double>(
145 auto const specific_storage =
147 .template value<double>(vars, pos, t, dt);
150 auto const porosity =
152 .template value<double>(vars, pos, t, dt);
156 .template value<double>(vars, pos, t, dt);
158 auto const solute_dispersivity_transverse =
160 .template value<double>(vars, pos, t, dt);
161 auto const solute_dispersivity_longitudinal =
163 .template value<double>(vars, pos, t, dt);
167 .template value<double>(vars, pos, t, dt);
170 .template value<double>(vars, pos, t, dt);
171 auto const pore_diffusion_coefficient =
172 MaterialPropertyLib::formEigenTensor<GlobalDim>(
174 .value(vars, pos, t, dt));
176 auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
179 vars[
static_cast<int>(
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 +
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;
232 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
233 std::vector<double>
const&
237 std::vector<GlobalVector*>
const& x,
238 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
239 std::vector<double>& cache)
const
241 unsigned const n_integration_points =
242 _integration_method.getNumberOfPoints();
244 constexpr
int process_id = 0;
247 assert(!indices.empty());
248 auto const local_x = x[process_id]->get(indices);
252 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
253 cache, GlobalDim, n_integration_points);
261 auto const& medium = *_process_data.media_map->getMedium(_element_id);
262 auto const& phase = medium.phase(
"AqueousLiquid");
264 auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
265 &local_x[ShapeFunction::NPOINTS], ShapeFunction::NPOINTS);
267 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
269 auto const& ip_data = _ip_data[ip];
270 auto const& N = ip_data.N;
271 auto const& dNdx = ip_data.dNdx;
273 pos.setIntegrationPoint(ip);
275 auto const dt = std::numeric_limits<double>::quiet_NaN();
276 auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
280 .template value<double>(vars, pos, t, dt);
282 double C_int_pt = 0.0;
283 double p_int_pt = 0.0;
287 vars[
static_cast<int>(
290 .template value<double>(vars, pos, t, dt);
292 vars[
static_cast<int>(
296 .template value<double>(vars, pos, t, dt);
298 cache_mat.col(ip).noalias() = -dNdx * p_nodal_values;
299 if (_process_data.has_gravity)
301 vars[
static_cast<int>(
303 vars[
static_cast<int>(
306 .template value<double>(vars, pos, t, dt);
307 auto const b = _process_data.specific_body_force;
309 cache_mat.col(ip).noalias() += rho_w * b;
311 cache_mat.col(ip).noalias() = k_rel / mu * (K * cache_mat.col(ip));
316 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
317 Eigen::Map<const Eigen::RowVectorXd>
319 const unsigned integration_point)
const
321 auto const& N = _ip_data[integration_point].N;
324 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
327 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
328 std::vector<double>
const&
332 std::vector<GlobalVector*>
const& x,
333 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
334 std::vector<double>& cache)
const
341 auto const& medium = *_process_data.media_map->getMedium(_element_id);
343 unsigned const n_integration_points =
344 _integration_method.getNumberOfPoints();
346 constexpr
int process_id = 0;
349 assert(!indices.empty());
350 auto const local_x = x[process_id]->get(indices);
353 auto cache_vec = MathLib::createZeroedVector<LocalVectorType>(
354 cache, n_integration_points);
356 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
358 auto const& ip_data = _ip_data[ip];
359 auto const& N = ip_data.N;
361 double C_int_pt = 0.0;
362 double p_int_pt = 0.0;
366 vars[
static_cast<int>(
368 auto const dt = std::numeric_limits<double>::quiet_NaN();
370 .template value<double>(vars, pos, t, dt);
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
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
std::vector< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType >, Eigen::aligned_allocator< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType > > > _ip_data
void assemble(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_xdot, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
LocalAssemblerData(MeshLib::Element const &element, std::size_t const local_matrix_size, bool is_axially_symmetric, unsigned const integration_order, RichardsComponentTransportProcessData const &process_data, ProcessVariable const &transport_process_variable)
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
IntegrationMethod const _integration_method
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
std::vector< double > const & getIntPtSaturation(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const 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< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
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)
const unsigned NUM_NODAL_DOF