From: Peter Munch <peterrmuench@gmail.com>
Date: Wed, 2 Dec 2020 22:16:46 +0000 (+0100)
Subject: MatrixFree: enable wedges (excluding faces integrals)
X-Git-Tag: v9.3.0-rc1~802^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11298%2Fhead;p=dealii.git

MatrixFree: enable wedges (excluding faces integrals)
---

diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h
index 444d5c1b8c..c2d822fff7 100644
--- a/include/deal.II/matrix_free/mapping_info.templates.h
+++ b/include/deal.II/matrix_free/mapping_info.templates.h
@@ -405,9 +405,17 @@ namespace internal
               if (flag == false)
                 {
 #ifdef DEAL_II_WITH_SIMPLEX_SUPPORT
-                  const auto quad_face = get_face_quadrature(quad[my_q][hpq]);
-                  face_data[my_q].descriptor[hpq].initialize(quad_face,
-                                                             update_default);
+                  try
+                    {
+                      const auto quad_face =
+                        get_face_quadrature(quad[my_q][hpq]);
+                      face_data[my_q].descriptor[hpq].initialize(
+                        quad_face, update_default);
+                    }
+                  catch (...)
+                    {
+                      // TODO: nothing to do for now for wedges and pyramids.
+                    }
 #else
                   Assert(false, ExcNotImplemented());
 #endif
@@ -431,9 +439,17 @@ namespace internal
               if (flag == false)
                 {
 #ifdef DEAL_II_WITH_SIMPLEX_SUPPORT
-                  const auto quad_face = get_face_quadrature(quad[my_q][hpq]);
-                  face_data_by_cells[my_q].descriptor[hpq].initialize(
-                    quad_face, update_default);
+                  try
+                    {
+                      const auto quad_face =
+                        get_face_quadrature(quad[my_q][hpq]);
+                      face_data_by_cells[my_q].descriptor[hpq].initialize(
+                        quad_face, update_default);
+                    }
+                  catch (...)
+                    {
+                      // TODO: nothing to do for now for wedges and pyramids.
+                    }
 #else
                   Assert(false, ExcNotImplemented());
 #endif
diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index 636127c5ba..f087691d63 100644
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -92,7 +92,8 @@ namespace internal
       const auto fe_poly = dynamic_cast<const FE_Poly<dim, dim> *>(fe);
 
       if (dynamic_cast<const Simplex::FE_P<dim, dim> *>(fe) != nullptr ||
-          dynamic_cast<const Simplex::FE_DGP<dim, dim> *>(fe) != nullptr)
+          dynamic_cast<const Simplex::FE_DGP<dim, dim> *>(fe) != nullptr ||
+          dynamic_cast<const Simplex::FE_WedgeP<dim, dim> *>(fe) != nullptr)
         {
           scalar_lexicographic.resize(fe->n_dofs_per_cell());
           for (unsigned int i = 0; i < scalar_lexicographic.size(); ++i)
@@ -204,7 +205,8 @@ namespace internal
               // Simplices are a special case since the polynomial family is not
               // indicative of their support
               if (dynamic_cast<const Simplex::FE_P<dim> *>(fe_poly_ptr) ||
-                  dynamic_cast<const Simplex::FE_DGP<dim> *>(fe_poly_ptr))
+                  dynamic_cast<const Simplex::FE_DGP<dim> *>(fe_poly_ptr) ||
+                  dynamic_cast<const Simplex::FE_WedgeP<dim> *>(fe_poly_ptr))
                 return true;
 #endif
 
@@ -243,6 +245,8 @@ namespace internal
           dynamic_cast<const Simplex::FE_P<dim> *>(
             &fe_in.base_element(base_element_number)) ||
           dynamic_cast<const Simplex::FE_DGP<dim> *>(
+            &fe_in.base_element(base_element_number)) ||
+          dynamic_cast<const Simplex::FE_WedgeP<dim> *>(
             &fe_in.base_element(base_element_number)))
         {
           // specialization for arbitrary finite elements and quadrature rules
@@ -306,10 +310,10 @@ namespace internal
                                   q] = grad[d];
               }
 
-          const auto reference_cell_type = ReferenceCell::get_simplex(dim);
-
           try
             {
+              const auto reference_cell_type = ReferenceCell::get_simplex(dim);
+
               const auto quad_face  = get_face_quadrature(quad);
               this->n_q_points_face = quad_face.size();