coordprovider.hh 13.3 KB
Newer Older
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
1
2
3
4
5
6
7
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_CRVSRF_COORDPROVIDER_HH
#define DUNE_CRVSRF_COORDPROVIDER_HH

#include <array>

8
9
10
#include <dune/localfunctions/lagrange/equidistantpoints.hh>
#include <dune/localfunctions/lagrange/lagrangelfecache.hh>

Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
namespace Dune
{

  namespace crvsrf
  {

    template< class Coordinate>
    Coordinate edgeCenter(Coordinate gc1, const Coordinate &gc2)
    {
      for(typename Coordinate::size_type i = 0, end = gc1.size(); i < end; ++i)
      {
        gc1[i] = (gc1[i] + gc2[i]) / 2.0;
      }
      return gc1;
    }

    template< class Coordinate>
    Coordinate edgeVertex(Coordinate gc1, const Coordinate &gc2, double factor)
    {
      for(typename Coordinate::size_type i = 0, end = gc1.size(); i < end; ++i)
      {
        gc1[i] = gc1[i] + (gc2[i] - gc1[i]) * factor;
      }
      return gc1;
    }

    template< class Coordinate>
    Coordinate faceVertex(Coordinate gc1, const Coordinate &gc2, const Coordinate &gc3,
                                            double coeff1, double coeff2, double coeff3)
    {
      for(typename Coordinate::size_type i = 0, end = gc1.size(); i < end; ++i)
      {
        gc1[i] = coeff1 * gc1[i] + coeff2 * gc2[i] + coeff3 * gc3[i];
      }
      return gc1;
    }

48
49


Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
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
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
    // InterpolatoryVerticesGenerator
    // ------------------------------

    template< class Coordinate, int mydim, int order >
    struct InterpolatoryVerticesGenerator;

    //Vertex
    template< class Coordinate, int order >
    struct InterpolatoryVerticesGenerator< Coordinate, 0, order >
    {
      static void generate(const std::array<Coordinate, 1> &corners, std::vector<Coordinate> &vertices)
      {
        vertices.resize(1);
        vertices[0] = corners[0];
      }
    };

    //Edge
    template< class Coordinate, int order >
    struct InterpolatoryVerticesGenerator< Coordinate, 1, order >
    {
      static void generate(const std::array<Coordinate, 2> &corners, std::vector<Coordinate> &vertices)
      {
        vertices.resize(order + 1);
        vertices[0] = corners[0];
        for(int i = 1; i < order; ++i)
        {
          vertices[i] = edgeVertex(corners[0], corners[1], (double)i / order);
        }
        vertices[order] = corners[1];
      }
    };

    //Triangle (order 1)
    template< class Coordinate >
    struct InterpolatoryVerticesGenerator< Coordinate, 2, 1 >
    {
      static void generate(const std::array<Coordinate, 4> &corners, std::vector<Coordinate> &vertices)
      {
        vertices.resize(3);
        vertices[0] = corners[0];
        vertices[1] = corners[1];
        vertices[2] = corners[2];
      }
    };

    //Triangle (order 2)
    template< class Coordinate >
    struct InterpolatoryVerticesGenerator< Coordinate, 2, 2 >
    {
      static void generate(const std::array<Coordinate, 4> &corners, std::vector<Coordinate> &vertices)
      {
        vertices.resize(6);
        vertices[0] = corners[0];
        vertices[1] = edgeCenter(corners[0], corners[1]);
        vertices[2] = corners[1];
        vertices[3] = edgeCenter(corners[0], corners[2]);
        vertices[4] = edgeCenter(corners[1], corners[2]);
        vertices[5] = corners[2];
      }
    };

    //Triangle (order 3)
    template< class Coordinate >
    struct InterpolatoryVerticesGenerator< Coordinate, 2, 3 >
    {
      static void generate(const std::array<Coordinate, 4> &corners, std::vector<Coordinate> &vertices)
      {
        static const double cf1 = 1.0 / 3.0;
        static const double cf2 = 2.0 / 3.0;
        vertices.resize(10);
        vertices[0] = corners[0];
        vertices[1] = edgeVertex(corners[0], corners[1], cf1);
        vertices[2] = edgeVertex(corners[0], corners[1], cf2);
        vertices[3] = corners[1];
        vertices[4] = edgeVertex(corners[0], corners[2], cf1);
        vertices[5] = faceVertex(corners[0], corners[1], corners[2], cf1, cf1, cf1);
        vertices[6] = edgeVertex(corners[1], corners[2], cf1);
        vertices[7] = edgeVertex(corners[0], corners[2], cf2);
        vertices[8] = edgeVertex(corners[1], corners[2], cf2);
        vertices[9] = corners[2];
      }
    };

