]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move compute_mapping_support_points() and friends to MappingQGeneric.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 13 Sep 2015 18:49:23 +0000 (13:49 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 17 Sep 2015 19:45:28 +0000 (14:45 -0500)
include/deal.II/fe/mapping_q.h
include/deal.II/fe/mapping_q_generic.h
source/fe/mapping_q.cc
source/fe/mapping_q_generic.cc

index 72e8678c09bd66f321daf12971737ed13404ecd0..81f16ea5bf4180bcf2c5dfb9297d8d60d01e38ef 100644 (file)
 
 
 #include <deal.II/base/config.h>
-#include <deal.II/base/table.h>
 #include <deal.II/fe/mapping_q1.h>
 #include <deal.II/fe/mapping_q_generic.h>
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/grid/manifold.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -266,94 +263,6 @@ protected:
 
 protected:
 
-  /**
-   * 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.
-   *
-   * The function works its way from the vertices (which it takes from
-   * the given cell) via the support points on the line (for which it
-   * calls the add_line_support_points() function) and the support
-   * points on the quad faces (in 3d, for which it calls the
-   * add_quad_support_points() function). It then adds interior
-   * support points that are either computed by interpolation from the
-   * surrounding points using weights computed by solving a Laplace
-   * equation, or if dim<spacedim, it asks the underlying manifold for
-   * the locations of interior points.
-   */
-  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 of the given cell to the
-   * vector @p a. Points located on the vertices of a line are not
-   * included.
-   *
-   * Needed by the @p compute_support_points() function. For
-   * <tt>dim=1</tt> this function is empty. The function uses the
-   * underlying manifold object of the line (or, if none is set, of
-   * the cell) for the location of the requested points.
-   *
-   * 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) of the given
-   * cell to the vector @p a. Points located on the vertices or lines
-   * of a quad are not included.
-   *
-   * Needed by the @p compute_support_points() function. For
-   * <tt>dim=1</tt> and <tt>dim=2</tt> this function is empty. The
-   * function uses the underlying manifold object of the quad (or, if
-   * none is set, of the cell) for the location of the requested
-   * points.
-   *
-   * 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;
-
-  /**
-   * Needed by the @p support_point_weights_on_quad function (for <tt>dim==2</tt>). Filled
-   * by the constructor.
-   *
-   * 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> support_point_weights_on_quad;
-
-  /**
-   * 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> support_point_weights_on_hex;
-
   /**
    * Exception.
    */
@@ -377,15 +286,6 @@ protected:
    */
   const bool use_mapping_q_on_all_cells;
 
-  /**
-   * An FE_Q object which is only needed in 3D, since it knows how to reorder
-   * shape functions/DoFs on non-standard faces. This is used to reorder
-   * support points in the same way. We could make this a pointer to prevent
-   * construction in 1D and 2D, but since memory and time requirements are not
-   * particularly high this seems unnecessary at the moment.
-   */
-  const FE_Q<dim> feq;
-
   /**
    * Pointer to a Q1 mapping. This mapping is used on interior cells unless
    * use_mapping_q_on_all_cells was set in the call to the
@@ -403,16 +303,11 @@ protected:
    */
   std_cxx11::unique_ptr<const MappingQ1<dim,spacedim> > q1_mapping;
 
-  /*
-   * The default line support points. These are used when computing
-   * the location in real space of the support points on lines and
-   * quads, which are asked to the Manifold<dim,spacedim> class.
-   *
-   * The number of quadrature points depends on the degree of this
-   * class, and it matches the number of degrees of freedom of an
-   * FE_Q<1>(this->degree).
-   */
-  QGaussLobatto<1> line_support_points;
+  //TODO: Remove again -- all the function does is bypass the inherited function from MappingQ1 and go back to the one in the MappingQGeneric base class
+  virtual
+  void
+  compute_mapping_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                                  std::vector<Point<spacedim> > &a) const;
 
   /**
    * Declare other MappingQ classes friends.
index f7b4ae5afdc726a0c9fdca6bcfc65d80d435fb62..46fa2b4b85d4a109b929298fa078046da01d67ac 100644 (file)
 #include <deal.II/base/derivative_form.h>
 #include <deal.II/base/config.h>
 #include <deal.II/base/table.h>
+#include <deal.II/base/quadrature_lib.h>
 #include <deal.II/base/qprojector.h>
 #include <deal.II/grid/tria_iterator.h>
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/fe/mapping.h>
+#include <deal.II/fe/fe_q.h>
 
 #include <cmath>
 
@@ -66,6 +68,11 @@ public:
    */
   MappingQGeneric (const unsigned int polynomial_degree);
 
+  /**
+   * Copy constructor.
+   */
+  MappingQGeneric (const MappingQGeneric<dim,spacedim> &mapping);
+
   /**
    * Return the degree of the mapping, i.e. the value which was passed to the
    * constructor.
@@ -453,6 +460,48 @@ protected:
    */
   const unsigned int polynomial_degree;
 
