This file is indexed.

/usr/include/dune/geometry/test/checkgeometry.hh is in libdune-geometry-dev 2.5.0-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