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();
282 const ctype tolerance = Traits::tolerance();
289 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 )
519 if( jacobianInverseTransposedComputed_ )
547 if( !integrationElementComputed_ )
550 integrationElementComputed_ =
true;
594 if( !jacobianInverseTransposedComputed_ )
597 jacobianInverseTransposedComputed_ =
true;
598 integrationElementComputed_ =
true;
613 mutable bool affine_ : 1;
615 mutable bool jacobianInverseTransposedComputed_ : 1;
616 mutable bool integrationElementComputed_ : 1;
624 template<
class ct,
int mydim,
int cdim,
class Traits >
625 inline const typename MultiLinearGeometry< ct, mydim, cdim, Traits >::JacobianInverseTransposed &
629 jacobianInverseTransposed_.
setup( jacobianTransposed( local ) );
630 return jacobianInverseTransposed_;
634 template<
class ct,
int mydim,
int cdim,
class Traits >
635 template<
bool add,
int dim >
641 const ctype xn = df*x[ dim-1 ];
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() || 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 ];
691 CornerIterator cit2( cit );
695 jacobianTransposed< add >( topologyId, integral_constant< int, dim-1 >(), cit2, df, x, rf*cxn, jt );
697 jacobianTransposed< true >( topologyId, integral_constant< int, dim-1 >(), cit2, df, x, rf*xn, jt );
699 global< add >( topologyId, integral_constant< int, dim-1 >(), cit, df, x, -rf, jt[ dim-1 ] );
700 global< true >( topologyId, integral_constant< int, dim-1 >(), cit, df, x, rf, jt[ dim-1 ] );
746 ctype dfcxn = (cxn > Traits::tolerance() || cxn < -Traits::tolerance()) ? df / cxn : 0;
751 global< add >( topologyId, integral_constant< int, dim-1 >(), cit, dfcxn, x, -rf, jt[ dim-1 ] );
753 jt[ dim-1 ].axpy( rf, *cit );
758 FieldMatrix<
ctype, dim-1, coorddimension > jt2;
760 jacobianTransposed< false >( topologyId, integral_constant< int, dim-1 >(), cit2, dfcxn, x, rf, jt2 );
764 for(
int j = 0; j < dim-1; ++j )
767 jt[ dim-1 ].axpy( dfcxn*x[ j ], jt2[ j ] );
773 jacobianTransposed< false >( topologyId, integral_constant< int, dim-1 >(), cit2, dfcxn, x, rf, jt );
775 for(
int j = 0; j < dim-1; ++j )
776 jt[ dim-1 ].axpy( dfcxn*x[ j ], jt[ j ] );
781 template<
class ct,
int mydim,
int cdim,
class Traits >
782 template<
bool add,
int rows >
786 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt )
793 template<
class ct,
int mydim,
int cdim,
class Traits >
799 if( !affine( topologyId, integral_constant< int, dim-1 >(), cit, jt ) )
806 if( !affine( topologyId, integral_constant< int, dim-1 >(), cit, jtTop ) )
811 for(
int i = 0; i < dim-1; ++i )
812 norm += (jtTop[ i ] - jt[ i ]).two_norm2();
813 if( norm >= Traits::tolerance() )
818 jt[ dim-1 ] = orgTop - orgBottom;
822 template<
class ct,
int mydim,
int cdim,
class Traits >
832 #endif // #ifndef DUNE_GEOMETRY_MULTILINEARGEOMETRY_HH