]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove simplex folders 11744/head
authorPeter Munch <peterrmuench@gmail.com>
Sat, 13 Feb 2021 20:50:35 +0000 (21:50 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Wed, 24 Feb 2021 10:56:34 +0000 (11:56 +0100)
103 files changed:
include/deal.II/base/polynomials_barycentric.h [moved from include/deal.II/simplex/barycentric_polynomials.h with 100% similarity]
include/deal.II/base/polynomials_pyramid.h [new file with mode: 0644]
include/deal.II/base/polynomials_wedge.h [moved from include/deal.II/simplex/polynomials.h with 54% similarity]
include/deal.II/base/quadrature_lib.h
include/deal.II/fe/fe_pyramid_p.h [new file with mode: 0644]
include/deal.II/fe/fe_simplex_p.h [new file with mode: 0644]
include/deal.II/fe/fe_simplex_p_bubbles.h [new file with mode: 0644]
include/deal.II/fe/fe_wedge_p.h [new file with mode: 0644]
include/deal.II/grid/grid_generator.h
include/deal.II/matrix_free/matrix_free.templates.h
include/deal.II/matrix_free/shape_info.templates.h
include/deal.II/matrix_free/util.h
include/deal.II/numerics/data_out_dof_data.templates.h
include/deal.II/simplex/fe_lib.h [deleted file]
include/deal.II/simplex/grid_generator.h [deleted file]
include/deal.II/simplex/quadrature_lib.h [deleted file]
source/CMakeLists.txt
source/base/CMakeLists.txt
source/base/polynomials_barycentric.cc [moved from source/simplex/barycentric_polynomials.cc with 99% similarity]
source/base/polynomials_pyramid.cc [moved from source/simplex/polynomials.cc with 53% similarity]
source/base/polynomials_wedge.cc [new file with mode: 0644]
source/base/qprojector.cc
source/base/quadrature_lib.cc
source/fe/CMakeLists.txt
source/fe/fe_dgq.cc
source/fe/fe_pyramid_p.cc [new file with mode: 0644]
source/fe/fe_pyramid_p.inst.in [moved from source/simplex/fe_lib.inst.in with 61% similarity]
source/fe/fe_q.cc
source/fe/fe_q_base.cc
source/fe/fe_simplex_p.cc [moved from source/simplex/fe_lib.cc with 52% similarity]
source/fe/fe_simplex_p.inst.in [moved from source/simplex/grid_generator.inst.in with 51% similarity]
source/fe/fe_simplex_p_bubbles.cc [new file with mode: 0644]
source/fe/fe_simplex_p_bubbles.inst.in [new file with mode: 0644]
source/fe/fe_wedge_p.cc [new file with mode: 0644]
source/fe/fe_wedge_p.inst.in [new file with mode: 0644]
source/grid/grid_generator.cc
source/grid/grid_generator.inst.in
source/grid/reference_cell.cc
source/simplex/CMakeLists.txt [deleted file]
source/simplex/grid_generator.cc [deleted file]
source/simplex/quadrature_lib.cc [deleted file]
tests/matrix_free/compute_diagonal_05.cc
tests/multigrid-global-coarsening/multigrid_util.h
tests/simplex/barycentric_01.cc
tests/simplex/cell_measure_01.cc
tests/simplex/compute_point_locations_01.cc
tests/simplex/compute_projection_matrices_01.cc
tests/simplex/conv_hex_to_simplex.cc
tests/simplex/data_out_write_hdf5_01.cc
tests/simplex/data_out_write_hdf5_02.cc
tests/simplex/data_out_write_vtk_01.cc
tests/simplex/data_out_write_vtk_02.cc
tests/simplex/data_out_write_vtu_01.cc
tests/simplex/extract_boundary_dofs.cc
tests/simplex/fe_lib_01.cc
tests/simplex/fe_lib_02.cc
tests/simplex/fe_p_bubbles_01.cc
tests/simplex/fe_p_bubbles_02.cc
tests/simplex/get_projection_matrix_01.cc
tests/simplex/hanging_nodes_01.cc
tests/simplex/hanging_nodes_02.cc
tests/simplex/hanging_nodes_03.cc
tests/simplex/hanging_nodes_hybrid_01.cc
tests/simplex/hanging_nodes_hybrid_02.cc
tests/simplex/hanging_nodes_hybrid_03.cc
tests/simplex/mapping_fe_01.cc
tests/simplex/mapping_fe_fields_01.cc
tests/simplex/mapping_fe_fields_02.cc
tests/simplex/mapping_transformations_01.cc
tests/simplex/matrix_free_01.cc
tests/simplex/matrix_free_02.cc
tests/simplex/matrix_free_03.cc
tests/simplex/matrix_free_04.cc
tests/simplex/matrix_free_shape_info_01.cc
tests/simplex/poisson_01.cc
tests/simplex/poisson_02.cc
tests/simplex/polynomials_01.cc
tests/simplex/q_projection_01.cc
tests/simplex/q_witherden_vincent_01.cc
tests/simplex/quadrature_lib_01.cc
tests/simplex/reference_cell_kind_01.cc
tests/simplex/refinement_01.cc
tests/simplex/refinement_02.cc
tests/simplex/refinement_03.cc
tests/simplex/step-01.cc
tests/simplex/step-02.cc
tests/simplex/step-03.cc
tests/simplex/step-04.cc
tests/simplex/step-06.cc
tests/simplex/step-07.cc
tests/simplex/step-08.cc
tests/simplex/step-12.cc
tests/simplex/step-12a.cc
tests/simplex/step-18.cc
tests/simplex/step-20.cc
tests/simplex/step-23.cc
tests/simplex/step-38.cc
tests/simplex/step-55.cc
tests/simplex/step-68.cc
tests/simplex/step-74.cc
tests/simplex/unit_tangential_vectors_01.cc
tests/simplex/variable_face_quadratures_03.cc
tests/simplex/wedge_01.cc

diff --git a/include/deal.II/base/polynomials_pyramid.h b/include/deal.II/base/polynomials_pyramid.h
new file mode 100644 (file)
index 0000000..31421e9
--- /dev/null
@@ -0,0 +1,151 @@
+// ---------------------------------------------------------------------
+//
+// 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
similarity index 54%
rename from include/deal.II/simplex/polynomials.h
rename to include/deal.II/base/polynomials_wedge.h
index 558e1f97a4e18635850dad9ab63961d9b7bae243..880660f27dad2dc369eaf6de128b88422d43437d 100644 (file)
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// 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
 
 /**
@@ -72,9 +71,6 @@ namespace Simplex
              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;
 
@@ -87,9 +83,6 @@ namespace Simplex
     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;
@@ -137,15 +130,9 @@ namespace Simplex
     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;
 
@@ -163,117 +150,6 @@ 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;
-
-    /**
-     * @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>
@@ -291,25 +167,6 @@ namespace Simplex
     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
index 3a832504d6ffe363d90ee3443e4e1fd01df30ce5..890440acd47a1ed989daa1f65e2cd16d4597bf68 100644 (file)
@@ -792,6 +792,98 @@ public:
   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 ------------- */
diff --git a/include/deal.II/fe/fe_pyramid_p.h b/include/deal.II/fe/fe_pyramid_p.h
new file mode 100644 (file)
index 0000000..0079600
--- /dev/null
@@ -0,0 +1,141 @@
+// ---------------------------------------------------------------------
+//
+// 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
diff --git a/include/deal.II/fe/fe_simplex_p.h b/include/deal.II/fe/fe_simplex_p.h
new file mode 100644 (file)
index 0000000..8105d22
--- /dev/null
@@ -0,0 +1,212 @@
+// ---------------------------------------------------------------------
+//
+// 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
diff --git a/include/deal.II/fe/fe_simplex_p_bubbles.h b/include/deal.II/fe/fe_simplex_p_bubbles.h
new file mode 100644 (file)
index 0000000..d6cdb8a
--- /dev/null
@@ -0,0 +1,119 @@
+// ---------------------------------------------------------------------
+//
+// 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
diff --git a/include/deal.II/fe/fe_wedge_p.h b/include/deal.II/fe/fe_wedge_p.h
new file mode 100644 (file)
index 0000000..3a01b2e
--- /dev/null
@@ -0,0 +1,141 @@
+// ---------------------------------------------------------------------
+//
+// 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
index f301aac359af0402ca47a904397086a9c2be12cd..aaf94235add5d915c144d95ec93aaf48f4e02fe9 100644 (file)
@@ -2146,8 +2146,6 @@ namespace GridGenerator
       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)
