lagrangepoints.hh 5.92 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
      // 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.
183
        static void barycentricIndex (int p, std::array<int,4>& bindex, int order);
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199

        // 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>