34 namespace RichardsFlow
36 template <
typename NodalRowVectorType,
typename GlobalDimNodalMatrixType,
37 typename NodalMatrixType>
41 GlobalDimNodalMatrixType
const& dNdx_,
42 double const& integration_weight_,
43 NodalMatrixType
const mass_operator_)
51 NodalRowVectorType
const N;
52 GlobalDimNodalMatrixType
const dNdx;
67 std::vector<GlobalVector*>
const& x,
68 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
69 std::vector<double>& cache)
const = 0;
73 std::vector<GlobalVector*>
const& x,
74 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
75 std::vector<double>& cache)
const = 0;
78 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
97 std::size_t
const local_matrix_size,
98 bool is_axially_symmetric,
99 unsigned const integration_order,
109 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
110 (void)local_matrix_size;
112 unsigned const n_integration_points =
114 _ip_data.reserve(n_integration_points);
116 auto const shape_matrices =
118 GlobalDim>(element, is_axially_symmetric,
121 for (
unsigned ip = 0; ip < n_integration_points; ip++)
123 auto const& sm = shape_matrices[ip];
124 const double integration_factor = sm.integralMeasure * sm.detJ;
129 sm.N.transpose() * sm.N * integration_factor *
135 std::vector<double>
const& local_x,
136 std::vector<double>
const& ,
137 std::vector<double>& local_M_data,
138 std::vector<double>& local_K_data,
139 std::vector<double>& local_b_data)
override
141 auto const local_matrix_size = local_x.size();
144 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
146 auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>(
147 local_M_data, local_matrix_size, local_matrix_size);
148 auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
149 local_K_data, local_matrix_size, local_matrix_size);
150 auto local_b = MathLib::createZeroedVector<NodalVectorType>(
151 local_b_data, local_matrix_size);
153 unsigned const n_integration_points =
160 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
166 .template value<double>(vars, pos, t, dt);
168 for (
unsigned ip = 0; ip < n_integration_points; ip++)
171 double p_int_pt = 0.0;
174 vars[
static_cast<int>(
176 vars[
static_cast<int>(
180 MaterialPropertyLib::formEigenTensor<GlobalDim>(
182 .value(vars, pos, t, dt));
184 auto const porosity =
186 .template value<double>(vars, pos, t, dt);
191 .template value<double>(vars, pos, t, dt);
193 vars[
static_cast<int>(
196 double const dSw_dpc =
199 .template dValue<double>(
203 auto const drhow_dp =
206 .template dValue<double>(
211 .template value<double>(vars, pos, t, dt);
212 double const mass_mat_coeff =
213 storage * Sw + porosity * Sw * drhow_dp - porosity * dSw_dpc;
215 local_M.noalias() += mass_mat_coeff *
_ip_data[ip].mass_operator;
221 .template value<double>(vars, pos, t, dt);
225 .template value<double>(vars, pos, t, dt);
228 _ip_data[ip].integration_weight * (k_rel / mu);
235 .template value<double>(vars, pos, t, dt);
237 assert(body_force.size() == GlobalDim);
241 local_b.noalias() += (k_rel / mu) * rho_w * gravity_operator;
246 for (
int idx_ml = 0; idx_ml < local_M.cols(); idx_ml++)
248 double const mass_lump_val = local_M.col(idx_ml).sum();
249 local_M.col(idx_ml).setZero();
250 local_M(idx_ml, idx_ml) = mass_lump_val;
256 const unsigned integration_point)
const override
258 auto const& N =
_ip_data[integration_point].N;
261 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
266 std::vector<GlobalVector*>
const& ,
267 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
268 std::vector<double>& )
const override
276 std::vector<GlobalVector*>
const& x,
277 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
278 std::vector<double>& cache)
const override
282 double const dt = std::numeric_limits<double>::quiet_NaN();
284 constexpr
int process_id = 0;
287 assert(!indices.empty());
288 auto const local_x = x[process_id]->get(indices);
295 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
302 .template value<double>(vars, pos, t, dt);
304 auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
305 &local_x[0], ShapeFunction::NPOINTS);
307 unsigned const n_integration_points =
312 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
313 cache, GlobalDim, n_integration_points);
315 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
317 double p_int_pt = 0.0;
319 vars[
static_cast<int>(
321 vars[
static_cast<int>(
327 .template value<double>(vars, pos, t, dt);
328 vars[
static_cast<int>(
332 MaterialPropertyLib::formEigenTensor<GlobalDim>(
334 .value(vars, pos, t, dt));
340 .template value<double>(vars, pos, t, dt);
343 .template value<double>(vars, pos, t, dt);
345 cache_vec.col(ip).noalias() =
346 -K_mat_coeff *
_ip_data[ip].dNdx * p_nodal_values;
352 .template value<double>(vars, pos, t, dt);
354 assert(body_force.size() == GlobalDim);
357 cache_vec.col(ip).noalias() += K_mat_coeff * rho_w * body_force;
Definition of the PiecewiseLinearInterpolation class.
virtual std::size_t getID() const final
Returns the ID of the element.
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
std::vector< double > _saturation
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::NodalVectorType NodalVectorType
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
IntegrationMethod const _integration_method
typename ShapeMatricesType::NodalMatrixType NodalMatrixType
RichardsFlowProcessData const & _process_data
LocalAssemblerData(MeshLib::Element const &element, std::size_t const local_matrix_size, bool is_axially_symmetric, unsigned const integration_order, RichardsFlowProcessData const &process_data)
MeshLib::Element const & _element
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
std::vector< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType >, Eigen::aligned_allocator< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType > > > _ip_data
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
typename ShapeMatricesType::GlobalDimNodalMatrixType GlobalDimNodalMatrixType
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
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
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 & getIntPtSaturation(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
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)
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
std::unique_ptr< MaterialPropertyLib::MaterialSpatialDistributionMap > media_map
bool const has_mass_lumping
Eigen::VectorXd const specific_body_force