]> https://gitweb.dealii.org/ - dealii.git/commitdiff
MatrixFree ShapeInfo: Move some code to UnivariateShapeData 15905/head
authorMartin Kronbichler <martin.kronbichler@uni-a.de>
Mon, 21 Aug 2023 21:23:54 +0000 (23:23 +0200)
committerMartin Kronbichler <martin.kronbichler@uni-a.de>
Tue, 22 Aug 2023 12:03:34 +0000 (14:03 +0200)
include/deal.II/matrix_free/shape_info.h
include/deal.II/matrix_free/shape_info.templates.h

index f280d555efa7f56b79cad5f3629c82b7efe2007e..258d8f23ceb00477ebdeb2d60202fcb8e322418d 100644 (file)
@@ -144,6 +144,35 @@ namespace internal
       std::size_t
       memory_consumption() const;
 
+      /**
+       * Evaluate the auxiliary polynomial space associated with the Lagrange
+       * polynomials in points of the given quadrature formula, filling the
+       * fields shape_[gradients,hessians]_collocation and related
+       * information.
+       */
+      template <int dim, int spacedim>
+      void
+      evaluate_collocation_space(
+        const FiniteElement<dim, spacedim> &fe,
+        const Quadrature<1> &               quad,
+        const std::vector<unsigned int> &   lexicographic);
+
+      /**
+       * Check whether we have symmetries in the shape values. In that case,
+       * also fill the shape_???_eo fields.
+       */
+      bool
+      check_and_set_shapes_symmetric();
+
+      /**
+       * Check whether symmetric 1d basis functions are such that the shape
+       * values form a diagonal matrix, i.e., the nodal points are collocated
+       * with the quadrature points. This allows for specialized algorithms
+       * that save some operations in the evaluation.
+       */
+      bool
+      check_shapes_collocation() const;
+
       /**
        * Encodes the type of element detected at construction. FEEvaluation
        * will select the most efficient algorithm based on the given element
@@ -578,25 +607,6 @@ namespace internal
        * quadrature points to represent the correct order.
        */
       dealii::Table<2, unsigned int> face_orientations_quad;
-
-    private:
-      /**
-       * Check whether we have symmetries in the shape values. In that case,
-       * also fill the shape_???_eo fields.
-       */
-      bool
-      check_1d_shapes_symmetric(
-        UnivariateShapeData<Number> &univariate_shape_data);
-
-      /**
-       * Check whether symmetric 1d basis functions are such that the shape
-       * values form a diagonal matrix, i.e., the nodal points are collocated
-       * with the quadrature points. This allows for specialized algorithms
-       * that save some operations in the evaluation.
-       */
-      bool
-      check_1d_shapes_collocation(
-        const UnivariateShapeData<Number> &univariate_shape_data) const;
     };
 
 
index 8f3c72cffa912b053b16b2baa0d187add461398d..c743f6831e2d82fcf635ab0b029ab4ab92a588a3 100644 (file)
@@ -286,7 +286,7 @@ namespace internal
             lex_normal.push_back(i);
           lex_normal.push_back(dofs_per_face_normal);
 