    //Triangle (order 4)
    template< class Coordinate >
    struct InterpolatoryVerticesGenerator< Coordinate, 2, 4 >
    {
      static void generate(const std::array<Coordinate, 4> &corners, std::vector<Coordinate> &vertices)
      {
        static const double cf1 = 1.0 / 4.0;
        static const double cf2 = 2.0 / 4.0;
        static const double cf3 = 3.0 / 4.0;
        vertices.resize(15);
        vertices[0] = corners[0];
        vertices[1] = edgeVertex(corners[0], corners[1], cf1);
        vertices[2] = edgeCenter(corners[0], corners[1]);
        vertices[3] = edgeVertex(corners[0], corners[1], cf3);
        vertices[4] = corners[1];
        vertices[5] = edgeVertex(corners[0], corners[2], cf1);
        vertices[6] = faceVertex(corners[0], corners[1], corners[2], cf2, cf1, cf1);
        vertices[7] = faceVertex(corners[0], corners[1], corners[2], cf1, cf2, cf1);
        vertices[8] = edgeVertex(corners[1], corners[2], cf1);
        vertices[9] = edgeCenter(corners[0], corners[2]);
        vertices[10] = faceVertex(corners[0], corners[1], corners[2], cf1, cf1, cf2);
        vertices[11] = edgeCenter(corners[1], corners[2]);
        vertices[12] = edgeVertex(corners[0], corners[2], cf3);
        vertices[13] = edgeVertex(corners[1], corners[2], cf3);
        vertices[14] = corners[2];
      }
    };

    //Triangle (order 5)
    template< class Coordinate >
    struct InterpolatoryVerticesGenerator< Coordinate, 2, 5 >
    {
      static void generate(const std::array<Coordinate, 4> &corners, std::vector<Coordinate> &vertices)
      {
        static const double cf1 = 1.0 / 5.0;
        static const double cf2 = 2.0 / 5.0;
        static const double cf3 = 3.0 / 5.0;
        static const double cf4 = 4.0 / 5.0;
        vertices.resize(21);
        vertices[0] = corners[0];
        vertices[1] = edgeVertex(corners[0], corners[1], cf1);
        vertices[2] = edgeVertex(corners[0], corners[1], cf2);
        vertices[3] = edgeVertex(corners[0], corners[1], cf3);
        vertices[4] = edgeVertex(corners[0], corners[1], cf4);
        vertices[5] = corners[1];
        vertices[6] = edgeVertex(corners[0], corners[2], cf1);
        vertices[7] = faceVertex(corners[0], corners[1], corners[2], cf3, cf1, cf1);
        vertices[8] = faceVertex(corners[0], corners[1], corners[2], cf2, cf2, cf1);
        vertices[9] = faceVertex(corners[0], corners[1], corners[2], cf1, cf3, cf1);
        vertices[10] = edgeVertex(corners[1], corners[2], cf1);
        vertices[11] = edgeVertex(corners[0], corners[2], cf2);
        vertices[12] = faceVertex(corners[0], corners[1], corners[2], cf2, cf1, cf2);
        vertices[13] = faceVertex(corners[0], corners[1], corners[2], cf1, cf2, cf2);
        vertices[14] = edgeVertex(corners[1], corners[2], cf2);
        vertices[15] = edgeVertex(corners[0], corners[2], cf3);
        vertices[16] = faceVertex(corners[0], corners[1], corners[2], cf1, cf1, cf3);
        vertices[17] = edgeVertex(corners[1], corners[2], cf3);
        vertices[18] = edgeVertex(corners[0], corners[2], cf4);
        vertices[19] = edgeVertex(corners[1], corners[2], cf4);
        vertices[20] = corners[2];
      }
    };



    // CoordProvider (adapted CoordVector from GeometryGrid)
    // -----------------------------------------------------

    template< int mydim, class Grid, bool fake >
    class CoordProvider;


    //CoordProvider (real)
    template< int mydim, class Grid >
    class CoordProvider< mydim, Grid, false >
    {
      typedef typename std::remove_const< Grid >::type::Traits Traits;

      typedef typename Traits::ctype ctype;

      static const int dimension = Traits::dimension;
      static const int mydimension = mydim;
      static const int codimension = dimension - mydimension;
      static const int dimensionworld = Traits::dimensionworld;

      typedef FieldVector< ctype, dimensionworld > Coordinate;

      typedef typename Traits::HostGrid HostGrid;
      typedef typename Traits::CoordFunction CoordFunction;

      typedef typename HostGrid::template Codim< codimension >::Entity HostEntity;

226
227
228
      using LocalFECache = LagrangeLocalFiniteElementCache<ctype, ctype, mydimension, Grid::order>;
      using LocalFiniteElement = typename LocalFECache::FiniteElementType;

Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
229
230
231
232
    public:
      CoordProvider ( const HostEntity &hostEntity,
                    const CoordFunction &coordFunction )
        : coordFunction_(coordFunction),
233
234
235
          hostEntity_(hostEntity),
          localFECache_(),
          localFE_(localFECache_.get(hostEntity.type()))
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
236
237
      {}

238
      void calculate ( std::vector< Coordinate > &vertices ) const
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
239
      {
240
241
242
243
244
245
246
247
248
        vertices.resize(localFE_.size());

        auto f = [this,geo = hostEntity_.geometry()](auto const& local) -> Coordinate {
          Coordinate res;
          this->coordFunction_.evaluate(geo.global(local), res);
          return res;
        };

        localFE_.localInterpolation().interpolate(f, vertices);
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
249
250
251
252
253
      }

