Dune Core Modules (2.5.0)

checkgeometry.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3
4#ifndef DUNE_CHECK_GEOMETRY_HH
5#define DUNE_CHECK_GEOMETRY_HH
6
7#include <limits>
8
12
13#include <dune/geometry/referenceelements.hh>
15#include <dune/geometry/referenceelements.hh>
16
17namespace Dune
18{
19
31 template <class TestGeometry>
32 bool checkGeometry ( const TestGeometry& geometry )
33 {
34 bool pass = true;
35
37 // Extract all static information
39
40 // dimension of the corresponding reference element
41 static const int mydim = TestGeometry::mydimension;
42
43 // dimension of the world space
44 static const int coorddim = TestGeometry::coorddimension;
45
46 // type used for coordinate coefficients
47 typedef typename TestGeometry::ctype ctype;
48
49 // vector type used for points in the domain
50 typedef typename TestGeometry::LocalCoordinate LocalCoordinate DUNE_UNUSED;
51
52 // vector type used for image points
53 typedef typename TestGeometry::GlobalCoordinate GlobalCoordinate DUNE_UNUSED;
54
55 // Matrix-like type for the return value of the jacobianTransposed method
56 typedef typename TestGeometry::JacobianTransposed JacobianTransposed;
57
58 // Matrix-like type for the return value of the jacobianInverseTransposed method
59 typedef typename TestGeometry::JacobianInverseTransposed JacobianInverseTransposed;
60
61 const ctype tolerance = std::sqrt( std::numeric_limits< ctype >::epsilon() );
62
64
65 const int corners = geometry.corners();
66 if( corners == 0 )
67 DUNE_THROW( Exception, "A geometry must have at least one corner." );
68
69 GlobalCoordinate cornerAvg ( 0 );
70 for( int i = 0; i < corners; ++i )
71 cornerAvg += geometry.corner( i );
72 cornerAvg /= ctype( corners );
73
74 const GlobalCoordinate center = geometry.center();
75
76 if( corners > 1 )
77 {
78 // if we have more than one corner, no corner may coincide with their average or the center
79 for( int i = 0; i < corners; ++i )
80 {
81 const GlobalCoordinate corner = geometry.corner( i );
82 if( (corner - center).two_norm() <= tolerance )
83 {
84 std::cerr << "Error: Geometry has multiple corners, but one corner coincides with the center." << std::endl;
85 pass = false;
86 }
87
88 if( (corner - cornerAvg).two_norm() <= tolerance )
89 {
90 std::cerr << "Error: Geometry has multiple corners, but one corner coincides with their average." << std::endl;
91 pass = false;
92 }
93 }
94 }
95 else
96 {
97 // the single corner must coincide with the center
98 if( (center - cornerAvg).two_norm() > tolerance )
99 {
100 std::cerr << "Error: Geometry has a single corner (" << cornerAvg << "), but it does not coincide with the center (" << center << ")." << std::endl;
101 pass = false;
102 }
103 }
104
105 const ctype volume = geometry.volume();
106 if( volume < tolerance )
107 {
108 std::cerr << "Error: Geometry has nearly vanishing volume (" << volume << ")" << std::endl;
109 pass = false;
110 }
111
113
114 const GeometryType type = geometry.type();
115 if( type.isNone() )
116 return pass;
117
119
120 // Test whether the return value of the method 'center' corresponds to the center of the
121 // reference element. That is the current definition of the method.
122 if( (center - geometry.global( refElement.position( 0, 0 ) )).two_norm() > tolerance )
123 DUNE_THROW( Exception, "center() is not consistent with global(refElem.position(0,0))." );
124
126 // Test whether the number and placement of the corners is consistent
127 // with the corners of the corresponding reference element.
129
130 if( refElement.size( mydim ) == corners )
131 {
132 for( int i = 0; i < geometry.corners(); ++i )
133 {
134 if( (geometry.corner( i ) - geometry.global( refElement.position( i, mydim ) )).two_norm() > tolerance )
135 {
136 std::cerr << "Error: Methods corner and global are inconsistent." << std::endl;
137 pass = false;
138 }
139 }
140 }
141 else
142 {
143 std::cerr << "Error: Incorrect number of corners (" << geometry.corners() << ", should be " << refElement.size( mydim ) << ")." << std::endl;
144 pass = false;
145 }
146
148 // Use a quadrature rule as a set of test points and loop over them
151 for (const auto& ip : quadrature)
152 {
153 const typename TestGeometry::LocalCoordinate &x = ip.position();
154
155 // Test whether the methods 'local' and 'global' are inverse to each other
156 if ( (x - geometry.local( geometry.global( x ) )).two_norm() > 1e-8 ) {
157 std::cerr << "Error: global and local are not inverse to each other." << std::endl;
158 pass = false;
159 }
160
161 // Test whether the methods 'jacobianTransposed' and 'jacobianInverseTransposed'
162 // return matrices that are inverse to each other.
163 const JacobianTransposed &jt = geometry.jacobianTransposed( x );
164 const JacobianInverseTransposed &jit = geometry.jacobianInverseTransposed( x );
165
166 // Transform to FieldMatrix, so we can have coefficent access and other goodies
167 // We need some black magic for the transformation, because there is no
168 // official easy way yet.
169 // The following code does the transformation by multiplying jt and jit from
170 // the right by identity matrices. That way, only the mv method is used.
172 for (int j=0; j<coorddim; j++) {
173
174 FieldVector<ctype,coorddim> idColumn(0);
175 idColumn[j] = 1;
176
178 jt.mv(idColumn,column);
179
180 for (int k=0; k<mydim; k++)
181 jtAsFieldMatrix[k][j] = column[k];
182
183 }
184
186 for (int j=0; j<mydim; j++) {
187
188 FieldVector<ctype,mydim> idColumn(0);
189 idColumn[j] = 1;
190
192 jit.mv(idColumn,column);
193
194 for (int k=0; k<coorddim; k++)
195 jitAsFieldMatrix[k][j] = column[k];
196
197 }
198
199
201 FMatrixHelp::multMatrix( jtAsFieldMatrix, jitAsFieldMatrix, id );
202 bool isId = true;
203 for( int j = 0; j < mydim; ++j )
204 for( int k = 0; k < mydim; ++k )
205 isId &= (std::abs( id[ j ][ k ] - (j == k ? 1 : 0) ) < 1e-8);
206 if( !isId)
207 {
208 std::cerr << "Error: jacobianTransposed and jacobianInverseTransposed are not inverse to each other." << std::endl;
209 std::cout << " id != [ ";
210 for( int j = 0; j < mydim; ++j )
211 std::cout << (j > 0 ? " | " : "") << id[ j ];
212 std::cout << " ]" << std::endl;
213 pass = false;
214 }
215
216 // Test whether integrationElement returns something nonnegative
217 if( geometry.integrationElement( x ) < 0 ) {
218 std::cerr << "Error: Negative integrationElement found." << std::endl;
219 pass = false;
220 }
221
223 for( int i = 0; i < mydim; ++i )
224 for( int j = 0; j < mydim; ++j )
225 for( int k = 0; k < coorddim; ++k )
226 jtj[ i ][ j ] += jtAsFieldMatrix[ i ][ k ] * jtAsFieldMatrix[ j ][ k ];
227
228 if( std::abs( std::sqrt( jtj.determinant() ) - geometry.integrationElement( x ) ) > 1e-8 ) {
229 std::cerr << "Error: integrationElement is not consistent with jacobianTransposed." << std::endl;
230 pass = false;
231 }
232 if (geometry.affine())
233 if( std::abs( geometry.volume() - refElement.volume()*geometry.integrationElement( x ) ) > 1e-8 ) {
234 std::cerr << "Error: volume is not consistent with jacobianTransposed." << std::endl;
235 pass = false;
236 }
237 }
238
239 return pass;
240 }
241
242}
243
244#endif // #ifndef DUNE_CHECK_GEOMETRY_HH
field_type determinant() const
calculates the determinant of this matrix
Base class for Dune-Exceptions.
Definition: exceptions.hh:94
vector space out of a tensor product of fields.
Definition: fvector.hh:93
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:268
bool isNone() const
Return true if entity is a singular of any dimension.
Definition: type.hh:560
generic geometry implementation based on corner coordinates
Definition: multilineargeometry.hh:190
static const int mydimension
geometry dimension
Definition: multilineargeometry.hh:198
Dune::GeometryType type() const
obtain the name of the reference element
Definition: multilineargeometry.hh:266
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian
Definition: multilineargeometry.hh:365
GlobalCoordinate global(const LocalCoordinate &local) const
evaluate the mapping
Definition: multilineargeometry.hh:287
GlobalCoordinate center() const
obtain the centroid of the mapping's image
Definition: multilineargeometry.hh:279
GlobalCoordinate corner(int i) const
obtain coordinates of the i-th corner
Definition: multilineargeometry.hh:272
ct ctype
coordinate type
Definition: multilineargeometry.hh:195
static const int coorddimension
coordinate dimension
Definition: multilineargeometry.hh:200
int corners() const
obtain number of corners of the corresponding reference element
Definition: multilineargeometry.hh:269
ctype volume() const
obtain the volume of the mapping's image
Definition: multilineargeometry.hh:351
ctype integrationElement(const LocalCoordinate &local) const
obtain the integration element
Definition: multilineargeometry.hh:338
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian's inverse
Definition: multilineargeometry.hh:671
bool affine() const
is this mapping affine?
Definition: multilineargeometry.hh:259
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:97
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:225
This class provides access to geometric and topological properties of a reference element.
Definition: referenceelements.hh:354
int size(int c) const
number of subentities of codimension c
Definition: referenceelements.hh:381
const FieldVector< ctype, dim > & position(int i, int c) const
position of the barycenter of entity (i,c)
Definition: referenceelements.hh:449
ctype volume() const
obtain the volume of the reference element
Definition: referenceelements.hh:486
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune namespace.
Definition: alignment.hh:11
bool checkGeometry(const TestGeometry &geometry)
Static and dynamic checks for all features of a Geometry.
Definition: checkgeometry.hh:32
static const ReferenceElement< ctype, dim > & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:757
Traits for type conversions and type information.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)