geometry.hh 21.5 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
20
#include <dune/curvedgeometry/curvedgeometry.hh>

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
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
90
91
92
93
94
#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>

namespace Dune
{
  namespace CGeo
  {
    namespace Impl
    {
      struct IdentityMatrix {};

      template <class K, int n, int m, int l>
      FieldMatrix<K,n,m> ABt(const FieldMatrix<K,n,l>& A, const FieldMatrix<K,m,l>& B)
      {
        FieldMatrix<K,n,m> ABt;
        for (int i = 0; i < n; ++i)
          for (int j = 0; j < m; ++j) {
            ABt[i][j] = 0;
            for (int k = 0; k < l; ++k)
              ABt[i][j] += A[i][k] * B[j][k];
          }

        return ABt;
      }

      template <class K, int n, int m>
      FieldMatrix<K,n,m> ABt(const DiagonalMatrix<K,n>& A, const FieldMatrix<K,m,n>& B)
      {
        FieldMatrix<K,n,m> ABt;
        for (int i = 0; i < n; ++i)
          for (int j = 0; j < m; ++j) {
            ABt[i][j] = A[i][i] * B[j][i];
          }

        return ABt;
      }


      template <class K, int n, int m, int l>
      FieldMatrix<K,n,m> AB(const FieldMatrix<K,n,l>& A, const FieldMatrix<K,l,m>& B)
      {
        FieldMatrix<K,n,m> ABt;
        for (int i = 0; i < n; ++i)
          for (int j = 0; j < m; ++j) {
            ABt[i][j] = 0;
            for (int k = 0; k < l; ++k)
              ABt[i][j] += A[i][k] * B[k][j];
          }

        return ABt;
      }

      template <class K, int n, int m>
      FieldMatrix<K,n,m> AB(const DiagonalMatrix<K,n>& A, const FieldMatrix<K,n,m>& B)
      {
        FieldMatrix<K,n,m> ABt;
        for (int i = 0; i < n; ++i)
          for (int j = 0; j < m; ++j) {
            ABt[i][j] = A[i][i] * B[i][j];
          }

        return ABt;
      }

      template <class K, int n, int m>
      const FieldMatrix<K,n,m>& AB(const IdentityMatrix& A, const FieldMatrix<K,n,m>& B)
      {
        return B;
      }


95
      // efficient implementation of a trivial local geometry for elements
96
97
98
99
100
101
102
103
104
105
106
107
108
109
      struct DefaultLocalGeometry
      {
        template <class LocalCoordinate>
        decltype(auto) global(LocalCoordinate&& local) const
        {
          return std::forward<LocalCoordinate>(local);
        }

        template <class LocalCoordinate>
        IdentityMatrix jacobianTransposed(LocalCoordinate&& local) const
        {
          return {};
        }
      };
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133

      // type-erased wrapper for local geometries
      template <class LG>
      struct LocalGeometryInterface
      {
        using LocalCoordinate = typename LG::LocalCoordinate;
        using GlobalCoordinate = typename LG::GlobalCoordinate;
        using JacobianTransposed = typename LG::JacobianTransposed;

        template <class LocalGeometry>
        LocalGeometryInterface(const LocalGeometry& lg)
          : global([lg](auto const& local) { return lg.global(local); })
          , jacobianTransposed([lg](auto const& local) { return lg.jacobianTransposed(local); })
        {}

        template <class GlobalFct, class JacobianTransposedFct>
        LocalGeometryInterface(GlobalFct&& g, JacobianTransposedFct&& jt)
          : global(std::forward<GlobalFct>(g))
          , jacobianTransposed(std::forward<JacobianTransposedFct>(jt))
        {}

        std::function<GlobalCoordinate(LocalCoordinate)> global;
        std::function<JacobianTransposed(LocalCoordinate)> jacobianTransposed;
      };
134
135
136
137
    }

    /// \brief Curved geometry implementation based on local basis function parametrization
    /**
138
     *  Parametrization of the geometry by a differentiable localFunction
139
140
141
142
     *
     *  \tparam  ct      coordinate type
     *  \tparam  mydim   geometry dimension
     *  \tparam  cdim    coordinate dimension
143
144
145
146
     *  \tparam  LF      localFunction parametrizing the geometry
     *  \tparam  LG      localGeometry for coordinate transform from subEntity to element,
     *                   see \ref Impl::DefaultLocalGeometry and \ref Impl::LocalGeometryInterface
     *  \tparam  order   Polynomial order of lagrange parametrization. If order < 0 use localFunction.
147
     */
148
    template <class ct, int mydim, int cdim, class LF, class LG, int order>
149
    class Geometry
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
        : public Dune::LagrangeCurvedGeometry<ct, mydim, cdim, (order > 0 ? order : 1)>
    {
      using Super = Dune::LagrangeCurvedGeometry<ct, mydim, cdim, (order > 0 ? order : 1)>;

      //! 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)...)
      {}
    };


