Store the correct polynomial degree in MappingQGeneric.

Currently, MappingQ derives from MappingQ1 which passes 1 down as polynomial degree.
This worked 'by accident' because the polynomial degree is only used by MappingQGeneric
in get_data() and friends, which are overloaded by MappingQ. However, the correct
place to store the polynomial degree is clearly in MappingQGeneric.

diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h
index 7f11705c25649e699cbf1a449c20289561e0f180..b954602e3d520582ac79f4d2997166d9262bf963 100644 (file)
--- a/include/deal.II/fe/mapping_q.h
+++ b/include/deal.II/fe/mapping_q.h
@@ -146,12 +146,6 @@ public:
              const typename Mapping<dim,spacedim>::InternalDataBase &internal,
              VectorSlice<std::vector<Tensor<3,spacedim> > >          output) const;
 
-  /**
-   * Return the degree of the mapping, i.e. the value which was passed to the
-   * constructor.
-   */
-  unsigned int get_degree () const;
-
   /**
    * Return a pointer to a copy of the present object. The caller of this copy
    * then assumes ownership of it.
@@ -421,12 +415,6 @@ protected:
                   int,
                   << "laplace_vector not set for degree=" << arg1 << ".");
 
-  /**
-   * Degree @p p of the polynomials used as shape functions for the Qp mapping
-   * of cells at the boundary.
-   */
-  const unsigned int degree;
-
   /**
    * Number of inner mapping shape functions.
    */

diff --git a/include/deal.II/fe/mapping_q1.h b/include/deal.II/fe/mapping_q1.h
index 92943025ace010f8f4309debfd041dbcc29d9d00..4ec4d7b928fe1e705af726dd90a926bdb8e990f7 100644 (file)
--- a/include/deal.II/fe/mapping_q1.h
+++ b/include/deal.II/fe/mapping_q1.h
@@ -98,14 +98,19 @@ public:
    */
   typedef typename MappingQGeneric<dim,spacedim>::InternalData InternalData;
 
-protected:
-
-
   /**
    * @}
    */
 
 protected:
+
+  /**
+   * Constructor. This constructor is for odd purposes: MappingQ is
+   * derived from this class (for historical reasons) and it needs a
+   * way to pass down the "true" polynomial degree of the mapping.
+   */
+  MappingQ1 (const unsigned int degree);
+
   /* Trick to templatize transform_real_to_unit_cell<dim, dim+1> */
   template<int dim_>
   Point<dim_>

diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h
index 5a1e97b2da7e6f59fb1306532667b0bf8885e1c7..3791c7da33d73d3035314f2a1fdc396599cf64fb 100644 (file)
--- a/include/deal.II/fe/mapping_q_generic.h
+++ b/include/deal.II/fe/mapping_q_generic.h
@@ -63,6 +63,12 @@ public:
    */
   MappingQGeneric (const unsigned int polynomial_degree);
 
+  /**
+   * Return the degree of the mapping, i.e. the value which was passed to the
+   * constructor.
+   */
+  unsigned int get_degree () const;
+
   /**
    * Always returns @p true because MappingQ1 preserves vertex locations.
    */

diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc
index 6635f30f90b62045c89c4128f1eddd95158504de..5eb71043b011d77aafa5336e4760ec8cbb41f3dc 100644 (file)
--- a/source/fe/mapping_q.cc
+++ b/source/fe/mapping_q.cc
@@ -54,10 +54,10 @@ MappingQ<dim,spacedim>::InternalData::memory_consumption () const
 
 
 template<int dim, int spacedim>
