3 #ifndef DUNE_GEOMETRY_MULTILINEARGEOMETRY_HH 4 #define DUNE_GEOMETRY_MULTILINEARGEOMETRY_HH 12 #include <dune/common/fmatrix.hh> 13 #include <dune/common/fvector.hh> 14 #include <dune/common/typetraits.hh> 26 template<
class ctype,
int dim >
27 class ReferenceElement;
29 template<
class ctype,
int dim >
30 struct ReferenceElements;
70 static ct
tolerance () {
return ct( 16 ) * std::numeric_limits< ct >::epsilon(); }
136 template<
int mydim,
int cdim >
139 typedef std::vector< FieldVector< ct, cdim > >
Type;
158 static const bool v =
false;
159 static const unsigned int topologyId = ~0u;
188 template<
class ct,
int mydim,
int cdim,
class Traits = MultiLinearGeometryTraits< ct > >
198 static const int mydimension= mydim;
200 static const int coorddimension = cdim;
221 typedef typename std::conditional< hasSingleGeometryType, std::integral_constant< unsigned int, Traits::template hasSingleGeometryType< mydimension >::topologyId >,
unsigned int >::type
TopologyId;
235 template<
class Corners >
237 const Corners &corners )
238 : refElement_( &refElement ),
251 template<
class Corners >
253 const Corners &corners )
254 : refElement_( &ReferenceElements::general( gt ) ),
261 JacobianTransposed jt;
269 int corners ()
const {
return refElement().size( mydimension ); }
274 assert( (i >= 0) && (i < corners()) );
275 return std::cref(corners_).get()[ i ];
279 GlobalCoordinate
center ()
const {
return global( refElement().position( 0, 0 ) ); }
287 GlobalCoordinate
global (
const LocalCoordinate &local )
const 291 auto cit = begin(std::cref(corners_).
get());
293 global< false >( topologyId(), std::integral_constant< int, mydimension >(), cit, ctype( 1 ), local, ctype( 1 ), y );
309 LocalCoordinate
local (
const GlobalCoordinate &globalCoord )
const 311 const ctype
tolerance = Traits::tolerance();
312 LocalCoordinate x = refElement().position( 0, 0 );
317 const GlobalCoordinate dglobal = (*this).global( x ) - globalCoord;
318 MatrixHelper::template xTRightInvA< mydimension, coorddimension >( jacobianTransposed( x ), dglobal, dx );
340 return MatrixHelper::template sqrtDetAAT< mydimension, coorddimension >( jacobianTransposed( local ) );
353 return integrationElement( refElement().position( 0, 0 ) ) * refElement().volume();
369 JacobianTransposed jt;
370 auto cit = begin(std::cref(corners_).
get());
371 jacobianTransposed< false >( topologyId(), std::integral_constant< int, mydimension >(), cit, ctype( 1 ), local, ctype( 1 ), jt );
386 const ReferenceElement &
refElement ()
const {
return *refElement_; }
390 return topologyId( std::integral_constant< bool, hasSingleGeometryType >() );
393 template<
bool add,
int dim,
class CornerIterator >
394 static void global ( TopologyId topologyId, std::integral_constant< int, dim >,
395 CornerIterator &cit,
const ctype &df,
const LocalCoordinate &x,
396 const ctype &rf, GlobalCoordinate &y );
397 template<
bool add,
class CornerIterator >
398 static void global ( TopologyId topologyId, std::integral_constant< int, 0 >,
399 CornerIterator &cit,
const ctype &df,
const LocalCoordinate &x,
400 const ctype &rf, GlobalCoordinate &y );
402 template<
bool add,
int rows,
int dim,
class CornerIterator >
403 static void jacobianTransposed ( TopologyId topologyId, std::integral_constant< int, dim >,
404 CornerIterator &cit,
const ctype &df,
const LocalCoordinate &x,
405 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt );
406 template<
bool add,
int rows,
class CornerIterator >
407 static void jacobianTransposed ( TopologyId topologyId, std::integral_constant< int, 0 >,
408 CornerIterator &cit,
const ctype &df,
const LocalCoordinate &x,
409 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt );
411 template<
int dim,
class CornerIterator >
412 static bool affine ( TopologyId topologyId, std::integral_constant< int, dim >, CornerIterator &cit, JacobianTransposed &jt );
413 template<
class CornerIterator >
414 static bool affine ( TopologyId topologyId, std::integral_constant< int, 0 >, CornerIterator &cit, JacobianTransposed &jt );
416 bool affine ( JacobianTransposed &jacobianT )
const 420 auto cit = begin(std::cref(corners_).
get());
421 return affine( topologyId(), std::integral_constant< int, mydimension >(), cit, jacobianT );
428 static unsigned int toUnsignedInt(
unsigned int i) {
return i; }
429 template<
unsigned int v>
430 static unsigned int toUnsignedInt(std::integral_constant<unsigned int,v> i) {
return v; }
431 TopologyId topologyId ( std::integral_constant< bool, true > )
const {
return TopologyId(); }
432 unsigned int topologyId ( std::integral_constant< bool, false > )
const {
return refElement().type().id(); }
434 const ReferenceElement *refElement_;
443 template<
class ct,
int mydim,
int cdim,
class Traits >
445 :
public FieldMatrix< ctype, coorddimension, mydimension >
447 typedef FieldMatrix< ctype, coorddimension, mydimension > Base;
452 detInv_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jt,
static_cast< Base &
>( *this ) );
457 detInv_ = MatrixHelper::template sqrtDetAAT< mydimension, coorddimension >( jt );
481 template<
class ct,
int mydim,
int cdim,
class Traits = MultiLinearGeometryTraits< ct > >
496 using Base::mydimension;
497 using Base::coorddimension;
505 template<
class CornerStorage >
507 : Base( referenceElement, cornerStorage ),
508 affine_( Base::affine( jacobianTransposed_ ) ),
509 jacobianInverseTransposedComputed_( false ),
510 integrationElementComputed_( false )
513 template<
class CornerStorage >
515 : Base( gt, cornerStorage ),
516 affine_( Base::affine( jacobianTransposed_ ) ),
517 jacobianInverseTransposedComputed_( false ),
518 integrationElementComputed_( false )
527 GlobalCoordinate
center ()
const {
return global( refElement().position( 0, 0 ) ); }
535 GlobalCoordinate
global (
const LocalCoordinate &local )
const 539 GlobalCoordinate global( corner( 0 ) );
540 jacobianTransposed_.umtv( local, global );
544 return Base::global( local );
559 LocalCoordinate
local (
const GlobalCoordinate &global )
const 563 LocalCoordinate local;
564 if( jacobianInverseTransposedComputed_ )
565 jacobianInverseTransposed_.mtv( global - corner( 0 ), local );
567 MatrixHelper::template xTRightInvA< mydimension, coorddimension >( jacobianTransposed_, global - corner( 0 ), local );
571 return Base::local( global );
592 if( !integrationElementComputed_ )
594 jacobianInverseTransposed_.setupDeterminant( jacobianTransposed_ );
595 integrationElementComputed_ =
true;
597 return jacobianInverseTransposed_.detInv();
600 return Base::integrationElement( local );
607 return integrationElement( refElement().position( 0, 0 ) ) * refElement().volume();
609 return Base::volume();
624 return jacobianTransposed_;
626 return Base::jacobianTransposed( local );
639 if( !jacobianInverseTransposedComputed_ )
641 jacobianInverseTransposed_.
setup( jacobianTransposed_ );
642 jacobianInverseTransposedComputed_ =
true;
643 integrationElementComputed_ =
true;
645 return jacobianInverseTransposed_;
648 return Base::jacobianInverseTransposed( local );
652 using Base::refElement;
655 mutable JacobianTransposed jacobianTransposed_;
656 mutable JacobianInverseTransposed jacobianInverseTransposed_;
658 mutable bool affine_ : 1;
660 mutable bool jacobianInverseTransposedComputed_ : 1;
661 mutable bool integrationElementComputed_ : 1;
669 template<
class ct,
int mydim,
int cdim,
class Traits >
674 jit.
setup( jacobianTransposed( local ) );
679 template<
class ct,
int mydim,
int cdim,
class Traits >
680 template<
bool add,
int dim,
class CornerIterator >
686 const ctype xn = df*x[ dim-1 ];
689 if(
Impl::isPrism( toUnsignedInt(topologyId), mydimension, mydimension-dim ) )
692 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, rf*cxn, y );
694 global< true >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, rf*xn, y );
698 assert(
Impl::isPyramid( toUnsignedInt(topologyId), mydimension, mydimension-dim ) );
700 if( cxn > Traits::tolerance() || cxn < -Traits::tolerance() )
701 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df/cxn, x, rf*cxn, y );
703 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x,
ctype( 0 ), y );
705 y.axpy( rf*xn, *cit );
710 template<
class ct,
int mydim,
int cdim,
class Traits >
711 template<
bool add,
class CornerIterator >
719 for(
int i = 0; i < coorddimension; ++i )
720 y[ i ] = (add ? y[ i ] + rf*origin[ i ] : rf*origin[ i ]);
724 template<
class ct,
int mydim,
int cdim,
class Traits >
725 template<
bool add,
int rows,
int dim,
class CornerIterator >
729 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt )
731 assert( rows >= dim );
733 const ctype xn = df*x[ dim-1 ];
737 if(
Impl::isPrism( toUnsignedInt(topologyId), mydimension, mydimension-dim ) )
740 jacobianTransposed< add >( topologyId, std::integral_constant< int, dim-1 >(), cit2, df, x, rf*cxn, jt );
742 jacobianTransposed< true >( topologyId, std::integral_constant< int, dim-1 >(), cit2, df, x, rf*xn, jt );
744 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, -rf, jt[ dim-1 ] );
745 global< true >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, rf, jt[ dim-1 ] );
749 assert(
Impl::isPyramid( toUnsignedInt(topologyId), mydimension, mydimension-dim ) );
791 ctype dfcxn = (cxn > Traits::tolerance() || cxn < -Traits::tolerance()) ?
ctype(df / cxn) :
ctype(0);
796 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, dfcxn, x, -rf, jt[ dim-1 ] );
798 jt[ dim-1 ].axpy( rf, *cit );
803 FieldMatrix<
ctype, dim-1, coorddimension > jt2;
805 jacobianTransposed< false >( topologyId, std::integral_constant< int, dim-1 >(), cit2, dfcxn, x, rf, jt2 );
809 for(
int j = 0; j < dim-1; ++j )
812 jt[ dim-1 ].axpy( dfcxn*x[ j ], jt2[ j ] );
818 jacobianTransposed< false >( topologyId, std::integral_constant< int, dim-1 >(), cit2, dfcxn, x, rf, jt );
820 for(
int j = 0; j < dim-1; ++j )
821 jt[ dim-1 ].axpy( dfcxn*x[ j ], jt[ j ] );
826 template<
class ct,
int mydim,
int cdim,
class Traits >
827 template<
bool add,
int rows,
class CornerIterator >
831 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt )
838 template<
class ct,
int mydim,
int cdim,
class Traits >
839 template<
int dim,
class CornerIterator >
844 if( !affine( topologyId, std::integral_constant< int, dim-1 >(), cit, jt ) )
848 if(
Impl::isPrism( toUnsignedInt(topologyId), mydimension, mydimension-dim ) )
851 if( !affine( topologyId, std::integral_constant< int, dim-1 >(), cit, jtTop ) )
856 for(
int i = 0; i < dim-1; ++i )
857 norm += (jtTop[ i ] - jt[ i ]).two_norm2();
858 if( norm >= Traits::tolerance() )
863 jt[ dim-1 ] = orgTop - orgBottom;
867 template<
class ct,
int mydim,
int cdim,
class Traits >
868 template<
class CornerIterator >
878 #endif // #ifndef DUNE_GEOMETRY_MULTILINEARGEOMETRY_HH JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian's inverse
Definition: multilineargeometry.hh:635
ctype integrationElement(const LocalCoordinate &local) const
obtain the integration element
Definition: multilineargeometry.hh:588
ct ctype
coordinate type
Definition: multilineargeometry.hh:195
Base::ctype ctype
Definition: multilineargeometry.hh:494
Class providing access to the singletons of the reference elements.
Definition: affinegeometry.hh:28
ctype detInv() const
Definition: multilineargeometry.hh:461
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:267
const ReferenceElement & refElement() const
Definition: multilineargeometry.hh:386
FieldVector< ctype, coorddimension > GlobalCoordinate
type of global coordinates
Definition: multilineargeometry.hh:205
Dune::ReferenceElements< ctype, mydimension > ReferenceElements
Definition: multilineargeometry.hh:223
static ct tolerance()
tolerance to numerical algorithms
Definition: multilineargeometry.hh:70
Base::ReferenceElement ReferenceElement
Definition: multilineargeometry.hh:492
GlobalCoordinate corner(int i) const
obtain coordinates of the i-th corner
Definition: multilineargeometry.hh:272
default traits class for MultiLinearGeometry
Definition: multilineargeometry.hh:47
ctype integrationElement(const LocalCoordinate &local) const
obtain the integration element
Definition: multilineargeometry.hh:338
bool affine() const
is this mapping affine?
Definition: multilineargeometry.hh:259
This class provides access to geometric and topological properties of a reference element...
Definition: affinegeometry.hh:25
ctype volume() const
obtain the volume of the mapping's image
Definition: multilineargeometry.hh:604
std::vector< FieldVector< ct, cdim > > Type
Definition: multilineargeometry.hh:139
Implement a MultiLinearGeometry with additional caching.
Definition: multilineargeometry.hh:482
Definition: affinegeometry.hh:18
will there be only one geometry type for a dimension?
Definition: multilineargeometry.hh:156
Base::JacobianTransposed JacobianTransposed
Definition: multilineargeometry.hh:502
A unique label for each type of element that can occur in a grid.
CachedMultiLinearGeometry(const ReferenceElement &referenceElement, const CornerStorage &cornerStorage)
Definition: multilineargeometry.hh:506
static bool isPyramid(unsigned int topologyId, int dim, int codim=0) noexcept
check whether a pyramid construction was used to create a given codimension
Definition: type.hh:109
Definition: multilineargeometry.hh:444
Dune::ReferenceElement< ctype, mydimension > ReferenceElement
type of reference element
Definition: multilineargeometry.hh:211
template specifying the storage for the corners
Definition: multilineargeometry.hh:137
Impl::FieldMatrixHelper< ct > MatrixHelper
helper structure containing some matrix routines
Definition: multilineargeometry.hh:67
int corners() const
obtain number of corners of the corresponding reference element
Definition: multilineargeometry.hh:269
An implementation of the Geometry interface for affine geometries.
std::conditional< hasSingleGeometryType, std::integral_constant< unsigned int, Traits::template hasSingleGeometryType< mydimension >::topologyId >, unsigned int >::type TopologyId
Definition: multilineargeometry.hh:221
GlobalCoordinate global(const LocalCoordinate &local) const
evaluate the mapping
Definition: multilineargeometry.hh:535
GlobalCoordinate center() const
obtain the centroid of the mapping's image
Definition: multilineargeometry.hh:279
bool affine(JacobianTransposed &jacobianT) const
Definition: multilineargeometry.hh:416
static const bool v
Definition: multilineargeometry.hh:158
Base::MatrixHelper MatrixHelper
Definition: multilineargeometry.hh:489
MultiLinearGeometry(Dune::GeometryType gt, const Corners &corners)
constructor
Definition: multilineargeometry.hh:252
void setup(const JacobianTransposed &jt)
Definition: multilineargeometry.hh:450
Base::JacobianInverseTransposed JacobianInverseTransposed
Definition: multilineargeometry.hh:503
LocalCoordinate local(const GlobalCoordinate &globalCoord) const
evaluate the inverse mapping
Definition: multilineargeometry.hh:309
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
type of jacobian transposed
Definition: multilineargeometry.hh:208
static bool isPrism(unsigned int topologyId, int dim, int codim=0) noexcept
check whether a prism construction was used to create a given codimension
Definition: type.hh:127
MultiLinearGeometry(const ReferenceElement &refElement, const Corners &corners)
constructor
Definition: multilineargeometry.hh:236
GlobalCoordinate center() const
obtain the centroid of the mapping's image
Definition: multilineargeometry.hh:527
CachedMultiLinearGeometry(Dune::GeometryType gt, const CornerStorage &cornerStorage)
Definition: multilineargeometry.hh:514
LocalCoordinate local(const GlobalCoordinate &global) const
evaluate the inverse mapping
Definition: multilineargeometry.hh:559
FieldVector< ctype, mydimension > LocalCoordinate
type of local coordinates
Definition: multilineargeometry.hh:203
Dune::GeometryType type() const
obtain the name of the reference element
Definition: multilineargeometry.hh:266
friend const ReferenceElement & referenceElement(const MultiLinearGeometry &geometry)
Definition: multilineargeometry.hh:383
Base::LocalCoordinate LocalCoordinate
Definition: multilineargeometry.hh:499
Definition: affinegeometry.hh:39
ctype det() const
Definition: multilineargeometry.hh:460
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian's inverse
Definition: multilineargeometry.hh:671
Base::GlobalCoordinate GlobalCoordinate
Definition: multilineargeometry.hh:500
Traits::MatrixHelper MatrixHelper
Definition: multilineargeometry.hh:220
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian
Definition: multilineargeometry.hh:621
void setupDeterminant(const JacobianTransposed &jt)
Definition: multilineargeometry.hh:455
ctype volume() const
obtain the volume of the mapping's image
Definition: multilineargeometry.hh:351
GlobalCoordinate global(const LocalCoordinate &local) const
evaluate the mapping
Definition: multilineargeometry.hh:287
bool affine() const
is this mapping affine?
Definition: multilineargeometry.hh:522
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian
Definition: multilineargeometry.hh:365
generic geometry implementation based on corner coordinates
Definition: multilineargeometry.hh:189
TopologyId topologyId() const
Definition: multilineargeometry.hh:388