From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Thu, 3 Sep 2015 13:15:26 +0000 (-0500)
Subject: Move the code that initializes the Mapping{Q,Q1}::InternalData into that structure.
X-Git-Tag: v8.4.0-rc2~474^2~9
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0faad4e08019d5f9909e70ebee390198ca978d32;p=dealii.git

Move the code that initializes the Mapping{Q,Q1}::InternalData into that structure.

This allows us to get rid of the two implementations of compute_shapes(). It also
significantly simplifies the logic which of the two functions actually need to
be called, since the InternalData object actually knows that itself.
---

diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h
index f6b588a1ec..2d4914eaaf 100644
--- a/include/deal.II/fe/mapping_q.h
+++ b/include/deal.II/fe/mapping_q.h
@@ -72,17 +72,10 @@ public:
             const bool use_mapping_q_on_all_cells = false);
 
   /**
-   * Copy constructor. Performs a deep copy, i.e. duplicates what #tensor_pols
-   * points to instead of simply copying the #tensor_pols pointer as done by a
-   * default copy constructor.
+   * Copy constructor.
    */
   MappingQ (const MappingQ<dim,spacedim> &mapping);
 
-  /**
-   * Destructor.
-   */
-  virtual ~MappingQ ();
-
   /**
    * Transforms the point @p p on the unit cell to the point @p p_real on the
    * real cell @p cell and returns @p p_real.
@@ -333,10 +326,7 @@ protected:
   /**
    * @}
    */
-public:
-  void
-  compute_shapes (const std::vector<Point<dim> > &unit_points,
-                  typename MappingQ1<dim,spacedim>::InternalData &data) const;
+
 protected:
 
   /**
@@ -462,12 +452,6 @@ protected:
    */
   const unsigned int n_shape_functions;
 
