1 #include <config.h>
2 #include <iostream>
3 #include <set>
4 
5 #include <dune/grid/common/gridinfo.hh>
6 
7 #include <dune/fem/gridpart/leafgridpart.hh>
8 #include <dune/fem/io/parameter.hh>
9 #include <dune/fem/quadrature/cachingquadrature.hh>
10 #include <dune/fem/quadrature/elementquadrature.hh>
11 
12 #include <dune/fem/quadrature/dunequadratures.hh>
13 #include <dune/fem/quadrature/lumpingquadrature.hh>
14 
15 #include "checkleafcodim1.hh"
16 
17 using namespace Dune;
18 using namespace Fem ;
19 
20 struct DuneQuadratures
21 {};
22 
23 template<class GridPartType, int codim, bool caching, class QuadratureTraits = int >
24 struct QuadratureChooser;
25 
26 template<class GridPartType, int codim>
27 struct QuadratureChooser<GridPartType,codim,false,int>
28 {
29   typedef ElementQuadrature< GridPartType, codim > Quadrature;
30 };
31 template<class GridPartType, int codim>
32 struct QuadratureChooser<GridPartType,codim,true,int>
33 {
34   typedef CachingQuadrature< GridPartType, codim > Quadrature;
35 };
36 
37 template<class GridPartType, int codim, class QuadratureTraits >
38 struct QuadratureChooser<GridPartType,codim,false, QuadratureTraits>
39 {
40   typedef ElementQuadrature< GridPartType, codim, Dune::Fem::DuneQuadratureTraits > Quadrature;
41 };
42 template<class GridPartType, int codim, class QuadratureTraits >
43 struct QuadratureChooser<GridPartType,codim,true, QuadratureTraits>
44 {
45   typedef CachingQuadrature< GridPartType, codim, Dune::Fem::DuneQuadratureTraits > Quadrature;
46 };
47 
48 template<class GridPartType, bool caching, class QuadratureTraits = int>
49 class TestCaching
50 {
51 private:
52   // types of iterators, entities and intersection
53   typedef typename GridPartType :: template Codim<0> :: IteratorType IteratorType;
54   typedef typename GridPartType :: IndexSetType  IndexSetType ;
55   typedef typename IteratorType :: Entity EntityType;
56   typedef typename EntityType :: Geometry  Geometry;
57   typedef typename GridPartType :: IntersectionIteratorType IntersectionIteratorType;
58   typedef typename IntersectionIteratorType :: Intersection IntersectionType;
59 
60   typedef typename QuadratureChooser<GridPartType, 0, caching, QuadratureTraits>::Quadrature VolumeQuadratureType;
61   typedef typename QuadratureChooser<GridPartType, 0, false,   QuadratureTraits>::Quadrature ElementQuadratureType;
62   typedef typename QuadratureChooser<GridPartType, 1, caching>::Quadrature FaceQuadratureType;
63 
64   typedef CachingLumpingQuadrature< GridPartType, 0> LumpingQuadratureType;
65 
66   GridPartType & gridPart_;
67   int order_;
68   const double eps_;
69   const bool skipFaces_;
70 
71 public:
72 
TestCaching(GridPartType & gridPart,int order,const double eps=1e-8)73   TestCaching(GridPartType & gridPart, int order, const double eps = 1e-8)
74     : gridPart_(gridPart), order_( std::abs(order) ), eps_(eps), skipFaces_( order < 0 )
75   { }
76 
runTest()77   void runTest()
78   {
79     testElementQuadratures();
80     if( skipFaces_ )
81     {
82       std::cout << "Skipping faces due to fem.skipfaces "<<std::endl;
83       return ;
84     }
85     testFaceQuadratures();
86   }
87 
testElementQuadratures()88   void testElementQuadratures()
89   {
90     IteratorType endit = gridPart_.template end<0>();
91     for (IteratorType it = gridPart_.template begin<0>(); it != endit; ++it)
92     {
93       const EntityType & entity = *it;
94       const Geometry& geo = entity.geometry();
95 
96       VolumeQuadratureType cacheQuad( entity, order_ );
97       ElementQuadratureType elemQuad( entity, order_ );
98       if( cacheQuad.nop() != elemQuad.nop() )
99       {
100         std::cout << " Error: nops not equal: cachequad = " << cacheQuad.nop()
101                   << "   elemQuad = " << elemQuad.nop() << std::endl;
102         return ;
103       }
104 
105       const int quadNop = cacheQuad.nop();
106       for (int qp = 0; qp < quadNop; ++qp)
107       {
108         Dune::FieldVector<typename GridPartType::ctype,GridPartType::dimensionworld> globalElem, globalCache;
109         globalCache  = geo.global( cacheQuad.point( qp ) );
110         globalElem   = geo.global( elemQuad.point( qp ) );
111         if( (globalCache-globalElem).two_norm() > eps_)
112         {
113           std::cout << " Error: x(cache) = " << globalCache << " != "
114                     << globalElem << " = x(elem) " << std::endl;
115           break;
116         }
117       }
118     }
119   }
120 
testFaceQuadratures()121   void testFaceQuadratures()
122   {
123     // just another face check
124     CachingQuadratureTest :: checkLeafsCodimOne( gridPart_, order_ );
125 
126     const IndexSetType& indexSet = gridPart_.indexSet();
127 
128     std::set< int > insideTwists;
129     std::set< int > outsideTwists;
130 
131     IteratorType endit = gridPart_.template end<0>();
132     for (IteratorType it = gridPart_.template begin<0>(); it != endit; ++it)
133     {
134       const EntityType & entity = *it;
135       const IntersectionIteratorType iend = gridPart_.iend(entity);
136       for( IntersectionIteratorType iit = gridPart_.ibegin(entity); iit != iend; ++iit )
137       {
138         const IntersectionType &intersection = *iit;
139         if (intersection.neighbor() && intersection.conforming())
140         {
141           EntityType inside = intersection.inside();
142           FaceQuadratureType faceQuadInner(gridPart_, intersection, order_,
143               FaceQuadratureType::INSIDE);
144 
145           EntityType outside = intersection.outside();
146           FaceQuadratureType faceQuadOuter(gridPart_, intersection, order_,
147               FaceQuadratureType::OUTSIDE);
148 
149           const int faceQuadInner_nop = faceQuadInner.nop();
150           for (int qp = 0; qp < faceQuadInner_nop; ++qp)
151           {
152             typedef typename GridPartType::TwistUtilityType TwistUtilityType;
153             Dune::FieldVector<typename GridPartType::ctype,GridPartType::dimensionworld> globalInside, globalOutside;
154             globalInside = inside.geometry().global(faceQuadInner.point(qp));
155             globalOutside = outside.geometry().global(faceQuadOuter.point(qp));
156             if( (globalInside-globalOutside).two_norm() > eps_)
157             {
158               const int twistInside  = TwistUtilityType::twistInSelf( gridPart_.grid(), intersection );
159               const int twistOutside = TwistUtilityType::twistInNeighbor( gridPart_.grid(), intersection );
160               insideTwists.insert( twistInside );
161               outsideTwists.insert( twistOutside );
162               //std::cout << "Inside  twist = " << twistInside << std::endl;
163               //std::cout << "Outside twist = " << twistOutside << std::endl;
164               std::cout << "On Element " << indexSet.index( entity ) << " with neighbor " << indexSet.index( outside ) << std::endl;
165               std::cout << " Error: x(inside) = " << globalInside << " != "
166                         << globalOutside << " = x(outside) "
167                         << " at intersection.indexInInside() = " << intersection.indexInInside()
168                         << ", intersection.indexInOutside() = " << intersection.indexInOutside() << std::endl;
169               break;
170             }
171           }
172         }
173       }
174     }
175 
176     typedef typename std::set<int> :: iterator iterator;
177     const iterator endin = insideTwists.end();
178     for( iterator it = insideTwists.begin(); it != endin; ++it )
179     {
180       std::cout << "Inside twst: " << (*it) << std::endl;
181     }
182     const iterator endout = outsideTwists.end();
183     for( iterator it = outsideTwists.begin(); it != endout; ++it )
184     {
185       std::cout << "Outside twst: " << (*it) << std::endl;
186     }
187   }
188 
testLumpingQuadrature()189   void testLumpingQuadrature()
190   {
191     IteratorType endit = gridPart_.template end<0>();
192     for (IteratorType it = gridPart_.template begin<0>(); it != endit; ++it)
193     {
194       const EntityType & entity = *it;
195       const Geometry& geo = entity.geometry();
196 
197       LumpingQuadratureType lumpQuad( entity, order_ );
198       const int quadNop = lumpQuad.nop();
199       if( quadNop != int(geo.corners()) )
200       {
201         std::cout << " Error: nops not equal: lumpQuad = " << quadNop
202                   << "   geo.corners = " << geo.corners() << std::endl;
203         return ;
204       }
205 
206       for (int qp = 0; qp < quadNop; ++qp)
207       {
208         Dune::FieldVector<typename GridPartType::ctype,GridPartType::dimensionworld> globalElem, globalCache;
209         globalCache  = geo.global( lumpQuad.point( qp ) );
210         globalElem   = geo.corner( qp );
211         if( (globalCache-globalElem).two_norm() > eps_)
212         {
213           std::cout << " Error: x(cache) = " << globalCache << " != "
214                     << globalElem << " = x(elem) " << std::endl;
215           break;
216         }
217       }
218     }
219   }
220 
221 };
222 
223 // main program
main(int argc,char ** argv)224 int main(int argc, char ** argv)
225 {
226   try
227   {
228     MPIManager :: initialize( argc, argv );
229     // *** Initialization
230     Parameter::append(argc, argv);
231     if (argc==2) Parameter::append(argv[1]);
232     else Parameter::append("parameter");
233 
234     typedef Dune::GridSelector::GridType GridType;
235     std::string filename;
236     std::stringstream fileKey ;
237     fileKey << "fem.gridfile" << GridType :: dimension ;
238     Parameter::get( fileKey.str(), filename);
239     Dune::GridPtr< GridType > gridptr( filename );
240     GridType &grid = *gridptr;
241 
242     int maxlevel = Parameter::getValue<int>("fem.maxlevel");
243     //grid.globalRefine( maxlevel );
244     //Dune::gridinfo(grid);
245 
246     int quadOrder = Parameter::getValue<int>("fem.quadorder");
247 
248     if ( Parameter::getValue<bool>("fem.skipfaces", false ) )
249       quadOrder = -quadOrder ;
250 
251     int nonConformOrigin = Parameter::getValue< int > ( "poisson.nonConformOrigin", 0 );
252     if ( nonConformOrigin )
253     {
254       const int refineelement = 1 ;
255       std::cout << "Create local refined grid" << std::endl;
256       for (int i=0;i<nonConformOrigin;++i)
257       {
258         if( grid.comm().rank() == 0)
259         {
260           typedef GridType::Codim<0>::LeafIterator IteratorType;
261           IteratorType endit = grid.leafend<0>();
262           for(IteratorType it = grid.leafbegin<0>(); it != endit ; ++it)
263           {
264             const IteratorType :: Entity & entity = *it ;
265             const IteratorType :: Entity :: Geometry& geo = entity.geometry();
266             if (geo.center().two_norm() < 0.5)
267             {
268               grid.mark(refineelement, entity );
269               std::cout << "mark" << std::endl;
270             }
271           }
272         }
273         grid.preAdapt();
274         grid.adapt();
275         grid.postAdapt();
276       }
277     }
278 
279     //typedef HierarchicGridPart< GridType > GridPartType;
280     typedef Dune::Fem::LeafGridPart< GridType > GridPartType;
281     GridPartType gridPart( grid );
282 
283     const double eps = 1e-8;
284 
285     for(int l=0; l<=maxlevel; ++l )
286     {
287       {
288         std::cout << "Testing default ElementQuadratures: " << std::endl;
289         TestCaching<GridPartType,false> testCaching(gridPart, quadOrder, eps);
290         testCaching.runTest();
291       }
292       {
293         std::cout << "Testing default CachingQuadratures: " << std::endl;
294         TestCaching<GridPartType,true> testCaching(gridPart, quadOrder, eps);
295         testCaching.runTest();
296       }
297 
298       {
299         std::cout << "Testing DUNE ElementQuadratures: " << std::endl;
300         TestCaching<GridPartType,false, DuneQuadratures> testCaching(gridPart, quadOrder, eps);
301         testCaching.runTest();
302       }
303       {
304         std::cout << "Testing DUNE CachingQuadratures: " << std::endl;
305         TestCaching<GridPartType,true, DuneQuadratures> testCaching(gridPart, quadOrder, eps);
306         testCaching.runTest();
307       }
308       {
309         std::cout << "Testing DUNE CachingLumpiongQuadrature: " << std::endl;
310         TestCaching<GridPartType,true, DuneQuadratures> testCaching(gridPart, quadOrder, eps);
311         testCaching.testLumpingQuadrature();
312       }
313       grid.globalRefine( 1 );
314     }
315     return 0;
316   }
317   catch( const Exception& e )
318   {
319     std :: cerr << e.what() << std :: endl;
320     return 1;
321   }
322 }
323