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 > >
200 class JacobianInverseTransposed;
212 static const bool hasSingleGeometryType = Traits::template hasSingleGeometryType< mydimension >::v;
216 typedef typename std::conditional< hasSingleGeometryType, std::integral_constant< unsigned int, Traits::template hasSingleGeometryType< mydimension >::topologyId >,
unsigned int >
::type TopologyId;
228 template<
class Corners >
244 template<
class Corners >
267 assert( (i >= 0) && (i <
corners()) );
268 return std::cref(corners_).get()[ i ];
284 auto cit = begin(std::cref(corners_).get());
286 global< false >(
topologyId(), std::integral_constant< int, mydimension >(), cit,
ctype( 1 ), local,
ctype( 1 ), y );
304 const ctype tolerance = Traits::tolerance();
311 MatrixHelper::template xTRightInvA< mydimension, coorddimension >(
jacobianTransposed( x ), dglobal, dx );
313 }
while( dx.two_norm2() > tolerance );
333 return MatrixHelper::template sqrtDetAAT< mydimension, coorddimension >(
jacobianTransposed( local ) );
363 auto cit = begin(std::cref(corners_).get());
364 jacobianTransposed< false >(
topologyId(), std::integral_constant< int, mydimension >(), cit,
ctype( 1 ), local,
ctype( 1 ), jt );
390 return topologyId( std::integral_constant< bool, hasSingleGeometryType >() );
393 template<
bool add,
int dim,
class CornerIterator >
397 template<
bool add,
class CornerIterator >
402 template<
bool add,
int rows,
int dim,
class CornerIterator >
405 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt );
406 template<
bool add,
int rows,
class CornerIterator >
409 const ctype &rf, FieldMatrix< ctype, rows, cdim > &jt );
411 template<
int dim,
class CornerIterator >
413 template<
class CornerIterator >
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; }
432 unsigned int topologyId ( std::integral_constant< bool, false > )
const {
return refElement().type().id(); }
435 typename Traits::template CornerStorage< mydimension, coorddimension >::Type corners_;
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 > >
505 template<
class 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 )
540 jacobianTransposed_.umtv( local,
global );
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_ )
595 integrationElementComputed_ =
true;
597 return jacobianInverseTransposed_.
detInv();
624 return jacobianTransposed_;
639 if( !jacobianInverseTransposedComputed_ )
641 jacobianInverseTransposed_.
setup( jacobianTransposed_ );
642 jacobianInverseTransposedComputed_ =
true;
643 integrationElementComputed_ =
true;
645 return 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 >
670 inline typename MultiLinearGeometry< ct, mydim, cdim, Traits >::JacobianInverseTransposed
673 JacobianInverseTransposed jit;
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