From: Wolfgang Bangerth Date: Wed, 29 Apr 1998 12:08:39 +0000 (+0000) Subject: Add functionality for the evaluation of finite element functions. X-Git-Tag: v8.0.0~23031 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=27ec2e2db58b60a7eb868f03f05bd96e5d415abf;p=dealii.git Add functionality for the evaluation of finite element functions. git-svn-id: https://svn.dealii.org/trunk@223 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index 0715c95ad8..5a2bba2448 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -61,6 +61,48 @@ template 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 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 &values) const; + /** * Return the gradient of the #i#th shape * function at the #j# quadrature point. @@ -137,6 +192,19 @@ class FEValues { */ const vector > > & 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 > &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::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 &values) const; + /** * Return the gradient of the #i#th shape * function at the #j# quadrature point. @@ -477,6 +574,19 @@ class FEFaceValues { */ const vector > > & 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 > &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::cell_iterator present_cell; + /** * Store the number of the face selected * last time the #reinit# function was diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index 6b76293f69..a51d5e85d8 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -59,6 +59,30 @@ double FEValues::shape_value (const unsigned int i, +template +void FEValues::get_function_values (const dVector &fe_function, + vector &values) const { + Assert (values.size() == n_quadrature_points, + ExcWrongVectorSize(values.size(), n_quadrature_points)); + + // get function values of dofs + // on this cell + vector 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 const Point & FEValues::shape_grad (const unsigned int i, @@ -72,6 +96,30 @@ FEValues::shape_grad (const unsigned int i, +template +void FEValues::get_function_grads (const dVector &fe_function, + vector > &gradients) const { + Assert (gradients.size() == n_quadrature_points, + ExcWrongVectorSize(gradients.size(), n_quadrature_points)); + + // get function values of dofs + // on this cell + vector 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()); + + // add up contributions of ansatz + // functions + for (unsigned int point=0; point const Point & FEValues::quadrature_point (const unsigned int i) const { Assert (i void FEValues::reinit (const typename DoFHandler::cell_iterator &cell, const FiniteElement &fe, const Boundary &boundary) { + present_cell = cell; // fill jacobi matrices and real // quadrature points if ((update_flags & update_jacobians) || @@ -209,13 +258,6 @@ FEFaceValues::FEFaceValues (const FiniteElement &fe, case 0: global_unit_quadrature_points[face][p] = Point(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(unit_quadrature_points[p](0),0) - << endl; break; case 1: global_unit_quadrature_points[face][p] @@ -249,7 +291,8 @@ FEFaceValues::FEFaceValues (const FiniteElement &fe, for (unsigned int i=0; i::shape_value (const unsigned int i, +template +void FEFaceValues::get_function_values (const dVector &fe_function, + vector &values) const { + Assert (values.size() == n_quadrature_points, + ExcWrongVectorSize(values.size(), n_quadrature_points)); + + // get function values of dofs + // on this cell + vector 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 const Point & FEFaceValues::shape_grad (const unsigned int i, @@ -278,17 +345,43 @@ FEFaceValues::shape_grad (const unsigned int i, ExcInvalidIndex (i, shape_values[selected_face].m())); Assert (j +void FEFaceValues::get_function_grads (const dVector &fe_function, + vector > &gradients) const { + Assert (gradients.size() == n_quadrature_points, + ExcWrongVectorSize(gradients.size(), n_quadrature_points)); + + // get function values of dofs + // on this cell + vector 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()); + + // add up contributions of ansatz + // functions + for (unsigned int point=0; point const Point & FEFaceValues::quadrature_point (const unsigned int i) const { Assert (i & FEFaceValues::quadrature_point (const unsigned int i) co template const Point & FEFaceValues::ansatz_point (const unsigned int i) const { Assert (i & FEFaceValues::ansatz_point (const unsigned int i) const template const Point & FEFaceValues::normal_vector (const unsigned int i) const { Assert (i & FEFaceValues::normal_vector (const unsigned int i) const template double FEFaceValues::JxW (const unsigned int i) const { Assert (i::reinit (const typename DoFHandler::cell_iterator &c const unsigned int face_no, const FiniteElement &fe, const Boundary &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::reinit (const typename DoFHandler::cell_iterator &c Assert (update_flags & update_jacobians, ExcCannotInitializeField()); for (unsigned int i=0; i(); + + for (unsigned int s=0; s::reinit (const typename DoFHandler::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; ivertex(i) << endl; - cout << "face no= " << face_no << endl; - - /* - cout << "----------------------------------------" << endl; - cout << "Unit face quadrature points:\n"; - for (unsigned i=0; i