]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Introduce FE_Wedge 11267/head
authorPeter Munch <peterrmuench@gmail.com>
Sat, 28 Nov 2020 08:17:39 +0000 (09:17 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Wed, 2 Dec 2020 11:55:18 +0000 (12:55 +0100)
include/deal.II/grid/reference_cell.h
include/deal.II/simplex/fe_lib.h
include/deal.II/simplex/polynomials.h
source/simplex/fe_lib.cc
source/simplex/fe_lib.inst.in
source/simplex/polynomials.cc
tests/simplex/wedge_01.cc [new file with mode: 0644]
tests/simplex/wedge_01.with_simplex_support=on.output [new file with mode: 0644]

index 32d1c49b729adb5c7db6c4e59071a3de71f839b6..5402c73267f2b23eeec4635cfd3aca7614f9accb 100644 (file)
@@ -822,8 +822,6 @@ namespace ReferenceCell
         standard_line_to_face_and_line_index(
           const unsigned int line) const override
         {
-          Assert(false, ExcNotImplemented());
-
           static const std::array<unsigned int, 2> table[6] = {
             {{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}};
 
@@ -986,10 +984,15 @@ namespace ReferenceCell
         standard_line_to_face_and_line_index(
           const unsigned int line) const override
         {
-          Assert(false, ExcNotImplemented());
-
-          static const std::array<unsigned int, 2> table[6] = {
-            {{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 1}}, {{1, 2}}, {{2, 1}}};
+          static const std::array<unsigned int, 2> table[9] = {{{0, 0}},
+                                                               {{0, 2}},
+                                                               {{0, 1}},
+                                                               {{1, 0}},
+                                                               {{1, 1}},
+                                                               {{1, 2}},
+                                                               {{2, 0}},
+                                                               {{2, 1}},
+                                                               {{3, 1}}};
 
           return table[line];
         }
@@ -1000,19 +1003,26 @@ namespace ReferenceCell
           const unsigned int  face,
           const unsigned char face_orientation) const override
         {
-          Assert(false, ExcNotImplemented());
-
-          (void)face;
-
-          static const std::array<std::array<unsigned int, 3>, 6> table = {
-            {{{2, 1, 0}},
-             {{0, 1, 2}},
-             {{1, 2, 0}},
-             {{0, 2, 1}},
-             {{1, 0, 2}},
-             {{2, 0, 1}}}};
+          if (face > 1) // QUAD
+            {
+              return GeometryInfo<3>::standard_to_real_face_line(
+                line,
+                get_bit(face_orientation, 0),
+                get_bit(face_orientation, 2),
+                get_bit(face_orientation, 1));
+            }
+          else // TRI
+            {
+              static const std::array<std::array<unsigned int, 3>, 6> table = {
+                {{{2, 1, 0}},
+                 {{0, 1, 2}},
+                 {{1, 2, 0}},
+                 {{0, 2, 1}},
+                 {{1, 0, 2}},
+                 {{2, 0, 1}}}};
 
-          return table[face_orientation][line];
+              return table[face_orientation][line];
+            }
         }
 
         bool
index f618ff8e5b3d90b95250a3261f8bafb6800fa5ca..e7c45a8ed97a858a3ac128d307114f509ea93706 100644 (file)
@@ -161,6 +161,67 @@ namespace Simplex
       const FiniteElement<dim, spacedim> &fe_other) const override;
   };
 
+
+  /**
+   * Implementation of a scalar Lagrange finite element on a wedge that yields
+   * the finite element space of continuous, piecewise polynomials of
+   * degree p.
+   *
+   * @ingroup simplex
+   */
+  template <int dim, int spacedim = dim>
+  class FE_WedgeP : public dealii::FE_Poly<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;
+  };
+
+
 } // namespace Simplex
 
 DEAL_II_NAMESPACE_CLOSE
index 7561d28a3afb43227ee6bf8b6c60eb96764f5f3d..a1c484aaa56c51ab3bb07a4d7b55e8c6246133cc 100644 (file)
@@ -145,6 +145,127 @@ namespace Simplex
     clone() const override;
   };
 
