]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Do not hardcode VectorizedArray in MatrixFreeFunctions::ShapeInfo. 4680/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 3 Aug 2017 21:16:35 +0000 (23:16 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 3 Aug 2017 21:16:35 +0000 (23:16 +0200)
Make class templated on generic type Number that can be VectorizedArray<double>.

include/deal.II/matrix_free/cuda_matrix_free.templates.h
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/matrix_free.h
include/deal.II/matrix_free/shape_info.h
include/deal.II/matrix_free/shape_info.templates.h
source/matrix_free/matrix_free.cc

index 8c9353162536b8122804c5a03d11d07785148517..b90a7bb3dbbdf2a66364432241109390536f9294 100644 (file)
@@ -461,7 +461,7 @@ namespace CUDAWrappers
     unsigned int size_shape_values = n_dofs_1d*n_q_points_1d*sizeof(Number);
 
     cudaError_t cuda_error = cudaMemcpyToSymbol(internal::global_shape_values,
-                                                &shape_info.shape_values_number[0],
+                                                &shape_info.shape_values[0],
                                                 size_shape_values,
                                                 0,
                                                 cudaMemcpyHostToDevice);
@@ -470,7 +470,7 @@ namespace CUDAWrappers
     if (update_flags & update_gradients)
       {
         cuda_error = cudaMemcpyToSymbol(internal::global_shape_gradients,
-                                        &shape_info.shape_gradient_number[0],
+                                        &shape_info.shape_gradient[0],
                                         size_shape_values,
                                         0,
                                         cudaMemcpyHostToDevice);
index 13de640678eee5b5814500e79d8a031520dd8562..d7cb55e0214f8ca6faf8807a50699148eafa1952 100644 (file)
@@ -100,7 +100,7 @@ namespace internal
   struct FEEvaluationImpl
   {
     static
-    void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+    void evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
                    VectorizedArray<Number> *values_dofs_actual[],
                    VectorizedArray<Number> *values_quad[],
                    VectorizedArray<Number> *gradients_quad[][dim],
@@ -111,7 +111,7 @@ namespace internal
                    const bool               evaluate_hessians);
 
     static
-    void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+    void integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
                     VectorizedArray<Number> *values_dofs_actual[],
                     VectorizedArray<Number> *values_quad[],
                     VectorizedArray<Number> *gradients_quad[][dim],
@@ -126,7 +126,7 @@ namespace internal
   inline
   void
   FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
-  ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+  ::evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
               VectorizedArray<Number> *values_dofs_actual[],
               VectorizedArray<Number> *values_quad[],
               VectorizedArray<Number> *gradients_quad[][dim],
@@ -351,7 +351,7 @@ namespace internal
   inline
   void
   FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
-  ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+  ::integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
                VectorizedArray<Number> *values_dofs_actual[],
                VectorizedArray<Number> *values_quad[],
                VectorizedArray<Number> *gradients_quad[][dim],
@@ -554,7 +554,7 @@ namespace internal
   struct FEEvaluationImplTransformToCollocation
   {
     static
-    void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+    void evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
                    VectorizedArray<Number> *values_dofs[],
                    VectorizedArray<Number> *values_quad[],
                    VectorizedArray<Number> *gradients_quad[][dim],
@@ -565,7 +565,7 @@ namespace internal
                    const bool               evaluate_hessians);
 
     static
-    void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+    void integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
                     VectorizedArray<Number> *values_dofs[],
                     VectorizedArray<Number> *values_quad[],
                     VectorizedArray<Number> *gradients_quad[][dim],
@@ -578,7 +578,7 @@ namespace internal
   inline
   void
   FEEvaluationImplTransformToCollocation<dim, fe_degree, n_components, Number>
-  ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+  ::evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
               VectorizedArray<Number> *values_dofs[],
               VectorizedArray<Number> *values_quad[],
               VectorizedArray<Number> *gradients_quad[][dim],
@@ -662,7 +662,7 @@ namespace internal
   inline
   void
   FEEvaluationImplTransformToCollocation<dim, fe_degree, n_components, Number>
-  ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+  ::integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
                VectorizedArray<Number> *values_dofs[],
                VectorizedArray<Number> *values_quad[],
                VectorizedArray<Number> *gradients_quad[][dim],
@@ -741,7 +741,7 @@ namespace internal
   struct FEEvaluationImplCollocation
   {
     static
-    void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+    void evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
                    VectorizedArray<Number> *values_dofs[],
                    VectorizedArray<Number> *values_quad[],
                    VectorizedArray<Number> *gradients_quad[][dim],
@@ -752,7 +752,7 @@ namespace internal
                    const bool               evaluate_hessians);
 
     static
-    void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+    void integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
                     VectorizedArray<Number> *values_dofs[],
                     VectorizedArray<Number> *values_quad[],
                     VectorizedArray<Number> *gradients_quad[][dim],
@@ -765,7 +765,7 @@ namespace internal
   inline
   void
   FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>
-  ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+  ::evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
               VectorizedArray<Number> *values_dofs[],
               VectorizedArray<Number> *values_quad[],
               VectorizedArray<Number> *gradients_quad[][dim],
@@ -829,7 +829,7 @@ namespace internal
   inline
   void
   FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>
-  ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+  ::integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
                VectorizedArray<Number> *values_dofs[],
                VectorizedArray<Number> *values_quad[],
                VectorizedArray<Number> *gradients_quad[][dim],
index 3058fb289145f93994376e14f3958bfd1ab047a3..eb000af7376bed902e6d4d46de7415ebfeed8709 100644 (file)
@@ -153,8 +153,8 @@ public:
   /**
    * Return a reference to the ShapeInfo object currently in use.
    */
-  const internal::MatrixFreeFunctions::ShapeInfo<Number> &
-  get_shape_info() const;
+  const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &
+      get_shape_info() const;
 
   /**
    * Fills the JxW values currently used.
@@ -808,7 +808,7 @@ protected:
    * product. Also contained in matrix_info, but it simplifies code if we
    * store a reference to it.
    */
-  const internal::MatrixFreeFunctions::ShapeInfo<Number> *data;
+  const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> *data;
 
   /**
    * A pointer to the Cartesian Jacobian information of the present cell. Only
@@ -2168,7 +2168,7 @@ FEEvaluationBase<dim,n_components_,Number>
   dof_info           (nullptr),
   mapping_info       (nullptr),
   // select the correct base element from the given FE component
-  data               (new internal::MatrixFreeFunctions::ShapeInfo<Number>(quadrature, fe, fe.component_to_base_index(first_selected_component).first)),
+  data               (new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>(quadrature, fe, fe.component_to_base_index(first_selected_component).first)),
   cartesian_data     (nullptr),
   jacobian           (nullptr),
   J_value            (nullptr),
@@ -2230,7 +2230,7 @@ FEEvaluationBase<dim,n_components_,Number>
   dof_info           (other.dof_info),
   mapping_info       (other.mapping_info),
   data               (other.matrix_info == nullptr ?
-                      new internal::MatrixFreeFunctions::ShapeInfo<Number>(*other.data) :
+                      new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>(*other.data) :
                       other.data),
   cartesian_data     (nullptr),
   jacobian           (nullptr),
@@ -2301,7 +2301,7 @@ FEEvaluationBase<dim,n_components_,Number>
   mapping_info = other.mapping_info;
   if (other.matrix_info == nullptr)
     {
-      data = new internal::MatrixFreeFunctions::ShapeInfo<Number>(*other.data);
+      data = new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>(*other.data);
       scratch_data_array = new AlignedVector<VectorizedArray<Number> >();
     }
   else
@@ -2552,8 +2552,8 @@ FEEvaluationBase<dim,n_components_,Number>::get_cell_type () const
 
 template <int dim, int n_components_, typename Number>
 inline
-const internal::MatrixFreeFunctions::ShapeInfo<Number> &
-FEEvaluationBase<dim,n_components_,Number>::get_shape_info() const
+const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &
+    FEEvaluationBase<dim,n_components_,Number>::get_shape_info() const
 {
   Assert(data != nullptr, ExcInternalError());
   return *data;
index 8a9d79b75290055c1572bdf6477d2c0fc8b6853e..7f1c9af072da9f0a0e956b78c65a5031c6cac124 100644 (file)
@@ -841,11 +841,11 @@ public:
   /**
    * Return the unit cell information for given hp index.
    */
-  const internal::MatrixFreeFunctions::ShapeInfo<Number> &
-  get_shape_info (const unsigned int fe_component = 0,
-                  const unsigned int quad_index   = 0,
-                  const unsigned int hp_active_fe_index = 0,
-                  const unsigned int hp_active_quad_index = 0) const;
+  const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &
+      get_shape_info (const unsigned int fe_component = 0,
+                      const unsigned int quad_index   = 0,
+                      const unsigned int hp_active_fe_index = 0,
+                      const unsigned int hp_active_quad_index = 0) const;
 
   /**
    * Obtains a scratch data object for internal use. Make sure to release it
@@ -980,7 +980,7 @@ private:
   /**
    * Contains shape value information on the unit cell.
    */
-  Table<4,internal::MatrixFreeFunctions::ShapeInfo<Number> > shape_info;
+  Table<4,internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>> shape_info;
 
   /**
    * Describes how the cells are gone through. With the cell level (first
@@ -1452,11 +1452,11 @@ MatrixFree<dim,Number>::get_ghost_set(const unsigned int dof_index) const
 
 template <int dim, typename Number>
 inline
-const internal::MatrixFreeFunctions::ShapeInfo<Number> &
-MatrixFree<dim,Number>::get_shape_info (const unsigned int index_fe,
-                                        const unsigned int index_quad,
-                                        const unsigned int active_fe_index,
-                                        const unsigned int active_quad_index) const
+const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &
+    MatrixFree<dim,Number>::get_shape_info (const unsigned int index_fe,
+                                            const unsigned int index_quad,
+                                            const unsigned int active_fe_index,
+                                            const unsigned int active_quad_index) const
 {
   AssertIndexRange (index_fe, shape_info.size(0));
   AssertIndexRange (index_quad, shape_info.size(1));
index e89de80f49e7fbb8d0cf6e714a6ccf52bb4e5a53..424ebd69b3641e063e46015704faa5b333aa6134 100644 (file)
@@ -20,7 +20,6 @@
 
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/vectorization.h>
 #include <deal.II/base/aligned_vector.h>
 #include <deal.II/fe/fe.h>
 
@@ -132,7 +131,7 @@ namespace internal
        * this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
        * points are the index running fastest.
        */
-      AlignedVector<VectorizedArray<Number> > shape_values;
+      AlignedVector<Number> shape_values;
 
       /**
        * Stores the shape gradients of the 1D finite element evaluated on all
@@ -141,7 +140,7 @@ namespace internal
        * this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
        * points are the index running fastest.
        */
-      AlignedVector<VectorizedArray<Number> > shape_gradients;
+      AlignedVector<Number> shape_gradients;
 
       /**
        * Stores the shape Hessians of the 1D finite element evaluated on all
@@ -150,67 +149,67 @@ namespace internal
        * this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
        * points are the index running fastest.
        */
-      AlignedVector<VectorizedArray<Number> > shape_hessians;
+      AlignedVector<Number> shape_hessians;
 
       /**
        * Stores the shape values in a different format, namely the so-called
        * even-odd scheme where the symmetries in shape_values are used for
        * faster evaluation.
        */
-      AlignedVector<VectorizedArray<Number> > shape_values_eo;
+      AlignedVector<Number> shape_values_eo;
 
       /**
        * Stores the shape gradients in a different format, namely the so-
        * called even-odd scheme where the symmetries in shape_gradients are
        * used for faster evaluation.
        */
-      AlignedVector<VectorizedArray<Number> > shape_gradients_eo;
+      AlignedVector<Number> shape_gradients_eo;
 
       /**
        * Stores the shape second derivatives in a different format, namely the
        * so-called even-odd scheme where the symmetries in shape_hessians are
        * used for faster evaluation.
        */
-      AlignedVector<VectorizedArray<Number> > shape_hessians_eo;
+      AlignedVector<Number> shape_hessians_eo;
 
       /**
        * Stores the shape gradients of the shape function space associated to
        * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). For
        * faster evaluation only the even-odd format is necessary.
        */
-      AlignedVector<VectorizedArray<Number> > shape_gradients_collocation_eo;
+      AlignedVector<Number> shape_gradients_collocation_eo;
 
       /**
        * Stores the shape hessians of the shape function space associated to
        * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). For
        * faster evaluation only the even-odd format is necessary.
        */
-      AlignedVector<VectorizedArray<Number> > shape_hessians_collocation_eo;
+      AlignedVector<Number> shape_hessians_collocation_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,
        * then gradients, then second derivatives.
        */
-      AlignedVector<VectorizedArray<Number> > shape_data_on_face[2];
+      AlignedVector<Number> shape_data_on_face[2];
 
       /**
        * Stores one-dimensional values of shape functions on subface. Since
        * there are two subfaces, store two variants.
        */
-      AlignedVector<VectorizedArray<Number> > values_within_subface[2];
+      AlignedVector<Number> values_within_subface[2];
 
       /**
        * Stores one-dimensional gradients of shape functions on subface. Since
        * there are two subfaces, store two variants.
        */
-      AlignedVector<VectorizedArray<Number> > gradients_within_subface[2];
+      AlignedVector<Number> gradients_within_subface[2];
 
       /**
        * Stores one-dimensional gradients of shape functions on subface. Since
        * there are two subfaces, store two variants.
        */
-      AlignedVector<VectorizedArray<Number> > hessians_within_subface[2];
+      AlignedVector<Number> hessians_within_subface[2];
 
       /**
        * Renumbering from deal.II's numbering of cell degrees of freedom to
@@ -259,31 +258,77 @@ namespace internal
       bool nodal_at_cell_boundaries;
 
       /**
-       * For nodal cells, we might get around by simply loading the indices to
-       * the degrees of freedom that act on a particular face, rather than the
-       * whole set of indices that is then interpolated down to the
-       * element. This array stores this indirect addressing.
+       * For nodal basis functions with nodes located at the boundary of the
+       * unit cell, face integrals that involve only the values of the shape
+       * functions (approximations of first derivatives in DG) do not need to
+       * load all degrees of freedom of the cell but rather only the degrees
+       * of freedom located on the face. While it would also be possible to
+       * compute these indices on the fly, we choose to simplify the code and
+       * store the indirect addressing in this class.
        *
        * The first table index runs through the faces of a cell, and the
        * second runs through the nodal degrees of freedom of the face, using
        * @p dofs_per_face entries.
        *
+       * The indices stored in this member variable are as follows. Consider
+       * for example a 2D element of degree 3 with the following degrees of
+       * freedom in lexicographic numbering:
+       * @code
+       * 12   13   14   15
+       * 8    9    10   11
+       * 4    5     6    7
+       * 0    1     2    3
+       * @endcode
+       *
+       * The first row stores the indices on the face with index 0, i.e., the
+       * numbers <code>0, 4, 8, 12</code>, the second row holds the indices
+       * <code>3, 7, 11, 15</code> for face 1, the third row holds the indices
+       * <code>0, 1, 2, 3</code> for face 2, and the last (fourth) row holds
+       * the indices <code>12, 13, 14, 15</code>. Similarly, the indices are
+       * stored in 3D. (Note that the y faces in 3D use indices reversed in
+       * terms of the lexicographic numbers due to the orientation of the
+       * coordinate system.)
+       *
        * @note This object is only filled in case @p nodal_at_cell_boundaries
        * evaluates to @p true.
        */
       dealii::Table<2,unsigned int> face_to_cell_index_nodal;
 
       /**
-       * For Hermite-type basis functions, the @p face_to_cell_index_nodal for
-       * the values on a face of the cell is used together with a respective
-       * slot in the derivative. In the lexicographic ordering, this index is
-       * in the next "layer" of degrees of freedom. This array stores the
-       * indirect addressing of both the values and the gradient.
+       * The @p face_to_cell_index_nodal provides a shortcut for the
+       * evaluation of values on the faces. For Hermite-type basis functions,
+       * one can go one step further and also use shortcuts to get derivatives
+       * more cheaply where only two layers of degrees of freedom contribute
+       * to the derivative on the face. In the lexicographic ordering, the
+       * respective indices is in the next "layer" of degrees of freedom as
+       * compared to the nodal values. This array stores the indirect
+       * addressing of both the values and the gradient.
        *
        * The first table index runs through the faces of a cell, and the
        * second runs through the pairs of the nodal degrees of freedom of the
        * face and the derivatives, using <code>2*dofs_per_face</code> entries.
        *
+       * The indices stored in this member variable are as follows. Consider
+       * for example a 2D element of degree 3 with the following degrees of
+       * freedom in lexicographic numbering:
+       * @code
+       * 20   21   22   23   24
+       * 15   16   17   18   19
+       * 10   11   12   13   14
+       * 5    6     7    8    9
+       * 0    1     2    3    4
+       * @endcode
+       *
+       * The first row stores the indices for values and gradients on the face
+       * with index 0, i.e., the numbers <code>0, 1, 5, 6, 10, 11, 15, 16, 20,
+       * 21</code>, the second row holds the indices <code>4, 3, 9, 8, 14, 13,
+       * 19, 18, 24, 23</code> for face 1, the third row holds the indices
+       * <code>0, 5, 1, 6, 2, 7, 3, 8, 4, 9</code> for face 2, and the last
+       * (fourth) row holds the indices <code>20, 15, 21, 16, 22, 17, 23, 18,
+       * 24, 19</code>. Similarly, the indices are stored in 3D. (Note that
+       * the y faces in 3D use indices reversed in terms of the lexicographic
+       * numbers due to the orientation of the coordinate system.)
+       *
        * @note This object is only filled in case @p element_type evaluates to
        * @p tensor_symmetric_hermite.
        */
index cefbd43e061a84f3b6ee0e072c5d10ca31f37d5c..9bc6bea4fe55fcf23fa6f5c6fe3087128eee2508 100644 (file)
@@ -40,6 +40,21 @@ namespace internal
 
     // ----------------- actual ShapeInfo functions --------------------
 
+    namespace
+    {
+      template <typename Number>
+      Number get_first_array_element(const Number a)
+      {
+        return a;
+      }
+
+      template <typename Number>
+      Number get_first_array_element(const VectorizedArray<Number> a)
+      {
+        return a[0];
+      }
+    }
+
     template <typename Number>
     ShapeInfo<Number>::ShapeInfo ()
       :
@@ -159,8 +174,8 @@ namespace internal
         if (fe->has_support_points())
           unit_point = fe->get_unit_support_points()[scalar_lexicographic[0]];
         Assert(fe->dofs_per_cell == 0 ||
-               std::fabs(fe->shape_value(scalar_lexicographic[0],
-                                         unit_point)-1) < 1e-13,
+               std::abs(fe->shape_value(scalar_lexicographic[0],
+                                        unit_point)-1) < 1e-13,
                ExcInternalError("Could not decode 1D shape functions for the "
                                 "element " + fe->get_name()));
       }
@@ -191,9 +206,6 @@ namespace internal
           const unsigned int my_i = scalar_lexicographic[i];
           for (unsigned int q=0; q<n_q_points_1d; ++q)
             {
-              // fill both vectors with
-              // VectorizedArray<Number>::n_array_elements copies for the
-              // shape information and non-vectorized fields
               Point<dim> q_point = unit_point;
               q_point[0] = quad.get_points()[q][0];
 
@@ -271,10 +283,10 @@ namespace internal
               // check if we are a Hermite type
               element_type = tensor_symmetric_hermite;
               for (unsigned int i=1; i<n_dofs_1d; ++i)
-                if (std::abs(this->shape_data_on_face[0][i][0]) > 1e-12)
+                if (std::abs(get_first_array_element(shape_data_on_face[0][i])) > 1e-12)
                   element_type = tensor_symmetric;
               for (unsigned int i=2; i<n_dofs_1d; ++i)
-                if (std::abs(this->shape_data_on_face[0][n_dofs_1d+i][0]) > 1e-12)
+                if (std::abs(get_first_array_element(shape_data_on_face[0][n_dofs_1d+i])) > 1e-12)
                   element_type = tensor_symmetric;
             }
         }