-          // 'direction' distingusishes between normal and tangential direction
+          // 'direction' distinguishes between normal and tangential direction
           for (unsigned int direction = 0; direction < 2; ++direction)
             {
               UnivariateShapeData<Number> &univariate_shape_data =
@@ -610,17 +610,10 @@ namespace internal
         return;
 
       // grant write access to common univariate shape data
-      auto &shape_values    = univariate_shape_data.shape_values;
-      auto &shape_gradients = univariate_shape_data.shape_gradients;
-      auto &shape_hessians  = univariate_shape_data.shape_hessians;
-      auto &shape_gradients_collocation =
-        univariate_shape_data.shape_gradients_collocation;
-      auto &shape_hessians_collocation =
-        univariate_shape_data.shape_hessians_collocation;
-      auto &inverse_shape_values = univariate_shape_data.inverse_shape_values;
-      auto &shape_data_on_face   = univariate_shape_data.shape_data_on_face;
-      auto &quadrature_data_on_face =
-        univariate_shape_data.quadrature_data_on_face;
+      auto &shape_values          = univariate_shape_data.shape_values;
+      auto &shape_gradients       = univariate_shape_data.shape_gradients;
+      auto &shape_hessians        = univariate_shape_data.shape_hessians;
+      auto &shape_data_on_face    = univariate_shape_data.shape_data_on_face;
       auto &values_within_subface = univariate_shape_data.values_within_subface;
       auto &gradients_within_subface =
         univariate_shape_data.gradients_within_subface;
@@ -737,26 +730,6 @@ namespace internal
             fe.shape_grad_grad(my_i, q_point)[0][0];
         }
 
-      if (n_q_points_1d < 200)
-        {
-          quadrature_data_on_face[0].resize(quad.size() * 3);
-          quadrature_data_on_face[1].resize(quad.size() * 3);
-
-          const std::vector<Polynomials::Polynomial<double>> poly_coll =
-            Polynomials::generate_complete_Lagrange_basis(quad.get_points());
-
-          for (unsigned int i = 0; i < quad.size(); ++i)
-            {
-              std::array<double, 3> values;
-              poly_coll[i].value(0.0, 2, values.data());
-              for (unsigned int d = 0; d < 3; ++d)
-                quadrature_data_on_face[0][i + d * quad.size()] = values[d];
-              poly_coll[i].value(1.0, 2, values.data());
-              for (unsigned int d = 0; d < 3; ++d)
-                quadrature_data_on_face[1][i + d * quad.size()] = values[d];
-            }
-        }
-
       if (dim > 1 && (dynamic_cast<const FE_Q<dim> *>(&fe) ||
                       dynamic_cast<const FE_Q_iso_Q1<dim> *>(&fe)))
         {
@@ -807,136 +780,17 @@ namespace internal
               }
         }
 
