From 2baea8f88fed3f461d42ef07229f7a29215223bc Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 4 May 1998 08:39:09 +0000 Subject: [PATCH] Next step towards restriction of finite elements to subfaces. git-svn-id: https://svn.dealii.org/trunk@239 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_values.h | 294 ++++++++++++------------- deal.II/deal.II/source/fe/fe_values.cc | 133 +++++------ 2 files changed, 194 insertions(+), 233 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index 97560206fd..81eb1c2738 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -20,7 +20,7 @@ template 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 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 @@ -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 &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 shape_values; + /** * Store the gradients of the shape * functions at the quadrature points. @@ -299,6 +418,16 @@ class FEValuesBase { */ vector 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 class FEValues : public FEValuesBase { @@ -419,36 +506,6 @@ class FEValues : public FEValuesBase { FEValues (const FiniteElement &, const Quadrature &, 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 &values) const; - /** * Reinitialize the gradients, Jacobi @@ -473,16 +530,6 @@ class FEValues : public FEValuesBase { const 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. - */ - dFMatrix shape_values; - /** * Store the gradients of the shape * functions at the quadrature points on @@ -570,7 +617,7 @@ class FEValues : public FEValuesBase { (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 class FEFaceValues : public FEValuesBase { @@ -595,35 +642,6 @@ class FEFaceValues : public FEValuesBase { const Quadrature &, 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 &values) const; - /** * Return the outward normal vector to * the cell at the #i#th quadrature @@ -664,18 +682,6 @@ class FEFaceValues : public FEValuesBase { const Boundary &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 { * in by the finite element class. */ vector > 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 { +template +inline +const dFMatrix & FEValuesBase::get_shape_values () const { + Assert (selected_dataset inline const vector > > & @@ -783,28 +792,9 @@ FEValuesBase::get_JxW_values () const { -/*------------------------ Inline functions: FEValues ----------------------------*/ - - -template -inline -const dFMatrix & FEValues::get_shape_values () const { - return shape_values; -}; - - - /*------------------------ Inline functions: FEFaceValues ------------------------*/ -template -inline -const dFMatrix & FEFaceValues::get_shape_values () const { - return shape_values[selected_face]; -}; - - - template inline const vector > & diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index 8cd9fe4f08..d10efec93a 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -17,20 +17,64 @@ template FEValuesBase::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 >(n_q_points)), weights (n_q_points, 0), JxW_values (n_q_points, 0), quadrature_points (n_q_points, Point()), ansatz_points (n_ansatz_points, Point()), jacobi_matrices (n_q_points, dFMatrix(dim,dim)), + selected_dataset (0), update_flags (update_flags) {}; +template +double FEValuesBase::shape_value (const unsigned int i, + const unsigned int j) const { + Assert (selected_dataset +void FEValuesBase::get_function_values (const dVector &fe_function, + vector &values) const { + Assert (selected_dataset 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 & FEValuesBase::shape_grad (const unsigned int i, @@ -110,8 +154,8 @@ FEValues::FEValues (const FiniteElement &fe, FEValuesBase (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 >(quadrature.n_quadrature_points)), unit_quadrature_points(quadrature.get_quad_points()) @@ -119,7 +163,7 @@ FEValues::FEValues (const FiniteElement &fe, for (unsigned int i=0; i::FEValues (const FiniteElement &fe, -template -double FEValues::shape_value (const unsigned int i, - const unsigned int j) const { - Assert (i -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::FEFaceValues (const FiniteElement &fe, FEValuesBase (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()), - selected_face(0) + normal_vectors (quadrature.n_quadrature_points,Point()) { 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 >(quadrature.n_quadrature_points)); + vector >(quadrature. + n_quadrature_points)); global_unit_quadrature_points[face].resize (quadrature.n_quadrature_points, Point()); }; @@ -295,10 +306,7 @@ FEFaceValues::FEFaceValues (const FiniteElement &fe, }; for (unsigned int i=0; i::FEFaceValues (const FiniteElement &fe, -template -double FEFaceValues::shape_value (const unsigned int i, - const unsigned int j) const { - Assert (i -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::normal_vector (const unsigned int i) const { Assert (i::reinit (const typename DoFHandler::cell_iterator &c const FiniteElement &fe, const Boundary &boundary) { present_cell = cell; - selected_face = face_no; + selected_dataset = face_no; // fill jacobi matrices and real // quadrature points if ((update_flags & update_jacobians) || -- 2.39.5