@@ -283,8 +295,8 @@ namespace internal
 
       nodal_at_cell_boundaries = true;
       for (unsigned int i=1; i<n_dofs_1d; ++i)
-        if (std::abs(this->shape_data_on_face[0][i][0]) > 1e-13 ||
-            std::abs(this->shape_data_on_face[1][i-1][0]) > 1e-13)
+        if (std::abs(get_first_array_element(shape_data_on_face[0][i])) > 1e-13 ||
+            std::abs(get_first_array_element(shape_data_on_face[1][i-1])) > 1e-13)
           nodal_at_cell_boundaries = false;
 
       if (nodal_at_cell_boundaries == true)
@@ -371,11 +383,10 @@ namespace internal
       const unsigned int n_dofs_1d = fe_degree + 1;
       for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
         for (unsigned int j=0; j<n_q_points_1d; ++j)
-          if (std::fabs(shape_values[i*n_q_points_1d+j][0] -
-                        shape_values[(n_dofs_1d-i)*n_q_points_1d
-                                     -j-1][0]) >
+          if (std::abs(get_first_array_element(shape_values[i*n_q_points_1d+j] -
+                                               shape_values[(n_dofs_1d-i)*n_q_points_1d-j-1])) >
               std::max(zero_tol, zero_tol*
-                       std::abs(shape_values[i*n_q_points_1d+j][0])))
+                       std::abs(get_first_array_element(shape_values[i*n_q_points_1d+j]))))
             return false;
 
       // shape values should be zero at x=0.5 for all basis functions except
