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