    private:
      const CoordFunction &coordFunction_;
      const HostEntity &hostEntity_;
254
255
256

      LocalFECache localFECache_;
      LocalFiniteElement const& localFE_;
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
    };

    //CoordProvider (fake)
    template< int mydim, class Grid >
    class CoordProvider< mydim, Grid, true >
    {
      typedef typename std::remove_const< Grid > :: type :: Traits Traits;

      typedef typename Traits::ctype ctype;

      static const int dimension = Traits::dimension;
      static const int mydimension = mydim;
      static const int codimension = dimension - mydimension;
      static const int dimensionworld = Traits::dimensionworld;

      typedef FieldVector< ctype, dimensionworld > Coordinate;

      typedef typename Traits::HostGrid HostGrid;
      typedef typename Traits::CoordFunction CoordFunction;

      typedef typename HostGrid::template Codim< 0 >::Entity HostElement;

    public:
      CoordProvider ( const HostElement &hostElement,
                    const unsigned int subEntity,
                    const CoordFunction &coordFunction )
        : coordFunction_(coordFunction),
284
          hostElement_(hostElement),
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
285
286
287
          subEntity_( subEntity )
      {}

288
      void calculate ( std::vector< Coordinate > &vertices ) const
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
289
      {
290
        throw "not implemented";
291
        const GeometryType type = hostElement_.geometry().type();
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
292
293
        auto refElement = referenceElement< ctype, dimension >( type );
        const std::size_t numCorners = refElement.size( subEntity_, codimension, dimension );
294
        std::array<Coordinate, (1 << mydim)> corners;
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
295
296
297
298
299
        for( std::size_t i = 0; i < numCorners; ++i )
        {
          const std::size_t j = refElement.subEntity( subEntity_, codimension, i, dimension );
          corners[i] = hostElement_.geometry().corner(j);
        }
300
301
        InterpolatoryVerticesGenerator<Coordinate, mydim, Grid::order>::generate(corners, vertices);
        for(auto &v : vertices) coordFunction_.evaluate(v, v);
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
302
303
304
305
      }

    private:
      const CoordFunction &coordFunction_;
306
      const HostElement &hostElement_;
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
      const unsigned int subEntity_;
    };



    // IntersectionCoordProvider
    // -----------------------

    template< class Grid >
    class IntersectionCoordProvider
    {
      typedef typename std::remove_const< Grid >::type::Traits Traits;

      typedef typename Traits::ctype ctype;

      static const int dimension = Traits::dimension;
      static const int codimension = 1;
      static const int mydimension = dimension-codimension;
      static const int dimensionworld = Traits::dimensionworld;

      typedef FieldVector< ctype, dimensionworld > Coordinate;

      typedef typename Traits::HostGrid HostGrid;
      typedef typename Traits::CoordFunction CoordFunction;

      typedef typename Traits::template Codim< 0 >::GeometryImpl ElementGeometryImpl;
      typedef typename Traits::template Codim< codimension >::LocalGeometry HostLocalGeometry;

    public:
      IntersectionCoordProvider ( const ElementGeometryImpl &elementGeometry,
337
                                  const HostLocalGeometry &hostLocalGeometry,
338
339
                                  const CoordFunction &coordFunction )
        : coordFunction_(coordFunction),
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
340
341
342
343
          elementGeometry_( elementGeometry ),
          hostLocalGeometry_( hostLocalGeometry )
      {}

344
      void calculate ( std::vector< Coordinate > &vertices ) const
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
345
      {
346
        throw "not implemented";
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
347
        const std::size_t numCorners = hostLocalGeometry_.corners();
348
349
        std::array<Coordinate, (1 << mydimension)> corners;
        // $flo: coordFunction is applied here which is in contrast to GeomentryGrid's behaviour!
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
350
351
        for( std::size_t i = 0; i < numCorners; ++i )
          corners[ i ] = elementGeometry_.global( hostLocalGeometry_.corner( i ) );
352
353
        InterpolatoryVerticesGenerator<Coordinate, mydimension, Grid::order>::generate(corners, vertices);
        for(auto &v : vertices) coordFunction_.evaluate(v, v);
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
354
355
356
357
358
359
360
361
362
      }

      template< unsigned int numCorners >
      void calculate ( Coordinate (&corners)[ numCorners ] ) const
      {
        assert( numCorners == hostLocalGeometry_.corners() );
      }

    private:
363
      const CoordFunction &coordFunction_;
Stenger, Florian's avatar
v0.1.0  
Stenger, Florian committed
364
365
366
367
368
369
370
371
372
      const ElementGeometryImpl &elementGeometry_;
      HostLocalGeometry hostLocalGeometry_;
    };

  } // namespace crvsrf

} // namespace Dune

#endif // #ifndef DUNE_CRVSRF_COORDPROVIDER_HH