@@ -383,11 +394,11 @@ namespace internal
       if (n_q_points_1d%2 == 1 && n_dofs_1d%2 == 1)
         {
           for (unsigned int i=0; i<n_dofs_1d/2; ++i)
-            if (std::fabs(shape_values[i*n_q_points_1d+
-                                       n_q_points_1d/2][0]) > zero_tol)
+            if (std::abs(get_first_array_element(shape_values[i*n_q_points_1d+
+                                                              n_q_points_1d/2])) > zero_tol)
               return false;
-          if (std::fabs(shape_values[(n_dofs_1d/2)*n_q_points_1d+
-                                     n_q_points_1d/2][0]-1.)> zero_tol)
+          if (std::abs(get_first_array_element(shape_values[(n_dofs_1d/2)*n_q_points_1d+
+                                                            n_q_points_1d/2])-1.)> zero_tol)
             return false;
         }
 
@@ -397,13 +408,14 @@ namespace internal
       const double zero_tol_gradient = zero_tol * std::sqrt(fe_degree+1.)*(fe_degree+1);
       for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
         for (unsigned int j=0; j<n_q_points_1d; ++j)
-          if (std::fabs(shape_gradients[i*n_q_points_1d+j][0] +
-                        shape_gradients[(n_dofs_1d-i)*n_q_points_1d-
-                                        j-1][0]) > zero_tol_gradient)
+          if (std::abs(get_first_array_element(shape_gradients[i*n_q_points_1d+j] +
+                                               shape_gradients[(n_dofs_1d-i)*n_q_points_1d-
+                                                               j-1])) > zero_tol_gradient)
             return false;
       if (n_dofs_1d%2 == 1 && n_q_points_1d%2 == 1)