+
+
+  /**
+   * Polynomials defined on wedge entities. This class is basis of
+   * Simplex::FE_WedgeP.
+   *
+   * The polynomials are created via a tensor product of a
+   * Simplex::ScalarPolynomial<2>(degree) and a
+   * Simplex::ScalarPolynomial<1>(degree), however, are re-numerated to better
+   * match the definition of FiniteElement.
+   */
+  template <int dim>
+  class ScalarWedgePolynomial : 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 (degree=1) and quadratic polynomials
+     *   (degree=2) are implemented.
+     */
+    ScalarWedgePolynomial(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()
+     *
+     * @note Not implemented yet.
+     */
+    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 functions
   template <int dim>
   template <int order>
@@ -162,6 +283,26 @@ namespace Simplex
 
     return der;
   }
+
+
+
+  template <int dim>
+  template <int order>
+  Tensor<order, dim>
+  ScalarWedgePolynomial<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 a034abd2344513a7008fb0ec65a018003e29a097..abd6e5904290d3147afa265f8ec017deb70211e4 100644 (file)
@@ -80,6 +80,41 @@ 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(const unsigned int degree)
+    {
+      internal::GenericDoFsPerObject dpo;
+
+      // all dofs are internal
+      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;
+    }
   } // namespace
 
 
@@ -371,6 +406,159 @@ namespace Simplex
 
 
 
+  template <int dim, int spacedim>
+  FE_WedgeP<dim, spacedim>::FE_WedgeP(const unsigned int degree)
+    : dealii::FE_Poly<dim, spacedim>(
+        Simplex::ScalarWedgePolynomial<dim>(degree),
+        FiniteElementData<dim>(get_dpo_vector_fe_wedge(degree),
+                               ReferenceCell::Type::Wedge,
+                               1,
+                               degree,
+                               FiniteElementData<dim>::L2),
+        std::vector<bool>(FiniteElementData<dim>(get_dpo_vector_fe_wedge(
+                                                   degree),
+                                                 ReferenceCell::Type::Wedge,
+                                                 1,
+                                                 degree)
+                            .dofs_per_cell,
+                          true),
+        std::vector<ComponentMask>(
+          FiniteElementData<dim>(get_dpo_vector_fe_wedge(degree),
+                                 ReferenceCell::Type::Wedge,
+                                 1,
+                                 degree)
+            .dofs_per_cell,
+          std::vector<bool>(1, true)))
+  {
+    AssertDimension(dim, 3);
+
+
+    if (degree == 1)
+      {
+        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, 0.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);
+        this->unit_support_points.emplace_back(0.0, 0.0, 1.0);
+
+        // TODO
+        this->unit_face_support_points[0].emplace_back(1.0, 0.0);
+        this->unit_face_support_points[0].emplace_back(0.0, 1.0);
+        this->unit_face_support_points[0].emplace_back(0.0, 0.0);
+      }
+  }
+
+
+
+  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
+  {
+    (void)fe_other;
+    (void)codim;
+
+    Assert((dynamic_cast<const Simplex::FE_P<dim, spacedim> *>(&fe_other)) ||
+             (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other)),
+           ExcNotImplemented());
+
+    return FiniteElementDomination::this_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;
+  }
+
+
+
 } // namespace Simplex
 
 // explicit instantiations
index 0f29f676d1957b033b34456da9d360e49d04b65c..160d876f950334e18f7a25e6c8f56939493a1a93 100644 (file)
@@ -21,5 +21,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS)
     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_WedgeP<deal_II_dimension,
+                                      deal_II_space_dimension>;
 #endif
   }
index 5c33c826dcf816a003150beee660d055c9d9f96a..bae98e5f953b83eab62bcabd502cbdfec50dedb2 100644 (file)
@@ -51,6 +51,23 @@ 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
 
 
@@ -485,9 +502,221 @@ namespace Simplex
     return std::make_unique<ScalarPolynomial<dim>>(*this);
   }
 