+  /*
+   * The default line support points. These are used when computing
+   * the location in real space of the support points on lines and
+   * quads, which are asked to the Manifold<dim,spacedim> class.
+   *
+   * The number of quadrature points depends on the degree of this
+   * class, and it matches the number of degrees of freedom of an
+   * FE_Q<1>(this->degree).
+   */
+  QGaussLobatto<1> line_support_points;
+
+  /**
+   * An FE_Q object which is only needed in 3D, since it knows how to reorder
+   * shape functions/DoFs on non-standard faces. This is used to reorder
+   * support points in the same way.
+   */
+  const std_cxx11::unique_ptr<FE_Q<dim> > fe_q;
+
+  /**
+   * A table of weights by which we multiply the locations of the
+   * support points on the perimeter of a quad to get the location of
+   * interior 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> support_point_weights_on_quad;
+
+  /**
+   * A table of weights by which we multiply the locations of the
+   * support points on the perimeter of a hex to get the location of
+   * interior support points.
+   *
+   * For the definition of this vector see equation (8) of the `mapping'
+   * report.
+   */
+  Table<2,double> support_point_weights_on_hex;
+
   /**
    * An interface that derived classes have to implement and that
    * computes the locations of support points for the mapping. For
@@ -462,11 +511,72 @@ protected:
    * the support points from the geometry of the current cell but
    * instead evaluating an externally given displacement field in
    * addition to the geometry of the cell.
+   *
+   * The default implementation of this function is appropriate for
+   * most cases. It takes the locations of support points on the
+   * boundary of the cell from the underlying manifold. Interior
+   * support points (ie. support points in quads for 2d, in hexes for
+   * 3d) are then 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.
+   *
+   * The function works its way from the vertices (which it takes from
+   * the given cell) via the support points on the line (for which it
+   * calls the add_line_support_points() function) and the support
+   * points on the quad faces (in 3d, for which it calls the
+   * add_quad_support_points() function). It then adds interior
+   * support points that are either computed by interpolation from the
+   * surrounding points using weights computed by solving a Laplace
+   * equation, or if dim<spacedim, it asks the underlying manifold for
+   * the locations of interior points.
    */
   virtual
   void
   compute_mapping_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
