]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Pass n_q_points_1D to Simplex::QGauss 11103/head
authorPeter Munch <peterrmuench@gmail.com>
Sun, 25 Oct 2020 18:05:00 +0000 (19:05 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Mon, 2 Nov 2020 19:33:17 +0000 (20:33 +0100)
include/deal.II/simplex/quadrature_lib.h
source/simplex/quadrature_lib.cc
tests/simplex/fe_lib_02.cc
tests/simplex/poisson_01.cc
tests/simplex/poisson_02.cc
tests/simplex/q_projection_01.cc
tests/simplex/quadrature_lib_01.cc
tests/simplex/quadrature_lib_01.output
tests/simplex/step-12.cc
tests/simplex/step-18.cc

index 41c02dc02fb2e989fe114ebc94e6e37fe6261548..09e2e97abda841651c06d5b78f7574f427261e12 100644 (file)
@@ -28,12 +28,15 @@ namespace Simplex
   /**
    * Integration rule for simplex entities.
    *
-   * Following number of quadrature points are currently supported for 2D and
-   * 3D:
+   * Users specify a number `n_points_1D` as an indication of what polynomial
+   * degree to be integrated exactly, similarly to the number of points in a
+   * QGauss quadrature object, even though the present quadrature formula is not
+   * a tensor product. The given value is translated for n_points_1D=1,2,3 to
+   * following number of quadrature points for 2D and 3D:
    *   - 2D: 1, 3, 7
    *   - 3D: 1, 4, 10
    *
-   * For 1D, the quadrature rule degenerates to a `QGauss<1>(n_points)`.
+   * For 1D, the quadrature rule degenerates to a `QGauss<1>(n_points_1D)`.
    *
    * @ingroup simplex
    */
@@ -42,9 +45,10 @@ namespace Simplex
   {
   public:
     /**
-     * Constructor taking the number of quadrature points @p n_points.
+     * Constructor taking the number of quadrature points in 1D direction
+     * @p n_points_1D.
      */
-    explicit QGauss(const unsigned int n_points);
+    explicit QGauss(const unsigned int n_points_1D);
   };
 } // namespace Simplex
 
