--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_polynomials_pyramid_h
+#define dealii_polynomials_pyramid_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/scalar_polynomials_base.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * A namespace for functions and classes that provide support for simplex
+ * reference cell entities, i.e., triangles and tetrahedrons.
+ *
+ * @ingroup simplex
+ */
+namespace Simplex
+{
+ /**
+ * Polynomials defined on pyramid entities. This class is basis of
+ * Simplex::FE_PyramidP.
+ */
+ template <int dim>
+ class ScalarPyramidPolynomial : public ScalarPolynomialsBase<dim>
+ {
+ public:
+ /**
+ * Make the dimension available to the outside.
+ */
+ static const unsigned int dimension = dim;
+
+ /*
+ * Constructor taking the polynomial @p degree as input.
+ *
+ * @note Currently, only linear polynomials (degree=1) are implemented.
+ */
+ ScalarPyramidPolynomial(const unsigned int degree);
+
+ /**
+ * @copydoc ScalarPolynomialsBase::evaluate()
+ *
+ * @note Currently, only the vectors @p values and @p grads are filled.
+ */
+ void
+ evaluate(const Point<dim> & unit_point,
+ std::vector<double> & values,
+ std::vector<Tensor<1, dim>> &grads,
+ std::vector<Tensor<2, dim>> &grad_grads,
+ std::vector<Tensor<3, dim>> &third_derivatives,
+ std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
+
+ double
+ compute_value(const unsigned int i, const Point<dim> &p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_derivative()
+ *
+ * @note Currently, only implemented for first derivative.
+ */
+ template <int order>
+ Tensor<order, dim>
+ compute_derivative(const unsigned int i, const Point<dim> &p) const;
+
+ Tensor<1, dim>
+ compute_1st_derivative(const unsigned int i,
+ const Point<dim> & p) const override;
+
+ Tensor<2, dim>
+ compute_2nd_derivative(const unsigned int i,
+ const Point<dim> & p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_3rd_derivative()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<3, dim>
+ compute_3rd_derivative(const unsigned int i,
+ const Point<dim> & p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_4th_derivative()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<4, dim>
+ compute_4th_derivative(const unsigned int i,
+ const Point<dim> & p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_grad()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<1, dim>
+ compute_grad(const unsigned int i, const Point<dim> &p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_grad_grad()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<2, dim>
+ compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
+
+ std::string
+ name() const override;
+
+ virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
+ clone() const override;
+ };
+
+
+
+ template <int dim>
+ template <int order>
+ Tensor<order, dim>
+ ScalarPyramidPolynomial<dim>::compute_derivative(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ Tensor<order, dim> der;
+
+ Assert(order == 1, ExcNotImplemented());
+ const auto grad = compute_grad(i, p);
+
+ for (unsigned int i = 0; i < dim; i++)
+ der[i] = grad[i];
+
+ return der;
+ }
+
+} // namespace Simplex
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
// ---------------------------------------------------------------------
//
-// Copyright (C) 2020 by the deal.II authors
+// Copyright (C) 2021 by the deal.II authors
//
// This file is part of the deal.II library.
//
// ---------------------------------------------------------------------
-#ifndef dealii_simplex_polynomials_h
-#define dealii_simplex_polynomials_h
+#ifndef dealii_base_polynomials_wedge_h
+#define dealii_base_polynomials_wedge_h
#include <deal.II/base/config.h>
+#include <deal.II/base/polynomials_barycentric.h>
#include <deal.II/base/scalar_polynomials_base.h>
-#include <deal.II/simplex/barycentric_polynomials.h>
-
DEAL_II_NAMESPACE_OPEN
/**
std::vector<Tensor<3, dim>> &third_derivatives,
std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
- /**
- * @copydoc ScalarPolynomialsBase::compute_value()
- */
double
compute_value(const unsigned int i, const Point<dim> &p) const override;
Tensor<order, dim>
compute_derivative(const unsigned int i, const Point<dim> &p) const;
- /**
- * @copydoc ScalarPolynomialsBase::compute_1st_derivative()
- */
Tensor<1, dim>
compute_1st_derivative(const unsigned int i,
const Point<dim> & p) const override;
Tensor<2, dim>
compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
- /**
- * @copydoc ScalarPolynomialsBase::name()
- */
std::string
name() const override;
- /**
- * @copydoc ScalarPolynomialsBase::clone()
- */
virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
clone() const override;
- /**
- * Polynomials defined on pyramid entities. This class is basis of
- * Simplex::FE_PyramidP.
- */
- template <int dim>
- class ScalarPyramidPolynomial : public ScalarPolynomialsBase<dim>
- {
- public:
- /**
- * Make the dimension available to the outside.
- */
- static const unsigned int dimension = dim;
-
- /*
- * Constructor taking the polynomial @p degree as input.
- *
- * @note Currently, only linear polynomials (degree=1) are implemented.
- */
- ScalarPyramidPolynomial(const unsigned int degree);
-
- /**
- * @copydoc ScalarPolynomialsBase::evaluate()
- *
- * @note Currently, only the vectors @p values and @p grads are filled.
- */
- void
- evaluate(const Point<dim> & unit_point,
- std::vector<double> & values,
- std::vector<Tensor<1, dim>> &grads,
- std::vector<Tensor<2, dim>> &grad_grads,
- std::vector<Tensor<3, dim>> &third_derivatives,
- std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
-
- /**
- * @copydoc ScalarPolynomialsBase::compute_value()
- */
- double
- compute_value(const unsigned int i, const Point<dim> &p) const override;
-
- /**
- * @copydoc ScalarPolynomialsBase::compute_derivative()
- *
- * @note Currently, only implemented for first derivative.
- */
- template <int order>
- Tensor<order, dim>
- compute_derivative(const unsigned int i, const Point<dim> &p) const;
-
- /**
- * @copydoc ScalarPolynomialsBase::compute_1st_derivative()
- */
- Tensor<1, dim>
- compute_1st_derivative(const unsigned int i,
- const Point<dim> & p) const override;
-
- /**
- * @copydoc ScalarPolynomialsBase::compute_2nd_derivative()
- */
- Tensor<2, dim>
- compute_2nd_derivative(const unsigned int i,
- const Point<dim> & p) const override;
-
- /**
- * @copydoc ScalarPolynomialsBase::compute_3rd_derivative()
- *
- * @note Not implemented yet.
- */
- Tensor<3, dim>
- compute_3rd_derivative(const unsigned int i,
- const Point<dim> & p) const override;
-
- /**
- * @copydoc ScalarPolynomialsBase::compute_4th_derivative()
- *
- * @note Not implemented yet.
- */
- Tensor<4, dim>
- compute_4th_derivative(const unsigned int i,
- const Point<dim> & p) const override;
-
- /**
- * @copydoc ScalarPolynomialsBase::compute_grad()
- *
- * @note Not implemented yet.
- */
- Tensor<1, dim>
- compute_grad(const unsigned int i, const Point<dim> &p) const override;
-
- /**
- * @copydoc ScalarPolynomialsBase::compute_grad_grad()
- *
- * @note Not implemented yet.
- */
- Tensor<2, dim>
- compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
-
- /**
- * @copydoc ScalarPolynomialsBase::name()
- */
- std::string
- name() const override;
-
- /**
- * @copydoc ScalarPolynomialsBase::clone()
- */
- virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
- clone() const override;
- };
-
-
-
template <int dim>
template <int order>
Tensor<order, dim>
return der;
}
-
-
- template <int dim>
- template <int order>
- Tensor<order, dim>
- ScalarPyramidPolynomial<dim>::compute_derivative(const unsigned int i,
- const Point<dim> & p) const
- {
- Tensor<order, dim> der;
-
- AssertDimension(order, 1);
- const auto grad = compute_grad(i, p);
-
- for (unsigned int i = 0; i < dim; i++)
- der[i] = grad[i];
-
- return der;
- }
-
} // namespace Simplex
DEAL_II_NAMESPACE_CLOSE
QSplit(const QSimplex<dim> &base, const Point<dim> &split_point);
};
+namespace Simplex
+{
+ /**
+ * Integration rule for simplex entities.
+ *
+ * Users specify a number `n_points_1D` as an indication of what polynomial
+ * degree to be integrated exactly, similarly to the number of points in a
+ * QGauss quadrature object, even though the present quadrature formula is not
+ * a tensor product. The given value is translated for n_points_1D=1,2,3,4 to
+ * following number of quadrature points for 2D and 3D:
+ * - 2D: 1, 3, 7, 15
+ * - 3D: 1, 4, 10, 35
+ *
+ * For 1D, the quadrature rule degenerates to a
+ * `dealii::QGauss<1>(n_points_1D)`.
+ *
+ * @ingroup simplex
+ */
+ template <int dim>
+ class QGauss : public QSimplex<dim>
+ {
+ public:
+ /**
+ * Constructor taking the number of quadrature points in 1D direction
+ * @p n_points_1D.
+ */
+ explicit QGauss(const unsigned int n_points_1D);
+ };
+
+ /**
+ * Witherden-Vincent rules for simplex entities.
+ *
+ * Like QGauss, users should specify a number `n_points_1D` as an indication
+ * of what polynomial degree to be integrated exactly (e.g., for $n$ points,
+ * the rule can integrate polynomials of degree $2 n - 1$ exactly). The given
+ * value for n_points_1D = 1, 2, 3, 4, 5 results in the following number of
+ * quadrature points in 2D and 3D:
+ * - 2D: 1, 6, 7, 15, 19
+ * - 3D: 1, 8, 14, 35, 59
+ *
+ * For 1D, the quadrature rule degenerates to a
+ * `dealii::QGauss<1>(n_points_1D)`.
+ *
+ * These rules match the ones listed for Witherden-Vincent in the quadpy
+ * @cite quadpy library and were first described in
+ * @cite witherden2015identification.
+ *
+ * @ingroup simplex
+ */
+ template <int dim>
+ class QWitherdenVincent : public QSimplex<dim>
+ {
+ public:
+ /**
+ * Constructor taking the number of quadrature points in 1D direction
+ * @p n_points_1D.
+ */
+ explicit QWitherdenVincent(const unsigned int n_points_1D);
+ };
+
+ /**
+ * Integration rule for wedge entities.
+ */
+ template <int dim>
+ class QGaussWedge : public Quadrature<dim>
+ {
+ public:
+ /**
+ * Users specify a number `n_points_1D` as an indication of what polynomial
+ * degree to be integrated exactly. For details, see the comments of
+ * Simplex::QGauss.
+ */
+ explicit QGaussWedge(const unsigned int n_points_1D);
+ };
+
+ /**
+ * Integration rule for pyramid entities.
+ */
+ template <int dim>
+ class QGaussPyramid : public Quadrature<dim>
+ {
+ public:
+ /**
+ * Users specify a number `n_points_1D` as an indication of what polynomial
+ * degree to be integrated exactly. For details, see the comments of
+ * Simplex::QGauss.
+ */
+ explicit QGaussPyramid(const unsigned int n_points_1D);
+ };
+
+} // namespace Simplex
+
/*@}*/
/* -------------- declaration of explicit specializations ------------- */
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_fe_fe_p_pyramids_h
+#define dealii_fe_fe_p_pyramids_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/polynomials_pyramid.h>
+
+#include <deal.II/fe/fe_poly.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Simplex
+{
+ /**
+ * Base class of FE_PyramidP and FE_PyramidDGP.
+ *
+ * @note Only implemented for 3D.
+ *
+ * @ingroup simplex
+ */
+ template <int dim, int spacedim = dim>
+ class FE_Pyramid : public dealii::FE_Poly<dim, spacedim>
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ FE_Pyramid(const unsigned int degree,
+ const internal::GenericDoFsPerObject & dpos,
+ const typename FiniteElementData<dim>::Conformity conformity);
+ };
+
+ /**
+ * Implementation of a scalar Lagrange finite element on a pyramid that yields
+ * the finite element space of continuous, piecewise polynomials of
+ * degree $k$.
+ *
+ * @ingroup simplex
+ */
+ template <int dim, int spacedim = dim>
+ class FE_PyramidP : public FE_Pyramid<dim, spacedim>
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ FE_PyramidP(const unsigned int degree);
+
+ /**
+ * @copydoc dealii::FiniteElement::clone()
+ */
+ std::unique_ptr<FiniteElement<dim, spacedim>>
+ clone() const override;
+
+ /**
+ * Return a string that uniquely identifies a finite element. This class
+ * returns <tt>Simplex::FE_PyramidP<dim>(degree)</tt>, with @p dim and @p degree
+ * replaced by appropriate values.
+ */
+ std::string
+ get_name() const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::compare_for_domination()
+ */
+ FiniteElementDomination::Domination
+ compare_for_domination(const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int codim) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_vertex_dof_identities()
+ */
+ std::vector<std::pair<unsigned int, unsigned int>>
+ hp_vertex_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_line_dof_identities()
+ */
+ std::vector<std::pair<unsigned int, unsigned int>>
+ hp_line_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_quad_dof_identities()
+ */
+ std::vector<std::pair<unsigned int, unsigned int>>
+ hp_quad_dof_identities(const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int face_no = 0) const override;
+ };
+
+ /**
+ * Implementation of a scalar Lagrange finite element on a pyramid that yields
+ * the finite element space of discontinuous, piecewise polynomials of
+ * degree $k$.
+ *
+ * @ingroup simplex
+ */
+ template <int dim, int spacedim = dim>
+ class FE_PyramidDGP : public FE_Pyramid<dim, spacedim>
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ FE_PyramidDGP(const unsigned int degree);
+
+ /**
+ * @copydoc dealii::FiniteElement::clone()
+ */
+ std::unique_ptr<FiniteElement<dim, spacedim>>
+ clone() const override;
+
+ /**
+ * Return a string that uniquely identifies a finite element. This class
+ * returns <tt>Simplex::FE_PyramidDGP<dim>(degree)</tt>, with @p dim and @p degree
+ * replaced by appropriate values.
+ */
+ std::string
+ get_name() const override;
+ };
+} // namespace Simplex
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_fe_fe_p_h
+#define dealii_fe_fe_p_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/polynomials_barycentric.h>
+
+#include <deal.II/fe/fe_poly.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Simplex
+{
+ /**
+ * Base class of FE_P and FE_DGP.
+ *
+ * @note Only implemented for 2D and 3D.
+ *
+ * @ingroup simplex
+ */
+ template <int dim, int spacedim = dim>
+ class FE_Poly : public dealii::FE_Poly<dim, spacedim>
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ FE_Poly(const unsigned int degree,
+ const std::vector<unsigned int> & dpo_vector,
+ const typename FiniteElementData<dim>::Conformity conformity);
+
+ /**
+ * Return a list of constant modes of the element. For this element, the
+ * list consists of true arguments for all components.
+ */
+ virtual std::pair<Table<2, bool>, std::vector<unsigned int>>
+ get_constant_modes() const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::get_prolongation_matrix()
+ *
+ * @note Only implemented for RefinementCase::isotropic_refinement.
+ */
+ virtual const FullMatrix<double> &
+ get_prolongation_matrix(
+ const unsigned int child,
+ const RefinementCase<dim> &refinement_case =
+ RefinementCase<dim>::isotropic_refinement) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::get_face_interpolation_matrix()
+ */
+ void
+ get_face_interpolation_matrix(const FiniteElement<dim, spacedim> &source_fe,
+ FullMatrix<double> &interpolation_matrix,
+ const unsigned int face_no) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::get_subface_interpolation_matrix()
+ */
+ void
+ get_subface_interpolation_matrix(
+ const FiniteElement<dim, spacedim> &x_source_fe,
+ const unsigned int subface,
+ FullMatrix<double> & interpolation_matrix,
+ const unsigned int face_no) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_constraints_are_implemented()
+ */
+ bool
+ hp_constraints_are_implemented() const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::convert_generalized_support_point_values_to_dof_values()
+ */
+ virtual void
+ convert_generalized_support_point_values_to_dof_values(
+ const std::vector<Vector<double>> &support_point_values,
+ std::vector<double> & nodal_values) const override;
+
+ mutable Threads::Mutex mutex;
+ };
+
+
+
+ /**
+ * Implementation of a scalar Lagrange finite element $P_k$ that yields
+ * the finite element space of continuous, piecewise polynomials of
+ * degree $k$.
+ *
+ * @ingroup simplex
+ */
+ template <int dim, int spacedim = dim>
+ class FE_P : public FE_Poly<dim, spacedim>
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ FE_P(const unsigned int degree);
+
+ /**
+ * @copydoc dealii::FiniteElement::clone()
+ */
+ std::unique_ptr<FiniteElement<dim, spacedim>>
+ clone() const override;
+
+ /**
+ * Return a string that uniquely identifies a finite element. This class
+ * returns <tt>Simplex::FE_P<dim>(degree)</tt>, with @p dim and @p degree
+ * replaced by appropriate values.
+ */
+ std::string
+ get_name() const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::compare_for_domination()
+ */
+ FiniteElementDomination::Domination
+ compare_for_domination(const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int codim) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_vertex_dof_identities()
+ */
+ std::vector<std::pair<unsigned int, unsigned int>>
+ hp_vertex_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_line_dof_identities()
+ */
+ std::vector<std::pair<unsigned int, unsigned int>>
+ hp_line_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const override;
+ };
+
+
+
+ /**
+ * Implementation of a scalar discontinuous Lagrange finite element
+ * $P_k$, sometimes denoted as $P_{-k}$, that yields the finite
+ * element space of discontinuous, piecewise polynomials of degree
+ * $k$.
+ *
+ * @ingroup simplex
+ */
+ template <int dim, int spacedim = dim>
+ class FE_DGP : public FE_Poly<dim, spacedim>
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ FE_DGP(const unsigned int degree);
+
+ /**
+ * @copydoc dealii::FiniteElement::clone()
+ */
+ std::unique_ptr<FiniteElement<dim, spacedim>>
+ clone() const override;
+
+ /**
+ * Return a string that uniquely identifies a finite element. This class
+ * returns <tt>Simplex::FE_DGP<dim>(degree)</tt>, with @p dim and @p degree
+ * replaced by appropriate values.
+ */
+ std::string
+ get_name() const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::compare_for_domination()
+ */
+ FiniteElementDomination::Domination
+ compare_for_domination(const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int codim) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_vertex_dof_identities()
+ */
+ std::vector<std::pair<unsigned int, unsigned int>>
+ hp_vertex_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_line_dof_identities()
+ */
+ std::vector<std::pair<unsigned int, unsigned int>>
+ hp_line_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const override;
+ };
+} // namespace Simplex
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_fe_fe_p_bubbles_h
+#define dealii_fe_fe_p_bubbles_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/polynomials_barycentric.h>
+
+#include <deal.II/fe/fe_poly.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Simplex
+{
+ /**
+ * @brief Enriched version of FE_P that can be used with nodal quadrature.
+ *
+ * Many explicit time integration schemes require solving a mass matrix at
+ * each time step. There are various ways around this requirement - for
+ * example, step-48 replaces the mass matrix with a diagonal approximation,
+ * which makes the solution step trivial. In step-48, and also commonly for
+ * tensor-product elements, this is done by computing the mass matrix with a
+ * lower-order quadrature point based on the nodes of the finite element
+ * (i.e., the nodal quadrature rule one obtains by using the shape functions
+ * as an interpolatory basis).
+ *
+ * A major drawback of standard simplex-based finite elements is that they
+ * cannot be used with nodal quadrature since some of the quadrature weights
+ * end up being either zero or negative, resulting in either an unsolvable or
+ * unstable approximation to the mass matrix. For example: the shape functions
+ * of FE_P<2>(2) with support points at vertices have mean values of zero so
+ * that element cannot be used with mass lumping.
+
+ * This element avoids this issue by replacing the shape functions of FE_P
+ * with an augmented space amendable to the construction of nodal quadrature
+ * rules. For example, on the triangle a single basis function is added
+ * corresponding to interpolation at the centroid (and all other basis
+ * functions are updated to preserve the partition of unity property). This
+ * results in shape functions with positive means (i.e., a valid nodal
+ * quadrature formula). Similarly, in 3D, the polynomial space of FE_P<3>(2)
+ * is enriched with five additional degrees of freedom (where four have
+ * support points at face centroids and one has a support point at the
+ * centroid) to enable construction of valid nodal quadrature rule.
+ *
+ * Since this FE space includes bubbles (i.e., extra functions which are
+ * nonzero only on element interiors), the polynomial degrees of the component
+ * basis functions are higher than the actual approximation degree of the
+ * element. For example, with a constructor argument <code>degree = 2</code>
+ * in 3D, the polynomials are in fact cubic (degree 3) but the order of the
+ * approximation is the same as if we were using quadratic (degree 2) finite
+ * elements.
+ *
+ * The 2D quadratic element was first described in @cite fried1975finite. The
+ * 3D quadratic element implemented here was first described in
+ * @cite Geevers_2018. Higher degree elements amendable to lumping exist but
+ * are not yet implemented in this class.
+ */
+ template <int dim, int spacedim = dim>
+ class FE_P_Bubbles : public dealii::FE_Poly<dim, spacedim>
+ {
+ public:
+ /**
+ * Constructor, taking the approximation degree as an argument. The
+ * polynomial space is typically one degree higher than the approximation
+ * space for this element: see the general documentation of this class for
+ * more information.
+ *
+ * @note For <code>degree == 1</code> this element is equivalent to FE_P(1).
+ */
+ FE_P_Bubbles(const unsigned int degree);
+
+ /**
+ * @copydoc dealii::FiniteElement::clone()
+ */
+ virtual std::unique_ptr<FiniteElement<dim, spacedim>>
+ clone() const override;
+
+ /**
+ * Return a string that uniquely identifies a finite element. This class
+ * returns <tt>Simplex::FE_F_Bubbles<dim,spacedim>(degree)</tt>, with
+ * @p dim, @p spacedim, and @p degree replaced by appropriate values. As
+ * usual, @p spacedim is omitted in the codimension zero case.
+ */
+ virtual std::string
+ get_name() const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::convert_generalized_support_point_values_to_dof_values()
+ */
+ virtual void
+ convert_generalized_support_point_values_to_dof_values(
+ const std::vector<Vector<double>> &support_point_values,
+ std::vector<double> & nodal_values) const override;
+
+ protected:
+ /**
+ * Degree of the approximation (i.e., the constructor argument).
+ */
+ unsigned int approximation_degree;
+ };
+} // namespace Simplex
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_fe_fe_p_wedge_h
+#define dealii_fe_fe_p_wedge_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/polynomials_wedge.h>
+
+#include <deal.II/fe/fe_poly.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Simplex
+{
+ /**
+ * Base class of FE_WedgeP and FE_WedgeDGP.
+ *
+ * @note Only implemented for 3D.
+ *
+ * @ingroup simplex
+ */
+ template <int dim, int spacedim = dim>
+ class FE_Wedge : public dealii::FE_Poly<dim, spacedim>
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ FE_Wedge(const unsigned int degree,
+ const internal::GenericDoFsPerObject & dpos,
+ const typename FiniteElementData<dim>::Conformity conformity);
+ };
+
+ /**
+ * Implementation of a scalar Lagrange finite element on a wedge that yields
+ * the finite element space of continuous, piecewise polynomials of
+ * degree $k$.
+ *
+ * @ingroup simplex
+ */
+ template <int dim, int spacedim = dim>
+ class FE_WedgeP : public FE_Wedge<dim, spacedim>
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ FE_WedgeP(const unsigned int degree);
+
+ /**
+ * @copydoc dealii::FiniteElement::clone()
+ */
+ std::unique_ptr<FiniteElement<dim, spacedim>>
+ clone() const override;
+
+ /**
+ * Return a string that uniquely identifies a finite element. This class
+ * returns <tt>Simplex::FE_WedgeP<dim>(degree)</tt>, with @p dim and @p degree
+ * replaced by appropriate values.
+ */
+ std::string
+ get_name() const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::compare_for_domination()
+ */
+ FiniteElementDomination::Domination
+ compare_for_domination(const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int codim) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_vertex_dof_identities()
+ */
+ std::vector<std::pair<unsigned int, unsigned int>>
+ hp_vertex_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_line_dof_identities()
+ */
+ std::vector<std::pair<unsigned int, unsigned int>>
+ hp_line_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const override;
+
+ /**
+ * @copydoc dealii::FiniteElement::hp_quad_dof_identities()
+ */
+ std::vector<std::pair<unsigned int, unsigned int>>
+ hp_quad_dof_identities(const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int face_no = 0) const override;
+ };
+
+ /**
+ * Implementation of a scalar Lagrange finite element on a wedge that yields
+ * the finite element space of discontinuous, piecewise polynomials of
+ * degree $k$.
+ *
+ * @ingroup simplex
+ */
+ template <int dim, int spacedim = dim>
+ class FE_WedgeDGP : public FE_Wedge<dim, spacedim>
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ FE_WedgeDGP(const unsigned int degree);
+
+ /**
+ * @copydoc dealii::FiniteElement::clone()
+ */
+ std::unique_ptr<FiniteElement<dim, spacedim>>
+ clone() const override;
+
+ /**
+ * Return a string that uniquely identifies a finite element. This class
+ * returns <tt>Simplex::FE_WedgeDGP<dim>(degree)</tt>, with @p dim and @p degree
+ * replaced by appropriate values.
+ */
+ std::string
+ get_name() const override;
+ };
+} // namespace Simplex
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
add_parameters(ParameterHandler &prm);
};
-
-
/**
* Initialize the given triangulation with a flow field around an airfoil,
* i.e., a mesh of C-Type approximating Joukowski or NACA (4 digit)
} // namespace Airfoil
+ /**
+ * Create a coordinate-parallel brick from the two diagonally opposite
+ * corner points @p p1 and @p p2. The number of vertices in coordinate
+ * direction @p i is given by <tt>repetitions[i]+1</tt>.
+ *
+ * @note This function connects internally 4/8 vertices to
+ * quadrilateral/hexahedral cells and subdivides these into 2/5
+ * triangular/tetrahedral cells.
+ *
+ * @note Currently, this function only works for `dim==spacedim`.
+ *
+ * @ingroup simplex
+ */
+ template <int dim, int spacedim>
+ void
+ subdivided_hyper_rectangle_with_simplices(
+ Triangulation<dim, spacedim> & tria,
+ const std::vector<unsigned int> &repetitions,
+ const Point<dim> & p1,
+ const Point<dim> & p2,
+ const bool colorize = false);
+
+ /**
+ * Initialize the given triangulation with a hypercube (square in 2D and
+ * cube in 3D) consisting of @p repetitions cells in each direction.
+ * The hypercube volume is the tensor product interval
+ * $[left,right]^{\text{dim}}$ in the present number of dimensions, where
+ * the limits are given as arguments. They default to zero and unity, then
+ * producing the unit hypercube.
+ *
+ * @note This function connects internally 4/8 vertices to
+ * quadrilateral/hexahedral cells and subdivides these into 2/5
+ * triangular/tetrahedral cells.
+ *
+ * @ingroup simplex
+ */
+ template <int dim, int spacedim>
+ void
+ subdivided_hyper_cube_with_simplices(Triangulation<dim, spacedim> &tria,
+ const unsigned int repetitions,
+ const double p1 = 0.0,
+ const double p2 = 1.0,
+ const bool colorize = false);
+
///@}
/**
#include <deal.II/base/mpi.h>
#include <deal.II/base/mpi_consensus_algorithms.h>
#include <deal.II/base/polynomials_piecewise.h>
+#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/tensor_product_polynomials.h>
#include <deal.II/base/utilities.h>
#include <deal.II/matrix_free/face_setup_internal.h>
#include <deal.II/matrix_free/matrix_free.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#ifdef DEAL_II_WITH_TBB
# include <deal.II/base/parallel.h>
#include <deal.II/fe/fe_dgp.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_poly.h>
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q_dg0.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/grid/reference_cell.h>
#include <deal.II/matrix_free/shape_info.h>
#include <deal.II/matrix_free/util.h>
-#include <deal.II/simplex/fe_lib.h>
-
DEAL_II_NAMESPACE_OPEN
#include <deal.II/base/config.h>
#include <deal.II/base/quadrature.h>
+#include <deal.II/base/quadrature_lib.h>
#include <deal.II/grid/reference_cell.h>
#include <deal.II/hp/q_collection.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
DEAL_II_NAMESPACE_OPEN
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_simplex_p.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping.h>
#include <deal.II/grid/tria.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/data_out_dof_data.h>
-#include <deal.II/simplex/fe_lib.h>
-
#include <memory>
#include <string>
#include <utility>
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2020 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-#ifndef dealii_simplex_fe_lib_h
-#define dealii_simplex_fe_lib_h
-
-#include <deal.II/base/config.h>
-
-#include <deal.II/fe/fe_poly.h>
-
-#include <deal.II/simplex/polynomials.h>
-
-DEAL_II_NAMESPACE_OPEN
-
-namespace Simplex
-{
- /**
- * Base class of FE_P and FE_DGP.
- *
- * @note Only implemented for 2D and 3D.
- *
- * @ingroup simplex
- */
- template <int dim, int spacedim = dim>
- class FE_Poly : public dealii::FE_Poly<dim, spacedim>
- {
- public:
- /**
- * Constructor.
- */
- FE_Poly(const unsigned int degree,
- const std::vector<unsigned int> & dpo_vector,
- const typename FiniteElementData<dim>::Conformity conformity);
-
- /**
- * Return a list of constant modes of the element. For this element, the
- * list consists of true arguments for all components.
- */
- virtual std::pair<Table<2, bool>, std::vector<unsigned int>>
- get_constant_modes() const override;
-
- /**
- * @copydoc dealii::FiniteElement::get_prolongation_matrix()
- *
- * @note Only implemented for RefinementCase::isotropic_refinement.
- */
- virtual const FullMatrix<double> &
- get_prolongation_matrix(
- const unsigned int child,
- const RefinementCase<dim> &refinement_case =
- RefinementCase<dim>::isotropic_refinement) const override;
-
- /**
- * @copydoc dealii::FiniteElement::get_face_interpolation_matrix()
- */
- void
- get_face_interpolation_matrix(const FiniteElement<dim, spacedim> &source_fe,
- FullMatrix<double> &interpolation_matrix,
- const unsigned int face_no) const override;
-
- /**
- * @copydoc dealii::FiniteElement::get_subface_interpolation_matrix()
- */
- void
- get_subface_interpolation_matrix(
- const FiniteElement<dim, spacedim> &x_source_fe,
- const unsigned int subface,
- FullMatrix<double> & interpolation_matrix,
- const unsigned int face_no) const override;
-
- /**
- * @copydoc dealii::FiniteElement::hp_constraints_are_implemented()
- */
- bool
- hp_constraints_are_implemented() const override;
-
- /**
- * @copydoc dealii::FiniteElement::convert_generalized_support_point_values_to_dof_values()
- */
- virtual void
- convert_generalized_support_point_values_to_dof_values(
- const std::vector<Vector<double>> &support_point_values,
- std::vector<double> & nodal_values) const override;
-
- mutable Threads::Mutex mutex;
- };
-
-
-
- /**
- * Implementation of a scalar Lagrange finite element $P_k$ that yields
- * the finite element space of continuous, piecewise polynomials of
- * degree $k$.
- *
- * @ingroup simplex
- */
- template <int dim, int spacedim = dim>
- class FE_P : public FE_Poly<dim, spacedim>
- {
- public:
- /**
- * Constructor.
- */
- FE_P(const unsigned int degree);
-
- /**
- * @copydoc dealii::FiniteElement::clone()
- */
- std::unique_ptr<FiniteElement<dim, spacedim>>
- clone() const override;
-
- /**
- * Return a string that uniquely identifies a finite element. This class
- * returns <tt>Simplex::FE_P<dim>(degree)</tt>, with @p dim and @p degree
- * replaced by appropriate values.
- */
- std::string
- get_name() const override;
-
- /**
- * @copydoc dealii::FiniteElement::compare_for_domination()
- */
- FiniteElementDomination::Domination
- compare_for_domination(const FiniteElement<dim, spacedim> &fe_other,
- const unsigned int codim) const override;
-
- /**
- * @copydoc dealii::FiniteElement::hp_vertex_dof_identities()
- */
- std::vector<std::pair<unsigned int, unsigned int>>
- hp_vertex_dof_identities(
- const FiniteElement<dim, spacedim> &fe_other) const override;
-
- /**
- * @copydoc dealii::FiniteElement::hp_line_dof_identities()
- */
- std::vector<std::pair<unsigned int, unsigned int>>
- hp_line_dof_identities(
- const FiniteElement<dim, spacedim> &fe_other) const override;
- };
-
-
-
- /**
- * Implementation of a scalar discontinuous Lagrange finite element
- * $P_k$, sometimes denoted as $P_{-k}$, that yields the finite
- * element space of discontinuous, piecewise polynomials of degree
- * $k$.
- *
- * @ingroup simplex
- */
- template <int dim, int spacedim = dim>
- class FE_DGP : public FE_Poly<dim, spacedim>
- {
- public:
- /**
- * Constructor.
- */
- FE_DGP(const unsigned int degree);
-
- /**
- * @copydoc dealii::FiniteElement::clone()
- */
- std::unique_ptr<FiniteElement<dim, spacedim>>
- clone() const override;
-
- /**
- * Return a string that uniquely identifies a finite element. This class
- * returns <tt>Simplex::FE_DGP<dim>(degree)</tt>, with @p dim and @p degree
- * replaced by appropriate values.
- */
- std::string
- get_name() const override;
-
- /**
- * @copydoc dealii::FiniteElement::compare_for_domination()
- */
- FiniteElementDomination::Domination
- compare_for_domination(const FiniteElement<dim, spacedim> &fe_other,
- const unsigned int codim) const override;
-
- /**
- * @copydoc dealii::FiniteElement::hp_vertex_dof_identities()
- */
- std::vector<std::pair<unsigned int, unsigned int>>
- hp_vertex_dof_identities(
- const FiniteElement<dim, spacedim> &fe_other) const override;
-
- /**
- * @copydoc dealii::FiniteElement::hp_line_dof_identities()
- */
- std::vector<std::pair<unsigned int, unsigned int>>
- hp_line_dof_identities(
- const FiniteElement<dim, spacedim> &fe_other) const override;
- };
-
- /**
- * Base class of FE_WedgeP and FE_WedgeDGP.
- *
- * @note Only implemented for 3D.
- *
- * @ingroup simplex
- */
- template <int dim, int spacedim = dim>
- class FE_Wedge : public dealii::FE_Poly<dim, spacedim>
- {
- public:
- /**
- * Constructor.
- */
- FE_Wedge(const unsigned int degree,
- const internal::GenericDoFsPerObject & dpos,
- const typename FiniteElementData<dim>::Conformity conformity);
- };
-
- /**
- * Implementation of a scalar Lagrange finite element on a wedge that yields
- * the finite element space of continuous, piecewise polynomials of
- * degree $k$.
- *
- * @ingroup simplex
- */
- template <int dim, int spacedim = dim>
- class FE_WedgeP : public FE_Wedge<dim, spacedim>
- {
- public:
- /**
- * Constructor.
- */
- FE_WedgeP(const unsigned int degree);
-
- /**
- * @copydoc dealii::FiniteElement::clone()
- */
- std::unique_ptr<FiniteElement<dim, spacedim>>
- clone() const override;
-
- /**
- * Return a string that uniquely identifies a finite element. This class
- * returns <tt>Simplex::FE_WedgeP<dim>(degree)</tt>, with @p dim and @p degree
- * replaced by appropriate values.
- */
- std::string
- get_name() const override;
-
- /**
- * @copydoc dealii::FiniteElement::compare_for_domination()
- */
- FiniteElementDomination::Domination
- compare_for_domination(const FiniteElement<dim, spacedim> &fe_other,
- const unsigned int codim) const override;
-
- /**
- * @copydoc dealii::FiniteElement::hp_vertex_dof_identities()
- */
- std::vector<std::pair<unsigned int, unsigned int>>
- hp_vertex_dof_identities(
- const FiniteElement<dim, spacedim> &fe_other) const override;
-
- /**
- * @copydoc dealii::FiniteElement::hp_line_dof_identities()
- */
- std::vector<std::pair<unsigned int, unsigned int>>
- hp_line_dof_identities(
- const FiniteElement<dim, spacedim> &fe_other) const override;
-
- /**
- * @copydoc dealii::FiniteElement::hp_quad_dof_identities()
- */
- std::vector<std::pair<unsigned int, unsigned int>>
- hp_quad_dof_identities(const FiniteElement<dim, spacedim> &fe_other,
- const unsigned int face_no = 0) const override;
- };
-
- /**
- * Implementation of a scalar Lagrange finite element on a wedge that yields
- * the finite element space of discontinuous, piecewise polynomials of
- * degree $k$.
- *
- * @ingroup simplex
- */
- template <int dim, int spacedim = dim>
- class FE_WedgeDGP : public FE_Wedge<dim, spacedim>
- {
- public:
- /**
- * Constructor.
- */
- FE_WedgeDGP(const unsigned int degree);
-
- /**
- * @copydoc dealii::FiniteElement::clone()
- */
- std::unique_ptr<FiniteElement<dim, spacedim>>
- clone() const override;
-
- /**
- * Return a string that uniquely identifies a finite element. This class
- * returns <tt>Simplex::FE_WedgeDGP<dim>(degree)</tt>, with @p dim and @p degree
- * replaced by appropriate values.
- */
- std::string
- get_name() const override;
- };
-
- /**
- * Base class of FE_PyramidP and FE_PyramidDGP.
- *
- * @note Only implemented for 3D.
- *
- * @ingroup simplex
- */
- template <int dim, int spacedim = dim>
- class FE_Pyramid : public dealii::FE_Poly<dim, spacedim>
- {
- public:
- /**
- * Constructor.
- */
- FE_Pyramid(const unsigned int degree,
- const internal::GenericDoFsPerObject & dpos,
- const typename FiniteElementData<dim>::Conformity conformity);
- };
-
- /**
- * Implementation of a scalar Lagrange finite element on a pyramid that yields
- * the finite element space of continuous, piecewise polynomials of
- * degree $k$.
- *
- * @ingroup simplex
- */
- template <int dim, int spacedim = dim>
- class FE_PyramidP : public FE_Pyramid<dim, spacedim>
- {
- public:
- /**
- * Constructor.
- */
- FE_PyramidP(const unsigned int degree);
-
- /**
- * @copydoc dealii::FiniteElement::clone()
- */
- std::unique_ptr<FiniteElement<dim, spacedim>>
- clone() const override;
-
- /**
- * Return a string that uniquely identifies a finite element. This class
- * returns <tt>Simplex::FE_PyramidP<dim>(degree)</tt>, with @p dim and @p degree
- * replaced by appropriate values.
- */
- std::string
- get_name() const override;
-
- /**
- * @copydoc dealii::FiniteElement::compare_for_domination()
- */
- FiniteElementDomination::Domination
- compare_for_domination(const FiniteElement<dim, spacedim> &fe_other,
- const unsigned int codim) const override;
-
- /**
- * @copydoc dealii::FiniteElement::hp_vertex_dof_identities()
- */
- std::vector<std::pair<unsigned int, unsigned int>>
- hp_vertex_dof_identities(
- const FiniteElement<dim, spacedim> &fe_other) const override;
-
- /**
- * @copydoc dealii::FiniteElement::hp_line_dof_identities()
- */
- std::vector<std::pair<unsigned int, unsigned int>>
- hp_line_dof_identities(
- const FiniteElement<dim, spacedim> &fe_other) const override;
-
- /**
- * @copydoc dealii::FiniteElement::hp_quad_dof_identities()
- */
- std::vector<std::pair<unsigned int, unsigned int>>
- hp_quad_dof_identities(const FiniteElement<dim, spacedim> &fe_other,
- const unsigned int face_no = 0) const override;
- };
-
- /**
- * Implementation of a scalar Lagrange finite element on a pyramid that yields
- * the finite element space of discontinuous, piecewise polynomials of
- * degree $k$.
- *
- * @ingroup simplex
- */
- template <int dim, int spacedim = dim>
- class FE_PyramidDGP : public FE_Pyramid<dim, spacedim>
- {
- public:
- /**
- * Constructor.
- */
- FE_PyramidDGP(const unsigned int degree);
-
- /**
- * @copydoc dealii::FiniteElement::clone()
- */
- std::unique_ptr<FiniteElement<dim, spacedim>>
- clone() const override;
-
- /**
- * Return a string that uniquely identifies a finite element. This class
- * returns <tt>Simplex::FE_PyramidDGP<dim>(degree)</tt>, with @p dim and @p degree
- * replaced by appropriate values.
- */
- std::string
- get_name() const override;
- };
-
- /**
- * @brief Enriched version of FE_P that can be used with nodal quadrature.
- *
- * Many explicit time integration schemes require solving a mass matrix at
- * each time step. There are various ways around this requirement - for
- * example, step-48 replaces the mass matrix with a diagonal approximation,
- * which makes the solution step trivial. In step-48, and also commonly for
- * tensor-product elements, this is done by computing the mass matrix with a
- * lower-order quadrature point based on the nodes of the finite element
- * (i.e., the nodal quadrature rule one obtains by using the shape functions
- * as an interpolatory basis).
- *
- * A major drawback of standard simplex-based finite elements is that they
- * cannot be used with nodal quadrature since some of the quadrature weights
- * end up being either zero or negative, resulting in either an unsolvable or
- * unstable approximation to the mass matrix. For example: the shape functions
- * of FE_P<2>(2) with support points at vertices have mean values of zero so
- * that element cannot be used with mass lumping.
-
- * This element avoids this issue by replacing the shape functions of FE_P
- * with an augmented space amendable to the construction of nodal quadrature
- * rules. For example, on the triangle a single basis function is added
- * corresponding to interpolation at the centroid (and all other basis
- * functions are updated to preserve the partition of unity property). This
- * results in shape functions with positive means (i.e., a valid nodal
- * quadrature formula). Similarly, in 3D, the polynomial space of FE_P<3>(2)
- * is enriched with five additional degrees of freedom (where four have
- * support points at face centroids and one has a support point at the
- * centroid) to enable construction of valid nodal quadrature rule.
- *
- * Since this FE space includes bubbles (i.e., extra functions which are
- * nonzero only on element interiors), the polynomial degrees of the component
- * basis functions are higher than the actual approximation degree of the
- * element. For example, with a constructor argument <code>degree = 2</code>
- * in 3D, the polynomials are in fact cubic (degree 3) but the order of the
- * approximation is the same as if we were using quadratic (degree 2) finite
- * elements.
- *
- * The 2D quadratic element was first described in @cite fried1975finite. The
- * 3D quadratic element implemented here was first described in
- * @cite Geevers_2018. Higher degree elements amendable to lumping exist but
- * are not yet implemented in this class.
- */
- template <int dim, int spacedim = dim>
- class FE_P_Bubbles : public dealii::FE_Poly<dim, spacedim>
- {
- public:
- /**
- * Constructor, taking the approximation degree as an argument. The
- * polynomial space is typically one degree higher than the approximation
- * space for this element: see the general documentation of this class for
- * more information.
- *
- * @note For <code>degree == 1</code> this element is equivalent to FE_P(1).
- */
- FE_P_Bubbles(const unsigned int degree);
-
- /**
- * @copydoc dealii::FiniteElement::clone()
- */
- virtual std::unique_ptr<FiniteElement<dim, spacedim>>
- clone() const override;
-
- /**
- * Return a string that uniquely identifies a finite element. This class
- * returns <tt>Simplex::FE_F_Bubbles<dim,spacedim>(degree)</tt>, with
- * @p dim, @p spacedim, and @p degree replaced by appropriate values. As
- * usual, @p spacedim is omitted in the codimension zero case.
- */
- virtual std::string
- get_name() const override;
-
- /**
- * @copydoc dealii::FiniteElement::convert_generalized_support_point_values_to_dof_values()
- */
- virtual void
- convert_generalized_support_point_values_to_dof_values(
- const std::vector<Vector<double>> &support_point_values,
- std::vector<double> & nodal_values) const override;
-
- protected:
- /**
- * Degree of the approximation (i.e., the constructor argument).
- */
- unsigned int approximation_degree;
- };
-} // namespace Simplex
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2020 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-#ifndef dealii_simplex_grid_generator_h
-#define dealii_simplex_grid_generator_h
-
-
-#include <deal.II/base/config.h>
-
-#include <deal.II/base/point.h>
-
-#include <deal.II/grid/tria.h>
-
-#include <vector>
-
-DEAL_II_NAMESPACE_OPEN
-
-namespace GridGenerator
-{
- /**
- * Create a coordinate-parallel brick from the two diagonally opposite
- * corner points @p p1 and @p p2. The number of vertices in coordinate
- * direction @p i is given by <tt>repetitions[i]+1</tt>.
- *
- * @note This function connects internally 4/8 vertices to
- * quadrilateral/hexahedral cells and subdivides these into 2/5
- * triangular/tetrahedral cells.
- *
- * @note Currently, this function only works for `dim==spacedim`.
- *
- * @ingroup simplex
- */
- template <int dim, int spacedim>
- void
- subdivided_hyper_rectangle_with_simplices(
- Triangulation<dim, spacedim> & tria,
- const std::vector<unsigned int> &repetitions,
- const Point<dim> & p1,
- const Point<dim> & p2,
- const bool colorize = false);
-
- /**
- * Initialize the given triangulation with a hypercube (square in 2D and
- * cube in 3D) consisting of @p repetitions cells in each direction.
- * The hypercube volume is the tensor product interval
- * $[left,right]^{\text{dim}}$ in the present number of dimensions, where
- * the limits are given as arguments. They default to zero and unity, then
- * producing the unit hypercube.
- *
- * @note This function connects internally 4/8 vertices to
- * quadrilateral/hexahedral cells and subdivides these into 2/5
- * triangular/tetrahedral cells.
- *
- * @ingroup simplex
- */
- template <int dim, int spacedim>
- void
- subdivided_hyper_cube_with_simplices(Triangulation<dim, spacedim> &tria,
- const unsigned int repetitions,
- const double p1 = 0.0,
- const double p2 = 1.0,
- const bool colorize = false);
-
-} // namespace GridGenerator
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2020 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-#ifndef dealii_tet_quadrature_lib_h
-#define dealii_tet_quadrature_lib_h
-
-
-#include <deal.II/base/config.h>
-
-#include <deal.II/base/quadrature_lib.h>
-
-DEAL_II_NAMESPACE_OPEN
-
-namespace Simplex
-{
- /**
- * Integration rule for simplex entities.
- *
- * Users specify a number `n_points_1D` as an indication of what polynomial
- * degree to be integrated exactly, similarly to the number of points in a
- * QGauss quadrature object, even though the present quadrature formula is not
- * a tensor product. The given value is translated for n_points_1D=1,2,3,4 to
- * following number of quadrature points for 2D and 3D:
- * - 2D: 1, 3, 7, 15
- * - 3D: 1, 4, 10, 35
- *
- * For 1D, the quadrature rule degenerates to a `QGauss<1>(n_points_1D)`.
- *
- * @ingroup simplex
- */
- template <int dim>
- class QGauss : public QSimplex<dim>
- {
- public:
- /**
- * Constructor taking the number of quadrature points in 1D direction
- * @p n_points_1D.
- */
- explicit QGauss(const unsigned int n_points_1D);
- };
-
- /**
- * Witherden-Vincent rules for simplex entities.
- *
- * Like QGauss, users should specify a number `n_points_1D` as an indication
- * of what polynomial degree to be integrated exactly (e.g., for $n$ points,
- * the rule can integrate polynomials of degree $2 n - 1$ exactly). The given
- * value for n_points_1D = 1, 2, 3, 4, 5 results in the following number of
- * quadrature points in 2D and 3D:
- * - 2D: 1, 6, 7, 15, 19
- * - 3D: 1, 8, 14, 35, 59
- *
- * For 1D, the quadrature rule degenerates to a `QGauss<1>(n_points_1D)`.
- *
- * These rules match the ones listed for Witherden-Vincent in the quadpy
- * @cite quadpy library and were first described in
- * @cite witherden2015identification.
- *
- * @ingroup simplex
- */
- template <int dim>
- class QWitherdenVincent : public QSimplex<dim>
- {
- public:
- /**
- * Constructor taking the number of quadrature points in 1D direction
- * @p n_points_1D.
- */
- explicit QWitherdenVincent(const unsigned int n_points_1D);
- };
-
- /**
- * Integration rule for wedge entities.
- */
- template <int dim>
- class QGaussWedge : public Quadrature<dim>
- {
- public:
- /**
- * Users specify a number `n_points_1D` as an indication of what polynomial
- * degree to be integrated exactly. For details, see the comments of
- * Simplex::QGauss.
- */
- explicit QGaussWedge(const unsigned int n_points_1D);
- };
-
- /**
- * Integration rule for pyramid entities.
- */
- template <int dim>
- class QGaussPyramid : public Quadrature<dim>
- {
- public:
- /**
- * Users specify a number `n_points_1D` as an indication of what polynomial
- * degree to be integrated exactly. For details, see the comments of
- * Simplex::QGauss.
- */
- explicit QGaussPyramid(const unsigned int n_points_1D);
- };
-
-} // namespace Simplex
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
ADD_SUBDIRECTORY(physics)
ADD_SUBDIRECTORY(optimization/rol)
ADD_SUBDIRECTORY(non_matching)
-ADD_SUBDIRECTORY(simplex)
ADD_SUBDIRECTORY(sundials)
FOREACH(build ${DEAL_II_BUILD_TYPES})
polynomials_abf.cc
polynomials_adini.cc
tensor_polynomials_base.cc
+ polynomials_barycentric.cc
polynomials_bernstein.cc
polynomials_bdm.cc
polynomials_bernardi_raugel.cc
polynomial_space.cc
polynomials_p.cc
polynomials_piecewise.cc
+ polynomials_pyramid.cc
polynomials_rannacher_turek.cc
polynomials_raviart_thomas.cc
polynomials_rt_bubbles.cc
+ polynomials_wedge.cc
process_grid.cc
qprojector.cc
quadrature.cc
// ---------------------------------------------------------------------
-#include <deal.II/simplex/barycentric_polynomials.h>
+#include <deal.II/base/polynomials_barycentric.h>
DEAL_II_NAMESPACE_OPEN
// ---------------------------------------------------------------------
-#include <deal.II/base/ndarray.h>
-
-#include <deal.II/simplex/barycentric_polynomials.h>
-#include <deal.II/simplex/polynomials.h>
+#include <deal.II/base/polynomials_barycentric.h>
+#include <deal.II/base/polynomials_pyramid.h>
DEAL_II_NAMESPACE_OPEN
return 0;
}
-
- unsigned int
- compute_n_polynomials_wedge(const unsigned int dim,
- const unsigned int degree)
- {
- if (dim == 3)
- {
- if (degree == 1)
- return 6;
- if (degree == 2)
- return 18;
- }
-
- Assert(false, ExcNotImplemented());
-
- return 0;
- }
} // namespace
- template <int dim>
- ScalarWedgePolynomial<dim>::ScalarWedgePolynomial(const unsigned int degree)
- : ScalarPolynomialsBase<dim>(degree,
- compute_n_polynomials_wedge(dim, degree))
- , poly_tri(BarycentricPolynomials<2>::get_fe_p_basis(degree))
- , poly_line(BarycentricPolynomials<1>::get_fe_p_basis(degree))
- {}
-
-
- namespace
- {
- /**
- * Decompose the shape-function index of a linear wedge into an index
- * to access the right shape function within the triangle and and within
- * the line.
- */
- static const constexpr ndarray<unsigned int, 6, 2> wedge_table_1{
- {{{0, 0}}, {{1, 0}}, {{2, 0}}, {{0, 1}}, {{1, 1}}, {{2, 1}}}};
-
- /**
- * Decompose the shape-function index of a quadratic wedge into an index
- * to access the right shape function within the triangle and and within
- * the line.
- */
- static const constexpr ndarray<unsigned int, 18, 2> wedge_table_2{
- {{{0, 0}},
- {{1, 0}},
- {{2, 0}},
- {{0, 1}},
- {{1, 1}},
- {{2, 1}},
- {{3, 0}},
- {{4, 0}},
- {{5, 0}},
- {{3, 1}},
- {{4, 1}},
- {{5, 1}},
- {{0, 2}},
- {{1, 2}},
- {{2, 2}},
- {{3, 2}},
- {{4, 2}},
- {{5, 2}}}};
- } // namespace
-
-
-
- template <int dim>
- double
- ScalarWedgePolynomial<dim>::compute_value(const unsigned int i,
- const Point<dim> & p) const
- {
- const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i];
-
- const Point<2> p_tri(p[0], p[1]);
- const auto v_tri = poly_tri.compute_value(pair[0], p_tri);
-
- const Point<1> p_line(p[2]);
- const auto v_line = poly_line.compute_value(pair[1], p_line);
-
- return v_tri * v_line;
- }
-
-
-
- template <int dim>
- Tensor<1, dim>
- ScalarWedgePolynomial<dim>::compute_grad(const unsigned int i,
- const Point<dim> & p) const
- {
- const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i];
-
- const Point<2> p_tri(p[0], p[1]);
- const auto v_tri = poly_tri.compute_value(pair[0], p_tri);
- const auto g_tri = poly_tri.compute_grad(pair[0], p_tri);
-
- const Point<1> p_line(p[2]);
- const auto v_line = poly_line.compute_value(pair[1], p_line);
- const auto g_line = poly_line.compute_grad(pair[1], p_line);
-
- Tensor<1, dim> grad;
- grad[0] = g_tri[0] * v_line;
- grad[1] = g_tri[1] * v_line;
- grad[2] = v_tri * g_line[0];
-
- return grad;
- }
-
-
-
- template <int dim>
- Tensor<2, dim>
- ScalarWedgePolynomial<dim>::compute_grad_grad(const unsigned int i,
- const Point<dim> & p) const
- {
- (void)i;
- (void)p;
-
- Assert(false, ExcNotImplemented());
- return Tensor<2, dim>();
- }
-
-
-
- template <int dim>
- void
- ScalarWedgePolynomial<dim>::evaluate(
- const Point<dim> & unit_point,
- std::vector<double> & values,
- std::vector<Tensor<1, dim>> &grads,
- std::vector<Tensor<2, dim>> &grad_grads,
- std::vector<Tensor<3, dim>> &third_derivatives,
- std::vector<Tensor<4, dim>> &fourth_derivatives) const
- {
- (void)grads;
- (void)grad_grads;
- (void)third_derivatives;
- (void)fourth_derivatives;
-
- if (values.size() == this->n())
- for (unsigned int i = 0; i < this->n(); i++)
- values[i] = compute_value(i, unit_point);
-
- if (grads.size() == this->n())
- for (unsigned int i = 0; i < this->n(); i++)
- grads[i] = compute_grad(i, unit_point);
- }
-
-
-
- template <int dim>
- Tensor<1, dim>
- ScalarWedgePolynomial<dim>::compute_1st_derivative(const unsigned int i,
- const Point<dim> & p) const
- {
- return compute_grad(i, p);
- }
-
-
-
- template <int dim>
- Tensor<2, dim>
- ScalarWedgePolynomial<dim>::compute_2nd_derivative(const unsigned int i,
- const Point<dim> & p) const
- {
- (void)i;
- (void)p;
-
- Assert(false, ExcNotImplemented());
-
- return {};
- }
-
-
-
- template <int dim>
- Tensor<3, dim>
- ScalarWedgePolynomial<dim>::compute_3rd_derivative(const unsigned int i,
- const Point<dim> & p) const
- {
- (void)i;
- (void)p;
-
- Assert(false, ExcNotImplemented());
-
- return {};
- }
-
-
-
- template <int dim>
- Tensor<4, dim>
- ScalarWedgePolynomial<dim>::compute_4th_derivative(const unsigned int i,
- const Point<dim> & p) const
- {
- (void)i;
- (void)p;
-
- Assert(false, ExcNotImplemented());
-
- return {};
- }
-
-
-
- template <int dim>
- std::string
- ScalarWedgePolynomial<dim>::name() const
- {
- return "ScalarWedgePolynomial";
- }
-
-
-
- template <int dim>
- std::unique_ptr<ScalarPolynomialsBase<dim>>
- ScalarWedgePolynomial<dim>::clone() const
- {
- return std::make_unique<ScalarWedgePolynomial<dim>>(*this);
- }
-
-
-
template <int dim>
ScalarPyramidPolynomial<dim>::ScalarPyramidPolynomial(
const unsigned int degree)
- template class ScalarWedgePolynomial<1>;
- template class ScalarWedgePolynomial<2>;
- template class ScalarWedgePolynomial<3>;
template class ScalarPyramidPolynomial<1>;
template class ScalarPyramidPolynomial<2>;
template class ScalarPyramidPolynomial<3>;
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/base/ndarray.h>
+#include <deal.II/base/polynomials_barycentric.h>
+#include <deal.II/base/polynomials_wedge.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Simplex
+{
+ namespace
+ {
+ unsigned int
+ compute_n_polynomials_wedge(const unsigned int dim,
+ const unsigned int degree)
+ {
+ if (dim == 3)
+ {
+ if (degree == 1)
+ return 6;
+ if (degree == 2)
+ return 18;
+ }
+
+ Assert(false, ExcNotImplemented());
+
+ return 0;
+ }
+ } // namespace
+
+
+
+ template <int dim>
+ ScalarWedgePolynomial<dim>::ScalarWedgePolynomial(const unsigned int degree)
+ : ScalarPolynomialsBase<dim>(degree,
+ compute_n_polynomials_wedge(dim, degree))
+ , poly_tri(BarycentricPolynomials<2>::get_fe_p_basis(degree))
+ , poly_line(BarycentricPolynomials<1>::get_fe_p_basis(degree))
+ {}
+
+
+ namespace
+ {
+ /**
+ * Decompose the shape-function index of a linear wedge into an index
+ * to access the right shape function within the triangle and and within
+ * the line.
+ */
+ static const constexpr ndarray<unsigned int, 6, 2> wedge_table_1{
+ {{{0, 0}}, {{1, 0}}, {{2, 0}}, {{0, 1}}, {{1, 1}}, {{2, 1}}}};
+
+ /**
+ * Decompose the shape-function index of a quadratic wedge into an index
+ * to access the right shape function within the triangle and and within
+ * the line.
+ */
+ static const constexpr ndarray<unsigned int, 18, 2> wedge_table_2{
+ {{{0, 0}},
+ {{1, 0}},
+ {{2, 0}},
+ {{0, 1}},
+ {{1, 1}},
+ {{2, 1}},
+ {{3, 0}},
+ {{4, 0}},
+ {{5, 0}},
+ {{3, 1}},
+ {{4, 1}},
+ {{5, 1}},
+ {{0, 2}},
+ {{1, 2}},
+ {{2, 2}},
+ {{3, 2}},
+ {{4, 2}},
+ {{5, 2}}}};
+ } // namespace
+
+
+
+ template <int dim>
+ double
+ ScalarWedgePolynomial<dim>::compute_value(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i];
+
+ const Point<2> p_tri(p[0], p[1]);
+ const auto v_tri = poly_tri.compute_value(pair[0], p_tri);
+
+ const Point<1> p_line(p[2]);
+ const auto v_line = poly_line.compute_value(pair[1], p_line);
+
+ return v_tri * v_line;
+ }
+
+
+
+ template <int dim>
+ Tensor<1, dim>
+ ScalarWedgePolynomial<dim>::compute_grad(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i];
+
+ const Point<2> p_tri(p[0], p[1]);
+ const auto v_tri = poly_tri.compute_value(pair[0], p_tri);
+ const auto g_tri = poly_tri.compute_grad(pair[0], p_tri);
+
+ const Point<1> p_line(p[2]);
+ const auto v_line = poly_line.compute_value(pair[1], p_line);
+ const auto g_line = poly_line.compute_grad(pair[1], p_line);
+
+ Tensor<1, dim> grad;
+ grad[0] = g_tri[0] * v_line;
+ grad[1] = g_tri[1] * v_line;
+ grad[2] = v_tri * g_line[0];
+
+ return grad;
+ }
+
+
+
+ template <int dim>
+ Tensor<2, dim>
+ ScalarWedgePolynomial<dim>::compute_grad_grad(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ (void)i;
+ (void)p;
+
+ Assert(false, ExcNotImplemented());
+ return Tensor<2, dim>();
+ }
+
+
+
+ template <int dim>
+ void
+ ScalarWedgePolynomial<dim>::evaluate(
+ const Point<dim> & unit_point,
+ std::vector<double> & values,
+ std::vector<Tensor<1, dim>> &grads,
+ std::vector<Tensor<2, dim>> &grad_grads,
+ std::vector<Tensor<3, dim>> &third_derivatives,
+ std::vector<Tensor<4, dim>> &fourth_derivatives) const
+ {
+ (void)grads;
+ (void)grad_grads;
+ (void)third_derivatives;
+ (void)fourth_derivatives;
+
+ if (values.size() == this->n())
+ for (unsigned int i = 0; i < this->n(); i++)
+ values[i] = compute_value(i, unit_point);
+
+ if (grads.size() == this->n())
+ for (unsigned int i = 0; i < this->n(); i++)
+ grads[i] = compute_grad(i, unit_point);
+ }
+
+
+
+ template <int dim>
+ Tensor<1, dim>
+ ScalarWedgePolynomial<dim>::compute_1st_derivative(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ return compute_grad(i, p);
+ }
+
+
+
+ template <int dim>
+ Tensor<2, dim>
+ ScalarWedgePolynomial<dim>::compute_2nd_derivative(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ (void)i;
+ (void)p;
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+
+
+ template <int dim>
+ Tensor<3, dim>
+ ScalarWedgePolynomial<dim>::compute_3rd_derivative(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ (void)i;
+ (void)p;
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+
+
+ template <int dim>
+ Tensor<4, dim>
+ ScalarWedgePolynomial<dim>::compute_4th_derivative(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ (void)i;
+ (void)p;
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+
+
+ template <int dim>
+ std::string
+ ScalarWedgePolynomial<dim>::name() const
+ {
+ return "ScalarWedgePolynomial";
+ }
+
+
+
+ template <int dim>
+ std::unique_ptr<ScalarPolynomialsBase<dim>>
+ ScalarWedgePolynomial<dim>::clone() const
+ {
+ return std::make_unique<ScalarWedgePolynomial<dim>>(*this);
+ }
+
+
+
+ template class ScalarWedgePolynomial<1>;
+ template class ScalarWedgePolynomial<2>;
+ template class ScalarWedgePolynomial<3>;
+} // namespace Simplex
+
+DEAL_II_NAMESPACE_CLOSE
#include <deal.II/base/derivative_form.h>
#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/polynomials_barycentric.h>
#include <deal.II/base/qprojector.h>
#include <deal.II/base/tensor_product_polynomials.h>
-#include <deal.II/simplex/barycentric_polynomials.h>
-
DEAL_II_NAMESPACE_OPEN
+namespace Simplex
+{
+ template <int dim>
+ QGauss<dim>::QGauss(const unsigned int n_points_1D)
+ : QSimplex<dim>(Quadrature<dim>())
+ {
+ // fill quadrature points and quadrature weights
+ if (dim == 1)
+ {
+ const dealii::QGauss<dim> quad(n_points_1D);
+
+ this->quadrature_points = quad.get_points();
+ this->weights = quad.get_weights();
+ }
+ else if (dim == 2)
+ {
+ if (n_points_1D == 1)
+ {
+ const double p = 1.0 / 3.0;
+ this->quadrature_points.emplace_back(p, p);
+ this->weights.emplace_back(0.5);
+ }
+ else if (n_points_1D == 2)
+ {
+ const double Q23 = 2.0 / 3.0;
+ const double Q16 = 1.0 / 6.0;
+
+ this->quadrature_points.emplace_back(Q23, Q16);
+ this->quadrature_points.emplace_back(Q16, Q23);
+ this->quadrature_points.emplace_back(Q16, Q16);
+ this->weights.emplace_back(Q16);
+ this->weights.emplace_back(Q16);
+ this->weights.emplace_back(Q16);
+ }
+ else if (n_points_1D == 3)
+ {
+ const double q12 = 0.5;
+
+ // clang-format off
+ this->quadrature_points.emplace_back(0.3333333333330, 0.3333333333330);
+ this->quadrature_points.emplace_back(0.7974269853530, 0.1012865073230);
+ this->quadrature_points.emplace_back(0.1012865073230, 0.7974269853530);
+ this->quadrature_points.emplace_back(0.1012865073230, 0.1012865073230);
+ this->quadrature_points.emplace_back(0.0597158717898, 0.4701420641050);
+ this->quadrature_points.emplace_back(0.4701420641050, 0.0597158717898);
+ this->quadrature_points.emplace_back(0.4701420641050, 0.4701420641050);
+ // clang-format on
+
+ this->weights.emplace_back(q12 * 0.225);
+ this->weights.emplace_back(q12 * 0.125939180545);
+ this->weights.emplace_back(q12 * 0.125939180545);
+ this->weights.emplace_back(q12 * 0.125939180545);
+ this->weights.emplace_back(q12 * 0.132394152789);
+ this->weights.emplace_back(q12 * 0.132394152789);
+ this->weights.emplace_back(q12 * 0.132394152789);
+ }
+ else if (n_points_1D == 4)
+ {
+ Quadrature<dim>::operator=(QWitherdenVincent<dim>(n_points_1D));
+ }
+ }
+ else if (dim == 3)
+ {
+ if (n_points_1D == 1)
+ {
+ const double Q14 = 1.0 / 4.0;
+ const double Q16 = 1.0 / 6.0;
+
+ this->quadrature_points.emplace_back(Q14, Q14, Q14);
+ this->weights.emplace_back(Q16);
+ }
+ else if (n_points_1D == 2)
+ {
+ const double Q124 = 1.0 / 6.0 / 4.0;
+
+ const double palpha = (5.0 + 3.0 * sqrt(5.0)) / 20.0;
+ const double pbeta = (5.0 - sqrt(5.0)) / 20.0;
+ this->quadrature_points.emplace_back(pbeta, pbeta, pbeta);
+ this->quadrature_points.emplace_back(palpha, pbeta, pbeta);
+ this->quadrature_points.emplace_back(pbeta, palpha, pbeta);
+ this->quadrature_points.emplace_back(pbeta, pbeta, palpha);
+ this->weights.emplace_back(Q124);
+ this->weights.emplace_back(Q124);
+ this->weights.emplace_back(Q124);
+ this->weights.emplace_back(Q124);
+ }
+ else if (n_points_1D == 3)
+ {
+ const double Q16 = 1.0 / 6.0;
+
+ // clang-format off
+ this->quadrature_points.emplace_back(0.5684305841968444, 0.1438564719343852, 0.1438564719343852);
+ this->quadrature_points.emplace_back(0.1438564719343852, 0.1438564719343852, 0.1438564719343852);
+ this->quadrature_points.emplace_back(0.1438564719343852, 0.1438564719343852, 0.5684305841968444);
+ this->quadrature_points.emplace_back(0.1438564719343852, 0.5684305841968444, 0.1438564719343852);
+ this->quadrature_points.emplace_back(0.0000000000000000, 0.5000000000000000, 0.5000000000000000);
+ this->quadrature_points.emplace_back(0.5000000000000000, 0.0000000000000000, 0.5000000000000000);
+ this->quadrature_points.emplace_back(0.5000000000000000, 0.5000000000000000, 0.0000000000000000);
+ this->quadrature_points.emplace_back(0.5000000000000000, 0.0000000000000000, 0.0000000000000000);
+ this->quadrature_points.emplace_back(0.0000000000000000, 0.5000000000000000, 0.0000000000000000);
+ this->quadrature_points.emplace_back(0.0000000000000000, 0.0000000000000000, 0.5000000000000000);
+ // clang-format on
+
+ this->weights.emplace_back(0.2177650698804054 * Q16);
+ this->weights.emplace_back(0.2177650698804054 * Q16);
+ this->weights.emplace_back(0.2177650698804054 * Q16);
+ this->weights.emplace_back(0.2177650698804054 * Q16);
+ this->weights.emplace_back(0.0214899534130631 * Q16);
+ this->weights.emplace_back(0.0214899534130631 * Q16);
+ this->weights.emplace_back(0.0214899534130631 * Q16);
+ this->weights.emplace_back(0.0214899534130631 * Q16);
+ this->weights.emplace_back(0.0214899534130631 * Q16);
+ this->weights.emplace_back(0.0214899534130631 * Q16);
+ }
+ else if (n_points_1D == 4)
+ {
+ Quadrature<dim>::operator=(QWitherdenVincent<dim>(n_points_1D));
+ }
+ }
+
+ AssertDimension(this->quadrature_points.size(), this->weights.size());
+ Assert(this->quadrature_points.size() > 0,
+ ExcNotImplemented(
+ "Simplex::QGauss is currently only implemented for "
+ "n_points_1D = 1, 2, 3, and 4 while you are asking for "
+ "n_points_1D = " +
+ Utilities::to_string(n_points_1D)));
+ }
+
+ namespace
+ {
+ template <std::size_t b_dim>
+ std::vector<std::array<double, b_dim>>
+ all_permutations(const std::array<double, b_dim> &b_point)
+ {
+ std::vector<std::array<double, b_dim>> output;
+
+ // We want all possible permutations of the barycentric coordinates.
+ // The easiest way to get all of them is to sort the input first and
+ // then use next_permutation to cycle through them all.
+ std::array<double, b_dim> temp = b_point;
+ std::sort(temp.begin(), temp.end());
+ do
+ {
+ output.push_back(temp);
+ }
+ while (std::next_permutation(temp.begin(), temp.end()));
+
+ return output;
+ }
+ } // namespace
+
+
+
+ template <int dim>
+ QWitherdenVincent<dim>::QWitherdenVincent(const unsigned int n_points_1D)
+ : QSimplex<dim>(Quadrature<dim>())
+ {
+ Assert(1 <= dim && dim <= 3, ExcNotImplemented());
+ // Just use Gauss in 1D: this is a high-order open rule so this is a
+ // reasonable equivalent for generic programming.
+ if (dim == 1)
+ {
+ Quadrature<dim>::operator=(dealii::QGauss<dim>(n_points_1D));
+ return;
+ }
+
+ std::array<double, dim + 1> centroid;
+ std::fill(centroid.begin(), centroid.end(), 1.0 / (dim + 1.0));
+ std::vector<std::vector<std::array<double, dim + 1>>> b_point_permutations;
+ std::vector<double> b_weights;
+
+ // We can simplify the implementation of these quadrature rules
+ // by quite a bit by exploiting symmetry - we do essentially the
+ // same thing for each barycentric coordinate, so we can express
+ // our quadrature rule as permutations of barycentric points
+ // instead of writing things out explicitly.
+
+ // Apply a Barycentric permutation where one point is different.
+ auto process_point_1 = [&](const double a, const double w) {
+ const double b = 1.0 - dim * a;
+ std::array<double, dim + 1> b_point;
+ std::fill(b_point.begin(), b_point.begin() + dim, a);
+ b_point[dim] = b;
+
+ b_weights.push_back(w);
+ b_point_permutations.push_back(all_permutations(b_point));
+ };
+
+ // Apply a Barycentric permutation where two points (in 3D) are different.
+ auto process_point_2 = [&](const double a, const double w) {
+ Assert(dim == 3, ExcInternalError());
+ const double b = (1.0 - 2.0 * a) / 2.0;
+ std::array<double, dim + 1> b_point;
+ std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
+ b_point[dim - 1] = b;
+ b_point[dim] = b;
+
+ b_weights.push_back(w);
+ b_point_permutations.push_back(all_permutations(b_point));
+ };
+
+ // Apply a Barycentric permutation where three (or four) points
+ // are different (since there are two inputs).
+ auto process_point_3 = [&](const double a, const double b, const double w) {
+ const double c = 1.0 - (dim - 1.0) * a - b;
+ std::array<double, dim + 1> b_point;
+ std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
+ b_point[dim - 1] = b;
+ b_point[dim] = c;
+
+ b_weights.push_back(w);
+ b_point_permutations.push_back(all_permutations(b_point));
+ };
+
+ if (n_points_1D == 1)
+ {
+ b_point_permutations.push_back({centroid});
+ b_weights.push_back(1.0);
+ }
+ else if (n_points_1D == 2)
+ {
+ // This is WV-4 in 2D and WV-3 in 3D
+ if (dim == 2)
+ {
+ process_point_1(9.1576213509770743e-02, 1.0995174365532187e-01);
+ process_point_1(4.4594849091596489e-01, 2.2338158967801147e-01);
+ }
+ else if (dim == 3)
+ {
+ process_point_1(3.281633025163817e-01, 1.362178425370874e-01);
+ process_point_1(1.080472498984286e-01, 1.137821574629126e-01);
+ }
+ }
+ else if (n_points_1D == 3)
+ {
+ // This is the WV-5 rule in both 2D and 3D
+ if (dim == 2)
+ {
+ b_weights.push_back(0.225);
+ b_point_permutations.push_back({centroid});
+
+ process_point_1(1.0128650732345634e-01, 1.2593918054482714e-01);
+ process_point_1(4.7014206410511511e-01, 1.3239415278850619e-01);
+ }
+ else if (dim == 3)
+ {
+ process_point_1(3.108859192633006e-01, 1.126879257180159e-01);
+ process_point_1(9.273525031089125e-02, 7.349304311636196e-02);
+
+ process_point_2(4.550370412564964e-02, 4.254602077708147e-02);
+ }
+ }
+ else if (n_points_1D == 4)
+ {
+ // This is the WV-7 rule in both 2D and 3D
+ if (dim == 2)
+ {
+ process_point_1(3.3730648554587850e-02, 1.6545050110792131e-02);
+ process_point_1(4.7430969250471822e-01, 7.7086646185986069e-02);
+ process_point_1(2.4157738259540357e-01, 1.2794417123015558e-01);
+ process_point_3(4.7036644652595216e-02,
+ 1.9868331479735168e-01,
+ 5.5878732903199779e-02);
+ }
+ else if (dim == 3)
+ {
+ b_point_permutations.push_back({centroid});
+ b_weights.push_back(9.548528946413085e-02);
+
+ process_point_1(3.157011497782028e-01, 4.232958120996703e-02);
+ process_point_2(5.048982259839635e-02, 3.189692783285758e-02);
+
+ process_point_3(1.888338310260010e-01,
+ 5.751716375870000e-01,
+ 3.720713072833462e-02);
+ process_point_3(2.126547254148314e-02,
+ 8.108302410985486e-01,
+ 8.110770829903342e-03);
+ }
+ }
+ else if (n_points_1D == 5)
+ {
+ // This is the WV-9 rule in both 2D and 3D
+ if (dim == 2)
+ {
+ b_point_permutations.push_back({centroid});
+ b_weights.push_back(9.7135796282798836e-02);
+
+ process_point_1(4.4729513394452691e-02, 2.5577675658698031e-02);
+ process_point_1(4.8968251919873762e-01, 3.1334700227139071e-02);
+ process_point_1(4.3708959149293664e-01, 7.7827541004774278e-02);
+ process_point_1(1.8820353561903275e-01, 7.9647738927210249e-02);
+
+ process_point_3(3.6838412054736258e-02,
+ 2.2196298916076568e-01,
+ 4.3283539377289376e-02);
+ }
+ else if (dim == 3)
+ {
+ b_point_permutations.push_back({centroid});
+ b_weights.push_back(5.801054891248025e-02);
+
+ process_point_1(6.198169755222693e-10, 6.431928175925639e-05);
+ process_point_1(1.607745353952616e-01, 2.317333846242546e-02);
+ process_point_1(3.222765218214210e-01, 2.956291233542929e-02);
+ process_point_1(4.510891834541358e-02, 8.063979979616182e-03);
+
+ process_point_2(1.122965460043761e-01, 3.813408010370246e-02);
+
+ process_point_3(4.588714487524592e-01,
+ 2.554579233041310e-03,
+ 8.384422198298552e-03);
+ process_point_3(3.377587068533860e-02,
+ 7.183503264420745e-01,
+ 1.023455935274533e-02);
+ process_point_3(1.836413698099279e-01,
+ 3.441591057817528e-02,
+ 2.052491596798814e-02);
+ }
+ }
+ else if (n_points_1D == 6)
+ {
+ // There is no WV-11 rule in 3D yet
+ if (dim == 2)
+ {
+ b_point_permutations.push_back({centroid});
+ b_weights.push_back(8.5761179732224219e-02);
+
+ process_point_1(2.8485417614371900e-02, 1.0431870512894697e-02);
+ process_point_1(4.9589190096589092e-01, 1.6606273054585369e-02);
+ process_point_1(1.0263548271224643e-01, 3.8630759237019321e-02);
+ process_point_1(4.3846592676435220e-01, 6.7316154079468296e-02);
+ process_point_1(2.1021995670317828e-01, 7.0515684111716576e-02);
+
+ process_point_3(7.3254276860644785e-03,
+ 1.4932478865208237e-01,
+ 1.0290289572953278e-02);
+ process_point_3(4.6010500165429957e-02,
+ 2.8958112563770588e-01,
+ 4.0332476640500554e-02);
+ }
+ else if (dim == 3)
+ {
+ Assert(false, ExcNotImplemented());
+ }
+ }
+ else
+ {
+ Assert(false, ExcNotImplemented());
+ }
+
+ Assert(b_point_permutations.size() == b_weights.size(), ExcInternalError());
+ for (unsigned int permutation_n = 0; permutation_n < b_weights.size();
+ ++permutation_n)
+ {
+ for (const std::array<double, dim + 1> &b_point :
+ b_point_permutations[permutation_n])
+ {
+ const double volume = (dim == 2 ? 1.0 / 2.0 : 1.0 / 6.0);
+ this->weights.emplace_back(volume * b_weights[permutation_n]);
+ Point<dim> c_point;
+ std::copy(b_point.begin(),
+ b_point.begin() + dim,
+ c_point.begin_raw());
+ this->quadrature_points.emplace_back(c_point);
+ }
+ }
+ }
+
+
+
+ template <int dim>
+ QGaussWedge<dim>::QGaussWedge(const unsigned int n_points)
+ : Quadrature<dim>()
+ {
+ AssertDimension(dim, 3);
+
+ Simplex::QGauss<2> quad_tri(n_points);
+ QGauss<1> quad_line(n_points);
+
+ for (unsigned int i = 0; i < quad_line.size(); ++i)
+ for (unsigned int j = 0; j < quad_tri.size(); ++j)
+ {
+ this->quadrature_points.emplace_back(quad_tri.point(j)[0],
+ quad_tri.point(j)[1],
+ quad_line.point(i)[0]);
+ this->weights.emplace_back(quad_tri.weight(j) * quad_line.weight(i));
+ }
+
+ AssertDimension(this->quadrature_points.size(), this->weights.size());
+ Assert(this->quadrature_points.size() > 0,
+ ExcMessage("No valid quadrature points!"));
+ }
+
+
+
+ template <int dim>
+ QGaussPyramid<dim>::QGaussPyramid(const unsigned int n_points_1D)
+ : Quadrature<dim>()
+ {
+ AssertDimension(dim, 3);
+
+ if (n_points_1D == 1)
+ {
+ const double Q14 = 1.0 / 4.0;
+ const double Q43 = 4.0 / 3.0;
+
+ this->quadrature_points.emplace_back(0, 0, Q14);
+ this->weights.emplace_back(Q43);
+ }
+ else if (n_points_1D == 2)
+ {
+ // clang-format off
+ this->quadrature_points.emplace_back(-0.26318405556971, -0.26318405556971, 0.54415184401122);
+ this->quadrature_points.emplace_back(-0.50661630334979, -0.50661630334979, 0.12251482265544);
+ this->quadrature_points.emplace_back(-0.26318405556971, +0.26318405556971, 0.54415184401122);
+ this->quadrature_points.emplace_back(-0.50661630334979, +0.50661630334979, 0.12251482265544);
+ this->quadrature_points.emplace_back(+0.26318405556971, -0.26318405556971, 0.54415184401122);
+ this->quadrature_points.emplace_back(+0.50661630334979, -0.50661630334979, 0.12251482265544);
+ this->quadrature_points.emplace_back(+0.26318405556971, +0.26318405556971, 0.54415184401122);
+ this->quadrature_points.emplace_back(+0.50661630334979, +0.50661630334979, 0.12251482265544);
+ // clang-format on
+
+ this->weights.emplace_back(0.10078588207983);
+ this->weights.emplace_back(0.23254745125351);
+ this->weights.emplace_back(0.10078588207983);
+ this->weights.emplace_back(0.23254745125351);
+ this->weights.emplace_back(0.10078588207983);
+ this->weights.emplace_back(0.23254745125351);
+ this->weights.emplace_back(0.10078588207983);
+ this->weights.emplace_back(0.23254745125351);
+ }
+
+ AssertDimension(this->quadrature_points.size(), this->weights.size());
+ Assert(this->quadrature_points.size() > 0,
+ ExcMessage("No valid quadrature points!"));
+ }
+
+} // namespace Simplex
+
+
+
// explicit specialization
// note that 1d formulae are specialized by implementation above
template class QGauss<2>;
template class QSplit<2>;
template class QSplit<3>;
+template class Simplex::QGauss<1>;
+template class Simplex::QGauss<2>;
+template class Simplex::QGauss<3>;
+template class Simplex::QGaussWedge<1>;
+template class Simplex::QGaussWedge<2>;
+template class Simplex::QGaussWedge<3>;
+template class Simplex::QGaussPyramid<1>;
+template class Simplex::QGaussPyramid<2>;
+template class Simplex::QGaussPyramid<3>;
+
+template class Simplex::QWitherdenVincent<1>;
+template class Simplex::QWitherdenVincent<2>;
+template class Simplex::QWitherdenVincent<3>;
+
DEAL_II_NAMESPACE_CLOSE
fe_dgp_nonparametric.cc
fe_dgq.cc
fe_dg_vector.cc
+ fe_enriched.cc
fe_face.cc
fe_nedelec.cc
fe_nedelec_sz.cc
fe_nothing.cc
fe_poly.cc
fe_poly_tensor.cc
+ fe_pyramid_p.cc
fe_p1nc.cc
fe_q_base.cc
fe_q.cc
fe_series_fourier.cc
fe_series_legendre.cc
fe_system.cc
- fe_enriched.cc
+ fe_simplex_p.cc
+ fe_simplex_p_bubbles.cc
fe_trace.cc
+ fe_wedge_p.cc
mapping_c1.cc
mapping_cartesian.cc
mapping.cc
fe_nothing.inst.in
fe_poly.inst.in
fe_poly_tensor.inst.in
+ fe_pyramid_p.inst.in
fe_q_base.inst.in
fe_q_bubbles.inst.in
fe_q_dg0.inst.in
fe_rt_bubbles.inst.in
fe_series_fourier.inst.in
fe_series_legendre.inst.in
+ fe_simplex_p.inst.in
+ fe_simplex_p_bubbles.inst.in
fe_system.inst.in
fe_enriched.inst.in
fe_tools.inst.in
fe_values.impl.1.inst.in
fe_values.impl.2.inst.in
fe_values.inst.in
+ fe_wedge_p.inst.in
mapping_c1.inst.in
mapping_cartesian.inst.in
mapping.inst.in
#include <deal.II/fe/fe_bernstein.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_q_bubbles.h>
#include <deal.II/fe/fe_q_dg0.h>
#include <deal.II/fe/fe_q_hierarchical.h>
#include <deal.II/fe/fe_q_iso_q1.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/lac/vector.h>
-#include <deal.II/simplex/fe_lib.h>
-
#include <iostream>
#include <memory>
#include <sstream>
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/polynomials_barycentric.h>
+#include <deal.II/base/qprojector.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_wedge_p.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Simplex
+{
+ namespace
+ {
+ /**
+ * Helper function to set up the dpo vector of FE_PyramidP for a given @p degree.
+ */
+ internal::GenericDoFsPerObject
+ get_dpo_vector_fe_pyramid_p(const unsigned int degree)
+ {
+ internal::GenericDoFsPerObject dpo;
+
+ if (degree == 1)
+ {
+ dpo.dofs_per_object_exclusive = {{1}, {0}, {0, 0, 0, 0, 0}, {0}};
+ dpo.dofs_per_object_inclusive = {{1}, {2}, {4, 3, 3, 3, 3}, {5}};
+ dpo.object_index = {{}, {5}, {5}, {5}};
+ dpo.first_object_index_on_face = {{},
+ {4, 3, 3, 3, 3},
+ {4, 3, 3, 3, 3}};
+ }
+ else
+ {
+ Assert(false, ExcNotImplemented());
+ }
+
+ return dpo;
+ }
+
+ /**
+ * Helper function to set up the dpo vector of FE_PyramidDGP for a given @p degree.
+ */
+ internal::GenericDoFsPerObject
+ get_dpo_vector_fe_pyramid_dgp(const unsigned int degree)
+ {
+ unsigned int n_dofs = 0;
+
+ if (degree == 1)
+ n_dofs = 5;
+ else
+ Assert(false, ExcNotImplemented());
+
+ return internal::expand(3, {{0, 0, 0, n_dofs}}, ReferenceCells::Pyramid);
+ }
+ } // namespace
+
+
+ template <int dim, int spacedim>
+ FE_Pyramid<dim, spacedim>::FE_Pyramid(
+ const unsigned int degree,
+ const internal::GenericDoFsPerObject & dpos,
+ const typename FiniteElementData<dim>::Conformity conformity)
+ : dealii::FE_Poly<dim, spacedim>(
+ Simplex::ScalarPyramidPolynomial<dim>(degree),
+ FiniteElementData<dim>(dpos,
+ ReferenceCells::Pyramid,
+ 1,
+ degree,
+ conformity),
+ std::vector<bool>(
+ FiniteElementData<dim>(dpos, ReferenceCells::Pyramid, 1, degree)
+ .dofs_per_cell,
+ true),
+ std::vector<ComponentMask>(
+ FiniteElementData<dim>(dpos, ReferenceCells::Pyramid, 1, degree)
+ .dofs_per_cell,
+ std::vector<bool>(1, true)))
+ {
+ AssertDimension(dim, 3);
+
+
+ if (degree == 1)
+ {
+ this->unit_support_points.emplace_back(-1.0, -1.0, 0.0);
+ this->unit_support_points.emplace_back(+1.0, -1.0, 0.0);
+ this->unit_support_points.emplace_back(-1.0, +1.0, 0.0);
+ this->unit_support_points.emplace_back(+1.0, +1.0, 0.0);
+ this->unit_support_points.emplace_back(+0.0, +0.0, 1.0);
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ FE_PyramidP<dim, spacedim>::FE_PyramidP(const unsigned int degree)
+ : FE_Pyramid<dim, spacedim>(degree,
+ get_dpo_vector_fe_pyramid_p(degree),
+ FiniteElementData<dim>::H1)
+ {}
+
+
+
+ template <int dim, int spacedim>
+ std::unique_ptr<FiniteElement<dim, spacedim>>
+ FE_PyramidP<dim, spacedim>::clone() const
+ {
+ return std::make_unique<FE_PyramidP<dim, spacedim>>(*this);
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::string
+ FE_PyramidP<dim, spacedim>::get_name() const
+ {
+ std::ostringstream namebuf;
+ namebuf << "FE_PyramidP<" << dim << ">(" << this->degree << ")";
+
+ return namebuf.str();
+ }
+
+
+
+ template <int dim, int spacedim>
+ FiniteElementDomination::Domination
+ FE_PyramidP<dim, spacedim>::compare_for_domination(
+ const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int codim) const
+ {
+ Assert(codim <= dim, ExcImpossibleInDim(dim));
+
+ // vertex/line/face domination
+ // (if fe_other is derived from FE_DGP)
+ // ------------------------------------
+ if (codim > 0)
+ if (dynamic_cast<const FE_DGP<dim, spacedim> *>(&fe_other) != nullptr)
+ // there are no requirements between continuous and discontinuous
+ // elements
+ return FiniteElementDomination::no_requirements;
+
+ // vertex/line/face domination
+ // (if fe_other is not derived from FE_DGP)
+ // & cell domination
+ // ----------------------------------------
+ if (const FE_PyramidP<dim, spacedim> *fe_pp_other =
+ dynamic_cast<const FE_PyramidP<dim, spacedim> *>(&fe_other))
+ {
+ if (this->degree < fe_pp_other->degree)
+ return FiniteElementDomination::this_element_dominates;
+ else if (this->degree == fe_pp_other->degree)
+ return FiniteElementDomination::either_element_can_dominate;
+ else
+ return FiniteElementDomination::other_element_dominates;
+ }
+ else if (const FE_P<dim, spacedim> *fe_p_other =
+ dynamic_cast<const FE_P<dim, spacedim> *>(&fe_other))
+ {
+ if (this->degree < fe_p_other->degree)
+ return FiniteElementDomination::this_element_dominates;
+ else if (this->degree == fe_p_other->degree)
+ return FiniteElementDomination::either_element_can_dominate;
+ else
+ return FiniteElementDomination::other_element_dominates;
+ }
+ else if (const FE_Q<dim, spacedim> *fe_q_other =
+ dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
+ {
+ if (this->degree < fe_q_other->degree)
+ return FiniteElementDomination::this_element_dominates;
+ else if (this->degree == fe_q_other->degree)
+ return FiniteElementDomination::either_element_can_dominate;
+ else
+ return FiniteElementDomination::other_element_dominates;
+ }
+ else if (const FE_Nothing<dim, spacedim> *fe_nothing =
+ dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
+ {
+ if (fe_nothing->is_dominating())
+ return FiniteElementDomination::other_element_dominates;
+ else
+ // the FE_Nothing has no degrees of freedom and it is typically used
+ // in a context where we don't require any continuity along the
+ // interface
+ return FiniteElementDomination::no_requirements;
+ }
+
+ Assert(false, ExcNotImplemented());
+ return FiniteElementDomination::neither_element_dominates;
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::vector<std::pair<unsigned int, unsigned int>>
+ FE_PyramidP<dim, spacedim>::hp_vertex_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const
+ {
+ (void)fe_other;
+
+ Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)) ||
+ (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+
+ return {{0, 0}};
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::vector<std::pair<unsigned int, unsigned int>>
+ FE_PyramidP<dim, spacedim>::hp_line_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const
+ {
+ (void)fe_other;
+
+ Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)) ||
+ (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+
+ std::vector<std::pair<unsigned int, unsigned int>> result;
+
+ for (unsigned int i = 0; i < this->degree - 1; ++i)
+ result.emplace_back(i, i);
+
+ return result;
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::vector<std::pair<unsigned int, unsigned int>>
+ FE_PyramidP<dim, spacedim>::hp_quad_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int face_no) const
+ {
+ (void)fe_other;
+
+
+ AssertIndexRange(face_no, 5);
+
+ if (face_no == 0)
+ {
+ Assert((dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+ }
+ else
+ {
+ Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+ }
+
+ std::vector<std::pair<unsigned int, unsigned int>> result;
+
+ for (unsigned int i = 0; i < this->n_dofs_per_quad(face_no); ++i)
+ result.emplace_back(i, i);
+
+ return result;
+ }
+
+
+
+ template <int dim, int spacedim>
+ FE_PyramidDGP<dim, spacedim>::FE_PyramidDGP(const unsigned int degree)
+ : FE_Pyramid<dim, spacedim>(degree,
+ get_dpo_vector_fe_pyramid_dgp(degree),
+ FiniteElementData<dim>::L2)
+ {}
+
+
+
+ template <int dim, int spacedim>
+ std::unique_ptr<FiniteElement<dim, spacedim>>
+ FE_PyramidDGP<dim, spacedim>::clone() const
+ {
+ return std::make_unique<FE_PyramidDGP<dim, spacedim>>(*this);
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::string
+ FE_PyramidDGP<dim, spacedim>::get_name() const
+ {
+ std::ostringstream namebuf;
+ namebuf << "FE_PyramidDGP<" << dim << ">(" << this->degree << ")";
+
+ return namebuf.str();
+ }
+} // namespace Simplex
+
+// explicit instantiations
+#include "fe_pyramid_p.inst"
+
+DEAL_II_NAMESPACE_CLOSE
for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS)
{
#if deal_II_dimension <= deal_II_space_dimension
- template class Simplex::FE_Poly<deal_II_dimension, deal_II_space_dimension>;
- template class Simplex::FE_P<deal_II_dimension, deal_II_space_dimension>;
- template class Simplex::FE_DGP<deal_II_dimension, deal_II_space_dimension>;
- template class Simplex::FE_Wedge<deal_II_dimension,
- deal_II_space_dimension>;
- template class Simplex::FE_WedgeP<deal_II_dimension,
- deal_II_space_dimension>;
- template class Simplex::FE_WedgeDGP<deal_II_dimension,
- deal_II_space_dimension>;
template class Simplex::FE_Pyramid<deal_II_dimension,
deal_II_space_dimension>;
template class Simplex::FE_PyramidP<deal_II_dimension,
deal_II_space_dimension>;
template class Simplex::FE_PyramidDGP<deal_II_dimension,
deal_II_space_dimension>;
- template class Simplex::FE_P_Bubbles<deal_II_dimension,
- deal_II_space_dimension>;
#endif
}
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/lac/vector.h>
-#include <deal.II/simplex/fe_lib.h>
-
#include <memory>
#include <sstream>
#include <vector>
#include <deal.II/fe/fe_dgp.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q_base.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_tools.h>
-
-#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <memory>
#include <sstream>
#include <deal.II/base/config.h>
+#include <deal.II/base/polynomials_barycentric.h>
#include <deal.II/base/qprojector.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
#include <deal.II/fe/fe_tools.h>
-#include <deal.II/simplex/barycentric_polynomials.h>
-#include <deal.II/simplex/fe_lib.h>
-
DEAL_II_NAMESPACE_OPEN
namespace Simplex
return dpo;
}
-
- /**
- * Helper function to set up the dpo vector of FE_WedgeP for a given @p degree.
- */
- internal::GenericDoFsPerObject
- get_dpo_vector_fe_wedge_p(const unsigned int degree)
- {
- internal::GenericDoFsPerObject dpo;
-
- if (degree == 1)
- {
- dpo.dofs_per_object_exclusive = {{1}, {0}, {0, 0, 0, 0, 0}, {0}};
- dpo.dofs_per_object_inclusive = {{1}, {2}, {3, 3, 4, 4, 4}, {6}};
- dpo.object_index = {{}, {6}, {6}, {6}};
- dpo.first_object_index_on_face = {{},
- {3, 3, 4, 4, 4},
- {3, 3, 4, 4, 4}};
- }
- else if (degree == 2)
- {
- dpo.dofs_per_object_exclusive = {{1}, {1}, {0, 0, 1, 1, 1}, {0}};
- dpo.dofs_per_object_inclusive = {{1}, {3}, {6, 6, 9, 9, 9}, {18}};
- dpo.object_index = {{}, {6}, {15, 15, 15, 16, 17}, {18}};
- dpo.first_object_index_on_face = {{},
- {3, 3, 4, 4, 4},
- {6, 6, 8, 8, 8}};
- }
- else
- {
- Assert(false, ExcNotImplemented());
- }
-
- return dpo;
- }
-
- /**
- * Helper function to set up the dpo vector of FE_WedgeDGP for a given @p degree.
- */
- internal::GenericDoFsPerObject
- get_dpo_vector_fe_wedge_dgp(const unsigned int degree)
- {
- unsigned int n_dofs = 0;
-
- if (degree == 1)
- n_dofs = 6;
- else if (degree == 2)
- n_dofs = 18;
- else
- Assert(false, ExcNotImplemented());
-
- return internal::expand(3, {{0, 0, 0, n_dofs}}, ReferenceCells::Wedge);
- }
-
- /**
- * Helper function to set up the dpo vector of FE_PyramidP for a given @p degree.
- */
- internal::GenericDoFsPerObject
- get_dpo_vector_fe_pyramid_p(const unsigned int degree)
- {
- internal::GenericDoFsPerObject dpo;
-
- if (degree == 1)
- {
- dpo.dofs_per_object_exclusive = {{1}, {0}, {0, 0, 0, 0, 0}, {0}};
- dpo.dofs_per_object_inclusive = {{1}, {2}, {4, 3, 3, 3, 3}, {5}};
- dpo.object_index = {{}, {5}, {5}, {5}};
- dpo.first_object_index_on_face = {{},
- {4, 3, 3, 3, 3},
- {4, 3, 3, 3, 3}};
- }
- else
- {
- Assert(false, ExcNotImplemented());
- }
-
- return dpo;
- }
-
- /**
- * Helper function to set up the dpo vector of FE_PyramidDGP for a given @p degree.
- */
- internal::GenericDoFsPerObject
- get_dpo_vector_fe_pyramid_dgp(const unsigned int degree)
- {
- unsigned int n_dofs = 0;
-
- if (degree == 1)
- n_dofs = 5;
- else
- Assert(false, ExcNotImplemented());
-
- return internal::expand(3, {{0, 0, 0, n_dofs}}, ReferenceCells::Pyramid);
- }
} // namespace
return {};
}
-
-
-
- template <int dim, int spacedim>
- FE_Wedge<dim, spacedim>::FE_Wedge(
- const unsigned int degree,
- const internal::GenericDoFsPerObject & dpos,
- const typename FiniteElementData<dim>::Conformity conformity)
- : dealii::FE_Poly<dim, spacedim>(
- Simplex::ScalarWedgePolynomial<dim>(degree),
- FiniteElementData<dim>(dpos,
- ReferenceCells::Wedge,
- 1,
- degree,
- conformity),
- std::vector<bool>(
- FiniteElementData<dim>(dpos, ReferenceCells::Wedge, 1, degree)
- .dofs_per_cell,
- true),
- std::vector<ComponentMask>(
- FiniteElementData<dim>(dpos, ReferenceCells::Wedge, 1, degree)
- .dofs_per_cell,
- std::vector<bool>(1, true)))
- {
- AssertDimension(dim, 3);
-
- if (degree == 1)
- {
- this->unit_support_points.emplace_back(0.0, 0.0, 0.0);
- this->unit_support_points.emplace_back(1.0, 0.0, 0.0);
- this->unit_support_points.emplace_back(0.0, 1.0, 0.0);
- this->unit_support_points.emplace_back(0.0, 0.0, 1.0);
- this->unit_support_points.emplace_back(1.0, 0.0, 1.0);
- this->unit_support_points.emplace_back(0.0, 1.0, 1.0);
- }
- }
-
-
-
- template <int dim, int spacedim>
- FE_WedgeP<dim, spacedim>::FE_WedgeP(const unsigned int degree)
- : FE_Wedge<dim, spacedim>(degree,
- get_dpo_vector_fe_wedge_p(degree),
- FiniteElementData<dim>::H1)
- {}
-
-
-
- template <int dim, int spacedim>
- std::unique_ptr<FiniteElement<dim, spacedim>>
- FE_WedgeP<dim, spacedim>::clone() const
- {
- return std::make_unique<FE_WedgeP<dim, spacedim>>(*this);
- }
-
-
-
- template <int dim, int spacedim>
- std::string
- FE_WedgeP<dim, spacedim>::get_name() const
- {
- std::ostringstream namebuf;
- namebuf << "FE_WedgeP<" << dim << ">(" << this->degree << ")";
-
- return namebuf.str();
- }
-
-
-
- template <int dim, int spacedim>
- FiniteElementDomination::Domination
- FE_WedgeP<dim, spacedim>::compare_for_domination(
- const FiniteElement<dim, spacedim> &fe_other,
- const unsigned int codim) const
- {
- Assert(codim <= dim, ExcImpossibleInDim(dim));
-
- // vertex/line/face domination
- // (if fe_other is derived from FE_DGP)
- // ------------------------------------
- if (codim > 0)
- if (dynamic_cast<const FE_DGP<dim, spacedim> *>(&fe_other) != nullptr)
- // there are no requirements between continuous and discontinuous
- // elements
- return FiniteElementDomination::no_requirements;
-
-
- // vertex/line/face domination
- // (if fe_other is not derived from FE_DGP)
- // & cell domination
- // ----------------------------------------
- if (const FE_WedgeP<dim, spacedim> *fe_wp_other =
- dynamic_cast<const FE_WedgeP<dim, spacedim> *>(&fe_other))
- {
- if (this->degree < fe_wp_other->degree)
- return FiniteElementDomination::this_element_dominates;
- else if (this->degree == fe_wp_other->degree)
- return FiniteElementDomination::either_element_can_dominate;
- else
- return FiniteElementDomination::other_element_dominates;
- }
- else if (const FE_P<dim, spacedim> *fe_p_other =
- dynamic_cast<const FE_P<dim, spacedim> *>(&fe_other))
- {
- if (this->degree < fe_p_other->degree)
- return FiniteElementDomination::this_element_dominates;
- else if (this->degree == fe_p_other->degree)
- return FiniteElementDomination::either_element_can_dominate;
- else
- return FiniteElementDomination::other_element_dominates;
- }
- else if (const FE_Q<dim, spacedim> *fe_q_other =
- dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
- {
- if (this->degree < fe_q_other->degree)
- return FiniteElementDomination::this_element_dominates;
- else if (this->degree == fe_q_other->degree)
- return FiniteElementDomination::either_element_can_dominate;
- else
- return FiniteElementDomination::other_element_dominates;
- }
- else if (const FE_Nothing<dim> *fe_nothing =
- dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
- {
- if (fe_nothing->is_dominating())
- return FiniteElementDomination::other_element_dominates;
- else
- // the FE_Nothing has no degrees of freedom and it is typically used
- // in a context where we don't require any continuity along the
- // interface
- return FiniteElementDomination::no_requirements;
- }
-
- Assert(false, ExcNotImplemented());
- return FiniteElementDomination::neither_element_dominates;
- }
-
-
-
- template <int dim, int spacedim>
- std::vector<std::pair<unsigned int, unsigned int>>
- FE_WedgeP<dim, spacedim>::hp_vertex_dof_identities(
- const FiniteElement<dim, spacedim> &fe_other) const
- {
- (void)fe_other;
-
- Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)) ||
- (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
- ExcNotImplemented());
-
- return {{0, 0}};
- }
-
-
-
- template <int dim, int spacedim>
- std::vector<std::pair<unsigned int, unsigned int>>
- FE_WedgeP<dim, spacedim>::hp_line_dof_identities(
- const FiniteElement<dim, spacedim> &fe_other) const
- {
- (void)fe_other;
-
- Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)) ||
- (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
- ExcNotImplemented());
-
- std::vector<std::pair<unsigned int, unsigned int>> result;
-
- for (unsigned int i = 0; i < this->degree - 1; ++i)
- result.emplace_back(i, i);
-
- return result;
- }
-
-
-
- template <int dim, int spacedim>
- std::vector<std::pair<unsigned int, unsigned int>>
- FE_WedgeP<dim, spacedim>::hp_quad_dof_identities(
- const FiniteElement<dim, spacedim> &fe_other,
- const unsigned int face_no) const
- {
- (void)fe_other;
-
- AssertIndexRange(face_no, 5);
-
- if (face_no < 2)
- {
- Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)),
- ExcNotImplemented());
- }
- else
- {
- Assert((dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
- ExcNotImplemented());
- }
-
- std::vector<std::pair<unsigned int, unsigned int>> result;
-
- for (unsigned int i = 0; i < this->n_dofs_per_quad(face_no); ++i)
- result.emplace_back(i, i);
-
- return result;
- }
-
-
-
- template <int dim, int spacedim>
- FE_WedgeDGP<dim, spacedim>::FE_WedgeDGP(const unsigned int degree)
- : FE_Wedge<dim, spacedim>(degree,
- get_dpo_vector_fe_wedge_dgp(degree),
- FiniteElementData<dim>::L2)
- {}
-
-
-
- template <int dim, int spacedim>
- std::unique_ptr<FiniteElement<dim, spacedim>>
- FE_WedgeDGP<dim, spacedim>::clone() const
- {
- return std::make_unique<FE_WedgeDGP<dim, spacedim>>(*this);
- }
-
-
-
- template <int dim, int spacedim>
- std::string
- FE_WedgeDGP<dim, spacedim>::get_name() const
- {
- std::ostringstream namebuf;
- namebuf << "FE_WedgeDGP<" << dim << ">(" << this->degree << ")";
-
- return namebuf.str();
- }
-
-
-
- template <int dim, int spacedim>
- FE_Pyramid<dim, spacedim>::FE_Pyramid(
- const unsigned int degree,
- const internal::GenericDoFsPerObject & dpos,
- const typename FiniteElementData<dim>::Conformity conformity)
- : dealii::FE_Poly<dim, spacedim>(
- Simplex::ScalarPyramidPolynomial<dim>(degree),
- FiniteElementData<dim>(dpos,
- ReferenceCells::Pyramid,
- 1,
- degree,
- conformity),
- std::vector<bool>(
- FiniteElementData<dim>(dpos, ReferenceCells::Pyramid, 1, degree)
- .dofs_per_cell,
- true),
- std::vector<ComponentMask>(
- FiniteElementData<dim>(dpos, ReferenceCells::Pyramid, 1, degree)
- .dofs_per_cell,
- std::vector<bool>(1, true)))
- {
- AssertDimension(dim, 3);
-
-
- if (degree == 1)
- {
- this->unit_support_points.emplace_back(-1.0, -1.0, 0.0);
- this->unit_support_points.emplace_back(+1.0, -1.0, 0.0);
- this->unit_support_points.emplace_back(-1.0, +1.0, 0.0);
- this->unit_support_points.emplace_back(+1.0, +1.0, 0.0);
- this->unit_support_points.emplace_back(+0.0, +0.0, 1.0);
- }
- }
-
-
-
- template <int dim, int spacedim>
- FE_PyramidP<dim, spacedim>::FE_PyramidP(const unsigned int degree)
- : FE_Pyramid<dim, spacedim>(degree,
- get_dpo_vector_fe_pyramid_p(degree),
- FiniteElementData<dim>::H1)
- {}
-
-
-
- template <int dim, int spacedim>
- std::unique_ptr<FiniteElement<dim, spacedim>>
- FE_PyramidP<dim, spacedim>::clone() const
- {
- return std::make_unique<FE_PyramidP<dim, spacedim>>(*this);
- }
-
-
-
- template <int dim, int spacedim>
- std::string
- FE_PyramidP<dim, spacedim>::get_name() const
- {
- std::ostringstream namebuf;
- namebuf << "FE_PyramidP<" << dim << ">(" << this->degree << ")";
-
- return namebuf.str();
- }
-
-
-
- template <int dim, int spacedim>
- FiniteElementDomination::Domination
- FE_PyramidP<dim, spacedim>::compare_for_domination(
- const FiniteElement<dim, spacedim> &fe_other,
- const unsigned int codim) const
- {
- Assert(codim <= dim, ExcImpossibleInDim(dim));
-
- // vertex/line/face domination
- // (if fe_other is derived from FE_DGP)
- // ------------------------------------
- if (codim > 0)
- if (dynamic_cast<const FE_DGP<dim, spacedim> *>(&fe_other) != nullptr)
- // there are no requirements between continuous and discontinuous
- // elements
- return FiniteElementDomination::no_requirements;
-
- // vertex/line/face domination
- // (if fe_other is not derived from FE_DGP)
- // & cell domination
- // ----------------------------------------
- if (const FE_PyramidP<dim, spacedim> *fe_pp_other =
- dynamic_cast<const FE_PyramidP<dim, spacedim> *>(&fe_other))
- {
- if (this->degree < fe_pp_other->degree)
- return FiniteElementDomination::this_element_dominates;
- else if (this->degree == fe_pp_other->degree)
- return FiniteElementDomination::either_element_can_dominate;
- else
- return FiniteElementDomination::other_element_dominates;
- }
- else if (const FE_P<dim, spacedim> *fe_p_other =
- dynamic_cast<const FE_P<dim, spacedim> *>(&fe_other))
- {
- if (this->degree < fe_p_other->degree)
- return FiniteElementDomination::this_element_dominates;
- else if (this->degree == fe_p_other->degree)
- return FiniteElementDomination::either_element_can_dominate;
- else
- return FiniteElementDomination::other_element_dominates;
- }
- else if (const FE_Q<dim, spacedim> *fe_q_other =
- dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
- {
- if (this->degree < fe_q_other->degree)
- return FiniteElementDomination::this_element_dominates;
- else if (this->degree == fe_q_other->degree)
- return FiniteElementDomination::either_element_can_dominate;
- else
- return FiniteElementDomination::other_element_dominates;
- }
- else if (const FE_Nothing<dim, spacedim> *fe_nothing =
- dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
- {
- if (fe_nothing->is_dominating())
- return FiniteElementDomination::other_element_dominates;
- else
- // the FE_Nothing has no degrees of freedom and it is typically used
- // in a context where we don't require any continuity along the
- // interface
- return FiniteElementDomination::no_requirements;
- }
-
- Assert(false, ExcNotImplemented());
- return FiniteElementDomination::neither_element_dominates;
- }
-
-
-
- template <int dim, int spacedim>
- std::vector<std::pair<unsigned int, unsigned int>>
- FE_PyramidP<dim, spacedim>::hp_vertex_dof_identities(
- const FiniteElement<dim, spacedim> &fe_other) const
- {
- (void)fe_other;
-
- Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)) ||
- (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
- ExcNotImplemented());
-
- return {{0, 0}};
- }
-
-
-
- template <int dim, int spacedim>
- std::vector<std::pair<unsigned int, unsigned int>>
- FE_PyramidP<dim, spacedim>::hp_line_dof_identities(
- const FiniteElement<dim, spacedim> &fe_other) const
- {
- (void)fe_other;
-
- Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)) ||
- (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
- ExcNotImplemented());
-
- std::vector<std::pair<unsigned int, unsigned int>> result;
-
- for (unsigned int i = 0; i < this->degree - 1; ++i)
- result.emplace_back(i, i);
-
- return result;
- }
-
-
-
- template <int dim, int spacedim>
- std::vector<std::pair<unsigned int, unsigned int>>
- FE_PyramidP<dim, spacedim>::hp_quad_dof_identities(
- const FiniteElement<dim, spacedim> &fe_other,
- const unsigned int face_no) const
- {
- (void)fe_other;
-
-
- AssertIndexRange(face_no, 5);
-
- if (face_no == 0)
- {
- Assert((dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
- ExcNotImplemented());
- }
- else
- {
- Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)),
- ExcNotImplemented());
- }
-
- std::vector<std::pair<unsigned int, unsigned int>> result;
-
- for (unsigned int i = 0; i < this->n_dofs_per_quad(face_no); ++i)
- result.emplace_back(i, i);
-
- return result;
- }
-
-
-
- template <int dim, int spacedim>
- FE_PyramidDGP<dim, spacedim>::FE_PyramidDGP(const unsigned int degree)
- : FE_Pyramid<dim, spacedim>(degree,
- get_dpo_vector_fe_pyramid_dgp(degree),
- FiniteElementData<dim>::L2)
- {}
-
-
-
- template <int dim, int spacedim>
- std::unique_ptr<FiniteElement<dim, spacedim>>
- FE_PyramidDGP<dim, spacedim>::clone() const
- {
- return std::make_unique<FE_PyramidDGP<dim, spacedim>>(*this);
- }
-
-
-
- template <int dim, int spacedim>
- std::string
- FE_PyramidDGP<dim, spacedim>::get_name() const
- {
- std::ostringstream namebuf;
- namebuf << "FE_PyramidDGP<" << dim << ">(" << this->degree << ")";
-
- return namebuf.str();
- }
-
-
- namespace FE_P_BubblesImplementation
- {
- template <int dim>
- std::vector<unsigned int>
- get_dpo_vector(const unsigned int degree)
- {
- std::vector<unsigned int> dpo(dim + 1);
- if (degree == 0)
- {
- dpo[dim] = 1; // single interior dof
- }
- else
- {
- Assert(degree == 1 || degree == 2, ExcNotImplemented());
- dpo[0] = 1; // vertex dofs
-
- if (degree == 2)
- {
- dpo[1] = 1; // line dofs
-
- if (dim > 1)
- dpo[dim] = 1; // the internal bubble function
- if (dim == 3)
- dpo[dim - 1] = 1; // face bubble functions
- }
- }
-
- return dpo;
- }
-
-
-
- template <int dim>
- std::vector<Point<dim>>
- unit_support_points(const unsigned int degree)
- {
- Assert(degree < 3, ExcNotImplemented());
- std::vector<Point<dim>> points = unit_support_points_fe_poly<dim>(degree);
-
- Point<dim> centroid;
- std::fill(centroid.begin_raw(),
- centroid.end_raw(),
- 1.0 / double(dim + 1));
-
- switch (dim)
- {
- case 1:
- // nothing more to do
- return points;
- case 2:
- {
- if (degree == 2)
- points.push_back(centroid);
- return points;
- }
- case 3:
- {
- if (degree == 2)
- {
- const double q13 = 1.0 / 3.0;
- points.emplace_back(q13, q13, 0.0);
- points.emplace_back(q13, 0.0, q13);
- points.emplace_back(0.0, q13, q13);
- points.emplace_back(q13, q13, q13);
- points.push_back(centroid);
- }
- return points;
- }
- default:
- Assert(false, ExcNotImplemented());
- }
- return points;
- }
-
-
-
- template <int dim>
- BarycentricPolynomials<dim>
- get_basis(const unsigned int degree)
- {
- Point<dim> centroid;
- std::fill(centroid.begin_raw(),
- centroid.end_raw(),
- 1.0 / double(dim + 1));
-
- auto M = [](const unsigned int d) {
- return BarycentricPolynomial<dim, double>::monomial(d);
- };
-
- switch (degree)
- {
- // we don't need to add bubbles to P0 or P1
- case 0:
- case 1:
- return BarycentricPolynomials<dim>::get_fe_p_basis(degree);
- case 2:
- {
- const auto fe_p =
- BarycentricPolynomials<dim>::get_fe_p_basis(degree);
- // no further work is needed in 1D
- if (dim == 1)
- return fe_p;
-
- // in 2D and 3D we add a centroid bubble function
- auto c_bubble = BarycentricPolynomial<dim>() + 1;
- for (unsigned int d = 0; d < dim + 1; ++d)
- c_bubble = c_bubble * M(d);
- c_bubble = c_bubble / c_bubble.value(centroid);
-
- std::vector<BarycentricPolynomial<dim>> bubble_functions;
- if (dim == 2)
- {
- bubble_functions.push_back(c_bubble);
- }
- else if (dim == 3)
- {
- // need 'face bubble' functions in addition to the centroid.
- // Furthermore we need to subtract them off from the other
- // functions so that we end up with an interpolatory basis
- auto b0 = 27 * M(0) * M(1) * M(2);
- bubble_functions.push_back(b0 -
- b0.value(centroid) * c_bubble);
- auto b1 = 27 * M(0) * M(1) * M(3);
- bubble_functions.push_back(b1 -
- b1.value(centroid) * c_bubble);
- auto b2 = 27 * M(0) * M(2) * M(3);
- bubble_functions.push_back(b2 -
- b2.value(centroid) * c_bubble);
- auto b3 = 27 * M(1) * M(2) * M(3);
- bubble_functions.push_back(b3 -
- b3.value(centroid) * c_bubble);
-
- bubble_functions.push_back(c_bubble);
- }
-
- // Extract out the support points for the extra bubble (both
- // volume and face) functions:
- const std::vector<Point<dim>> support_points =
- unit_support_points<dim>(degree);
- const std::vector<Point<dim>> bubble_support_points(
- support_points.begin() + fe_p.n(), support_points.end());
- Assert(bubble_support_points.size() == bubble_functions.size(),
- ExcInternalError());
- const unsigned int n_bubbles = bubble_support_points.size();
-
- // Assemble the final basis:
- std::vector<BarycentricPolynomial<dim>> lump_polys;
- for (unsigned int i = 0; i < fe_p.n(); ++i)
- {
- BarycentricPolynomial<dim> p = fe_p[i];
-
- for (unsigned int j = 0; j < n_bubbles; ++j)
- {
- p = p - p.value(bubble_support_points[j]) *
- bubble_functions[j];
- }
-
- lump_polys.push_back(p);
- }
-
- for (auto &p : bubble_functions)
- lump_polys.push_back(std::move(p));
-
- // Sanity check:
-#ifdef DEBUG
- BarycentricPolynomial<dim> unity;
- for (const auto &p : lump_polys)
- unity = unity + p;
-
- Point<dim> test;
- for (unsigned int d = 0; d < dim; ++d)
- test[d] = 2.0;
- Assert(std::abs(unity.value(test) - 1.0) < 1e-10,
- ExcInternalError());
-#endif
-
- return BarycentricPolynomials<dim>(lump_polys);
- }
- default:
- Assert(degree < 3, ExcNotImplemented());
- }
-
- Assert(degree < 3, ExcNotImplemented());
- // bogus return to placate compilers
- return BarycentricPolynomials<dim>::get_fe_p_basis(degree);
- }
-
-
-
- template <int dim>
- FiniteElementData<dim>
- get_fe_data(const unsigned int degree)
- {
- // It's not efficient, but delegate computation of the degree of the
- // finite element (which is different from the input argument) to the
- // basis.
- const auto polys = get_basis<dim>(degree);
- return FiniteElementData<dim>(get_dpo_vector<dim>(degree),
- ReferenceCells::get_simplex<dim>(),
- 1, // n_components
- polys.degree(),
- FiniteElementData<dim>::H1);
- }
- } // namespace FE_P_BubblesImplementation
-
-
-
- template <int dim, int spacedim>
- FE_P_Bubbles<dim, spacedim>::FE_P_Bubbles(const unsigned int degree)
- : dealii::FE_Poly<dim, spacedim>(
- FE_P_BubblesImplementation::get_basis<dim>(degree),
- FE_P_BubblesImplementation::get_fe_data<dim>(degree),
- std::vector<bool>(
- FE_P_BubblesImplementation::get_fe_data<dim>(degree).dofs_per_cell,
- true),
- std::vector<ComponentMask>(
- FE_P_BubblesImplementation::get_fe_data<dim>(degree).dofs_per_cell,
- std::vector<bool>(1, true)))
- , approximation_degree(degree)
- {
- this->unit_support_points =
- FE_P_BubblesImplementation::unit_support_points<dim>(degree);
-
- // TODO
- // this->unit_face_support_points =
- // unit_face_support_points_fe_poly<dim>(degree);
- }
-
-
-
- template <int dim, int spacedim>
- std::string
- FE_P_Bubbles<dim, spacedim>::get_name() const
- {
- return "Simplex::FE_P_Bubbles<" + Utilities::dim_string(dim, spacedim) +
- ">" + "(" + std::to_string(approximation_degree) + ")";
- }
-
-
-
- template <int dim, int spacedim>
- void
- FE_P_Bubbles<dim, spacedim>::
- convert_generalized_support_point_values_to_dof_values(
- const std::vector<Vector<double>> &support_point_values,
- std::vector<double> & nodal_values) const
- {
- AssertDimension(support_point_values.size(),
- this->get_unit_support_points().size());
- AssertDimension(support_point_values.size(), nodal_values.size());
- AssertDimension(this->dofs_per_cell, nodal_values.size());
-
- for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
- {
- AssertDimension(support_point_values[i].size(), 1);
-
- nodal_values[i] = support_point_values[i](0);
- }
- }
-
-
-
- template <int dim, int spacedim>
- std::unique_ptr<FiniteElement<dim, spacedim>>
- FE_P_Bubbles<dim, spacedim>::clone() const
- {
- return std::make_unique<FE_P_Bubbles<dim, spacedim>>(*this);
- }
} // namespace Simplex
// explicit instantiations
-#include "fe_lib.inst"
+#include "fe_simplex_p.inst"
DEAL_II_NAMESPACE_CLOSE
-for (deal_II_dimension : DIMENSIONS)
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS)
{
- template void GridGenerator::subdivided_hyper_rectangle_with_simplices(
- Triangulation<deal_II_dimension> & tria,
- const std::vector<unsigned int> &repetitions,
- const Point<deal_II_dimension> & p1,
- const Point<deal_II_dimension> & p2,
- const bool colorize);
-
- template void GridGenerator::subdivided_hyper_cube_with_simplices(
- Triangulation<deal_II_dimension> & tria,
- const unsigned int repetitions,
- const double p1,
- const double p2,
- const bool colorize);
+#if deal_II_dimension <= deal_II_space_dimension
+ template class Simplex::FE_Poly<deal_II_dimension, deal_II_space_dimension>;
+ template class Simplex::FE_P<deal_II_dimension, deal_II_space_dimension>;
+ template class Simplex::FE_DGP<deal_II_dimension, deal_II_space_dimension>;
+#endif
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/polynomials_barycentric.h>
+#include <deal.II/base/qprojector.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_tools.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Simplex
+{
+ namespace
+ {
+ /**
+ * Set up a vector that contains the unit (reference) cell support points
+ * for FE_Poly and sufficiently similar elements.
+ */
+ template <int dim>
+ std::vector<Point<dim>>
+ unit_support_points_fe_poly_bubbles(const unsigned int degree)
+ {
+ std::vector<Point<dim>> unit_points;
+
+ // Piecewise constants are a special case: use a support point at the
+ // centroid and only the centroid
+ if (degree == 0)
+ {
+ Point<dim> centroid;
+ std::fill(centroid.begin_raw(),
+ centroid.end_raw(),
+ 1.0 / double(dim + 1));
+ unit_points.emplace_back(centroid);
+ return unit_points;
+ }
+
+ if (dim == 1)
+ {
+ // We don't really have dim = 1 support for simplex elements yet, but
+ // its convenient for populating the face array
+ Assert(degree <= 2, ExcNotImplemented());
+ if (degree >= 1)
+ {
+ unit_points.emplace_back(0.0);
+ unit_points.emplace_back(1.0);
+
+ if (degree == 2)
+ unit_points.emplace_back(0.5);
+ }
+ }
+ else if (dim == 2)
+ {
+ Assert(degree <= 2, ExcNotImplemented());
+ if (degree >= 1)
+ {
+ unit_points.emplace_back(0.0, 0.0);
+ unit_points.emplace_back(1.0, 0.0);
+ unit_points.emplace_back(0.0, 1.0);
+
+ if (degree == 2)
+ {
+ unit_points.emplace_back(0.5, 0.0);
+ unit_points.emplace_back(0.5, 0.5);
+ unit_points.emplace_back(0.0, 0.5);
+ }
+ }
+ }
+ else if (dim == 3)
+ {
+ Assert(degree <= 2, ExcNotImplemented());
+ if (degree >= 1)
+ {
+ unit_points.emplace_back(0.0, 0.0, 0.0);
+ unit_points.emplace_back(1.0, 0.0, 0.0);
+ unit_points.emplace_back(0.0, 1.0, 0.0);
+ unit_points.emplace_back(0.0, 0.0, 1.0);
+
+ if (degree == 2)
+ {
+ unit_points.emplace_back(0.5, 0.0, 0.0);
+ unit_points.emplace_back(0.5, 0.5, 0.0);
+ unit_points.emplace_back(0.0, 0.5, 0.0);
+ unit_points.emplace_back(0.0, 0.0, 0.5);
+ unit_points.emplace_back(0.5, 0.0, 0.5);
+ unit_points.emplace_back(0.0, 0.5, 0.5);
+ }
+ }
+ }
+ else
+ {
+ Assert(false, ExcNotImplemented());
+ }
+
+ return unit_points;
+ }
+ } // namespace
+
+ namespace FE_P_BubblesImplementation
+ {
+ template <int dim>
+ std::vector<unsigned int>
+ get_dpo_vector(const unsigned int degree)
+ {
+ std::vector<unsigned int> dpo(dim + 1);
+ if (degree == 0)
+ {
+ dpo[dim] = 1; // single interior dof
+ }
+ else
+ {
+ Assert(degree == 1 || degree == 2, ExcNotImplemented());
+ dpo[0] = 1; // vertex dofs
+
+ if (degree == 2)
+ {
+ dpo[1] = 1; // line dofs
+
+ if (dim > 1)
+ dpo[dim] = 1; // the internal bubble function
+ if (dim == 3)
+ dpo[dim - 1] = 1; // face bubble functions
+ }
+ }
+
+ return dpo;
+ }
+
+
+
+ template <int dim>
+ std::vector<Point<dim>>
+ unit_support_points(const unsigned int degree)
+ {
+ Assert(degree < 3, ExcNotImplemented());
+ std::vector<Point<dim>> points =
+ unit_support_points_fe_poly_bubbles<dim>(degree);
+
+ Point<dim> centroid;
+ std::fill(centroid.begin_raw(),
+ centroid.end_raw(),
+ 1.0 / double(dim + 1));
+
+ switch (dim)
+ {
+ case 1:
+ // nothing more to do
+ return points;
+ case 2:
+ {
+ if (degree == 2)
+ points.push_back(centroid);
+ return points;
+ }
+ case 3:
+ {
+ if (degree == 2)
+ {
+ const double q13 = 1.0 / 3.0;
+ points.emplace_back(q13, q13, 0.0);
+ points.emplace_back(q13, 0.0, q13);
+ points.emplace_back(0.0, q13, q13);
+ points.emplace_back(q13, q13, q13);
+ points.push_back(centroid);
+ }
+ return points;
+ }
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return points;
+ }
+
+
+
+ template <int dim>
+ BarycentricPolynomials<dim>
+ get_basis(const unsigned int degree)
+ {
+ Point<dim> centroid;
+ std::fill(centroid.begin_raw(),
+ centroid.end_raw(),
+ 1.0 / double(dim + 1));
+
+ auto M = [](const unsigned int d) {
+ return BarycentricPolynomial<dim, double>::monomial(d);
+ };
+
+ switch (degree)
+ {
+ // we don't need to add bubbles to P0 or P1
+ case 0:
+ case 1:
+ return BarycentricPolynomials<dim>::get_fe_p_basis(degree);
+ case 2:
+ {
+ const auto fe_p =
+ BarycentricPolynomials<dim>::get_fe_p_basis(degree);
+ // no further work is needed in 1D
+ if (dim == 1)
+ return fe_p;
+
+ // in 2D and 3D we add a centroid bubble function
+ auto c_bubble = BarycentricPolynomial<dim>() + 1;
+ for (unsigned int d = 0; d < dim + 1; ++d)
+ c_bubble = c_bubble * M(d);
+ c_bubble = c_bubble / c_bubble.value(centroid);
+
+ std::vector<BarycentricPolynomial<dim>> bubble_functions;
+ if (dim == 2)
+ {
+ bubble_functions.push_back(c_bubble);
+ }
+ else if (dim == 3)
+ {
+ // need 'face bubble' functions in addition to the centroid.
+ // Furthermore we need to subtract them off from the other
+ // functions so that we end up with an interpolatory basis
+ auto b0 = 27 * M(0) * M(1) * M(2);
+ bubble_functions.push_back(b0 -
+ b0.value(centroid) * c_bubble);
+ auto b1 = 27 * M(0) * M(1) * M(3);
+ bubble_functions.push_back(b1 -
+ b1.value(centroid) * c_bubble);
+ auto b2 = 27 * M(0) * M(2) * M(3);
+ bubble_functions.push_back(b2 -
+ b2.value(centroid) * c_bubble);
+ auto b3 = 27 * M(1) * M(2) * M(3);
+ bubble_functions.push_back(b3 -
+ b3.value(centroid) * c_bubble);
+
+ bubble_functions.push_back(c_bubble);
+ }
+
+ // Extract out the support points for the extra bubble (both
+ // volume and face) functions:
+ const std::vector<Point<dim>> support_points =
+ unit_support_points<dim>(degree);
+ const std::vector<Point<dim>> bubble_support_points(
+ support_points.begin() + fe_p.n(), support_points.end());
+ Assert(bubble_support_points.size() == bubble_functions.size(),
+ ExcInternalError());
+ const unsigned int n_bubbles = bubble_support_points.size();
+
+ // Assemble the final basis:
+ std::vector<BarycentricPolynomial<dim>> lump_polys;
+ for (unsigned int i = 0; i < fe_p.n(); ++i)
+ {
+ BarycentricPolynomial<dim> p = fe_p[i];
+
+ for (unsigned int j = 0; j < n_bubbles; ++j)
+ {
+ p = p - p.value(bubble_support_points[j]) *
+ bubble_functions[j];
+ }
+
+ lump_polys.push_back(p);
+ }
+
+ for (auto &p : bubble_functions)
+ lump_polys.push_back(std::move(p));
+
+ // Sanity check:
+#ifdef DEBUG
+ BarycentricPolynomial<dim> unity;
+ for (const auto &p : lump_polys)
+ unity = unity + p;
+
+ Point<dim> test;
+ for (unsigned int d = 0; d < dim; ++d)
+ test[d] = 2.0;
+ Assert(std::abs(unity.value(test) - 1.0) < 1e-10,
+ ExcInternalError());
+#endif
+
+ return BarycentricPolynomials<dim>(lump_polys);
+ }
+ default:
+ Assert(degree < 3, ExcNotImplemented());
+ }
+
+ Assert(degree < 3, ExcNotImplemented());
+ // bogus return to placate compilers
+ return BarycentricPolynomials<dim>::get_fe_p_basis(degree);
+ }
+
+
+
+ template <int dim>
+ FiniteElementData<dim>
+ get_fe_data(const unsigned int degree)
+ {
+ // It's not efficient, but delegate computation of the degree of the
+ // finite element (which is different from the input argument) to the
+ // basis.
+ const auto polys = get_basis<dim>(degree);
+ return FiniteElementData<dim>(get_dpo_vector<dim>(degree),
+ ReferenceCells::get_simplex<dim>(),
+ 1, // n_components
+ polys.degree(),
+ FiniteElementData<dim>::H1);
+ }
+ } // namespace FE_P_BubblesImplementation
+
+
+
+ template <int dim, int spacedim>
+ FE_P_Bubbles<dim, spacedim>::FE_P_Bubbles(const unsigned int degree)
+ : dealii::FE_Poly<dim, spacedim>(
+ FE_P_BubblesImplementation::get_basis<dim>(degree),
+ FE_P_BubblesImplementation::get_fe_data<dim>(degree),
+ std::vector<bool>(
+ FE_P_BubblesImplementation::get_fe_data<dim>(degree).dofs_per_cell,
+ true),
+ std::vector<ComponentMask>(
+ FE_P_BubblesImplementation::get_fe_data<dim>(degree).dofs_per_cell,
+ std::vector<bool>(1, true)))
+ , approximation_degree(degree)
+ {
+ this->unit_support_points =
+ FE_P_BubblesImplementation::unit_support_points<dim>(degree);
+
+ // TODO
+ // this->unit_face_support_points =
+ // unit_face_support_points_fe_poly<dim>(degree);
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::string
+ FE_P_Bubbles<dim, spacedim>::get_name() const
+ {
+ return "Simplex::FE_P_Bubbles<" + Utilities::dim_string(dim, spacedim) +
+ ">" + "(" + std::to_string(approximation_degree) + ")";
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ FE_P_Bubbles<dim, spacedim>::
+ convert_generalized_support_point_values_to_dof_values(
+ const std::vector<Vector<double>> &support_point_values,
+ std::vector<double> & nodal_values) const
+ {
+ AssertDimension(support_point_values.size(),
+ this->get_unit_support_points().size());
+ AssertDimension(support_point_values.size(), nodal_values.size());
+ AssertDimension(this->dofs_per_cell, nodal_values.size());
+
+ for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
+ {
+ AssertDimension(support_point_values[i].size(), 1);
+
+ nodal_values[i] = support_point_values[i](0);
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::unique_ptr<FiniteElement<dim, spacedim>>
+ FE_P_Bubbles<dim, spacedim>::clone() const
+ {
+ return std::make_unique<FE_P_Bubbles<dim, spacedim>>(*this);
+ }
+} // namespace Simplex
+
+// explicit instantiations
+#include "fe_simplex_p_bubbles.inst"
+
+DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS)
+ {
+#if deal_II_dimension <= deal_II_space_dimension
+ template class Simplex::FE_P_Bubbles<deal_II_dimension,
+ deal_II_space_dimension>;
+#endif
+ }
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/polynomials_barycentric.h>
+#include <deal.II/base/qprojector.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_wedge_p.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Simplex
+{
+ namespace
+ {
+ /**
+ * Helper function to set up the dpo vector of FE_WedgeP for a given @p degree.
+ */
+ internal::GenericDoFsPerObject
+ get_dpo_vector_fe_wedge_p(const unsigned int degree)
+ {
+ internal::GenericDoFsPerObject dpo;
+
+ if (degree == 1)
+ {
+ dpo.dofs_per_object_exclusive = {{1}, {0}, {0, 0, 0, 0, 0}, {0}};
+ dpo.dofs_per_object_inclusive = {{1}, {2}, {3, 3, 4, 4, 4}, {6}};
+ dpo.object_index = {{}, {6}, {6}, {6}};
+ dpo.first_object_index_on_face = {{},
+ {3, 3, 4, 4, 4},
+ {3, 3, 4, 4, 4}};
+ }
+ else if (degree == 2)
+ {
+ dpo.dofs_per_object_exclusive = {{1}, {1}, {0, 0, 1, 1, 1}, {0}};
+ dpo.dofs_per_object_inclusive = {{1}, {3}, {6, 6, 9, 9, 9}, {18}};
+ dpo.object_index = {{}, {6}, {15, 15, 15, 16, 17}, {18}};
+ dpo.first_object_index_on_face = {{},
+ {3, 3, 4, 4, 4},
+ {6, 6, 8, 8, 8}};
+ }
+ else
+ {
+ Assert(false, ExcNotImplemented());
+ }
+
+ return dpo;
+ }
+
+ /**
+ * Helper function to set up the dpo vector of FE_WedgeDGP for a given @p degree.
+ */
+ internal::GenericDoFsPerObject
+ get_dpo_vector_fe_wedge_dgp(const unsigned int degree)
+ {
+ unsigned int n_dofs = 0;
+
+ if (degree == 1)
+ n_dofs = 6;
+ else if (degree == 2)
+ n_dofs = 18;
+ else
+ Assert(false, ExcNotImplemented());
+
+ return internal::expand(3, {{0, 0, 0, n_dofs}}, ReferenceCells::Wedge);
+ }
+ } // namespace
+
+ template <int dim, int spacedim>
+ FE_Wedge<dim, spacedim>::FE_Wedge(
+ const unsigned int degree,
+ const internal::GenericDoFsPerObject & dpos,
+ const typename FiniteElementData<dim>::Conformity conformity)
+ : dealii::FE_Poly<dim, spacedim>(
+ Simplex::ScalarWedgePolynomial<dim>(degree),
+ FiniteElementData<dim>(dpos,
+ ReferenceCells::Wedge,
+ 1,
+ degree,
+ conformity),
+ std::vector<bool>(
+ FiniteElementData<dim>(dpos, ReferenceCells::Wedge, 1, degree)
+ .dofs_per_cell,
+ true),
+ std::vector<ComponentMask>(
+ FiniteElementData<dim>(dpos, ReferenceCells::Wedge, 1, degree)
+ .dofs_per_cell,
+ std::vector<bool>(1, true)))
+ {
+ AssertDimension(dim, 3);
+
+ if (degree == 1)
+ {
+ this->unit_support_points.emplace_back(0.0, 0.0, 0.0);
+ this->unit_support_points.emplace_back(1.0, 0.0, 0.0);
+ this->unit_support_points.emplace_back(0.0, 1.0, 0.0);
+ this->unit_support_points.emplace_back(0.0, 0.0, 1.0);
+ this->unit_support_points.emplace_back(1.0, 0.0, 1.0);
+ this->unit_support_points.emplace_back(0.0, 1.0, 1.0);
+ }
+ }
+
+
+
+ template <int dim, int spacedim>
+ FE_WedgeP<dim, spacedim>::FE_WedgeP(const unsigned int degree)
+ : FE_Wedge<dim, spacedim>(degree,
+ get_dpo_vector_fe_wedge_p(degree),
+ FiniteElementData<dim>::H1)
+ {}
+
+
+
+ template <int dim, int spacedim>
+ std::unique_ptr<FiniteElement<dim, spacedim>>
+ FE_WedgeP<dim, spacedim>::clone() const
+ {
+ return std::make_unique<FE_WedgeP<dim, spacedim>>(*this);
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::string
+ FE_WedgeP<dim, spacedim>::get_name() const
+ {
+ std::ostringstream namebuf;
+ namebuf << "FE_WedgeP<" << dim << ">(" << this->degree << ")";
+
+ return namebuf.str();
+ }
+
+
+
+ template <int dim, int spacedim>
+ FiniteElementDomination::Domination
+ FE_WedgeP<dim, spacedim>::compare_for_domination(
+ const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int codim) const
+ {
+ Assert(codim <= dim, ExcImpossibleInDim(dim));
+
+ // vertex/line/face domination
+ // (if fe_other is derived from FE_DGP)
+ // ------------------------------------
+ if (codim > 0)
+ if (dynamic_cast<const FE_DGP<dim, spacedim> *>(&fe_other) != nullptr)
+ // there are no requirements between continuous and discontinuous
+ // elements
+ return FiniteElementDomination::no_requirements;
+
+
+ // vertex/line/face domination
+ // (if fe_other is not derived from FE_DGP)
+ // & cell domination
+ // ----------------------------------------
+ if (const FE_WedgeP<dim, spacedim> *fe_wp_other =
+ dynamic_cast<const FE_WedgeP<dim, spacedim> *>(&fe_other))
+ {
+ if (this->degree < fe_wp_other->degree)
+ return FiniteElementDomination::this_element_dominates;
+ else if (this->degree == fe_wp_other->degree)
+ return FiniteElementDomination::either_element_can_dominate;
+ else
+ return FiniteElementDomination::other_element_dominates;
+ }
+ else if (const FE_P<dim, spacedim> *fe_p_other =
+ dynamic_cast<const FE_P<dim, spacedim> *>(&fe_other))
+ {
+ if (this->degree < fe_p_other->degree)
+ return FiniteElementDomination::this_element_dominates;
+ else if (this->degree == fe_p_other->degree)
+ return FiniteElementDomination::either_element_can_dominate;
+ else
+ return FiniteElementDomination::other_element_dominates;
+ }
+ else if (const FE_Q<dim, spacedim> *fe_q_other =
+ dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
+ {
+ if (this->degree < fe_q_other->degree)
+ return FiniteElementDomination::this_element_dominates;
+ else if (this->degree == fe_q_other->degree)
+ return FiniteElementDomination::either_element_can_dominate;
+ else
+ return FiniteElementDomination::other_element_dominates;
+ }
+ else if (const FE_Nothing<dim> *fe_nothing =
+ dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
+ {
+ if (fe_nothing->is_dominating())
+ return FiniteElementDomination::other_element_dominates;
+ else
+ // the FE_Nothing has no degrees of freedom and it is typically used
+ // in a context where we don't require any continuity along the
+ // interface
+ return FiniteElementDomination::no_requirements;
+ }
+
+ Assert(false, ExcNotImplemented());
+ return FiniteElementDomination::neither_element_dominates;
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::vector<std::pair<unsigned int, unsigned int>>
+ FE_WedgeP<dim, spacedim>::hp_vertex_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const
+ {
+ (void)fe_other;
+
+ Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)) ||
+ (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+
+ return {{0, 0}};
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::vector<std::pair<unsigned int, unsigned int>>
+ FE_WedgeP<dim, spacedim>::hp_line_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other) const
+ {
+ (void)fe_other;
+
+ Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)) ||
+ (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+
+ std::vector<std::pair<unsigned int, unsigned int>> result;
+
+ for (unsigned int i = 0; i < this->degree - 1; ++i)
+ result.emplace_back(i, i);
+
+ return result;
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::vector<std::pair<unsigned int, unsigned int>>
+ FE_WedgeP<dim, spacedim>::hp_quad_dof_identities(
+ const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int face_no) const
+ {
+ (void)fe_other;
+
+ AssertIndexRange(face_no, 5);
+
+ if (face_no < 2)
+ {
+ Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+ }
+ else
+ {
+ Assert((dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
+ ExcNotImplemented());
+ }
+
+ std::vector<std::pair<unsigned int, unsigned int>> result;
+
+ for (unsigned int i = 0; i < this->n_dofs_per_quad(face_no); ++i)
+ result.emplace_back(i, i);
+
+ return result;
+ }
+
+
+
+ template <int dim, int spacedim>
+ FE_WedgeDGP<dim, spacedim>::FE_WedgeDGP(const unsigned int degree)
+ : FE_Wedge<dim, spacedim>(degree,
+ get_dpo_vector_fe_wedge_dgp(degree),
+ FiniteElementData<dim>::L2)
+ {}
+
+
+
+ template <int dim, int spacedim>
+ std::unique_ptr<FiniteElement<dim, spacedim>>
+ FE_WedgeDGP<dim, spacedim>::clone() const
+ {
+ return std::make_unique<FE_WedgeDGP<dim, spacedim>>(*this);
+ }
+
+
+
+ template <int dim, int spacedim>
+ std::string
+ FE_WedgeDGP<dim, spacedim>::get_name() const
+ {
+ std::ostringstream namebuf;
+ namebuf << "FE_WedgeDGP<" << dim << ">(" << this->degree << ")";
+
+ return namebuf.str();
+ }
+} // namespace Simplex
+
+// explicit instantiations
+#include "fe_wedge_p.inst"
+
+DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS)
+ {
+#if deal_II_dimension <= deal_II_space_dimension
+ template class Simplex::FE_Wedge<deal_II_dimension,
+ deal_II_space_dimension>;
+ template class Simplex::FE_WedgeP<deal_II_dimension,
+ deal_II_space_dimension>;
+ template class Simplex::FE_WedgeDGP<deal_II_dimension,
+ deal_II_space_dimension>;
+#endif
+ }
return surface_to_volume_mapping;
}
+
+
+ template <int dim, int spacedim>
+ void
+ subdivided_hyper_rectangle_with_simplices(
+ Triangulation<dim, spacedim> & tria,
+ const std::vector<unsigned int> &repetitions,
+ const Point<dim> & p1,
+ const Point<dim> & p2,
+ const bool colorize)
+ {
+# ifndef DEAL_II_WITH_SIMPLEX_SUPPORT
+ Assert(false, ExcNeedsSimplexSupport());
+# endif
+ AssertDimension(dim, spacedim);
+
+ AssertThrow(colorize == false, ExcNotImplemented());
+
+ std::vector<Point<spacedim>> vertices;
+ std::vector<CellData<dim>> cells;
+
+ if (dim == 2)
+ {
+ // determine cell sizes
+ const Point<dim> dx((p2[0] - p1[0]) / repetitions[0],
+ (p2[1] - p1[1]) / repetitions[1]);
+
+ // create vertices
+ for (unsigned int j = 0; j <= repetitions[1]; ++j)
+ for (unsigned int i = 0; i <= repetitions[0]; ++i)
+ vertices.push_back(
+ Point<spacedim>(p1[0] + dx[0] * i, p1[1] + dx[1] * j));
+
+ // create cells
+ for (unsigned int j = 0; j < repetitions[1]; ++j)
+ for (unsigned int i = 0; i < repetitions[0]; ++i)
+ {
+ // create reference QUAD cell
+ std::array<unsigned int, 4> quad{{
+ (j + 0) * (repetitions[0] + 1) + i + 0, //
+ (j + 0) * (repetitions[0] + 1) + i + 1, //
+ (j + 1) * (repetitions[0] + 1) + i + 0, //
+ (j + 1) * (repetitions[0] + 1) + i + 1 //
+ }}; //
+
+ // TRI cell 0
+ {
+ CellData<dim> tri;
+ tri.vertices = {quad[0], quad[1], quad[2]};
+ cells.push_back(tri);
+ }
+
+ // TRI cell 1
+ {
+ CellData<dim> tri;
+ tri.vertices = {quad[3], quad[2], quad[1]};
+ cells.push_back(tri);
+ }
+ }
+ }
+ else if (dim == 3)
+ {
+ // determine cell sizes
+ const Point<dim> dx((p2[0] - p1[0]) / repetitions[0],
+ (p2[1] - p1[1]) / repetitions[1],
+ (p2[2] - p1[2]) / repetitions[2]);
+
+ // create vertices
+ for (unsigned int k = 0; k <= repetitions[2]; ++k)
+ for (unsigned int j = 0; j <= repetitions[1]; ++j)
+ for (unsigned int i = 0; i <= repetitions[0]; ++i)
+ vertices.push_back(Point<spacedim>(p1[0] + dx[0] * i,
+ p1[1] + dx[1] * j,
+ p1[2] + dx[2] * k));
+
+ // create cells
+ for (unsigned int k = 0; k < repetitions[2]; ++k)
+ for (unsigned int j = 0; j < repetitions[1]; ++j)
+ for (unsigned int i = 0; i < repetitions[0]; ++i)
+ {
+ // create reference HEX cell
+ std::array<unsigned int, 8> quad{
+ {(k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+ (j + 0) * (repetitions[0] + 1) + i + 0,
+ (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+ (j + 0) * (repetitions[0] + 1) + i + 1,
+ (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+ (j + 1) * (repetitions[0] + 1) + i + 0,
+ (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+ (j + 1) * (repetitions[0] + 1) + i + 1,
+ (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+ (j + 0) * (repetitions[0] + 1) + i + 0,
+ (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+ (j + 0) * (repetitions[0] + 1) + i + 1,
+ (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+ (j + 1) * (repetitions[0] + 1) + i + 0,
+ (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
+ (j + 1) * (repetitions[0] + 1) + i + 1}};
+
+ // TET cell 0
+ {
+ CellData<dim> cell;
+ if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
+ cell.vertices = {{quad[0], quad[1], quad[2], quad[4]}};
+ else
+ cell.vertices = {{quad[0], quad[1], quad[3], quad[5]}};
+
+ cells.push_back(cell);
+ }
+
+ // TET cell 1
+ {
+ CellData<dim> cell;
+ if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
+ cell.vertices = {{quad[2], quad[1], quad[3], quad[7]}};
+ else
+ cell.vertices = {{quad[0], quad[3], quad[2], quad[6]}};
+ cells.push_back(cell);
+ }
+
+ // TET cell 2
+ {
+ CellData<dim> cell;
+ if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
+ cell.vertices = {{quad[1], quad[4], quad[5], quad[7]}};
+ else
+ cell.vertices = {{quad[0], quad[4], quad[5], quad[6]}};
+ cells.push_back(cell);
+ }
+
+ // TET cell 3
+ {
+ CellData<dim> cell;
+ if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
+ cell.vertices = {{quad[2], quad[4], quad[7], quad[6]}};
+ else
+ cell.vertices = {{quad[3], quad[5], quad[7], quad[6]}};
+ cells.push_back(cell);
+ }
+
+ // TET cell 4
+ {
+ CellData<dim> cell;
+ if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
+ cell.vertices = {{quad[1], quad[2], quad[4], quad[7]}};
+ else
+ cell.vertices = {{quad[0], quad[3], quad[6], quad[5]}};
+ cells.push_back(cell);
+ }
+ }
+ }
+ else
+ {
+ AssertThrow(colorize == false, ExcNotImplemented());
+ }
+
+ // actually create triangulation
+ tria.create_triangulation(vertices, cells, SubCellData());
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ subdivided_hyper_cube_with_simplices(Triangulation<dim, spacedim> &tria,
+ const unsigned int repetitions,
+ const double p1,
+ const double p2,
+ const bool colorize)
+ {
+ if (dim == 2)
+ {
+ subdivided_hyper_rectangle_with_simplices(
+ tria, {{repetitions, repetitions}}, {p1, p1}, {p2, p2}, colorize);
+ }
+ else if (dim == 3)
+ {
+ subdivided_hyper_rectangle_with_simplices(
+ tria,
+ {{repetitions, repetitions, repetitions}},
+ {p1, p1, p1},
+ {p2, p2, p2},
+ colorize);
+ }
+ else
+ {
+ AssertThrow(false, ExcNotImplemented())
+ }
+ }
+
} // namespace GridGenerator
// explicit instantiations
#endif
\}
}
+
+for (deal_II_dimension : DIMENSIONS)
+ {
+ template void GridGenerator::subdivided_hyper_rectangle_with_simplices(
+ Triangulation<deal_II_dimension> & tria,
+ const std::vector<unsigned int> &repetitions,
+ const Point<deal_II_dimension> & p1,
+ const Point<deal_II_dimension> & p2,
+ const bool colorize);
+
+ template void GridGenerator::subdivided_hyper_cube_with_simplices(
+ Triangulation<deal_II_dimension> & tria,
+ const unsigned int repetitions,
+ const double p1,
+ const double p2,
+ const bool colorize);
+ }
// ---------------------------------------------------------------------
#include <deal.II/base/polynomial.h>
+#include <deal.II/base/polynomials_barycentric.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/tensor_product_polynomials.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/fe/mapping_q_generic.h>
#include <deal.II/grid/reference_cell.h>
#include <deal.II/grid/tria.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/polynomials.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include <memory>
DEAL_II_NAMESPACE_OPEN
+++ /dev/null
-## ---------------------------------------------------------------------
-##
-## Copyright (C) 2020 by the deal.II authors
-##
-## This file is part of the deal.II library.
-##
-## The deal.II library is free software; you can use it, redistribute
-## it, and/or modify it under the terms of the GNU Lesser General
-## Public License as published by the Free Software Foundation; either
-## version 2.1 of the License, or (at your option) any later version.
-## The full text of the license can be found in the file LICENSE.md at
-## the top level directory of deal.II.
-##
-## ---------------------------------------------------------------------
-
-INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
-
-SET(_unity_include_src
- fe_lib.cc
- grid_generator.cc
- barycentric_polynomials.cc
- polynomials.cc
- quadrature_lib.cc
- )
-
-SET(_separate_src
- )
-
-# include all files in the unity file
-SET(_n_includes_per_unity_file 20)
-
-SETUP_SOURCE_LIST("${_unity_include_src}"
- "${_separate_src}"
- ${_n_includes_per_unity_file}
- _src
- )
-
-SET(_inst
- fe_lib.inst.in
- grid_generator.inst.in
-)
-
-FILE(GLOB _header
- ${CMAKE_SOURCE_DIR}/include/deal.II/simplex/*.h
- )
-
-DEAL_II_ADD_LIBRARY(obj_simplex OBJECT ${_src} ${_header} ${_inst})
-EXPAND_INSTANTIATIONS(obj_simplex "${_inst}")
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2020 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-#include <deal.II/base/config.h>
-
-#include <deal.II/simplex/grid_generator.h>
-
-DEAL_II_NAMESPACE_OPEN
-
-namespace GridGenerator
-{
- template <int dim, int spacedim>
- void
- subdivided_hyper_rectangle_with_simplices(
- Triangulation<dim, spacedim> & tria,
- const std::vector<unsigned int> &repetitions,
- const Point<dim> & p1,
- const Point<dim> & p2,
- const bool colorize)
- {
-#ifndef DEAL_II_WITH_SIMPLEX_SUPPORT
- Assert(false, ExcNeedsSimplexSupport());
-#endif
- AssertDimension(dim, spacedim);
-
- AssertThrow(colorize == false, ExcNotImplemented());
-
- std::vector<Point<spacedim>> vertices;
- std::vector<CellData<dim>> cells;
-
- if (dim == 2)
- {
- // determine cell sizes
- const Point<dim> dx((p2[0] - p1[0]) / repetitions[0],
- (p2[1] - p1[1]) / repetitions[1]);
-
- // create vertices
- for (unsigned int j = 0; j <= repetitions[1]; ++j)
- for (unsigned int i = 0; i <= repetitions[0]; ++i)
- vertices.push_back(
- Point<spacedim>(p1[0] + dx[0] * i, p1[1] + dx[1] * j));
-
- // create cells
- for (unsigned int j = 0; j < repetitions[1]; ++j)
- for (unsigned int i = 0; i < repetitions[0]; ++i)
- {
- // create reference QUAD cell
- std::array<unsigned int, 4> quad{{
- (j + 0) * (repetitions[0] + 1) + i + 0, //
- (j + 0) * (repetitions[0] + 1) + i + 1, //
- (j + 1) * (repetitions[0] + 1) + i + 0, //
- (j + 1) * (repetitions[0] + 1) + i + 1 //
- }}; //
-
- // TRI cell 0
- {
- CellData<dim> tri;
- tri.vertices = {quad[0], quad[1], quad[2]};
- cells.push_back(tri);
- }
-
- // TRI cell 1
- {
- CellData<dim> tri;
- tri.vertices = {quad[3], quad[2], quad[1]};
- cells.push_back(tri);
- }
- }
- }
- else if (dim == 3)
- {
- // determine cell sizes
- const Point<dim> dx((p2[0] - p1[0]) / repetitions[0],
- (p2[1] - p1[1]) / repetitions[1],
- (p2[2] - p1[2]) / repetitions[2]);
-
- // create vertices
- for (unsigned int k = 0; k <= repetitions[2]; ++k)
- for (unsigned int j = 0; j <= repetitions[1]; ++j)
- for (unsigned int i = 0; i <= repetitions[0]; ++i)
- vertices.push_back(Point<spacedim>(p1[0] + dx[0] * i,
- p1[1] + dx[1] * j,
- p1[2] + dx[2] * k));
-
- // create cells
- for (unsigned int k = 0; k < repetitions[2]; ++k)
- for (unsigned int j = 0; j < repetitions[1]; ++j)
- for (unsigned int i = 0; i < repetitions[0]; ++i)
- {
- // create reference HEX cell
- std::array<unsigned int, 8> quad{
- {(k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
- (j + 0) * (repetitions[0] + 1) + i + 0,
- (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
- (j + 0) * (repetitions[0] + 1) + i + 1,
- (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
- (j + 1) * (repetitions[0] + 1) + i + 0,
- (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
- (j + 1) * (repetitions[0] + 1) + i + 1,
- (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
- (j + 0) * (repetitions[0] + 1) + i + 0,
- (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
- (j + 0) * (repetitions[0] + 1) + i + 1,
- (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
- (j + 1) * (repetitions[0] + 1) + i + 0,
- (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
- (j + 1) * (repetitions[0] + 1) + i + 1}};
-
- // TET cell 0
- {
- CellData<dim> cell;
- if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
- cell.vertices = {{quad[0], quad[1], quad[2], quad[4]}};
- else
- cell.vertices = {{quad[0], quad[1], quad[3], quad[5]}};
-
- cells.push_back(cell);
- }
-
- // TET cell 1
- {
- CellData<dim> cell;
- if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
- cell.vertices = {{quad[2], quad[1], quad[3], quad[7]}};
- else
- cell.vertices = {{quad[0], quad[3], quad[2], quad[6]}};
- cells.push_back(cell);
- }
-
- // TET cell 2
- {
- CellData<dim> cell;
- if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
- cell.vertices = {{quad[1], quad[4], quad[5], quad[7]}};
- else
- cell.vertices = {{quad[0], quad[4], quad[5], quad[6]}};
- cells.push_back(cell);
- }
-
- // TET cell 3
- {
- CellData<dim> cell;
- if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
- cell.vertices = {{quad[2], quad[4], quad[7], quad[6]}};
- else
- cell.vertices = {{quad[3], quad[5], quad[7], quad[6]}};
- cells.push_back(cell);
- }
-
- // TET cell 4
- {
- CellData<dim> cell;
- if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
- cell.vertices = {{quad[1], quad[2], quad[4], quad[7]}};
- else
- cell.vertices = {{quad[0], quad[3], quad[6], quad[5]}};
- cells.push_back(cell);
- }
- }
- }
- else
- {
- AssertThrow(colorize == false, ExcNotImplemented());
- }
-
- // actually create triangulation
- tria.create_triangulation(vertices, cells, SubCellData());
- }
-
- template <int dim, int spacedim>
- void
- subdivided_hyper_cube_with_simplices(Triangulation<dim, spacedim> &tria,
- const unsigned int repetitions,
- const double p1,
- const double p2,
- const bool colorize)
- {
- if (dim == 2)
- {
- subdivided_hyper_rectangle_with_simplices(
- tria, {{repetitions, repetitions}}, {p1, p1}, {p2, p2}, colorize);
- }
- else if (dim == 3)
- {
- subdivided_hyper_rectangle_with_simplices(
- tria,
- {{repetitions, repetitions, repetitions}},
- {p1, p1, p1},
- {p2, p2, p2},
- colorize);
- }
- else
- {
- AssertThrow(false, ExcNotImplemented())
- }
- }
-} // namespace GridGenerator
-
-// explicit instantiations
-#include "grid_generator.inst"
-
-DEAL_II_NAMESPACE_CLOSE
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2020 by the deal.II authors
-//
-// This file is part of the deal.II library.
-//
-// The deal.II library is free software; you can use it, redistribute
-// it, and/or modify it under the terms of the GNU Lesser General
-// Public License as published by the Free Software Foundation; either
-// version 2.1 of the License, or (at your option) any later version.
-// The full text of the license can be found in the file LICENSE.md at
-// the top level directory of deal.II.
-//
-// ---------------------------------------------------------------------
-
-#include <deal.II/base/geometry_info.h>
-#include <deal.II/base/polynomial.h>
-
-#include <deal.II/simplex/quadrature_lib.h>
-
-#include <algorithm>
-#include <cmath>
-#include <functional>
-#include <limits>
-
-
-DEAL_II_NAMESPACE_OPEN
-
-namespace Simplex
-{
- template <int dim>
- QGauss<dim>::QGauss(const unsigned int n_points_1D)
- : QSimplex<dim>(Quadrature<dim>())
- {
- // fill quadrature points and quadrature weights
- if (dim == 1)
- {
- const dealii::QGauss<dim> quad(n_points_1D);
-
- this->quadrature_points = quad.get_points();
- this->weights = quad.get_weights();
- }
- else if (dim == 2)
- {
- if (n_points_1D == 1)
- {
- const double p = 1.0 / 3.0;
- this->quadrature_points.emplace_back(p, p);
- this->weights.emplace_back(0.5);
- }
- else if (n_points_1D == 2)
- {
- const double Q23 = 2.0 / 3.0;
- const double Q16 = 1.0 / 6.0;
-
- this->quadrature_points.emplace_back(Q23, Q16);
- this->quadrature_points.emplace_back(Q16, Q23);
- this->quadrature_points.emplace_back(Q16, Q16);
- this->weights.emplace_back(Q16);
- this->weights.emplace_back(Q16);
- this->weights.emplace_back(Q16);
- }
- else if (n_points_1D == 3)
- {
- const double q12 = 0.5;
-
- // clang-format off
- this->quadrature_points.emplace_back(0.3333333333330, 0.3333333333330);
- this->quadrature_points.emplace_back(0.7974269853530, 0.1012865073230);
- this->quadrature_points.emplace_back(0.1012865073230, 0.7974269853530);
- this->quadrature_points.emplace_back(0.1012865073230, 0.1012865073230);
- this->quadrature_points.emplace_back(0.0597158717898, 0.4701420641050);
- this->quadrature_points.emplace_back(0.4701420641050, 0.0597158717898);
- this->quadrature_points.emplace_back(0.4701420641050, 0.4701420641050);
- // clang-format on
-
- this->weights.emplace_back(q12 * 0.225);
- this->weights.emplace_back(q12 * 0.125939180545);
- this->weights.emplace_back(q12 * 0.125939180545);
- this->weights.emplace_back(q12 * 0.125939180545);
- this->weights.emplace_back(q12 * 0.132394152789);
- this->weights.emplace_back(q12 * 0.132394152789);
- this->weights.emplace_back(q12 * 0.132394152789);
- }
- else if (n_points_1D == 4)
- {
- Quadrature<dim>::operator=(QWitherdenVincent<dim>(n_points_1D));
- }
- }
- else if (dim == 3)
- {
- if (n_points_1D == 1)
- {
- const double Q14 = 1.0 / 4.0;
- const double Q16 = 1.0 / 6.0;
-
- this->quadrature_points.emplace_back(Q14, Q14, Q14);
- this->weights.emplace_back(Q16);
- }
- else if (n_points_1D == 2)
- {
- const double Q124 = 1.0 / 6.0 / 4.0;
-
- const double palpha = (5.0 + 3.0 * sqrt(5.0)) / 20.0;
- const double pbeta = (5.0 - sqrt(5.0)) / 20.0;
- this->quadrature_points.emplace_back(pbeta, pbeta, pbeta);
- this->quadrature_points.emplace_back(palpha, pbeta, pbeta);
- this->quadrature_points.emplace_back(pbeta, palpha, pbeta);
- this->quadrature_points.emplace_back(pbeta, pbeta, palpha);
- this->weights.emplace_back(Q124);
- this->weights.emplace_back(Q124);
- this->weights.emplace_back(Q124);
- this->weights.emplace_back(Q124);
- }
- else if (n_points_1D == 3)
- {
- const double Q16 = 1.0 / 6.0;
-
- // clang-format off
- this->quadrature_points.emplace_back(0.5684305841968444, 0.1438564719343852, 0.1438564719343852);
- this->quadrature_points.emplace_back(0.1438564719343852, 0.1438564719343852, 0.1438564719343852);
- this->quadrature_points.emplace_back(0.1438564719343852, 0.1438564719343852, 0.5684305841968444);
- this->quadrature_points.emplace_back(0.1438564719343852, 0.5684305841968444, 0.1438564719343852);
- this->quadrature_points.emplace_back(0.0000000000000000, 0.5000000000000000, 0.5000000000000000);
- this->quadrature_points.emplace_back(0.5000000000000000, 0.0000000000000000, 0.5000000000000000);
- this->quadrature_points.emplace_back(0.5000000000000000, 0.5000000000000000, 0.0000000000000000);
- this->quadrature_points.emplace_back(0.5000000000000000, 0.0000000000000000, 0.0000000000000000);
- this->quadrature_points.emplace_back(0.0000000000000000, 0.5000000000000000, 0.0000000000000000);
- this->quadrature_points.emplace_back(0.0000000000000000, 0.0000000000000000, 0.5000000000000000);
- // clang-format on
-
- this->weights.emplace_back(0.2177650698804054 * Q16);
- this->weights.emplace_back(0.2177650698804054 * Q16);
- this->weights.emplace_back(0.2177650698804054 * Q16);
- this->weights.emplace_back(0.2177650698804054 * Q16);
- this->weights.emplace_back(0.0214899534130631 * Q16);
- this->weights.emplace_back(0.0214899534130631 * Q16);
- this->weights.emplace_back(0.0214899534130631 * Q16);
- this->weights.emplace_back(0.0214899534130631 * Q16);
- this->weights.emplace_back(0.0214899534130631 * Q16);
- this->weights.emplace_back(0.0214899534130631 * Q16);
- }
- else if (n_points_1D == 4)
- {
- Quadrature<dim>::operator=(QWitherdenVincent<dim>(n_points_1D));
- }
- }
-
- AssertDimension(this->quadrature_points.size(), this->weights.size());
- Assert(this->quadrature_points.size() > 0,
- ExcNotImplemented(
- "Simplex::QGauss is currently only implemented for "
- "n_points_1D = 1, 2, 3, and 4 while you are asking for "
- "n_points_1D = " +
- Utilities::to_string(n_points_1D)));
- }
-
- namespace
- {
- template <std::size_t b_dim>
- std::vector<std::array<double, b_dim>>
- all_permutations(const std::array<double, b_dim> &b_point)
- {
- std::vector<std::array<double, b_dim>> output;
-
- // We want all possible permutations of the barycentric coordinates.
- // The easiest way to get all of them is to sort the input first and
- // then use next_permutation to cycle through them all.
- std::array<double, b_dim> temp = b_point;
- std::sort(temp.begin(), temp.end());
- do
- {
- output.push_back(temp);
- }
- while (std::next_permutation(temp.begin(), temp.end()));
-
- return output;
- }
- } // namespace
-
-
-
- template <int dim>
- QWitherdenVincent<dim>::QWitherdenVincent(const unsigned int n_points_1D)
- : QSimplex<dim>(Quadrature<dim>())
- {
- Assert(1 <= dim && dim <= 3, ExcNotImplemented());
- // Just use Gauss in 1D: this is a high-order open rule so this is a
- // reasonable equivalent for generic programming.
- if (dim == 1)
- {
- Quadrature<dim>::operator=(dealii::QGauss<dim>(n_points_1D));
- return;
- }
-
- std::array<double, dim + 1> centroid;
- std::fill(centroid.begin(), centroid.end(), 1.0 / (dim + 1.0));
- std::vector<std::vector<std::array<double, dim + 1>>> b_point_permutations;
- std::vector<double> b_weights;
-
- // We can simplify the implementation of these quadrature rules
- // by quite a bit by exploiting symmetry - we do essentially the
- // same thing for each barycentric coordinate, so we can express
- // our quadrature rule as permutations of barycentric points
- // instead of writing things out explicitly.
-
- // Apply a Barycentric permutation where one point is different.
- auto process_point_1 = [&](const double a, const double w) {
- const double b = 1.0 - dim * a;
- std::array<double, dim + 1> b_point;
- std::fill(b_point.begin(), b_point.begin() + dim, a);
- b_point[dim] = b;
-
- b_weights.push_back(w);
- b_point_permutations.push_back(all_permutations(b_point));
- };
-
- // Apply a Barycentric permutation where two points (in 3D) are different.
- auto process_point_2 = [&](const double a, const double w) {
- Assert(dim == 3, ExcInternalError());
- const double b = (1.0 - 2.0 * a) / 2.0;
- std::array<double, dim + 1> b_point;
- std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
- b_point[dim - 1] = b;
- b_point[dim] = b;
-
- b_weights.push_back(w);
- b_point_permutations.push_back(all_permutations(b_point));
- };
-
- // Apply a Barycentric permutation where three (or four) points
- // are different (since there are two inputs).
- auto process_point_3 = [&](const double a, const double b, const double w) {
- const double c = 1.0 - (dim - 1.0) * a - b;
- std::array<double, dim + 1> b_point;
- std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
- b_point[dim - 1] = b;
- b_point[dim] = c;
-
- b_weights.push_back(w);
- b_point_permutations.push_back(all_permutations(b_point));
- };
-
- if (n_points_1D == 1)
- {
- b_point_permutations.push_back({centroid});
- b_weights.push_back(1.0);
- }
- else if (n_points_1D == 2)
- {
- // This is WV-4 in 2D and WV-3 in 3D
- if (dim == 2)
- {
- process_point_1(9.1576213509770743e-02, 1.0995174365532187e-01);
- process_point_1(4.4594849091596489e-01, 2.2338158967801147e-01);
- }
- else if (dim == 3)
- {
- process_point_1(3.281633025163817e-01, 1.362178425370874e-01);
- process_point_1(1.080472498984286e-01, 1.137821574629126e-01);
- }
- }
- else if (n_points_1D == 3)
- {
- // This is the WV-5 rule in both 2D and 3D
- if (dim == 2)
- {
- b_weights.push_back(0.225);
- b_point_permutations.push_back({centroid});
-
- process_point_1(1.0128650732345634e-01, 1.2593918054482714e-01);
- process_point_1(4.7014206410511511e-01, 1.3239415278850619e-01);
- }
- else if (dim == 3)
- {
- process_point_1(3.108859192633006e-01, 1.126879257180159e-01);
- process_point_1(9.273525031089125e-02, 7.349304311636196e-02);
-
- process_point_2(4.550370412564964e-02, 4.254602077708147e-02);
- }
- }
- else if (n_points_1D == 4)
- {
- // This is the WV-7 rule in both 2D and 3D
- if (dim == 2)
- {
- process_point_1(3.3730648554587850e-02, 1.6545050110792131e-02);
- process_point_1(4.7430969250471822e-01, 7.7086646185986069e-02);
- process_point_1(2.4157738259540357e-01, 1.2794417123015558e-01);
- process_point_3(4.7036644652595216e-02,
- 1.9868331479735168e-01,
- 5.5878732903199779e-02);
- }
- else if (dim == 3)
- {
- b_point_permutations.push_back({centroid});
- b_weights.push_back(9.548528946413085e-02);
-
- process_point_1(3.157011497782028e-01, 4.232958120996703e-02);
- process_point_2(5.048982259839635e-02, 3.189692783285758e-02);
-
- process_point_3(1.888338310260010e-01,
- 5.751716375870000e-01,
- 3.720713072833462e-02);
- process_point_3(2.126547254148314e-02,
- 8.108302410985486e-01,
- 8.110770829903342e-03);
- }
- }
- else if (n_points_1D == 5)
- {
- // This is the WV-9 rule in both 2D and 3D
- if (dim == 2)
- {
- b_point_permutations.push_back({centroid});
- b_weights.push_back(9.7135796282798836e-02);
-
- process_point_1(4.4729513394452691e-02, 2.5577675658698031e-02);
- process_point_1(4.8968251919873762e-01, 3.1334700227139071e-02);
- process_point_1(4.3708959149293664e-01, 7.7827541004774278e-02);
- process_point_1(1.8820353561903275e-01, 7.9647738927210249e-02);
-
- process_point_3(3.6838412054736258e-02,
- 2.2196298916076568e-01,
- 4.3283539377289376e-02);
- }
- else if (dim == 3)
- {
- b_point_permutations.push_back({centroid});
- b_weights.push_back(5.801054891248025e-02);
-
- process_point_1(6.198169755222693e-10, 6.431928175925639e-05);
- process_point_1(1.607745353952616e-01, 2.317333846242546e-02);
- process_point_1(3.222765218214210e-01, 2.956291233542929e-02);
- process_point_1(4.510891834541358e-02, 8.063979979616182e-03);
-
- process_point_2(1.122965460043761e-01, 3.813408010370246e-02);
-
- process_point_3(4.588714487524592e-01,
- 2.554579233041310e-03,
- 8.384422198298552e-03);
- process_point_3(3.377587068533860e-02,
- 7.183503264420745e-01,
- 1.023455935274533e-02);
- process_point_3(1.836413698099279e-01,
- 3.441591057817528e-02,
- 2.052491596798814e-02);
- }
- }
- else if (n_points_1D == 6)
- {
- // There is no WV-11 rule in 3D yet
- if (dim == 2)
- {
- b_point_permutations.push_back({centroid});
- b_weights.push_back(8.5761179732224219e-02);
-
- process_point_1(2.8485417614371900e-02, 1.0431870512894697e-02);
- process_point_1(4.9589190096589092e-01, 1.6606273054585369e-02);
- process_point_1(1.0263548271224643e-01, 3.8630759237019321e-02);
- process_point_1(4.3846592676435220e-01, 6.7316154079468296e-02);
- process_point_1(2.1021995670317828e-01, 7.0515684111716576e-02);
-
- process_point_3(7.3254276860644785e-03,
- 1.4932478865208237e-01,
- 1.0290289572953278e-02);
- process_point_3(4.6010500165429957e-02,
- 2.8958112563770588e-01,
- 4.0332476640500554e-02);
- }
- else if (dim == 3)
- {
- Assert(false, ExcNotImplemented());
- }
- }
- else
- {
- Assert(false, ExcNotImplemented());
- }
-
- Assert(b_point_permutations.size() == b_weights.size(), ExcInternalError());
- for (unsigned int permutation_n = 0; permutation_n < b_weights.size();
- ++permutation_n)
- {
- for (const std::array<double, dim + 1> &b_point :
- b_point_permutations[permutation_n])
- {
- const double volume = (dim == 2 ? 1.0 / 2.0 : 1.0 / 6.0);
- this->weights.emplace_back(volume * b_weights[permutation_n]);
- Point<dim> c_point;
- std::copy(b_point.begin(),
- b_point.begin() + dim,
- c_point.begin_raw());
- this->quadrature_points.emplace_back(c_point);
- }
- }
- }
-
-
-
- template <int dim>
- QGaussWedge<dim>::QGaussWedge(const unsigned int n_points)
- : Quadrature<dim>()
- {
- AssertDimension(dim, 3);
-
- Simplex::QGauss<2> quad_tri(n_points);
- QGauss<1> quad_line(n_points);
-
- for (unsigned int i = 0; i < quad_line.size(); ++i)
- for (unsigned int j = 0; j < quad_tri.size(); ++j)
- {
- this->quadrature_points.emplace_back(quad_tri.point(j)[0],
- quad_tri.point(j)[1],
- quad_line.point(i)[0]);
- this->weights.emplace_back(quad_tri.weight(j) * quad_line.weight(i));
- }
-
- AssertDimension(this->quadrature_points.size(), this->weights.size());
- Assert(this->quadrature_points.size() > 0,
- ExcMessage("No valid quadrature points!"));
- }
-
-
-
- template <int dim>
- QGaussPyramid<dim>::QGaussPyramid(const unsigned int n_points_1D)
- : Quadrature<dim>()
- {
- AssertDimension(dim, 3);
-
- if (n_points_1D == 1)
- {
- const double Q14 = 1.0 / 4.0;
- const double Q43 = 4.0 / 3.0;
-
- this->quadrature_points.emplace_back(0, 0, Q14);
- this->weights.emplace_back(Q43);
- }
- else if (n_points_1D == 2)
- {
- // clang-format off
- this->quadrature_points.emplace_back(-0.26318405556971, -0.26318405556971, 0.54415184401122);
- this->quadrature_points.emplace_back(-0.50661630334979, -0.50661630334979, 0.12251482265544);
- this->quadrature_points.emplace_back(-0.26318405556971, +0.26318405556971, 0.54415184401122);
- this->quadrature_points.emplace_back(-0.50661630334979, +0.50661630334979, 0.12251482265544);
- this->quadrature_points.emplace_back(+0.26318405556971, -0.26318405556971, 0.54415184401122);
- this->quadrature_points.emplace_back(+0.50661630334979, -0.50661630334979, 0.12251482265544);
- this->quadrature_points.emplace_back(+0.26318405556971, +0.26318405556971, 0.54415184401122);
- this->quadrature_points.emplace_back(+0.50661630334979, +0.50661630334979, 0.12251482265544);
- // clang-format on
-
- this->weights.emplace_back(0.10078588207983);
- this->weights.emplace_back(0.23254745125351);
- this->weights.emplace_back(0.10078588207983);
- this->weights.emplace_back(0.23254745125351);
- this->weights.emplace_back(0.10078588207983);
- this->weights.emplace_back(0.23254745125351);
- this->weights.emplace_back(0.10078588207983);
- this->weights.emplace_back(0.23254745125351);
- }
-
- AssertDimension(this->quadrature_points.size(), this->weights.size());
- Assert(this->quadrature_points.size() > 0,
- ExcMessage("No valid quadrature points!"));
- }
-
-} // namespace Simplex
-
-
-template class Simplex::QGauss<1>;
-template class Simplex::QGauss<2>;
-template class Simplex::QGauss<3>;
-template class Simplex::QGaussWedge<1>;
-template class Simplex::QGaussWedge<2>;
-template class Simplex::QGaussWedge<3>;
-template class Simplex::QGaussPyramid<1>;
-template class Simplex::QGaussPyramid<2>;
-template class Simplex::QGaussPyramid<3>;
-
-template class Simplex::QWitherdenVincent<1>;
-template class Simplex::QWitherdenVincent<2>;
-template class Simplex::QWitherdenVincent<3>;
-
-DEAL_II_NAMESPACE_CLOSE
// Test MatrixFreeTools::compute_diagonal() for a Laplace operator on a simplex
// mesh.
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_p.h>
#include <deal.II/fe/mapping_fe.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
+#include <deal.II/grid/grid_generator.h>
#include "compute_diagonal_util.h"
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
template <int dim_, typename Number>
// Test Simplex::BarycentricPolynomial and Simplex::BarycentricPolynomials.
#include <deal.II/base/point.h>
+#include <deal.II/base/polynomials_barycentric.h>
#include <deal.II/base/table.h>
-#include <deal.II/simplex/barycentric_polynomials.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/polynomials.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include "../tests.h"
// Test TriaAccessor::measure() and TriaAccessor::diameter().
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
-#include <deal.II/simplex/fe_lib.h>
-
#include "../tests.h"
template <int dim, int spacedim>
//
// ---------------------------------------------------------------------
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_tools_cache.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
// Create a simplex mesh in the unit cube. Check points distribute to each cell
// (and indirectly FETools::compute_embedding_matrices() for simplices).
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/grid/grid_generator.h>
+
#include <deal.II/lac/householder.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
using namespace dealii;
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/simplex/grid_generator.h>
-
#include <fstream>
#include <iostream>
#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-
#include "../tests.h"
using namespace dealii;
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-
#include "../tests.h"
using namespace dealii;
// Test DataOut::write_vtk() for simplex meshes.
+#include <deal.II/base/quadrature_lib.h>
+
#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
using namespace dealii;
// Test DataOut::write_vtk() for simplex meshes.
+#include <deal.II/base/quadrature_lib.h>
+
#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
#include "./simplex_grids.h"
#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-
#include "../tests.h"
using namespace dealii;
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
+
+#include <deal.II/grid/grid_generator.h>
#include "../tests.h"
// Test n_dofs_per-methods of Simplex::FE_P and Simplex::FE_DGP.
-#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include "../tests.h"
// Evaluate Simplex::FE_P and Simplex::FE_DGP at quadrature points.
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/quadrature_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include "../tests.h"
// Evaluate Simplex::FE_P_Bubbles.
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/quadrature_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include "../tests.h"
// convergence rate.
#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/numerics/vector_tools_integrate_difference.h>
#include <deal.II/numerics/vector_tools_interpolate.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
using namespace dealii;
// Test FETools::get_projection_matrix for simplices.
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/grid/grid_generator.h>
+
#include <deal.II/lac/householder.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
using namespace dealii;
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
+
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
#include <deal.II/lac/affine_constraints.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-
#include "../tests.h"
#include "hanging_nodes.h"
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
+
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
#include <deal.II/lac/affine_constraints.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-
#include "../tests.h"
#include "hanging_nodes.h"
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
+
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
#include <deal.II/lac/affine_constraints.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-
#include "../tests.h"
#include "hanging_nodes.h"
*/
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/hp/fe_collection.h>
-#include <deal.II/simplex/fe_lib.h>
-
#include <vector>
#include "../tests.h"
*/
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/hp/fe_collection.h>
-#include <deal.II/simplex/fe_lib.h>
-
#include <vector>
#include "../tests.h"
// verify hanging node constraints on locally hp-refined hybrid meshes
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/hp/fe_collection.h>
-#include <deal.II/simplex/fe_lib.h>
-
#include <vector>
#include "../tests.h"
// Distribute Simplex::FE_Wedge on a DoFHandler.
+#include <deal.II/base/quadrature_lib.h>
+
#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/numerics/data_out.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
using namespace dealii;
// Test MappingFEField and VectorTools::get_position_vector() for simplex
// meshes.
+#include <deal.II/base/quadrature_lib.h>
+
#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/fe/mapping_fe_field.h>
+#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
using namespace dealii;
// Like mapping_fe_fields_01 but for deformed meshes tested for linear and
// quadratic mapping.
+#include <deal.II/base/quadrature_lib.h>
+
#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/fe/mapping_fe_field.h>
+#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
using namespace dealii;
#include <deal.II/base/patterns.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/simplex/fe_lib.h>
-
#include <fstream>
#include "../tests.h"
// continuous elements and compare results between matrix-free and matrix-based
// implementations.
+#include <deal.II/base/quadrature_lib.h>
+
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
#include "./simplex_grids.h"
// Same as matrix_free_01 but testing mixed meshes (and also pure simplex and
// hypercube mesh as special case of mixed meshs).
+#include <deal.II/base/quadrature_lib.h>
+
#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
#include "./simplex_grids.h"
// continuous elements and compare results between matrix-free and matrix-based
// implementations.
+#include <deal.II/base/quadrature_lib.h>
+
#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
using namespace dealii;
// Solve Poisson problem problem on a mixed mesh with DG and MatrixFree.
+#include <deal.II/base/quadrature_lib.h>
+
#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
#include "./simplex_grids.h"
// Test ShapeData for Simplex::FE_P and Simplex::QGauss
+#include <deal.II/base/quadrature_lib.h>
+
#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
using namespace dealii;
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/mpi.h>
+#include <deal.II/base/quadrature_lib.h>
#include <deal.II/distributed/fully_distributed_tria.h>
#include <deal.II/distributed/shared_tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
#include "simplex_grids.h"
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include <fstream>
#include <iostream>
// quadrature rule.
+#include <deal.II/base/polynomials_barycentric.h>
#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/simplex/polynomials.h>
-
#include "../tests.h"
using namespace dealii;
#include <deal.II/base/qprojector.h>
-
-#include <deal.II/simplex/quadrature_lib.h>
+#include <deal.II/base/quadrature_lib.h>
#include "../tests.h"
// ---------------------------------------------------------------------
#include <deal.II/base/function_lib.h>
-
-#include <deal.II/simplex/quadrature_lib.h>
+#include <deal.II/base/quadrature_lib.h>
#include "../tests.h"
// Test Simplex::QGauss: output its quadrature points and weights.
-#include <deal.II/simplex/quadrature_lib.h>
+#include <deal.II/base/quadrature_lib.h>
#include "../tests.h"
// Test ReferenceCell::Kind::faces_for_given_vertex().
-#include <deal.II/grid/reference_cell.h>
+#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/simplex/quadrature_lib.h>
+#include <deal.II/grid/reference_cell.h>
#include "../tests.h"
// Test refinement of simplex mesh.
+#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
-#include <deal.II/simplex/grid_generator.h>
-
#include "../tests.h"
#include "./simplex_grids.h"
// Test the correct setup of neighbors during refinement of simplex mesh.
+#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
-#include <deal.II/simplex/grid_generator.h>
-
#include "../tests.h"
#include "./simplex_grids.h"
// Verify CellAccessor::neighbor_child_on_subface for a triangle mesh.
+#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
-#include <deal.II/simplex/grid_generator.h>
-
#include "../tests.h"
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/simplex/grid_generator.h>
-
#include <cmath>
#include <fstream>
#include <iostream>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
+
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/simplex/fe_lib.h>
-
#include <fstream>
#include "../tests.h"
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include <fstream>
#include <iostream>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include <fstream>
#include <iostream>
#ifdef USE_SIMPLEX
+# include <deal.II/base/quadrature_lib.h>
# include <deal.II/base/types.h>
-# include <deal.II/simplex/fe_lib.h>
-# include <deal.II/simplex/quadrature_lib.h>
+# include <deal.II/fe/fe_pyramid_p.h>
+# include <deal.II/fe/fe_simplex_p.h>
+# include <deal.II/fe/fe_simplex_p_bubbles.h>
+# include <deal.II/fe/fe_wedge_p.h>
#else
# include <deal.II/base/quadrature_lib.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include <array>
#include <fstream>
#include <iostream>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include <fstream>
#include <iostream>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/numerics/derivative_approximation.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include <fstream>
#include <iostream>
// Finally, the new include file for using the mesh_loop from the MeshWorker
// framework
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/grid/grid_generator.h>
+
#include <deal.II/lac/precondition.h>
#include <deal.II/meshworker/mesh_loop.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include <fstream>
#include <iostream>
#include "../tests.h"
// simplex
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
+#include <deal.II/grid/grid_generator.h>
//#define HEX
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include <fstream>
#include <iostream>
#include "../tests.h"
// simplex
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
+#include <deal.II/grid/grid_generator.h>
// #define HEX
#include <deal.II/numerics/vector_tools.h>
// simplex
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
+
+#include <deal.II/grid/grid_generator.h>
#include <fstream>
#include <iostream>
#ifdef HEX
# include <deal.II/distributed/tria.h>
#else
+# include <deal.II/base/quadrature_lib.h>
+
# include <deal.II/distributed/fully_distributed_tria.h>
+# include <deal.II/fe/fe_pyramid_p.h>
+# include <deal.II/fe/fe_simplex_p.h>
+# include <deal.II/fe/fe_simplex_p_bubbles.h>
+# include <deal.II/fe/fe_wedge_p.h>
# include <deal.II/fe/mapping_fe.h>
-
-# include <deal.II/simplex/fe_lib.h>
-# include <deal.II/simplex/quadrature_lib.h>
#endif
#include <deal.II/base/discrete_time.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/parameter_acceptor.h>
+#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/timer.h>
#include <deal.II/distributed/cell_weights.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/particles/generators.h>
#include <deal.II/particles/particle_handler.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include <cmath>
#include <iostream>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/numerics/derivative_approximation.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
#define USE_SIMPLEX
// for all reference-cell types.
-#include <deal.II/grid/reference_cell.h>
+#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/simplex/quadrature_lib.h>
+#include <deal.II/grid/reference_cell.h>
#include "../tests.h"
#include <deal.II/base/qprojector.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/grid/grid_generator.h>
+
#include <deal.II/hp/q_collection.h>
#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/simplex/fe_lib.h>
-#include <deal.II/simplex/grid_generator.h>
-#include <deal.II/simplex/quadrature_lib.h>
-
#include "../tests.h"
using namespace dealii;
#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
+
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria.h>
-#include <deal.II/simplex/fe_lib.h>
-
#include "../tests.h"
using namespace dealii;
projects / dealii.git / commitdiff