]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Change MappingQ support points from equispaced points into Gauss-Lobatto points....
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 23 Aug 2013 15:48:29 +0000 (15:48 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 23 Aug 2013 15:48:29 +0000 (15:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@30450 0785d39b-7218-0410-832d-ea1e28bc413d

15 files changed:
deal.II/include/deal.II/fe/mapping_c1.h
deal.II/include/deal.II/fe/mapping_q.h
deal.II/include/deal.II/fe/mapping_q_eulerian.h
deal.II/include/deal.II/grid/tria.h
deal.II/include/deal.II/grid/tria_boundary.h
deal.II/source/fe/mapping_c1.cc
deal.II/source/fe/mapping_q.cc
deal.II/source/fe/mapping_q_eulerian.cc
deal.II/source/grid/tria_boundary.cc
deal.II/source/grid/tria_boundary_lib.cc
tests/fe/jacobians/cmp/generic
tests/fe/mapping/cmp/generic
tests/fe/mapping_q_eulerian/cmp/generic
tests/fe/mapping_real_to_unit_q4_sphere.cc
tests/fe/mapping_real_to_unit_q4_sphere/cmp/generic

index 85f900bf28023d8804255ec866eca9c5031cf9e2..62bbafecdaf5b49f89e7b5b5ffb5b2daba663ed4 100644 (file)
@@ -37,7 +37,7 @@ DEAL_II_NAMESPACE_OPEN
  *
  * To use this class, make sure that the
  * Boundary::@p get_normals_at_vertices function is implemented
- * for the users boundary object.
+ * for the user's boundary object.
  *
  * For more information about the <tt>spacedim</tt> template parameter
  * check the documentation of FiniteElement or the one of
@@ -50,71 +50,46 @@ class MappingC1 : public MappingQ<dim,spacedim>
 {
 public:
   /**
-   * Constructor. Pass the fixed
-   * degree @p 3 down to the base
-   * class, as a cubic mapping
-   * suffices to generate a
-   * continuous mapping of the
-   * boundary.
+   * Constructor. Pass the fixed degree @p 3 down to the base class, as a
+   * cubic mapping suffices to generate a continuous mapping of the boundary.
    */
   MappingC1 ();
 
   /**
-   * Return a pointer to a copy of the
-   * present object. The caller of this
-   * copy then assumes ownership of it.
+   * Return a pointer to a copy of the present object. The caller of this copy
+   * then assumes ownership of it.
    */
   virtual
   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 on vertices are not
-   * included.
+   * 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 on vertices are not included.
    *
-   * Needed by the
-   * <tt>compute_support_points_simple(laplace)</tt>
-   * functions. For <tt>dim=1</tt> this
-   * function is empty.
+   * Needed by the <tt>compute_support_points_simple(laplace)</tt>
+   * functions. For <tt>dim=1</tt> this function is empty.
    *
-   * This function chooses the
-   * respective points not such
-   * that they are interpolating
-   * the boundary (as does the base
-   * class), but rather such that
-   * the resulting cubic mapping is
-   * a continuous one.
+   * This function chooses the respective points not such that they are
+   * interpolating the boundary (as does the base class), but rather such that
+   * the resulting cubic mapping is a continuous one.
    */
   virtual void
   add_line_support_points (const typename Triangulation<dim>::cell_iterator &cell,
                            std::vector<Point<dim> > &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 line but on vertices
-   * are not included.
+   * 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 line but on vertices are not included.
    *
-   * Needed by the
-   * @p compute_support_points_laplace
-   * function. For <tt>dim=1</tt> and 2
-   * this function is empty.
+   * Needed by the @p compute_support_points_laplace function. For
+   * <tt>dim=1</tt> and 2 this function is empty.
    *
-   * This function chooses the
-   * respective points not such
-   * that they are interpolating
-   * the boundary (as does the base
-   * class), but rather such that
-   * the resulting cubic mapping is
-   * a continuous one.
+   * This function chooses the respective points not such that they are
+   * interpolating the boundary (as does the base class), but rather such that
+   * the resulting cubic mapping is a continuous one.
    */
   virtual void
   add_quad_support_points(const typename Triangulation<dim>::cell_iterator &cell,
index 484963727b0c16cd8ee65da5836b5d8e7237d942..28e4aefcbf1e75d334bc88c11bcfc2198c4b87d2 100644 (file)
@@ -32,9 +32,11 @@ template <int dim, typename POLY> class TensorProductPolynomials;
 /*@{*/
 
 /**
- * Mapping class that uses Qp-mappings on boundary cells. The mapping
- * shape functions make use of tensor product polynomials with
- * equidistant (on the unit cell) support points.
+ * Mapping class that uses Qp-mappings on boundary cells. The mapping shape
+ * functions make use of tensor product polynomials with unit cell support
+ * points equal to the points of the Gauss-Lobatto quadrature formula. These
+ * points give a well-conditioned interpolation also for very high orders and
+ * are therefore preferred over equidistant support points.
  *
  * For more details about Qp-mappings, see the `mapping' report at
  * <tt>deal.II/doc/reports/mapping_q/index.html</tt> in the `Reports'
@@ -44,6 +46,10 @@ template <int dim, typename POLY> class TensorProductPolynomials;
  * check the documentation of FiniteElement or the one of
  * Triangulation.
  *
+ * @note Since the boundary description is closely tied to the unit cell
+ * support points, new boundary descriptions need to explicitly use the
+ * Gauss-Lobatto points.
+ *
  * @author Ralf Hartmann, 2000, 2001, 2005; Guido Kanschat 2000, 2001
  */
 template <int dim, int spacedim=dim>
@@ -51,41 +57,25 @@ class MappingQ : public MappingQ1<dim,spacedim>
 {
 public:
   /**
-   * Constructor.  @p p gives the
-   * degree of mapping polynomials
-   * on boundary cells.
+   * Constructor.  @p p gives the degree of mapping polynomials on boundary
+   * cells.
    *
-   * The second argument determines
-   * whether the higher order
-   * mapping should also be used on
-   * interior cells. If its value
-   * is <code>false</code> (the
-   * default), the a lower-order
-   * mapping is used in the
-   * interior. This is sufficient
-   * for most cases where higher
-   * order mappings are only used
-   * to better approximate the
-   * boundary. In that case, cells
-   * bounded by straight lines are
-   * acceptable in the
-   * interior. However, there are
-   * cases where one would also
-   * like to use a higher order
-   * mapping in the interior. The
-   * MappingQEulerian class is one
-   * such case.
+   * The second argument determines whether the higher order mapping should
+   * also be used on interior cells. If its value is <code>false</code> (the
+   * default), the a lower-order mapping is used in the interior. This is
+   * sufficient for most cases where higher order mappings are only used to
+   * better approximate the boundary. In that case, cells bounded by straight
+   * lines are acceptable in the interior. However, there are cases where one
+   * would also like to use a higher order mapping in the interior. The
+   * MappingQEulerian class is one such case.
    */
   MappingQ (const unsigned int p,
             const bool use_mapping_q_on_all_cells = false);
 
   /**
-   * Copy constructor. Performs a
-   * deep copy, i.e. duplicates
-   * what #tensor_pols points to
-   * instead of simply copying the
-   * #tensor_pols pointer as done
-   * by a default copy constructor.
+   * Copy constructor. Performs a deep copy, i.e. duplicates what #tensor_pols
+   * points to instead of simply copying the #tensor_pols pointer as done by a
+   * default copy constructor.
    */
   MappingQ (const MappingQ<dim,spacedim> &mapping);
 
@@ -95,10 +85,8 @@ public:
   virtual ~MappingQ ();
 
   /**
-   * Transforms the point @p p on
-   * the unit cell to the point
-   * @p p_real on the real cell
-   * @p cell and returns @p p_real.
+   * Transforms the point @p p on the unit cell to the point @p p_real on the
+   * real cell @p cell and returns @p p_real.
    */
   virtual Point<spacedim>
   transform_unit_to_real_cell (
@@ -106,46 +94,25 @@ public:
     const Point<dim>                                 &p) const;
 
   /**
-   * Transforms the point @p p on
-   * the real cell to the point
-   * @p p_unit on the unit cell
-   * @p cell and returns @p p_unit.
+   * Transforms the point @p p on the real cell to the point @p p_unit on the
+   * unit cell @p cell and returns @p p_unit.
    *
-   * Uses Newton iteration and the
-   * @p transform_unit_to_real_cell
-   * function.
+   * Uses Newton iteration and the @p transform_unit_to_real_cell function.
    *
-   * In the codimension one case,
-   * this function returns the
-   * normal projection of the real
-   * point @p p on the curve or
-   * surface identified by the @p
-   * cell.
+   * In the codimension one case, this function returns the normal projection
+   * of the real point @p p on the curve or surface identified by the @p cell.
    *
-   * @note Polynomial mappings from
-   * the reference (unit) cell coordinates
-   * to the coordinate system of a real
-   * cell are not always invertible if
-   * the point for which the inverse
-   * mapping is to be computed lies
-   * outside the cell's boundaries.
-   * In such cases, the current function
-   * may fail to compute a point on
-   * the reference cell whose image
-   * under the mapping equals the given
-   * point @p p.  If this is the case
-   * then this function throws an
-   * exception of type
-   * Mapping::ExcTransformationFailed .
-   * Whether the given point @p p lies
-   * outside the cell can therefore be
-   * determined by checking whether the
-   * return reference coordinates lie
-   * inside of outside the reference
-   * cell (e.g., using
-   * GeometryInfo::is_inside_unit_cell)
-   * or whether the exception mentioned
-   * above has been thrown.
+   * @note Polynomial mappings from the reference (unit) cell coordinates to
+   * the coordinate system of a real cell are not always invertible if the
+   * point for which the inverse mapping is to be computed lies outside the
+   * cell's boundaries.  In such cases, the current function may fail to
+   * compute a point on the reference cell whose image under the mapping
+   * equals the given point @p p.  If this is the case then this function
+   * throws an exception of type Mapping::ExcTransformationFailed .  Whether
+   * the given point @p p lies outside the cell can therefore be determined by
+   * checking whether the return reference coordinates lie inside of outside
+   * the reference cell (e.g., using GeometryInfo::is_inside_unit_cell) or
+   * whether the exception mentioned above has been thrown.
    */
   virtual Point<dim>
   transform_real_to_unit_cell (
@@ -172,23 +139,20 @@ public:
              const MappingType type) const;
 
   /**
-   * Return the degree of the
-   * mapping, i.e. the value which
-   * was passed to the constructor.
+   * Return the degree of the mapping, i.e. the value which was passed to the
+   * constructor.
    */
   unsigned int get_degree () const;
 
   /**
-   * Return a pointer to a copy of the
-   * present object. The caller of this
-   * copy then assumes ownership of it.
+   * Return a pointer to a copy of the present object. The caller of this copy
+   * then assumes ownership of it.
    */
   virtual
   Mapping<dim,spacedim> *clone () const;
 
   /**
-   * Storage for internal data of
-   * Q_degree transformation.
+   * Storage for internal data of Q_degree transformation.
    */
   class InternalData : public MappingQ1<dim,spacedim>::InternalData
   {
@@ -200,50 +164,35 @@ public:
 
 
     /**
-     * Return an estimate (in
-     * bytes) or the memory
-     * consumption of this
-     * object.
+     * Return an estimate (in bytes) or the memory consumption of this object.
      */
     virtual std::size_t memory_consumption () const;
 
     /**
-     * Unit normal vectors. Used
-     * for the alternative
-     * computation of the normal
-     * vectors. See doc of the
-     * @p alternative_normals_computation
-     * flag.
+     * Unit normal vectors. Used for the alternative computation of the normal
+     * vectors. See doc of the @p alternative_normals_computation flag.
      *
-     * Filled (hardcoded) once in
-     * @p get_face_data.
+     * Filled (hardcoded) once in @p get_face_data.
      */
     std::vector<std::vector<Point<dim> > > unit_normals;
 
     /**
-     * Flag that is set by the
-     * <tt>fill_fe_[[sub]face]_values</tt>
-     * function.
+     * Flag that is set by the <tt>fill_fe_[[sub]face]_values</tt> function.
      *
-     * If this flag is @p true
-     * we are on an interior cell
-     * and the
-     * @p mapping_q1_data is
-     * used.
+     * If this flag is @p true we are on an interior cell and the @p
+     * mapping_q1_data is used.
      */
     bool use_mapping_q1_on_current_cell;
 
     /**
-     * On interior cells
-     * @p MappingQ1 is used.
+     * On interior cells @p MappingQ1 is used.
      */
     typename MappingQ1<dim,spacedim>::InternalData mapping_q1_data;
   };
 
 protected:
   /**
-   * Implementation of the interface in
-   * Mapping.
+   * Implementation of the interface in Mapping.
    */
   virtual void
   fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
@@ -258,8 +207,7 @@ protected:
                   CellSimilarity::Similarity                           &cell_similarity) const ;
 
   /**
-   * Implementation of the interface in
-   * Mapping.
+   * Implementation of the interface in Mapping.
    */
   virtual void
   fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
@@ -272,8 +220,7 @@ protected:
                        typename std::vector<Point<spacedim> >        &normal_vectors) const ;
 
   /**
-   * Implementation of the interface in
-   * Mapping.
+   * Implementation of the interface in Mapping.
    */
   virtual void
   fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
@@ -287,55 +234,32 @@ protected:
                           typename std::vector<Point<spacedim> >        &normal_vectors) 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.
+   * 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.
+   * 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.
+   * 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.
+   * 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.
+   * 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.
+   * 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,
@@ -359,145 +283,99 @@ private:
                     const Quadrature<dim-1>& quadrature) const;
 
   /**
-   * Compute shape values and/or
-   * derivatives.
+   * Compute shape values and/or derivatives.
    */
   virtual void
   compute_shapes_virtual (const std::vector<Point<dim> > &unit_points,
                           typename MappingQ1<dim,spacedim>::InternalData &data) const;
 
   /**
-   * This function is needed by the
-   * constructor of <tt>MappingQ<dim,spacedim></tt>
-   * for <tt>dim=</tt> 2 and 3.
+   * 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
-   * the hardcoded data. For
-   * <tt>degree>=4</tt> and MappingQ<2>
-   * this vector is computed.
+   * For <tt>degree<4</tt> this function sets the @p laplace_on_quad_vector 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 equation (8)
-   * of the `mapping' report.
+   * For the definition of the @p laplace_on_quad_vector please refer to
+   * equation (8) of the `mapping' report.
    */
   void
   set_laplace_on_quad_vector(Table<2,double> &loqvs) const;
 
   /**
-   * This function is needed by the
-   * constructor of <tt>MappingQ<3></tt>.
+   * 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
-   * the hardcoded data. For
-   * <tt>degree>2</tt> this vector is
-   * computed.
+   * For <tt>degree==2</tt> this function sets the @p laplace_on_hex_vector 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 equation (8)
-   * of the `mapping' report.
+   * For the definition of the @p laplace_on_hex_vector please refer to
+   * equation (8) of the `mapping' report.
    */
   void set_laplace_on_hex_vector(Table<2,double> &lohvs) const;
 
   /**
-   * Computes the
-   * <tt>laplace_on_quad(hex)_vector</tt>.
+   * Computes the <tt>laplace_on_quad(hex)_vector</tt>.
    *
-   * Called by the
-   * <tt>set_laplace_on_quad(hex)_vector</tt>
-   * functions if the data is not
-   * yet hardcoded.
+   * Called by the <tt>set_laplace_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 refer to equation (8)
-   * of the `mapping' report.
+   * For the definition of the <tt>laplace_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.
+   * 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.
+   * 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.
+   * See equation (7) of the `mapping' report.
    */
   void apply_laplace_vector(const Table<2,double>   &lvs,
                             std::vector<Point<spacedim> > &a) const;
 
   /**
-   * Computes the support points of
-   * the mapping.
+   * 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;
 
   /**
-   * 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.
+   * 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.
    */
   void compute_support_points_laplace(
     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 by the
-   * constructor.
+   * Needed by the @p laplace_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, i.e.
-   *   unit_support_points on the
-   *   boundary of the quad
+   * laplace_on_quad_vector.size()= number of inner unit_support_points
+   * laplace_on_quad_vector[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.
+   * For the definition of this vector see equation (8) of the `mapping'
+   * report.
    */
   Table<2,double> laplace_on_quad_vector;
 
   /**
-   * Needed by the
-   * @p laplace_on_hex function
-   * (for <tt>dim==3</tt>). Filled by the
-   * constructor.
+   * Needed by the @p laplace_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.
+   * For the definition of this vector see equation (8) of the `mapping'
+   * report.
    */
   Table<2,double> laplace_on_hex_vector;
 
@@ -509,69 +387,57 @@ private:
                   << "laplace_vector not set for degree=" << arg1 << ".");
 
   /**
-   * Degree @p p of the
-   * polynomials used as shape
-   * functions for the Qp mapping
+   * Degree @p p of the polynomials used as shape functions for the Qp mapping
    * of cells at the boundary.
    */
   const unsigned int degree;
 
   /**
-   * Number of inner mapping shape
-   * functions.
+   * Number of inner mapping shape functions.
    */
   const unsigned int n_inner;
 
   /**
-   * Number of mapping shape
-   * functions on the boundary.
+   * Number of mapping shape functions on the boundary.
    */
   const unsigned int n_outer;
 
   /**
-   * Pointer to the
-   * @p dim-dimensional tensor
-   * product polynomials used as
-   * shape functions for the Qp
-   * mapping of cells at the
-   * boundary.
+   * Pointer to the @p dim-dimensional tensor product polynomials used as
+   * shape functions for the Qp mapping of cells at the boundary.
    */
   const TensorProductPolynomials<dim> *tensor_pols;
 
   /**
-   * Number of the Qp tensor
-   * product shape functions.
+   * Number of the Qp tensor product shape functions.
    */
   const unsigned int n_shape_functions;
 
   /**
-   * Mapping from lexicographic to
-   * to the Qp shape function
-   * numbering. Its size is
-   * @p dofs_per_cell.
+   * Mapping from lexicographic to to the Qp shape function numbering. Its
+   * size is @p dofs_per_cell.
    */
   const std::vector<unsigned int> renumber;
 
   /**
-   * If this flag is set @p true
-   * then @p MappingQ is used on
-   * all cells, not only on
-   * boundary cells.
+   * If this flag is set @p true then @p MappingQ is used on all cells, not
+   * only on boundary cells.
    */
   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.
+   * 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;
+
+  /**
+   * Declare other MappingQ classes friends.
+   */
+  template <int,int> friend class MappingQ;
 };
 
 /*@}*/
index 4d4ba9cc9cfdcc0ea58d0a83562922b964f56555..f0a79f1066f4c30b4026c06e082e0b6cdead19cf 100644 (file)
@@ -95,25 +95,14 @@ class MappingQEulerian : public MappingQ<dim, spacedim>
 {
 public:
   /**
-   * Constructor. The first argument is
-   * the polynomical degree of the desired
-   * Qp mapping.  It then takes a
-   * <tt>Vector<double> &</tt> to specify the
-   * transformation of the domain
-   * from the reference to
-   * the current configuration.
-   * The organization of the
-   * elements in the @p Vector
-   * must follow the concept how
-   * deal.II stores solutions that
-   * are associated to a
-   * triangulation.  This is
-   * automatically the case if the
-   * @p Vector represents the
-   * solution of the previous step
-   * of a nonlinear problem.
-   * Alternatively, the @p Vector
-   * can be initialized by
+   * Constructor. The first argument is the polynomical degree of the desired
+   * Qp mapping.  It then takes a <tt>Vector<double> &</tt> to specify the
+   * transformation of the domain from the reference to the current
+   * configuration.  The organization of the elements in the @p Vector must
+   * follow the concept how deal.II stores solutions that are associated to a
+   * triangulation.  This is automatically the case if the @p Vector
+   * represents the solution of the previous step of a nonlinear problem.
+   * Alternatively, the @p Vector can be initialized by
    * <tt>DoFAccessor::set_dof_values()</tt>.
    */
 
@@ -122,20 +111,16 @@ public:
                     const DoFHandler<dim,spacedim>  &euler_dof_handler);
 
   /**
-   * Return a pointer to a copy of the
-   * present object. The caller of this
-   * copy then assumes ownership of it.
+   * Return a pointer to a copy of the present object. The caller of this copy
+   * then assumes ownership of it.
    */
   virtual
   Mapping<dim,spacedim> *clone () const;
 
   /**
-   * Always returns @p false because
-   * MappingQ1Eulerian does not in general
-   * preserve vertex locations (unless the
-   * translation vector happens to provide
-   * for zero displacements at vertex
-   * locations).
+   * Always returns @p false because MappingQ1Eulerian does not in general
+   * preserve vertex locations (unless the translation vector happens to
+   * provide for zero displacements at vertex locations).
    */
   bool preserves_vertex_locations () const;
 
@@ -146,10 +131,8 @@ public:
 
 protected:
   /**
-   * Implementation of the interface in
-   * MappingQ. Overrides the function in
-   * the base class, since we cannot use
-   * any cell similarity for this class.
+   * Implementation of the interface in MappingQ. Overrides the function in
+   * the base class, since we cannot use any cell similarity for this class.
    */
   virtual void
   fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
@@ -164,36 +147,29 @@ protected:
                   CellSimilarity::Similarity                           &cell_similarity) const;
 
   /**
-   * Reference to the vector of
-   * shifts.
+   * Reference to the vector of shifts.
    */
 
   SmartPointer<const VECTOR, MappingQEulerian<dim,VECTOR,spacedim> > euler_vector;
 
   /**
-   * Pointer to the DoFHandler to
-   * which the mapping vector is
-   * associated.
+   * Pointer to the DoFHandler to which the mapping vector is associated.
    */
-
   SmartPointer<const DoFHandler<dim,spacedim>,MappingQEulerian<dim,VECTOR,spacedim> > euler_dof_handler;
 
 
 private:
 
   /**
-   * Special quadrature rule used
-   * to define the support points
-   * in the reference configuration.
+   * Special quadrature rule used to define the support points in the
+   * reference configuration.
    */
 
   class SupportQuadrature : public Quadrature<dim>
   {
   public:
     /**
-     * Constructor, with an argument
-     * defining the desired polynomial
-     * degree.
+     * Constructor, with an argument defining the desired polynomial degree.
      */
 
     SupportQuadrature (const unsigned int map_degree);
@@ -201,36 +177,26 @@ private:
   };
 
   /**
-   * A member variable holding the
-   * quadrature points in the right
-   * order.
+   * A member variable holding the quadrature points in the right order.
    */
   const SupportQuadrature support_quadrature;
 
   /**
-   * FEValues object used to query the
-   * the given finite element field
-   * at the support points in the
-   * reference configuration.
+   * FEValues object used to query the the given finite element field at the
+   * support points in the reference configuration.
    *
-   * The variable is marked as
-   * mutable since we have to call
-   * FEValues::reinit from
-   * compute_mapping_support_points,
-   * a function that is 'const'.
+   * The variable is marked as mutable since we have to call FEValues::reinit
+   * from compute_mapping_support_points, a function that is 'const'.
    */
   mutable FEValues<dim,spacedim> fe_values;
 
   /**
-   * A variable to guard access to
-   * the fe_values variable.
+   * A variable to guard access to the fe_values variable.
    */
   mutable Threads::Mutex fe_values_mutex;
 
   /**
-   * Compute the positions of the
-   * support points in the current
-   * configuration
+   * Compute the positions of the support points in the current configuration
    */
   virtual void compute_mapping_support_points(
     const typename Triangulation<dim,spacedim>::cell_iterator &cell,
index 72b71db88f3b1686ab5ece69a927d98ff12488ce..b478da200952aa36a33704af2ba3d699cfe518ed 100644 (file)
@@ -1566,7 +1566,7 @@ public:
     virtual
     void
     create_notification (const Triangulation<dim, spacedim> &tria);
-  } DEAL_II_DEPRECATED;
+  };
 
   /**
    * A structure that is used as an
index ba830652009fd7ca89fa0cb4ff92d06d50448ad7..ca948550fc7d75638a145bea0bcb8ed527f39ab9 100644 (file)
@@ -22,6 +22,8 @@
 
 #include <deal.II/base/config.h>
 #include <deal.II/base/subscriptor.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/thread_management.h>
 #include <deal.II/base/point.h>
 #include <deal.II/grid/tria.h>
 
@@ -85,120 +87,75 @@ class Boundary : public Subscriptor
 public:
 
   /**
-   * Type keeping information about
-   * the normals at the vertices of
-   * a face of a cell. Thus, there
-   * are
-   * <tt>GeometryInfo<dim>::vertices_per_face</tt>
-   * normal vectors, that define
-   * the tangent spaces of the
-   * boundary at the vertices. Note
-   * that the vectors stored in
-   * this object are not required
-   * to be normalized, nor to
-   * actually point outward, as one
-   * often will only want to check
-   * for orthogonality to define
-   * the tangent plane; if a
-   * function requires the normals
-   * to be normalized, then it must
-   * do so itself.
+   * Type keeping information about the normals at the vertices of a face of a
+   * cell. Thus, there are <tt>GeometryInfo<dim>::vertices_per_face</tt>
+   * normal vectors, that define the tangent spaces of the boundary at the
+   * vertices. Note that the vectors stored in this object are not required to
+   * be normalized, nor to actually point outward, as one often will only want
+   * to check for orthogonality to define the tangent plane; if a function
+   * requires the normals to be normalized, then it must do so itself.
    *
-   * For obvious reasons, this
-   * type is not useful in 1d.
+   * For obvious reasons, this type is not useful in 1d.
    */
   typedef Tensor<1,spacedim> FaceVertexNormals[GeometryInfo<dim>::vertices_per_face];
 
   /**
-   * Destructor. Does nothing here, but
-   * needs to be declared to make it
+   * Destructor. Does nothing here, but needs to be declared to make it
    * virtual.
    */
   virtual ~Boundary ();
 
   /**
-   * Return the point which shall
-   * become the new middle vertex
-   * of the two children of a
-   * regular line. In 2D, this line
-   * is a line at the boundary,
-   * while in 3d, it is bounding a
-   * face at the boundary (the
-   * lines therefore is also on the
-   * boundary).
+   * Return the point which shall become the new middle vertex of the two
+   * children of a regular line. In 2D, this line is a line at the boundary,
+   * while in 3d, it is bounding a face at the boundary (the lines therefore
+   * is also on the boundary).
    */
   virtual
   Point<spacedim>
   get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const = 0;
 
   /**
-   * Return the point which shall
-   * become the common point of the
-   * four children of a quad at the
-   * boundary in three or more
-   * spatial dimensions. This
-   * function therefore is only
-   * useful in at least three
-   * dimensions and should not be
-   * called for lower dimensions.
+   * Return the point which shall become the common point of the four children
+   * of a quad at the boundary in three or more spatial dimensions. This
+   * function therefore is only useful in at least three dimensions and should
+   * not be called for lower dimensions.
    *
-   * This function is called after
-   * the four lines bounding the
-   * given @p quad are refined, so
-   * you may want to use the
-   * information provided by
-   * <tt>quad->line(i)->child(j)</tt>,
-   * <tt>i=0...3</tt>, <tt>j=0,1</tt>.
+   * This function is called after the four lines bounding the given @p quad
+   * are refined, so you may want to use the information provided by
+   * <tt>quad->line(i)->child(j)</tt>, <tt>i=0...3</tt>, <tt>j=0,1</tt>.
    *
-   * Because in 2D, this function
-   * is not needed, it is not made
-   * pure virtual, to avoid the
-   * need to overload it.  The
-   * default implementation throws
-   * an error in any case, however.
+   * Because in 2D, this function is not needed, it is not made pure virtual,
+   * to avoid the need to overload it.  The default implementation throws an
+   * error in any case, however.
    */
   virtual
   Point<spacedim>
   get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
 
   /**
-   * Depending on <tt>dim=2</tt> or
-   * <tt>dim=3</tt> this function
-   * calls the
-   * get_new_point_on_line or the
-   * get_new_point_on_quad
-   * function. It throws an
-   * exception for
-   * <tt>dim=1</tt>. This wrapper
-   * allows dimension independent
+   * Depending on <tt>dim=2</tt> or <tt>dim=3</tt> this function calls the
+   * get_new_point_on_line or the get_new_point_on_quad function. It throws an
+   * exception for <tt>dim=1</tt>. This wrapper allows dimension independent
    * programming.
    */
   Point<spacedim>
   get_new_point_on_face (const typename Triangulation<dim,spacedim>::face_iterator &face) const;
 
   /**
-   * Return equally spaced
-   * intermediate points on a line.
+   * Return intermediate points on a line spaced according to the interior
+   * support points of the 1D Gauss-Lobatto quadrature formula.
    *
-   * The number of points requested
-   * is given by the size of the
-   * vector @p points. It is the
-   * task of the derived classes to
-   * arrange the points in
+   * The number of points requested is given by the size of the vector @p
+   * points. It is the task of the derived classes to arrange the points in
    * approximately equal distances.
    *
-   * This function is called by the
-   * @p MappingQ class. This
-   * happens on each face line of a
-   * cells that has got at least
-   * one boundary line.
+   * This function is called by the @p MappingQ class. This happens on each
+   * face line of a cells that has got at least one boundary line.
    *
-   * As this function is not needed
-   * for @p MappingQ1, it is not
-   * made pure virtual, to avoid
-   * the need to overload it.  The
-   * default implementation throws
-   * an error in any case, however.
+   * As this function is not needed for @p MappingQ1, it is not made pure
+   * virtual, to avoid the need to overload it.  The default implementation
+   * throws an error in any case, however.
    */
   virtual
   void
