]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add functionality for the evaluation of finite element functions.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 29 Apr 1998 12:08:39 +0000 (12:08 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 29 Apr 1998 12:08:39 +0000 (12:08 +0000)
git-svn-id: https://svn.dealii.org/trunk@223 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 0715c95ad82a49de2dff042e9d4cbd748053cc51..5a2bba2448ac462ca568eaa9d9364d1f5d3f1cd5 100644 (file)
@@ -61,6 +61,48 @@ template <int dim> class Quadrature;
   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}
   */
 template <int dim>
 class FEValues {
@@ -112,6 +154,19 @@ class FEValues {
                                      */
     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 gradient of the #i#th shape
                                      * function at the #j# quadrature point.
@@ -137,6 +192,19 @@ class FEValues {
                                      */
     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.
@@ -225,6 +293,13 @@ class FEValues {
                                      * Exception
                                      */
     DeclException0 (ExcCannotInitializeField);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcWrongVectorSize,
+                   int, int,
+                   << "Vector has wrong size " << arg1
+                   << ", expected size " << arg2);
     
   private:
                                     /**
@@ -319,6 +394,15 @@ class FEValues {
                                      * 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;
 };
 
 
@@ -451,7 +535,20 @@ class FEFaceValues {
                                      * 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 gradient of the #i#th shape
                                      * function at the #j# quadrature point.
@@ -477,6 +574,19 @@ class FEFaceValues {
                                      */
     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.
@@ -604,7 +714,13 @@ class FEFaceValues {
                                      * Exception
                                      */
     DeclException0 (ExcNotImplemented);
-    
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcWrongVectorSize,
+                   int, int,
+                   << "Vector has wrong size " << arg1
+                   << ", expected size " << arg2);
   private:
                                     /**
                                      * Store the values of the shape functions
@@ -744,6 +860,15 @@ class FEFaceValues {
                                      */
     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
index 6b76293f695be8d3d495bee2daf56cc8fa1f6026..a51d5e85d8c697314853a0190e2da81c6a49335c 100644 (file)
@@ -59,6 +59,30 @@ double FEValues<dim>::shape_value (const unsigned int i,
 
 
 
+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>
 const Point<dim> &
 FEValues<dim>::shape_grad (const unsigned int i,
@@ -72,6 +96,30 @@ 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 {
+  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> & FEValues<dim>::quadrature_point (const unsigned int i) const {
   Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
@@ -106,6 +154,7 @@ template <int dim>
 void FEValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &cell,
                            const FiniteElement<dim>                      &fe,
                            const Boundary<dim>                           &boundary) {
+  present_cell = cell;
                                   // fill jacobi matrices and real
                                   // quadrature points
   if ((update_flags & update_jacobians) ||
@@ -209,13 +258,6 @@ FEFaceValues<dim>::FEFaceValues (const FiniteElement<dim> &fe,
                case 0:
                      global_unit_quadrature_points[face][p]
                        = Point<dim>(unit_quadrature_points[p](0),0);
-                     cout << "  face=" << face << ", p=" << p
-                          << ",  uqp[p]=" << unit_quadrature_points[p]
-                          << ", guqp[face][p]="
-                          << global_unit_quadrature_points[face][p]
-                          << ", ??="
-                          << Point<dim>(unit_quadrature_points[p](0),0)
-                          << endl;
                      break;       
                case 1:
                      global_unit_quadrature_points[face][p]
@@ -249,7 +291,8 @@ FEFaceValues<dim>::FEFaceValues (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[face](i,j) = fe.shape_value(i, global_unit_quadrature_points[face][j]);
+         shape_values[face](i,j)
+           = fe.shape_value(i, global_unit_quadrature_points[face][j]);
          unit_shape_gradients[face][i][j]
            = fe.shape_grad(i, global_unit_quadrature_points[face][j]);
        };
@@ -270,6 +313,30 @@ double FEFaceValues<dim>::shape_value (const unsigned int i,
 
 
 
+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>::shape_grad (const unsigned int i,
@@ -278,17 +345,43 @@ FEFaceValues<dim>::shape_grad (const unsigned int i,
          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());
+  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());
+  Assert (update_flags & update_q_points,
+         ExcAccessToUninitializedField());
   
   return quadrature_points[i];
 };
@@ -298,7 +391,8 @@ const Point<dim> & FEFaceValues<dim>::quadrature_point (const unsigned int i) co
 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());
+  Assert (update_flags & update_ansatz_points,
+         ExcAccessToUninitializedField());
   
   return ansatz_points[i];
 };
@@ -308,7 +402,8 @@ const Point<dim> & FEFaceValues<dim>::ansatz_point (const unsigned int i) const
 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()));
-  Assert (update_flags & update_normal_vectors, ExcAccessToUninitializedField());
+  Assert (update_flags & update_normal_vectors,
+         ExcAccessToUninitializedField());
   
   return normal_vectors[i];
 };
@@ -318,7 +413,8 @@ 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());
+  Assert (update_flags & update_JxW_values,
+         ExcAccessToUninitializedField());
   
   return JxW_values[i];
 };
@@ -331,8 +427,7 @@ void FEFaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &c
                                const unsigned int                             face_no,
                                const FiniteElement<dim>                      &fe,
                                const Boundary<dim>                           &boundary) {
-  cout << "ATTENTION: update_flags all set to " << (update_flags=255) << endl;
-  
+  present_cell  = cell;
   selected_face = face_no;
                                   // fill jacobi matrices and real
                                   // quadrature points
@@ -363,19 +458,19 @@ void FEFaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &c
       Assert (update_flags & update_jacobians, ExcCannotInitializeField());
       
       for (unsigned int i=0; i<fe.total_dofs; ++i)
-       for (unsigned int j=0; j<n_quadrature_points; ++j)
-         for (unsigned int s=0; s<dim; ++s)
-           {
-             shape_gradients[i][j](s) = 0;
-             
+       for (unsigned int j=0; j<n_quadrature_points; ++j) 
+         {
+           shape_gradients[i][j] = Point<dim>();
+           
+           for (unsigned int s=0; s<dim; ++s)
                                               // (grad psi)_s =
                                               // (grad_{\xi\eta})_b J_{bs}
                                               // with J_{bs}=(d\xi_b)/(dx_s)
              for (unsigned int b=0; b<dim; ++b)
                shape_gradients[i][j](s)
-                 +=
-                 unit_shape_gradients[face_no][i][j](b) * jacobi_matrices[j](b,s);
-           };
+                 += (unit_shape_gradients[face_no][i][j](b) *
+                     jacobi_matrices[j](b,s));
+         };
     };
   
   
@@ -386,45 +481,11 @@ void FEFaceValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &c
                                   // determinant
   if (update_flags & update_JxW_values) 
     {
-      Assert (update_flags & update_jacobians, ExcCannotInitializeField());
+      Assert (update_flags & update_jacobians,
+             ExcCannotInitializeField());
       for (unsigned int i=0; i<n_quadrature_points; ++i)
        JxW_values[i] = weights[i] * face_jacobi_determinants[i];
     };
-
-  cout << "----------------------------------------" << endl;
-  cout << "cell vertices: \n";
-  for (unsigned int i=0; i<2*dim; ++i)
-    cout << "  " << cell->vertex(i) << endl;
-  cout << "face no= " << face_no << endl;
-
-                                  /*
-  cout << "----------------------------------------" << endl;
-  cout << "Unit face quadrature points:\n";
-  for (unsigned i=0; i<unit_quadrature_points.size(); ++i)
-    cout << "  " << unit_quadrature_points[i];
-  cout << endl;
-                                  */
-  cout << "----------------------------------------" << endl;
-  cout << "Unit global quadrature points:\n";
-  for (unsigned int face=0; face<2*dim; ++face)
-    {
-      cout << "  face=" << face << endl;
-      for (unsigned int i=0; i<unit_quadrature_points.size(); ++i)
-       cout << global_unit_quadrature_points[face][i] << endl;
-    };
-                                  /*
-  cout << "----------------------------------------" << endl;
-  cout << "shape values: \n";
-  for (unsigned int i=0; i<2*dim; ++i) 
-    {
-      cout << "face=" << i << endl;
-      shape_values[i].print_formatted(cout);
-    };
-  cout << endl;
-                                  */
-  
-  cout << "----------------------------------------" << endl;
-  Assert (false,ExcAccessToUninitializedField());
 };
 
 

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.