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