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_MULTILINEARGEOMETRY_HH 4 #define DUNE_GEOMETRY_MULTILINEARGEOMETRY_HH 5 6 #include <cassert> 7 #include <functional> 8 #include <iterator> 9 #include <limits> 10 #include <vector> 11 12 #include <dune/common/fmatrix.hh> 13 #include <dune/common/fvector.hh> 14 #include <dune/common/typetraits.hh> 15 16 #include <dune/geometry/affinegeometry.hh> 17 #include <dune/geometry/referenceelements.hh> 18 #include <dune/geometry/type.hh> 19 20 namespace Dune 21 { 22 23 // MultiLinearGeometryTraits 24 // ------------------------- 25 26 /** \brief default traits class for MultiLinearGeometry 27 * 28 * The MultiLinearGeometry (and CachedMultiLinearGeometry) allow tweaking 29 * some implementation details through a traits class. 30 * 31 * This structure provides the default values. 32 * 33 * \tparam ct coordinate type 34 */ 35 template< class ct > 36 struct MultiLinearGeometryTraits 37 { 38 /** \brief helper structure containing some matrix routines 39 * 40 * This helper allows exchanging the matrix inversion algorithms. 41 * It must provide the following static methods: 42 * \code 43 * template< int m, int n > 44 * static ctype sqrtDetAAT ( const FieldMatrix< ctype, m, n > &A ); 45 * 46 * template< int m, int n > 47 * static ctype rightInvA ( const FieldMatrix< ctype, m, n > &A, 48 * FieldMatrix< ctype, n, m > &ret ); 49 * 50 * template< int m, int n > 51 * static void xTRightInvA ( const FieldMatrix< ctype, m, n > &A, 52 * const FieldVector< ctype, n > &x, 53 * FieldVector< ctype, m > &y ); 54 * \endcode 55 */ 56 typedef Impl::FieldMatrixHelper< ct > MatrixHelper; 57 58 /** \brief tolerance to numerical algorithms */ toleranceDune::MultiLinearGeometryTraits59 static ct tolerance () { return ct( 16 ) * std::numeric_limits< ct >::epsilon(); } 60 61 /** \brief template specifying the storage for the corners 62 * 63 * Internally, the MultiLinearGeometry needs to store the corners of the 64 * geometry. 65 * 66 * The corner storage may be chosen depending on geometry dimension and 67 * coordinate dimension. It is required to contain a type named Type, e.g., 68 * \code 69 * template< int mydim, int cdim > 70 * struct CornerStorage 71 * { 72 * typedef std::vector< FieldVector< ctype, cdim > > Type; 73 * }; 74 * \endcode 75 * By default, a std::vector of FieldVector is used. 76 * 77 * Apart from being copy constructable and assignable, an \c const corner 78 * storage object \c corners must support the expressions \c 79 * begin(corners), \c end(corners), and subscription \c corners[i]. \c 80 * begin() and \c end() are looked up via ADL and in namespace \c std: 81 * \code 82 * using std::begin; 83 * using std::end; 84 * // it is a const_iterator over the corners in Dune-ordering 85 * auto it = begin(corners); 86 * FieldVector<ctype, cdim> c0 = *it; 87 * auto itend = end(corners); 88 * while(it != itend) { 89 * //... 90 * } 91 * 92 * // elements must be accessible by subscription, indexed in 93 * // Dune-ordering 94 * FieldVector<ctype, cdim> c1 = corners[1]; 95 * \endcode 96 * This means that all of the following qualify: \c 97 * FieldVector<ctype,cdim>[1<<mydim], \c 98 * std::array<FieldVector<ctype,cdim>,(1<<mydim)>, \c 99 * std::vector<FieldVector<ctype,cdim>>. 100 * 101 * \note The expression \c end(corners) isn't actually used by the 102 * implementation currently, but we require it anyway so we can add 103 * runtime checks for the container size when we feel like it. 104 * 105 * It is also possible to use a \c std::reference_wrapper of a suitable 106 * container as the type for the corner storage. The implementation 107 * automatically calls \c corners.get() on internally stored \c 108 * std::reference_wrapper objects before applying \c begin(), \c end(), 109 * or subscription in that case. 110 * 111 * \note Using \c std::reference_wrapper of some container as the corner 112 * storage means that the geometry has no control over the lifetime 113 * of or the access to that container. When the lifetime of the 114 * container ends, or the container itself or its elements are 115 * modified, any geometry object that still references that 116 * container becomes invalid. The only valid operation on invalid 117 * geometry objects are destruction and assignment from another 118 * geometry. If invalidation happens concurrently with some 119 * operation (other than destruction or assignment) on the 120 * geometry, that is a race. 121 * 122 * \tparam mydim geometry dimension 123 * \tparam cdim coordinate dimension 124 */ 125 template< int mydim, int cdim > 126 struct CornerStorage 127 { 128 typedef std::vector< FieldVector< ct, cdim > > Type; 129 }; 130 131 /** \brief will there be only one geometry type for a dimension? 132 * 133 * If there is only a single geometry type for a certain dimension, 134 * <em>hasSingleGeometryType::v</em> can be set to true. 135 * Supporting only one geometry type might yield a gain in performance. 136 * 137 * If <em>hasSingleGeometryType::v</em> is set to true, an additional 138 * parameter <em>topologyId</em> is required. 139 * Here's an example: 140 * \code 141 * static const unsigned int topologyId = GeometryTypes::simplex(dim).id(); 142 * \endcode 143 */ 144 template< int dim > 145 struct hasSingleGeometryType 146 { 147 static const bool v = false; 148 static const unsigned int topologyId = ~0u; 149 }; 150 }; 151 152 153 154 // MultiLinearGeometry 155 // ------------------- 156 157 /** \brief generic geometry implementation based on corner coordinates 158 * 159 * Based on the recursive definition of the reference elements, the 160 * MultiLinearGeometry provides a generic implementation of a geometry given 161 * the corner coordinates. 162 * 163 * The geometric mapping is multilinear in the classical sense only in the 164 * case of cubes; for simplices it is linear. 165 * The name is still justified, because the mapping satisfies the important 166 * property of begin linear along edges. 167 * 168 * \tparam ct coordinate type 169 * \tparam mydim geometry dimension 170 * \tparam cdim coordinate dimension 171 * \tparam Traits traits allowing to tweak some implementation details 172 * (optional) 173 * 174 * The requirements on the traits are documented along with their default, 175 * MultiLinearGeometryTraits. 176 */ 177 template< class ct, int mydim, int cdim, class Traits = MultiLinearGeometryTraits< ct > > 178 class MultiLinearGeometry 179 { 180 typedef MultiLinearGeometry< ct, mydim, cdim, Traits > This; 181 182 public: 183 //! coordinate type 184 typedef ct ctype; 185 186 //! geometry dimension 187 static const int mydimension= mydim; 188 //! coordinate dimension 189 static const int coorddimension = cdim; 190 191 //! type of local coordinates 192 typedef FieldVector< ctype, mydimension > LocalCoordinate; 193 //! type of global coordinates 194 typedef FieldVector< ctype, coorddimension > GlobalCoordinate; 195 //! type of volume 196 typedef ctype Volume; 197 198 //! type of jacobian transposed 199 typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed; 200 201 //! type of jacobian inverse transposed 202 class JacobianInverseTransposed; 203 204 protected: 205 206 typedef Dune::ReferenceElements< ctype, mydimension > ReferenceElements; 207 208 public: 209 210 //! type of reference element 211 typedef typename ReferenceElements::ReferenceElement ReferenceElement; 212 213 private: 214 static const bool hasSingleGeometryType = Traits::template hasSingleGeometryType< mydimension >::v; 215 216 protected: 217 typedef typename Traits::MatrixHelper MatrixHelper; 218 typedef typename std::conditional< hasSingleGeometryType, std::integral_constant< unsigned int, Traits::template hasSingleGeometryType< mydimension >::topologyId >, unsigned int >::type TopologyId; 219 220 public: 221 /** \brief constructor 222 * 223 * \param[in] refElement reference element for the geometry 224 * \param[in] corners corners to store internally 225 * 226 * \note The type of corners is actually a template argument. 227 * It is only required that the internal corner storage can be 228 * constructed from this object. 229 */ 230 template< class Corners > MultiLinearGeometry(const ReferenceElement & refElement,const Corners & corners)231 MultiLinearGeometry ( const ReferenceElement &refElement, 232 const Corners &corners ) 233 : refElement_( refElement ), 234 corners_( corners ) 235 {} 236 237 /** \brief constructor 238 * 239 * \param[in] gt geometry type 240 * \param[in] corners corners to store internally 241 * 242 * \note The type of corners is actually a template argument. 243 * It is only required that the internal corner storage can be 244 * constructed from this object. 245 */ 246 template< class Corners > MultiLinearGeometry(Dune::GeometryType gt,const Corners & corners)247 MultiLinearGeometry ( Dune::GeometryType gt, 248 const Corners &corners ) 249 : refElement_( ReferenceElements::general( gt ) ), 250 corners_( corners ) 251 {} 252 253 /** \brief is this mapping affine? */ affine() const254 bool affine () const 255 { 256 JacobianTransposed jt; 257 return affine( jt ); 258 } 259 260 /** \brief obtain the name of the reference element */ type() const261 Dune::GeometryType type () const { return GeometryType( toUnsignedInt(topologyId()), mydimension ); } 262 263 /** \brief obtain number of corners of the corresponding reference element */ corners() const264 int corners () const { return refElement().size( mydimension ); } 265 266 /** \brief obtain coordinates of the i-th corner */ corner(int i) const267 GlobalCoordinate corner ( int i ) const 268 { 269 assert( (i >= 0) && (i < corners()) ); 270 return std::cref(corners_).get()[ i ]; 271 } 272 273 /** \brief obtain the centroid of the mapping's image */ center() const274 GlobalCoordinate center () const { return global( refElement().position( 0, 0 ) ); } 275 276 /** \brief evaluate the mapping 277 * 278 * \param[in] local local coordinate to map 279 * 280 * \returns corresponding global coordinate 281 */ global(const LocalCoordinate & local) const282 GlobalCoordinate global ( const LocalCoordinate &local ) const 283 { 284 using std::begin; 285 286 auto cit = begin(std::cref(corners_).get()); 287 GlobalCoordinate y; 288 global< false >( topologyId(), std::integral_constant< int, mydimension >(), cit, ctype( 1 ), local, ctype( 1 ), y ); 289 return y; 290 } 291 292 /** \brief evaluate the inverse mapping 293 * 294 * \param[in] globalCoord global coordinate to map 295 * 296 * \return corresponding local coordinate 297 * 298 * \note For given global coordinate y the returned local coordinate x that minimizes 299 * the following function over the local coordinate space spanned by the reference element. 300 * \code 301 * (global( x ) - y).two_norm() 302 * \endcode 303 */ local(const GlobalCoordinate & globalCoord) const304 LocalCoordinate local ( const GlobalCoordinate &globalCoord ) const 305 { 306 const ctype tolerance = Traits::tolerance(); 307 LocalCoordinate x = refElement().position( 0, 0 ); 308 LocalCoordinate dx; 309 const bool affineMapping = this->affine(); 310 do 311 { 312 // Newton's method: DF^n dx^n = F^n, x^{n+1} -= dx^n 313 const GlobalCoordinate dglobal = (*this).global( x ) - globalCoord; 314 const bool invertible = 315 MatrixHelper::template xTRightInvA< mydimension, coorddimension >( jacobianTransposed( x ), dglobal, dx ); 316 if( ! invertible ) 317 return LocalCoordinate( std::numeric_limits< ctype > :: max() ); 318 319 // update x with correction 320 x -= dx; 321 322 // for affine mappings only one iteration is needed 323 if ( affineMapping ) break; 324 } while( dx.two_norm2() > tolerance ); 325 return x; 326 } 327 328 /** \brief obtain the integration element 329 * 330 * If the Jacobian of the mapping is denoted by $J(x)$, the integration 331 * integration element \f$\mu(x)\f$ is given by 332 * \f[ \mu(x) = \sqrt{|\det (J^T(x) J(x))|}.\f] 333 * 334 * \param[in] local local coordinate to evaluate the integration element in 335 * 336 * \returns the integration element \f$\mu(x)\f$. 337 * 338 * \note For affine mappings, it is more efficient to call 339 * jacobianInverseTransposed before integrationElement, if both 340 * are required. 341 */ integrationElement(const LocalCoordinate & local) const342 ctype integrationElement ( const LocalCoordinate &local ) const 343 { 344 return MatrixHelper::template sqrtDetAAT< mydimension, coorddimension >( jacobianTransposed( local ) ); 345 } 346 347 /** \brief obtain the volume of the mapping's image 348 * 349 * \note The current implementation just returns 350 * \code 351 * integrationElement( refElement().position( 0, 0 ) ) * refElement().volume() 352 * \endcode 353 * which is wrong for n-linear surface maps and other nonlinear maps. 354 */ volume() const355 Volume volume () const 356 { 357 return integrationElement( refElement().position( 0, 0 ) ) * refElement().volume(); 358 } 359 360 /** \brief obtain the transposed of the Jacobian 361 * 362 * \param[in] local local coordinate to evaluate Jacobian in 363 * 364 * \returns a reference to the transposed of the Jacobian 365 * 366 * \note The returned reference is reused on the next call to 367 * JacobianTransposed, destroying the previous value. 368 */ jacobianTransposed(const LocalCoordinate & local) const369 JacobianTransposed jacobianTransposed ( const LocalCoordinate &local ) const 370 { 371 using std::begin; 372 373 JacobianTransposed jt; 374 auto cit = begin(std::cref(corners_).get()); 375 jacobianTransposed< false >( topologyId(), std::integral_constant< int, mydimension >(), cit, ctype( 1 ), local, ctype( 1 ), jt ); 376 return jt; 377 } 378 379 /** \brief obtain the transposed of the Jacobian's inverse 380 * 381 * The Jacobian's inverse is defined as a pseudo-inverse. If we denote 382 * the Jacobian by \f$J(x)\f$, the following condition holds: 383 * \f[J^{-1}(x) J(x) = I.\f] 384 */ 385 JacobianInverseTransposed jacobianInverseTransposed ( const LocalCoordinate &local ) const; 386 referenceElement(const MultiLinearGeometry & geometry)387 friend ReferenceElement referenceElement ( const MultiLinearGeometry &geometry ) 388 { 389 return geometry.refElement(); 390 } 391 392 protected: 393 refElement() const394 ReferenceElement refElement () const 395 { 396 return refElement_; 397 } 398 topologyId() const399 TopologyId topologyId () const 400 { 401 return topologyId( std::integral_constant< bool, hasSingleGeometryType >() ); 402 } 403 404 template< bool add, int dim, class CornerIterator > 405 static void global ( TopologyId topologyId, std::integral_constant< int, dim >, 406 CornerIterator &cit, const ctype &df, const LocalCoordinate &x, 407 const ctype &rf, GlobalCoordinate &y ); 408 template< bool add, class CornerIterator > 409 static void global ( TopologyId topologyId, std::integral_constant< int, 0 >, 410 CornerIterator &cit, const ctype &df, const LocalCoordinate &x, 411 const ctype &rf, GlobalCoordinate &y ); 412 413 template< bool add, int rows, int dim, class CornerIterator > 414 static void jacobianTransposed ( TopologyId topologyId, std::integral_constant< int, dim >, 415 CornerIterator &cit, const ctype &df, const LocalCoordinate &x, 416 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt ); 417 template< bool add, int rows, class CornerIterator > 418 static void jacobianTransposed ( TopologyId topologyId, std::integral_constant< int, 0 >, 419 CornerIterator &cit, const ctype &df, const LocalCoordinate &x, 420 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt ); 421 422 template< int dim, class CornerIterator > 423 static bool affine ( TopologyId topologyId, std::integral_constant< int, dim >, CornerIterator &cit, JacobianTransposed &jt ); 424 template< class CornerIterator > 425 static bool affine ( TopologyId topologyId, std::integral_constant< int, 0 >, CornerIterator &cit, JacobianTransposed &jt ); 426 affine(JacobianTransposed & jacobianT) const427 bool affine ( JacobianTransposed &jacobianT ) const 428 { 429 using std::begin; 430 431 auto cit = begin(std::cref(corners_).get()); 432 return affine( topologyId(), std::integral_constant< int, mydimension >(), cit, jacobianT ); 433 } 434 435 private: 436 // The following methods are needed to convert the return type of topologyId to 437 // unsigned int with g++-4.4. It has problems casting integral_constant to the 438 // integral type. toUnsignedInt(unsigned int i)439 static unsigned int toUnsignedInt(unsigned int i) { return i; } 440 template<unsigned int v> toUnsignedInt(std::integral_constant<unsigned int,v>)441 static unsigned int toUnsignedInt(std::integral_constant<unsigned int,v> ) { return v; } topologyId(std::integral_constant<bool,true>) const442 TopologyId topologyId ( std::integral_constant< bool, true > ) const { return TopologyId(); } topologyId(std::integral_constant<bool,false>) const443 unsigned int topologyId ( std::integral_constant< bool, false > ) const { return refElement().type().id(); } 444 445 ReferenceElement refElement_; 446 typename Traits::template CornerStorage< mydimension, coorddimension >::Type corners_; 447 }; 448 449 450 451 // MultiLinearGeometry::JacobianInverseTransposed 452 // ---------------------------------------------- 453 454 template< class ct, int mydim, int cdim, class Traits > 455 class MultiLinearGeometry< ct, mydim, cdim, Traits >::JacobianInverseTransposed 456 : public FieldMatrix< ctype, coorddimension, mydimension > 457 { 458 typedef FieldMatrix< ctype, coorddimension, mydimension > Base; 459 460 public: setup(const JacobianTransposed & jt)461 void setup ( const JacobianTransposed &jt ) 462 { 463 detInv_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jt, static_cast< Base & >( *this ) ); 464 } 465 setupDeterminant(const JacobianTransposed & jt)466 void setupDeterminant ( const JacobianTransposed &jt ) 467 { 468 detInv_ = MatrixHelper::template sqrtDetAAT< mydimension, coorddimension >( jt ); 469 } 470 det() const471 ctype det () const { return ctype( 1 ) / detInv_; } detInv() const472 ctype detInv () const { return detInv_; } 473 474 private: 475 ctype detInv_; 476 }; 477 478 479 480 /** \brief Implement a MultiLinearGeometry with additional caching 481 * 482 * This class implements the same interface and functionality as MultiLinearGeometry. 483 * However, it additionally implements caching for various results. 484 * 485 * \tparam ct coordinate type 486 * \tparam mydim geometry dimension 487 * \tparam cdim coordinate dimension 488 * \tparam Traits traits allowing to tweak some implementation details 489 * (optional) 490 * 491 */ 492 template< class ct, int mydim, int cdim, class Traits = MultiLinearGeometryTraits< ct > > 493 class CachedMultiLinearGeometry 494 : public MultiLinearGeometry< ct, mydim, cdim, Traits > 495 { 496 typedef CachedMultiLinearGeometry< ct, mydim, cdim, Traits > This; 497 typedef MultiLinearGeometry< ct, mydim, cdim, Traits > Base; 498 499 protected: 500 typedef typename Base::MatrixHelper MatrixHelper; 501 502 public: 503 typedef typename Base::ReferenceElement ReferenceElement; 504 505 typedef typename Base::ctype ctype; 506 507 using Base::mydimension; 508 using Base::coorddimension; 509 510 typedef typename Base::LocalCoordinate LocalCoordinate; 511 typedef typename Base::GlobalCoordinate GlobalCoordinate; 512 typedef typename Base::Volume Volume; 513 514 typedef typename Base::JacobianTransposed JacobianTransposed; 515 typedef typename Base::JacobianInverseTransposed JacobianInverseTransposed; 516 517 template< class CornerStorage > CachedMultiLinearGeometry(const ReferenceElement & referenceElement,const CornerStorage & cornerStorage)518 CachedMultiLinearGeometry ( const ReferenceElement &referenceElement, const CornerStorage &cornerStorage ) 519 : Base( referenceElement, cornerStorage ), 520 affine_( Base::affine( jacobianTransposed_ ) ), 521 jacobianInverseTransposedComputed_( false ), 522 integrationElementComputed_( false ) 523 {} 524 525 template< class CornerStorage > CachedMultiLinearGeometry(Dune::GeometryType gt,const CornerStorage & cornerStorage)526 CachedMultiLinearGeometry ( Dune::GeometryType gt, const CornerStorage &cornerStorage ) 527 : Base( gt, cornerStorage ), 528 affine_( Base::affine( jacobianTransposed_ ) ), 529 jacobianInverseTransposedComputed_( false ), 530 integrationElementComputed_( false ) 531 {} 532 533 /** \brief is this mapping affine? */ affine() const534 bool affine () const { return affine_; } 535 536 using Base::corner; 537 538 /** \brief obtain the centroid of the mapping's image */ center() const539 GlobalCoordinate center () const { return global( refElement().position( 0, 0 ) ); } 540 541 /** \brief evaluate the mapping 542 * 543 * \param[in] local local coordinate to map 544 * 545 * \returns corresponding global coordinate 546 */ global(const LocalCoordinate & local) const547 GlobalCoordinate global ( const LocalCoordinate &local ) const 548 { 549 if( affine() ) 550 { 551 GlobalCoordinate global( corner( 0 ) ); 552 jacobianTransposed_.umtv( local, global ); 553 return global; 554 } 555 else 556 return Base::global( local ); 557 } 558 559 /** \brief evaluate the inverse mapping 560 * 561 * \param[in] global global coordinate to map 562 * 563 * \return corresponding local coordinate 564 * 565 * \note For given global coordinate y the returned local coordinate x that minimizes 566 * the following function over the local coordinate space spanned by the reference element. 567 * \code 568 * (global( x ) - y).two_norm() 569 * \endcode 570 */ local(const GlobalCoordinate & global) const571 LocalCoordinate local ( const GlobalCoordinate &global ) const 572 { 573 if( affine() ) 574 { 575 LocalCoordinate local; 576 if( jacobianInverseTransposedComputed_ ) 577 jacobianInverseTransposed_.mtv( global - corner( 0 ), local ); 578 else 579 MatrixHelper::template xTRightInvA< mydimension, coorddimension >( jacobianTransposed_, global - corner( 0 ), local ); 580 return local; 581 } 582 else 583 return Base::local( global ); 584 } 585 586 /** \brief obtain the integration element 587 * 588 * If the Jacobian of the mapping is denoted by $J(x)$, the integration 589 * integration element \f$\mu(x)\f$ is given by 590 * \f[ \mu(x) = \sqrt{|\det (J^T(x) J(x))|}.\f] 591 * 592 * \param[in] local local coordinate to evaluate the integration element in 593 * 594 * \returns the integration element \f$\mu(x)\f$. 595 * 596 * \note For affine mappings, it is more efficient to call 597 * jacobianInverseTransposed before integrationElement, if both 598 * are required. 599 */ integrationElement(const LocalCoordinate & local) const600 ctype integrationElement ( const LocalCoordinate &local ) const 601 { 602 if( affine() ) 603 { 604 if( !integrationElementComputed_ ) 605 { 606 jacobianInverseTransposed_.setupDeterminant( jacobianTransposed_ ); 607 integrationElementComputed_ = true; 608 } 609 return jacobianInverseTransposed_.detInv(); 610 } 611 else 612 return Base::integrationElement( local ); 613 } 614 615 /** \brief obtain the volume of the mapping's image */ volume() const616 Volume volume () const 617 { 618 if( affine() ) 619 return integrationElement( refElement().position( 0, 0 ) ) * refElement().volume(); 620 else 621 return Base::volume(); 622 } 623 624 /** \brief obtain the transposed of the Jacobian 625 * 626 * \param[in] local local coordinate to evaluate Jacobian in 627 * 628 * \returns a reference to the transposed of the Jacobian 629 * 630 * \note The returned reference is reused on the next call to 631 * JacobianTransposed, destroying the previous value. 632 */ jacobianTransposed(const LocalCoordinate & local) const633 JacobianTransposed jacobianTransposed ( const LocalCoordinate &local ) const 634 { 635 if( affine() ) 636 return jacobianTransposed_; 637 else 638 return Base::jacobianTransposed( local ); 639 } 640 641 /** \brief obtain the transposed of the Jacobian's inverse 642 * 643 * The Jacobian's inverse is defined as a pseudo-inverse. If we denote 644 * the Jacobian by \f$J(x)\f$, the following condition holds: 645 * \f[J^{-1}(x) J(x) = I.\f] 646 */ jacobianInverseTransposed(const LocalCoordinate & local) const647 JacobianInverseTransposed jacobianInverseTransposed ( const LocalCoordinate &local ) const 648 { 649 if( affine() ) 650 { 651 if( !jacobianInverseTransposedComputed_ ) 652 { 653 jacobianInverseTransposed_.setup( jacobianTransposed_ ); 654 jacobianInverseTransposedComputed_ = true; 655 integrationElementComputed_ = true; 656 } 657 return jacobianInverseTransposed_; 658 } 659 else 660 return Base::jacobianInverseTransposed( local ); 661 } 662 663 protected: 664 using Base::refElement; 665 666 private: 667 mutable JacobianTransposed jacobianTransposed_; 668 mutable JacobianInverseTransposed jacobianInverseTransposed_; 669 670 mutable bool affine_ : 1; 671 672 mutable bool jacobianInverseTransposedComputed_ : 1; 673 mutable bool integrationElementComputed_ : 1; 674 }; 675 676 677 678 // Implementation of MultiLinearGeometry 679 // ------------------------------------- 680 681 template< class ct, int mydim, int cdim, class Traits > 682 inline typename MultiLinearGeometry< ct, mydim, cdim, Traits >::JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate & local) const683 MultiLinearGeometry< ct, mydim, cdim, Traits >::jacobianInverseTransposed ( const LocalCoordinate &local ) const 684 { 685 JacobianInverseTransposed jit; 686 jit.setup( jacobianTransposed( local ) ); 687 return jit; 688 } 689 690 691 template< class ct, int mydim, int cdim, class Traits > 692 template< bool add, int dim, class CornerIterator > 693 inline void MultiLinearGeometry< ct, mydim, cdim, Traits > global(TopologyId topologyId,std::integral_constant<int,dim>,CornerIterator & cit,const ctype & df,const LocalCoordinate & x,const ctype & rf,GlobalCoordinate & y)694 ::global ( TopologyId topologyId, std::integral_constant< int, dim >, 695 CornerIterator &cit, const ctype &df, const LocalCoordinate &x, 696 const ctype &rf, GlobalCoordinate &y ) 697 { 698 const ctype xn = df*x[ dim-1 ]; 699 const ctype cxn = ctype( 1 ) - xn; 700 701 if( Impl::isPrism( toUnsignedInt(topologyId), mydimension, mydimension-dim ) ) 702 { 703 // apply (1-xn) times mapping for bottom 704 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, rf*cxn, y ); 705 // apply xn times mapping for top 706 global< true >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, rf*xn, y ); 707 } 708 else 709 { 710 assert( Impl::isPyramid( toUnsignedInt(topologyId), mydimension, mydimension-dim ) ); 711 // apply (1-xn) times mapping for bottom (with argument x/(1-xn)) 712 if( cxn > Traits::tolerance() || cxn < -Traits::tolerance() ) 713 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df/cxn, x, rf*cxn, y ); 714 else 715 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, ctype( 0 ), y ); 716 // apply xn times the tip 717 y.axpy( rf*xn, *cit ); 718 ++cit; 719 } 720 } 721 722 template< class ct, int mydim, int cdim, class Traits > 723 template< bool add, class CornerIterator > 724 inline void MultiLinearGeometry< ct, mydim, cdim, Traits > global(TopologyId,std::integral_constant<int,0>,CornerIterator & cit,const ctype &,const LocalCoordinate &,const ctype & rf,GlobalCoordinate & y)725 ::global ( TopologyId, std::integral_constant< int, 0 >, 726 CornerIterator &cit, const ctype &, const LocalCoordinate &, 727 const ctype &rf, GlobalCoordinate &y ) 728 { 729 const GlobalCoordinate &origin = *cit; 730 ++cit; 731 for( int i = 0; i < coorddimension; ++i ) 732 y[ i ] = (add ? y[ i ] + rf*origin[ i ] : rf*origin[ i ]); 733 } 734 735 736 template< class ct, int mydim, int cdim, class Traits > 737 template< bool add, int rows, int dim, class CornerIterator > 738 inline void MultiLinearGeometry< ct, mydim, cdim, Traits > jacobianTransposed(TopologyId topologyId,std::integral_constant<int,dim>,CornerIterator & cit,const ctype & df,const LocalCoordinate & x,const ctype & rf,FieldMatrix<ctype,rows,cdim> & jt)739 ::jacobianTransposed ( TopologyId topologyId, std::integral_constant< int, dim >, 740 CornerIterator &cit, const ctype &df, const LocalCoordinate &x, 741 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt ) 742 { 743 assert( rows >= dim ); 744 745 const ctype xn = df*x[ dim-1 ]; 746 const ctype cxn = ctype( 1 ) - xn; 747 748 auto cit2( cit ); 749 if( Impl::isPrism( toUnsignedInt(topologyId), mydimension, mydimension-dim ) ) 750 { 751 // apply (1-xn) times Jacobian for bottom 752 jacobianTransposed< add >( topologyId, std::integral_constant< int, dim-1 >(), cit2, df, x, rf*cxn, jt ); 753 // apply xn times Jacobian for top 754 jacobianTransposed< true >( topologyId, std::integral_constant< int, dim-1 >(), cit2, df, x, rf*xn, jt ); 755 // compute last row as difference between top value and bottom value 756 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, -rf, jt[ dim-1 ] ); 757 global< true >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, rf, jt[ dim-1 ] ); 758 } 759 else 760 { 761 assert( Impl::isPyramid( toUnsignedInt(topologyId), mydimension, mydimension-dim ) ); 762 /* 763 * In the pyramid case, we need a transformation Tb: B -> R^n for the 764 * base B \subset R^{n-1}. The pyramid transformation is then defined as 765 * T: P \subset R^n -> R^n 766 * (x, xn) |-> (1-xn) Tb(x*) + xn t (x \in R^{n-1}, xn \in R) 767 * with the tip of the pyramid mapped to t and x* = x/(1-xn) 768 * the projection of (x,xn) onto the base. 769 * 770 * For the Jacobi matrix DT we get 771 * DT = ( A | b ) 772 * with A = DTb(x*) (n x n-1 matrix) 773 * and b = dT/dxn (n-dim column vector). 774 * Furthermore 775 * b = -Tb(x*) + t + \sum_i dTb/dx_i(x^*) x_i/(1-xn) 776 * 777 * Note that both A and b are not defined in the pyramid tip (x=0, xn=1)! 778 * Indeed for B the unit square, Tb mapping B to the quadrilateral given 779 * by the vertices (0,0,0), (2,0,0), (0,1,0), (1,1,0) and t=(0,0,1), we get 780 * 781 * T(x,y,xn) = ( x(2-y/(1-xn)), y, xn ) 782 * / 2-y/(1-xn) -x 0 \ 783 * DT(x,y,xn) = | 0 1 0 | 784 * \ 0 0 1 / 785 * which is not continuous for xn -> 1, choose for example 786 * x=0, y=1-xn, xn -> 1 --> DT -> diag(1,1,1) 787 * x=1-xn, y=0, xn -> 1 --> DT -> diag(2,1,1) 788 * 789 * However, for Tb affine-linear, Tb(y) = My + y0, DTb = M: 790 * A = M 791 * b = -M x* - y0 + t + \sum_i M_i x_i/(1-xn) 792 * = -M x* - y0 + t + M x* 793 * = -y0 + t 794 * which is continuous for xn -> 1. Note that this b is also given by 795 * b = -Tb(0) + t + \sum_i dTb/dx_i(0) x_i/1 796 * that is replacing x* by 1 and 1-xn by 1 in the formular above. 797 * 798 * For xn -> 1, we can thus set x*=0, "1-xn"=1 (or anything != 0) and get 799 * the right result in case Tb is affine-linear. 800 */ 801 802 /* The second case effectively results in x* = 0 */ 803 ctype dfcxn = (cxn > Traits::tolerance() || cxn < -Traits::tolerance()) ? ctype(df / cxn) : ctype(0); 804 805 // initialize last row 806 // b = -Tb(x*) 807 // (b = -Tb(0) = -y0 in case xn -> 1 and Tb affine-linear) 808 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, dfcxn, x, -rf, jt[ dim-1 ] ); 809 // b += t 810 jt[ dim-1 ].axpy( rf, *cit ); 811 ++cit; 812 // apply Jacobian for bottom (with argument x/(1-xn)) and correct last row 813 if( add ) 814 { 815 FieldMatrix< ctype, dim-1, coorddimension > jt2; 816 // jt2 = dTb/dx_i(x*) 817 jacobianTransposed< false >( topologyId, std::integral_constant< int, dim-1 >(), cit2, dfcxn, x, rf, jt2 ); 818 // A = dTb/dx_i(x*) (jt[j], j=0..dim-1) 819 // b += \sum_i dTb/dx_i(x*) x_i/(1-xn) (jt[dim-1]) 820 // (b += 0 in case xn -> 1) 821 for( int j = 0; j < dim-1; ++j ) 822 { 823 jt[ j ] += jt2[ j ]; 824 jt[ dim-1 ].axpy( dfcxn*x[ j ], jt2[ j ] ); 825 } 826 } 827 else 828 { 829 // jt = dTb/dx_i(x*) 830 jacobianTransposed< false >( topologyId, std::integral_constant< int, dim-1 >(), cit2, dfcxn, x, rf, jt ); 831 // b += \sum_i dTb/dx_i(x*) x_i/(1-xn) 832 for( int j = 0; j < dim-1; ++j ) 833 jt[ dim-1 ].axpy( dfcxn*x[ j ], jt[ j ] ); 834 } 835 } 836 } 837 838 template< class ct, int mydim, int cdim, class Traits > 839 template< bool add, int rows, class CornerIterator > 840 inline void MultiLinearGeometry< ct, mydim, cdim, Traits > jacobianTransposed(TopologyId,std::integral_constant<int,0>,CornerIterator & cit,const ctype &,const LocalCoordinate &,const ctype &,FieldMatrix<ctype,rows,cdim> &)841 ::jacobianTransposed ( TopologyId, std::integral_constant< int, 0 >, 842 CornerIterator &cit, const ctype &, const LocalCoordinate &, 843 const ctype &, FieldMatrix< ctype, rows, cdim > & ) 844 { 845 ++cit; 846 } 847 848 849 850 template< class ct, int mydim, int cdim, class Traits > 851 template< int dim, class CornerIterator > 852 inline bool MultiLinearGeometry< ct, mydim, cdim, Traits > affine(TopologyId topologyId,std::integral_constant<int,dim>,CornerIterator & cit,JacobianTransposed & jt)853 ::affine ( TopologyId topologyId, std::integral_constant< int, dim >, CornerIterator &cit, JacobianTransposed &jt ) 854 { 855 const GlobalCoordinate &orgBottom = *cit; 856 if( !affine( topologyId, std::integral_constant< int, dim-1 >(), cit, jt ) ) 857 return false; 858 const GlobalCoordinate &orgTop = *cit; 859 860 if( Impl::isPrism( toUnsignedInt(topologyId), mydimension, mydimension-dim ) ) 861 { 862 JacobianTransposed jtTop; 863 if( !affine( topologyId, std::integral_constant< int, dim-1 >(), cit, jtTop ) ) 864 return false; 865 866 // check whether both jacobians are identical 867 ctype norm( 0 ); 868 for( int i = 0; i < dim-1; ++i ) 869 norm += (jtTop[ i ] - jt[ i ]).two_norm2(); 870 if( norm >= Traits::tolerance() ) 871 return false; 872 } 873 else 874 ++cit; 875 jt[ dim-1 ] = orgTop - orgBottom; 876 return true; 877 } 878 879 template< class ct, int mydim, int cdim, class Traits > 880 template< class CornerIterator > 881 inline bool MultiLinearGeometry< ct, mydim, cdim, Traits > affine(TopologyId,std::integral_constant<int,0>,CornerIterator & cit,JacobianTransposed &)882 ::affine ( TopologyId, std::integral_constant< int, 0 >, CornerIterator &cit, JacobianTransposed & ) 883 { 884 ++cit; 885 return true; 886 } 887 888 } // namespace Dune 889 890 #endif // #ifndef DUNE_GEOMETRY_MULTILINEARGEOMETRY_HH 891