]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move several functions from the public interface of MappingQ into the .cc file.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 12 Sep 2015 01:07:14 +0000 (20:07 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sat, 12 Sep 2015 10:51:51 +0000 (05:51 -0500)
include/deal.II/fe/mapping_q.h
source/fe/mapping_q.cc

index 96bce39faff817c84c518c39ede1041c41b8d332..19af543546fe875cc40a7699587705ddbe079e70 100644 (file)
@@ -163,68 +163,6 @@ public:
   Mapping<dim,spacedim> *clone () const;
 
 
-protected:
-
-  /**
-   * For <tt>dim=2,3</tt>. Append the support points of all shape functions
-   * located on bounding lines to the vector @p a. Points located on the line
-   * but not on vertices are not included.
-   *
-   * Needed by the @p compute_support_points_laplace function . For
-   * <tt>dim=1</tt> this function is empty.
-   *
-   * This function is made virtual in order to allow derived classes to choose
-   * shape function support points differently than the present class, which
-   * chooses the points as interpolation points on the boundary.
-   */
-  virtual void
-  add_line_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-                           std::vector<Point<spacedim> > &a) const;
-
-  /**
-   * For <tt>dim=3</tt>. Append the support points of all shape functions
-   * located on bounding faces (quads in 3d) to the vector @p a. Points
-   * located on the quad but not on vertices are not included.
-   *
-   * Needed by the @p compute_support_points_laplace function. For
-   * <tt>dim=1</tt> and <tt>dim=2</tt> this function is empty.
-   *
-   * This function is made virtual in order to allow derived classes to choose
-   * shape function support points differently than the present class, which
-   * chooses the points as interpolation points on the boundary.
-   */
-  virtual void
-  add_quad_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-                          std::vector<Point<spacedim> > &a) const;
-
-
-private:
-  /**
-   * Ask the manifold descriptor to return intermediate points on lines or
-   * faces. The function needs to return one or multiple points (depending on
-   * the number of elements in the output vector @p points that lie inside a
-   * line, quad or hex). Whether it is a line, quad or hex doesn't really
-   * matter to this function but it can be inferred from the number of input
-   * points in the @p surrounding_points vector.
-   */
-  void get_intermediate_points(const Manifold<dim, spacedim> &manifold,
-                               const std::vector<Point<spacedim> > &surrounding_points,
-                               std::vector<Point<spacedim> > &points) const;
-
-
-  /**
-   * Ask the manifold descriptor to return intermediate points on the object
-   * pointed to by the TriaIterator @p iter. This function tries to be
-   * backward compatible with respect to the differences between
-   * Boundary<dim,spacedim> and Manifold<dim,spacedim>, querying the first
-   * whenever the passed @p manifold can be upgraded to a
-   * Boundary<dim,spacedim>.
-   */
-  template <class TriaIterator>
-  void get_intermediate_points_on_object(const Manifold<dim, spacedim> &manifold,
-                                         const TriaIterator &iter,
-                                         std::vector<Point<spacedim> > &points) const;
-
   /**
    * @name Interface with FEValues
    * @{
@@ -332,90 +270,106 @@ protected:
    * This function is needed by the constructor of
    * <tt>MappingQ<dim,spacedim></tt> for <tt>dim=</tt> 2 and 3.
    *
-   * For <tt>degree<4</tt> this function sets the @p laplace_on_quad_vector to
+   * For <tt>degree<4</tt> this function sets the @p support_point_weights_on_quad to
    * the hardcoded data. For <tt>degree>=4</tt> and MappingQ<2> this vector is
    * computed.
    *
-   * For the definition of the @p laplace_on_quad_vector please refer to
+   * For the definition of the @p support_point_weights_on_quad please refer to
    * equation (8) of the `mapping' report.
    */
   void
-  set_laplace_on_quad_vector(Table<2,double> &loqvs) const;
+  set_support_point_weights_on_quad(Table<2,double> &loqvs) const;
 
   /**
    * This function is needed by the constructor of <tt>MappingQ<3></tt>.
    *
-   * For <tt>degree==2</tt> this function sets the @p laplace_on_hex_vector to
+   * For <tt>degree==2</tt> this function sets the @p support_point_weights_on_hex to
    * the hardcoded data. For <tt>degree>2</tt> this vector is computed.
    *
-   * For the definition of the @p laplace_on_hex_vector please refer to
+   * For the definition of the @p support_point_weights_on_hex please refer to
    * equation (8) of the `mapping' report.
    */
