Clean up some functions in FEFieldFunction.

diff --git a/include/deal.II/numerics/fe_field_function.templates.h b/include/deal.II/numerics/fe_field_function.templates.h
index b80d1bffd92ea069cebda3c802f9e726ea8a603c..0f3b8f614e40d13d19bc9a328f4db0875d5dd19f 100644 (file)
--- a/include/deal.II/numerics/fe_field_function.templates.h
+++ b/include/deal.II/numerics/fe_field_function.templates.h
@@ -273,10 +273,9 @@ namespace Functions
 
     const unsigned int n_cells =
       compute_point_locations(points, cells, qpoints, maps);
-    hp::MappingCollection<dim>   mapping_collection(mapping);
-    const hp::FECollection<dim> &fe_collection = dh->get_fe_collection();
-    hp::QCollection<dim>         quadrature_collection;
+
     // Create quadrature collection
+    hp::QCollection<dim> quadrature_collection;
     for (unsigned int i = 0; i < n_cells; ++i)
       {
         // Number of quadrature points on this cell
@@ -286,21 +285,25 @@ namespace Functions
 
         quadrature_collection.push_back(Quadrature<dim>(qpoints[i], ww));
       }
-    // Get a function value object
-    hp::FEValues<dim> fe_v(mapping_collection,
-                           fe_collection,
+
+    // Now gather all the information we need
+    hp::MappingCollection<dim> mapping_collection(mapping);
+    hp::FEValues<dim>          fe_v(mapping_collection,
+                           dh->get_fe_collection(),
                            quadrature_collection,
                            update_values);
-    // Now gather all the information we need
     for (unsigned int i = 0; i < n_cells; ++i)
       {
         AssertThrow(!cells[i]->is_artificial(),
                     VectorTools::ExcPointNotAvailableHere());
+
         fe_v.reinit(cells[i], i, 0);
         const unsigned int nq = qpoints[i].size();
+
         std::vector<Vector<typename VectorType::value_type>> vvalues(
           nq, Vector<typename VectorType::value_type>(this->n_components));
         fe_v.get_present_fe_values().get_function_values(data_vector, vvalues);
+
         for (unsigned int q = 0; q < nq; ++q)
           values[maps[i][q]] = vvalues[q];
       }
@@ -317,10 +320,17 @@ namespace Functions
   {
     Assert(points.size() == values.size(),
            ExcDimensionMismatch(points.size(), values.size()));
+
+    // Simply forward everything to the vector_value_list()
+    // function. This requires a temporary object, but everything we
+    // do here is so expensive that that really doesn't make any
+    // difference any more.
     std::vector<Vector<typename VectorType::value_type>> vvalues(
       points.size(),
       Vector<typename VectorType::value_type>(this->n_components));
+
     vector_value_list(points, vvalues);
+
     for (unsigned int q = 0; q < points.size(); ++q)
       values[q] = vvalues[q](component);
   }
@@ -343,10 +353,9 @@ namespace Functions
 
     const unsigned int n_cells =
       compute_point_locations(points, cells, qpoints, maps);
-    hp::MappingCollection<dim>   mapping_collection(mapping);
-    const hp::FECollection<dim> &fe_collection = dh->get_fe_collection();
-    hp::QCollection<dim>         quadrature_collection;
+
     // Create quadrature collection
+    hp::QCollection<dim> quadrature_collection;
     for (unsigned int i = 0; i < n_cells; ++i)
       {
         // Number of quadrature points on this cell
@@ -356,17 +365,20 @@ namespace Functions
 
         quadrature_collection.push_back(Quadrature<dim>(qpoints[i], ww));
       }
-    // Get a function value object
-    hp::FEValues<dim> fe_v(mapping_collection,
-                           fe_collection,
+
+    // Now gather all the information we need
+    hp::MappingCollection<dim> mapping_collection(mapping);
+    hp::FEValues<dim>          fe_v(mapping_collection,
+                           dh->get_fe_collection(),
                            quadrature_collection,
                            update_gradients);
-    // Now gather all the information we need
     for (unsigned int i = 0; i < n_cells; ++i)
       {
         AssertThrow(!cells[i]->is_artificial(),
                     VectorTools::ExcPointNotAvailableHere());
+
         fe_v.reinit(cells[i], i, 0);
+
         const unsigned int nq = qpoints[i].size();
         std::vector<
           std::vector<Tensor<1, dim, typename VectorType::value_type>>>
@@ -375,6 +387,7 @@ namespace Functions
                    this->n_components));
         fe_v.get_present_fe_values().get_function_gradients(data_vector,
                                                             vgrads);
+
         for (unsigned int q = 0; q < nq; ++q)
           {
             const unsigned int s = vgrads[q].size();
@@ -385,6 +398,8 @@ namespace Functions
       }
   }
 
