From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Sun, 6 Sep 2015 21:47:01 +0000 (-0500)
Subject: Let MappingQ have a MappingQ1 member.
X-Git-Tag: v8.4.0-rc2~439^2~8
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=306a739e62608bdf3d0aaad0b6e2a4ecd936905e;p=dealii.git

Let MappingQ have a MappingQ1 member.

This mostly breaks the derivation of MappingQ from MappingQ1. There are
a few hold-outs left, though.
---

diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h
index b954602e3d..b2c49db3e4 100644
--- a/include/deal.II/fe/mapping_q.h
+++ b/include/deal.II/fe/mapping_q.h
@@ -62,16 +62,25 @@ public:
    *
    * The second argument determines whether the higher order mapping should
    * also be used on interior cells. If its value is <code>false</code> (the
-   * default), the a lower-order mapping is used in the interior. This is
+   * default), then a lower order mapping is used in the interior. This is
    * sufficient for most cases where higher order mappings are only used to
    * better approximate the boundary. In that case, cells bounded by straight
    * lines are acceptable in the interior. However, there are cases where one
    * would also like to use a higher order mapping in the interior. The
    * MappingQEulerian class is one such case.
+   *
+   * The value of @p use_mapping_q_on_all_cells is ignore if @p dim is not
+   * equal to @p spacedim, i.e., if we are considering meshes on surfaces
+   * embedded into higher dimensional spaces.
    */
   MappingQ (const unsigned int polynomial_degree,
             const bool use_mapping_q_on_all_cells = false);
 
+  /**
+   * Copy constructor.
+   */
+  MappingQ (const MappingQ<dim,spacedim> &mapping);
+
   /**
    * 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.
@@ -277,13 +286,13 @@ protected:
 
   // documentation can be found in Mapping::get_face_data()
   virtual
-  typename Mapping<dim,spacedim>::InternalDataBase *
+  InternalData *
   get_face_data (const UpdateFlags flags,
                  const Quadrature<dim-1>& quadrature) const;
 
   // documentation can be found in Mapping::get_subface_data()
   virtual
-  typename Mapping<dim,spacedim>::InternalDataBase *
+  InternalData *
   get_subface_data (const UpdateFlags flags,
                     const Quadrature<dim-1>& quadrature) const;
 
@@ -440,6 +449,12 @@ protected:
    */
   const FE_Q<dim> feq;
 
+  /**
+   * Pointer to a Q1 mapping to be used on interior cells unless
+   * use_mapping_q_on_all_cells was set in the call to the
+   * constructor.
+   */
+  std_cxx11::unique_ptr<const MappingQ1<dim,spacedim> > q1_mapping;
 
   /*
    * The default line support points. These are used when computing
diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h
index 3791c7da33..cb5df7fd60 100644
--- a/include/deal.II/fe/mapping_q_generic.h
+++ b/include/deal.II/fe/mapping_q_generic.h
@@ -29,6 +29,9 @@
 
 DEAL_II_NAMESPACE_OPEN
 
+template <int,int> class MappingQ;
+
+
 /*!@addtogroup mapping */
 /*@{*/
 
@@ -374,7 +377,6 @@ public:
     mutable std::vector<double> volume_elements;
   };
 
-protected:
 
   // documentation can be found in Mapping::requires_update_flags()
   virtual
@@ -389,13 +391,13 @@ protected:
 
   // documentation can be found in Mapping::get_face_data()
   virtual
-  typename Mapping<dim,spacedim>::InternalDataBase *
+  InternalData *
   get_face_data (const UpdateFlags flags,
                  const Quadrature<dim-1>& quadrature) const;
 
   // documentation can be found in Mapping::get_subface_data()
   virtual
-  typename Mapping<dim,spacedim>::InternalDataBase *
+  InternalData *
   get_subface_data (const UpdateFlags flags,
                     const Quadrature<dim-1>& quadrature) const;
 
@@ -429,6 +431,7 @@ protected:
    * @}
    */
 
