]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Move unit support points on cells and faces to base class. Initialize them in derived...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 8 Mar 2001 15:18:12 +0000 (15:18 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 8 Mar 2001 15:18:12 +0000 (15:18 +0000)
git-svn-id: https://svn.dealii.org/trunk@4173 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/fe/fe_base.h
deal.II/deal.II/include/fe/fe_dgq.h
deal.II/deal.II/include/fe/fe_q.h
deal.II/deal.II/include/fe/fe_system.h
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_dgq.cc
deal.II/deal.II/source/fe/fe_q.cc
deal.II/deal.II/source/fe/fe_system.cc
deal.II/deal.II/source/fe/fe_tools.cc
deal.II/deal.II/source/numerics/vectors.cc

index 2a86519a752fd424b2d9bbdff6f3d9b67a2308f3..6e4edfc9f72be16019844536a9a7399512b72366 100644 (file)
@@ -67,42 +67,6 @@ class FiniteElement : public FiniteElementBase<dim>
                                      * are deleted properly.
                                      */
     virtual ~FiniteElement ();
-
-                                    /**
-                                     * Return the support points of the
-                                     * trial functions on the unit cell.
-                                     *
-                                     * The order of points in the
-                                     * array matches that returned by
-                                     * the #cell->get_dof_indices#
-                                     * function, but:
-                                     *
-                                     * If the shape functions are not
-                                     * Lagrangian interpolants at
-                                     * some points, the size of the
-                                     * array will be zero after
-                                     * calling this function. This is
-                                     * the standard behavior, if the
-                                     * function is not overloaded.
-                                     */
-    virtual void get_unit_support_points (std::vector<Point<dim> > &points) const;    
-    
-                                    /**
-                                     * Return the support points of
-                                     * the trial functions on the
-                                     * first face of the unit cell.
-                                     *
-                                     * The order of points in
-                                     * the array matches that returned by
-                                     * the #cell->get_dof_indices# function.
-                                     *
-                                     * If the shape functions are not
-                                     * Lagrangian interpolants at some
-                                     * points, the size of the array
-                                     * will be zero. This is the standard behavior,
-                                     * if the function is not overloaded.
-                                     */
-    virtual void get_unit_face_support_points (std::vector<Point<dim-1> > &points) const;    
     
                                     /**
                                      * Number of base elements in a mixed
index c8cca72bf2da0e224a412a8b793608aa4fb8d9f4..518c86c78111dd384d8735a5947e1a086b9c8be4 100644 (file)
@@ -375,18 +375,6 @@ class FiniteElementBase : public Subscriptor,
                                      */
     FiniteElementBase (const FiniteElementData<dim> &fe_data,
                       const std::vector<bool> &restriction_is_additive_flags);
-    
-                                     /**
-                                     * Compute second derivatives by
-                                     * finite differences of
-                                     * gradients.
-                                     */
-    void compute_2nd (const Mapping<dim>                      &mapping,
-                     const DoFHandler<dim>::cell_iterator    &cell,
-                     const unsigned int                       offset,
-                     typename Mapping<dim>::InternalDataBase &mapping_internal,
-                     InternalDataBase                        &fe_internal,
-                     FEValuesData<dim>                       &data) const;
 
                                     /**
                                      * Projection from a fine grid
@@ -503,8 +491,8 @@ class FiniteElementBase : public Subscriptor,
                                     /**
                                      * Compute system index from components.
                                      */
-    unsigned int component_to_system_index (unsigned int component,
-                                           unsigned int component_index) const;
+    unsigned int component_to_system_index (const unsigned int component,
+                                           const unsigned int component_index) const;
   
                                     /**
                                      * Compute component and index from
@@ -515,13 +503,13 @@ class FiniteElementBase : public Subscriptor,
                                      * component.
                                      */
     std::pair<unsigned int,unsigned int>
-    system_to_component_index (unsigned int index) const; 
+    system_to_component_index (const unsigned int index) const; 
     
                                     /**
                                      * Compute system index from components on a face.
                                      */
-    unsigned int face_component_to_system_index (unsigned int component,
-                                                unsigned int component_index) const;
+    unsigned int face_component_to_system_index (const unsigned int component,
+                                                const unsigned int component_index) const;
   
                                     /**
                                      * Compute component and index from system
@@ -532,7 +520,7 @@ class FiniteElementBase : public Subscriptor,
                                      * component.
                                      */
     std::pair<unsigned int,unsigned int>
-    face_system_to_component_index (unsigned int index) const;
+    face_system_to_component_index (const unsigned int index) const;
     
                                     /**
                                      * The base element establishing a
@@ -556,6 +544,109 @@ class FiniteElementBase : public Subscriptor,
                                      */
     bool restriction_is_additive (const unsigned int component) const;
 
+                                    /**
+                                     * Return the support points of
+                                     * the trial functions on the
+                                     * unit cell, if the derived
+                                     * finite element defines some.
+                                     * Finite elements that allow
+                                     * some kind of interpolation
+                                     * operation usually have support
+                                     * points. On the other hand,
+                                     * elements that define their
+                                     * degrees of freedom by, for
+                                     * example, moments on faces, or
+                                     * as derivatives, don't have
+                                     * support points. In that case,
+                                     * the returned field is empty.
+                                     *
+                                     * If the finite element defines
+                                     * support points, then their
+                                     * number equals the number of
+                                     * degrees of freedom of the
+                                     * element.  The order of points
+                                     * in the array matches that
+                                     * returned by the
+                                     * @p{cell->get_dof_indices}
+                                     * function.
+                                     */
+    const std::vector<Point<dim> > & get_unit_support_points () const;    
+
+                                    /**
+                                     * Return whether a finite
+                                     * element has defined support
+                                     * points. If the result is true,
+                                     * then a call to the
+                                     * @p{get_unit_support_points}
+                                     * yields a non-empty array.
+                                     */
+    bool has_support_points () const;
+    
+                                    /**
+                                     * Return the support points of
+                                     * the trial functions on the
+                                     * unit face, if the derived
+                                     * finite element defines some.
+                                     * Finite elements that allow
+                                     * some kind of interpolation
+                                     * operation usually have support
+                                     * points. On the other hand,
+                                     * elements that define their
+                                     * degrees of freedom by, for
+                                     * example, moments on faces, or
+                                     * as derivatives, don't have
+                                     * support points. In that case,
+                                     * the returned field is empty
+                                     *
+                                     * Note that elements that have
+                                     * support points need not
+                                     * necessarily have some on the
+                                     * faces, even if the
+                                     * interpolation points are
+                                     * located physically on a
+                                     * face. For example, the
+                                     * discontinuous elements have
+                                     * interpolation points on the
+                                     * vertices, and for higher
+                                     * degree elements also on the
+                                     * faces, but they are not
+                                     * defined to be on faces since
+                                     * in that case degrees of
+                                     * freedom from both sides of a
+                                     * face (or from all adjacent
+                                     * elements to a vertex) would be
+                                     * identified with each other,
+                                     * which is not what we would
+                                     * like to have). Logically,
+                                     * these degrees of freedom are
+                                     * therefore defined to belong to
+                                     * the cell, rather than the face
+                                     * or vertex. In that case, the
+                                     * returned element would
+                                     * therefore have length zero.
+                                     *
+                                     * If the finite element defines
+                                     * support points, then their
+                                     * number equals the number of
+                                     * degrees of freedom on the face
+                                     * (@p{dofs_per_face}). The order
+                                     * of points in the array matches
+                                     * that returned by the
+                                     * @p{cell->get_dof_indices}
+                                     * function.
+                                     */
+    const std::vector<Point<dim-1> > & get_unit_face_support_points () const;    
+
+                                    /**
+                                     * Return whether a finite
+                                     * element has defined support
+                                     * points on faces. If the result
+                                     * is true, then a call to the
+                                     * @p{get_unit_support_points}
+                                     * yields a non-empty array.
+                                     */
+    bool has_face_support_points () const;
+    
                                     /**
                                      * Determine an estimate for the
                                      * memory consumption (in bytes)
@@ -721,6 +812,49 @@ class FiniteElementBase : public Subscriptor,
                                      */
     const std::vector<bool> restriction_is_additive_flags;
 
+                                    /**
+                                     * List of support points on the
+                                     * unit cell, in case the finite
+                                     * element has any. The
+                                     * constructor leaves this field
+                                     * empty, derived classes may
+                                     * write in some contents.
+                                     *
+                                     * Finite elements that allow
+                                     * some kind of interpolation
+                                     * operation usually have support
+                                     * points. On the other hand,
+                                     * elements that define their
+                                     * degrees of freedom by, for
+                                     * example, moments on faces, or
+                                     * as derivatives, don't have
+                                     * support points. In that case,
+                                     * this field remains empty.
+                                     */
+    std::vector<Point<dim> > unit_support_points;
+
+                                    /**
+                                     * Same for the faces. See the
+                                     * description of the
+                                     * @p{get_unit_face_support_points}
+                                     * function for a discussion of
+                                     * what contributes a face
+                                     * support point.
+                                     */
+    std::vector<Point<dim-1> > unit_face_support_points;
+
+                                     /**
+                                     * Compute second derivatives by
+                                     * finite differences of
+                                     * gradients.
+                                     */
+    void compute_2nd (const Mapping<dim>                      &mapping,
+                     const DoFHandler<dim>::cell_iterator    &cell,
+                     const unsigned int                       offset,
+                     typename Mapping<dim>::InternalDataBase &mapping_internal,
+                     InternalDataBase                        &fe_internal,
+                     FEValuesData<dim>                       &data) const;
+    
                                     /**
                                      * Allow the FESystem class to
                                      * access the restriction and
index aa3ba2f9226e579115d6abbd9892056c5d72d151..dedbb0b76d26c30e15a07cac4f9a0e089385d8fb 100644 (file)
@@ -38,40 +38,6 @@ class FE_DGQ : public FiniteElement<dim>
                                      */
     ~FE_DGQ ();
 
-                                    /**
-                                     * Return the support points of the
-                                     * trial functions on the unit cell.
-                                     *
-                                     * The order of points in
-                                     * the array matches that returned by
-                                     * the #cell->get_dof_indices# function, but:
-                                     *
-                                     * If the shape functions are not
-                                     * Lagrangian interpolants at some
-                                     * points, the size of the array
-                                     * will be zero. This is the standard behavior,
-                                     * if the function is not overloaded.
-                                     */
-    virtual void get_unit_support_points (std::vector<Point<dim> > &) const;    
-    
-                                    /**
-                                     * Return the support points of
-                                     * the trial functions on the
-                                     * first face of the unit cell.
-                                     *
-                                     * The order of points in
-                                     * the array matches that returned by
-                                     * the #cell->get_dof_indices# function.
-                                     *
-                                     * If the shape functions are not
-                                     * Lagrangian interpolants at some
-                                     * points, the size of the array
-                                     * will be zero. This is the standard behavior,
-                                     * if the function is not overloaded.
-                                     */
-    virtual void get_unit_face_support_points (std::vector<Point<dim-1> > &) const;
-
-
                                     /**
                                      * Return the polynomial degree
                                      * of this finite element,
@@ -227,13 +193,6 @@ class FE_DGQ : public FiniteElement<dim>
                                      */
     virtual UpdateFlags update_each (UpdateFlags flags) const;
   
-                                    /**
-                                     * Compute support points, only
-                                     * for @p{degree>0}.
-                                     */
-    static void compute_support_points (std::vector<Point<dim> > &support_points,
-                                       const                     unsigned int degree);
-
                                     /**
                                      * Compute renumbering for rotation
                                      * of degrees of freedom.
index 237ad8a6412e7f9bac5c593ad53a1cc59b98f3cb..04206bf33a4fd901134709287f0b2cd6b7b8d79b 100644 (file)
@@ -40,44 +40,6 @@ class FE_Q : public FiniteElement<dim>
                                      */
     ~FE_Q ();
 
-                                    /**
-                                     * Return the support points of
-                                     * the trial functions on the
-                                     * unit cell.
-                                     *
-                                     * The order of points in the
-                                     * array matches that returned by
-                                     * the @p{cell->get_dof_indices}
-                                     * function, but:
-                                     *
-                                     * If the shape functions are not
-                                     * Lagrangian interpolants at
-                                     * some points, the size of the
-                                     * array will be zero. This is
-                                     * the standard behavior, if the
-                                     * function is not overloaded.
-                                     */
-    virtual void get_unit_support_points (std::vector<Point<dim> > &) const;    
-    
-                                    /**
-                                     * Return the support points of
-                                     * the trial functions on the
-                                     * first face of the unit cell.
-                                     *
-                                     * The order of points in the
-                                     * array matches that returned by
-                                     * the @p{cell->get_dof_indices}
-                                     * function.
-                                     *
-                                     * If the shape functions are not
-                                     * Lagrangian interpolants at
-                                     * some points, the size of the
-                                     * array will be zero. This is
-                                     * the standard behavior, if the
-                                     * function is not overloaded.
-                                     */
-    virtual void get_unit_face_support_points (std::vector<Point<dim-1> > &) const;
-
 //TODO: make get_renumber private or move it some other place    
                                     /**
                                      * Read-only access to the
@@ -304,13 +266,23 @@ class FE_Q : public FiniteElement<dim>
                                        std::vector<unsigned int>      &numbering);
 
                                     /**
-                                     * Compute support points.
+                                     * Initialize the
+                                     * @p{unit_support_points} field
+                                     * of the @ref{FiniteElementBase}
+                                     * class. Called from the
+                                     * constructor.
                                      */
-    static void compute_support_points (const unsigned int               degree,
-                                       const std::vector<unsigned int> &renumber,
-                                       std::vector<Point<dim> >        &support_points);
-
+    void initialize_unit_support_points ();
 
+                                    /**
+                                     * Initialize the
+                                     * @p{unit_face_support_points} field
+                                     * of the @ref{FiniteElementBase}
+                                     * class. Called from the
+                                     * constructor.
+                                     */
+    void initialize_unit_face_support_points ();
+    
                                     /**
                                      * Compute flags for initial
                                      * update only.
index a287fab7007b7af3594954792e52f8bd32a6aaa0..a6246becc20fb28c7bacac774077ffede0a2ccc8 100644 (file)
@@ -170,47 +170,7 @@ class FESystem : public FiniteElement<dim>
                                      * pointers to their base class,
                                      * rather than the class itself.
                                      */
-    virtual unsigned int memory_consumption () const;
-    
-                                    /**
-                                     * Return the support points of
-                                     * the trial functions on the
-                                     * unit cell.
-                                     *
-                                     * The order of points in the
-                                     * array matches that returned by
-                                     * the @p{cell->get_dof_indices}
-                                     * function, but:
-                                     *
-                                     * If one of the base elements
-                                     * has no support points, then it
-                                     * makes no sense to define
-                                     * support points for the
-                                     * composed element, so return an
-                                     * empty array to demonstrate
-                                     * that fact.
-                                     */
-    virtual void get_unit_support_points (typename std::vector<Point<dim> > &) const;    
-
-                                    /**
-                                     * Return the support points of
-                                     * the trial functions on the
-                                     * first face of the unit cell.
-                                     *
-                                     * The order of points in the
-                                     * array matches that returned by
-                                     * the @p{cell->get_dof_indices}
-                                     * function, but:
-                                     *
-                                     * If one of the base elements
-                                     * has no support points, then it
-                                     * makes no sense to define
-                                     * support points for the
-                                     * composed element, so return an
-                                     * empty array to demonstrate
-                                     * that fact.
-                                     */
-    virtual void get_unit_face_support_points (typename std::vector<Point<dim-1> > &) const;    
+    virtual unsigned int memory_consumption () const;    
   
   protected:
                                     /**
@@ -340,6 +300,24 @@ class FESystem : public FiniteElement<dim>
     typename std::vector<ElementPair> base_elements;
 
 
+                                    /**
+                                     * Initialize the
+                                     * @p{unit_support_points} field
+                                     * of the @ref{FiniteElementBase}
+                                     * class. Called from the
+                                     * constructor.
+                                     */
+    void initialize_unit_support_points ();
+
+                                    /**
+                                     * Initialize the
+                                     * @p{unit_face_support_points} field
+                                     * of the @ref{FiniteElementBase}
+                                     * class. Called from the
+                                     * constructor.
+                                     */
+    void initialize_unit_face_support_points ();
+
                                     /**
                                      * Helper function used in the constructor:
                                      * take a @p{FiniteElementData} object
index 01eda10290284aa4777871616527c175bf883a52..67abddfe975c697e143b616c53f7e2c964fd7943 100644 (file)
@@ -220,6 +220,56 @@ bool FiniteElementBase<dim>::operator == (const FiniteElementBase<dim> &f) const
 
 
 
+template <int dim>
+const std::vector<Point<dim> > &
+FiniteElementBase<dim>::get_unit_support_points () const
+{
+                                  // a finite element may define
+                                  // support points, but only if
+                                  // there are as many as there are
+                                  // degrees of freedom
+  Assert ((unit_support_points.size() == 0) ||
+         (unit_support_points.size() == dofs_per_cell),
+         ExcInternalError());
+  return unit_support_points;
+};
+
+
+
+template <int dim>
+bool
+FiniteElementBase<dim>::has_support_points () const
+{
+  return (unit_support_points.size() != 0);
+};
+
+
+
+template <int dim>
+const std::vector<Point<dim-1> > &
+FiniteElementBase<dim>::get_unit_face_support_points () const
+{
+                                  // a finite element may define
+                                  // support points, but only if
+                                  // there are as many as there are
+                                  // degrees of freedom on a face
+  Assert ((unit_face_support_points.size() == 0) ||
+         (unit_face_support_points.size() == dofs_per_face),
+         ExcInternalError());
+  return unit_face_support_points;
+};
+
+
+
+template <int dim>
+bool
+FiniteElementBase<dim>::has_face_support_points () const
+{
+  return (unit_face_support_points.size() != 0);
+};
+
+
+
 template <int dim>
 unsigned int
 FiniteElementBase<dim>::memory_consumption () const
@@ -333,26 +383,6 @@ FiniteElement<dim>::~FiniteElement ()
 
 
 
-template <int dim>
-void
-FiniteElement<dim>::get_unit_support_points (std::vector<Point<dim> > &points) const
-{
-                                  // default implementation
-  points.resize(0);
-}
-
-
-    
-template <int dim>
-void
-FiniteElement<dim>::get_unit_face_support_points (std::vector<Point<dim-1> > &points) const
-{
-                                  // default implementation
-  points.resize(0);
-}
-
-    
-
 template <int dim>
 Mapping<dim>::InternalDataBase*
 FiniteElement<dim>::get_face_data (const UpdateFlags       flags,
index e9d98c7a3ff529af00e1ae0fec1770132425e02b..b9da2f48ef2789f124822ec27adc8ac05815a5e9 100644 (file)
@@ -161,6 +161,46 @@ FE_DGQ<dim>::FE_DGQ (unsigned int degree)
                                     // matrix undefined, set size to zero
     for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
       restriction[i].reinit(0);
+
+  
+                                  // finally fill in support points
+  if (degree == 0)
+    {
+                                      // constant elements, take
+                                      // midpoint
+      unit_support_points.resize(1);
+      for (unsigned int i=0; i<dim; ++i)
+       unit_support_points[0](i) = 0.5;
+    }
+  else
+    {
+                                      // number of points: (degree+1)^dim
+      unsigned int n = degree+1;
+      for (unsigned int i=1; i<dim; ++i)
+       n *= degree+1;
+      
+      unit_support_points.resize(n);
+      
+      const double step = 1./degree;
+      Point<dim> p;
+      
+      unsigned int k=0;
+      for (unsigned int iz=0; iz <= ((dim>2) ? degree : 0) ; ++iz)
+       for (unsigned int iy=0; iy <= ((dim>1) ? degree : 0) ; ++iy)
+         for (unsigned int ix=0; ix<=degree; ++ix)
+           {
+             p(0) = ix * step;
+             if (dim>1)
+               p(1) = iy * step;
+             if (dim>2)
+               p(2) = iz * step;
+             
+             unit_support_points[k++] = p;
+           };
+    };
+
+                                  // note: no face support points for
+                                  // DG elements
 };
 
 
@@ -185,31 +225,6 @@ FE_DGQ<dim>::clone() const
 
 
 
-template <int dim>
-void
-FE_DGQ<dim>::get_unit_support_points (std::vector<Point<dim> > &points) const
-{
-  if (degree>0)
-    compute_support_points (points, degree);
-  else
-    {
-                                      // constant elements, take
-                                      // midpoint
-      points.resize(1);
-      for (unsigned int i=0; i<dim; ++i)
-       points[0](i)=0.5;
-    }
-}
-
-    
-template <int dim>
-void
-FE_DGQ<dim>::get_unit_face_support_points (std::vector<Point<dim-1> > &points) const
-{
-                                  // no face support points
-  points.resize(0);
-}
-
 
 //----------------------------------------------------------------------
 // Auxiliary functions
@@ -257,37 +272,6 @@ FE_DGQ<dim>::update_each (const UpdateFlags flags) const
 }
 
 
-template <int dim>
-void
-FE_DGQ<dim>::compute_support_points (std::vector<Point<dim> > &support_points,
-                                    const unsigned int degree)
-{
-  Assert(degree>0, ExcInternalError());
-                                  // number of points: (degree+1)^dim
-  unsigned int n = degree+1;
-  for (unsigned int i=1;i<dim;++i)
-    n *= degree+1;
-
-  support_points.resize(n);
-  
-  const double step = 1./degree;
-  Point<dim> p;
-  
-  unsigned int k=0;
-  for (unsigned int iz=0;iz <= ((dim>2) ? degree : 0) ; ++iz)
-    for (unsigned int iy=0;iy <= ((dim>1) ? degree : 0) ; ++iy)
-      for (unsigned int ix=0;ix<=degree;++ix)
-        {
-         p(0) = ix * step;
-         if (dim>1)
-           p(1) = iy * step;
-         if (dim>2)
-           p(2) = iz * step;
-         
-         support_points[k++] = p;
-       }
-}
-
 
 template <int dim>
 void
index 90b17638d29e817b196368e6a89e6fd6db58c69a..24a8ab12f8e78a82e8536ca147b41001a4deddd8 100644 (file)
@@ -570,7 +570,13 @@ FE_Q<dim>::FE_Q (const unsigned int degree)
       default:
            Assert (false,ExcNotImplemented());
     }
-}
+
+                                  // finally fill in support points
+                                  // on cell and face
+  initialize_unit_support_points ();
+  initialize_unit_face_support_points ();
+};
+
 
 
 template <int dim>
@@ -593,71 +599,82 @@ FE_Q<dim>::clone() const
 
 
 
-template <int dim>
-void
-FE_Q<dim>::get_unit_support_points (std::vector<Point<dim> > &points) const
-{
-  compute_support_points (degree, renumber, points);
-}
+//----------------------------------------------------------------------
+// Auxiliary functions
+//----------------------------------------------------------------------
+
 
 
-    
 template <int dim>
-void
-FE_Q<dim>::get_unit_face_support_points (std::vector<Point<dim-1> > &points) const
+void FE_Q<dim>::initialize_unit_support_points ()
 {
-  FE_Q<dim-1>::compute_support_points (degree, face_renumber, points);
-}
+                                  // number of points: (degree+1)^dim
+  unsigned int n = degree+1;
+  for (unsigned int i=1; i<dim; ++i)
+    n *= degree+1;
+  
+  unit_support_points.resize(n);
+  
+  const double step = 1./degree;
+  Point<dim> p;
+  
+  unsigned int k=0;
+  for (unsigned int iz=0; iz <= ((dim>2) ? degree : 0) ; ++iz)
+    for (unsigned int iy=0; iy <= ((dim>1) ? degree : 0) ; ++iy)
+      for (unsigned int ix=0; ix<=degree; ++ix)
+       {
+         p(0) = ix * step;
+         if (dim>1)
+           p(1) = iy * step;
+         if (dim>2)
+           p(2) = iz * step;
+         
+         unit_support_points[renumber[k++]] = p;
+       };
+};
 
 
 #if deal_II_dimension == 1
 
 template <>
-void
-FE_Q<1>::get_unit_face_support_points (std::vector<Point<0> > &points) const
+void FE_Q<1>::initialize_unit_face_support_points ()
 {
-  points.resize(0);
-}
+                                  // no faces in 1d, so nothing to do
+};
 
 #endif
 
 
-
-//----------------------------------------------------------------------
-// Auxiliary functions
-//----------------------------------------------------------------------
-
-
 template <int dim>
-void
-FE_Q<dim>::compute_support_points (const unsigned int               degree,
-                                  const std::vector<unsigned int> &renumber,
-                                  std::vector<Point<dim> >        &support_points)
+void FE_Q<dim>::initialize_unit_face_support_points ()
 {
-                                  // number of points: (degree+1)^dim
+  const unsigned int codim = dim-1;
+  
+                                  // number of points: (degree+1)^codim
   unsigned int n = degree+1;
-  for (unsigned int i=1; i<dim; ++i)
+  for (unsigned int i=1; i<codim; ++i)
     n *= degree+1;
-
-  support_points.resize(n);
+  
+  unit_face_support_points.resize(n);
   
   const double step = 1./degree;
-  Point<dim> p;
+  Point<codim> p;
   
   unsigned int k=0;
-  for (unsigned int iz=0; iz <= ((dim>2) ? degree : 0) ; ++iz)
-    for (unsigned int iy=0; iy <= ((dim>1) ? degree : 0) ; ++iy)
+  for (unsigned int iz=0; iz <= ((codim>2) ? degree : 0) ; ++iz)
+    for (unsigned int iy=0; iy <= ((codim>1) ? degree : 0) ; ++iy)
       for (unsigned int ix=0; ix<=degree; ++ix)
-        {
+       {
          p(0) = ix * step;
-         if (dim>1)
+         if (codim>1)
            p(1) = iy * step;
-         if (dim>2)
+         if (codim>2)
            p(2) = iz * step;
          
-         support_points[renumber[k++]] = p;
-       }
-}
+         unit_face_support_points[face_renumber[k++]] = p;
+       };
+};
+
 
 
 template <int dim>
index 140c00d0ab1724af368847f1741531544b42cbe0..e7f923b1451a70b84c4dcdcf89bf902757c7a8a4 100644 (file)
@@ -230,104 +230,6 @@ FESystem<dim>::clone() const
 
 
 
-template <int dim>
-void
-FESystem<dim>::
-get_unit_support_points (typename std::vector<Point<dim> > &unit_support_points) const
-{
-                                  // generate unit support points
-                                  // from unit support points of sub
-                                  // elements
-  unit_support_points.resize(dofs_per_cell);
-  
-  typename std::vector<Point<dim> > base_unit_support_points (base_element(0).dofs_per_cell);
-  unsigned int comp = 0;
-  for (unsigned int base_el=0; base_el<n_base_elements(); ++base_el)
-    {
-      const unsigned int
-       base_element_dofs_per_cell = base_element(base_el).dofs_per_cell;
-      base_element(base_el).get_unit_support_points (base_unit_support_points);
-      
-                                      // if one of the base elements
-                                      // has no support points, then
-                                      // it makes no sense to define
-                                      // support points for the
-                                      // composed element, so return
-                                      // an empty array to
-                                      // demonstrate that fact
-      if (base_unit_support_points.size()==0)
-       {
-         unit_support_points.resize(0);
-         return;
-       }
-
-                                      // otherwise distribute the
-                                      // support points of this base
-                                      // element to all degrees of
-                                      // freedom contributed by this
-                                      // base element
-      Assert(base_unit_support_points.size()==base_element_dofs_per_cell,
-            ExcInternalError());
-      for (unsigned int n=0; n<element_multiplicity(base_el); ++n, ++comp)
-       for (unsigned int i=0; i<base_element_dofs_per_cell; ++i)
-         unit_support_points[component_to_system_index(comp,i)]
-           = base_unit_support_points[i];
-    }
-}
-
-
-
-template <int dim>
-void
-FESystem<dim>::
-get_unit_face_support_points (typename std::vector<Point<dim-1> > &unit_support_points) const
-{
-                                  // generate unit face support points
-                                  // from unit support points of sub
-                                  // elements
-  unit_support_points.resize(dofs_per_face);
-  
-  typename std::vector<Point<dim-1> >
-    base_unit_support_points (base_element(0).dofs_per_cell);
-  
-  unsigned int comp = 0;
-  for (unsigned int base_el=0 ; base_el<n_base_elements(); ++base_el)
-    {
-      const unsigned int
-       base_element_dofs_per_face = base_element(base_el).dofs_per_face;
-      base_element(base_el).get_unit_face_support_points (base_unit_support_points);
-      
-                                      // if one of the base elements
-                                      // has no support points, then
-                                      // it makes no sense to define
-                                      // support points for the
-                                      // composed element, so return
-                                      // an empty array to
-                                      // demonstrate that fact
-      if (base_unit_support_points.size()==0)
-       {
-         unit_support_points.resize(0);
-         return;
-       }
-
-                                      // otherwise distribute the
-                                      // support points of this base
-                                      // element to all degrees of
-                                      // freedom contributed by this
-                                      // base element
-      Assert(base_unit_support_points.size()==base_element_dofs_per_face,
-            ExcNotImplemented());
-      for (unsigned int n=0; n<element_multiplicity(base_el); ++n, ++comp)
-       for (unsigned int i=0; i<base_element_dofs_per_face; ++i)
-         unit_support_points[face_component_to_system_index(comp,i)]
-           = base_unit_support_points[i];
-    }
-}
-
-
-
 template <int dim>
 UpdateFlags
 FESystem<dim>::update_once (const UpdateFlags flags) const
@@ -1181,6 +1083,11 @@ void FESystem<dim>::initialize ()
 //FESystem(FE_Q<dim> (3), 2) and for FESystem<dim>(FE_Q<dim> (1), 2,
 //FE_Q<dim> (3), 1) and for FESystem<dim>(FE_Q<dim> (4), 2))
   build_interface_constraints ();
+
+                                  // finally fill in support points
+                                  // on cell and face
+  initialize_unit_support_points ();
+  initialize_unit_face_support_points ();
 };
 
 
@@ -1201,6 +1108,144 @@ FESystem<dim>::multiply_dof_numbers (const FiniteElementData<dim> &fe_data,
 
 
 
+
+template <int dim>
+void
+FESystem<dim>::
+initialize_unit_support_points ()
+{
+                                      // if one of the base elements
+                                      // has no support points, then
+                                      // it makes no sense to define
+                                      // support points for the
+                                      // composed element, so return
+                                      // an empty array to
+                                      // demonstrate that
+                                      // fact
+  for (unsigned int base_el=0 ; base_el<n_base_elements(); ++base_el)
+    if (!base_element(base_el).has_support_points())
+      {
+       unit_support_points.resize(0);
+       return;
+      };
+
+                                  // generate unit support points
+                                  // from unit support points of sub
+                                  // elements
+  unit_support_points.resize(dofs_per_cell);
+  
+  unsigned int comp = 0;
+  for (unsigned int base_el=0; base_el<n_base_elements(); ++base_el)
+    {
+                                      // we know that there are
+                                      // support points on the cell,
+                                      // so collect them
+      const unsigned int
+       base_element_dofs_per_cell = base_element(base_el).dofs_per_cell;
+      const typename std::vector<Point<dim> >
+       & base_unit_support_points = base_element(base_el).get_unit_support_points ();
+      
+                                      // otherwise distribute the
+                                      // support points of this base
+                                      // element to all degrees of
+                                      // freedom contributed by this
+                                      // base element
+      Assert(base_unit_support_points.size()==base_element_dofs_per_cell,
+            ExcInternalError());
+      for (unsigned int n=0; n<element_multiplicity(base_el); ++n, ++comp)
+       for (unsigned int i=0; i<base_element_dofs_per_cell; ++i)
+         unit_support_points[component_to_system_index(comp,i)]
+           = base_unit_support_points[i];
+    }
+}
+
+
+#if deal_II_dimension == 1
+
+template <>
+void
+FESystem<1>::
+initialize_unit_face_support_points ()
+{
+                                  // no faces no work
+};
+
+#endif
+
+
+template <int dim>
+void
+FESystem<dim>::
+initialize_unit_face_support_points ()
+{
+                                      // if one of the base elements
+                                      // has no support points, then
+                                      // it makes no sense to define
+                                      // support points for the
+                                      // composed element, so return
+                                      // an empty array to
+                                      // demonstrate that fact (note
+                                      // that we ask whether the base
+                                      // element has no support
+                                      // points at all, not only none
+                                      // on the face!)
+  for (unsigned int base_el=0 ; base_el<n_base_elements(); ++base_el)
+    if (!base_element(base_el).has_support_points())
+      {
+       unit_face_support_points.resize(0);
+       return;
+      };
+  
+
+                                  // generate unit face support points
+                                  // from unit support points of sub
+                                  // elements
+  unit_face_support_points.resize(dofs_per_face);
+  
+  unsigned int comp = 0;
+  for (unsigned int base_el=0 ; base_el<n_base_elements(); ++base_el)
+    {
+                                      // in some cases, finite
+                                      // elements have support points
+                                      // (we have made sure that they
+                                      // have above) but don't have
+                                      // any on the face (e.g. DG
+                                      // elements). in that case,
+                                      // don't even bother with this
+                                      // base element and directly go
+                                      // to the next one:
+      if (!base_element(base_el).has_face_support_points())
+       {
+         comp += element_multiplicity(base_el);
+         continue;
+       };
+
+                                      // otherwise, we know that
+                                      // there are support points on
+                                      // the face, so collect them
+      const unsigned int
+       base_element_dofs_per_face = base_element(base_el).dofs_per_face;
+      const typename std::vector<Point<dim-1> > &
+       base_unit_support_points = base_element(base_el).get_unit_face_support_points ();
+      
+                                      // distribute the support
+                                      // points of this base element
+                                      // to all degrees of freedom
+                                      // contributed by this base
+                                      // element
+      Assert(base_unit_support_points.size()==base_element_dofs_per_face,
+            ExcNotImplemented());
+      for (unsigned int n=0; n<element_multiplicity(base_el); ++n, ++comp)
+       for (unsigned int i=0; i<base_element_dofs_per_face; ++i)
+         unit_face_support_points[face_component_to_system_index(comp,i)]
+           = base_unit_support_points[i];
+    }
+}
+
+
+
 template <int dim>
 FiniteElementData<dim>
 FESystem<dim>::multiply_dof_numbers (const FiniteElementData<dim> &fe1,
index 4d7ce886396380ccc929faa77530eaff5bf0e230..09b897d8a7beaeb0a57a757b3a46b34d995b85c2 100644 (file)
@@ -50,8 +50,8 @@ void FETools::get_interpolation_matrix(const FiniteElement<dim> &fe1,
                                   // the unit support points of the
                                   // fe2 element.
   std::vector<double> phantom_weights(fe2.dofs_per_cell,1.);
-  std::vector<Point<dim> > fe2_support_points (fe2.dofs_per_cell);
-  fe2.get_unit_support_points (fe2_support_points);
+  const typename std::vector<Point<dim> > &
+    fe2_support_points = fe2.get_unit_support_points ();
   Quadrature<dim> fe2_support_points_quadrature(fe2_support_points,
                                                phantom_weights);
 
index 69f15b5e67e2811a34a4c1ea7dd06974130f3955..98864bfb46b4527f863fd4dfebdd9720d5af4ada 100644 (file)
@@ -91,8 +91,8 @@ void VectorTools::interpolate (const DoFHandler<dim> &dof,
                                   // avoided to evaluate
                                   // the vectorfunction multiply at
                                   // the same point on a cell.
-  std::vector<Point<dim> > unit_support_points;
-  fe.get_unit_support_points(unit_support_points);
+  const typename std::vector<Point<dim> > &
+    unit_support_points = fe.get_unit_support_points();
   Assert (unit_support_points.size() != 0,
          ExcNonInterpolatingFE());
 
@@ -630,8 +630,9 @@ VectorTools::interpolate_boundary_values (const DoFHandler<dim>    &dof,
                                   // support points. this wil be used
                                   // to obtain the quadrature points
                                   // on the real cell's face
-  typename std::vector<Point<dim-1> > unit_support_points;
-  fe.get_unit_face_support_points(unit_support_points);
+  const typename std::vector<Point<dim-1> >
+    & unit_support_points = fe.get_unit_face_support_points();
+  
                                   // check whether there are support
                                   // points on the face, if not, then
                                   // this FE does not allow to be

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.