+
+
+  template <int dim>
+  ScalarWedgePolynomial<dim>::ScalarWedgePolynomial(const unsigned int degree)
+    : ScalarPolynomialsBase<dim>(degree,
+                                 compute_n_polynomials_wedge(dim, degree))
+  {}
+
+
+  namespace
+  {
+    /**
+     * TODO
+     */
+    static const constexpr std::array<std::array<unsigned int, 2>, 6>
+      wedge_table_1{
+        {{{0, 0}}, {{1, 0}}, {{2, 0}}, {{0, 1}}, {{1, 1}}, {{2, 1}}}};
+
+    /**
+     * TODO
+     */
+    static const constexpr std::array<std::array<unsigned int, 2>, 18>
+      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
+  {
+    AssertDimension(dim, 3);
+    AssertIndexRange(this->degree(), 3);
+
+    const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i];
+
+    const ScalarPolynomial<2> poly_tri(this->degree());
+    const Point<2>            p_tri(p[0], p[1]);
+    const auto                v_tri = poly_tri.compute_value(pair[0], p_tri);
+
+    const ScalarPolynomial<1> poly_line(this->degree());
+    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
+  {
+    AssertDimension(dim, 3);
+    AssertIndexRange(this->degree(), 3);
+
+    const auto pair = this->degree() == 1 ? wedge_table_1[i] : wedge_table_2[i];
+
+    const ScalarPolynomial<2> poly_tri(this->degree());
+    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 ScalarPolynomial<1> poly_line(this->degree());
+    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 ScalarPolynomial<1>;
   template class ScalarPolynomial<2>;
   template class ScalarPolynomial<3>;
+  template class ScalarWedgePolynomial<1>;
+  template class ScalarWedgePolynomial<2>;
+  template class ScalarWedgePolynomial<3>;
 
 } // namespace Simplex
 
diff --git a/tests/simplex/wedge_01.cc b/tests/simplex/wedge_01.cc
new file mode 100644 (file)
index 0000000..5629248
--- /dev/null
@@ -0,0 +1,113 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Distribute Simplex::FE_Wedge on a DoFHandler.
+
+#include <deal.II/dofs/dof_handler.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;
+
+void
+test_3()
+{
+  const int dim      = 3;
+  const int spacedim = 3;
+
+  std::vector<Point<spacedim>> vertices;
+  std::vector<CellData<dim>>   cells;
+
+#if true
+  vertices.emplace_back(0.0, 0.0, 0.0);
+  vertices.emplace_back(1.0, 0.0, 0.0);
+  vertices.emplace_back(0.0, 1.0, 0.0);
+  vertices.emplace_back(0.0, 0.0, 1.0);
+  vertices.emplace_back(1.0, 0.0, 1.0);
+  vertices.emplace_back(0.0, 1.0, 1.0);
+
+  {
+    CellData<dim> cell;
+    cell.vertices = {0, 1, 2, 3, 4, 5};
+    cells.push_back(cell);
+  }
+#else
+  vertices.emplace_back(0.0, 0.0, 0.0);
+  vertices.emplace_back(1.0, 0.0, 0.0);
+  vertices.emplace_back(0.0, 1.0, 0.0);
+  vertices.emplace_back(1.0, 1.0, 0.0);
+  vertices.emplace_back(0.0, 0.0, 1.0);
+  vertices.emplace_back(1.0, 0.0, 1.0);
+  vertices.emplace_back(0.0, 1.0, 1.0);
+  vertices.emplace_back(1.0, 1.0, 1.0);
+
+  {
+    CellData<dim> cell;
+    cell.vertices = {0, 1, 2, 4, 5, 6};
+    cells.push_back(cell);
+  }
+  {
+    CellData<dim> cell;
+    cell.vertices = {2, 1, 3, 6, 5, 7};
+    cells.push_back(cell);
+  }
+#endif
+
+  Triangulation<dim, spacedim> tria;
+  tria.create_triangulation(vertices, cells, SubCellData());
+
+  GridOut       grid_out;
+  std::ofstream out("mesh.vtk");
+  grid_out.write_vtk(tria, out);
+
+  Simplex::FE_WedgeP<dim, spacedim> fe(2);
+
+  DoFHandler<dim, spacedim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  deallog << dof_handler.n_dofs() << std::endl;
+
+  std::vector<types::global_dof_index> dof_indices;
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      for (const auto face : cell->face_indices())
+        {
+#if false
+          dof_indices.resize(face <= 1 ? 3 : 4);
+#else
+          dof_indices.resize(fe.n_dofs_per_face(face));
+#endif
+          cell->face(face)->get_dof_indices(dof_indices);
+
+          for (const auto i : dof_indices)
+            deallog << i << " ";
+          deallog << std::endl;
+        }
+    }
+}
+
+int
+main()
+{
+  initlog();
+
+  test_3();
+}
diff --git a/tests/simplex/wedge_01.with_simplex_support=on.output b/tests/simplex/wedge_01.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..ad3705a
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL::18
+DEAL::1 0 2 6 8 7 
+DEAL::3 4 5 9 10 11 
+DEAL::0 1 3 4 12 13 6 9 15 
+DEAL::1 2 4 5 13 14 7 10 16 
+DEAL::2 0 5 3 14 12 11 8 17 
\ No newline at end of file

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.