-  /**
-   * Mapping from lexicographic to to the Qp shape function numbering. Its
-   * size is @p dofs_per_cell.
-   */
-  const std::vector<unsigned int> renumber;
-
   /**
    * If this flag is set @p true then @p MappingQ is used on all cells, not
    * only on boundary cells.
diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h
index 6a32f40ceb..06d4beff9b 100644
--- a/include/deal.II/fe/mapping_q1.h
+++ b/include/deal.II/fe/mapping_q1.h
@@ -184,6 +184,18 @@ public:
                      const Quadrature<dim> &quadrature,
                      const unsigned int     n_original_q_points);
 
+    /**
+     * Compute the values and/or derivatives of the shape functions
+     * used for the mapping.
+     *
+     * 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.
+     */
+    void compute_shape_function_values (const std::vector<Point<dim> > &unit_points);
+
+
     /**
      * Shape function at quadrature point. Shape functions are in tensor
      * product order, so vertices must be reordered to obtain transformation.
@@ -465,16 +477,6 @@ protected:
   DataSetDescriptor;
 
 
-  /**
-   * Compute shape values and/or derivatives.
-   *
-   * Calls either the @p compute_shapes_virtual of this class or that of the
-   * derived class, depending on whether <tt>data.polynomial_degree</tt>
-   * is one or is greater than that.
-   */
-  void compute_shapes (const std::vector<Point<dim> > &unit_points,
-                       InternalData &data) const;
-
   /**
    * Transforms a point @p p on the unit cell to the point @p p_real on the
    * real cell @p cell and returns @p p_real.
diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc
index e2d6ea4737..158437471b 100644
--- a/source/fe/mapping_q.cc
+++ b/source/fe/mapping_q.cc
@@ -24,7 +24,6 @@
 #include <deal.II/grid/tria_iterator.h>
 #include <deal.II/grid/tria_boundary.h>
 #include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/fe/fe_tools.h>
 #include <deal.II/fe/mapping_q.h>
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/fe_values.h>
@@ -55,22 +54,6 @@ MappingQ<dim,spacedim>::InternalData::memory_consumption () const
 
 
 
-namespace
-{
-  template <int dim>
-  std::vector<unsigned int>
-  get_dpo_vector (const unsigned int degree)
-  {
-    std::vector<unsigned int> dpo(dim+1, 1U);
-    for (unsigned int i=1; i<dpo.size(); ++i)
-      dpo[i]=dpo[i-1]*(degree-1);
-    return dpo;
-  }
-}
-
-
-
-
 template<int dim, int spacedim>
 MappingQ<dim,spacedim>::MappingQ (const unsigned int p,
                                   const bool use_mapping_q_on_all_cells)
@@ -81,23 +64,12 @@ MappingQ<dim,spacedim>::MappingQ (const unsigned int p,
           ((dim==2) ?
            4+4*(degree-1) :
            8+12*(degree-1)+6*(degree-1)*(degree-1))),
-  tensor_pols(0),
   n_shape_functions(Utilities::fixed_power<dim>(degree+1)),
-  renumber(FETools::
-           lexicographic_to_hierarchic_numbering (
-             FiniteElementData<dim> (get_dpo_vector<dim>(degree), 1,
-                                     degree))),
   use_mapping_q_on_all_cells (use_mapping_q_on_all_cells
                               || (dim != spacedim)),
   feq(degree),
   line_support_points(degree+1)
 {
-  // Construct the tensor product polynomials used as shape functions for the
-  // Qp mapping of cells at the boundary.
-  tensor_pols = new TensorProductPolynomials<dim>
-  (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points()));
-  Assert (n_shape_functions==tensor_pols->n(),
-          ExcInternalError());
   Assert(n_inner+n_outer==n_shape_functions, ExcInternalError());
 
   // build laplace_on_quad_vector
@@ -118,105 +90,16 @@ MappingQ<dim,spacedim>::MappingQ (const MappingQ<dim,spacedim> &mapping)
   degree(mapping.degree),
   n_inner(mapping.n_inner),
   n_outer(mapping.n_outer),
-  tensor_pols(0),
   n_shape_functions(mapping.n_shape_functions),
-  renumber(mapping.renumber),
   use_mapping_q_on_all_cells (mapping.use_mapping_q_on_all_cells),
   feq(degree),
   line_support_points(degree+1)
 {
-  tensor_pols=new TensorProductPolynomials<dim> (*mapping.tensor_pols);
   laplace_on_quad_vector=mapping.laplace_on_quad_vector;
   laplace_on_hex_vector=mapping.laplace_on_hex_vector;
 }
 
 
-template<int dim, int spacedim>
-MappingQ<dim,spacedim>::~MappingQ ()
-{
-  delete tensor_pols;
-}
-
-
-template<int dim, int spacedim>
-void
-MappingQ<dim,spacedim>::compute_shapes (const std::vector<Point<dim> > &unit_points,
-                                        typename MappingQ1<dim,spacedim>::InternalData &data) const
-{
-  const unsigned int n_points=unit_points.size();
-  std::vector<double> values;
-  std::vector<Tensor<1,dim> > grads;
-  if (data.shape_values.size()!=0)
-    {
-      Assert(data.shape_values.size()==n_shape_functions*n_points,
-             ExcInternalError());
-      values.resize(n_shape_functions);
-    }
-  if (data.shape_derivatives.size()!=0)
-    {
-      Assert(data.shape_derivatives.size()==n_shape_functions*n_points,
-             ExcInternalError());
-      grads.resize(n_shape_functions);
-    }
-
-//                                 // dummy variable of size 0
-  std::vector<Tensor<2,dim> > grad2;
-  if (data.shape_second_derivatives.size()!=0)
-    {
-      Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points,
-             ExcInternalError());
-      grad2.resize(n_shape_functions);
-    }
-
-  std::vector<Tensor<3,dim> > grad3;
-  if (data.shape_third_derivatives.size()!=0)
-    {
-      Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points,
-             ExcInternalError());
-      grad3.resize(n_shape_functions);
-    }
-
-  std::vector<Tensor<4,dim> > grad4;
-  if (data.shape_fourth_derivatives.size()!=0)
-    {
-      Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points,
-             ExcInternalError());
-      grad4.resize(n_shape_functions);
-    }
-
-
-  if (data.shape_values.size()!=0 ||
-      data.shape_derivatives.size()!=0 ||
-      data.shape_second_derivatives.size()!=0 ||
-      data.shape_third_derivatives.size()!=0 ||
-      data.shape_fourth_derivatives.size()!=0 )
-    for (unsigned int point=0; point<n_points; ++point)
-      {
-        tensor_pols->compute(unit_points[point], values, grads, grad2, grad3, grad4);
-
-        if (data.shape_values.size()!=0)
-          for (unsigned int i=0; i<n_shape_functions; ++i)
-            data.shape(point,renumber[i]) = values[i];
-
-        if (data.shape_derivatives.size()!=0)
-          for (unsigned int i=0; i<n_shape_functions; ++i)
-            data.derivative(point,renumber[i]) = grads[i];
-
-        if (data.shape_second_derivatives.size()!=0)
-          for (unsigned int i=0; i<n_shape_functions; ++i)
-            data.second_derivative(point,renumber[i]) = grad2[i];
-
-        if (data.shape_third_derivatives.size()!=0)
-          for (unsigned int i=0; i<n_shape_functions; ++i)
-            data.third_derivative(point,renumber[i]) = grad3[i];
-
-        if (data.shape_fourth_derivatives.size()!=0)
-          for (unsigned int i=0; i<n_shape_functions; ++i)
-            data.fourth_derivative(point,renumber[i]) = grad4[i];
-      }
-}
-
-
 
 template<int dim, int spacedim>
 typename MappingQ<dim,spacedim>::InternalData *
@@ -243,9 +126,9 @@ MappingQ<dim,spacedim>::get_data (const UpdateFlags update_flags,
   tasks.join_all ();
 
   // TODO: parallelize this as well
-  this->compute_shapes (quadrature.get_points(), *data);
+  data->compute_shape_function_values (quadrature.get_points());
   if (!use_mapping_q_on_all_cells)
-    this->MappingQ1<dim,spacedim>::compute_shapes (quadrature.get_points(), data->mapping_q1_data);
+    data->mapping_q1_data.compute_shape_function_values (quadrature.get_points());
 
 
   return data;
@@ -279,9 +162,9 @@ MappingQ<dim,spacedim>::get_face_data (const UpdateFlags update_flags,
   tasks.join_all ();
 
   // TODO: parallelize this as well
-  this->compute_shapes (q.get_points(), *data);
+  data->compute_shape_function_values (q.get_points());
   if (!use_mapping_q_on_all_cells)
-    this->MappingQ1<dim,spacedim>::compute_shapes (q.get_points(), data->mapping_q1_data);
+    data->mapping_q1_data.compute_shape_function_values (q.get_points());
 
   return data;
 }
@@ -314,9 +197,9 @@ MappingQ<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
   tasks.join_all ();
 
   // TODO: parallelize this as well
-  this->compute_shapes (q.get_points(), *data);
+  data->compute_shape_function_values (q.get_points());
   if (!use_mapping_q_on_all_cells)
-    this->MappingQ1<dim,spacedim>::compute_shapes (q.get_points(), data->mapping_q1_data);
+    data->mapping_q1_data.compute_shape_function_values (q.get_points());
 
 
   return data;
@@ -616,7 +499,7 @@ MappingQ<dim,spacedim>::compute_laplace_vector(Table<2,double> &lvs) const
 
   InternalData quadrature_data(degree);
   quadrature_data.shape_derivatives.resize(n_shape_functions * n_q_points);
-  this->compute_shapes(quadrature.get_points(), quadrature_data);
+  quadrature_data.compute_shape_function_values(quadrature.get_points());
 
   // Compute the stiffness matrix of the inner dofs
   FullMatrix<long double> S(n_inner);
diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc
index 69ce4b5d8f..ca19ad1b8e 100644
--- a/source/fe/mapping_q1.cc
+++ b/source/fe/mapping_q1.cc
@@ -24,6 +24,8 @@
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/tria_iterator.h>
 #include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/mapping_q1.h>
 #include <deal.II/fe/mapping_q.h>
@@ -326,9 +328,9 @@ namespace internal
 
     template <int spacedim>
     void
-    compute_shapes (const unsigned int            n_shape_functions,
-                    const std::vector<Point<1> > &unit_points,
-                    typename dealii::MappingQ1<1,spacedim>::InternalData &data)
+    compute_shape_function_values (const unsigned int            n_shape_functions,
+                                   const std::vector<Point<1> > &unit_points,
+                                   typename dealii::MappingQ1<1,spacedim>::InternalData &data)
     {
       (void)n_shape_functions;
       const unsigned int n_points=unit_points.size();
@@ -391,9 +393,9 @@ namespace internal
 
     template <int spacedim>
     void
-    compute_shapes (const unsigned int            n_shape_functions,
-                    const std::vector<Point<2> > &unit_points,
-                    typename dealii::MappingQ1<2,spacedim>::InternalData &data)
+    compute_shape_function_values (const unsigned int            n_shape_functions,
+                                   const std::vector<Point<2> > &unit_points,
+                                   typename dealii::MappingQ1<2,spacedim>::InternalData &data)
     {
       (void)n_shape_functions;
       const unsigned int n_points=unit_points.size();
@@ -469,9 +471,9 @@ namespace internal
 
     template <int spacedim>
     void
-    compute_shapes (const unsigned int            n_shape_functions,
-                    const std::vector<Point<3> > &unit_points,
-                    typename dealii::MappingQ1<3,spacedim>::InternalData &data)
+    compute_shape_function_values (const unsigned int            n_shape_functions,
+                                   const std::vector<Point<3> > &unit_points,
+                                   typename dealii::MappingQ1<3,spacedim>::InternalData &data)
     {
       (void)n_shape_functions;
       const unsigned int n_points=unit_points.size();
@@ -653,14 +655,123 @@ namespace internal
 }
 
 
+namespace
+{
+  template <int dim>
+  std::vector<unsigned int>
+  get_dpo_vector (const unsigned int degree)
+  {
+    std::vector<unsigned int> dpo(dim+1, 1U);
+    for (unsigned int i=1; i<dpo.size(); ++i)
+      dpo[i]=dpo[i-1]*(degree-1);
+    return dpo;
+  }
+}
+
+
+
 template<int dim, int spacedim>
 void
-MappingQ1<dim,spacedim>::compute_shapes (const std::vector<Point<dim> > &unit_points,
-                                         InternalData &data) const
+MappingQ1<dim,spacedim>::InternalData::
+compute_shape_function_values (const std::vector<Point<dim> > &unit_points)
 {
-  internal::MappingQ1::
-  compute_shapes<spacedim> (n_shape_functions,
-                            unit_points, data);
+  // if the polynomial degree is one, then we can simplify code a bit
+  // by using hard-coded shape functions.
+  if ((polynomial_degree == 1)
+      &&
+      (dim == spacedim))
+    internal::MappingQ1::compute_shape_function_values<spacedim> (n_shape_functions,
+        unit_points, *this);
+  else
+    // otherwise ask an object that describes the polynomial space
+    {
+      const unsigned int n_points=unit_points.size();
+
+      // Construct the tensor product polynomials used as shape functions for the
+      // Qp mapping of cells at the boundary.
+      const QGaussLobatto<1> line_support_points (polynomial_degree + 1);
+      const TensorProductPolynomials<dim>
+      tensor_pols (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points()));
+      Assert (n_shape_functions==tensor_pols.n(),
+              ExcInternalError());
+
+      // then also construct the mapping from lexicographic to the Qp shape function numbering
+      const std::vector<unsigned int>
+      renumber (FETools::
+                lexicographic_to_hierarchic_numbering (
+                  FiniteElementData<dim> (get_dpo_vector<dim>(polynomial_degree), 1,
+                                          polynomial_degree)));
+
+      std::vector<double> values;
+      std::vector<Tensor<1,dim> > grads;
+      if (shape_values.size()!=0)
+        {
+          Assert(shape_values.size()==n_shape_functions*n_points,
+                 ExcInternalError());
+          values.resize(n_shape_functions);
+        }
+      if (shape_derivatives.size()!=0)
+        {
+          Assert(shape_derivatives.size()==n_shape_functions*n_points,
+                 ExcInternalError());
+          grads.resize(n_shape_functions);
+        }
+
+      std::vector<Tensor<2,dim> > grad2;
+      if (shape_second_derivatives.size()!=0)
+        {
+          Assert(shape_second_derivatives.size()==n_shape_functions*n_points,
+                 ExcInternalError());
+          grad2.resize(n_shape_functions);
+        }
+
+      std::vector<Tensor<3,dim> > grad3;
+      if (shape_third_derivatives.size()!=0)
+        {
+          Assert(shape_third_derivatives.size()==n_shape_functions*n_points,
+                 ExcInternalError());
+          grad3.resize(n_shape_functions);
+        }
+
+      std::vector<Tensor<4,dim> > grad4;
+      if (shape_fourth_derivatives.size()!=0)
+        {
+          Assert(shape_fourth_derivatives.size()==n_shape_functions*n_points,
+                 ExcInternalError());
+          grad4.resize(n_shape_functions);
+        }
+
+
+      if (shape_values.size()!=0 ||
+          shape_derivatives.size()!=0 ||
+          shape_second_derivatives.size()!=0 ||
+          shape_third_derivatives.size()!=0 ||
+          shape_fourth_derivatives.size()!=0 )
+        for (unsigned int point=0; point<n_points; ++point)
+          {
+            tensor_pols.compute(unit_points[point], values, grads, grad2, grad3, grad4);
+
+            if (shape_values.size()!=0)
+              for (unsigned int i=0; i<n_shape_functions; ++i)
+                shape(point,renumber[i]) = values[i];
+
+            if (shape_derivatives.size()!=0)
+              for (unsigned int i=0; i<n_shape_functions; ++i)
+                derivative(point,renumber[i]) = grads[i];
+
+            if (shape_second_derivatives.size()!=0)
+              for (unsigned int i=0; i<n_shape_functions; ++i)
+                second_derivative(point,renumber[i]) = grad2[i];
+
+            if (shape_third_derivatives.size()!=0)
+              for (unsigned int i=0; i<n_shape_functions; ++i)
+                third_derivative(point,renumber[i]) = grad3[i];
+
+            if (shape_fourth_derivatives.size()!=0)
+              for (unsigned int i=0; i<n_shape_functions; ++i)
+                fourth_derivative(point,renumber[i]) = grad4[i];
+          }
+    }
 }
 
 
@@ -732,7 +843,7 @@ MappingQ1<dim,spacedim>::get_data (const UpdateFlags update_flags,
 {
   InternalData *data = new InternalData(1);
   data->initialize (requires_update_flags(update_flags), q, q.size());
-  compute_shapes (q.get_points(), *data);
+  data->compute_shape_function_values (q.get_points());
 
   return data;
 }
@@ -748,8 +859,7 @@ 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());
-  compute_shapes (QProjector<dim>::project_to_all_faces(quadrature).get_points(),
-                  *data);
+  data->compute_shape_function_values (QProjector<dim>::project_to_all_faces(quadrature).get_points());
 
   return data;
 }
@@ -765,8 +875,7 @@ 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());
-  compute_shapes (QProjector<dim>::project_to_all_subfaces(quadrature).get_points(),
-                  *data);
+  data->compute_shape_function_values (QProjector<dim>::project_to_all_subfaces(quadrature).get_points());
 
 
   return data;
@@ -2513,10 +2622,7 @@ transform_real_to_unit_cell_internal
 
   Point<dim> p_unit = initial_p_unit;
 
-  if (dynamic_cast<typename MappingQ<dim,spacedim>::InternalData *>(&mdata) == 0)
-    this->compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata);
-  else
-    dynamic_cast<const MappingQ<dim,spacedim>*>(this)->compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata);
+  mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit));
 
   Point<spacedim> p_real = transform_unit_to_real_cell_internal(mdata);
   Tensor<1,spacedim> f = p_real-p;
@@ -2604,11 +2710,7 @@ transform_real_to_unit_cell_internal
 
           // shape values and derivatives
           // at new p_unit point
-          if (dynamic_cast<typename MappingQ<dim,spacedim>::InternalData *>(&mdata) == 0)
-            this->compute_shapes(std::vector<Point<dim> > (1, p_unit_trial), mdata);
-          else
-            dynamic_cast<const MappingQ<dim,spacedim>*>(this)->compute_shapes(std::vector<Point<dim> > (1, p_unit_trial), mdata);
-
+          mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit_trial));
 
           // f(x)
           Point<spacedim> p_real_trial = transform_unit_to_real_cell_internal(mdata);
@@ -2750,10 +2852,8 @@ transform_real_to_unit_cell_internal_codim1
   Tensor<2,dim1>  df;
 
   //Evaluate first and second derivatives
-  if (dynamic_cast<typename MappingQ<dim,spacedim>::InternalData *>(&mdata) == 0)
-    this->compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata);
-  else
-    dynamic_cast<const MappingQ<dim,spacedim>*>(this)->compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata);
+  mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit));
+
   for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
     {
       const Tensor<1,dim1>   &grad_phi_k = mdata.derivative(0,k);
@@ -2805,10 +2905,8 @@ transform_real_to_unit_cell_internal_codim1
             D2F[j][l].clear();
         }
 
-      if (dynamic_cast<typename MappingQ<dim,spacedim>::InternalData *>(&mdata) == 0)
-        this->compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata);
-      else
-        dynamic_cast<const MappingQ<dim,spacedim>*>(this)->compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata);
+      mdata.compute_shape_function_values(std::vector<Point<dim> > (1, p_unit));
+
       for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
         {
           const Tensor<1,dim1>   &grad_phi_k = mdata.derivative(0,k);