Fix MatrixFree + FE_Q_iso_Q1(2)

include/deal.II/matrix_free/evaluation_kernels.h

@@ -2303,7 +2303,8 @@ namespace internal
                               values_dofs,
                               fe_eval);
         }
-      else if (fe_degree >= 0 && element_type <= ElementType::tensor_symmetric)
+      else if (fe_degree >= 0 &&
+               element_type <= ElementType::tensor_symmetric_no_collocation)
         {
           FEEvaluationImpl<ElementType::tensor_symmetric,
                            dim,
@@ -2438,7 +2439,8 @@ namespace internal
                                fe_eval,
                                sum_into_values_array);
         }
-      else if (fe_degree >= 0 && element_type <= ElementType::tensor_symmetric)
+      else if (fe_degree >= 0 &&
+               element_type <= ElementType::tensor_symmetric_no_collocation)
         {
           FEEvaluationImpl<ElementType::tensor_symmetric,
                            dim,
@@ -5174,7 +5176,7 @@ namespace internal
     {
       Assert(fe_degree > -1, ExcInternalError());
       Assert(fe_eval.get_shape_info().element_type <=
-               MatrixFreeFunctions::tensor_symmetric,
+               MatrixFreeFunctions::tensor_symmetric_no_collocation,
              ExcInternalError());
 
       const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
@@ -5293,7 +5295,7 @@ namespace internal
           (evaluation_flag & EvaluationFlags::hessians) ||
           vector_ptr == nullptr ||
           shape_info.data.front().element_type >
-            MatrixFreeFunctions::tensor_symmetric ||
+            MatrixFreeFunctions::tensor_symmetric_no_collocation ||
           storage <
             MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous)
         return false;
@@ -5394,7 +5396,7 @@ namespace internal
     {
       Assert(fe_degree > -1, ExcInternalError());
       Assert(fe_eval.get_shape_info().element_type <=
-               MatrixFreeFunctions::tensor_symmetric,
+               MatrixFreeFunctions::tensor_symmetric_no_collocation,
              ExcInternalError());
 
       const unsigned int dofs_per_face = Utilities::pow(fe_degree + 1, dim - 1);
@@ -5601,7 +5603,7 @@ namespace internal
 
       Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
       Assert(fe_eval.get_shape_info().element_type <=
-               MatrixFreeFunctions::tensor_symmetric,
+               MatrixFreeFunctions::tensor_symmetric_no_collocation,
              ExcNotImplemented());
 
       EvaluatorTensorProduct<evaluate_evenodd,

include/deal.II/matrix_free/shape_info.h

@@ -81,35 +81,41 @@ namespace internal
        */
       tensor_symmetric = 2,
 
+      /**
+       * For function spaces that are not equivalent to a polynom of degree p,
+       * which the assumption of collocation is.
+       */
+      tensor_symmetric_no_collocation = 3,
+
       /**
        * Tensor product shape functions without further particular properties
        */
-      tensor_general = 3,
+      tensor_general = 4,
 
       /**
        * Polynomials of complete degree rather than tensor degree which can be
        * described by a truncated tensor product
        */
-      truncated_tensor = 4,
+      truncated_tensor = 5,
 
       /**
        * Tensor product shape functions that are symmetric about the midpoint
        * of the unit interval 0.5 that additionally add a constant shape
        * function according to FE_Q_DG0.
        */
-      tensor_symmetric_plus_dg0 = 5,
+      tensor_symmetric_plus_dg0 = 6,
 
       /**
        * Special case of the FE_RaviartThomasNodal element with anisotropic
        * tensor product shape functions, i.e degree (k + 1) in normal direction,
        * and k in tangential direction.
        */
-      tensor_raviart_thomas = 6,
+      tensor_raviart_thomas = 7,
 
       /**
        * Shape functions without a tensor product properties.
        */
-      tensor_none = 7
+      tensor_none = 8
 
 
     };

include/deal.II/matrix_free/shape_info.templates.h

@@ -34,6 +34,7 @@
 #include <deal.II/fe/fe_pyramid_p.h>
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/fe_q_dg0.h>
+#include <deal.II/fe/fe_q_iso_q1.h>
 #include <deal.II/fe/fe_raviart_thomas.h>
 #include <deal.II/fe/fe_simplex_p.h>
 #include <deal.II/fe/fe_simplex_p_bubbles.h>
@@ -913,10 +914,14 @@ namespace internal
       if (element_type == tensor_general &&
           check_1d_shapes_symmetric(univariate_shape_data))
         {
-          if (check_1d_shapes_collocation(univariate_shape_data))
+          if (dynamic_cast<const FE_Q_iso_Q1<dim> *>(&fe) &&
+              fe.tensor_degree() > 1)
+            element_type = tensor_symmetric_no_collocation;
+          else if (check_1d_shapes_collocation(univariate_shape_data))
             element_type = tensor_symmetric_collocation;
           else
             element_type = tensor_symmetric;
+
           if (n_dofs_1d > 2 && element_type == tensor_symmetric)
             {
               // check if we are a Hermite type

tests/matrix_free/fe_q_iso_q1_01.cc

@@ -0,0 +1,98 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test the evaluated values and gradients at the quadrature points for
+// FE_Q_iso_Q1.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_q_iso_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+test(const unsigned int n_subdivisions)
+{
+  const FE_Q_iso_Q1<dim> fe(n_subdivisions);
+  const QIterated<1>     quad(QGauss<1>(2), n_subdivisions);
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+
+  AffineConstraints<double> constraints;
+  constraints.close();
+
+  using number = double;
+  MatrixFree<dim, number>                          matrix_free;
+  typename MatrixFree<dim, number>::AdditionalData data;
+  data.mapping_update_flags = update_values | update_gradients;
+  matrix_free.reinit(MappingQ1<dim>(), dof, constraints, quad, data);
+
+  FEEvaluation<dim, -1, 0, 1, number> fe_eval(matrix_free);
+  fe_eval.reinit(0);
+
+  for (unsigned int i = 0; i <= dim; ++i)
+    {
+      for (unsigned int j = 0; j < fe_eval.dofs_per_cell; ++j)
+        {
+          for (unsigned int k = 0; k < fe_eval.dofs_per_cell; ++k)
+            fe_eval.begin_dof_values()[k] = static_cast<number>(j == k);
+
+          fe_eval.evaluate(EvaluationFlags::values |
+                           EvaluationFlags::gradients);
+
+          for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
+            {
+              const auto temp = (i == 0) ? fe_eval.get_value(q) :
+                                           fe_eval.get_gradient(q)[i - 1];
+              deallog << ((std::abs(temp[0]) > 1e-8) ? 1 : 0) << " ";
+            }
+
+          deallog << std::endl;
+        }
+
+      deallog << std::endl;
+    }
+}
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+  MPILogInitAll                    log;
+
+  {
+    deallog.push("2d");
+    for (unsigned int i = 1; i <= 4; ++i)
+      test<2>(i);
+    deallog.pop();
+    deallog.push("3d");
+    for (unsigned int i = 1; i <= 4; ++i)
+      test<2>(i);
+    deallog.pop();
+  }
+}