Get rid of the awkward dispatch mechanism in MappingQ*::compute_shapes().

In particular, remove the manual dispatch to either the Q1 or Qp functions.

diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h
index 317eb7953a886253929062723358bd21c7c2b8f7..f6b588a1ecdf3d3d86bc83b0ccf77fdaa6b0e1db 100644 (file)
--- a/include/deal.II/fe/mapping_q.h
+++ b/include/deal.II/fe/mapping_q.h
@@ -261,8 +261,7 @@ public:
     /**
      * Constructor.
      */
-    InternalData (const unsigned int polynomial_degree,
-                  const unsigned int n_shape_functions);
+    InternalData (const unsigned int polynomial_degree);
 
 
     /**
@@ -334,13 +333,11 @@ protected:
   /**
    * @}
    */
-
-  /**
-   * Compute shape values and/or derivatives.
-   */
-  virtual void
-  compute_shapes_virtual (const std::vector<Point<dim> > &unit_points,
-                          typename MappingQ1<dim,spacedim>::InternalData &data) const;
+public:
+  void
+  compute_shapes (const std::vector<Point<dim> > &unit_points,
+                  typename MappingQ1<dim,spacedim>::InternalData &data) const;
+protected:
 
   /**
    * This function is needed by the constructor of

diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h
index 082f5e74669f084c65a96cbd2cac40e25c2449a5..6a32f40cebdae2541feeed928e0915fb0878fb2d 100644 (file)
--- a/include/deal.II/fe/mapping_q1.h
+++ b/include/deal.II/fe/mapping_q1.h
@@ -161,10 +161,10 @@ public:
   {
   public:
     /**
-     * Constructor. Pass the number of shape functions.
+     * Constructor. The argument denotes the polynomial degree of
+     * the mapping to which this object will correspond.
      */
-    InternalData(const unsigned int polynomial_degree,
-                 const unsigned int n_shape_functions);
+    InternalData(const unsigned int polynomial_degree);
 
     /**
      * Initialize the object's member variables related to cell data
@@ -315,8 +315,11 @@ public:
      * the number of vertices per cell. However, since also derived classes
      * use this class (e.g. the Mapping_Q() class), the number of shape
      * functions may also be different.
+     *
+     * In general, it is $(p+1)^\text{dim}$, where $p$ is the
+     * polynomial degree of the mapping.
      */
-    unsigned int n_shape_functions;
+    const unsigned int n_shape_functions;
 
     /**
      * Tensors of covariant transformation at each of the quadrature points.
@@ -472,12 +475,6 @@ protected:
   void compute_shapes (const std::vector<Point<dim> > &unit_points,
                        InternalData &data) const;
 
-  /**
-   * Compute shape values and/or derivatives.
-   */
-  virtual void compute_shapes_virtual (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 08068cc01ca36689750b323fc171dc291e0cd5af..e2d6ea4737f77d6d5c9c8a71d3955cc688f3b2f9 100644 (file)
--- a/source/fe/mapping_q.cc
+++ b/source/fe/mapping_q.cc
@@ -35,14 +35,11 @@ DEAL_II_NAMESPACE_OPEN
 
 
 template<int dim, int spacedim>
-MappingQ<dim,spacedim>::InternalData::InternalData (const unsigned int polynomial_degree,
-                                                    const unsigned int n_shape_functions)
+MappingQ<dim,spacedim>::InternalData::InternalData (const unsigned int polynomial_degree)
   :
-  MappingQ1<dim,spacedim>::InternalData(polynomial_degree,
-                                        n_shape_functions),
+  MappingQ1<dim,spacedim>::InternalData(polynomial_degree),
   use_mapping_q1_on_current_cell(false),
-  mapping_q1_data(1,
-                  GeometryInfo<dim>::vertices_per_cell)
+  mapping_q1_data(1)
 {}
 
 
@@ -141,11 +138,10 @@ MappingQ<dim,spacedim>::~MappingQ ()
 }
 
 
-
 template<int dim, int spacedim>
 void
