From: Peter Munch <peterrmuench@gmail.com>
Date: Sat, 5 Dec 2020 12:00:33 +0000 (+0100)
Subject: Implement ScalarPyramidPolynomial and FE_PyramidP
X-Git-Tag: v9.3.0-rc1~781^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=953765f7af66e5634452301bee731591d014b493;p=dealii.git

Implement ScalarPyramidPolynomial and FE_PyramidP
---

diff --git a/include/deal.II/simplex/fe_lib.h b/include/deal.II/simplex/fe_lib.h
index e7c45a8ed9..030b89d526 100644
--- a/include/deal.II/simplex/fe_lib.h
+++ b/include/deal.II/simplex/fe_lib.h
@@ -222,6 +222,66 @@ namespace Simplex
   };
 
 
+  /**
+   * Implementation of a scalar Lagrange finite element on a pyramid that yields
+   * the finite element space of continuous, piecewise polynomials of
+   * degree p.
+   *
+   * @ingroup simplex
+   */
+  template <int dim, int spacedim = dim>
+  class FE_PyramidP : public dealii::FE_Poly<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;
+  };
+
+
 } // namespace Simplex
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/include/deal.II/simplex/polynomials.h b/include/deal.II/simplex/polynomials.h
index a1c484aaa5..051877bdc7 100644
--- a/include/deal.II/simplex/polynomials.h
+++ b/include/deal.II/simplex/polynomials.h
@@ -266,6 +266,119 @@ 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()
+     *
+     * @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>
@@ -303,6 +416,25 @@ 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
diff --git a/source/simplex/fe_lib.cc b/source/simplex/fe_lib.cc
index abd6e59042..9312656e3c 100644
--- a/source/simplex/fe_lib.cc
+++ b/source/simplex/fe_lib.cc
@@ -115,6 +115,32 @@ 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_pyramid(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}, {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;
+    }
   } // namespace
 
 
@@ -559,6 +585,153 @@ namespace Simplex
 
 
 
+  template <int dim, int spacedim>
+  FE_PyramidP<dim, spacedim>::FE_PyramidP(const unsigned int degree)
+    : dealii::FE_Poly<dim, spacedim>(
+        Simplex::ScalarPyramidPolynomial<dim>(degree),
+        FiniteElementData<dim>(get_dpo_vector_fe_pyramid(degree),
+                               ReferenceCell::Type::Pyramid,
+                               1,
+                               degree,
+                               FiniteElementData<dim>::H1),
+        std::vector<bool>(FiniteElementData<dim>(get_dpo_vector_fe_pyramid(
+                                                   degree),
+                                                 ReferenceCell::Type::Pyramid,
+                                                 1,
+                                                 degree)
+                            .dofs_per_cell,
+                          true),
+        std::vector<ComponentMask>(
+          FiniteElementData<dim>(get_dpo_vector_fe_pyramid(degree),
+                                 ReferenceCell::Type::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>
+  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
+  {
+    (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_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;
+  }
+
+
+
 } // namespace Simplex
 
 // explicit instantiations
diff --git a/source/simplex/fe_lib.inst.in b/source/simplex/fe_lib.inst.in
index 160d876f95..36b7279cd3 100644
--- a/source/simplex/fe_lib.inst.in
+++ b/source/simplex/fe_lib.inst.in
@@ -23,5 +23,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS)
     template class Simplex::FE_DGP<deal_II_dimension, deal_II_space_dimension>;
     template class Simplex::FE_WedgeP<deal_II_dimension,
                                       deal_II_space_dimension>;
+    template class Simplex::FE_PyramidP<deal_II_dimension,
+                                        deal_II_space_dimension>;
 #endif
   }
diff --git a/source/simplex/polynomials.cc b/source/simplex/polynomials.cc
index bae98e5f95..dc9ec63a43 100644
--- a/source/simplex/polynomials.cc
+++ b/source/simplex/polynomials.cc
@@ -52,6 +52,21 @@ namespace Simplex
       return 0;
     }
 
+    unsigned int
+    compute_n_polynomials_pyramid(const unsigned int dim,
+                                  const unsigned int degree)
+    {
+      if (dim == 3)
+        {
+          if (degree == 1)
+            return 5;
+        }
+
+      Assert(false, ExcNotImplemented());
+
+      return 0;
+    }
+
     unsigned int
     compute_n_polynomials_wedge(const unsigned int dim,
                                 const unsigned int degree)
@@ -711,12 +726,255 @@ namespace Simplex
 
 
 
