]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Major reorganisation to allow code sharing and make restriction of finite elements...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 30 Apr 1998 14:14:51 +0000 (14:14 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 30 Apr 1998 14:14:51 +0000 (14:14 +0000)
git-svn-id: https://svn.dealii.org/trunk@228 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_values.h
deal.II/deal.II/source/fe/fe_values.cc

index 5a2bba2448ac462ca568eaa9d9364d1f5d3f1cd5..97560206fd3842f5bff47bc70a2a6bf0725b309f 100644 (file)
@@ -19,155 +19,69 @@ template <int dim> class Quadrature;
 
 
 
-
-
 /**
-  Represent a finite element evaluated with a specific quadrature rule.
-  This class is an optimization which avoids evaluating the shape functions
-  at the quadrature points each time a quadrature takes place. Rather, the
-  values and gradients (and possibly higher order derivatives in future
-  versions of this library) are evaluated once and for all on the unit
-  cell before doing the quadrature itself. Only the Jacobian matrix of
-  the transformation from the unit cell to the real cell and the integration
-  points in real space are calculated each time we move on to a new cell.
-
-  The unit cell is defined to be the tensor product of the interval $[0,1]$
-  in the present number of dimensions. In part of the literature, the convention
-  is used that the unit cell be the tensor product of the interval $[-1,1]$,
-  which is to distinguished properly.
-
-  Objects of this class store a multitude of different values needed to
-  do the assemblage steps on real cells rather than on the unit cell. Among
-  these values are the values and gradients of the shape functions at the
-  quadrature points on the real and the unit cell, the location of the
-  quadrature points on the real and on the unit cell, the weights of the
-  quadrature points, the Jacobian matrices of the mapping from the unit to
-  the real cell at the quadrature points and so on.
-
-  The Jacobian matrix is defined to be
-  $$ J_{ij} = {d\xi_i \over dx_j} $$
-  where the $\xi_i$ are the coordinates on the unit cell and the $x_i$ are
-  the coordinates on the real cell.
-  This is the form needed to compute the gradient on the real cell from
-  the gradient on the unit cell. If we want to transform the area element
-  $dx dy$ from the real to the unit cell, we have to take the determinant of
-  the inverse matrix, which is the reciprocal value of the determinant of the
-  matrix defined above.
-
-  The #FEValues# object keeps track of those fields which really need to
-  be computed, since the computation of the gradients of the ansatz functions
-  on each real cell can be quite an expensive thing if it is not needed. The
-  object knows about which fields are needed by the #UpdateFlags# object
-  passed through the constructor. In debug mode, the accessor functions, which
-  return values from the different fields, check whether the required field
-  was initialized, thus avoiding use of unitialized data.
-
-
-  {\bf Member functions}
-
-  The functions of this class fall into different cathegories:
-  \begin{itemize}
-  \item #shape_value#, #shape_grad#, etc: return one of the values
-    of this object at a time. In many cases you will want to get
-    a whole bunch at a time for performance or convenience reasons,
-    then use the #get_*# functions.
-    
-  \item #get_shape_values#, #get_shape_grads#, etc: these return
-    a reference to a whole field. Usually these fields contain
-    the values of all ansatz functions at all quadrature points.
-
-  \item #get_function_values#, #get_function_gradients#: these
-    two functions offer a simple way to avoid the detour of the
-    ansatz functions, if you have a finite solution (resp. the
-    vector of values associated with the different ansatz functions.)
-    Then you may want to get information from the restriction of
-    the finite element function to a certain cell, e.g. the values
-    of the function at the quadrature points or the values of its
-    gradient. These two functions provide the information needed:
-    you pass it a vector holding the finite element solution and the
-    functions return the values or gradients of the finite element
-    function restricted to the cell which was given last time the
-    #reinit# function was given.
-    
-    Though possible in principle, these functions do not call the
-    #reinit# function, you have to do so yourself beforehand. On the
-    other hand, a copy of the cell iterator is stored which was used
-    last time the #reinit# function was called. This frees us from
-    the need to pass the cell iterator again to these two functions,
-    which guarantees that the cell used here is in sync with that used
-    for the #reinit# function. You should, however, make sure that
-    nothing substantial happens to the #DoFHandler# object or any
-    other involved instance between the #reinit# and the #get_function_*#
-    functions are called.
-
-  \item #reinit#: initialize the #FEValues# object for a certain cell.
-    See above for more information.
-  \end{itemize}
-  */
+   This class offers a multitude of arrays and other fields which are used by
+   the derived classes #FEValues# and #FEFaceValues#. In principle, it is the
+   back end of the front end for the unification of a certain finite element
+   and a quadrature formula which evaluates certain aspects of the finite
+   element at quadrature points.
+   
+   This class is an optimization which avoids evaluating the shape functions
+   at the quadrature points each time a quadrature takes place. Rather, the
+   values and gradients (and possibly higher order derivatives in future
+   versions of this library) are evaluated once and for all on the unit
+   cell or face before doing the quadrature itself. Only the Jacobian matrix of
+   the transformation from the unit cell or face to the real cell or face and
+   the integration points in real space are calculated each time we move on
+   to a new face.
+
+   Actually, this class does none of the evaluations at startup itself; this is
+   all done by the derived classes. It only offers the basic functionality,
+   like providing those fields that are common to the derived classes and
+   access to these fields. Any computations are in the derived classes. See there
+   for more information.
+
+   @author Wolfgang Bangerth, 1998
+*/
 template <int dim>
-class FEValues {
+class FEValuesBase {
   public:
-
-    
-    
+        
                                     /**
                                      * Number of quadrature points.
                                      */
     const unsigned int n_quadrature_points;
 
                                     /**
-                                     * Total number of shape functions.
+                                     * Total number of shape functions
+                                     * per cell. If we use this base class
+                                     * to evaluate a finite element on
+                                     * faces of cells, this is still the
+                                     * number of degrees of freedom per
+                                     * cell, not per face.
                                      */
     const unsigned int total_dofs;
-    
-                                    /**
-                                     * Constructor. Fill all arrays with the
-                                     * values of the shape functions of the
-                                     * specified finite element using the
-                                     * quadrature points of the given
-                                     * quadrature rule.
-                                     *
-                                     * This function actually only fills
-                                     * the fields related to the unit face,
-                                     * the fields related to a real face (like
-                                     * gradients, true quadrature points, etc.)
-                                     * need to be initialized using the
-                                     * #reinit# function.
-                                     */
-    FEValues (const FiniteElement<dim> &,
-             const Quadrature<dim> &,
-             const UpdateFlags);
-
-                                    /**
-                                     * Return the value of the #i#th shape
-                                     * function at the #j# quadrature point.
-                                     */
-    double shape_value (const unsigned int i,
-                       const unsigned int j) const;
 
                                     /**
-                                     * Return a pointer to the matrix holding
-                                     * all values of shape functions at all
-                                     * integration points, on the present cell.
-                                     * For the format of this matrix, see the
-                                     * documentation for the matrix itself.
-                                     */
-    const dFMatrix & get_shape_values () const;
+                                     * Constructor. Set up the array sizes
+                                     * with #n_q_points# quadrature points,
+                                     * #n_ansatz_points# ansatz points (on
+                                     * the cell or face), #n_dof# ansatz
+                                     * functions per cell and with the
+                                     * given pattern to update the fields
+                                     * when the #reinit# function of the
+                                     * derived classes is called. The
+                                     * fields themselves are not set up,
+                                     * this must happen in the derived
+                                     * class's constructor, only the sizes
+                                     * are set correctly.
+                                     */
+    FEValuesBase (const unsigned int n_q_points,
+                 const unsigned int n_ansatz_points,
+                 const unsigned int n_dofs,
+                 const UpdateFlags update_flags);
     
-                                    /**
-                                     * Return the values of the finite
-                                     * element function characterized
-                                     * by #fe_function# restricted to
-                                     * #cell# at the quadrature points.
-                                     *
-                                     * The function assumes that the
-                                     * #values# object already has the
-                                     * right size.
-                                     */
-    void get_function_values (const dVector      &fe_function,
-                             vector<double>     &values) const;
-
-                                    /**
+                                    /**
                                      * Return the gradient of the #i#th shape
                                      * function at the #j# quadrature point.
                                      * If you want to get the derivative in
@@ -208,6 +122,15 @@ class FEValues {
                                     /**
                                      * Return the position of the #i#th
                                      * quadrature point in real space.
+                                     *
+                                     * If this object is used to evaluate
+                                     * finite elements on faces of cells,
+                                     * and for curved boundary cells, using
+                                     * biquadratic or higher mappings
+                                     * of the unit cell to the real cell,
+                                     * these points may not be on the
+                                     * plane submannifold on which the
+                                     * vertices of the face lie.
                                      */
     const Point<dim> & quadrature_point (const unsigned int i) const;
 
@@ -232,6 +155,11 @@ class FEValues {
                                      * we have to take the continuous
                                      * function's value at the ansatz function
                                      * locations.
+                                     *
+                                     * For the evaluation of finite elements on
+                                     * faces of cells, #i# is the number
+                                     * of the ansatz function on the face, not
+                                     * on the cell.
                                      */
     const Point<dim> & ansatz_point (const unsigned int i) const;
 
@@ -246,6 +174,16 @@ class FEValues {
                                      * Return the Jacobi determinant times
                                      * the weight of the #i#th quadrature
                                      * point.
+                                     *
+                                     * If faces of cells are concerned,
+                                     * the jacobi determinant is that of the
+                                     * transformation of the unit face to
+                                     * the present face, not of the unit
+                                     * cell to the real cell (unlike for
+                                     * the #jacobi_matrix# array of the
+                                     * derived classes which store the cell
+                                     * transformation's Jacobi matrix in
+                                     * all cases).
                                      */
     double JxW (const unsigned int i) const;
 
@@ -255,28 +193,6 @@ class FEValues {
                                      * quadrature points.
                                      */
     const vector<double> & get_JxW_values () const;
-    
-                                    /**
-                                     * Reinitialize the gradients, Jacobi
-                                     * determinants, etc for the given cell
-                                     * and the given finite element.
-                                     *
-                                     * This function needs a boundary object
-                                     * passed, since this class needs to know
-                                     * how to handle faces which are located
-                                     * on the boundary of the domain. In that
-                                     * case, faces may be curved and the
-                                     * calculation of quadrature points,
-                                     * gradients and the like may need
-                                     * additional effort, depending on the
-                                     * mapping from the unit to the real cell
-                                     * (linear mappings use straight boundary
-                                     * segments, but higher order elements
-                                     * may use other ways.)
-                                     */
-    void reinit (const typename DoFHandler<dim>::cell_iterator &,
-                const FiniteElement<dim> &,
-                const Boundary<dim> &);
 
                                     /**
                                      * Exception
@@ -300,19 +216,17 @@ class FEValues {
                    int, int,
                    << "Vector has wrong size " << arg1
                    << ", expected size " << arg2);
-    
-  private:
                                     /**
-                                     * Store the values of the shape functions
-                                     * at the quadrature points. Rows in this
-                                     * matrix denote the values of a single
-                                     * shape function at the different points,
-                                     * columns are for a single point with the
-                                     * different shape functions.
+                                     * Exception
                                      */
-    dFMatrix             shape_values;
-
+    DeclException0 (ExcInternalError);
                                     /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcNotImplemented);
+    
+  protected:
+                                    /**
                                      * Store the gradients of the shape
                                      * functions at the quadrature points.
                                      * Since unfortunately the full matrix
@@ -327,20 +241,16 @@ class FEValues {
                                      */
     vector<vector<Point<dim> > >  shape_gradients;
 
-                                    /**
-                                     * Store the gradients of the shape
-                                     * functions at the quadrature points on
-                                     * the unit cell.
-                                     * This field is set up upon construction
-                                     * of the object and contains the gradients
-                                     * on the reference element.
-                                     */
-    vector<vector<Point<dim> > >   unit_shape_gradients;
-    
                                     /**
                                      * Store an array of the weights of the
                                      * quadrature points. This array is
                                      * set up upon construction.
+                                     *
+                                     * If faces rather than cells are
+                                     * considered, the weights are stored
+                                     * only once still, since they are
+                                     * not transformed and are thus the same
+                                     * for all faces.
                                      */
     vector<double>       weights;
 
@@ -365,15 +275,7 @@ class FEValues {
                                      */
     vector<Point<dim> >  quadrature_points;
 
-                                    /**
-                                     * Array of quadrature points in the unit
-                                     * cell. This array is set up upon
-                                     * construction and contains the quadrature
-                                     * points on the reference element.
-                                     */
-    vector<Point<dim> >  unit_quadrature_points;
-    
-                                    /**
+                                    /**
                                      * Array of points denoting the off-point
                                      * of the ansatz functions. In real space
                                      * (no-one seems to need the off-point
@@ -386,6 +288,14 @@ class FEValues {
                                      * Store the jacobi matrices at the
                                      * different quadrature points. This field
                                      * is set each time #reinit# is called.
+                                     *
+                                     * If faces rather than cells are considered
+                                     * this is the Jacobi matrix of the
+                                     * transformation of the unit cell to the
+                                     * real cell, not of the unit face to the
+                                     * face. We need this full matrix for the
+                                     * transformation of the gradients to the
+                                     * real cell.
                                      */
     vector<dFMatrix>     jacobi_matrices;
 
@@ -407,16 +317,199 @@ class FEValues {
 
 
 
+/**
+  Represent a finite element evaluated with a specific quadrature rule on
+  a cell.
+
+  The unit cell is defined to be the tensor product of the interval $[0,1]$
+  in the present number of dimensions. In part of the literature, the convention
+  is used that the unit cell be the tensor product of the interval $[-1,1]$,
+  which is to distinguished properly.
+
+  Objects of this class store a multitude of different values needed to
+  do the assemblage steps on real cells rather than on the unit cell. Among
+  these values are the values and gradients of the shape functions at the
+  quadrature points on the real and the unit cell, the location of the
+  quadrature points on the real and on the unit cell, the weights of the
+  quadrature points, the Jacobian matrices of the mapping from the unit to
+  the real cell at the quadrature points and so on.
+
+  The Jacobian matrix is defined to be
+  $$ J_{ij} = {d\xi_i \over dx_j} $$
+  where the $\xi_i$ are the coordinates on the unit cell and the $x_i$ are
+  the coordinates on the real cell.
+  This is the form needed to compute the gradient on the real cell from
+  the gradient on the unit cell. If we want to transform the area element
+  $dx dy$ from the real to the unit cell, we have to take the determinant of
+  the inverse matrix, which is the reciprocal value of the determinant of the
+  matrix defined above.
+
+  The #FEValues# object keeps track of those fields which really need to
+  be computed, since the computation of the gradients of the ansatz functions
+  on each real cell can be quite an expensive thing if it is not needed. The
+  object knows about which fields are needed by the #UpdateFlags# object
+  passed through the constructor. In debug mode, the accessor functions, which
+  return values from the different fields, check whether the required field
+  was initialized, thus avoiding use of unitialized data.
+
+
+  {\bf Member functions}
+
+  The functions of this class fall into different cathegories:
+  \begin{itemize}
+  \item #shape_value#, #shape_grad#, etc: return one of the values
+    of this object at a time. In many cases you will want to get
+    a whole bunch at a time for performance or convenience reasons,
+    then use the #get_*# functions.
+    
+  \item #get_shape_values#, #get_shape_grads#, etc: these return
+    a reference to a whole field. Usually these fields contain
+    the values of all ansatz functions at all quadrature points.
+
+  \item #get_function_values#, #get_function_gradients#: these
+    two functions offer a simple way to avoid the detour of the
+    ansatz functions, if you have a finite solution (resp. the
+    vector of values associated with the different ansatz functions.)
+    Then you may want to get information from the restriction of
+    the finite element function to a certain cell, e.g. the values
+    of the function at the quadrature points or the values of its
+    gradient. These two functions provide the information needed:
+    you pass it a vector holding the finite element solution and the
+    functions return the values or gradients of the finite element
+    function restricted to the cell which was given last time the
+    #reinit# function was given.
+    
+    Though possible in principle, these functions do not call the
+    #reinit# function, you have to do so yourself beforehand. On the
+    other hand, a copy of the cell iterator is stored which was used
+    last time the #reinit# function was called. This frees us from
+    the need to pass the cell iterator again to these two functions,
+    which guarantees that the cell used here is in sync with that used
+    for the #reinit# function. You should, however, make sure that
+    nothing substantial happens to the #DoFHandler# object or any
+    other involved instance between the #reinit# and the #get_function_*#
+    functions are called.
+
+  \item #reinit#: initialize the #FEValues# object for a certain cell.
+    See above for more information.
+  \end{itemize}
+
+   @author Wolfgang Bangerth, 1998  
+  */
+template <int dim>
+class FEValues : public FEValuesBase<dim> {
+  public:
+
+    
+    
+                                    /**
+                                     * Constructor. Fill all arrays with the
+                                     * values of the shape functions of the
+                                     * specified finite element using the
+                                     * quadrature points of the given
+                                     * quadrature rule.
+                                     *
+                                     * This function actually only fills
+                                     * the fields related to the unit face,
+                                     * the fields related to a real face (like
+                                     * gradients, true quadrature points, etc.)
+                                     * need to be initialized using the
+                                     * #reinit# function.
+                                     */
+    FEValues (const FiniteElement<dim> &,
+             const Quadrature<dim> &,
+             const UpdateFlags);
+
+                                    /**
+                                     * Return the value of the #i#th shape
+                                     * function at the #j# quadrature point.
+                                     */
+    double shape_value (const unsigned int i,
+                       const unsigned int j) const;
+
+                                    /**
+                                     * Return a pointer to the matrix holding
+                                     * all values of shape functions at all
+                                     * integration points, on the present cell.
+                                     * For the format of this matrix, see the
+                                     * documentation for the matrix itself.
+                                     */
+    const dFMatrix & get_shape_values () const;
+    
+                                    /**
+                                     * Return the values of the finite
+                                     * element function characterized
+                                     * by #fe_function# restricted to
+                                     * #cell# at the quadrature points.
+                                     *
+                                     * The function assumes that the
+                                     * #values# object already has the
+                                     * right size.
+                                     */
+    void get_function_values (const dVector      &fe_function,
+                             vector<double>     &values) const;
+
+    
+                                    /**
+                                     * Reinitialize the gradients, Jacobi
+                                     * determinants, etc for the given cell
+                                     * and the given finite element.
+                                     *
+                                     * This function needs a boundary object
+                                     * passed, since this class needs to know
+                                     * how to handle faces which are located
+                                     * on the boundary of the domain. In that
+                                     * case, faces may be curved and the
+                                     * calculation of quadrature points,
+                                     * gradients and the like may need
+                                     * additional effort, depending on the
+                                     * mapping from the unit to the real cell
+                                     * (linear mappings use straight boundary
+                                     * segments, but higher order elements
+                                     * may use other ways.)
+                                     */
+    void reinit (const typename DoFHandler<dim>::cell_iterator &,
+                const FiniteElement<dim> &,
+                const Boundary<dim> &);
+
+  private:
+                                    /**
+                                     * Store the values of the shape functions
+                                     * at the quadrature points. Rows in this
+                                     * matrix denote the values of a single
+                                     * shape function at the different points,
+                                     * columns are for a single point with the
+                                     * different shape functions.
+                                     */
+    dFMatrix             shape_values;
+
+                                    /**
+                                     * Store the gradients of the shape
+                                     * functions at the quadrature points on
+                                     * the unit cell.
+                                     * This field is set up upon construction
+                                     * of the object and contains the gradients
+                                     * on the reference element.
+                                     */
+    vector<vector<Point<dim> > >   unit_shape_gradients;
+    
+
+                                    /**
+                                     * Array of quadrature points in the unit
+                                     * cell. This array is set up upon
+                                     * construction and contains the quadrature
+                                     * points on the reference element.
+                                     */
+    vector<Point<dim> >  unit_quadrature_points;
+    
+};
+
+
+
 
 /**
-  Represent a finite element evaluated with a specific quadrature rule.
-  This class is an optimization which avoids evaluating the shape functions
-  at the quadrature points each time a quadrature takes place. Rather, the
-  values and gradients (and possibly higher order derivatives in future
-  versions of this library) are evaluated once and for all on the unit
-  face before doing the quadrature itself. Only the Jacobian matrix of
-  the transformation from the unit face to the real face and the integration
-  points in real space are calculated each time we move on to a new face.
+  Represent a finite element evaluated with a specific quadrature rule on
+  the face of a cell.
 
   The unit face is defined to be the tensor product of the interval $[0,1]$
   in the present number of dimensions minus one. In part of the literature,
@@ -473,33 +566,15 @@ class FEValues {
 
   Finally, we will often need the outward normal to a cell at the quadrature
   points. While this could in principle be easily done using the Jacobi
-  matrices at the quadrature points and the normal vectors to the unit cell
-  (also easily derived, since they have an appealingly easy form for the unit
-  cell ;-), it is more efficiently done by the finite element class itself.
-  For example for (bi-, tri-)linear mappings the normal vector is readily
-  available without compicated matrix-vector-multiplications.
-  */
-template <int dim>
-class FEFaceValues {
-  public:
-
-    
-    
-                                    /**
-                                     * Number of quadrature points on
-                                     * the face.
-                                     */
-    const unsigned int n_quadrature_points;
-
-                                    /**
-                                     * Total number of shape functions
-                                     * on the cell adjacent to this face.
-                                     * This number is not the same as the
-                                     * number of shape functions of which
-                                     * the center is located on the face.
-                                     */
-    const unsigned int total_dofs;
-    
+  matrices at the quadrature points and the normal vectors to the unit cell
+  (also easily derived, since they have an appealingly easy form for the unit
+  cell ;-), it is more efficiently done by the finite element class itself.
+  For example for (bi-, tri-)linear mappings the normal vector is readily
+  available without compicated matrix-vector-multiplications.
+  */
+template <int dim>
+class FEFaceValues : public FEValuesBase<dim> {
+  public:
                                     /**
                                      * Constructor. Fill all arrays with the
                                      * values of the shape functions of the
@@ -549,109 +624,6 @@ class FEFaceValues {
     void get_function_values (const dVector      &fe_function,
                              vector<double>     &values) const;
 
-                                    /**
-                                     * Return the gradient of the #i#th shape
-                                     * function at the #j# quadrature point.
-                                     * If you want to get the derivative in
-                                     * one of the coordinate directions, use
-                                     * the appropriate function of the #Point#
-                                     * class to extract one component. Since
-                                     * only a reference to the gradient's value
-                                     * is returned, there should be no major
-                                     * performance drawback.
-                                     * The function returns the gradient on the
-                                     * real element, not the reference element.
-                                     */
-    const Point<dim> & shape_grad (const unsigned int i,
-                                  const unsigned int j) const;
-
-                                    /** 
-                                     * Return a pointer to the matrix holding
-                                     * all gradients of shape functions at all
-                                     * integration points, on the present cell.
-                                     * For the format of this matrix, see the
-                                     * documentation for the matrix itself.
-                                     */
-    const vector<vector<Point<dim> > > & get_shape_grads () const;
-    
-                                    /**
-                                     * Return the gradients of the finite
-                                     * element function characterized
-                                     * by #fe_function# restricted to
-                                     * #cell# at the quadrature points.
-                                     *
-                                     * The function assumes that the
-                                     * #gradients# object already has the
-                                     * right size.
-                                     */
-    void get_function_grads (const dVector       &fe_function,
-                            vector<Point<dim> > &gradients) const;
-
-                                    /**
-                                     * Return the position of the #i#th
-                                     * quadrature point in real space.
-                                     *
-                                     * For curved boundary cells, using
-                                     * biquadratic or higher mappings
-                                     * of the unit cell to the real cell,
-                                     * these points may not be on the
-                                     * plane submannifold on which the
-                                     * vertices of the face lie.
-                                     */
-    const Point<dim> & quadrature_point (const unsigned int i) const;
-
-                                    /**
-                                     * Return a pointer to the vector of
-                                     * quadrature points.
-                                     */
-    const vector<Point<dim> > & get_quadrature_points () const;
-
-                                    /**
-                                     * Return the point in real space where
-                                     * the #i#th ansatz function is located
-                                     * (location is in the sense of where it
-                                     * assumes its nominal properties, e.g. at
-                                     * the vertex of a cell, at the center of
-                                     * a line, etc).
-                                     *
-                                     * This function is needed for the
-                                     * interpolation problem: if we want to
-                                     * transfer a continuous function to a
-                                     * finite element function by interpolation
-                                     * we have to take the continuous
-                                     * function's value at the ansatz function
-                                     * locations.
-                                     */
-    const Point<dim> & ansatz_point (const unsigned int i) const;
-
-                                    /**
-                                     * Return a pointer to the vector of points
-                                     * denoting the location of the ansatz
-                                     * functions.
-                                     */
-    const vector<Point<dim> > & get_ansatz_points () const;
-    
-                                    /**
-                                     * Return the Jacobi determinant times
-                                     * the weight of the #i#th quadrature
-                                     * point. The Jacobi determinant is that
-                                     * of the transformation of the unit
-                                     * face to the real face, not of the
-                                     * alike cells.
-                                     */
-    double JxW (const unsigned int i) const;
-
-                                    /**
-                                     * Return a pointer to the array holding
-                                     * the JxW values at the different
-                                     * quadrature points. The Jacobi
-                                     * determinant is that
-                                     * of the transformation of the unit
-                                     * face to the real face, not of the
-                                     * alike cells.
-                                     */
-    const vector<double> & get_JxW_values () const;
-
                                     /**
                                      * Return the outward normal vector to
                                      * the cell at the #i#th quadrature
@@ -691,36 +663,6 @@ class FEFaceValues {
                 const FiniteElement<dim>             &fe,
                 const Boundary<dim>                  &boundary);
 
-                                    /**
-                                     * Exception
-                                     */
-    DeclException2 (ExcInvalidIndex,
-                   int, int,
-                   << "The index " << arg1
-                   << " is out of range, it should be less than " << arg2);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcAccessToUninitializedField);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcCannotInitializeField);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcInternalError);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcNotImplemented);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException2 (ExcWrongVectorSize,
-                   int, int,
-                   << "Vector has wrong size " << arg1
-                   << ", expected size " << arg2);
   private:
                                     /**
                                      * Store the values of the shape functions
@@ -734,25 +676,6 @@ class FEFaceValues {
                                      */
     dFMatrix             shape_values[2*dim];
 
-                                    /**
-                                     * Store the gradients of the shape
-                                     * functions at the quadrature points.
-                                     * Since unfortunately the full matrix
-                                     * classes of DEAL are not templated,
-                                     * we have to store them in an
-                                     * archetypic style.
-                                     *
-                                     * This field is reset each time
-                                     * #reinit# is called and contains the
-                                     * gradients on the real element, rather
-                                     * than on the reference element. This
-                                     * function does the transformation from
-                                     * the unit cell to the real cell using
-                                     * the #unit_shape_gradients# for the
-                                     * selected face.
-                                     */
-    vector<vector<Point<dim> > >  shape_gradients;
-
                                     /**
                                      * Store the gradients of the shape
                                      * functions at the quadrature points on
@@ -765,36 +688,6 @@ class FEFaceValues {
                                      */
     vector<vector<Point<dim> > >   unit_shape_gradients[2*dim];
     
-                                    /**
-                                     * Store an array of the weights of the
-                                     * quadrature points. This array is
-                                     * set up upon construction.
-                                     *
-                                     * Since these weights are not transformed
-                                     * they are the same for all faces.
-                                     */
-    vector<double>       weights;
-
-                                    /**
-                                     * Store an array of weights times the
-                                     * Jacobi determinant at the quadrature
-                                     * points. This function is reset each time
-                                     * #reinit# is called. The Jacobi determinant
-                                     * is actually the reciprocal value of the
-                                     * Jacobi matrices stored in this class,
-                                     * see the general documentation of this
-                                     * class for more information.
-                                     */
-    vector<double>       JxW_values;
-
-                                    /**
-                                     * Array of quadrature points. This array
-                                     * is set up upon calling #reinit# and
-                                     * contains the quadrature points on the
-                                     * real element, rather than on the
-                                     * reference element.
-                                     */
-    vector<Point<dim> >  quadrature_points;
 
                                     /**
                                      * Array of quadrature points on the
@@ -817,32 +710,10 @@ class FEFaceValues {
                                      */
     vector<Point<dim> >  global_unit_quadrature_points[2*dim];
     
-                                    /**
-                                     * Array of points denoting the off-point
-                                     * of the ansatz functions. In real space
-                                     * (no-one seems to need the off-point
-                                     * on the unit cell, so no function is
-                                     * provided for this).
-                                     */
-    vector<Point<dim> >  ansatz_points;
-    
-                                    /**
-                                     * Store the jacobi matrices at the
-                                     * different quadrature points. This field
-                                     * is set each time #reinit# is called.
-                                     * This is the Jacobi matrix of the
-                                     * transformation of the unit cell to the
-                                     * real cell, not of the unit face to the
-                                     * face. We need this full matrix for the
-                                     * transformation of the gradients to the
-                                     * real cell.
-                                     */
-    vector<dFMatrix>     jacobi_matrices;
-
                                     /**
                                      * List of values denoting the determinant
                                      * of the transformation from the unit face
-                                     * to the real face. Neede to actually
+                                     * to the real face. Needed to actually
                                      * compute the JxW values.
                                      */
     vector<double>       face_jacobi_determinants;
@@ -854,21 +725,6 @@ class FEFaceValues {
                                      */
     vector<Point<dim> >  normal_vectors;
     
-                                    /**
-                                     * Store which fields are to be updated by
-                                     * the reinit function.
-                                     */
-    UpdateFlags          update_flags;
-
-                                    /**
-                                     * Store the cell selected last time
-                                     * the #reinit# function was called
-                                     * to make access
-                                     * to the #get_function_*# functions
-                                     * safer.
-                                     */
-    DoFHandler<dim>::cell_iterator present_cell;
-
                                     /**
                                      * Store the number of the face selected
                                      * last time the #reinit# function was
@@ -881,15 +737,7 @@ class FEFaceValues {
 
 
 
-/*------------------------ Inline functions: FEValues ----------------------------*/
-
-
-
-template <int dim>
-inline
-const dFMatrix & FEValues<dim>::get_shape_values () const {
-  return shape_values;
-};
+/*------------------------ Inline functions: FEValuesBase ------------------------*/
 
 
 
@@ -897,7 +745,7 @@ const dFMatrix & FEValues<dim>::get_shape_values () const {
 template <int dim>
 inline
 const vector<vector<Point<dim> > > &
-FEValues<dim>::get_shape_grads () const {
+FEValuesBase<dim>::get_shape_grads () const {
   Assert (update_flags & update_gradients, ExcAccessToUninitializedField());
   return shape_gradients;
 };
@@ -907,7 +755,7 @@ FEValues<dim>::get_shape_grads () const {
 template <int dim>
 inline
 const vector<Point<dim> > &
-FEValues<dim>::get_quadrature_points () const {
+FEValuesBase<dim>::get_quadrature_points () const {
   Assert (update_flags & update_q_points, ExcAccessToUninitializedField());
   return quadrature_points;
 };
@@ -917,7 +765,7 @@ FEValues<dim>::get_quadrature_points () const {
 template <int dim>
 inline
 const vector<Point<dim> > &
-FEValues<dim>::get_ansatz_points () const {
+FEValuesBase<dim>::get_ansatz_points () const {
   Assert (update_flags & update_ansatz_points, ExcAccessToUninitializedField());
   return ansatz_points;
 };
@@ -927,7 +775,7 @@ FEValues<dim>::get_ansatz_points () const {
 template <int dim>
 inline
 const vector<double> &
-FEValues<dim>::get_JxW_values () const {
+FEValuesBase<dim>::get_JxW_values () const {
   Assert (update_flags & update_JxW_values, ExcAccessToUninitializedField());
   return JxW_values;
 };
@@ -935,55 +783,24 @@ FEValues<dim>::get_JxW_values () const {
 
 
 
-
-/*------------------------ Inline functions: FEFaceValues ------------------------*/
-
-
-template <int dim>
-inline
-const dFMatrix & FEFaceValues<dim>::get_shape_values () const {
-  return shape_values[selected_face];
-};
-
-
-
-
-template <int dim>
-inline
-const vector<vector<Point<dim> > > &
-FEFaceValues<dim>::get_shape_grads () const {
-  Assert (update_flags & update_gradients, ExcAccessToUninitializedField());
-  return shape_gradients;
-};
-
+/*------------------------ Inline functions: FEValues ----------------------------*/
 
 
 template <int dim>
 inline
-const vector<Point<dim> > &
-FEFaceValues<dim>::get_quadrature_points () const {
-  Assert (update_flags & update_q_points, ExcAccessToUninitializedField());
-  return quadrature_points;
+const dFMatrix & FEValues<dim>::get_shape_values () const {
+  return shape_values;
 };
 
 
 
-template <int dim>
-inline
-const vector<Point<dim> > &
-FEFaceValues<dim>::get_ansatz_points () const {
-  Assert (update_flags & update_ansatz_points, ExcAccessToUninitializedField());
-  return ansatz_points;
-};
-
+/*------------------------ Inline functions: FEFaceValues ------------------------*/
 
 
 template <int dim>
 inline
-const vector<double> &
-FEFaceValues<dim>::get_JxW_values () const {
-  Assert (update_flags & update_JxW_values, ExcAccessToUninitializedField());
-  return JxW_values;
+const dFMatrix & FEFaceValues<dim>::get_shape_values () const {
+  return shape_values[selected_face];
 };
 
 
index a51d5e85d8c697314853a0190e2da81c6a49335c..8cd9fe4f089c8aeffb409fd3b2749237dc25b7c2 100644 (file)
@@ -6,89 +6,39 @@
 #include <grid/tria_iterator.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_boundary.h>
+#include <grid/dof_accessor.h>
 
 
 
-/*------------------------------- FEValues -------------------------------*/
-
-
-template <int dim>
-FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
-                        const Quadrature<dim>    &quadrature,
-                        const UpdateFlags         update_flags) :
-               n_quadrature_points(quadrature.n_quadrature_points),
-               total_dofs(fe.total_dofs),
-               shape_values(fe.total_dofs, quadrature.n_quadrature_points),
-               shape_gradients(fe.total_dofs,
-                               vector<Point<dim> >(quadrature.n_quadrature_points)),
-               unit_shape_gradients(fe.total_dofs,
-                                    vector<Point<dim> >(quadrature.n_quadrature_points)),
-               weights(quadrature.n_quadrature_points, 0),
-               JxW_values(quadrature.n_quadrature_points, 0),
-               quadrature_points(quadrature.n_quadrature_points, Point<dim>()),
-               unit_quadrature_points(quadrature.get_quad_points()),
-               ansatz_points (fe.total_dofs, Point<dim>()),
-               jacobi_matrices (quadrature.n_quadrature_points,
-                                dFMatrix(dim,dim)),
-               update_flags (update_flags)
-{
-  for (unsigned int i=0; i<fe.total_dofs; ++i)
-    for (unsigned int j=0; j<n_quadrature_points; ++j) 
-      {
-       shape_values(i,j) = fe.shape_value(i, quadrature.quad_point(j));
-       unit_shape_gradients[i][j]
-         = fe.shape_grad(i, quadrature.quad_point(j));
-      };
-
-  for (unsigned int i=0; i<n_quadrature_points; ++i) 
-    {
-      weights[i] = quadrature.weight(i);
-    };
-};
-
-
-
-template <int dim>
-double FEValues<dim>::shape_value (const unsigned int i,
-                                  const unsigned int j) const {
-  Assert (i<shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
-  Assert (j<shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
-
-  return shape_values(i,j);
-};
-
+/*------------------------------- FEValuesBase ---------------------------*/
 
 
 template <int dim>
-void FEValues<dim>::get_function_values (const dVector  &fe_function,
-                                        vector<double> &values) const {
-  Assert (values.size() == n_quadrature_points,
-         ExcWrongVectorSize(values.size(), n_quadrature_points));
-
-                                  // get function values of dofs
-                                  // on this cell
-  vector<double> dof_values (total_dofs, 0);
-  present_cell->get_dof_values (fe_function, dof_values);
-
-                                  // initialize with zero
-  fill_n (values.begin(), n_quadrature_points, 0);
+FEValuesBase<dim>::FEValuesBase (const unsigned int n_q_points,
+                                const unsigned int n_ansatz_points,
+                                const unsigned int n_dofs,
+                                const UpdateFlags update_flags) :
+               n_quadrature_points (n_q_points),
+               total_dofs (n_dofs),
+               shape_gradients (n_dofs, vector<Point<dim> >(n_q_points)),
+               weights (n_q_points, 0),
+               JxW_values (n_q_points, 0),
+               quadrature_points (n_q_points, Point<dim>()),
+               ansatz_points (n_ansatz_points, Point<dim>()),
+               jacobi_matrices (n_q_points, dFMatrix(dim,dim)),
+               update_flags (update_flags) {};
 
-                                  // add up contributions of ansatz
-                                  // functions
-  for (unsigned int point=0; point<n_quadrature_points; ++point)
-    for (unsigned int shape_func=0; shape_func<total_dofs; ++shape_func)
-      values[point] += (dof_values[shape_func] *
-                       shape_values(shape_func, point));
-};
 
 
 
 template <int dim>
 const Point<dim> &
-FEValues<dim>::shape_grad (const unsigned int i,
-                          const unsigned int j) const {
-  Assert (i<shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
-  Assert (j<shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
+FEValuesBase<dim>::shape_grad (const unsigned int i,
+                              const unsigned int j) const {
+  Assert (i<shape_gradients.size(),
+         ExcInvalidIndex (i, shape_gradients.size()));
+  Assert (j<shape_gradients[i].size(),
+         ExcInvalidIndex (j, shape_gradients[i].size()));
   Assert (update_flags & update_gradients, ExcAccessToUninitializedField());
 
   return shape_gradients[i][j];
@@ -97,8 +47,8 @@ FEValues<dim>::shape_grad (const unsigned int i,
 
 
 template <int dim>
-void FEValues<dim>::get_function_grads (const dVector       &fe_function,
-                                       vector<Point<dim> > &gradients) const {
+void FEValuesBase<dim>::get_function_grads (const dVector       &fe_function,
+                                           vector<Point<dim> > &gradients) const {
   Assert (gradients.size() == n_quadrature_points,
          ExcWrongVectorSize(gradients.size(), n_quadrature_points));
 
@@ -121,7 +71,7 @@ void FEValues<dim>::get_function_grads (const dVector       &fe_function,
 
 
 template <int dim>
-const Point<dim> & FEValues<dim>::quadrature_point (const unsigned int i) const {
+const Point<dim> & FEValuesBase<dim>::quadrature_point (const unsigned int i) const {
   Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
   Assert (update_flags & update_q_points, ExcAccessToUninitializedField());
   
@@ -131,7 +81,7 @@ const Point<dim> & FEValues<dim>::quadrature_point (const unsigned int i) const
 
 
 template <int dim>
-const Point<dim> & FEValues<dim>::ansatz_point (const unsigned int i) const {
+const Point<dim> & FEValuesBase<dim>::ansatz_point (const unsigned int i) const {
   Assert (i<ansatz_points.size(), ExcInvalidIndex(i, ansatz_points.size()));
   Assert (update_flags & update_ansatz_points, ExcAccessToUninitializedField());
   
@@ -141,7 +91,7 @@ const Point<dim> & FEValues<dim>::ansatz_point (const unsigned int i) const {
 
 
 template <int dim>
-double FEValues<dim>::JxW (const unsigned int i) const {
+double FEValuesBase<dim>::JxW (const unsigned int i) const {
   Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
   Assert (update_flags & update_JxW_values, ExcAccessToUninitializedField());
   
@@ -150,6 +100,75 @@ double FEValues<dim>::JxW (const unsigned int i) const {
 
 
 
+
+/*------------------------------- FEValues -------------------------------*/
+
+template <int dim>
+FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
+                        const Quadrature<dim>    &quadrature,
+                        const UpdateFlags         update_flags) :
+               FEValuesBase<dim> (quadrature.n_quadrature_points,
+                                  fe.total_dofs,
+                                  fe.total_dofs,
+                                  update_flags),
+               shape_values(fe.total_dofs, quadrature.n_quadrature_points),
+               unit_shape_gradients(fe.total_dofs,
+                                    vector<Point<dim> >(quadrature.n_quadrature_points)),
+               unit_quadrature_points(quadrature.get_quad_points())
+{
+  for (unsigned int i=0; i<fe.total_dofs; ++i)
+    for (unsigned int j=0; j<n_quadrature_points; ++j) 
+      {
+       shape_values(i,j) = fe.shape_value(i, quadrature.quad_point(j));
+       unit_shape_gradients[i][j]
+         = fe.shape_grad(i, quadrature.quad_point(j));
+      };
+
+  for (unsigned int i=0; i<n_quadrature_points; ++i) 
+    {
+      weights[i] = quadrature.weight(i);
+    };
+};
+
+
+
+template <int dim>
+double FEValues<dim>::shape_value (const unsigned int i,
+                                  const unsigned int j) const {
+  Assert (i<shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
+  Assert (j<shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
+
+  return shape_values(i,j);
+};
+
+
+
+template <int dim>
+void FEValues<dim>::get_function_values (const dVector  &fe_function,
+                                        vector<double> &values) const {
+  Assert (values.size() == n_quadrature_points,
+         ExcWrongVectorSize(values.size(), n_quadrature_points));
+
+                                  // get function values of dofs
+                                  // on this cell
+  vector<double> dof_values (total_dofs, 0);
+  present_cell->get_dof_values (fe_function, dof_values);
+
+                                  // initialize with zero
+  fill_n (values.begin(), n_quadrature_points, 0);
+
+                                  // add up contributions of ansatz
+                                  // functions
+  for (unsigned int point=0; point<n_quadrature_points; ++point)
+    for (unsigned int shape_func=0; shape_func<total_dofs; ++shape_func)
+      values[point] += (dof_values[shape_func] *
+                       shape_values(shape_func, point));
+};
+
+
+
+
+
 template <int dim>
 void FEValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &cell,
                            const FiniteElement<dim>                      &fe,
@@ -217,19 +236,13 @@ template <int dim>
 FEFaceValues<dim>::FEFaceValues (const FiniteElement<dim> &fe,
                                 const Quadrature<dim-1>  &quadrature,
                                 const UpdateFlags         update_flags) :
-               n_quadrature_points(quadrature.n_quadrature_points),
-               total_dofs(fe.total_dofs),
-               shape_gradients(fe.total_dofs,
-                               vector<Point<dim> >(quadrature.n_quadrature_points)),
-               weights(quadrature.n_quadrature_points, 0),
-               JxW_values(quadrature.n_quadrature_points, 0),
-               quadrature_points(quadrature.n_quadrature_points, Point<dim>()),
+               FEValuesBase<dim> (quadrature.n_quadrature_points,
+                                  fe.dofs_per_face,
+                                  fe.total_dofs,
+                                  update_flags),
                unit_quadrature_points(quadrature.get_quad_points()),
-               ansatz_points (fe.dofs_per_face, Point<dim>()),
-               jacobi_matrices (quadrature.n_quadrature_points,dFMatrix(dim,dim)),
                face_jacobi_determinants (quadrature.n_quadrature_points,0),
                normal_vectors (quadrature.n_quadrature_points,Point<dim>()),
-               update_flags (update_flags),
                selected_face(0)
 {
   for (unsigned int face=0; face<2*dim; ++face)
@@ -337,68 +350,6 @@ void FEFaceValues<dim>::get_function_values (const dVector  &fe_function,
 
 
 
-template <int dim>
-const Point<dim> &
-FEFaceValues<dim>::shape_grad (const unsigned int i,
-                              const unsigned int j) const {
-  Assert (i<shape_values[selected_face].m(),
-         ExcInvalidIndex (i, shape_values[selected_face].m()));
-  Assert (j<shape_values[selected_face].n(),
-         ExcInvalidIndex (j, shape_values[selected_face].n()));
-  Assert (update_flags & update_gradients,
-         ExcAccessToUninitializedField());
-
-  return shape_gradients[i][j];
-};
-
-
-
-template <int dim>
-void FEFaceValues<dim>::get_function_grads (const dVector       &fe_function,
-                                           vector<Point<dim> > &gradients) const {
-  Assert (gradients.size() == n_quadrature_points,
-         ExcWrongVectorSize(gradients.size(), n_quadrature_points));
-
-                                  // get function values of dofs
-                                  // on this cell
-  vector<double> dof_values (total_dofs, 0);
-  present_cell->get_dof_values (fe_function, dof_values);
-
-                                  // initialize with zero
-  fill_n (gradients.begin(), n_quadrature_points, Point<dim>());
-
-                                  // add up contributions of ansatz
-                                  // functions
-  for (unsigned int point=0; point<n_quadrature_points; ++point)
-    for (unsigned int shape_func=0; shape_func<total_dofs; ++shape_func)
-      gradients[point] += (dof_values[shape_func] *
-                          shape_gradients[shape_func][point]);
-};
-
-
-
-template <int dim>
-const Point<dim> & FEFaceValues<dim>::quadrature_point (const unsigned int i) const {
-  Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
-  Assert (update_flags & update_q_points,
-         ExcAccessToUninitializedField());
-  
-  return quadrature_points[i];
-};
-
-
-
-template <int dim>
-const Point<dim> & FEFaceValues<dim>::ansatz_point (const unsigned int i) const {
-  Assert (i<ansatz_points.size(), ExcInvalidIndex(i, ansatz_points.size()));
-  Assert (update_flags & update_ansatz_points,
-         ExcAccessToUninitializedField());
-  
-  return ansatz_points[i];
-};
-
-
-
 template <int dim>
 const Point<dim> & FEFaceValues<dim>::normal_vector (const unsigned int i) const {
   Assert (i<normal_vectors.size(), ExcInvalidIndex(i, normal_vectors.size()));
@@ -410,18 +361,6 @@ const Point<dim> & FEFaceValues<dim>::normal_vector (const unsigned int i) const
 
 
 
-template <int dim>
-double FEFaceValues<dim>::JxW (const unsigned int i) const {
-  Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
-  Assert (update_flags & update_JxW_values,
-         ExcAccessToUninitializedField());
-  
-  return JxW_values[i];
-};
-
-
-#include <grid/dof_accessor.h>
-
 template <int dim>
 void FEFaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &cell,
                                const unsigned int                             face_no,
@@ -494,6 +433,9 @@ void FEFaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &c
 
 /*------------------------------- Explicit Instantiations -------------*/
 
+template class FEValuesBase<1>;
+template class FEValuesBase<2>;
+
 template class FEValues<1>;
 template class FEValues<2>;
 

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.