-  void set_laplace_on_hex_vector(Table<2,double> &lohvs) const;
+  void set_support_point_weights_on_hex(Table<2,double> &lohvs) const;
 
   /**
-   * Computes the <tt>laplace_on_quad(hex)_vector</tt>.
+   * Compute the <tt>support_point_weights_on_quad(hex)_vector</tt>.
    *
-   * Called by the <tt>set_laplace_on_quad(hex)_vector</tt> functions if the
+   * Called by the <tt>set_support_point_weights_on_quad(hex)_vector</tt> functions if the
    * data is not yet hardcoded.
    *
-   * For the definition of the <tt>laplace_on_quad(hex)_vector</tt> please
+   * For the definition of the <tt>support_point_weights_on_quad(hex)_vector</tt> please
    * refer to equation (8) of the `mapping' report.
    */
   void compute_laplace_vector(Table<2,double> &lvs) const;
 
   /**
-   * Takes a <tt>laplace_on_hex(quad)_vector</tt> and applies it to the vector
-   * @p a to compute the inner support points as a linear combination of the
-   * exterior points.
-   *
-   * The vector @p a initially contains the locations of the @p n_outer
-   * points, the @p n_inner computed inner points are appended.
-   *
-   * See equation (7) of the `mapping' report.
+   * Compute the support points of the mapping. Interior support
+   * points (ie. support points in quads for 2d, in hexes for 3d) are
+   * computed using the solution of a Laplace equation with the
+   * position of the outer support points as boundary values, in order
+   * to make the transformation as smooth as possible.
    */
-  void apply_laplace_vector(const Table<2,double>   &lvs,
-                            std::vector<Point<spacedim> > &a) const;
+  virtual
+  void
+  compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                                 std::vector<Point<spacedim> > &a) const;
+
 
   /**
-   * Computes the support points of the mapping.
-   */
-  virtual void compute_mapping_support_points(
-    const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-    std::vector<Point<spacedim> > &a) const;
+  * For <tt>dim=2,3</tt>. Append the support points of all shape functions
+  * located on bounding lines to the vector @p a. Points located on the line
+  * but not on vertices are not included.
+  *
+  * Needed by the @p compute_support_points_laplace function . For
+  * <tt>dim=1</tt> this function is empty.
+  *
+  * This function is made virtual in order to allow derived classes to choose
+  * shape function support points differently than the present class, which
+  * chooses the points as interpolation points on the boundary.
+  */
+  virtual
+  void
+  add_line_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                           std::vector<Point<spacedim> > &a) const;
 
   /**
-   * Computes all support points of the mapping shape functions. The inner
-   * support points (ie. support points in quads for 2d, in hexes for 3d) are
-   * computed using the solution of a Laplace equation with the position of
-   * the outer support points as boundary values, in order to make the
-   * transformation as smooth as possible.
+   * For <tt>dim=3</tt>. Append the support points of all shape functions
+   * located on bounding faces (quads in 3d) to the vector @p a. Points
+   * located on the quad but not on vertices are not included.
+   *
+   * Needed by the @p compute_support_points_laplace function. For
+   * <tt>dim=1</tt> and <tt>dim=2</tt> this function is empty.
+   *
+   * This function is made virtual in order to allow derived classes to choose
+   * shape function support points differently than the present class, which
+   * chooses the points as interpolation points on the boundary.
    */
-  void compute_support_points_laplace(
-    const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-    std::vector<Point<spacedim> > &a) const;
+  virtual
+  void
+  add_quad_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                          std::vector<Point<spacedim> > &a) const;
 
   /**
-   * Needed by the @p laplace_on_quad function (for <tt>dim==2</tt>). Filled
+   * Needed by the @p support_point_weights_on_quad function (for <tt>dim==2</tt>). Filled
    * by the constructor.
    *
-   * Sizes: laplace_on_quad_vector.size()= number of inner unit_support_points
-   * laplace_on_quad_vector[i].size()= number of outer unit_support_points,
+   * Sizes: support_point_weights_on_quad.size()= number of inner unit_support_points
+   * support_point_weights_on_quad[i].size()= number of outer unit_support_points,
    * i.e.  unit_support_points on the boundary of the quad
    *
    * For the definition of this vector see equation (8) of the `mapping'
    * report.
    */