@@ -2207,6 +2205,50 @@ namespace GridGenerator
 
   } // 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);
+
   ///@}
 
   /**
index df3f382117b05f57569a63202c37b6bab213470f..1e0d715dd4cf68c29ca7e632454fd79e792e7ecd 100644 (file)
@@ -23,6 +23,7 @@
 #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>
 
@@ -41,8 +42,6 @@
 #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>
 
index 4298f7a9bae05cb4cf5638e63e930f23996ce61e..6cf1cb2125e950cd2e9d920d33202c82a63ddd8f 100644 (file)
 #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>
 
@@ -38,8 +42,6 @@
 #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
 
index 495c6a6e6dc7840c556ec697f6e0bb6ab62f33b8..dec9542a16e9bf243385c0683f5efa76806a564a 100644 (file)
 #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
 
 
index bf8962e28220f8fe9553f12aee75fc8b926656f9..4222b7f11218fc87ed76b66c8f03b4cf6d15a861 100644 (file)
@@ -30,7 +30,9 @@
 #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>
@@ -46,8 +48,6 @@
 #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>
diff --git a/include/deal.II/simplex/fe_lib.h b/include/deal.II/simplex/fe_lib.h
deleted file mode 100644 (file)
index 7cc7800..0000000
+++ /dev/null
@@ -1,517 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// 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
diff --git a/include/deal.II/simplex/grid_generator.h b/include/deal.II/simplex/grid_generator.h
deleted file mode 100644 (file)
index 3eaba78..0000000
+++ /dev/null
@@ -1,80 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// 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
diff --git a/include/deal.II/simplex/quadrature_lib.h b/include/deal.II/simplex/quadrature_lib.h
deleted file mode 100644 (file)
index ba55488..0000000
+++ /dev/null
@@ -1,118 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// 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
index a6d85d8b68025eb3a16929cf0e20b7eb5dab8c50..7173d25b85b730b911fa2218b57e3ca4da96e513 100644 (file)
@@ -71,7 +71,6 @@ ADD_SUBDIRECTORY(differentiation)
 ADD_SUBDIRECTORY(physics)
 ADD_SUBDIRECTORY(optimization/rol)
 ADD_SUBDIRECTORY(non_matching)
-ADD_SUBDIRECTORY(simplex)
 ADD_SUBDIRECTORY(sundials)
 
 FOREACH(build ${DEAL_II_BUILD_TYPES})
index 2a924e7fb68479277e52afa56d91aaf0e76fc5d3..c48dce79cc3eb44f03b06ea9bff667174df52a67 100644 (file)
@@ -64,6 +64,7 @@ SET(_unity_include_src
   polynomials_abf.cc
   polynomials_adini.cc
   tensor_polynomials_base.cc
+  polynomials_barycentric.cc
   polynomials_bernstein.cc
   polynomials_bdm.cc
   polynomials_bernardi_raugel.cc
@@ -72,9 +73,11 @@ SET(_unity_include_src
   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
similarity index 99%
rename from source/simplex/barycentric_polynomials.cc
rename to source/base/polynomials_barycentric.cc
index 8e4e5416a5c1664b8b3f8231940ebd5e60c77350..d916f9e9441a0426b7f6cf5cb63fd68eed9324d9 100644 (file)
@@ -14,7 +14,7 @@
 // ---------------------------------------------------------------------
 
 
-#include <deal.II/simplex/barycentric_polynomials.h>
+#include <deal.II/base/polynomials_barycentric.h>
 
 DEAL_II_NAMESPACE_OPEN
 
similarity index 53%
rename from source/simplex/polynomials.cc
rename to source/base/polynomials_pyramid.cc
index 4808d455849c29d1c58f7e6e3d758c9c78468809..9321dfb3a96e0e254c53ed6db411fb168946b9c4 100644 (file)
 // ---------------------------------------------------------------------
 
 
-#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
 
@@ -39,230 +37,10 @@ namespace Simplex
 
       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)
@@ -503,9 +281,6 @@ namespace Simplex
 
 
 
-  template class ScalarWedgePolynomial<1>;
-  template class ScalarWedgePolynomial<2>;
-  template class ScalarWedgePolynomial<3>;
   template class ScalarPyramidPolynomial<1>;
   template class ScalarPyramidPolynomial<2>;
   template class ScalarPyramidPolynomial<3>;
diff --git a/source/base/polynomials_wedge.cc b/source/base/polynomials_wedge.cc
new file mode 100644 (file)
index 0000000..8ac20bb
--- /dev/null
@@ -0,0 +1,255 @@
+// ---------------------------------------------------------------------
+//
+// 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
index 0a310707ce00a4b389b7a3f9bd37b2f95bf884f4..f64d67cc5a3f60689475cf21089d1360a9784e43 100644 (file)
 
 #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
 
 
index 9be07c25bda21ce0c24e03f9a850f1469f0564f7..f84f13fb6cd0f0081c959af46b38cdff12aa52b7 100644 (file)
@@ -1350,6 +1350,449 @@ QSplit<dim>::QSplit(const QSimplex<dim> &base, const Point<dim> &split_point)
 
 
 
+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>;
@@ -1396,4 +1839,18 @@ template class QSplit<1>;
 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
index 666999b0a5b0985de05fe1526267904af725f143..9d9fc097033c6e47b7ae2ad985986a6f1c78518e 100644 (file)
@@ -29,12 +29,14 @@ SET(_unity_include_src
   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
@@ -50,8 +52,10 @@ SET(_unity_include_src
   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
@@ -105,6 +109,7 @@ SET(_inst
   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
@@ -117,6 +122,8 @@ SET(_inst
   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
@@ -128,6 +135,7 @@ SET(_inst
   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
index f75719806e9d030b654bb3c2cca3b6bff1f9a2c0..73acc4e8d59537c5a3fcc4dbc4c0da86adff742b 100644 (file)
 #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>
diff --git a/source/fe/fe_pyramid_p.cc b/source/fe/fe_pyramid_p.cc
new file mode 100644 (file)
index 0000000..fe286a5
--- /dev/null
@@ -0,0 +1,315 @@
+// ---------------------------------------------------------------------
+//
+// 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
similarity index 61%
rename from source/simplex/fe_lib.inst.in
rename to source/fe/fe_pyramid_p.inst.in
index f40343af26fe67244b8371319e5210e533cbf3b4..cce82dcff2be484d81adbe8d713466e2885229ac 100644 (file)
 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
   }
index 9be753fa616c0fbc92154312b93195436a13be37..39db78bb5ceb94c09ad7c91b14a0697f5b6cc353 100644 (file)
 
 #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>
index 64eff9f42951ac0c37cb25a8ec954cfeade58cae..c663f75efe469539a5994f05dfbc10b58d8638ff 100644 (file)
 #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>
similarity index 52%
rename from source/simplex/fe_lib.cc
rename to source/fe/fe_simplex_p.cc
index fd51f0402436e1f7e4236ce3c235d83da22c8dd3..9adbe33f9e320ccd54d01c915fe192a481651602 100644 (file)
 
 #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
@@ -260,99 +259,6 @@ 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
 
 
@@ -949,747 +855,9 @@ namespace Simplex
 
     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
similarity index 51%
rename from source/simplex/grid_generator.inst.in
rename to source/fe/fe_simplex_p.inst.in
index b65ed758bbbf2930258b53a90170cc8e58890bad..0f29f676d1957b033b34456da9d360e49d04b65c 100644 (file)
 
 
 
-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
   }
diff --git a/source/fe/fe_simplex_p_bubbles.cc b/source/fe/fe_simplex_p_bubbles.cc
new file mode 100644 (file)
index 0000000..8c2ee43
--- /dev/null
@@ -0,0 +1,390 @@
+// ---------------------------------------------------------------------
+//
+// 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
diff --git a/source/fe/fe_simplex_p_bubbles.inst.in b/source/fe/fe_simplex_p_bubbles.inst.in
new file mode 100644 (file)
index 0000000..ca82953
--- /dev/null
@@ -0,0 +1,24 @@
+// ---------------------------------------------------------------------
+//
+// 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
+  }
diff --git a/source/fe/fe_wedge_p.cc b/source/fe/fe_wedge_p.cc
new file mode 100644 (file)
index 0000000..1a1d1fd
--- /dev/null
@@ -0,0 +1,325 @@
+// ---------------------------------------------------------------------
+//
+// 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
diff --git a/source/fe/fe_wedge_p.inst.in b/source/fe/fe_wedge_p.inst.in
new file mode 100644 (file)
index 0000000..9b6ee61
--- /dev/null
@@ -0,0 +1,28 @@
+// ---------------------------------------------------------------------
+//
+// 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
+  }
index ce59fc0d2c3f117a43f95653c025264b7808eb92..c845f2d454974ef59c35aa82e673c3925e4f8b7e 100644 (file)
@@ -8220,6 +8220,196 @@ namespace GridGenerator
     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
index 17b7d02ad0204610783f2e974dcf2d2dfea24621..3f063e77d091d963f082eb7da6018725a5f98c94 100644 (file)
@@ -281,3 +281,20 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
 #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);
+  }
index 41b656fc6f854485fc0a19055c7eef4d927675f4..a9bc30f535abd4009acdd1ced876a36d563bb1b6 100644 (file)
 // ---------------------------------------------------------------------
 
 #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
diff --git a/source/simplex/CMakeLists.txt b/source/simplex/CMakeLists.txt
deleted file mode 100644 (file)
index 6814164..0000000
+++ /dev/null
@@ -1,48 +0,0 @@
-## ---------------------------------------------------------------------
-##
-## 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}")
diff --git a/source/simplex/grid_generator.cc b/source/simplex/grid_generator.cc
deleted file mode 100644 (file)
index aaed037..0000000
+++ /dev/null
@@ -1,214 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// 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
diff --git a/source/simplex/quadrature_lib.cc b/source/simplex/quadrature_lib.cc
deleted file mode 100644 (file)
index cbbd1c9..0000000
+++ /dev/null
@@ -1,485 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// 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
index 79a5ec2fbe539277bdcb67ea0006212e9d1948a7..5daa7c3d059db5b67b1cf82aac6912b3647cd830 100644 (file)
 // 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"
 
index d52957f3bba237072ebd2cb3a9f6fe9e99285d01..bc59afd3cf33f262999654b1bcb7422badc60932 100644 (file)
@@ -27,6 +27,7 @@
 
 #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>
index 6bd73fbff96b6a36de57cc4e2b2e2a9a992c2028..277d7220be0cf42d7ef261d0505514ca15a059b1 100644 (file)
 // 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"
 
index a21e5837c8312d57ca884ace8aaa17e6fe117c1f..12488df21f48a61b54252d59ffcacd9b213719ad 100644 (file)
 
 // 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>
index 66dc56cb39e37b50b7c04af65f6f1cccf1debc31..7f01defc34718f9f9f842dcc5489fceff7debcf5 100644 (file)
 //
 // ---------------------------------------------------------------------
 
+#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
index 02d85f17b5d4607b76859d4cd92ee365260126ca..0eb96bc768582a392117bc2a3c9642ec78f159c4 100644 (file)
 // (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;
index 10dab721d22fffbe24b35561b8bf4b9448cb5365..975c8a8e55073562e185c7302afbe7b7595e81f4 100644 (file)
@@ -27,8 +27,6 @@
 #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>
 
index 2c9b1b93c71e5d082e9f78b81500c17ef958619f..bcae95fde259650c8d32f924af4513876559c278 100644 (file)
 
 #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>
@@ -31,9 +36,6 @@
 #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;
index 9b1735310663cef278f51923a96f218655bd0947..038c15220bb609fb4250fe1214f31416f0145bdb 100644 (file)
 #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>
@@ -34,9 +39,6 @@
 #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;
index 0886b0fc426c7983e414a4cda9ebce92b423335c..b2d6b364804c7875f13315b491466a01150bbd54 100644 (file)
 
 // 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;
index 2e9ddbd96ce8edd2482c599b825df9a6d8825ac6..d88279770dcb43a7a4a4b8372cdf05b98c70d839 100644 (file)
 
 // 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"
index aa197b2b209d11166c8b6d676e1b2ffd06ef297a..a58cbf23072bf3b780b3db225edf5d253a21aba6 100644 (file)
 
 #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>
@@ -29,9 +34,6 @@
 #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;
index 711f9445248000d70583a1fe3387e2b1a920789a..9329aca91da7c62b900946e8add0849252d26552 100644 (file)
 #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"
 
index 9c46e8cbacb572e9fed4ddb620e902a9c38ffd07..50460a2f09f75e4ffccbe3e7bfbc74c11bec0d80 100644 (file)
 // 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"
 
index f04b0bd6603c12653520695b093734f2c7197ea9..bb96eef5f3110c29f7f1c7b339f022d2c3eb98e8 100644 (file)
 // 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"
 
index 4c90162db741111347721c172222e12c25e2238d..c79de2e615557642af3846a00dd7fdd8d2c6de6d 100644 (file)
 // 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"
 
index 09e869623f0cb02c49f58d3734009ee2456e5b13..ca12e668510580d3439b34e6befb4f4653388bc9 100644 (file)
 // 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>
@@ -30,9 +35,6 @@
 #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;
index 095fee6f5ef7ed8be01e6e24bb393bb5a0f6b07b..327b17589dc43dc05bbb785d940a2b61449d6e37 100644 (file)
 // 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;
index ac96c5eb37bd8c6ba78e226dcd5efc8e80b2cdc1..cd10837641837ec054bfe23677713aba4b5e521a 100644 (file)
 #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"
index 838092d0f8127580753e455825bdc310534d9294..3fce2d9bf69c8815fbdbd9914079af543520b0b9 100644 (file)
 #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>
@@ -39,9 +44,6 @@
 
 #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"
index c0e802b5545182a43472c1dce3278f972c2bfcd2..120aecb5a56776b7a37782b6387030e2a82819fe 100644 (file)
 #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>
@@ -39,9 +44,6 @@
 
 #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"
index 6a98030110df8091100fec2739c3e260d2cce7e4..27368b8f356fb32ceedd415fb353bac182378643 100644 (file)
  */
 
 
