This file is indexed.

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

    ////////////////////////////////////////////////////////////////////////

    GeometryType type = geometry.type();

    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.
    const FieldVector<ctype, coorddim> center = geometry.global(refElement.position(0,0));
    if( std::abs( (geometry.center() - center).two_norm() ) > 1e-8 )
      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 ) == geometry.corners() )
    {
      for( int i = 0; i < geometry.corners(); ++i )
      {
        if( (geometry.corner( i ) - geometry.global( refElement.position( i, mydim ) )).two_norm() > 1e-8 ) {
          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 (size_t i = 0; i < quadrature.size(); ++i)
    {
      const typename TestGeometry::LocalCoordinate &x = quadrature[i].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