Convert internal vector<vector<T>> to Table<2,T>.

This should improve data locality and simplify memory management.

diff --git a/include/deal.II/fe/fe_poly.templates.h b/include/deal.II/fe/fe_poly.templates.h
index 945bf378f5bfc13de4e9097860b5116296c265e9..30df2b076e87f17222b871ab9c46167305024400 100644 (file)
--- a/include/deal.II/fe/fe_poly.templates.h
+++ b/include/deal.II/fe/fe_poly.templates.h
@@ -262,7 +262,7 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &,
       mapping.transform (make_array_view(fe_data.shape_gradients[k]),
                          mapping_covariant,
                          mapping_internal,
-                         make_array_view(output_data.shape_gradients[k]));
+                         make_array_view(output_data.shape_gradients, k));
 
   if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
     {
@@ -270,7 +270,7 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &,
         mapping.transform (make_array_view(fe_data.shape_hessians[k]),
                            mapping_covariant_gradient,
                            mapping_internal,
-                           make_array_view(output_data.shape_hessians[k]));
+                           make_array_view(output_data.shape_hessians, k));
 
       for (unsigned int k=0; k<this->dofs_per_cell; ++k)
         for (unsigned int i=0; i<quadrature.size(); ++i)
@@ -286,7 +286,7 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &,
         mapping.transform (make_array_view(fe_data.shape_3rd_derivatives[k]),
                            mapping_covariant_hessian,
                            mapping_internal,
-                           make_array_view(output_data.shape_3rd_derivatives[k]));
+                           make_array_view(output_data.shape_3rd_derivatives, k));
 
       for (unsigned int k=0; k<this->dofs_per_cell; ++k)
         correct_third_derivatives(output_data, mapping_data, quadrature.size(), k);
@@ -340,7 +340,7 @@ fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator
       mapping.transform (ArrayView<const Tensor<1,dim> >(&fe_data.shape_gradients[k][offset], quadrature.size()),
                          mapping_covariant,
                          mapping_internal,
-                         make_array_view(output_data.shape_gradients[k]));
+                         make_array_view(output_data.shape_gradients, k));
 
   if (flags & update_hessians)
     {
@@ -349,7 +349,7 @@ fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator
                                                            quadrature.size()),
                            mapping_covariant_gradient,
                            mapping_internal,
-                           make_array_view(output_data.shape_hessians[k]));
+                           make_array_view(output_data.shape_hessians, k));
 
       for (unsigned int k=0; k<this->dofs_per_cell; ++k)
         for (unsigned int i=0; i<quadrature.size(); ++i)
@@ -366,7 +366,7 @@ fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator
                                                            quadrature.size()),
                            mapping_covariant_hessian,
                            mapping_internal,
-                           make_array_view(output_data.shape_3rd_derivatives[k]));
+                           make_array_view(output_data.shape_3rd_derivatives, k));
 
       for (unsigned int k=0; k<this->dofs_per_cell; ++k)
         correct_third_derivatives(output_data, mapping_data, quadrature.size(), k);
@@ -423,7 +423,7 @@ fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterato
                                                          quadrature.size()),
                          mapping_covariant,
                          mapping_internal,
-                         make_array_view(output_data.shape_gradients[k]));
+                         make_array_view(output_data.shape_gradients, k));
 
   if (flags & update_hessians)
     {
@@ -432,7 +432,7 @@ fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterato
                                                            quadrature.size()),
                            mapping_covariant_gradient,
                            mapping_internal,
-                           make_array_view(output_data.shape_hessians[k]));
+                           make_array_view(output_data.shape_hessians, k));
 
       for (unsigned int k=0; k<this->dofs_per_cell; ++k)
         for (unsigned int i=0; i<quadrature.size(); ++i)
@@ -449,7 +449,7 @@ fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterato
                                                            quadrature.size()),
                            mapping_covariant_hessian,
                            mapping_internal,
-                           make_array_view(output_data.shape_3rd_derivatives[k]));
+                           make_array_view(output_data.shape_3rd_derivatives, k));
 
       for (unsigned int k=0; k<this->dofs_per_cell; ++k)
         correct_third_derivatives(output_data, mapping_data, quadrature.size(), k);