-  Table<2,double> laplace_on_quad_vector;
+  Table<2,double> support_point_weights_on_quad;
 
   /**
-   * Needed by the @p laplace_on_hex function (for <tt>dim==3</tt>). Filled by
+   * Needed by the @p support_point_weights_on_hex function (for <tt>dim==3</tt>). Filled by
    * the constructor.
    *
    * For the definition of this vector see equation (8) of the `mapping'
    * report.
    */
-  Table<2,double> laplace_on_hex_vector;
+  Table<2,double> support_point_weights_on_hex;
 
   /**
    * Exception.
@@ -490,63 +444,15 @@ protected:
 #ifndef DOXYGEN
 
 template <>
-void MappingQ<1>::set_laplace_on_quad_vector(Table<2,double> &) const;
+void MappingQ<1>::set_support_point_weights_on_quad(Table<2,double> &) const;
 
 template <>
-void MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const;
+void MappingQ<3>::set_support_point_weights_on_hex(Table<2,double> &lohvs) const;
 
 template <>
 void MappingQ<1>::compute_laplace_vector(Table<2,double> &) const;
 
 
-template<>
-void MappingQ<3>::add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
-                                          std::vector<Point<3> >                &a) const;
-
-// ---- Templated functions ---- //
-template <int dim, int spacedim>
-template <class TriaIterator>
-void
-MappingQ<dim,spacedim>::get_intermediate_points_on_object(const Manifold<dim, spacedim> &manifold,
-                                                          const TriaIterator &iter,
-                                                          std::vector<Point<spacedim> > &points) const
-{
-  const unsigned int structdim = TriaIterator::AccessorType::structure_dimension;
-
-  // Try backward compatibility option.
-  if (const Boundary<dim,spacedim> *boundary
-      = dynamic_cast<const Boundary<dim,spacedim> *>(&manifold))
-    // This is actually a boundary. Call old methods.
-    {
-      switch (structdim)
-        {
-        case 1:
-        {
-          const typename Triangulation<dim,spacedim>::line_iterator line = iter;
-          boundary->get_intermediate_points_on_line(line, points);
-          return;
-        }
-        case 2:
-        {
-          const typename Triangulation<dim,spacedim>::quad_iterator quad = iter;
-          boundary->get_intermediate_points_on_quad(quad, points);
-          return;
-        }
-        default:
-          Assert(false, ExcInternalError());
-          return;
-        }
-    }
-  else
-    {
-      std::vector<Point<spacedim> > sp(GeometryInfo<structdim>::vertices_per_cell);
-      for (unsigned int i=0; i<sp.size(); ++i)
-        sp[i] = iter->vertex(i);
-      get_intermediate_points(manifold, sp, points);
-    }
-}
-
-
 #endif // DOXYGEN
 
 DEAL_II_NAMESPACE_CLOSE
index 13443f4338c871fb3cdffc6b2633273a142fbf51..9db5d010b29c8668ac80f409efbfa669b98811fe 100644 (file)
@@ -83,13 +83,13 @@ MappingQ<dim,spacedim>::MappingQ (const unsigned int degree,
   Assert(n_inner+n_outer==Utilities::fixed_power<dim>(degree+1),
          ExcInternalError());
 
-  // build laplace_on_quad_vector
+  // build support_point_weights_on_quad
   if (degree>1)
     {
       if (dim >= 2)
-        set_laplace_on_quad_vector(laplace_on_quad_vector);
+        set_support_point_weights_on_quad(support_point_weights_on_quad);
       if (dim >= 3)
-        set_laplace_on_hex_vector(laplace_on_hex_vector);
+        set_support_point_weights_on_hex(support_point_weights_on_hex);
     }
 }
 
@@ -111,13 +111,13 @@ MappingQ<dim,spacedim>::MappingQ (const MappingQ<dim,spacedim> &mapping)
   Assert(n_inner+n_outer==Utilities::fixed_power<dim>(this->polynomial_degree+1),
          ExcInternalError());
 
-  // build laplace_on_quad_vector
+  // build support_point_weights_on_quad
   if (this->polynomial_degree>1)
     {
       if (dim >= 2)
-        set_laplace_on_quad_vector(laplace_on_quad_vector);
+        set_support_point_weights_on_quad(support_point_weights_on_quad);
       if (dim >= 3)
-        set_laplace_on_hex_vector(laplace_on_hex_vector);
+        set_support_point_weights_on_hex(support_point_weights_on_hex);
     }
 }
 
@@ -341,7 +341,7 @@ fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterato
 
 template <>
 void
-MappingQ<1>::set_laplace_on_quad_vector(Table<2,double> &) const
+MappingQ<1>::set_support_point_weights_on_quad(Table<2,double> &) const
 {
   Assert(false, ExcInternalError());
 }
@@ -350,7 +350,7 @@ MappingQ<1>::set_laplace_on_quad_vector(Table<2,double> &) const
 
 template<int dim, int spacedim>
 void
-MappingQ<dim,spacedim>::set_laplace_on_quad_vector(Table<2,double> &loqvs) const
+MappingQ<dim,spacedim>::set_support_point_weights_on_quad(Table<2,double> &loqvs) const
 {
   Assert(this->polynomial_degree>1, ExcInternalError());
   const unsigned int n_inner_2d=(this->polynomial_degree-1)*(this->polynomial_degree-1);
@@ -400,7 +400,7 @@ MappingQ<dim,spacedim>::set_laplace_on_quad_vector(Table<2,double> &loqvs) const
       else if (dim == 3)
         {
           MappingQ<2,2> mapping_2d(this->polynomial_degree);
-          loqvs = mapping_2d.laplace_on_quad_vector;
+          loqvs = mapping_2d.support_point_weights_on_quad;
         }
     }
 
@@ -416,7 +416,7 @@ MappingQ<dim,spacedim>::set_laplace_on_quad_vector(Table<2,double> &loqvs) const
 
 template <>
 void
-MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const
+MappingQ<3>::set_support_point_weights_on_hex(Table<2,double> &lohvs) const
 {
   Assert(this->polynomial_degree>1, ExcInternalError());
 
@@ -459,7 +459,7 @@ MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const
 
 template<int dim, int spacedim>
 void
-MappingQ<dim,spacedim>::set_laplace_on_hex_vector(Table<2,double> &) const
+MappingQ<dim,spacedim>::set_support_point_weights_on_hex(Table<2,double> &) const
 {
   Assert(false, ExcInternalError());
 }
@@ -541,107 +541,176 @@ MappingQ<dim,spacedim>::compute_laplace_vector(Table<2,double> &lvs) const
 
 
 
-template<int dim, int spacedim>
-void
-MappingQ<dim,spacedim>::apply_laplace_vector(const Table<2,double> &lvs,
-                                             std::vector<Point<spacedim> > &a) const
-{
-  // check whether the data we need is really available. if you fail here and
-  // if lvs==laplace_on_quad_vector in the calling function, then we didn't
-  // 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(this->polynomial_degree));
-
-  const unsigned int n_inner_apply=lvs.n_rows();
-  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,
-         ExcDimensionMismatch(a.size(), n_outer_apply));
-
-  // compute each inner point as linear combination of the outer points. the
-  // weights are given by the lvs entries, the outer points are the first
-  // (existing) elements of a
-  for (unsigned int unit_point=0; unit_point<n_inner_apply; ++unit_point)
-    {
-      Assert(lvs.n_cols()==n_outer_apply, ExcInternalError());
-      Point<spacedim> p;
-      for (unsigned int k=0; k<n_outer_apply; ++k)
-        p+=lvs[unit_point][k]*a[k];
-
-      a.push_back(p);
-    }
-}
-
-
-template<int dim, int spacedim>
-void
-MappingQ<dim,spacedim>::compute_mapping_support_points(
-  const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-  std::vector<Point<spacedim> > &a) const
-{
-  // if this is a cell for which we want to compute the full mapping, then get
-  // them from the following function
-  if (use_mapping_q_on_all_cells || cell->has_boundary_lines())
-    compute_support_points_laplace(cell, a);
-  else
-    // otherwise: use a Q1 mapping for which the mapping shape function
-    // support points are simply the vertices of the cell
-    {
-      a.resize(GeometryInfo<dim>::vertices_per_cell);
-
-      for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-        a[i] = cell->vertex(i);
-    }
-}
-
-
-template<int dim, int spacedim>
-void
-MappingQ<dim,spacedim>::compute_support_points_laplace(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-                                                       std::vector<Point<spacedim> > &a) const
+namespace
 {
-  // in any case, we need the vertices first
-  a.resize(GeometryInfo<dim>::vertices_per_cell);
-  for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-    a[i] = cell->vertex(i);
-
-  if (this->polynomial_degree>1)
-    switch (dim)
+  /**
+   * Ask the manifold descriptor to return intermediate points on lines or
+   * faces. The function needs to return one or multiple points (depending on
+   * the number of elements in the output vector @p points that lie inside a
+   * line, quad or hex). Whether it is a line, quad or hex doesn't really
+   * matter to this function but it can be inferred from the number of input
+   * points in the @p surrounding_points vector.
+   */
+  template<int dim, int spacedim>
+  void
+  get_intermediate_points (const Manifold<dim, spacedim> &manifold,
+                           const QGaussLobatto<1>        &line_support_points,
+                           const std::vector<Point<spacedim> > &surrounding_points,
+                           std::vector<Point<spacedim> > &points)
+  {
+    Assert(surrounding_points.size() >= 2, ExcMessage("At least 2 surrounding points are required"));
+    const unsigned int n=points.size();
+    Assert(n>0, ExcMessage("You can't ask for 0 intermediate points."));
+    std::vector<double> w(surrounding_points.size());
+
+    switch (surrounding_points.size())
       {
-      case 1:
-        add_line_support_points(cell, a);
-        break;
       case 2:
-        // in 2d, add the points on the four bounding lines to the exterior
-        // (outer) points
-        add_line_support_points (cell, a);
-        if (dim != spacedim)
-          add_quad_support_points(cell, a);
-        else
-          apply_laplace_vector (laplace_on_quad_vector,a);
+      {
+        // If two points are passed, these are the two vertices, and
+        // we can only compute degree-1 intermediate points.
+        for (unsigned int i=0; i<n; ++i)
+          {
+            const double x = line_support_points.point(i+1)[0];
+            w[1] = x;
+            w[0] = (1-x);
+            Quadrature<spacedim> quadrature(surrounding_points, w);
+            points[i] = manifold.get_new_point(quadrature);
+          }
         break;
+      }
 
-      case 3:
+      case 4:
       {
-        // in 3d also add the points located on the boundary faces
-        add_line_support_points (cell, a);
-        add_quad_support_points (cell, a);
-        apply_laplace_vector (laplace_on_hex_vector, a);
+        Assert(spacedim >= 2, ExcImpossibleInDim(spacedim));
+        const unsigned m=
+          static_cast<unsigned int>(std::sqrt(static_cast<double>(n)));
+        // is n a square number
+        Assert(m*m==n, ExcInternalError());
+
+        // If four points are passed, these are the two vertices, and
+        // we can only compute (degree-1)*(degree-1) intermediate
+        // points.
+        for (unsigned int i=0; i<m; ++i)
+          {
+            const double y=line_support_points.point(1+i)[0];
+            for (unsigned int j=0; j<m; ++j)
+              {
+                const double x=line_support_points.point(1+j)[0];
+
+                w[0] = (1-x)*(1-y);
+                w[1] =     x*(1-y);
+                w[2] = (1-x)*y    ;
+                w[3] =     x*y    ;
+                Quadrature<spacedim> quadrature(surrounding_points, w);
+                points[i*m+j]=manifold.get_new_point(quadrature);
+              }
+          }
         break;
       }
-      default:
+
+      case 8:
         Assert(false, ExcNotImplemented());
         break;
+      default:
+        Assert(false, ExcInternalError());
+        break;
+      }
+  }
+
+
+
+
+  /**
+   * Ask the manifold descriptor to return intermediate points on the object
+   * pointed to by the TriaIterator @p iter. This function tries to be
+   * backward compatible with respect to the differences between
+   * Boundary<dim,spacedim> and Manifold<dim,spacedim>, querying the first
+   * whenever the passed @p manifold can be upgraded to a
+   * Boundary<dim,spacedim>.
+   */
+  template <int dim, int spacedim, class TriaIterator>
+  void get_intermediate_points_on_object(const Manifold<dim, spacedim> &manifold,
+                                         const QGaussLobatto<1>        &line_support_points,
+                                         const TriaIterator &iter,
+                                         std::vector<Point<spacedim> > &points)
+  {
+    const unsigned int structdim = TriaIterator::AccessorType::structure_dimension;
+
+    // Try backward compatibility option.
+    if (const Boundary<dim,spacedim> *boundary
+        = dynamic_cast<const Boundary<dim,spacedim> *>(&manifold))
+      // This is actually a boundary. Call old methods.
+      {
+        switch (structdim)
+          {
+          case 1:
+          {
+            const typename Triangulation<dim,spacedim>::line_iterator line = iter;
+            boundary->get_intermediate_points_on_line(line, points);
+            return;
+          }
+          case 2:
+          {
+            const typename Triangulation<dim,spacedim>::quad_iterator quad = iter;
+            boundary->get_intermediate_points_on_quad(quad, points);
+            return;
+          }
+          default:
+            Assert(false, ExcInternalError());
+            return;
+          }
       }
+    else
+      {
+        std::vector<Point<spacedim> > sp(GeometryInfo<structdim>::vertices_per_cell);
+        for (unsigned int i=0; i<sp.size(); ++i)
+          sp[i] = iter->vertex(i);
+        get_intermediate_points(manifold, line_support_points, sp, points);
+      }
+  }
+
+
+  /**
+   * Take a <tt>support_point_weights_on_hex(quad)</tt> and apply it to the vector
+   * @p a to compute the inner support points as a linear combination of the
+   * exterior points.
+   *
+   * The vector @p a initially contains the locations of the @p n_outer
+   * points, the @p n_inner computed inner points are appended.
+   *
+   * See equation (7) of the `mapping' report.
+   */
+  template <int spacedim>
+  void add_weighted_interior_points(const Table<2,double>   &lvs,
+                                    std::vector<Point<spacedim> > &a)
+  {
+    const unsigned int n_inner_apply=lvs.n_rows();
+    const unsigned int n_outer_apply=lvs.n_cols();
+    Assert(a.size()==n_outer_apply,
+           ExcDimensionMismatch(a.size(), n_outer_apply));
+
+    // compute each inner point as linear combination of the outer points. the
+    // weights are given by the lvs entries, the outer points are the first
+    // (existing) elements of a
+    for (unsigned int unit_point=0; unit_point<n_inner_apply; ++unit_point)
+      {
+        Assert(lvs.n_cols()==n_outer_apply, ExcInternalError());
+        Point<spacedim> p;
+        for (unsigned int k=0; k<n_outer_apply; ++k)
+          p+=lvs[unit_point][k]*a[k];
+
+        a.push_back(p);
+      }
+  }
 }
 
 
