affinegeometry.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_GEOMETRY_AFFINEGEOMETRY_HH
#define DUNE_GEOMETRY_AFFINEGEOMETRY_HH
 
#include <cmath>
 
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
 
#include <dune/geometry/type.hh>
 
namespace Dune
{
 
  // External Forward Declarations
  // -----------------------------
 
  template< class ctype, int dim >
  class ReferenceElement;
 
  template< class ctype, int dim >
  struct ReferenceElements;
 
 
 
  namespace Impl
  {
 
    // FieldMatrixHelper
    // -----------------
 
    template< class ct >
    struct FieldMatrixHelper
    {
      typedef ct ctype;
 
      template< int m, int n >
      static void Ax ( const FieldMatrix< ctype, m, n > &A, const FieldVector< ctype, n > &x, FieldVector< ctype, m > &ret )
      {
        for( int i = 0; i < m; ++i )
        {
          ret[ i ] = ctype( 0 );
          for( int j = 0; j < n; ++j )
            ret[ i ] += A[ i ][ j ] * x[ j ];
        }
      }
 
      template< int m, int n >
      static void ATx ( const FieldMatrix< ctype, m, n > &A, const FieldVector< ctype, m > &x, FieldVector< ctype, n > &ret )
      {
        for( int i = 0; i < n; ++i )
        {
          ret[ i ] = ctype( 0 );
          for( int j = 0; j < m; ++j )
            ret[ i ] += A[ j ][ i ] * x[ j ];
        }
      }
 
      template< int m, int n, int p >
      static void AB ( const FieldMatrix< ctype, m, n > &A, const FieldMatrix< ctype, n, p > &B, FieldMatrix< ctype, m, p > &ret )
      {
        for( int i = 0; i < m; ++i )
        {
          for( int j = 0; j < p; ++j )
          {
            ret[ i ][ j ] = ctype( 0 );
            for( int k = 0; k < n; ++k )
              ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
          }
        }
      }
 
      template< int m, int n, int p >
      static void ATBT ( const FieldMatrix< ctype, m, n > &A, const FieldMatrix< ctype, p, m > &B, FieldMatrix< ctype, n, p > &ret )
      {
        for( int i = 0; i < n; ++i )
        {
          for( int j = 0; j < p; ++j )
          {
            ret[ i ][ j ] = ctype( 0 );
            for( int k = 0; k < m; ++k )
              ret[ i ][ j ] += A[ k ][ i ] * B[ j ][ k ];
          }
        }
      }
 
      template< int m, int n >
      static void ATA_L ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, n > &ret )
      {
        for( int i = 0; i < n; ++i )
        {
          for( int j = 0; j <= i; ++j )
          {
            ret[ i ][ j ] = ctype( 0 );
            for( int k = 0; k < m; ++k )
              ret[ i ][ j ] += A[ k ][ i ] * A[ k ][ j ];
          }
        }
      }
 
      template< int m, int n >
      static void ATA ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, n > &ret )
      {
        for( int i = 0; i < n; ++i )
        {
          for( int j = 0; j <= i; ++j )
          {
            ret[ i ][ j ] = ctype( 0 );
            for( int k = 0; k < m; ++k )
              ret[ i ][ j ] += A[ k ][ i ] * A[ k ][ j ];
            ret[ j ][ i ] = ret[ i ][ j ];
          }
 
          ret[ i ][ i ] = ctype( 0 );
          for( int k = 0; k < m; ++k )
            ret[ i ][ i ] += A[ k ][ i ] * A[ k ][ i ];
        }
      }
 