-                                  std::vector<Point<spacedim> > &a) const = 0;
+                                  std::vector<Point<spacedim> > &a) const;
+
+  /**
+   * For <tt>dim=2,3</tt>. Append the support points of all shape
+   * functions located on bounding lines of the given cell to the
+   * vector @p a. Points located on the vertices of a line are not
+   * included.
+   *
+   * Needed by the @p compute_support_points() function. For
+   * <tt>dim=1</tt> this function is empty. The function uses the
+   * underlying manifold object of the line (or, if none is set, of
+   * the cell) for the location of the requested points.
+   *
+   * 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) of the given
+   * cell to the vector @p a. Points located on the vertices or lines
+   * of a quad are not included.
+   *
+   * Needed by the @p compute_support_points() function. For
+   * <tt>dim=1</tt> and <tt>dim=2</tt> this function is empty. The
+   * function uses the underlying manifold object of the quad (or, if
+   * none is set, of the cell) for the location of the requested
+   * points.
+   *
+   * 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;
 
   /**
    * Make MappingQ a friend since it needs to call the
index 62b02f8197459bc3cbcf9b99207d06c32ce5c219..6ac2a071d5e146ea8590d52d6b148b2c3e9aec9d 100644 (file)
@@ -22,7 +22,6 @@
 #include <deal.II/base/std_cxx11/unique_ptr.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_boundary.h>
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/fe/mapping_q.h>
 #include <deal.II/fe/fe_q.h>
@@ -53,222 +52,6 @@ MappingQ<dim,spacedim>::InternalData::memory_consumption () const
 
 
 
-namespace
-{
-  /**
-   * Compute the <tt>support_point_weights_on_quad(hex)</tt> arrays.
-   *
-   * Called by the <tt>compute_support_point_weights_on_quad(hex)</tt> functions if the
-   * data is not yet hardcoded.
-   *
-   * For the definition of the <tt>support_point_weights_on_quad(hex)</tt> please
-   * refer to equation (8) of the `mapping' report.
-   */
-  template<int dim>
-  Table<2,double>
-  compute_laplace_vector(const unsigned int polynomial_degree)
-  {
-    Table<2,double> lvs;
-
-    Assert(lvs.n_rows()==0, ExcInternalError());
-    Assert(dim==2 || dim==3, ExcNotImplemented());
-
-    // for degree==1, we shouldn't have to compute any support points, since all
-    // of them are on the vertices
-    Assert(polynomial_degree>1, ExcInternalError());
-
-    const unsigned int n_inner = Utilities::fixed_power<dim>(polynomial_degree-1);
-    const unsigned int n_outer = (dim==1) ? 2 :
-                                 ((dim==2) ?
-                                  4+4*(polynomial_degree-1) :
-                                  8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1));
-
-
-    // compute the shape gradients at the quadrature points on the unit cell
-    const QGauss<dim> quadrature(polynomial_degree+1);
-    const unsigned int n_q_points=quadrature.size();
-
-    typename MappingQGeneric<dim>::InternalData quadrature_data(polynomial_degree);
-    quadrature_data.shape_derivatives.resize(quadrature_data.n_shape_functions *
-                                             n_q_points);
-    quadrature_data.compute_shape_function_values(quadrature.get_points());
-
-    // Compute the stiffness matrix of the inner dofs
-    FullMatrix<long double> S(n_inner);
-    for (unsigned int point=0; point<n_q_points; ++point)
-      for (unsigned int i=0; i<n_inner; ++i)
-        for (unsigned int j=0; j<n_inner; ++j)
-          {
-            long double res = 0.;
-            for (unsigned int l=0; l<dim; ++l)
-              res += (long double)quadrature_data.derivative(point, n_outer+i)[l] *
-                     (long double)quadrature_data.derivative(point, n_outer+j)[l];
-
-            S(i,j) += res * (long double)quadrature.weight(point);
-          }
-
-    // Compute the components of T to be the product of gradients of inner and
-    // outer shape functions.
-    FullMatrix<long double> T(n_inner, n_outer);
-    for (unsigned int point=0; point<n_q_points; ++point)
-      for (unsigned int i=0; i<n_inner; ++i)
-        for (unsigned int k=0; k<n_outer; ++k)
-          {
-            long double res = 0.;
-            for (unsigned int l=0; l<dim; ++l)
-              res += (long double)quadrature_data.derivative(point, n_outer+i)[l] *
-                     (long double)quadrature_data.derivative(point, k)[l];
-
-            T(i,k) += res *(long double)quadrature.weight(point);
-          }
-
-    FullMatrix<long double> S_1(n_inner);
-    S_1.invert(S);
-
-    FullMatrix<long double> S_1_T(n_inner, n_outer);
-
-    // S:=S_1*T
-    S_1.mmult(S_1_T,T);
-
-    // Resize and initialize the lvs
-    lvs.reinit (n_inner, n_outer);
-    for (unsigned int i=0; i<n_inner; ++i)
-      for (unsigned int k=0; k<n_outer; ++k)
-        lvs(i,k) = -S_1_T(i,k);
-
-    return lvs;
-  }
-
-
-  /**
-   * 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 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 support_point_weights_on_quad please refer to
-   * equation (8) of the `mapping' report.
-   */
-  template<int dim>
-  Table<2,double>
-  compute_support_point_weights_on_quad(const unsigned int polynomial_degree)
-  {
-    Table<2,double> loqvs;
-
-    // in 1d, there are no quads, so return an empty object
-    if (dim == 1)
-      return loqvs;
-
-    // we are asked to compute weights for interior support points, but
-    // there are no interior points if degree==1
-    if (polynomial_degree == 1)
-      return loqvs;
-
-    const unsigned int n_inner_2d=(polynomial_degree-1)*(polynomial_degree-1);
-    const unsigned int n_outer_2d=4+4*(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
-    if (polynomial_degree == 2)
-      {
-        // (checked these values against the output of compute_laplace_vector
-        // again, and found they're indeed right -- just in case someone wonders
-        // where they come from -- WB)
-        static const double loqv2[1*8]
-          = {1/16., 1/16., 1/16., 1/16., 3/16., 3/16., 3/16., 3/16.};
-        Assert (sizeof(loqv2)/sizeof(loqv2[0]) ==
-                n_inner_2d * n_outer_2d,
-                ExcInternalError());
-
-        // copy and return
-        loqvs.reinit(n_inner_2d, n_outer_2d);
-        for (unsigned int unit_point=0; unit_point<n_inner_2d; ++unit_point)
-          for (unsigned int k=0; k<n_outer_2d; ++k)
-            loqvs[unit_point][k] = loqv2[unit_point*n_outer_2d+k];
-      }
-    else
-      {
-        // not precomputed, then do so now
-        loqvs = compute_laplace_vector<2>(polynomial_degree);
-      }
-
-    // the sum of weights of the points at the outer rim should be one. check
-    // 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*polynomial_degree,
-             ExcInternalError());
-
-    return loqvs;
-  }
-
-
-
-  /**
-   * This function is needed by the constructor of <tt>MappingQ<3></tt>.
-   *
-   * 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 support_point_weights_on_hex please refer to
-   * equation (8) of the `mapping' report.
-   */
-  template <int dim>
-  Table<2,double>
-  compute_support_point_weights_on_hex(const unsigned int polynomial_degree)
-  {
-    Table<2,double> lohvs;
-
-    // in 1d and 2d, there are no hexes, so return an empty object
-    if (dim < 3)
-      return lohvs;
-
-    // we are asked to compute weights for interior support points, but
-    // there are no interior points if degree==1
-    if (polynomial_degree == 1)
-      return lohvs;
-
-    const unsigned int n_inner = Utilities::fixed_power<dim>(polynomial_degree-1);
-    const unsigned int n_outer = 8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(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
-    if (polynomial_degree == 2)
-      {
-        static const double lohv2[26]
-          = {1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128.,
-             7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192.,
-             7/192., 7/192., 7/192., 7/192.,
-             1/12., 1/12., 1/12., 1/12., 1/12., 1/12.
-            };
-
-        // copy and return
-        lohvs.reinit(n_inner, n_outer);
-        for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
-          for (unsigned int k=0; k<n_outer; ++k)
-            lohvs[unit_point][k] = lohv2[unit_point*n_outer+k];
-      }
-    else
-      {
-        // not precomputed, then do so now
-        lohvs = compute_laplace_vector<dim>(polynomial_degree);
-      }
-
-    // the sum of weights of the points at the outer rim should be one. check
-    // 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*polynomial_degree,
-             ExcInternalError());
-
-    return lohvs;
-  }
-}
-
-
-
 template<int dim, int spacedim>
 MappingQ<dim,spacedim>::MappingQ (const unsigned int degree,
                                   const bool use_mapping_q_on_all_cells)
