25 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
27 assemble(
double const t,
double const dt,
28 std::vector<double>
const& local_x,
29 std::vector<double>
const& ,
30 std::vector<double>& local_M_data,
31 std::vector<double>& local_K_data,
32 std::vector<double>& local_b_data)
37 auto const& medium = *_process_data.media_map->getMedium(_element.getID());
41 .
template value<double>(vars, pos, t, dt);
43 std::numeric_limits<double>::quiet_NaN();
45 MaterialPropertyLib::formEigenTensor<ShapeFunction::DIM>(
56 assembleMatrixAndVector<IsotropicCalculator>(
57 t, dt, local_x, local_M_data, local_K_data, local_b_data);
61 assembleMatrixAndVector<AnisotropicCalculator>(
62 t, dt, local_x, local_M_data, local_K_data, local_b_data);
66 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
70 std::vector<double>
const& local_x)
const
74 double const dt = std::numeric_limits<double>::quiet_NaN();
78 auto const shape_matrices =
83 std::array{p_local_coords})[0];
89 auto const& medium = *_process_data.media_map->getMedium(_element.getID());
90 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
94 double pressure = 0.0;
100 MaterialPropertyLib::formEigenTensor<ShapeFunction::DIM>(
105 .template value<double>(vars, pos, t, dt);
107 Eigen::Vector3d flux(0.0, 0.0, 0.0);
108 auto const flux_local =
109 (-intrinsic_permeability /
viscosity * shape_matrices.dNdx *
110 Eigen::Map<const NodalVectorType>(local_x.data(), local_x.size()))
113 flux.head<GlobalDim>() =
114 _process_data.element_rotation_matrices[_element.getID()] * flux_local;
119 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
120 template <
typename LaplacianGravityVelocityCalculator>
123 std::vector<double>
const& local_x,
124 std::vector<double>& local_M_data,
125 std::vector<double>& local_K_data,
126 std::vector<double>& local_b_data)
128 auto const local_matrix_size = local_x.size();
129 assert(local_matrix_size == ShapeFunction::NPOINTS);
131 auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>(
132 local_M_data, local_matrix_size, local_matrix_size);
133 auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
134 local_K_data, local_matrix_size, local_matrix_size);
135 auto local_b = MathLib::createZeroedVector<NodalVectorType>(
136 local_b_data, local_matrix_size);
138 unsigned const n_integration_points =
139 _integration_method.getNumberOfPoints();
144 auto const& medium = *_process_data.media_map->getMedium(_element.getID());
145 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
150 .
template value<double>(vars, pos, t, dt);
153 _process_data.element_rotation_matrices[_element.getID()].transpose() *
154 _process_data.specific_body_force;
156 for (
unsigned ip = 0; ip < n_integration_points; ip++)
158 auto const& ip_data = _ip_data[ip];
168 .template value<double>(vars, pos, t, dt);
170 auto const ddensity_dpressure =
172 .template dValue<double>(
176 auto const porosity =
178 .template value<double>(vars, pos, t, dt);
180 .template value<double>(vars, pos, t, dt);
185 ip_data.N.transpose() * ip_data.N * ip_data.integration_weight;
190 .template value<double>(vars, pos, t, dt);
194 MaterialPropertyLib::formEigenTensor<ShapeFunction::DIM>(
199 LaplacianGravityVelocityCalculator::calculateLaplacianAndGravityTerm(
201 projected_body_force_vector, _process_data.has_gravity);
205 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
206 template <
typename VelocityCacheType>
209 const double dt, std::vector<double>
const& local_x,
211 VelocityCacheType& darcy_velocity_at_ips)
const
213 if (is_scalar_permeability)
215 computeProjectedDarcyVelocity<IsotropicCalculator>(
216 t, dt, local_x, pos, darcy_velocity_at_ips);
220 computeProjectedDarcyVelocity<AnisotropicCalculator>(
221 t, dt, local_x, pos, darcy_velocity_at_ips);
225 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
226 std::vector<double>
const&
230 std::vector<GlobalVector*>
const& x,
231 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
232 std::vector<double>& velocity_cache)
const
236 double const dt = std::numeric_limits<double>::quiet_NaN();
238 constexpr
int process_id = 0;
241 auto const local_x = x[process_id]->get(indices);
242 auto const n_integration_points = _integration_method.getNumberOfPoints();
243 velocity_cache.clear();
248 auto const& medium = *_process_data.media_map->getMedium(_element.getID());
252 .
template value<double>(vars, pos, t, dt);
254 std::numeric_limits<double>::quiet_NaN();
257 MaterialPropertyLib::formEigenTensor<ShapeFunction::DIM>(
261 assert(
permeability.rows() == _element.getDimension() ||
264 bool const is_scalar_permeability = (
permeability.size() == 1);
267 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
268 velocity_cache, GlobalDim, n_integration_points);
270 computeDarcyVelocity(is_scalar_permeability, t, dt, local_x, pos,
271 velocity_cache_vectors);
273 return velocity_cache;
276 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
277 template <
typename LaplacianGravityVelocityCalculator,
278 typename VelocityCacheType>
281 const double t,
const double dt, std::vector<double>
const& local_x,
283 VelocityCacheType& darcy_velocity_at_ips)
const
285 auto const local_matrix_size = local_x.size();
286 assert(local_matrix_size == ShapeFunction::NPOINTS);
288 const auto local_p_vec =
289 MathLib::toVector<NodalVectorType>(local_x, local_matrix_size);
291 unsigned const n_integration_points =
292 _integration_method.getNumberOfPoints();
294 auto const& medium = *_process_data.media_map->getMedium(_element.getID());
295 auto const& liquid_phase = medium.phase(
"AqueousLiquid");
300 .
template value<double>(vars, pos, t, dt);
303 _process_data.element_rotation_matrices[_element.getID()].transpose() *
304 _process_data.specific_body_force;
306 for (
unsigned ip = 0; ip < n_integration_points; ip++)
308 auto const& ip_data = _ip_data[ip];
317 .template value<double>(vars, pos, t, dt);
321 .template value<double>(vars, pos, t, dt);
324 MaterialPropertyLib::formEigenTensor<ShapeFunction::DIM>(
328 Eigen::VectorXd
const velocity_at_ip =
329 LaplacianGravityVelocityCalculator::calculateVelocity(
331 projected_body_force_vector, _process_data.has_gravity);
333 Eigen::MatrixXd
const& rotation_matrix =
334 _process_data.element_rotation_matrices[_element.getID()];
335 darcy_velocity_at_ips.col(ip) = rotation_matrix * velocity_at_ip;
339 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
342 Eigen::Map<NodalMatrixType>& local_K,
343 Eigen::Map<NodalVectorType>& local_b,
347 DimVectorType const& specific_body_force,
bool const has_gravity)
350 const double fac = K * ip_data.integration_weight;
351 local_K.noalias() += fac * ip_data.dNdx.transpose() * ip_data.dNdx;
356 (fac * rho_L) * ip_data.dNdx.transpose() * specific_body_force;
360 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
361 Eigen::Matrix<double, ShapeFunction::DIM, 1>
364 Eigen::Map<const NodalVectorType>
const& local_p,
368 DimVectorType const& specific_body_force,
bool const has_gravity)
376 velocity += (K * rho_L) * specific_body_force;
381 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
384 Eigen::Map<NodalMatrixType>& local_K,
385 Eigen::Map<NodalVectorType>& local_b,
389 DimVectorType const& specific_body_force,
bool const has_gravity)
391 const double fac = ip_data.integration_weight / mu;
393 fac * ip_data.dNdx.transpose() *
permeability * ip_data.dNdx;
397 local_b.noalias() += (fac * rho_L) * ip_data.dNdx.transpose() *
402 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
403 Eigen::Matrix<double, ShapeFunction::DIM, 1>
406 Eigen::Map<const NodalVectorType>
const& local_p,
410 DimVectorType const& specific_body_force,
bool const has_gravity)
418 velocity += (rho_L / mu) *
permeability * specific_body_force;
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
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 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::DimVectorType DimVectorType
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
typename ShapeMatricesType::NodalRowVectorType NodalRowVectorType
typename ShapeMatricesType::GlobalDimNodalMatrixType GlobalDimNodalMatrixType
Eigen::Vector3d getFlux(MathLib::Point3d const &p_local_coords, double const t, std::vector< double > const &local_x) const override
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::DimMatrixType DimMatrixType
void computeProjectedDarcyVelocity(const double t, const double dt, std::vector< double > const &local_x, ParameterLib::SpatialPosition const &pos, VelocityCacheType &darcy_velocity_at_ips) const
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)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > computeShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, PointContainer const &points)
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)
double fluid_density(const double p, const double T, const double x)
static void calculateLaplacianAndGravityTerm(Eigen::Map< NodalMatrixType > &local_K, Eigen::Map< NodalVectorType > &local_b, IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType > const &ip_data, DimMatrixType const &permeability, double const mu, double const rho_L, DimVectorType const &specific_body_force, bool const has_gravity)
static Eigen::Matrix< double, ShapeFunction::DIM, 1 > calculateVelocity(Eigen::Map< const NodalVectorType > const &local_p, IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType > const &ip_data, DimMatrixType const &permeability, double const mu, double const rho_L, DimVectorType const &specific_body_force, bool const has_gravity)
static void calculateLaplacianAndGravityTerm(Eigen::Map< NodalMatrixType > &local_K, Eigen::Map< NodalVectorType > &local_b, IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType > const &ip_data, DimMatrixType const &permeability, double const mu, double const rho_L, DimVectorType const &specific_body_force, bool const has_gravity)
static Eigen::Matrix< double, ShapeFunction::DIM, 1 > calculateVelocity(Eigen::Map< const NodalVectorType > const &local_p, IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType > const &ip_data, DimMatrixType const &permeability, double const mu, double const rho_L, DimVectorType const &specific_body_force, bool const has_gravity)