From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Wed, 2 Sep 2015 19:46:45 +0000 (-0500)
Subject: Attempt to disentangle MappingQ1 and MappingQ.
X-Git-Tag: v8.4.0-rc2~474^2~11
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a332d1db2803d48c3720c45b63b0bd067138b5f5;p=dealii.git

Attempt to disentangle MappingQ1 and MappingQ.

Store the polynomial degree in MappingQ1::InternalData, rather than
which kind of mapping created it.
---

diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h
index ec29075344..317eb7953a 100644
--- a/include/deal.II/fe/mapping_q.h
+++ b/include/deal.II/fe/mapping_q.h
@@ -238,7 +238,7 @@ private:
    * @{
    */
 
-protected:
+public:
 
   /**
    * Storage for internal data of this mapping. See Mapping::InternalDataBase
@@ -261,7 +261,8 @@ protected:
     /**
      * Constructor.
      */
-    InternalData (const unsigned int n_shape_functions);
+    InternalData (const unsigned int polynomial_degree,
+                  const unsigned int n_shape_functions);
 
 
     /**
@@ -284,6 +285,8 @@ protected:
     typename MappingQ1<dim,spacedim>::InternalData mapping_q1_data;
   };
 
+protected:
+
   // documentation can be found in Mapping::get_data()
   virtual
   InternalData *
diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h
index 1b4c05bc77..082f5e7466 100644
--- a/include/deal.II/fe/mapping_q1.h
+++ b/include/deal.II/fe/mapping_q1.h
@@ -163,7 +163,8 @@ public:
     /**
      * Constructor. Pass the number of shape functions.
      */
-    InternalData(const unsigned int n_shape_functions);
+    InternalData(const unsigned int polynomial_degree,
+                 const unsigned int n_shape_functions);
 
     /**
      * Initialize the object's member variables related to cell data
@@ -303,12 +304,11 @@ public:
     std::vector<std::vector<Tensor<1,dim> > > unit_tangentials;
 
     /**
-     * Default value of this flag is @p true. If <tt>*this</tt> is an object
-     * of a derived class, this flag is set to @p false. (This is, for example,
-     * the case for MappingQ, which derives MappingQ::InternalData from the
-     * current MappingQ1::InternalData.)
+     * The polynomial degree of the mapping. Since the objects here
+     * are also used (with minor adjustments) by MappingQ, we need to
+     * store this.
      */
-    bool is_mapping_q1_data;
+    unsigned int polynomial_degree;
 
     /**
      * Number of shape functions. If this is a Q1 mapping, then it is simply
@@ -466,8 +466,8 @@ protected:
    * 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.is_mapping_q1_data</tt>
-   * equals @p true or @p false.
+   * 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;
diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc
index 4f811729ed..08068cc01c 100644
--- a/source/fe/mapping_q.cc
+++ b/source/fe/mapping_q.cc
@@ -35,14 +35,15 @@ DEAL_II_NAMESPACE_OPEN
 
 
 template<int dim, int spacedim>
-MappingQ<dim,spacedim>::InternalData::InternalData (const unsigned int n_shape_functions)
+MappingQ<dim,spacedim>::InternalData::InternalData (const unsigned int polynomial_degree,
+                                                    const unsigned int n_shape_functions)
   :
-  MappingQ1<dim,spacedim>::InternalData(n_shape_functions),
+  MappingQ1<dim,spacedim>::InternalData(polynomial_degree,
+                                        n_shape_functions),
   use_mapping_q1_on_current_cell(false),
-  mapping_q1_data(1 << dim)
-{
-  this->is_mapping_q1_data=false;
-}
+  mapping_q1_data(1,
+                  GeometryInfo<dim>::vertices_per_cell)
+{}
 
 
 
@@ -226,7 +227,7 @@ typename MappingQ<dim,spacedim>::InternalData *
 MappingQ<dim,spacedim>::get_data (const UpdateFlags update_flags,
                                   const Quadrature<dim> &quadrature) const
 {
-  InternalData *data = new InternalData(n_shape_functions);
+  InternalData *data = new InternalData(degree, n_shape_functions);
 
   // fill the data of both the Q_p and the Q_1 objects in parallel
   Threads::TaskGroup<> tasks;
@@ -261,7 +262,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(n_shape_functions);
+  InternalData *data = new InternalData(degree, n_shape_functions);
   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
@@ -296,7 +297,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(n_shape_functions);
+  InternalData *data = new InternalData(degree, n_shape_functions);
   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
@@ -617,7 +618,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(n_shape_functions);
+  InternalData quadrature_data(degree, n_shape_functions);
   quadrature_data.shape_derivatives.resize(n_shape_functions * n_q_points);
   this->compute_shapes(quadrature.get_points(), quadrature_data);
 
@@ -1009,15 +1010,13 @@ transform (const VectorSlice<const std::vector<Tensor<1,dim> > >   input,
   Assert(q1_data!=0, ExcInternalError());
 
   // If it is a genuine MappingQ::InternalData, we have to test further
-  if (!q1_data->is_mapping_q1_data)
+  if (const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data))
     {
-      Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
-              ExcInternalError());
-      const InternalData &data = static_cast<const InternalData &>(mapping_data);
       // If we only use the Q1-portion, we have to extract that data object
-      if (data.use_mapping_q1_on_current_cell)
-        q1_data = &data.mapping_q1_data;
+      if (data->use_mapping_q1_on_current_cell)
+        q1_data = &data->mapping_q1_data;
     }
+
   // Now, q1_data should have the right tensors in it and we call the base
   // classes transform function
   MappingQ1<dim,spacedim>::transform(input, mapping_type, *q1_data, output);
@@ -1041,15 +1040,13 @@ transform (const VectorSlice<const std::vector<DerivativeForm<1, dim ,spacedim>
   Assert(q1_data!=0, ExcInternalError());
 
   // If it is a genuine MappingQ::InternalData, we have to test further
-  if (!q1_data->is_mapping_q1_data)
+  if (const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data))
     {
-      Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
-              ExcInternalError());
-      const InternalData &data = static_cast<const InternalData &>(mapping_data);
       // If we only use the Q1-portion, we have to extract that data object
-      if (data.use_mapping_q1_on_current_cell)
-        q1_data = &data.mapping_q1_data;
+      if (data->use_mapping_q1_on_current_cell)
+        q1_data = &data->mapping_q1_data;
     }
+
   // Now, q1_data should have the right tensors in it and we call the base
   // classes transform function
   MappingQ1<dim,spacedim>::transform(input, mapping_type, *q1_data, output);
@@ -1071,15 +1068,13 @@ transform (const VectorSlice<const std::vector<Tensor<2, dim> > >  input,
   Assert(q1_data!=0, ExcInternalError());
 
   // If it is a genuine MappingQ::InternalData, we have to test further
-  if (!q1_data->is_mapping_q1_data)
+  if (const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data))
     {
-      Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
-              ExcInternalError());
-      const InternalData &data = static_cast<const InternalData &>(mapping_data);
       // If we only use the Q1-portion, we have to extract that data object
-      if (data.use_mapping_q1_on_current_cell)
-        q1_data = &data.mapping_q1_data;
+      if (data->use_mapping_q1_on_current_cell)
+        q1_data = &data->mapping_q1_data;
     }
+
   // Now, q1_data should have the right tensors in it and we call the base
   // classes transform function
   MappingQ1<dim,spacedim>::transform(input, mapping_type, *q1_data, output);
@@ -1103,15 +1098,13 @@ transform (const VectorSlice<const std::vector<DerivativeForm<2, dim ,spacedim>
   Assert(q1_data!=0, ExcInternalError());
 
   // If it is a genuine MappingQ::InternalData, we have to test further
-  if (!q1_data->is_mapping_q1_data)
+  if (const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data))
     {
-      Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
-              ExcInternalError());
-      const InternalData &data = static_cast<const InternalData &>(mapping_data);
       // If we only use the Q1-portion, we have to extract that data object
-      if (data.use_mapping_q1_on_current_cell)
-        q1_data = &data.mapping_q1_data;
+      if (data->use_mapping_q1_on_current_cell)
+        q1_data = &data->mapping_q1_data;
     }
+
   // Now, q1_data should have the right tensors in it and we call the base
   // classes transform function
   MappingQ1<dim,spacedim>::transform(input, mapping_type, *q1_data, output);
@@ -1133,15 +1126,13 @@ transform (const VectorSlice<const std::vector<Tensor<3, dim> > >  input,
   Assert(q1_data!=0, ExcInternalError());
 
   // If it is a genuine MappingQ::InternalData, we have to test further
-  if (!q1_data->is_mapping_q1_data)
+  if (const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data))
     {
-      Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
-              ExcInternalError());
-      const InternalData &data = static_cast<const InternalData &>(mapping_data);
       // If we only use the Q1-portion, we have to extract that data object
-      if (data.use_mapping_q1_on_current_cell)
-        q1_data = &data.mapping_q1_data;
+      if (data->use_mapping_q1_on_current_cell)
+        q1_data = &data->mapping_q1_data;
     }
+
   // Now, q1_data should have the right tensors in it and we call the base
   // classes transform function
   MappingQ1<dim,spacedim>::transform(input, mapping_type, *q1_data, output);
diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc
index 56e83e29fa..e6a176f48a 100644
--- a/source/fe/mapping_q1.cc
+++ b/source/fe/mapping_q1.cc
@@ -26,6 +26,7 @@
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/mapping_q1.h>
+#include <deal.II/fe/mapping_q.h>
 #include <deal.II/fe/mapping_q1_eulerian.h>
 
 #include <cmath>
@@ -41,11 +42,16 @@ const unsigned int MappingQ1<dim,spacedim>::n_shape_functions;
 
 
 template<int dim, int spacedim>
-MappingQ1<dim,spacedim>::InternalData::InternalData (const unsigned int n_shape_functions)
+MappingQ1<dim,spacedim>::InternalData::InternalData (const unsigned int polynomial_degree,
+                                                     const unsigned int n_shape_functions)
   :
-  is_mapping_q1_data(true),
+  polynomial_degree (polynomial_degree),
   n_shape_functions (n_shape_functions)
-{}
+{
+  Assert (n_shape_functions ==
+          Utilities::fixed_power<dim>(polynomial_degree+1),
+          ExcInternalError());
+}
 
 
 
@@ -63,7 +69,7 @@ MappingQ1<dim,spacedim>::InternalData::memory_consumption () const
           MemoryConsumption::memory_consumption (mapping_support_points) +
           MemoryConsumption::memory_consumption (cell_of_current_support_points) +
           MemoryConsumption::memory_consumption (volume_elements) +
-          MemoryConsumption::memory_consumption (is_mapping_q1_data) +
+          MemoryConsumption::memory_consumption (polynomial_degree) +
           MemoryConsumption::memory_consumption (n_shape_functions));
 }
 
@@ -204,10 +210,9 @@ 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 (data.is_mapping_q1_data)
+  // 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);
@@ -745,7 +750,7 @@ typename MappingQ1<dim,spacedim>::InternalData *
 MappingQ1<dim,spacedim>::get_data (const UpdateFlags update_flags,
                                    const Quadrature<dim> &q) const
 {
-  InternalData *data = new InternalData(n_shape_functions);
+  InternalData *data = new InternalData(1, n_shape_functions);
   data->initialize (requires_update_flags(update_flags), q, q.size());
   compute_shapes (q.get_points(), *data);
 
@@ -759,7 +764,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(n_shape_functions);
+  InternalData *data = new InternalData(1, n_shape_functions);
   data->initialize_face (requires_update_flags(update_flags),
                          QProjector<dim>::project_to_all_faces(quadrature),
                          quadrature.size());
@@ -776,7 +781,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(n_shape_functions);
+  InternalData *data = new InternalData(1, n_shape_functions);
   data->initialize_face (requires_update_flags(update_flags),
                          QProjector<dim>::project_to_all_subfaces(quadrature),
                          quadrature.size());