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>
59 static ct
tolerance () {
return ct( 16 ) * std::numeric_limits< ct >::epsilon(); }
125 template<
int mydim,
int cdim >
128 typedef std::vector< FieldVector< ct, cdim > >
Type;
147 static const bool v =
false;
177 template<
class ct,
int mydim,
int cdim,
class Traits = MultiLinearGeometryTraits< ct > >
214 static const bool hasSingleGeometryType = Traits::template hasSingleGeometryType< mydimension >::v;
218 typedef typename std::conditional< hasSingleGeometryType, std::integral_constant< unsigned int, Traits::template hasSingleGeometryType< mydimension >::topologyId >,
unsigned int >
::type TopologyId;
230 template<
class Corners >
246 template<
class Corners >
269 assert( (i >= 0) && (i <
corners()) );
270 return std::cref(corners_).get()[ i ];
286 auto cit = begin(std::cref(corners_).get());
288 global< false >(
topologyId(), std::integral_constant< int, mydimension >(), cit,
ctype( 1 ), local,
ctype( 1 ), y );
306 const ctype tolerance = Traits::tolerance();
309 const bool affineMapping = this->
affine();
313 const GlobalCoordinate dglobal = (*this).global( x ) - globalCoord;
314 const bool invertible =
315 MatrixHelper::template xTRightInvA< mydimension, coorddimension >(
jacobianTransposed( x ), dglobal, dx );
323 if ( affineMapping )
break;
324 }
while( dx.two_norm2() > tolerance );
344 return MatrixHelper::template sqrtDetAAT< mydimension, coorddimension >(
jacobianTransposed( local ) );
374 auto cit = begin(std::cref(corners_).get());
375 jacobianTransposed< false >(
topologyId(), std::integral_constant< int, mydimension >(), cit,
ctype( 1 ), local,
ctype( 1 ), jt );
401 return topologyId( std::integral_constant< bool, hasSingleGeometryType >() );
404 template<
bool add,
int dim,
class CornerIterator >
408 template<
bool add,
class CornerIterator >
413 template<
bool add,
int rows,
int dim,
class CornerIterator >
416 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt );
417 template<
bool add,
int rows,
class CornerIterator >
420 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt );
422 template<
int dim,
class CornerIterator >
424 template<
class CornerIterator >
431 auto cit = begin(std::cref(corners_).get());
432 return affine(
topologyId(), std::integral_constant< int, mydimension >(), cit, jacobianT );
439 static unsigned int toUnsignedInt(
unsigned int i) {
return i; }
440 template<
unsigned int v>
441 static unsigned int toUnsignedInt(std::integral_constant<unsigned int,v> ) {
return v; }
443 unsigned int topologyId ( std::integral_constant< bool, false > )
const {
return refElement().type().id(); }
446 typename Traits::template CornerStorage< mydimension, coorddimension >::Type corners_;
454 template<
class ct,
int mydim,
int cdim,
class Traits >
456 :
public FieldMatrix< ctype, coorddimension, mydimension >
458 typedef FieldMatrix< ctype, coorddimension, mydimension > Base;
463 detInv_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jt,
static_cast< Base &
>( *
this ) );
468 detInv_ = MatrixHelper::template sqrtDetAAT< mydimension, coorddimension >( jt );
492 template<
class ct,
int mydim,
int cdim,
class Traits = MultiLinearGeometryTraits< ct > >
517 template<
class CornerStorage >
520 affine_(
Base::
affine( jacobianTransposed_ ) ),
521 jacobianInverseTransposedComputed_( false ),
522 integrationElementComputed_( false )
525 template<
class CornerStorage >
527 :
Base( gt, cornerStorage ),
528 affine_(
Base::
affine( jacobianTransposed_ ) ),
529 jacobianInverseTransposedComputed_( false ),
530 integrationElementComputed_( false )
552 jacobianTransposed_.umtv( local,
global );
576 if( jacobianInverseTransposedComputed_ )
577 jacobianInverseTransposed_.mtv(
global -
corner( 0 ), local );
579 MatrixHelper::template xTRightInvA< mydimension, coorddimension >( jacobianTransposed_,
global -
corner( 0 ), local );
583 return Base::local(
global );
604 if( !integrationElementComputed_ )
607 integrationElementComputed_ =
true;
609 return jacobianInverseTransposed_.
detInv();
636 return jacobianTransposed_;
651 if( !jacobianInverseTransposedComputed_ )
653 jacobianInverseTransposed_.
setup( jacobianTransposed_ );
654 jacobianInverseTransposedComputed_ =
true;
655 integrationElementComputed_ =
true;
657 return jacobianInverseTransposed_;
670 mutable bool affine_ : 1;
672 mutable bool jacobianInverseTransposedComputed_ : 1;
673 mutable bool integrationElementComputed_ : 1;
681 template<
class ct,
int mydim,
int cdim,
class Traits >
682 inline typename MultiLinearGeometry< ct, mydim, cdim, Traits >::JacobianInverseTransposed
686 jit.
setup( jacobianTransposed( local ) );
691 template<
class ct,
int mydim,
int cdim,
class Traits >
692 template<
bool add,
int dim,
class CornerIterator >
698 const ctype xn = df*x[ dim-1 ];
701 if( Impl::isPrism( toUnsignedInt(topologyId), mydimension, mydimension-dim ) )
704 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, rf*cxn, y );
706 global< true >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, rf*xn, y );
710 assert( Impl::isPyramid( toUnsignedInt(topologyId), mydimension, mydimension-dim ) );
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 );
715 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x,
ctype( 0 ), y );
717 y.axpy( rf*xn, *cit );
722 template<
class ct,
int mydim,
int cdim,
class Traits >
723 template<
bool add,
class CornerIterator >
731 for(
int i = 0; i < coorddimension; ++i )
732 y[ i ] = (add ? y[ i ] + rf*origin[ i ] : rf*origin[ i ]);
736 template<
class ct,
int mydim,
int cdim,
class Traits >
737 template<
bool add,
int rows,
int dim,
class CornerIterator >
741 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt )
743 assert( rows >= dim );
745 const ctype xn = df*x[ dim-1 ];
749 if( Impl::isPrism( toUnsignedInt(topologyId), mydimension, mydimension-dim ) )
752 jacobianTransposed< add >( topologyId, std::integral_constant< int, dim-1 >(), cit2, df, x, rf*cxn, jt );
754 jacobianTransposed< true >( topologyId, std::integral_constant< int, dim-1 >(), cit2, df, x, rf*xn, jt );
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 ] );
761 assert( Impl::isPyramid( toUnsignedInt(topologyId), mydimension, mydimension-dim ) );
803 ctype dfcxn = (cxn > Traits::tolerance() || cxn < -Traits::tolerance()) ?
ctype(df / cxn) :
ctype(0);
808 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, dfcxn, x, -rf, jt[ dim-1 ] );
810 jt[ dim-1 ].axpy( rf, *cit );
815 FieldMatrix<
ctype, dim-1, coorddimension > jt2;
817 jacobianTransposed< false >( topologyId, std::integral_constant< int, dim-1 >(), cit2, dfcxn, x, rf, jt2 );
821 for(
int j = 0; j < dim-1; ++j )
824 jt[ dim-1 ].axpy( dfcxn*x[ j ], jt2[ j ] );
830 jacobianTransposed< false >( topologyId, std::integral_constant< int, dim-1 >(), cit2, dfcxn, x, rf, jt );
832 for(
int j = 0; j < dim-1; ++j )
833 jt[ dim-1 ].axpy( dfcxn*x[ j ], jt[ j ] );
838 template<
class ct,
int mydim,
int cdim,
class Traits >
839 template<
bool add,
int rows,
class CornerIterator >
843 const ctype &, FieldMatrix< ctype, rows, cdim > & )
850 template<
class ct,
int mydim,
int cdim,
class Traits >
851 template<
int dim,
class CornerIterator >
856 if( !affine( topologyId, std::integral_constant< int, dim-1 >(), cit, jt ) )
860 if( Impl::isPrism( toUnsignedInt(topologyId), mydimension, mydimension-dim ) )
863 if( !affine( topologyId, std::integral_constant< int, dim-1 >(), cit, jtTop ) )
868 for(
int i = 0; i < dim-1; ++i )
869 norm += (jtTop[ i ] - jt[ i ]).two_norm2();
870 if( norm >= Traits::tolerance() )
875 jt[ dim-1 ] = orgTop - orgBottom;
879 template<
class ct,
int mydim,
int cdim,
class Traits >
880 template<
class CornerIterator >
An implementation of the Geometry interface for affine geometries.
A unique label for each type of element that can occur in a grid.
Definition: affinegeometry.hh:19
Class providing access to the singletons of the reference elements.
Definition: referenceelements.hh:168
typename Container::ReferenceElement ReferenceElement
The reference element type.
Definition: referenceelements.hh:186
default traits class for MultiLinearGeometry
Definition: multilineargeometry.hh:37
Impl::FieldMatrixHelper< ct > MatrixHelper
helper structure containing some matrix routines
Definition: multilineargeometry.hh:56
static ct tolerance()
tolerance to numerical algorithms
Definition: multilineargeometry.hh:59
template specifying the storage for the corners
Definition: multilineargeometry.hh:127
std::vector< FieldVector< ct, cdim > > Type
Definition: multilineargeometry.hh:128
will there be only one geometry type for a dimension?
Definition: multilineargeometry.hh:146
static const unsigned int topologyId
Definition: multilineargeometry.hh:148
static const bool v
Definition: multilineargeometry.hh:147
generic geometry implementation based on corner coordinates
Definition: multilineargeometry.hh:179
static void global(TopologyId topologyId, std::integral_constant< int, dim >, CornerIterator &cit, const ctype &df, const LocalCoordinate &x, const ctype &rf, GlobalCoordinate &y)
Definition: multilineargeometry.hh:694
static const int mydimension
geometry dimension
Definition: multilineargeometry.hh:187
Dune::GeometryType type() const
obtain the name of the reference element
Definition: multilineargeometry.hh:261
Traits::MatrixHelper MatrixHelper
Definition: multilineargeometry.hh:217
FieldVector< ctype, coorddimension > GlobalCoordinate
type of global coordinates
Definition: multilineargeometry.hh:194
static void jacobianTransposed(TopologyId topologyId, std::integral_constant< int, 0 >, CornerIterator &cit, const ctype &df, const LocalCoordinate &x, const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt)
Definition: multilineargeometry.hh:841
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian
Definition: multilineargeometry.hh:369
ReferenceElement refElement() const
Definition: multilineargeometry.hh:394
GlobalCoordinate global(const LocalCoordinate &local) const
evaluate the mapping
Definition: multilineargeometry.hh:282
GlobalCoordinate center() const
obtain the centroid of the mapping's image
Definition: multilineargeometry.hh:274
GlobalCoordinate corner(int i) const
obtain coordinates of the i-th corner
Definition: multilineargeometry.hh:267
Dune::ReferenceElements< ctype, mydimension > ReferenceElements
Definition: multilineargeometry.hh:202
ct ctype
coordinate type
Definition: multilineargeometry.hh:184
static void jacobianTransposed(TopologyId topologyId, std::integral_constant< int, dim >, CornerIterator &cit, const ctype &df, const LocalCoordinate &x, const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt)
Definition: multilineargeometry.hh:739
static const int coorddimension
coordinate dimension
Definition: multilineargeometry.hh:189
int corners() const
obtain number of corners of the corresponding reference element
Definition: multilineargeometry.hh:264
TopologyId topologyId() const
Definition: multilineargeometry.hh:399
friend ReferenceElement referenceElement(const MultiLinearGeometry &geometry)
Definition: multilineargeometry.hh:387
static void global(TopologyId topologyId, std::integral_constant< int, 0 >, CornerIterator &cit, const ctype &df, const LocalCoordinate &x, const ctype &rf, GlobalCoordinate &y)
Definition: multilineargeometry.hh:725
Volume volume() const
obtain the volume of the mapping's image
Definition: multilineargeometry.hh:355
FieldVector< ctype, mydimension > LocalCoordinate
type of local coordinates
Definition: multilineargeometry.hh:192
MultiLinearGeometry(const ReferenceElement &refElement, const Corners &corners)
constructor
Definition: multilineargeometry.hh:231
static bool affine(TopologyId topologyId, std::integral_constant< int, 0 >, CornerIterator &cit, JacobianTransposed &jt)
Definition: multilineargeometry.hh:882
std::conditional< hasSingleGeometryType, std::integral_constant< unsignedint, Traits::templatehasSingleGeometryType< mydimension >::topologyId >, unsignedint >::type TopologyId
Definition: multilineargeometry.hh:218
ctype Volume
type of volume
Definition: multilineargeometry.hh:196
static bool affine(TopologyId topologyId, std::integral_constant< int, dim >, CornerIterator &cit, JacobianTransposed &jt)
Definition: multilineargeometry.hh:853
MultiLinearGeometry(Dune::GeometryType gt, const Corners &corners)
constructor
Definition: multilineargeometry.hh:247
ctype integrationElement(const LocalCoordinate &local) const
obtain the integration element
Definition: multilineargeometry.hh:342
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
type of jacobian transposed
Definition: multilineargeometry.hh:199
ReferenceElements::ReferenceElement ReferenceElement
type of reference element
Definition: multilineargeometry.hh:211
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian's inverse
Definition: multilineargeometry.hh:683
bool affine() const
is this mapping affine?
Definition: multilineargeometry.hh:254
bool affine(JacobianTransposed &jacobianT) const
Definition: multilineargeometry.hh:427
Definition: multilineargeometry.hh:457
void setup(const JacobianTransposed &jt)
Definition: multilineargeometry.hh:461
ctype det() const
Definition: multilineargeometry.hh:471
ctype detInv() const
Definition: multilineargeometry.hh:472
void setupDeterminant(const JacobianTransposed &jt)
Definition: multilineargeometry.hh:466
Implement a MultiLinearGeometry with additional caching.
Definition: multilineargeometry.hh:495
GlobalCoordinate global(const LocalCoordinate &local) const
evaluate the mapping
Definition: multilineargeometry.hh:547
Base::ReferenceElement ReferenceElement
Definition: multilineargeometry.hh:503
bool affine() const
is this mapping affine?
Definition: multilineargeometry.hh:534
CachedMultiLinearGeometry(const ReferenceElement &referenceElement, const CornerStorage &cornerStorage)
Definition: multilineargeometry.hh:518
ReferenceElement refElement() const
Definition: multilineargeometry.hh:394
Base::MatrixHelper MatrixHelper
Definition: multilineargeometry.hh:500
Base::LocalCoordinate LocalCoordinate
Definition: multilineargeometry.hh:510
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian
Definition: multilineargeometry.hh:633
GlobalCoordinate corner(int i) const
obtain coordinates of the i-th corner
Definition: multilineargeometry.hh:267
Volume volume() const
obtain the volume of the mapping's image
Definition: multilineargeometry.hh:616
CachedMultiLinearGeometry(Dune::GeometryType gt, const CornerStorage &cornerStorage)
Definition: multilineargeometry.hh:526
ctype integrationElement(const LocalCoordinate &local) const
obtain the integration element
Definition: multilineargeometry.hh:600
Base::ctype ctype
Definition: multilineargeometry.hh:505
Base::JacobianInverseTransposed JacobianInverseTransposed
Definition: multilineargeometry.hh:515
Base::JacobianTransposed JacobianTransposed
Definition: multilineargeometry.hh:514
Base::Volume Volume
Definition: multilineargeometry.hh:512
Base::GlobalCoordinate GlobalCoordinate
Definition: multilineargeometry.hh:511
GlobalCoordinate center() const
obtain the centroid of the mapping's image
Definition: multilineargeometry.hh:539
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian's inverse
Definition: multilineargeometry.hh:647
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:123