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>
27 template<
class ctype,
int dim >
28 class ReferenceElement;
30 template<
class ctype,
int dim >
31 struct ReferenceElements;
71 static ct
tolerance () {
return ct( 16 ) * std::numeric_limits< ct >::epsilon(); }
137 template<
int mydim,
int cdim >
140 typedef std::vector< FieldVector< ct, cdim > >
Type;
159 static const bool v =
false;
189 template<
class ct,
int mydim,
int cdim,
class Traits = MultiLinearGeometryTraits< ct > >
218 static const bool hasSingleGeometryType = Traits::template hasSingleGeometryType< mydimension >::v;
222 typedef typename std::conditional< hasSingleGeometryType, std::integral_constant< unsigned int, Traits::template hasSingleGeometryType< mydimension >::topologyId >,
unsigned int >
::type TopologyId;
236 template<
class Corners >
252 template<
class Corners >
275 assert( (i >= 0) && (i <
corners()) );
276 return std::cref(corners_).get()[ i ];
292 auto cit = begin(std::cref(corners_).get());
312 const ctype tolerance = Traits::tolerance();
319 MatrixHelper::template xTRightInvA< mydimension, coorddimension >(
jacobianTransposed( x ), dglobal, dx );
321 }
while( dx.two_norm2() > tolerance );
371 auto cit = begin(std::cref(corners_).get());
372 jacobianTransposed< false >(
topologyId(), std::integral_constant< int, mydimension >(), cit,
ctype( 1 ),
local,
ctype( 1 ), jt );
391 return topologyId( std::integral_constant< bool, hasSingleGeometryType >() );
394 template<
bool add,
int dim,
class CornerIterator >
398 template<
bool add,
class CornerIterator >
403 template<
bool add,
int rows,
int dim,
class CornerIterator >
406 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt );
407 template<
bool add,
int rows,
class CornerIterator >
410 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt );
412 template<
int dim,
class CornerIterator >
414 template<
class CornerIterator >
421 auto cit = begin(std::cref(corners_).get());
422 return affine(
topologyId(), std::integral_constant< int, mydimension >(), cit, jacobianT );
429 static unsigned int toUnsignedInt(
unsigned int i) {
return i; }
430 template<
unsigned int v>
431 static unsigned int toUnsignedInt(std::integral_constant<unsigned int,v> i) {
return v; }
436 typename Traits::template CornerStorage< mydimension, coorddimension >::Type corners_;
444 template<
class ct,
int mydim,
int cdim,
class Traits >
446 :
public FieldMatrix< ctype, coorddimension, mydimension >
448 typedef FieldMatrix< ctype, coorddimension, mydimension > Base;
453 detInv_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jt,
static_cast< Base &
>( *
this ) );
458 detInv_ = MatrixHelper::template sqrtDetAAT< mydimension, coorddimension >( jt );
482 template<
class ct,
int mydim,
int cdim,
class Traits = MultiLinearGeometryTraits< ct > >
506 template<
class CornerStorage >
509 affine_(
Base::
affine( jacobianTransposed_ ) ),
510 jacobianInverseTransposedComputed_( false ),
511 integrationElementComputed_( false )
514 template<
class CornerStorage >
516 :
Base( gt, cornerStorage ),
517 affine_(
Base::
affine( jacobianTransposed_ ) ),
518 jacobianInverseTransposedComputed_( false ),
519 integrationElementComputed_( false )
565 if( jacobianInverseTransposedComputed_ )
568 MatrixHelper::template xTRightInvA< mydimension, coorddimension >( jacobianTransposed_,
global -
corner( 0 ),
local );
593 if( !integrationElementComputed_ )
596 integrationElementComputed_ =
true;
598 return jacobianInverseTransposed_.
detInv();
625 return jacobianTransposed_;
640 if( !jacobianInverseTransposedComputed_ )
642 jacobianInverseTransposed_.
setup( jacobianTransposed_ );
643 jacobianInverseTransposedComputed_ =
true;
644 integrationElementComputed_ =
true;
646 return jacobianInverseTransposed_;
659 mutable bool affine_ : 1;
661 mutable bool jacobianInverseTransposedComputed_ : 1;
662 mutable bool integrationElementComputed_ : 1;
670 template<
class ct,
int mydim,
int cdim,
class Traits >
671 inline typename MultiLinearGeometry< ct, mydim, cdim, Traits >::JacobianInverseTransposed
675 jit.
setup( jacobianTransposed( local ) );
680 template<
class ct,
int mydim,
int cdim,
class Traits >
681 template<
bool add,
int dim,
class CornerIterator >
687 const ctype xn = df*x[ dim-1 ];
693 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, rf*cxn, y );
695 global< true >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, rf*xn, y );
701 if( cxn > Traits::tolerance() || cxn < -Traits::tolerance() )
702 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df/cxn, x, rf*cxn, y );
704 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x,
ctype( 0 ), y );
706 y.axpy( rf*xn, *cit );
711 template<
class ct,
int mydim,
int cdim,
class Traits >
712 template<
bool add,
class CornerIterator >
720 for(
int i = 0; i < coorddimension; ++i )
721 y[ i ] = (add ? y[ i ] + rf*origin[ i ] : rf*origin[ i ]);
725 template<
class ct,
int mydim,
int cdim,
class Traits >
726 template<
bool add,
int rows,
int dim,
class CornerIterator >
730 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt )
732 assert( rows >= dim );
734 const ctype xn = df*x[ dim-1 ];
741 jacobianTransposed< add >( topologyId, std::integral_constant< int, dim-1 >(), cit2, df, x, rf*cxn, jt );
743 jacobianTransposed< true >( topologyId, std::integral_constant< int, dim-1 >(), cit2, df, x, rf*xn, jt );
745 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, -rf, jt[ dim-1 ] );
746 global< true >( topologyId, std::integral_constant< int, dim-1 >(), cit, df, x, rf, jt[ dim-1 ] );
792 ctype dfcxn = (cxn > Traits::tolerance() || cxn < -Traits::tolerance()) ?
ctype(df / cxn) :
ctype(0);
797 global< add >( topologyId, std::integral_constant< int, dim-1 >(), cit, dfcxn, x, -rf, jt[ dim-1 ] );
799 jt[ dim-1 ].axpy( rf, *cit );
804 FieldMatrix<
ctype, dim-1, coorddimension > jt2;
806 jacobianTransposed< false >( topologyId, std::integral_constant< int, dim-1 >(), cit2, dfcxn, x, rf, jt2 );
810 for(
int j = 0; j < dim-1; ++j )
813 jt[ dim-1 ].axpy( dfcxn*x[ j ], jt2[ j ] );
819 jacobianTransposed< false >( topologyId, std::integral_constant< int, dim-1 >(), cit2, dfcxn, x, rf, jt );
821 for(
int j = 0; j < dim-1; ++j )
822 jt[ dim-1 ].axpy( dfcxn*x[ j ], jt[ j ] );
827 template<
class ct,
int mydim,
int cdim,
class Traits >
828 template<
bool add,
int rows,
class CornerIterator >
832 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt )
839 template<
class ct,
int mydim,
int cdim,
class Traits >
840 template<
int dim,
class CornerIterator >
845 if( !affine( topologyId, std::integral_constant< int, dim-1 >(), cit, jt ) )
852 if( !affine( topologyId, std::integral_constant< int, dim-1 >(), cit, jtTop ) )
857 for(
int i = 0; i < dim-1; ++i )
858 norm += (jtTop[ i ] - jt[ i ]).two_norm2();
859 if( norm >= Traits::tolerance() )
864 jt[ dim-1 ] = orgTop - orgBottom;
868 template<
class ct,
int mydim,
int cdim,
class Traits >
869 template<
class CornerIterator >
A unique label for each type of element that can occur in a grid.
Definition affinegeometry.hh:19
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
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
This class provides access to geometric and topological properties of a reference element.
Definition referenceelements.hh:53
int size(int c) const
number of subentities of codimension c
Definition referenceelements.hh:80
const FieldVector< ctype, dim > & position(int i, int c) const
position of the barycenter of entity (i,c)
Definition referenceelements.hh:148
ctype volume() const
obtain the volume of the reference element
Definition referenceelements.hh:186
const GeometryType & type(int i, int c) const
obtain the type of subentity (i,c)
Definition referenceelements.hh:130
Class providing access to the singletons of the reference elements.
Definition referenceelements.hh:455
Definition matrixhelper.hh:36
default traits class for MultiLinearGeometry
Definition multilineargeometry.hh:49
static ct tolerance()
tolerance to numerical algorithms
Definition multilineargeometry.hh:71
GenericGeometry::MatrixHelper< GenericGeometry::DuneCoordTraits< ct > > MatrixHelper
helper structure containing some matrix routines
Definition multilineargeometry.hh:68
template specifying the storage for the corners
Definition multilineargeometry.hh:139
std::vector< FieldVector< ct, cdim > > Type
Definition multilineargeometry.hh:140
will there be only one geometry type for a dimension?
Definition multilineargeometry.hh:158
static const unsigned int topologyId
Definition multilineargeometry.hh:160
static const bool v
Definition multilineargeometry.hh:159
generic geometry implementation based on corner coordinates
Definition multilineargeometry.hh:191
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:683
static const int mydimension
geometry dimension
Definition multilineargeometry.hh:199
Dune::ReferenceElement< ctype, mydimension > ReferenceElement
type of reference element
Definition multilineargeometry.hh:215
Dune::GeometryType type() const
obtain the name of the reference element
Definition multilineargeometry.hh:267
Traits::MatrixHelper MatrixHelper
Definition multilineargeometry.hh:221
FieldVector< ctype, coorddimension > GlobalCoordinate
type of global coordinates
Definition multilineargeometry.hh:206
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:830
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian
Definition multilineargeometry.hh:366
GlobalCoordinate global(const LocalCoordinate &local) const
evaluate the mapping
Definition multilineargeometry.hh:288
GlobalCoordinate center() const
obtain the centroid of the mapping's image
Definition multilineargeometry.hh:280
GlobalCoordinate corner(int i) const
obtain coordinates of the i-th corner
Definition multilineargeometry.hh:273
Dune::ReferenceElements< ctype, mydimension > ReferenceElements
Definition multilineargeometry.hh:224
ct ctype
coordinate type
Definition multilineargeometry.hh:196
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:728
static const int coorddimension
coordinate dimension
Definition multilineargeometry.hh:201
int corners() const
obtain number of corners of the corresponding reference element
Definition multilineargeometry.hh:270
TopologyId topologyId() const
Definition multilineargeometry.hh:389
LocalCoordinate local(const GlobalCoordinate &globalCoord) const
evaluate the inverse mapping
Definition multilineargeometry.hh:310
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:714
FieldVector< ctype, mydimension > LocalCoordinate
type of local coordinates
Definition multilineargeometry.hh:204
ctype volume() const
obtain the volume of the mapping's image
Definition multilineargeometry.hh:352
MultiLinearGeometry(const ReferenceElement &refElement, const Corners &corners)
constructor
Definition multilineargeometry.hh:237
static bool affine(TopologyId topologyId, std::integral_constant< int, 0 >, CornerIterator &cit, JacobianTransposed &jt)
Definition multilineargeometry.hh:871
std::conditional< hasSingleGeometryType, std::integral_constant< unsignedint, Traits::templatehasSingleGeometryType< mydimension >::topologyId >, unsignedint >::type TopologyId
Definition multilineargeometry.hh:222
static bool affine(TopologyId topologyId, std::integral_constant< int, dim >, CornerIterator &cit, JacobianTransposed &jt)
Definition multilineargeometry.hh:842
MultiLinearGeometry(Dune::GeometryType gt, const Corners &corners)
constructor
Definition multilineargeometry.hh:253
ctype integrationElement(const LocalCoordinate &local) const
obtain the integration element
Definition multilineargeometry.hh:339
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
type of jacobian transposed
Definition multilineargeometry.hh:209
friend const ReferenceElement & referenceElement(const MultiLinearGeometry &geometry)
Definition multilineargeometry.hh:384
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian's inverse
Definition multilineargeometry.hh:672
bool affine() const
is this mapping affine?
Definition multilineargeometry.hh:260
const ReferenceElement & refElement() const
Definition multilineargeometry.hh:387
bool affine(JacobianTransposed &jacobianT) const
Definition multilineargeometry.hh:417
Definition multilineargeometry.hh:447
void setup(const JacobianTransposed &jt)
Definition multilineargeometry.hh:451
ctype det() const
Definition multilineargeometry.hh:461
ctype detInv() const
Definition multilineargeometry.hh:462
void setupDeterminant(const JacobianTransposed &jt)
Definition multilineargeometry.hh:456
Implement a MultiLinearGeometry with additional caching.
Definition multilineargeometry.hh:485
GlobalCoordinate global(const LocalCoordinate &local) const
evaluate the mapping
Definition multilineargeometry.hh:536
Base::ReferenceElement ReferenceElement
Definition multilineargeometry.hh:493
bool affine() const
is this mapping affine?
Definition multilineargeometry.hh:523
CachedMultiLinearGeometry(const ReferenceElement &referenceElement, const CornerStorage &cornerStorage)
Definition multilineargeometry.hh:507
LocalCoordinate local(const GlobalCoordinate &global) const
evaluate the inverse mapping
Definition multilineargeometry.hh:560
Base::MatrixHelper MatrixHelper
Definition multilineargeometry.hh:490
Base::LocalCoordinate LocalCoordinate
Definition multilineargeometry.hh:500
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian
Definition multilineargeometry.hh:622
GlobalCoordinate corner(int i) const
obtain coordinates of the i-th corner
Definition multilineargeometry.hh:273
CachedMultiLinearGeometry(Dune::GeometryType gt, const CornerStorage &cornerStorage)
Definition multilineargeometry.hh:515
ctype volume() const
obtain the volume of the mapping's image
Definition multilineargeometry.hh:605
ctype integrationElement(const LocalCoordinate &local) const
obtain the integration element
Definition multilineargeometry.hh:589
Base::ctype ctype
Definition multilineargeometry.hh:495
Base::JacobianInverseTransposed JacobianInverseTransposed
Definition multilineargeometry.hh:504
Base::JacobianTransposed JacobianTransposed
Definition multilineargeometry.hh:503
const ReferenceElement & refElement() const
Definition multilineargeometry.hh:387
Base::GlobalCoordinate GlobalCoordinate
Definition multilineargeometry.hh:501
GlobalCoordinate center() const
obtain the centroid of the mapping's image
Definition multilineargeometry.hh:528
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian's inverse
Definition multilineargeometry.hh:636
Unique label for each type of entities that can occur in DUNE grids.
Definition type.hh:25
unsigned int id() const
Return the topology id the type.
Definition type.hh:326