-        if (std::fabs(shape_gradients[(n_dofs_1d/2)*n_q_points_1d+
-                                      (n_q_points_1d/2)][0]) > zero_tol_gradient)
+        if (std::abs(get_first_array_element(shape_gradients[(n_dofs_1d/2)*n_q_points_1d+
+                                                             (n_q_points_1d/2)]))
+            > zero_tol_gradient)
           return false;
 
       // symmetry for Hessian. Multiply tolerance by degree^3 of the element
@@ -411,9 +423,9 @@ namespace internal
       const double zero_tol_hessian = zero_tol * (fe_degree+1)*(fe_degree+1)*(fe_degree+1);
       for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
         for (unsigned int j=0; j<n_q_points_1d; ++j)
-          if (std::fabs(shape_hessians[i*n_q_points_1d+j][0] -
-                        shape_hessians[(n_dofs_1d-i)*n_q_points_1d-
-                                       j-1][0]) > zero_tol_hessian)
+          if (std::abs(get_first_array_element(shape_hessians[i*n_q_points_1d+j] -
+                                               shape_hessians[(n_dofs_1d-i)*n_q_points_1d-
+                                                              j-1])) > zero_tol_hessian)
             return false;
 
       const unsigned int stride = (n_q_points_1d+1)/2;
@@ -424,25 +436,25 @@ namespace internal
         for (unsigned int q=0; q<stride; ++q)
           {
             shape_values_eo[i*stride+q] =
-              Number(0.5) * (shape_values[i*n_q_points_1d+q] +
-                             shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
+              0.5 * (shape_values[i*n_q_points_1d+q] +
+                     shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
             shape_values_eo[(fe_degree-i)*stride+q] =
-              Number(0.5) * (shape_values[i*n_q_points_1d+q] -
-                             shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
+              0.5 * (shape_values[i*n_q_points_1d+q] -
+                     shape_values[i*n_q_points_1d+n_q_points_1d-1-q]);
 
             shape_gradients_eo[i*stride+q] =
-              Number(0.5) * (shape_gradients[i*n_q_points_1d+q] +
-                             shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]);
+              0.5 * (shape_gradients[i*n_q_points_1d+q] +
+                     shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]);
             shape_gradients_eo[(fe_degree-i)*stride+q] =
-              Number(0.5) * (shape_gradients[i*n_q_points_1d+q] -
-                             shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]);
+              0.5 * (shape_gradients[i*n_q_points_1d+q] -
+                     shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]);
 
             shape_hessians_eo[i*stride+q] =