+  template <int dim>
+  ScalarPyramidPolynomial<dim>::ScalarPyramidPolynomial(
+    const unsigned int degree)
+    : ScalarPolynomialsBase<dim>(degree,
+                                 compute_n_polynomials_pyramid(dim, degree))
+  {}
+
+
+  template <int dim>
+  double
+  ScalarPyramidPolynomial<dim>::compute_value(const unsigned int i,
+                                              const Point<dim> & p) const
+  {
+    AssertDimension(dim, 3);
+    AssertIndexRange(this->degree(), 2);
+
+    const double Q14 = 0.25;
+    double       ration;
+
+    const double r = p[0];
+    const double s = p[1];
+    const double t = p[2];
+
+    if (fabs(t - 1.0) > 1.0e-14)
+      {
+        ration = (r * s * t) / (1.0 - t);
+      }
+    else
+      {
+        ration = 0.0;
+      }
+
+    if (i == 0)
+      return Q14 * ((1.0 - r) * (1.0 - s) - t + ration);
+    if (i == 1)
+      return Q14 * ((1.0 + r) * (1.0 - s) - t - ration);
+    if (i == 2)
+      return Q14 * ((1.0 - r) * (1.0 + s) - t - ration);
+    if (i == 3)
+      return Q14 * ((1.0 + r) * (1.0 + s) - t + ration);
+    else
+      return t;
+  }
+
+
+
+  template <int dim>
+  Tensor<1, dim>
+  ScalarPyramidPolynomial<dim>::compute_grad(const unsigned int i,
+                                             const Point<dim> & p) const
+  {
+    AssertDimension(dim, 3);
+    AssertIndexRange(this->degree(), 4);
+
+    Tensor<1, dim> grad;
+
+    if (this->degree() == 1)
+      {
+        const double Q14 = 0.25;
+
+        const double r = p[0];
+        const double s = p[1];
+        const double t = p[2];
+
+        double rationdr;
+        double rationds;
+        double rationdt;
+
+        if (fabs(t - 1.0) > 1.0e-14)
+          {
+            rationdr = s * t / (1.0 - t);
+            rationds = r * t / (1.0 - t);
+            rationdt = r * s / ((1.0 - t) * (1.0 - t));
+          }
+        else
+          {
+            rationdr = 1.0;
+            rationds = 1.0;
+            rationdt = 1.0;
+          }
+
+
+        if (i == 0)
+          {
+            grad[0] = Q14 * (-1.0 * (1.0 - s) + rationdr);
+            grad[1] = Q14 * (-1.0 * (1.0 - r) + rationds);
+            grad[2] = Q14 * (rationdt - 1.0);
+          }
+        else if (i == 1)
+          {
+            grad[0] = Q14 * (1.0 * (1.0 - s) - rationdr);
+            grad[1] = Q14 * (-1.0 * (1.0 + r) - rationds);
+            grad[2] = Q14 * (-1.0 * rationdt - 1.0);
+          }
+        else if (i == 2)
+          {
+            grad[0] = Q14 * (-1.0 * (1.0 + s) - rationdr);
+            grad[1] = Q14 * (1.0 * (1.0 - r) - rationds);
+            grad[2] = Q14 * (-1.0 * rationdt - 1.0);
+          }
+        else if (i == 3)
+          {
+            grad[0] = Q14 * (1.0 * (1.0 + s) + rationdr);
+            grad[1] = Q14 * (1.0 * (1.0 + r) + rationds);
+            grad[2] = Q14 * (rationdt - 1.0);
+          }
+        else if (i == 4)
+          {
+            grad[0] = 0.0;
+            grad[1] = 0.0;
+            grad[2] = 1.0;
+          }
+        else
+          {
+            Assert(false, ExcNotImplemented());
+          }
+      }
+
+    return grad;
+  }
+
+
+
+  template <int dim>
+  Tensor<2, dim>
+  ScalarPyramidPolynomial<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
+  ScalarPyramidPolynomial<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>
+  ScalarPyramidPolynomial<dim>::compute_1st_derivative(
+    const unsigned int i,
+    const Point<dim> & p) const
+  {
+    return compute_grad(i, p);
+  }
+
+
+
+  template <int dim>
+  Tensor<2, dim>
+  ScalarPyramidPolynomial<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>
+  ScalarPyramidPolynomial<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>
+  ScalarPyramidPolynomial<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
+  ScalarPyramidPolynomial<dim>::name() const
+  {
+    return "ScalarPyramidPolynomial";
+  }
+
+
+
+  template <int dim>
+  std::unique_ptr<ScalarPolynomialsBase<dim>>
+  ScalarPyramidPolynomial<dim>::clone() const
+  {
+    return std::make_unique<ScalarPyramidPolynomial<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>;
+  template class ScalarPyramidPolynomial<1>;
+  template class ScalarPyramidPolynomial<2>;
+  template class ScalarPyramidPolynomial<3>;
 
 } // namespace Simplex