@@ -290,17 +73,12 @@ MappingQ<dim,spacedim>::MappingQ (const unsigned int degree,
                               use_mapping_q_on_all_cells
                               ||
                               (dim != spacedim)),
-  feq(degree),
   // create a Q1 mapping for use on interior cells (if necessary)
   // or to create a good initial guess in transform_real_to_unit_cell()
-  q1_mapping (new MappingQ1<dim,spacedim>()),
-  line_support_points(degree+1)
+  q1_mapping (new MappingQ1<dim,spacedim>())
 {
   Assert(n_inner+n_outer==Utilities::fixed_power<dim>(degree+1),
          ExcInternalError());
-
-  support_point_weights_on_quad = compute_support_point_weights_on_quad<dim>(this->polynomial_degree);
-  support_point_weights_on_hex = compute_support_point_weights_on_hex<dim>(this->polynomial_degree);
 }
 
 
@@ -312,17 +90,12 @@ MappingQ<dim,spacedim>::MappingQ (const MappingQ<dim,spacedim> &mapping)
   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()),
   // clone the Q1 mapping for use on interior cells (if necessary)
   // or to create a good initial guess in transform_real_to_unit_cell()
-  q1_mapping (dynamic_cast<MappingQ1<dim,spacedim>*>(mapping.q1_mapping->clone())),
-  line_support_points(mapping.line_support_points)
+  q1_mapping (dynamic_cast<MappingQ1<dim,spacedim>*>(mapping.q1_mapping->clone()))
 {
   Assert(n_inner+n_outer==Utilities::fixed_power<dim>(this->polynomial_degree+1),
          ExcInternalError());
-
-  support_point_weights_on_quad = compute_support_point_weights_on_quad<dim>(this->polynomial_degree);
-  support_point_weights_on_hex = compute_support_point_weights_on_hex<dim>(this->polynomial_degree);
 }
 
 
@@ -542,452 +315,13 @@ fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterato
 }
 
 