index 683c9ab7cdb150d810c8d63a7568d6885f523c75..bb60716e6c6ed7e67790fbead015319ba4af3a2f 100644 (file)
@@ -29,26 +29,26 @@ DEAL_II_NAMESPACE_OPEN
 namespace Simplex
 {
   template <int dim>
-  QGauss<dim>::QGauss(const unsigned int n_points)
+  QGauss<dim>::QGauss(const unsigned int n_points_1D)
     : QSimplex<dim>(Quadrature<dim>())
   {
     // fill quadrature points and quadrature weights
     if (dim == 1)
       {
-        const dealii::QGauss<dim> quad(n_points);
+        const dealii::QGauss<dim> quad(n_points_1D);
 
         this->quadrature_points = quad.get_points();
         this->weights           = quad.get_weights();
       }
     else if (dim == 2)
       {
-        if (n_points == 1)
+        if (n_points_1D == 1)
           {
             const double p = 1.0 / 3.0;
             this->quadrature_points.emplace_back(p, p);
             this->weights.emplace_back(1.0);
           }
-        else if (n_points == 3)
+        else if (n_points_1D == 2)
           {
             const double Q23 = 2.0 / 3.0;
             const double Q16 = 1.0 / 6.0;
@@ -60,7 +60,7 @@ namespace Simplex
             this->weights.emplace_back(Q16);
             this->weights.emplace_back(Q16);
           }
-        else if (n_points == 7)
+        else if (n_points_1D == 3)
           {
             const double q12 = 0.5;
 
@@ -85,7 +85,7 @@ namespace Simplex
       }
     else if (dim == 3)
       {
-        if (n_points == 1)
+        if (n_points_1D == 1)
           {
             const double Q14 = 1.0 / 4.0;
             const double Q16 = 1.0 / 6.0;
@@ -93,7 +93,7 @@ namespace Simplex
             this->quadrature_points.emplace_back(Q14, Q14, Q14);
             this->weights.emplace_back(Q16);
           }
-        else if (n_points == 4)
+        else if (n_points_1D == 2)
           {
             const double Q124 = 1.0 / 6.0 / 4.0;
 
@@ -108,7 +108,7 @@ namespace Simplex
             this->weights.emplace_back(Q124);
             this->weights.emplace_back(Q124);
           }
-        else if (n_points == 10)
+        else if (n_points_1D == 3)
           {
             const double Q16 = 1.0 / 6.0;
 
index 8e4c141a9434e9c94cffe57d9507d0d40f8dd7d4..feee438b32636c6918d74b49d199dee8fe7d606d 100644 (file)
@@ -49,12 +49,12 @@ main()
 {
   initlog();
 
-  test(Simplex::FE_P<2>(1), Simplex::QGauss<2>(3));
-  test(Simplex::FE_P<2>(2), Simplex::QGauss<2>(7));
-  test(Simplex::FE_P<3>(1), Simplex::QGauss<3>(4));
-  test(Simplex::FE_P<3>(2), Simplex::QGauss<3>(10));
-  test(Simplex::FE_DGP<2>(1), Simplex::QGauss<2>(3));
-  test(Simplex::FE_DGP<2>(2), Simplex::QGauss<2>(7));
-  test(Simplex::FE_DGP<3>(1), Simplex::QGauss<3>(4));
-  test(Simplex::FE_DGP<3>(2), Simplex::QGauss<3>(10));
+  test(Simplex::FE_P<2>(1), Simplex::QGauss<2>(2));
+  test(Simplex::FE_P<2>(2), Simplex::QGauss<2>(3));
+  test(Simplex::FE_P<3>(1), Simplex::QGauss<3>(2));
+  test(Simplex::FE_P<3>(2), Simplex::QGauss<3>(3));
+  test(Simplex::FE_DGP<2>(1), Simplex::QGauss<2>(2));
+  test(Simplex::FE_DGP<2>(2), Simplex::QGauss<2>(3));
+  test(Simplex::FE_DGP<3>(1), Simplex::QGauss<3>(2));
+  test(Simplex::FE_DGP<3>(2), Simplex::QGauss<3>(3));
 }
index e920c6e5397e76b1e11ab5576069dfa1748b81eb..e020592b7df493bbc234e80e73565f74ad3f3f83 100644 (file)
@@ -336,11 +336,9 @@ test_tet(const MPI_Comm &comm, const Parameters<dim> &params)
   // 3) Select components
   Simplex::FE_P<dim> fe(params.degree);
 
-  Simplex::QGauss<dim> quad(dim == 2 ? (params.degree == 1 ? 3 : 7) :
-                                       (params.degree == 1 ? 4 : 10));
+  Simplex::QGauss<dim> quad(params.degree + 1);
 
-  Simplex::QGauss<dim - 1> face_quad(dim == 2 ? (params.degree == 1 ? 2 : 3) :
-                                                (params.degree == 1 ? 3 : 7));
+  Simplex::QGauss<dim - 1> face_quad(params.degree + 1);
 
   Simplex::FE_P<dim> fe_mapping(1);
   MappingFE<dim>     mapping(fe_mapping);
index 494dee8715912ee65fa8c9670fd662f87373e600..c117671b3d73aab9c9ce6e9e4bb0ea9a1ac85ca6 100644 (file)
@@ -189,10 +189,8 @@ public:
       false,
       new Simplex::FE_DGP<dim>(degree),
       new MappingFE<dim>(Simplex::FE_P<dim>(1)),
-      new Simplex::QGauss<dim>(dim == 2 ? (degree == 1 ? 3 : 7) :
-                                          (degree == 1 ? 4 : 10)),
-      new Simplex::QGauss<dim - 1>(dim == 2 ? (degree == 1 ? 2 : 3) :
-                                              (degree == 1 ? 3 : 7)),
+      new Simplex::QGauss<dim>(degree + 1),
+      new Simplex::QGauss<dim - 1>(degree + 1),
       initial_refinement,
       number_refinement);
   }
index 324cdc6080f14d23db4774ee385dbc7b4216a567..7639b3a0c7719afe40a7fa31754970757d31fb70 100644 (file)
@@ -56,8 +56,8 @@ test<2>(const unsigned int n_points)
                                                       face_orientation,
                                                       false,
                                                       false,
-                                                      n_points);
-         q < n_points;
+                                                      quad_ref.size());
+         q < quad_ref.size();
          ++q, ++i)
       {
         deallog << quad.point(i) << " ";
@@ -102,8 +102,8 @@ test<3>(const unsigned int n_points)
                                                       face_orientation,
                                                       face_flip,
                                                       face_rotation,
-                                                      n_points);
-         q < n_points;
+                                                      quad_ref.size());
+         q < quad_ref.size();
          ++q, ++i)
       {
         deallog << quad.point(i) << " ";
@@ -152,12 +152,12 @@ main()
   }
   {
     deallog.push("3d-4");
-    test<3>(3);
+    test<3>(2);
     deallog.pop();
   }
   {
     deallog.push("3d-10");
-    test<3>(7);
+    test<3>(3);
     deallog.pop();
   }
 }
