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