+#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"
index cf7dc743d0d3ec64ff7a6726e949af1971bcd204..7e2bd1d6f10adc390aea9724ec12d51395dc4010 100644 (file)
  */
 
 
+#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"
index 63440a780134a2f16ef2ada54dfb1517950d886c..4c56033805bdbfad85089162d7553be93dd04746 100644 (file)
 // 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"
index 6bcb8ce79fa094ab212feec1a73736eea204a00c..dc23a7bc0cd9fecb3c3d8ff2de2c01603858e0c6 100644 (file)
 
 // 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>
@@ -26,9 +32,6 @@
 
 #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;
index 96577c03feb0a13ae830a9ace17bfe5c13095546..4bc83fa77b192e0d10e5c96b3fb0e344b1d47a9e 100644 (file)
 // 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;
index 9580d402a717ed4087e2b38a82221b627dacfaef..3e9051d4d4d10890dc03f627bd4a5553155cd94e 100644 (file)
 // 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;
index 830fa8a4e1404fbe22cfcedb652046555dd7c883..970dce5548c43ef096aca3b8235e389aa85c5e28 100644 (file)
 
 #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>
@@ -25,8 +29,6 @@
 #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"
index daeeb53115599067e14ec2322f22df85a0e28bb4..1e607dca78c67ee3d8b1e4e061c404e729028759 100644 (file)
 // 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"
