lagrangepoints.hh 6.27 KB
Newer Older
1
2
3
4
5
6
#pragma once

#include <cassert>
#include <array>

#include <dune/common/exceptions.hh>
7
#include <dune/common/version.hh>
8
9
10
#include <dune/geometry/type.hh>
#include <dune/localfunctions/lagrange/equidistantpoints.hh>

11
12
namespace Dune
{
13
  namespace Gmsh4
14
  {
15
    namespace Impl
16
    {
17
18
19
      // forward declaration
      template <class K, unsigned int dim>
      class LagrangePointSetBuilder;
20
21
22
    }


23
24
25
26
27
28
29
30
    /// \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>
31
    {
32
      using Super = EmptyPointSet<K, dim>;
33

34
35
    public:
      static const unsigned int dimension = dim;
36

37
38
39
40
      LagrangePointSet (std::size_t order)
        : Super(order)
      {
        assert(order > 0);
41
42
      }

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

50
      /// Fill the lagrange points for the given topology type `Topology`
51
#if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
52
53
54
55
56
57
      template <class Topology>
      bool build ()
      {
        build(GeometryType(Topology{}));
        return true;
      }
58
59
60
61
62
63
64
65
#else
      template <GeometryType::Id geometryId>
      bool build ()
      {
        build(GeometryType(geometryId));
        return true;
      }
#endif
66

67
68
      /// Returns whether the point set support the given topology type `Topology` and can
      /// generate point for the given order.
69
#if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
70
      template <class Topology>
71
72
73
#else
      template <GeometryType::Id geometryId>
#endif
74
75
76
77
      static bool supports (std::size_t order)
      {
        return true;
      }
78

79
      using Super::order;
80

81
82
83
84
    private:
      using Super::points_;
      Impl::LagrangePointSetBuilder<K,dim> builder_;
    };
85
86


87
    namespace Impl
88
    {
89
90
91
92
93
94
95
96
97
98
99
100
101
      // 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.");
        }
      };
102
103


104
105
106
107
108
109
110
111
112
113
114
115
116
      // 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;
      };
117
118


119
120
121
122
123
124
125
126
      // 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;
127

128
129
130
131
      public:
        template <class Points>
        void operator()(GeometryType gt, int order, Points& points) const;
      };
132
133


134
135
136
      // Lagrange points on 2d geometries
      template <class K>
      class LagrangePointSetBuilder<K,2>
137
      {
138
139
140
141
        static constexpr int dim = 2;
        using LP = LagrangePoint<K,dim>;
        using Vec = typename LP::Vector;
        using Key = LocalKey;
142

143
        friend class LagrangePointSetBuilder<K,3>;
144

145
146
147
      public:
        template <class Points>
        void operator()(GeometryType gt, int order, Points& points) const;
148

149
      private: // implementation details
150

151
152
153
154
        // 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;
155

156
157
158
159
        // "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);
160

161
162
163
164
165
        // 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;
166

167
168
169
170
        // 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);
      };
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
200
201
202
203
204
205
206
207
208
209
210
211
212
213
      // 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 Gmsh4
} // end namespace Dune

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