DUNE-FEM (unstable)

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// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5
6#ifndef DUNE_CHECK_GEOMETRY_HH
7#define DUNE_CHECK_GEOMETRY_HH
8
9#include <limits>
10
14
15#include <dune/geometry/referenceelements.hh>
17#include <dune/geometry/referenceelements.hh>
18
19namespace Dune
20{
21
22 namespace Impl
23 {
24
25 // Helper function to convert a matrix to a FieldMatrix
26 // using only the mv method.
27 template<class ctype, int rows, int cols, class M>
28 Dune::FieldMatrix<ctype, rows, cols> toFieldMatrix(const M& m){
30 for (int j=0; j<cols; j++) {
31 FieldVector<ctype,cols> idColumn(0);
32 idColumn[j] = 1;
33 FieldVector<ctype,rows> column;
34 m.mv(idColumn,column);
35 for (int k=0; k<rows; k++)
36 result[k][j] = column[k];
37 }
38 return result;
39 }
40
41 // Helper function to check if FieldMatrix is an identity matrix
42 template<class ctype, int rows, int cols>
43 bool isIdentity(const Dune::FieldMatrix<ctype, rows, cols>& m, double tolerance){
44 if (rows != cols)
45 return false;
46 for(int i=0; i<rows; ++i)
47 for(int j=0; j<rows; ++j)
48 if (std::abs(m[i][j] - (i == j ? 1 : 0)) >= tolerance)
49 return false;
50 return true;
51 }
52
53 }
54
66 template <class TestGeometry>
67 bool checkGeometry ( const TestGeometry& geometry )
68 {
69 using std::sqrt;
70 using std::abs;
71 bool pass = true;
72
74 // Extract all static information
76
77 // dimension of the corresponding reference element
78 static const int mydim = TestGeometry::mydimension;
79
80 // dimension of the world space
81 static const int coorddim = TestGeometry::coorddimension;
82
83 // type used for coordinate coefficients
84 typedef typename TestGeometry::ctype ctype;
85
86 // vector type used for points in the domain
87 [[maybe_unused]] typedef typename TestGeometry::LocalCoordinate LocalCoordinate;
88
89 // vector type used for image points
90 [[maybe_unused]] typedef typename TestGeometry::GlobalCoordinate GlobalCoordinate;
91
92 // Matrix-like type for the return value of the jacobianTransposed method
93 typedef typename TestGeometry::JacobianTransposed JacobianTransposed;
94
95 // Matrix-like type for the return value of the jacobianInverseTransposed method
96 typedef typename TestGeometry::JacobianInverseTransposed JacobianInverseTransposed;
97
98 // Matrix-like type for the return value of the jacobian method
99 typedef typename TestGeometry::Jacobian Jacobian;
100
101 // Matrix-like type for the return value of the jacobianInverse method
102 typedef typename TestGeometry::JacobianInverse JacobianInverse;
103
104 const ctype tolerance = std::sqrt( std::numeric_limits< ctype >::epsilon() );
105
107
108 const int corners = geometry.corners();
109 if( corners == 0 )
110 DUNE_THROW( Exception, "A geometry must have at least one corner." );
111
112 GlobalCoordinate cornerAvg ( 0 );
113 for( int i = 0; i < corners; ++i )
114 cornerAvg += geometry.corner( i );
115 cornerAvg /= ctype( corners );
116
117 const GlobalCoordinate center = geometry.center();
118
119 if( corners > 1 )
120 {
121 // if we have more than one corner, no corner may coincide with their average or the center
122 for( int i = 0; i < corners; ++i )
123 {
124 const GlobalCoordinate corner = geometry.corner( i );
125 if( (corner - center).two_norm() <= tolerance )
126 {
127 std::cerr << "Error: Geometry has multiple corners, but one corner coincides with the center." << std::endl;
128 pass = false;
129 }
130
131 if( (corner - cornerAvg).two_norm() <= tolerance )
132 {
133 std::cerr << "Error: Geometry has multiple corners, but one corner coincides with their average." << std::endl;
134 pass = false;
135 }
136 }
137 }
138 else
139 {
140 // the single corner must coincide with the center
141 if( (center - cornerAvg).two_norm() > tolerance )
142 {
143 std::cerr << "Error: Geometry has a single corner (" << cornerAvg << "), but it does not coincide with the center (" << center << ")." << std::endl;
144 pass = false;
145 }
146 }
147
148 const ctype volume = geometry.volume();
149 if( volume < tolerance )
150 {
151 std::cerr << "Error: Geometry has nearly vanishing volume (" << volume << ")" << std::endl;
152 pass = false;
153 }
154
156
157 const GeometryType type = geometry.type();
158 if( type.isNone() )
159 return pass;
160
161 // make sure the reference element type lookup works
162 ReferenceElement< TestGeometry > refElement = referenceElement( geometry );
163
164 // Test whether the return value of the method 'center' corresponds to the center of the
165 // reference element. That is the current definition of the method.
166 if( (center - geometry.global( refElement.position( 0, 0 ) )).two_norm() > tolerance )
167 DUNE_THROW( Exception, "center() is not consistent with global(refElem.position(0,0))." );
168
170 // Test whether the number and placement of the corners is consistent
171 // with the corners of the corresponding reference element.
173
174 if( refElement.size( mydim ) == corners )
175 {
176 for( int i = 0; i < geometry.corners(); ++i )
177 {
178 if( (geometry.corner( i ) - geometry.global( refElement.position( i, mydim ) )).two_norm() > tolerance )
179 {
180 std::cerr << "Error: Methods corner and global are inconsistent." << std::endl;
181 pass = false;
182 }
183 }
184 }
185 else
186 {
187 std::cerr << "Error: Incorrect number of corners (" << geometry.corners() << ", should be " << refElement.size( mydim ) << ")." << std::endl;
188 pass = false;
189 }
190
192 // Use a quadrature rule as a set of test points and loop over them
195 for (const auto& ip : quadrature)
196 {
197 const typename TestGeometry::LocalCoordinate &x = ip.position();
198
199 // Test whether the methods 'local' and 'global' are inverse to each other
200 if ( (x - geometry.local( geometry.global( x ) )).two_norm() > tolerance ) {
201 std::cerr << "Error: global and local are not inverse to each other." << std::endl;
202 pass = false;
203 }
204
205 const JacobianTransposed &Jt = geometry.jacobianTransposed( x );
206 const JacobianInverseTransposed &Jit = geometry.jacobianInverseTransposed( x );
207 const Jacobian &J = geometry.jacobian( x );
208 const JacobianInverse &Ji = geometry.jacobianInverse( x );
209
210 // Transform to FieldMatrix, so we can have coefficient access and other goodies
211 auto JtAsFieldMatrix = Impl::toFieldMatrix< ctype, mydim, coorddim >(Jt);
212 auto JitAsFieldMatrix = Impl::toFieldMatrix< ctype, coorddim, mydim >(Jit);
213 auto JAsFieldMatrix = Impl::toFieldMatrix< ctype, coorddim, mydim >(J);
214 auto JiAsFieldMatrix = Impl::toFieldMatrix< ctype, mydim, coorddim >(Ji);
215
216
217 // Test whether the methods 'jacobianTransposed' and 'jacobianInverseTransposed'
218 // return matrices that are inverse to each other.
219 {
220 FieldMatrix< ctype, mydim, mydim > id = JtAsFieldMatrix * JitAsFieldMatrix;
221 if(not Impl::isIdentity(id, tolerance))
222 {
223 std::cerr << "Error: jacobianTransposed and jacobianInverseTransposed are not inverse to each other." << std::endl;
224 std::cout << " id != [ ";
225 for( int j = 0; j < mydim; ++j )
226 std::cout << (j > 0 ? " | " : "") << id[ j ];
227 std::cout << " ]" << std::endl;
228 pass = false;
229 }
230 }
231
232 // Test whether the methods 'jacobian' and 'jacobianInverse'
233 // return matrices that are inverse to each other.
234 {
235 FieldMatrix< ctype, mydim, mydim > id = JiAsFieldMatrix * JAsFieldMatrix;
236 if(not Impl::isIdentity(id, tolerance))
237 {
238 std::cerr << "Error: jacobian and jacobianInverse are not inverse to each other." << std::endl;
239 std::cout << " id != [ ";
240 for( int j = 0; j < mydim; ++j )
241 std::cout << (j > 0 ? " | " : "") << id[ j ];
242 std::cout << " ]" << std::endl;
243 pass = false;
244 }
245 }
246
247 // Test whether the methods 'jacobianTransposed' and 'jacobianInverseTransposed'
248 // are the transposed of 'jacobian' and 'jacobianInverse', respectively.
249 {
250 if( (JtAsFieldMatrix - JAsFieldMatrix.transposed()).infinity_norm() != 0 )
251 {
252 std::cerr << "Error: jacobian and jacobianTransposed are not transposed to each other." << std::endl;
253 pass = false;
254 }
255 }
256 {
257 if( (JitAsFieldMatrix - JiAsFieldMatrix.transposed()).infinity_norm() != 0 )
258 {
259 std::cerr << "Error: jacobianInverse and jacobianInverseTransposed are not transposed to each other." << std::endl;
260 pass = false;
261 }
262 }
263
264
265 // Test whether integrationElement returns something nonnegative
266 if( geometry.integrationElement( x ) < 0 ) {
267 std::cerr << "Error: Negative integrationElement found." << std::endl;
268 pass = false;
269 }
270
272 for( int i = 0; i < mydim; ++i )
273 for( int j = 0; j < mydim; ++j )
274 for( int k = 0; k < coorddim; ++k )
275 jtj[ i ][ j ] += JtAsFieldMatrix[ i ][ k ] * JtAsFieldMatrix[ j ][ k ];
276
277 if( abs( sqrt( jtj.determinant() ) - geometry.integrationElement( x ) ) > tolerance ) {
278 std::cerr << "Error: integrationElement is not consistent with jacobianTransposed." << std::endl;
279 pass = false;
280 }
281 if (geometry.affine())
282 if( abs( geometry.volume() - refElement.volume()*geometry.integrationElement( x ) ) > tolerance ) {
283 std::cerr << "Error: volume is not consistent with jacobianTransposed." << std::endl;
284 pass = false;
285 }
286 }
287
288 return pass;
289 }
290
291}
292
293#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:96
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
constexpr bool isNone() const
Return true if entity is a singular of any dimension.
Definition: type.hh:355
generic geometry implementation based on corner coordinates
Definition: multilineargeometry.hh:181
static const int mydimension
geometry dimension
Definition: multilineargeometry.hh:189
Dune::GeometryType type() const
obtain the name of the reference element
Definition: multilineargeometry.hh:269
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian
Definition: multilineargeometry.hh:377
GlobalCoordinate global(const LocalCoordinate &local) const
evaluate the mapping
Definition: multilineargeometry.hh:290
GlobalCoordinate center() const
obtain the centroid of the mapping's image
Definition: multilineargeometry.hh:282
GlobalCoordinate corner(int i) const
obtain coordinates of the i-th corner
Definition: multilineargeometry.hh:275
ct ctype
coordinate type
Definition: multilineargeometry.hh:186
static const int coorddimension
coordinate dimension
Definition: multilineargeometry.hh:191
int corners() const
obtain number of corners of the corresponding reference element
Definition: multilineargeometry.hh:272
Volume volume() const
obtain the volume of the mapping's image
Definition: multilineargeometry.hh:363
JacobianInverse jacobianInverse(const LocalCoordinate &local) const
Obtain the Jacobian's inverse.
Definition: multilineargeometry.hh:418
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
obtain the transposed of the Jacobian's inverse
Definition: multilineargeometry.hh:738
bool affine() const
is this mapping affine?
Definition: multilineargeometry.hh:262
Volume integrationElement(const LocalCoordinate &local) const
obtain the integration element
Definition: multilineargeometry.hh:350
Jacobian jacobian(const LocalCoordinate &local) const
Obtain the Jacobian.
Definition: multilineargeometry.hh:407
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:214
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:326
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:218
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
Dune namespace.
Definition: alignedallocator.hh:13
bool checkGeometry(const TestGeometry &geometry)
Static and dynamic checks for all features of a Geometry.
Definition: checkgeometry.hh:67
Traits for type conversions and type information.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)