-              Number(0.5) * (shape_hessians[i*n_q_points_1d+q] +
-                             shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]);
+              0.5 * (shape_hessians[i*n_q_points_1d+q] +
+                     shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]);
             shape_hessians_eo[(fe_degree-i)*stride+q] =
-              Number(0.5) * (shape_hessians[i*n_q_points_1d+q] -
-                             shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]);
+              0.5 * (shape_hessians[i*n_q_points_1d+q] -
+                     shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]);
           }
       if (fe_degree % 2 == 0)
         for (unsigned int q=0; q<stride; ++q)
@@ -475,13 +487,12 @@ namespace internal
         for (unsigned int j=0; j<n_points_1d; ++j)
           if (i!=j)
             {
-              if (std::fabs(shape_values[i*n_points_1d+j][0])>zero_tol)
+              if (std::abs(get_first_array_element(shape_values[i*n_points_1d+j]))>zero_tol)
                 return false;
             }
           else
             {
-              if (std::fabs(shape_values[i*n_points_1d+
-                                         j][0]-1.)>zero_tol)
+              if (std::abs(get_first_array_element(shape_values[i*n_points_1d+j])-1.)>zero_tol)
                 return false;
             }
       return true;
index d6a07ec0f43c696b0fd21f6d4f0da12a6f7679e6..ea2c2e421d044f3a67aa6b043a64a6aa24b70dc9 100644 (file)
@@ -16,6 +16,7 @@
 
 #include <deal.II/matrix_free/matrix_free.templates.h>
 #include <deal.II/base/utilities.h>
+#include <deal.II/base/vectorization.h>
 #include <deal.II/base/conditional_ostream.h>
 
 #include <iostream>
@@ -26,5 +27,7 @@ DEAL_II_NAMESPACE_OPEN
 
 template struct internal::MatrixFreeFunctions::ShapeInfo<double>;
 template struct internal::MatrixFreeFunctions::ShapeInfo<float>;
+template struct internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<double>>;
+template struct internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<float>>;
 
 DEAL_II_NAMESPACE_CLOSE

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.