geometry.hh 8.23 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_CURVED_SURFACE_GRID_GEOMETRY_HH
#define DUNE_CURVED_SURFACE_GRID_GEOMETRY_HH

#include <cassert>
#include <functional>
#include <iterator>
#include <limits>
#include <vector>

#include <dune/common/diagonalmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/std/optional.hh>
#include <dune/common/std/type_traits.hh>

19
#include <dune/curvedgeometry/curvedgeometry.hh>
20
#include <dune/curvedgeometry/localfunctiongeometry.hh>
21

22
23
24
25
26
27
#include <dune/geometry/affinegeometry.hh>
#include <dune/geometry/multilineargeometry.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>

28
29
#include "localgeometrywrapper.hh"

30
31
32
33
34
35
namespace Dune
{
  namespace CGeo
  {
    /// \brief Curved geometry implementation based on local basis function parametrization
    /**
36
     *  Parametrization of the geometry by a differentiable localFunction
37
     *
38
39
40
     *  \tparam  LF      localFunction parametrizing the geometry
     *  \tparam  LG      localGeometry for coordinate transform from subEntity to element,
     *                   see \ref Impl::DefaultLocalGeometry and \ref Impl::LocalGeometryInterface
41
     *  \tparam  ORDER   Polynomial order of lagrange parametrization. If order < 0 use localFunction.
42
     */
43
    template <class ct, int mydim, int cdim, class LF, class LG, int ORDER = -1>
44
    class Geometry
45
        : public LagrangeCurvedGeometry<ct, mydim, cdim, (ORDER > 0 ? ORDER : 1)>
46
    {
47
      using Super = LagrangeCurvedGeometry<ct, mydim, cdim, (ORDER > 0 ? ORDER : 1)>;
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89

      //! type of the localFunction representation the geometry parametrization
      using LocalFunction = LF;

      //! type of coordinate transformation for subEntities to codim=0 entities
      using LocalGeometry = LG;

      /// type of reference element
      using ReferenceElements = Dune::ReferenceElements<ct, mydim>;
      using ReferenceElement = typename ReferenceElements::ReferenceElement;

    public:
      /// \brief constructor from localFunction to parametrize the geometry
      /**
       *  \param[in]  refElement     reference element for the geometry
       *  \param[in]  localFunction  localFunction for the parametrization of the geometry
       *                             (stored by value)
       *  \param[in]  args...        argument to construct the local geometry
       *
       */
      template <class... Args>
      Geometry (const ReferenceElement& refElement, const LocalFunction& localFunction, Args&&... args)
        : Super(refElement, [localFunction, localGeometry=LocalGeometry{std::forward<Args>(args)...}](auto const& local)
          {
            return localFunction(localGeometry.global(local));
          })
      {}

      /// \brief constructor, forwarding to the other constructor that take a reference-element
      /**
       *  \param[in]  gt             geometry type
       *  \param[in]  localFunction  localFunction for the parametrization of the geometry
       *                             (stored by value)
       *  \param[in]  args...        argument to construct the local geometry
       */
      template <class... Args>
      Geometry (Dune::GeometryType gt, const LocalFunction& localFunction, Args&&... args)
        : Geometry(ReferenceElements::general(gt), localFunction, std::forward<Args>(args)...)
      {}
    };


90
    // Specialization for the case that ORDER < 0: Use LocalFunction `LF` directly as parametrization
91
92
    template <class ct, int mydim, int cdim, class LF, class LG>
    class Geometry<ct,mydim,cdim,LF,LG,-1>
93
        : public LocalFunctionGeometry<LF,LG>
94
    {
95
      using Super = LocalFunctionGeometry<LF,LG>;
96
97
98
99

      //! type of the localFunction representation the geometry parametrization
      using LocalFunction = LF;

100
101
      //! type of coordinate transformation for subEntities to codim=0 entities
      using LocalGeometry = LG;
102

103
104
105
      /// type of reference element
      using ReferenceElements = Dune::ReferenceElements<ct, mydim>;
      using ReferenceElement = typename ReferenceElements::ReferenceElement;
106
107

    public:
108
      /// \brief constructor from localFunction to parametrize the geometry
109
      /**
110
111
112
113
       *  \param[in]  refElement     reference element for the geometry
       *  \param[in]  localFunction  localFunction for the parametrization of the geometry
       *                             (stored by value)
       *  \param[in]  args...        argument to construct the local geometry
114
115
       *
       */
116
117
      template <class... Args>
      Geometry (const ReferenceElement& refElement, const LocalFunction& localFunction, Args&&... args)
118
        : Super(refElement, localFunction, LocalGeometry{std::forward<Args>(args)...})
119
120
121
122
123
      {}

      /// \brief constructor, forwarding to the other constructor that take a reference-element
      /**
       *  \param[in]  gt             geometry type
124
125
126
       *  \param[in]  localFunction  localFunction for the parametrization of the geometry
       *                             (stored by value)
       *  \param[in]  args...        argument to construct the local geometry
127
       */
128
129
      template <class... Args>
      Geometry (Dune::GeometryType gt, const LocalFunction& localFunction, Args&&... args)
130
        : Super(gt, localFunction, LocalGeometry{std::forward<Args>(args)...})
131
132
133
134
135
136
137
      {}
    };