@@ -206,34 +163,23 @@ public:
                                    std::vector<Point<spacedim> > &points) const;
 
   /**
-   * Return equally spaced
-   * intermediate points on a
-   * boundary quad.
+   * Return intermediate points on a line spaced according to the tensor
+   * product of the interior support points of the 1D Gauss-Lobatto quadrature
+   * formula.
    *
-   * The number of points requested
-   * is given by the size of the
-   * vector @p points. It is
-   * required that this number is a
-   * square of another integer,
-   * i.e. <tt>n=points.size()=m*m</tt>. It
-   * is the task of the derived
-   * classes to arrange the points
-   * such they split the quad into
-   * <tt>(m+1)(m+1)</tt> approximately
-   * equal-sized subquads.
+   * The number of points requested is given by the size of the vector @p
+   * points. It is required that this number is a square of another integer,
+   * i.e. <tt>n=points.size()=m*m</tt>. It is the task of the derived classes
+   * to arrange the points such they split the quad into <tt>(m+1)(m+1)</tt>
+   * approximately equal-sized subquads.
    *
-   * This function is called by the
-   * <tt>MappingQ<3></tt> class. This
-   * happens each face quad of
-   * cells in 3d that has got at
-   * least one boundary face quad.
+   * This function is called by the <tt>MappingQ<3></tt> class. This happens
+   * each face quad of cells in 3d that has got at least one boundary face
+   * quad.
    *
-   * As this function is not needed
-   * for @p MappingQ1, it is not
-   * made pure virtual, to avoid
-   * the need to overload it.  The
-   * default implementation throws
-   * an error in any case, however.
+   * As this function is not needed for @p MappingQ1, it is not made pure
+   * virtual, to avoid the need to overload it.  The default implementation
+   * throws an error in any case, however.
    */
   virtual
   void
@@ -241,56 +187,36 @@ public:
                                    std::vector<Point<spacedim> > &points) const;
 
   /**
-   * Depending on <tt>dim=2</tt> or
-   * <tt>dim=3</tt> this function
-   * calls the
-   * get_intermediate_points_on_line
-   * or the
-   * get_intermediate_points_on_quad
-   * function. It throws an
-   * exception for
-   * <tt>dim=1</tt>. This wrapper
-   * allows dimension independent
-   * programming.
+   * Depending on <tt>dim=2</tt> or <tt>dim=3</tt> this function calls the
+   * get_intermediate_points_on_line or the get_intermediate_points_on_quad
+   * function. It throws an exception for <tt>dim=1</tt>. This wrapper allows
+   * dimension independent programming.
    */
   void
   get_intermediate_points_on_face (const typename Triangulation<dim,spacedim>::face_iterator &face,
                                    std::vector<Point<spacedim> > &points) const;
 
   /**
-   * Return the normal vector to the surface
-   * at the point p. If p is not in fact
-   * on the surface, but only closeby,
-   * try to return something reasonable,
-   * for example the normal vector
-   * at the surface point closest to p.
-   * (The point p will in fact not normally
-   * lie on the actual surface, but rather
-   * be a quadrature point mapped by some
-   * polynomial mapping; the mapped surface,
-   * however, will not usually coincide with
-   * the actual surface.)
+   * Return the normal vector to the surface at the point p. If p is not in
+   * fact on the surface, but only close-by, try to return something
+   * reasonable, for example the normal vector at the surface point closest to
+   * p.  (The point p will in fact not normally lie on the actual surface, but
+   * rather be a quadrature point mapped by some polynomial mapping; the
+   * mapped surface, however, will not usually coincide with the actual
+   * surface.)
    *
-   * The face iterator gives an indication
-   * which face this function is supposed
-   * to compute the normal vector for.
-   * This is useful if the boundary of
-   * the domain is composed of different
-   * nondifferential pieces (for example
-   * when using the StraightBoundary class
-   * to approximate a geometry that is
-   * completely described by the coarse mesh,
-   * with piecewise (bi-)linear components
-   * between the vertices, but where the
-   * boundary may have a kink at the vertices
-   * itself).
+   * The face iterator gives an indication which face this function is
+   * supposed to compute the normal vector for.  This is useful if the
+   * boundary of the domain is composed of different nondifferential pieces
+   * (for example when using the StraightBoundary class to approximate a
+   * geometry that is completely described by the coarse mesh, with piecewise
+   * (bi-)linear components between the vertices, but where the boundary may
+   * have a kink at the vertices itself).
    *
-   * @note Implementations of this function
-   * should be able to assume that the point p
-   * lies within or close to the face described by the
-   * first argument. In turn, callers of this
-   * function should ensure that this is
-   * in fact the case.
+   * @note Implementations of this function should be able to assume that the
+   * point p lies within or close to the face described by the first
+   * argument. In turn, callers of this function should ensure that this is in
+   * fact the case.
    */
   virtual
   Tensor<1,spacedim>
