25template <
typename ShapeFunction,
int GlobalDim>
27 double const t,
double const dt, std::vector<double>
const& local_x,
28 std::vector<double>
const& ,
29 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
30 std::vector<double>& local_b_data)
35 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
39 .template value<double>(vars, pos, t, dt);
48 assert(permeability.rows() == GlobalDim || permeability.rows() == 1);
50 if (permeability.size() == 1)
52 assembleMatrixAndVector<IsotropicCalculator>(
53 t, dt, local_x, local_M_data, local_K_data, local_b_data);
57 assembleMatrixAndVector<AnisotropicCalculator>(
58 t, dt, local_x, local_M_data, local_K_data, local_b_data);
62template <
typename ShapeFunction,
int GlobalDim>
65 std::vector<double>
const& local_x)
const
69 double const dt = std::numeric_limits<double>::quiet_NaN();
73 auto const shape_matrices =
77 std::array{p_local_coords})[0];
83 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
84 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
88 double pressure = 0.0;
96 auto const viscosity =
98 .template value<double>(vars, pos, t, dt);
100 Eigen::Vector3d flux(0.0, 0.0, 0.0);
101 flux.head<GlobalDim>() =
102 -intrinsic_permeability / viscosity * shape_matrices.dNdx *
103 Eigen::Map<const NodalVectorType>(local_x.data(), local_x.size());
108template <
typename ShapeFunction,
int GlobalDim>
109template <
typename LaplacianGravityVelocityCalculator>
112 std::vector<double>
const& local_x,
113 std::vector<double>& local_M_data,
114 std::vector<double>& local_K_data,
115 std::vector<double>& local_b_data)
117 auto const local_matrix_size = local_x.size();
118 assert(local_matrix_size == ShapeFunction::NPOINTS);
121 local_M_data, local_matrix_size, local_matrix_size);
123 local_K_data, local_matrix_size, local_matrix_size);
125 local_b_data, local_matrix_size);
127 unsigned const n_integration_points =
128 _integration_method.getNumberOfPoints();
133 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
134 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
139 .template value<double>(vars, pos, t, dt);
142 _process_data.element_rotation_matrices[_element.getID()] *
143 _process_data.element_rotation_matrices[_element.getID()].transpose() *
144 _process_data.specific_body_force;
146 auto const& Ns = _process_data.shape_matrix_cache
147 .NsHigherOrder<
typename ShapeFunction::MeshElement>();
149 for (
unsigned ip = 0; ip < n_integration_points; ip++)
151 auto const& ip_data = _ip_data[ip];
152 auto const& N = Ns[ip];
159 auto const fluid_density =
161 .template value<double>(vars, pos, t, dt);
162 assert(fluid_density > 0.);
165 auto const ddensity_dpressure =
167 .template dValue<double>(
171 auto const porosity =
173 .template value<double>(vars, pos, t, dt);
175 .template value<double>(vars, pos, t, dt);
179 (porosity * ddensity_dpressure / fluid_density + storage) *
180 N.transpose() * N * ip_data.integration_weight;
183 auto const viscosity =
185 .template value<double>(vars, pos, t, dt);
194 LaplacianGravityVelocityCalculator::calculateLaplacianAndGravityTerm(
195 local_K, local_b, ip_data, permeability, viscosity, fluid_density,
196 projected_body_force_vector, _process_data.has_gravity);
200template <
typename ShapeFunction,
int GlobalDim>
201template <
typename VelocityCacheType>
203 bool const is_scalar_permeability,
const double t,
const double dt,
204 std::vector<double>
const& local_x,
206 VelocityCacheType& darcy_velocity_at_ips)
const
208 if (is_scalar_permeability)
210 computeProjectedDarcyVelocity<IsotropicCalculator>(
211 t, dt, local_x, pos, darcy_velocity_at_ips);
215 computeProjectedDarcyVelocity<AnisotropicCalculator>(
216 t, dt, local_x, pos, darcy_velocity_at_ips);
220template <
typename ShapeFunction,
int GlobalDim>
221std::vector<double>
const&
224 std::vector<GlobalVector*>
const& x,
225 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
226 std::vector<double>& velocity_cache)
const
230 double const dt = std::numeric_limits<double>::quiet_NaN();
232 constexpr int process_id = 0;
235 auto const local_x = x[process_id]->get(indices);
236 auto const n_integration_points = _integration_method.getNumberOfPoints();
237 velocity_cache.clear();
242 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
246 .template value<double>(vars, pos, t, dt);
254 assert(permeability.rows() == GlobalDim || permeability.rows() == 1);
256 bool const is_scalar_permeability = (permeability.size() == 1);
259 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
260 velocity_cache, GlobalDim, n_integration_points);
262 computeDarcyVelocity(is_scalar_permeability, t, dt, local_x, pos,
263 velocity_cache_vectors);
265 return velocity_cache;
268template <
typename ShapeFunction,
int GlobalDim>
269template <
typename LaplacianGravityVelocityCalculator,
270 typename VelocityCacheType>
273 const double t,
const double dt, std::vector<double>
const& local_x,
275 VelocityCacheType& darcy_velocity_at_ips)
const
277 auto const local_matrix_size = local_x.size();
278 assert(local_matrix_size == ShapeFunction::NPOINTS);
280 const auto local_p_vec =
283 unsigned const n_integration_points =
284 _integration_method.getNumberOfPoints();
286 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
287 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
292 .template value<double>(vars, pos, t, dt);
295 _process_data.element_rotation_matrices[_element.getID()] *
296 _process_data.element_rotation_matrices[_element.getID()].transpose() *
297 _process_data.specific_body_force;
299 auto const& Ns = _process_data.shape_matrix_cache
300 .NsHigherOrder<
typename ShapeFunction::MeshElement>();
302 for (
unsigned ip = 0; ip < n_integration_points; ip++)
304 auto const& ip_data = _ip_data[ip];
305 auto const& N = Ns[ip];
311 auto const fluid_density =
313 .template value<double>(vars, pos, t, dt);
317 auto const viscosity =
319 .template value<double>(vars, pos, t, dt);
326 darcy_velocity_at_ips.col(ip) =
327 LaplacianGravityVelocityCalculator::calculateVelocity(
328 local_p_vec, ip_data, permeability, viscosity, fluid_density,
329 projected_body_force_vector, _process_data.has_gravity);
333template <
typename ShapeFunction,
int GlobalDim>
336 Eigen::Map<NodalMatrixType>& local_K,
337 Eigen::Map<NodalVectorType>& local_b,
341 bool const has_gravity)
343 const double K = permeability(0, 0) / mu;
345 local_K.noalias() += fac * ip_data.
dNdx.transpose() * ip_data.
dNdx;
350 (fac * rho_L) * ip_data.
dNdx.transpose() * specific_body_force;
354template <
typename ShapeFunction,
int GlobalDim>
355Eigen::Matrix<double, GlobalDim, 1>
358 Eigen::Map<const NodalVectorType>
const& local_p,
362 bool const has_gravity)
364 const double K = permeability(0, 0) / mu;
370 velocity += (K * rho_L) * specific_body_force;
375template <
typename ShapeFunction,
int GlobalDim>
378 Eigen::Map<NodalMatrixType>& local_K,
379 Eigen::Map<NodalVectorType>& local_b,
383 bool const has_gravity)
387 fac * ip_data.
dNdx.transpose() * permeability * ip_data.
dNdx;
391 local_b.noalias() += (fac * rho_L) * ip_data.
dNdx.transpose() *
392 permeability * specific_body_force;
396template <
typename ShapeFunction,
int GlobalDim>
397Eigen::Matrix<double, GlobalDim, 1>
400 Eigen::Map<const NodalVectorType>
const& local_p,
404 bool const has_gravity)
412 velocity += (rho_L / mu) * permeability * specific_body_force;
double liquid_phase_pressure
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
Eigen::Vector3d getFlux(MathLib::Point3d const &p_local_coords, double const t, std::vector< double > const &local_x) const override
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
void computeProjectedDarcyVelocity(const double t, const double dt, std::vector< double > const &local_x, ParameterLib::SpatialPosition const &pos, VelocityCacheType &darcy_velocity_at_ips) const
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 assembleMatrixAndVector(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
void computeDarcyVelocity(bool const is_scalar_permeability, const double t, const double dt, std::vector< double > const &local_x, ParameterLib::SpatialPosition const &pos, VelocityCacheType &darcy_velocity_at_ips) const
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 > &velocity_cache) const override
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
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)
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 > > computeShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, PointContainer const &points)
GlobalDimNodalMatrixType const dNdx
double const integration_weight
static void calculateLaplacianAndGravityTerm(Eigen::Map< NodalMatrixType > &local_K, Eigen::Map< NodalVectorType > &local_b, IntegrationPointData< GlobalDimNodalMatrixType > const &ip_data, GlobalDimMatrixType const &permeability, double const mu, double const rho_L, GlobalDimVectorType const &specific_body_force, bool const has_gravity)
static Eigen::Matrix< double, GlobalDim, 1 > calculateVelocity(Eigen::Map< const NodalVectorType > const &local_p, IntegrationPointData< GlobalDimNodalMatrixType > const &ip_data, GlobalDimMatrixType const &permeability, double const mu, double const rho_L, GlobalDimVectorType const &specific_body_force, bool const has_gravity)
static void calculateLaplacianAndGravityTerm(Eigen::Map< NodalMatrixType > &local_K, Eigen::Map< NodalVectorType > &local_b, IntegrationPointData< GlobalDimNodalMatrixType > const &ip_data, GlobalDimMatrixType const &permeability, double const mu, double const rho_L, GlobalDimVectorType const &specific_body_force, bool const has_gravity)
static Eigen::Matrix< double, GlobalDim, 1 > calculateVelocity(Eigen::Map< const NodalVectorType > const &local_p, IntegrationPointData< GlobalDimNodalMatrixType > const &ip_data, GlobalDimMatrixType const &permeability, double const mu, double const rho_L, GlobalDimVectorType const &specific_body_force, bool const has_gravity)