Commit 76b1d4bf authored by Christian Engwer's avatar Christian Engwer
Browse files

Merge branch 'feature/remove-generic-geometry' into 'master'

remove generic geometry subdirectory

On the last developer meeting, we agree on the following things:
- The `genericgeometry` subdirectory should vanish in the next DUNE release.
- All content of the namespace GenericGeometry is considered an implemenation detail.

This merge implements these decisions. All functionality that is sitll required is migrated into the appropriate headers in `dune/geometry` and put into the namespace `Impl`. 

**Note**: There is *no transition phase* for the content of the `genericgeometry` subdirectory.

See merge request !34
parents 448af15f 2a0f48db
Pipeline #1102 passed with stage
in 12 minutes and 51 seconds
add_subdirectory("quadraturerules")
add_subdirectory("genericgeometry")
add_subdirectory("refinement")
add_subdirectory("test")
......@@ -23,3 +22,7 @@ install(FILES
# install the header as done for the auto-tools
install(FILES test/checkgeometry.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/geometry/test)
dune_add_library(geometry OBJECT
referenceelements.cc
)
......@@ -8,12 +8,12 @@
* \author Martin Nolte
*/
#include <cmath>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/geometry/type.hh>
#include <dune/geometry/genericgeometry/geometrytraits.hh>
#include <dune/geometry/genericgeometry/matrixhelper.hh>
namespace Dune
{
......@@ -29,6 +29,428 @@ namespace Dune
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 );
}
/** \brief Compute the square root of the determinant of A times A transposed
*
* This is the volume element for an embedded submanifold and needed to
* implement the method integrationElement().
*/
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 );
}
/** \brief Compute right pseudo-inverse of matrix A */
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
/** \brief Implementation of the Geometry interface for affine geometries
* \tparam ct Type used for coordinates
* \tparam mydim Dimension of the geometry
......@@ -67,7 +489,7 @@ namespace Dune
typedef Dune::ReferenceElements< ctype, mydimension > ReferenceElements;
// Helper class to compute a matrix pseudo inverse
typedef GenericGeometry::MatrixHelper< GenericGeometry::DuneCoordTraits< ct > > MatrixHelper;
typedef Impl::FieldMatrixHelper< ct > MatrixHelper;
public:
/** \brief Create affine geometry from reference element, one vertex, and the Jacobian matrix */
......
add_subdirectory("test")
install(FILES
codimtable.hh
conversion.hh
geometrytraits.hh
matrixhelper.hh
referencedomain.hh
referenceelements.hh
subtopologies.hh
topologytypes.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/geometry/genericgeometry)
dune_add_library(genericgeometry OBJECT
referencedomain.cc
subtopologies.cc
)
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_GEOMETRY_GENERICGEOMETRY_CODIMTABLE_HH
#define DUNE_GEOMETRY_GENERICGEOMETRY_CODIMTABLE_HH
#include <dune/common/typetraits.hh>
#include <dune/common/tupleutility.hh>
#include <dune/common/unused.hh>
namespace Dune
{
namespace GenericGeometry
{
template< template< int > class Element, int dim >
class CodimTable
{
friend class CodimTable< Element, dim+1 >;
typedef typename PushBackTuple<
typename CodimTable< Element, dim-1 >::ElementTuple,
Element< dim > >::type ElementTuple;
ElementTuple map_;
public:
template< int codim >
const Element< codim > &
operator[] ( const std::integral_constant< int, codim > codimVariable ) const
{
DUNE_UNUSED_PARAMETER(codimVariable);
return std::get<codim>(map_);
}
template< int codim >
Element< codim > &
operator[] ( const std::integral_constant< int, codim > codimVariable )
{
DUNE_UNUSED_PARAMETER(codimVariable);
return std::get<codim>(map_);
}
};
template< template< int > class Element>
class CodimTable< Element, -1 >
{
friend class CodimTable< Element, 0 >;
typedef typename std::tuple<> ElementTuple;
};
}
}
#endif // #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_CODIMTABLE_HH
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_GEOMETRY_GENERICGEOMETRY_CONVERSION_HH
#define DUNE_GEOMETRY_GENERICGEOMETRY_CONVERSION_HH
#include <dune/common/visibility.hh>
#include <dune/geometry/type.hh>
#include <dune/geometry/genericgeometry/topologytypes.hh>
#include <dune/geometry/genericgeometry/subtopologies.hh>
namespace Dune
{
namespace GenericGeometry
{
// DuneGeometryType
// ----------------
/** \class DuneGeometryType
* \ingroup GenericGeometry
* \brief statically convert a generic topology type into a GeometryType
*
* \tparam Topology topology type to be converted
* \tparam linetype basic geometry type to assign to a line
* (either simplex or cube)
*/
template< class Topology, GeometryType::BasicType linetype >
class DuneGeometryType;
template< GeometryType::BasicType linetype >
class DuneGeometryType< Point, linetype >
{
static_assert((linetype == GeometryType::simplex)
|| (linetype == GeometryType::cube),
"Parameter linetype may only be a simplex or a cube.");
public:
static const unsigned int dimension = 0;
static const GeometryType::BasicType basicType = linetype;
};