diff --git a/include/deal.II/fe/fe_update_flags.h b/include/deal.II/fe/fe_update_flags.h
index 84d473ddfffa89c682ecadbe2ec7fb3636cf2cf0..cd322cbfa041e3d2ed25c87e07137506ec43f22a 100644 (file)
--- a/include/deal.II/fe/fe_update_flags.h
+++ b/include/deal.II/fe/fe_update_flags.h
@@ -528,17 +528,17 @@ namespace internal
        * Storage type for gradients. The layout of data is the same as for the
        * #ShapeVector data type.
        */
-      typedef std::vector<std::vector<Tensor<1,spacedim> > > GradientVector;
+      typedef dealii::Table<2,Tensor<1,spacedim> > GradientVector;
 
       /**
        * Likewise for second order derivatives.
        */
-      typedef std::vector<std::vector<Tensor<2,spacedim> > > HessianVector;
+      typedef dealii::Table<2,Tensor<2,spacedim> > HessianVector;
 
       /**
        * And the same also applies to the third order derivatives.
        */
-      typedef std::vector<std::vector<Tensor<3,spacedim> > > ThirdDerivativeVector;
+      typedef dealii::Table<2,Tensor<3,spacedim> > ThirdDerivativeVector;
 
       /**
        * Store the values of the shape functions at the quadrature points. See the

diff --git a/include/deal.II/fe/fe_values.h b/include/deal.II/fe/fe_values.h
index 02fcf605200767cbebb090710551bebe7b556488..495822167b354cf4390967b79b36f16202c99899 100644 (file)
--- a/include/deal.II/fe/fe_values.h
+++ b/include/deal.II/fe/fe_values.h
@@ -4140,10 +4140,6 @@ FEValuesBase<dim,spacedim>::shape_grad (const unsigned int i,
           ExcAccessToUninitializedField("update_gradients"));
   Assert (fe->is_primitive (i),
           ExcShapeFunctionNotPrimitive(i));
-  Assert (i<this->finite_element_output.shape_gradients.size(),
-          ExcIndexRange (i, 0, this->finite_element_output.shape_gradients.size()));
-  Assert (j<this->finite_element_output.shape_gradients[0].size(),
-          ExcIndexRange (j, 0, this->finite_element_output.shape_gradients[0].size()));
 
   // if the entire FE is primitive,
   // then we can take a short-cut:
@@ -4209,10 +4205,6 @@ FEValuesBase<dim,spacedim>::shape_hessian (const unsigned int i,
           ExcAccessToUninitializedField("update_hessians"));
   Assert (fe->is_primitive (i),
           ExcShapeFunctionNotPrimitive(i));
-  Assert (i<this->finite_element_output.shape_hessians.size(),
-          ExcIndexRange (i, 0, this->finite_element_output.shape_hessians.size()));
-  Assert (j<this->finite_element_output.shape_hessians[0].size(),
-          ExcIndexRange (j, 0, this->finite_element_output.shape_hessians[0].size()));
 
   // if the entire FE is primitive,
   // then we can take a short-cut:
@@ -4278,10 +4270,6 @@ FEValuesBase<dim,spacedim>::shape_3rd_derivative (const unsigned int i,
           ExcAccessToUninitializedField("update_3rd_derivatives"));
   Assert (fe->is_primitive (i),
           ExcShapeFunctionNotPrimitive(i));
-  Assert (i<this->finite_element_output.shape_3rd_derivatives.size(),
-          ExcIndexRange (i, 0, this->finite_element_output.shape_3rd_derivatives.size()));
-  Assert (j<this->finite_element_output.shape_3rd_derivatives[0].size(),
-          ExcIndexRange (j, 0, this->finite_element_output.shape_3rd_derivatives[0].size()));
 
   // if the entire FE is primitive,
   // then we can take a short-cut:

diff --git a/source/fe/fe_poly.cc b/source/fe/fe_poly.cc
index c3c4c091a0587a68bf8a6f456f9e426254a64b30..e9c78aedcb38cd29dc40b8d9fba0c687fc95f0a6 100644 (file)
--- a/source/fe/fe_poly.cc
+++ b/source/fe/fe_poly.cc
@@ -55,7 +55,7 @@ fill_fe_values (const Triangulation<1,2>::cell_iterator &,
       mapping.transform (make_array_view(fe_data.shape_gradients[k]),
                          mapping_covariant,
                          mapping_internal,
-                         make_array_view(output_data.shape_gradients[k]));
+                         make_array_view(output_data.shape_gradients, k));
 
   if (fe_data.update_each & update_hessians && cell_similarity != CellSimilarity::translation)
     {
@@ -63,7 +63,7 @@ fill_fe_values (const Triangulation<1,2>::cell_iterator &,
         mapping.transform (make_array_view(fe_data.shape_hessians[k]),
                            mapping_covariant_gradient,
                            mapping_internal,
-                           make_array_view(output_data.shape_hessians[k]));
+                           make_array_view(output_data.shape_hessians, k));
 
       for (unsigned int k=0; k<this->dofs_per_cell; ++k)
         for (unsigned int i=0; i<quadrature.size(); ++i)
@@ -79,7 +79,7 @@ fill_fe_values (const Triangulation<1,2>::cell_iterator &,
         mapping.transform (make_array_view(fe_data.shape_3rd_derivatives[k]),
                            mapping_covariant_hessian,
                            mapping_internal,
-                           make_array_view(output_data.shape_3rd_derivatives[k]));
+                           make_array_view(output_data.shape_3rd_derivatives, k));
 
       for (unsigned int k=0; k<this->dofs_per_cell; ++k)
         correct_third_derivatives(output_data, mapping_data, quadrature.size(), k);
@@ -114,7 +114,7 @@ fill_fe_values (const Triangulation<2,3>::cell_iterator &,
       mapping.transform (make_array_view(fe_data.shape_gradients[k]),
                          mapping_covariant,
                          mapping_internal,
-                         make_array_view(output_data.shape_gradients[k]));
+                         make_array_view(output_data.shape_gradients, k));
 
   if (fe_data.update_each & update_hessians && cell_similarity != CellSimilarity::translation)
     {
@@ -122,7 +122,7 @@ fill_fe_values (const Triangulation<2,3>::cell_iterator &,
         mapping.transform (make_array_view(fe_data.shape_hessians[k]),
                            mapping_covariant_gradient,
                            mapping_internal,
-                           make_array_view(output_data.shape_hessians[k]));
+                           make_array_view(output_data.shape_hessians, k));
 
       for (unsigned int k=0; k<this->dofs_per_cell; ++k)
         for (unsigned int i=0; i<quadrature.size(); ++i)
@@ -138,7 +138,7 @@ fill_fe_values (const Triangulation<2,3>::cell_iterator &,
         mapping.transform (make_array_view(fe_data.shape_3rd_derivatives[k]),
                            mapping_covariant_hessian,
                            mapping_internal,
-                           make_array_view(output_data.shape_3rd_derivatives[k]));
+                           make_array_view(output_data.shape_3rd_derivatives, k));
 
       for (unsigned int k=0; k<this->dofs_per_cell; ++k)
         correct_third_derivatives(output_data, mapping_data, quadrature.size(), k);
@@ -174,7 +174,7 @@ fill_fe_values (const Triangulation<1,2>::cell_iterator &,
       mapping.transform (make_array_view(fe_data.shape_gradients[k]),
                          mapping_covariant,
                          mapping_internal,
-                         make_array_view(output_data.shape_gradients[k]));
+                         make_array_view(output_data.shape_gradients, k));
 
   if (fe_data.update_each & update_hessians && cell_similarity != CellSimilarity::translation)
     {
@@ -182,7 +182,7 @@ fill_fe_values (const Triangulation<1,2>::cell_iterator &,
         mapping.transform (make_array_view(fe_data.shape_hessians[k]),
                            mapping_covariant_gradient,
                            mapping_internal,
-                           make_array_view(output_data.shape_hessians[k]));
+                           make_array_view(output_data.shape_hessians, k));
 
       for (unsigned int k=0; k<this->dofs_per_cell; ++k)
         for (unsigned int i=0; i<quadrature.size(); ++i)
@@ -198,7 +198,7 @@ fill_fe_values (const Triangulation<1,2>::cell_iterator &,
         mapping.transform (make_array_view(fe_data.shape_3rd_derivatives[k]),
                            mapping_covariant_hessian,
                            mapping_internal,
-                           make_array_view(output_data.shape_3rd_derivatives[k]));
+                           make_array_view(output_data.shape_3rd_derivatives, k));
 
       for (unsigned int k=0; k<this->dofs_per_cell; ++k)
         correct_third_derivatives(output_data, mapping_data, quadrature.size(), k);
@@ -229,7 +229,7 @@ fill_fe_values (const Triangulation<2,3>::cell_iterator &,
       mapping.transform (make_array_view(fe_data.shape_gradients[k]),
                          mapping_covariant,
                          mapping_internal,
-                         make_array_view(output_data.shape_gradients[k]));
+                         make_array_view(output_data.shape_gradients, k));
 
   if (fe_data.update_each & update_hessians && cell_similarity != CellSimilarity::translation)
     {
@@ -237,7 +237,7 @@ fill_fe_values (const Triangulation<2,3>::cell_iterator &,
         mapping.transform (make_array_view(fe_data.shape_hessians[k]),
                            mapping_covariant_gradient,
                            mapping_internal,
-                           make_array_view(output_data.shape_hessians[k]));
+                           make_array_view(output_data.shape_hessians, k));
 
       for (unsigned int k=0; k<this->dofs_per_cell; ++k)
         for (unsigned int i=0; i<quadrature.size(); ++i)
@@ -253,7 +253,7 @@ fill_fe_values (const Triangulation<2,3>::cell_iterator &,
         mapping.transform (make_array_view(fe_data.shape_3rd_derivatives[k]),
                            mapping_covariant_hessian,
                            mapping_internal,
-                           make_array_view(output_data.shape_3rd_derivatives[k]));
+                           make_array_view(output_data.shape_3rd_derivatives, k));
 
       for (unsigned int k=0; k<this->dofs_per_cell; ++k)
         correct_third_derivatives(output_data, mapping_data, quadrature.size(), k);

diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc
index 17165b72266216414818121cf656785f13c66031..4392910aea41c75552aec34a0338bfa5a1f664b8 100644 (file)
--- a/source/fe/fe_values.cc
+++ b/source/fe/fe_values.cc
@@ -487,7 +487,7 @@ namespace FEValuesViews
     template <int order, int dim, int spacedim, typename Number>
     void
     do_function_derivatives (const ::dealii::Vector<Number> &dof_values,
-                             const std::vector<std::vector<dealii::Tensor<order,spacedim> > > &shape_derivatives,
+                             const Table<2,dealii::Tensor<order,spacedim> > &shape_derivatives,
                              const std::vector<typename Scalar<dim,spacedim>::ShapeFunctionData> &shape_function_data,
                              std::vector<typename ProductType<Number,dealii::Tensor<order,spacedim> >::type> &derivatives)
     {
@@ -520,7 +520,7 @@ namespace FEValuesViews
     template <int dim, int spacedim, typename Number>
     void
     do_function_laplacians (const ::dealii::Vector<Number> &dof_values,
-                            const std::vector<std::vector<dealii::Tensor<2,spacedim> > > &shape_hessians,
+                            const Table<2,dealii::Tensor<2,spacedim> > &shape_hessians,
                             const std::vector<typename Scalar<dim,spacedim>::ShapeFunctionData> &shape_function_data,
                             std::vector<typename ProductType<Number,double>::type>           &laplacians)
     {
@@ -601,7 +601,7 @@ namespace FEValuesViews
     template <int order, int dim, int spacedim, typename Number>
     void
     do_function_derivatives (const ::dealii::Vector<Number> &dof_values,
-                             const std::vector<std::vector<dealii::Tensor<order,spacedim> > > &shape_derivatives,
+                             const Table<2,dealii::Tensor<order,spacedim> > &shape_derivatives,
                              const std::vector<typename Vector<dim,spacedim>::ShapeFunctionData> &shape_function_data,
                              std::vector<typename ProductType<Number,dealii::Tensor<order+1,spacedim> >::type> &derivatives)
     {
@@ -655,7 +655,7 @@ namespace FEValuesViews
     template <int dim, int spacedim, typename Number>
     void
     do_function_symmetric_gradients (const ::dealii::Vector<Number> &dof_values,
-                                     const std::vector<std::vector<dealii::Tensor<1,spacedim> > > &shape_gradients,
+                                     const Table<2,dealii::Tensor<1,spacedim> > &shape_gradients,
                                      const std::vector<typename Vector<dim,spacedim>::ShapeFunctionData> &shape_function_data,
                                      std::vector<typename ProductType<Number,dealii::SymmetricTensor<2,spacedim> >::type> &symmetric_gradients)
     {
@@ -708,7 +708,7 @@ namespace FEValuesViews
     template <int dim, int spacedim, typename Number>
     void
     do_function_divergences (const ::dealii::Vector<Number> &dof_values,
-                             const std::vector<std::vector<dealii::Tensor<1,spacedim> > > &shape_gradients,
+                             const Table<2,dealii::Tensor<1,spacedim> > &shape_gradients,
                              const std::vector<typename Vector<dim,spacedim>::ShapeFunctionData> &shape_function_data,
                              std::vector<typename ProductType<Number,double>::type> &divergences)
     {
@@ -758,7 +758,7 @@ namespace FEValuesViews
     template <int dim, int spacedim, typename Number>
     void
     do_function_curls (const ::dealii::Vector<Number> &dof_values,
-                       const std::vector<std::vector<dealii::Tensor<1,spacedim> > > &shape_gradients,
+                       const Table<2,dealii::Tensor<1,spacedim> > &shape_gradients,
                        const std::vector<typename Vector<dim,spacedim>::ShapeFunctionData> &shape_function_data,
                        std::vector<typename ProductType<Number,typename dealii::internal::CurlType<spacedim>::type>::type> &curls)
     {
@@ -950,7 +950,7 @@ namespace FEValuesViews
     template <int dim, int spacedim, typename Number>
     void
     do_function_laplacians (const ::dealii::Vector<Number> &dof_values,
-                            const std::vector<std::vector<dealii::Tensor<2,spacedim> > > &shape_hessians,
+                            const Table<2,dealii::Tensor<2,spacedim> > &shape_hessians,
                             const std::vector<typename Vector<dim,spacedim>::ShapeFunctionData> &shape_function_data,
                             std::vector<typename ProductType<Number,dealii::Tensor<1,spacedim> >::type> &laplacians)
     {
@@ -1058,7 +1058,7 @@ namespace FEValuesViews
     template <int dim, int spacedim, typename Number>
     void
     do_function_divergences (const ::dealii::Vector<Number> &dof_values,
-                             const std::vector<std::vector<dealii::Tensor<1,spacedim> > > &shape_gradients,
+                             const Table<2,dealii::Tensor<1,spacedim> > &shape_gradients,
                              const std::vector<typename SymmetricTensor<2,dim,spacedim>::ShapeFunctionData> &shape_function_data,
                              std::vector<typename ProductType<Number,dealii::Tensor<1,spacedim> >::type> &divergences)
     {
@@ -1202,7 +1202,7 @@ namespace FEValuesViews
     template <int dim, int spacedim, typename Number>
     void
     do_function_divergences (const ::dealii::Vector<Number> &dof_values,
-                             const std::vector<std::vector<dealii::Tensor<1,spacedim> > > &shape_gradients,
+                             const Table<2,dealii::Tensor<1,spacedim> > &shape_gradients,
                              const std::vector<typename Tensor<2,dim,spacedim>::ShapeFunctionData> &shape_function_data,
                              std::vector<typename ProductType<Number,dealii::Tensor<1,spacedim> >::type> &divergences)
     {
@@ -2230,19 +2230,25 @@ namespace internal
         }
 
       if (flags & update_gradients)
-        this->shape_gradients.resize (n_nonzero_shape_components,
-                                      std::vector<Tensor<1,spacedim> > (n_quadrature_points,
-                                          numbers::signaling_nan<Tensor<1,spacedim> >()));
+        {
+          this->shape_gradients.reinit(n_nonzero_shape_components,
+                                       n_quadrature_points);
+          this->shape_gradients.fill (numbers::signaling_nan<Tensor<1,spacedim> >());
+        }
 
       if (flags & update_hessians)
-        this->shape_hessians.resize (n_nonzero_shape_components,
-                                     std::vector<Tensor<2,spacedim> > (n_quadrature_points,
-                                         numbers::signaling_nan<Tensor<2,spacedim> >()));
+        {
+          this->shape_hessians.reinit(n_nonzero_shape_components,
+                                      n_quadrature_points);
+          this->shape_hessians.fill (numbers::signaling_nan<Tensor<2,spacedim> >());
+        }
 
       if (flags & update_3rd_derivatives)
-        this->shape_3rd_derivatives.resize (n_nonzero_shape_components,
-                                            std::vector<Tensor<3,spacedim> > (n_quadrature_points,
-                                                numbers::signaling_nan<Tensor<3,spacedim> >()));
+        {
+          this->shape_3rd_derivatives.reinit(n_nonzero_shape_components,
+                                             n_quadrature_points);
+          this->shape_3rd_derivatives.fill (numbers::signaling_nan<Tensor<3,spacedim> >());
+        }
     }
 
 
@@ -2440,10 +2446,10 @@ namespace internal
   template <int order, int spacedim, typename Number>
   void
   do_function_derivatives (const Number                     *dof_values_ptr,
-                           const std::vector<std::vector<Tensor<order,spacedim> > > &shape_derivatives,
+                           const dealii::Table<2,Tensor<order,spacedim> > &shape_derivatives,
                            std::vector<Tensor<order,spacedim,Number> > &derivatives)
   {
-    const unsigned int dofs_per_cell = shape_derivatives.size();
+    const unsigned int dofs_per_cell = shape_derivatives.size()[0];
     const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
                                              shape_derivatives[0].size() : derivatives.size();
     AssertDimension(derivatives.size(), n_quadrature_points);
@@ -2475,7 +2481,7 @@ namespace internal
   template <int order, int dim, int spacedim, typename Number>
   void
   do_function_derivatives (const Number                      *dof_values_ptr,
-                           const std::vector<std::vector<Tensor<order,spacedim> > > &shape_derivatives,
+                           const dealii::Table<2,Tensor<order,spacedim> > &shape_derivatives,
                            const FiniteElement<dim,spacedim> &fe,
                            const std::vector<unsigned int> &shape_function_to_row_table,
                            VectorSlice<std::vector<std::vector<Tensor<order,spacedim,Number> > > > &derivatives,
@@ -2571,10 +2577,10 @@ namespace internal
   template <int spacedim, typename Number, typename Number2>
   void
   do_function_laplacians (const Number2        *dof_values_ptr,
-                          const std::vector<std::vector<Tensor<2,spacedim> > > &shape_hessians,
+                          const dealii::Table<2,Tensor<2,spacedim> > &shape_hessians,
                           std::vector<Number> &laplacians)
   {
-    const unsigned int dofs_per_cell = shape_hessians.size();
+    const unsigned int dofs_per_cell = shape_hessians.size()[0];
     const unsigned int n_quadrature_points = dofs_per_cell > 0 ?
                                              shape_hessians[0].size() : laplacians.size();
     AssertDimension(laplacians.size(), n_quadrature_points);
@@ -2601,7 +2607,7 @@ namespace internal
   template <int dim, int spacedim, typename VectorType, typename Number>
   void
   do_function_laplacians (const Number                    *dof_values_ptr,
-                          const std::vector<std::vector<Tensor<2,spacedim> > > &shape_hessians,
+                          const dealii::Table<2,Tensor<2,spacedim> > &shape_hessians,
                           const FiniteElement<dim,spacedim> &fe,
                           const std::vector<unsigned int> &shape_function_to_row_table,
                           std::vector<VectorType>         &laplacians,