From 9ed4f2c8ae505bc13fb39c1fd2303f237c09eb45 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Thu, 3 Sep 2015 10:00:04 -0500
Subject: [PATCH] Move the computation of the Mapping values into initialize().

Rather than having to do it yourself right after initialize(), let initialize()
call the function that does that. This greatly simplifies a good deal of code,
and also leads to better parallelization.
---
 include/deal.II/fe/mapping_q1.h | 18 +++++++++++++++---
 source/fe/mapping_q.cc          | 17 -----------------
 source/fe/mapping_q1.cc         |  8 ++++----
 3 files changed, 19 insertions(+), 24 deletions(-)

diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h
index 06d4beff9b..0bd58ab147 100644
--- a/include/deal.II/fe/mapping_q1.h
+++ b/include/deal.II/fe/mapping_q1.h
@@ -169,6 +169,10 @@ public:
     /**
      * Initialize the object's member variables related to cell data
      * based on the given arguments.
+     *
+     * The function also calls compute_shape_function_values() to
+     * actually set the member variables related to the values and
+     * derivatives of the mapping shape functions.
      */
     void
     initialize (const UpdateFlags      update_flags,
@@ -177,7 +181,8 @@ public:
 
     /**
      * Initialize the object's member variables related to cell and
-     * face data based on the given arguments.
+     * face data based on the given arguments. In order to initialize
+     * cell data, this function calls initialize().
      */
     void
     initialize_face (const UpdateFlags      update_flags,
@@ -190,8 +195,15 @@ public:
      *
      * Which values, derivatives, or higher order derivatives are
      * computed is determined by which of the member arrays have
-     * nonzero sizes. They is typically set to their appropriate sizes
-     * by the initialize() and initialize_face() functions.
+     * nonzero sizes. They are typically set to their appropriate
+     * sizes by the initialize() and initialize_face() functions,
+     * which indeed call this function internally. However, it is
+     * possible (and at times useful) to do the resizing by hand and
+     * then call this function directly. An example is in a Newton
+     * iteration where we update the location of a quadrature point
+     * (e.g., in MappingQ::transform_real_to_uni_cell()) and need to
+     * re-compute the mapping and its derivatives at this location,
+     * but have already sized all internal arrays correctly.
      */
     void compute_shape_function_values (const std::vector<Point<dim> > &unit_points);
 
diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc
index 158437471b..0a957cc02d 100644
--- a/source/fe/mapping_q.cc
+++ b/source/fe/mapping_q.cc
@@ -125,12 +125,6 @@ MappingQ<dim,spacedim>::get_data (const UpdateFlags update_flags,
                                                  quadrature.size())));
   tasks.join_all ();
 
-  // TODO: parallelize this as well
-  data->compute_shape_function_values (quadrature.get_points());
-  if (!use_mapping_q_on_all_cells)
-    data->mapping_q1_data.compute_shape_function_values (quadrature.get_points());
-
-
   return data;
 }
 
@@ -161,11 +155,6 @@ MappingQ<dim,spacedim>::get_face_data (const UpdateFlags update_flags,
                                                  quadrature.size())));
   tasks.join_all ();
 
-  // TODO: parallelize this as well
-  data->compute_shape_function_values (q.get_points());
-  if (!use_mapping_q_on_all_cells)
-    data->mapping_q1_data.compute_shape_function_values (q.get_points());
-
   return data;
 }
 
@@ -196,12 +185,6 @@ MappingQ<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
                                                  quadrature.size())));
   tasks.join_all ();
 
-  // TODO: parallelize this as well
-  data->compute_shape_function_values (q.get_points());
-  if (!use_mapping_q_on_all_cells)
-    data->mapping_q1_data.compute_shape_function_values (q.get_points());
-
-
   return data;
 }
 
diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc
index ca19ad1b8e..96ab2065de 100644
--- a/source/fe/mapping_q1.cc
+++ b/source/fe/mapping_q1.cc
@@ -124,9 +124,13 @@ initialize (const UpdateFlags      update_flags,
   if (this->update_each &
       (update_jacobian_3rd_derivatives | update_jacobian_pushed_forward_3rd_derivatives) )
     shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
+
+  // now also fill the various fields with their correct values
+  compute_shape_function_values (q.get_points());
 }
 
 
+
 template <int dim, int spacedim>
 void
 MappingQ1<dim,spacedim>::InternalData::
@@ -843,7 +847,6 @@ MappingQ1<dim,spacedim>::get_data (const UpdateFlags update_flags,
 {
   InternalData *data = new InternalData(1);
   data->initialize (requires_update_flags(update_flags), q, q.size());
-  data->compute_shape_function_values (q.get_points());
 
   return data;
 }
@@ -859,7 +862,6 @@ MappingQ1<dim,spacedim>::get_face_data (const UpdateFlags        update_flags,
   data->initialize_face (requires_update_flags(update_flags),
                          QProjector<dim>::project_to_all_faces(quadrature),
                          quadrature.size());
-  data->compute_shape_function_values (QProjector<dim>::project_to_all_faces(quadrature).get_points());
 
   return data;
 }
@@ -875,8 +877,6 @@ MappingQ1<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
   data->initialize_face (requires_update_flags(update_flags),
                          QProjector<dim>::project_to_all_subfaces(quadrature),
                          quadrature.size());
-  data->compute_shape_function_values (QProjector<dim>::project_to_all_subfaces(quadrature).get_points());
-
 
   return data;
 }
-- 
2.39.5