]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Enable FEFaceValues for wedges and pyramids 11361/head
authorPeter Munch <peterrmuench@gmail.com>
Sat, 12 Dec 2020 14:22:33 +0000 (15:22 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Sat, 19 Dec 2020 16:45:01 +0000 (17:45 +0100)
source/base/qprojector.cc
source/fe/mapping_fe.cc
tests/simplex/poisson_01.cc
tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output
tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output

index 5515b5465207b103b09221f8e548040818d13f84..b0cd657d780b6a9d0512ce5b68a5068eb9eb0f8c 100644 (file)
@@ -16,6 +16,7 @@
 #include <deal.II/base/derivative_form.h>
 #include <deal.II/base/geometry_info.h>
 #include <deal.II/base/qprojector.h>
+#include <deal.II/base/tensor_product_polynomials.h>
 
 #include <deal.II/simplex/polynomials.h>
 
@@ -642,77 +643,88 @@ QProjector<3>::project_to_all_faces(
   const ReferenceCell::Type reference_cell_type,
   const hp::QCollection<2> &quadrature)
 {
-  if (reference_cell_type == ReferenceCell::Type::Tet)
-    {
-      // reference faces (defined by its support points and its area)
-      // note: the area is later not used as a scaling factor but recomputed
-      const std::array<std::pair<std::array<Point<3>, 3>, double>, 4> faces = {
-        {{{{Point<3>(0.0, 0.0, 0.0),
-            Point<3>(1.0, 0.0, 0.0),
-            Point<3>(0.0, 1.0, 0.0)}},
-          0.5},
-         {{{Point<3>(1.0, 0.0, 0.0),
-            Point<3>(0.0, 0.0, 0.0),
-            Point<3>(0.0, 0.0, 1.0)}},
-          0.5},
-         {{{Point<3>(0.0, 0.0, 0.0),
-            Point<3>(0.0, 1.0, 0.0),
-            Point<3>(0.0, 0.0, 1.0)}},
-          0.5},
-         {{{Point<3>(0.0, 1.0, 0.0),
-            Point<3>(1.0, 0.0, 0.0),
-            Point<3>(0.0, 0.0, 1.0)}},
-          0.5 * sqrt(3.0) /*equilateral triangle*/}}};
+  const auto support_points_tri =
+    [](const auto &face, const auto &orientation) -> std::vector<Point<3>> {
+    // determine support point of the current line with the correct
+    // orientation
+    switch (orientation)
+      {
+        case 1:
+          return {{face.first[0], face.first[1], face.first[2]}};
+        case 3:
+          return {{face.first[1], face.first[0], face.first[2]}};
+        case 5:
+          return {{face.first[2], face.first[0], face.first[1]}};
+        case 0:
+          return {{face.first[0], face.first[2], face.first[1]}};
+        case 2:
+          return {{face.first[1], face.first[2], face.first[0]}};
+        case 4:
+          return {{face.first[2], face.first[1], face.first[0]}};
+        default:
+          Assert(false, ExcNotImplemented());
+          return {{}};
+      }
+  };
 
-      // linear polynomial to map the reference quadrature points correctly
-      // on faces
-      const Simplex::ScalarPolynomial<2> poly(1);
+  const auto support_points_quad =
+    [](const auto &face, const auto &orientation) -> std::vector<Point<3>> {
+    switch (orientation)
+      {
+        case 1:
+          return {{face.first[0], face.first[1], face.first[2], face.first[3]}};
+        case 3:
+          return {{face.first[1], face.first[3], face.first[0], face.first[2]}};
+        case 5:
+          return {{face.first[3], face.first[2], face.first[1], face.first[0]}};
+        case 7:
+          return {{face.first[2], face.first[0], face.first[3], face.first[1]}};
+        case 0:
+          return {{face.first[0], face.first[2], face.first[1], face.first[3]}};
+        case 2:
+          return {{face.first[2], face.first[3], face.first[0], face.first[1]}};
+        case 4:
+          return {{face.first[3], face.first[1], face.first[2], face.first[0]}};
+        case 6:
+          return {{face.first[1], face.first[0], face.first[3], face.first[2]}};
+        default:
+          Assert(false, ExcNotImplemented());
+          return {{}};
+      }
+  };
 
-      // new (projected) quadrature points and weights
-      std::vector<Point<3>> points;
-      std::vector<double>   weights;
+  const auto process = [&](const auto faces) {
+    // new (projected) quadrature points and weights
+    std::vector<Point<3>> points;
+    std::vector<double>   weights;
+
+    const Simplex::ScalarPolynomial<2> poly_tri(1);
+    const TensorProductPolynomials<2>  poly_quad(
+      Polynomials::generate_complete_Lagrange_basis(
+        {Point<1>(0.0), Point<1>(1.0)}));
+
+    // loop over all faces (triangles) ...
+    for (unsigned int face_no = 0; face_no < faces.size(); ++face_no)
+      {
+        // linear polynomial to map the reference quadrature points correctly
+        // on faces
+        const unsigned int n_shape_functions = faces[face_no].first.size();
+
+        const auto &poly =
+          n_shape_functions == 3 ?
+            static_cast<const ScalarPolynomialsBase<2> &>(poly_tri) :
+            static_cast<const ScalarPolynomialsBase<2> &>(poly_quad);
 
-      // loop over all faces (triangles) ...
-      for (unsigned int face_no = 0; face_no < faces.size(); ++face_no)
         // ... and over all possible orientations
-        for (unsigned int orientation = 0; orientation < 6; ++orientation)
+        for (unsigned int orientation = 0;
+             orientation < (n_shape_functions * 2);
+             ++orientation)
           {
             const auto &face = faces[face_no];
 
-            std::array<Point<3>, 3> support_points;
-
-            // determine support point of the current line with the correct
-            // orientation
-            switch (orientation)
-              {
-                case 1:
-                  support_points = {
-                    {face.first[0], face.first[1], face.first[2]}};
-                  break;
-                case 3:
-                  support_points = {
-                    {face.first[1], face.first[0], face.first[2]}};
-                  break;
-                case 5:
-                  support_points = {
-                    {face.first[2], face.first[0], face.first[1]}};
-                  break;
-                case 0:
-                  support_points = {
-                    {face.first[0], face.first[2], face.first[1]}};
-                  break;
-                case 2:
-                  support_points = {
-                    {face.first[1], face.first[2], face.first[0]}};
-                  break;
-                case 4:
-                  support_points = {
-                    {face.first[2], face.first[1], face.first[0]}};
-                  break;
-                default:
-                  Assert(false, ExcNotImplemented());
-              }
-
+            const auto support_points =
+              n_shape_functions == 3 ? support_points_tri(face, orientation) :
+                                       support_points_quad(face, orientation);
 
             // the quadrature rule to be projected ...
             const auto &sub_quadrature_points =
@@ -726,7 +738,7 @@ QProjector<3>::project_to_all_faces(
                 Point<3> mapped_point;
 
                 // map reference quadrature point
-                for (unsigned int i = 0; i < 3; ++i)
+                for (unsigned int i = 0; i < n_shape_functions; ++i)
                   mapped_point +=
                     support_points[i] *
                     poly.compute_value(i, sub_quadrature_points[j]);
@@ -735,18 +747,16 @@ QProjector<3>::project_to_all_faces(
 
                 // scale quadrature weight
                 const double scaling = [&]() {
-                  const auto &supp_pts = support_points;
-
-                  const unsigned int n_shape_functions = 3;
-                  const unsigned int dim_              = 2;
-                  const unsigned int spacedim          = 3;
+                  const auto &       supp_pts = support_points;
+                  const unsigned int dim_     = 2;
+                  const unsigned int spacedim = 3;
 
                   double result[spacedim][dim_];
 
                   std::vector<Tensor<1, dim_>> shape_derivatives(
                     n_shape_functions);
 
-                  for (unsigned int i = 0; i < 3; ++i)
+                  for (unsigned int i = 0; i < n_shape_functions; ++i)
                     shape_derivatives[i] =
                       poly.compute_1st_derivative(i, sub_quadrature_points[j]);
 
@@ -782,9 +792,91 @@ QProjector<3>::project_to_all_faces(
                 weights.push_back(sub_quadrature_weights[j] * scaling);
               }
           }
+      }
 
-      // construct new quadrature rule
-      return {points, weights};
+    // construct new quadrature rule
+    return Quadrature<3>(points, weights);
+  };
+
+  if (reference_cell_type == ReferenceCell::Type::Tet)
+    {
+      // reference faces (defined by its support points and its area)
+      // note: the area is later not used as a scaling factor but recomputed
+      const std::vector<std::pair<std::vector<Point<3>>, double>> faces = {
+        {{{{Point<3>(0.0, 0.0, 0.0),
+            Point<3>(1.0, 0.0, 0.0),
+            Point<3>(0.0, 1.0, 0.0)}},
+          0.5},
+         {{{Point<3>(1.0, 0.0, 0.0),
+            Point<3>(0.0, 0.0, 0.0),
+            Point<3>(0.0, 0.0, 1.0)}},
+          0.5},
+         {{{Point<3>(0.0, 0.0, 0.0),
+            Point<3>(0.0, 1.0, 0.0),
+            Point<3>(0.0, 0.0, 1.0)}},
+          0.5},
+         {{{Point<3>(0.0, 1.0, 0.0),
+            Point<3>(1.0, 0.0, 0.0),
+            Point<3>(0.0, 0.0, 1.0)}},
+          0.5 * sqrt(3.0) /*equilateral triangle*/}}};
+
+      return process(faces);
+    }
+  else if (reference_cell_type == ReferenceCell::Type::Wedge)
+    {
+      const std::vector<std::pair<std::vector<Point<3>>, double>> faces = {
+        {{{{Point<3>(0.0, 1.0, 0.0),
+            Point<3>(1.0, 0.0, 0.0),
+            Point<3>(0.0, 0.0, 0.0)}},
+          0.5},
+         {{{Point<3>(1.0, 0.0, 1.0),
+            Point<3>(0.0, 1.0, 1.0),
+            Point<3>(0.0, 0.0, 1.0)}},
+          0.5},
+         {{{Point<3>(1.0, 0.0, 0.0),
+            Point<3>(0.0, 1.0, 0.0),
+            Point<3>(1.0, 0.0, 1.0),
+            Point<3>(0.0, 1.0, 1.0)}},
+          std::sqrt(2.0)},
+         {{{Point<3>(0.0, 1.0, 0.0),
+            Point<3>(0.0, 0.0, 0.0),
+            Point<3>(0.0, 1.0, 1.0),
+            Point<3>(0.0, 0.0, 1.0)}},
+          1.0},
+         {{{Point<3>(0.0, 0.0, 0.0),
+            Point<3>(1.0, 0.0, 0.0),
+            Point<3>(0.0, 0.0, 1.0),
+            Point<3>(1.0, 0.0, 1.0)}},
+          1.0}}};
+
+      return process(faces);
+    }
+  else if (reference_cell_type == ReferenceCell::Type::Pyramid)
+    {
+      const std::vector<std::pair<std::vector<Point<3>>, double>> faces = {
+        {{{{Point<3>(-1.0, -1.0, 0.0),
+            Point<3>(+1.0, -1.0, 0.0),
+            Point<3>(-1.0, +1.0, 0.0),
+            Point<3>(+1.0, +1.0, 0.0)}},
+          4.0},
+         {{{Point<3>(-1.0, -1.0, 0.0),
+            Point<3>(-1.0, +1.0, 0.0),
+            Point<3>(+0.0, +0.0, 1.0)}},
+          std::sqrt(2.0)},
+         {{{Point<3>(+1.0, +1.0, 0.0),
+            Point<3>(+1.0, -1.0, 0.0),
+            Point<3>(+0.0, +0.0, 1.0)}},
+          std::sqrt(2.0)},
+         {{{Point<3>(+1.0, -1.0, 0.0),
+            Point<3>(-1.0, -1.0, 0.0),
+            Point<3>(+0.0, +0.0, 1.0)}},
+          std::sqrt(2.0)},
+         {{{Point<3>(-1.0, +1.0, 0.0),
+            Point<3>(+1.0, +1.0, 0.0),
+            Point<3>(+0.0, +0.0, 1.0)}},
+          std::sqrt(2.0)}}};
+
+      return process(faces);
     }
 
 
@@ -1331,18 +1423,36 @@ QProjector<dim>::DataSetDescriptor::face(
   const hp::QCollection<dim - 1> &quadrature)
 {
   if (reference_cell_type == ReferenceCell::Type::Tri ||
-      reference_cell_type == ReferenceCell::Type::Tet)
+      reference_cell_type == ReferenceCell::Type::Tet ||
+      reference_cell_type == ReferenceCell::Type::Wedge ||
+      reference_cell_type == ReferenceCell::Type::Pyramid)
     {
       unsigned int offset = 0;
 
+      static const unsigned int X = numbers::invalid_unsigned_int;
+      static const std::array<unsigned int, 5> scale_tri   = {{2, 2, 2, X, X}};
+      static const std::array<unsigned int, 5> scale_tet   = {{6, 6, 6, 6, X}};
+      static const std::array<unsigned int, 5> scale_wedge = {{6, 6, 8, 8, 8}};
+      static const std::array<unsigned int, 5> scale_pyramid = {
+        {8, 6, 6, 6, 6}};
+
+      const auto &scale =
+        (reference_cell_type == ReferenceCell::Type::Tri) ?
+          scale_tri :
+          ((reference_cell_type == ReferenceCell::Type::Tet) ?
+             scale_tet :
+             ((reference_cell_type == ReferenceCell::Type::Wedge) ?
+                scale_wedge :
+                scale_pyramid));
+
       if (quadrature.size() == 1)
-        offset = quadrature[0].size() * face_no;
+        offset = scale[0] * quadrature[0].size() * face_no;
       else
         for (unsigned int i = 0; i < face_no; ++i)
-          offset += quadrature[i].size();
+          offset += scale[i] * quadrature[i].size();
 
       if (dim == 2)
-        return {2 * offset +
+        return {offset +
                 face_orientation *
                   quadrature[quadrature.size() == 1 ? 0 : face_no].size()};
       else if (dim == 3)
@@ -1350,7 +1460,7 @@ QProjector<dim>::DataSetDescriptor::face(
           const unsigned int orientation =
             (face_flip * 2 + face_rotation) * 2 + face_orientation;
 
-          return {6 * offset +
+          return {offset +
                   orientation *
                     quadrature[quadrature.size() == 1 ? 0 : face_no].size()};
         }
index a5c428ea4ceebce61fcb5d8d9a00e9eb84ce5849..169f944bf561a0cc6e018656f64299311451647a 100644 (file)
@@ -246,6 +246,160 @@ MappingFE<dim, spacedim>::InternalData::initialize_face(
           for (unsigned int i = 0; i < n_original_q_points; i++)
             unit_tangentials[7].emplace_back(t1);
         }
+      else if (this->fe.reference_cell_type() == ReferenceCell::Type::Wedge)
+        {
+          Tensor<1, dim> t1;
+          constexpr int  d0 = 0;
+          constexpr int  d1 = 1 % dim;
+          constexpr int  d2 = 2 % dim;
+
+          // face 0
+          t1[d0] = 0;
+          t1[d1] = 1;
+          t1[d2] = 0;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[0].emplace_back(t1);
+
+          // face 0
+          t1[d0] = 1;
+          t1[d1] = 0;
+          t1[d2] = 0;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[5].emplace_back(t1);
+
+          // face 1
+          t1[d0] = 1;
+          t1[d1] = 0;
+          t1[d2] = 0;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[1].emplace_back(t1);
+
+          // face 1
+          t1[d0] = 0;
+          t1[d1] = 1;
+          t1[d2] = 0;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[6].emplace_back(t1);
+
+          // face 2
+          t1[d0] = -1 / std::sqrt(2.0);
+          t1[d1] = +1 / std::sqrt(2.0);
+          t1[d2] = 0;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[2].emplace_back(t1);
+
+          // face 2
+          t1[d0] = 0;
+          t1[d1] = 0;
+          t1[d2] = 1;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[7].emplace_back(t1);
+
+          // face 3
+          t1[d0] = +0;
+          t1[d1] = +0;
+          t1[d2] = +1;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[3].emplace_back(t1);
+
+          // face 3
+          t1[d0] = +0;
+          t1[d1] = +1;
+          t1[d2] = +0;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[8].emplace_back(t1);
+
+          // face 4
+          t1[d0] = 1;
+          t1[d1] = 0;
+          t1[d2] = 0;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[4].emplace_back(t1);
+
+          // face 4
+          t1[d0] = 0;
+          t1[d1] = 0;
+          t1[d2] = 1;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[9].emplace_back(t1);
+        }
+      else if (this->fe.reference_cell_type() == ReferenceCell::Type::Pyramid)
+        {
+          Tensor<1, dim> t1;
+          constexpr int  d0 = 0;
+          constexpr int  d1 = 1 % dim;
+          constexpr int  d2 = 2 % dim;
+
+          // face 0
+          t1[d0] = 0;
+          t1[d1] = 1;
+          t1[d2] = 0;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[0].emplace_back(t1);
+
+          // face 0
+          t1[d0] = 1;
+          t1[d1] = 0;
+          t1[d2] = 0;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[5].emplace_back(t1);
+
+          // face 1
+          t1[d0] = 1.0 / sqrt(2.0);
+          t1[d1] = 0;
+          t1[d2] = 1.0 / sqrt(2.0);
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[1].emplace_back(t1);
+
+          // face 1
+          t1[d0] = 0;
+          t1[d1] = 1;
+          t1[d2] = 0;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[6].emplace_back(t1);
+
+          // face 2
+          t1[d0] = 1.0 / sqrt(2.0);
+          t1[d1] = 0;
+          t1[d2] = -1.0 / sqrt(2.0);
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[2].emplace_back(t1);
+
+          // face 2
+          t1[d0] = 0;
+          t1[d1] = 1;
+          t1[d2] = 0;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[7].emplace_back(t1);
+
+          // face 3
+          t1[d0] = 1;
+          t1[d1] = 0;
+          t1[d2] = 0;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[3].emplace_back(t1);
+
+          // face 3
+          t1[d0] = 0;
+          t1[d1] = 1.0 / sqrt(2.0);
+          t1[d2] = 1.0 / sqrt(2.0);
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[8].emplace_back(t1);
+
+          // face 4
+          t1[d0] = 1;
+          t1[d1] = 0;
+          t1[d2] = 0;
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[4].emplace_back(t1);
+
+          // face 4
+          t1[d0] = 0;
+          t1[d1] = +1.0 / sqrt(2.0);
+          t1[d2] = -1.0 / sqrt(2.0);
+          for (unsigned int i = 0; i < n_original_q_points; i++)
+            unit_tangentials[9].emplace_back(t1);
+        }
       else
         {
           Assert(false, ExcNotImplemented());
index 44e9f0580049d4f5cded51d3bd6e99f47c53a889..56a2020e412966ec6e7c2214b074ad292392d12e 100644 (file)
@@ -38,6 +38,8 @@
 #include <deal.II/grid/grid_tools.h>
 #include <deal.II/grid/tria.h>
 
+#include <deal.II/hp/q_collection.h>
+
 #include <deal.II/lac/affine_constraints.h>
 #include <deal.II/lac/dynamic_sparsity_pattern.h>
 #include <deal.II/lac/full_matrix.h>
@@ -96,7 +98,7 @@ void
 test(const Triangulation<dim, spacedim> &tria,
      const FiniteElement<dim, spacedim> &fe,
      const Quadrature<dim> &             quad,
-     const Quadrature<dim - 1> &         face_quad,
+     const hp::QCollection<dim - 1> &    face_quad,
      const Mapping<dim, spacedim> &      mapping,
      const double                        r_boundary,
      const bool                          do_use_fe_face_values = true)
@@ -230,7 +232,7 @@ test(const Triangulation<dim, spacedim> &tria,
           if (face->at_boundary() && (face->boundary_id() == 1))
             {
               fe_face_values->reinit(cell, face);
-              for (unsigned int q = 0; q < face_quad.size(); ++q)
+              for (const auto q : fe_face_values->quadrature_point_indices())
                 for (unsigned int i = 0; i < dofs_per_cell; ++i)
                   cell_rhs(i) +=
                     (1.0 *                               // 1.0
@@ -289,7 +291,9 @@ test(const Triangulation<dim, spacedim> &tria,
 
       data_out.build_patches(mapping);
 
-      std::ofstream output("result.vtk");
+      static unsigned int counter = 0;
+
+      std::ofstream output("result" + std::to_string(counter++) + ".vtk");
       data_out.write_vtk(output);
     }
 }
@@ -372,7 +376,8 @@ test_tet(const MPI_Comm &comm, const Parameters<dim> &params)
 
   Simplex::QGauss<dim> quad(params.degree + 1);
 
-  Simplex::QGauss<dim - 1> face_quad(params.degree + 1);
+  hp::QCollection<dim - 1> face_quad{
+    Simplex::QGauss<dim - 1>(params.degree + 1)};
 
   Simplex::FE_P<dim> fe_mapping(1);
   MappingFE<dim>     mapping(fe_mapping);
@@ -415,12 +420,12 @@ test_hex(const MPI_Comm &comm, const Parameters<dim> &params)
 
   QGauss<dim> quad(params.degree + 1);
 
-  QGauss<dim - 1> quad_face(params.degree + 1);
+  hp::QCollection<dim - 1> face_quad{QGauss<dim - 1>(params.degree + 1)};
 
   MappingQ<dim, spacedim> mapping(1);
 
   // 4) Perform test (independent of mesh type)
-  test(tria, fe, quad, quad_face, mapping, params.p2[0]);
+  test(tria, fe, quad, face_quad, mapping, params.p2[0]);
 }
 
 template <int dim, int spacedim = dim>
@@ -501,13 +506,18 @@ test_wedge(const MPI_Comm &comm, const Parameters<dim> &params)
 
   Simplex::QGaussWedge<dim> quad(params.degree + 1);
 
-  Quadrature<dim - 1> face_quad; // not needed
+  hp::QCollection<dim - 1> face_quad{
+    Simplex::QGauss<dim - 1>(params.degree + 1),
+    Simplex::QGauss<dim - 1>(params.degree + 1),
+    QGauss<dim - 1>(params.degree + 1),
+    QGauss<dim - 1>(params.degree + 1),
+    QGauss<dim - 1>(params.degree + 1)};
 
   Simplex::FE_WedgeP<dim> fe_mapping(1);
   MappingFE<dim>          mapping(fe_mapping);
 
   // 4) Perform test (independent of mesh type)
-  test(*tria, fe, quad, face_quad, mapping, params.p2[0], false);
+  test(*tria, fe, quad, face_quad, mapping, params.p2[0], true);
 }
 
 template <int dim, int spacedim = dim>
@@ -588,13 +598,18 @@ test_pyramid(const MPI_Comm &comm, const Parameters<dim> &params)
 
   Simplex::QGaussPyramid<dim> quad(params.degree + 1);
 
-  Quadrature<dim - 1> face_quad; // not needed
+  hp::QCollection<dim - 1> face_quad{
+    QGauss<dim - 1>(params.degree + 1),
+    Simplex::QGauss<dim - 1>(params.degree + 1),
+    Simplex::QGauss<dim - 1>(params.degree + 1),
+    Simplex::QGauss<dim - 1>(params.degree + 1),
+    Simplex::QGauss<dim - 1>(params.degree + 1)};
 
   Simplex::FE_PyramidP<dim> fe_mapping(1);
   MappingFE<dim>            mapping(fe_mapping);
 
   // 4) Perform test (independent of mesh type)
-  test(*tria, fe, quad, face_quad, mapping, params.p2[0], false);
+  test(*tria, fe, quad, face_quad, mapping, params.p2[0], true);
 }
 
 int
index f239fef1ac4483ae415f7cc8238837ee4b5a4e4e..1df99dcc927b94a9a568934e453e6a3cd1307f25 100644 (file)
@@ -25,13 +25,13 @@ DEAL::   with 134 CG iterations needed to obtain convergence
 DEAL::
 DEAL::Solve problem on WEDGE mesh:
 DEAL::   on parallel::fullydistributed::Triangulation
-DEAL:cg::Starting value 0.0132550
-DEAL:cg::Convergence step 186 value 8.72381e-13
-DEAL::   with 186 CG iterations needed to obtain convergence
+DEAL:cg::Starting value 0.0561953
+DEAL:cg::Convergence step 197 value 8.61398e-13
+DEAL::   with 197 CG iterations needed to obtain convergence
 DEAL::
 DEAL::Solve problem on PYRAMID mesh:
 DEAL::   on parallel::fullydistributed::Triangulation
-DEAL:cg::Starting value 0.0230669
-DEAL:cg::Convergence step 79 value 7.70260e-13
-DEAL::   with 79 CG iterations needed to obtain convergence
+DEAL:cg::Starting value 0.101163
+DEAL:cg::Convergence step 83 value 8.81456e-13
+DEAL::   with 83 CG iterations needed to obtain convergence
 DEAL::
index 01f8dc29c4cf0541a69069e856f7c42f0fccc385..6c0ab9a4a121705efac77969cc460450dd57ccf2 100644 (file)
@@ -25,13 +25,13 @@ DEAL::   with 134 CG iterations needed to obtain convergence
 DEAL::
 DEAL::Solve problem on WEDGE mesh:
 DEAL::   on parallel::fullydistributed::Triangulation
-DEAL:cg::Starting value 0.0132550
-DEAL:cg::Convergence step 186 value 8.72414e-13
-DEAL::   with 186 CG iterations needed to obtain convergence
+DEAL:cg::Starting value 0.0561953
+DEAL:cg::Convergence step 197 value 8.60684e-13
+DEAL::   with 197 CG iterations needed to obtain convergence
 DEAL::
 DEAL::Solve problem on PYRAMID mesh:
 DEAL::   on parallel::fullydistributed::Triangulation
-DEAL:cg::Starting value 0.0230669
-DEAL:cg::Convergence step 79 value 7.71392e-13
-DEAL::   with 79 CG iterations needed to obtain convergence
+DEAL:cg::Starting value 0.101163
+DEAL:cg::Convergence step 83 value 9.01991e-13
+DEAL::   with 83 CG iterations needed to obtain convergence
 DEAL::

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.