+protected:
 
   /**
    * The degree of the polynomials used as shape functions for the mapping
@@ -450,6 +453,12 @@ protected:
   void
   compute_mapping_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                                   std::vector<Point<spacedim> > &a) const = 0;
+
+  /**
+   * Make MappingQ a friend since it needs to call the
+   * fill_fe_values() functions on its MappingQ1 sub-object.
+   */
+  template <int, int> friend class MappingQ;
 };
 
 
diff --git a/source/fe/mapping_c1.cc b/source/fe/mapping_c1.cc
index 52e53a9b33..9e69706fdd 100644
--- a/source/fe/mapping_c1.cc
+++ b/source/fe/mapping_c1.cc
@@ -193,7 +193,7 @@ template<int dim, int spacedim>
 Mapping<dim, spacedim> *
 MappingC1<dim,spacedim>::clone () const
 {
-  return new MappingC1<dim,spacedim>(*this);
+  return new MappingC1<dim,spacedim>();
 }
 
 
diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc
index 5eb71043b0..ef90792511 100644
--- a/source/fe/mapping_q.cc
+++ b/source/fe/mapping_q.cc
@@ -36,7 +36,7 @@ DEAL_II_NAMESPACE_OPEN
 template<int dim, int spacedim>
 MappingQ<dim,spacedim>::InternalData::InternalData (const unsigned int polynomial_degree)
   :
-  MappingQ1<dim,spacedim>::InternalData(polynomial_degree),
+  MappingQGeneric<dim,spacedim>::InternalData(polynomial_degree),
   use_mapping_q1_on_current_cell(false)
 {}
 