      template< int m, int n >
      static void AAT_L ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, m, m > &ret )
      {
        /*
           if (m==2) {
           ret[0][0] = A[0]*A[0];
           ret[1][1] = A[1]*A[1];
           ret[1][0] = A[0]*A[1];
           }
           else
         */
        for( int i = 0; i < m; ++i )
        {
          for( int j = 0; j <= i; ++j )
          {
            ctype &retij = ret[ i ][ j ];
            retij = A[ i ][ 0 ] * A[ j ][ 0 ];
            for( int k = 1; k < n; ++k )
              retij += A[ i ][ k ] * A[ j ][ k ];
          }
        }
      }
 
      template< int m, int n >
      static void AAT ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, m, m > &ret )
      {
        for( int i = 0; i < m; ++i )
        {
          for( int j = 0; j < i; ++j )
          {
            ret[ i ][ j ] = ctype( 0 );
            for( int k = 0; k < n; ++k )
              ret[ i ][ j ] += A[ i ][ k ] * A[ j ][ k ];
            ret[ j ][ i ] = ret[ i ][ j ];
          }
          ret[ i ][ i ] = ctype( 0 );
          for( int k = 0; k < n; ++k )
            ret[ i ][ i ] += A[ i ][ k ] * A[ i ][ k ];
        }
      }
 
      template< int n >
      static void Lx ( const FieldMatrix< ctype, n, n > &L, const FieldVector< ctype, n > &x, FieldVector< ctype, n > &ret )
      {
        for( int i = 0; i < n; ++i )
        {
          ret[ i ] = ctype( 0 );
          for( int j = 0; j <= i; ++j )
            ret[ i ] += L[ i ][ j ] * x[ j ];
        }
      }
 
      template< int n >
      static void LTx ( const FieldMatrix< ctype, n, n > &L, const FieldVector< ctype, n > &x, FieldVector< ctype, n > &ret )
      {
        for( int i = 0; i < n; ++i )
        {
          ret[ i ] = ctype( 0 );
          for( int j = i; j < n; ++j )
            ret[ i ] += L[ j ][ i ] * x[ j ];
        }
      }
 
      template< int n >
      static void LTL ( const FieldMatrix< ctype, n, n > &L, FieldMatrix< ctype, n, n > &ret )
      {
        for( int i = 0; i < n; ++i )
        {
          for( int j = 0; j < i; ++j )
          {
            ret[ i ][ j ] = ctype( 0 );
            for( int k = i; k < n; ++k )
              ret[ i ][ j ] += L[ k ][ i ] * L[ k ][ j ];
            ret[ j ][ i ] = ret[ i ][ j ];
          }
          ret[ i ][ i ] = ctype( 0 );
          for( int k = i; k < n; ++k )
            ret[ i ][ i ] += L[ k ][ i ] * L[ k ][ i ];
        }
      }
 
      template< int n >
      static void LLT ( const FieldMatrix< ctype, n, n > &L, FieldMatrix< ctype, n, n > &ret )
      {
        for( int i = 0; i < n; ++i )
        {
          for( int j = 0; j < i; ++j )
          {
            ret[ i ][ j ] = ctype( 0 );
            for( int k = 0; k <= j; ++k )
              ret[ i ][ j ] += L[ i ][ k ] * L[ j ][ k ];
            ret[ j ][ i ] = ret[ i ][ j ];
          }
          ret[ i ][ i ] = ctype( 0 );
          for( int k = 0; k <= i; ++k )
            ret[ i ][ i ] += L[ i ][ k ] * L[ i ][ k ];
        }
      }
 
      template< int n >
      static void cholesky_L ( const FieldMatrix< ctype, n, n > &A, FieldMatrix< ctype, n, n > &ret )
      {
        for( int i = 0; i < n; ++i )
        {
          ctype &rii = ret[ i ][ i ];
 
          ctype xDiag = A[ i ][ i ];
          for( int j = 0; j < i; ++j )
            xDiag -= ret[ i ][ j ] * ret[ i ][ j ];
          assert( xDiag > ctype( 0 ) );
          rii = sqrt( xDiag );
 
          ctype invrii = ctype( 1 ) / rii;
          for( int k = i+1; k < n; ++k )
          {
            ctype x = A[ k ][ i ];
            for( int j = 0; j < i; ++j )
              x -= ret[ i ][ j ] * ret[ k ][ j ];
            ret[ k ][ i ] = invrii * x;
          }
        }
      }
 
      template< int n >
      static ctype detL ( const FieldMatrix< ctype, n, n > &L )
      {
        ctype det( 1 );
        for( int i = 0; i < n; ++i )
          det *= L[ i ][ i ];
        return det;
      }
 
      template< int n >
      static ctype invL ( FieldMatrix< ctype, n, n > &L )
      {
        ctype det( 1 );
        for( int i = 0; i < n; ++i )
        {
          ctype &lii = L[ i ][ i ];
          det *= lii;
          lii = ctype( 1 ) / lii;
          for( int j = 0; j < i; ++j )
          {
            ctype &lij = L[ i ][ j ];
            ctype x = lij * L[ j ][ j ];
            for( int k = j+1; k < i; ++k )
              x += L[ i ][ k ] * L[ k ][ j ];
            lij = (-lii) * x;
          }
        }
        return det;
      }
 
      // calculates x := L^{-1} x
      template< int n >
      static void invLx ( FieldMatrix< ctype, n, n > &L, FieldVector< ctype, n > &x )
      {
        for( int i = 0; i < n; ++i )
        {
          for( int j = 0; j < i; ++j )
            x[ i ] -= L[ i ][ j ] * x[ j ];
          x[ i ] /= L[ i ][ i ];
        }
      }
 
      // calculates x := L^{-T} x
      template< int n >
      static void invLTx ( FieldMatrix< ctype, n, n > &L, FieldVector< ctype, n > &x )
      {
        for( int i = n; i > 0; --i )
        {
          for( int j = i; j < n; ++j )
            x[ i-1 ] -= L[ j ][ i-1 ] * x[ j ];
          x[ i-1 ] /= L[ i-1 ][ i-1 ];
        }
      }
 
      template< int n >
      static ctype spdDetA ( const FieldMatrix< ctype, n, n > &A )
      {
        // return A[0][0]*A[1][1]-A[1][0]*A[1][0];
        FieldMatrix< ctype, n, n > L;
        cholesky_L( A, L );
        return detL( L );
      }
 
      template< int n >
      static ctype spdInvA ( FieldMatrix< ctype, n, n > &A )
      {
        FieldMatrix< ctype, n, n > L;
        cholesky_L( A, L );
        const ctype det = invL( L );
        LTL( L, A );
        return det;
      }
 
      // calculate x := A^{-1} x
      template< int n >
      static void spdInvAx ( FieldMatrix< ctype, n, n > &A, FieldVector< ctype, n > &x )
      {
        FieldMatrix< ctype, n, n > L;
        cholesky_L( A, L );
        invLx( L, x );
        invLTx( L, x );
      }
 
      template< int m, int n >
      static ctype detATA ( const FieldMatrix< ctype, m, n > &A )
      {
        if( m >= n )
        {
          FieldMatrix< ctype, n, n > ata;
          ATA_L( A, ata );
          return spdDetA( ata );
        }
        else
          return ctype( 0 );
      }
 
      template< int m, int n >
      static ctype sqrtDetAAT ( const FieldMatrix< ctype, m, n > &A )
      {
        using std::abs;
        using std::sqrt;
        // These special cases are here not only for speed reasons:
        // The general implementation aborts if the matrix is almost singular,
        // and the special implementation provide a stable way to handle that case.
        if( (n == 2) && (m == 2) )
        {
          // Special implementation for 2x2 matrices: faster and more stable
          return abs( A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ] );
        }
        else if( (n == 3) && (m == 3) )
        {
          // Special implementation for 3x3 matrices
          const ctype v0 = A[ 0 ][ 1 ] * A[ 1 ][ 2 ] - A[ 1 ][ 1 ] * A[ 0 ][ 2 ];
          const ctype v1 = A[ 0 ][ 2 ] * A[ 1 ][ 0 ] - A[ 1 ][ 2 ] * A[ 0 ][ 0 ];
          const ctype v2 = A[ 0 ][ 0 ] * A[ 1 ][ 1 ] - A[ 1 ][ 0 ] * A[ 0 ][ 1 ];
          return abs( v0 * A[ 2 ][ 0 ] + v1 * A[ 2 ][ 1 ] + v2 * A[ 2 ][ 2 ] );
        }
        else if ( (n == 3) && (m == 2) )
        {
          // Special implementation for 2x3 matrices
          const ctype v0 = A[ 0 ][ 0 ] * A[ 1 ][ 1 ] - A[ 0 ][ 1 ] * A[ 1 ][ 0 ];
          const ctype v1 = A[ 0 ][ 0 ] * A[ 1 ][ 2 ] - A[ 1 ][ 0 ] * A[ 0 ][ 2 ];
          const ctype v2 = A[ 0 ][ 1 ] * A[ 1 ][ 2 ] - A[ 0 ][ 2 ] * A[ 1 ][ 1 ];
          return sqrt( v0*v0 + v1*v1 + v2*v2);
        }
        else if( n >= m )
        {
          // General case
          FieldMatrix< ctype, m, m > aat;
          AAT_L( A, aat );
          return spdDetA( aat );
        }
        else
          return ctype( 0 );
      }
 
      // A^{-1}_L = (A^T A)^{-1} A^T
      // => A^{-1}_L A = I
      template< int m, int n >
      static ctype leftInvA ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, m > &ret )
      {
        static_assert((m >= n), "Matrix has no left inverse.");
        FieldMatrix< ctype, n, n > ata;
        ATA_L( A, ata );
        const ctype det = spdInvA( ata );
        ATBT( ata, A, ret );
        return det;
      }
 
      template< int m, int n >
      static void leftInvAx ( const FieldMatrix< ctype, m, n > &A, const FieldVector< ctype, m > &x, FieldVector< ctype, n > &y )
      {
        static_assert((m >= n), "Matrix has no left inverse.");
        FieldMatrix< ctype, n, n > ata;
        ATx( A, x, y );
        ATA_L( A, ata );
        spdInvAx( ata, y );
      }
 
      template< int m, int n >
      static ctype rightInvA ( const FieldMatrix< ctype, m, n > &A, FieldMatrix< ctype, n, m > &ret )
      {
        static_assert((n >= m), "Matrix has no right inverse.");
        using std::abs;
        if( (n == 2) && (m == 2) )
        {
          const ctype det = (A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ]);
          const ctype detInv = ctype( 1 ) / det;
          ret[ 0 ][ 0 ] = A[ 1 ][ 1 ] * detInv;
          ret[ 1 ][ 1 ] = A[ 0 ][ 0 ] * detInv;
          ret[ 1 ][ 0 ] = -A[ 1 ][ 0 ] * detInv;
          ret[ 0 ][ 1 ] = -A[ 0 ][ 1 ] * detInv;
          return abs( det );
        }
        else
        {
          FieldMatrix< ctype, m , m > aat;
          AAT_L( A, aat );
          const ctype det = spdInvA( aat );
          ATBT( A , aat , ret );
          return det;
        }
      }
 
      template< int m, int n >
      static void xTRightInvA ( const FieldMatrix< ctype, m, n > &A, const FieldVector< ctype, n > &x, FieldVector< ctype, m > &y )
      {
        static_assert((n >= m), "Matrix has no right inverse.");
        FieldMatrix< ctype, m, m > aat;
        Ax( A, x, y );
        AAT_L( A, aat );
        spdInvAx( aat, y );
      }
    };
 
  } // namespace Impl
 
 
 
  template< class ct, int mydim, int cdim>
  class AffineGeometry
  {
  public:
 
    typedef ct ctype;
 
    static const int mydimension= mydim;
 
    static const int coorddimension = cdim;
 
    typedef FieldVector< ctype, mydimension > LocalCoordinate;
 
    typedef FieldVector< ctype, coorddimension > GlobalCoordinate;
 
    typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
 
    typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;
 
  private:
    typedef Dune::ReferenceElement< ctype, mydimension > ReferenceElement;
 
    typedef Dune::ReferenceElements< ctype, mydimension > ReferenceElements;
 
    // Helper class to compute a matrix pseudo inverse
    typedef Impl::FieldMatrixHelper< ct > MatrixHelper;
 
  public:
    AffineGeometry ( const ReferenceElement &refElement, const GlobalCoordinate &origin,
                     const JacobianTransposed &jt )
      : refElement_(&refElement), origin_(origin), jacobianTransposed_(jt)
    {
      integrationElement_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jacobianTransposed_, jacobianInverseTransposed_ );
    }
 
    AffineGeometry ( Dune::GeometryType gt, const GlobalCoordinate &origin,
                     const JacobianTransposed &jt )
      : refElement_( &ReferenceElements::general( gt ) ), origin_(origin), jacobianTransposed_( jt )
    {
      integrationElement_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jacobianTransposed_, jacobianInverseTransposed_ );
    }
 
    template< class CoordVector >
    AffineGeometry ( const ReferenceElement &refElement, const CoordVector &coordVector )
      : refElement_(&refElement), origin_(coordVector[0])
    {
      for( int i = 0; i < mydimension; ++i )
        jacobianTransposed_[ i ] = coordVector[ i+1 ] - origin_;
      integrationElement_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jacobianTransposed_, jacobianInverseTransposed_ );
    }
 
    template< class CoordVector >
    AffineGeometry ( Dune::GeometryType gt, const CoordVector &coordVector )
      : refElement_(&ReferenceElements::general( gt )), origin_(coordVector[0] )
    {
      for( int i = 0; i < mydimension; ++i )
        jacobianTransposed_[ i ] = coordVector[ i+1 ] - origin_;
      integrationElement_ = MatrixHelper::template rightInvA< mydimension, coorddimension >( jacobianTransposed_, jacobianInverseTransposed_ );
    }
 
    bool affine () const { return true; }
 
    Dune::GeometryType type () const { return refElement_->type(); }
 
    int corners () const { return refElement_->size( mydimension ); }
 
    GlobalCoordinate corner ( int i ) const
    {
      return global( refElement_->position( i, mydimension ) );
    }
 
    GlobalCoordinate center () const { return global( refElement_->position( 0, 0 ) ); }
 
    GlobalCoordinate global ( const LocalCoordinate &local ) const
    {
      GlobalCoordinate global( origin_ );
      jacobianTransposed_.umtv( local, global );
      return global;
    }
 
    LocalCoordinate local ( const GlobalCoordinate &global ) const
    {
      LocalCoordinate local;
      jacobianInverseTransposed_.mtv( global - origin_, local );
      return local;
    }
 
    ctype integrationElement ( const LocalCoordinate &local ) const
    {
      DUNE_UNUSED_PARAMETER(local);
      return integrationElement_;
    }
 
    ctype volume () const
    {
      return integrationElement_ * refElement_->volume();
    }
 
    const JacobianTransposed &jacobianTransposed ( const LocalCoordinate &local ) const
    {
      DUNE_UNUSED_PARAMETER(local);
      return jacobianTransposed_;
    }
 
    const JacobianInverseTransposed &jacobianInverseTransposed ( const LocalCoordinate &local ) const
    {
      DUNE_UNUSED_PARAMETER(local);
      return jacobianInverseTransposed_;
    }
 
    friend const ReferenceElement &referenceElement ( const AffineGeometry &geometry ) { return *geometry.refElement_; }
 
  private:
    const ReferenceElement* refElement_;
    GlobalCoordinate origin_;
    JacobianTransposed jacobianTransposed_;
    JacobianInverseTransposed jacobianInverseTransposed_;
    ctype integrationElement_;
  };
 
} // namespace Dune
 
#endif // #ifndef DUNE_GEOMETRY_AFFINEGEOMETRY_HH