lagrangepoints.hh

#pragma once
 
#include <cassert>
#include <array>
 
#include <dune/common/exceptions.hh>
#include <dune/common/version.hh>
#include <dune/geometry/type.hh>
#include <dune/localfunctions/lagrange/equidistantpoints.hh>
 
namespace Dune::Vtk
{
  namespace Impl
  {
    // forward declaration
    template <class K, unsigned int dim>
    class LagrangePointSetBuilder;
  }
 
 
 
  template <class K, unsigned int dim>
  class LagrangePointSet
      : public EmptyPointSet<K, dim>
  {
    using Super = EmptyPointSet<K, dim>;
 
  public:
    static const unsigned int dimension = dim;
 
    LagrangePointSet (std::size_t order)
      : Super(order)
    {
      assert(order > 0);
    }
 
    void build (GeometryType gt)
    {
      assert(gt.dim() == dimension);
      builder_(gt, order(), points_);
    }
 
    template <GeometryType::Id geometryId>
    bool build ()
    {
      build(GeometryType(geometryId));
      return true;
    }
 
    template <GeometryType::Id geometryId>
    static bool supports (std::size_t /*order*/)
    {
      return true;
    }
 
    using Super::order;
 
  private:
    using Super::points_;
    Impl::LagrangePointSetBuilder<K,dim> builder_;
  };
 
 
  namespace Impl
  {
    // Build for lagrange point sets in different dimensions
    // Specialized for dim=1,2,3
    template <class K, unsigned int dim>
    class LagrangePointSetBuilder
    {
    public:
      template <class Points>
      void operator()(GeometryType, unsigned int, Points& points) const
      {
        DUNE_THROW(Dune::NotImplemented,
          "Lagrange points not yet implemented for this GeometryType.");
      }
    };
 
 
    // Lagrange points on point geometries
    template <class K>
    class LagrangePointSetBuilder<K,0>
    {
      static constexpr int dim = 0;
      using LP = LagrangePoint<K,dim>;
      using Vec = typename LP::Vector;
      using Key = LocalKey;
 
    public:
      template <class Points>
      void operator()(GeometryType gt, int /*order*/, Points& points) const;
    };
 
 
    // Lagrange points on line geometries
    template <class K>
    class LagrangePointSetBuilder<K,1>
    {
      static constexpr int dim = 1;
      using LP = LagrangePoint<K,dim>;
      using Vec = typename LP::Vector;
      using Key = LocalKey;
 
    public:
      template <class Points>
      void operator()(GeometryType gt, int order, Points& points) const;
    };
 
 
    // Lagrange points on 2d geometries
    template <class K>
    class LagrangePointSetBuilder<K,2>
    {
      static constexpr int dim = 2;
      using LP = LagrangePoint<K,dim>;
      using Vec = typename LP::Vector;
      using Key = LocalKey;
 
      friend class LagrangePointSetBuilder<K,3>;
 
    public:
      template <class Points>
      void operator()(GeometryType gt, int order, Points& points) const;
 
      // implementation details
      // ----------------------
 
      // Construct the point set in a triangle element.
      // Loop from the outside to the inside
      template <class Points>
      void buildTriangle (int order, Points& points) const;
 
      // "Barycentric index" is a triplet of integers, each running from 0 to
      // <Order>. It is the index of a point on the triangle in barycentric
      // coordinates.
      static void barycentricIndex (int index, std::array<int,3>& bindex, int order);
 
      // Construct the point set in the quad element
      // 1. build equispaced points with index tuple (i,j)
      // 2. map index tuple to DOF index and LocalKey
      template <class Points>
      void buildQuad(int orderx, int ordery, Points& points) const;
 
      // Obtain the VTK DOF index of the node (i,j) in the quad element
      // and construct a LocalKey
      static std::pair<int,Key> calcQuadKey (int i, int j, std::array<int,2> order);
    };
 
 
    // Lagrange points on 3d geometries
    template <class K>
    class LagrangePointSetBuilder<K,3>
    {
      static constexpr int dim = 3;
      using LP = LagrangePoint<K,dim>;
      using Vec = typename LP::Vector;
      using Key = LocalKey;
 
    public:
      template <class Points>
      void operator() (GeometryType gt, unsigned int order, Points& points) const;
 
      // implementation details
      // ----------------------
 
      // Construct the point set in the tetrahedron element
      // 1. construct barycentric (index) coordinates
      // 2. obtains the DOF index, LocalKey and actual coordinate from barycentric index
      template <class Points>
      void buildTetra (int order, Points& points) const;
 
      // "Barycentric index" is a set of 4 integers, each running from 0 to
      // <Order>. It is the index of a point in the tetrahedron in barycentric
      // coordinates.
      static void barycentricIndex (int p, std::array<int,4>& bindex, int order);
 
      // Construct the point set in the heyhedral element
      // 1. build equispaced points with index tuple (i,j,k)
      // 2. map index tuple to DOF index and LocalKey
      template <class Points>
      void buildHex (int order, Points& points) const;
 
      // Obtain the VTK DOF index of the node (i,j,k) in the hexahedral element
      static std::pair<int,Key> calcHexKey (int i, int j, int k, std::array<int,3> order);
 
      // Construct the point set in the prism element
      // with polynomial order `order_tri` in the triangular faces and along the
      // z-direction with `order_quad` polynomial order.
      template <class Points>
      void buildWedge (int order_tri, int order_quad, Points& points) const;
    };
 
    template <class K, unsigned int dim>
    class LexicographicPointSetBuilder
    {
    public:
      // equidistant points along the line with `a+1` points.
      template <class Points>
      static void buildLine(int a, Points& points);
 
      // equidistant points in the quad with `a+1` points along the x-axis and
      // `b+1` points along the y-axis, lexicographically ordered.
      template <class Points>
      static void buildQuad(int a, int b, Points& points);
 
      // equidistant points in the hexahedron with `a+1` points along the x-axis,
      // `b+1` points along the y-axis, and `c+1` points along the z-axis,
      // lexicographically ordered.
      template <class Points>
      static void buildHex(int a, int b, int c, Points& points);
    };
 
  } // end namespace Impl
} // end namespace Dune::Vtk
 
#include <dune/vtk/utility/lagrangepoints.impl.hh>