index 9f128749ec55d5b1b1fabeaa0b5ddfaeda26b6a0..d73069d0a0e419a4e0dcf2d047fb377470738635 100644 (file)
 // 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"
index 5eee6625c60a7ef3f71b1cdfa09e367d619a0eab..5e9e00f13718063304f2c6c3e31101de443a0881 100644 (file)
@@ -18,6 +18,8 @@
 // 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;
index 5d613014aeb7577c2d0c573384eeca7bc2ec5f3c..52e750f4efe8e3bf71ad0773734bc1e516839ffb 100644 (file)
@@ -16,6 +16,8 @@
 
 // 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"
index cfe1180855b3696b56ddfa2f1ec0fe54dcedb8b9..126faf65d39af1f08fc7ca37c9971ee586adbc7a 100644 (file)
 
 // 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;
index 56a2020e412966ec6e7c2214b074ad292392d12e..2828aa18827fbbea8266f410d2935c81fcf67361 100644 (file)
@@ -19,6 +19,7 @@
 
 #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"
index c117671b3d73aab9c9ce6e9e4bb0ea9a1ac85ca6..f2906722bd4a377b508bca321b3ff082d3c34446 100644 (file)
 
 #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>
 
index c501094a14876cbf49e6441629815c8d79ce2eac..96ffeb06c0d050a53ff650ef9e95d055405bb249 100644 (file)
 // 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;