@@ -66,6 +66,11 @@ MappingQ<dim,spacedim>::MappingQ (const unsigned int degree,
   use_mapping_q_on_all_cells (use_mapping_q_on_all_cells
                               || (dim != spacedim)),
   feq(degree),
+  q1_mapping (this->use_mapping_q_on_all_cells
+              ?
+              0
+              :
+              new MappingQ1<dim,spacedim>()),
   line_support_points(degree+1)
 {
   Assert(n_inner+n_outer==Utilities::fixed_power<dim>(degree+1),
@@ -83,6 +88,36 @@ MappingQ<dim,spacedim>::MappingQ (const unsigned int degree,
 
 
 
+template<int dim, int spacedim>
+MappingQ<dim,spacedim>::MappingQ (const MappingQ<dim,spacedim> &mapping)
+  :
+  MappingQ1<dim,spacedim>(mapping.get_degree()),
+  n_inner(mapping.n_inner),
+  n_outer(mapping.n_outer),
+  use_mapping_q_on_all_cells (mapping.use_mapping_q_on_all_cells),
+  feq(mapping.get_degree()),
+  q1_mapping (mapping.q1_mapping != 0
+              ?
+              dynamic_cast<MappingQ1<dim,spacedim>*>(mapping.q1_mapping->clone())
+              :
+              0),
+  line_support_points(mapping.line_support_points)
+{
+  Assert(n_inner+n_outer==Utilities::fixed_power<dim>(this->polynomial_degree+1),
+         ExcInternalError());
+
+  // build laplace_on_quad_vector
+  if (this->polynomial_degree>1)
+    {
+      if (dim >= 2)
+        set_laplace_on_quad_vector(laplace_on_quad_vector);
+      if (dim >= 3)
+        set_laplace_on_hex_vector(laplace_on_hex_vector);
+    }
+}
+
+
+
 template<int dim, int spacedim>
 typename MappingQ<dim,spacedim>::InternalData *
 MappingQ<dim,spacedim>::get_data (const UpdateFlags update_flags,
@@ -91,32 +126,24 @@ MappingQ<dim,spacedim>::get_data (const UpdateFlags update_flags,
   InternalData *data = new InternalData(this->polynomial_degree);
 
   // fill the data of both the Q_p and the Q_1 objects in parallel
-  Threads::TaskGroup<> tasks;
-  tasks += Threads::new_task (std_cxx11::function<void()>
-                              (std_cxx11::bind(&MappingQ1<dim,spacedim>::InternalData::initialize,
-                                               data,
-                                               update_flags,
-                                               std_cxx11::cref(quadrature),
-                                               quadrature.size())));
+  Threads::Task<>
+  task = Threads::new_task (std_cxx11::function<void()>
+                            (std_cxx11::bind(&MappingQ1<dim,spacedim>::InternalData::initialize,
+                                             data,
+                                             update_flags,
+                                             std_cxx11::cref(quadrature),
+                                             quadrature.size())));
   if (!use_mapping_q_on_all_cells)
-    {
-      data->mapping_q1_data.reset (new typename MappingQ1<dim,spacedim>::InternalData(1));
-      tasks += Threads::new_task (std_cxx11::function<void()>
-                                  (std_cxx11::bind(&MappingQ1<dim,spacedim>::InternalData::initialize,
-                                                   data->mapping_q1_data.get(),
-                                                   update_flags,
-                                                   std_cxx11::cref(quadrature),
-                                                   quadrature.size())));
-    }
+    data->mapping_q1_data.reset (q1_mapping->get_data (update_flags, quadrature));
 
-  tasks.join_all ();
+  task.join ();
   return data;
 }
 
 
 
 template<int dim, int spacedim>
-typename Mapping<dim,spacedim>::InternalDataBase *
+typename MappingQ<dim,spacedim>::InternalData *
 MappingQ<dim,spacedim>::get_face_data (const UpdateFlags update_flags,
                                        const Quadrature<dim-1>& quadrature) const
 {
@@ -124,32 +151,24 @@ MappingQ<dim,spacedim>::get_face_data (const UpdateFlags update_flags,
   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
-  Threads::TaskGroup<> tasks;
-  tasks += Threads::new_task (std_cxx11::function<void()>
-                              (std_cxx11::bind(&MappingQ1<dim,spacedim>::InternalData::initialize_face,
-                                               data,
-                                               update_flags,
-                                               std_cxx11::cref(q),
-                                               quadrature.size())));
+  Threads::Task<>
+  task = Threads::new_task (std_cxx11::function<void()>
+                            (std_cxx11::bind(&MappingQ1<dim,spacedim>::InternalData::initialize_face,
+                                             data,
+                                             update_flags,
+                                             std_cxx11::cref(q),
+                                             quadrature.size())));
   if (!use_mapping_q_on_all_cells)
-    {
-      data->mapping_q1_data.reset (new typename MappingQ1<dim,spacedim>::InternalData(1));
-      tasks += Threads::new_task (std_cxx11::function<void()>
-                                  (std_cxx11::bind(&MappingQ1<dim,spacedim>::InternalData::initialize_face,
-                                                   data->mapping_q1_data.get(),
-                                                   update_flags,
-                                                   std_cxx11::cref(q),
-                                                   quadrature.size())));
-    }
+    data->mapping_q1_data.reset (q1_mapping->get_face_data (update_flags, quadrature));
 
-  tasks.join_all ();
+  task.join ();
   return data;
 }
 
 
 
 template<int dim, int spacedim>
-typename Mapping<dim,spacedim>::InternalDataBase *
+typename MappingQ<dim,spacedim>::InternalData *
 MappingQ<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
                                           const Quadrature<dim-1>& quadrature) const
 {
@@ -157,25 +176,17 @@ MappingQ<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
   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
-  Threads::TaskGroup<> tasks;
-  tasks += Threads::new_task (std_cxx11::function<void()>
-                              (std_cxx11::bind(&MappingQ1<dim,spacedim>::InternalData::initialize_face,
-                                               data,
-                                               update_flags,
-                                               std_cxx11::cref(q),
-                                               quadrature.size())));
+  Threads::Task<>
+  task = Threads::new_task (std_cxx11::function<void()>
+                            (std_cxx11::bind(&MappingQ1<dim,spacedim>::InternalData::initialize_face,
+                                             data,
+                                             update_flags,
+                                             std_cxx11::cref(q),
+                                             quadrature.size())));
   if (!use_mapping_q_on_all_cells)
-    {
-      data->mapping_q1_data.reset (new typename MappingQ1<dim,spacedim>::InternalData(1));
-      tasks += Threads::new_task (std_cxx11::function<void()>
-                                  (std_cxx11::bind(&MappingQ1<dim,spacedim>::InternalData::initialize_face,
-                                                   data->mapping_q1_data.get(),
-                                                   update_flags,
-                                                   std_cxx11::cref(q),
-                                                   quadrature.size())));
-    }
+    data->mapping_q1_data.reset (q1_mapping->get_subface_data (update_flags, quadrature));
 
-  tasks.join_all ();
+  task.join ();
   return data;
 }
 
@@ -201,14 +212,6 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
   data.use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells
                                           || cell->has_boundary_lines());
 
-  // depending on this result, use this or the other data object for the
-  // mapping.
-  const typename MappingQ1<dim,spacedim>::InternalData *
-  p_data = (data.use_mapping_q1_on_current_cell
-            ?
-            data.mapping_q1_data.get()
-            :
-            &data);
 
   // call the base class. we need to ensure that the flag indicating whether
   // we can use some similarity has to be modified - for a general MappingQ,
@@ -224,11 +227,21 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
        CellSimilarity::invalid_next_cell
        :
        cell_similarity);
-  MappingQ1<dim,spacedim>::fill_fe_values(cell,
-                                          updated_cell_similarity,
-                                          quadrature,
-                                          *p_data,
-                                          output_data);
+
+  // depending on the results above, decide whether the Q1 mapping or
+  // the Qp mapping needs to handle this cell
+  if (data.use_mapping_q1_on_current_cell)
+    q1_mapping->fill_fe_values (cell,
+                                updated_cell_similarity,
+                                quadrature,
+                                *data.mapping_q1_data,
+                                output_data);
+  else
+    MappingQ1<dim,spacedim>::fill_fe_values(cell,
+                                            updated_cell_similarity,
+                                            quadrature,
+                                            data,
+                                            output_data);
 
   return updated_cell_similarity;
 }
@@ -259,20 +272,20 @@ fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator &
   data.use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells
                                           || cell->has_boundary_lines());
 
-  // depending on this result, use this or the other data object for the
-  // mapping and call the function in the base class with this
-  // data object to do the work
-  const typename MappingQ1<dim,spacedim>::InternalData *p_data
-    = (data.use_mapping_q1_on_current_cell
-       ?
-       data.mapping_q1_data.get()
-       :
-       &data);
-  MappingQ1<dim,spacedim>::fill_fe_face_values(cell,
-                                               face_no,
-                                               quadrature,
-                                               *p_data,
-                                               output_data);
+  // depending on the results above, decide whether the Q1 mapping or
+  // the Qp mapping needs to handle this cell
+  if (data.use_mapping_q1_on_current_cell)
+    q1_mapping->fill_fe_face_values (cell,
+                                     face_no,
+                                     quadrature,
+                                     *data.mapping_q1_data,
+                                     output_data);
+  else
+    MappingQ1<dim,spacedim>::fill_fe_face_values(cell,
+                                                 face_no,
+                                                 quadrature,
+                                                 data,
+                                                 output_data);
 }
 
 
