37template <
typename NodalRowVectorType,
typename GlobalDimNodalMatrixType,
38 typename NodalMatrixType>
42 GlobalDimNodalMatrixType
const& dNdx_,
43 double const& integration_weight_,
44 NodalMatrixType
const mass_operator_)
52 NodalRowVectorType
const N;
53 GlobalDimNodalMatrixType
const dNdx;
68 std::vector<GlobalVector*>
const& x,
69 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
70 std::vector<double>& cache)
const = 0;
74 std::vector<GlobalVector*>
const& x,
75 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
76 std::vector<double>& cache)
const = 0;
79template <
typename ShapeFunction,
int GlobalDim>
99 std::size_t
const local_matrix_size,
101 bool is_axially_symmetric,
111 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
112 (void)local_matrix_size;
114 unsigned const n_integration_points =
116 _ip_data.reserve(n_integration_points);
118 auto const shape_matrices =
120 GlobalDim>(element, is_axially_symmetric,
123 for (
unsigned ip = 0; ip < n_integration_points; ip++)
125 auto const& sm = shape_matrices[ip];
126 const double integration_factor = sm.integralMeasure * sm.detJ;
131 sm.N.transpose() * sm.N * integration_factor *
137 std::vector<double>
const& local_x,
138 std::vector<double>
const& ,
139 std::vector<double>& local_M_data,
140 std::vector<double>& local_K_data,
141 std::vector<double>& local_b_data)
override
143 auto const local_matrix_size = local_x.size();
146 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
148 auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>(
149 local_M_data, local_matrix_size, local_matrix_size);
150 auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
151 local_K_data, local_matrix_size, local_matrix_size);
152 auto local_b = MathLib::createZeroedVector<NodalVectorType>(
153 local_b_data, local_matrix_size);
155 unsigned const n_integration_points =
162 auto const& liquid_phase = medium.
phase(
"AqueousLiquid");
168 .template value<double>(vars, pos, t, dt);
170 for (
unsigned ip = 0; ip < n_integration_points; ip++)
173 double p_int_pt = 0.0;
179 auto const permeability =
180 MaterialPropertyLib::formEigenTensor<GlobalDim>(
182 .value(vars, pos, t, dt));
184 auto const porosity =
186 .template value<double>(vars, pos, t, dt);
190 .template value<double>(vars, pos, t, dt);
194 double const dSw_dpc =
196 .template dValue<double>(
200 auto const drhow_dp =
203 .template dValue<double>(
209 .template value<double>(vars, pos, t, dt);
210 double const mass_mat_coeff =
211 storage * Sw + porosity * Sw * drhow_dp - porosity * dSw_dpc;
213 local_M.noalias() += mass_mat_coeff *
_ip_data[ip].mass_operator;
218 relative_permeability)
219 .template value<double>(vars, pos, t, dt);
223 .template value<double>(vars, pos, t, dt);
224 local_K.noalias() +=
_ip_data[ip].dNdx.transpose() * permeability *
226 _ip_data[ip].integration_weight * (k_rel / mu);
233 .template value<double>(vars, pos, t, dt);
235 assert(body_force.size() == GlobalDim);
237 _ip_data[ip].dNdx.transpose() * permeability * body_force *
239 local_b.noalias() += (k_rel / mu) * rho_w * gravity_operator;
244 for (
int idx_ml = 0; idx_ml < local_M.cols(); idx_ml++)
246 double const mass_lump_val = local_M.col(idx_ml).sum();
247 local_M.col(idx_ml).setZero();
248 local_M(idx_ml, idx_ml) = mass_lump_val;
254 const unsigned integration_point)
const override
256 auto const& N =
_ip_data[integration_point].N;
259 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
264 std::vector<GlobalVector*>
const& ,
265 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
266 std::vector<double>& )
const override
274 std::vector<GlobalVector*>
const& x,
275 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
276 std::vector<double>& cache)
const override
280 double const dt = std::numeric_limits<double>::quiet_NaN();
282 constexpr int process_id = 0;
285 assert(!indices.empty());
286 auto const local_x = x[process_id]->get(indices);
293 auto const& liquid_phase = medium.
phase(
"AqueousLiquid");
300 .template value<double>(vars, pos, t, dt);
302 auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
303 &local_x[0], ShapeFunction::NPOINTS);
305 unsigned const n_integration_points =
310 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
311 cache, GlobalDim, n_integration_points);
313 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
315 double p_int_pt = 0.0;
322 .template value<double>(vars, pos, t, dt);
325 auto const permeability =
326 MaterialPropertyLib::formEigenTensor<GlobalDim>(
328 .value(vars, pos, t, dt));
333 relative_permeability)
334 .template value<double>(vars, pos, t, dt);
337 .template value<double>(vars, pos, t, dt);
338 auto const K_mat_coeff = permeability * (k_rel / mu);
339 cache_vec.col(ip).noalias() =
340 -K_mat_coeff *
_ip_data[ip].dNdx * p_nodal_values;
346 .template value<double>(vars, pos, t, dt);
348 assert(body_force.size() == GlobalDim);
351 cache_vec.col(ip).noalias() += K_mat_coeff * rho_w * body_force;
Definition of the PiecewiseLinearInterpolation class.
Medium * getMedium(std::size_t element_id)
Phase const & phase(std::size_t index) const
double capillary_pressure
double liquid_phase_pressure
virtual std::size_t getID() const final
Returns the ID of the element.
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
unsigned getNumberOfPoints() const
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
NumLib::GenericIntegrationMethod const & _integration_method
std::vector< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType >, Eigen::aligned_allocator< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType > > > _ip_data
typename ShapeMatricesType::GlobalDimNodalMatrixType GlobalDimNodalMatrixType
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
std::vector< double > _saturation
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< double > const & getIntPtSaturation(const double, std::vector< GlobalVector * > const &, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &) const override
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
RichardsFlowProcessData const & _process_data
LocalAssemblerData(MeshLib::Element const &element, std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const &integration_method, bool is_axially_symmetric, RichardsFlowProcessData const &process_data)
typename ShapeMatricesType::NodalRowVectorType NodalRowVectorType
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
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
typename ShapeMatricesType::NodalMatrixType NodalMatrixType
MeshLib::Element const & _element
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
typename ShapeMatricesType::NodalVectorType NodalVectorType
virtual 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 =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
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
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)
const unsigned NUM_NODAL_DOF
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
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
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
GlobalDimNodalMatrixType const dNdx
double const integration_weight
NodalRowVectorType const N
IntegrationPointData(NodalRowVectorType const &N_, GlobalDimNodalMatrixType const &dNdx_, double const &integration_weight_, NodalMatrixType const mass_operator_)
NodalMatrixType const mass_operator
MaterialPropertyLib::MaterialSpatialDistributionMap media_map
bool const has_mass_lumping
Eigen::VectorXd const specific_body_force