From e533b4a00a6a4b2430602c57ba0f723f5551e48f Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 24 Apr 1998 14:30:35 +0000 Subject: [PATCH] File reorganization git-svn-id: https://svn.dealii.org/trunk@192 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe.h | 711 +----------------------- deal.II/deal.II/include/fe/fe_values.h | 731 +++++++++++++++++++++++++ deal.II/deal.II/source/fe/fe.cc | 156 ------ deal.II/deal.II/source/fe/fe_values.cc | 168 ++++++ 4 files changed, 900 insertions(+), 866 deletions(-) create mode 100644 deal.II/deal.II/include/fe/fe_values.h create mode 100644 deal.II/deal.II/source/fe/fe_values.cc diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 8a81163cbb..187bf90e63 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -4,667 +4,12 @@ #define __fe_H /*---------------------------- fe.h ---------------------------*/ -#include -#include #include #include +#include -// forward declarations -template class Boundary; -template class FiniteElement; -template class Quadrature; - - -/** - Provide a set of flags which tells the #FEValues<>::reinit# function, which - fields are to be updated for each cell. E.g. if you do not need the - gradients since you want to assemble the mass matrix, you can switch that - off. By default, all flags are off, i.e. no reinitialization will be done. - - A variable of this type has to be passed to the constructor of the - #FEValues# object. You can select more than one flag by concatenation - using the #|# (bitwise #or#) operator. - */ -enum UpdateFields { - /** - * Default: update nothing. - */ - update_default = 0, - /** - * Compute quadrature points in real - * space (not on unit cell). - */ - update_q_points = 1, - /** - * Transform gradients on unit cell to - * gradients on real cell. - */ - update_gradients = 2, - /** - * Compute jacobian matrices of the - * transform between unit and real cell - * in the evaluation points. - */ - update_jacobians = 4, - /** - * Compute the JxW values (Jacobian - * determinant at the quadrature point - * times the weight of this point). - */ - update_JxW_values = 8, - /** - * Compute the points on the real cell - * on which the ansatz functions are - * located. - */ - update_ansatz_points = 16 -}; - - - - -/** - 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 #UpdateFields# 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. - */ -template -class FEValues { - public: - - - - /** - * Number of quadrature points. - */ - const unsigned int n_quadrature_points; - - /** - * Total number of shape functions. - */ - 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 UpdateFields); - - /** - * 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 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 position of the #i#th - * quadrature point in real space. - */ - 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. - */ - double JxW (const unsigned int i) const; - - /** - * Return a pointer to the array holding - * the JxW values at the different - * 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 Triangulation::cell_iterator &, - const FiniteElement &, - const 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); - - 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. - * 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. - */ - 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. - */ - 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 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 - * 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. - */ - vector jacobi_matrices; - - /** - * Store which fields are to be updated by - * the reinit function. - */ - UpdateFields update_flags; -}; - - - - -/** - 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. - - 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, - the convention is used that the unit cell be the tensor product of the - interval $[-1,1]$, which is to distinguished properly. - - This class is very similar to the #FEValues# class; see there for more - documentation. - */ -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; - - /** - * 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 for the face, which - * has a dimension one less than the - * cell. - * - * 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. - */ - FEFaceValues (const FiniteElement &, - const Quadrature &, - const UpdateFields); - - /** - * 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 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 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. - */ - double JxW (const unsigned int i) const; - - /** - * Return a pointer to the array holding - * the JxW values at the different - * quadrature points. - */ - const vector & get_JxW_values () const; - - /** - * Reinitialize the gradients, Jacobi - * determinants, etc for the given cell - * and the given finite element. - * - * The constructor 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 Triangulation::face_iterator &, - const FiniteElement &, - const 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); - - 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. - * 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. - */ - 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. - */ - 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 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 - * 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. - */ - vector jacobi_matrices; - - /** - * Store a pointer to the object describing - * the boundary of the domain. - */ - const Boundary &boundary; - - /** - * Store which fields are to be updated by - * the reinit function. - */ - UpdateFields update_flags; -}; - @@ -1219,60 +564,6 @@ class FiniteElement<2> : public FiniteElementBase<2> { - -/*------------------------ Inline functions -----------------------------------*/ - - - -template -inline -const dFMatrix & FEValues::get_shape_values () const { - return shape_values; -}; - - - - -template -inline -const vector > > & -FEValues::get_shape_grads () const { - Assert (update_flags | update_gradients, ExcAccessToUninitializedField()); - return shape_gradients; -}; - - - -template -inline -const vector > & -FEValues::get_quadrature_points () const { - Assert (update_flags | update_q_points, ExcAccessToUninitializedField()); - return quadrature_points; -}; - - - -template -inline -const vector > & -FEValues::get_ansatz_points () const { - Assert (update_flags | update_ansatz_points, ExcAccessToUninitializedField()); - return ansatz_points; -}; - - - -template -inline -const vector & -FEValues::get_JxW_values () const { - Assert (update_flags | update_JxW_values, ExcAccessToUninitializedField()); - return JxW_values; -}; - - - /*---------------------------- fe.h ---------------------------*/ /* end of #ifndef __fe_H */ diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h new file mode 100644 index 0000000000..ef030e3b61 --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -0,0 +1,731 @@ +/*---------------------------- fe_values.h ---------------------------*/ +/* $Id$ */ +#ifndef __fe_values_H +#define __fe_values_H +/*---------------------------- fe_values.h ---------------------------*/ + + +#include +#include +#include +#include + + + +// forward declarations +template class Boundary; +template class FiniteElement; +template class Quadrature; + + +/** + Provide a set of flags which tells the #FEValues<>::reinit# function, which + fields are to be updated for each cell. E.g. if you do not need the + gradients since you want to assemble the mass matrix, you can switch that + off. By default, all flags are off, i.e. no reinitialization will be done. + + A variable of this type has to be passed to the constructor of the + #FEValues# object. You can select more than one flag by concatenation + using the #|# (bitwise #or#) operator. + */ +enum UpdateFields { + /** + * Default: update nothing. + */ + update_default = 0, + /** + * Compute quadrature points in real + * space (not on unit cell). + */ + update_q_points = 1, + /** + * Transform gradients on unit cell to + * gradients on real cell. + */ + update_gradients = 2, + /** + * Compute jacobian matrices of the + * transform between unit and real cell + * in the evaluation points. + */ + update_jacobians = 4, + /** + * Compute the JxW values (Jacobian + * determinant at the quadrature point + * times the weight of this point). + */ + update_JxW_values = 8, + /** + * Compute the points on the real cell + * on which the ansatz functions are + * located. + */ + update_ansatz_points = 16 +}; + + + + +/** + 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 #UpdateFields# 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. + */ +template +class FEValues { + public: + + + + /** + * Number of quadrature points. + */ + const unsigned int n_quadrature_points; + + /** + * Total number of shape functions. + */ + 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 UpdateFields); + + /** + * 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 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 position of the #i#th + * quadrature point in real space. + */ + 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. + */ + double JxW (const unsigned int i) const; + + /** + * Return a pointer to the array holding + * the JxW values at the different + * 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 Triangulation::cell_iterator &, + const FiniteElement &, + const 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); + + 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. + * 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. + */ + 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. + */ + 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 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 + * 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. + */ + vector jacobi_matrices; + + /** + * Store which fields are to be updated by + * the reinit function. + */ + UpdateFields update_flags; +}; + + + + +/** + 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. + + 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, + the convention is used that the unit cell be the tensor product of the + interval $[-1,1]$, which is to distinguished properly. + + This class is very similar to the #FEValues# class; see there for more + documentation. + */ +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; + + /** + * 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 for the face, which + * has a dimension one less than the + * cell. + * + * 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. + */ + FEFaceValues (const FiniteElement &, + const Quadrature &, + const UpdateFields); + + /** + * 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 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 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. + */ + double JxW (const unsigned int i) const; + + /** + * Return a pointer to the array holding + * the JxW values at the different + * quadrature points. + */ + const vector & get_JxW_values () const; + + /** + * Reinitialize the gradients, Jacobi + * determinants, etc for the given cell + * and the given finite element. + * + * The constructor 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 Triangulation::face_iterator &, + const FiniteElement &, + const 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); + + 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. + * 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. + */ + 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. + */ + 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 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 + * 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. + */ + vector jacobi_matrices; + + /** + * Store a pointer to the object describing + * the boundary of the domain. + */ + const Boundary &boundary; + + /** + * Store which fields are to be updated by + * the reinit function. + */ + UpdateFields update_flags; +}; + + + + + +/*------------------------ Inline functions -----------------------------------*/ + + + +template +inline +const dFMatrix & FEValues::get_shape_values () const { + return shape_values; +}; + + + + +template +inline +const vector > > & +FEValues::get_shape_grads () const { + Assert (update_flags | update_gradients, ExcAccessToUninitializedField()); + return shape_gradients; +}; + + + +template +inline +const vector > & +FEValues::get_quadrature_points () const { + Assert (update_flags | update_q_points, ExcAccessToUninitializedField()); + return quadrature_points; +}; + + + +template +inline +const vector > & +FEValues::get_ansatz_points () const { + Assert (update_flags | update_ansatz_points, ExcAccessToUninitializedField()); + return ansatz_points; +}; + + + +template +inline +const vector & +FEValues::get_JxW_values () const { + Assert (update_flags | update_JxW_values, ExcAccessToUninitializedField()); + return JxW_values; +}; + + + + + +/*---------------------------- fe_values.h ---------------------------*/ +/* end of #ifndef __fe_values_H */ +#endif +/*---------------------------- fe_values.h ---------------------------*/ diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 911f77d438..f4a962c847 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -9,159 +9,6 @@ -/*------------------------------- FEValues -------------------------------*/ - - -template -FEValues::FEValues (const FiniteElement &fe, - const Quadrature &quadrature, - const UpdateFields 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.n_quadrature_points, Point()), - 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 -const Point & -FEValues::shape_grad (const unsigned int i, - const unsigned int j) const { - Assert (i -const Point & FEValues::quadrature_point (const unsigned int i) const { - Assert (i -const Point & FEValues::ansatz_point (const unsigned int i) const { - Assert (i -double FEValues::JxW (const unsigned int i) const { - Assert (i -void FEValues::reinit (const typename Triangulation::cell_iterator &cell, - const FiniteElement &fe, - const Boundary &boundary) { - // fill jacobi matrices and real - // quadrature points - if ((update_flags | update_jacobians) || - (update_flags | update_q_points)) - fe.fill_fe_values (cell, - unit_quadrature_points, - jacobi_matrices, - update_flags | update_jacobians, - ansatz_points, - update_flags | update_ansatz_points, - quadrature_points, - update_flags | update_q_points, - boundary); - - // compute gradients on real element if - // requested - if (update_flags | update_gradients) - { - Assert (update_flags | update_jacobians, ExcCannotInitializeField()); - - for (unsigned int i=0; i::face_ansatz_points (const typename Triangulation<2>::face /*------------------------------- Explicit Instantiations -------------*/ -template class FEValues<1>; -template class FEValues<2>; - template class FiniteElementBase<1>; template class FiniteElementBase<2>; diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc new file mode 100644 index 0000000000..9f127da495 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -0,0 +1,168 @@ +/* $Id$ */ + +#include +#include +#include +#include +#include +#include + + + +/*------------------------------- FEValues -------------------------------*/ + + +template +FEValues::FEValues (const FiniteElement &fe, + const Quadrature &quadrature, + const UpdateFields 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.n_quadrature_points, Point()), + 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 +const Point & +FEValues::shape_grad (const unsigned int i, + const unsigned int j) const { + Assert (i +const Point & FEValues::quadrature_point (const unsigned int i) const { + Assert (i +const Point & FEValues::ansatz_point (const unsigned int i) const { + Assert (i +double FEValues::JxW (const unsigned int i) const { + Assert (i +void FEValues::reinit (const typename Triangulation::cell_iterator &cell, + const FiniteElement &fe, + const Boundary &boundary) { + // fill jacobi matrices and real + // quadrature points + if ((update_flags | update_jacobians) || + (update_flags | update_q_points)) + fe.fill_fe_values (cell, + unit_quadrature_points, + jacobi_matrices, + update_flags | update_jacobians, + ansatz_points, + update_flags | update_ansatz_points, + quadrature_points, + update_flags | update_q_points, + boundary); + + // compute gradients on real element if + // requested + if (update_flags | update_gradients) + { + Assert (update_flags | update_jacobians, ExcCannotInitializeField()); + + for (unsigned int i=0; i; +template class FEValues<2>; + -- 2.39.5