From: wolf Date: Thu, 30 Apr 1998 14:14:51 +0000 (+0000) Subject: Major reorganisation to allow code sharing and make restriction of finite elements... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a3bbcb633069699f88869c3c2f5b614d3c04a533;p=dealii-svn.git Major reorganisation to allow code sharing and make restriction of finite elements to subfaces possible. git-svn-id: https://svn.dealii.org/trunk@228 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 5a2bba2448..97560206fd 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -19,155 +19,69 @@ template 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 -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 &, - 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; + * 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 &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 & 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 & 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 & 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::cell_iterator &, - const FiniteElement &, - const Boundary &); /** * 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 > > 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 > > 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 weights; @@ -365,15 +275,7 @@ class FEValues { */ vector > 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 > 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 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 +class FEValues : public FEValuesBase { + 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 &, + 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 + * 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::cell_iterator &, + const FiniteElement &, + 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 + * the unit cell. + * This field is set up upon construction + * of the object and contains the gradients + * on the reference element. + */ + vector > > 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 > 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 -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 +class FEFaceValues : public FEValuesBase { + 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 &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 & 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 > > & 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. - * - * 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 & quadrature_point (const unsigned int i) const; - - /** - * Return a pointer to the vector of - * quadrature points. - */ - const vector > & 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 & ansatz_point (const unsigned int i) const; - - /** - * Return a pointer to the vector of points - * denoting the location of the ansatz - * functions. - */ - const vector > & 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 & 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 &fe, const Boundary &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 > > shape_gradients; - /** * Store the gradients of the shape * functions at the quadrature points on @@ -765,36 +688,6 @@ class FEFaceValues { */ vector > > 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 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 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 > quadrature_points; /** * Array of quadrature points on the @@ -817,32 +710,10 @@ class FEFaceValues { */ vector > 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 > 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 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 face_jacobi_determinants; @@ -854,21 +725,6 @@ class FEFaceValues { */ vector > 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::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 -inline -const dFMatrix & FEValues::get_shape_values () const { - return shape_values; -}; +/*------------------------ Inline functions: FEValuesBase ------------------------*/ @@ -897,7 +745,7 @@ const dFMatrix & FEValues::get_shape_values () const { template inline const vector > > & -FEValues::get_shape_grads () const { +FEValuesBase::get_shape_grads () const { Assert (update_flags & update_gradients, ExcAccessToUninitializedField()); return shape_gradients; }; @@ -907,7 +755,7 @@ FEValues::get_shape_grads () const { template inline const vector > & -FEValues::get_quadrature_points () const { +FEValuesBase::get_quadrature_points () const { Assert (update_flags & update_q_points, ExcAccessToUninitializedField()); return quadrature_points; }; @@ -917,7 +765,7 @@ FEValues::get_quadrature_points () const { template inline const vector > & -FEValues::get_ansatz_points () const { +FEValuesBase::get_ansatz_points () const { Assert (update_flags & update_ansatz_points, ExcAccessToUninitializedField()); return ansatz_points; }; @@ -927,7 +775,7 @@ FEValues::get_ansatz_points () const { template inline const vector & -FEValues::get_JxW_values () const { +FEValuesBase::get_JxW_values () const { Assert (update_flags & update_JxW_values, ExcAccessToUninitializedField()); return JxW_values; }; @@ -935,55 +783,24 @@ FEValues::get_JxW_values () const { - -/*------------------------ Inline functions: FEFaceValues ------------------------*/ - - -template -inline -const dFMatrix & FEFaceValues::get_shape_values () const { - return shape_values[selected_face]; -}; - - - - -template -inline -const vector > > & -FEFaceValues::get_shape_grads () const { - Assert (update_flags & update_gradients, ExcAccessToUninitializedField()); - return shape_gradients; -}; - +/*------------------------ Inline functions: FEValues ----------------------------*/ template inline -const vector > & -FEFaceValues::get_quadrature_points () const { - Assert (update_flags & update_q_points, ExcAccessToUninitializedField()); - return quadrature_points; +const dFMatrix & FEValues::get_shape_values () const { + return shape_values; }; -template -inline -const vector > & -FEFaceValues::get_ansatz_points () const { - Assert (update_flags & update_ansatz_points, ExcAccessToUninitializedField()); - return ansatz_points; -}; - +/*------------------------ Inline functions: FEFaceValues ------------------------*/ template inline -const vector & -FEFaceValues::get_JxW_values () const { - Assert (update_flags & update_JxW_values, ExcAccessToUninitializedField()); - return JxW_values; +const dFMatrix & FEFaceValues::get_shape_values () const { + return shape_values[selected_face]; }; diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index a51d5e85d8..8cd9fe4f08 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -6,89 +6,39 @@ #include #include #include +#include -/*------------------------------- FEValues -------------------------------*/ - - -template -FEValues::FEValues (const FiniteElement &fe, - const Quadrature &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 >(quadrature.n_quadrature_points)), - unit_shape_gradients(fe.total_dofs, - vector >(quadrature.n_quadrature_points)), - weights(quadrature.n_quadrature_points, 0), - JxW_values(quadrature.n_quadrature_points, 0), - quadrature_points(quadrature.n_quadrature_points, Point()), - unit_quadrature_points(quadrature.get_quad_points()), - ansatz_points (fe.total_dofs, Point()), - jacobi_matrices (quadrature.n_quadrature_points, - dFMatrix(dim,dim)), - update_flags (update_flags) -{ - for (unsigned int i=0; i -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); +FEValuesBase::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 >(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)), + update_flags (update_flags) {}; - // add up contributions of ansatz - // functions - for (unsigned int point=0; point const Point & -FEValues::shape_grad (const unsigned int i, - const unsigned int j) const { - Assert (i::shape_grad (const unsigned int i, + const unsigned int j) const { + Assert (i::shape_grad (const unsigned int i, template -void FEValues::get_function_grads (const dVector &fe_function, - vector > &gradients) const { +void FEValuesBase::get_function_grads (const dVector &fe_function, + vector > &gradients) const { Assert (gradients.size() == n_quadrature_points, ExcWrongVectorSize(gradients.size(), n_quadrature_points)); @@ -121,7 +71,7 @@ void FEValues::get_function_grads (const dVector &fe_function, template -const Point & FEValues::quadrature_point (const unsigned int i) const { +const Point & FEValuesBase::quadrature_point (const unsigned int i) const { Assert (i & FEValues::quadrature_point (const unsigned int i) const template -const Point & FEValues::ansatz_point (const unsigned int i) const { +const Point & FEValuesBase::ansatz_point (const unsigned int i) const { Assert (i & FEValues::ansatz_point (const unsigned int i) const { template -double FEValues::JxW (const unsigned int i) const { +double FEValuesBase::JxW (const unsigned int i) const { Assert (i::JxW (const unsigned int i) const { + +/*------------------------------- FEValues -------------------------------*/ + +template +FEValues::FEValues (const FiniteElement &fe, + const Quadrature &quadrature, + const UpdateFlags update_flags) : + FEValuesBase (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 >(quadrature.n_quadrature_points)), + unit_quadrature_points(quadrature.get_quad_points()) +{ + for (unsigned int i=0; i +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 void FEValues::reinit (const typename DoFHandler::cell_iterator &cell, const FiniteElement &fe, @@ -217,19 +236,13 @@ template FEFaceValues::FEFaceValues (const FiniteElement &fe, const Quadrature &quadrature, const UpdateFlags update_flags) : - n_quadrature_points(quadrature.n_quadrature_points), - total_dofs(fe.total_dofs), - shape_gradients(fe.total_dofs, - vector >(quadrature.n_quadrature_points)), - weights(quadrature.n_quadrature_points, 0), - JxW_values(quadrature.n_quadrature_points, 0), - quadrature_points(quadrature.n_quadrature_points, Point()), + FEValuesBase (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()), - jacobi_matrices (quadrature.n_quadrature_points,dFMatrix(dim,dim)), face_jacobi_determinants (quadrature.n_quadrature_points,0), normal_vectors (quadrature.n_quadrature_points,Point()), - update_flags (update_flags), selected_face(0) { for (unsigned int face=0; face<2*dim; ++face) @@ -337,68 +350,6 @@ void FEFaceValues::get_function_values (const dVector &fe_function, -template -const Point & -FEFaceValues::shape_grad (const unsigned int i, - const unsigned int j) const { - Assert (i -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 -const Point & FEFaceValues::ansatz_point (const unsigned int i) const { - Assert (i 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 - template void FEFaceValues::reinit (const typename DoFHandler::cell_iterator &cell, const unsigned int face_no, @@ -494,6 +433,9 @@ void FEFaceValues::reinit (const typename DoFHandler::cell_iterator &c /*------------------------------- Explicit Instantiations -------------*/ +template class FEValuesBase<1>; +template class FEValuesBase<2>; + template class FEValues<1>; template class FEValues<2>;