-      // get gradient and Hessian transformation matrix for the polynomial
-      // space associated with the quadrature rule (collocation space). We
-      // need to avoid the case with more than a few hundreds of quadrature
-      // points when the Lagrange polynomials might underflow. Note that 200
-      // is not an exact value, as different quadrature formulas behave
-      // slightly differently, but 200 has been observed to be low enough for
-      // all common quadrature formula types. For QGauss, the actual limit is
-      // 517 points, for example.
-      if (n_q_points_1d < 200)
-        {
-          shape_gradients_collocation.resize(n_q_points_1d * n_q_points_1d);
-          shape_hessians_collocation.resize(n_q_points_1d * n_q_points_1d);
-          const std::vector<Polynomials::Polynomial<double>> poly_coll =
-            Polynomials::generate_complete_Lagrange_basis(quad.get_points());
-          std::array<double, 3> values;
-          for (unsigned int i = 0; i < n_q_points_1d; ++i)
-            for (unsigned int q = 0; q < n_q_points_1d; ++q)
-              {
-                poly_coll[i].value(quad.get_points()[q][0], 2, values.data());
-                shape_gradients_collocation[i * n_q_points_1d + q] = values[1];
-                shape_hessians_collocation[i * n_q_points_1d + q]  = values[2];
-              }
-
-          // compute the inverse shape functions in three steps: we first
-          // change from the given quadrature formula and the associated
-          // Lagrange polynomials to the Lagrange polynomials at quadrature
-          // points. in this basis, we can then perform the second step, which
-          // is the computation of a projection matrix from the potentially
-          // higher polynomial degree associated to the quadrature points to a
-          // polynomial space of degree equal to the degree of the given
-          // elements. in the third step, we change from the Lagrange
-          // polynomials in the Gauss quadrature points to the polynomial
-          // space of the given element
-
-          // step 1: change basis from the Lagrange polynomials at the given
-          // quadrature points to the Lagrange basis at Gauss quadrature
-          // points. this is often the identity operation as we often compute
-          // with Gaussian quadrature, but not necessarily so
-          QGauss<1>          quad_gauss(n_q_points_1d);
-          FullMatrix<double> transform_to_gauss(n_q_points_1d, n_q_points_1d);
-          for (unsigned int i = 0; i < n_q_points_1d; ++i)
-            for (unsigned int j = 0; j < n_q_points_1d; ++j)
-              transform_to_gauss(i, j) =
-                poly_coll[j].value(quad_gauss.point(i)[0]);
-
-          // step 2: computation for the projection (in reference coordinates)
-          // from higher to lower polynomial degree
-          //
-          // loop over quadrature points, multiply by q-weight on high degree
-          // integrate loop going from high degree to low degree loop over new
-          // points, multiply by inverse q-weight on low degree
-          //
-          // This projection step is for the special case of Lagrange
-          // polynomials where most of the interpolation matrices are unit
-          // matrices when applying the inverse mass matrix, so we do not need
-          // to compute much.
-          QGauss<1> quad_project(n_dofs_1d);
-          const std::vector<Polynomials::Polynomial<double>> poly_project =
-            Polynomials::generate_complete_Lagrange_basis(
-              quad_project.get_points());
-
-          FullMatrix<double> project_gauss(n_dofs_1d, n_q_points_1d);
-
-          for (unsigned int i = 0; i < n_dofs_1d; ++i)
-            for (unsigned int q = 0; q < n_q_points_1d; ++q)
-              project_gauss(i, q) =
-                poly_project[i].value(quad_gauss.get_points()[q][0]) *
-                (quad_gauss.weight(q) / quad_project.weight(i));
-          FullMatrix<double> project_to_dof_space(n_dofs_1d, n_q_points_1d);
-          project_gauss.mmult(project_to_dof_space, transform_to_gauss);
-
-          // step 3: change the basis back to the given finite element
-          // space. we can use a shortcut for elements that define support
-          // points, in which case we can evaluate the Lagrange polynomials of
-          // the Gauss quadrature in those points. this will give more
-          // accurate results than the inversion of a matrix. for more general
-          // polynomial spaces, we must invert a matrix of a Vandermonde type,
-          // which we do by a Householder transformation to keep roundoff
-          // errors low.
-          inverse_shape_values.resize_fast(array_size);
-          FullMatrix<double> transform_from_gauss(n_dofs_1d, n_dofs_1d);
-          if (fe.has_support_points())
-            {
-              for (unsigned int i = 0; i < n_dofs_1d; ++i)
-                for (unsigned int j = 0; j < n_dofs_1d; ++j)
-                  transform_from_gauss(i, j) = poly_project[j].value(
-                    fe.get_unit_support_points()[scalar_lexicographic[i]][0]);
-              FullMatrix<double> result(n_dofs_1d, n_q_points_1d);
-              transform_from_gauss.mmult(result, project_to_dof_space);
-
-              // set very small entries to zero - we are in reference space
-              // with normalized numbers, so this is straight-forward to check
-              // here
-              for (unsigned int i = 0; i < n_dofs_1d; ++i)
-                for (unsigned int q = 0; q < n_q_points_1d; ++q)
-                  inverse_shape_values[i * n_q_points_1d + q] =
-                    std::abs(result(i, q)) < 1e-15 ? 0 : result(i, q);
-            }
-          else
-            {
-              for (unsigned int i = 0; i < n_dofs_1d; ++i)
-                for (unsigned int j = 0; j < n_dofs_1d; ++j)
-                  {
-                    Point<dim> q_point = unit_point;
-                    q_point[0]         = quad_project.point(i)[0];
-
-                    transform_from_gauss(i, j) =
-                      fe.shape_value(scalar_lexicographic[j], q_point);
-                  }
-              Householder<double> H(transform_from_gauss);
-              Vector<double>      in(n_dofs_1d), out(n_dofs_1d);
-              for (unsigned int q = 0; q < n_q_points_1d; ++q)
-                {
-                  for (unsigned int i = 0; i < n_dofs_1d; ++i)
-                    in(i) = project_to_dof_space(i, q);
-                  H.least_squares(out, in);
-                  for (unsigned int i = 0; i < n_dofs_1d; ++i)
-                    inverse_shape_values[i * n_q_points_1d + q] =
-                      std::abs(out(i)) < 1e-15 ? 0. : out(i);
-                }
-            }
-        }
+      univariate_shape_data.evaluate_collocation_space(fe,
+                                                       quad,
+                                                       scalar_lexicographic);
 
       if (element_type == tensor_general &&
-          check_1d_shapes_symmetric(univariate_shape_data))
+          univariate_shape_data.check_and_set_shapes_symmetric())
         {
           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))