  #ifndef DOXYGEN

    // Specialization for vertex geometries
138
139
    template <class ct, int cdim, class LF, class LG>
    class VertexGeometry
140
141
142
143
        : public AffineGeometry<ct,0,cdim>
    {
      using Super = AffineGeometry<ct,0,cdim>;

144
145
146
      using LocalFunction = LF;
      using LocalGeometry = LG;

147
148
149
      //! tolerance to numerical algorithms
      static ct tolerance () { return ct(16) * std::numeric_limits<ct>::epsilon(); }

150
151
      struct Construct {};

152
153
154
155
156
157
158
159
160
161
162
163
    public:
      /// type of reference element
      using ReferenceElements = Dune::ReferenceElements<ct, 0>;
      using ReferenceElement = typename ReferenceElements::ReferenceElement;

      /// type of local coordinates
      using LocalCoordinate = FieldVector<ct, 0>;

      /// type of global coordinates
      using GlobalCoordinate = FieldVector<ct, cdim>;

    protected:
164
      VertexGeometry (Construct, const ReferenceElement& refElement, const LocalFunction& localFunction,
165
                      const LocalGeometry& lg)
166
        : Super(refElement, std::vector<GlobalCoordinate>{localFunction(lg.global(refElement.position(0,0)))})
167
168
169
      {}

    public:
170
171
172
      template <class... Args>
      VertexGeometry (const ReferenceElement& refElement, const LocalFunction& lf, Args&&... args)
        : VertexGeometry(Construct{}, refElement, lf, LocalGeometry(std::forward<Args>(args)...))
173
174
      {}

175
176
177
178
      template <class... Args>
      VertexGeometry (Dune::GeometryType gt, const LocalFunction& lf, Args&&... args)
        : VertexGeometry(Construct{}, ReferenceElements::general(gt), lf,
                   LocalGeometry(std::forward<Args>(args)...))
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
      {}

      Std::optional<LocalCoordinate> checkedLocal (const GlobalCoordinate& globalCoord) const
      {
        auto localCoord = Super::local(globalCoord);
        if ((globalCoord - Super::global(localCoord)).two_norm2() > tolerance())
          return {};

        return localCoord;
      }

      GlobalCoordinate normal (const LocalCoordinate& local) const
      {
        DUNE_THROW(Dune::NotImplemented,
          "ERROR: normal() method only defined for edges in 2D and faces in 3D");
        return GlobalCoordinate(0);
      }

197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
    private:
      VertexGeometry const& flatGeometry () const { return *this; }
    };


    // Specialization for vertex geometries
    template <class ct, int cdim, class LF, class LG, int order>
    class Geometry<ct,0,cdim,LF,LG,order>
        : public VertexGeometry<ct,cdim,LF,LG>
    {
      using Super = VertexGeometry<ct,cdim,LF,LG>;
    public:
      using Super::Super;
    };

    template <class ct, int cdim, class LF, class LG>
    class Geometry<ct,0,cdim,LF,LG,-1>
        : public VertexGeometry<ct,cdim,LF,LG>
    {
      using Super = VertexGeometry<ct,cdim,LF,LG>;
217
    public:
218
      using Super::Super;
219
220
221
222
223
224
225
226
    };

#endif // DOXYGEN

  } // end namespace CGeo
} // namespace Dune

#endif // DUNE_CURVED_SURFACE_GRID_GEOMETRY_HH