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.

diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h
index 06d4beff9b77d2773688e67d052e0efd68f6fc38..0bd58ab14721576e70a452bbd06e7121cbe1e3e2 100644 (file)
--- 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 158437471be32b86efb94e954e9550a4d4fd790b..0a957cc02d09620bab93d4cc33c4ed1c929c046e 100644 (file)
--- 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 ca19ad1b8e98d55faeefe45f19df7100302195d7..96ab2065de8837a141fd6cb5070edda7765deef2 100644 (file)
--- 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;
 }