+          else if (univariate_shape_data.check_shapes_collocation())
             element_type = tensor_symmetric_collocation;
           else
             element_type = tensor_symmetric;
@@ -956,7 +810,7 @@ namespace internal
             }
         }
       else if (element_type == tensor_symmetric_plus_dg0)
-        check_1d_shapes_symmetric(univariate_shape_data);
+        univariate_shape_data.check_and_set_shapes_symmetric();
 
       nodal_at_cell_boundaries = true;
       for (unsigned int i = 1; i < n_dofs_1d; ++i)
@@ -1057,34 +911,163 @@ namespace internal
 
 
 
+    template <typename Number>
+    template <int dim, int spacedim>
+    void
+    UnivariateShapeData<Number>::evaluate_collocation_space(
+      const FiniteElement<dim, spacedim> &fe,
+      const Quadrature<1> &               quad,
+      const std::vector<unsigned int> &   lexicographic)
+    {
+      const unsigned int n_q_points_1d = quad.size();
+      const unsigned int n_dofs_1d =
+        std::min(fe.n_dofs_per_cell(), fe_degree + 1);
+
+      // get gradient and Hessian transformation matrix for the polynomial
+      // space associated with the quadrature rule (collocation space). We
+      // need to avoid the case with more than a few hundreds of quadrature
+      // points when the Lagrange polynomials might underflow. Note that 200
+      // is not an exact value, as different quadrature formulas behave
+      // slightly differently, but 200 has been observed to be low enough for
+      // all common quadrature formula types. For QGauss, the actual limit is
+      // 517 points, for example.
+      if (n_q_points_1d >= 200)
+        return;
+
+      shape_gradients_collocation.resize(n_q_points_1d * n_q_points_1d);
+      shape_hessians_collocation.resize(n_q_points_1d * n_q_points_1d);
+      const std::vector<Polynomials::Polynomial<double>> poly_coll =
+        Polynomials::generate_complete_Lagrange_basis(quad.get_points());
+      std::array<double, 3> values;
+      for (unsigned int i = 0; i < n_q_points_1d; ++i)
+        for (unsigned int q = 0; q < n_q_points_1d; ++q)
+          {
+            poly_coll[i].value(quad.get_points()[q][0], 2, values.data());
+            shape_gradients_collocation[i * n_q_points_1d + q] = values[1];
+            shape_hessians_collocation[i * n_q_points_1d + q]  = values[2];
+          }
+
+      // compute the inverse shape functions in three steps: we first
+      // change from the given quadrature formula and the associated
+      // Lagrange polynomials to the Lagrange polynomials at quadrature
+      // points. in this basis, we can then perform the second step, which
+      // is the computation of a projection matrix from the potentially
+      // higher polynomial degree associated to the quadrature points to a
+      // polynomial space of degree equal to the degree of the given
+      // elements. in the third step, we change from the Lagrange
+      // polynomials in the Gauss quadrature points to the polynomial
+      // space of the given element
+
+      // step 1: change basis from the Lagrange polynomials at the given
+      // quadrature points to the Lagrange basis at Gauss quadrature
+      // points. this is often the identity operation as we often compute
+      // with Gaussian quadrature, but not necessarily so
+      QGauss<1>          quad_gauss(n_q_points_1d);
+      FullMatrix<double> transform_to_gauss(n_q_points_1d, n_q_points_1d);
+      for (unsigned int i = 0; i < n_q_points_1d; ++i)
+        for (unsigned int j = 0; j < n_q_points_1d; ++j)
+          transform_to_gauss(i, j) = poly_coll[j].value(quad_gauss.point(i)[0]);
+
+      // step 2: computation for the projection (in reference coordinates)
+      // from higher to lower polynomial degree
+      //
+      // loop over quadrature points, multiply by q-weight on high degree
+      // integrate loop going from high degree to low degree loop over new
+      // points, multiply by inverse q-weight on low degree
+      //
+      // This projection step is for the special case of Lagrange
+      // polynomials where most of the interpolation matrices are unit
+      // matrices when applying the inverse mass matrix, so we do not need
+      // to compute much.
+      QGauss<1> quad_project(n_dofs_1d);
+      const std::vector<Polynomials::Polynomial<double>> poly_project =
+        Polynomials::generate_complete_Lagrange_basis(
+          quad_project.get_points());
+
+      FullMatrix<double> project_gauss(n_dofs_1d, n_q_points_1d);
+
+      for (unsigned int i = 0; i < n_dofs_1d; ++i)
+        for (unsigned int q = 0; q < n_q_points_1d; ++q)
+          project_gauss(i, q) =
+            poly_project[i].value(quad_gauss.get_points()[q][0]) *
+            (quad_gauss.weight(q) / quad_project.weight(i));
+      FullMatrix<double> project_to_dof_space(n_dofs_1d, n_q_points_1d);
+      project_gauss.mmult(project_to_dof_space, transform_to_gauss);
+
+      // step 3: change the basis back to the given finite element
+      // space. we can use a shortcut for elements that define support
+      // points, in which case we can evaluate the Lagrange polynomials of
+      // the Gauss quadrature in those points. this will give more
+      // accurate results than the inversion of a matrix. for more general
+      // polynomial spaces, we must invert a matrix of a Vandermonde type,
+      // which we do by a Householder transformation to keep roundoff
+      // errors low.
+      inverse_shape_values.resize_fast(shape_values.size());
+      FullMatrix<double> transform_from_gauss(n_dofs_1d, n_dofs_1d);
+      if (fe.has_support_points())
+        {
+          for (unsigned int i = 0; i < n_dofs_1d; ++i)
+            for (unsigned int j = 0; j < n_dofs_1d; ++j)
+              transform_from_gauss(i, j) = poly_project[j].value(
+                fe.get_unit_support_points()[lexicographic[i]][0]);
+          FullMatrix<double> result(n_dofs_1d, n_q_points_1d);
+          transform_from_gauss.mmult(result, project_to_dof_space);
+
+          // set very small entries to zero - we are in reference space
+          // with normalized numbers, so this is straight-forward to check
+          // here
+          for (unsigned int i = 0; i < n_dofs_1d; ++i)
+            for (unsigned int q = 0; q < n_q_points_1d; ++q)
+              inverse_shape_values[i * n_q_points_1d + q] =
+                std::abs(result(i, q)) < 1e-15 ? 0 : result(i, q);
+        }
+      else
+        {
+          for (unsigned int i = 0; i < n_dofs_1d; ++i)
+            for (unsigned int j = 0; j < n_dofs_1d; ++j)
+              {
+                Point<dim> q_point;
+                q_point[0] = quad_project.point(i)[0];
+
+                transform_from_gauss(i, j) =
+                  fe.shape_value(lexicographic[j], q_point);
+              }
+          Householder<double> H(transform_from_gauss);
+          Vector<double>      in(n_dofs_1d), out(n_dofs_1d);
+          for (unsigned int q = 0; q < n_q_points_1d; ++q)
+            {
+              for (unsigned int i = 0; i < n_dofs_1d; ++i)
+                in(i) = project_to_dof_space(i, q);
+              H.least_squares(out, in);
+              for (unsigned int i = 0; i < n_dofs_1d; ++i)
+                inverse_shape_values[i * n_q_points_1d + q] =
+                  std::abs(out(i)) < 1e-15 ? 0. : out(i);
+            }
+        }
+      quadrature_data_on_face[0].resize(quad.size() * 3);
+      quadrature_data_on_face[1].resize(quad.size() * 3);
+
+      for (unsigned int i = 0; i < quad.size(); ++i)
+        {
+          std::array<double, 3> values;
+          poly_coll[i].value(0.0, 2, values.data());
+          for (unsigned int d = 0; d < 3; ++d)
+            quadrature_data_on_face[0][i + d * quad.size()] = values[d];
+          poly_coll[i].value(1.0, 2, values.data());
+          for (unsigned int d = 0; d < 3; ++d)
+            quadrature_data_on_face[1][i + d * quad.size()] = values[d];
+        }
+    }
+
+
+
     template <typename Number>
     bool