-namespace
-{
-  /**
-   * 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 2:
-      {
-        // 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 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.
-        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;
-      }
-  }
-
-
-
-
-  /**
-   * 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>
-void
-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)
-    {
-      for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
-        {
-          const typename Triangulation<dim,spacedim>::line_iterator line =
-            (dim == 1  ?
-             static_cast<typename Triangulation<dim,spacedim>::line_iterator>(cell) :
-             cell->line(line_no));
-
-          const Manifold<dim,spacedim> &manifold =
-            ( ( line->manifold_id() == numbers::invalid_manifold_id ) &&
-              ( 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
-    // from the boundary description
-    {
-      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)
-        {
-          const typename Triangulation<dim,spacedim>::line_iterator
-          line = (dim == 1
-                  ?
-                  static_cast<typename Triangulation<dim,spacedim>::line_iterator>(cell)
-                  :
-                  cell->line(line_no));
-
-          const Manifold<dim,spacedim> &manifold =
-            ( ( line->manifold_id() == numbers::invalid_manifold_id ) &&
-              ( dim < spacedim )
-              ?
-              cell->get_manifold() :
-              line->get_manifold() );
-
-          get_intermediate_points_on_object (manifold, line_support_points, line, line_points);
-
-          if (dim==3)
-            {
-              // in 3D, lines might be in wrong orientation. if so, reverse
-              // the vector
-              if (cell->line_orientation(line_no))
-                a.insert (a.end(), line_points.begin(), line_points.end());
-              else
-                a.insert (a.end(), line_points.rbegin(), line_points.rend());
-            }
-          else
-            // in 2D, lines always have the correct orientation. simply append
-            // all points
-            a.insert (a.end(), line_points.begin(), line_points.end());
-        }
-    }
-}
-
-
-
-template <>
-void
-MappingQ<3,3>::
-add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
-                        std::vector<Point<3> >                &a) const
-{
-  const unsigned int faces_per_cell    = GeometryInfo<3>::faces_per_cell,
-                     vertices_per_face = GeometryInfo<3>::vertices_per_face,
-                     lines_per_face    = GeometryInfo<3>::lines_per_face,
-                     vertices_per_cell = GeometryInfo<3>::vertices_per_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 ((polynomial_degree-1)*(polynomial_degree-1));
-  // used if only one line of face quad is at boundary
-  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.
-  std::vector<Point<3> > vertices(4);
-
-  // loop over all faces and collect points on them
-  for (unsigned int face_no=0; face_no<faces_per_cell; ++face_no)
-    {
-      const Triangulation<3>::face_iterator face = cell->face(face_no);
-
-      // select the correct mappings for the present face
-      const bool face_orientation = cell->face_orientation(face_no),
-                 face_flip        = cell->face_flip       (face_no),
-                 face_rotation    = cell->face_rotation   (face_no);
-
-#ifdef DEBUG
-      // some sanity checks up front
-      for (unsigned int i=0; i<vertices_per_face; ++i)
-        Assert(face->vertex_index(i)==cell->vertex_index(
-                 GeometryInfo<3>::face_to_cell_vertices(face_no, i,
-                                                        face_orientation,
-                                                        face_flip,
-                                                        face_rotation)),
-               ExcInternalError());
-
-      // indices of the lines that bound a face are given by GeometryInfo<3>::
-      // face_to_cell_lines
-      for (unsigned int i=0; i<lines_per_face; ++i)
-        Assert(face->line(i)==cell->line(GeometryInfo<3>::face_to_cell_lines(
-                                           face_no, i, face_orientation, face_flip, face_rotation)),
-               ExcInternalError());
-#endif
-
-      // if face at boundary, then ask boundary object to return intermediate
-      // points on it
-      if (face->at_boundary())
-        {
-          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 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,
-                                    face_flip,
-                                    face_rotation)]);
-        }
-      else
-        {
-          // face is not at boundary, but maybe some of its lines are. count
-          // them
-          unsigned int lines_at_boundary=0;
-          for (unsigned int i=0; i<lines_per_face; ++i)
-            if (face->line(i)->at_boundary())
-              ++lines_at_boundary;
-
-          Assert(lines_at_boundary<=lines_per_face, ExcInternalError());
-
-          // if at least one of the lines bounding this quad is at the
-          // boundary, then collect points separately
-          if (lines_at_boundary>0)
-            {
-              // 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*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*(polynomial_degree-1),
-                     ExcDimensionMismatch(b.size(),
-                                          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
-              // orientation. Doing so, the calculated points will be in
-              // standard orientation as well.
-              for (unsigned int i=0; i<vertices_per_face; ++i)
-                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<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)*(polynomial_degree-1)+j];
-
-              // Now b includes the support points on the quad and we can
-              // apply the laplace vector
-              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<(polynomial_degree-1)*(polynomial_degree-1); ++i)
-                a.push_back(b[4*polynomial_degree+i]);
-            }
-          else
-            {
-              // face is entirely in the interior. get intermediate
-              // points from the relevant manifold object.
-              vertices.resize(4);
-              for (unsigned int i=0; i<4; ++i)
-                vertices[i] = face->vertex(i);
-              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
-              // uses the same shape function as an FEQ, we can ask a FEQ 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,
-                                        face_flip,
-                                        face_rotation)]);
-            }
-        }
-    }
-}
-
-
-
-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 ((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>
-void
-MappingQ<dim,spacedim>::
-add_quad_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &,
-                        std::vector<Point<spacedim> > &) const
-{
-  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;
-      }
+  MappingQGeneric<dim,spacedim>::compute_mapping_support_points (cell, a);
 }
 
 
index 155e6a23a938eae8d052cb833d2af9b8e22c83aa..adaac27fe753d85c48a5896f4e5eec994c6729c4 100644 (file)
@@ -25,6 +25,7 @@
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary.h>
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/fe/fe_tools.h>
 #include <deal.II/fe/fe.h>
@@ -767,11 +768,243 @@ compute_shape_function_values (const std::vector<Point<dim> > &unit_points)
 }
 
 
+namespace
+{
+  /**
+   * Compute the <tt>support_point_weights_on_quad(hex)</tt> arrays.
+   *
+   * Called by the <tt>compute_support_point_weights_on_quad(hex)</tt> functions if the
+   * data is not yet hardcoded.
+   *
+   * For the definition of the <tt>support_point_weights_on_quad(hex)</tt> please
+   * refer to equation (8) of the `mapping' report.
+   */
+  template<int dim>
+  Table<2,double>
+  compute_laplace_vector(const unsigned int polynomial_degree)
+  {
+    Table<2,double> lvs;
+
+    Assert(lvs.n_rows()==0, ExcInternalError());
+    Assert(dim==2 || dim==3, ExcNotImplemented());
+
+    // for degree==1, we shouldn't have to compute any support points, since all
+    // of them are on the vertices
+    Assert(polynomial_degree>1, ExcInternalError());
+
+    const unsigned int n_inner = Utilities::fixed_power<dim>(polynomial_degree-1);
+    const unsigned int n_outer = (dim==1) ? 2 :
+                                 ((dim==2) ?
+                                  4+4*(polynomial_degree-1) :
+                                  8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1));
+
+
+    // compute the shape gradients at the quadrature points on the unit cell
+    const QGauss<dim> quadrature(polynomial_degree+1);
+    const unsigned int n_q_points=quadrature.size();
+
+    typename MappingQGeneric<dim>::InternalData quadrature_data(polynomial_degree);
+    quadrature_data.shape_derivatives.resize(quadrature_data.n_shape_functions *
+                                             n_q_points);
+    quadrature_data.compute_shape_function_values(quadrature.get_points());
+
+    // Compute the stiffness matrix of the inner dofs
+    FullMatrix<long double> S(n_inner);
+    for (unsigned int point=0; point<n_q_points; ++point)
+      for (unsigned int i=0; i<n_inner; ++i)
+        for (unsigned int j=0; j<n_inner; ++j)
+          {
+            long double res = 0.;
+            for (unsigned int l=0; l<dim; ++l)
+              res += (long double)quadrature_data.derivative(point, n_outer+i)[l] *
+                     (long double)quadrature_data.derivative(point, n_outer+j)[l];
+
+            S(i,j) += res * (long double)quadrature.weight(point);
+          }
+
+    // Compute the components of T to be the product of gradients of inner and
+    // outer shape functions.
+    FullMatrix<long double> T(n_inner, n_outer);
+    for (unsigned int point=0; point<n_q_points; ++point)
+      for (unsigned int i=0; i<n_inner; ++i)
+        for (unsigned int k=0; k<n_outer; ++k)
+          {
+            long double res = 0.;
+            for (unsigned int l=0; l<dim; ++l)
+              res += (long double)quadrature_data.derivative(point, n_outer+i)[l] *
+                     (long double)quadrature_data.derivative(point, k)[l];
+
+            T(i,k) += res *(long double)quadrature.weight(point);
+          }
+
+    FullMatrix<long double> S_1(n_inner);
+    S_1.invert(S);
+
+    FullMatrix<long double> S_1_T(n_inner, n_outer);
+
+    // S:=S_1*T
+    S_1.mmult(S_1_T,T);
+
+    // Resize and initialize the lvs
+    lvs.reinit (n_inner, n_outer);
+    for (unsigned int i=0; i<n_inner; ++i)
+      for (unsigned int k=0; k<n_outer; ++k)
+        lvs(i,k) = -S_1_T(i,k);
+
+    return lvs;
+  }
+
+
+  /**
+   * 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 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 support_point_weights_on_quad please refer to
+   * equation (8) of the `mapping' report.
+   */
+  template<int dim>
+  Table<2,double>
+  compute_support_point_weights_on_quad(const unsigned int polynomial_degree)
+  {
+    Table<2,double> loqvs;
+
+    // in 1d, there are no quads, so return an empty object
+    if (dim == 1)
+      return loqvs;
+
+    // we are asked to compute weights for interior support points, but
+    // there are no interior points if degree==1
+    if (polynomial_degree == 1)
+      return loqvs;
+
+    const unsigned int n_inner_2d=(polynomial_degree-1)*(polynomial_degree-1);
+    const unsigned int n_outer_2d=4+4*(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
+    if (polynomial_degree == 2)
+      {
+        // (checked these values against the output of compute_laplace_vector
+        // again, and found they're indeed right -- just in case someone wonders
+        // where they come from -- WB)
+        static const double loqv2[1*8]
+          = {1/16., 1/16., 1/16., 1/16., 3/16., 3/16., 3/16., 3/16.};
+        Assert (sizeof(loqv2)/sizeof(loqv2[0]) ==
+                n_inner_2d * n_outer_2d,
+                ExcInternalError());
+
+        // copy and return
+        loqvs.reinit(n_inner_2d, n_outer_2d);
+        for (unsigned int unit_point=0; unit_point<n_inner_2d; ++unit_point)
+          for (unsigned int k=0; k<n_outer_2d; ++k)
+            loqvs[unit_point][k] = loqv2[unit_point*n_outer_2d+k];
+      }
+    else
+      {
+        // not precomputed, then do so now
+        loqvs = compute_laplace_vector<2>(polynomial_degree);
+      }
+
+    // the sum of weights of the points at the outer rim should be one. check
+    // 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*polynomial_degree,
+             ExcInternalError());
+
+    return loqvs;
+  }
+
+
+
+  /**
+   * This function is needed by the constructor of <tt>MappingQ<3></tt>.
+   *
+   * 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 support_point_weights_on_hex please refer to
+   * equation (8) of the `mapping' report.
+   */
+  template <int dim>
+  Table<2,double>
+  compute_support_point_weights_on_hex(const unsigned int polynomial_degree)
+  {
+    Table<2,double> lohvs;
+
+    // in 1d and 2d, there are no hexes, so return an empty object
+    if (dim < 3)
+      return lohvs;
+
+    // we are asked to compute weights for interior support points, but
+    // there are no interior points if degree==1
+    if (polynomial_degree == 1)
+      return lohvs;
+
+    const unsigned int n_inner = Utilities::fixed_power<dim>(polynomial_degree-1);
+    const unsigned int n_outer = 8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(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
+    if (polynomial_degree == 2)
+      {
+        static const double lohv2[26]
+          = {1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128.,
+             7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192.,
+             7/192., 7/192., 7/192., 7/192.,
+             1/12., 1/12., 1/12., 1/12., 1/12., 1/12.
+            };
+
+        // copy and return
+        lohvs.reinit(n_inner, n_outer);
+        for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
+          for (unsigned int k=0; k<n_outer; ++k)
+            lohvs[unit_point][k] = lohv2[unit_point*n_outer+k];
+      }
+    else
+      {
+        // not precomputed, then do so now
+        lohvs = compute_laplace_vector<dim>(polynomial_degree);
+      }
+
+    // the sum of weights of the points at the outer rim should be one. check
+    // 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*polynomial_degree,
+             ExcInternalError());
+
+    return lohvs;
+  }
+}
+
+
+
 
 template<int dim, int spacedim>
 MappingQGeneric<dim,spacedim>::MappingQGeneric (const unsigned int p)
   :
