1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_GEOMETRY_AFFINEGEOMETRY_HH
4 #define DUNE_GEOMETRY_AFFINEGEOMETRY_HH
5 
6 /** \file
7  *  \brief An implementation of the Geometry interface for affine geometries
8  *  \author Martin Nolte
9  */
10 
11 #include <cmath>
12 
13 #include <dune/common/fmatrix.hh>
14 #include <dune/common/fvector.hh>
15 
16 #include <dune/geometry/type.hh>
17 
18 namespace Dune
19 {
20 
21   // External Forward Declarations
22   // -----------------------------
23 
24   namespace Geo
25   {
26 
27     template< typename Implementation >
28     class ReferenceElement;
29 
30     template< class ctype, int dim >
31     class ReferenceElementImplementation;
32 
33     template< class ctype, int dim >
34     struct ReferenceElements;
35 
36   }
37 
38 
39   namespace Impl
40   {
41 
42     // FieldMatrixHelper
43     // -----------------
44 
45     template< class ct >
46     struct FieldMatrixHelper
47     {
48       typedef ct ctype;
49 
50       template< int m, int n >
AxDune::Impl::FieldMatrixHelper51       static void Ax ( const FieldMatrix< ctype, m, n > &A, const FieldVector< ctype, n > &x, FieldVector< ctype, m > &ret )
52       {
53         for( int i = 0; i < m; ++i )
54         {
55           ret[ i ] = ctype( 0 );
56           for( int j = 0; j < n; ++j )
57             ret[ i ] += A[ i ][ j ] * x[ j ];
58         }
59       }
60 
61       template< int m, int n >
ATxDune::Impl::FieldMatrixHelper62       static void ATx ( const FieldMatrix< ctype, m, n > &A, const FieldVector< ctype, m > &x, FieldVector< ctype, n > &ret )
63       {
64         for( int i = 0; i < n; ++i )
65         {
66           ret[ i ] = ctype( 0 );
67           for( int j = 0; j < m; ++j )
68             ret[ i ] += A[ j ][ i ] * x[ j ];
69         }
70       }
71 
72       template< int m, int n, int p >
ABDune::Impl::FieldMatrixHelper73       static void AB ( const FieldMatrix< ctype, m, n > &A, const FieldMatrix< ctype, n, p > &B, FieldMatrix< ctype, m, p > &ret )
74       {
75         for( int i = 0; i < m; ++i )
76         {
77           for( int j = 0; j < p; ++j )
78           {
79             ret[ i ][ j ] = ctype( 0 );
80             for( int k = 0; k < n; ++k )
81               ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
82           }
83         }
84       }
85 
86       template< int m, int n, int p >
ATBTDune::Impl::FieldMatrixHelper87       static void ATBT ( const FieldMatrix< ctype, m, n > &A, const FieldMatrix< ctype, p, m > &B, FieldMatrix< ctype, n, p > &ret )
88       {
89         for( int i = 0; i < n; ++i )
90         {
91           for( int j = 0; j < p; ++j )
92           {
93             ret[ i ][ j ] = ctype( 0 );
94             for( int k = 0; k < m; ++k )
95               ret[ i ][ j ] += A[ k ][ i ] * B[ j ][ k ];
96           }
97         }
98       }
99 
100       template< int m, int n >
ATA_LDune::Impl::FieldMatrixHelper101       static void ATA_L ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, n > &ret )
102       {
103         for( int i = 0; i < n; ++i )
104         {
105           for( int j = 0; j <= i; ++j )
106           {
107             ret[ i ][ j ] = ctype( 0 );
108             for( int k = 0; k < m; ++k )
109               ret[ i ][ j ] += A[ k ][ i ] * A[ k ][ j ];
110           }
111         }
112       }
113 
114       template< int m, int n >
ATADune::Impl::FieldMatrixHelper115       static void ATA ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, n > &ret )
116       {
117         for( int i = 0; i < n; ++i )
118         {
119           for( int j = 0; j <= i; ++j )
120           {
121             ret[ i ][ j ] = ctype( 0 );
122             for( int k = 0; k < m; ++k )
123               ret[ i ][ j ] += A[ k ][ i ] * A[ k ][ j ];
124             ret[ j ][ i ] = ret[ i ][ j ];
125           }
126 
127           ret[ i ][ i ] = ctype( 0 );
128           for( int k = 0; k < m; ++k )
129             ret[ i ][ i ] += A[ k ][ i ] * A[ k ][ i ];
130         }
131       }
132 
133       template< int m, int n >
AAT_LDune::Impl::FieldMatrixHelper134       static void AAT_L ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, m, m > &ret )
135       {
136         /*
137            if (m==2) {
138            ret[0][0] = A[0]*A[0];
139            ret[1][1] = A[1]*A[1];
140            ret[1][0] = A[0]*A[1];
141            }
142            else
143          */
144         for( int i = 0; i < m; ++i )
145         {
146           for( int j = 0; j <= i; ++j )
147           {
148             ctype &retij = ret[ i ][ j ];
149             retij = A[ i ][ 0 ] * A[ j ][ 0 ];
150             for( int k = 1; k < n; ++k )
151               retij += A[ i ][ k ] * A[ j ][ k ];
152           }
153         }
154       }
155 
156       template< int m, int n >
AATDune::Impl::FieldMatrixHelper157       static void AAT ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, m, m > &ret )
158       {
159         for( int i = 0; i < m; ++i )
160         {
161           for( int j = 0; j < i; ++j )
162           {
163             ret[ i ][ j ] = ctype( 0 );
164             for( int k = 0; k < n; ++k )
165               ret[ i ][ j ] += A[ i ][ k ] * A[ j ][ k ];
166             ret[ j ][ i ] = ret[ i ][ j ];
167           }
168           ret[ i ][ i ] = ctype( 0 );
169           for( int k = 0; k < n; ++k )
170             ret[ i ][ i ] += A[ i ][ k ] * A[ i ][ k ];
171         }
172       }
173 
174       template< int n >
LxDune::Impl::FieldMatrixHelper175       static void Lx ( const FieldMatrix< ctype, n, n > &L, const FieldVector< ctype, n > &x, FieldVector< ctype, n > &ret )
176       {
177         for( int i = 0; i < n; ++i )
178         {
179           ret[ i ] = ctype( 0 );
180           for( int j = 0; j <= i; ++j )
181             ret[ i ] += L[ i ][ j ] * x[ j ];
182         }
183       }
184 
185       template< int n >
LTxDune::Impl::FieldMatrixHelper186       static void LTx ( const FieldMatrix< ctype, n, n > &L, const FieldVector< ctype, n > &x, FieldVector< ctype, n > &ret )
187       {
188         for( int i = 0; i < n; ++i )
189         {
190           ret[ i ] = ctype( 0 );
191           for( int j = i; j < n; ++j )
192             ret[ i ] += L[ j ][ i ] * x[ j ];
193         }
194       }
195 
196       template< int n >
LTLDune::Impl::FieldMatrixHelper197       static void LTL ( const FieldMatrix< ctype, n, n > &L, FieldMatrix< ctype, n, n > &ret )
198       {
199         for( int i = 0; i < n; ++i )
200         {
201           for( int j = 0; j < i; ++j )
202           {
203             ret[ i ][ j ] = ctype( 0 );
204             for( int k = i; k < n; ++k )
205               ret[ i ][ j ] += L[ k ][ i ] * L[ k ][ j ];
206             ret[ j ][ i ] = ret[ i ][ j ];
207           }
208           ret[ i ][ i ] = ctype( 0 );
209           for( int k = i; k < n; ++k )
210             ret[ i ][ i ] += L[ k ][ i ] * L[ k ][ i ];
211         }
212       }
213 
214       template< int n >
LLTDune::Impl::FieldMatrixHelper215       static void LLT ( const FieldMatrix< ctype, n, n > &L, FieldMatrix< ctype, n, n > &ret )
216       {
217         for( int i = 0; i < n; ++i )
218         {
219           for( int j = 0; j < i; ++j )
220           {
221             ret[ i ][ j ] = ctype( 0 );
222             for( int k = 0; k <= j; ++k )
223               ret[ i ][ j ] += L[ i ][ k ] * L[ j ][ k ];
224             ret[ j ][ i ] = ret[ i ][ j ];
225           }
226           ret[ i ][ i ] = ctype( 0 );
227           for( int k = 0; k <= i; ++k )
228             ret[ i ][ i ] += L[ i ][ k ] * L[ i ][ k ];
229         }
230       }
231 
232       template< int n >
cholesky_LDune::Impl::FieldMatrixHelper233       static bool cholesky_L ( const FieldMatrix< ctype, n, n > &A, FieldMatrix< ctype, n, n > &ret, const bool checkSingular = false )
234       {
235         using std::sqrt;
236         for( int i = 0; i < n; ++i )
237         {
238           ctype &rii = ret[ i ][ i ];
239 
240           ctype xDiag = A[ i ][ i ];
241           for( int j = 0; j < i; ++j )
242             xDiag -= ret[ i ][ j ] * ret[ i ][ j ];
243 
244           // in some cases A can be singular, e.g. when checking local for
245           // outside points during checkInside
246           if( checkSingular && ! ( xDiag > ctype( 0 )) )
247             return false ;
248 
249           // otherwise this should be true always
250           assert( xDiag > ctype( 0 ) );
251           rii = sqrt( xDiag );
252 
253           ctype invrii = ctype( 1 ) / rii;
254           for( int k = i+1; k < n; ++k )
255           {
256             ctype x = A[ k ][ i ];
257             for( int j = 0; j < i; ++j )
258               x -= ret[ i ][ j ] * ret[ k ][ j ];
259             ret[ k ][ i ] = invrii * x;
260           }
261         }
262 
263         // return true for meaning A is non-singular
264         return true;
265       }
266 
267       template< int n >
detLDune::Impl::FieldMatrixHelper268       static ctype detL ( const FieldMatrix< ctype, n, n > &L )
269       {
270         ctype det( 1 );
271         for( int i = 0; i < n; ++i )
272           det *= L[ i ][ i ];
273         return det;
274       }
275 
276       template< int n >
invLDune::Impl::FieldMatrixHelper277       static ctype invL ( FieldMatrix< ctype, n, n > &L )
278       {
279         ctype det( 1 );
280         for( int i = 0; i < n; ++i )
281         {
282           ctype &lii = L[ i ][ i ];
283           det *= lii;
284           lii = ctype( 1 ) / lii;
285           for( int j = 0; j < i; ++j )
286           {
287             ctype &lij = L[ i ][ j ];
288             ctype x = lij * L[ j ][ j ];
289             for( int k = j+1; k < i; ++k )
290               x += L[ i ][ k ] * L[ k ][ j ];
291             lij = (-lii) * x;
292           }
293         }
294         return det;
295       }
296 
297       // calculates x := L^{-1} x
298       template< int n >
invLxDune::Impl::FieldMatrixHelper299       static void invLx ( FieldMatrix< ctype, n, n > &L, FieldVector< ctype, n > &x )
300       {
301         for( int i = 0; i < n; ++i )
302         {
303           for( int j = 0; j < i; ++j )
304             x[ i ] -= L[ i ][ j ] * x[ j ];
305           x[ i ] /= L[ i ][ i ];
306         }
307       }
308 
309       // calculates x := L^{-T} x
310       template< int n >
invLTxDune::Impl::FieldMatrixHelper311       static void invLTx ( FieldMatrix< ctype, n, n > &L, FieldVector< ctype, n > &x )
312       {
313         for( int i = n; i > 0; --i )
314         {
315           for( int j = i; j < n; ++j )
316             x[ i-1 ] -= L[ j ][ i-1 ] * x[ j ];
317           x[ i-1 ] /= L[ i-1 ][ i-1 ];
318         }
319       }
320 
321       template< int n >
spdDetADune::Impl::FieldMatrixHelper322       static ctype spdDetA ( const FieldMatrix< ctype, n, n > &A )
323       {
324         // return A[0][0]*A[1][1]-A[1][0]*A[1][0];
325         FieldMatrix< ctype, n, n > L;
326         cholesky_L( A, L );
327         return detL( L );
328       }
329 
330       template< int n >
spdInvADune::Impl::FieldMatrixHelper331       static ctype spdInvA ( FieldMatrix< ctype, n, n > &A )
332       {
333         FieldMatrix< ctype, n, n > L;
334         cholesky_L( A, L );
335         const ctype det = invL( L );
336         LTL( L, A );
337         return det;
338       }
339 
340       // calculate x := A^{-1} x
341       template< int n >
spdInvAxDune::Impl::FieldMatrixHelper342       static bool spdInvAx ( FieldMatrix< ctype, n, n > &A, FieldVector< ctype, n > &x, const bool checkSingular = false )
343       {
344         FieldMatrix< ctype, n, n > L;
345         const bool invertible = cholesky_L( A, L, checkSingular );
346         if( ! invertible ) return invertible ;
347         invLx( L, x );
348         invLTx( L, x );
349         return invertible;
350       }
351 
352       template< int m, int n >
detATADune::Impl::FieldMatrixHelper353       static ctype detATA ( const FieldMatrix< ctype, m, n > &A )
354       {
355         if( m >= n )
356         {
357           FieldMatrix< ctype, n, n > ata;
358           ATA_L( A, ata );
359           return spdDetA( ata );
360         }
361         else
362           return ctype( 0 );
363       }
364 
365       /** \brief Compute the square root of the determinant of A times A transposed
366        *
367        *  This is the volume element for an embedded submanifold and needed to
368        *  implement the method integrationElement().
369        */
370       template< int m, int n >
sqrtDetAATDune::Impl::FieldMatrixHelper371       static ctype sqrtDetAAT ( const FieldMatrix< ctype, m, n > &A )
372       {
373         using std::abs;
374         using std::sqrt;
375         // These special cases are here not only for speed reasons:
376         // The general implementation aborts if the matrix is almost singular,
377         // and the special implementation provide a stable way to handle that case.
378         if( (n == 2) && (m == 2) )
379         {
380           // Special implementation for 2x2 matrices: faster and more stable
381           return abs( A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ] );
382         }
383         else if( (n == 3) && (m == 3) )
384         {
385           // Special implementation for 3x3 matrices
386           const ctype v0 = A[ 0 ][ 1 ] * A[ 1 ][ 2 ] - A[ 1 ][ 1 ] * A[ 0 ][ 2 ];
387           const ctype v1 = A[ 0 ][ 2 ] * A[ 1 ][ 0 ] - A[ 1 ][ 2 ] * A[ 0 ][ 0 ];
388           const ctype v2 = A[ 0 ][ 0 ] * A[ 1 ][ 1 ] - A[ 1 ][ 0 ] * A[ 0 ][ 1 ];
389           return abs( v0 * A[ 2 ][ 0 ] + v1 * A[ 2 ][ 1 ] + v2 * A[ 2 ][ 2 ] );
390         }
391         else if ( (n == 3) && (m == 2) )
392         {
393           // Special implementation for 2x3 matrices
394           const ctype v0 = A[ 0 ][ 0 ] * A[ 1 ][ 1 ] - A[ 0 ][ 1 ] * A[ 1 ][ 0 ];
395           const ctype v1 = A[ 0 ][ 0 ] * A[ 1 ][ 2 ] - A[ 1 ][ 0 ] * A[ 0 ][ 2 ];
396           const ctype v2 = A[ 0 ][ 1 ] * A[ 1 ][ 2 ] - A[ 0 ][ 2 ] * A[ 1 ][ 1 ];
397           return sqrt( v0*v0 + v1*v1 + v2*v2);
398         }
399         else if( n >= m )
400         {
401           // General case
402           FieldMatrix< ctype, m, m > aat;
403           AAT_L( A, aat );
404           return spdDetA( aat );
405         }
406         else
407           return ctype( 0 );
408       }
409 
410       // A^{-1}_L = (A^T A)^{-1} A^T
411       // => A^{-1}_L A = I
412       template< int m, int n >
leftInvADune::Impl::FieldMatrixHelper413       static ctype leftInvA ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, m > &ret )
414       {
415         static_assert((m >= n), "Matrix has no left inverse.");
416         FieldMatrix< ctype, n, n > ata;
417         ATA_L( A, ata );
418         const ctype det = spdInvA( ata );
419         ATBT( ata, A, ret );
420         return det;
421       }
422 
423       template< int m, int n >
leftInvAxDune::Impl::FieldMatrixHelper424       static void leftInvAx ( const FieldMatrix< ctype, m, n > &A, const FieldVector< ctype, m > &x, FieldVector< ctype, n > &y )
425       {
426         static_assert((m >= n), "Matrix has no left inverse.");
427         FieldMatrix< ctype, n, n > ata;
428         ATx( A, x, y );
429         ATA_L( A, ata );
430         spdInvAx( ata, y );
431       }
432 
433       /** \brief Compute right pseudo-inverse of matrix A */
434       template< int m, int n >
rightInvADune::Impl::FieldMatrixHelper435       static ctype rightInvA ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, m > &ret )
436       {
437         static_assert((n >= m), "Matrix has no right inverse.");
438         using std::abs;
439         if( (n == 2) && (m == 2) )
440         {
441           const ctype det = (A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ]);
442           const ctype detInv = ctype( 1 ) / det;
443           ret[ 0 ][ 0 ] = A[ 1 ][ 1 ] * detInv;
444           ret[ 1 ][ 1 ] = A[ 0 ][ 0 ] * detInv;
445           ret[ 1 ][ 0 ] = -A[ 1 ][ 0 ] * detInv;
446           ret[ 0 ][ 1 ] = -A[ 0 ][ 1 ] * detInv;
447           return abs( det );
448         }
449         else
450         {
451           FieldMatrix< ctype, m , m > aat;
452           AAT_L( A, aat );
453           const ctype det = spdInvA( aat );
454           ATBT( A , aat , ret );
455           return det;
456         }
457       }
458 
459       template< int m, int n >
xTRightInvADune::Impl::FieldMatrixHelper460       static bool xTRightInvA ( const FieldMatrix< ctype, m, n > &A, const FieldVector< ctype, n > &x, FieldVector< ctype, m > &y )
461       {
462         static_assert((n >= m), "Matrix has no right inverse.");
463         FieldMatrix< ctype, m, m > aat;
464         Ax( A, x, y );
465         AAT_L( A, aat );
466         // check whether aat is singular and return true if non-singular
467         return spdInvAx( aat, y, true );
468       }
469     };
470 
471   } // namespace Impl
472 
473 
474 
475   /** \brief Implementation of the Geometry interface for affine geometries
476    * \tparam ct Type used for coordinates
477    * \tparam mydim Dimension of the geometry
478    * \tparam cdim Dimension of the world space
479    */
480   template< class ct, int mydim, int cdim>
481   class AffineGeometry
482   {
483   public:
484 
485     /** \brief Type used for coordinates */
486     typedef ct ctype;
487 
488     /** \brief Dimension of the geometry */
489     static const int mydimension= mydim;
490 
491     /** \brief Dimension of the world space */
492     static const int coorddimension = cdim;
493 
494     /** \brief Type for local coordinate vector */
495     typedef FieldVector< ctype, mydimension > LocalCoordinate;
496 
497     /** \brief Type for coordinate vector in world space */
498     typedef FieldVector< ctype, coorddimension > GlobalCoordinate;
499 
500     /** \brief Type used for volume */
501     typedef ctype Volume;
502 
503     /** \brief Type for the transposed Jacobian matrix */
504     typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
505 
506     /** \brief Type for the transposed inverse Jacobian matrix */
507     typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;
508 
509   private:
510     //! type of reference element
511     typedef Geo::ReferenceElement< Geo::ReferenceElementImplementation< ctype, mydimension > > ReferenceElement;
512 
513     typedef Geo::ReferenceElements< ctype, mydimension > ReferenceElements;
514 
515     // Helper class to compute a matrix pseudo inverse
516     typedef Impl::FieldMatrixHelper< ct > MatrixHelper;
517 
518   public:
519     /** \brief Create affine geometry from reference element, one vertex, and the Jacobian matrix */
AffineGeometry(const ReferenceElement & refElement,const GlobalCoordinate & origin,const JacobianTransposed & jt)520     AffineGeometry ( const ReferenceElement &refElement, const GlobalCoordinate &origin,
521                      const JacobianTransposed &jt )
522       : refElement_(refElement), origin_(origin), jacobianTransposed_(jt)
523     {
524       integrationElement_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jacobianTransposed_, jacobianInverseTransposed_ );
525     }
526 
527     /** \brief Create affine geometry from GeometryType, one vertex, and the Jacobian matrix */
AffineGeometry(Dune::GeometryType gt,const GlobalCoordinate & origin,const JacobianTransposed & jt)528     AffineGeometry ( Dune::GeometryType gt, const GlobalCoordinate &origin,
529                      const JacobianTransposed &jt )
530       : AffineGeometry(ReferenceElements::general( gt ), origin, jt)
531     { }
532 
533     /** \brief Create affine geometry from reference element and a vector of vertex coordinates */
534     template< class CoordVector >
AffineGeometry(const ReferenceElement & refElement,const CoordVector & coordVector)535     AffineGeometry ( const ReferenceElement &refElement, const CoordVector &coordVector )
536       : refElement_(refElement), origin_(coordVector[0])
537     {
538       for( int i = 0; i < mydimension; ++i )
539         jacobianTransposed_[ i ] = coordVector[ i+1 ] - origin_;
540       integrationElement_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jacobianTransposed_, jacobianInverseTransposed_ );
541     }
542 
543     /** \brief Create affine geometry from GeometryType and a vector of vertex coordinates */
544     template< class CoordVector >
AffineGeometry(Dune::GeometryType gt,const CoordVector & coordVector)545     AffineGeometry ( Dune::GeometryType gt, const CoordVector &coordVector )
546       : AffineGeometry(ReferenceElements::general( gt ), coordVector)
547     { }
548 
549     /** \brief Always true: this is an affine geometry */
affine() const550     bool affine () const { return true; }
551 
552     /** \brief Obtain the type of the reference element */
type() const553     Dune::GeometryType type () const { return refElement_.type(); }
554 
555     /** \brief Obtain number of corners of the corresponding reference element */
corners() const556     int corners () const { return refElement_.size( mydimension ); }
557 
558     /** \brief Obtain coordinates of the i-th corner */
corner(int i) const559     GlobalCoordinate corner ( int i ) const
560     {
561       return global( refElement_.position( i, mydimension ) );
562     }
563 
564     /** \brief Obtain the centroid of the mapping's image */
center() const565     GlobalCoordinate center () const { return global( refElement_.position( 0, 0 ) ); }
566 
567     /** \brief Evaluate the mapping
568      *
569      *  \param[in]  local  local coordinate to map
570      *
571      *  \returns corresponding global coordinate
572      */
global(const LocalCoordinate & local) const573     GlobalCoordinate global ( const LocalCoordinate &local ) const
574     {
575       GlobalCoordinate global( origin_ );
576       jacobianTransposed_.umtv( local, global );
577       return global;
578     }
579 
580     /** \brief Evaluate the inverse mapping
581      *
582      *  \param[in]  global  global coordinate to map
583      *
584      *  \return corresponding local coordinate
585      *
586      *  The returned local coordinate y minimizes
587      *  \code
588      *  (global( y ) - x).two_norm()
589      *  \endcode
590      *  on the entire affine hull of the reference element.  This degenerates
591      *  to the inverse map if the argument y is in the range of the map.
592      */
local(const GlobalCoordinate & global) const593     LocalCoordinate local ( const GlobalCoordinate &global ) const
594     {
595       LocalCoordinate local;
596       jacobianInverseTransposed_.mtv( global - origin_, local );
597       return local;
598     }
599 
600     /** \brief Obtain the integration element
601      *
602      *  If the Jacobian of the mapping is denoted by $J(x)$, the integration
603      *  integration element \f$\mu(x)\f$ is given by
604      *  \f[ \mu(x) = \sqrt{|\det (J^T(x) J(x))|}.\f]
605      *
606      *  \param[in]  local  local coordinate to evaluate the integration element in
607      *
608      *  \returns the integration element \f$\mu(x)\f$.
609      */
integrationElement(const LocalCoordinate & local) const610     ctype integrationElement ([[maybe_unused]] const LocalCoordinate &local) const
611     {
612       return integrationElement_;
613     }
614 
615     /** \brief Obtain the volume of the element */
volume() const616     Volume volume () const
617     {
618       return integrationElement_ * refElement_.volume();
619     }
620 
621     /** \brief Obtain the transposed of the Jacobian
622      *
623      *  \param[in]  local  local coordinate to evaluate Jacobian in
624      *
625      *  \returns a reference to the transposed of the Jacobian
626      */
jacobianTransposed(const LocalCoordinate & local) const627     const JacobianTransposed &jacobianTransposed ([[maybe_unused]] const LocalCoordinate &local) const
628     {
629       return jacobianTransposed_;
630     }
631 
632     /** \brief Obtain the transposed of the Jacobian's inverse
633      *
634      *  The Jacobian's inverse is defined as a pseudo-inverse. If we denote
635      *  the Jacobian by \f$J(x)\f$, the following condition holds:
636      *  \f[J^{-1}(x) J(x) = I.\f]
637      */
jacobianInverseTransposed(const LocalCoordinate & local) const638     const JacobianInverseTransposed &jacobianInverseTransposed ([[maybe_unused]] const LocalCoordinate &local) const
639     {
640       return jacobianInverseTransposed_;
641     }
642 
referenceElement(const AffineGeometry & geometry)643     friend ReferenceElement referenceElement ( const AffineGeometry &geometry )
644     {
645       return geometry.refElement_;
646     }
647 
648   private:
649     ReferenceElement refElement_;
650     GlobalCoordinate origin_;
651     JacobianTransposed jacobianTransposed_;
652     JacobianInverseTransposed jacobianInverseTransposed_;
653     ctype integrationElement_;
654   };
655 
656 } // namespace Dune
657 
658 #endif // #ifndef DUNE_GEOMETRY_AFFINEGEOMETRY_HH
659