@@ -298,32 +224,17 @@ public:
                  const Point<spacedim> &p) const;
 
   /**
-   * Compute the normal vectors to
-   * the boundary at each vertex of
-   * the given face. It is not
-   * required that the normal
-   * vectors be normed
-   * somehow. Neither is it
-   * required that the normals
-   * actually point outward.
+   * Compute the normal vectors to the boundary at each vertex of the given
+   * face. It is not required that the normal vectors be normed
+   * somehow. Neither is it required that the normals actually point outward.
    *
-   * This function is
-   * needed to compute data for C1
-   * mappings. The default
-   * implementation is to throw an
-   * error, so you need not
-   * overload this function in case
-   * you do not intend to use C1
-   * mappings.
+   * This function is needed to compute data for C1 mappings. The default
+   * implementation is to throw an error, so you need not overload this
+   * function in case you do not intend to use C1 mappings.
    *
-   * Note that when computing
-   * normal vectors at a vertex
-   * where the boundary is not
-   * differentiable, you have to
-   * make sure that you compute the
-   * one-sided limits, i.e. limit
-   * with respect to points inside
-   * the given face.
+   * Note that when computing normal vectors at a vertex where the boundary is
+   * not differentiable, you have to make sure that you compute the one-sided
+   * limits, i.e. limit with respect to points inside the given face.
    */
   virtual
   void
@@ -331,35 +242,20 @@ public:
                            FaceVertexNormals &face_vertex_normals) const;
 
   /**
-   * Given a candidate point and a
-   * line segment characterized by
-   * the iterator, return a point
-   * that lies on the surface
-   * described by this object. This
-   * function is used in some mesh
-   * smoothing algorithms that try
-   * to move around points in order
-   * to improve the mesh quality
-   * but need to ensure that points
-   * that were on the boundary
-   * remain on the boundary.
+   * Given a candidate point and a line segment characterized by the iterator,
+   * return a point that lies on the surface described by this object. This
+   * function is used in some mesh smoothing algorithms that try to move
+   * around points in order to improve the mesh quality but need to ensure
+   * that points that were on the boundary remain on the boundary.
    *
-   * If spacedim==1, then the line
-   * represented by the line
-   * iterator is the entire space
-   * (i.e. it is a cell, not a part
-   * of the boundary), and the
-   * returned point equals the
-   * given input point.
+   * If spacedim==1, then the line represented by the line iterator is the
+   * entire space (i.e. it is a cell, not a part of the boundary), and the
+   * returned point equals the given input point.
    *
-   * Derived classes do not need to
-   * implement this function unless
-   * mesh smoothing algorithms are
-   * used with a particular
-   * boundary object. The default
-   * implementation of this
-   * function throws an exception
-   * of type ExcPureFunctionCalled.
+   * Derived classes do not need to implement this function unless mesh
+   * smoothing algorithms are used with a particular boundary object. The
+   * default implementation of this function throws an exception of type
+   * ExcPureFunctionCalled.
    */
   virtual
   Point<spacedim>
@@ -367,19 +263,12 @@ public:
                       const Point<spacedim> &candidate) const;
 
   /**
-   * Same function as above but for
-   * a point that is to be
-   * projected onto the area
-   * characterized by the given
-   * quad.
+   * Same function as above but for a point that is to be projected onto the
+   * area characterized by the given quad.
    *
-   * If spacedim<=2, then the surface
-   * represented by the quad
-   * iterator is the entire space
-   * (i.e. it is a cell, not a part
-   * of the boundary), and the
-   * returned point equals the
-   * given input point.
+   * If spacedim<=2, then the surface represented by the quad iterator is the
+   * entire space (i.e. it is a cell, not a part of the boundary), and the
+   * returned point equals the given input point.
    */
   virtual
   Point<spacedim>
@@ -387,24 +276,40 @@ public:
                       const Point<spacedim> &candidate) const;
 
   /**
-   * Same function as above but for
-   * a point that is to be
-   * projected onto the area
-   * characterized by the given
-   * quad.
+   * Same function as above but for a point that is to be projected onto the
+   * area characterized by the given quad.
    *
-   * If spacedim<=3, then the manifold
-   * represented by the hex
-   * iterator is the entire space
-   * (i.e. it is a cell, not a part
-   * of the boundary), and the
-   * returned point equals the
-   * given input point.
+   * If spacedim<=3, then the manifold represented by the hex iterator is the
+   * entire space (i.e. it is a cell, not a part of the boundary), and the
+   * returned point equals the given input point.
    */
   virtual
   Point<spacedim>
   project_to_surface (const typename Triangulation<dim,spacedim>::hex_iterator &hex,
                       const Point<spacedim> &candidate) const;
+
+protected:
+  /**
+   * Returns the support points of the Gauss-Lobatto quadrature formula used
+   * for intermediate points.
+   *
+   * @note Since the boundary description is closely tied to the unit cell
+   * support points of MappingQ, new boundary descriptions need to explicitly
+   * use these Gauss-Lobatto points and not equidistant points.
+   */
+  const std::vector<Point<1> > &
+  get_line_support_points (const unsigned int n_intermediate_points) const;
+
+private:
+  /**
+   * Point generator for the intermediate points on a boundary.
+   */
+  mutable std::vector<std_cxx1x::shared_ptr<QGaussLobatto<1> > > points;
+
+  /**
+   * Mutex for protecting the points array.
+   */
+  mutable Threads::Mutex mutex;
 };
 
 
@@ -429,56 +334,37 @@ class StraightBoundary : public Boundary<dim,spacedim>
 {
 public:
   /**
-   * Default constructor. Some
-   * compilers require this for
-   * some reasons.
+   * Default constructor. Some compilers require this for some reasons.
    */
   StraightBoundary ();
 
   /**
-   * Let the new point be the
-   * arithmetic mean of the two
-   * vertices of the line.
+   * Let the new point be the arithmetic mean of the two vertices of the line.
    *
-   * Refer to the general
-   * documentation of this class
-   * and the documentation of the
-   * base class for more
-   * information.
+   * Refer to the general documentation of this class and the documentation of
+   * the base class for more information.
    */
   virtual Point<spacedim>
   get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const;
 
   /**
-   * Let the new point be the
-   * arithmetic mean of the four
-   * vertices of this quad and the
-   * four midpoints of the lines,
-   * which are already created at
-   * the time of calling this
-   * function.
+   * Let the new point be the arithmetic mean of the four vertices of this
+   * quad and the four midpoints of the lines, which are already created at
+   * the time of calling this function.
    *
-   * Refer to the general
-   * documentation of this class
-   * and the documentation of the
-   * base class for more
-   * information.
+   * Refer to the general documentation of this class and the documentation of
+   * the base class for more information.
    */
   virtual
   Point<spacedim>
   get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
 
   /**
-   * Gives <tt>n=points.size()</tt>
-   * points that splits the
-   * StraightBoundary line into
-   * $n+1$ partitions of equal
-   * lengths.
+   * Gives <tt>n=points.size()</tt> points that splits the StraightBoundary
+   * line into $n+1$ partitions of equal lengths.
    *
-   * Refer to the general
-   * documentation of this class
-   * and the documentation of the
-   * base class.
+   * Refer to the general documentation of this class and the documentation of
+   * the base class.
    */
   virtual
   void
@@ -486,16 +372,11 @@ public:
                                    std::vector<Point<spacedim> > &points) const;
 
   /**
-   * Gives <tt>n=points.size()=m*m</tt>
-   * points that splits the
-   * StraightBoundary quad into
-   * $(m+1)(m+1)$ subquads of equal
-   * size.
+   * Gives <tt>n=points.size()=m*m</tt> points that splits the
+   * StraightBoundary quad into $(m+1)(m+1)$ subquads of equal size.
    *
-   * Refer to the general
-   * documentation of this class
-   * and the documentation of the
-   * base class.
+   * Refer to the general documentation of this class and the documentation of
+   * the base class.
    */
   virtual
   void
@@ -503,13 +384,10 @@ public:
                                    std::vector<Point<spacedim> > &points) const;
 
   /**
-   * Implementation of the function
-   * declared in the base class.
+   * Implementation of the function declared in the base class.
    *
-   * Refer to the general
-   * documentation of this class
-   * and the documentation of the
-   * base class.
+   * Refer to the general documentation of this class and the documentation of
+   * the base class.
    */
   virtual
   Tensor<1,spacedim>
@@ -517,14 +395,10 @@ public:
                  const Point<spacedim> &p) const;
 
   /**
-   * Compute the normals to the
-   * boundary at the vertices of
-   * the given face.
+   * Compute the normals to the boundary at the vertices of the given face.
    *
-   * Refer to the general
-   * documentation of this class
-   * and the documentation of the
-   * base class.
+   * Refer to the general documentation of this class and the documentation of
+   * the base class.
    */
   virtual
   void
@@ -532,32 +406,18 @@ public:
                            typename Boundary<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const;
 
   /**
-   * Given a candidate point and a
-   * line segment characterized by
-   * the iterator, return a point
-   * that lies on the surface
-   * described by this object. This
-   * function is used in some mesh
-   * smoothing algorithms that try
-   * to move around points in order
-   * to improve the mesh quality
-   * but need to ensure that points
-   * that were on the boundary
-   * remain on the boundary.
+   * Given a candidate point and a line segment characterized by the iterator,
+   * return a point that lies on the surface described by this object. This
+   * function is used in some mesh smoothing algorithms that try to move
+   * around points in order to improve the mesh quality but need to ensure
+   * that points that were on the boundary remain on the boundary.
    *
-   * The point returned is the
-   * projection of the candidate
-   * point onto the line through
-   * the two vertices of the given
-   * line iterator.
+   * The point returned is the projection of the candidate point onto the line
+   * through the two vertices of the given line iterator.
    *
-   * If spacedim==1, then the line
-   * represented by the line
-   * iterator is the entire space
-   * (i.e. it is a cell, not a part
-   * of the boundary), and the
-   * returned point equals the
-   * given input point.
+   * If spacedim==1, then the line represented by the line iterator is the
+   * entire space (i.e. it is a cell, not a part of the boundary), and the
+   * returned point equals the given input point.
    */
   virtual
   Point<spacedim>
@@ -565,26 +425,15 @@ public:
                       const Point<spacedim> &candidate) const;
 
   /**
-   * Same function as above but for
-   * a point that is to be
-   * projected onto the area
-   * characterized by the given
-   * quad.
+   * Same function as above but for a point that is to be projected onto the
+   * area characterized by the given quad.
    *
-   * The point returned is the
-   * projection of the candidate
-   * point onto the bilinear
-   * surface spanned by the four
-   * vertices of the given quad
-   * iterator.
+   * The point returned is the projection of the candidate point onto the
+   * bilinear surface spanned by the four vertices of the given quad iterator.
    *
-   * If spacedim<=2, then the surface
-   * represented by the quad
-   * iterator is the entire space
-   * (i.e. it is a cell, not a part
-   * of the boundary), and the
-   * returned point equals the
-   * given input point.
+   * If spacedim<=2, then the surface represented by the quad iterator is the
+   * entire space (i.e. it is a cell, not a part of the boundary), and the
+   * returned point equals the given input point.
    */
   virtual
   Point<spacedim>
@@ -592,26 +441,16 @@ public:
                       const Point<spacedim> &candidate) const;
 
   /**
-   * Same function as above but for
-   * a point that is to be
-   * projected onto the area
-   * characterized by the given
-   * quad.
+   * Same function as above but for a point that is to be projected onto the
+   * area characterized by the given quad.
    *
-   * The point returned is the
-   * projection of the candidate
-   * point onto the trilinear
-   * manifold spanned by the eight
-   * vertices of the given hex
+   * The point returned is the projection of the candidate point onto the
+   * trilinear manifold spanned by the eight vertices of the given hex
    * iterator.
    *
-   * If spacedim<=3, then the manifold
-   * represented by the hex
-   * iterator is the entire space
-   * (i.e. it is a cell, not a part
-   * of the boundary), and the
-   * returned point equals the
-   * given input point.
+   * If spacedim<=3, then the manifold represented by the hex iterator is the
+   * entire space (i.e. it is a cell, not a part of the boundary), and the
+   * returned point equals the given input point.
    */
   virtual
   Point<spacedim>
index 1873d9c737122c81a3bb2dbdef839006b6dbf461..5a6115220a5fadf5960ce6d4fa60c5c5422872e6 100644 (file)
@@ -53,88 +53,46 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce
   const unsigned int dim = 2;
   std::vector<Point<dim> > line_points (2);
 
-  // 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. if not at the boundary,
-  // get the respective points
-  // from another function
+  // 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. if not at
+  // the boundary, get the respective points from another function
   for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
     {
       const Triangulation<dim>::line_iterator line = cell->line(line_no);
 
       if (line->at_boundary())
         {
-          // first get the normal
-          // vectors at the two
-          // vertices of this line
-          // from the boundary
-          // description
+          // first get the normal vectors at the two vertices of this line
+          // from the boundary description
           const Boundary<dim> &boundary
             = line->get_triangulation().get_boundary(line->boundary_indicator());
 
           Boundary<dim>::FaceVertexNormals face_vertex_normals;
           boundary.get_normals_at_vertices (line, face_vertex_normals);
 
-          // then transform them into
-          // interpolation points for
-          // a cubic polynomial
+          // then transform them into interpolation points for a cubic
+          // polynomial
           //
-          // for this, note that if
-          // we describe the boundary
-          // curve as a polynomial in
-          // tangential coordinate
-          // @p{t=0..1} (along the
-          // line) and @p{s} in
-          // normal direction, then
-          // the cubic mapping is
-          // such that @p{s = a*t**3
-          // + b*t**2 + c*t + d}, and
-          // we want to determine the
-          // interpolation points at
-          // @p{t=1/3} and
-          // @p{t=2/3}. Since at
-          // @p{t=0,1} we want a
-          // vertex which is actually
-          // at the boundary, we know
-          // that @p{d=0} and
-          // @p{a=-b-c}. As
-          // side-conditions, we want
-          // that the derivatives at
-          // @p{t=0} and @p{t=1},
-          // i.e. at the vertices
-          // match those returned by
-          // the boundary. We then
-          // have that
-          // @p{s(1/3)=1/27(2b+8c)}
-          // and
-          // @p{s(2/3)=4/27b+10/27c}.
+          // for this, note that if we describe the boundary curve as a
+          // polynomial in tangential coordinate @p{t=0..1} (along the line)
+          // and @p{s} in normal direction, then the cubic mapping is such
+          // that @p{s = a*t**3 + b*t**2 + c*t + d}, and we want to determine
+          // the interpolation points at @p{t=0.276} and @p{t=0.724}
+          // (Gauss-Lobatto points). Since at @p{t=0,1} we want a vertex which
+          // is actually at the boundary, we know that @p{d=0} and @p{a=-b-c},
+          // which gives @p{s(0.276)} and @{s(0.726)} in terms of @p{b,c}. As
+          // side-conditions, we want that the derivatives at @p{t=0} and
+          // @p{t=1}, i.e. at the vertices match those returned by the
+          // boundary.
           //
-          // The task is then first
-          // to determine the
-          // coefficients from the
-          // tangentials. for that,
-          // first rotate the
-          // tangents of @p{s(t)}
-          // into the global
-          // coordinate system. they
-          // are @p{A (1,c)} and @p{A
-          // (1,-b-2c)} with @p{A} the
-          // rotation matrix, since
-          // the tangentials in the
-          // coordinate system
-          // relative to the line are
-          // @p{(1,c)} and @p{(1,-b-2c)}
-          // at the two vertices,
-          // respectively. We then
-          // have to make sure by
-          // matching @p{b,c} that
-          // these tangentials are
-          // orthogonal to the normals
-          // returned by the boundary
-          // object
+          // The task is then first to determine the coefficients from the
+          // tangentials. for that, first rotate the tangents of @p{s(t)} into
+          // the global coordinate system. they are @p{A (1,c)} and @p{A
+          // (1,-b-2c)} with @p{A} the rotation matrix, since the tangentials
+          // in the coordinate system relative to the line are @p{(1,c)} and
+          // @p{(1,-b-2c)} at the two vertices, respectively. We then have to
+          // make sure by matching @p{b,c} that these tangentials are
+          // orthogonal to the normals returned by the boundary object
           const Tensor<1,2> coordinate_vector = line->vertex(1) - line->vertex(0);
           const double      h                 = std::sqrt(coordinate_vector * coordinate_vector);
           Tensor<1,2> coordinate_axis = coordinate_vector;
@@ -151,18 +109,18 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce
                              -face_vertex_normals[1][0] * std::sin(alpha)))
                            -2*c;
 