-MappingQ<dim,spacedim>::compute_shapes_virtual (const std::vector<Point<dim> > &unit_points,
-                                                typename MappingQ1<dim,spacedim>::InternalData &data) const
+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;
@@ -227,7 +223,7 @@ typename MappingQ<dim,spacedim>::InternalData *
 MappingQ<dim,spacedim>::get_data (const UpdateFlags update_flags,
                                   const Quadrature<dim> &quadrature) const
 {
-  InternalData *data = new InternalData(degree, n_shape_functions);
+  InternalData *data = new InternalData(degree);
 
   // fill the data of both the Q_p and the Q_1 objects in parallel
   Threads::TaskGroup<> tasks;
@@ -249,7 +245,7 @@ MappingQ<dim,spacedim>::get_data (const UpdateFlags update_flags,
   // TODO: parallelize this as well
   this->compute_shapes (quadrature.get_points(), *data);
   if (!use_mapping_q_on_all_cells)
-    this->compute_shapes (quadrature.get_points(), data->mapping_q1_data);
+    this->MappingQ1<dim,spacedim>::compute_shapes (quadrature.get_points(), data->mapping_q1_data);
 
 
   return data;
@@ -262,7 +258,7 @@ typename Mapping<dim,spacedim>::InternalDataBase *
 MappingQ<dim,spacedim>::get_face_data (const UpdateFlags update_flags,
                                        const Quadrature<dim-1>& quadrature) const
 {
-  InternalData *data = new InternalData(degree, n_shape_functions);
+  InternalData *data = new InternalData(degree);
   const Quadrature<dim> q (QProjector<dim>::project_to_all_faces(quadrature));
 
   // fill the data of both the Q_p and the Q_1 objects in parallel
@@ -285,7 +281,7 @@ MappingQ<dim,spacedim>::get_face_data (const UpdateFlags update_flags,
   // TODO: parallelize this as well
   this->compute_shapes (q.get_points(), *data);
   if (!use_mapping_q_on_all_cells)
-    this->compute_shapes (q.get_points(), data->mapping_q1_data);
+    this->MappingQ1<dim,spacedim>::compute_shapes (q.get_points(), data->mapping_q1_data);
 
   return data;
 }
@@ -297,7 +293,7 @@ typename Mapping<dim,spacedim>::InternalDataBase *
 MappingQ<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
                                           const Quadrature<dim-1>& quadrature) const
 {
-  InternalData *data = new InternalData(degree, n_shape_functions);
+  InternalData *data = new InternalData(degree);
   const Quadrature<dim> q (QProjector<dim>::project_to_all_subfaces(quadrature));
 
   // fill the data of both the Q_p and the Q_1 objects in parallel
@@ -320,7 +316,7 @@ MappingQ<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
   // TODO: parallelize this as well
   this->compute_shapes (q.get_points(), *data);
   if (!use_mapping_q_on_all_cells)
-    this->compute_shapes (q.get_points(), data->mapping_q1_data);
+    this->MappingQ1<dim,spacedim>::compute_shapes (q.get_points(), data->mapping_q1_data);
 
 
   return data;
@@ -618,7 +614,7 @@ MappingQ<dim,spacedim>::compute_laplace_vector(Table<2,double> &lvs) const
   const QGauss<dim> quadrature(degree+1);
   const unsigned int n_q_points=quadrature.size();
 
-  InternalData quadrature_data(degree, n_shape_functions);
+  InternalData quadrature_data(degree);
   quadrature_data.shape_derivatives.resize(n_shape_functions * n_q_points);
   this->compute_shapes(quadrature.get_points(), quadrature_data);
 

diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc
index e6a176f48ab0033d3d7108a3aacf806aff17daf6..69ce4b5d8f8c65df270fd2e7e0951298d36bbcb7 100644 (file)
--- a/source/fe/mapping_q1.cc
+++ b/source/fe/mapping_q1.cc
@@ -42,16 +42,11 @@ const unsigned int MappingQ1<dim,spacedim>::n_shape_functions;
 
 
 template<int dim, int spacedim>
-MappingQ1<dim,spacedim>::InternalData::InternalData (const unsigned int polynomial_degree,
-                                                     const unsigned int n_shape_functions)
+MappingQ1<dim,spacedim>::InternalData::InternalData (const unsigned int polynomial_degree)
   :
   polynomial_degree (polynomial_degree),
-  n_shape_functions (n_shape_functions)
-{
-  Assert (n_shape_functions ==
-          Utilities::fixed_power<dim>(polynomial_degree+1),
-          ExcInternalError());
-}
+  n_shape_functions (Utilities::fixed_power<dim>(polynomial_degree+1))
+{}
 
 
 
@@ -205,20 +200,6 @@ MappingQ1<dim,spacedim>::MappingQ1 ()
 
 
 
-template<int dim, int spacedim>
-void
-MappingQ1<dim,spacedim>::compute_shapes (const std::vector<Point<dim> > &unit_points,
-                                         InternalData &data) const
-{
-  // choose either the function implemented in this class, or whatever
-  // a virtual function call resolves to
-  if (dynamic_cast<typename MappingQ<dim,spacedim>::InternalData *>(&data) == 0)
-    MappingQ1<dim,spacedim>::compute_shapes_virtual(unit_points, data);
-  else
-    compute_shapes_virtual(unit_points, data);
-}
-
-
 namespace internal
 {
   namespace MappingQ1
@@ -345,9 +326,9 @@ namespace internal
 
     template <int spacedim>
     void
-    compute_shapes_virtual (const unsigned int            n_shape_functions,
-                            const std::vector<Point<1> > &unit_points,
-                            typename dealii::MappingQ1<1,spacedim>::InternalData &data)
+    compute_shapes (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();
@@ -410,9 +391,9 @@ namespace internal
 
     template <int spacedim>
     void
-    compute_shapes_virtual (const unsigned int            n_shape_functions,
-                            const std::vector<Point<2> > &unit_points,
-                            typename dealii::MappingQ1<2,spacedim>::InternalData &data)
+    compute_shapes (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();
@@ -488,9 +469,9 @@ namespace internal
 
     template <int spacedim>
     void
-    compute_shapes_virtual (const unsigned int            n_shape_functions,
-                            const std::vector<Point<3> > &unit_points,
-                            typename dealii::MappingQ1<3,spacedim>::InternalData &data)
+    compute_shapes (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();
@@ -674,13 +655,12 @@ namespace internal
 
 template<int dim, int spacedim>
 void
-MappingQ1<dim, spacedim>::
-compute_shapes_virtual (const std::vector<Point<dim> > &unit_points,
-                        InternalData &data) const
+MappingQ1<dim,spacedim>::compute_shapes (const std::vector<Point<dim> > &unit_points,
+                                         InternalData &data) const
 {
   internal::MappingQ1::
-  compute_shapes_virtual<spacedim> (n_shape_functions,
-                                    unit_points, data);
+  compute_shapes<spacedim> (n_shape_functions,
+                            unit_points, data);
 }
 
 
@@ -750,7 +730,7 @@ typename MappingQ1<dim,spacedim>::InternalData *
 MappingQ1<dim,spacedim>::get_data (const UpdateFlags update_flags,
                                    const Quadrature<dim> &q) const
 {
-  InternalData *data = new InternalData(1, n_shape_functions);
+  InternalData *data = new InternalData(1);
   data->initialize (requires_update_flags(update_flags), q, q.size());
   compute_shapes (q.get_points(), *data);
 
@@ -764,7 +744,7 @@ typename Mapping<dim,spacedim>::InternalDataBase *
 MappingQ1<dim,spacedim>::get_face_data (const UpdateFlags        update_flags,
                                         const Quadrature<dim-1> &quadrature) const
 {
-  InternalData *data = new InternalData(1, n_shape_functions);
+  InternalData *data = new InternalData(1);
   data->initialize_face (requires_update_flags(update_flags),
                          QProjector<dim>::project_to_all_faces(quadrature),
                          quadrature.size());
@@ -781,7 +761,7 @@ typename Mapping<dim,spacedim>::InternalDataBase *
 MappingQ1<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
                                            const Quadrature<dim-1>& quadrature) const
 {
-  InternalData *data = new InternalData(1, n_shape_functions);
+  InternalData *data = new InternalData(1);
   data->initialize_face (requires_update_flags(update_flags),
                          QProjector<dim>::project_to_all_subfaces(quadrature),
                          quadrature.size());
@@ -2533,7 +2513,11 @@ transform_real_to_unit_cell_internal
 
   Point<dim> p_unit = initial_p_unit;
 
-  compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata);
+  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);
+
   Point<spacedim> p_real = transform_unit_to_real_cell_internal(mdata);
   Tensor<1,spacedim> f = p_real-p;
 
@@ -2620,7 +2604,11 @@ transform_real_to_unit_cell_internal
 
           // shape values and derivatives
           // at new p_unit point
-          compute_shapes(std::vector<Point<dim> > (1, p_unit_trial), mdata);
+          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);
+
 
           // f(x)
           Point<spacedim> p_real_trial = transform_unit_to_real_cell_internal(mdata);
@@ -2762,7 +2750,10 @@ transform_real_to_unit_cell_internal_codim1
   Tensor<2,dim1>  df;
 
   //Evaluate first and second derivatives
-  compute_shapes(std::vector<Point<dim1> > (1, p_unit), mdata);
+  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);
   for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
     {
       const Tensor<1,dim1>   &grad_phi_k = mdata.derivative(0,k);
@@ -2814,7 +2805,10 @@ transform_real_to_unit_cell_internal_codim1
             D2F[j][l].clear();
         }
 
-      compute_shapes(std::vector<Point<dim1> > (1, p_unit), mdata);
+      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);
       for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
         {
           const Tensor<1,dim1>   &grad_phi_k = mdata.derivative(0,k);