index 7d6f6a1065423268d17edc6569ccabf4c9bec793..60f4cf89a5fb7552131a9ff07389db9e3a0f3ca7 100644 (file)
@@ -18,8 +18,7 @@
 
 
 #include <deal.II/base/qprojector.h>
-
-#include <deal.II/simplex/quadrature_lib.h>
+#include <deal.II/base/quadrature_lib.h>
 
 #include "../tests.h"
 
index 67010b0a31e226b03211b58e2363349eb48bb942..b521142c98cec87c146aac6e58a4529900c59f41 100644 (file)
@@ -14,8 +14,7 @@
 // ---------------------------------------------------------------------
 
 #include <deal.II/base/function_lib.h>
-
-#include <deal.II/simplex/quadrature_lib.h>
+#include <deal.II/base/quadrature_lib.h>
 
 #include "../tests.h"
 
index fa05a42fd0946d1a05693c8c4407c13811cbef17..a3e417b5a7eabf5cab896e9ad1a3900e97e47337 100644 (file)
@@ -17,7 +17,7 @@
 // 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"
 
index 0d3e7fe611dc182096a072ba5cbc58ff41f920d1..7885dc9b343f517361602884efe3a931242f01b0 100644 (file)
@@ -17,9 +17,9 @@
 // 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"
 