index 85d6b4a0aa60ab569b59c31cfb35f42dfb2ba0df..fa05a42fd0946d1a05693c8c4407c13811cbef17 100644 (file)
@@ -54,7 +54,7 @@ main()
   }
   {
     deallog.push("1d-3");
-    test<1>(2);
+    test<1>(3);
     deallog.pop();
   }
   {
@@ -64,12 +64,12 @@ main()
   }
   {
     deallog.push("2d-3");
-    test<2>(3);
+    test<2>(2);
     deallog.pop();
   }
   {
     deallog.push("2d-7");
-    test<2>(7);
+    test<2>(3);
     deallog.pop();
   }
 
@@ -80,12 +80,12 @@ main()
   }
   {
     deallog.push("3d-4");
-    test<3>(4);
+    test<3>(2);
     deallog.pop();
   }
   {
     deallog.push("3d-10");
-    test<3>(10);
+    test<3>(3);
     deallog.pop();
   }
 }
index 24a8420bbd77a909e7c998c3754d651f2e397d13..94888c477fefdf5c416719a198b2e1355972176b 100644 (file)
@@ -2,8 +2,9 @@
 DEAL:1d-1::0.500000 1.00000 
 DEAL:1d-2::0.211325 0.500000 
 DEAL:1d-2::0.788675 0.500000 
-DEAL:1d-3::0.211325 0.500000 
-DEAL:1d-3::0.788675 0.500000 
+DEAL:1d-3::0.112702 0.277778 
+DEAL:1d-3::0.500000 0.444444 
+DEAL:1d-3::0.887298 0.277778 
 DEAL:2d-1::0.333333 0.333333 1.00000 
 DEAL:2d-3::0.666667 0.166667 0.166667 
 DEAL:2d-3::0.166667 0.666667 0.166667 
index 4c7fe7e2c1df53c823dc78d2c9ee833b10e0908e..29484a3aa275befba6168b0a641b17f5e1686417 100644 (file)
@@ -391,11 +391,9 @@ namespace Step12
 
     QGauss<dim - 1> face_quad(degree + 1);
 #else
-    Simplex::QGauss<dim> quad(dim == 2 ? (degree == 1 ? 3 : 7) :
-                                         (degree == 1 ? 4 : 10));
+    Simplex::QGauss<dim> quad(degree + 1);
 
-    Simplex::QGauss<dim - 1> face_quad(dim == 2 ? (degree == 1 ? 2 : 3) :
-                                                  (degree == 1 ? 3 : 7));
+    Simplex::QGauss<dim - 1> face_quad(degree + 1);
 #endif
 
     ScratchData<dim> scratch_data(mapping, fe, quad, face_quad);
index 2bac10c4d140a17d2f311e3bc18ed8087dc3c9b7..3e2196005ce03e05bc5ebedc4736978cf08582fa 100644 (file)
@@ -419,7 +419,7 @@ namespace Step18
     : triangulation()
     , fe(Simplex::FE_P<dim>(degree), dim)
     , dof_handler(triangulation)
-    , quadrature_formula(Simplex::QGauss<dim>(fe.degree == 1 ? 4 : 10))
+    , quadrature_formula(Simplex::QGauss<dim>(fe.degree + 1))
     , mapping(Simplex::FE_P<dim>(1))
     , present_time(0.0)
     , present_timestep(1.0)

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.