1 #ifndef RTN_PK_PK_PK_HPP
2 #define RTN_PK_PK_PK_HPP
4 #include <boost/fusion/include/map.hpp>
14 class RTNPkPkPk:
public HYPRE {
16 typedef boost::fusion::map<boost::fusion::pair<RolyCellType, std::string>,
36 return m_cell_dofs_map;
50 Eigen::VectorXd _interpolate_velocity(
55 void _construct_local_operators(
size_t iT);
56 void _construct_stabilization(
size_t iT);
57 std::pair<Eigen::MatrixXd, Eigen::VectorXd> _velocity_pressure_coupling(
60 const Eigen::VectorXd & uT_old,
61 const Eigen::VectorXd & pT_old
63 Eigen::MatrixXd _forcing_terms(
71 std::unique_ptr<DiscreteSpace> m_pkpo_pk;
77 Eigen::VectorXd RTNPkPkPk::_interpolate_velocity(
79 const InterpolateParameters & parameters
85 int doe_cell = (parameters.doe_cell > 0 ? doe_cell : 2 *
m_degree + 1);
86 int doe_edge = (parameters.doe_edge > 0 ? doe_edge : 2 *
m_degree + 3);
90 = [
this,
v, &vh, &doe_edge](
size_t start,
size_t end)->
void
94 for (
size_t iE = start; iE <
end; iE++) {
96 const auto & edge_dofs_E = edge_dofs[E.global_index()];
99 vh.segment(this->
m_velocity_space->globalOffset(E), this->m_velocity_space->numberOfLocalEdgeDofs())
107 = [
this, &vh,
v, &doe_cell, &doe_edge](
size_t start,
size_t end)->
void
118 size_t dim_cell_dofs = dim_roly_cell_dofs + dim_roly_compl_cell_dofs;
120 for (
size_t iT = start; iT <
end; iT++) {
124 Eigen::MatrixXd MIT = Eigen::MatrixXd::Zero(dim_cell_dofs, dim_cell_dofs);
125 Eigen::VectorXd bIT = Eigen::VectorXd(dim_cell_dofs);
130 if (dim_rtn_cell_dofs > 0) {
133 MIT.block(m_rtn_dof_table.
localOffset(
T), 0, dim_rtn_cell_dofs, dim_roly_cell_dofs)
134 =
GramMatrix(
T, *rtn_cell_dofs, *roly_cell_dofs[iT], int_mono_cell);
135 MIT.block(m_rtn_dof_table.
localOffset(
T), dim_roly_cell_dofs, dim_rtn_cell_dofs, dim_roly_compl_cell_dofs)
136 =
GramMatrix(
T, *rtn_cell_dofs, *roly_compl_cell_dofs[iT], int_mono_cell);
139 bIT.segment(m_rtn_dof_table.
localOffset(
T), dim_rtn_cell_dofs)
143 for (
size_t iE = 0; iE <
T.n_edges(); iE++) {
144 const Edge & E = *
T.edge(iE);
145 auto nTE =
T.edge_normal(iE);
147 const auto & rtn_edge_dofs_E = rtn_edge_dofs[E.global_index()];
149 auto normal_component = [&nTE](
const VectorRd & w)->
double {
return w.dot(nTE); };
152 auto roly_cell_dofs_quad
155 auto roly_compl_cell_dofs_quad
162 = [
v, &nTE](
const VectorRd & x) {
return v(x).dot(nTE); };
164 MIT.block(m_rtn_dof_table.
localOffset(
T, E), 0, dim_rtn_edge_dofs, dim_roly_cell_dofs)
166 MIT.block(m_rtn_dof_table.
localOffset(
T, E), dim_roly_cell_dofs, dim_rtn_edge_dofs, dim_roly_compl_cell_dofs)
169 bIT.segment(m_rtn_dof_table.
localOffset(
T, E), dim_rtn_edge_dofs)
170 =
l2_projection(evaluate_normal_component, *rtn_edge_dofs_E, quad_E, rtn_edge_dofs_quad,
GramMatrix(E, *rtn_edge_dofs_E));
173 vh.segment(this->
m_velocity_space->globalOffset(T), this->m_velocity_space->numberOfLocalCellDofs())
174 = MIT.fullPivLu().solve(bIT);
The BoundaryConditions class provides definition of boundary conditions.
Definition: BoundaryConditions.hpp:45
Definition: discrete-space.hpp:226
std::function< VectorRd(const VectorRd &)> VelocityType
Definition: hypre.hpp:35
size_t degree() const
Returns the degree.
Definition: hypre.hpp:62
const Mesh & mesh() const
Returns the mesh.
Definition: hypre.hpp:164
std::unique_ptr< DiscreteSpace > m_velocity_space
Definition: hypre.hpp:260
size_t m_degree
Definition: hypre.hpp:244
std::function< double(const VectorRd &)> PressureType
Definition: hypre.hpp:37
bool m_use_threads
Definition: hypre.hpp:247
Eigen::VectorXd extendCellDofs(const Eigen::VectorXd &vh) const
Extend cell DOFs from to .
const DiscreteSpace & extendedDofsSpace() const
Return the discrete space corresponding to the extended DOFs.
Definition: rtn-pk-pk-pk.hpp:43
RTNPkPkPk(const Mesh &mesh, size_t degree, const BoundaryConditions &bc, const HYPREParameters ¶meters)
Constructor.
const CellDofsMapType & cellDofsMap() const
Return cell DOFs map.
Definition: rtn-pk-pk-pk.hpp:34
boost::fusion::map< boost::fusion::pair< RolyCellType, std::string >, boost::fusion::pair< RolyComplCellType, std::string > > CellDofsMapType
Definition: rtn-pk-pk-pk.hpp:17
Eigen::VectorXd interpolate(const VelocityType &u, const PressureType &p) const
Interpolate both velocity and pressure.
Base class for global DOF spaces. Provides functions to manipulate global DOFs (the local version bei...
Definition: globaldofspace.hpp:16
Definition: Mesh2D.hpp:26
function[maxeig, mineig, cond, mat]
Definition: compute_eigs.m:1
end
Definition: compute_eigs.m:12
Eigen::Vector2d VectorRd
Definition: basis.hpp:55
Eigen::MatrixXd compute_gram_matrix(const boost::multi_array< VectorRd, 2 > &B1, const boost::multi_array< double, 2 > &B2, const QuadratureRule &qr)
Compute the Gram-like matrix given a family of vector-valued and one of scalar-valued functions by te...
Definition: basis.cpp:225
static boost::multi_array< typename detail::basis_evaluation_traits< BasisType, BasisFunction >::ReturnValue, 2 > compute(const BasisType &basis, const QuadratureRule &quad)
Generic basis evaluation.
Definition: basis.hpp:2288
Eigen::VectorXd l2_projection(const std::function< typename BasisType::FunctionValue(const VectorRd &)> &f, const BasisType &basis, QuadratureRule &quad, const boost::multi_array< typename BasisType::FunctionValue, 2 > &basis_quad, const Eigen::MatrixXd &mass_basis=Eigen::MatrixXd::Zero(1, 1))
Compute the L2-projection of a function.
Definition: basis.hpp:2886
size_t localOffset(const Edge &E, const Vertex &V) const
Returns the local offset of the vertex V with respect to the edge E.
Definition: localdofspace.hpp:143
static void parallel_for(unsigned nb_elements, std::function< void(size_t start, size_t end)> functor, bool use_threads=true, unsigned nb_threads_max=1e9)
Generic function to execute threaded processes.
Definition: parallel_for.hpp:42
Polytope< DIMENSION > Cell
Definition: Polytope2D.hpp:151
Polytope< 1 > Edge
A Face is a Polytope with object_dim = DIMENSION - 1.
Definition: Polytope2D.hpp:147
std::vector< QuadratureNode > QuadratureRule
Definition: quadraturerule.hpp:55
std::unordered_map< VectorZd, double, VecHash > MonomialCellIntegralsType
Type for list of integrals of monomials.
Definition: GMpoly_cell.hpp:53
MonomialCellIntegralsType IntegrateCellMonomials(const Cell &T, const size_t maxdeg)
Compute all integrals on a cell of monomials up to a total degree, using vertex values.
Definition: GMpoly_cell.cpp:7
Eigen::MatrixXd GramMatrix(const Cell &T, const MonomialScalarBasisCell &basis1, const MonomialScalarBasisCell &basis2, MonomialCellIntegralsType mono_int_map={})
Computes the Gram Matrix of a pair of local scalar monomial bases.
Definition: GMpoly_cell.cpp:86
QuadratureRule generate_quadrature_rule(const Cell &T, const int doe, const bool force_split)
Generate quadrature rule on mesh element.
Definition: quadraturerule.cpp:10
if(strcmp(field, 'real')) % real valued entries T
Definition: mmread.m:93
const std::vector< std::unique_ptr< T > > & get(const DiscreteSpace &discrete_space, const std::string &name)
TensorizedVectorFamily< PolyCellType, dimspace > PolynCellType
Definition: discrete-space.hpp:22
TensorizedVectorFamily< Family< MonomialScalarBasisEdge >, dimspace > PolynEdgeType
Definition: discrete-space.hpp:23
Family< RolyComplBasisCell > RolyComplCellType
Definition: discrete-space.hpp:29
Family< CurlBasis< ShiftedBasis< MonomialScalarBasisCell > > > RolyCellType
Definition: discrete-space.hpp:28
Family< MonomialScalarBasisEdge > PolyEdgeType
Definition: discrete-space.hpp:21
Definition: ddr-klplate.hpp:27
static auto v
Definition: ddrcore-test.hpp:32
Definition: hho-interpolate.hpp:15
Definition: ns-solutions.hpp:15
Evaluate a basis at quadrature nodes. 'BasisFunction' (=Function, Gradient, Curl, Divergence,...
Definition: basis.hpp:2284