@@ -301,21 +314,22 @@ fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterato
   data.use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells
                                           || cell->has_boundary_lines());
 
-  // depending on this result, use this or the other data object for the
-  // mapping and call the function in the base class with this
-  // data object to do the work
-  const typename MappingQ1<dim,spacedim>::InternalData *p_data
-    = (data.use_mapping_q1_on_current_cell
-       ?
-       data.mapping_q1_data.get()
-       :
-       &data);
-  MappingQ1<dim,spacedim>::fill_fe_subface_values(cell,
-                                                  face_no,
-                                                  subface_no,
-                                                  quadrature,
-                                                  *p_data,
-                                                  output_data);
+  // depending on the results above, decide whether the Q1 mapping or
+  // the Qp mapping needs to handle this cell
+  if (data.use_mapping_q1_on_current_cell)
+    q1_mapping->fill_fe_subface_values (cell,
+                                        face_no,
+                                        subface_no,
+                                        quadrature,
+                                        *data.mapping_q1_data,
+                                        output_data);
+  else
+    MappingQ1<dim,spacedim>::fill_fe_subface_values(cell,
+                                                    face_no,
+                                                    subface_no,
+                                                    quadrature,
+                                                    data,
+                                                    output_data);
 }
 
 
@@ -865,7 +879,7 @@ transform (const VectorSlice<const std::vector<Tensor<1,dim> > >   input,
   // whether we should in fact work on the Q1 portion of it
   if (const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data))
     if (data->use_mapping_q1_on_current_cell)
