From 9bf90fb8fbd5eeb5d5e09a7d88fa7948b9315254 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Fri, 24 May 2013 09:53:12 +0000 Subject: [PATCH] Initialize restriction and prolongation matrix in FESystem and FE_Q only upon request. This significantly reduces memory requirements for higher order elements as most of these matrices are never needed, e.g. anisotropic ones without anisotropic refinement. git-svn-id: https://svn.dealii.org/trunk@29561 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/fe/fe.h | 2367 ++++++++---------------- deal.II/include/deal.II/fe/fe_q_base.h | 73 +- deal.II/include/deal.II/fe/fe_system.h | 745 +++----- deal.II/source/fe/fe.cc | 34 +- deal.II/source/fe/fe_q_base.cc | 503 ++--- deal.II/source/fe/fe_system.cc | 1072 +++++------ 6 files changed, 1786 insertions(+), 3008 deletions(-) diff --git a/deal.II/include/deal.II/fe/fe.h b/deal.II/include/deal.II/fe/fe.h index 97ba56192f..fbf458b9fe 100644 --- a/deal.II/include/deal.II/fe/fe.h +++ b/deal.II/include/deal.II/fe/fe.h @@ -343,14 +343,11 @@ class FiniteElement : public Subscriptor, { public: /** - * Base class for internal data. - * Adds data for second derivatives to + * Base class for internal data. Adds data for second derivatives to * Mapping::InternalDataBase() * - * For information about the - * general purpose of this class, - * see the documentation of the - * base class. + * For information about the general purpose of this class, see the + * documentation of the base class. * * @author Guido Kanschat, 2001 */ @@ -358,35 +355,24 @@ public: { public: /** - * Destructor. Needed to - * avoid memory leaks with - * difference quotients. + * Destructor. Needed to avoid memory leaks with difference quotients. */ virtual ~InternalDataBase (); /** - * Initialize some pointers - * used in the computation of - * second derivatives by - * finite differencing of - * gradients. + * Initialize some pointers used in the computation of second derivatives + * by finite differencing of gradients. */ void initialize_2nd (const FiniteElement *element, const Mapping &mapping, const Quadrature &quadrature); /** - * Storage for FEValues - * objects needed to - * approximate second - * derivatives. + * Storage for FEValues objects needed to approximate second derivatives. * - * The ordering is p+hx, - * p+hy, p+hz, - * p-hx, p-hy, - * p-hz, where unused - * entries in lower dimensions - * are missing. + * The ordering is p+hx, p+hy, p+hz, p-hx, + * p-hy, p-hz, where unused entries in lower dimensions are + * missing. */ std::vector*> differences; }; @@ -400,74 +386,43 @@ public: const std::vector &nonzero_components); /** - * Virtual destructor. Makes sure - * that pointers to this class - * are deleted properly. + * Virtual destructor. Makes sure that pointers to this class are deleted + * properly. */ virtual ~FiniteElement (); /** - * Return a string that uniquely - * identifies a finite - * element. The general - * convention is that this is the - * class name, followed by the - * dimension in angle - * brackets, and the polynomial - * degree and whatever else is - * necessary in parentheses. For - * example, FE_Q<2>(3) is the - * value returned for a cubic - * element in 2d. + * Return a string that uniquely identifies a finite element. The general + * convention is that this is the class name, followed by the dimension in + * angle brackets, and the polynomial degree and whatever else is necessary + * in parentheses. For example, FE_Q<2>(3) is the value returned + * for a cubic element in 2d. * - * Systems of elements have their - * own naming convention, see the - * FESystem class. + * Systems of elements have their own naming convention, see the FESystem + * class. */ virtual std::string get_name () const = 0; /** - * This operator returns a - * reference to the present - * object if the argument given - * equals to zero. While this - * does not seem particularly - * useful, it is helpful in - * writing code that works with - * both ::DoFHandler and the hp - * version hp::DoFHandler, since - * one can then write code like - * this: + * This operator returns a reference to the present object if the argument + * given equals to zero. While this does not seem particularly useful, it is + * helpful in writing code that works with both ::DoFHandler and the hp + * version hp::DoFHandler, since one can then write code like this: * @code * dofs_per_cell * = dof_handler->get_fe()[cell->active_fe_index()].dofs_per_cell; * @endcode * - * This code doesn't work in both - * situations without the present - * operator because - * DoFHandler::get_fe() returns a - * finite element, whereas - * hp::DoFHandler::get_fe() - * returns a collection of finite - * elements that doesn't offer a - * dofs_per_cell - * member variable: one first has - * to select which finite element - * to work on, which is done - * using the - * operator[]. Fortunately, - * cell-@>active_fe_index() - * also works for non-hp classes - * and simply returns zero in - * that case. The present - * operator[] accepts this zero - * argument, by returning the - * finite element with index zero - * within its collection (that, - * of course, consists only of - * the present finite element - * anyway). + * This code doesn't work in both situations without the present operator + * because DoFHandler::get_fe() returns a finite element, whereas + * hp::DoFHandler::get_fe() returns a collection of finite elements that + * doesn't offer a dofs_per_cell member variable: one first has + * to select which finite element to work on, which is done using the + * operator[]. Fortunately, cell-@>active_fe_index() also works + * for non-hp classes and simply returns zero in that case. The present + * operator[] accepts this zero argument, by returning the finite element + * with index zero within its collection (that, of course, consists only of + * the present finite element anyway). */ const FiniteElement &operator[] (const unsigned int fe_index) const; @@ -477,139 +432,75 @@ public: */ /** - * Return the value of the - * @p ith shape function at the - * point @p p. @p p is a point - * on the reference element. If - * the finite element is - * vector-valued, then return the - * value of the only non-zero - * component of the vector value - * of this shape function. If the - * shape function has more than - * one non-zero component (which - * we refer to with the term - * non-primitive), then derived - * classes implementing this - * function should throw an - * exception of type - * ExcShapeFunctionNotPrimitive. In - * that case, use the - * shape_value_component() - * function. + * Return the value of the @p ith shape function at the point @p p. @p p is + * a point on the reference element. If the finite element is vector-valued, + * then return the value of the only non-zero component of the vector value + * of this shape function. If the shape function has more than one non-zero + * component (which we refer to with the term non-primitive), then derived + * classes implementing this function should throw an exception of type + * ExcShapeFunctionNotPrimitive. In that case, use the + * shape_value_component() function. * - * An - * ExcUnitShapeValuesDoNotExist - * is thrown if the shape values - * of the FiniteElement under - * consideration depends on the - * shape of the cell in real - * space. + * An ExcUnitShapeValuesDoNotExist is thrown if the shape values of the + * FiniteElement under consideration depends on the shape of the cell in + * real space. */ virtual double shape_value (const unsigned int i, const Point &p) const; /** - * Just like for shape_value(), - * but this function will be - * called when the shape function - * has more than one non-zero - * vector component. In that - * case, this function should - * return the value of the - * @p component-th vector - * component of the @p ith shape - * function at point @p p. + * Just like for shape_value(), but this function will be called when the + * shape function has more than one non-zero vector component. In that case, + * this function should return the value of the @p component-th vector + * component of the @p ith shape function at point @p p. */ virtual double shape_value_component (const unsigned int i, const Point &p, const unsigned int component) const; /** - * Return the gradient of the - * @p ith shape function at the - * point @p p. @p p is a point - * on the reference element, and - * likewise the gradient is the - * gradient on the unit cell with - * respect to unit cell - * coordinates. If - * the finite element is - * vector-valued, then return the - * value of the only non-zero - * component of the vector value - * of this shape function. If the - * shape function has more than - * one non-zero component (which - * we refer to with the term - * non-primitive), then derived - * classes implementing this - * function should throw an - * exception of type - * ExcShapeFunctionNotPrimitive. In - * that case, use the - * shape_grad_component() - * function. + * Return the gradient of the @p ith shape function at the point @p p. @p p + * is a point on the reference element, and likewise the gradient is the + * gradient on the unit cell with respect to unit cell coordinates. If the + * finite element is vector-valued, then return the value of the only + * non-zero component of the vector value of this shape function. If the + * shape function has more than one non-zero component (which we refer to + * with the term non-primitive), then derived classes implementing this + * function should throw an exception of type + * ExcShapeFunctionNotPrimitive. In that case, use the + * shape_grad_component() function. * - * An - * ExcUnitShapeValuesDoNotExist - * is thrown if the shape values - * of the FiniteElement under - * consideration depends on the - * shape of the cell in real - * space. + * An ExcUnitShapeValuesDoNotExist is thrown if the shape values of the + * FiniteElement under consideration depends on the shape of the cell in + * real space. */ virtual Tensor<1,dim> shape_grad (const unsigned int i, const Point &p) const; /** - * Just like for shape_grad(), - * but this function will be - * called when the shape function - * has more than one non-zero - * vector component. In that - * case, this function should - * return the gradient of the - * @p component-th vector - * component of the @p ith shape - * function at point @p p. + * Just like for shape_grad(), but this function will be called when the + * shape function has more than one non-zero vector component. In that case, + * this function should return the gradient of the @p component-th vector + * component of the @p ith shape function at point @p p. */ virtual Tensor<1,dim> shape_grad_component (const unsigned int i, const Point &p, const unsigned int component) const; /** - * Return the tensor of second - * derivatives of the @p ith - * shape function at point @p p - * on the unit cell. The - * derivatives are derivatives on - * the unit cell with respect to - * unit cell coordinates. If - * the finite element is - * vector-valued, then return the - * value of the only non-zero - * component of the vector value - * of this shape function. If the - * shape function has more than - * one non-zero component (which - * we refer to with the term - * non-primitive), then derived - * classes implementing this - * function should throw an - * exception of type - * ExcShapeFunctionNotPrimitive. In - * that case, use the - * shape_grad_grad_component() - * function. + * Return the tensor of second derivatives of the @p ith shape function at + * point @p p on the unit cell. The derivatives are derivatives on the unit + * cell with respect to unit cell coordinates. If the finite element is + * vector-valued, then return the value of the only non-zero component of + * the vector value of this shape function. If the shape function has more + * than one non-zero component (which we refer to with the term + * non-primitive), then derived classes implementing this function should + * throw an exception of type ExcShapeFunctionNotPrimitive. In that case, + * use the shape_grad_grad_component() function. * - * An - * ExcUnitShapeValuesDoNotExist - * is thrown if the shape values - * of the FiniteElement under - * consideration depends on the - * shape of the cell in real - * space. + * An ExcUnitShapeValuesDoNotExist is thrown if the shape values of the + * FiniteElement under consideration depends on the shape of the cell in + * real space. */ virtual Tensor<2,dim> shape_grad_grad (const unsigned int i, const Point &p) const; @@ -630,21 +521,14 @@ public: const Point &p, const unsigned int component) const; /** - * Check for non-zero values on a - * face in order to optimize out - * matrix elements. + * Check for non-zero values on a face in order to optimize out matrix + * elements. * - * This function returns - * @p true, if the shape - * function @p shape_index has - * non-zero values on the face - * @p face_index. + * This function returns @p true, if the shape function @p shape_index has + * non-zero values on the face @p face_index. * - * A default implementation is - * provided in this basis class - * which always returns @p - * true. This is the safe way to - * go. + * A default implementation is provided in this basis class which always + * returns @p true. This is the safe way to go. */ virtual bool has_support_on_face (const unsigned int shape_index, const unsigned int face_index) const; @@ -656,374 +540,214 @@ public: */ /** - * Projection from a fine grid - * space onto a coarse grid - * space. If this projection - * operator is associated with a - * matrix @p P, then the - * restriction of this matrix - * @p P_i to a single child cell - * is returned here. - * - * The matrix @p P is the - * concatenation or the sum of - * the cell matrices @p P_i, - * depending on the - * #restriction_is_additive_flags. This - * distinguishes interpolation - * (concatenation) and projection - * with respect to scalar + * Projection from a fine grid space onto a coarse grid space. If this + * projection operator is associated with a matrix @p P, then the + * restriction of this matrix @p P_i to a single child cell is returned + * here. + * + * The matrix @p P is the concatenation or the sum of the cell matrices @p + * P_i, depending on the #restriction_is_additive_flags. This distinguishes + * interpolation (concatenation) and projection with respect to scalar * products (summation). * - * Row and column indices are - * related to coarse grid and - * fine grid spaces, - * respectively, consistent with - * the definition of the - * associated operator. - * - * If projection matrices are not - * implemented in the derived - * finite element class, this - * function aborts with - * ExcProjectionVoid. You can - * check whether this is the case - * by calling the - * restriction_is_implemented() - * or the - * isotropic_restriction_is_implemented() - * function. + * Row and column indices are related to coarse grid and fine grid spaces, + * respectively, consistent with the definition of the associated operator. + * + * If projection matrices are not implemented in the derived finite element + * class, this function aborts with ExcProjectionVoid. You can check whether + * this is the case by calling the restriction_is_implemented() or the + * isotropic_restriction_is_implemented() function. */ - const FullMatrix & + virtual const FullMatrix & get_restriction_matrix (const unsigned int child, const RefinementCase &refinement_case=RefinementCase::isotropic_refinement) const; /** * Embedding matrix between grids. * - * The identity operator from a - * coarse grid space into a fine - * grid space is associated with - * a matrix @p P. The - * restriction of this matrix @p P_i to - * a single child cell is - * returned here. - * - * The matrix @p P is the - * concatenation, not the sum of - * the cell matrices - * @p P_i. That is, if the same - * non-zero entry j,k exists - * in in two different child - * matrices @p P_i, the value - * should be the same in both - * matrices and it is copied into - * the matrix @p P only once. - * - * Row and column indices are - * related to fine grid and - * coarse grid spaces, - * respectively, consistent with - * the definition of the - * associated operator. - * - * These matrices are used by - * routines assembling the - * prolongation matrix for - * multi-level methods. Upon - * assembling the transfer matrix - * between cells using this - * matrix array, zero elements in - * the prolongation matrix are - * discarded and will not fill up - * the transfer matrix. - * - * If projection matrices are not - * implemented in the derived - * finite element class, this - * function aborts with - * ExcEmbeddingVoid. You can - * check whether this is the case - * by calling the - * prolongation_is_implemented() - * or the - * isotropic_prolongation_is_implemented() - * function. + * The identity operator from a coarse grid space into a fine grid space is + * associated with a matrix @p P. The restriction of this matrix @p P_i to a + * single child cell is returned here. + * + * The matrix @p P is the concatenation, not the sum of the cell matrices @p + * P_i. That is, if the same non-zero entry j,k exists in in two + * different child matrices @p P_i, the value should be the same in both + * matrices and it is copied into the matrix @p P only once. + * + * Row and column indices are related to fine grid and coarse grid spaces, + * respectively, consistent with the definition of the associated operator. + * + * These matrices are used by routines assembling the prolongation matrix + * for multi-level methods. Upon assembling the transfer matrix between + * cells using this matrix array, zero elements in the prolongation matrix + * are discarded and will not fill up the transfer matrix. + * + * If projection matrices are not implemented in the derived finite element + * class, this function aborts with ExcEmbeddingVoid. You can check whether + * this is the case by calling the prolongation_is_implemented() or the + * isotropic_prolongation_is_implemented() function. */ - const FullMatrix & + virtual const FullMatrix & get_prolongation_matrix (const unsigned int child, const RefinementCase &refinement_case=RefinementCase::isotropic_refinement) const; /** - * Return whether this element implements - * its prolongation matrices. The return - * value also indicates whether a call to - * the get_prolongation_matrix() - * function will generate an error or - * not. - * - * Note, that this function - * returns true only - * if the prolongation matrices of - * the isotropic and all - * anisotropic refinement cases - * are implemented. If you are - * interested in the prolongation - * matrices for isotropic - * refinement only, use the - * isotropic_prolongation_is_implemented - * function instead. - * - * This function is mostly here in order - * to allow us to write more efficient - * test programs which we run on all - * kinds of weird elements, and for which - * we simply need to exclude certain - * tests in case something is not - * implemented. It will in general - * probably not be a great help in - * applications, since there is not much - * one can do if one needs these features - * and they are not implemented. This - * function could be used to check - * whether a call to - * get_prolongation_matrix() will - * succeed; however, one then still needs - * to cope with the lack of information - * this just expresses. + * Return whether this element implements its prolongation matrices. The + * return value also indicates whether a call to the + * get_prolongation_matrix() function will generate an error or not. + * + * Note, that this function returns true only if the + * prolongation matrices of the isotropic and all anisotropic refinement + * cases are implemented. If you are interested in the prolongation matrices + * for isotropic refinement only, use the + * isotropic_prolongation_is_implemented function instead. + * + * This function is mostly here in order to allow us to write more efficient + * test programs which we run on all kinds of weird elements, and for which + * we simply need to exclude certain tests in case something is not + * implemented. It will in general probably not be a great help in + * applications, since there is not much one can do if one needs these + * features and they are not implemented. This function could be used to + * check whether a call to get_prolongation_matrix() will succeed; + * however, one then still needs to cope with the lack of information this + * just expresses. */ bool prolongation_is_implemented () const; /** - * Return whether this element implements - * its prolongation matrices for isotropic - * children. The return value also - * indicates whether a call to the @p - * get_prolongation_matrix function will - * generate an error or not. - * - * This function is mostly here in order - * to allow us to write more efficient - * test programs which we run on all - * kinds of weird elements, and for which - * we simply need to exclude certain - * tests in case something is not - * implemented. It will in general - * probably not be a great help in - * applications, since there is not much - * one can do if one needs these features - * and they are not implemented. This - * function could be used to check - * whether a call to - * get_prolongation_matrix() will - * succeed; however, one then still needs - * to cope with the lack of information - * this just expresses. + * Return whether this element implements its prolongation matrices for + * isotropic children. The return value also indicates whether a call to the + * @p get_prolongation_matrix function will generate an error or not. + * + * This function is mostly here in order to allow us to write more efficient + * test programs which we run on all kinds of weird elements, and for which + * we simply need to exclude certain tests in case something is not + * implemented. It will in general probably not be a great help in + * applications, since there is not much one can do if one needs these + * features and they are not implemented. This function could be used to + * check whether a call to get_prolongation_matrix() will succeed; + * however, one then still needs to cope with the lack of information this + * just expresses. */ bool isotropic_prolongation_is_implemented () const; /** - * Return whether this element implements - * its restriction matrices. The return - * value also indicates whether a call to - * the get_restriction_matrix() - * function will generate an error or - * not. - * - * Note, that this function - * returns true only - * if the restriction matrices of - * the isotropic and all - * anisotropic refinement cases - * are implemented. If you are - * interested in the restriction - * matrices for isotropic - * refinement only, use the - * isotropic_restriction_is_implemented - * function instead. - * - * This function is mostly here in order - * to allow us to write more efficient - * test programs which we run on all - * kinds of weird elements, and for which - * we simply need to exclude certain - * tests in case something is not - * implemented. It will in general - * probably not be a great help in - * applications, since there is not much - * one can do if one needs these features - * and they are not implemented. This - * function could be used to check - * whether a call to - * get_restriction_matrix() will - * succeed; however, one then still needs - * to cope with the lack of information - * this just expresses. + * Return whether this element implements its restriction matrices. The + * return value also indicates whether a call to the + * get_restriction_matrix() function will generate an error or not. + * + * Note, that this function returns true only if the + * restriction matrices of the isotropic and all anisotropic refinement + * cases are implemented. If you are interested in the restriction matrices + * for isotropic refinement only, use the + * isotropic_restriction_is_implemented function instead. + * + * This function is mostly here in order to allow us to write more efficient + * test programs which we run on all kinds of weird elements, and for which + * we simply need to exclude certain tests in case something is not + * implemented. It will in general probably not be a great help in + * applications, since there is not much one can do if one needs these + * features and they are not implemented. This function could be used to + * check whether a call to get_restriction_matrix() will succeed; + * however, one then still needs to cope with the lack of information this + * just expresses. */ bool restriction_is_implemented () const; /** - * Return whether this element implements - * its restriction matrices for isotropic - * children. The return value also - * indicates whether a call to the @p - * get_restriction_matrix function will - * generate an error or not. - * - * This function is mostly here in order - * to allow us to write more efficient - * test programs which we run on all - * kinds of weird elements, and for which - * we simply need to exclude certain - * tests in case something is not - * implemented. It will in general - * probably not be a great help in - * applications, since there is not much - * one can do if one needs these features - * and they are not implemented. This - * function could be used to check - * whether a call to - * get_restriction_matrix() will - * succeed; however, one then still needs - * to cope with the lack of information - * this just expresses. + * Return whether this element implements its restriction matrices for + * isotropic children. The return value also indicates whether a call to the + * @p get_restriction_matrix function will generate an error or not. + * + * This function is mostly here in order to allow us to write more efficient + * test programs which we run on all kinds of weird elements, and for which + * we simply need to exclude certain tests in case something is not + * implemented. It will in general probably not be a great help in + * applications, since there is not much one can do if one needs these + * features and they are not implemented. This function could be used to + * check whether a call to get_restriction_matrix() will succeed; + * however, one then still needs to cope with the lack of information this + * just expresses. */ bool isotropic_restriction_is_implemented () const; /** - * Access the - * #restriction_is_additive_flags - * field. See there for more + * Access the #restriction_is_additive_flags field. See there for more * information on its contents. * - * The index must be between zero - * and the number of shape - * functions of this element. + * The index must be between zero and the number of shape functions of this + * element. */ bool restriction_is_additive (const unsigned int index) const; /** - * Return a read only reference to - * the matrix which describes the - * constraints at the interface - * between a refined and an - * unrefined cell. - * - * The matrix is obviously empty - * in only one dimension, - * since there are no constraints - * then. - * - * Note that some finite elements - * do not (yet) implement hanging - * node constraints. If this is - * the case, then this function - * will generate an exception, - * since no useful return value - * can be generated. If you - * should have a way to live with - * this, then you might want to - * use the - * constraints_are_implemented() - * function to check up front - * whether this function will - * succeed or generate the - * exception. + * Return a read only reference to the matrix which describes the + * constraints at the interface between a refined and an unrefined cell. + * + * The matrix is obviously empty in only one dimension, since there are no + * constraints then. + * + * Note that some finite elements do not (yet) implement hanging node + * constraints. If this is the case, then this function will generate an + * exception, since no useful return value can be generated. If you should + * have a way to live with this, then you might want to use the + * constraints_are_implemented() function to check up front whether this + * function will succeed or generate the exception. */ const FullMatrix &constraints (const dealii::internal::SubfaceCase &subface_case=dealii::internal::SubfaceCase::case_isotropic) const; /** - * Return whether this element - * implements its hanging node - * constraints. The return value - * also indicates whether a call - * to the constraints() function + * Return whether this element implements its hanging node constraints. The + * return value also indicates whether a call to the constraints() function * will generate an error or not. * - * This function is mostly here - * in order to allow us to write - * more efficient test programs - * which we run on all kinds of - * weird elements, and for which - * we simply need to exclude - * certain tests in case hanging - * node constraints are not - * implemented. It will in - * general probably not be a - * great help in applications, - * since there is not much one - * can do if one needs hanging - * node constraints and they are - * not implemented. This function - * could be used to check whether - * a call to constraints() - * will succeed; however, one - * then still needs to cope with - * the lack of information this + * This function is mostly here in order to allow us to write more efficient + * test programs which we run on all kinds of weird elements, and for which + * we simply need to exclude certain tests in case hanging node constraints + * are not implemented. It will in general probably not be a great help in + * applications, since there is not much one can do if one needs hanging + * node constraints and they are not implemented. This function could be + * used to check whether a call to constraints() will succeed; + * however, one then still needs to cope with the lack of information this * just expresses. */ bool constraints_are_implemented (const dealii::internal::SubfaceCase &subface_case=dealii::internal::SubfaceCase::case_isotropic) const; /** - * Return whether this element - * implements its hanging node - * constraints in the new way, - * which has to be used to make - * elements "hp compatible". - * That means, the element properly - * implements the - * get_face_interpolation_matrix - * and get_subface_interpolation_matrix - * methods. Therefore the return - * value also indicates whether a call - * to the get_face_interpolation_matrix() - * method and the get_subface_interpolation_matrix() - * method will generate an error or not. - * - * Currently the main purpose of this - * function is to allow the - * make_hanging_node_constraints method - * to decide whether the new procedures, - * which are supposed to work in the hp - * framework can be used, or if the old - * well verified but not hp capable - * functions should be used. Once the - * transition to the new scheme for - * computing the interface constraints is - * complete, this function will be - * superfluous and will probably go away. - * - * Derived classes should implement this - * function accordingly. The default - * assumption is that a finite element - * does not provide hp capable face - * interpolation, and the default - * implementation therefore returns @p - * false. + * Return whether this element implements its hanging node constraints in + * the new way, which has to be used to make elements "hp compatible". That + * means, the element properly implements the get_face_interpolation_matrix + * and get_subface_interpolation_matrix methods. Therefore the return value + * also indicates whether a call to the get_face_interpolation_matrix() + * method and the get_subface_interpolation_matrix() method will generate an + * error or not. + * + * Currently the main purpose of this function is to allow the + * make_hanging_node_constraints method to decide whether the new + * procedures, which are supposed to work in the hp framework can be used, + * or if the old well verified but not hp capable functions should be used. + * Once the transition to the new scheme for computing the interface + * constraints is complete, this function will be superfluous and will + * probably go away. + * + * Derived classes should implement this function accordingly. The default + * assumption is that a finite element does not provide hp capable face + * interpolation, and the default implementation therefore returns @p false. */ virtual bool hp_constraints_are_implemented () const; /** - * Return the matrix - * interpolating from the given - * finite element to the present - * one. The size of the matrix is - * then #dofs_per_cell times + * Return the matrix interpolating from the given finite element to the + * present one. The size of the matrix is then #dofs_per_cell times * source.#dofs_per_cell. * - * Derived elements will have to - * implement this function. They - * may only provide interpolation - * matrices for certain source - * finite elements, for example - * those from the same family. If - * they don't implement - * interpolation from a given - * element, then they must throw - * an exception of type + * Derived elements will have to implement this function. They may only + * provide interpolation matrices for certain source finite elements, for + * example those from the same family. If they don't implement interpolation + * from a given element, then they must throw an exception of type * ExcInterpolationNotImplemented. */ virtual void @@ -1038,24 +762,14 @@ public: /** - * Return the matrix - * interpolating from a face of - * of one element to the face of - * the neighboring element. - * The size of the matrix is - * then source.#dofs_per_face times - * this->#dofs_per_face. + * Return the matrix interpolating from a face of of one element to the face + * of the neighboring element. The size of the matrix is then + * source.#dofs_per_face times this->#dofs_per_face. * - * Derived elements will have to - * implement this function. They - * may only provide interpolation - * matrices for certain source - * finite elements, for example - * those from the same family. If - * they don't implement - * interpolation from a given - * element, then they must throw - * an exception of type + * Derived elements will have to implement this function. They may only + * provide interpolation matrices for certain source finite elements, for + * example those from the same family. If they don't implement interpolation + * from a given element, then they must throw an exception of type * ExcInterpolationNotImplemented. */ virtual void @@ -1064,24 +778,14 @@ public: /** - * Return the matrix - * interpolating from a face of - * of one element to the subface of - * the neighboring element. - * The size of the matrix is - * then source.#dofs_per_face times - * this->#dofs_per_face. - * - * Derived elements will have to - * implement this function. They - * may only provide interpolation - * matrices for certain source - * finite elements, for example - * those from the same family. If - * they don't implement - * interpolation from a given - * element, then they must throw - * an exception of type + * Return the matrix interpolating from a face of of one element to the + * subface of the neighboring element. The size of the matrix is then + * source.#dofs_per_face times this->#dofs_per_face. + * + * Derived elements will have to implement this function. They may only + * provide interpolation matrices for certain source finite elements, for + * example those from the same family. If they don't implement interpolation + * from a given element, then they must throw an exception of type * ExcInterpolationNotImplemented. */ virtual void @@ -1134,15 +838,11 @@ public: hp_quad_dof_identities (const FiniteElement &fe_other) const; /** - * Return whether this element dominates - * the one given as argument when they - * meet at a common face, - * whether it is the other way around, - * whether neither dominates, or if - * either could dominate. + * Return whether this element dominates the one given as argument when they + * meet at a common face, whether it is the other way around, whether + * neither dominates, or if either could dominate. * - * For a definition of domination, see - * FiniteElementBase::Domination and in + * For a definition of domination, see FiniteElementBase::Domination and in * particular the @ref hp_paper "hp paper". */ virtual @@ -1152,21 +852,13 @@ public: //@} /** - * Comparison operator. We also - * check for equality of the - * constraint matrix, which is - * quite an expensive operation. - * Do therefore use this function - * with care, if possible only - * for debugging purposes. + * Comparison operator. We also check for equality of the constraint matrix, + * which is quite an expensive operation. Do therefore use this function + * with care, if possible only for debugging purposes. * - * Since this function is not - * that important, we avoid an - * implementational question - * about comparing arrays and do - * not compare the matrix arrays - * #restriction and - * #prolongation. + * Since this function is not that important, we avoid an implementational + * question about comparing arrays and do not compare the matrix arrays + * #restriction and #prolongation. */ bool operator == (const FiniteElement &) const; @@ -1175,116 +867,65 @@ public: * @{ */ /** - * Compute vector component and - * index of this shape function - * within the shape functions - * corresponding to this - * component from the index of a - * shape function within this - * finite element. + * Compute vector component and index of this shape function within the + * shape functions corresponding to this component from the index of a shape + * function within this finite element. * - * If the element is scalar, then - * the component is always zero, - * and the index within this - * component is equal to the - * overall index. - * - * If the shape function - * referenced has more than one - * non-zero component, then it - * cannot be associated with one - * vector component, and an - * exception of type - * ExcShapeFunctionNotPrimitive - * will be raised. - * - * Note that if the element is - * composed of other (base) - * elements, and a base element - * has more than one component - * but all its shape functions - * are primitive (i.e. are - * non-zero in only one - * component), then this mapping - * contains valid - * information. However, the - * index of a shape function of - * this element within one - * component (i.e. the second - * number of the respective entry - * of this array) does not - * indicate the index of the - * respective shape function - * within the base element (since - * that has more than one - * vector-component). For this - * information, refer to the - * #system_to_base_table field - * and the - * system_to_base_index() - * function. + * If the element is scalar, then the component is always zero, and the + * index within this component is equal to the overall index. + * + * If the shape function referenced has more than one non-zero component, + * then it cannot be associated with one vector component, and an exception + * of type ExcShapeFunctionNotPrimitive will be raised. + * + * Note that if the element is composed of other (base) elements, and a base + * element has more than one component but all its shape functions are + * primitive (i.e. are non-zero in only one component), then this mapping + * contains valid information. However, the index of a shape function of + * this element within one component (i.e. the second number of the + * respective entry of this array) does not indicate the index of the + * respective shape function within the base element (since that has more + * than one vector-component). For this information, refer to the + * #system_to_base_table field and the system_to_base_index() function. * - * The use of this function is - * explained extensively in the - * step-8 and @ref - * step_20 "step-20" tutorial - * programs as well as in the - * @ref vector_valued module. + * The use of this function is explained extensively in the step-8 and @ref + * step_20 "step-20" tutorial programs as well as in the @ref vector_valued + * module. */ std::pair system_to_component_index (const unsigned int index) const; /** - * Compute the shape function for - * the given vector component and - * index. + * Compute the shape function for the given vector component and index. * - * If the element is scalar, then - * the component must be zero, - * and the index within this - * component is equal to the - * overall index. + * If the element is scalar, then the component must be zero, and the index + * within this component is equal to the overall index. * - * This is the opposite operation - * from the system_to_component_index() + * This is the opposite operation from the system_to_component_index() * function. */ unsigned int component_to_system_index(const unsigned int component, const unsigned int index) const; /** - * Same as - * system_to_component_index(), - * but do it for shape functions - * and their indices on a - * face. The range of allowed - * indices is therefore + * Same as system_to_component_index(), but do it for shape functions and + * their indices on a face. The range of allowed indices is therefore * 0..#dofs_per_face. * - * You will rarely need this - * function in application - * programs, since almost all - * application codes only need to - * deal with cell indices, not - * face indices. The function is - * mainly there for use inside - * the library. + * You will rarely need this function in application programs, since almost + * all application codes only need to deal with cell indices, not face + * indices. The function is mainly there for use inside the library. */ std::pair face_system_to_component_index (const unsigned int index) const; /** - * For faces with non-standard - * face_orientation in 3D, the dofs on - * faces (quads) have to be permuted in - * order to be combined with the correct - * shape functions. Given a local dof @p - * index on a quad, return the local index, - * if the face has non-standard - * face_orientation, face_flip or - * face_rotation. In 2D and 1D there is no - * need for permutation and consequently - * an exception is thrown. + * For faces with non-standard face_orientation in 3D, the dofs on faces + * (quads) have to be permuted in order to be combined with the correct + * shape functions. Given a local dof @p index on a quad, return the local + * index, if the face has non-standard face_orientation, face_flip or + * face_rotation. In 2D and 1D there is no need for permutation and + * consequently an exception is thrown. */ unsigned int adjust_quad_dof_index_for_face_orientation (const unsigned int index, const bool face_orientation, @@ -1292,268 +933,166 @@ public: const bool face_rotation) const; /** - * For lines with non-standard - * line_orientation in 3D, the dofs on - * lines have to be permuted in order to be - * combined with the correct shape - * functions. Given a local dof @p index on - * a line, return the local index, if the - * line has non-standard - * line_orientation. In 2D and 1D there is - * no need for permutation, so the given - * index is simply returned. + * For lines with non-standard line_orientation in 3D, the dofs on lines + * have to be permuted in order to be combined with the correct shape + * functions. Given a local dof @p index on a line, return the local index, + * if the line has non-standard line_orientation. In 2D and 1D there is no + * need for permutation, so the given index is simply returned. */ unsigned int adjust_line_dof_index_for_line_orientation (const unsigned int index, const bool line_orientation) const; /** - * Return in which of the vector - * components of this finite - * element the @p ith shape - * function is non-zero. The - * length of the returned array - * is equal to the number of - * vector components of this - * element. + * Return in which of the vector components of this finite element the @p + * ith shape function is non-zero. The length of the returned array is equal + * to the number of vector components of this element. + * + * For most finite element spaces, the result of this function will be a + * vector with exactly one element being @p true, since for most spaces the + * individual vector components are independent. In that case, the component + * with the single zero is also the first element of what + * system_to_component_index() returns. * - * For most finite element - * spaces, the result of this - * function will be a vector with - * exactly one element being - * @p true, since for most - * spaces the individual vector - * components are independent. In - * that case, the component with - * the single zero is also the - * first element of what - * system_to_component_index() - * returns. - * - * Only for those spaces that couple the - * components, for example to make a - * shape function divergence free, will - * there be more than one @p true entry. - * Elements for which this is true are - * called non-primitive (see - * @ref GlossPrimitive). + * Only for those spaces that couple the components, for example to make a + * shape function divergence free, will there be more than one @p true + * entry. Elements for which this is true are called non-primitive (see + * @ref GlossPrimitive). */ const ComponentMask & get_nonzero_components (const unsigned int i) const; /** - * Return in how many vector - * components the @p ith shape - * function is non-zero. This - * value equals the number of - * entries equal to @p true in - * the result of the - * get_nonzero_components() - * function. + * Return in how many vector components the @p ith shape function is + * non-zero. This value equals the number of entries equal to @p true in the + * result of the get_nonzero_components() function. * - * For most finite element - * spaces, the result will be - * equal to one. It is not equal - * to one only for those ansatz - * spaces for which vector-valued - * shape functions couple the - * individual components, for - * example in order to make them - * divergence-free. + * For most finite element spaces, the result will be equal to one. It is + * not equal to one only for those ansatz spaces for which vector-valued + * shape functions couple the individual components, for example in order to + * make them divergence-free. */ unsigned int n_nonzero_components (const unsigned int i) const; /** - * Return whether the @p ith - * shape function is primitive in - * the sense that the shape - * function is non-zero in only - * one vector - * component. Non-primitive shape - * functions would then, for - * example, be those of - * divergence free ansatz spaces, - * in which the individual vector + * Return whether the @p ith shape function is primitive in the sense that + * the shape function is non-zero in only one vector + * component. Non-primitive shape functions would then, for example, be + * those of divergence free ansatz spaces, in which the individual vector * components are coupled. * - * The result of the function is - * @p true if and only if the - * result of - * n_nonzero_components(i) is - * equal to one. + * The result of the function is @p true if and only if the result of + * n_nonzero_components(i) is equal to one. */ bool is_primitive (const unsigned int i) const; /** - * Import function that is overloaded - * by the one above and would otherwise + * Import function that is overloaded by the one above and would otherwise * be hidden. */ using FiniteElementData::is_primitive; /** - * Number of base elements in a - * mixed discretization. - * - * Note that even for vector - * valued finite elements, the - * number of components needs not - * coincide with the number of - * base elements, since they may - * be reused. For example, if you - * create a FESystem with - * three identical finite element - * classes by using the - * constructor that takes one - * finite element and a - * multiplicity, then the number - * of base elements is still one, - * although the number of - * components of the finite - * element is equal to the - * multiplicity. + * Number of base elements in a mixed discretization. + * + * Note that even for vector valued finite elements, the number of + * components needs not coincide with the number of base elements, since + * they may be reused. For example, if you create a FESystem with three + * identical finite element classes by using the constructor that takes one + * finite element and a multiplicity, then the number of base elements is + * still one, although the number of components of the finite element is + * equal to the multiplicity. */ unsigned int n_base_elements () const; /** - * Access to base element - * objects. If the element is - * atomic, then - * base_element(0) is - * @p this. + * Access to base element objects. If the element is atomic, then + * base_element(0) is @p this. */ virtual const FiniteElement & base_element (const unsigned int index) const; /** - * This index denotes how often - * the base element @p index is - * used in a composed element. If - * the element is atomic, then - * the result is always equal to - * one. See the documentation for - * the n_base_elements() - * function for more details. + * This index denotes how often the base element @p index is used in a + * composed element. If the element is atomic, then the result is always + * equal to one. See the documentation for the n_base_elements() function + * for more details. */ unsigned int element_multiplicity (const unsigned int index) const; /** - * Return for shape function - * @p index the base element it - * belongs to, the number of the - * copy of this base element - * (which is between zero and the - * multiplicity of this element), - * and the index of this shape - * function within this base - * element. + * Return for shape function @p index the base element it belongs to, the + * number of the copy of this base element (which is between zero and the + * multiplicity of this element), and the index of this shape function + * within this base element. * - * If the element is not composed of - * others, then base and instance - * are always zero, and the index - * is equal to the number of the - * shape function. If the element - * is composed of single - * instances of other elements - * (i.e. all with multiplicity - * one) all of which are scalar, - * then base values and dof - * indices within this element - * are equal to the - * #system_to_component_table. It - * differs only in case the - * element is composed of other - * elements and at least one of - * them is vector-valued itself. - * - * This function returns valid - * values also in the case of - * vector-valued - * (i.e. non-primitive) shape - * functions, in contrast to the - * system_to_component_index() - * function. + * If the element is not composed of others, then base and instance are + * always zero, and the index is equal to the number of the shape + * function. If the element is composed of single instances of other + * elements (i.e. all with multiplicity one) all of which are scalar, then + * base values and dof indices within this element are equal to the + * #system_to_component_table. It differs only in case the element is + * composed of other elements and at least one of them is vector-valued + * itself. + * + * This function returns valid values also in the case of vector-valued + * (i.e. non-primitive) shape functions, in contrast to the + * system_to_component_index() function. */ std::pair, unsigned int> system_to_base_index (const unsigned int index) const; /** - * Same as - * system_to_base_index(), but - * for degrees of freedom located - * on a face. The range of allowed - * indices is therefore - * 0..#dofs_per_face. + * Same as system_to_base_index(), but for degrees of freedom located on a + * face. The range of allowed indices is therefore 0..#dofs_per_face. * - * You will rarely need this - * function in application - * programs, since almost all - * application codes only need to - * deal with cell indices, not - * face indices. The function is - * mainly there for use inside - * the library. + * You will rarely need this function in application programs, since almost + * all application codes only need to deal with cell indices, not face + * indices. The function is mainly there for use inside the library. */ std::pair, unsigned int> face_system_to_base_index (const unsigned int index) const; /** - * Given a base element number, - * return the first block of a - * BlockVector it would generate. + * Given a base element number, return the first block of a BlockVector it + * would generate. */ unsigned int first_block_of_base (const unsigned int b) const; /** - * For each vector component, - * return which base - * element implements this - * component and which vector - * component in this base element - * this is. This information is - * only of interest for - * vector-valued finite elements - * which are composed of several - * sub-elements. In that case, - * one may want to obtain - * information about the element - * implementing a certain vector - * component, which can be done - * using this function and the - * FESystem::base_element() - * function. + * For each vector component, return which base element implements this + * component and which vector component in this base element this is. This + * information is only of interest for vector-valued finite elements which + * are composed of several sub-elements. In that case, one may want to + * obtain information about the element implementing a certain vector + * component, which can be done using this function and the + * FESystem::base_element() function. * - * If this is a scalar finite - * element, then the return value - * is always equal to a pair of - * zeros. + * If this is a scalar finite element, then the return value is always equal + * to a pair of zeros. */ std::pair component_to_base_index (const unsigned int component) const; /** - * Return the base element for - * this block and the number of - * the copy of the base element. + * Return the base element for this block and the number of the copy of the + * base element. */ std::pair block_to_base_index (const unsigned int block) const; /** - * The vector block and the index - * inside the block for this - * shape function. + * The vector block and the index inside the block for this shape function. */ std::pair system_to_block_index (const unsigned int component) const; /** - * The vector block for this - * component. + * The vector block for this component. */ unsigned int component_to_block_index (const unsigned int component) const; @@ -1566,57 +1105,53 @@ public: */ /** - * Return a component mask with as many elements as this - * object has vector components and of which exactly the - * one component is true that corresponds to the given - * argument. See @ref GlossComponentMask "the glossary" - * for more information. + * Return a component mask with as many elements as this object has vector + * components and of which exactly the one component is true that + * corresponds to the given argument. See @ref GlossComponentMask "the + * glossary" for more information. * - * @param scalar An object that represents a single scalar - * vector component of this finite element. - * @return A component mask that is false in all components - * except for the one that corresponds to the argument. + * @param scalar An object that represents a single scalar vector component + * of this finite element. + * @return A component mask that is false in all components except for the + * one that corresponds to the argument. */ ComponentMask component_mask (const FEValuesExtractors::Scalar &scalar) const; /** - * Return a component mask with as many elements as this - * object has vector components and of which exactly the - * dim components are true that correspond to the given - * argument. See @ref GlossComponentMask "the glossary" - * for more information. + * Return a component mask with as many elements as this object has vector + * components and of which exactly the dim components are true + * that correspond to the given argument. See @ref GlossComponentMask "the + * glossary" for more information. * - * @param vector An object that represents dim - * vector components of this finite element. - * @return A component mask that is false in all components - * except for the ones that corresponds to the argument. + * @param vector An object that represents dim vector components of this + * finite element. + * @return A component mask that is false in all components except for the + * ones that corresponds to the argument. */ ComponentMask component_mask (const FEValuesExtractors::Vector &vector) const; /** - * Return a component mask with as many elements as this - * object has vector components and of which exactly the - * dim*(dim+1)/2 components are true that - * correspond to the given argument. See @ref GlossComponentMask "the glossary" - * for more information. + * Return a component mask with as many elements as this object has vector + * components and of which exactly the dim*(dim+1)/2 components + * are true that correspond to the given argument. See @ref + * GlossComponentMask "the glossary" for more information. * - * @param sym_tensor An object that represents dim*(dim+1)/2 - * components of this finite element that are jointly to be - * interpreted as forming a symmetric tensor. - * @return A component mask that is false in all components - * except for the ones that corresponds to the argument. + * @param sym_tensor An object that represents dim*(dim+1)/2 components of + * this finite element that are jointly to be interpreted as forming a + * symmetric tensor. @return A component mask that is false in all + * components except for the ones that corresponds to the argument. */ ComponentMask component_mask (const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const; /** * Given a block mask (see @ref GlossBlockMask "this glossary entry"), - * produce a component mask (see @ref GlossComponentMask "this glossary entry") - * that represents the components that correspond to the blocks selected in - * the input argument. This is essentially a conversion operator from - * BlockMask to ComponentMask. + * produce a component mask (see @ref GlossComponentMask "this glossary + * entry") that represents the components that correspond to the blocks + * selected in the input argument. This is essentially a conversion operator + * from BlockMask to ComponentMask. * * @param block_mask The mask that selects individual blocks of the finite * element @@ -1627,84 +1162,80 @@ public: component_mask (const BlockMask &block_mask) const; /** - * Return a block mask with as many elements as this - * object has blocks and of which exactly the - * one component is true that corresponds to the given - * argument. See @ref GlossBlockMask "the glossary" - * for more information. + * Return a block mask with as many elements as this object has blocks and + * of which exactly the one component is true that corresponds to the given + * argument. See @ref GlossBlockMask "the glossary" for more information. * - * @note This function will only succeed if the scalar referenced - * by the argument encompasses a complete block. In other words, - * if, for example, you pass an extractor for the single - * $x$ velocity and this object represents an FE_RaviartThomas - * object, then the single scalar object you selected is part - * of a larger block and consequently there is no block mask that - * would represent it. The function will then produce an exception. + * @note This function will only succeed if the scalar referenced by the + * argument encompasses a complete block. In other words, if, for example, + * you pass an extractor for the single $x$ velocity and this object + * represents an FE_RaviartThomas object, then the single scalar object you + * selected is part of a larger block and consequently there is no block + * mask that would represent it. The function will then produce an + * exception. * - * @param scalar An object that represents a single scalar - * vector component of this finite element. - * @return A component mask that is false in all components - * except for the one that corresponds to the argument. + * @param scalar An object that represents a single scalar vector component + * of this finite element. + * @return A component mask that is false in all components except for the + * one that corresponds to the argument. */ BlockMask block_mask (const FEValuesExtractors::Scalar &scalar) const; /** - * Return a component mask with as many elements as this - * object has vector components and of which exactly the - * dim components are true that correspond to the given - * argument. See @ref GlossBlockMask "the glossary" - * for more information. + * Return a component mask with as many elements as this object has vector + * components and of which exactly the dim components are true + * that correspond to the given argument. See @ref GlossBlockMask "the + * glossary" for more information. * * @note The same caveat applies as to the version of the function above: - * The extractor object passed as argument must be so that it corresponds - * to full blocks and does not split blocks of this element. + * The extractor object passed as argument must be so that it corresponds to + * full blocks and does not split blocks of this element. * - * @param vector An object that represents dim - * vector components of this finite element. - * @return A component mask that is false in all components - * except for the ones that corresponds to the argument. + * @param vector An object that represents dim vector components of this + * finite element. + * @return A component mask that is false in all components except for the + * ones that corresponds to the argument. */ BlockMask block_mask (const FEValuesExtractors::Vector &vector) const; /** - * Return a component mask with as many elements as this - * object has vector components and of which exactly the - * dim*(dim+1)/2 components are true that - * correspond to the given argument. See @ref GlossBlockMask "the glossary" - * for more information. + * Return a component mask with as many elements as this object has vector + * components and of which exactly the dim*(dim+1)/2 components + * are true that correspond to the given argument. See @ref GlossBlockMask + * "the glossary" for more information. * * @note The same caveat applies as to the version of the function above: - * The extractor object passed as argument must be so that it corresponds - * to full blocks and does not split blocks of this element. + * The extractor object passed as argument must be so that it corresponds to + * full blocks and does not split blocks of this element. * - * @param sym_tensor An object that represents dim*(dim+1)/2 - * components of this finite element that are jointly to be - * interpreted as forming a symmetric tensor. - * @return A component mask that is false in all components - * except for the ones that corresponds to the argument. + * @param sym_tensor An object that represents dim*(dim+1)/2 components of + * this finite element that are jointly to be interpreted as forming a + * symmetric tensor. + * @return A component mask that is false in all components except for the + * ones that corresponds to the argument. */ BlockMask block_mask (const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const; /** - * Given a component mask (see @ref GlossComponentMask "this glossary entry"), - * produce a block mask (see @ref GlossBlockMask "this glossary entry") - * that represents the blocks that correspond to the components selected in - * the input argument. This is essentially a conversion operator from - * ComponentMask to BlockMask. + * Given a component mask (see @ref GlossComponentMask "this glossary + * entry"), produce a block mask (see @ref GlossBlockMask "this glossary + * entry") that represents the blocks that correspond to the components + * selected in the input argument. This is essentially a conversion operator + * from ComponentMask to BlockMask. * - * @note This function will only succeed if the components referenced - * by the argument encompasses complete blocks. In other words, - * if, for example, you pass an component mask for the single - * $x$ velocity and this object represents an FE_RaviartThomas - * object, then the single component you selected is part - * of a larger block and consequently there is no block mask that - * would represent it. The function will then produce an exception. + * @note This function will only succeed if the components referenced by the + * argument encompasses complete blocks. In other words, if, for example, + * you pass an component mask for the single $x$ velocity and this object + * represents an FE_RaviartThomas object, then the single component you + * selected is part of a larger block and consequently there is no block + * mask that would represent it. The function will then produce an + * exception. * - * @param component_mask The mask that selects individual components of the finite - * element + * @param component_mask The mask that selects individual components of the + * finite element * @return A mask that selects those blocks corresponding to the selected * blocks of the input argument. */ @@ -1719,199 +1250,118 @@ public: */ /** - * Return the support points of - * the trial functions on the - * unit cell, if the derived - * finite element defines some. - * Finite elements that allow - * some kind of interpolation - * operation usually have support - * points. On the other hand, - * elements that define their - * degrees of freedom by, for - * example, moments on faces, or - * as derivatives, don't have - * support points. In that case, - * the returned field is empty. - * - * If the finite element defines - * support points, then their - * number equals the number of - * degrees of freedom of the - * element. The order of points - * in the array matches that - * returned by the - * cell->get_dof_indices + * Return the support points of the trial functions on the unit cell, if the + * derived finite element defines some. Finite elements that allow some + * kind of interpolation operation usually have support points. On the other + * hand, elements that define their degrees of freedom by, for example, + * moments on faces, or as derivatives, don't have support points. In that + * case, the returned field is empty. + * + * If the finite element defines support points, then their number equals + * the number of degrees of freedom of the element. The order of points in + * the array matches that returned by the cell->get_dof_indices * function. * - * See the class documentation - * for details on support points. + * See the class documentation for details on support points. */ const std::vector > & get_unit_support_points () const; /** - * Return whether a finite - * element has defined support - * points. If the result is true, - * then a call to the - * get_unit_support_points() - * yields a non-empty array. + * Return whether a finite element has defined support points. If the result + * is true, then a call to the get_unit_support_points() yields a non-empty + * array. * - * The result may be false if an - * element is not defined by - * interpolating shape functions, - * for example by P-elements on - * quadrilaterals. It will - * usually only be true if the - * element constructs its shape - * functions by the requirement - * that they be one at a certain - * point and zero at all the - * points associated with the - * other shape functions. - * - * In composed elements (i.e. for - * the FESystem class, the - * result will be true if all all - * the base elements have defined - * support points. + * The result may be false if an element is not defined by interpolating + * shape functions, for example by P-elements on quadrilaterals. It will + * usually only be true if the element constructs its shape functions by the + * requirement that they be one at a certain point and zero at all the + * points associated with the other shape functions. + * + * In composed elements (i.e. for the FESystem class, the result will be + * true if all all the base elements have defined support points. */ bool has_support_points () const; /** - * Return the position of the - * support point of the - * @p indexth shape function. If - * it does not exist, raise an - * exception. + * Return the position of the support point of the @p indexth shape + * function. If it does not exist, raise an exception. * - * The default implementation - * simply returns the respective - * element from the array you get - * from - * get_unit_support_points(), - * but derived elements may - * overload this function. In - * particular, note that the - * FESystem class overloads - * it so that it can return the - * support points of individual - * base elements, if not all the - * base elements define support - * points. In this way, you can - * still ask for certain support - * points, even if - * get_unit_support_points() - * only returns an empty array. + * The default implementation simply returns the respective element from the + * array you get from get_unit_support_points(), but derived elements may + * overload this function. In particular, note that the FESystem class + * overloads it so that it can return the support points of individual base + * elements, if not all the base elements define support points. In this + * way, you can still ask for certain support points, even if + * get_unit_support_points() only returns an empty array. */ virtual Point unit_support_point (const unsigned int index) const; /** - * Return the support points of - * the trial functions on the - * unit face, if the derived - * finite element defines some. - * Finite elements that allow - * some kind of interpolation - * operation usually have support - * points. On the other hand, - * elements that define their - * degrees of freedom by, for - * example, moments on faces, or - * as derivatives, don't have - * support points. In that case, - * the returned field is empty - * - * Note that elements that have - * support points need not - * necessarily have some on the - * faces, even if the - * interpolation points are - * located physically on a - * face. For example, the - * discontinuous elements have - * interpolation points on the - * vertices, and for higher - * degree elements also on the - * faces, but they are not - * defined to be on faces since - * in that case degrees of - * freedom from both sides of a - * face (or from all adjacent - * elements to a vertex) would be - * identified with each other, - * which is not what we would - * like to have). Logically, - * these degrees of freedom are - * therefore defined to belong to - * the cell, rather than the face - * or vertex. In that case, the - * returned element would - * therefore have length zero. - * - * If the finite element defines - * support points, then their - * number equals the number of - * degrees of freedom on the face - * (#dofs_per_face). The order - * of points in the array matches - * that returned by the - * cell->get_dof_indices - * function. - * - * See the class documentation - * for details on support points. + * Return the support points of the trial functions on the unit face, if the + * derived finite element defines some. Finite elements that allow some + * kind of interpolation operation usually have support points. On the other + * hand, elements that define their degrees of freedom by, for example, + * moments on faces, or as derivatives, don't have support points. In that + * case, the returned field is empty + * + * Note that elements that have support points need not necessarily have + * some on the faces, even if the interpolation points are located + * physically on a face. For example, the discontinuous elements have + * interpolation points on the vertices, and for higher degree elements also + * on the faces, but they are not defined to be on faces since in that case + * degrees of freedom from both sides of a face (or from all adjacent + * elements to a vertex) would be identified with each other, which is not + * what we would like to have). Logically, these degrees of freedom are + * therefore defined to belong to the cell, rather than the face or + * vertex. In that case, the returned element would therefore have length + * zero. + * + * If the finite element defines support points, then their number equals + * the number of degrees of freedom on the face (#dofs_per_face). The order + * of points in the array matches that returned by the + * cell->get_dof_indices function. + * + * See the class documentation for details on support points. */ const std::vector > & get_unit_face_support_points () const; /** - * Return whether a finite - * element has defined support - * points on faces. If the result - * is true, then a call to the - * get_unit_face_support_points() + * Return whether a finite element has defined support points on faces. If + * the result is true, then a call to the get_unit_face_support_points() * yields a non-empty array. * - * For more information, see the - * documentation for the - * has_support_points() + * For more information, see the documentation for the has_support_points() * function. */ bool has_face_support_points () const; /** - * The function corresponding to - * the unit_support_point() - * function, but for faces. See - * there for more information. + * The function corresponding to the unit_support_point() function, but for + * faces. See there for more information. */ virtual Point unit_face_support_point (const unsigned int index) const; /** - * Return a support point vector - * for generalized interpolation. - * - * See the @ref GlossGeneralizedSupport "glossary entry on generalized support points" - * for more information. + * Return a support point vector for generalized interpolation. + * + * See the @ref GlossGeneralizedSupport "glossary entry on generalized + * support points" for more information. */ const std::vector > & get_generalized_support_points () const; /** - * Returns true if the - * class provides nonempty - * vectors either from - * get_unit_support_points() or - * get_generalized_support_points(). - * - * See the @ref GlossGeneralizedSupport "glossary entry on generalized support points" - * for more information. + * Returns true if the class provides nonempty vectors either from + * get_unit_support_points() or get_generalized_support_points(). + * + * See the @ref GlossGeneralizedSupport "glossary entry on generalized + * support points" for more information. */ bool has_generalized_support_points () const; @@ -1922,61 +1372,43 @@ public: get_generalized_face_support_points () const; /** - * Return whether a finite - * element has defined - * generalized support - * points on faces. If the result - * is true, then a call to the - * get_generalized_face_support_points - * yields a non-empty array. + * Return whether a finite element has defined generalized support points on + * faces. If the result is true, then a call to the + * get_generalized_face_support_points yields a non-empty array. * - * For more information, see the - * documentation for the - * has_support_points() + * For more information, see the documentation for the has_support_points() * function. */ bool has_generalized_face_support_points () const; /** - * Interpolate a set of scalar - * values, computed in the - * generalized support points. + * Interpolate a set of scalar values, computed in the generalized support + * points. * - * @note This function is - * implemented in - * FiniteElement for the case - * that the element has support - * points. In this case, the - * resulting coefficients are - * just the values in the suport - * points. All other elements - * must reimplement it. + * @note This function is implemented in FiniteElement for the case that the + * element has support points. In this case, the resulting coefficients are + * just the values in the suport points. All other elements must reimplement + * it. */ virtual void interpolate(std::vector &local_dofs, const std::vector &values) const; /** - * Interpolate a set of vector - * values, computed in the - * generalized support points. + * Interpolate a set of vector values, computed in the generalized support + * points. * - * Since a finite element often - * only interpolates part of a - * vector, offset is - * used to determine the first - * component of the vector to be - * interpolated. Maybe consider - * changing your data structures - * to use the next function. + * Since a finite element often only interpolates part of a vector, + * offset is used to determine the first component of the vector to + * be interpolated. Maybe consider changing your data structures to use the + * next function. */ virtual void interpolate(std::vector &local_dofs, const std::vector > &values, unsigned int offset = 0) const; /** - * Interpolate a set of vector - * values, computed in the - * generalized support points. + * Interpolate a set of vector values, computed in the generalized support + * points. */ virtual void interpolate( std::vector &local_dofs, @@ -1985,15 +1417,12 @@ public: //@} /** - * Determine an estimate for the - * memory consumption (in bytes) - * of this object. + * Determine an estimate for the memory consumption (in bytes) of this + * object. * - * This function is made virtual, - * since finite element objects - * are usually accessed through - * pointers to their base class, - * rather than the class itself. + * This function is made virtual, since finite element objects are usually + * accessed through pointers to their base class, rather than the class + * itself. */ virtual std::size_t memory_consumption () const; /** @@ -2023,29 +1452,24 @@ public: DeclException0 (ExcUnitShapeValuesDoNotExist); /** - * Attempt to access support - * points of a finite element - * which is not Lagrangian. + * Attempt to access support points of a finite element which is not + * Lagrangian. * * @ingroup Exceptions */ DeclException0 (ExcFEHasNoSupportPoints); /** - * Attempt to access embedding - * matrices of a finite element - * which did not implement these - * matrices. + * Attempt to access embedding matrices of a finite element which did not + * implement these matrices. * * @ingroup Exceptions */ DeclException0 (ExcEmbeddingVoid); /** - * Attempt to access restriction - * matrices of a finite element - * which did not implement these - * matrices. + * Attempt to access restriction matrices of a finite element which did not + * implement these matrices. * * Exception * @ingroup Exceptions @@ -2053,10 +1477,8 @@ public: DeclException0 (ExcProjectionVoid); /** - * Attempt to access constraint - * matrices of a finite element - * which did not implement these - * matrices. + * Attempt to access constraint matrices of a finite element which did not + * implement these matrices. * * Exception * @ingroup Exceptions @@ -2102,235 +1524,140 @@ public: protected: /** - * Reinit the vectors of - * restriction and prolongation - * matrices to the right sizes: - * For every refinement case, - * except for - * RefinementCase::no_refinement, - * and for every child of that - * refinement case the space of - * one restriction and - * prolongation matrix is - * allocated, see the - * documentation of the - * restriction and prolongation - * vectors for more detail on the - * actual vector sizes. - * - * @param - * isotropic_restriction_only - * only the restriction matrices - * required for isotropic - * refinement are reinited to the - * right size. - * @param - * isotropic_prolongation_only - * only the prolongation matrices - * required for isotropic - * refinement are reinited to the - * right size. + * Reinit the vectors of restriction and prolongation matrices to the right + * sizes: For every refinement case, except for + * RefinementCase::no_refinement, and for every child of that refinement + * case the space of one restriction and prolongation matrix is allocated, + * see the documentation of the restriction and prolongation vectors for + * more detail on the actual vector sizes. + * + * @param isotropic_restriction_only only the restriction matrices required + * for isotropic refinement are reinited to the right size. + * @param isotropic_prolongation_only only the prolongation matrices + * required for isotropic refinement are reinited to the right size. */ void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false); /** - * Vector of projection - * matrices. See - * get_restriction_matrix() - * above. The constructor - * initializes these matrices to - * zero dimensions, which can be - * changed by derived classes - * implementing them. - * - * Note, that - * restriction[refinement_case-1][child] - * includes the restriction - * matrix of child - * child for the - * RefinementCase - * refinement_case. Here, - * we use - * refinement_case-1 - * instead of - * refinement_case - * as for - * RefinementCase::no_refinement(=0) - * there are no restriction - * matrices available. + * Vector of projection matrices. See get_restriction_matrix() above. The + * constructor initializes these matrices to zero dimensions, which can be + * changed by derived classes implementing them. + * + * Note, that restriction[refinement_case-1][child] includes + * the restriction matrix of child child for the RefinementCase + * refinement_case. Here, we use refinement_case-1 + * instead of refinement_case as for + * RefinementCase::no_refinement(=0) there are no restriction matrices + * available. */ std::vector > > restriction; /** - * Vector of embedding - * matrices. See - * get_prolongation_matrix() - * above. The constructor - * initializes these matrices to - * zero dimensions, which can be - * changed by derived classes - * implementing them. - * - * Note, that - * prolongation[refinement_case-1][child] - * includes the prolongation - * matrix of child - * child for the - * RefinementCase - * refinement_case. Here, - * we use - * refinement_case-1 - * instead of - * refinement_case - * as for - * RefinementCase::no_refinement(=0) - * there are no prolongation - * matrices available. + * Vector of embedding matrices. See get_prolongation_matrix() + * above. The constructor initializes these matrices to zero dimensions, + * which can be changed by derived classes implementing them. + * + * Note, that prolongation[refinement_case-1][child] includes + * the prolongation matrix of child child for the + * RefinementCase refinement_case. Here, we use + * refinement_case-1 instead of refinement_case as + * for RefinementCase::no_refinement(=0) there are no prolongation matrices + * available. */ std::vector > > prolongation; /** - * Specify the constraints which - * the dofs on the two sides of a - * cell interface underly if the - * line connects two cells of - * which one is refined once. + * Specify the constraints which the dofs on the two sides of a cell + * interface underly if the line connects two cells of which one is refined + * once. * - * For further details see the - * general description of the - * derived class. + * For further details see the general description of the derived class. * - * This field is obviously - * useless in one dimension - * and has there a zero size. + * This field is obviously useless in one dimension and has there a zero + * size. */ FullMatrix interface_constraints; /** - * List of support points on the - * unit cell, in case the finite - * element has any. The - * constructor leaves this field - * empty, derived classes may - * write in some contents. - * - * Finite elements that allow - * some kind of interpolation - * operation usually have support - * points. On the other hand, - * elements that define their - * degrees of freedom by, for - * example, moments on faces, or - * as derivatives, don't have - * support points. In that case, - * this field remains empty. + * List of support points on the unit cell, in case the finite element has + * any. The constructor leaves this field empty, derived classes may write + * in some contents. + * + * Finite elements that allow some kind of interpolation operation usually + * have support points. On the other hand, elements that define their + * degrees of freedom by, for example, moments on faces, or as derivatives, + * don't have support points. In that case, this field remains empty. */ std::vector > unit_support_points; /** - * Same for the faces. See the - * description of the - * get_unit_face_support_points() - * function for a discussion of - * what contributes a face - * support point. + * Same for the faces. See the description of the + * get_unit_face_support_points() function for a discussion of what + * contributes a face support point. */ std::vector > unit_face_support_points; /** - * Support points used for - * interpolation functions of - * non-Lagrangian elements. + * Support points used for interpolation functions of non-Lagrangian + * elements. */ std::vector > generalized_support_points; /** - * Face support points used for - * interpolation functions of - * non-Lagrangian elements. + * Face support points used for interpolation functions of non-Lagrangian + * elements. */ std::vector > generalized_face_support_points; /** - * For faces with non-standard - * face_orientation in 3D, the dofs on - * faces (quads) have to be permuted in - * order to be combined with the correct - * shape functions. Given a local dof @p - * index on a quad, return the shift in the - * local index, if the face has - * non-standard face_orientation, - * i.e. old_index + shift = - * new_index. In 2D and 1D there is - * no need for permutation so the vector is - * empty. In 3D it has the size of - * #dofs_per_quad * 8 , where 8 is - * the number of orientations, a face can - * be in (all combinations of the three - * bool flags face_orientation, face_flip - * and face_rotation). - * - * The standard implementation fills this - * with zeros, i.e. no permuatation at - * all. Derived finite element classes have - * to fill this Table with the correct - * values. + * For faces with non-standard face_orientation in 3D, the dofs on faces + * (quads) have to be permuted in order to be combined with the correct + * shape functions. Given a local dof @p index on a quad, return the shift + * in the local index, if the face has non-standard face_orientation, + * i.e. old_index + shift = new_index. In 2D and 1D there is no + * need for permutation so the vector is empty. In 3D it has the size of + * #dofs_per_quad * 8 , where 8 is the number of orientations, + * a face can be in (all combinations of the three bool flags + * face_orientation, face_flip and face_rotation). + * + * The standard implementation fills this with zeros, i.e. no permuatation + * at all. Derived finite element classes have to fill this Table with the + * correct values. */ Table<2,int> adjust_quad_dof_index_for_face_orientation_table; /** - * For lines with non-standard - * line_orientation in 3D, the dofs on - * lines have to be permuted in - * order to be combined with the correct - * shape functions. Given a local dof @p - * index on a line, return the shift in the - * local index, if the line has - * non-standard line_orientation, - * i.e. old_index + shift = - * new_index. In 2D and 1D there is - * no need for permutation so the vector is - * empty. In 3D it has the size of + * For lines with non-standard line_orientation in 3D, the dofs on lines + * have to be permuted in order to be combined with the correct shape + * functions. Given a local dof @p index on a line, return the shift in the + * local index, if the line has non-standard line_orientation, + * i.e. old_index + shift = new_index. In 2D and 1D there is no + * need for permutation so the vector is empty. In 3D it has the size of * #dofs_per_line. * - * The standard implementation fills this - * with zeros, i.e. no permutation at - * all. Derived finite element classes have - * to fill this vector with the correct - * values. + * The standard implementation fills this with zeros, i.e. no permutation at + * all. Derived finite element classes have to fill this vector with the + * correct values. */ std::vector adjust_line_dof_index_for_line_orientation_table; /** - * Return the size of interface - * constraint matrices. Since - * this is needed in every - * derived finite element class - * when initializing their size, - * it is placed into this - * function, to avoid having to - * recompute the - * dimension-dependent size of - * these matrices each time. - * - * Note that some elements do not - * implement the interface - * constraints for certain - * polynomial degrees. In this - * case, this function still - * returns the size these - * matrices should have when - * implemented, but the actual - * matrices are empty. + * Return the size of interface constraint matrices. Since this is needed in + * every derived finite element class when initializing their size, it is + * placed into this function, to avoid having to recompute the + * dimension-dependent size of these matrices each time. + * + * Note that some elements do not implement the interface constraints for + * certain polynomial degrees. In this case, this function still returns the + * size these matrices should have when implemented, but the actual matrices + * are empty. */ TableIndices<2> interface_constraints_size () const; /** - * Compute second derivatives by - * finite differences of - * gradients. + * Compute second derivatives by finite differences of gradients. */ void compute_2nd (const Mapping &mapping, const typename Triangulation::cell_iterator &cell, @@ -2340,174 +1667,107 @@ protected: FEValuesData &data) const; /** - * Given the pattern of nonzero - * components for each shape - * function, compute for each - * entry how many components are - * non-zero for each shape - * function. This function is - * used in the constructor of - * this class. + * Given the pattern of nonzero components for each shape function, compute + * for each entry how many components are non-zero for each shape + * function. This function is used in the constructor of this class. */ static std::vector compute_n_nonzero_components (const std::vector &nonzero_components); /** - * Determine the values a finite - * element should compute on - * initialization of data for - * FEValues. - * - * Given a set of flags - * indicating what quantities are - * requested from a FEValues - * object, update_once() and - * update_each() compute which - * values must really be - * computed. Then, the - * fill_*_values functions - * are called with the result of - * these. - * - * Furthermore, values must be - * computed either on the unit - * cell or on the physical - * cell. For instance, the - * function values of FE_Q do - * only depend on the quadrature - * points on the unit - * cell. Therefore, this flags - * will be returned by - * update_once(). The gradients - * require computation of the - * covariant transformation - * matrix. Therefore, - * @p update_covariant_transformation - * and @p update_gradients will - * be returned by - * update_each(). - * - * For an example see the same - * function in the derived class - * FE_Q. + * Determine the values a finite element should compute on initialization of + * data for FEValues. + * + * Given a set of flags indicating what quantities are requested from a + * FEValues object, update_once() and update_each() compute which values + * must really be computed. Then, the fill_*_values functions are + * called with the result of these. + * + * Furthermore, values must be computed either on the unit cell or on the + * physical cell. For instance, the function values of FE_Q do only depend + * on the quadrature points on the unit cell. Therefore, this flags will be + * returned by update_once(). The gradients require computation of the + * covariant transformation matrix. Therefore, @p + * update_covariant_transformation and @p update_gradients will be returned + * by update_each(). + * + * For an example see the same function in the derived class FE_Q. */ virtual UpdateFlags update_once (const UpdateFlags flags) const = 0; /** - * Complementary function for - * update_once(). + * Complementary function for update_once(). * - * While update_once() returns - * the values to be computed on - * the unit cell for yielding the - * required data, this function - * determines the values that - * must be recomputed on each - * cell. + * While update_once() returns the values to be computed on the unit cell + * for yielding the required data, this function determines the values that + * must be recomputed on each cell. * - * Refer to update_once() for - * more details. + * Refer to update_once() for more details. */ virtual UpdateFlags update_each (const UpdateFlags flags) const = 0; /** - * A sort of virtual copy - * constructor. Some places in - * the library, for example the - * constructors of FESystem as - * well as the hp::FECollection - * class, need to make copies of - * finite elements without - * knowing their exact type. They - * do so through this function. + * A sort of virtual copy constructor. Some places in the library, for + * example the constructors of FESystem as well as the hp::FECollection + * class, need to make copies of finite elements without knowing their exact + * type. They do so through this function. */ virtual FiniteElement *clone() const = 0; private: /** - * Store what - * system_to_component_index() - * will return. + * Store what system_to_component_index() will return. */ std::vector< std::pair > system_to_component_table; /** - * Map between linear dofs and - * component dofs on face. This - * is filled with default values - * in the constructor, but - * derived classes will have to - * overwrite the information if - * necessary. + * Map between linear dofs and component dofs on face. This is filled with + * default values in the constructor, but derived classes will have to + * overwrite the information if necessary. * - * By component, we mean the - * vector component, not the base - * element. The information thus - * makes only sense if a shape - * function is non-zero in only + * By component, we mean the vector component, not the base element. The + * information thus makes only sense if a shape function is non-zero in only * one component. */ std::vector< std::pair > face_system_to_component_table; /** - * For each shape function, store - * to which base element and - * which instance of this base - * element (in case its - * multiplicity is greater than - * one) it belongs, and its index - * within this base element. If - * the element is not composed of - * others, then base and instance - * are always zero, and the index - * is equal to the number of the - * shape function. If the element - * is composed of single - * instances of other elements - * (i.e. all with multiplicity - * one) all of which are scalar, - * then base values and dof - * indices within this element - * are equal to the - * #system_to_component_table. It - * differs only in case the - * element is composed of other - * elements and at least one of - * them is vector-valued itself. - * - * This array has valid values - * also in the case of - * vector-valued - * (i.e. non-primitive) shape - * functions, in contrast to the + * For each shape function, store to which base element and which instance + * of this base element (in case its multiplicity is greater than one) it + * belongs, and its index within this base element. If the element is not + * composed of others, then base and instance are always zero, and the index + * is equal to the number of the shape function. If the element is composed + * of single instances of other elements (i.e. all with multiplicity one) + * all of which are scalar, then base values and dof indices within this + * element are equal to the #system_to_component_table. It differs only in + * case the element is composed of other elements and at least one of them + * is vector-valued itself. + * + * This array has valid values also in the case of vector-valued + * (i.e. non-primitive) shape functions, in contrast to the * #system_to_component_table. */ std::vector,unsigned int> > system_to_base_table; /** - * Likewise for the indices on - * faces. + * Likewise for the indices on faces. */ std::vector,unsigned int> > face_system_to_base_table; /** - * For each base element, store - * the number of blocks generated - * by the base and the first block in a block - * vector it will generate. + * For each base element, store the number of blocks generated by the base + * and the first block in a block vector it will generate. */ BlockIndices base_to_block_indices; /** - * The base element establishing - * a component. + * The base element establishing a component. * - * For each component number - * c, the entries have - * the following meaning: + * For each component number c, the entries have the following + * meaning: *
*
table[c].first.first
*
Number of the base element for c.
@@ -2517,147 +1777,82 @@ private: *
Multiple of the base element for c.
*
* - * This variable is set to the - * correct size by the - * constructor of this class, but - * needs to be initialized by - * derived classes, unless its - * size is one and the only entry - * is a zero, which is the case - * for scalar elements. In that - * case, the initialization by - * the base class is sufficient. + * This variable is set to the correct size by the constructor of this + * class, but needs to be initialized by derived classes, unless its size is + * one and the only entry is a zero, which is the case for scalar + * elements. In that case, the initialization by the base class is + * sufficient. */ std::vector, unsigned int> > component_to_base_table; /** - * Projection matrices are - * concatenated or summed up. - * - * This flags decides on how the - * projection matrices of the - * children of the same father - * are put together to one - * operator. The possible modes - * are concatenation and - * summation. - * - * If the projection is defined - * by an interpolation operator, - * the child matrices are - * concatenated, i.e. values - * belonging to the same node - * functional are identified and - * enter the interpolated value - * only once. In this case, the - * flag must be @p false. - * - * For projections with respect - * to scalar products, the child - * matrices must be summed up to - * build the complete matrix. The - * flag should be @p true. - * - * For examples of use of these - * flags, see the places in the - * library where it is queried. - * - * There is one flag per shape - * function, indicating whether - * it belongs to the class of - * shape functions that are - * additive in the restriction or - * not. - * - * Note that in previous versions - * of the library, there was one - * flag per vector component of - * the element. This is based on - * the fact that all the shape - * functions that belong to the - * same vector component must - * necessarily behave in the same - * way, to make things - * reasonable. However, the - * problem is that it is - * sometimes impossible to query - * this flag in the vector-valued - * case: this used to be done - * with the - * #system_to_component_index - * function that returns which - * vector component a shape - * function is associated - * with. The point is that since - * we now support shape functions - * that are associated with more - * than one vector component (for - * example the shape functions of - * Raviart-Thomas, or Nedelec - * elements), that function can - * no more be used, so it can be - * difficult to find out which - * for vector component we would - * like to query the + * Projection matrices are concatenated or summed up. + * + * This flags decides on how the projection matrices of the children of the + * same father are put together to one operator. The possible modes are + * concatenation and summation. + * + * If the projection is defined by an interpolation operator, the child + * matrices are concatenated, i.e. values belonging to the same node + * functional are identified and enter the interpolated value only once. In + * this case, the flag must be @p false. + * + * For projections with respect to scalar products, the child matrices must + * be summed up to build the complete matrix. The flag should be @p true. + * + * For examples of use of these flags, see the places in the library where + * it is queried. + * + * There is one flag per shape function, indicating whether it belongs to + * the class of shape functions that are additive in the restriction or not. + * + * Note that in previous versions of the library, there was one flag per + * vector component of the element. This is based on the fact that all the + * shape functions that belong to the same vector component must necessarily + * behave in the same way, to make things reasonable. However, the problem + * is that it is sometimes impossible to query this flag in the + * vector-valued case: this used to be done with the + * #system_to_component_index function that returns which vector component a + * shape function is associated with. The point is that since we now support + * shape functions that are associated with more than one vector component + * (for example the shape functions of Raviart-Thomas, or Nedelec elements), + * that function can no more be used, so it can be difficult to find out + * which for vector component we would like to query the * restriction-is-additive flags. */ const std::vector restriction_is_additive_flags; /** - * For each shape function, give - * a vector of bools (with size - * equal to the number of vector - * components which this finite - * element has) indicating in - * which component each of these - * shape functions is non-zero. + * For each shape function, give a vector of bools (with size equal to the + * number of vector components which this finite element has) indicating in + * which component each of these shape functions is non-zero. * - * For primitive elements, there - * is only one non-zero - * component. + * For primitive elements, there is only one non-zero component. */ const std::vector nonzero_components; /** - * This array holds how many - * values in the respective entry - * of the #nonzero_components - * element are non-zero. The - * array is thus a short-cut to - * allow faster access to this - * information than if we had to - * count the non-zero entries - * upon each request for this - * information. The field is - * initialized in the constructor - * of this class. + * This array holds how many values in the respective entry of the + * #nonzero_components element are non-zero. The array is thus a short-cut + * to allow faster access to this information than if we had to count the + * non-zero entries upon each request for this information. The field is + * initialized in the constructor of this class. */ const std::vector n_nonzero_components_table; /** - * Second derivatives of shapes - * functions are not computed - * analytically, but by finite - * differences of the - * gradients. This static - * variable denotes the step - * length to be used for - * that. It's value is set to - * 1e-6. + * Second derivatives of shapes functions are not computed analytically, but + * by finite differences of the gradients. This static variable denotes the + * step length to be used for that. It's value is set to 1e-6. */ static const double fd_step_length; /** - * Prepare internal data - * structures and fill in values - * independent of the - * cell. Returns a pointer to an - * object of which the caller of - * this function then has to - * assume ownership (which - * includes destruction when it - * is no more needed). + * Prepare internal data structures and fill in values independent of the + * cell. Returns a pointer to an object of which the caller of this function + * then has to assume ownership (which includes destruction when it is no + * more needed). */ virtual typename Mapping::InternalDataBase * get_data (const UpdateFlags flags, @@ -2665,16 +1860,10 @@ private: const Quadrature &quadrature) const = 0; /** - * Prepare internal data - * structure for transformation - * of faces and fill in values - * independent of the - * cell. Returns a pointer to an - * object of which the caller of - * this function then has to - * assume ownership (which - * includes destruction when it - * is no more needed). + * Prepare internal data structure for transformation of faces and fill in + * values independent of the cell. Returns a pointer to an object of which + * the caller of this function then has to assume ownership (which includes + * destruction when it is no more needed). */ virtual typename Mapping::InternalDataBase * get_face_data (const UpdateFlags flags, @@ -2682,16 +1871,10 @@ private: const Quadrature &quadrature) const; /** - * Prepare internal data - * structure for transformation - * of children of faces and fill - * in values independent of the - * cell. Returns a pointer to an - * object of which the caller of - * this function then has to - * assume ownership (which - * includes destruction when it - * is no more needed). + * Prepare internal data structure for transformation of children of faces + * and fill in values independent of the cell. Returns a pointer to an + * object of which the caller of this function then has to assume ownership + * (which includes destruction when it is no more needed). */ virtual typename Mapping::InternalDataBase * get_subface_data (const UpdateFlags flags, @@ -2699,15 +1882,11 @@ private: const Quadrature &quadrature) const; /** - * Fill the fields of - * FEValues. This function - * performs all the operations - * needed to compute the data of an - * FEValues object. + * Fill the fields of FEValues. This function performs all the operations + * needed to compute the data of an FEValues object. * - * The same function in - * @p mapping must have been - * called for the same cell first! + * The same function in @p mapping must have been called for the same cell + * first! */ virtual void fill_fe_values (const Mapping &mapping, @@ -2719,15 +1898,11 @@ private: CellSimilarity::Similarity &cell_similarity) const = 0; /** - * Fill the fields of - * FEFaceValues. This function - * performs all the operations - * needed to compute the data of an - * FEFaceValues object. + * Fill the fields of FEFaceValues. This function performs all the + * operations needed to compute the data of an FEFaceValues object. * - * The same function in - * @p mapping must have been - * called for the same cell first! + * The same function in @p mapping must have been called for the same cell + * first! */ virtual void fill_fe_face_values (const Mapping &mapping, @@ -2739,15 +1914,11 @@ private: FEValuesData &data) const = 0; /** - * Fill the fields of - * FESubfaceValues. This function - * performs all the operations - * needed to compute the data of an - * FESubfaceValues object. + * Fill the fields of FESubfaceValues. This function performs all the + * operations needed to compute the data of an FESubfaceValues object. * - * The same function in - * @p mapping must have been - * called for the same cell first! + * The same function in @p mapping must have been called for the same cell + * first! */ virtual void fill_fe_subface_values (const Mapping &mapping, diff --git a/deal.II/include/deal.II/fe/fe_q_base.h b/deal.II/include/deal.II/fe/fe_q_base.h index 9ecbfd0eb9..e3fe259c17 100644 --- a/deal.II/include/deal.II/fe/fe_q_base.h +++ b/deal.II/include/deal.II/fe/fe_q_base.h @@ -14,6 +14,7 @@ #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -95,6 +96,62 @@ public: virtual bool has_support_on_face (const unsigned int shape_index, const unsigned int face_index) const; + /** + * Projection from a fine grid space onto a coarse grid space. Overrides the + * respective method in FiniteElement, implementing lazy evaluation + * (initialize when requested). + * + * If this projection operator is associated with a matrix @p P, then the + * restriction of this matrix @p P_i to a single child cell is returned + * here. + * + * The matrix @p P is the concatenation or the sum of the cell matrices @p + * P_i, depending on the #restriction_is_additive_flags. This distinguishes + * interpolation (concatenation) and projection with respect to scalar + * products (summation). + * + * Row and column indices are related to coarse grid and fine grid spaces, + * respectively, consistent with the definition of the associated operator. + * + * If projection matrices are not implemented in the derived finite element + * class, this function aborts with ExcProjectionVoid. You can check whether + * this is the case by calling the restriction_is_implemented() or the + * isotropic_restriction_is_implemented() function. + */ + virtual const FullMatrix & + get_restriction_matrix (const unsigned int child, + const RefinementCase &refinement_case=RefinementCase::isotropic_refinement) const; + + /** + * Embedding matrix between grids. Overrides the respective method in + * FiniteElement, implementing lazy evaluation (initialize when queried). + * + * The identity operator from a coarse grid space into a fine grid space is + * associated with a matrix @p P. The restriction of this matrix @p P_i to a + * single child cell is returned here. + * + * The matrix @p P is the concatenation, not the sum of the cell matrices @p + * P_i. That is, if the same non-zero entry j,k exists in in two + * different child matrices @p P_i, the value should be the same in both + * matrices and it is copied into the matrix @p P only once. + * + * Row and column indices are related to fine grid and coarse grid spaces, + * respectively, consistent with the definition of the associated operator. + * + * These matrices are used by routines assembling the prolongation matrix + * for multi-level methods. Upon assembling the transfer matrix between + * cells using this matrix array, zero elements in the prolongation matrix + * are discarded and will not fill up the transfer matrix. + * + * If projection matrices are not implemented in the derived finite element + * class, this function aborts with ExcEmbeddingVoid. You can check whether + * this is the case by calling the prolongation_is_implemented() or the + * isotropic_prolongation_is_implemented() function. + */ + virtual const FullMatrix & + get_prolongation_matrix (const unsigned int child, + const RefinementCase &refinement_case=RefinementCase::isotropic_refinement) const; + /** * @name Functions to support hp * @{ @@ -182,16 +239,6 @@ protected: */ void initialize_constraints (const std::vector > &points); - /** - * Initialize the embedding matrices. Called from initialize(). - */ - void initialize_embedding (); - - /** - * Initialize the restriction matrices. Called from initialize(). - */ - void initialize_restriction (); - /** * Initialize the @p unit_support_points field of the FiniteElement * class. Called from initialize(). @@ -222,6 +269,12 @@ protected: * Declare implementation friend. */ friend struct FE_Q_Base::Implementation; + +private: + /* + * Mutex for protecting initialization of restriction and embedding matrix. + */ + mutable Threads::Mutex mutex; }; diff --git a/deal.II/include/deal.II/fe/fe_system.h b/deal.II/include/deal.II/fe/fe_system.h index e5f408b544..650dbc5a43 100644 --- a/deal.II/include/deal.II/fe/fe_system.h +++ b/deal.II/include/deal.II/fe/fe_system.h @@ -17,6 +17,7 @@ #include +#include #include #include #include @@ -152,32 +153,23 @@ class FESystem : public FiniteElement public: /** - * Constructor. Take a finite element type - * and the number of elements you want to - * group together using this class. + * Constructor. Take a finite element type and the number of elements you + * want to group together using this class. * - * In fact, the object @p fe is not used, - * apart from getting the number of dofs - * per vertex, line, etc for that finite - * element class. The objects creates its - * own copy of the finite element object - * at construction time (but after - * the initialization of the base class - * @p FiniteElement, which is why we need - * a valid finite element object passed - * to the constructor). + * In fact, the object @p fe is not used, apart from getting the number of + * dofs per vertex, line, etc for that finite element class. The objects + * creates its own copy of the finite element object at construction time + * (but after the initialization of the base class @p FiniteElement, which + * is why we need a valid finite element object passed to the constructor). * - * Obviously, the template finite element - * class needs to be of the same dimension - * as is this object. + * Obviously, the template finite element class needs to be of the same + * dimension as is this object. */ FESystem (const FiniteElement &fe, const unsigned int n_elements); /** - * Constructor for mixed - * discretizations with two - * base elements. + * Constructor for mixed discretizations with two base elements. * * See the other constructor. */ @@ -185,9 +177,7 @@ public: const FiniteElement &fe2, const unsigned int n2); /** - * Constructor for mixed - * discretizations with three - * base elements. + * Constructor for mixed discretizations with three base elements. * * See the other constructor. */ @@ -196,9 +186,7 @@ public: const FiniteElement &fe3, const unsigned int n3); /** - * Constructor for mixed - * discretizations with four - * base elements. + * Constructor for mixed discretizations with four base elements. * * See the other constructor. */ @@ -208,9 +196,7 @@ public: const FiniteElement &fe4, const unsigned int n4); /** - * Constructor for mixed - * discretizations with five - * base elements. + * Constructor for mixed discretizations with five base elements. * * See the other constructor. */ @@ -221,15 +207,9 @@ public: const FiniteElement &fe5, const unsigned int n5); /** - * Same as above but for any - * number of base - * elements. Pointers to the base - * elements and their - * multiplicities are passed as - * vectors to this - * constructor. The length of - * these vectors is assumed to be - * equal. + * Same as above but for any number of base elements. Pointers to the base + * elements and their multiplicities are passed as vectors to this + * constructor. The length of these vectors is assumed to be equal. */ FESystem (const std::vector*> &fes, @@ -241,172 +221,95 @@ public: virtual ~FESystem (); /** - * Return a string that uniquely - * identifies a finite - * element. This element returns - * a string that is composed of - * the strings - * @p name1...@p nameN returned - * by the basis elements. From - * these, we create a sequence - * FESystem[name1^m1-name2^m2-...-nameN^mN], - * where @p mi are the - * multiplicities of the basis - * elements. If a multiplicity is - * equal to one, then the - * superscript is omitted. + * Return a string that uniquely identifies a finite element. This element + * returns a string that is composed of the strings @p name1...@p nameN + * returned by the basis elements. From these, we create a sequence + * FESystem[name1^m1-name2^m2-...-nameN^mN], where @p mi are + * the multiplicities of the basis elements. If a multiplicity is equal to + * one, then the superscript is omitted. */ virtual std::string get_name () const; /** - * Return the value of the - * @p ith shape function at the - * point @p p. @p p is a point - * on the reference element. Since - * this finite element is always - * vector-valued, we return the - * value of the only non-zero - * component of the vector value - * of this shape function. If the - * shape function has more than - * one non-zero component (which - * we refer to with the term - * non-primitive), then throw an - * exception of type - * @p ExcShapeFunctionNotPrimitive. + * Return the value of the @p ith shape function at the point @p p. @p p is + * a point on the reference element. Since this finite element is always + * vector-valued, we return the value of the only non-zero component of the + * vector value of this shape function. If the shape function has more than + * one non-zero component (which we refer to with the term non-primitive), + * then throw an exception of type @p ExcShapeFunctionNotPrimitive. * - * An - * @p ExcUnitShapeValuesDoNotExist - * is thrown if the shape values - * of the @p FiniteElement - * (corresponding to the @p ith - * shape function) depend on the - * shape of the cell in real - * space. + * An @p ExcUnitShapeValuesDoNotExist is thrown if the shape values of the + * @p FiniteElement (corresponding to the @p ith shape function) depend on + * the shape of the cell in real space. */ virtual double shape_value (const unsigned int i, const Point &p) const; /** - * Return the value of the - * @p componentth vector - * component of the @p ith shape - * function at the point - * @p p. See the - * FiniteElement base - * class for more information - * about the semantics of this - * function. + * Return the value of the @p componentth vector component of the @p ith + * shape function at the point @p p. See the FiniteElement base class for + * more information about the semantics of this function. * - * Since this element is vector - * valued in general, it relays - * the computation of these - * values to the base elements. + * Since this element is vector valued in general, it relays the computation + * of these values to the base elements. */ virtual double shape_value_component (const unsigned int i, const Point &p, const unsigned int component) const; /** - * Return the gradient of the - * @p ith shape function at the - * point @p p. @p p is a point - * on the reference element, and - * likewise the gradient is the - * gradient on the unit cell with - * respect to unit cell - * coordinates. Since - * this finite element is always - * vector-valued, we return the - * value of the only non-zero - * component of the vector value - * of this shape function. If the - * shape function has more than - * one non-zero component (which - * we refer to with the term - * non-primitive), then throw an - * exception of type - * @p ExcShapeFunctionNotPrimitive. + * Return the gradient of the @p ith shape function at the point @p p. @p p + * is a point on the reference element, and likewise the gradient is the + * gradient on the unit cell with respect to unit cell coordinates. Since + * this finite element is always vector-valued, we return the value of the + * only non-zero component of the vector value of this shape function. If + * the shape function has more than one non-zero component (which we refer + * to with the term non-primitive), then throw an exception of type @p + * ExcShapeFunctionNotPrimitive. * - * An - * @p ExcUnitShapeValuesDoNotExist - * is thrown if the shape values - * of the @p FiniteElement - * (corresponding to the @p ith - * shape function) depend on the - * shape of the cell in real - * space. + * An @p ExcUnitShapeValuesDoNotExist is thrown if the shape values of the + * @p FiniteElement (corresponding to the @p ith shape function) depend on + * the shape of the cell in real space. */ virtual Tensor<1,dim> shape_grad (const unsigned int i, const Point &p) const; /** - * Return the gradient of the - * @p componentth vector - * component of the @p ith shape - * function at the point - * @p p. See the - * FiniteElement base - * class for more information - * about the semantics of this - * function. + * Return the gradient of the @p componentth vector component of the @p ith + * shape function at the point @p p. See the FiniteElement base class for + * more information about the semantics of this function. * - * Since this element is vector - * valued in general, it relays - * the computation of these - * values to the base elements. + * Since this element is vector valued in general, it relays the computation + * of these values to the base elements. */ virtual Tensor<1,dim> shape_grad_component (const unsigned int i, const Point &p, const unsigned int component) const; /** - * Return the tensor of second - * derivatives of the @p ith - * shape function at point @p p - * on the unit cell. The - * derivatives are derivatives on - * the unit cell with respect to - * unit cell coordinates. Since - * this finite element is always - * vector-valued, we return the - * value of the only non-zero - * component of the vector value - * of this shape function. If the - * shape function has more than - * one non-zero component (which - * we refer to with the term - * non-primitive), then throw an - * exception of type - * @p ExcShapeFunctionNotPrimitive. + * Return the tensor of second derivatives of the @p ith shape function at + * point @p p on the unit cell. The derivatives are derivatives on the unit + * cell with respect to unit cell coordinates. Since this finite element is + * always vector-valued, we return the value of the only non-zero component + * of the vector value of this shape function. If the shape function has + * more than one non-zero component (which we refer to with the term + * non-primitive), then throw an exception of type @p + * ExcShapeFunctionNotPrimitive. * - * An - * @p ExcUnitShapeValuesDoNotExist - * is thrown if the shape values - * of the @p FiniteElement - * (corresponding to the @p ith - * shape function) depend on the - * shape of the cell in real - * space. + * An @p ExcUnitShapeValuesDoNotExist is thrown if the shape values of the + * @p FiniteElement (corresponding to the @p ith shape function) depend on + * the shape of the cell in real space. */ virtual Tensor<2,dim> shape_grad_grad (const unsigned int i, const Point &p) const; /** - * Return the second derivatives - * of the @p componentth vector - * component of the @p ith shape - * function at the point - * @p p. See the - * FiniteElement base - * class for more information - * about the semantics of this - * function. + * Return the second derivatives of the @p componentth vector component of + * the @p ith shape function at the point @p p. See the FiniteElement base + * class for more information about the semantics of this function. * - * Since this element is vector - * valued in general, it relays - * the computation of these - * values to the base elements. + * Since this element is vector valued in general, it relays the computation + * of these values to the base elements. */ virtual Tensor<2,dim> @@ -415,76 +318,105 @@ public: const unsigned int component) const; /** - * Return the matrix - * interpolating from the given - * finite element to the present - * one. The size of the matrix is - * then @p dofs_per_cell times + * Return the matrix interpolating from the given finite element to the + * present one. The size of the matrix is then @p dofs_per_cell times * source.dofs_per_cell. * - * These matrices are available - * if source and destination - * element are both @p FESystem - * elements, have the same number - * of base elements with same - * element multiplicity, and if - * these base elements also - * implement their - * @p get_interpolation_matrix - * functions. Otherwise, an - * exception of type - * FiniteElement::ExcInterpolationNotImplemented - * is thrown. + * These matrices are available if source and destination element are both + * @p FESystem elements, have the same number of base elements with same + * element multiplicity, and if these base elements also implement their @p + * get_interpolation_matrix functions. Otherwise, an exception of type + * FiniteElement::ExcInterpolationNotImplemented is thrown. */ virtual void get_interpolation_matrix (const FiniteElement &source, FullMatrix &matrix) const; /** - * Access to a composing - * element. The index needs to be - * smaller than the number of - * base elements. Note that the - * number of base elements may in - * turn be smaller than the - * number of components of the - * system element, if the - * multiplicities are greater - * than one. + * Access to a composing element. The index needs to be smaller than the + * number of base elements. Note that the number of base elements may in + * turn be smaller than the number of components of the system element, if + * the multiplicities are greater than one. */ virtual const FiniteElement & base_element (const unsigned int index) const; /** - * Check for non-zero values on a - * face. + * Check for non-zero values on a face. * - * This function returns - * @p true, if the shape - * function @p shape_index has - * non-zero values on the face - * @p face_index. + * This function returns @p true, if the shape function @p shape_index has + * non-zero values on the face @p face_index. * - * Implementation of the - * interface in - * FiniteElement + * Implementation of the interface in FiniteElement */ virtual bool has_support_on_face (const unsigned int shape_index, const unsigned int face_index) const; /** - * Implementation of the - * respective function in the - * base class. + * Projection from a fine grid space onto a coarse grid space. Overrides the + * respective method in FiniteElement, implementing lazy evaluation + * (initialize when requested). + * + * If this projection operator is associated with a matrix @p P, then the + * restriction of this matrix @p P_i to a single child cell is returned + * here. + * + * The matrix @p P is the concatenation or the sum of the cell matrices @p + * P_i, depending on the #restriction_is_additive_flags. This distinguishes + * interpolation (concatenation) and projection with respect to scalar + * products (summation). + * + * Row and column indices are related to coarse grid and fine grid spaces, + * respectively, consistent with the definition of the associated operator. + * + * If projection matrices are not implemented in the derived finite element + * class, this function aborts with ExcProjectionVoid. You can check whether + * this is the case by calling the restriction_is_implemented() or the + * isotropic_restriction_is_implemented() function. + */ + virtual const FullMatrix & + get_restriction_matrix (const unsigned int child, + const RefinementCase &refinement_case=RefinementCase::isotropic_refinement) const; + + /** + * Embedding matrix between grids. Overrides the respective method in + * FiniteElement, implementing lazy evaluation (initialize when queried). + * + * The identity operator from a coarse grid space into a fine grid space is + * associated with a matrix @p P. The restriction of this matrix @p P_i to a + * single child cell is returned here. + * + * The matrix @p P is the concatenation, not the sum of the cell matrices @p + * P_i. That is, if the same non-zero entry j,k exists in in two + * different child matrices @p P_i, the value should be the same in both + * matrices and it is copied into the matrix @p P only once. + * + * Row and column indices are related to fine grid and coarse grid spaces, + * respectively, consistent with the definition of the associated operator. + * + * These matrices are used by routines assembling the prolongation matrix + * for multi-level methods. Upon assembling the transfer matrix between + * cells using this matrix array, zero elements in the prolongation matrix + * are discarded and will not fill up the transfer matrix. + * + * If projection matrices are not implemented in the derived finite element + * class, this function aborts with ExcEmbeddingVoid. You can check whether + * this is the case by calling the prolongation_is_implemented() or the + * isotropic_prolongation_is_implemented() function. + */ + virtual const FullMatrix & + get_prolongation_matrix (const unsigned int child, + const RefinementCase &refinement_case=RefinementCase::isotropic_refinement) const; + + /** + * Implementation of the respective function in the base class. */ virtual Point unit_support_point (const unsigned int index) const; /** - * Implementation of the - * respective function in the - * base class. + * Implementation of the respective function in the base class. */ virtual Point @@ -496,39 +428,25 @@ public: */ /** - * Return whether this element - * implements its hanging node - * constraints in the new way, - * which has to be used to make - * elements "hp compatible". + * Return whether this element implements its hanging node constraints in + * the new way, which has to be used to make elements "hp compatible". * - * This function returns @p true iff all - * its base elements return @p true for - * this function. + * This function returns @p true iff all its base elements return @p true + * for this function. */ virtual bool hp_constraints_are_implemented () const; /** - * Return the matrix - * interpolating from a face of - * of one element to the face of - * the neighboring element. - * The size of the matrix is - * then source.dofs_per_face times - * this->dofs_per_face. + * Return the matrix interpolating from a face of of one element to the face + * of the neighboring element. The size of the matrix is then + * source.dofs_per_face times this->dofs_per_face. * - * Base elements of this element will - * have to implement this function. They - * may only provide interpolation - * matrices for certain source finite - * elements, for example those from the - * same family. If they don't implement - * interpolation from a given element, - * then they must throw an exception of - * type - * FiniteElement::ExcInterpolationNotImplemented, - * which will get propagated out from - * this element. + * Base elements of this element will have to implement this function. They + * may only provide interpolation matrices for certain source finite + * elements, for example those from the same family. If they don't implement + * interpolation from a given element, then they must throw an exception of + * type FiniteElement::ExcInterpolationNotImplemented, which + * will get propagated out from this element. */ virtual void get_face_interpolation_matrix (const FiniteElement &source, @@ -536,26 +454,16 @@ public: /** - * Return the matrix - * interpolating from a face of - * of one element to the subface of - * the neighboring element. - * The size of the matrix is - * then source.dofs_per_face times - * this->dofs_per_face. + * Return the matrix interpolating from a face of of one element to the + * subface of the neighboring element. The size of the matrix is then + * source.dofs_per_face times this->dofs_per_face. * - * Base elements of this element will - * have to implement this function. They - * may only provide interpolation - * matrices for certain source finite - * elements, for example those from the - * same family. If they don't implement - * interpolation from a given element, - * then they must throw an exception of - * type - * FiniteElement::ExcInterpolationNotImplemented, - * which will get propagated out from - * this element. + * Base elements of this element will have to implement this function. They + * may only provide interpolation matrices for certain source finite + * elements, for example those from the same family. If they don't implement + * interpolation from a given element, then they must throw an exception of + * type FiniteElement::ExcInterpolationNotImplemented, which + * will get propagated out from this element. */ virtual void get_subface_interpolation_matrix (const FiniteElement &source, @@ -600,15 +508,11 @@ public: hp_quad_dof_identities (const FiniteElement &fe_other) const; /** - * Return whether this element dominates - * the one given as argument when they - * meet at a common face, - * whether it is the other way around, - * whether neither dominates, or if - * either could dominate. + * Return whether this element dominates the one given as argument when they + * meet at a common face, whether it is the other way around, whether + * neither dominates, or if either could dominate. * - * For a definition of domination, see - * FiniteElementBase::Domination and in + * For a definition of domination, see FiniteElementBase::Domination and in * particular the @ref hp_paper "hp paper". */ virtual @@ -617,37 +521,30 @@ public: //@} /** - * Determine an estimate for the - * memory consumption (in bytes) - * of this object. + * Determine an estimate for the memory consumption (in bytes) of this + * object. * - * This function is made virtual, - * since finite element objects - * are usually accessed through - * pointers to their base class, - * rather than the class itself. + * This function is made virtual, since finite element objects are usually + * accessed through pointers to their base class, rather than the class + * itself. */ virtual std::size_t memory_consumption () const; protected: /** - * Compute flags for initial - * update only. + * Compute flags for initial update only. */ virtual UpdateFlags update_once (const UpdateFlags flags) const; /** - * Compute flags for update on - * each cell. + * Compute flags for update on each cell. */ virtual UpdateFlags update_each (const UpdateFlags flags) const; /** - * @p clone function instead of - * a copy constructor. + * @p clone function instead of a copy constructor. * - * This function is needed by the - * constructors of @p FESystem. + * This function is needed by the constructors of @p FESystem. */ virtual FiniteElement *clone() const; @@ -667,15 +564,10 @@ protected: const Quadrature &quadrature) const ; /** - * Implementation of the same - * function in - * FiniteElement. + * Implementation of the same function in FiniteElement. * - * Passes on control to - * @p compute_fill that does the - * work for all three - * fill_fe*_values - * functions. + * Passes on control to @p compute_fill that does the work for all three + * fill_fe*_values functions. */ virtual void fill_fe_values (const Mapping &mapping, @@ -687,13 +579,9 @@ protected: CellSimilarity::Similarity &cell_similarity) const; /** - * Implementation of the same - * function in - * FiniteElement. + * Implementation of the same function in FiniteElement. * - * Passes on control to - * @p compute_fill that does the - * work for all three + * Passes on control to @p compute_fill that does the work for all three * fill_fe*_values functions. */ virtual void @@ -706,13 +594,9 @@ protected: FEValuesData &data) const ; /** - * Implementation of the same - * function in - * FiniteElement. + * Implementation of the same function in FiniteElement. * - * Passes on control to - * @p compute_fill that does the - * work for all three + * Passes on control to @p compute_fill that does the work for all three * fill_fe*_values functions. */ virtual void @@ -727,26 +611,14 @@ protected: /** - * Do the work for the three - * fill_fe*_values functions. + * Do the work for the three fill_fe*_values functions. * - * Calls (among other things) - * fill_fe_([sub]face)_values - * of the base elements. Calls - * @p fill_fe_values if - * face_no==invalid_face_no - * and - * sub_no==invalid_face_no; - * calls @p fill_fe_face_values - * if - * face_no==invalid_face_no - * and - * sub_no!=invalid_face_no; - * and calls - * @p fill_fe_subface_values if - * face_no!=invalid_face_no - * and - * sub_no!=invalid_face_no. + * Calls (among other things) fill_fe_([sub]face)_values of the + * base elements. Calls @p fill_fe_values if + * face_no==invalid_face_no and sub_no==invalid_face_no; + * calls @p fill_fe_face_values if face_no==invalid_face_no and + * sub_no!=invalid_face_no; and calls @p fill_fe_subface_values if + * face_no!=invalid_face_no and sub_no!=invalid_face_no. */ template void compute_fill (const Mapping &mapping, @@ -762,9 +634,7 @@ protected: private: /** - * Value to indicate that a given - * face or subface number is - * invalid. + * Value to indicate that a given face or subface number is invalid. */ static const unsigned int invalid_face_number = numbers::invalid_unsigned_int; @@ -788,42 +658,28 @@ private: /** - * Initialize the - * @p unit_support_points field - * of the FiniteElement - * class. Called from the - * constructor. + * Initialize the @p unit_support_points field of the FiniteElement + * class. Called from the constructor. */ void initialize_unit_support_points (); /** - * Initialize the - * @p unit_face_support_points field - * of the FiniteElement - * class. Called from the - * constructor. + * Initialize the @p unit_face_support_points field of the FiniteElement + * class. Called from the constructor. */ void initialize_unit_face_support_points (); /** - * Initialize the - * @p adjust_quad_dof_index_for_face_orientation_table field - * of the FiniteElement - * class. Called from the - * constructor. + * Initialize the @p adjust_quad_dof_index_for_face_orientation_table field + * of the FiniteElement class. Called from the constructor. */ void initialize_quad_dof_index_permutation (); /** - * Helper function used in the constructor: - * take a @p FiniteElementData object - * and return an object of the same type - * with the number of degrees of - * freedom per vertex, line, etc. - * multiplied by @p n. Don't touch the - * number of functions for the - * transformation from unit to real - * cell. + * Helper function used in the constructor: take a @p FiniteElementData + * object and return an object of the same type with the number of degrees + * of freedom per vertex, line, etc. multiplied by @p n. Don't touch the + * number of functions for the transformation from unit to real cell. */ static FiniteElementData multiply_dof_numbers (const FiniteElement *fe1, @@ -838,8 +694,7 @@ private: const unsigned int N5=0); /** - * Same as above but for - * any number of sub-elements. + * Same as above but for any number of sub-elements. */ static FiniteElementData multiply_dof_numbers (const std::vector*> &fes, @@ -848,15 +703,10 @@ private: /** - * Helper function used in the - * constructor: takes a - * @p FiniteElement object and - * returns an boolean vector - * including the - * @p restriction_is_additive_flags - * of the mixed element - * consisting of @p N elements - * of the sub-element @p fe. + * Helper function used in the constructor: takes a @p FiniteElement object + * and returns an boolean vector including the @p + * restriction_is_additive_flags of the mixed element consisting of @p N + * elements of the sub-element @p fe. */ static std::vector compute_restriction_is_additive_flags ( @@ -872,11 +722,8 @@ private: const unsigned int N5=0); /** - * Compute the named flags for a - * list of finite elements with - * multiplicities given in the - * second argument. This function - * is called from all the above + * Compute the named flags for a list of finite elements with multiplicities + * given in the second argument. This function is called from all the above * functions. */ static std::vector @@ -886,9 +733,7 @@ private: /** - * Compute the non-zero vector - * components of a composed - * finite element. + * Compute the non-zero vector components of a composed finite element. */ static std::vector compute_nonzero_components (const FiniteElement *fe1, @@ -903,26 +748,18 @@ private: const unsigned int N5=0); /** - * Compute the nonzero components - * of a list of finite elements - * with multiplicities given in - * the second argument. This - * function is called from all - * the above functions. + * Compute the nonzero components of a list of finite elements with + * multiplicities given in the second argument. This function is called from + * all the above functions. */ static std::vector compute_nonzero_components (const std::vector*> &fes, const std::vector &multiplicities); /** - * This function is simply - * singled out of the - * constructors since there are - * several of them. It sets up - * the index table for the system - * as well as @p restriction and - * @p prolongation - * matrices. + * This function is simply singled out of the constructors since there are + * several of them. It sets up the index table for the system as well as @p + * restriction and @p prolongation matrices. */ void initialize (const std::vector*> &fes, const std::vector &multiplicities); @@ -943,176 +780,116 @@ private: void build_interface_constraints (); /** - * A function that computes the - * hp_vertex_dof_identities(), - * hp_line_dof_identities(), or - * hp_quad_dof_identities(), depending on - * the value of the template parameter. + * A function that computes the hp_vertex_dof_identities(), + * hp_line_dof_identities(), or hp_quad_dof_identities(), depending on the + * value of the template parameter. */ template std::vector > hp_object_dof_identities (const FiniteElement &fe_other) const; /** - * Usually: Fields of - * cell-independent data. + * Usually: Fields of cell-independent data. * - * However, here, this class does - * not itself store the data but - * only pointers to - * @p InternalData objects for - * each of the base elements. + * However, here, this class does not itself store the data but only + * pointers to @p InternalData objects for each of the base elements. */ class InternalData : public FiniteElement::InternalDataBase { public: /** - * Constructor. Is called by - * the @p get_data - * function. Sets the size of - * the @p base_fe_datas - * vector to - * @p n_base_elements. + * Constructor. Is called by the @p get_data function. Sets the size of + * the @p base_fe_datas vector to @p n_base_elements. */ InternalData (const unsigned int n_base_elements); /** - * Destructor. Deletes all - * @p InternalDatas whose - * pointers are stored by the - * @p base_fe_datas - * vector. + * Destructor. Deletes all @p InternalDatas whose pointers are stored by + * the @p base_fe_datas vector. */ ~InternalData(); /** - * Flag indicating whether - * second derivatives shall - * be computed. + * Flag indicating whether second derivatives shall be computed. */ bool compute_hessians; /** - * Gives write-access to the - * pointer to a - * @p InternalData of the - * @p base_noth base - * element. + * Gives write-access to the pointer to a @p InternalData of the @p + * base_noth base element. */ void set_fe_data(const unsigned int base_no, typename FiniteElement::InternalDataBase *); /** - * Gives read-access to the - * pointer to a - * @p InternalData of the - * @p base_noth base element. + * Gives read-access to the pointer to a @p InternalData of the @p + * base_noth base element. */ typename FiniteElement::InternalDataBase & get_fe_data (const unsigned int base_no) const; /** - * Gives write-access to the - * pointer to a - * @p FEValuesData for the - * @p base_noth base - * element. + * Gives write-access to the pointer to a @p FEValuesData for the @p + * base_noth base element. */ void set_fe_values_data (const unsigned int base_no, FEValuesData *); /** - * Gives read-access to the - * pointer to a - * @p FEValuesData for the - * @p base_noth base element. + * Gives read-access to the pointer to a @p FEValuesData for the @p + * base_noth base element. */ FEValuesData &get_fe_values_data (const unsigned int base_no) const; /** - * Deletes the - * @p FEValuesData the - * fe_datas[base_no] - * pointer is pointing - * to. Sets - * fe_datas[base_no] to - * zero. + * Deletes the @p FEValuesData the fe_datas[base_no] pointer is + * pointing to. Sets fe_datas[base_no] to zero. * - * This function is used to - * delete @p FEValuesData - * that are needed only on - * the first cell but not any - * more afterwards. This is - * the case for - * e.g. Lagrangian elements - * (see e.g. @p FE_Q - * classes). + * This function is used to delete @p FEValuesData that are needed only on + * the first cell but not any more afterwards. This is the case for + * e.g. Lagrangian elements (see e.g. @p FE_Q classes). */ void delete_fe_values_data (const unsigned int base_no); /** - * Set the @p first_cell - * flag to @p false. Used by - * the @p FEValues class to - * indicate that we have - * already done the work on - * the first cell. + * Set the @p first_cell flag to @p false. Used by the @p FEValues class + * to indicate that we have already done the work on the first cell. * - * In addition to calling the - * respective function of the - * base class, this function - * also calls the functions - * of the sub-data objects. + * In addition to calling the respective function of the base class, this + * function also calls the functions of the sub-data objects. */ virtual void clear_first_cell (); private: /** - * Pointers to - * @p InternalData objects - * for each of the base - * elements. They are - * accessed to by the - * @p set_ and - * @p get_fe_data functions. + * Pointers to @p InternalData objects for each of the base elements. They + * are accessed to by the @p set_ and @p get_fe_data functions. * - * The size of this vector is - * set to @p n_base_elements - * by the InternalData - * constructor. It is filled - * by the @p get_data - * function. Note that since - * the data for each instance - * of a base class is - * necessarily the same, we - * only need as many of these - * objects as there are base - * elements, irrespective of - * their multiplicity. + * The size of this vector is set to @p n_base_elements by the + * InternalData constructor. It is filled by the @p get_data + * function. Note that since the data for each instance of a base class is + * necessarily the same, we only need as many of these objects as there + * are base elements, irrespective of their multiplicity. */ typename std::vector::InternalDataBase *> base_fe_datas; /** - * Pointers to the - * @p FEValuesData objects - * that are given to the - * @p fill_fe_values - * function of the base - * elements. They are - * accessed to by the - * @p set_ and - * @p get_fe_values_data - * functions. + * Pointers to the @p FEValuesData objects that are given to the @p + * fill_fe_values function of the base elements. They are accessed to by + * the @p set_ and @p get_fe_values_data functions. * - * The size of this vector is - * set to @p n_base_elements - * by the InternalData - * constructor. + * The size of this vector is set to @p n_base_elements by the + * InternalData constructor. */ std::vector *> base_fe_values_datas; }; + + /* + * Mutex for protecting initialization of restriction and embedding matrix. + */ + mutable Threads::Mutex mutex; }; diff --git a/deal.II/source/fe/fe.cc b/deal.II/source/fe/fe.cc index c5504e0cb8..1a9e9f9fa4 100644 --- a/deal.II/source/fe/fe.cc +++ b/deal.II/source/fe/fe.cc @@ -180,11 +180,9 @@ FiniteElement::FiniteElement ( ExcInternalError()); }; - // initialize some tables in the - // default way, i.e. if there is - // only one (vector-)component; if - // the element is not primitive, - // leave these tables empty. + // initialize some tables in the default way, i.e. if there is only one + // (vector-)component; if the element is not primitive, leave these tables + // empty. if (this->is_primitive()) { system_to_component_table.resize(this->dofs_per_cell); @@ -200,15 +198,11 @@ FiniteElement::FiniteElement ( face_system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j); } } - // Fill with default value; may be - // changed by constructor of - // derived class. + // Fill with default value; may be changed by constructor of derived class. base_to_block_indices.reinit(1,1); - // initialize the restriction and - // prolongation matrices. the default - // contructur of FullMatrix initializes - // them with size zero + // initialize the restriction and prolongation matrices. the default + // contructur of FullMatrix initializes them with size zero prolongation.resize(RefinementCase::isotropic_refinement); restriction.resize(RefinementCase::isotropic_refinement); for (unsigned int ref=RefinementCase::cut_x; @@ -337,11 +331,9 @@ FiniteElement::get_restriction_matrix (const unsigned int child, ExcMessage("Restriction matrices are only available for refined cells!")); Assert (child::n_children(RefinementCase(refinement_case)), ExcIndexRange(child,0,GeometryInfo::n_children(RefinementCase(refinement_case)))); - // we use refinement_case-1 here. the -1 takes care - // of the origin of the vector, as for - // RefinementCase::no_refinement (=0) there is no - // data available and so the vector indices - // are shifted + // we use refinement_case-1 here. the -1 takes care of the origin of the + // vector, as for RefinementCase::no_refinement (=0) there is no data + // available and so the vector indices are shifted Assert (restriction[refinement_case-1][child].n() == this->dofs_per_cell, ExcProjectionVoid()); return restriction[refinement_case-1][child]; } @@ -614,6 +606,8 @@ FiniteElement::prolongation_is_implemented () const for (unsigned int c=0; c::n_children(RefinementCase(ref_case)); ++c) { + // make sure also the lazily initialized matrices are created + get_prolongation_matrix(c, RefinementCase(ref_case)); Assert ((prolongation[ref_case-1][c].m() == this->dofs_per_cell) || (prolongation[ref_case-1][c].m() == 0), ExcInternalError()); @@ -638,6 +632,8 @@ FiniteElement::restriction_is_implemented () const for (unsigned int c=0; c::n_children(RefinementCase(ref_case)); ++c) { + // make sure also the lazily initialized matrices are created + get_restriction_matrix(c, RefinementCase(ref_case)); Assert ((restriction[ref_case-1][c].m() == this->dofs_per_cell) || (restriction[ref_case-1][c].m() == 0), ExcInternalError()); @@ -662,6 +658,8 @@ FiniteElement::isotropic_prolongation_is_implemented () const for (unsigned int c=0; c::n_children(RefinementCase(ref_case)); ++c) { + // make sure also the lazily initialized matrices are created + get_prolongation_matrix(c, RefinementCase(ref_case)); Assert ((prolongation[ref_case-1][c].m() == this->dofs_per_cell) || (prolongation[ref_case-1][c].m() == 0), ExcInternalError()); @@ -686,6 +684,8 @@ FiniteElement::isotropic_restriction_is_implemented () const for (unsigned int c=0; c::n_children(RefinementCase(ref_case)); ++c) { + // make sure also the lazily initialized matrices are created + get_restriction_matrix(c, RefinementCase(ref_case)); Assert ((restriction[ref_case-1][c].m() == this->dofs_per_cell) || (restriction[ref_case-1][c].m() == 0), ExcInternalError()); diff --git a/deal.II/source/fe/fe_q_base.cc b/deal.II/source/fe/fe_q_base.cc index b99c8b5d8a..9705535b40 100644 --- a/deal.II/source/fe/fe_q_base.cc +++ b/deal.II/source/fe/fe_q_base.cc @@ -58,13 +58,13 @@ namespace FE_Q_Helper - // in initialize_embedding() and initialize_restriction(), want to undo + // in get_restriction_matrix() and get_prolongation_matrix(), want to undo // tensorization on inner loops for performance reasons. this clears a // dim-array template inline void - zero_indices (unsigned int indices[dim]) + zero_indices (unsigned int (&indices)[dim]) { for (unsigned int d=0; d inline void - increment_indices (unsigned int indices[dim], - const unsigned int dofs1d) + increment_indices (unsigned int (&indices)[dim], + const unsigned int dofs1d) { ++indices[0]; for (int d=0; d::initialize (const std::vector > &points) // reinit constraints initialize_constraints (points); - // Reinit the vectors of restriction and prolongation matrices to the right - // sizes and compute the matrices - this->reinit_restriction_and_prolongation_matrices(); - initialize_embedding(); - initialize_restriction(); + // do not initialize embedding and restriction here. these matrices are + // initialized on demand in get_restriction_matrix and + // get_prolongation_matrix this->initialize_quad_dof_index_permutation(); } @@ -1009,267 +1007,292 @@ FE_Q_Base::initialize_constraints (const std::vector template -void -FE_Q_Base::initialize_embedding () +const FullMatrix& +FE_Q_Base +::get_prolongation_matrix (const unsigned int child, + const RefinementCase &refinement_case) const { - // distinguish q/q_dg0 case: only treat Q dofs first - const unsigned int q_dofs_per_cell = Utilities::fixed_power(this->degree+1); + Assert (refinement_case::isotropic_refinement+1, + ExcIndexRange(refinement_case,0,RefinementCase::isotropic_refinement+1)); + Assert (refinement_case!=RefinementCase::no_refinement, + ExcMessage("Prolongation matrices are only available for refined cells!")); + Assert (child::n_children(refinement_case), + ExcIndexRange(child,0,GeometryInfo::n_children(refinement_case))); + + // needs initialization? + Threads::Mutex::ScopedLock(this->mutex); + if (this->prolongation[refinement_case-1][child].n() == 0) + { + // distinguish q/q_dg0 case: only treat Q dofs first + const unsigned int q_dofs_per_cell = Utilities::fixed_power(this->degree+1); - // compute the interpolation matrices in much the same way as we do for the - // constraints. it's actually simpler here, since we don't have this weird - // renumbering stuff going on. The trick is again that we the interpolation - // matrix is formed by a permutation of the indices of the cell matrix. The - // value eps is used a threshold to decide when certain evaluations of the - // Lagrange polynomials are zero or one. - const double eps = 1e-15*this->degree*dim; + // compute the interpolation matrices in much the same way as we do for + // the constraints. it's actually simpler here, since we don't have this + // weird renumbering stuff going on. The trick is again that we the + // interpolation matrix is formed by a permutation of the indices of the + // cell matrix. The value eps is used a threshold to decide when certain + // evaluations of the Lagrange polynomials are zero or one. + const double eps = 1e-15*this->degree*dim; #ifdef DEBUG - // in DEBUG mode, check that the evaluation of support points in the current - // numbering gives the identity operation - for (unsigned int i=0; ipoly_space.compute_value - (i, this->unit_support_points[i])) < eps, - ExcInternalError()); - for (unsigned int j=0; jpoly_space.compute_value - (i, this->unit_support_points[j])) < eps, + // in DEBUG mode, check that the evaluation of support points in the + // current numbering gives the identity operation + for (unsigned int i=0; ipoly_space.compute_value + (i, this->unit_support_points[i])) < eps, ExcInternalError()); - } + for (unsigned int j=0; jpoly_space.compute_value + (i, this->unit_support_points[j])) < eps, + ExcInternalError()); + } #endif - // to efficiently evaluate the polynomial at the subcell, make use of the - // tensor product structure of this element and only evaluate 1D information - // from the polynomial. This makes the cost of this function almost - // negligible also for high order elements - const unsigned int dofs1d = this->degree+1; - std::vector > - subcell_evaluations (dim, Table<2,double>(dofs1d, dofs1d)); - const std::vector &index_map_inverse = - this->poly_space.get_numbering_inverse(); - - // helper value: step size how to walk through diagonal and how many points - // we have left apart from the first dimension - unsigned int step_size_diag = 0; - { - unsigned int factor = 1; - for (unsigned int d=0; ddegree+1; + std::vector > + subcell_evaluations (dim, Table<2,double>(dofs1d, dofs1d)); + const std::vector &index_map_inverse = + this->poly_space.get_numbering_inverse(); - // next evaluate the functions for the different refinement cases. - for (unsigned int ref=0; ref::isotropic_refinement; ++ref) - for (unsigned int child=0; child::n_children(RefinementCase(ref+1)); ++child) + // helper value: step size how to walk through diagonal and how many + // points we have left apart from the first dimension + unsigned int step_size_diag = 0; { - // go through the points in diagonal to capture variation in all - // directions simultaneously - for (unsigned int j=0; j p_subcell = this->unit_support_points[diag_comp]; - const Point p_cell = - GeometryInfo::child_to_cell_coordinates (p_subcell, child, - RefinementCase(ref+1)); - for (unsigned int i=0; i point; - point[0] = p_cell[d]; - const double cell_value = - this->poly_space.compute_value(index_map_inverse[i], point); - - // cut off values that are too small. note that we have here - // Lagrange interpolation functions, so they should be zero - // at almost all points, and one at the others, at least on - // the subcells. so set them to their exact values - // - // the actual cut-off value is somewhat fuzzy, but it works - // for 2e-13*degree*dim (see above), which is kind of - // reasonable given that we compute the values of the - // polynomials via an degree-step recursion and then - // multiply the 1d-values. this gives us a linear growth in - // degree*dim, times a small constant. - // - // the embedding matrix is given by applying the inverse of - // the subcell matrix on the cell_interpolation - // matrix. since the subcell matrix is actually only a - // permutation vector, all we need to do is to switch the - // rows we write the data into. moreover, cut off very small - // values here - if (std::fabs(cell_value) < eps) - subcell_evaluations[d](j,i) = 0; - else - subcell_evaluations[d](j,i) = cell_value; - } + step_size_diag += factor; + factor *= dofs1d; } + } - // now expand from 1D info. block innermost dimension (x_0) in order - // to avoid difficult checks at innermost loop - unsigned int j_indices[dim]; - FE_Q_Helper::zero_indices (j_indices); - for (unsigned int j=0; j (i_indices); - for (unsigned int i=0; iprolongation[ref][child](j,i) = - // this->poly_space.compute_value (i, p_cell); - for (unsigned int jj=0; jjprolongation[ref][child](j_ind,index_map_inverse[i+ii]) - = val_extra_dim * subcell_evaluations[0](jj,ii); - } + FullMatrix& prolongate = const_cast&> + (this->prolongation[refinement_case-1][child]); + prolongate.reinit(this->dofs_per_cell, this->dofs_per_cell); - // update indices that denote the tensor product position. a - // bit fuzzy and therefore not done for innermost x_0 - // direction - FE_Q_Helper::increment_indices (i_indices, dofs1d); + // go through the points in diagonal to capture variation in all + // directions simultaneously + for (unsigned int j=0; j p_subcell = this->unit_support_points[diag_comp]; + const Point p_cell = + GeometryInfo::child_to_cell_coordinates (p_subcell, child, + refinement_case); + for (unsigned int i=0; i point; + point[0] = p_cell[d]; + const double cell_value = + this->poly_space.compute_value(index_map_inverse[i], point); + + // cut off values that are too small. note that we have here + // Lagrange interpolation functions, so they should be zero at + // almost all points, and one at the others, at least on the + // subcells. so set them to their exact values + // + // the actual cut-off value is somewhat fuzzy, but it works + // for 2e-13*degree*dim (see above), which is kind of + // reasonable given that we compute the values of the + // polynomials via an degree-step recursion and then multiply + // the 1d-values. this gives us a linear growth in degree*dim, + // times a small constant. + // + // the embedding matrix is given by applying the inverse of + // the subcell matrix on the cell_interpolation matrix. since + // the subcell matrix is actually only a permutation vector, + // all we need to do is to switch the rows we write the data + // into. moreover, cut off very small values here + if (std::fabs(cell_value) < eps) + subcell_evaluations[d](j,i) = 0; + else + subcell_evaluations[d](j,i) = cell_value; } - Assert (i_indices[dim-1] == 1, ExcInternalError()); - FE_Q_Helper::increment_indices (j_indices, dofs1d); - } + } + + // now expand from 1D info. block innermost dimension (x_0) in order to + // avoid difficult checks at innermost loop + unsigned int j_indices[dim]; + FE_Q_Helper::zero_indices (j_indices); + for (unsigned int j=0; j (i_indices); + for (unsigned int i=0; ipoly_space.compute_value (i, p_cell) + for (unsigned int jj=0; jj (i_indices, dofs1d); + } + Assert (i_indices[dim-1] == 1, ExcInternalError()); + FE_Q_Helper::increment_indices (j_indices, dofs1d); + } - // the discontinuous node is simply mapped on the discontinuous node - // on the child element - if (q_dofs_per_cell < this->dofs_per_cell) - this->prolongation[ref][child](q_dofs_per_cell,q_dofs_per_cell) = 1.; + // the discontinuous node is simply mapped on the discontinuous node on + // the child element + if (q_dofs_per_cell < this->dofs_per_cell) + prolongate(q_dofs_per_cell,q_dofs_per_cell) = 1.; - // and make sure that the row sum is 1. this must be so since for this - // element, the shape functions add up to one + // and make sure that the row sum is 1. this must be so since for this + // element, the shape functions add up to one #ifdef DEBUG - for (unsigned int row=0; rowdofs_per_cell; ++row) - { - double sum = 0; - for (unsigned int col=0; coldofs_per_cell; ++col) - sum += this->prolongation[ref][child](row,col); - Assert (std::fabs(sum-1.) < eps, ExcInternalError()); - } + for (unsigned int row=0; rowdofs_per_cell; ++row) + { + double sum = 0; + for (unsigned int col=0; coldofs_per_cell; ++col) + sum += prolongate(row,col); + Assert (std::fabs(sum-1.) < eps, ExcInternalError()); + } #endif - } + } + + // finally return the matrix + return this->prolongation[refinement_case-1][child]; } template -void -FE_Q_Base::initialize_restriction () +const FullMatrix& +FE_Q_Base +::get_restriction_matrix (const unsigned int child, + const RefinementCase &refinement_case) const { - // distinguish q/q_dg0 case - const unsigned int q_dofs_per_cell = Utilities::fixed_power(this->degree+1); + Assert (refinement_case::isotropic_refinement+1, + ExcIndexRange(refinement_case,0,RefinementCase::isotropic_refinement+1)); + Assert (refinement_case!=RefinementCase::no_refinement, + ExcMessage("Restriction matrices are only available for refined cells!")); + Assert (child::n_children(refinement_case), + ExcIndexRange(child,0,GeometryInfo::n_children(refinement_case))); + + // needs initialization? + Threads::Mutex::ScopedLock(this->mutex); + if (this->restriction[refinement_case-1][child].n() == 0) + { + FullMatrix &restriction = + const_cast &>(this->restriction[refinement_case-1][child]); + + // distinguish q/q_dg0 case + const unsigned int q_dofs_per_cell = Utilities::fixed_power(this->degree+1); + + // for these Lagrange interpolation polynomials, construction of the + // restriction matrices is relatively simple. the reason is that the + // interpolation points on the mother cell are (except for the case with + // arbitrary nonequidistant nodes) always also interpolation points for + // some shape function on one or the other child, because we have chosen + // equidistant Lagrange interpolation points for the polynomials + // + // so the only thing we have to find out is: for each shape function on + // the mother cell, which is the child cell (possibly more than one) on + // which it is located, and which is the corresponding shape function + // there. rather than doing it for the shape functions on the mother + // cell, we take the interpolation points there + // + // note that the interpolation point of a shape function can be on the + // boundary between subcells. in that case, restriction from children to + // mother may produce two or more entries for a dof on the mother + // cell. however, this doesn't hurt: since the element is continuous, + // the contribution from each child should yield the same result, and + // since the element is non-additive we just overwrite one value + // (compute on one child) by the same value (compute on a later child), + // so we don't have to care about this + + const double eps = 1e-15*this->degree*dim; + const std::vector &index_map_inverse = + this->poly_space.get_numbering_inverse(); - // for these Lagrange interpolation polynomials, construction of the - // restriction matrices is relatively simple. the reason is that the - // interpolation points on the mother cell are (except for the case with - // arbitrary nonequidistant nodes) always also interpolation points for some - // shape function on one or the other child, because we have chosen - // equidistant Lagrange interpolation points for the polynomials - // - // so the only thing we have to find out is: for each shape function on the - // mother cell, which is the child cell (possibly more than one) on which it - // is located, and which is the corresponding shape function there. rather - // than doing it for the shape functions on the mother cell, we take the - // interpolation points there - // - // note that the interpolation point of a shape function can be on the - // boundary between subcells. in that case, restriction from children to - // mother may produce two or more entries for a dof on the mother - // cell. however, this doesn't hurt: since the element is continuous, the - // contribution from each child should yield the same result, and since the - // element is non-additive we just overwrite one value (compute on one - // child) by the same value (compute on a later child), so we don't have to - // care about this - - const double eps = 1e-15*this->degree*dim; - const std::vector &index_map_inverse = - this->poly_space.get_numbering_inverse(); + const unsigned int dofs1d = this->degree+1; + std::vector > evaluations1d (dofs1d); - // recreate 1D polynomials for faster evaluation of polynomial - const unsigned int dofs1d = this->degree+1; - std::vector > evaluations1d (dofs1d); + restriction.reinit(this->dofs_per_cell, this->dofs_per_cell); - for (unsigned int i=0; i p_cell = this->unit_support_points[mother_dof]; - - // then find the children on which the interpolation point is located - for (unsigned int ref=RefinementCase::cut_x; ref<=RefinementCase::isotropic_refinement; ++ref) - for (unsigned int child=0; child::n_children(RefinementCase(ref)); ++child) - { - // check whether this interpolation point is inside this child - // cell - const Point p_subcell - = GeometryInfo::cell_to_child_coordinates (p_cell, child, RefinementCase(ref)); - if (GeometryInfo::is_inside_unit_cell (p_subcell)) - { - // same logic as in initialize_embedding to evaluate the - // polynomial faster than from the tensor product: since we - // evaluate all polynomials, it is much faster to just compute - // the 1D values for all polynomials before and then get the - // dim-data. - for (unsigned int j=0; j point; - point[0] = p_subcell[d]; - evaluations1d[j][d] = - this->poly_space.compute_value(index_map_inverse[j], point); - } - unsigned int j_indices[dim]; - FE_Q_Helper::zero_indices (j_indices); - double sum_check = 0; - for (unsigned int j = 0; j p_cell = this->unit_support_points[mother_dof]; + + // check whether this interpolation point is inside this child cell + const Point p_subcell + = GeometryInfo::cell_to_child_coordinates (p_cell, child, + refinement_case); + if (GeometryInfo::is_inside_unit_cell (p_subcell)) + { + // same logic as in initialize_embedding to evaluate the + // polynomial faster than from the tensor product: since we + // evaluate all polynomials, it is much faster to just compute + // the 1D values for all polynomials before and then get the + // dim-data. + for (unsigned int j=0; jrestriction[ref-1][child](mother_dof,child_dof)=1.; - else if (std::fabs(val) > eps) - this->restriction[ref-1][child](mother_dof,child_dof)=val; - sum_check += val; - } - FE_Q_Helper::increment_indices (j_indices, dofs1d); + Point point; + point[0] = p_subcell[d]; + evaluations1d[j][d] = + this->poly_space.compute_value(index_map_inverse[j], point); } - Assert (std::fabs(sum_check-1) < eps, - ExcInternalError()); - } + unsigned int j_indices[dim]; + FE_Q_Helper::zero_indices (j_indices); + double sum_check = 0; + for (unsigned int j = 0; jdofs_per_cell) - this->restriction[ref-1][child] - (this->dofs_per_cell-1,this->dofs_per_cell-1) = - 1./(double)GeometryInfo::n_children - (RefinementCase(ref)); - } + // find the child shape function(s) corresponding to + // this point. Usually this is just one function; + // however, when we use FE_Q on arbitrary nodes a parent + // support point might not be a child support point, and + // then we will get more than one nonzero value per + // row. Still, the values should sum up to 1 + const double val + = val_extra_dim * evaluations1d[jj][0]; + const unsigned int child_dof = + index_map_inverse[j+jj]; + if (std::fabs (val-1.) < eps) + restriction(mother_dof,child_dof)=1.; + else if (std::fabs(val) > eps) + restriction(mother_dof,child_dof)=val; + sum_check += val; + } + FE_Q_Helper::increment_indices (j_indices, dofs1d); + } + Assert (std::fabs(sum_check-1) < eps, + ExcInternalError()); + } + + // part for FE_Q_DG0 + if (q_dofs_per_cell < this->dofs_per_cell) + restriction(this->dofs_per_cell-1,this->dofs_per_cell-1) = + 1./GeometryInfo::n_children(RefinementCase(refinement_case)); + } } + return this->restriction[refinement_case-1][child]; } diff --git a/deal.II/source/fe/fe_system.cc b/deal.II/source/fe/fe_system.cc index 6bb7de8489..c3ecf56590 100644 --- a/deal.II/source/fe/fe_system.cc +++ b/deal.II/source/fe/fe_system.cc @@ -54,8 +54,7 @@ FESystem::InternalData::InternalData(const unsigned int n_base_ele template FESystem::InternalData::~InternalData() { - // delete pointers and set them to - // zero to avoid inadvertent use + // delete pointers and set them to zero to avoid inadvertent use for (unsigned int i=0; i::shape_grad_component (const unsigned int i, Assert (component < this->n_components(), ExcIndexRange (component, 0, this->n_components())); - // if this value is supposed to be - // zero, then return right away... + // if this value is supposed to be zero, then return right away... if (this->nonzero_components[i][component] == false) return Tensor<1,dim>(); - // ...otherwise: first find out to - // which of the base elements this - // desired component belongs, and - // which component within this base - // element it is + // ...otherwise: first find out to which of the base elements this desired + // component belongs, and which component within this base element it is const unsigned int base = this->component_to_base_index(component).first; const unsigned int component_in_base = this->component_to_base_index(component).second; - // then get value from base - // element. note that that will - // throw an error should the - // respective shape function not be - // primitive; thus, there is no + // then get value from base element. note that that will throw an error + // should the respective shape function not be primitive; thus, there is no // need to check this here return (base_element(base). shape_grad_component(this->system_to_base_table[i].second, @@ -505,24 +497,17 @@ FESystem::shape_grad_grad_component (const unsigned int i, Assert (component < this->n_components(), ExcIndexRange (component, 0, this->n_components())); - // if this value is supposed to be - // zero, then return right away... + // if this value is supposed to be zero, then return right away... if (this->nonzero_components[i][component] == false) return Tensor<2,dim>(); - // ...otherwise: first find out to - // which of the base elements this - // desired component belongs, and - // which component within this base - // element it is + // ...otherwise: first find out to which of the base elements this desired + // component belongs, and which component within this base element it is const unsigned int base = this->component_to_base_index(component).first; const unsigned int component_in_base = this->component_to_base_index(component).second; - // then get value from base - // element. note that that will - // throw an error should the - // respective shape function not be - // primitive; thus, there is no + // then get value from base element. note that that will throw an error + // should the respective shape function not be primitive; thus, there is no // need to check this here return (base_element(base). shape_grad_grad_component(this->system_to_base_table[i].second, @@ -545,12 +530,10 @@ FESystem::get_interpolation_matrix ( ExcDimensionMismatch (interpolation_matrix.m(), x_source_fe.dofs_per_cell)); - // there are certain conditions - // that the two elements have to - // satisfy so that this can work. + // there are certain conditions that the two elements have to satisfy so + // that this can work. // - // condition 1: the other element - // must also be a system element + // condition 1: the other element must also be a system element typedef FiniteElement FEL; AssertThrow ((x_source_fe.get_name().find ("FESystem<") == 0) @@ -559,38 +542,29 @@ FESystem::get_interpolation_matrix ( typename FEL:: ExcInterpolationNotImplemented()); - // ok, source is a system element, - // so we may be able to do the work + // ok, source is a system element, so we may be able to do the work const FESystem &source_fe = dynamic_cast&>(x_source_fe); - // condition 2: same number of - // basis elements + // condition 2: same number of basis elements AssertThrow (this->n_base_elements() == source_fe.n_base_elements(), typename FEL:: ExcInterpolationNotImplemented()); - // condition 3: same number of - // basis elements + // condition 3: same number of basis elements for (unsigned int i=0; i< this->n_base_elements(); ++i) AssertThrow (this->element_multiplicity(i) == source_fe.element_multiplicity(i), typename FEL:: ExcInterpolationNotImplemented()); - // ok, so let's try whether it - // works: - - // first let's see whether all the - // basis elements actually generate - // their interpolation matrices. if - // we get past the following loop, - // then apparently none of the - // called base elements threw an - // exception, so we're fine - // continuing and assembling the - // one big matrix from the small - // ones of the base elements + // ok, so let's try whether it works: + + // first let's see whether all the basis elements actually generate their + // interpolation matrices. if we get past the following loop, then + // apparently none of the called base elements threw an exception, so we're + // fine continuing and assembling the one big matrix from the small ones of + // the base elements std::vector > base_matrices (this->n_base_elements()); for (unsigned int i=0; in_base_elements(); ++i) { @@ -600,11 +574,9 @@ FESystem::get_interpolation_matrix ( base_matrices[i]); } - // first clear big matrix, to make - // sure that entries that would - // couple different bases (or - // multiplicity indices) are really - // zero. then assign entries + // first clear big matrix, to make sure that entries that would couple + // different bases (or multiplicity indices) are really zero. then assign + // entries interpolation_matrix = 0; for (unsigned int i=0; idofs_per_cell; ++i) for (unsigned int j=0; j::get_interpolation_matrix ( +template +const FullMatrix& +FESystem +::get_restriction_matrix (const unsigned int child, + const RefinementCase &refinement_case) const +{ + Assert (refinement_case::isotropic_refinement+1, + ExcIndexRange(refinement_case,0,RefinementCase::isotropic_refinement+1)); + Assert (refinement_case!=RefinementCase::no_refinement, + ExcMessage("Restriction matrices are only available for refined cells!")); + Assert (child::n_children(refinement_case), + ExcIndexRange(child,0,GeometryInfo::n_children(refinement_case))); + + // needs initialization? + Threads::Mutex::ScopedLock(this->mutex); + if (this->restriction[refinement_case-1][child].n() == 0) + { + // Check if some of the matrices of the base elements are void. + bool do_restriction = true; + + // shortcut for accessing local restrictions further down + std::vector *> + base_matrices(this->n_base_elements()); + + for (unsigned int i=0; in_base_elements(); ++i) + { + base_matrices[i] = + &base_element(i).get_restriction_matrix(child, refinement_case); + if (base_matrices[i]->n() != base_element(i).dofs_per_cell) + do_restriction = false; + } + Assert(do_restriction, + (typename FiniteElement::ExcProjectionVoid())); + + // if we did not encounter void matrices, initialize the matrix sizes + if (do_restriction) + { + FullMatrix &restriction = + const_cast&>(this->restriction[refinement_case-1] + [child]); + restriction.reinit(this->dofs_per_cell, this->dofs_per_cell); + + // distribute the matrices of the base finite elements to the + // matrices of this object. for this, loop over all degrees of + // freedom and take the respective entry of the underlying base + // element. + // + // note that we by definition of a base element, they are + // independent, i.e. do not couple. only DoFs that belong to the + // same instance of a base element may couple + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int j=0; jdofs_per_cell; ++j) + { + // first find out to which base element indices i and j + // belong, and which instance thereof in case the base element + // has a multiplicity greater than one. if they should not + // happen to belong to the same instance of a base element, + // then they cannot couple, so go on with the next index + if (this->system_to_base_table[i].first != + this->system_to_base_table[j].first) + continue; + + // so get the common base element and the indices therein: + const unsigned int + base = this->system_to_base_table[i].first.first; + + const unsigned int + base_index_i = this->system_to_base_table[i].second, + base_index_j = this->system_to_base_table[j].second; + + // if we are sure that DoFs i and j may couple, then copy + // entries of the matrices: + restriction(i,j) = (*base_matrices[base])(base_index_i,base_index_j); + } + } + } + + return this->restriction[refinement_case-1][child]; +} + + + +template +const FullMatrix& +FESystem +::get_prolongation_matrix (const unsigned int child, + const RefinementCase &refinement_case) const +{ + Assert (refinement_case::isotropic_refinement+1, + ExcIndexRange(refinement_case,0,RefinementCase::isotropic_refinement+1)); + Assert (refinement_case!=RefinementCase::no_refinement, + ExcMessage("Restriction matrices are only available for refined cells!")); + Assert (child::n_children(refinement_case), + ExcIndexRange(child,0,GeometryInfo::n_children(refinement_case))); + + // needs initialization? + Threads::Mutex::ScopedLock(this->mutex); + if (this->prolongation[refinement_case-1][child].n() == 0) + { + // Check if some of the matrices of the base elements are void. + bool do_prolongation = true; + + // shortcut for accessing local prolongations further down + std::vector *> + base_matrices(this->n_base_elements()); + + for (unsigned int i=0; in_base_elements(); ++i) + { + base_matrices[i] = + &base_element(i).get_prolongation_matrix(child, refinement_case); + if (base_matrices[i]->n() != base_element(i).dofs_per_cell) + do_prolongation = false; + } + Assert(do_prolongation, + (typename FiniteElement::ExcEmbeddingVoid())); + + + if (do_prolongation) + { + FullMatrix &prolongate = + const_cast &>(this->prolongation[refinement_case-1][child]); + prolongate.reinit(this->dofs_per_cell, this->dofs_per_cell); + + // distribute the matrices of the base finite elements to the + // matrices of this object. for this, loop over all degrees of + // freedom and take the respective entry of the underlying base + // element. + // + // note that we by definition of a base element, they are + // independent, i.e. do not couple. only DoFs that belong to the + // same instance of a base element may couple + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int j=0; jdofs_per_cell; ++j) + { + // first find out to which base element indices i and j + // belong, and which instance thereof in case the base element + // has a multiplicity greater than one. if they should not + // happen to belong to the same instance of a base element, + // then they cannot couple, so go on with the next index + if (this->system_to_base_table[i].first != + this->system_to_base_table[j].first) + continue; + + // so get the common base element and the indices therein: + const unsigned int + base = this->system_to_base_table[i].first.first; + + const unsigned int + base_index_i = this->system_to_base_table[i].second, + base_index_j = this->system_to_base_table[j].second; + + // if we are sure that DoFs i and j may couple, then copy + // entries of the matrices: + prolongate(i,j) = (*base_matrices[base])(base_index_i,base_index_j); + } + } + } + + return this->prolongation[refinement_case-1][child]; +} + + + //--------------------------------------------------------------------------- // Data field initialization //--------------------------------------------------------------------------- @@ -941,25 +1076,14 @@ FESystem::compute_fill ( { if (fe_data.is_first_cell()) { - // Initialize the - // FEValuesDatas for the - // base elements. - // Originally this was the - // task of - // FEValues::FEValues() but - // the latter initializes - // the FEValuesDatas only - // of the FESystem, not of - // the FEValuesDatas needed - // by the base elements - // (and: how should it know - // of their existence, - // after all). + // Initialize the FEValuesDatas for the base elements. Originally + // this was the task of FEValues::FEValues() but the latter + // initializes the FEValuesDatas only of the FESystem, not of the + // FEValuesDatas needed by the base elements (and: how should it + // know of their existence, after all). for (unsigned int base_no=0; base_non_base_elements(); ++base_no) { - // Pointer needed to get - // the update flags of the - // base element + // Pointer needed to get the update flags of the base element typename Mapping::InternalDataBase & base_fe_data = fe_data.get_fe_data(base_no); @@ -967,28 +1091,20 @@ FESystem::compute_fill ( const UpdateFlags base_update_flags = mapping_data.update_flags | base_fe_data.update_flags; - // Initialize the FEValuesDatas - // for the base elements. + // Initialize the FEValuesDatas for the base elements. FEValuesData &base_data=fe_data.get_fe_values_data(base_no); const FiniteElement &base_fe=base_element(base_no); base_data.initialize (n_q_points, base_fe, base_update_flags); } } - // fill_fe_face_values needs - // argument Quadrature - // for both cases - // dim_1==dim-1 and - // dim_1=dim. Hence the - // following workaround + // fill_fe_face_values needs argument Quadrature for both cases + // dim_1==dim-1 and dim_1=dim. Hence the following workaround const Quadrature *cell_quadrature = 0; const Quadrature *face_quadrature = 0; - // static cast to the - // common base class of - // quadrature being either - // Quadrature or - // Quadrature: + // static cast to the common base class of quadrature being either + // Quadrature or Quadrature: const Subscriptor *quadrature_base_pointer = &quadrature; if (face_no==invalid_face_number) @@ -1010,8 +1126,7 @@ FESystem::compute_fill ( = static_cast *>(quadrature_base_pointer); } - // let base elements update the - // necessary data + // let base elements update the necessary data for (unsigned int base_no=0; base_non_base_elements(); ++base_no) { const FiniteElement & @@ -1021,34 +1136,21 @@ FESystem::compute_fill ( FEValuesData & base_data = fe_data.get_fe_values_data(base_no); - //TODO: Think about a smarter alternative - // Copy quadrature points. These - // are required for computing the - // determinant in the FEPolyTensor - // class. The determinant is - // one ingredient of the Piola - // transformation, which is applied - // to correctly map the RT space from - // the reference element to the global - // coordinate system. + //TODO: Think about a smarter alternative Copy quadrature + // points. These are required for computing the determinant in the + // FEPolyTensor class. The determinant is one ingredient of the + // Piola transformation, which is applied to correctly map the RT + // space from the reference element to the global coordinate system. if (cell_similarity != CellSimilarity::translation) base_data.JxW_values = data.JxW_values; - // Make sure that in the - // case of fill_fe_values - // the data is only copied - // from base_data to data - // if base_data is - // changed. therefore use - // fe_fe_data.current_update_flags() + // Make sure that in the case of fill_fe_values the data is only + // copied from base_data to data if base_data is changed. therefore + // use fe_fe_data.current_update_flags() // - // for the case of - // fill_fe_(sub)face_values - // the data needs to be - // copied from base_data to - // data on each face, - // therefore use + // for the case of fill_fe_(sub)face_values the data needs to be + // copied from base_data to data on each face, therefore use // base_fe_data.update_flags. if (face_no==invalid_face_number) base_fe.fill_fe_values(mapping, cell, *cell_quadrature, mapping_data, @@ -1060,45 +1162,25 @@ FESystem::compute_fill ( base_fe.fill_fe_subface_values(mapping, cell, face_no, sub_no, *face_quadrature, mapping_data, base_fe_data, base_data); - // now data has been - // generated, so copy - // it. we used to work by - // looping over all base - // elements (i.e. this - // outer loop), then over - // multiplicity, then over - // the shape functions from - // that base element, but - // that requires that we - // can infer the global - // number of a shape - // function from its number - // in the base element. for - // that we had the - // component_to_system_table. + // now data has been generated, so copy it. we used to work by + // looping over all base elements (i.e. this outer loop), then over + // multiplicity, then over the shape functions from that base + // element, but that requires that we can infer the global number of + // a shape function from its number in the base element. for that we + // had the component_to_system_table. // - // however, this does of - // course no longer work - // since we have - // non-primitive - // elements. so we go the - // other way round: loop - // over all shape functions - // of the composed element, - // and here only treat - // those shape functions - // that belong to a given - // base element + // however, this does of course no longer work since we have + // non-primitive elements. so we go the other way round: loop over + // all shape functions of the composed element, and here only treat + // those shape functions that belong to a given base element //TODO: Introduce the needed table and loop only over base element shape functions. This here is not efficient at all AND very bad style const UpdateFlags base_flags(dim_1==dim ? base_fe_data.current_update_flags() : base_fe_data.update_flags); - // if the current cell is just a - // translation of the previous - // one, the underlying data has - // not changed, and we don't even - // need to enter this section + // if the current cell is just a translation of the previous one, + // the underlying data has not changed, and we don't even need to + // enter this section if (cell_similarity != CellSimilarity::translation) for (unsigned int system_index=0; system_indexdofs_per_cell; ++system_index) @@ -1108,22 +1190,12 @@ FESystem::compute_fill ( base_index = this->system_to_base_table[system_index].second; Assert (base_indexn_nonzero_components(i); @@ -1131,9 +1203,7 @@ FESystem::compute_fill ( for (unsigned int i=0; in_nonzero_components(system_index) == base_fe.n_nonzero_components(base_index), ExcInternalError()); @@ -1149,15 +1219,8 @@ FESystem::compute_fill ( data.shape_gradients[out_index+s][q]= base_data.shape_gradients[in_index+s][q]; - // _we_ handle - // computation of - // second - // derivatives, - // so the base - // elements - // should not - // have computed - // them! + // _we_ handle computation of second derivatives, so the + // base elements should not have computed them! Assert (!(base_flags & update_hessians), ExcInternalError()); }; @@ -1197,9 +1260,8 @@ template void FESystem::build_cell_tables() { - // If the system is not primitive, - // these have not been initialized - // by FiniteElement + // If the system is not primitive, these have not been initialized by + // FiniteElement this->system_to_component_table.resize(this->dofs_per_cell); this->face_system_to_component_table.resize(this->dofs_per_face); @@ -1217,20 +1279,16 @@ FESystem::build_cell_tables() Assert (total_index == this->component_to_base_table.size(), ExcInternalError()); - // Initialize index tables. - // Multi-component base elements - // have to be thought of. For - // non-primitive shape functions, - // have a special invalid index. + // Initialize index tables. Multi-component base elements have to be + // thought of. For non-primitive shape functions, have a special invalid + // index. const std::pair non_primitive_index (numbers::invalid_unsigned_int, numbers::invalid_unsigned_int); - // First enumerate vertex indices, - // where we first enumerate all - // indices on the first vertex in - // the order of the base elements, - // then of the second vertex, etc + // First enumerate vertex indices, where we first enumerate all indices on + // the first vertex in the order of the base elements, then of the second + // vertex, etc total_index = 0; for (unsigned int vertex_number=0; vertex_number::vertices_per_cell; @@ -1388,14 +1446,10 @@ template void FESystem::build_face_tables() { - // Initialize index tables. do this - // in the same way as done for the - // cell tables, except that we now - // loop over the objects of faces - - // For non-primitive shape - // functions, have a special - // invalid index + // Initialize index tables. do this in the same way as done for the cell + // tables, except that we now loop over the objects of faces + + // For non-primitive shape functions, have a special invalid index const std::pair non_primitive_index (numbers::invalid_unsigned_int, numbers::invalid_unsigned_int); @@ -1414,30 +1468,13 @@ FESystem::build_face_tables() local_index < base_element(base).dofs_per_vertex; ++local_index, ++total_index) { - // get (cell) index of - // this shape function - // inside the base - // element to see - // whether the shape - // function is - // primitive (assume - // that all shape - // functions on - // vertices share the - // same primitivity - // property; assume - // likewise for all - // shape functions - // located on lines, - // quads, etc. this - // way, we can ask for - // primitivity of only - // _one_ shape - // function, which is - // taken as - // representative for - // all others located - // on the same type of + // get (cell) index of this shape function inside the base + // element to see whether the shape function is primitive + // (assume that all shape functions on vertices share the same + // primitivity property; assume likewise for all shape functions + // located on lines, quads, etc. this way, we can ask for + // primitivity of only _one_ shape function, which is taken as + // representative for all others located on the same type of // object): const unsigned int index_in_base = (base_element(base).dofs_per_vertex*vertex_number + @@ -1480,9 +1517,7 @@ FESystem::build_face_tables() local_index < base_element(base).dofs_per_line; ++local_index, ++total_index) { - // do everything - // alike for this - // type of object + // do everything alike for this type of object const unsigned int index_in_base = (base_element(base).dofs_per_line*line_number + local_index + @@ -1526,9 +1561,7 @@ FESystem::build_face_tables() local_index < base_element(base).dofs_per_quad; ++local_index, ++total_index) { - // do everything - // alike for this - // type of object + // do everything alike for this type of object const unsigned int index_in_base = (base_element(base).dofs_per_quad*quad_number + local_index + @@ -1569,14 +1602,10 @@ FESystem::build_face_tables() template void FESystem::build_interface_constraints () { - // check whether all base elements - // implement their interface - // constraint matrices. if this is - // not the case, then leave the - // interface costraints of this - // composed element empty as well; - // however, the rest of the element - // is usable + // check whether all base elements implement their interface constraint + // matrices. if this is not the case, then leave the interface costraints of + // this composed element empty as well; however, the rest of the element is + // usable for (unsigned int base=0; basen_base_elements(); ++base) if (base_element(base).constraints_are_implemented() == false) return; @@ -1584,58 +1613,35 @@ void FESystem::build_interface_constraints () this->interface_constraints. TableBase<2,double>::reinit (this->interface_constraints_size()); - // the layout of the constraints - // matrix is described in the - // FiniteElement class. you may - // want to look there first before - // trying to understand the - // following, especially the - // mapping of the @p{m} index. + // the layout of the constraints matrix is described in the FiniteElement + // class. you may want to look there first before trying to understand the + // following, especially the mapping of the @p{m} index. // - // in order to map it to the - // fe-system class, we have to know - // which base element a degree of - // freedom within a vertex, line, - // etc belongs to. this can be - // accomplished by the - // system_to_component_index - // function in conjunction with the - // numbers - // first_{line,quad,...}_index + // in order to map it to the fe-system class, we have to know which base + // element a degree of freedom within a vertex, line, etc belongs to. this + // can be accomplished by the system_to_component_index function in + // conjunction with the numbers first_{line,quad,...}_index for (unsigned int n=0; ninterface_constraints.n(); ++n) for (unsigned int m=0; minterface_constraints.m(); ++m) { - // for the pair (n,m) find - // out which base element - // they belong to and the - // number therein + // for the pair (n,m) find out which base element they belong to and + // the number therein // - // first for the n - // index. this is simple - // since the n indices are in - // the same order as they are - // usually on a face. note - // that for the data type, - // first value in pair is - // (base element,instance of - // base element), second is - // index within this instance + // first for the n index. this is simple since the n indices are in + // the same order as they are usually on a face. note that for the + // data type, first value in pair is (base element,instance of base + // element), second is index within this instance const std::pair, unsigned int> n_index = this->face_system_to_base_table[n]; - // likewise for the m - // index. this is more - // complicated due to the - // strange ordering we have - // for the dofs on the - // refined faces. + // likewise for the m index. this is more complicated due to the + // strange ordering we have for the dofs on the refined faces. std::pair, unsigned int> m_index; switch (dim) { case 1: { - // we should never get here! - // (in 1d, the constraints matrix + // we should never get here! (in 1d, the constraints matrix // should be of size zero) Assert (false, ExcInternalError()); break; @@ -1643,21 +1649,14 @@ void FESystem::build_interface_constraints () case 2: { - // the indices m=0..d_v-1 are - // from the center vertex. - // their order is the same - // as for the first vertex - // of the whole cell, so we - // can use the - // system_to_base_table - // variable (using the - // face_s_t_base_t function would - // yield the same) + // the indices m=0..d_v-1 are from the center vertex. their order + // is the same as for the first vertex of the whole cell, so we + // can use the system_to_base_table variable (using the + // face_s_t_base_t function would yield the same) if (m < this->dofs_per_vertex) m_index = this->system_to_base_table[m]; else - // then come the two sets of - // line indices + // then come the two sets of line indices { const unsigned int index_in_line = (m-this->dofs_per_vertex) % this->dofs_per_line; @@ -1665,49 +1664,21 @@ void FESystem::build_interface_constraints () = (m-this->dofs_per_vertex) / this->dofs_per_line; Assert (sub_line < 2, ExcInternalError()); - // from this - // information, try - // to get base - // element and - // instance of base - // element. we do - // so by - // constructing the - // corresponding - // face index of m - // in the present - // element, then - // use + // from this information, try to get base element and instance + // of base element. we do so by constructing the corresponding + // face index of m in the present element, then use // face_system_to_base_table const unsigned int tmp1 = 2*this->dofs_per_vertex+index_in_line; m_index.first = this->face_system_to_base_table[tmp1].first; - // what we are - // still missing is - // the index of m - // within the base - // elements - // interface_constraints - // table + // what we are still missing is the index of m within the base + // elements interface_constraints table // - // here, the second - // value of - // face_system_to_base_table - // can help: it - // denotes the face - // index of that - // shape function - // within the base - // element. since - // we know that it - // is a line dof, - // we can construct - // the rest: tmp2 - // will denote the - // index of this - // shape function - // among the line - // shape functions: + // here, the second value of face_system_to_base_table can + // help: it denotes the face index of that shape function + // within the base element. since we know that it is a line + // dof, we can construct the rest: tmp2 will denote the index + // of this shape function among the line shape functions: Assert (this->face_system_to_base_table[tmp1].second >= 2*base_element(m_index.first.first).dofs_per_vertex, ExcInternalError()); @@ -1724,31 +1695,18 @@ void FESystem::build_interface_constraints () case 3: { - // same way as above, - // although a little - // more complicated... - - // the indices - // m=0..5*d_v-1 are - // from the center and - // the four subline - // vertices. their - // order is the same as - // for the first vertex - // of the whole cell, - // so we can use the - // simple arithmetic + // same way as above, although a little more complicated... + + // the indices m=0..5*d_v-1 are from the center and the four + // subline vertices. their order is the same as for the first + // vertex of the whole cell, so we can use the simple arithmetic if (m < 5*this->dofs_per_vertex) m_index = this->system_to_base_table[m]; else - // then come the 12 sets of - // line indices + // then come the 12 sets of line indices if (m < 5*this->dofs_per_vertex + 12*this->dofs_per_line) { - // for the - // meaning of all - // this, see the - // 2d part + // for the meaning of all this, see the 2d part const unsigned int index_in_line = (m-5*this->dofs_per_vertex) % this->dofs_per_line; const unsigned int sub_line @@ -1772,10 +1730,7 @@ void FESystem::build_interface_constraints () else // on one of the four sub-quads { - // for the - // meaning of all - // this, see the - // 2d part + // for the meaning of all this, see the 2d part const unsigned int index_in_quad = (m-5*this->dofs_per_vertex-12*this->dofs_per_line) % this->dofs_per_quad; @@ -1815,14 +1770,9 @@ void FESystem::build_interface_constraints () Assert (false, ExcNotImplemented()); }; - // now that we gathered all - // information: use it to - // build the matrix. note - // that if n and m belong to - // different base elements or - // instances, then there - // definitely will be no - // coupling + // now that we gathered all information: use it to build the + // matrix. note that if n and m belong to different base elements or + // instances, then there definitely will be no coupling if (n_index.first == m_index.first) this->interface_constraints(m,n) = (base_element(n_index.first.first).constraints()(m_index.second, @@ -1869,106 +1819,13 @@ void FESystem::initialize (const std::vector::cut_x; - ref_case::isotropic_refinement+1; - ++ref_case) - { - // Check if some of the matrices of - // the base elements are void. - // repeat this check for each RefineCase - bool do_restriction = true; - bool do_prolongation = true; - - for (unsigned int i=0; in_base_elements(); ++i) - { - if (base_element(i).restriction[ref_case-1][0].n() != base_element(i).dofs_per_cell) - do_restriction = false; - if (base_element(i).prolongation[ref_case-1][0].n() != base_element(i).dofs_per_cell) - do_prolongation = false; - } + // restriction and prolongation matrices are build on demand - // if we did not encounter void - // matrices, initialize the - // respective matrix sizes - if (do_restriction) - for (unsigned int i=0; i::n_children(RefinementCase(ref_case)); ++i) - this->restriction[ref_case-1][i].reinit(this->dofs_per_cell, - this->dofs_per_cell); - if (do_prolongation) - for (unsigned int i=0; i::n_children(RefinementCase(ref_case)); ++i) - this->prolongation[ref_case-1][i].reinit(this->dofs_per_cell, - this->dofs_per_cell); - - // distribute the matrices of the - // base finite elements to the - // matrices of this object. for - // this, loop over all degrees of - // freedom and take the respective - // entry of the underlying base - // element. - // - // note that we by definition of a - // base element, they are - // independent, i.e. do not - // couple. only DoFs that belong to - // the same instance of a base - // element may couple - for (unsigned int i=0; idofs_per_cell; ++i) - for (unsigned int j=0; jdofs_per_cell; ++j) - { - // first find out to which - // base element indices i and - // j belong, and which - // instance thereof in case - // the base element has a - // multiplicity greater than - // one. if they should not - // happen to belong to the - // same instance of a base - // element, then they cannot - // couple, so go on with the - // next index - if (this->system_to_base_table[i].first != - this->system_to_base_table[j].first) - continue; - - // so get the common base - // element and the indices - // therein: - const unsigned int - base = this->system_to_base_table[i].first.first; - - const unsigned int - base_index_i = this->system_to_base_table[i].second, - base_index_j = this->system_to_base_table[j].second; - - // if we are sure that DoFs i - // and j may couple, then - // copy entries of the - // matrices: - for (unsigned int child=0; child::n_children(RefinementCase(ref_case)); ++child) - { - if (do_restriction) - this->restriction[ref_case-1][child] (i,j) - = (base_element(base) - .get_restriction_matrix(child, RefinementCase(ref_case))( - base_index_i,base_index_j)); - - if (do_prolongation) - this->prolongation[ref_case-1][child] (i,j) - = (base_element(base) - .get_prolongation_matrix(child, RefinementCase(ref_case))( - base_index_i,base_index_j)); - }; - }; - } - // now set up the interface constraints. - // this is kind'o hairy, so don't try + // now set up the interface constraints. this is kind'o hairy, so don't try // to do it dimension independent build_interface_constraints (); - // finally fill in support points - // on cell and face + // finally fill in support points on cell and face initialize_unit_support_points (); initialize_unit_face_support_points (); @@ -1985,14 +1842,9 @@ void FESystem:: initialize_unit_support_points () { - // if one of the base elements - // has no support points, then - // it makes no sense to define - // support points for the - // composed element, so return - // an empty array to - // demonstrate that - // fact + // if one of the base elements has no support points, then it makes no sense + // to define support points for the composed element, so return an empty + // array to demonstrate that fact for (unsigned int base_el=0; base_eln_base_elements(); ++base_el) if (!base_element(base_el).has_support_points()) { @@ -2000,9 +1852,7 @@ initialize_unit_support_points () return; }; - // generate unit support points - // from unit support points of sub - // elements + // generate unit support points from unit support points of sub elements this->unit_support_points.resize(this->dofs_per_cell); for (unsigned int i=0; idofs_per_cell; ++i) @@ -2028,22 +1878,15 @@ initialize_unit_face_support_points () if (dim == 1) return; - // if one of the base elements has no - // support points, then it makes no sense - // to define support points for the - // composed element. In that case, return - // an empty array to demonstrate that fact - // (note that we ask whether the base - // element has no support points at all, - // not only none on the face!) + // if one of the base elements has no support points, then it makes no sense + // to define support points for the composed element. In that case, return + // an empty array to demonstrate that fact (note that we ask whether the + // base element has no support points at all, not only none on the face!) // - // on the other hand, if there is an - // element that simply has no degrees of - // freedom on the face at all, then we - // don't care whether it has support points - // or not. this is, for example, the case - // for the stable Stokes element Q(p)^dim - // \times DGP(p-1). + // on the other hand, if there is an element that simply has no degrees of + // freedom on the face at all, then we don't care whether it has support + // points or not. this is, for example, the case for the stable Stokes + // element Q(p)^dim \times DGP(p-1). for (unsigned int base_el=0; base_eln_base_elements(); ++base_el) if (!base_element(base_el).has_support_points() && @@ -2054,8 +1897,7 @@ initialize_unit_face_support_points () } - // generate unit face support points - // from unit support points of sub + // generate unit face support points from unit support points of sub // elements this->unit_face_support_points.resize(this->dofs_per_face); @@ -2082,14 +1924,13 @@ FESystem::initialize_quad_dof_index_permutation () if (dim < 3) return; - // the array into which we want to write - // should have the correct size already. + // the array into which we want to write should have the correct size + // already. Assert (this->adjust_quad_dof_index_for_face_orientation_table.n_elements()== 8*this->dofs_per_quad, ExcInternalError()); - // to obtain the shifts for this composed - // element, copy the shift information of - // the base elements + // to obtain the shifts for this composed element, copy the shift + // information of the base elements unsigned int index = 0; for (unsigned int b=0; bn_base_elements(); ++b) { @@ -2107,8 +1948,7 @@ FESystem::initialize_quad_dof_index_permutation () Assert (index == this->dofs_per_quad, ExcInternalError()); - // aditionally compose the permutation - // information for lines + // aditionally compose the permutation information for lines Assert (this->adjust_line_dof_index_for_line_orientation_table.size()== this->dofs_per_line, ExcInternalError()); index = 0; @@ -2164,31 +2004,23 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, ExcDimensionMismatch (interpolation_matrix.m(), x_source_fe.dofs_per_face)); - // since dofs for each base are - // independent, we only have to stack - // things up from base element to base - // element + // since dofs for each base are independent, we only have to stack things up + // from base element to base element // - // the problem is that we have to work with - // two FEs (this and fe_other). only deal - // with the case that both are FESystems - // and that they both have the same number - // of bases (counting multiplicity) each of - // which match in their number of - // components. this covers - // FESystem(FE_Q(p),1,FE_Q(q),2) vs - // FESystem(FE_Q(r),2,FE_Q(s),1), but not + // the problem is that we have to work with two FEs (this and + // fe_other). only deal with the case that both are FESystems and that they + // both have the same number of bases (counting multiplicity) each of which + // match in their number of components. this covers + // FESystem(FE_Q(p),1,FE_Q(q),2) vs FESystem(FE_Q(r),2,FE_Q(s),1), but not // FESystem(FE_Q(p),1,FE_Q(q),2) vs // FESystem(FESystem(FE_Q(r),2),1,FE_Q(s),1) const FESystem *fe_other_system = dynamic_cast*>(&x_source_fe); - // clear matrix, since we will not get to - // set all elements + // clear matrix, since we will not get to set all elements interpolation_matrix = 0; - // loop over all the base elements of - // this and the other element, counting + // loop over all the base elements of this and the other element, counting // their multiplicities unsigned int base_index = 0, base_index_other = 0; @@ -2212,11 +2044,9 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, base.get_face_interpolation_matrix (base_other, base_to_base_interpolation); - // now translate entries. we'd like to - // have something like - // face_base_to_system_index, but that - // doesn't exist. rather, all we have - // is the reverse. well, use that then + // now translate entries. we'd like to have something like + // face_base_to_system_index, but that doesn't exist. rather, all we + // have is the reverse. well, use that then for (unsigned int i=0; idofs_per_face; ++i) if (this->face_system_to_base_index(i).first == @@ -2229,12 +2059,9 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, = base_to_base_interpolation(fe_other_system->face_system_to_base_index(j).second, this->face_system_to_base_index(i).second); - // advance to the next base element - // for this and the other - // fe_system; see if we can simply - // advance the multiplicity by one, - // or if have to move on to the - // next base element + // advance to the next base element for this and the other fe_system; + // see if we can simply advance the multiplicity by one, or if have to + // move on to the next base element ++multiplicity; if (multiplicity == this->element_multiplicity(base_index)) { @@ -2249,10 +2076,8 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, ++base_index_other; } - // see if we have reached the end - // of the present element. if so, - // we should have reached the end - // of the other one as well + // see if we have reached the end of the present element. if so, we + // should have reached the end of the other one as well if (base_index == this->n_base_elements()) { Assert (base_index_other == fe_other_system->n_base_elements(), @@ -2260,10 +2085,8 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, break; } - // if we haven't reached the end of - // this element, we shouldn't have - // reached the end of the other one - // either + // if we haven't reached the end of this element, we shouldn't have + // reached the end of the other one either Assert (base_index_other != fe_other_system->n_base_elements(), ExcInternalError()); } @@ -2292,31 +2115,23 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe ExcDimensionMismatch (interpolation_matrix.m(), x_source_fe.dofs_per_face)); - // since dofs for each base are - // independent, we only have to stack - // things up from base element to base - // element + // since dofs for each base are independent, we only have to stack things up + // from base element to base element // - // the problem is that we have to work with - // two FEs (this and fe_other). only deal - // with the case that both are FESystems - // and that they both have the same number - // of bases (counting multiplicity) each of - // which match in their number of - // components. this covers - // FESystem(FE_Q(p),1,FE_Q(q),2) vs - // FESystem(FE_Q(r),2,FE_Q(s),1), but not + // the problem is that we have to work with two FEs (this and + // fe_other). only deal with the case that both are FESystems and that they + // both have the same number of bases (counting multiplicity) each of which + // match in their number of components. this covers + // FESystem(FE_Q(p),1,FE_Q(q),2) vs FESystem(FE_Q(r),2,FE_Q(s),1), but not // FESystem(FE_Q(p),1,FE_Q(q),2) vs // FESystem(FESystem(FE_Q(r),2),1,FE_Q(s),1) const FESystem *fe_other_system = dynamic_cast*>(&x_source_fe); - // clear matrix, since we will not get to - // set all elements + // clear matrix, since we will not get to set all elements interpolation_matrix = 0; - // loop over all the base elements of - // this and the other element, counting + // loop over all the base elements of this and the other element, counting // their multiplicities unsigned int base_index = 0, base_index_other = 0; @@ -2341,11 +2156,9 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe subface, base_to_base_interpolation); - // now translate entries. we'd like to - // have something like - // face_base_to_system_index, but that - // doesn't exist. rather, all we have - // is the reverse. well, use that then + // now translate entries. we'd like to have something like + // face_base_to_system_index, but that doesn't exist. rather, all we + // have is the reverse. well, use that then for (unsigned int i=0; idofs_per_face; ++i) if (this->face_system_to_base_index(i).first == @@ -2358,12 +2171,9 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe = base_to_base_interpolation(fe_other_system->face_system_to_base_index(j).second, this->face_system_to_base_index(i).second); - // advance to the next base element - // for this and the other - // fe_system; see if we can simply - // advance the multiplicity by one, - // or if have to move on to the - // next base element + // advance to the next base element for this and the other fe_system; + // see if we can simply advance the multiplicity by one, or if have to + // move on to the next base element ++multiplicity; if (multiplicity == this->element_multiplicity(base_index)) { @@ -2378,10 +2188,8 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe ++base_index_other; } - // see if we have reached the end - // of the present element. if so, - // we should have reached the end - // of the other one as well + // see if we have reached the end of the present element. if so, we + // should have reached the end of the other one as well if (base_index == this->n_base_elements()) { Assert (base_index_other == fe_other_system->n_base_elements(), @@ -2389,10 +2197,8 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe break; } - // if we haven't reached the end of - // this element, we shouldn't have - // reached the end of the other one - // either + // if we haven't reached the end of this element, we shouldn't have + // reached the end of the other one either Assert (base_index_other != fe_other_system->n_base_elements(), ExcInternalError()); } @@ -2405,38 +2211,29 @@ template std::vector > FESystem::hp_object_dof_identities (const FiniteElement &fe_other) const { - // since dofs on each subobject (vertex, - // line, ...) are ordered such that first - // come all from the first base element all - // multiplicities, then second base element - // all multiplicities, etc., we simply have - // to stack all the identities after each - // other + // since dofs on each subobject (vertex, line, ...) are ordered such that + // first come all from the first base element all multiplicities, then + // second base element all multiplicities, etc., we simply have to stack all + // the identities after each other // - // the problem is that we have to work with - // two FEs (this and fe_other). only deal - // with the case that both are FESystems - // and that they both have the same number - // of bases (counting multiplicity) each of - // which match in their number of - // components. this covers - // FESystem(FE_Q(p),1,FE_Q(q),2) vs - // FESystem(FE_Q(r),2,FE_Q(s),1), but not + // the problem is that we have to work with two FEs (this and + // fe_other). only deal with the case that both are FESystems and that they + // both have the same number of bases (counting multiplicity) each of which + // match in their number of components. this covers + // FESystem(FE_Q(p),1,FE_Q(q),2) vs FESystem(FE_Q(r),2,FE_Q(s),1), but not // FESystem(FE_Q(p),1,FE_Q(q),2) vs // FESystem(FESystem(FE_Q(r),2),1,FE_Q(s),1) if (const FESystem *fe_other_system = dynamic_cast*>(&fe_other)) { - // loop over all the base elements of - // this and the other element, counting - // their multiplicities + // loop over all the base elements of this and the other element, + // counting their multiplicities unsigned int base_index = 0, base_index_other = 0; unsigned int multiplicity = 0, multiplicity_other = 0; - // we also need to keep track of the - // number of dofs already treated for + // we also need to keep track of the number of dofs already treated for // each of the elements unsigned int dof_offset = 0, dof_offset_other = 0; @@ -2452,10 +2249,8 @@ FESystem::hp_object_dof_identities (const FiniteElement > base_identities; switch (structdim) { @@ -2477,17 +2272,13 @@ FESystem::hp_object_dof_identities (const FiniteElement(); dof_offset_other += base_other.template n_dofs_per_object(); - // advance to the next base element - // for this and the other - // fe_system; see if we can simply - // advance the multiplicity by one, - // or if have to move on to the - // next base element + // advance to the next base element for this and the other + // fe_system; see if we can simply advance the multiplicity by one, + // or if have to move on to the next base element ++multiplicity; if (multiplicity == this->element_multiplicity(base_index)) { @@ -2502,10 +2293,8 @@ FESystem::hp_object_dof_identities (const FiniteElementn_base_elements()) { Assert (base_index_other == fe_other_system->n_base_elements(), @@ -2513,10 +2302,8 @@ FESystem::hp_object_dof_identities (const FiniteElementn_base_elements(), ExcInternalError()); } @@ -2562,9 +2349,8 @@ FiniteElementDomination::Domination FESystem:: compare_for_face_domination (const FiniteElement &fe_other) const { - // at present all we can do is to compare - // with other FESystems that have the same - // number of components and bases + // at present all we can do is to compare with other FESystems that have the + // same number of components and bases if (const FESystem *fe_sys_other = dynamic_cast*>(&fe_other)) { @@ -2576,8 +2362,7 @@ compare_for_face_domination (const FiniteElement &fe_other) const FiniteElementDomination::Domination domination = FiniteElementDomination::no_requirements; - // loop over all base elements and do - // some sanity checks + // loop over all base elements and do some sanity checks for (unsigned int b=0; bn_base_elements(); ++b) { Assert (this->base_element(b).n_components() == @@ -2587,8 +2372,7 @@ compare_for_face_domination (const FiniteElement &fe_other) const fe_sys_other->element_multiplicity(b), ExcNotImplemented()); - // for this pair of base elements, - // check who dominates and combine + // for this pair of base elements, check who dominates and combine // with previous result const FiniteElementDomination::Domination base_domination = (this->base_element(b) @@ -2596,8 +2380,7 @@ compare_for_face_domination (const FiniteElement &fe_other) const domination = domination & base_domination; } - // if we've gotten here, then we've - // either found a winner or either + // if we've gotten here, then we've either found a winner or either // element is fine being dominated Assert (domination != FiniteElementDomination::neither_element_dominates, @@ -2756,34 +2539,23 @@ compute_restriction_is_additive_flags (const std::vector0) // check needed as fe might be NULL n_shape_functions += fes[i]->dofs_per_cell * multiplicities[i]; - // generate the array that will - // hold the output + // generate the array that will hold the output std::vector retval (n_shape_functions, false); - // finally go through all the shape - // functions of the base elements, - // and copy their flags. this - // somehow copies the code in - // build_cell_table, which is not - // nice as it uses too much - // implicit knowledge about the - // layout of the individual bases - // in the composed FE, but there - // seems no way around... + // finally go through all the shape functions of the base elements, and copy + // their flags. this somehow copies the code in build_cell_table, which is + // not nice as it uses too much implicit knowledge about the layout of the + // individual bases in the composed FE, but there seems no way around... // - // for each shape function, copy - // the flags from the base element - // to this one, taking into account - // multiplicities, and other - // complications + // for each shape function, copy the flags from the base element to this + // one, taking into account multiplicities, and other complications unsigned int total_index = 0; for (unsigned int vertex_number=0; vertex_number::vertices_per_cell; @@ -2925,9 +2697,8 @@ compute_nonzero_components (const std::vector* { Assert (fes.size() == multiplicities.size(), ExcInternalError()); - // first count the number of dofs - // and components that will emerge - // from the given FEs + // first count the number of dofs and components that will emerge from the + // given FEs unsigned int n_shape_functions = 0; for (unsigned int i=0; i0) //needed because fe might be NULL @@ -2938,28 +2709,18 @@ compute_nonzero_components (const std::vector* if (multiplicities[i]>0) //needed because fe might be NULL n_components += fes[i]->n_components() * multiplicities[i]; - // generate the array that will - // hold the output + // generate the array that will hold the output std::vector > retval (n_shape_functions, std::vector (n_components, false)); - // finally go through all the shape - // functions of the base elements, - // and copy their flags. this - // somehow copies the code in - // build_cell_table, which is not - // nice as it uses too much - // implicit knowledge about the - // layout of the individual bases - // in the composed FE, but there - // seems no way around... + // finally go through all the shape functions of the base elements, and copy + // their flags. this somehow copies the code in build_cell_table, which is + // not nice as it uses too much implicit knowledge about the layout of the + // individual bases in the composed FE, but there seems no way around... // - // for each shape function, copy - // the non-zero flags from the base - // element to this one, taking into - // account multiplicities, multiple - // components in base elements, and - // other complications + // for each shape function, copy the non-zero flags from the base element to + // this one, taking into account multiplicities, multiple components in base + // elements, and other complications unsigned int total_index = 0; for (unsigned int vertex_number=0; vertex_number::vertices_per_cell; @@ -3143,14 +2904,12 @@ FESystem::unit_support_point (const unsigned int index) const (this->unit_support_points.size() == 0), typename FEL::ExcFEHasNoSupportPoints ()); - // let's see whether we have the - // information pre-computed + // let's see whether we have the information pre-computed if (this->unit_support_points.size() != 0) return this->unit_support_points[index]; else - // no. ask the base element - // whether it would like to - // provide this information + // no. ask the base element whether it would like to provide this + // information return (base_element(this->system_to_base_index(index).first.first) .unit_support_point(this->system_to_base_index(index).second)); } @@ -3168,14 +2927,12 @@ FESystem::unit_face_support_point (const unsigned int index) const (this->unit_face_support_points.size() == 0), typename FEL::ExcFEHasNoSupportPoints ()); - // let's see whether we have the - // information pre-computed + // let's see whether we have the information pre-computed if (this->unit_face_support_points.size() != 0) return this->unit_face_support_points[index]; else - // no. ask the base element - // whether it would like to - // provide this information + // no. ask the base element whether it would like to provide this + // information return (base_element(this->face_system_to_base_index(index).first.first) .unit_face_support_point(this->face_system_to_base_index(index).second)); } @@ -3186,11 +2943,8 @@ template std::size_t FESystem::memory_consumption () const { - // neglect size of data stored in - // @p{base_elements} due to some - // problems with the - // compiler. should be neglectable - // after all, considering the size + // neglect size of data stored in @p{base_elements} due to some problems + // with the compiler. should be neglectable after all, considering the size // of the data of the subelements std::size_t mem = (FiniteElement::memory_consumption () + sizeof (base_elements)); -- 2.39.5