-  polynomial_degree(p)
+  polynomial_degree(p),
+  line_support_points(this->polynomial_degree+1),
+  fe_q(dim == 3 ? new FE_Q<dim>(this->polynomial_degree) : 0),
+  support_point_weights_on_quad (compute_support_point_weights_on_quad<dim>(this->polynomial_degree)),
+  support_point_weights_on_hex (compute_support_point_weights_on_hex<dim>(this->polynomial_degree))
+{}
+
+
+
+template<int dim, int spacedim>
+MappingQGeneric<dim,spacedim>::MappingQGeneric (const MappingQGeneric<dim,spacedim> &mapping)
+  :
+  polynomial_degree(mapping.polynomial_degree),
+  line_support_points(mapping.line_support_points),
+  fe_q(dim == 3 ? new FE_Q<dim>(*mapping.fe_q) : 0),
+  support_point_weights_on_quad (mapping.support_point_weights_on_quad),
+  support_point_weights_on_hex (mapping.support_point_weights_on_hex)
 {}
 
 
@@ -2281,7 +2514,6 @@ transform (const VectorSlice<const std::vector< Tensor<3,dim> > >   input,
            const typename Mapping<dim,spacedim>::InternalDataBase  &mapping_data,
            VectorSlice<std::vector<Tensor<3,spacedim> > >           output) const
 {
-
   switch (mapping_type)
     {
     case mapping_piola_hessian:
@@ -2292,7 +2524,454 @@ transform (const VectorSlice<const std::vector< Tensor<3,dim> > >   input,
     default:
       Assert(false, ExcNotImplemented());
     }
+}
+
+
+
+namespace
+{
+  /**
+   * 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 2:
+      {
+        // 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 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.
+        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;
+      }
+  }
+
+
+
+
+  /**
+   * 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>
+void
+MappingQGeneric<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)
+    {
+      for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
+        {
+          const typename Triangulation<dim,spacedim>::line_iterator line =
+            (dim == 1  ?
+             static_cast<typename Triangulation<dim,spacedim>::line_iterator>(cell) :
+             cell->line(line_no));
+
+          const Manifold<dim,spacedim> &manifold =
+            ( ( line->manifold_id() == numbers::invalid_manifold_id ) &&
+              ( 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
+    // from the boundary description
+    {
+      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)
+        {
+          const typename Triangulation<dim,spacedim>::line_iterator
+          line = (dim == 1
+                  ?
+                  static_cast<typename Triangulation<dim,spacedim>::line_iterator>(cell)
+                  :
+                  cell->line(line_no));
+
+          const Manifold<dim,spacedim> &manifold =
+            ( ( line->manifold_id() == numbers::invalid_manifold_id ) &&
+              ( dim < spacedim )
+              ?
+              cell->get_manifold() :
+              line->get_manifold() );
+
+          get_intermediate_points_on_object (manifold, line_support_points, line, line_points);
+
+          if (dim==3)
+            {
+              // in 3D, lines might be in wrong orientation. if so, reverse
+              // the vector
+              if (cell->line_orientation(line_no))
+                a.insert (a.end(), line_points.begin(), line_points.end());
+              else
+                a.insert (a.end(), line_points.rbegin(), line_points.rend());
+            }
+          else
+            // in 2D, lines always have the correct orientation. simply append
+            // all points
+            a.insert (a.end(), line_points.begin(), line_points.end());
+        }
+    }
+}
+
+
+
+template <>
+void
+MappingQGeneric<3,3>::
+add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
+                        std::vector<Point<3> >                &a) const
+{
+  const unsigned int faces_per_cell    = GeometryInfo<3>::faces_per_cell,
+                     vertices_per_face = GeometryInfo<3>::vertices_per_face,
+                     lines_per_face    = GeometryInfo<3>::lines_per_face,
+                     vertices_per_cell = GeometryInfo<3>::vertices_per_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 ((polynomial_degree-1)*(polynomial_degree-1));
+  // used if only one line of face quad is at boundary
+  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.
+  std::vector<Point<3> > vertices(4);
+
+  // loop over all faces and collect points on them
+  for (unsigned int face_no=0; face_no<faces_per_cell; ++face_no)
+    {
+      const Triangulation<3>::face_iterator face = cell->face(face_no);
+
+      // select the correct mappings for the present face
+      const bool face_orientation = cell->face_orientation(face_no),
+                 face_flip        = cell->face_flip       (face_no),
+                 face_rotation    = cell->face_rotation   (face_no);
+
+#ifdef DEBUG
+      // some sanity checks up front
+      for (unsigned int i=0; i<vertices_per_face; ++i)
+        Assert(face->vertex_index(i)==cell->vertex_index(
+                 GeometryInfo<3>::face_to_cell_vertices(face_no, i,
+                                                        face_orientation,
+                                                        face_flip,
+                                                        face_rotation)),
+               ExcInternalError());
+
+      // indices of the lines that bound a face are given by GeometryInfo<3>::
+      // face_to_cell_lines
+      for (unsigned int i=0; i<lines_per_face; ++i)
+        Assert(face->line(i)==cell->line(GeometryInfo<3>::face_to_cell_lines(
+                                           face_no, i, face_orientation, face_flip, face_rotation)),
+               ExcInternalError());
+#endif
+
+      // if face at boundary, then ask boundary object to return intermediate
+      // points on it
+      if (face->at_boundary())
+        {
+          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 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[fe_q->adjust_quad_dof_index_for_face_orientation(i,
+                                    face_orientation,
+                                    face_flip,
+                                    face_rotation)]);
+        }
+      else
+        {
+          // face is not at boundary, but maybe some of its lines are. count
+          // them
+          unsigned int lines_at_boundary=0;
+          for (unsigned int i=0; i<lines_per_face; ++i)
+            if (face->line(i)->at_boundary())
+              ++lines_at_boundary;
+
+          Assert(lines_at_boundary<=lines_per_face, ExcInternalError());
+
+          // if at least one of the lines bounding this quad is at the
+          // boundary, then collect points separately
+          if (lines_at_boundary>0)
+            {
+              // 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*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*(polynomial_degree-1),
+                     ExcDimensionMismatch(b.size(),
+                                          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
+              // orientation. Doing so, the calculated points will be in
+              // standard orientation as well.
+              for (unsigned int i=0; i<vertices_per_face; ++i)
+                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<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)*(polynomial_degree-1)+j];
+
+              // Now b includes the support points on the quad and we can
+              // apply the laplace vector
+              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<(polynomial_degree-1)*(polynomial_degree-1); ++i)
+                a.push_back(b[4*polynomial_degree+i]);
+            }
+          else
+            {
+              // face is entirely in the interior. get intermediate
+              // points from the relevant manifold object.
+              vertices.resize(4);
+              for (unsigned int i=0; i<4; ++i)
+                vertices[i] = face->vertex(i);
+              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
+              // uses the same shape 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[fe_q->adjust_quad_dof_index_for_face_orientation(i,
+                                        face_orientation,
+                                        face_flip,
+                                        face_rotation)]);
+            }
+        }
+    }
+}
+
+
+
+template <>
+void
+MappingQGeneric<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 ((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>
+void
+MappingQGeneric<dim,spacedim>::
+add_quad_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &,
+                        std::vector<Point<spacedim> > &) const
+{
+  Assert (false, ExcInternalError());
+}
+
+
+
+template<int dim, int spacedim>
+void
+MappingQGeneric<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;
+      }
 }
 
 

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.