+
+
   template <int dim, typename DoFHandlerType, typename VectorType>
   void
   FEFieldFunction<dim, DoFHandlerType, VectorType>::gradient_list(
@@ -394,16 +409,24 @@ namespace Functions
   {
     Assert(points.size() == values.size(),
            ExcDimensionMismatch(points.size(), values.size()));
+
+    // Simply forward everything to the vector_gradient_list()
+    // function. This requires a temporary object, but everything we
+    // do here is so expensive that that really doesn't make any
+    // difference any more.
     std::vector<std::vector<Tensor<1, dim, typename VectorType::value_type>>>
       vvalues(points.size(),
               std::vector<Tensor<1, dim, typename VectorType::value_type>>(
                 this->n_components));
+
     vector_gradient_list(points, vvalues);
+
     for (unsigned int q = 0; q < points.size(); ++q)
       values[q] = vvalues[q][component];
   }
 
 
+
   template <int dim, typename DoFHandlerType, typename VectorType>
   void
   FEFieldFunction<dim, DoFHandlerType, VectorType>::vector_laplacian_list(
@@ -419,10 +442,9 @@ namespace Functions
 
     const unsigned int n_cells =
       compute_point_locations(points, cells, qpoints, maps);
-    hp::MappingCollection<dim>   mapping_collection(mapping);
-    const hp::FECollection<dim> &fe_collection = dh->get_fe_collection();
-    hp::QCollection<dim>         quadrature_collection;
+
     // Create quadrature collection
+    hp::QCollection<dim> quadrature_collection;
     for (unsigned int i = 0; i < n_cells; ++i)
       {
         // Number of quadrature points on this cell
@@ -432,9 +454,11 @@ namespace Functions
 
         quadrature_collection.push_back(Quadrature<dim>(qpoints[i], ww));
       }
-    // Get a function value object
-    hp::FEValues<dim> fe_v(mapping_collection,
-                           fe_collection,
+
+    // Now gather all the information we need
+    hp::MappingCollection<dim> mapping_collection(mapping);
+    hp::FEValues<dim>          fe_v(mapping_collection,
+                           dh->get_fe_collection(),
                            quadrature_collection,
                            update_hessians);
     // Now gather all the information we need
@@ -442,17 +466,22 @@ namespace Functions
       {
         AssertThrow(!cells[i]->is_artificial(),
                     VectorTools::ExcPointNotAvailableHere());
+
         fe_v.reinit(cells[i], i, 0);
+
         const unsigned int nq = qpoints[i].size();
         std::vector<Vector<typename VectorType::value_type>> vvalues(
           nq, Vector<typename VectorType::value_type>(this->n_components));
         fe_v.get_present_fe_values().get_function_laplacians(data_vector,
                                                              vvalues);
+
         for (unsigned int q = 0; q < nq; ++q)
           values[maps[i][q]] = vvalues[q];
       }
   }
 
+
+
   template <int dim, typename DoFHandlerType, typename VectorType>
   void
   FEFieldFunction<dim, DoFHandlerType, VectorType>::laplacian_list(
@@ -462,10 +491,17 @@ namespace Functions
   {
     Assert(points.size() == values.size(),
            ExcDimensionMismatch(points.size(), values.size()));
+
+    // Simply forward everything to the vector_gradient_list()
+    // function. This requires a temporary object, but everything we
+    // do here is so expensive that that really doesn't make any
+    // difference any more.
     std::vector<Vector<typename VectorType::value_type>> vvalues(
       points.size(),
       Vector<typename VectorType::value_type>(this->n_components));
+
     vector_laplacian_list(points, vvalues);
+
     for (unsigned int q = 0; q < points.size(); ++q)
       values[q] = vvalues[q](component);
   }