index 5346fe9a191ed22862839064d023d50738dfbfbd..aa52255570236890d69de5b12ef33efacc093b5b 100644 (file)
 
 // 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"
index d74bfe73391768d3e00639236594df12ba25c8a8..984c499cd46751854a7b31af790cf6c7b26c6702 100644 (file)
 
 // 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"
index 5236adfd48e948a348593606b1241275d7b53eaa..b72cd179471765782d7edce8064fe3257b8b43da 100644 (file)
 // 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"
 
 
index b9a41e241824f3daee6be7c21590ce0494212d32..5ce2637c7b0007f12aa2889e896269f714b3d434 100644 (file)
@@ -27,8 +27,6 @@
 #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>
index 32cf8295d457a667d469d1efb1abdb24bedc7faa..7798f2e45c6373d454d8dcd7e3465fa099a680d8 100644 (file)
 #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>
@@ -34,8 +39,6 @@
 #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"
index ae8a03e29d26c9623a92c854159854e4d2ceb95a..6924433f29ae1cb787633da9adf96ae4c24ace38 100644 (file)
 #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>
 
index 892c1efb0449cea9482fa9ed08c97a6f2ca42b71..60d0b14681bf2b599584eb09e79de571725d7b80 100644 (file)
 #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>
 
index d3559515918b9e48d645401219964ab0cc59f897..d13a24d4d919ca7bb4234a296cf03b38fe1e8357 100644 (file)
 
 
 #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>
 