-    ShapeInfo<Number>::check_1d_shapes_symmetric(
-      UnivariateShapeData<Number> &univariate_shape_data)
+    UnivariateShapeData<Number>::check_and_set_shapes_symmetric()
     {
-      if (dofs_per_component_on_cell == 0)
+      if (fe_degree == 0)
         return false;
 
-      const auto n_q_points_1d   = univariate_shape_data.n_q_points_1d;
-      const auto fe_degree       = univariate_shape_data.fe_degree;
-      auto &     shape_values    = univariate_shape_data.shape_values;
-      auto &     shape_gradients = univariate_shape_data.shape_gradients;
-      auto &     shape_hessians  = univariate_shape_data.shape_hessians;
-      auto &     shape_gradients_collocation =
-        univariate_shape_data.shape_gradients_collocation;
-      auto &shape_hessians_collocation =
-        univariate_shape_data.shape_hessians_collocation;
-      auto &shape_values_eo    = univariate_shape_data.shape_values_eo;
-      auto &shape_gradients_eo = univariate_shape_data.shape_gradients_eo;
-      auto &shape_hessians_eo  = univariate_shape_data.shape_hessians_eo;
-      auto &shape_gradients_collocation_eo =
-        univariate_shape_data.shape_gradients_collocation_eo;
-      auto &shape_hessians_collocation_eo =
-        univariate_shape_data.shape_hessians_collocation_eo;
-      auto &inverse_shape_values = univariate_shape_data.inverse_shape_values;
-      auto &inverse_shape_values_eo =
-        univariate_shape_data.inverse_shape_values_eo;
-
       const double zero_tol =
         std::is_same_v<Number, double> == true ? 1e-12 : 1e-7;
       // symmetry for values
@@ -1196,15 +1179,11 @@ namespace internal
 
     template <typename Number>
     bool
-    ShapeInfo<Number>::check_1d_shapes_collocation(
-      const UnivariateShapeData<Number> &univariate_shape_data) const
+    UnivariateShapeData<Number>::check_shapes_collocation() const
     {
-      if (dofs_per_component_on_cell != n_q_points)
+      if (fe_degree + 1 != n_q_points_1d)
         return false;
 
-      const auto fe_degree    = univariate_shape_data.fe_degree;
-      auto &     shape_values = univariate_shape_data.shape_values;
-
       const double zero_tol =
         std::is_same_v<Number, double> == true ? 1e-12 : 1e-7;
       // check: identity operation for shape values

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.