-      return MappingQ1<dim,spacedim>::transform (input, mapping_type, *data->mapping_q1_data, output);
+      return q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output);
 
   // otherwise just stick with what we already had
   MappingQ1<dim,spacedim>::transform(input, mapping_type, mapping_data, output);
@@ -890,7 +904,7 @@ transform (const VectorSlice<const std::vector<DerivativeForm<1, dim ,spacedim>
   // whether we should in fact work on the Q1 portion of it
   if (const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data))
     if (data->use_mapping_q1_on_current_cell)
-      return MappingQ1<dim,spacedim>::transform (input, mapping_type, *data->mapping_q1_data, output);
+      return q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output);
 
   // otherwise just stick with what we already had
   MappingQ1<dim,spacedim>::transform(input, mapping_type, mapping_data, output);
@@ -913,7 +927,7 @@ transform (const VectorSlice<const std::vector<Tensor<2, dim> > >  input,
   // whether we should in fact work on the Q1 portion of it
   if (const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data))
     if (data->use_mapping_q1_on_current_cell)
-      return MappingQ1<dim,spacedim>::transform (input, mapping_type, *data->mapping_q1_data, output);
+      return q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output);
 
   // otherwise just stick with what we already had
   MappingQ1<dim,spacedim>::transform(input, mapping_type, mapping_data, output);
@@ -938,7 +952,7 @@ transform (const VectorSlice<const std::vector<DerivativeForm<2, dim ,spacedim>
   // whether we should in fact work on the Q1 portion of it
   if (const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data))
     if (data->use_mapping_q1_on_current_cell)
-      return MappingQ1<dim,spacedim>::transform (input, mapping_type, *data->mapping_q1_data, output);
+      return q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output);
 
   // otherwise just stick with what we already had
   MappingQ1<dim,spacedim>::transform(input, mapping_type, mapping_data, output);
@@ -961,7 +975,7 @@ transform (const VectorSlice<const std::vector<Tensor<3, dim> > >  input,
   // whether we should in fact work on the Q1 portion of it
   if (const InternalData *data = dynamic_cast<const InternalData *>(&mapping_data))
     if (data->use_mapping_q1_on_current_cell)
-      return MappingQ1<dim,spacedim>::transform (input, mapping_type, *data->mapping_q1_data, output);
+      return q1_mapping->transform (input, mapping_type, *data->mapping_q1_data, output);
 
   // otherwise just stick with what we already had
   MappingQ1<dim,spacedim>::transform(input, mapping_type, mapping_data, output);
diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc
index 5fe49778e6..06a42dd06f 100644
--- a/source/fe/mapping_q_generic.cc
+++ b/source/fe/mapping_q_generic.cc
@@ -859,7 +859,7 @@ MappingQGeneric<dim,spacedim>::get_data (const UpdateFlags update_flags,
 
 
 template<int dim, int spacedim>
-typename Mapping<dim,spacedim>::InternalDataBase *
+typename MappingQGeneric<dim,spacedim>::InternalData *
 MappingQGeneric<dim,spacedim>::get_face_data (const UpdateFlags        update_flags,
                                               const Quadrature<dim-1> &quadrature) const
 {
@@ -874,7 +874,7 @@ MappingQGeneric<dim,spacedim>::get_face_data (const UpdateFlags        update_fl
 
 
 template<int dim, int spacedim>
-typename Mapping<dim,spacedim>::InternalDataBase *
+typename MappingQGeneric<dim,spacedim>::InternalData *
 MappingQGeneric<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
                                                  const Quadrature<dim-1>& quadrature) const
 {