-template<int dim, int spacedim>
+template <int dim, int spacedim>
 void
-MappingQ<dim,spacedim>::add_line_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-                                                 std::vector<Point<spacedim> > &a) const
+MappingQ<dim,spacedim>::
+add_line_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                         std::vector<Point<spacedim> > &a) const
 {
   // if we only need the midpoint, then ask for it.
   if (this->polynomial_degree==2)
@@ -655,10 +724,13 @@ MappingQ<dim,spacedim>::add_line_support_points (const typename Triangulation<di
 
           const Manifold<dim,spacedim> &manifold =
             ( ( line->manifold_id() == numbers::invalid_manifold_id ) &&
-              ( dim < spacedim ) ? cell->get_manifold() :
+              ( dim < spacedim )
+              ?
+              cell->get_manifold()
+              :
               line->get_manifold() );
           a.push_back(manifold.get_new_point_on_line(line));
-        };
+        }
     }
   else
     // otherwise call the more complicated functions and ask for inner points
@@ -683,7 +755,7 @@ MappingQ<dim,spacedim>::add_line_support_points (const typename Triangulation<di
               cell->get_manifold() :
               line->get_manifold() );
 
-          get_intermediate_points_on_object (manifold, line, line_points);
+          get_intermediate_points_on_object (manifold, line_support_points, line, line_points);
 
           if (dim==3)
             {
@@ -698,16 +770,15 @@ MappingQ<dim,spacedim>::add_line_support_points (const typename Triangulation<di
             // in 2D, lines always have the correct orientation. simply append
             // all points
             a.insert (a.end(), line_points.begin(), line_points.end());
-
         }
     }
 }
 
 
 
