1 #ifndef DUNE_FEM_GRIDPART_TEST_CHECKINTERSECTIONS_HH 2 #define DUNE_FEM_GRIDPART_TEST_CHECKINTERSECTIONS_HH 3 4 //- dune-common includes 5 #include <dune/common/exceptions.hh> 6 #include <dune/common/stdstreams.hh> 7 8 //- dune-geometry includes 9 #include <dune/geometry/quadraturerules.hh> 10 #include <dune/geometry/referenceelements.hh> 11 12 //- dune-grid includes 13 #include <dune/grid/test/checkgeometry.hh> 14 #include <dune/grid/test/checkintersectionit.hh> 15 16 //- dune-fem includes 17 #include <dune/fem/gridpart/test/failure.hh> 18 #include <dune/fem/quadrature/cachingquadrature.hh> 19 #include <dune/fem/quadrature/intersectionquadrature.hh> 20 21 namespace Dune 22 { 23 namespace Fem 24 { 25 template< class GridPartType, class FailureHandler > 26 struct CheckIntersections 27 { 28 /** \brief type of intersection iterator */ 29 typedef typename GridPartType::IntersectionIteratorType IntersectionIteratorType; 30 /** \brief type of intersection */ 31 typedef typename GridPartType::IntersectionType IntersectionType; 32 /** \brief loal geometry type */ 33 typedef typename IntersectionType::Geometry GeometryType; 34 /** \brief loal geometry type */ 35 typedef typename IntersectionType::LocalGeometry LocalGeometryType; 36 /** \brief entity type */ 37 typedef typename IntersectionType::Entity EntityType; 38 39 /** \brief failure type */ 40 struct AssignmentFailure; 41 /** \brief failure type */ 42 struct SumNormalFailure; 43 44 private: 45 // static tests 46 static_assert( ( std::is_same< typename IntersectionType::ctype, 47 typename GridPartType::ctype 48 >::value 49 ), "Conflicting types for Intersection" ); 50 51 static_assert( ( static_cast< int>( IntersectionType::dimensionworld ) 52 == static_cast< int >( GridPartType::dimensionworld ) 53 ), "IntersectionIterator has wrong dimensionworld" ); 54 55 // check intersection iterator assignment 56 template< class IntersectionIterator > 57 static void checkIntersectionIteratorAssignment( const IntersectionIterator &it, 58 FailureHandler &failureHandler ); 59 60 // check intersection iterator assignment 61 static void checkOuterNormals ( const EntityType &entity, 62 const GridPartType &gridPart, 63 FailureHandler &failureHandler ); 64 65 // check local geometries 66 static void checkLocalGeometries ( const IntersectionType &intersection, 67 const EntityType &entity, 68 FailureHandler &failureHandler ); 69 70 public: checkDune::Fem::CheckIntersections71 static void check ( const GridPartType &gridPart, FailureHandler &failureHandler ) 72 { 73 for( const auto& entity : elements( gridPart ) ) 74 { 75 for( auto it = gridPart.ibegin( entity ); it != gridPart.iend( entity ); ++it ) 76 { 77 const auto& intersection = *it; 78 79 checkIntersection( intersection, false ); 80 81 // create intersection caching quadrature to check twists 82 IntersectionQuadrature< CachingQuadrature< GridPartType, 1 >, true > 83 inter( gridPart, intersection, 2 ); 84 85 checkIntersectionIteratorAssignment( it, failureHandler ); 86 checkLocalGeometries( intersection, entity, failureHandler ); 87 } 88 89 checkOuterNormals( entity, gridPart, failureHandler ); 90 if( SumNormalFailure::failed() ) 91 failureHandler( SumNormalFailure::instance() ); 92 } 93 } 94 }; 95 96 97 98 // Implementation of CheckIntersections 99 // ------------------------------------ 100 101 template< class GridPartType, class FailureHandler > 102 template< class IntersectionIterator > 103 inline void CheckIntersections< GridPartType, FailureHandler > checkIntersectionIteratorAssignment(const IntersectionIterator & it,FailureHandler & failureHandler)104 ::checkIntersectionIteratorAssignment( const IntersectionIterator &it, FailureHandler &failureHandler ) 105 { 106 // type of intersection iterator 107 typedef IntersectionIterator IntersectionIteratorType; 108 109 AssignmentFailure failure; 110 // assert assignment 111 IntersectionIteratorType it2 = it; 112 if( it != it2 ) 113 failureHandler( failure ); 114 if( it->inside() != it2->inside() ) 115 failureHandler( failure ); 116 if( it->neighbor() ) 117 if( !it->neighbor() || it->outside() != it2->outside() ) 118 failureHandler( failure ); 119 120 // now increment second iterator 121 ++it2; 122 if( it == it2 ) 123 failureHandler( failure ); 124 } 125 126 127 template< class GridPartType, class FailureHandler > 128 inline void CheckIntersections< GridPartType, FailureHandler > checkOuterNormals(const EntityType & entity,const GridPartType & gridPart,FailureHandler & failureHandler)129 ::checkOuterNormals ( const EntityType &entity, 130 const GridPartType &gridPart, 131 FailureHandler &failureHandler ) 132 { 133 if (!entity.geometry().affine() && (int)entity.geometry().coorddimension>(int)entity.dimension) 134 // this test is wrong in the case of a non affine embedded surface 135 return; 136 // global coordinate type 137 typedef typename GeometryType::GlobalCoordinate GlobalCoordinateType; 138 // single coordinate type 139 typedef typename GridPartType::ctype ctype; 140 141 // compute sum of all integration outer normals 142 GlobalCoordinateType sum( 0 ); 143 for( const auto& intersection : intersections( gridPart, entity ) ) 144 { 145 // get Gauss quadrature points 146 const auto& quadrature = Dune::QuadratureRules< ctype, LocalGeometryType::mydimension >::rule( intersection.type(), 3 ); 147 148 for( const auto& qp : quadrature ) 149 { 150 auto integrationOuterNormal = intersection.integrationOuterNormal( qp.position() ); 151 sum.axpy( qp.weight(), integrationOuterNormal ); 152 } 153 } 154 155 // check wheter error is within fixed tolerance 156 double error = sum.two_norm(); 157 const double tolerance = SumNormalFailure::tolerance; 158 if( (sum.two_norm() > tolerance ) && ( entity.partitionType() != Dune::GhostEntity) ) 159 SumNormalFailure::add( entity, error ); 160 } 161 162 163 template< class GridPartType, class FailureHandler > 164 inline void CheckIntersections< GridPartType, FailureHandler > checkLocalGeometries(const IntersectionType & intersection,const EntityType & entity,FailureHandler & failureHandler)165 ::checkLocalGeometries ( const IntersectionType &intersection, 166 const EntityType &entity, 167 FailureHandler &failureHandler ) 168 { 169 auto geometryInInside = intersection.geometryInInside(); 170 geometryInInside.type(); 171 Dune::checkLocalGeometry( geometryInInside, entity.type(), "geometryInInside" ); 172 173 if( intersection.neighbor() ) 174 { 175 LocalGeometryType geometryInOutside = intersection.geometryInOutside(); 176 Dune::checkLocalGeometry( geometryInOutside, entity.type(), "geometryInOutside" ); 177 } 178 } 179 180 181 182 // Implementation of Failures used in CheckIntersections 183 // ----------------------------------------------------- 184 185 template< class GridPartType, class FailureHandler > 186 struct CheckIntersections< GridPartType, FailureHandler >::AssignmentFailure 187 : public Failure 188 { writeToDune::Fem::CheckIntersections::AssignmentFailure189 virtual void writeTo ( std::ostream &out ) const 190 { 191 out << __FILE__ 192 << ":" << __LINE__ << ": Failure :" 193 << "Assignment failure"; 194 } 195 }; 196 197 198 199 template< class GridPartType, class FailureHandler > 200 struct CheckIntersections< GridPartType, FailureHandler >::SumNormalFailure 201 : public Failure 202 { 203 /** \brief entity type */ 204 typedef typename GridPartType::template Codim< 0 >::EntityType EntityType; 205 206 /** \brief tolerance */ 207 static const double tolerance; 208 209 protected: 210 /** \brief constructor */ SumNormalFailureDune::Fem::CheckIntersections::SumNormalFailure211 SumNormalFailure () 212 : fails_( 0 ) 213 {} 214 215 public: 216 /** \brief return instance of singleton */ instanceDune::Fem::CheckIntersections::SumNormalFailure217 static SumNormalFailure &instance () 218 { 219 static SumNormalFailure instance_; 220 return instance_; 221 } 222 223 /** \brief return true, if actual failure */ failedDune::Fem::CheckIntersections::SumNormalFailure224 static bool failed () 225 { 226 return ( instance().fails_ > 0 ); 227 } 228 229 /** \brief add failed entity */ addDune::Fem::CheckIntersections::SumNormalFailure230 static void add ( const EntityType &entity, double error ) 231 { 232 instance().fails_ += 1; 233 } 234 235 /** \brief write message to stream */ writeToDune::Fem::CheckIntersections::SumNormalFailure236 virtual void writeTo ( std::ostream &out ) const 237 { 238 assert( failed() ); 239 out << __FILE__ 240 << ":" << __LINE__ << ": Failure :" 241 << "Summation of normals fails for " 242 << instance().fails_ << " entities in grid part"; 243 } 244 245 private: 246 int fails_; 247 }; 248 249 250 251 template< class GridPartType, class FailureHandler > 252 const double CheckIntersections< GridPartType, FailureHandler >::SumNormalFailure::tolerance = 1e-8; 253 254 } // end namespace Fem 255 256 } // end namespace Dune 257 258 #endif // #ifndef DUNE_FEM_GRIDPART_TEST_CHECKINTERSECTIONS_HH 259