Dune Core Modules (2.8.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 using std::sqrt;
35 using std::abs;
36 bool pass = true;
37
39 // Extract all static information
41
42 // dimension of the corresponding reference element
43 static const int mydim = TestGeometry::mydimension;
44
45 // dimension of the world space
46 static const int coorddim = TestGeometry::coorddimension;
47
48 // type used for coordinate coefficients
49 typedef typename TestGeometry::ctype ctype;
50
51 // vector type used for points in the domain
52 [[maybe_unused]] typedef typename TestGeometry::LocalCoordinate LocalCoordinate;
53
54 // vector type used for image points
55 [[maybe_unused]] typedef typename TestGeometry::GlobalCoordinate GlobalCoordinate;
56
57 // Matrix-like type for the return value of the jacobianTransposed method
58 typedef typename TestGeometry::JacobianTransposed JacobianTransposed;
59
60 // Matrix-like type for the return value of the jacobianInverseTransposed method
61 typedef typename TestGeometry::JacobianInverseTransposed JacobianInverseTransposed;
62
63 const ctype tolerance = std::sqrt( std::numeric_limits< ctype >::epsilon() );
64
66
67 const int corners = geometry.corners();
68 if( corners == 0 )
69 DUNE_THROW( Exception, "A geometry must have at least one corner." );
70
71 GlobalCoordinate cornerAvg ( 0 );
72 for( int i = 0; i < corners; ++i )
73 cornerAvg += geometry.corner( i );
74 cornerAvg /= ctype( corners );
75
76 const GlobalCoordinate center = geometry.center();
77
78 if( corners > 1 )
79 {
80 // if we have more than one corner, no corner may coincide with their average or the center
81 for( int i = 0; i < corners; ++i )
82 {
83 const GlobalCoordinate corner = geometry.corner( i );
84 if( (corner - center).two_norm() <= tolerance )
85 {
86 std::cerr << "Error: Geometry has multiple corners, but one corner coincides with the center." << std::endl;
87 pass = false;
88 }
89
90 if( (corner - cornerAvg).two_norm() <= tolerance )
91 {
92 std::cerr << "Error: Geometry has multiple corners, but one corner coincides with their average." << std::endl;
93 pass = false;
94 }
95 }
96 }
97 else
98 {
99 // the single corner must coincide with the center
100 if( (center - cornerAvg).two_norm() > tolerance )
101 {
102 std::cerr << "Error: Geometry has a single corner (" << cornerAvg << "), but it does not coincide with the center (" << center << ")." << std::endl;
103 pass = false;
104 }
105 }
106
107 const ctype volume = geometry.volume();
108 if( volume < tolerance )
109 {
110 std::cerr << "Error: Geometry has nearly vanishing volume (" << volume << ")" << std::endl;
111 pass = false;
112 }
113
115
116 const GeometryType type = geometry.type();
117 if( type.isNone() )
118 return pass;
119
120 // make sure the reference element type lookup works
121 ReferenceElement< TestGeometry > refElement = referenceElement( geometry );
122
123 // Test whether the return value of the method 'center' corresponds to the center of the
124 // reference element. That is the current definition of the method.
125 if( (center - geometry.global( refElement.position( 0, 0 ) )).two_norm() > tolerance )
126 DUNE_THROW( Exception, "center() is not consistent with global(refElem.position(0,0))." );
127
129 // Test whether the number and placement of the corners is consistent
130 // with the corners of the corresponding reference element.
132
133 if( refElement.size( mydim ) == corners )
134 {
135 for( int i = 0; i < geometry.corners(); ++i )
136 {
137 if( (geometry.corner( i ) - geometry.global( refElement.position( i, mydim ) )).two_norm() > tolerance )
138 {
139 std::cerr << "Error: Methods corner and global are inconsistent." << std::endl;
140 pass = false;
141 }
142 }
143 }
144 else
145 {
146 std::cerr << "Error: Incorrect number of corners (" << geometry.corners() << ", should be " << refElement.size( mydim ) << ")." << std::endl;
147 pass = false;
148 }
149
151 // Use a quadrature rule as a set of test points and loop over them
154 for (const auto& ip : quadrature)
155 {
156 const typename TestGeometry::LocalCoordinate &x = ip.position();
157
158 // Test whether the methods 'local' and 'global' are inverse to each other
159 if ( (x - geometry.local( geometry.global( x ) )).two_norm() > tolerance ) {
160 std::cerr << "Error: global and local are not inverse to each other." << std::endl;
161 pass = false;
162 }
163
164 // Test whether the methods 'jacobianTransposed' and 'jacobianInverseTransposed'
165 // return matrices that are inverse to each other.
166 const JacobianTransposed &jt = geometry.jacobianTransposed( x );
167 const JacobianInverseTransposed &jit = geometry.jacobianInverseTransposed( x );
168
169 // Transform to FieldMatrix, so we can have coefficent access and other goodies
170 // We need some black magic for the transformation, because there is no
171 // official easy way yet.
172 // The following code does the transformation by multiplying jt and jit from
173 // the right by identity matrices. That way, only the mv method is used.
175 for (int j=0; j<coorddim; j++) {
176
177 FieldVector<ctype,coorddim> idColumn(0);
178 idColumn[j] = 1;
179
181 jt.mv(idColumn,column);
182
183 for (int k=0; k<mydim; k++)
184 jtAsFieldMatrix[k][j] = column[k];
185
186 }
187
189 for (int j=0; j<mydim; j++) {
190
191 FieldVector<ctype,mydim> idColumn(0);
192 idColumn[j] = 1;
193
195 jit.mv(idColumn,column);
196
197 for (int k=0; k<coorddim; k++)
198 jitAsFieldMatrix[k][j] = column[k];
199
200 }
201
202
204 FMatrixHelp::multMatrix( jtAsFieldMatrix, jitAsFieldMatrix, id );
205 bool isId = true;
206 for( int j = 0; j < mydim; ++j )
207 for( int k = 0; k < mydim; ++k )
208 isId &= (std::abs( id[ j ][ k ] - (j == k ? 1 : 0) ) < tolerance);
209 if( !isId)
210 {
211 std::cerr << "Error: jacobianTransposed and jacobianInverseTransposed are not inverse to each other." << std::endl;
212 std::cout << " id != [ ";
213 for( int j = 0; j < mydim; ++j )
214 std::cout << (j > 0 ? " | " : "") << id[ j ];
215 std::cout << " ]" << std::endl;
216 pass = false;
217 }
218
219 // Test whether integrationElement returns something nonnegative
220 if( geometry.integrationElement( x ) < 0 ) {
221 std::cerr << "Error: Negative integrationElement found." << std::endl;
222 pass = false;
223 }
224
226 for( int i = 0; i < mydim; ++i )
227 for( int j = 0; j < mydim; ++j )
228 for( int k = 0; k < coorddim; ++k )
229 jtj[ i ][ j ] += jtAsFieldMatrix[ i ][ k ] * jtAsFieldMatrix[ j ][ k ];
230
231 if( abs( sqrt( jtj.determinant() ) - geometry.integrationElement( x ) ) > tolerance ) {
232 std::cerr << "Error: integrationElement is not consistent with jacobianTransposed." << std::endl;
233 pass = false;
234 }
235 if (geometry.affine())
236 if( abs( geometry.volume() - refElement.volume()*geometry.integrationElement( x ) ) > tolerance ) {
237 std::cerr << "Error: volume is not consistent with jacobianTransposed." << std::endl;
238 pass = false;
239 }
240 }
241
242 return pass;
243 }
244
245}
246
247#endif // #ifndef DUNE_CHECK_GEOMETRY_HH
field_type determinant(bool doPivoting=true) 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:95
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:123
constexpr bool isNone() const
Return true if entity is a singular of any dimension.
Definition: type.hh:364
generic geometry implementation based on corner coordinates
Definition: multilineargeometry.hh:179
static const int mydimension
geometry dimension
Definition: multilineargeometry.hh:187
Dune::GeometryType type() const
obtain the name of the reference element
Definition: multilineargeometry.hh:261
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian
Definition: multilineargeometry.hh:369
GlobalCoordinate global(const LocalCoordinate &local) const
evaluate the mapping
Definition: multilineargeometry.hh:282
GlobalCoordinate center() const
obtain the centroid of the mapping's image
Definition: multilineargeometry.hh:274
GlobalCoordinate corner(int i) const
obtain coordinates of the i-th corner
Definition: multilineargeometry.hh:267
ct ctype
coordinate type
Definition: multilineargeometry.hh:184
static const int coorddimension
coordinate dimension
Definition: multilineargeometry.hh:189
int corners() const
obtain number of corners of the corresponding reference element
Definition: multilineargeometry.hh:264
Volume volume() const
obtain the volume of the mapping's image
Definition: multilineargeometry.hh:355
ctype integrationElement(const LocalCoordinate &local) const
obtain the integration element
Definition: multilineargeometry.hh:342
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian's inverse
Definition: multilineargeometry.hh:683
bool affine() const
is this mapping affine?
Definition: multilineargeometry.hh:254
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:152
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:280
Traits for type conversions and type information.
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
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
unspecified-type ReferenceElement
Returns the type of reference element for the argument type T.
Definition: referenceelements.hh:495
Dune namespace.
Definition: alignedallocator.hh:11
bool checkGeometry(const TestGeometry &geometry)
Static and dynamic checks for all features of a Geometry.
Definition: checkgeometry.hh:32
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 22, 23:30, 2024)