lagrangepoints.hh 5.93 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1
2
3
4
5
6
7
8
9
#pragma once

#include <cassert>
#include <array>

#include <dune/common/exceptions.hh>
#include <dune/geometry/type.hh>
#include <dune/localfunctions/lagrange/equidistantpoints.hh>

Stenger, Florian's avatar
Stenger, Florian committed
10
11
namespace Dune
{
12
  namespace Vtk
Praetorius, Simon's avatar
Praetorius, Simon committed
13
  {
14
    namespace Impl
Praetorius, Simon's avatar
Praetorius, Simon committed
15
    {
16
17
18
      // forward declaration
      template <class K, unsigned int dim>
      class LagrangePointSetBuilder;
Praetorius, Simon's avatar
Praetorius, Simon committed
19
20
21
    }


22
23
24
25
26
27
28
29
    /// \brief A set of lagrange points compatible with the numbering of VTK and Gmsh
    /**
     * \tparam K    Field-type for the coordinates
     * \tparam dim  Dimension of the coordinates
     **/
    template <class K, unsigned int dim>
    class LagrangePointSet
        : public EmptyPointSet<K, dim>
Praetorius, Simon's avatar
Praetorius, Simon committed
30
    {
31
      using Super = EmptyPointSet<K, dim>;
Praetorius, Simon's avatar
Praetorius, Simon committed
32

33
34
    public:
      static const unsigned int dimension = dim;
Praetorius, Simon's avatar
Praetorius, Simon committed
35

36
37
38
39
      LagrangePointSet (std::size_t order)
        : Super(order)
      {
        assert(order > 0);
Praetorius, Simon's avatar
Praetorius, Simon committed
40
      }
Stenger, Florian's avatar
Stenger, Florian committed
41

42
43
44
45
46
      /// Fill the lagrange points for the given geometry type
      void build (GeometryType gt)
      {
        assert(gt.dim() == dimension);
        builder_(gt, order(), points_);
Praetorius, Simon's avatar
Praetorius, Simon committed
47
48
      }

49
50
51
52
53
54
55
      /// Fill the lagrange points for the given topology type `Topology`
      template <class Topology>
      bool build ()
      {
        build(GeometryType(Topology{}));
        return true;
      }
Praetorius, Simon's avatar
Praetorius, Simon committed
56

57
58
59
60
61
62
63
      /// Returns whether the point set support the given topology type `Topology` and can
      /// generate point for the given order.
      template <class Topology>
      static bool supports (std::size_t order)
      {
        return true;
      }
Praetorius, Simon's avatar
Praetorius, Simon committed
64

65
      using Super::order;
Praetorius, Simon's avatar
Praetorius, Simon committed
66

67
68
69
70
    private:
      using Super::points_;
      Impl::LagrangePointSetBuilder<K,dim> builder_;
    };
Praetorius, Simon's avatar
Praetorius, Simon committed
71
72


73
    namespace Impl
Praetorius, Simon's avatar
Praetorius, Simon committed
74
    {
75
76
77
78
79
80
81
82
83
84
85
86
87
      // 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.");
        }
      };
Praetorius, Simon's avatar
Praetorius, Simon committed
88

Stenger, Florian's avatar
Stenger, Florian committed
89

90
91
92
93
94
95
96
97
98
99
100
101
102
      // 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;
      };
Praetorius, Simon's avatar
Praetorius, Simon committed
103
104


105
106
107
108
109
110
111
112
      // 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;
Praetorius, Simon's avatar
Praetorius, Simon committed
113

114
115
116
117
      public:
        template <class Points>
        void operator()(GeometryType gt, int order, Points& points) const;
      };
Praetorius, Simon's avatar
Praetorius, Simon committed
118
119


120
121
122
      // Lagrange points on 2d geometries
      template <class K>
      class LagrangePointSetBuilder<K,2>
Praetorius, Simon's avatar
Praetorius, Simon committed
123
      {
124
125
126
127
        static constexpr int dim = 2;
        using LP = LagrangePoint<K,dim>;
        using Vec = typename LP::Vector;
        using Key = LocalKey;
Praetorius, Simon's avatar
Praetorius, Simon committed
128

129
        friend class LagrangePointSetBuilder<K,3>;
Stenger, Florian's avatar
Stenger, Florian committed
130

131
132
133
      public:
        template <class Points>
        void operator()(GeometryType gt, int order, Points& points) const;
Praetorius, Simon's avatar
Praetorius, Simon committed
134

135
      private: // implementation details
Praetorius, Simon's avatar
Praetorius, Simon committed
136

137
138
139
140
        // Construct the point set in a triangle element.
        // Loop from the outside to the inside
        template <class Points>
        void buildTriangle (std::size_t nPoints, int order, Points& points) const;
Praetorius, Simon's avatar
Praetorius, Simon committed
141

142
143
144
145
        // "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);
Praetorius, Simon's avatar
Praetorius, Simon committed
146

147
148
149
150
151
        // 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(std::size_t nPoints, int order, Points& points) const;
Praetorius, Simon's avatar
Praetorius, Simon committed
152

153
154
155
156
        // 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);
      };
Praetorius, Simon's avatar
Praetorius, Simon committed
157
158


159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
      // 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;

      private: // 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 (std::size_t nPoints, 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 (std::size_t 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 (std::size_t nPoints, 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);
      };

    } // end namespace Impl
  } // end namespace Vtk
} // end namespace Dune

#include <dune/vtk/utility/lagrangepoints.impl.hh>