1 #include <config.h>
2 #include <set>
3 
4 #include "quad_test.hh"
5 
6 namespace Dune {
7   namespace Fem {
8 
9   namespace {
doTest(const double & a,const double & b)10     static void doTest( const double& a, const double& b )
11     {
12       if( std::abs( a - b ) > 1e-12 )
13       {
14         assert( false );
15         std::abort();
16       }
17     }
18 
doTest(const bool value)19     static void doTest( const bool value )
20     {
21       if( ! value )
22       {
23         assert( false );
24         std::abort();
25       }
26     }
27   }
28 
run()29   void Quad_Test::run() {
30 
31     // 1d types
32     GeometryType line ( GeometryType::cube, 1);
33     GeometryType cube1 = line;
34     GeometryType simplex1 = line;
35 
36 
37     // 2d types
38     GeometryType quadrilateral ( GeometryType::cube, 2);
39     GeometryType cube2 = quadrilateral;
40     GeometryType triangle ( GeometryType::simplex, 2);
41     GeometryType simplex2 = triangle;
42 
43     // 3d types
44     GeometryType hexahedron ( GeometryType::cube, 3);
45     GeometryType cube3 = hexahedron;
46     GeometryType tetrahedron ( GeometryType::simplex, 3);
47     GeometryType simplex3 = tetrahedron;
48     GeometryType prism ( GeometryType::prism, 3);
49     GeometryType pyramid ( GeometryType::pyramid, 3);
50 
51     // 1D Quadratures
52     Quadrature<double, 1> quadCube1d1(cube1, 1);
53     Quadrature<double, 1> quadLine1d1(line, 1);
54     Quadrature<double, 1> quadSimplex1d1(simplex1, 1);
55     Quadrature<double, 1> quadLine1d3(line, 3);
56     Quadrature<double, 1> quadLine1d10(line, 10);
57 
58     // 2D Quadratures
59     // - Cubes
60     Quadrature<double, 2> quadCube2d1(cube2, 1);
61     Quadrature<double, 2> quadQuad2d1(quadrilateral, 1);
62     Quadrature<double, 2> quadCube2d3(cube2, 3);
63     Quadrature<double, 2> quadCube2d9(cube2, 9);
64     // - Simplices
65     Quadrature<double, 2> quadSimplex2d1(simplex2, 1);
66     Quadrature<double, 2> quadTriangle2d1(triangle, 1);
67     Quadrature<double, 2> quadSimplex2d4(simplex2, 4);
68     Quadrature<double, 2> quadSimplex2d11(simplex2, 11);
69 
70     // 3D Quadratures
71     // - Cubes
72     Quadrature<double, 3> quadCube3d1(cube3, 1);
73     Quadrature<double, 3> quadHexa3d1(hexahedron, 1);
74     Quadrature<double, 3> quadCube3d2(cube3, 2);
75     Quadrature<double, 3> quadCube3d11(cube3, 11);
76     // - Simplices
77     Quadrature<double, 3> quadSimplex3d1(simplex3, 1);
78     Quadrature<double, 3> quadTetra3d1(tetrahedron, 1);
79     Quadrature<double, 3> quadSimplex3d4(simplex3, 4);
80     Quadrature<double, 3> quadSimplex3d7(simplex3, 7);
81     // - Prism and pyramids
82     Quadrature<double, 3> quadPrism3d1(prism, 1);
83     Quadrature<double, 3> quadPyramid3d1(pyramid, 1);
84 
85     /*
86     // FixedOrder quads
87     typedef FieldVector<double, 2> Coord2;
88     typedef FieldVector<double, 3> Coord3;
89     FixedOrderQuad<double, Coord2, 1> fixedSimplex2d1(simplex);
90     FixedOrderQuad<double, Coord2, 4> fixedSimplex2d4(simplex);
91     FixedOrderQuad<double, Coord2, 11> fixedSimplex2d11(simplex);
92 
93     FixedOrderQuad<double, Coord3, 1> fixedSimplex3d1(simplex);
94     FixedOrderQuad<double, Coord3, 4> fixedSimplex3d4(simplex);
95     FixedOrderQuad<double, Coord3, 7> fixedSimplex3d7(simplex);
96 
97     //- Test fixed order for simplex
98     #ifndef HAVE_ALBERTA
99     // if Alberta is present, the new quadratures use the Alberta quadratures
100     // while the old ones stick to the UG quadratures
101     fixedOrderComparisonExec(quadSimplex2d1, fixedSimplex2d1);
102     fixedOrderComparisonExec(quadSimplex2d4, fixedSimplex2d4);
103     fixedOrderComparisonExec(quadSimplex2d11, fixedSimplex2d11);
104     fixedOrderComparisonExec(quadSimplex3d1, fixedSimplex3d1);
105     fixedOrderComparisonExec(quadSimplex3d4, fixedSimplex3d4);
106     fixedOrderComparisonExec(quadSimplex3d7, fixedSimplex3d7);
107     #endif
108     */
109 
110     //- Test weight summation for all
111     weightSummationExec(quadCube1d1);
112     weightSummationExec(quadLine1d1);
113     weightSummationExec(quadSimplex1d1);
114     weightSummationExec(quadLine1d3);
115     weightSummationExec(quadLine1d10);
116     weightSummationExec(quadCube2d1);
117     weightSummationExec(quadQuad2d1);
118     weightSummationExec(quadCube2d3);
119     weightSummationExec(quadCube2d9);
120     weightSummationExec(quadSimplex2d1);
121     weightSummationExec(quadTriangle2d1);
122     weightSummationExec(quadSimplex2d4);
123     weightSummationExec(quadSimplex2d11);
124     weightSummationExec(quadCube3d1);
125     weightSummationExec(quadTetra3d1);
126     weightSummationExec(quadSimplex3d4);
127     weightSummationExec(quadSimplex3d7);
128     weightSummationExec(quadPrism3d1);
129     weightSummationExec(quadPyramid3d1);
130 
131     //- Simple integration test for simplex and cube
132     integrationExec(quadCube1d1);
133     integrationExec(quadLine1d10);
134     integrationExec(quadCube2d3);
135     integrationExec(quadCube2d9);
136     integrationExec(quadSimplex2d4);
137     integrationExec(quadSimplex2d11);
138     integrationExec(quadCube3d11);
139     integrationExec(quadSimplex3d4);
140     integrationExec(quadSimplex3d7);
141 
142     //- Test same geometry type
143     sameGeometryExec(quadCube1d1, quadLine1d1);
144     sameGeometryExec(quadCube2d1, quadQuad2d1);
145     sameGeometryExec(quadCube3d1, quadHexa3d1);
146     sameGeometryExec(quadSimplex1d1, quadLine1d1);
147     sameGeometryExec(quadSimplex2d1, quadTriangle2d1);
148     sameGeometryExec(quadSimplex3d1, quadTetra3d1);
149 
150     //- Test order
151     orderExec(quadCube1d1, 1);
152     orderExec(quadLine1d1, 1);
153     orderExec(quadSimplex1d1, 1);
154     orderExec(quadLine1d3, 3);
155     orderExec(quadLine1d10, 10);
156     orderExec(quadCube2d1, 1);
157     orderExec(quadQuad2d1, 1);
158     orderExec(quadCube2d3, 3);
159     orderExec(quadCube2d9, 9);
160     orderExec(quadSimplex2d1, 1);
161     orderExec(quadTriangle2d1, 1);
162     orderExec(quadSimplex2d4, 4);
163     orderExec(quadSimplex2d11, 11);
164     orderExec(quadCube3d1, 1);
165     orderExec(quadTetra3d1, 1);
166     orderExec(quadSimplex3d4, 4);
167     orderExec(quadSimplex3d7, 7); // This one fails, but it's because of ugquadratures.cc
168 
169     //- Test uniqueness of indices
170     indicesTest();
171 
172     std::set<int> ids;
173     doTest(ids.insert(quadCube1d1.id()).second);
174     doTest(ids.insert(quadLine1d3.id()).second);
175     doTest(ids.insert(quadLine1d10.id()).second);
176     doTest(ids.insert(quadCube2d1.id()).second);
177     doTest(ids.insert(quadCube2d3.id()).second);
178     doTest(ids.insert(quadCube2d9.id()).second);
179     doTest(ids.insert(quadSimplex2d1.id()).second);
180     doTest(ids.insert(quadSimplex2d4.id()).second);
181     doTest(ids.insert(quadSimplex2d11.id()).second);
182     doTest(ids.insert(quadCube3d1.id()).second);
183     doTest(ids.insert(quadCube3d2.id()).second);
184     doTest(ids.insert(quadCube3d11.id()).second);
185     doTest(ids.insert(quadSimplex3d1.id()).second);
186     doTest(ids.insert(quadSimplex3d4.id()).second);
187     doTest(ids.insert(quadSimplex3d7.id()).second);
188     doTest(ids.insert(quadPrism3d1.id()).second);
189     doTest(ids.insert(quadPyramid3d1.id()).second);
190 
191   }
192 
193   template <class Quad, class Fixed>
fixedOrderComparisonExec(Quad & quad,Fixed & fixed)194   void Quad_Test::fixedOrderComparisonExec(Quad& quad, Fixed& fixed)
195   {
196     for (int i = 0; i < quad.nop(); ++i) {
197       for (int j = 0; j < Quad::dimension; ++j) {
198         doTest(quad.point(i)[j], fixed.point(i)[j]);
199       }
200       doTest(quad.weight(i), fixed.weight(i));
201     }
202   }
203 
204   template <class Quad>
weightSummationExec(Quad & quad)205   void Quad_Test::weightSummationExec(Quad& quad)
206   {
207     double sum = 0.;
208     for (int i = 0; i < quad.nop(); ++i) {
209       sum += quad.weight(i);
210     }
211 
212     const GeometryType geom = quad.geometry();
213     if ( geom.isCube())
214     {
215       doTest(sum, 1.0);
216     }
217     else if ( geom.isSimplex())
218     {
219       switch(Quad::dimension) {
220       case 1:
221         doTest(sum, 1.0);
222         break;
223       case 2:
224         doTest(sum, 0.5);
225         break;
226       case 3:
227         doTestTol(sum, 0.16666666666666666667, 1e-04);
228         break;
229       }
230     }
231     else if( geom.isPrism() )
232     {
233       doTest(sum, 0.5);
234     }
235     else if ( geom.isPyramid() )
236     {
237       doTestTol(sum, 0.3333333333333333, 1e-04);
238     }
239     else
240     {
241       DUNE_THROW(NotImplemented, "What geometry type ey?");
242     }
243   }
244 
245   template <class Quad>
integrationExec(Quad & quad)246   void Quad_Test::integrationExec(Quad& quad)
247   {
248     Func f;
249     double result = 0.;
250     for (int i = 0; i < quad.nop(); ++i) {
251       result += f(quad.point(i))*quad.weight(i);
252     }
253 
254     switch (Quad::dimension) {
255     case 1:
256       doTest(result, -0.25);
257       break;
258     case 2:
259       if (quad.geometry().isCube() ) {
260         doTest(result, -0.2708333333);
261       }
262       else {
263         doTest(result, -0.0333333333);
264       }
265       break;
266     case 3:
267       if (quad.geometry().isCube()) {
268         doTest(result, -0.03854166666);
269       }
270       else {
271         doTestTol(result, -0.00551388888, 1e-05);
272       }
273       break;
274     default:
275       _fail("should not get here");
276     }
277   }
278 
279   template <class Quad>
orderExec(Quad & quad,int order)280   void Quad_Test::orderExec(Quad& quad, int order)
281   {
282     doTest(quad.order() >= order);
283   }
284 
285   template <class Quad>
sameGeometryExec(Quad & quad1,Quad & quad2)286   void Quad_Test::sameGeometryExec(Quad& quad1, Quad& quad2)
287   {
288     doTest(quad1.id() == quad2.id());
289   }
290 
indicesTest()291   void Quad_Test::indicesTest()
292   {
293     const GeometryType line ( GeometryType::cube, 1 );
294     size_t id;
295     {
296       Quadrature<double, 1> quadTemp(line, 5);
297       id = quadTemp.id();
298     }
299     Quadrature<double, 1> quadTemp(line, 5);
300     doTest(quadTemp.id() == id);
301   }
302   } // end namespace Fem
303 } // end namespace Dune
304