    template <class ct, int mydim, int cdim, class LF, class LG>
    class Geometry<ct,mydim,cdim,LF,LG,-1>
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
226
227
228
229
230
231
    {
    public:
      /// coordinate type
      using ctype = ct;

      /// geometry dimension
      static const int mydimension = mydim;

      /// coordinate dimension
      static const int coorddimension = cdim;

      /// type of local coordinates
      using LocalCoordinate = FieldVector<ctype, mydimension>;

      /// type of global coordinates
      using GlobalCoordinate = FieldVector<ctype, coorddimension>;

      /// type of volume
      using Volume = ctype;

      /// type of jacobian transposed
      using JacobianTransposed = FieldMatrix<ctype, mydimension, coorddimension>;

      /// type of jacobian inverse transposed
      using JacobianInverseTransposed = FieldMatrix<ctype, coorddimension, mydimension>;

    public:
      /// type of reference element
      using ReferenceElements = Dune::ReferenceElements<ctype, mydimension>;
      using ReferenceElement = typename ReferenceElements::ReferenceElement;

    protected:
      //! helper structure containing some matrix routines. See affinegeometry.hh
      using MatrixHelper = Dune::Impl::FieldMatrixHelper<ct>;

232
233
234
235
236
237
      //! type of coordinate transformation for subEntities to codim=0 entities
      using LocalGeometry = LG;

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

238
      //! type of the LocalFunction representing the derivative
239
      using DerivativeLocalFunction = std::decay_t<decltype(derivative(std::declval<LF>()))>;
240
241
242
243
244
245
246
247

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

      //! maximal number of Newton iteration in `geometry.local(global)`
      static int maxIteration () { return 100; }

    public:
248
      /// \brief constructor from localFunction to parametrize the geometry
249
      /**
250
251
252
253
       *  \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
254
255
       *
       */
256
257
      template <class... Args>
      Geometry (const ReferenceElement& refElement, const LocalFunction& localFunction, Args&&... args)
258
259
        : refElement_(refElement)
        , localFunction_(localFunction)
260
        , localGeometry_(std::forward<Args>(args)...)
261
262
263
264
265
      {}

      /// \brief constructor, forwarding to the other constructor that take a reference-element
      /**
       *  \param[in]  gt             geometry type
266
267
268
       *  \param[in]  localFunction  localFunction for the parametrization of the geometry
       *                             (stored by value)
       *  \param[in]  args...        argument to construct the local geometry
269
       */
270
271
272
      template <class... Args>
      Geometry (Dune::GeometryType gt, const LocalFunction& localFunction, Args&&... args)
        : Geometry(ReferenceElements::general(gt), localFunction, std::forward<Args>(args)...)
273
274
275
276
277
278
279
280
      {}

      /// \brief is this mapping affine? No!
      bool affine () const
      {
        return false;
      }

281
      /// \brief obtain the element type from the reference element
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
      Dune::GeometryType type () const
      {
        return refElement_.type();
      }

      /// \brief obtain number of corners of the corresponding reference element
      int corners () const
      {
        return refElement_.size(mydimension);
      }

      /// \brief obtain coordinates of the i-th corner
      GlobalCoordinate corner (int i) const
      {
        assert( (i >= 0) && (i < corners()) );
        return global(refElement_.position(i, mydimension));
      }

      /// \brief obtain the centroid of the mapping's image
      GlobalCoordinate center () const
      {
        return global(refElement_.position(0, 0));
      }

      /// \brief Evaluate the coordinate mapping
      /**
       *  Evaluate the local function in local coordinates
       *
       *  \param[in] local  local coordinates
       *
       *  \returns corresponding global coordinate
       */
      GlobalCoordinate global (const LocalCoordinate& local) const
      {
        return localFunction()(localGeometry().global(local));
      }

      /// \brief Evaluate the inverse coordinate mapping
      /**
       *  \param[in] globalCoord global coordinate to map
       *
       *  \return corresponding local coordinate
324
325
326
       *  \throws in case of an error indicating that the local coordinate can not be obtained,
       *          an exception is thrown. See \ref checkedLocal for a variant that returns
       *          an optional instead.
327
328
329
330
331
332
333
334
335
336
337
       *
       *  \note For given global coordinate `y` the returned local coordinate `x` that minimizes
       *  the following function over the local coordinate space spanned by the reference element.
       *  \code
       *  (global( x ) - y).two_norm()
       *  \endcode
       */
      LocalCoordinate local (const GlobalCoordinate& globalCoord) const
      {
        auto localCoord = checkedLocal(globalCoord);
        if (!localCoord)
338
339
          DUNE_THROW(Dune::Exception,
            "local coordinate can not be recovered from given global coordinate " << globalCoord);
340
341
342
343

        return *localCoord;
      }

344
345
346
347
348
349
350
      /// \brief Evaluate the inverse coordinate mapping
      /**
       *  \param[in] globalCoord global coordinate to map
       *
       *  \return corresponding local coordinate or nothing, in case the local coordinate
       *          could not be calculated. The return value is wrapped in an optional.
       **/
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
      Std::optional<LocalCoordinate> checkedLocal (const GlobalCoordinate& globalCoord) const
      {
        LocalCoordinate x = flatGeometry().local(globalCoord);
        LocalCoordinate dx;

        for (int i = 0; i < maxIteration(); ++i)
        {
          // Newton's method: DF^n dx^n = F^n, x^{n+1} -= dx^n
          const GlobalCoordinate dglobal = global(x) - globalCoord;
          const bool invertible = MatrixHelper::xTRightInvA(jacobianTransposed(x), dglobal, dx);

          // break if jacobian is not invertible
          if (!invertible)
            return {};

          // update x with correction
          x -= dx;

          // break if tolerance is reached.
          if (dx.two_norm2() < tolerance())
            break;
        }

        if (dx.two_norm2() > tolerance())
          return {};

        return x;
      }

      /// \brief Construct a global coordinate normal of the curvilinear element evaluated at
      /// a given local coordinate
      /**
       * \note Implemented for codim=1 entities only, i.e. edges in 2D and faces in 3D
       **/
      GlobalCoordinate normal (const LocalCoordinate& local) const
      {
        if ((mydim == 1) && (cdim == 2)) { return normal1D(local, jacobianTransposed(local)); }
        if ((mydim == 2) && (cdim == 3)) { return normal2D(local, jacobianTransposed(local)); }

        DUNE_THROW(Dune::NotImplemented,
          "ERROR: normal() method only defined for edges in 2D and faces in 3D");

        return GlobalCoordinate(0);
      }

      ///  \brief Obtain the integration element
      /**
       *  If the Jacobian of the mapping is denoted by $J(x)$, the integration
       *  integration element \f$\mu(x)\f$ is given by
       *  \f[ \mu(x) = \sqrt{|\det (J^T(x) J(x))|}.\f]
       *
       *  \param[in]  local  local coordinate to evaluate the integration element in
       *
       *  \returns the integration element \f$\mu(x)\f$.
       */
      ctype integrationElement (const LocalCoordinate& local) const
      {
        return MatrixHelper::sqrtDetAAT(jacobianTransposed(local));
      }

      /// \brief Obtain the volume of the mapping's image
      /**
       * Calculates the volume of the entity by numerical integration
       *
       * \todo Check the quadrature order of the volume integral.
       */
      Volume volume () const
      {
        const int p = 10; // TODO: use adaptive quadrature rule
        const auto& quadRule = QuadratureRules<ctype, mydimension>::rule(type(), p);
        Volume vol(0);
        for (auto const& qp : quadRule)
          vol += integrationElement(qp.position()) * qp.weight();

        return vol;
      }

      /// \brief Obtain the transposed of the Jacobian
      /**
       *  \param[in]  local  local coordinate to evaluate Jacobian in
       *
       *  \returns the matrix corresponding to the transposed of the Jacobian
       */
      JacobianTransposed jacobianTransposed (const LocalCoordinate& local) const
      {
        if (!derivativeLocalFunction_) {
          derivativeLocalFunction_.emplace(derivative(localFunction()));
          derivativeLocalFunction_->bind(localFunction().localContext());
        }

441
442
443
444
445
446
447
        // coordinate in the localContext of the localFunction
        auto&& elementLocal = localGeometry().global(local);

        auto jLocal = localGeometry().jacobianTransposed(local);
        auto jGlobal = hostGeometry().jacobianTransposed(elementLocal);
        auto jacobian = (*derivativeLocalFunction_)(elementLocal);
        return Impl::AB(jLocal,Impl::ABt(jGlobal,jacobian));
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
      }

      /// \brief obtain the transposed of the Jacobian's inverse
      /**
       *  The Jacobian's inverse is defined as a pseudo-inverse. If we denote
       *  the Jacobian by \f$J(x)\f$, the following condition holds:
       *  \f[J^{-1}(x) J(x) = I.\f]
       */
      JacobianInverseTransposed jacobianInverseTransposed (const LocalCoordinate& local) const
      {
        JacobianInverseTransposed out;
        MatrixHelper::rightInvA(jacobianTransposed(local), out);
        return out;
      }

      /// \brief obtain the reference-element related to this geometry
      friend ReferenceElement referenceElement (const Geometry& geometry)
      {
        return geometry.refElement();
      }

    protected:
      // the internal stored reference element
      const ReferenceElement& refElement () const
      {
        return refElement_;
      }

      // normal vector to an edge line-element
      GlobalCoordinate normal1D (const LocalCoordinate& /*local*/, const JacobianTransposed& J) const
      {
        GlobalCoordinate res{
           J[0][1],
          -J[0][0]};

        return res /= res.two_norm();
      }

      // normal vector to a triangle or quad face element
      GlobalCoordinate normal2D (const LocalCoordinate& /*local*/, const JacobianTransposed& J) const
      {
        GlobalCoordinate res{
          J[0][1] * J[1][2] - J[0][2] * J[1][1],
          J[0][2] * J[1][0] - J[0][0] * J[1][2],
          J[0][0] * J[1][1] - J[0][1] * J[1][0]};

        return res /= res.two_norm();
      }

497
    private:
498
499
500
501
502
503
504
      const LocalFunction& localFunction () const
      {
        return localFunction_;
      }

      const LocalGeometry& localGeometry () const
      {
505
        return localGeometry_;
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
      }

      using FlatGeometry = MultiLinearGeometry<ctype, mydim, cdim>;
      const FlatGeometry& flatGeometry () const
      {
        if (!flatGeometry_) {
          std::vector<GlobalCoordinate> corners;
          corners.reserve(refElement_.size(mydimension));
          for (int i = 0; i < refElement_.size(mydimension); ++i)
            corners.push_back(global(refElement_.position(i, mydimension)));

          flatGeometry_.emplace(refElement_, std::move(corners));
        }

        return *flatGeometry_;
      }

523
      using HostGeometry = std::decay_t<decltype(std::declval<LF>().localContext().geometry())>;
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
      const HostGeometry& hostGeometry () const
      {
        if (!hostGeometry_)
          hostGeometry_.emplace(localFunction().localContext().geometry());

        return *hostGeometry_;
      }

    private:
      //! Reference of the geometry
      ReferenceElement refElement_;

      //! local parametrization of the element
      LocalFunction localFunction_;

      //! transformation of local coordinates to element-local coordinates
540
      LocalGeometry localGeometry_;
541
542
543
544
545
546
547
548
549
550
551

      // some data optionally provided
      mutable Std::optional<DerivativeLocalFunction> derivativeLocalFunction_;
      mutable Std::optional<FlatGeometry> flatGeometry_;
      mutable Std::optional<HostGeometry> hostGeometry_;
    };