index 4cbbd6a9fafc47be5a1b873da421b76b679e2a22..2ffa48a28c60911fb35b325e6143e1485319cb04 100644 (file)
 #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>
@@ -58,9 +62,6 @@
 #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>
index eececbf1b465124944932f33769939587381ccb6..483a891611cd357cc97bf42d05f8ba9f8f483d02 100644 (file)
 #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>
 
index 94d64a085605e9751810a326d3fc70ed0baddebc..5e8936a356ad0dc67758154634e924255a3ee8fd 100644 (file)
 
 #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>
 
index 29484a3aa275befba6168b0a641b17f5e1686417..40df59beefbee64df79f1d35a6f8e0053a7d178a 100644 (file)
 
 // 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>
 
index f2459a1000870424b56b3550a5d4e4153fa3f9bd..25cd89769981298bd428876b557e3487554cbbde 100644 (file)
 #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
 
index 34e482fe7c870b013fb10d1017af6e5d7b954507..c16551e69fa9892cc1b93f24294336a09a5d747d 100644 (file)
 #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>
 
index 1a8dee444d255d681624611db3df5f06fe503929..52af8814de587e0a89463f2f2bda12d64c61a9ef 100644 (file)
 #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
 
index cda26794bb4d3979a26afc940afc6516120f7ae9..e47803bce1c2034ac0beff760b50d2d9f8781fc9 100644 (file)
 #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>
index 76561d44bc2522a243f8598a737a4e08b1c40412..c657613c0cafac517ec04dd5395497b0b1b0707d 100644 (file)
@@ -86,12 +86,15 @@ namespace LA
 #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
 
 
index 8bfbc11d637427c1854805afa7baed33a7938d46..861b224e97c3f32c0015bc629dac0e8e2cf731ff 100644 (file)
@@ -30,6 +30,7 @@
 #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>
 
index 7c9ce0553c3ed08ace6c86cb313aed71073b5e8d..340cd3c42a0a9aa52699fbbfce5d50f81e563e80 100644 (file)
 
 #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>
 
@@ -56,9 +60,6 @@
 #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
index 61fec47d1d2a325da7a760817c2ed88f6ba73ab2..d72e26a48a85c500cc1da69c7f8addb1b06919d0 100644 (file)
@@ -18,9 +18,9 @@
 // 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"
 
index 95855659f9a7996f17387b7c1254ba295fb709fa..6560f1515a428880f4eb766994d41f4050582662 100644 (file)
 
 
 #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;
index 5629248bee807e20a7df63317364043258db6f8f..62ecb8ddfb6565fc913b915ac464f8b0ca7b5a3b 100644 (file)
 
 #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;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.