]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Next step towards restriction of finite elements to subfaces.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 4 May 1998 08:39:09 +0000 (08:39 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 4 May 1998 08:39:09 +0000 (08:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@239 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 97560206fd3842f5bff47bc70a2a6bf0725b309f..81eb1c27386683d5948ec02b9c1becd3b9ec86de 100644 (file)
@@ -20,7 +20,7 @@ template <int dim> class Quadrature;
 
 
 /**
-   This class offers a multitude of arrays and other fields which are used by
+  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
@@ -41,6 +41,65 @@ template <int dim> class Quadrature;
    access to these fields. Any computations are in the derived classes. See there
    for more information.
 
+   It has support for the restriction of finite elements to faces of cells or
+   even to subfaces (i.e. refined faces). For this purpose, it offers an array
+   of matrices of ansatz function values, rather than one. Since the value of
+   a function at a quadrature point is an invariant under the transformation
+   from the unit cell to the real cell, it is only evaluated once upon startup.
+   However, when considering the restriction of a finite element to a face of
+   a cell (using a given quadrature rule), we may be tempted to compute the
+   restriction to all faces at startup (thus ending in four array of ansatz
+   function values in two dimensions, one per face, and even more in higher
+   dimensions) and let the respective #reinit# function of the derived classes
+   set a number which of the fields is to be taken when the user requests the
+   function values. This is done through the #selected_dataset# variable. See
+   the derived classes and the #get_values# function for the exact usage of
+   this variable.
+
+   
+   {\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.
+     This function is not in the present class but only in the derived
+     classes and has a variable call syntax. 
+     See the docs for the derived classes for more information.
+  \end{itemize}
+
    @author Wolfgang Bangerth, 1998
 */
 template <int dim>
@@ -79,8 +138,48 @@ class FEValuesBase {
     FEValuesBase (const unsigned int n_q_points,
                  const unsigned int n_ansatz_points,
                  const unsigned int n_dofs,
+                 const unsigned int n_values_array,
                  const UpdateFlags update_flags);
     
+
+                                    /**
+                                     * Return the value of the #i#th shape
+                                     * function at the #j# quadrature point
+                                     * on the cell, face or subface selected
+                                     * the last time the #reinit# function
+                                     * of the derived class was called.
+                                     */
+    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,
+                                     * face or subface selected
+                                     * the last time the #reinit# function
+                                     * of the derived class was called.
+                                     * 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
+                                     * the cell, face or subface selected
+                                     * the last time the #reinit# function
+                                     * of the derived class was called,
+                                     * 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.
@@ -226,6 +325,26 @@ class FEValuesBase {
     DeclException0 (ExcNotImplemented);
     
   protected:
+                                    /**
+                                     * Store the values of the shape functions
+                                     * at the quadrature points. Rows in the
+                                     * matrices denote the values of a single
+                                     * shape function at the different points,
+                                     * columns are for a single point with the
+                                     * different shape functions.
+                                     *
+                                     * For cell values, the vector contains
+                                     * only one entry, representing the
+                                     * restriction of the finite element ansatz
+                                     * space to a cell. For face values, the
+                                     * vector contains as many elements as
+                                     * there are faces, for subfaces the same
+                                     * applies. Which of the matrices is active
+                                     * is determined by the #selected_dataset#
+                                     * variable.
+                                     */
+    vector<dFMatrix>     shape_values;
+
                                     /**
                                      * Store the gradients of the shape
                                      * functions at the quadrature points.
@@ -299,6 +418,16 @@ class FEValuesBase {
                                      */
     vector<dFMatrix>     jacobi_matrices;
 
+                                    /**
+                                     * Store which of the data sets in the
+                                     * #shape_values# array is presently
+                                     * active. This variable is set by the
+                                     * #reinit# functions of the derived
+                                     * classes. For the exact meaning see
+                                     * there and in the doc for this class.
+                                     */
+    unsigned int         selected_dataset;
+    
                                     /**
                                      * Store which fields are to be updated by
                                      * the reinit function.
@@ -352,49 +481,7 @@ class FEValuesBase {
   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  
+  @author Wolfgang Bangerth, 1998  
   */
 template <int dim>
 class FEValues : public FEValuesBase<dim> {
@@ -419,36 +506,6 @@ class FEValues : public FEValuesBase<dim> {
     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
@@ -473,16 +530,6 @@ class FEValues : public FEValuesBase<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
@@ -570,7 +617,7 @@ class FEValues : public FEValuesBase<dim> {
   (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.
+  available without complicated matrix-vector-multiplications.
   */
 template <int dim>
 class FEFaceValues : public FEValuesBase<dim> {
@@ -595,35 +642,6 @@ class FEFaceValues : public FEValuesBase<dim> {
                  const Quadrature<dim-1> &,
                  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;
-
                                     /**
                                      * Return the outward normal vector to
                                      * the cell at the #i#th quadrature
@@ -664,18 +682,6 @@ class FEFaceValues : public FEValuesBase<dim> {
                 const Boundary<dim>                  &boundary);
 
   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.
-                                     *
-                                     * There is one matrix for each face.
-                                     */
-    dFMatrix             shape_values[2*dim];
-
                                     /**
                                      * Store the gradients of the shape
                                      * functions at the quadrature points on
@@ -724,13 +730,6 @@ class FEFaceValues : public FEValuesBase<dim> {
                                      * in by the finite element class.
                                      */
     vector<Point<dim> >  normal_vectors;
-    
-                                    /**
-                                     * Store the number of the face selected
-                                     * last time the #reinit# function was
-                                     * called.
-                                     */
-    unsigned int         selected_face;
 };
 
 
@@ -742,6 +741,16 @@ class FEFaceValues : public FEValuesBase<dim> {
 
 
 
+template <int dim>
+inline
+const dFMatrix & FEValuesBase<dim>::get_shape_values () const {
+  Assert (selected_dataset<shape_values.size(),
+         ExcInvalidIndex (selected_dataset, shape_values.size()));
+  return shape_values[selected_dataset];
+};
+
+
+
 template <int dim>
 inline
 const vector<vector<Point<dim> > > &
@@ -783,28 +792,9 @@ FEValuesBase<dim>::get_JxW_values () const {
 
 
 
-/*------------------------ Inline functions: FEValues ----------------------------*/
-
-
-template <int dim>
-inline
-const dFMatrix & FEValues<dim>::get_shape_values () const {
-  return shape_values;
-};
-
-
-
 /*------------------------ 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<Point<dim> > &
index 8cd9fe4f089c8aeffb409fd3b2749237dc25b7c2..d10efec93ad946b11cd5b544edfee72430d6ee30 100644 (file)
@@ -17,20 +17,64 @@ template <int dim>
 FEValuesBase<dim>::FEValuesBase (const unsigned int n_q_points,
                                 const unsigned int n_ansatz_points,
                                 const unsigned int n_dofs,
+                                const unsigned int n_values_arrays,
                                 const UpdateFlags update_flags) :
                n_quadrature_points (n_q_points),
                total_dofs (n_dofs),
+               shape_values (n_values_arrays, dFMatrix(n_dofs, n_q_points)),
                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)),
+               selected_dataset (0),
                update_flags (update_flags) {};
 
 
 
 
+template <int dim>
+double FEValuesBase<dim>::shape_value (const unsigned int i,
+                                      const unsigned int j) const {
+  Assert (selected_dataset<shape_values.size(),
+         ExcInvalidIndex (selected_dataset, shape_values.size()));
+  Assert (i<shape_values[selected_dataset].m(),
+         ExcInvalidIndex (i, shape_values[selected_dataset].m()));
+  Assert (j<shape_values[selected_dataset].n(),
+         ExcInvalidIndex (j, shape_values[selected_dataset].n()));
+
+  return shape_values[selected_dataset](i,j);
+};
+
+
+
+template <int dim>
+void FEValuesBase<dim>::get_function_values (const dVector  &fe_function,
+                                            vector<double> &values) const {
+  Assert (selected_dataset<shape_values.size(),
+         ExcInvalidIndex (selected_dataset, shape_values.size()));
+  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[selected_dataset](shape_func, point));
+};
+
+
+
 template <int dim>
 const Point<dim> &
 FEValuesBase<dim>::shape_grad (const unsigned int i,
@@ -110,8 +154,8 @@ FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
                FEValuesBase<dim> (quadrature.n_quadrature_points,
                                   fe.total_dofs,
                                   fe.total_dofs,
+                                  1,
                                   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())
@@ -119,7 +163,7 @@ FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
   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));
+       shape_values[0](i,j) = fe.shape_value(i, quadrature.quad_point(j));
        unit_shape_gradients[i][j]
          = fe.shape_grad(i, quadrature.quad_point(j));
       };
@@ -132,39 +176,6 @@ FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
 
 
 
-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));
-};
-
 
 
 
@@ -239,17 +250,17 @@ FEFaceValues<dim>::FEFaceValues (const FiniteElement<dim> &fe,
                FEValuesBase<dim> (quadrature.n_quadrature_points,
                                   fe.dofs_per_face,
                                   fe.total_dofs,
+                                  2*dim,
                                   update_flags),
                unit_quadrature_points(quadrature.get_quad_points()),
                face_jacobi_determinants (quadrature.n_quadrature_points,0),
-               normal_vectors (quadrature.n_quadrature_points,Point<dim>()),
-               selected_face(0)
+               normal_vectors (quadrature.n_quadrature_points,Point<dim>())
 {
   for (unsigned int face=0; face<2*dim; ++face)
     {
-      shape_values[face].reinit(fe.total_dofs, quadrature.n_quadrature_points);
       unit_shape_gradients[face].resize (fe.total_dofs,
-                                        vector<Point<dim> >(quadrature.n_quadrature_points));
+                                        vector<Point<dim> >(quadrature.
+                                                            n_quadrature_points));
       global_unit_quadrature_points[face].resize (quadrature.n_quadrature_points,
                                                  Point<dim>());
     };
@@ -295,10 +306,7 @@ FEFaceValues<dim>::FEFaceValues (const FiniteElement<dim> &fe,
        };
 
   for (unsigned int i=0; i<n_quadrature_points; ++i) 
-    {
-      weights[i] = quadrature.weight(i);
-      unit_quadrature_points[i] = quadrature.quad_point(i);
-    };
+    weights[i] = quadrature.weight(i);
 
   for (unsigned int face=0; face<2*dim; ++face)
     for (unsigned int i=0; i<fe.total_dofs; ++i)
@@ -313,43 +321,6 @@ FEFaceValues<dim>::FEFaceValues (const FiniteElement<dim> &fe,
 
 
 
-template <int dim>
-double FEFaceValues<dim>::shape_value (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()));
-
-  return shape_values[selected_face](i,j);
-};
-
-
-
-template <int dim>
-void FEFaceValues<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[selected_face](shape_func, point));
-};
-
-
-
 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()));
@@ -367,7 +338,7 @@ void FEFaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &c
                                const FiniteElement<dim>                      &fe,
                                const Boundary<dim>                           &boundary) {
   present_cell  = cell;
-  selected_face = face_no;
+  selected_dataset = face_no;
                                   // fill jacobi matrices and real
                                   // quadrature points
   if ((update_flags & update_jacobians) ||

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.