13 #include <Eigen/Dense>
31 template <
typename ShapeFunctionLiquidVelocity,
typename ShapeFunctionPressure,
32 typename IntegrationMethod,
int GlobalDim>
37 ShapeFunctionLiquidVelocity::NPOINTS * GlobalDim;
40 ShapeFunctionLiquidVelocity::NPOINTS * GlobalDim;
53 bool const is_axially_symmetric,
54 unsigned const integration_order,
61 unsigned const n_integration_points =
63 _ip_data.reserve(n_integration_points);
65 auto const shape_matrices_v =
70 auto const shape_matrices_p =
75 for (
unsigned ip = 0; ip < n_integration_points; ip++)
78 shape_matrices_v[ip].N, shape_matrices_v[ip].dNdx,
79 shape_matrices_p[ip].N, shape_matrices_p[ip].dNdx,
81 shape_matrices_v[ip].integralMeasure *
82 shape_matrices_v[ip].detJ);
86 ip_data.N_v_op = ShapeMatrixTypeLiquidVelocity::template MatrixType<
91 ShapeMatrixTypeLiquidVelocity::template MatrixType<
92 GlobalDim * GlobalDim,
96 for (
int i = 0; i < GlobalDim; ++i)
101 .noalias() = shape_matrices_v[ip].N;
104 .template block<GlobalDim,
107 .noalias() = shape_matrices_v[ip].dNdx;
113 std::vector<double>
const& local_x,
114 std::vector<double>
const& ,
115 std::vector<double>& ,
116 std::vector<double>& local_K_data,
117 std::vector<double>& local_b_data)
override
120 assert(local_x.size() == local_matrix_size);
122 auto local_p = Eigen::Map<const NodalVectorType>(
126 typename ShapeMatrixTypeLiquidVelocity::template MatrixType<
127 local_matrix_size, local_matrix_size>>(
128 local_K_data, local_matrix_size, local_matrix_size);
130 typename ShapeMatrixTypeLiquidVelocity::template VectorType<
131 local_matrix_size>>(local_b_data, local_matrix_size);
136 local_K.template block<liquid_velocity_size, liquid_velocity_size>(
139 auto Kvp = local_K.template block<liquid_velocity_size, pressure_size>(
142 auto Kpv = local_K.template block<pressure_size, liquid_velocity_size>(
145 auto bv = local_b.template segment<liquid_velocity_size>(
148 unsigned const n_integration_points =
159 assert(GlobalDim == 2);
160 auto const identity_mat =
161 Eigen::Matrix<double, GlobalDim, GlobalDim>::Identity().eval();
162 auto const unit_vector = Eigen::Map<const Eigen::RowVectorXd>(
163 identity_mat.data(), identity_mat.size())
169 auto const& phase = medium.phase(
"AqueousLiquid");
171 for (
unsigned ip(0); ip < n_integration_points; ++ip)
177 auto const& N_p = ip_data.N_p;
178 auto const& N_v_op = ip_data.N_v_op;
180 auto const& dNdx_v = ip_data.dNdx_v;
181 auto const& dNdx_v_op = ip_data.dNdx_v_op;
182 auto const& dNdx_p = ip_data.dNdx_p;
184 auto const& w = ip_data.integration_weight;
186 double p_int_pt = 0.0;
190 vars[
static_cast<int>(
195 .template value<double>(vars, pos, t, dt);
201 auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
203 .value(vars, pos, t, dt));
206 w * N_v_op.transpose() * mu * K.inverse() * N_v_op;
209 for (
int i = 0; i < GlobalDim; ++i)
215 .noalias() += w * dNdx_v.transpose() * mu * dNdx_v;
219 Kvp.noalias() += w * N_v_op.transpose() * dNdx_p;
222 Kpv.noalias() += w * N_p.transpose() * unit_vector * dNdx_v_op;
225 bv.noalias() += w * N_v_op.transpose() * b;
230 const unsigned integration_point)
const override
232 auto const& N_p =
_ip_data[integration_point].N_p;
235 return Eigen::Map<const Eigen::RowVectorXd>(N_p.data(), N_p.size());
241 Eigen::VectorXd
const& local_x,
242 Eigen::VectorXd
const& )
override
248 ShapeFunctionPressure,
249 typename ShapeFunctionLiquidVelocity::MeshElement, GlobalDim>(
262 ShapeFunctionLiquidVelocity::NPOINTS>,
265 GlobalDim, ShapeFunctionLiquidVelocity::NPOINTS>>>
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< IntegrationPointData< ShapeMatrixTypeLiquidVelocity, ShapeMatrixTypePressure, GlobalDim, ShapeFunctionLiquidVelocity::NPOINTS >, Eigen::aligned_allocator< IntegrationPointData< ShapeMatrixTypeLiquidVelocity, ShapeMatrixTypePressure, GlobalDim, ShapeFunctionLiquidVelocity::NPOINTS > > > _ip_data
void assemble(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &, std::vector< double > &, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
void computeSecondaryVariableConcrete(double const, double const, Eigen::VectorXd const &local_x, Eigen::VectorXd const &) override
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
IntegrationMethod const _integration_method
StokesFlowProcessData const & _process_data
MeshLib::Element const & _element
static const int liquid_velocity_size
bool const _is_axially_symmetric
static const int liquid_velocity_index
LocalAssemblerData(MeshLib::Element const &element, std::size_t const, bool const is_axially_symmetric, unsigned const integration_order, StokesFlowProcessData const &process_data)
static const int pressure_size
ShapeMatrixPolicyType< ShapeFunctionLiquidVelocity, GlobalDim > ShapeMatrixTypeLiquidVelocity
static const int pressure_index
typename ShapeMatrixTypePressure::NodalVectorType NodalVectorType
ShapeMatrixPolicyType< ShapeFunctionPressure, GlobalDim > ShapeMatrixTypePressure
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)
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
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)
VectorType< ShapeFunction::NPOINTS > NodalVectorType
Eigen::VectorXd const specific_body_force
an external force that applies in the bulk of the fluid, like gravity.
std::unique_ptr< MaterialPropertyLib::MaterialSpatialDistributionMap > media_map
MeshLib::PropertyVector< double > * pressure_interpolated
bool const use_stokes_brinkman_form