  #ifndef DOXYGEN

    // Specialization for vertex geometries
552
553
    template <class ct, int cdim, class LF, class LG>
    class VertexGeometry
554
555
556
557
        : public AffineGeometry<ct,0,cdim>
    {
      using Super = AffineGeometry<ct,0,cdim>;

558
559
560
      using LocalFunction = LF;
      using LocalGeometry = LG;

561
562
563
      //! tolerance to numerical algorithms
      static ct tolerance () { return ct(16) * std::numeric_limits<ct>::epsilon(); }

564
565
      struct Construct {};

566
567
568
569
570
571
572
573
574
575
576
577
    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:
578
      VertexGeometry (Construct, const ReferenceElement& refElement, const LocalFunction& localFunction,
579
                      const LocalGeometry& lg)
580
        : Super(refElement, std::vector<GlobalCoordinate>{localFunction(lg.global(refElement.position(0,0)))})
581
582
583
      {}

    public:
584
585
586
      template <class... Args>
      VertexGeometry (const ReferenceElement& refElement, const LocalFunction& lf, Args&&... args)
        : VertexGeometry(Construct{}, refElement, lf, LocalGeometry(std::forward<Args>(args)...))
587
588
      {}

589
590
591
592
      template <class... Args>
      VertexGeometry (Dune::GeometryType gt, const LocalFunction& lf, Args&&... args)
        : VertexGeometry(Construct{}, ReferenceElements::general(gt), lf,
                   LocalGeometry(std::forward<Args>(args)...))
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
      {}

      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);
      }

611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
    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>;
631
    public:
632
      using Super::Super;
633
634
635
636
637
638
639
640
    };

#endif // DOXYGEN

  } // end namespace CGeo
} // namespace Dune

#endif // DUNE_CURVED_SURFACE_GRID_GEOMETRY_HH