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