-template<>
+template <>
 void
-MappingQ<3>::
+MappingQ<3,3>::
 add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
                         std::vector<Point<3> >                &a) const
 {
@@ -718,9 +789,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 ((this->polynomial_degree-1)*(this->polynomial_degree-1));
+  std::vector<Point<3> > quad_points ((polynomial_degree-1)*(polynomial_degree-1));
   // used if only one line of face quad is at boundary
-  std::vector<Point<3> > b(4*this->polynomial_degree);
+  std::vector<Point<3> > b(4*polynomial_degree);
 
   // Used by the new Manifold interface. This vector collects the
   // vertices used to compute the intermediate points.
@@ -758,12 +829,12 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
       // points on it
       if (face->at_boundary())
         {
-          get_intermediate_points_on_object(face->get_manifold(), face, quad_points);
+          get_intermediate_points_on_object(face->get_manifold(), line_support_points, face, quad_points);
 
           // in 3D, the orientation, flip and rotation of the face might not
           // match what we expect here, namely the standard orientation. thus
           // reorder points accordingly. since a Mapping uses the same shape
-          // function as an FEQ, we can ask a FEQ to do the reordering for us.
+          // function as an FE_Q, we can ask a FE_Q to do the reordering for us.
           for (unsigned int i=0; i<quad_points.size(); ++i)
             a.push_back(quad_points[feq.adjust_quad_dof_index_for_face_orientation(i,
                                     face_orientation,
@@ -785,15 +856,15 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
           // boundary, then collect points separately
           if (lines_at_boundary>0)
             {
-              // call of function apply_laplace_vector increases size of b
+              // call of function add_weighted_interior_points increases size of b
               // about 1. There resize b for the case the mentioned function
               // was already called.
-              b.resize(4*this->polynomial_degree);
+              b.resize(4*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*(this->polynomial_degree-1),
+              Assert(b.size()==vertices_per_face+lines_per_face*(polynomial_degree-1),
                      ExcDimensionMismatch(b.size(),
-                                          vertices_per_face+lines_per_face*(this->polynomial_degree-1)));
+                                          vertices_per_face+lines_per_face*(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
@@ -803,19 +874,20 @@ 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<this->polynomial_degree-1; ++j)
-                  b[vertices_per_face+i*(this->polynomial_degree-1)+j]=
+                for (unsigned int j=0; j<polynomial_degree-1; ++j)
+                  b[vertices_per_face+i*(polynomial_degree-1)+j]=
                     a[vertices_per_cell + GeometryInfo<3>::face_to_cell_lines(
-                        face_no, i)*(this->polynomial_degree-1)+j];
+                        face_no, i)*(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*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)));
+              add_weighted_interior_points (support_point_weights_on_quad, b);
+              AssertDimension (b.size(),
+                               4*this->polynomial_degree +
+                               (this->polynomial_degree-1)*(this->polynomial_degree-1));
 
-              for (unsigned int i=0; i<(this->polynomial_degree-1)*(this->polynomial_degree-1); ++i)
-                a.push_back(b[4*this->polynomial_degree+i]);
+              for (unsigned int i=0; i<(polynomial_degree-1)*(polynomial_degree-1); ++i)
+                a.push_back(b[4*polynomial_degree+i]);
             }
           else
             {
@@ -824,7 +896,7 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
               vertices.resize(4);
               for (unsigned int i=0; i<4; ++i)
                 vertices[i] = face->vertex(i);
-              get_intermediate_points (face->get_manifold(), vertices, quad_points);
+              get_intermediate_points (face->get_manifold(), line_support_points, vertices, quad_points);
               // in 3D, the orientation, flip and rotation of the face might
               // not match what we expect here, namely the standard
               // orientation. thus reorder points accordingly. since a Mapping
@@ -842,27 +914,79 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
 
 
 
-template<>
+template <>
 void
 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 ((this->polynomial_degree-1)*(this->polynomial_degree-1));
-  get_intermediate_points_on_object (cell->get_manifold(), cell, quad_points);
+  std::vector<Point<3> > quad_points ((polynomial_degree-1)*(polynomial_degree-1));
+  get_intermediate_points_on_object (cell->get_manifold(), line_support_points,
+                                     cell, quad_points);
   for (unsigned int i=0; i<quad_points.size(); ++i)
     a.push_back(quad_points[i]);
 }
 
 
 
-template<int dim, int spacedim>
+template <int dim, int spacedim>
 void
 MappingQ<dim,spacedim>::
 add_quad_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &,
                         std::vector<Point<spacedim> > &) const
 {
-  Assert (dim > 2, ExcImpossibleInDim(dim));
+  Assert (false, ExcInternalError());
+}
+
+
+
+
+
+template<int dim, int spacedim>
+void
+MappingQ<dim,spacedim>::
+compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                               std::vector<Point<spacedim> > &a) const
+{
+  // get the vertices first
+  a.resize(GeometryInfo<dim>::vertices_per_cell);
+  for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+    a[i] = cell->vertex(i);
+
+  if (this->polynomial_degree>1)
+    switch (dim)
+      {
+      case 1:
+        add_line_support_points(cell, a);
+        break;
+      case 2:
+        // in 2d, add the points on the four bounding lines to the exterior
+        // (outer) points
+        add_line_support_points(cell, a);
+
+        // then get the support points on the quad if we are on a
+        // manifold, otherwise compute them from the points around it
+        if (dim != spacedim)
+          add_quad_support_points(cell, a);
+        else
+          add_weighted_interior_points (support_point_weights_on_quad, a);
+        break;
+
+      case 3:
+      {
+        // in 3d also add the points located on the boundary faces
+        add_line_support_points (cell, a);
+        add_quad_support_points (cell, a);
+
+        // then compute the interior points
+        add_weighted_interior_points (support_point_weights_on_hex, a);
+        break;
+      }
+
+      default:
+        Assert(false, ExcNotImplemented());
+        break;
+      }
 }
 
 
@@ -1093,76 +1217,6 @@ MappingQ<dim,spacedim>::clone () const
 }
 
 
