]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Compute inverse of shape values as part of MF::ShapeInfo
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sun, 8 Dec 2019 19:09:03 +0000 (20:09 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 9 Dec 2019 12:05:14 +0000 (13:05 +0100)
include/deal.II/matrix_free/shape_info.h
include/deal.II/matrix_free/shape_info.templates.h
include/deal.II/matrix_free/tensor_product_kernels.h
tests/matrix_free/shape_info.output

index 5473845a58b69b938080b5dfc835c398461e8748..afc38492bed80d06fd5ce3766895d16d8e9b3ab0 100644 (file)
@@ -223,6 +223,27 @@ namespace internal
        */
       AlignedVector<Number> shape_hessians_collocation_eo;
 
+      /**
+       * Stores the inverse transformation from the data at quadrature points
+       * to the basis defined by the shape_values fields. The data at
+       * quadrature points is interpreted either implicitly by its polynomial
+       * interpolation, or explicitly in terms of separate polynomials such as
+       * with the `_collocation` fields. The size of the array equals the
+       * layout of the `shape_values` array, and it is combined with the shape
+       * values array such that this matrix is the pseudo inverse of
+       * shape_values. In case the number of 1D quadrature points equals the
+       * size of the basis, this array is exactly the inverse of the
+       * shape_values array. The length of this array is <tt>n_dofs_1d *
+       * n_q_points_1d</tt> and quadrature points are the index running
+       * fastest.
+       */
+      AlignedVector<Number> inverse_shape_values;
+
+      /**
+       * Stores the even-odd variant of the `inverse_shape_values` field.
+       */
+      AlignedVector<Number> inverse_shape_values_eo;
+
       /**
        * Collects all data of 1D shape values evaluated at the point 0 and 1
        * (the vertices) in one data structure. Sorting is first the values,
index 61eb965e0bcfe4563f50e71754d5abcb045a5ae1..3be0eea9b96608614aef2570294616df0eadd4c9 100644 (file)
@@ -30,6 +30,8 @@
 #include <deal.II/fe/fe_poly.h>
 #include <deal.II/fe/fe_q_dg0.h>
 
+#include <deal.II/lac/householder.h>
+
 #include <deal.II/matrix_free/shape_info.h>
 
 
@@ -277,15 +279,115 @@ namespace internal
         {
           shape_gradients_collocation.resize(n_q_points_1d * n_q_points_1d);
           shape_hessians_collocation.resize(n_q_points_1d * n_q_points_1d);
-          FE_DGQArbitraryNodes<1> fe(quad.get_points());
+          FE_DGQArbitraryNodes<1> fe_coll(quad.get_points());
           for (unsigned int i = 0; i < n_q_points_1d; ++i)
             for (unsigned int q = 0; q < n_q_points_1d; ++q)
               {
                 shape_gradients_collocation[i * n_q_points_1d + q] =
-                  fe.shape_grad(i, quad.get_points()[q])[0];
+                  fe_coll.shape_grad(i, quad.get_points()[q])[0];
                 shape_hessians_collocation[i * n_q_points_1d + q] =
-                  fe.shape_grad_grad(i, quad.get_points()[q])[0][0];
+                  fe_coll.shape_grad_grad(i, quad.get_points()[q])[0][0];
               }
+
+          // 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) =
+                fe_coll.shape_value(j, quad_gauss.point(i));
+
+          // 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);
+          FE_DGQArbitraryNodes<1> fe_project(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) =
+                fe_project.shape_value(i, quad_gauss.get_points()[q]) *
+                (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) = fe_project.shape_value(
+                    j,
+                    Point<1>(
+                      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);
+                }
+            }
         }
 
       if (element_type == tensor_general &&
@@ -419,8 +521,8 @@ namespace internal
             return false;
 
       // shape values should be zero at x=0.5 for all basis functions except
-      // for the middle one
-      if (n_q_points_1d % 2 == 1 && n_dofs_1d % 2 == 1)
+      // for the middle one for degrees of 4 and higher
+      if (n_dofs_1d > 3 && n_q_points_1d % 2 == 1 && n_dofs_1d % 2 == 1)
         {
           for (unsigned int i = 0; i < n_dofs_1d / 2; ++i)
             if (std::abs(get_first_array_element(
@@ -460,36 +562,36 @@ namespace internal
             return false;
 
       auto convert_to_eo = [](const AlignedVector<Number> &array,
-                              const unsigned               n_points_dst,
-                              const unsigned               n_points_src) {
-        const unsigned int    stride = (n_points_dst + 1) / 2;
-        AlignedVector<Number> array_eo(n_points_src * stride);
-        for (unsigned int i = 0; i < n_points_src / 2; ++i)
+                              const unsigned               n_rows,
+                              const unsigned               n_cols) {
+        const unsigned int    stride = (n_cols + 1) / 2;
+        AlignedVector<Number> array_eo(n_rows * stride);
+        for (unsigned int i = 0; i < n_rows / 2; ++i)
           for (unsigned int q = 0; q < stride; ++q)
             {
               array_eo[i * stride + q] =
-                0.5 * (array[i * n_points_dst + q] +
-                       array[i * n_points_dst + n_points_dst - 1 - q]);
-              array_eo[(n_points_src - 1 - i) * stride + q] =
-                0.5 * (array[i * n_points_dst + q] -
-                       array[i * n_points_dst + n_points_dst - 1 - q]);
+                0.5 *
+                (array[i * n_cols + q] + array[i * n_cols + n_cols - 1 - q]);
+              array_eo[(n_rows - 1 - i) * stride + q] =
+                0.5 *
+                (array[i * n_cols + q] - array[i * n_cols + n_cols - 1 - q]);
             }
-        if ((n_points_src - 1) % 2 == 0)
+        if ((n_rows - 1) % 2 == 0)
           for (unsigned int q = 0; q < stride; ++q)
             {
-              array_eo[(n_points_src - 1) / 2 * stride + q] =
-                array[((n_points_src - 1) / 2) * n_points_dst + q];
+              array_eo[(n_rows - 1) / 2 * stride + q] =
+                array[((n_rows - 1) / 2) * n_cols + q];
             }
 
         return array_eo;
       };
 
       shape_values_eo =
-        convert_to_eo(shape_values, n_q_points_1d, fe_degree + 1);
+        convert_to_eo(shape_values, fe_degree + 1, n_q_points_1d);
       shape_gradients_eo =
-        convert_to_eo(shape_gradients, n_q_points_1d, fe_degree + 1);
+        convert_to_eo(shape_gradients, fe_degree + 1, n_q_points_1d);
       shape_hessians_eo =
-        convert_to_eo(shape_hessians, n_q_points_1d, fe_degree + 1);
+        convert_to_eo(shape_hessians, fe_degree + 1, n_q_points_1d);
 
       // FE_DGQArbitraryNodes underflow (see also above where
       // shape_gradients_collocation and shape_hessians_collocation is set up).
@@ -503,6 +605,8 @@ namespace internal
             convert_to_eo(shape_hessians_collocation,
                           n_q_points_1d,
                           n_q_points_1d);
+          inverse_shape_values_eo =
+            convert_to_eo(inverse_shape_values, fe_degree + 1, n_q_points_1d);
         }
 
       return true;
@@ -572,10 +676,9 @@ namespace internal
       return memory;
     }
 
-    // end of functions for ShapeInfo
+  } // namespace MatrixFreeFunctions
 
-  } // end of namespace MatrixFreeFunctions
-} // end of namespace internal
+} // namespace internal
 
 
 DEAL_II_NAMESPACE_CLOSE
index 7c799fc226357bfa61acef82cfb84fb70e948955..98434e847b522996d1bfbb83e2e40098aa9667cf 100644 (file)
@@ -1716,7 +1716,7 @@ namespace internal
                     else
                       r0 += shapes[mid * offset + col] * xmid;
                   }
-                else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0))
+                else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0 || mm == 3))
                   r0 += shapes[col * offset + mid] * xmid;
 
                 if (add == false)
@@ -1737,7 +1737,7 @@ namespace internal
                   }
               }
             if (type == 0 && contract_over_rows == true && nn % 2 == 1 &&
-                mm % 2 == 1)
+                mm % 2 == 1 && mm > 3)
               {
                 if (add == false)
                   out[stride * n_cols] = shapes[mid * offset + n_cols] * xmid;
index 87b8dbfbd584185534925fc615fae1ccdf8be821..186fa91f1a5367fca2f01be4e9a16093acbc03d3 100644 (file)
@@ -31,7 +31,7 @@ DEAL::Detected shape info type for FE_DGQ<2>(17): 2
 DEAL::Detected shape info type for FE_DGQ<2>(25): 2
 DEAL::Detected shape info type for FE_DGQHermite<2>(1): 2
 DEAL::Detected shape info type for FE_DGQHermite<2>(1): 2
-DEAL::Detected shape info type for FE_DGQHermite<2>(2): 3
+DEAL::Detected shape info type for FE_DGQHermite<2>(2): 1
 DEAL::Detected shape info type for FE_DGQHermite<2>(2): 1
 DEAL::Detected shape info type for FE_DGQHermite<2>(3): 1
 DEAL::Detected shape info type for FE_DGQHermite<2>(3): 1

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.