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;
149 template<
class ct,
int mydim,
int cdim,
class Traits = MultiLinearGeometryTraits< ct > >
187 static const bool hasSingleGeometryType = Traits::template hasSingleGeometryType< mydimension >::v;
191 typedef typename conditional< hasSingleGeometryType, integral_constant< unsigned int, Traits::template hasSingleGeometryType< mydimension >::topologyId >,
unsigned int >
::type TopologyId;
196 typedef typename Traits::template CornerStorage< mydimension, coorddimension >::Type::const_iterator CornerIterator;
208 template<
class Corners >
211 : refElement_( &refElement ),
224 template<
class Corners >
234 CornerIterator cit = corners_.begin();
247 assert( (i >= 0) && (i <
corners()) );
248 return corners_[ i ];
262 CornerIterator cit = corners_.begin();
281 const ctype tolerance = Traits::tolerance();
288 MatrixHelper::template xTRightInvA< mydimension, coorddimension >(
jacobianTransposed( x ), dglobal, dx );
291 }
while( dx.two_norm2() > tolerance );
311 return MatrixHelper::template sqrtDetAAT< mydimension, coorddimension >(
jacobianTransposed( local ) );
338 CornerIterator cit = corners_.begin();
356 return topologyId( integral_constant< bool, hasSingleGeometryType >() );
359 template<
bool add,
int dim >
368 template<
bool add,
int rows,
int dim >
371 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt );
372 template<
bool add,
int rows >
375 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt );
390 typename Traits::template CornerStorage< mydimension, coorddimension >::Type corners_;
398 template<
class ct,
int mydim,
int cdim,
class Traits >
400 :
public FieldMatrix< ctype, coorddimension, mydimension >
402 typedef FieldMatrix< ctype, coorddimension, mydimension > Base;
407 detInv_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jt,
static_cast< Base &
>( *this ) );
412 detInv_ = MatrixHelper::template sqrtDetAAT< mydimension, coorddimension >( jt );
436 template<
class ct,
int mydim,
int cdim,
class Traits = MultiLinearGeometryTraits< ct > >
460 template<
class CornerStorage >
462 :
Base( refElement, cornerStorage ),
464 jacobianInverseTransposedComputed_( false ),
465 integrationElementComputed_( false )
468 template<
class CornerStorage >
470 :
Base( gt, cornerStorage ),
472 jacobianInverseTransposedComputed_( false ),
473 integrationElementComputed_( false )
518 if( jacobianInverseTransposedComputed_ )
546 if( !integrationElementComputed_ )
549 integrationElementComputed_ =
true;
593 if( !jacobianInverseTransposedComputed_ )
596 jacobianInverseTransposedComputed_ =
true;
597 integrationElementComputed_ =
true;
612 mutable bool affine_ : 1;
614 mutable bool jacobianInverseTransposedComputed_ : 1;
615 mutable bool integrationElementComputed_ : 1;
623 template<
class ct,
int mydim,
int cdim,
class Traits >
624 inline const typename MultiLinearGeometry< ct, mydim, cdim, Traits >::JacobianInverseTransposed &
628 jacobianInverseTransposed_.
setup( jacobianTransposed( local ) );
629 return jacobianInverseTransposed_;
633 template<
class ct,
int mydim,
int cdim,
class Traits >
634 template<
bool add,
int dim >
640 const ctype xn = df*x[ dim-1 ];
642 assert( (xn > -Traits::tolerance()) && (cxn > -Traits::tolerance()) );
647 global< add >( topologyId, integral_constant< int, dim-1 >(), cit, df, x, rf*cxn, y );
649 global< true >( topologyId, integral_constant< int, dim-1 >(), cit, df, x, rf*xn, y );
655 if( cxn > Traits::tolerance() )
656 global< add >( topologyId, integral_constant< int, dim-1 >(), cit, df/cxn, x, rf*cxn, y );
658 global< add >( topologyId, integral_constant< int, dim-1 >(), cit, df, x,
ctype( 0 ), y );
660 y.axpy( rf*xn, *cit );
665 template<
class ct,
int mydim,
int cdim,
class Traits >
674 for(
int i = 0; i < coorddimension; ++i )
675 y[ i ] = (add ? y[ i ] + rf*origin[ i ] : rf*origin[ i ]);
679 template<
class ct,
int mydim,
int cdim,
class Traits >
680 template<
bool add,
int rows,
int dim >
684 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt )
686 assert( rows >= dim );
688 const ctype xn = df*x[ dim-1 ];
690 assert( (xn > -Traits::tolerance()) && (cxn > -Traits::tolerance()) );
692 CornerIterator cit2( cit );
696 jacobianTransposed< add >( topologyId, integral_constant< int, dim-1 >(), cit2, df, x, rf*cxn, jt );
698 jacobianTransposed< true >( topologyId, integral_constant< int, dim-1 >(), cit2, df, x, rf*xn, jt );
700 global< add >( topologyId, integral_constant< int, dim-1 >(), cit, df, x, -rf, jt[ dim-1 ] );
701 global< true >( topologyId, integral_constant< int, dim-1 >(), cit, df, x, rf, jt[ dim-1 ] );
707 global< add >( topologyId, integral_constant< int, dim-1 >(), cit, df/cxn, x, -rf, jt[ dim-1 ] );
708 jt[ dim-1 ].axpy( rf, *cit );
713 FieldMatrix<
ctype, dim-1, coorddimension > jt2;
714 jacobianTransposed< false >( topologyId, integral_constant< int, dim-1 >(), cit2, df/cxn, x, rf, jt2 );
715 for(
int j = 0; j < dim-1; ++j )
718 jt[ dim-1 ].axpy( (df/cxn)*x[ j ], jt2[ j ] );
723 jacobianTransposed< false >( topologyId, integral_constant< int, dim-1 >(), cit2, df/cxn, x, rf, jt );
724 for(
int j = 0; j < dim-1; ++j )
725 jt[ dim-1 ].axpy( (df/cxn)*x[ j ], jt[ j ] );
730 template<
class ct,
int mydim,
int cdim,
class Traits >
731 template<
bool add,
int rows >
735 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt )
742 template<
class ct,
int mydim,
int cdim,
class Traits >
748 if( !affine( topologyId, integral_constant< int, dim-1 >(), cit, jt ) )
755 if( !affine( topologyId, integral_constant< int, dim-1 >(), cit, jtTop ) )
760 for(
int i = 0; i < dim-1; ++i )
761 norm += (jtTop[ i ] - jt[ i ]).two_norm2();
762 if( norm >= Traits::tolerance() )
767 jt[ dim-1 ] = orgTop - orgBottom;
771 template<
class ct,
int mydim,
int cdim,
class Traits >
781 #endif // #ifndef DUNE_GEOMETRY_MULTILINEARGEOMETRY_HH
Traits::MatrixHelper MatrixHelper
Definition: multilineargeometry.hh:190
CachedMultiLinearGeometry(const ReferenceElement &refElement, const CornerStorage &cornerStorage)
Definition: multilineargeometry.hh:461
Dune::GeometryType type() const
obtain the name of the reference element
Definition: multilineargeometry.hh:239
JacobianInverseTransposed jacobianInverseTransposed_
Definition: multilineargeometry.hh:386
static const bool v
Definition: multilineargeometry.hh:116
will there be only one geometry type for a dimension?
Definition: multilineargeometry.hh:114
TopologyId topologyId() const
Definition: multilineargeometry.hh:354
MultiLinearGeometry(const ReferenceElement &refElement, const Corners &corners)
constructor
Definition: multilineargeometry.hh:209
Definition: matrixhelper.hh:33
LocalCoordinate local(const GlobalCoordinate &global) const
evaluate the inverse mapping
Definition: multilineargeometry.hh:513
Implement a MultiLinearGeometry with additional caching.
Definition: multilineargeometry.hh:437
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
type of jacobian transposed
Definition: multilineargeometry.hh:173
Base::ctype ctype
Definition: multilineargeometry.hh:449
GlobalCoordinate center() const
obtain the centroid of the mapping's image
Definition: multilineargeometry.hh:252
bool affine() const
is this mapping affine?
Definition: multilineargeometry.hh:477
unsigned int topologyId(integral_constant< bool, false >) const
Definition: multilineargeometry.hh:383
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:169
default traits class for MultiLinearGeometry
Definition: multilineargeometry.hh:46
int size(int c) const
number of subentities of codimension c
Definition: referenceelements.hh:86
Dune::ReferenceElement< ctype, mydimension > ReferenceElement
type of reference element
Definition: multilineargeometry.hh:184
A unique label for each type of element that can occur in a grid.
ctype volume() const
obtain the volume of the mapping's image
Definition: multilineargeometry.hh:558
ctype volume() const
obtain the volume of the reference element
Definition: referenceelements.hh:280
void setupDeterminant(const JacobianTransposed &jt)
Definition: multilineargeometry.hh:410
const ReferenceElement & refElement() const
Definition: multilineargeometry.hh:352
FieldVector< ctype, mydimension > LocalCoordinate
type of user data
Definition: multilineargeometry.hh:168
static const unsigned int topologyId
Definition: multilineargeometry.hh:117
This class provides access to geometric and topological properties of a reference element...
Definition: affinegeometry.hh:25
const FieldVector< ctype, dim > & position(int i, int c) const
position of the barycenter of entity (i,c)
Definition: referenceelements.hh:154
static const int coorddimension
coordinate dimension
Definition: multilineargeometry.hh:161
Base::JacobianTransposed JacobianTransposed
Definition: multilineargeometry.hh:457
static ct tolerance()
tolerance to numerical algorithms
Definition: multilineargeometry.hh:69
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:151
void setup(const JacobianTransposed &jt)
Definition: multilineargeometry.hh:405
ct ctype
coordinate type
Definition: multilineargeometry.hh:156
CachedMultiLinearGeometry(Dune::GeometryType gt, const CornerStorage &cornerStorage)
Definition: multilineargeometry.hh:469
const JacobianTransposed & jacobianTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian
Definition: multilineargeometry.hh:336
const JacobianInverseTransposed & jacobianInverseTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian's inverse
Definition: multilineargeometry.hh:626
Base::MatrixHelper MatrixHelper
Definition: multilineargeometry.hh:444
GlobalCoordinate corner(int i) const
obtain coordinates of the i-th corner
Definition: multilineargeometry.hh:245
GenericGeometry::MatrixHelper< GenericGeometry::DuneCoordTraits< ct > > MatrixHelper
helper structure containing some matrix routines
Definition: multilineargeometry.hh:66
unsigned int id() const
Return the topology id the type.
Definition: type.hh:327
Class providing access to the singletons of the reference elements. Special methods are available for...
Definition: affinegeometry.hh:28
GlobalCoordinate center() const
obtain the centroid of the mapping's image
Definition: multilineargeometry.hh:482
MultiLinearGeometry(Dune::GeometryType gt, const Corners &corners)
constructor
Definition: multilineargeometry.hh:225
Dune::ReferenceElements< ctype, mydimension > ReferenceElements
Definition: multilineargeometry.hh:193
GlobalCoordinate global(const LocalCoordinate &local) const
evaluate the mapping
Definition: multilineargeometry.hh:490
int corners() const
obtain number of corners of the corresponding reference element
Definition: multilineargeometry.hh:242
FieldVector< ctype, coorddimension > GlobalCoordinate
type of global coordinates
Definition: multilineargeometry.hh:170
ctype det() const
Definition: multilineargeometry.hh:415
TopologyId topologyId(integral_constant< bool, true >) const
Definition: multilineargeometry.hh:382
Base::ReferenceElement ReferenceElement
Definition: multilineargeometry.hh:447
Base::LocalCoordinate LocalCoordinate
Definition: multilineargeometry.hh:454
LocalCoordinate local(const GlobalCoordinate &global) const
evaluate the inverse mapping
Definition: multilineargeometry.hh:279
Definition: multilineargeometry.hh:399
ctype integrationElement(const LocalCoordinate &local) const
obtain the integration element
Definition: multilineargeometry.hh:542
bool checkInside(unsigned int topologyId, int dim, const FieldVector< ct, cdim > &x, ct tolerance, ct factor=ct(1))
Definition: referencedomain.hh:284
Base::GlobalCoordinate GlobalCoordinate
Definition: multilineargeometry.hh:455
Base::JacobianInverseTransposed JacobianInverseTransposed
Definition: multilineargeometry.hh:458
GlobalCoordinate global(const LocalCoordinate &local) const
evaluate the mapping
Definition: multilineargeometry.hh:260
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:24
JacobianInverseTransposed Jacobian
For backward-compatibility, export the type JacobianInverseTransposed as Jacobian.
Definition: multilineargeometry.hh:176
ctype integrationElement(const LocalCoordinate &local) const
obtain the integration element
Definition: multilineargeometry.hh:309
ctype detInv() const
Definition: multilineargeometry.hh:416
const JacobianInverseTransposed & jacobianInverseTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian's inverse
Definition: multilineargeometry.hh:589
const GeometryType & type(int i, int c) const
obtain the type of subentity (i,c)
Definition: referenceelements.hh:136
template specifying the storage for the corners
Definition: multilineargeometry.hh:95
conditional< hasSingleGeometryType, integral_constant< unsigned int, Traits::template hasSingleGeometryType< mydimension >::topologyId >, unsigned int >::type TopologyId
Definition: multilineargeometry.hh:191
const JacobianTransposed & jacobianTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian
Definition: multilineargeometry.hh:575
generic geometry implementation based on corner coordinates
Definition: multilineargeometry.hh:150
bool affine() const
is this mapping affine?
Definition: multilineargeometry.hh:232
static const int mydimension
geometry dimension
Definition: multilineargeometry.hh:159
ctype volume() const
obtain the volume of the mapping's image
Definition: multilineargeometry.hh:322
std::vector< FieldVector< ct, cdim > > Type
Definition: multilineargeometry.hh:97
JacobianTransposed jacobianTransposed_
Definition: multilineargeometry.hh:385