/usr/include/dune/geometry/test/checkgeometry.hh is in libdune-geometry-dev 2.5.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 | // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_CHECK_GEOMETRY_HH
#define DUNE_CHECK_GEOMETRY_HH
#include <limits>
#include <dune/common/typetraits.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/referenceelements.hh>
namespace Dune
{
/**
* \brief Static and dynamic checks for all features of a Geometry
*
* This excludes anything related to being part of a grid.
*
* \tparam TestGeometry The type of the geometry to be tested
*
* \param geometry The TestGeometry object to be tested
*
* \returns true if check passed
*/
template <class TestGeometry>
bool checkGeometry ( const TestGeometry& geometry )
{
bool pass = true;
////////////////////////////////////////////////////////////////
// Extract all static information
////////////////////////////////////////////////////////////////
// dimension of the corresponding reference element
static const int mydim = TestGeometry::mydimension;
// dimension of the world space
static const int coorddim = TestGeometry::coorddimension;
// type used for coordinate coefficients
typedef typename TestGeometry::ctype ctype;
// vector type used for points in the domain
typedef typename TestGeometry::LocalCoordinate LocalCoordinate DUNE_UNUSED;
// vector type used for image points
typedef typename TestGeometry::GlobalCoordinate GlobalCoordinate DUNE_UNUSED;
// Matrix-like type for the return value of the jacobianTransposed method
typedef typename TestGeometry::JacobianTransposed JacobianTransposed;
// Matrix-like type for the return value of the jacobianInverseTransposed method
typedef typename TestGeometry::JacobianInverseTransposed JacobianInverseTransposed;
const ctype tolerance = std::sqrt( std::numeric_limits< ctype >::epsilon() );
////////////////////////////////////////////////////////////////////////
const int corners = geometry.corners();
if( corners == 0 )
DUNE_THROW( Exception, "A geometry must have at least one corner." );
GlobalCoordinate cornerAvg ( 0 );
for( int i = 0; i < corners; ++i )
cornerAvg += geometry.corner( i );
cornerAvg /= ctype( corners );
const GlobalCoordinate center = geometry.center();
if( corners > 1 )
{
// if we have more than one corner, no corner may coincide with their average or the center
for( int i = 0; i < corners; ++i )
{
const GlobalCoordinate corner = geometry.corner( i );
if( (corner - center).two_norm() <= tolerance )
{
std::cerr << "Error: Geometry has multiple corners, but one corner coincides with the center." << std::endl;
pass = false;
}
if( (corner - cornerAvg).two_norm() <= tolerance )
{
std::cerr << "Error: Geometry has multiple corners, but one corner coincides with their average." << std::endl;
pass = false;
}
}
}
else
{
// the single corner must coincide with the center
if( (center - cornerAvg).two_norm() > tolerance )
{
std::cerr << "Error: Geometry has a single corner (" << cornerAvg << "), but it does not coincide with the center (" << center << ")." << std::endl;
pass = false;
}
}
const ctype volume = geometry.volume();
if( volume < tolerance )
{
std::cerr << "Error: Geometry has nearly vanishing volume (" << volume << ")" << std::endl;
pass = false;
}
////////////////////////////////////////////////////////////////////////
const GeometryType type = geometry.type();
if( type.isNone() )
return pass;
const ReferenceElement< ctype, mydim > &refElement = ReferenceElements< ctype, mydim >::general(type);
// Test whether the return value of the method 'center' corresponds to the center of the
// reference element. That is the current definition of the method.
if( (center - geometry.global( refElement.position( 0, 0 ) )).two_norm() > tolerance )
DUNE_THROW( Exception, "center() is not consistent with global(refElem.position(0,0))." );
////////////////////////////////////////////////////////////////////////
// Test whether the number and placement of the corners is consistent
// with the corners of the corresponding reference element.
////////////////////////////////////////////////////////////////////////
if( refElement.size( mydim ) == corners )
{
for( int i = 0; i < geometry.corners(); ++i )
{
if( (geometry.corner( i ) - geometry.global( refElement.position( i, mydim ) )).two_norm() > tolerance )
{
std::cerr << "Error: Methods corner and global are inconsistent." << std::endl;
pass = false;
}
}
}
else
{
std::cerr << "Error: Incorrect number of corners (" << geometry.corners() << ", should be " << refElement.size( mydim ) << ")." << std::endl;
pass = false;
}
///////////////////////////////////////////////////////////////////////////////
// Use a quadrature rule as a set of test points and loop over them
///////////////////////////////////////////////////////////////////////////////
const Dune::QuadratureRule<ctype, mydim> & quadrature = Dune::QuadratureRules<ctype, mydim>::rule(geometry.type(), 2);
for (const auto& ip : quadrature)
{
const typename TestGeometry::LocalCoordinate &x = ip.position();
// Test whether the methods 'local' and 'global' are inverse to each other
if ( (x - geometry.local( geometry.global( x ) )).two_norm() > 1e-8 ) {
std::cerr << "Error: global and local are not inverse to each other." << std::endl;
pass = false;
}
// Test whether the methods 'jacobianTransposed' and 'jacobianInverseTransposed'
// return matrices that are inverse to each other.
const JacobianTransposed &jt = geometry.jacobianTransposed( x );
const JacobianInverseTransposed &jit = geometry.jacobianInverseTransposed( x );
// Transform to FieldMatrix, so we can have coefficent access and other goodies
// We need some black magic for the transformation, because there is no
// official easy way yet.
// The following code does the transformation by multiplying jt and jit from
// the right by identity matrices. That way, only the mv method is used.
FieldMatrix< ctype, mydim, coorddim > jtAsFieldMatrix;
for (int j=0; j<coorddim; j++) {
FieldVector<ctype,coorddim> idColumn(0);
idColumn[j] = 1;
FieldVector<ctype,mydim> column;
jt.mv(idColumn,column);
for (int k=0; k<mydim; k++)
jtAsFieldMatrix[k][j] = column[k];
}
FieldMatrix< ctype, coorddim, mydim > jitAsFieldMatrix;
for (int j=0; j<mydim; j++) {
FieldVector<ctype,mydim> idColumn(0);
idColumn[j] = 1;
FieldVector<ctype,coorddim> column;
jit.mv(idColumn,column);
for (int k=0; k<coorddim; k++)
jitAsFieldMatrix[k][j] = column[k];
}
FieldMatrix< ctype, mydim, mydim > id;
FMatrixHelp::multMatrix( jtAsFieldMatrix, jitAsFieldMatrix, id );
bool isId = true;
for( int j = 0; j < mydim; ++j )
for( int k = 0; k < mydim; ++k )
isId &= (std::abs( id[ j ][ k ] - (j == k ? 1 : 0) ) < 1e-8);
if( !isId)
{
std::cerr << "Error: jacobianTransposed and jacobianInverseTransposed are not inverse to each other." << std::endl;
std::cout << " id != [ ";
for( int j = 0; j < mydim; ++j )
std::cout << (j > 0 ? " | " : "") << id[ j ];
std::cout << " ]" << std::endl;
pass = false;
}
// Test whether integrationElement returns something nonnegative
if( geometry.integrationElement( x ) < 0 ) {
std::cerr << "Error: Negative integrationElement found." << std::endl;
pass = false;
}
FieldMatrix< ctype, mydim, mydim > jtj( 0 );
for( int i = 0; i < mydim; ++i )
for( int j = 0; j < mydim; ++j )
for( int k = 0; k < coorddim; ++k )
jtj[ i ][ j ] += jtAsFieldMatrix[ i ][ k ] * jtAsFieldMatrix[ j ][ k ];
if( std::abs( std::sqrt( jtj.determinant() ) - geometry.integrationElement( x ) ) > 1e-8 ) {
std::cerr << "Error: integrationElement is not consistent with jacobianTransposed." << std::endl;
pass = false;
}
if (geometry.affine())
if( std::abs( geometry.volume() - refElement.volume()*geometry.integrationElement( x ) ) > 1e-8 ) {
std::cerr << "Error: volume is not consistent with jacobianTransposed." << std::endl;
pass = false;
}
}
return pass;
}
}
#endif // #ifndef DUNE_CHECK_GEOMETRY_HH
|