-template<int dim, int spacedim>
-void
-MappingQ<dim,spacedim>::get_intermediate_points (const Manifold<dim, spacedim> &manifold,
-                                                 const std::vector<Point<spacedim> > &surrounding_points,
-                                                 std::vector<Point<spacedim> > &points) const
-{
-  Assert(surrounding_points.size() >= 2, ExcMessage("At least 2 surrounding points are required"));
-  const unsigned int n=points.size();
-  Assert(n>0, ExcMessage("You can't ask for 0 intermediate points."));
-  std::vector<double> w(surrounding_points.size());
-
-  switch (surrounding_points.size())
-    {
-    case 2:
-    {
-      // If two points are passed, these are the two vertices, and
-      // we can only compute degree-1 intermediate points.
-      AssertDimension(n, this->polynomial_degree-1);
-      for (unsigned int i=0; i<n; ++i)
-        {
-          const double x = line_support_points.point(i+1)[0];
-          w[1] = x;
-          w[0] = (1-x);
-          Quadrature<spacedim> quadrature(surrounding_points, w);
-          points[i] = manifold.get_new_point(quadrature);
-        }
-      break;
-    }
-
-    case 4:
-    {
-      Assert(spacedim >= 2, ExcImpossibleInDim(spacedim));
-      const unsigned m=
-        static_cast<unsigned int>(std::sqrt(static_cast<double>(n)));
-      // is n a square number
-      Assert(m*m==n, ExcInternalError());
-
-      // If four points are passed, these are the two vertices, and
-      // we can only compute (degree-1)*(degree-1) intermediate
-      // points.
-      AssertDimension(m, this->polynomial_degree-1);
-
-      for (unsigned int i=0; i<m; ++i)
-        {
-          const double y=line_support_points.point(1+i)[0];
-          for (unsigned int j=0; j<m; ++j)
-            {
-              const double x=line_support_points.point(1+j)[0];
-
-              w[0] = (1-x)*(1-y);
-              w[1] =     x*(1-y);
-              w[2] = (1-x)*y    ;
-              w[3] =     x*y    ;
-              Quadrature<spacedim> quadrature(surrounding_points, w);
-              points[i*m+j]=manifold.get_new_point(quadrature);
-            }
-        }
-      break;
-    }
-
-    case 8:
-      Assert(false, ExcNotImplemented());
-      break;
-    default:
-      Assert(false, ExcInternalError());
-      break;
-    }
-}
-
-
 
 // explicit instantiations
 #include "mapping_q.inst"

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.