+          QGaussLobatto<1> quad_points(4);
+          const double t1 = quad_points.point(1)[0];
+          const double t2 = quad_points.point(2)[0];
+          const double s_t1 = (((-b-c)*t1+b)*t1+c)*t1;
+          const double s_t2 = (((-b-c)*t2+b)*t2+c)*t2;
 
-          // next evaluate the so
-          // determined cubic
-          // polynomial at the points
-          // 1/3 and 2/3, first in
-          // unit coordinates
-          const Point<2> new_unit_points[2] = { Point<2>(1./3., 1./27.*(2*b+8*c)),
-                                                Point<2>(2./3., 4./27.*b+10./27.*c)
+          // next evaluate the so determined cubic polynomial at the points
+          // 1/3 and 2/3, first in unit coordinates
+          const Point<2> new_unit_points[2] = { Point<2>(t1, s_t1),
+                                                Point<2>(t2, s_t2)
                                               };
-          // then transform these
-          // points to real
-          // coordinates by rotating,
+          // then transform these points to real coordinates by rotating,
           // scaling and shifting
           for (unsigned int i=0; i<2; ++i)
             {
index 3269194d14875b619c395b77d5f273d6ceeda602..4ee1d4374c86afca511c663c7f6ffba075b4554d 100644 (file)
@@ -134,7 +134,7 @@ MappingQ<dim,spacedim>::MappingQ (const unsigned int p,
   // Qp mapping of cells at the boundary.
   const QGaussLobatto<1> points(degree+1);
   tensor_pols = new TensorProductPolynomials<dim>
-    (Polynomials::LagrangeEquidistant::generate_complete_basis(degree));
+    (Polynomials::generate_complete_Lagrange_basis(points.get_points()));
   Assert (n_shape_functions==tensor_pols->n(),
           ExcInternalError());
   Assert(n_inner+n_outer==n_shape_functions, ExcInternalError());
@@ -290,10 +290,8 @@ MappingQ<dim,spacedim>::get_subface_data (const UpdateFlags update_flags,
 }
 
 
-// Note that the CellSimilarity flag is
-// modifyable, since MappingQ can need to
-// recalculate data even when cells are
-// similar.
+// Note that the CellSimilarity flag is modifyable, since MappingQ can need to
+// recalculate data even when cells are similar.
 template<int dim, int spacedim>
 void
 MappingQ<dim,spacedim>::fill_fe_values (
@@ -308,35 +306,23 @@ MappingQ<dim,spacedim>::fill_fe_values (
   std::vector<Point<spacedim> >                             &normal_vectors,
   CellSimilarity::Similarity                           &cell_similarity) const
 {
-  // convert data object to internal
-  // data for this class. fails with
-  // an exception if that is not
-  // possible
+  // convert data object to internal data for this class. fails with an
+  // exception if that is not possible
   Assert (dynamic_cast<InternalData *> (&mapping_data) != 0, ExcInternalError());
   InternalData &data = static_cast<InternalData &> (mapping_data);
 
-  // check whether this cell needs
-  // the full mapping or can be
-  // treated by a reduced Q1 mapping,
-  // e.g. if the cell is in the
-  // interior of the domain
+  // check whether this cell needs the full mapping or can be treated by a
+  // reduced Q1 mapping, e.g. if the cell is in the interior of the domain
   data.use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells
                                           || cell->has_boundary_lines());
 
-  // depending on this result, use this or
-  // the other data object for the
-  // mapping. furthermore, we need to
-  // ensure that the flag indicating
-  // whether we can use some similarity has
-  // to be modified - for a general
-  // MappingQ, the data needs to be
-  // recomputed anyway since then the
-  // mapping changes the data. this needs
-  // to be known also for later operations,
-  // so modify the variable here. this also
-  // affects the calculation of the next
-  // cell -- if we use Q1 data on the next
-  // cell, the data will still be invalid.
+  // depending on this result, use this or the other data object for the
+  // mapping. furthermore, we need to ensure that the flag indicating whether
+  // we can use some similarity has to be modified - for a general MappingQ,
+  // the data needs to be recomputed anyway since then the mapping changes the
+  // data. this needs to be known also for later operations, so modify the
+  // variable here. this also affects the calculation of the next cell -- if
+  // we use Q1 data on the next cell, the data will still be invalid.
   typename MappingQ1<dim,spacedim>::InternalData *p_data=0;
   if (data.use_mapping_q1_on_current_cell)
     p_data=&data.mapping_q1_data;
@@ -368,34 +354,23 @@ MappingQ<dim,spacedim>::fill_fe_face_values (
   std::vector<Tensor<1,spacedim> >  &exterior_forms,
   std::vector<Point<spacedim> >     &normal_vectors) const
 {
-  // convert data object to internal
-  // data for this class. fails with
-  // an exception if that is not
-  // possible
+  // convert data object to internal data for this class. fails with an
+  // exception if that is not possible
   Assert (dynamic_cast<InternalData *> (&mapping_data) != 0,
           ExcInternalError());
   InternalData &data = static_cast<InternalData &> (mapping_data);
 
-  // check whether this cell needs
-  // the full mapping or can be
-  // treated by a reduced Q1 mapping,
-  // e.g. if the cell is entirely in
-  // the interior of the domain. note
-  // that it is not sufficient to ask
-  // whether the present _face_ is in
-  // the interior, as the mapping on
-  // the face depends on the mapping
-  // of the cell, which in turn
-  // depends on the fact whether
-  // _any_ of the faces of this cell
-  // is at the boundary, not only the
-  // present face
+  // check whether this cell needs the full mapping or can be treated by a
+  // reduced Q1 mapping, e.g. if the cell is entirely in the interior of the
+  // domain. note that it is not sufficient to ask whether the present _face_
+  // is in the interior, as the mapping on the face depends on the mapping of
+  // the cell, which in turn depends on the fact whether _any_ of the faces of
+  // this cell is at the boundary, not only the present face
   data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells
                                         || cell->has_boundary_lines());
 
-  // depending on this result, use
-  // this or the other data object
-  // for the mapping
+  // depending on this result, use this or the other data object for the
+  // mapping
   typename MappingQ1<dim,spacedim>::InternalData *p_data=0;
   if (data.use_mapping_q1_on_current_cell)
     p_data=&data.mapping_q1_data;
@@ -430,34 +405,23 @@ MappingQ<dim,spacedim>::fill_fe_subface_values (const typename Triangulation<dim
                                                 std::vector<Tensor<1,spacedim> >  &exterior_forms,
                                                 std::vector<Point<spacedim> >     &normal_vectors) const
 {
-  // convert data object to internal
-  // data for this class. fails with
-  // an exception if that is not
-  // possible
+  // convert data object to internal data for this class. fails with an
+  // exception if that is not possible
   Assert (dynamic_cast<InternalData *> (&mapping_data) != 0,
           ExcInternalError());
   InternalData &data = static_cast<InternalData &> (mapping_data);
 
-  // check whether this cell needs
-  // the full mapping or can be
-  // treated by a reduced Q1 mapping,
-  // e.g. if the cell is entirely in
-  // the interior of the domain. note
-  // that it is not sufficient to ask
-  // whether the present _face_ is in
-  // the interior, as the mapping on
-  // the face depends on the mapping
-  // of the cell, which in turn
-  // depends on the fact whether
-  // _any_ of the faces of this cell
-  // is at the boundary, not only the
-  // present face
+  // check whether this cell needs the full mapping or can be treated by a
+  // reduced Q1 mapping, e.g. if the cell is entirely in the interior of the
+  // domain. note that it is not sufficient to ask whether the present _face_
+  // is in the interior, as the mapping on the face depends on the mapping of
+  // the cell, which in turn depends on the fact whether _any_ of the faces of
+  // this cell is at the boundary, not only the present face
   data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells
                                         || cell->has_boundary_lines());
 
-  // depending on this result, use
-  // this or the other data object
-  // for the mapping
+  // depending on this result, use this or the other data object for the
+  // mapping
   typename MappingQ1<dim,spacedim>::InternalData *p_data=0;
   if (data.use_mapping_q1_on_current_cell)
     p_data=&data.mapping_q1_data;
@@ -499,25 +463,19 @@ MappingQ<dim,spacedim>::set_laplace_on_quad_vector(Table<2,double> &loqvs) const
   const unsigned int n_inner_2d=(degree-1)*(degree-1);
   const unsigned int n_outer_2d=4+4*(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
+  // first check whether we have precomputed the values for some polynomial
+  // degree; the sizes of arrays is n_inner_2d*n_outer_2d
   double const *loqv_ptr=0;
   switch (degree)
     {
-      // for degree==1, we shouldn't have to
-      // compute any support points, since
+      // for degree==1, we shouldn't have to compute any support points, since
       // all of them are on the vertices
 
     case 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)
+      // (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.};
       loqv_ptr=&loqv2[0];
@@ -528,133 +486,14 @@ MappingQ<dim,spacedim>::set_laplace_on_quad_vector(Table<2,double> &loqvs) const
       break;
     }
 
-    case 3:
-    {
-      // (same as above)
-      static const double loqv3[4*12]
-        = {80/1053., 1/81., 1/81., 11/1053., 25/117., 44/351.,
-           7/117., 16/351., 25/117., 44/351., 7/117., 16/351.,
-           1/81., 80/1053., 11/1053., 1/81., 7/117., 16/351.,
-           25/117., 44/351., 44/351., 25/117., 16/351., 7/117.,
-           1/81., 11/1053., 80/1053., 1/81., 44/351., 25/117.,
-           16/351., 7/117., 7/117., 16/351., 25/117., 44/351.,
-           11/1053., 1/81., 1/81., 80/1053., 16/351., 7/117.,
-           44/351., 25/117., 16/351., 7/117., 44/351., 25/117.
-          };
-      Assert (sizeof(loqv3)/sizeof(loqv3[0]) ==
-              n_inner_2d * n_outer_2d,
-              ExcInternalError());
-
-      loqv_ptr=&loqv3[0];
-
-      break;
-    }
-
-
-    case 4:
-    {
-      static const double loqv4[9*16]
-        = {0.07405921850311571, -0.001075744628905992,
-           -0.001075744628906007, 0.001914292239071463,
-           0.2231273865431892, 0.1346851306015187,
-           0.03812914216116724, 0.02913160002633252,
-           0.02200737428129396, 0.01600835564431224,
-           0.2231273865431891, 0.1346851306015187,
-           0.03812914216116723, 0.02913160002633253,
-           0.02200737428129391, 0.01600835564431222,
-
-           0.00664803151334206, 0.006648031513342719,
-           0.002873452861657458, 0.002873452861657626,
-           0.07903572682584378, 0.05969238281250031,
-           0.03619864817415824, 0.07903572682584187,
-           0.0596923828124999, 0.03619864817415815,
-           0.1527716818820237, 0.2348152760709273,
-           0.152771681882024, 0.02496269311797778,
-           0.04081948955407129, 0.02496269311797789,
-
-           -0.001075744628906923, 0.07405921850311589,
-           0.001914292239071339, -0.001075744628905884,
-           0.02913160002633509, 0.02200737428129395,
-           0.01600835564431229, 0.2231273865431878,
-           0.1346851306015183, 0.0381291421611672,
-           0.03812914216116729, 0.1346851306015185,
-           0.2231273865431898, 0.01600835564431217,
-           0.02200737428129394, 0.02913160002633262,
-
-           0.006648031513342073, 0.002873452861657473,
-           0.006648031513342726, 0.002873452861657636,
-           0.1527716818820238, 0.2348152760709273,
-           0.152771681882024, 0.02496269311797779,
-           0.04081948955407131, 0.0249626931179779,
-           0.07903572682584376, 0.05969238281250026,
-           0.03619864817415824, 0.07903572682584187,
-           0.0596923828124998, 0.0361986481741581,
-
-           0.01106770833333302, 0.01106770833333336,
-           0.01106770833333337, 0.01106770833333374,
-           0.06770833333333424, 0.1035156250000011,
-           0.0677083333333344, 0.06770833333333376,
-           0.103515624999999, 0.06770833333333399,
-           0.06770833333333422, 0.1035156250000009,
-           0.06770833333333436, 0.0677083333333337,
-           0.1035156249999988, 0.0677083333333339,
-
-           0.002873452861657185, 0.006648031513342362,
-           0.002873452861657334, 0.006648031513343038,
-           0.02496269311797779, 0.04081948955407401,
-           0.02496269311797788, 0.1527716818820234,
-           0.234815276070926, 0.1527716818820237,
-           0.03619864817415819, 0.05969238281250028,
-           0.07903572682584407, 0.03619864817415804,
-           0.05969238281249986, 0.0790357268258422,
-
-           -0.001075744628906913, 0.00191429223907134,
-           0.07405921850311592, -0.001075744628905865,
-           0.03812914216116729, 0.1346851306015185,
-           0.2231273865431899, 0.01600835564431217,
-           0.02200737428129396, 0.02913160002633264,
-           0.02913160002633509, 0.02200737428129391,
-           0.01600835564431228, 0.2231273865431878,
-           0.1346851306015183, 0.03812914216116718,
-
-           0.002873452861657176, 0.002873452861657321,
-           0.006648031513342374, 0.006648031513343037,
-           0.03619864817415817, 0.05969238281250032,
-           0.07903572682584409, 0.03619864817415805,
-           0.05969238281249992, 0.07903572682584221,
-           0.02496269311797776, 0.04081948955407392,
-           0.02496269311797785, 0.1527716818820233,
-           0.2348152760709258, 0.1527716818820236,
-
-           0.001914292239071237, -0.001075744628906803,
-           -0.001075744628906778, 0.07405921850311617,
-           0.01600835564431228, 0.02200737428129401,
-           0.02913160002633524, 0.03812914216116726,
-           0.1346851306015182, 0.2231273865431886,
-           0.01600835564431228, 0.02200737428129397,
-           0.02913160002633523, 0.03812914216116726,
-           0.1346851306015181, 0.2231273865431886,
-          };
-
-      Assert (sizeof(loqv4)/sizeof(loqv4[0]) ==
-              n_inner_2d * n_outer_2d,
-              ExcInternalError());
-
-      loqv_ptr=&loqv4[0];
-
-      break;
-    }
-
-    // no other cases implemented,
-    // so simply fall through
+    // no other cases implemented, so simply fall through
     default:
       break;
     }
 
   if (loqv_ptr!=0)
     {
-      // precomputed. copy values to
-      // the loqvs array
+      // precomputed. copy values to the loqvs array
       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)
@@ -663,33 +502,17 @@ MappingQ<dim,spacedim>::set_laplace_on_quad_vector(Table<2,double> &loqvs) const
   else
     {
       // not precomputed, then do so now
-      if (dim==2)
+      if (dim == 2)
         compute_laplace_vector(loqvs);
-      else
-        // computing the Laplace vector for
-        // faces is not supported in 3d at
-        // present. presumably, doing so
-        // would not be so hard: we would
-        // only have to call the function in
-        // 2d, i.e. the quad(=face) values in
-        // 3d are equal to the quad(=cell)
-        // values in 2d. however, that would
-        // require us to link in the 2d
-        // library, which is kind of awkward
-        // (note that compute_laplace_vector
-        // really makes use of a lot of 2d
-        // stuff, such as FEValues etc). an
-        // alternative would be to precompute
-        // the values of this array for a
-        // couple of higher mapping orders,
-        // pin down their values and insert
-        // them into the array above.
-        Assert (false, ExcNotImplemented());
+      else if (dim == 3)
+        {
+          MappingQ<2,2> mapping_2d(this->degree);
+          loqvs = mapping_2d.laplace_on_quad_vector;
+        }
     }
 
-  // the sum of weights of the points
-  // at the outer rim should be
-  // one. check this
+  // 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*this->degree,
@@ -704,9 +527,8 @@ MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const
 {
   Assert(degree>1, ExcInternalError());
 
-  // first check whether we have
-  // precomputed the values for some
-  // polynomial degree
+  // first check whether we have precomputed the values for some polynomial
+  // degree
   double const *lohv_ptr=0;
   if (degree==2)
     {
@@ -722,8 +544,7 @@ MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const
 
   if (lohv_ptr!=0)
     {
-      // precomputed. copy values to
-      // the lohvs array
+      // precomputed. copy values to the lohvs array
       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)
@@ -733,9 +554,8 @@ MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const
     // not precomputed, then do so now
     compute_laplace_vector(lohvs);
 
-  // the sum of weights of the points
-  // at the outer rim should be
-  // one. check this
+  // 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-12*this->degree*this->degree,
@@ -769,14 +589,11 @@ MappingQ<dim,spacedim>::compute_laplace_vector(Table<2,double> &lvs) const
   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
+  // for degree==1, we shouldn't have to compute any support points, since all
+  // of them are on the vertices
   Assert(degree>1, ExcInternalError());
 
-  // compute the shape
-  // gradients at the quadrature
-  // points on the unit cell
+  // compute the shape gradients at the quadrature points on the unit cell
   const QGauss<dim> quadrature(degree+1);
   const unsigned int n_q_points=quadrature.size();
 
@@ -784,8 +601,7 @@ MappingQ<dim,spacedim>::compute_laplace_vector(Table<2,double> &lvs) const
   quadrature_data.shape_derivatives.resize(n_shape_functions * n_q_points);
   this->compute_shapes(quadrature.get_points(), quadrature_data);
 
-  // Compute the stiffness matrix of
-  // the inner dofs
+  // 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)
@@ -799,8 +615,7 @@ MappingQ<dim,spacedim>::compute_laplace_vector(Table<2,double> &lvs) const
          S(i,j) += res * (long double)quadrature.weight(point);
        }
 
-  // Compute the components of T to be the
-  // product of gradients of inner and
+  // 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)
@@ -823,8 +638,7 @@ MappingQ<dim,spacedim>::compute_laplace_vector(Table<2,double> &lvs) const
   // S:=S_1*T
   S_1.mmult(S_1_T,T);
 
-  // Resize and initialize the
-  // lvs
+  // 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)
@@ -838,16 +652,10 @@ 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
+  // 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(degree));
 
@@ -858,12 +666,9 @@ MappingQ<dim,spacedim>::apply_laplace_vector(const Table<2,double> &lvs,
   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
+  // 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());
@@ -882,18 +687,13 @@ 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 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
+    // 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);
 
@@ -908,8 +708,7 @@ void
 MappingQ<dim,spacedim>::compute_support_points_laplace(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                                                        std::vector<Point<spacedim> > &a) const
 {
-  // in any case, we need the
-  // vertices first
+  // 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);
@@ -921,11 +720,8 @@ MappingQ<dim,spacedim>::compute_support_points_laplace(const typename Triangulat
         add_line_support_points(cell, a);
         break;
       case 2:
-        // in 2d, add the
-        // points on the four
-        // bounding lines to
-        // the exterior (outer)
-        // points
+        // 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);
@@ -935,9 +731,7 @@ MappingQ<dim,spacedim>::compute_support_points_laplace(const typename Triangulat
 
       case 3:
       {
-        // in 3d also add the
-        // points located on
-        // the boundary faces
+        // 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);
@@ -956,8 +750,7 @@ void
 MappingQ<1>::add_line_support_points (const Triangulation<1>::cell_iterator &,
                                       std::vector<Point<1> > &) const
 {
-  // there are no points on bounding
-  // lines which are to be added
+  // there are no points on bounding lines which are to be added
   const unsigned int dim=1;
   Assert (dim > 1, ExcImpossibleInDim(dim));
 }
@@ -971,8 +764,7 @@ MappingQ<1,2>::add_line_support_points (const Triangulation<1,2>::cell_iterator
 {
   const unsigned int dim      = 1;
   const unsigned int spacedim = 2;
-  // Ask for the mid point, if that's
-  // the only thing we need.
+  // Ask for the mid point, if that's the only thing we need.
   if (degree==2)
     {
       const Boundary<dim,spacedim> *const boundary
@@ -980,10 +772,8 @@ MappingQ<1,2>::add_line_support_points (const Triangulation<1,2>::cell_iterator
       a.push_back(boundary->get_new_point_on_line(cell));
     }
   else
-    // otherwise call the more
-    // complicated functions and ask
-    // for inner points from the
-    // boundary description
+    // otherwise call the more complicated functions and ask for inner points
+    // from the boundary description
     {
       std::vector<Point<spacedim> > line_points (degree-1);
 
@@ -1005,8 +795,7 @@ MappingQ<1,3>::add_line_support_points (const Triangulation<1,3>::cell_iterator
 {
   const unsigned int dim      = 1;
   const unsigned int spacedim = 3;
-  // Ask for the mid point, if that's
-  // the only thing we need.
+  // Ask for the mid point, if that's the only thing we need.
   if (degree==2)
     {
       const Boundary<dim,spacedim> *const boundary
@@ -1014,10 +803,8 @@ MappingQ<1,3>::add_line_support_points (const Triangulation<1,3>::cell_iterator
       a.push_back(boundary->get_new_point_on_line(cell));
     }
   else
-    // otherwise call the more
-    // complicated functions and ask
-    // for inner points from the
-    // boundary description
+    // otherwise call the more complicated functions and ask for inner points
+    // from the boundary description
     {
       std::vector<Point<spacedim> > line_points (degree-1);
 
@@ -1038,8 +825,7 @@ MappingQ<dim,spacedim>::add_line_support_points (const typename Triangulation<di
                                                  std::vector<Point<spacedim> > &a) const
 {
   static const StraightBoundary<dim,spacedim> straight_boundary;
-  // if we only need the midpoint,
-  // then ask for it.
+  // if we only need the midpoint, then ask for it.
   if (degree==2)
     {
       for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
@@ -1056,19 +842,13 @@ MappingQ<dim,spacedim>::add_line_support_points (const typename Triangulation<di
         };
     }
   else
-    // otherwise call the more
-    // complicated functions and ask
-    // for inner points from the
-    // boundary description
+    // otherwise call the more complicated functions and ask for inner points
+    // from the boundary description
     {
       std::vector<Point<spacedim> > line_points (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
+      // 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 = cell->line(line_no);
@@ -1083,8 +863,7 @@ MappingQ<dim,spacedim>::add_line_support_points (const typename Triangulation<di
           boundary->get_intermediate_points_on_line (line, line_points);
           if (dim==3)
             {
-              // in 3D, lines might be in wrong
-              // orientation. if so, reverse
+              // 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());
@@ -1092,9 +871,8 @@ MappingQ<dim,spacedim>::add_line_support_points (const typename Triangulation<di
                 a.insert (a.end(), line_points.rbegin(), line_points.rend());
             }
           else
-            // in 2D, lines always have the
-            // correct orientation. simply
-            // append all points
+            // in 2D, lines always have the correct orientation. simply append
+            // all points
             a.insert (a.end(), line_points.begin(), line_points.end());
 
         }
@@ -1115,23 +893,18 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
                      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
+  // used if face quad at boundary or entirely in the interior of the domain
   std::vector<Point<3> > quad_points ((degree-1)*(degree-1));
-  // used if only one line of face
-  // quad is at boundary
+  // used if only one line of face quad is at boundary
   std::vector<Point<3> > b(4*degree);
 
 
-  // loop over all faces and collect
-  // points on them
+  // 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
+      // 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);
@@ -1146,9 +919,7 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
                                                         face_rotation)),
                ExcInternalError());
 
-      // indices of the lines that
-      // bound a face are given by
-      // GeometryInfo<3>::
+      // 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(
@@ -1156,22 +927,16 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
                ExcInternalError());
 #endif
 
-      // if face at boundary, then
-      // ask boundary object to
-      // return intermediate points
-      // on it
+      // if face at boundary, then ask boundary object to return intermediate
+      // points on it
       if (face->at_boundary())
         {
           face->get_triangulation().get_boundary(face->boundary_indicator())
           .get_intermediate_points_on_quad (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.
+          // 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,
@@ -1180,9 +945,8 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
         }
       else
         {
-          // face is not at boundary,
-          // but maybe some of its
-          // lines are. count them
+          // 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())
@@ -1190,37 +954,23 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
 
           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 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
-              // apply_laplace_vector
-              // increases size of b
-              // about 1. There
-              // resize b for the
-              // case the mentioned
-              // function was already
-              // called.
+              // call of function apply_laplace_vector increases size of b
+              // about 1. There resize b for the case the mentioned function
+              // was already called.
               b.resize(4*degree);
 
-              // b is of size
-              // 4*degree, make sure
-              // that this is the
-              // right size
+              // b is of size 4*degree, make sure that this is the right size
               Assert(b.size()==vertices_per_face+lines_per_face*(degree-1),
                      ExcDimensionMismatch(b.size(),
                                           vertices_per_face+lines_per_face*(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
+              // 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)];
@@ -1231,8 +981,7 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
                     a[vertices_per_cell + GeometryInfo<3>::face_to_cell_lines(
                         face_no, i)*(degree-1)+j];
 
-              // Now b includes the support
-              // points on the quad and we can
+              // Now b includes the support points on the quad and we can
               // apply the laplace vector
               apply_laplace_vector(laplace_on_quad_vector, b);
               Assert(b.size()==4*degree+(degree-1)*(degree-1),
@@ -1243,22 +992,14 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
             }
           else
             {
-              // face is entirely in
-              // the interior. get
-              // intermediate points
-              // from a straight
-              // boundary object
+              // face is entirely in the interior. get intermediate points
+              // from a straight boundary object
               straight_boundary.get_intermediate_points_on_quad (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.
+              // 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,
@@ -1307,31 +1048,24 @@ MappingQ<dim,spacedim>::transform (
   const MappingType mapping_type) const
 {
   AssertDimension (input.size(), output.size());
-  // The data object may be just a
-  // MappingQ1::InternalData, so we
-  // have to test for this first.
+  // The data object may be just a MappingQ1::InternalData, so we have to test
+  // for this first.
   const typename MappingQ1<dim,spacedim>::InternalData *q1_data =
     dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *> (&mapping_data);
   Assert(q1_data!=0, ExcInternalError());
 
-  // If it is a genuine
-  // MappingQ::InternalData, we have
-  // to test further
+  // If it is a genuine MappingQ::InternalData, we have to test further
   if (!q1_data->is_mapping_q1_data)
     {
       Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
               ExcInternalError());
       const InternalData &data = static_cast<const InternalData &>(mapping_data);
-      // If we only use the
-      // Q1-portion, we have to
-      // extract that data object
+      // If we only use the Q1-portion, we have to extract that data object
       if (data.use_mapping_q1_on_current_cell)
         q1_data = &data.mapping_q1_data;
     }
-  // Now, q1_data should have the
-  // right tensors in it and we call
-  // the base classes transform
-  // function
+  // Now, q1_data should have the right tensors in it and we call the base
+  // classes transform function
   MappingQ1<dim,spacedim>::transform(input, output, *q1_data, mapping_type);
 }
 
@@ -1346,31 +1080,24 @@ MappingQ<dim,spacedim>::transform (
   const MappingType mapping_type) const
 {
   AssertDimension (input.size(), output.size());
-  // The data object may be just a
-  // MappingQ1::InternalData, so we
-  // have to test for this first.
+  // The data object may be just a MappingQ1::InternalData, so we have to test
+  // for this first.
   const typename MappingQ1<dim,spacedim>::InternalData *q1_data =
     dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *> (&mapping_data);
   Assert(q1_data!=0, ExcInternalError());
 
-  // If it is a genuine
-  // MappingQ::InternalData, we have
-  // to test further
+  // If it is a genuine MappingQ::InternalData, we have to test further
   if (!q1_data->is_mapping_q1_data)
     {
       Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
               ExcInternalError());
       const InternalData &data = static_cast<const InternalData &>(mapping_data);
-      // If we only use the
-      // Q1-portion, we have to
-      // extract that data object
+      // If we only use the Q1-portion, we have to extract that data object
       if (data.use_mapping_q1_on_current_cell)
         q1_data = &data.mapping_q1_data;
     }
-  // Now, q1_data should have the
-  // right tensors in it and we call
-  // the base classes transform
-  // function
+  // Now, q1_data should have the right tensors in it and we call the base
+  // classes transform function
   MappingQ1<dim,spacedim>::transform(input, output, *q1_data, mapping_type);
 }
 
@@ -1383,31 +1110,24 @@ void MappingQ<dim,spacedim>::transform
  const MappingType mapping_type) const
 {
   AssertDimension (input.size(), output.size());
-  // The data object may be just a
-  // MappingQ1::InternalData, so we
-  // have to test for this first.
+  // The data object may be just a MappingQ1::InternalData, so we have to test
+  // for this first.
   const typename MappingQ1<dim,spacedim>::InternalData *q1_data =
     dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *> (&mapping_data);
   Assert(q1_data!=0, ExcInternalError());
 
-  // If it is a genuine
-  // MappingQ::InternalData, we have
-  // to test further
+  // If it is a genuine MappingQ::InternalData, we have to test further
   if (!q1_data->is_mapping_q1_data)
     {
       Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
               ExcInternalError());
       const InternalData &data = static_cast<const InternalData &>(mapping_data);
-      // If we only use the
-      // Q1-portion, we have to
-      // extract that data object
+      // If we only use the Q1-portion, we have to extract that data object
       if (data.use_mapping_q1_on_current_cell)
         q1_data = &data.mapping_q1_data;
     }
-  // Now, q1_data should have the
-  // right tensors in it and we call
-  // the base classes transform
-  // function
+  // Now, q1_data should have the right tensors in it and we call the base
+  // classes transform function
   MappingQ1<dim,spacedim>::transform(input, output, *q1_data, mapping_type);
 }
 
@@ -1418,11 +1138,9 @@ MappingQ<dim,spacedim>::
 transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                              const Point<dim>                                 &p) const
 {
-  // Use the get_data function to
-  // create an InternalData with data
-  // vectors of the right size and
-  // transformation shape values
-  // already computed at point p.
+  // Use the get_data function to create an InternalData with data vectors of
+  // the right size and transformation shape values already computed at point
+  // p.
   const Quadrature<dim> point_quadrature(p);
   std::auto_ptr<InternalData>
   mdata (dynamic_cast<InternalData *> (
@@ -1437,8 +1155,7 @@ transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_it
              &*mdata);
 
   compute_mapping_support_points(cell, p_data->mapping_support_points);
-  // If this should be Q1, ignore all
-  // other support points.
+  // If this should be Q1, ignore all other support points.
   if (p_data->shape_values.size()<p_data->mapping_support_points.size())
     p_data->mapping_support_points.resize
     (GeometryInfo<dim>::vertices_per_cell);
@@ -1455,24 +1172,15 @@ MappingQ<dim,spacedim>::
 transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                              const Point<spacedim>                            &p) const
 {
-  // first a Newton iteration based
-  // on a Q1 mapping to get a good
-  // starting point, the idea being
-  // that this is cheaper than trying
-  // to start with the real mapping
-  // and likely also more robust.
+  // first a Newton iteration based on a Q1 mapping to get a good starting
+  // point, the idea being that this is cheaper than trying to start with the
+  // real mapping and likely also more robust.
   //
-  // that said, this doesn't always
-  // work: there are cases where the
-  // point is outside the cell and
-  // the inverse mapping doesn't
-  // converge. in that case, use the
-  // center point of the cell as a
-  // starting point if we are to go
-  // on using the full mapping, or
-  // just propagate up the exception
-  // if we had no intention of
-  // continuing with the full mapping
+  // that said, this doesn't always work: there are cases where the point is
+  // outside the cell and the inverse mapping doesn't converge. in that case,
+  // use the center point of the cell as a starting point if we are to go on
+  // using the full mapping, or just propagate up the exception if we had no
+  // intention of continuing with the full mapping
   Point<dim> initial_p_unit;
   try
     {
@@ -1481,11 +1189,8 @@ transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_it
     }
   catch (const typename Mapping<dim,spacedim>::ExcTransformationFailed &)
     {
-      // mirror the conditions of the
-      // code below to determine if
-      // we need to use an arbitrary
-      // starting point or if we just
-      // need to rethrow the
+      // mirror the conditions of the code below to determine if we need to
+      // use an arbitrary starting point or if we just need to rethrow the
       // exception
       if (cell->has_boundary_lines()
           ||
@@ -1500,30 +1205,22 @@ transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_it
         throw;
     }
 
-  // then a Newton iteration based on the
-  // full MappingQ if we need this. note that
-  // for interior cells with dim==spacedim,
-  // the mapping used is in fact a Q1
-  // mapping, so there is nothing we need to
-  // do unless the iteration above failed
+  // then a Newton iteration based on the full MappingQ if we need this. note
+  // that for interior cells with dim==spacedim, the mapping used is in fact a
+  // Q1 mapping, so there is nothing we need to do unless the iteration above
+  // failed
   if (cell->has_boundary_lines()
       ||
       use_mapping_q_on_all_cells
       ||
       (dim!=spacedim) )
     {
-      // use the full mapping. in case the
-      // function above should have given us
-      // something back that lies outside the
-      // unit cell (that might happen because
-      // we may have given a point 'p' that
-      // lies inside the cell with the higher
-      // order mapping, but outside the
-      // Q1-mapped reference cell), then
-      // project it back into the reference
-      // cell in hopes that this gives a
-      // better starting point to the
-      // following iteration
+      // use the full mapping. in case the function above should have given us
+      // something back that lies outside the unit cell (that might happen
+      // because we may have given a point 'p' that lies inside the cell with
+      // the higher order mapping, but outside the Q1-mapped reference cell),
+      // then project it back into the reference cell in hopes that this gives
+      // a better starting point to the following iteration
       initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
 
       const Quadrature<dim> point_quadrature(initial_p_unit);
@@ -1539,13 +1236,11 @@ transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_it
 
       compute_mapping_support_points (cell, mdata->mapping_support_points);
 
-      // If this is a q1 mapping,
-      // then only use the support
-      // points on the vertices.
+      // If this is a q1 mapping, then only use the support points on the
+      // vertices.
       if (mdata->shape_values.size() < mdata->mapping_support_points.size())
         mdata->mapping_support_points.resize(GeometryInfo<dim>::vertices_per_cell);
 
-
       return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit, *mdata);
     }
   else
index c8ba7374746cc74fef9a4a8b1b01953d25fd7886..902c6cd3592a3d6cd24ad237a43ce77ad841aaa5 100644 (file)
@@ -71,20 +71,15 @@ SupportQuadrature (const unsigned int map_degree)
   :
   Quadrature<dim>(Utilities::fixed_power<dim>(map_degree+1))
 {
-  // first we determine the support points
-  // on the unit cell in lexicographic order.
-  // for this purpose we can use an interated
-  // trapezoidal quadrature rule.
-  const QTrapez<1> q1d;
-  const QIterated<dim> q_iterated(q1d,map_degree);
+  // first we determine the support points on the unit cell in lexicographic
+  // order, which are (in accordance with MappingQ) the support points of
+  // QGaussLobatto.
+  const QGaussLobatto<dim> q_iterated(map_degree+1);
   const unsigned int n_q_points = q_iterated.size();
 
-  // we then need to define a renumbering
-  // vector that allows us to go from a
-  // lexicographic numbering scheme to a hierarchic
-  // one.  this fragment is taking almost verbatim
-  // from the MappingQ class.
-
+  // we then need to define a renumbering vector that allows us to go from a
+  // lexicographic numbering scheme to a hierarchic one.  this fragment is
+  // taking almost verbatim from the MappingQ class.
   std::vector<unsigned int> renumber(n_q_points);
   std::vector<unsigned int> dpo(dim+1, 1U);
   for (unsigned int i=1; i<dpo.size(); ++i)
@@ -93,9 +88,7 @@ SupportQuadrature (const unsigned int map_degree)
   FETools::lexicographic_to_hierarchic_numbering (
     FiniteElementData<dim> (dpo, 1, map_degree), renumber);
 
-  // finally we assign the quadrature points in the
-  // required order.
-
+  // finally we assign the quadrature points in the required order.
   for (unsigned int q=0; q<n_q_points; ++q)
     this->quadrature_points[renumber[q]] = q_iterated.point(q);
 }
@@ -112,43 +105,30 @@ compute_mapping_support_points
  std::vector<Point<spacedim> > &a) const
 {
 
-  // first, basic assertion
-  // with respect to vector size,
+  // first, basic assertion with respect to vector size,
 
   const types::global_dof_index n_dofs  = euler_dof_handler->n_dofs();
   const types::global_dof_index vector_size = euler_vector->size();
 
   AssertDimension(vector_size,n_dofs);
 
-  // we then transform our tria iterator
-  // into a dof iterator so we can
-  // access data not associated with
-  // triangulations
+  // we then transform our tria iterator into a dof iterator so we can access
+  // data not associated with triangulations
   typename DoFHandler<dim,spacedim>::cell_iterator dof_cell(*cell, euler_dof_handler);
 
   Assert (dof_cell->active() == true, ExcInactiveCell());
 
-  // our quadrature rule is chosen
-  // so that each quadrature point
-  // corresponds to a support point
-  // in the undeformed configuration.
-  // we can then query the given
-  // displacement field at these points
-  // to determine the shift vector that
-  // maps the support points to the
-  // deformed configuration.
-
-  // we assume that the given field contains
-  // dim displacement components, but
-  // that there may be other solution
-  // components as well (e.g. pressures).
-  // this class therefore assumes that the
-  // first dim components represent the
-  // actual shift vector we need, and simply
-  // ignores any components after that.
-  // this implies that the user should order
-  // components appropriately, or create a
-  // separate dof handler for the displacements.
+  // our quadrature rule is chosen so that each quadrature point corresponds
+  // to a support point in the undeformed configuration.  we can then query
+  // the given displacement field at these points to determine the shift
+  // vector that maps the support points to the deformed configuration.
+
+  // we assume that the given field contains dim displacement components, but
+  // that there may be other solution components as well (e.g. pressures).
+  // this class therefore assumes that the first dim components represent the
+  // actual shift vector we need, and simply ignores any components after
+  // that.  this implies that the user should order components appropriately,
+  // or create a separate dof handler for the displacements.
 
   const unsigned int n_support_pts = support_quadrature.size();
   const unsigned int n_components  = euler_dof_handler->get_fe().n_components();
@@ -157,20 +137,15 @@ compute_mapping_support_points
 
   std::vector<Vector<double> > shift_vector(n_support_pts,Vector<double>(n_components));
 
-  // fill shift vector for each
-  // support point using an fe_values
-  // object. make sure that the
-  // fe_values variable isn't used
-  // simulatenously from different
+  // fill shift vector for each support point using an fe_values object. make
+  // sure that the fe_values variable isn't used simulatenously from different
   // threads
   Threads::Mutex::ScopedLock lock(fe_values_mutex);
   fe_values.reinit(dof_cell);
   fe_values.get_function_values(*euler_vector, shift_vector);
 
-  // and finally compute the positions of the
-  // support points in the deformed
+  // and finally compute the positions of the support points in the deformed
   // configuration.
-
   a.resize(n_support_pts);
   for (unsigned int q=0; q<n_support_pts; ++q)
     {
@@ -196,8 +171,7 @@ MappingQEulerian<dim,EulerVectorType,spacedim>::fill_fe_values (
   std::vector<Point<spacedim> >                             &normal_vectors,
   CellSimilarity::Similarity                           &cell_similarity) const
 {
-  // disable any previously detected
-  // similarity and hand on to the respective
+  // disable any previously detected similarity and hand on to the respective
   // function of the base class.
   cell_similarity = CellSimilarity::invalid_next_cell;
   MappingQ<dim,spacedim>::fill_fe_values (cell, q, mapping_data,
index f956c1d16c3a4886733683590c5c18461597a794..65b3bafe7bb4e26c1320f961a9332e7f73e08ae8 100644 (file)
@@ -247,6 +247,33 @@ project_to_surface (const typename Triangulation<dim, spacedim>::hex_iterator &,
 }
 
 
+
+template <int dim, int spacedim>
+const std::vector<Point<1> > &
+Boundary<dim,spacedim>::
+get_line_support_points (const unsigned int n_intermediate_points) const
+{
+  if (points.size() <= n_intermediate_points ||
+      points[n_intermediate_points].get() == 0)
+    {
+      Threads::Mutex::ScopedLock lock(mutex);
+      if (points.size() <= n_intermediate_points)
+        points.resize(n_intermediate_points+1);
+
+      // another thread might have created points in the meantime
+      if (points[n_intermediate_points].get() == 0)
+        {
+          std_cxx1x::shared_ptr<QGaussLobatto<1> >
+            quadrature (new QGaussLobatto<1>(n_intermediate_points+2));
+          points[n_intermediate_points] = quadrature;
+        }
+    }
+  return points[n_intermediate_points]->get_points();
+}
+
+
+
+
 /* -------------------------- StraightBoundary --------------------- */
 
 
@@ -394,15 +421,18 @@ get_intermediate_points_on_line (const Triangulation<1, 2>::line_iterator &line,
   const unsigned int n=points.size();
   Assert(n>0, ExcInternalError());
 
-  const double dx=1./(n+1);
-  double x=dx;
-
+  // Use interior points of QGaussLobatto quadrature formula support points
+  // for consistency with MappingQ
+  const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
   const Point<spacedim> vertices[2] = { line->vertex(0),
                                         line->vertex(1)
                                       };
 
-  for (unsigned int i=0; i<n; ++i, x+=dx)
-    points[i] = (1-x)*vertices[0] + x*vertices[1];
+  for (unsigned int i=0; i<n; ++i)
+    {
+      const double x = line_points[i+1][0];
+      points[i] = (1-x)*vertices[0] + x*vertices[1];
+    }
 }
 
 
@@ -417,15 +447,19 @@ get_intermediate_points_on_line (const typename Triangulation<dim, spacedim>::li
   const unsigned int n=points.size();
   Assert(n>0, ExcInternalError());
 
-  const double dx=1./(n+1);
-  double x=dx;
+  // Use interior points of QGaussLobatto quadrature formula support points
+  // for consistency with MappingQ
+  const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
 
   const Point<spacedim> vertices[2] = { line->vertex(0),
                                         line->vertex(1)
                                       };
 
-  for (unsigned int i=0; i<n; ++i, x+=dx)
-    points[i] = (1-x)*vertices[0] + x*vertices[1];
+  for (unsigned int i=0; i<n; ++i)
+    {
+      const double x = line_points[1+i][0];
+      points[i] = (1-x)*vertices[0] + x*vertices[1];
+    }
 }
 
 
@@ -455,8 +489,7 @@ get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad,
   // is n a square number
   Assert(m*m==n, ExcInternalError());
 
-  const double ds=1./(m+1);
-  double y=ds;
+  const std::vector<Point<1> > &line_points = this->get_line_support_points(m);
 
   const Point<spacedim> vertices[4] = { quad->vertex(0),
                                         quad->vertex(1),
@@ -464,14 +497,17 @@ get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad,
                                         quad->vertex(3)
                                       };
 
-  for (unsigned int i=0; i<m; ++i, y+=ds)
+  for (unsigned int i=0; i<m; ++i)
     {
-      double x=ds;
-      for (unsigned int j=0; j<m; ++j, x+=ds)
-        points[i*m+j]=((1-x) * vertices[0] +
-                       x     * vertices[1]) * (1-y) +
-                      ((1-x) * vertices[2] +
-                       x     * vertices[3]) * y;
+      const double y=line_points[1+i][0];
+      for (unsigned int j=0; j<m; ++j)
+        {
+          const double x=line_points[1+j][0];
+          points[i*m+j]=((1-x) * vertices[0] +
+                         x     * vertices[1]) * (1-y) +
+                        ((1-x) * vertices[2] +
+                         x     * vertices[3]) * y;
+        }
     }
 }
 
@@ -490,8 +526,7 @@ get_intermediate_points_on_quad (const Triangulation<2,3>::quad_iterator &quad,
   // is n a square number
   Assert(m*m==n, ExcInternalError());
 
-  const double ds=1./(m+1);
-  double y=ds;
+  const std::vector<Point<1> > &line_points = this->get_line_support_points(m);
 
   const Point<spacedim> vertices[4] = { quad->vertex(0),
                                         quad->vertex(1),
@@ -499,14 +534,17 @@ get_intermediate_points_on_quad (const Triangulation<2,3>::quad_iterator &quad,
                                         quad->vertex(3)
                                       };
 
-  for (unsigned int i=0; i<m; ++i, y+=ds)
+  for (unsigned int i=0; i<m; ++i)
     {
-      double x=ds;
-      for (unsigned int j=0; j<m; ++j, x+=ds)
-        points[i*m+j]=((1-x) * vertices[0] +
-                       x     * vertices[1]) * (1-y) +
-                      ((1-x) * vertices[2] +
-                       x     * vertices[3]) * y;
+      const double y=line_points[1+i][0];
+      for (unsigned int j=0; j<m; ++j)
+        {
+          const double x=line_points[1+j][0];
+          points[i*m+j]=((1-x) * vertices[0] +
+                         x     * vertices[1]) * (1-y) +
+                        ((1-x) * vertices[2] +
+                         x     * vertices[3]) * y;
+        }
     }
 }
 
index df1f4e991126968a1602937a0df630f9d5b536c4..1cbd41c77a3848531d5927b8cebe859081c39d2b 100644 (file)
@@ -165,14 +165,14 @@ CylinderBoundary<dim,spacedim>::get_intermediate_points_between_points (
   const unsigned int n=points.size();
   Assert(n>0, ExcInternalError());
 
-  // Do a simple linear interpolation
-  // followed by projection, using
-  // the same algorithm as above
-  const Point<spacedim> ds = (v1-v0) / (n+1);
+  // Do a simple linear interpolation followed by projection, using the same
+  // algorithm as above
+  const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
 
   for (unsigned int i=0; i<n; ++i)
     {
-      const Point<spacedim> middle = v0 + (i+1)*ds;
+      const double x = line_points[i+1][0];
+      const Point<spacedim> middle = (1-x)*v0 + x*v1;
 
       const Point<spacedim> vector_from_axis = (middle-point_on_axis) -
                                                ((middle-point_on_axis) * direction) * direction;
@@ -307,24 +307,23 @@ get_intermediate_points_between_points (const Point<dim> &p0,
 {
   const unsigned int n = points.size ();
   const Point<dim> axis = x_1 - x_0;
-  const Point<dim> dx = (p1 - p0) / (n + 1);
 
   Assert (n > 0, ExcInternalError ());
 
-  for (unsigned int i = 0; i < n; ++i)
+  const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
+
+  for (unsigned int i=0; i<n; ++i)
     {
+      const double x = line_points[i+1][0];
+
       // Compute the current point.
-      const Point<dim> x_i = p0 + (i + 1) * dx;
-      // To project this point on the
-      // boundary of the cone we first
-      // compute the orthogonal
-      // projection of this point onto
-      // the axis of the cone.
+      const Point<dim> x_i = (1-x)*p0 + x*p1;
+      // To project this point on the boundary of the cone we first compute
+      // the orthogonal projection of this point onto the axis of the cone.
       const double c = (x_i - x_0) * axis / axis.square ();
       const Point<dim> x_ip = x_0 + c * axis;
-      // Compute the projection of
-      // the middle point on the
-      // boundary of the cone.
+      // Compute the projection of the middle point on the boundary of the
+      // cone.
       points[i] = x_ip + get_radius (x_ip) *  (x_i - x_ip) / (x_i - x_ip).norm ();
     }
 }
@@ -337,16 +336,11 @@ get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) c
   const Point<dim> axis = x_1 - x_0;
   // Compute the middle point of the line.
   const Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
-  // To project it on the boundary of
-  // the cone we first compute the
-  // orthogonal projection of the
-  // middle point onto the axis of
-  // the cone.
+  // To project it on the boundary of the cone we first compute the orthogonal
+  // projection of the middle point onto the axis of the cone.
   const double c = (middle - x_0) * axis / axis.square ();
   const Point<dim> middle_p = x_0 + c * axis;
-  // Compute the projection of the
-  // middle point on the boundary
-  // of the cone.
+  // Compute the projection of the middle point on the boundary of the cone.
   return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
 }
 
@@ -360,20 +354,14 @@ get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
   const int dim = 3;
 
   const Point<dim> axis = x_1 - x_0;
-  // Compute the middle point of the
-  // quad.
+  // Compute the middle point of the quad.
   const Point<dim> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
-  // Same algorithm as above: To
-  // project it on the boundary of
-  // the cone we first compute the
-  // orthogonal projection of the
-  // middle point onto the axis of
-  // the cone.
+  // Same algorithm as above: To project it on the boundary of the cone we
+  // first compute the orthogonal projection of the middle point onto the axis
+  // of the cone.
   const double c = (middle - x_0) * axis / axis.square ();
   const Point<dim> middle_p = x_0 + c * axis;
-  // Compute the projection of the
-  // middle point on the boundary
-  // of the cone.
+  // Compute the projection of the middle point on the boundary of the cone.
   return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
 }
 
@@ -475,15 +463,12 @@ get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
 
   for (unsigned int vertex = 0; vertex < GeometryInfo<dim>::vertices_per_cell; ++vertex)
     {
-      // Compute the orthogonal
-      // projection of the vertex onto
-      // the axis of the cone.
+      // Compute the orthogonal projection of the vertex onto the axis of the
+      // cone.
       const double c = (face->vertex (vertex) - x_0) * axis / axis.square ();
       const Point<dim> vertex_p = x_0 + c * axis;
-      // Then compute the vector
-      // pointing from the point
-      // <tt>vertex_p</tt> on the axis
-      // to the vertex.
+      // Then compute the vector pointing from the point <tt>vertex_p</tt> on
+      // the axis to the vertex.
       const Point<dim> axis_to_vertex = face->vertex (vertex) - vertex_p;
 
       face_vertex_normals[vertex] = axis_to_vertex / axis_to_vertex.norm ();
@@ -629,49 +614,28 @@ HyperBallBoundary<dim,spacedim>::get_intermediate_points_between_points (
   Assert(std::fabs(v1.square()-r2)<eps*r2, ExcInternalError());
 
   const double alpha=std::acos((v0*v1)/std::sqrt(v0.square()*v1.square()));
-  const double d_alpha=alpha/(n+1);
   const Point<spacedim> pm=0.5*(v0+v1);
 
   const double h=std::sqrt(pm.square());
 
-  double beta=0;
-  unsigned int left_index=0, right_index=0;
-  if ((n+1)%2==0)
-    {
-      // if the number of
-      // parts is even insert
-      // the midpoint
-      left_index=(n-1)/2;
-      right_index=left_index;
-      points[left_index]=pm;
-      --left_index;
-      ++right_index;
-      beta=d_alpha;
-    }
-  else
-    {
-      beta=(n/2+1)*d_alpha-alpha/2.;
-      right_index=n/2;
-      left_index=n/2-1;
-    }
-
   // n even:  m=n/2,
   // n odd:   m=(n-1)/2
+  const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
   const unsigned int m=n/2;
-  for (unsigned int i=0; i<m ; ++i, ++right_index, --left_index, beta+=d_alpha)
+  for (unsigned int i=0; i<m ; ++i)
     {
-      Assert(right_index<n, ExcInternalError());
-      Assert(left_index+1>=1, ExcInternalError());
-
-      double d=h*std::tan(beta);
-      points[right_index]=pm+d/length*(v1-v0);
-      points[left_index]=pm-d/length*(v1-v0);
+      const double beta = alpha * (line_points[i+1][0]-0.5);
+      const double d=h*std::tan(beta);
+      points[i]=pm+d/length*(v1-v0);
+      points[n-1-i]=pm-d/length*(v1-v0);
     }
 
+  if ((n+1)%2==0)
+    // if the number of parts is even insert the midpoint
+    points[(n-1)/2]=pm;
 
-  // project the points from the
-  // straight line to the
-  // HyperBallBoundary
+
+  // project the points from the straight line to the HyperBallBoundary
   for (unsigned int i=0; i<n; ++i)
     {
       points[i] *= r / std::sqrt(points[i].square());
@@ -885,10 +849,8 @@ HalfHyperBallBoundary<dim>::
 get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
                                  std::vector<Point<dim> > &points) const
 {
-  // check whether center of object is
-  // at x==0, since then it belongs
-  // to the plane part of the
-  // boundary
+  // check whether center of object is at x==0, since then it belongs to the
+  // plane part of the boundary
   const Point<dim> line_center = line->center();
   if (line_center(0) == this->center(0))
     return StraightBoundary<dim>::get_intermediate_points_on_line (line, points);
@@ -908,10 +870,8 @@ get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterato
     points[0]=get_new_point_on_quad(quad);
   else
     {
-      // check whether center of
-      // object is at x==0, since
-      // then it belongs to the plane
-      // part of the boundary
+      // check whether center of object is at x==0, since then it belongs to
+      // the plane part of the boundary
       const Point<dim> quad_center = quad->center();
       if (quad_center(0) == this->center(0))
         StraightBoundary<dim>::get_intermediate_points_on_quad (quad, points);
@@ -950,10 +910,8 @@ HalfHyperBallBoundary<dim>::
 get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
                          typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
 {
-  // check whether center of object is
-  // at x==0, since then it belongs
-  // to the plane part of the
-  // boundary
+  // check whether center of object is at x==0, since then it belongs to the
+  // plane part of the boundary
   const Point<dim> quad_center = face->center();
   if (quad_center(0) == this->center(0))
     StraightBoundary<dim>::get_normals_at_vertices (face, face_vertex_normals);
@@ -1005,10 +963,8 @@ get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) c
 {
   switch (dim)
     {
-      // in 2d, first check whether the two
-      // end points of the line are on the
-      // axis of symmetry. if so, then return
-      // the mid point
+      // in 2d, first check whether the two end points of the line are on the
+      // axis of symmetry. if so, then return the mid point
     case 2:
     {
       if ((line->vertex(0)(0) == this->center(0))
@@ -1016,17 +972,13 @@ get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) c
           (line->vertex(1)(0) == this->center(0)))
         return (line->vertex(0) + line->vertex(1))/2;
       else
-        // otherwise we are on the outer or
-        // inner part of the shell. proceed
+        // otherwise we are on the outer or inner part of the shell. proceed
         // as in the base class
         return HyperShellBoundary<dim>::get_new_point_on_line (line);
     }
 
-    // in 3d, a line is a straight
-    // line if it is on the symmetry
-    // plane and if not both of its
-    // end points are on either the
-    // inner or outer sphere
+    // in 3d, a line is a straight line if it is on the symmetry plane and if
+    // not both of its end points are on either the inner or outer sphere
     case 3:
     {
 
@@ -1047,8 +999,7 @@ get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) c
                          - outer_radius) < 1e-12 * outer_radius))))
         return (line->vertex(0) + line->vertex(1))/2;
       else
-        // otherwise we are on the outer or
-        // inner part of the shell. proceed
+        // otherwise we are on the outer or inner part of the shell. proceed
         // as in the base class
         return HyperShellBoundary<dim>::get_new_point_on_line (line);
     }
@@ -1079,10 +1030,8 @@ Point<dim>
 HalfHyperShellBoundary<dim>::
 get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
 {
-  // if this quad is on the symmetry plane,
-  // take the center point and project it
-  // outward to the same radius as the
-  // centers of the two radial lines
+  // if this quad is on the symmetry plane, take the center point and project
+  // it outward to the same radius as the centers of the two radial lines
   if ((quad->vertex(0)(0) == this->center(0)) &&
       (quad->vertex(1)(0) == this->center(0)) &&
       (quad->vertex(2)(0) == this->center(0)) &&
@@ -1119,9 +1068,8 @@ get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) c
         Assert (false, ExcInternalError());
     }
 
-  // otherwise we are on the outer or
-  // inner part of the shell. proceed
-  // as in the base class
+  // otherwise we are on the outer or inner part of the shell. proceed as in
+  // the base class
   return HyperShellBoundary<dim>::get_new_point_on_quad (quad);
 }
 
@@ -1135,10 +1083,8 @@ get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterato
 {
   switch (dim)
     {
-      // in 2d, first check whether the two
-      // end points of the line are on the
-      // axis of symmetry. if so, then return
-      // the mid point
+      // in 2d, first check whether the two end points of the line are on the
+      // axis of symmetry. if so, then return the mid point
     case 2:
     {
       if ((line->vertex(0)(0) == this->center(0))
@@ -1146,18 +1092,14 @@ get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterato
           (line->vertex(1)(0) == this->center(0)))
         StraightBoundary<dim>::get_intermediate_points_on_line (line, points);
       else
-        // otherwise we are on the outer or
-        // inner part of the shell. proceed
+        // otherwise we are on the outer or inner part of the shell. proceed
         // as in the base class
         HyperShellBoundary<dim>::get_intermediate_points_on_line (line, points);
       break;
     }
 
-    // in 3d, a line is a straight
-    // line if it is on the symmetry
-    // plane and if not both of its
-    // end points are on either the
-    // inner or outer sphere
+    // in 3d, a line is a straight line if it is on the symmetry plane and if
+    // not both of its end points are on either the inner or outer sphere
     case 3:
     {
       if (((line->vertex(0)(0) == this->center(0))
@@ -1177,8 +1119,7 @@ get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterato
                          - outer_radius) < 1e-12 * outer_radius))))
         StraightBoundary<dim>::get_intermediate_points_on_line (line, points);
       else
-        // otherwise we are on the outer or
-        // inner part of the shell. proceed
+        // otherwise we are on the outer or inner part of the shell. proceed
         // as in the base class
         HyperShellBoundary<dim>::get_intermediate_points_on_line (line, points);
 
@@ -1200,10 +1141,8 @@ get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterato
 {
   Assert (dim < 3, ExcNotImplemented());
 
-  // check whether center of object is
-  // at x==0, since then it belongs
-  // to the plane part of the
-  // boundary
+  // check whether center of object is at x==0, since then it belongs to the
+  // plane part of the boundary
   const Point<dim> quad_center = quad->center();
   if (quad_center(0) == this->center(0))
     StraightBoundary<dim>::get_intermediate_points_on_quad (quad, points);
@@ -1350,9 +1289,9 @@ TorusBoundary<2,3>::get_new_point_on_line (const Triangulation<2,3>::line_iterat
 
   Point<2>  middle(0,0);
 
-  //Take care for periodic conditions,
-  //For instance phi0= 0, phi1= 3/2*Pi  middle has to be 7/4*Pi not 3/4*Pi
-  //This also works for -Pi/2 + Pi, middle is 5/4*Pi
+  //Take care for periodic conditions, For instance phi0= 0, phi1= 3/2*Pi
+  //middle has to be 7/4*Pi not 3/4*Pi. This also works for -Pi/2 + Pi, middle
+  //is 5/4*Pi
   for (unsigned int i=0; i<2; i++)
     if (std::abs(p0(i)-p1(i))> numbers::PI)
       middle(i)=2*numbers::PI;
@@ -1440,7 +1379,7 @@ TorusBoundary<2,3>::
 get_intermediate_points_on_line (const Triangulation<2, 3>::line_iterator   &line,
                                  std::vector< Point< 3 > > &points) const
 {
-  //Almost the same implementation as  StraightBoundary<2,3>
+  //Almost the same implementation as StraightBoundary<2,3>
   unsigned int npoints=points.size();
   if (npoints==0) return;
 
@@ -1453,10 +1392,10 @@ get_intermediate_points_on_line (const Triangulation<2, 3>::line_iterator   &lin
   offset[0]=0;
   offset[1]=0;
 
-  //Take care for periodic conditions & negative angles,
-  //see get_new_point_on_line() above
-  //Because we dont have a symmetric interpolation (just the middle) we need to
-  //add 2*Pi to each almost zero and negative angles.
+  //Take care for periodic conditions & negative angles, see
+  //get_new_point_on_line() above. Because we dont have a symmetric
+  //interpolation (just the middle) we need to add 2*Pi to each almost zero
+  //and negative angles.
   for (unsigned int i=0; i<2; i++)
     for (unsigned int j=1; j<2; j++)
       {
@@ -1473,11 +1412,11 @@ get_intermediate_points_on_line (const Triangulation<2, 3>::line_iterator   &lin
         p[j](i)+=2*numbers::PI*offset[i];
 
 
-  double dx=1.0/(npoints+1);
-  double x=dx;
   Point<2>  target;
-  for (unsigned int i=0; i<npoints; i++,x+=dx)
+  const std::vector<Point<1> > &line_points = this->get_line_support_points(npoints);
+  for (unsigned int i=0; i<npoints; i++)
     {
+      const double x = line_points[i+1][0];
       target=  (1-x)*p[0] + x*p[1];
       points[i]=get_real_coord(target);
     }
@@ -1497,9 +1436,6 @@ get_intermediate_points_on_quad (const Triangulation< 2, 3 >::quad_iterator &qua
   // is n a square number
   Assert(m*m==n, ExcInternalError());
 
-  const double ds=1./(m+1);
-  double y=ds;
-
   Point<2>  p[4];
 
   for (unsigned int i=0; i<4; i++)
@@ -1510,10 +1446,10 @@ get_intermediate_points_on_quad (const Triangulation< 2, 3 >::quad_iterator &qua
   offset[0]=0;
   offset[1]=0;
 
-  //Take care for periodic conditions & negative angles,
-  //see get_new_point_on_line() above
-  //Because we dont have a symmetric interpolation (just the middle) we need to
-  //add 2*Pi to each almost zero and negative angles.
+  //Take care for periodic conditions & negative angles, see
+  //get_new_point_on_line() above.  Because we dont have a symmetric
+  //interpolation (just the middle) we need to add 2*Pi to each almost zero
+  //and negative angles.
   for (unsigned int i=0; i<2; i++)
     for (unsigned int j=1; j<4; j++)
       {
@@ -1529,11 +1465,13 @@ get_intermediate_points_on_quad (const Triangulation< 2, 3 >::quad_iterator &qua
       if (p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles
         p[j](i)+=2*numbers::PI*offset[i];
 
-  for (unsigned int i=0; i<m; ++i, y+=ds)
+  const std::vector<Point<1> > &line_points = this->get_line_support_points(m);
+  for (unsigned int i=0; i<m; ++i)
     {
-      double x=ds;
-      for (unsigned int j=0; j<m; ++j, x+=ds)
+      const double y=line_points[i+1][0];
+      for (unsigned int j=0; j<m; ++j)
         {
+          const double x=line_points[j+1][0];
           target=((1-x) * p[0] +
                   x     * p[1]) * (1-y) +
                  ((1-x) * p[2] +
index dc931c1b9c7682f3cb510d2d56184e4663abc2ee..7a8e9489addf83873a65f595d1b35ce8f38aa3c8 100644 (file)
@@ -1,33 +1,33 @@
 
 DEAL::2d Jacobians:
-DEAL::0.9742 0.0664 -0.0314 0.5822 
-DEAL::0.6722 0.0174 -0.0273 1.0693 
+DEAL::0.9742 0.0664 -0.0313 0.5814 
+DEAL::0.6710 0.0173 -0.0272 1.0693 
 DEAL::0.5858 0 0 0.5858 
-DEAL::0.0314 -0.5822 0.9742 0.0664 
-DEAL::-0.0273 1.0693 -0.6722 -0.0174 
+DEAL::0.0313 -0.5814 0.9742 0.0664 
+DEAL::-0.0272 1.0693 -0.6710 -0.0173 
 DEAL::
 DEAL::2d inverse Jacobians:
-DEAL::1.0228 -0.1167 0.0552 1.7113 
-DEAL::1.4866 -0.0242 0.0379 0.9345 
+DEAL::1.0228 -0.1169 0.0550 1.7137 
+DEAL::1.4893 -0.0241 0.0380 0.9345 
 DEAL::1.7071 0 0 1.7071 
-DEAL::0.1167 1.0228 -1.7113 0.0552 
-DEAL::-0.0242 -1.4866 0.9345 -0.0379 
+DEAL::0.1169 1.0228 -1.7137 0.0550 
+DEAL::-0.0241 -1.4893 0.9345 -0.0380 
 DEAL::
 DEAL::3d Jacobians:
 DEAL::0.4226 0 0 0 0.4226 0 0 0 0.4226 
-DEAL::0.7166 0.0006 0.0633 0.0006 0.7163 -0.0221 -0.0143 0.0053 0.481
-DEAL::0.0204 -0.5582 -0.0364 0.8017 0.0685 0.0034 0.0038 -0.1227 0.8075 
-DEAL::0.8017 0.0685 0.0034 0.0038 -0.1227 0.8075 0.0204 -0.5582 -0.0364 
-DEAL::0.6915 0.0121 0.0583 -0.0304 0.9077 -0.0021 -0.1531 -0.0021 0.9115 
-DEAL::0.8017 0.0685 0.0034 -0.0204 0.5582 0.0364 0.0038 -0.1227 0.8075 
-DEAL::-0.0304 0.9077 -0.0021 -0.6915 -0.0121 -0.0583 -0.1531 -0.0021 0.9115 
+DEAL::0.7165 0.0006 0.0634 0.0006 0.7162 -0.0222 -0.0141 0.0052 0.480
+DEAL::0.0201 -0.5556 -0.0358 0.8018 0.0685 0.0034 0.0038 -0.1227 0.8074 
+DEAL::0.8018 0.0685 0.0034 0.0038 -0.1227 0.8074 0.0201 -0.5556 -0.0358 
+DEAL::0.6859 0.0119 0.0572 -0.0304 0.9075 -0.0021 -0.1524 -0.0021 0.9114 
+DEAL::0.8018 0.0685 0.0034 -0.0201 0.5556 0.0358 0.0038 -0.1227 0.8074 
+DEAL::-0.0304 0.9075 -0.0021 -0.6859 -0.0119 -0.0572 -0.1524 -0.0021 0.9114 
 DEAL::
 DEAL::3d inverse Jacobians:
 DEAL::2.3660 0 0 0 2.3660 0 0 0 2.3660 
-DEAL::1.3918 0.0002 -0.1831 0.0002 1.3957 0.0642 0.0414 -0.0155 2.0720 
-DEAL::0.1523 1.2434 0.0016 -1.7683 0.0455 -0.0800 -0.2694 0.0011 1.2263 
-DEAL::1.2434 0.0016 0.1523 0.0455 -0.0800 -1.7683 0.0011 1.2263 -0.2694 
-DEAL::1.4251 -0.0193 -0.0912 0.0483 1.1010 -0.0005 0.2394 -0.0007 1.0818 
-DEAL::1.2434 -0.1523 0.0016 0.0455 1.7683 -0.0800 0.0011 0.2694 1.2263 
-DEAL::-0.0193 -1.4251 -0.0912 1.1010 -0.0483 -0.0005 -0.0007 -0.2394 1.0818 
+DEAL::1.3921 0.0002 -0.1837 0.0001 1.3959 0.0646 0.0408 -0.0152 2.0764 
+DEAL::0.1531 1.2433 0.0015 -1.7770 0.0449 -0.0790 -0.2707 0.0010 1.2265 
+DEAL::1.2433 0.0015 0.1531 0.0449 -0.0790 -1.7770 0.0010 1.2265 -0.2707 
+DEAL::1.4370 -0.0191 -0.0902 0.0487 1.1012 -0.0005 0.2404 -0.0006 1.0822 
+DEAL::1.2433 -0.1531 0.0015 0.0449 1.7770 -0.0790 0.0010 0.2707 1.2265 
+DEAL::-0.0191 -1.4370 -0.0902 1.1012 -0.0487 -0.0005 -0.0006 -0.2404 1.0822 
 DEAL::
index eef2917ed6c9ef39e4b6de67f8d23c66d262ef10..3cc05ead1ca0dc2ff4deec08c1e744937aa3aab6 100644 (file)
@@ -575,7 +575,7 @@ DEAL:Mapping2d-2-Q2::
 DEAL::  area=6.8
 DEAL::MappingFace2d-2-Q2
 DEAL:MappingFace2d-2-Q2::1.0 1.1       -1.0 0.0        0.35
-DEAL:MappingFace2d-2-Q2::1.0 1.7       -1.0 -1.1e-16   0.65
+DEAL:MappingFace2d-2-Q2::1.0 1.7       -1.0 0.0        0.65
 DEAL:MappingFace2d-2-Q2::1.0 2.3       -1.0 0.0        0.65
 DEAL:MappingFace2d-2-Q2::1.0 2.9       -1.0 0.0        0.35
 DEAL:MappingFace2d-2-Q2::
@@ -587,7 +587,7 @@ DEAL:MappingFace2d-2-Q2::3.2 3.0    0.16 0.99       0.60
 DEAL:MappingFace2d-2-Q2::
 DEAL:MappingFace2d-2-Q2::
 DEAL:MappingFace2d-2-Q2::1.3 1.0       0.0 -1.0        0.70
-DEAL:MappingFace2d-2-Q2::2.3 1.0       -5.6e-17 -1.0   1.3
+DEAL:MappingFace2d-2-Q2::2.3 1.0       0.0 -1.0        1.3
 DEAL:MappingFace2d-2-Q2::3.7 1.0       0.0 -1.0        1.3
 DEAL:MappingFace2d-2-Q2::4.7 1.0       0.0 -1.0        0.70
 DEAL:MappingFace2d-2-Q2::
@@ -599,10 +599,10 @@ DEAL:MappingFace2d-2-Q2::2.9 2.9  -0.38 0.93      0.38
 DEAL:MappingFace2d-2-Q2::
 DEAL:MappingFace2d-2-Q2::
 DEAL::MappingSubface2d-2-Q2
-DEAL:MappingSubface2d-2-Q2::1.0 1.1    -1.0 -2.2e-16
+DEAL:MappingSubface2d-2-Q2::1.0 1.1    -1.0 0.0
 DEAL:MappingSubface2d-2-Q2::1.0 1.3    -1.0 0.0
 DEAL:MappingSubface2d-2-Q2::1.0 1.7    -1.0 0.0
-DEAL:MappingSubface2d-2-Q2::1.0 1.9    -1.0 -1.1e-16
+DEAL:MappingSubface2d-2-Q2::1.0 1.9    -1.0 0.0
 DEAL:MappingSubface2d-2-Q2::
 DEAL:MappingSubface2d-2-Q2::
 DEAL:MappingSubface2d-2-Q2::1.0 2.1    -1.0 0.0
@@ -623,10 +623,10 @@ DEAL:MappingSubface2d-2-Q2::3.6 2.9       0.28 0.96
 DEAL:MappingSubface2d-2-Q2::3.1 3.0    0.13 0.99
 DEAL:MappingSubface2d-2-Q2::
 DEAL:MappingSubface2d-2-Q2::
-DEAL:MappingSubface2d-2-Q2::1.1 1.0    -1.1e-16 -1.0
+DEAL:MappingSubface2d-2-Q2::1.1 1.0    0.0 -1.0
 DEAL:MappingSubface2d-2-Q2::1.7 1.0    0.0 -1.0
 DEAL:MappingSubface2d-2-Q2::2.3 1.0    0.0 -1.0
-DEAL:MappingSubface2d-2-Q2::2.9 1.0    -5.6e-17 -1.0
+DEAL:MappingSubface2d-2-Q2::2.9 1.0    0.0 -1.0
 DEAL:MappingSubface2d-2-Q2::
 DEAL:MappingSubface2d-2-Q2::
 DEAL:MappingSubface2d-2-Q2::3.1 1.0    0.0 -1.0
@@ -672,11 +672,11 @@ DEAL:Mapping2d-2-Q3::1.5 1.6 7.1
 DEAL:Mapping2d-2-Q3::2.0 1.6 7.1
 DEAL:Mapping2d-2-Q3::2.5 1.6 7.4
 DEAL:Mapping2d-2-Q3::3.0 1.6 8.1
-DEAL:Mapping2d-2-Q3::3.6 1.7 9.4
+DEAL:Mapping2d-2-Q3::3.6 1.7 9.3
 DEAL:Mapping2d-2-Q3::4.2 1.8 11.
 DEAL:Mapping2d-2-Q3::4.8 1.9 14.
 DEAL:Mapping2d-2-Q3::
-DEAL:Mapping2d-2-Q3::1.0 1.9 6.7
+DEAL:Mapping2d-2-Q3::1.0 1.9 6.6
 DEAL:Mapping2d-2-Q3::1.5 1.9 6.4
 DEAL:Mapping2d-2-Q3::1.9 1.9 6.3
 DEAL:Mapping2d-2-Q3::2.4 1.9 6.6
@@ -691,7 +691,7 @@ DEAL:Mapping2d-2-Q3::1.9 2.1 5.4
 DEAL:Mapping2d-2-Q3::2.3 2.1 5.6
 DEAL:Mapping2d-2-Q3::2.7 2.2 6.1
 DEAL:Mapping2d-2-Q3::3.2 2.3 7.1
-DEAL:Mapping2d-2-Q3::3.7 2.4 8.8
+DEAL:Mapping2d-2-Q3::3.7 2.4 8.7
 DEAL:Mapping2d-2-Q3::4.2 2.6 11.
 DEAL:Mapping2d-2-Q3::
 DEAL:Mapping2d-2-Q3::1.0 2.4 5.3
@@ -707,7 +707,7 @@ DEAL:Mapping2d-2-Q3::1.0 2.7 4.6
 DEAL:Mapping2d-2-Q3::1.3 2.7 3.8
 DEAL:Mapping2d-2-Q3::1.7 2.6 3.4
 DEAL:Mapping2d-2-Q3::2.0 2.6 3.3
-DEAL:Mapping2d-2-Q3::2.4 2.6 3.4
+DEAL:Mapping2d-2-Q3::2.4 2.6 3.5
 DEAL:Mapping2d-2-Q3::2.7 2.7 3.9
 DEAL:Mapping2d-2-Q3::3.1 2.8 4.7
 DEAL:Mapping2d-2-Q3::3.4 3.0 6.1
@@ -724,21 +724,21 @@ DEAL:Mapping2d-2-Q3::
 DEAL:Mapping2d-2-Q3::
 DEAL::  area=6.8
 DEAL::MappingFace2d-2-Q3
-DEAL:MappingFace2d-2-Q3::1.0 1.1       -1.0 4.4e-16    0.35
+DEAL:MappingFace2d-2-Q3::1.0 1.1       -1.0 -7.8e-16   0.35
 DEAL:MappingFace2d-2-Q3::1.0 1.7       -1.0 0.0        0.65
-DEAL:MappingFace2d-2-Q3::1.0 2.3       -1.0 -6.7e-16   0.65
+DEAL:MappingFace2d-2-Q3::1.0 2.3       -1.0 5.6e-17    0.65
 DEAL:MappingFace2d-2-Q3::1.0 2.9       -1.0 0.0        0.35
 DEAL:MappingFace2d-2-Q3::
 DEAL:MappingFace2d-2-Q3::
-DEAL:MappingFace2d-2-Q3::5.0 1.2       1.0 0.098       0.55
+DEAL:MappingFace2d-2-Q3::5.0 1.2       0.99 0.10       0.55
 DEAL:MappingFace2d-2-Q3::4.7 2.0       0.86 0.51       1.0
 DEAL:MappingFace2d-2-Q3::4.0 2.7       0.51 0.86       1.0
-DEAL:MappingFace2d-2-Q3::3.2 3.0       0.098 1.0       0.55
+DEAL:MappingFace2d-2-Q3::3.2 3.0       0.10 0.99       0.55
 DEAL:MappingFace2d-2-Q3::
 DEAL:MappingFace2d-2-Q3::
-DEAL:MappingFace2d-2-Q3::1.3 1.0       2.2e-16 -1.0    0.70
+DEAL:MappingFace2d-2-Q3::1.3 1.0       -3.9e-16 -1.0   0.70
 DEAL:MappingFace2d-2-Q3::2.3 1.0       0.0 -1.0        1.3
-DEAL:MappingFace2d-2-Q3::3.7 1.0       -3.3e-16 -1.0   1.3
+DEAL:MappingFace2d-2-Q3::3.7 1.0       2.8e-17 -1.0    1.3
 DEAL:MappingFace2d-2-Q3::4.7 1.0       0.0 -1.0        0.70
 DEAL:MappingFace2d-2-Q3::
 DEAL:MappingFace2d-2-Q3::
@@ -749,40 +749,40 @@ DEAL:MappingFace2d-2-Q3::2.9 2.9  -0.39 0.92      0.36
 DEAL:MappingFace2d-2-Q3::
 DEAL:MappingFace2d-2-Q3::
 DEAL::MappingSubface2d-2-Q3
-DEAL:MappingSubface2d-2-Q3::1.0 1.1    -1.0 -4.4e-16
-DEAL:MappingSubface2d-2-Q3::1.0 1.3    -1.0 -2.2e-16
-DEAL:MappingSubface2d-2-Q3::1.0 1.7    -1.0 -4.4e-16
-DEAL:MappingSubface2d-2-Q3::1.0 1.9    -1.0 2.2e-16
+DEAL:MappingSubface2d-2-Q3::1.0 1.1    -1.0 -8.9e-16
+DEAL:MappingSubface2d-2-Q3::1.0 1.3    -1.0 -5.0e-16
+DEAL:MappingSubface2d-2-Q3::1.0 1.7    -1.0 0.0
+DEAL:MappingSubface2d-2-Q3::1.0 1.9    -1.0 -2.2e-16
 DEAL:MappingSubface2d-2-Q3::
 DEAL:MappingSubface2d-2-Q3::
-DEAL:MappingSubface2d-2-Q3::1.0 2.1    -1.0 -8.9e-16
-DEAL:MappingSubface2d-2-Q3::1.0 2.3    -1.0 -1.1e-16
+DEAL:MappingSubface2d-2-Q3::1.0 2.1    -1.0 2.2e-16
+DEAL:MappingSubface2d-2-Q3::1.0 2.3    -1.0 0.0
 DEAL:MappingSubface2d-2-Q3::1.0 2.7    -1.0 0.0
-DEAL:MappingSubface2d-2-Q3::1.0 2.9    -1.0 -4.4e-16
+DEAL:MappingSubface2d-2-Q3::1.0 2.9    -1.0 4.4e-16
 DEAL:MappingSubface2d-2-Q3::
 DEAL:MappingSubface2d-2-Q3::
-DEAL:MappingSubface2d-2-Q3::5.0 1.1    1.0 0.036
-DEAL:MappingSubface2d-2-Q3::4.9 1.5    0.97 0.26
+DEAL:MappingSubface2d-2-Q3::5.0 1.1    1.0 0.039
+DEAL:MappingSubface2d-2-Q3::4.9 1.5    0.96 0.26
 DEAL:MappingSubface2d-2-Q3::4.7 2.0    0.86 0.51
 DEAL:MappingSubface2d-2-Q3::4.5 2.3    0.74 0.67
 DEAL:MappingSubface2d-2-Q3::
 DEAL:MappingSubface2d-2-Q3::
 DEAL:MappingSubface2d-2-Q3::4.3 2.5    0.67 0.74
 DEAL:MappingSubface2d-2-Q3::4.0 2.7    0.51 0.86
-DEAL:MappingSubface2d-2-Q3::3.5 2.9    0.26 0.97
-DEAL:MappingSubface2d-2-Q3::3.1 3.0    0.036 1.0
+DEAL:MappingSubface2d-2-Q3::3.5 2.9    0.26 0.96
+DEAL:MappingSubface2d-2-Q3::3.1 3.0    0.039 1.0
 DEAL:MappingSubface2d-2-Q3::
 DEAL:MappingSubface2d-2-Q3::
-DEAL:MappingSubface2d-2-Q3::1.1 1.0    -2.2e-16 -1.0
-DEAL:MappingSubface2d-2-Q3::1.7 1.0    -1.1e-16 -1.0
-DEAL:MappingSubface2d-2-Q3::2.3 1.0    -2.2e-16 -1.0
-DEAL:MappingSubface2d-2-Q3::2.9 1.0    1.1e-16 -1.0
+DEAL:MappingSubface2d-2-Q3::1.1 1.0    -4.4e-16 -1.0
+DEAL:MappingSubface2d-2-Q3::1.7 1.0    -2.5e-16 -1.0
+DEAL:MappingSubface2d-2-Q3::2.3 1.0    0.0 -1.0
+DEAL:MappingSubface2d-2-Q3::2.9 1.0    -1.1e-16 -1.0
 DEAL:MappingSubface2d-2-Q3::
 DEAL:MappingSubface2d-2-Q3::
-DEAL:MappingSubface2d-2-Q3::3.1 1.0    -4.4e-16 -1.0
-DEAL:MappingSubface2d-2-Q3::3.7 1.0    -5.6e-17 -1.0
+DEAL:MappingSubface2d-2-Q3::3.1 1.0    1.1e-16 -1.0
+DEAL:MappingSubface2d-2-Q3::3.7 1.0    0.0 -1.0
 DEAL:MappingSubface2d-2-Q3::4.3 1.0    0.0 -1.0
-DEAL:MappingSubface2d-2-Q3::4.9 1.0    -2.2e-16 -1.0
+DEAL:MappingSubface2d-2-Q3::4.9 1.0    2.2e-16 -1.0
 DEAL:MappingSubface2d-2-Q3::
 DEAL:MappingSubface2d-2-Q3::
 DEAL:MappingSubface2d-2-Q3::1.1 3.0    0.42 0.91
index ee6d13f6ea79b441170809c35c01268bd7388727..a0152980555335487bf869a886d80296adf9279b 100644 (file)
@@ -29,5 +29,5 @@ cells dofs    area         error
 4     162  3.14159639 1.1900e-06 5.92 
 16    578  3.14159271 1.8860e-08 5.98 
 64    2178 3.14159265 2.9574e-10 5.99 
-256   8450 3.14159265 4.6257e-12 6.00 
+256   8450 3.14159265 4.6356e-12 6.00 
 DEAL::
index 8a46c0b5aa28dfeb573febd2713d824f98ea6bf6..1c2d37f4a521550ebbf54ca835d0b0a40eac0ef3 100644 (file)
@@ -24,6 +24,9 @@
 // initial guess for the inverse Q4 mapping by running an inverse Q1
 // mapping, and we failed to find a solution in the latter.
 //
+// todo: this test used to fail before r30450 but the new mapping support
+// points seem to fix this
+//
 // a redux of this test is in the mapping_real_to_unit_q4_sphere_*
 // tests
 
index d193d0581ee13dac85bc89508803949c51e807c1..8b137891791fe96927ad78e64b0aad7bded08bdc 100644 (file)
@@ -1,2 +1 @@
 
-DEAL::Point is outside!

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.