-MappingQ<dim,spacedim>::MappingQ (const unsigned int p,
+MappingQ<dim,spacedim>::MappingQ (const unsigned int degree,
                                   const bool use_mapping_q_on_all_cells)
   :
-  degree(p),
+  MappingQ1<dim,spacedim>(degree),
   n_inner(Utilities::fixed_power<dim>(degree-1)),
   n_outer((dim==1) ? 2 :
           ((dim==2) ?
@@ -88,7 +88,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);
+  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;
@@ -120,7 +120,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);
+  InternalData *data = new InternalData(this->polynomial_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
@@ -153,7 +153,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);
+  InternalData *data = new InternalData(this->polynomial_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
@@ -219,7 +219,7 @@ fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
   const CellSimilarity::Similarity updated_cell_similarity
     = ((data.use_mapping_q1_on_current_cell == false)
        &&
-       (get_degree() > 1)
+       (this->polynomial_degree > 1)
        ?
        CellSimilarity::invalid_next_cell
        :
@@ -333,14 +333,14 @@ template<int dim, int spacedim>
 void
 MappingQ<dim,spacedim>::set_laplace_on_quad_vector(Table<2,double> &loqvs) const
 {
-  Assert(degree>1, ExcInternalError());
-  const unsigned int n_inner_2d=(degree-1)*(degree-1);
-  const unsigned int n_outer_2d=4+4*(degree-1);
+  Assert(this->polynomial_degree>1, ExcInternalError());
+  const unsigned int n_inner_2d=(this->polynomial_degree-1)*(this->polynomial_degree-1);
+  const unsigned int n_outer_2d=4+4*(this->polynomial_degree-1);
 
   // first check whether we have precomputed the values for some polynomial
   // degree; the sizes of arrays is n_inner_2d*n_outer_2d
   double const *loqv_ptr=0;
-  switch (degree)
+  switch (this->polynomial_degree)
     {
     // for degree==1, we shouldn't have to compute any support points, since
     // all of them are on the vertices
@@ -380,7 +380,7 @@ MappingQ<dim,spacedim>::set_laplace_on_quad_vector(Table<2,double> &loqvs) const
         compute_laplace_vector(loqvs);
       else if (dim == 3)
         {
-          MappingQ<2,2> mapping_2d(this->degree);
+          MappingQ<2,2> mapping_2d(this->polynomial_degree);
           loqvs = mapping_2d.laplace_on_quad_vector;
         }
     }
@@ -389,7 +389,7 @@ MappingQ<dim,spacedim>::set_laplace_on_quad_vector(Table<2,double> &loqvs) const
   // this
   for (unsigned int unit_point=0; unit_point<loqvs.n_rows(); ++unit_point)
     Assert(std::fabs(std::accumulate(loqvs[unit_point].begin(),
-                                     loqvs[unit_point].end(),0.)-1)<1e-13*this->degree,
+                                     loqvs[unit_point].end(),0.)-1)<1e-13*this->polynomial_degree,
            ExcInternalError());
 }
 
@@ -399,12 +399,12 @@ template <>
 void
 MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const
 {
-  Assert(degree>1, ExcInternalError());
+  Assert(this->polynomial_degree>1, ExcInternalError());
 
   // first check whether we have precomputed the values for some polynomial
   // degree
   double const *lohv_ptr=0;
-  if (degree==2)
+  if (this->polynomial_degree==2)
     {
       static const double loqv2[26]
         = {1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128.,
@@ -432,7 +432,7 @@ MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const
   // this
   for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
     Assert(std::fabs(std::accumulate(lohvs[unit_point].begin(),
-                                     lohvs[unit_point].end(),0.) - 1)<1e-13*this->degree,
+                                     lohvs[unit_point].end(),0.) - 1)<1e-13*this->polynomial_degree,
            ExcInternalError());
 }
 
@@ -465,13 +465,13 @@ MappingQ<dim,spacedim>::compute_laplace_vector(Table<2,double> &lvs) const
 
   // for degree==1, we shouldn't have to compute any support points, since all
   // of them are on the vertices
-  Assert(degree>1, ExcInternalError());
+  Assert(this->polynomial_degree>1, ExcInternalError());
 
   // compute the shape gradients at the quadrature points on the unit cell
-  const QGauss<dim> quadrature(degree+1);
+  const QGauss<dim> quadrature(this->polynomial_degree+1);
   const unsigned int n_q_points=quadrature.size();
 
-  InternalData quadrature_data(degree);
+  InternalData quadrature_data(this->polynomial_degree);
   quadrature_data.shape_derivatives.resize(quadrature_data.n_shape_functions *
                                            n_q_points);
   quadrature_data.compute_shape_function_values(quadrature.get_points());
@@ -532,10 +532,10 @@ MappingQ<dim,spacedim>::apply_laplace_vector(const Table<2,double> &lvs,
   // compute the quad laplace vector. this is mentioned in the constructor of
   // this class, although I don't understand the reason for not aborting there
   // any more [WB]
-  Assert(lvs.n_rows()!=0, ExcLaplaceVectorNotSet(degree));
+  Assert(lvs.n_rows()!=0, ExcLaplaceVectorNotSet(this->polynomial_degree));
 
   const unsigned int n_inner_apply=lvs.n_rows();
-  Assert(n_inner_apply==n_inner || n_inner_apply==(degree-1)*(degree-1),
+  Assert(n_inner_apply==n_inner || n_inner_apply==(this->polynomial_degree-1)*(this->polynomial_degree-1),
          ExcInternalError());
   const unsigned int n_outer_apply=lvs.n_cols();
   Assert(a.size()==n_outer_apply,
@@ -588,7 +588,7 @@ MappingQ<dim,spacedim>::compute_support_points_laplace(const typename Triangulat
   for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
     a[i] = cell->vertex(i);
 
-  if (degree>1)
+  if (this->polynomial_degree>1)
     switch (dim)
       {
       case 1:
@@ -625,7 +625,7 @@ MappingQ<dim,spacedim>::add_line_support_points (const typename Triangulation<di
                                                  std::vector<Point<spacedim> > &a) const
 {
   // if we only need the midpoint, then ask for it.
-  if (degree==2)
+  if (this->polynomial_degree==2)
     {
       for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
         {
@@ -645,7 +645,7 @@ MappingQ<dim,spacedim>::add_line_support_points (const typename Triangulation<di
     // otherwise call the more complicated functions and ask for inner points
     // from the boundary description
     {
-      std::vector<Point<spacedim> > line_points (degree-1);
+      std::vector<Point<spacedim> > line_points (this->polynomial_degree-1);
       // loop over each of the lines, and if it is at the boundary, then first
       // get the boundary description and second compute the points on it
       for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
@@ -699,9 +699,9 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
 
   static const StraightBoundary<3> straight_boundary;
   // used if face quad at boundary or entirely in the interior of the domain
-  std::vector<Point<3> > quad_points ((degree-1)*(degree-1));
+  std::vector<Point<3> > quad_points ((this->polynomial_degree-1)*(this->polynomial_degree-1));
   // used if only one line of face quad is at boundary
-  std::vector<Point<3> > b(4*degree);
+  std::vector<Point<3> > b(4*this->polynomial_degree);
 
   // Used by the new Manifold interface. This vector collects the
   // vertices used to compute the intermediate points.
@@ -769,12 +769,12 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
               // call of function apply_laplace_vector increases size of b
               // about 1. There resize b for the case the mentioned function
               // was already called.
-              b.resize(4*degree);
+              b.resize(4*this->polynomial_degree);
 
               // b is of size 4*degree, make sure that this is the right size
-              Assert(b.size()==vertices_per_face+lines_per_face*(degree-1),
+              Assert(b.size()==vertices_per_face+lines_per_face*(this->polynomial_degree-1),
                      ExcDimensionMismatch(b.size(),
-                                          vertices_per_face+lines_per_face*(degree-1)));
+                                          vertices_per_face+lines_per_face*(this->polynomial_degree-1)));
 
               // sort the points into b. We used access from the cell (not
               // from the face) to fill b, so we can assume a standard face
@@ -784,19 +784,19 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
                 b[i]=a[GeometryInfo<3>::face_to_cell_vertices(face_no, i)];
 
               for (unsigned int i=0; i<lines_per_face; ++i)
-                for (unsigned int j=0; j<degree-1; ++j)
-                  b[vertices_per_face+i*(degree-1)+j]=
+                for (unsigned int j=0; j<this->polynomial_degree-1; ++j)
+                  b[vertices_per_face+i*(this->polynomial_degree-1)+j]=
                     a[vertices_per_cell + GeometryInfo<3>::face_to_cell_lines(
-                        face_no, i)*(degree-1)+j];
+                        face_no, i)*(this->polynomial_degree-1)+j];
 
               // Now b includes the support points on the quad and we can
               // apply the laplace vector
               apply_laplace_vector(laplace_on_quad_vector, b);
-              Assert(b.size()==4*degree+(degree-1)*(degree-1),
-                     ExcDimensionMismatch(b.size(), 4*degree+(degree-1)*(degree-1)));
+              Assert(b.size()==4*this->polynomial_degree+(this->polynomial_degree-1)*(this->polynomial_degree-1),
+                     ExcDimensionMismatch(b.size(), 4*this->polynomial_degree+(this->polynomial_degree-1)*(this->polynomial_degree-1)));
 
-              for (unsigned int i=0; i<(degree-1)*(degree-1); ++i)
-                a.push_back(b[4*degree+i]);
+              for (unsigned int i=0; i<(this->polynomial_degree-1)*(this->polynomial_degree-1); ++i)
+                a.push_back(b[4*this->polynomial_degree+i]);
             }
           else
             {
@@ -829,7 +829,7 @@ MappingQ<2,3>::
 add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell,
                         std::vector<Point<3> >                &a) const
 {
-  std::vector<Point<3> > quad_points ((degree-1)*(degree-1));
+  std::vector<Point<3> > quad_points ((this->polynomial_degree-1)*(this->polynomial_degree-1));
   get_intermediate_points_on_object (cell->get_manifold(), cell, quad_points);
   for (unsigned int i=0; i<quad_points.size(); ++i)
     a.push_back(quad_points[i]);
@@ -1083,20 +1083,11 @@ transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_it
 
 
 
-template<int dim, int spacedim>
-unsigned int
-MappingQ<dim,spacedim>::get_degree() const
-{
-  return degree;
-}
-
-
-
 template<int dim, int spacedim>
 Mapping<dim,spacedim> *
 MappingQ<dim,spacedim>::clone () const
 {
-  return new MappingQ<dim,spacedim>(degree);
+  return new MappingQ<dim,spacedim>(this->polynomial_degree);
 }
 
 
@@ -1117,7 +1108,7 @@ MappingQ<dim,spacedim>::get_intermediate_points (const Manifold<dim, spacedim> &
     {
       // If two points are passed, these are the two vertices, and
       // we can only compute degree-1 intermediate points.
-      AssertDimension(n, degree-1);
+      AssertDimension(n, this->polynomial_degree-1);
       for (unsigned int i=0; i<n; ++i)
         {
           const double x = line_support_points.point(i+1)[0];
@@ -1140,7 +1131,7 @@ MappingQ<dim,spacedim>::get_intermediate_points (const Manifold<dim, spacedim> &
       // If four points are passed, these are the two vertices, and
       // we can only compute (degree-1)*(degree-1) intermediate
       // points.
-      AssertDimension(m, degree-1);
+      AssertDimension(m, this->polynomial_degree-1);
 
       for (unsigned int i=0; i<m; ++i)
         {

diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc
index 6bed941c46c63ff9e85d034ef81ae2c3a5ece1f4..98287bf2a9eb41c010039e669027f5a107045fc6 100644 (file)
--- a/source/fe/mapping_q1.cc
+++ b/source/fe/mapping_q1.cc
@@ -48,6 +48,14 @@ MappingQ1<dim,spacedim>::MappingQ1 ()
 
 
 
+template<int dim, int spacedim>
+MappingQ1<dim,spacedim>::MappingQ1 (const unsigned int p)
+  :
+  MappingQGeneric<dim,spacedim> (p)
+{}
+
+
+
 namespace internal
 {
   namespace MappingQ1
@@ -197,13 +205,14 @@ MappingQ1<dim,spacedim>::transform_unit_to_real_cell (
   const typename Triangulation<dim,spacedim>::cell_iterator &cell,
   const Point<dim> &p) const
 {
-  // Use the get_data function to create an InternalData with data
-  // vectors of the right size and transformation shape values already
-  // computed at point p.
   const Quadrature<dim> point_quadrature(p);
 
-  std_cxx11::unique_ptr<InternalData> mdata (get_data(update_quadrature_points | update_jacobians,
-                                                      point_quadrature));
+  //TODO: Use get_data() here once MappingQ is no longer derived from
+  //MappingQ1. this doesn't currently work because we here really need
+  //a Q1 InternalData, but MappingQGeneric produces one with the
+  //polynomial degree of the MappingQ
+  std_cxx11::unique_ptr<InternalData> mdata (new InternalData(1));
+  mdata->initialize (this->requires_update_flags (update_quadrature_points | update_jacobians), point_quadrature, 1);
 
   // compute the mapping support
   // points
@@ -476,20 +485,18 @@ transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_it
 // it shouldn't really make any difference...
 //      initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
 
-      // Use the get_data function to
-      // create an InternalData with data
-      // vectors of the right size and
-      // transformation shape values and
-      // derivatives already computed at
-      // point initial_p_unit.
       const Quadrature<dim> point_quadrature(initial_p_unit);
 
       UpdateFlags update_flags = update_quadrature_points | update_jacobians;
       if (spacedim>dim)
         update_flags |= update_jacobian_grads;
 
-      std_cxx11::unique_ptr<InternalData> mdata(MappingQ1<dim,spacedim>::get_data(update_flags,
-                                                point_quadrature));
+      //TODO: Use get_data() here once MappingQ is no longer derived from
+      //MappingQ1. this doesn't currently work because we here really need
+      //a Q1 InternalData, but MappingQGeneric produces one with the
+      //polynomial degree of the MappingQ
+      std_cxx11::unique_ptr<InternalData> mdata (new InternalData(1));
+      mdata->initialize (this->requires_update_flags (update_flags), point_quadrature, 1);
 
       compute_mapping_support_points (cell, mdata->mapping_support_points);
       // The support points have to be at

diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc
index 3d8df89ddffaba0133bd1b053264cb4d38e0ec1d..53534ffb3bc1e8dbf09e7d1f3165582a0242528f 100644 (file)
--- a/source/fe/mapping_q_generic.cc
+++ b/source/fe/mapping_q_generic.cc
@@ -776,6 +776,15 @@ MappingQGeneric<dim,spacedim>::MappingQGeneric (const unsigned int p)
 
 
 
+template<int dim, int spacedim>
+unsigned int
+MappingQGeneric<dim,spacedim>::get_degree() const
+{
+  return polynomial_degree;
+}
+
+
+
 template<int dim, int spacedim>
 UpdateFlags
 MappingQGeneric<dim,spacedim>::requires_update_flags (const UpdateFlags in) const