3 #ifndef DUNE_GEOMETRY_MULTILINEARGEOMETRY_HH
4 #define DUNE_GEOMETRY_MULTILINEARGEOMETRY_HH
10 #include <dune/common/fmatrix.hh>
11 #include <dune/common/fvector.hh>
12 #include <dune/common/typetraits.hh>
25 template<
class ctype,
int dim >
26 class ReferenceElement;
28 template<
class ctype,
int dim >
29 struct ReferenceElements;
69 static ct
tolerance () {
return ct( 16 ) * std::numeric_limits< ct >::epsilon(); }
94 template<
int mydim,
int cdim >
97 typedef std::vector< FieldVector< ct, cdim > >
Type;
116 static const bool v =
false;
146 template<
class ct,
int mydim,
int cdim,
class Traits = MultiLinearGeometryTraits< ct > >
175 static const bool hasSingleGeometryType = Traits::template hasSingleGeometryType< mydimension >::v;
179 typedef typename conditional< hasSingleGeometryType, integral_constant< unsigned int, Traits::template hasSingleGeometryType< mydimension >::topologyId >,
unsigned int >
::type TopologyId;
184 typedef typename Traits::template CornerStorage< mydimension, coorddimension >::Type::const_iterator CornerIterator;
196 template<
class Corners >
199 : refElement_( &refElement ),
212 template<
class Corners >
215 : refElement_( &ReferenceElements::general( gt ) ),
222 JacobianTransposed jt;
235 assert( (i >= 0) && (i <
corners()) );
236 return corners_[ i ];
250 CornerIterator cit = corners_.begin();
270 const ctype tolerance = Traits::tolerance();
276 const GlobalCoordinate dglobal = (*this).global( x ) -
global;
277 MatrixHelper::template xTRightInvA< mydimension, coorddimension >(
jacobianTransposed( x ), dglobal, dx );
279 }
while( dx.two_norm2() > tolerance );
299 return MatrixHelper::template sqrtDetAAT< mydimension, coorddimension >(
jacobianTransposed( local ) );
326 JacobianTransposed jt;
327 CornerIterator cit = corners_.begin();
328 jacobianTransposed< false >(
topologyId(), integral_constant< int, mydimension >(), cit,
ctype( 1 ),
local,
ctype( 1 ), jt );
341 const ReferenceElement &
refElement ()
const {
return *refElement_; }
345 return topologyId( integral_constant< bool, hasSingleGeometryType >() );
348 template<
bool add,
int dim >
350 CornerIterator &cit,
const ctype &df,
const LocalCoordinate &x,
351 const ctype &rf, GlobalCoordinate &y );
354 CornerIterator &cit,
const ctype &df,
const LocalCoordinate &x,
355 const ctype &rf, GlobalCoordinate &y );
357 template<
bool add,
int rows,
int dim >
359 CornerIterator &cit,
const ctype &df,
const LocalCoordinate &x,
360 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt );
361 template<
bool add,
int rows >
363 CornerIterator &cit,
const ctype &df,
const LocalCoordinate &x,
364 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt );
367 static bool affine ( TopologyId
topologyId, integral_constant< int, dim >, CornerIterator &cit, JacobianTransposed &jt );
368 static bool affine ( TopologyId
topologyId, integral_constant< int, 0 >, CornerIterator &cit, JacobianTransposed &jt );
372 CornerIterator cit = corners_.begin();
373 return affine(
topologyId(), integral_constant< int, mydimension >(), cit, jacobianTransposed );
380 static unsigned int toUnsignedInt(
unsigned int i) {
return i; }
381 template<
unsigned int v>
382 static unsigned int toUnsignedInt(std::integral_constant<unsigned int,v> i) {
return v; }
386 const ReferenceElement *refElement_;
387 typename Traits::template CornerStorage< mydimension, coorddimension >::Type corners_;
395 template<
class ct,
int mydim,
int cdim,
class Traits >
397 :
public FieldMatrix< ctype, coorddimension, mydimension >
399 typedef FieldMatrix< ctype, coorddimension, mydimension > Base;
402 void setup (
const JacobianTransposed &jt )
404 detInv_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jt,
static_cast< Base &
>( *this ) );
409 detInv_ = MatrixHelper::template sqrtDetAAT< mydimension, coorddimension >( jt );
412 ctype
det ()
const {
return ctype( 1 ) / detInv_; }
413 ctype
detInv ()
const {
return detInv_; }
433 template<
class ct,
int mydim,
int cdim,
class Traits = MultiLinearGeometryTraits< ct > >
457 template<
class CornerStorage >
459 : Base( refElement, cornerStorage ),
460 affine_( Base::
affine( jacobianTransposed_ ) ),
461 jacobianInverseTransposedComputed_( false ),
462 integrationElementComputed_( false )
465 template<
class CornerStorage >
467 : Base( gt, cornerStorage ),
468 affine_( Base::
affine( jacobianTransposed_ ) ),
469 jacobianInverseTransposedComputed_( false ),
470 integrationElementComputed_( false )
492 jacobianTransposed_.umtv( local, global );
515 LocalCoordinate
local;
516 if( jacobianInverseTransposedComputed_ )
517 jacobianInverseTransposed_.mtv( global -
corner( 0 ), local );
519 MatrixHelper::template xTRightInvA< mydimension, coorddimension >( jacobianTransposed_, global -
corner( 0 ),
local );
544 if( !integrationElementComputed_ )
547 integrationElementComputed_ =
true;
549 return jacobianInverseTransposed_.
detInv();
576 return jacobianTransposed_;
591 if( !jacobianInverseTransposedComputed_ )
593 jacobianInverseTransposed_.
setup( jacobianTransposed_ );
594 jacobianInverseTransposedComputed_ =
true;
595 integrationElementComputed_ =
true;
597 return jacobianInverseTransposed_;
607 mutable JacobianTransposed jacobianTransposed_;
608 mutable JacobianInverseTransposed jacobianInverseTransposed_;
610 mutable bool affine_ : 1;
612 mutable bool jacobianInverseTransposedComputed_ : 1;
613 mutable bool integrationElementComputed_ : 1;
621 template<
class ct,
int mydim,
int cdim,
class Traits >
622 inline typename MultiLinearGeometry< ct, mydim, cdim, Traits >::JacobianInverseTransposed
626 jit.
setup( jacobianTransposed( local ) );
631 template<
class ct,
int mydim,
int cdim,
class Traits >
632 template<
bool add,
int dim >
634 ::global ( TopologyId topologyId, integral_constant< int, dim >,
635 CornerIterator &cit,
const ctype &df,
const LocalCoordinate &x,
636 const ctype &rf, GlobalCoordinate &y )
638 const ctype xn = df*x[ dim-1 ];
639 const ctype cxn = ctype( 1 ) - xn;
644 global< add >( topologyId, integral_constant< int, dim-1 >(), cit, df, x, rf*cxn, y );
646 global< true >( topologyId, integral_constant< int, dim-1 >(), cit, df, x, rf*xn, y );
652 if( cxn > Traits::tolerance() || cxn < -Traits::tolerance() )
653 global< add >( topologyId, integral_constant< int, dim-1 >(), cit, df/cxn, x, rf*cxn, y );
655 global< add >( topologyId, integral_constant< int, dim-1 >(), cit, df, x, ctype( 0 ), y );
657 y.axpy( rf*xn, *cit );
662 template<
class ct,
int mydim,
int cdim,
class Traits >
665 ::global ( TopologyId topologyId, integral_constant< int, 0 >,
666 CornerIterator &cit,
const ctype &df,
const LocalCoordinate &x,
667 const ctype &rf, GlobalCoordinate &y )
669 const GlobalCoordinate &origin = *cit;
671 for(
int i = 0; i < coorddimension; ++i )
672 y[ i ] = (add ? y[ i ] + rf*origin[ i ] : rf*origin[ i ]);
676 template<
class ct,
int mydim,
int cdim,
class Traits >
677 template<
bool add,
int rows,
int dim >
680 CornerIterator &cit,
const ctype &df,
const LocalCoordinate &x,
681 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt )
683 assert( rows >= dim );
685 const ctype xn = df*x[ dim-1 ];
686 const ctype cxn = ctype( 1 ) - xn;
688 CornerIterator cit2( cit );
692 jacobianTransposed< add >( topologyId, integral_constant< int, dim-1 >(), cit2, df, x, rf*cxn, jt );
694 jacobianTransposed< true >( topologyId, integral_constant< int, dim-1 >(), cit2, df, x, rf*xn, jt );
696 global< add >( topologyId, integral_constant< int, dim-1 >(), cit, df, x, -rf, jt[ dim-1 ] );
697 global< true >( topologyId, integral_constant< int, dim-1 >(), cit, df, x, rf, jt[ dim-1 ] );
743 ctype dfcxn = (cxn > Traits::tolerance() || cxn < -Traits::tolerance()) ? ctype(df / cxn) : ctype(0);
748 global< add >( topologyId, integral_constant< int, dim-1 >(), cit, dfcxn, x, -rf, jt[ dim-1 ] );
750 jt[ dim-1 ].axpy( rf, *cit );
755 FieldMatrix< ctype, dim-1, coorddimension > jt2;
757 jacobianTransposed< false >( topologyId, integral_constant< int, dim-1 >(), cit2, dfcxn, x, rf, jt2 );
761 for(
int j = 0; j < dim-1; ++j )
764 jt[ dim-1 ].axpy( dfcxn*x[ j ], jt2[ j ] );
770 jacobianTransposed< false >( topologyId, integral_constant< int, dim-1 >(), cit2, dfcxn, x, rf, jt );
772 for(
int j = 0; j < dim-1; ++j )
773 jt[ dim-1 ].axpy( dfcxn*x[ j ], jt[ j ] );
778 template<
class ct,
int mydim,
int cdim,
class Traits >
779 template<
bool add,
int rows >
782 CornerIterator &cit,
const ctype &df,
const LocalCoordinate &x,
783 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt )
790 template<
class ct,
int mydim,
int cdim,
class Traits >
793 ::affine ( TopologyId topologyId, integral_constant< int, dim >, CornerIterator &cit, JacobianTransposed &jt )
795 const GlobalCoordinate &orgBottom = *cit;
796 if( !affine( topologyId, integral_constant< int, dim-1 >(), cit, jt ) )
798 const GlobalCoordinate &orgTop = *cit;
802 JacobianTransposed jtTop;
803 if( !affine( topologyId, integral_constant< int, dim-1 >(), cit, jtTop ) )
808 for(
int i = 0; i < dim-1; ++i )
809 norm += (jtTop[ i ] - jt[ i ]).two_norm2();
810 if( norm >= Traits::tolerance() )
815 jt[ dim-1 ] = orgTop - orgBottom;
819 template<
class ct,
int mydim,
int cdim,
class Traits >
821 ::affine ( TopologyId topologyId, integral_constant< int, 0 >, CornerIterator &cit, JacobianTransposed &jt )
829 #endif // #ifndef DUNE_GEOMETRY_MULTILINEARGEOMETRY_HH
void setup(const JacobianTransposed &jt)
Definition: multilineargeometry.hh:402
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
type of jacobian transposed
Definition: multilineargeometry.hh:166
bool affine() const
is this mapping affine?
Definition: multilineargeometry.hh:474
unsigned int id() const
Return the topology id the type.
Definition: type.hh:326
ctype integrationElement(const LocalCoordinate &local) const
obtain the integration element
Definition: multilineargeometry.hh:540
int size(int c) const
number of subentities of codimension c
Definition: referenceelements.hh:82
Dune::GeometryType type() const
obtain the name of the reference element
Definition: multilineargeometry.hh:227
const ReferenceElement & refElement() const
Definition: multilineargeometry.hh:341
Definition: affinegeometry.hh:18
int corners() const
obtain number of corners of the corresponding reference element
Definition: multilineargeometry.hh:230
Class providing access to the singletons of the reference elements.
Definition: affinegeometry.hh:28
Base::MatrixHelper MatrixHelper
Definition: multilineargeometry.hh:441
MultiLinearGeometry(const ReferenceElement &refElement, const Corners &corners)
constructor
Definition: multilineargeometry.hh:197
generic geometry implementation based on corner coordinates
Definition: multilineargeometry.hh:147
static const unsigned int topologyId
Definition: multilineargeometry.hh:117
GlobalCoordinate center() const
obtain the centroid of the mapping's image
Definition: multilineargeometry.hh:479
GlobalCoordinate global(const LocalCoordinate &local) const
evaluate the mapping
Definition: multilineargeometry.hh:487
Base::JacobianTransposed JacobianTransposed
Definition: multilineargeometry.hh:454
ct ctype
coordinate type
Definition: multilineargeometry.hh:153
will there be only one geometry type for a dimension?
Definition: multilineargeometry.hh:114
FieldVector< ctype, coorddimension > GlobalCoordinate
type of global coordinates
Definition: multilineargeometry.hh:163
bool isPyramid(unsigned int topologyId, int dim, int codim=0)
check whether a pyramid construction was used to create a given codimension
Definition: topologytypes.hh:150
Definition: multilineargeometry.hh:396
ctype integrationElement(const LocalCoordinate &local) const
obtain the integration element
Definition: multilineargeometry.hh:297
LocalCoordinate local(const GlobalCoordinate &global) const
evaluate the inverse mapping
Definition: multilineargeometry.hh:268
GlobalCoordinate global(const LocalCoordinate &local) const
evaluate the mapping
Definition: multilineargeometry.hh:248
Definition: matrixhelper.hh:32
bool affine(JacobianTransposed &jacobianTransposed) const
Definition: multilineargeometry.hh:370
GlobalCoordinate center() const
obtain the centroid of the mapping's image
Definition: multilineargeometry.hh:240
Dune::ReferenceElement< ctype, mydimension > ReferenceElement
type of reference element
Definition: multilineargeometry.hh:169
This class provides access to geometric and topological properties of a reference element...
Definition: affinegeometry.hh:25
std::vector< FieldVector< ct, cdim > > Type
Definition: multilineargeometry.hh:97
static ct tolerance()
tolerance to numerical algorithms
Definition: multilineargeometry.hh:69
static const bool v
Definition: multilineargeometry.hh:116
TopologyId topologyId() const
Definition: multilineargeometry.hh:343
ctype volume() const
obtain the volume of the mapping's image
Definition: multilineargeometry.hh:556
bool isPrism(unsigned int topologyId, int dim, int codim=0)
check whether a prism construction was used to create a given codimension
Definition: topologytypes.hh:168
Base::ReferenceElement ReferenceElement
Definition: multilineargeometry.hh:444
ctype det() const
Definition: multilineargeometry.hh:412
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian's inverse
Definition: multilineargeometry.hh:587
template specifying the storage for the corners
Definition: multilineargeometry.hh:95
FieldVector< ctype, mydimension > LocalCoordinate
type of local coordinates
Definition: multilineargeometry.hh:161
Implement a MultiLinearGeometry with additional caching.
Definition: multilineargeometry.hh:434
GenericGeometry::MatrixHelper< GenericGeometry::DuneCoordTraits< ct > > MatrixHelper
helper structure containing some matrix routines
Definition: multilineargeometry.hh:66
conditional< hasSingleGeometryType, integral_constant< unsigned int, Traits::template hasSingleGeometryType< mydimension >::topologyId >, unsigned int >::type TopologyId
Definition: multilineargeometry.hh:179
A unique label for each type of element that can occur in a grid.
Base::JacobianInverseTransposed JacobianInverseTransposed
Definition: multilineargeometry.hh:455
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian's inverse
Definition: multilineargeometry.hh:623
GlobalCoordinate corner(int i) const
obtain coordinates of the i-th corner
Definition: multilineargeometry.hh:233
MultiLinearGeometry(Dune::GeometryType gt, const Corners &corners)
constructor
Definition: multilineargeometry.hh:213
LocalCoordinate local(const GlobalCoordinate &global) const
evaluate the inverse mapping
Definition: multilineargeometry.hh:511
Base::LocalCoordinate LocalCoordinate
Definition: multilineargeometry.hh:451
bool affine() const
is this mapping affine?
Definition: multilineargeometry.hh:220
Base::ctype ctype
Definition: multilineargeometry.hh:446
CachedMultiLinearGeometry(Dune::GeometryType gt, const CornerStorage &cornerStorage)
Definition: multilineargeometry.hh:466
const FieldVector< ctype, dim > & position(int i, int c) const
position of the barycenter of entity (i,c)
Definition: referenceelements.hh:150
ctype volume() const
obtain the volume of the mapping's image
Definition: multilineargeometry.hh:310
Traits::MatrixHelper MatrixHelper
Definition: multilineargeometry.hh:178
default traits class for MultiLinearGeometry
Definition: multilineargeometry.hh:46
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian
Definition: multilineargeometry.hh:573
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:24
Base::GlobalCoordinate GlobalCoordinate
Definition: multilineargeometry.hh:452
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian
Definition: multilineargeometry.hh:324
static const int coorddimension
coordinate dimension
Definition: multilineargeometry.hh:158
ctype volume() const
obtain the volume of the reference element
Definition: referenceelements.hh:210
const GeometryType & type(int i, int c) const
obtain the type of subentity (i,c)
Definition: referenceelements.hh:132
CachedMultiLinearGeometry(const ReferenceElement &refElement, const CornerStorage &cornerStorage)
Definition: multilineargeometry.hh:458
ctype detInv() const
Definition: multilineargeometry.hh:413
static const int mydimension
geometry dimension
Definition: multilineargeometry.hh:156
void setupDeterminant(const JacobianTransposed &jt)
Definition: multilineargeometry.hh:407
Dune::ReferenceElements< ctype, mydimension > ReferenceElements
Definition: multilineargeometry.hh:181