From 51cd5aeb60442dc4c0d2410d0257d460ce7c1aa9 Mon Sep 17 00:00:00 2001 From: guido Date: Thu, 25 Aug 2005 10:45:31 +0000 Subject: [PATCH] FiniteElementBase removed git-svn-id: https://svn.dealii.org/trunk@11326 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/dofs/dof_tools.h | 2 +- deal.II/deal.II/include/fe/fe.h | 1723 +++++++++++++++- deal.II/deal.II/include/fe/fe_base.h | 1746 +---------------- deal.II/deal.II/include/fe/fe_dgp_monomial.h | 4 +- .../deal.II/include/fe/fe_dgp_nonparametric.h | 14 +- deal.II/deal.II/include/fe/fe_dgq.h | 4 +- deal.II/deal.II/include/fe/fe_nedelec.h | 12 +- deal.II/deal.II/include/fe/fe_poly.h | 16 +- deal.II/deal.II/include/fe/fe_poly_tensor.h | 2 +- deal.II/deal.II/include/fe/fe_q.h | 8 +- .../deal.II/include/fe/fe_q_hierarchical.h | 4 +- .../deal.II/include/fe/fe_raviart_thomas.h | 22 +- deal.II/deal.II/include/fe/fe_system.h | 22 +- deal.II/deal.II/include/fe/fe_tools.h | 4 +- deal.II/deal.II/source/dofs/dof_tools.cc | 4 +- deal.II/deal.II/source/fe/fe.cc | 143 +- deal.II/deal.II/source/fe/fe_dgp_monomial.cc | 2 +- deal.II/deal.II/source/fe/fe_dgq.cc | 4 +- deal.II/deal.II/source/fe/fe_poly_tensor.cc | 6 +- deal.II/deal.II/source/fe/fe_q.cc | 8 +- .../deal.II/source/fe/fe_raviart_thomas.cc | 6 +- deal.II/deal.II/source/fe/fe_system.cc | 30 +- deal.II/deal.II/source/fe/fe_tools.cc | 4 +- deal.II/doc/news/changes.html | 7 + 24 files changed, 1868 insertions(+), 1929 deletions(-) diff --git a/deal.II/deal.II/include/dofs/dof_tools.h b/deal.II/deal.II/include/dofs/dof_tools.h index 520b77e004..0e5245033c 100644 --- a/deal.II/deal.II/include/dofs/dof_tools.h +++ b/deal.II/deal.II/include/dofs/dof_tools.h @@ -130,7 +130,7 @@ template class Mapping; * this is the point where it has the function value 1), is * located on the boundary. We do not check this directly, the * criterion is rather defined through the information the finite - * element class gives: the FiniteElementBase class defines the + * element class gives: the FiniteElement class defines the * numbers of basis functions per vertex, per line, and so on and the * basis functions are numbered after this information; a basis * function is to be considered to be on the face of a cell (and thus diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 6c5e66119b..817f1f9178 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -30,9 +30,7 @@ template class FECollection; /** * Common interface of all finite elements. Here, the functions to - * fill the data fields of FEValues are declared. While - * FiniteElementBase provides implementation of common - * functionality, this class only serves as an abstract base class. + * fill the data fields of FEValues are declared. * * The interface of this class is very restrictive. The reason is that * finite element values should be accessed only by use of @@ -54,8 +52,59 @@ template class FECollection; * @author Wolfgang Bangerth, Guido Kanschat, Ralf Hartmann, 1998, 2000, 2001 */ template -class FiniteElement : public FiniteElementBase +class FiniteElement : public Subscriptor, + public FiniteElementData { + public: + /** + * 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. + * + * @author Guido Kanschat, 2001 + */ + class InternalDataBase : public Mapping::InternalDataBase + { + public: + /** + * 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. + */ + void initialize_2nd (const FiniteElement *element, + const Mapping &mapping, + const Quadrature &quadrature); + + /** + * Storage for @p FEValues + * objects needed to + * approximate second + * derivatives. + * + * The ordering is p+hx, + * p+hy, p+hz, + * @p p-hx, @p p-hy, + * @p p-hz, where unused + * entries in lower dimensions + * are missing. + */ + std::vector*> differences; + }; + public: /** * Copy constructor. This one is declared @@ -87,6 +136,694 @@ class FiniteElement : public FiniteElementBase */ virtual ~FiniteElement (); + /** + * Return a string that uniquely + * identifies a finite + * element. The general + * convention is that this is the + * class name, followed by the + * space 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. + */ + virtual std::string get_name () const = 0; + + /** + * @name Shape function access + * @{ + */ + + /** + * 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 + * @p ExcShapeFunctionNotPrimitive. In + * that case, use the + * shape_value_component() + * function. + * + * An + * @p ExcUnitShapeValuesDoNotExist + * is thrown if the shape values + * of the @p 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 @p 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 + * @p ExcShapeFunctionNotPrimitive. In + * that case, use the + * shape_grad_component() + * function. + * + * An + * @p ExcUnitShapeValuesDoNotExist + * is thrown if the shape values + * of the @p 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 @p 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 + * @p ExcShapeFunctionNotPrimitive. In + * that case, use the + * shape_grad_grad_component() + * function. + * + * An + * @p ExcUnitShapeValuesDoNotExist + * is thrown if the shape values + * of the @p 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; + + /** + * Just like for @p shape_grad_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<2,dim> shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + /** + * 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. + */ + virtual bool has_support_on_face (const unsigned int shape_index, + const unsigned int face_index) const = 0; + + //@} + /** + * @name Transfer and constraint matrices + * @{ + */ + + /** + * 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 + * @p restriction_is_additive_flags + * given to the constructor. 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 + * @p ExcProjectionVoid. + */ + const FullMatrix & + get_restriction_matrix (const unsigned int child) 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 + * @p ExcEmbeddingVoid. You can + * check whether this is the case + * by calling the + * prolongation_is_implemented(). + */ + const FullMatrix & + get_prolongation_matrix (const unsigned int child) const; + + /** + * Return whether this element implements + * its prolongation matrices. 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 prolongation_is_implemented () const; + + /** + * Return whether this element implements + * its restriction matrices. 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 restriction_is_implemented () const; + + /** + * Access the + * @p restriction_is_additive_flag + * field. See there for more + * information on its contents. + * + * 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 readonly 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 space 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 + * @p constraints_are_implemented + * function to check up front + * whethehr this function will + * succeed or generate the + * exception. + */ + const FullMatrix & constraints () const; + + /** + * Return whether this element + * implements its hanging node + * constraints. The return value + * also indicates whether a call + * to the @p constraint 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 + * just expresses. + */ + bool constraints_are_implemented () 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 + * 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 + * FiniteElement::ExcInterpolationNotImplemented. + */ + virtual void + get_interpolation_matrix (const FiniteElement &source, + FullMatrix &matrix) const; + //@} + + /** + * 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 + * @p restriction and + * @p prolongation. + */ + bool operator == (const FiniteElement &) const; + + /** + * @name Index computations + * @{ + */ + /** + * 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 + * @p 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 + * @p system_to_base_table field + * and the + * @p system_to_base_index + * function. + */ + std::pair + system_to_component_index (const unsigned int index) const; + + /** + * 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. + * + * This is the opposite operation + * from the @p system_to_component_index + * function. + */ + unsigned int component_to_system_index(const unsigned int component, + const unsigned int index) const; + + /** + * Same as above, but do it for + * shape functions and their + * indices on a face. + */ + std::pair + face_system_to_component_index (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. + * + * 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 + * @p 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 + * @p system_to_component_index + * function. + */ + std::pair, unsigned int> + system_to_base_index (const unsigned int index) const; + + /** + * Same as + * @p system_to_base_index, but + * for degrees of freedom located + * on a face. + */ + std::pair, unsigned int> + face_system_to_base_index (const unsigned int index) const; + + /** + * Return in which of the vector + * components of this finite + * element the @p ithe 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(i) + * 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. + */ + const std::vector & + 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 + * @p 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. + */ + 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 + * 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. + */ + bool + is_primitive (const unsigned int i) const; + + /** + * Return whether the entire + * finite element is primitive, + * in the sense that all its + * shape functions are + * primitive. If the finite + * element is scalar, then this + * is always the case. + * + * Since this is an extremely + * common operation, the result + * is cached in the + * @p cached_primitivity + * variable which is computed in + * the constructor. + */ + bool + is_primitive () const; + /** * Number of base elements in a * mixed discretization. @@ -110,7 +847,7 @@ class FiniteElement : public FiniteElementBase * multiplicity. */ virtual unsigned int n_base_elements () const = 0; - + /** * Access to base element * objects. If the element is @@ -136,17 +873,293 @@ class FiniteElement : public FiniteElementBase unsigned int element_multiplicity (const unsigned int index) const = 0; + /** + * Given a vector component, + * return an index 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::@p base_element + * function. + * + * If this is a scalar finite + * element, then the return value + * is always equal to a pair of + * zeros. + */ + std::pair + component_to_base (const unsigned int component) const; + //@} + /** - * Check for non-zero values on a face. + * @name Support points and interpolation + * @{ + */ + + /** + * 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. * - * This function returns - * @p true, if the shape - * function @p shape_index has - * non-zero values on the face - * @p face_index. + * 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. */ - virtual bool has_support_on_face (const unsigned int shape_index, - const unsigned int face_index) const = 0; + 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 + * @p 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. + */ + 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. + * + * 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, of not all the + * base elements define support + * points. In this way, you can + * still ask for certain support + * points, even if + * @p 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 + * (@p 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 + * @p get_unit_support_points + * yields a non-empty array. + * + * 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. + */ + virtual + Point + unit_face_support_point (const unsigned int index) const; + + /** + * Return a support point vector + * for generalized interpolation. + */ + const std::vector > & + get_generalized_support_points () const; + + /** + * + */ + bool has_generalized_support_points () const; + + /** + * + */ + const std::vector > & + get_generalized_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 + * @p get_unit_support_points + * yields a non-empty array. + * + * 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. + * + * @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. + * + * 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. + */ + virtual void interpolate( + std::vector& local_dofs, + const VectorSlice > >& values) const; + + //@} /** * Determine an estimate for the @@ -160,6 +1173,474 @@ class FiniteElement : public FiniteElementBase * rather than the class itself. */ virtual unsigned int memory_consumption () const; + /** + * Exception + * + * @ingroup Exceptions + */ + DeclException1 (ExcShapeFunctionNotPrimitive, + int, + << "The shape function with index " << arg1 + << " is not primitive, i.e. it is vector-valued and " + << "has more than one non-zero vector component. This " + << "function cannot be called for these shape functions. " + << "Maybe you want to use the same function with the " + << "_component suffix?"); + /** + * Exception + * + * @ingroup Exceptions + */ + DeclException0 (ExcFENotPrimitive); + /** + * Exception + * + * @ingroup Exceptions + */ + DeclException0 (ExcUnitShapeValuesDoNotExist); + + /** + * 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. + * + * @ingroup Exceptions + */ + DeclException0 (ExcEmbeddingVoid); + + /** + * Attempt to access restriction + * matrices of a finite element + * which did not implement these + * matrices. + * + * Exception + * @ingroup Exceptions + */ + DeclException0 (ExcProjectionVoid); + + /** + * Attempt to access constraint + * matrices of a finite element + * which did not implement these + * matrices. + * + * Exception + * @ingroup Exceptions + */ + DeclException0 (ExcConstraintsVoid); + + /** + * Exception + * @ingroup Exceptions + */ + DeclException2 (ExcWrongInterfaceMatrixSize, + int, int, + << "The interface matrix has a size of " << arg1 + << "x" << arg2 + << ", which is not reasonable in the present dimension."); + /** + * Exception + * @ingroup Exceptions + */ + DeclException2 (ExcComponentIndexInvalid, + int, int, + << "The component-index pair (" << arg1 << ", " << arg2 + << ") is invalid, i.e. non-existent"); + /** + * Exception + * @ingroup Exceptions + */ + DeclException0 (ExcInterpolationNotImplemented); + + protected: + /** + * Array of projection matrices. See + * get_restriction_matrix() above. + * + * Matrices in this array are + * automatically initialized to + * correct size. If the derived + * finite element class does not + * implement these matrices, they + * should be resized to zero + * size. + */ + FullMatrix restriction[GeometryInfo::children_per_cell]; + + /** + * Array of embedding matrices. See + * get_prolongation_matrix() above. + * + * Matrices in this array are + * automatically initialized to + * correct size. If the derived + * finite element class does not + * implement these matrices, they + * should be resized to zero + * size. + */ + FullMatrix prolongation[GeometryInfo::children_per_cell]; + + /** + * 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. + * + * This field is obviously + * useless in one space dimension + * and has there a zero size. + */ + FullMatrix interface_constraints; + + /** + * 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; + + /** + * Store what + * @p 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. + * + * 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 + * @p 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 + * @p system_to_component_table. + */ + std::vector,unsigned int> > + system_to_base_table; + + /** + * Likewise for the indices on + * faces. + */ + std::vector,unsigned int> > + face_system_to_base_table; + + /** + * The base element establishing + * a component. + * + * This table converts a + * component number to a pair + * consisting of the + * @p base_element number, and + * the component within this base + * element. While component + * information contains + * multiplicity of base elements, + * the result allows access to + * shape functions of the base + * element. + * + * 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 > 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 + * @p 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; + + /** + * 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 + * @p 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. + */ + std::vector > generalized_support_points; + + /** + * Face support points used for + * interpolation functions of + * non-Lagrangian elements. + */ + std::vector > generalized_face_support_points; + + /** + * 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. + */ + const std::vector > nonzero_components; + + /** + * This array holds how many + * values in the respective entry + * of the @p 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; + + /** + * Store whether all shape + * functions are primitive. Since + * finding this out is a very + * common operation, we cache the + * result, i.e. compute the value + * in the constructor for simpler + * access. + */ + const bool cached_primitivity; + + /** + * Compute second derivatives by + * finite differences of + * gradients. + */ + void compute_2nd (const Mapping &mapping, + const typename Triangulation::cell_iterator &cell, + const unsigned int offset, + typename Mapping::InternalDataBase &mapping_internal, + InternalDataBase &fe_internal, + 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. + */ + static + std::vector + compute_n_nonzero_components (const std::vector > &nonzero_components); + + /** + * Allow the FESystem class to access the + * restriction and prolongation matrices + * directly. Hence, FESystem has the + * possibility to see if these matrices + * are initialized without accessing + * these matrices through the + * @p get_restriction_matrix and + * @p get_prolongation_matrix + * functions. This is important as these + * functions include assertions that + * throw if the matrices are not already + * initialized. + */ + template friend class FESystem; + + /** + * Make the inner class a + * friend. This is not strictly + * necessary, but the Intel + * compiler seems to want this. + */ + friend class InternalDataBase; + /** * Exception @@ -244,6 +1725,20 @@ class FiniteElement : public FiniteElementBase */ virtual FiniteElement *clone() const = 0; + private: + /** + * 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 @@ -367,5 +1862,207 @@ class FiniteElement : public FiniteElementBase }; /*@}*/ +//----------------------------------------------------------------------// + +template +inline +std::pair +FiniteElement::system_to_component_index (const unsigned int index) const +{ + Assert (index < system_to_component_table.size(), + ExcIndexRange(index, 0, system_to_component_table.size())); + Assert (is_primitive (index), + typename FiniteElement::ExcShapeFunctionNotPrimitive(index)); + return system_to_component_table[index]; +} + +template +inline +unsigned int +FiniteElement::component_to_system_index (const unsigned int component, + const unsigned int index) const +{ + std::vector< std::pair >::const_iterator + it = std::find(system_to_component_table.begin(), system_to_component_table.end(), + std::pair(component, index)); + + Assert(it != system_to_component_table.end(), ExcComponentIndexInvalid(component, index)); + return std::distance(system_to_component_table.begin(), it); +} + + + +template +inline +std::pair +FiniteElement::face_system_to_component_index (const unsigned int index) const +{ + Assert(index < face_system_to_component_table.size(), + ExcIndexRange(index, 0, face_system_to_component_table.size())); + + // in debug mode, check whether the + // function is primitive, since + // otherwise the result may have no + // meaning + // + // since the primitivity tables are + // all geared towards cell dof + // indices, rather than face dof + // indices, we have to work a + // little bit... + // + // in 1d, the face index is equal + // to the cell index + Assert (((dim == 1) && is_primitive(index)) + || + // in 2d, construct it like + // this: + ((dim == 2) && + is_primitive (index < (GeometryInfo<2>::vertices_per_face * + this->dofs_per_vertex) + ? + index + : + GeometryInfo<2>::vertices_per_cell * + this->dofs_per_vertex + + (index - + GeometryInfo<2>::vertices_per_face * + this->dofs_per_vertex))) + || + // likewise in 3d, but more + // complicated + ((dim == 3) && + is_primitive (index < (GeometryInfo<3>::vertices_per_face * + this->dofs_per_vertex) + ? + index + : + (index < (GeometryInfo<3>::vertices_per_face * + this->dofs_per_vertex + + + GeometryInfo<3>::lines_per_face * + this->dofs_per_line) + ? + GeometryInfo<3>::vertices_per_cell * + this->dofs_per_vertex + + (index - + GeometryInfo<3>::vertices_per_face * + this->dofs_per_vertex) + : + GeometryInfo<3>::vertices_per_cell * + this->dofs_per_vertex + + GeometryInfo<3>::lines_per_cell * + this->dofs_per_line + + (index - + GeometryInfo<3>::vertices_per_face * + this->dofs_per_vertex + - + GeometryInfo<3>::lines_per_face * + this->dofs_per_line)))), + typename FiniteElement::ExcShapeFunctionNotPrimitive(index)); + + return face_system_to_component_table[index]; +} + + + +template +inline +std::pair,unsigned int> +FiniteElement::system_to_base_index (const unsigned int index) const +{ + Assert (index < system_to_base_table.size(), + ExcIndexRange(index, 0, system_to_base_table.size())); + return system_to_base_table[index]; +} + + + + +template +inline +std::pair,unsigned int> +FiniteElement::face_system_to_base_index (const unsigned int index) const +{ + Assert(index < face_system_to_base_table.size(), + ExcIndexRange(index, 0, face_system_to_base_table.size())); + return face_system_to_base_table[index]; +} + + + +template +inline +std::pair +FiniteElement::component_to_base (const unsigned int index) const +{ + Assert(index < component_to_base_table.size(), + ExcIndexRange(index, 0, component_to_base_table.size())); + + return component_to_base_table[index]; +} + + +template +inline +bool +FiniteElement::restriction_is_additive (const unsigned int index) const +{ + Assert(index < this->dofs_per_cell, + ExcIndexRange(index, 0, this->dofs_per_cell)); + return restriction_is_additive_flags[index]; +} + + +template +inline +const std::vector & +FiniteElement::get_nonzero_components (const unsigned int i) const +{ + Assert (i < this->dofs_per_cell, ExcIndexRange (i, 0, this->dofs_per_cell)); + return nonzero_components[i]; +} + + + +template +inline +unsigned int +FiniteElement::n_nonzero_components (const unsigned int i) const +{ + Assert (i < this->dofs_per_cell, ExcIndexRange (i, 0, this->dofs_per_cell)); + return n_nonzero_components_table[i]; +} + + + +template +inline +bool +FiniteElement::is_primitive (const unsigned int i) const +{ + Assert (i < this->dofs_per_cell, ExcIndexRange (i, 0, this->dofs_per_cell)); + + // return primitivity of a shape + // function by checking whether it + // has more than one non-zero + // component or not. we could cache + // this value in an array of bools, + // but accessing a bit-vector (as + // std::vector is) is + // probably more expensive than + // just comparing against 1 + return (n_nonzero_components_table[i] == 1); +} + + +template +inline +bool +FiniteElement::is_primitive () const +{ + return cached_primitivity; +} + #endif diff --git a/deal.II/deal.II/include/fe/fe_base.h b/deal.II/deal.II/include/fe/fe_base.h index 1e5c1ac93b..33837eefe3 100644 --- a/deal.II/deal.II/include/fe/fe_base.h +++ b/deal.II/deal.II/include/fe/fe_base.h @@ -35,14 +35,14 @@ template class FESystem; /** * Dimension independent data for finite elements. See the derived - * class FiniteElementBase class for information on its use. All + * class FiniteElement class for information on its use. All * its data are available to the implementation in a concrete finite * element class. * * Remark on a change in implementation: it is now wrong to cast a * pointer to FiniteElement to a pointer to FiniteElementData and * delete it. The virtual destructor has been moved up. In a later - * version, FiniteElementData and FiniteElementBase should be private + * version, FiniteElementData and FiniteElement should be private * base classes of FiniteElement. * * @author Wolfgang Bangerth, Guido Kanschat, 1998, 1999, 2000, 2001, 2003, 2005 @@ -560,1546 +560,6 @@ class FiniteElementData * * @author Wolfgang Bangerth, 1998, 2002; Ralf Hartmann, Guido Kanschat, 2001 */ -template -class FiniteElementBase : public Subscriptor, - public FiniteElementData -{ - public: - /** - * 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. - * - * @author Guido Kanschat, 2001 - */ - class InternalDataBase : public Mapping::InternalDataBase - { - public: - /** - * 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. - */ - void initialize_2nd (const FiniteElement *element, - const Mapping &mapping, - const Quadrature &quadrature); - - /** - * Storage for @p FEValues - * objects needed to - * approximate second - * derivatives. - * - * The ordering is p+hx, - * p+hy, p+hz, - * @p p-hx, @p p-hy, - * @p p-hz, where unused - * entries in lower dimensions - * are missing. - */ - std::vector*> differences; - }; - - /** - * Construct an object of this - * type. You have to set some - * member variables, for example - * some matrices, explicitly - * after calling this base class' - * constructor. For this see the - * existing finite element - * classes. For the second and - * third parameter of this - * constructor, see the documentation - * of the respective member - * variables. - * - * @note Both vector parameters - * should have length - * dofs_per_cell. Nevertheless, - * it is allowed to use vectors - * of length one. In this case, - * the vector is resized to the - * correct length and filled with - * the entry value. - */ - FiniteElementBase (const FiniteElementData &fe_data, - const std::vector &restriction_is_additive_flags, - const std::vector > &nonzero_components); - - /** - * Return a string that uniquely - * identifies a finite - * element. The general - * convention is that this is the - * class name, followed by the - * space 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. - */ - virtual std::string get_name () const = 0; - - /** - * @name Shape function access - * @{ - */ - - /** - * 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 - * @p ExcShapeFunctionNotPrimitive. In - * that case, use the - * shape_value_component() - * function. - * - * An - * @p ExcUnitShapeValuesDoNotExist - * is thrown if the shape values - * of the @p 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 @p 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 - * @p ExcShapeFunctionNotPrimitive. In - * that case, use the - * shape_grad_component() - * function. - * - * An - * @p ExcUnitShapeValuesDoNotExist - * is thrown if the shape values - * of the @p 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 @p 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 - * @p ExcShapeFunctionNotPrimitive. In - * that case, use the - * shape_grad_grad_component() - * function. - * - * An - * @p ExcUnitShapeValuesDoNotExist - * is thrown if the shape values - * of the @p 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; - - /** - * Just like for @p shape_grad_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<2,dim> shape_grad_grad_component (const unsigned int i, - const Point &p, - const unsigned int component) const; - //@} - /** - * @name Transfer and constraint matrices - * @{ - */ - - /** - * 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 - * @p restriction_is_additive_flags - * given to the constructor. 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 - * @p ExcProjectionVoid. - */ - const FullMatrix & - get_restriction_matrix (const unsigned int child) 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 - * @p ExcEmbeddingVoid. You can - * check whether this is the case - * by calling the - * prolongation_is_implemented(). - */ - const FullMatrix & - get_prolongation_matrix (const unsigned int child) const; - - /** - * Return whether this element implements - * its prolongation matrices. 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 prolongation_is_implemented () const; - - /** - * Return whether this element implements - * its restriction matrices. 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 restriction_is_implemented () const; - - /** - * Access the - * @p restriction_is_additive_flag - * field. See there for more - * information on its contents. - * - * 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 readonly 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 space 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 - * @p constraints_are_implemented - * function to check up front - * whethehr this function will - * succeed or generate the - * exception. - */ - const FullMatrix & constraints () const; - - /** - * Return whether this element - * implements its hanging node - * constraints. The return value - * also indicates whether a call - * to the @p constraint 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 - * just expresses. - */ - bool constraints_are_implemented () 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 - * 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 - * FiniteElementBase::ExcInterpolationNotImplemented. - */ - virtual void - get_interpolation_matrix (const FiniteElementBase &source, - FullMatrix &matrix) const; - //@} - - /** - * 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 - * @p restriction and - * @p prolongation. - */ - bool operator == (const FiniteElementBase &) const; - - /** - * @name Index computations - * @{ - */ - /** - * 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 - * @p 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 - * @p system_to_base_table field - * and the - * @p system_to_base_index - * function. - */ - std::pair - system_to_component_index (const unsigned int index) const; - - /** - * 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. - * - * This is the opposite operation - * from the @p system_to_component_index - * function. - */ - unsigned int component_to_system_index(const unsigned int component, - const unsigned int index) const; - - /** - * Same as above, but do it for - * shape functions and their - * indices on a face. - */ - std::pair - face_system_to_component_index (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. - * - * 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 - * @p 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 - * @p system_to_component_index - * function. - */ - std::pair, unsigned int> - system_to_base_index (const unsigned int index) const; - - /** - * Same as - * @p system_to_base_index, but - * for degrees of freedom located - * on a face. - */ - std::pair, unsigned int> - face_system_to_base_index (const unsigned int index) const; - - /** - * Given a vector component, - * return an index 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::@p base_element - * function. - * - * If this is a scalar finite - * element, then the return value - * is always equal to a pair of - * zeros. - */ - std::pair - component_to_base (const unsigned int component) const; - //@} - - /** - * @name Support points and interpolation - * @{ - */ - - /** - * 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. - */ - 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 - * @p 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. - */ - 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. - * - * 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, of not all the - * base elements define support - * points. In this way, you can - * still ask for certain support - * points, even if - * @p 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 - * (@p 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 - * @p get_unit_support_points - * yields a non-empty array. - * - * 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. - */ - virtual - Point - unit_face_support_point (const unsigned int index) const; - - /** - * Return a support point vector - * for generalized interpolation. - */ - const std::vector > & - get_generalized_support_points () const; - - /** - * - */ - bool has_generalized_support_points () const; - - /** - * - */ - const std::vector > & - get_generalized_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 - * @p get_unit_support_points - * yields a non-empty array. - * - * 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. - * - * @note This function is - * implemented in - * FiniteElementBase 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. - * - * 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. - */ - virtual void interpolate( - std::vector& local_dofs, - const VectorSlice > >& values) const; - - //@} - - /** - * Return in which of the vector - * components of this finite - * element the @p ithe 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(i) - * 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. - */ - const std::vector & - 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 - * @p 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. - */ - 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 - * 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. - */ - bool - is_primitive (const unsigned int i) const; - - /** - * Return whether the entire - * finite element is primitive, - * in the sense that all its - * shape functions are - * primitive. If the finite - * element is scalar, then this - * is always the case. - * - * Since this is an extremely - * common operation, the result - * is cached in the - * @p cached_primitivity - * variable which is computed in - * the constructor. - */ - bool - is_primitive () const; - - /** - * Determine an estimate for the - * memory consumption (in bytes) - * of this object. - * - * This function is not - * virtual. Use a - * FiniteElement object to - * get the actual size of a - * concrete element. - */ - unsigned int memory_consumption () const; - - /** - * Exception - * - * @ingroup Exceptions - */ - DeclException1 (ExcShapeFunctionNotPrimitive, - int, - << "The shape function with index " << arg1 - << " is not primitive, i.e. it is vector-valued and " - << "has more than one non-zero vector component. This " - << "function cannot be called for these shape functions. " - << "Maybe you want to use the same function with the " - << "_component suffix?"); - /** - * Exception - * - * @ingroup Exceptions - */ - DeclException0 (ExcFENotPrimitive); - /** - * Exception - * - * @ingroup Exceptions - */ - DeclException0 (ExcUnitShapeValuesDoNotExist); - - /** - * 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. - * - * @ingroup Exceptions - */ - DeclException0 (ExcEmbeddingVoid); - - /** - * Attempt to access restriction - * matrices of a finite element - * which did not implement these - * matrices. - * - * Exception - * @ingroup Exceptions - */ - DeclException0 (ExcProjectionVoid); - - /** - * Attempt to access constraint - * matrices of a finite element - * which did not implement these - * matrices. - * - * Exception - * @ingroup Exceptions - */ - DeclException0 (ExcConstraintsVoid); - - /** - * Exception - * @ingroup Exceptions - */ - DeclException2 (ExcWrongInterfaceMatrixSize, - int, int, - << "The interface matrix has a size of " << arg1 - << "x" << arg2 - << ", which is not reasonable in the present dimension."); - /** - * Exception - * @ingroup Exceptions - */ - DeclException2 (ExcComponentIndexInvalid, - int, int, - << "The component-index pair (" << arg1 << ", " << arg2 - << ") is invalid, i.e. non-existent"); - /** - * Exception - * @ingroup Exceptions - */ - DeclException0 (ExcInterpolationNotImplemented); - - protected: - /** - * Array of projection matrices. See - * get_restriction_matrix() above. - * - * Matrices in this array are - * automatically initialized to - * correct size. If the derived - * finite element class does not - * implement these matrices, they - * should be resized to zero - * size. - */ - FullMatrix restriction[GeometryInfo::children_per_cell]; - - /** - * Array of embedding matrices. See - * get_prolongation_matrix() above. - * - * Matrices in this array are - * automatically initialized to - * correct size. If the derived - * finite element class does not - * implement these matrices, they - * should be resized to zero - * size. - */ - FullMatrix prolongation[GeometryInfo::children_per_cell]; - - /** - * 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. - * - * This field is obviously - * useless in one space dimension - * and has there a zero size. - */ - FullMatrix interface_constraints; - - /** - * 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; - - /** - * Store what - * @p 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. - * - * 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 - * @p 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 - * @p system_to_component_table. - */ - std::vector,unsigned int> > - system_to_base_table; - - /** - * Likewise for the indices on - * faces. - */ - std::vector,unsigned int> > - face_system_to_base_table; - - /** - * The base element establishing - * a component. - * - * This table converts a - * component number to a pair - * consisting of the - * @p base_element number, and - * the component within this base - * element. While component - * information contains - * multiplicity of base elements, - * the result allows access to - * shape functions of the base - * element. - * - * 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 > 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 - * @p 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; - - /** - * 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 - * @p 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. - */ - std::vector > generalized_support_points; - - /** - * Face support points used for - * interpolation functions of - * non-Lagrangian elements. - */ - std::vector > generalized_face_support_points; - - /** - * 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. - */ - const std::vector > nonzero_components; - - /** - * This array holds how many - * values in the respective entry - * of the @p 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; - - /** - * Store whether all shape - * functions are primitive. Since - * finding this out is a very - * common operation, we cache the - * result, i.e. compute the value - * in the constructor for simpler - * access. - */ - const bool cached_primitivity; - - /** - * Compute second derivatives by - * finite differences of - * gradients. - */ - void compute_2nd (const Mapping &mapping, - const typename Triangulation::cell_iterator &cell, - const unsigned int offset, - typename Mapping::InternalDataBase &mapping_internal, - InternalDataBase &fe_internal, - FEValuesData &data) const; - - private: - /** - * 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; - - /** - * 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); - - /** - * Allow the FESystem class to access the - * restriction and prolongation matrices - * directly. Hence, FESystem has the - * possibility to see if these matrices - * are initialized without accessing - * these matrices through the - * @p get_restriction_matrix and - * @p get_prolongation_matrix - * functions. This is important as these - * functions include assertions that - * throw if the matrices are not already - * initialized. - */ - template friend class FESystem; - - /** - * Make the inner class a - * friend. This is not strictly - * necessary, but the Intel - * compiler seems to want this. - */ - friend class InternalDataBase; -}; /*@}*/ @@ -2186,208 +646,6 @@ FiniteElementData::conforms (Conformity space) const return ((space & conforming_space) != 0); } -//----------------------------------------------------------------------// - -template -inline -std::pair -FiniteElementBase::system_to_component_index (const unsigned int index) const -{ - Assert (index < system_to_component_table.size(), - ExcIndexRange(index, 0, system_to_component_table.size())); - Assert (is_primitive (index), - typename FiniteElementBase::ExcShapeFunctionNotPrimitive(index)); - return system_to_component_table[index]; -} - -template -inline -unsigned int -FiniteElementBase::component_to_system_index (const unsigned int component, - const unsigned int index) const -{ - std::vector< std::pair >::const_iterator - it = std::find(system_to_component_table.begin(), system_to_component_table.end(), - std::pair(component, index)); - - Assert(it != system_to_component_table.end(), ExcComponentIndexInvalid(component, index)); - return std::distance(system_to_component_table.begin(), it); -} - - - -template -inline -std::pair -FiniteElementBase::face_system_to_component_index (const unsigned int index) const -{ - Assert(index < face_system_to_component_table.size(), - ExcIndexRange(index, 0, face_system_to_component_table.size())); - - // in debug mode, check whether the - // function is primitive, since - // otherwise the result may have no - // meaning - // - // since the primitivity tables are - // all geared towards cell dof - // indices, rather than face dof - // indices, we have to work a - // little bit... - // - // in 1d, the face index is equal - // to the cell index - Assert (((dim == 1) && is_primitive(index)) - || - // in 2d, construct it like - // this: - ((dim == 2) && - is_primitive (index < (GeometryInfo<2>::vertices_per_face * - this->dofs_per_vertex) - ? - index - : - GeometryInfo<2>::vertices_per_cell * - this->dofs_per_vertex + - (index - - GeometryInfo<2>::vertices_per_face * - this->dofs_per_vertex))) - || - // likewise in 3d, but more - // complicated - ((dim == 3) && - is_primitive (index < (GeometryInfo<3>::vertices_per_face * - this->dofs_per_vertex) - ? - index - : - (index < (GeometryInfo<3>::vertices_per_face * - this->dofs_per_vertex - + - GeometryInfo<3>::lines_per_face * - this->dofs_per_line) - ? - GeometryInfo<3>::vertices_per_cell * - this->dofs_per_vertex + - (index - - GeometryInfo<3>::vertices_per_face * - this->dofs_per_vertex) - : - GeometryInfo<3>::vertices_per_cell * - this->dofs_per_vertex + - GeometryInfo<3>::lines_per_cell * - this->dofs_per_line + - (index - - GeometryInfo<3>::vertices_per_face * - this->dofs_per_vertex - - - GeometryInfo<3>::lines_per_face * - this->dofs_per_line)))), - typename FiniteElementBase::ExcShapeFunctionNotPrimitive(index)); - - return face_system_to_component_table[index]; -} - - - -template -inline -std::pair,unsigned int> -FiniteElementBase::system_to_base_index (const unsigned int index) const -{ - Assert (index < system_to_base_table.size(), - ExcIndexRange(index, 0, system_to_base_table.size())); - return system_to_base_table[index]; -} - - - - -template -inline -std::pair,unsigned int> -FiniteElementBase::face_system_to_base_index (const unsigned int index) const -{ - Assert(index < face_system_to_base_table.size(), - ExcIndexRange(index, 0, face_system_to_base_table.size())); - return face_system_to_base_table[index]; -} - - - -template -inline -std::pair -FiniteElementBase::component_to_base (const unsigned int index) const -{ - Assert(index < component_to_base_table.size(), - ExcIndexRange(index, 0, component_to_base_table.size())); - - return component_to_base_table[index]; -} - - -template -inline -bool -FiniteElementBase::restriction_is_additive (const unsigned int index) const -{ - Assert(index < this->dofs_per_cell, - ExcIndexRange(index, 0, this->dofs_per_cell)); - return restriction_is_additive_flags[index]; -} - - -template -inline -const std::vector & -FiniteElementBase::get_nonzero_components (const unsigned int i) const -{ - Assert (i < this->dofs_per_cell, ExcIndexRange (i, 0, this->dofs_per_cell)); - return nonzero_components[i]; -} - - - -template -inline -unsigned int -FiniteElementBase::n_nonzero_components (const unsigned int i) const -{ - Assert (i < this->dofs_per_cell, ExcIndexRange (i, 0, this->dofs_per_cell)); - return n_nonzero_components_table[i]; -} - - - -template -inline -bool -FiniteElementBase::is_primitive (const unsigned int i) const -{ - Assert (i < this->dofs_per_cell, ExcIndexRange (i, 0, this->dofs_per_cell)); - - // return primitivity of a shape - // function by checking whether it - // has more than one non-zero - // component or not. we could cache - // this value in an array of bools, - // but accessing a bit-vector (as - // std::vector is) is - // probably more expensive than - // just comparing against 1 - return (n_nonzero_components_table[i] == 1); -} - - -template -inline -bool -FiniteElementBase::is_primitive () const -{ - return cached_primitivity; -} - diff --git a/deal.II/deal.II/include/fe/fe_dgp_monomial.h b/deal.II/deal.II/include/fe/fe_dgp_monomial.h index 2409122a6e..5698ce39fb 100644 --- a/deal.II/deal.II/include/fe/fe_dgp_monomial.h +++ b/deal.II/deal.II/include/fe/fe_dgp_monomial.h @@ -68,11 +68,11 @@ class FE_DGPMonomial : public FE_Poly,dim> * element is also a @p FE_Q * element. Otherwise, an * exception of type - * FiniteElementBase::ExcInterpolationNotImplemented + * FiniteElement::ExcInterpolationNotImplemented * is thrown. */ virtual void - get_interpolation_matrix (const FiniteElementBase &source, + get_interpolation_matrix (const FiniteElement &source, FullMatrix &matrix) const; /** diff --git a/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h b/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h index 6f6274308a..0b7fab8400 100644 --- a/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h +++ b/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h @@ -70,7 +70,7 @@ class FE_DGPNonparametric : public FiniteElement * Return the value of the * @p ith shape function at the * point @p p. See the - * FiniteElementBase base + * FiniteElement base * class for more information * about the semantics of this * function. @@ -84,7 +84,7 @@ class FE_DGPNonparametric : public FiniteElement * component of the @p ith shape * function at the point * @p p. See the - * FiniteElementBase base + * FiniteElement base * class for more information * about the semantics of this * function. @@ -104,7 +104,7 @@ class FE_DGPNonparametric : public FiniteElement * Return the gradient of the * @p ith shape function at the * point @p p. See the - * FiniteElementBase base + * FiniteElement base * class for more information * about the semantics of this * function. @@ -118,7 +118,7 @@ class FE_DGPNonparametric : public FiniteElement * component of the @p ith shape * function at the point * @p p. See the - * FiniteElementBase base + * FiniteElement base * class for more information * about the semantics of this * function. @@ -139,7 +139,7 @@ class FE_DGPNonparametric : public FiniteElement * derivatives of the @p ith * shape function at point @p p * on the unit cell. See the - * FiniteElementBase base + * FiniteElement base * class for more information * about the semantics of this * function. @@ -153,7 +153,7 @@ class FE_DGPNonparametric : public FiniteElement * component of the @p ith shape * function at the point * @p p. See the - * FiniteElementBase base + * FiniteElement base * class for more information * about the semantics of this * function. @@ -418,7 +418,7 @@ class FE_DGPNonparametric : public FiniteElement * see the documentation of the * base class. */ - class InternalData : public FiniteElementBase::InternalDataBase + class InternalData : public FiniteElement::InternalDataBase { public: // have some scratch arrays diff --git a/deal.II/deal.II/include/fe/fe_dgq.h b/deal.II/deal.II/include/fe/fe_dgq.h index 7980415915..27fff613f0 100644 --- a/deal.II/deal.II/include/fe/fe_dgq.h +++ b/deal.II/deal.II/include/fe/fe_dgq.h @@ -115,11 +115,11 @@ class FE_DGQ : public FE_Poly,dim> * element is also a @p FE_DGQ * element. Otherwise, an * exception of type - * FiniteElementBase::ExcInterpolationNotImplemented + * FiniteElement::ExcInterpolationNotImplemented * is thrown. */ virtual void - get_interpolation_matrix (const FiniteElementBase &source, + get_interpolation_matrix (const FiniteElement &source, FullMatrix &matrix) const; /** diff --git a/deal.II/deal.II/include/fe/fe_nedelec.h b/deal.II/deal.II/include/fe/fe_nedelec.h index f5b06b2935..e4f8b12f20 100644 --- a/deal.II/deal.II/include/fe/fe_nedelec.h +++ b/deal.II/deal.II/include/fe/fe_nedelec.h @@ -205,7 +205,7 @@ class FE_Nedelec : public FiniteElement * component of the @p ith shape * function at the point * @p p. See the - * FiniteElementBase base + * FiniteElement base * class for more information * about the semantics of this * function. @@ -220,7 +220,7 @@ class FE_Nedelec : public FiniteElement * component of the @p ith shape * function at the point * @p p. See the - * FiniteElementBase base + * FiniteElement base * class for more information * about the semantics of this * function. @@ -235,7 +235,7 @@ class FE_Nedelec : public FiniteElement * component of the @p ith shape * function at the point * @p p. See the - * FiniteElementBase base + * FiniteElement base * class for more information * about the semantics of this * function. @@ -482,7 +482,7 @@ class FE_Nedelec : public FiniteElement /** * Initialize the * @p unit_support_points field - * of the FiniteElementBase + * of the FiniteElement * class. Called from the * constructor. */ @@ -491,7 +491,7 @@ class FE_Nedelec : public FiniteElement /** * Initialize the * @p unit_face_support_points field - * of the FiniteElementBase + * of the FiniteElement * class. Called from the * constructor. */ @@ -549,7 +549,7 @@ class FE_Nedelec : public FiniteElement * see the documentation of the * base class. */ - class InternalData : public FiniteElementBase::InternalDataBase + class InternalData : public FiniteElement::InternalDataBase { public: /** diff --git a/deal.II/deal.II/include/fe/fe_poly.h b/deal.II/deal.II/include/fe/fe_poly.h index 4b15ec7b2f..bffcdc241d 100644 --- a/deal.II/deal.II/include/fe/fe_poly.h +++ b/deal.II/deal.II/include/fe/fe_poly.h @@ -42,7 +42,7 @@ * * This class is not a fully implemented FiniteElement class. Instead * there are several pure virtual functions declared in the - * FiniteElement and FiniteElementBase classes which cannot + * FiniteElement and FiniteElement classes which cannot * implemented by this class but are left for implementation in * derived classes. * @@ -83,7 +83,7 @@ class FE_Poly : public FiniteElement * Return the value of the * ith shape function at * the point p. See the - * FiniteElementBase base class + * FiniteElement base class * for more information about the * semantics of this function. */ @@ -96,7 +96,7 @@ class FE_Poly : public FiniteElement * component of the ith * shape function at the point * p. See the - * FiniteElementBase base class + * FiniteElement base class * for more information about the * semantics of this function. * @@ -115,7 +115,7 @@ class FE_Poly : public FiniteElement * Return the gradient of the * ith shape function at * the point p. See the - * FiniteElementBase base class + * FiniteElement base class * for more information about the * semantics of this function. */ @@ -128,7 +128,7 @@ class FE_Poly : public FiniteElement * component of the ith * shape function at the point * p. See the - * FiniteElementBase base class + * FiniteElement base class * for more information about the * semantics of this function. * @@ -149,7 +149,7 @@ class FE_Poly : public FiniteElement * ith shape function at * point p on the unit * cell. See the - * FiniteElementBase base class + * FiniteElement base class * for more information about the * semantics of this function. */ @@ -162,7 +162,7 @@ class FE_Poly : public FiniteElement * vector component of the * ith shape function at * the point p. See the - * FiniteElementBase base class + * FiniteElement base class * for more information about the * semantics of this function. * @@ -336,7 +336,7 @@ class FE_Poly : public FiniteElement * see the documentation of the * base class. */ - class InternalData : public FiniteElementBase::InternalDataBase + class InternalData : public FiniteElement::InternalDataBase { public: /** diff --git a/deal.II/deal.II/include/fe/fe_poly_tensor.h b/deal.II/deal.II/include/fe/fe_poly_tensor.h index 7378a58d63..c0cf51b345 100644 --- a/deal.II/deal.II/include/fe/fe_poly_tensor.h +++ b/deal.II/deal.II/include/fe/fe_poly_tensor.h @@ -318,7 +318,7 @@ class FE_PolyTensor : public FiniteElement * accessed by indices * (i,k). */ - class InternalData : public FiniteElementBase::InternalDataBase + class InternalData : public FiniteElement::InternalDataBase { public: /** diff --git a/deal.II/deal.II/include/fe/fe_q.h b/deal.II/deal.II/include/fe/fe_q.h index 8ca83af781..21588ebe37 100644 --- a/deal.II/deal.II/include/fe/fe_q.h +++ b/deal.II/deal.II/include/fe/fe_q.h @@ -269,11 +269,11 @@ class FE_Q : public FE_Poly,dim> * element is also a @p FE_Q * element. Otherwise, an * exception of type - * FiniteElementBase::ExcInterpolationNotImplemented + * FiniteElement::ExcInterpolationNotImplemented * is thrown. */ virtual void - get_interpolation_matrix (const FiniteElementBase &source, + get_interpolation_matrix (const FiniteElement &source, FullMatrix &matrix) const; /** @@ -364,7 +364,7 @@ class FE_Q : public FE_Poly,dim> /** * Initialize the * @p unit_support_points field - * of the FiniteElementBase + * of the FiniteElement * class. Called from the * constructor. */ @@ -373,7 +373,7 @@ class FE_Q : public FE_Poly,dim> /** * Initialize the * @p unit_face_support_points field - * of the FiniteElementBase + * of the FiniteElement * class. Called from the * constructor. */ diff --git a/deal.II/deal.II/include/fe/fe_q_hierarchical.h b/deal.II/deal.II/include/fe/fe_q_hierarchical.h index 3ba126f3ad..1c8aa5e768 100644 --- a/deal.II/deal.II/include/fe/fe_q_hierarchical.h +++ b/deal.II/deal.II/include/fe/fe_q_hierarchical.h @@ -409,7 +409,7 @@ class FE_Q_Hierarchical : public FE_Poly,dim> /** * Initialize the * @p unit_support_points field - * of the FiniteElementBase + * of the FiniteElement * class. Called from the * constructor. */ @@ -418,7 +418,7 @@ class FE_Q_Hierarchical : public FE_Poly,dim> /** * Initialize the * @p unit_face_support_points field - * of the FiniteElementBase + * of the FiniteElement * class. Called from the * constructor. */ diff --git a/deal.II/deal.II/include/fe/fe_raviart_thomas.h b/deal.II/deal.II/include/fe/fe_raviart_thomas.h index 0cea7b8eb6..2f7f3f07a0 100644 --- a/deal.II/deal.II/include/fe/fe_raviart_thomas.h +++ b/deal.II/deal.II/include/fe/fe_raviart_thomas.h @@ -163,7 +163,7 @@ class FE_RaviartThomas : public FiniteElement * component of the @p ith shape * function at the point * @p p. See the - * FiniteElementBase base + * FiniteElement base * class for more information * about the semantics of this * function. @@ -178,7 +178,7 @@ class FE_RaviartThomas : public FiniteElement * component of the @p ith shape * function at the point * @p p. See the - * FiniteElementBase base + * FiniteElement base * class for more information * about the semantics of this * function. @@ -193,7 +193,7 @@ class FE_RaviartThomas : public FiniteElement * component of the @p ith shape * function at the point * @p p. See the - * FiniteElementBase base + * FiniteElement base * class for more information * about the semantics of this * function. @@ -237,11 +237,11 @@ class FE_RaviartThomas : public FiniteElement * element is also a Raviart * Thomas element. Otherwise, an * exception of type - * FiniteElementBase::ExcInterpolationNotImplemented + * FiniteElement::ExcInterpolationNotImplemented * is thrown. */ virtual void - get_interpolation_matrix (const FiniteElementBase &source, + get_interpolation_matrix (const FiniteElement &source, FullMatrix &matrix) const; /** @@ -439,7 +439,7 @@ class FE_RaviartThomas : public FiniteElement /** * Initialize the * @p unit_support_points field - * of the FiniteElementBase + * of the FiniteElement * class. Called from the * constructor. */ @@ -448,7 +448,7 @@ class FE_RaviartThomas : public FiniteElement /** * Initialize the * @p unit_face_support_points field - * of the FiniteElementBase + * of the FiniteElement * class. Called from the * constructor. */ @@ -501,7 +501,7 @@ class FE_RaviartThomas : public FiniteElement * see the documentation of the * base class. */ - class InternalData : public FiniteElementBase::InternalDataBase + class InternalData : public FiniteElement::InternalDataBase { public: /** @@ -655,8 +655,8 @@ class FE_RaviartThomasNodal get_ria_vector (const unsigned int degree); /** * Initialize the - * FiniteElementBase::unit_support_points - * and FiniteElementBase::unit_face_support_points + * FiniteElement::unit_support_points + * and FiniteElement::unit_face_support_points * fields. Called from the * constructor. */ @@ -739,7 +739,7 @@ FE_RaviartThomas<3>::initialize_restriction (); template <> void FE_RaviartThomas<1>:: -get_interpolation_matrix (const FiniteElementBase<1> &, +get_interpolation_matrix (const FiniteElement<1> &, FullMatrix &) const; /// @endif diff --git a/deal.II/deal.II/include/fe/fe_system.h b/deal.II/deal.II/include/fe/fe_system.h index 1928683a2a..64f34b0539 100644 --- a/deal.II/deal.II/include/fe/fe_system.h +++ b/deal.II/deal.II/include/fe/fe_system.h @@ -174,7 +174,7 @@ class FESystem : public FiniteElement * component of the @p ith shape * function at the point * @p p. See the - * FiniteElementBase base + * FiniteElement base * class for more information * about the semantics of this * function. @@ -227,7 +227,7 @@ class FESystem : public FiniteElement * component of the @p ith shape * function at the point * @p p. See the - * FiniteElementBase base + * FiniteElement base * class for more information * about the semantics of this * function. @@ -279,7 +279,7 @@ class FESystem : public FiniteElement * component of the @p ith shape * function at the point * @p p. See the - * FiniteElementBase base + * FiniteElement base * class for more information * about the semantics of this * function. @@ -314,11 +314,11 @@ class FESystem : public FiniteElement * @p get_interpolation_matrix * functions. Otherwise, an * exception of type - * FiniteElementBase::ExcInterpolationNotImplemented + * FiniteElement::ExcInterpolationNotImplemented * is thrown. */ virtual void - get_interpolation_matrix (const FiniteElementBase &source, + get_interpolation_matrix (const FiniteElement &source, FullMatrix &matrix) const; /** @@ -558,7 +558,7 @@ class FESystem : public FiniteElement /** * Initialize the * @p unit_support_points field - * of the FiniteElementBase + * of the FiniteElement * class. Called from the * constructor. */ @@ -567,7 +567,7 @@ class FESystem : public FiniteElement /** * Initialize the * @p unit_face_support_points field - * of the FiniteElementBase + * of the FiniteElement * class. Called from the * constructor. */ @@ -743,7 +743,7 @@ class FESystem : public FiniteElement * @p InternalData objects for * each of the base elements. */ - class InternalData : public FiniteElementBase::InternalDataBase + class InternalData : public FiniteElement::InternalDataBase { public: /** @@ -780,7 +780,7 @@ class FESystem : public FiniteElement * element. */ void set_fe_data(const unsigned int base_no, - typename FiniteElementBase::InternalDataBase *); + typename FiniteElement::InternalDataBase *); /** * Gives read-access to the @@ -788,7 +788,7 @@ class FESystem : public FiniteElement * @p InternalData of the * @p base_noth base element. */ - typename FiniteElementBase::InternalDataBase & + typename FiniteElement::InternalDataBase & get_fe_data (const unsigned int base_no) const; @@ -872,7 +872,7 @@ class FESystem : public FiniteElement * elements, irrespective of * their multiplicity. */ - typename std::vector::InternalDataBase *> base_fe_datas; + typename std::vector::InternalDataBase *> base_fe_datas; /** * Pointers to the diff --git a/deal.II/deal.II/include/fe/fe_tools.h b/deal.II/deal.II/include/fe/fe_tools.h index 5bdb67003e..5c88bc49fc 100644 --- a/deal.II/deal.II/include/fe/fe_tools.h +++ b/deal.II/deal.II/include/fe/fe_tools.h @@ -195,7 +195,7 @@ class FETools * @param matrices A pointer to * 2dim FullMatrix * objects. This is the format - * used in FiniteElementBase, + * used in FiniteElement, * where we want to use ths * function mostly. */ @@ -241,7 +241,7 @@ class FETools * matrices. @arg matrices A pointer to * GeometryInfo::children_per_cell=2dim * FullMatrix objects. This is the format - * used in FiniteElementBase, where we + * used in FiniteElement, where we * want to use this function mostly. */ template diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 3e8a53a2f6..bc82fc01c4 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -92,7 +92,7 @@ DoFTools::compute_row_length_vector( const Table<2,Coupling>& flux_couplings) { const FiniteElement& fe = dofs.get_fe(); - Assert (fe.is_primitive(), typename FiniteElementBase::ExcFENotPrimitive()); + Assert (fe.is_primitive(), typename FiniteElement::ExcFENotPrimitive()); Assert (row_lengths.size() == dofs.n_dofs(), ExcDimensionMismatch(row_lengths.size(), dofs.n_dofs())); Assert (couplings.n_rows()==fe.n_components(), @@ -329,7 +329,7 @@ DoFTools::compute_row_length_vector( const Table<2,Coupling>& flux_couplings) { const FiniteElement& fe = dofs.get_fe(); - Assert (fe.is_primitive(), typename FiniteElementBase::ExcFENotPrimitive()); + Assert (fe.is_primitive(), typename FiniteElement::ExcFENotPrimitive()); Assert (row_lengths.size() == dofs.n_dofs(), ExcDimensionMismatch(row_lengths.size(), dofs.n_dofs())); Assert (couplings.n_rows()==fe.n_components(), diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 1463402d32..6747e01832 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -29,16 +29,16 @@ -/*------------------------------- FiniteElementBase ----------------------*/ +/*------------------------------- FiniteElement ----------------------*/ template -const double FiniteElementBase::fd_step_length = 1.0e-6; +const double FiniteElement::fd_step_length = 1.0e-6; template void -FiniteElementBase:: +FiniteElement:: InternalDataBase::initialize_2nd (const FiniteElement *element, const Mapping &mapping, const Quadrature &quadrature) @@ -89,7 +89,7 @@ InternalDataBase::initialize_2nd (const FiniteElement *element, template -FiniteElementBase::InternalDataBase::~InternalDataBase () +FiniteElement::InternalDataBase::~InternalDataBase () { for (unsigned int i=0; i::InternalDataBase::~InternalDataBase () template -FiniteElementBase::FiniteElementBase ( +FiniteElement::FiniteElement ( const FiniteElementData &fe_data, const std::vector &r_i_a_f, const std::vector > &nonzero_c) @@ -199,9 +199,31 @@ FiniteElementBase::FiniteElementBase ( } +template +FiniteElement::FiniteElement (const FiniteElement &) + : + Subscriptor(), + FiniteElementData() +{ + Assert (false, + ExcMessage ("Finite element objects don't support copying " + "semantics through the copy constructor. If " + "you want to copy a finite element, use the " + "clone() function.")); +} + + + +template +FiniteElement::~FiniteElement () +{} + + + + template double -FiniteElementBase::shape_value (const unsigned int, +FiniteElement::shape_value (const unsigned int, const Point &) const { AssertThrow(false, ExcUnitShapeValuesDoNotExist()); @@ -212,7 +234,7 @@ FiniteElementBase::shape_value (const unsigned int, template double -FiniteElementBase::shape_value_component (const unsigned int, +FiniteElement::shape_value_component (const unsigned int, const Point &, const unsigned int) const { @@ -224,7 +246,7 @@ FiniteElementBase::shape_value_component (const unsigned int, template Tensor<1,dim> -FiniteElementBase::shape_grad (const unsigned int, +FiniteElement::shape_grad (const unsigned int, const Point &) const { AssertThrow(false, ExcUnitShapeValuesDoNotExist()); @@ -235,7 +257,7 @@ FiniteElementBase::shape_grad (const unsigned int, template Tensor<1,dim> -FiniteElementBase::shape_grad_component (const unsigned int, +FiniteElement::shape_grad_component (const unsigned int, const Point &, const unsigned int) const { @@ -247,7 +269,7 @@ FiniteElementBase::shape_grad_component (const unsigned int, template Tensor<2,dim> -FiniteElementBase::shape_grad_grad (const unsigned int, +FiniteElement::shape_grad_grad (const unsigned int, const Point &) const { AssertThrow(false, ExcUnitShapeValuesDoNotExist()); @@ -258,7 +280,7 @@ FiniteElementBase::shape_grad_grad (const unsigned int, template Tensor<2,dim> -FiniteElementBase::shape_grad_grad_component (const unsigned int, +FiniteElement::shape_grad_grad_component (const unsigned int, const Point &, const unsigned int) const { @@ -269,7 +291,7 @@ FiniteElementBase::shape_grad_grad_component (const unsigned int, template const FullMatrix & -FiniteElementBase::get_restriction_matrix (const unsigned int child) const +FiniteElement::get_restriction_matrix (const unsigned int child) const { Assert (child::children_per_cell, ExcIndexRange(child, 0, GeometryInfo::children_per_cell)); @@ -281,7 +303,7 @@ FiniteElementBase::get_restriction_matrix (const unsigned int child) const template const FullMatrix & -FiniteElementBase::get_prolongation_matrix (const unsigned int child) const +FiniteElement::get_prolongation_matrix (const unsigned int child) const { Assert (child::children_per_cell, ExcIndexRange(child, 0, GeometryInfo::children_per_cell)); @@ -292,7 +314,7 @@ FiniteElementBase::get_prolongation_matrix (const unsigned int child) const template bool -FiniteElementBase::prolongation_is_implemented () const +FiniteElement::prolongation_is_implemented () const { for (unsigned int c=0; c::children_per_cell; ++c) { @@ -313,7 +335,7 @@ FiniteElementBase::prolongation_is_implemented () const template bool -FiniteElementBase::restriction_is_implemented () const +FiniteElement::restriction_is_implemented () const { for (unsigned int c=0; c::children_per_cell; ++c) { @@ -334,7 +356,7 @@ FiniteElementBase::restriction_is_implemented () const template bool -FiniteElementBase::constraints_are_implemented () const +FiniteElement::constraints_are_implemented () const { return (this->dofs_per_face == 0) || (interface_constraints.m() != 0); } @@ -343,7 +365,7 @@ FiniteElementBase::constraints_are_implemented () const template const FullMatrix & -FiniteElementBase::constraints () const +FiniteElement::constraints () const { Assert ((this->dofs_per_face == 0) || (interface_constraints.m() != 0), ExcConstraintsVoid()); @@ -360,7 +382,7 @@ FiniteElementBase::constraints () const template TableIndices<2> -FiniteElementBase::interface_constraints_size () const +FiniteElement::interface_constraints_size () const { switch (dim) { @@ -386,15 +408,15 @@ FiniteElementBase::interface_constraints_size () const template void -FiniteElementBase:: -get_interpolation_matrix (const FiniteElementBase &, +FiniteElement:: +get_interpolation_matrix (const FiniteElement &, FullMatrix &) const { // by default, no interpolation // implemented. so throw exception, // as documentation says AssertThrow (false, - typename FiniteElementBase:: + typename FiniteElement:: ExcInterpolationNotImplemented()); } @@ -403,7 +425,7 @@ get_interpolation_matrix (const FiniteElementBase &, template bool -FiniteElementBase::operator == (const FiniteElementBase &f) const +FiniteElement::operator == (const FiniteElement &f) const { return ((static_cast&>(*this) == static_cast&>(f)) && @@ -414,7 +436,7 @@ FiniteElementBase::operator == (const FiniteElementBase &f) const template const std::vector > & -FiniteElementBase::get_unit_support_points () const +FiniteElement::get_unit_support_points () const { // a finite element may define // support points, but only if @@ -430,7 +452,7 @@ FiniteElementBase::get_unit_support_points () const template bool -FiniteElementBase::has_support_points () const +FiniteElement::has_support_points () const { return (unit_support_points.size() != 0); } @@ -439,7 +461,7 @@ FiniteElementBase::has_support_points () const template const std::vector > & -FiniteElementBase::get_generalized_support_points () const +FiniteElement::get_generalized_support_points () const { // a finite element may define // support points, but only if @@ -455,7 +477,7 @@ FiniteElementBase::get_generalized_support_points () const template bool -FiniteElementBase::has_generalized_support_points () const +FiniteElement::has_generalized_support_points () const { return (generalized_support_points.size() != 0); } @@ -464,7 +486,7 @@ FiniteElementBase::has_generalized_support_points () const template Point -FiniteElementBase::unit_support_point (const unsigned index) const +FiniteElement::unit_support_point (const unsigned index) const { Assert (index < this->dofs_per_cell, ExcIndexRange (index, 0, this->dofs_per_cell)); @@ -477,7 +499,7 @@ FiniteElementBase::unit_support_point (const unsigned index) const template const std::vector > & -FiniteElementBase::get_unit_face_support_points () const +FiniteElement::get_unit_face_support_points () const { // a finite element may define // support points, but only if @@ -493,7 +515,7 @@ FiniteElementBase::get_unit_face_support_points () const template bool -FiniteElementBase::has_face_support_points () const +FiniteElement::has_face_support_points () const { return (unit_face_support_points.size() != 0); } @@ -502,7 +524,7 @@ FiniteElementBase::has_face_support_points () const template const std::vector > & -FiniteElementBase::get_generalized_face_support_points () const +FiniteElement::get_generalized_face_support_points () const { // a finite element may define // support points, but only if @@ -518,7 +540,7 @@ FiniteElementBase::get_generalized_face_support_points () const template bool -FiniteElementBase::has_generalized_face_support_points () const +FiniteElement::has_generalized_face_support_points () const { return (generalized_face_support_points.size() != 0); } @@ -527,7 +549,7 @@ FiniteElementBase::has_generalized_face_support_points () const template Point -FiniteElementBase::unit_face_support_point (const unsigned index) const +FiniteElement::unit_face_support_point (const unsigned index) const { Assert (index < this->dofs_per_face, ExcIndexRange (index, 0, this->dofs_per_face)); @@ -540,7 +562,7 @@ FiniteElementBase::unit_face_support_point (const unsigned index) const template void -FiniteElementBase::interpolate( +FiniteElement::interpolate( std::vector& local_dofs, const std::vector& values) const { @@ -562,7 +584,7 @@ FiniteElementBase::interpolate( template void -FiniteElementBase::interpolate( +FiniteElement::interpolate( std::vector& local_dofs, const std::vector >& values, unsigned int offset) const @@ -588,7 +610,7 @@ FiniteElementBase::interpolate( template void -FiniteElementBase::interpolate( +FiniteElement::interpolate( std::vector& local_dofs, const VectorSlice > >& values) const { @@ -613,7 +635,7 @@ FiniteElementBase::interpolate( template unsigned int -FiniteElementBase::memory_consumption () const +FiniteElement::memory_consumption () const { return (sizeof(FiniteElementData) + MemoryConsumption:: @@ -637,7 +659,7 @@ FiniteElementBase::memory_consumption () const template void -FiniteElementBase:: +FiniteElement:: compute_2nd (const Mapping &mapping, const typename Triangulation::cell_iterator &cell, const unsigned int, @@ -780,7 +802,7 @@ compute_2nd (const Mapping &mapping, template std::vector -FiniteElementBase::compute_n_nonzero_components ( +FiniteElement::compute_n_nonzero_components ( const std::vector > &nonzero_components) { std::vector retval (nonzero_components.size()); @@ -795,40 +817,6 @@ FiniteElementBase::compute_n_nonzero_components ( /*------------------------------- FiniteElement ----------------------*/ -template -FiniteElement::FiniteElement (const FiniteElementData &fe_data, - const std::vector &restriction_is_additive_flags, - const std::vector > &nonzero_components) - : - FiniteElementBase (fe_data, - restriction_is_additive_flags, - nonzero_components) -{} - - - -template -FiniteElement::FiniteElement (const FiniteElement &) - : - FiniteElementBase (FiniteElementData(), - std::vector (), - std::vector >()) -{ - Assert (false, - ExcMessage ("Finite element objects don't support copying " - "semantics through the copy constructor. If " - "you want to copy a finite element, use the " - "clone() function.")); -} - - - -template -FiniteElement::~FiniteElement () -{} - - - template typename Mapping::InternalDataBase * FiniteElement::get_face_data (const UpdateFlags flags, @@ -853,16 +841,6 @@ FiniteElement::get_subface_data (const UpdateFlags flags, -template -unsigned int -FiniteElement::memory_consumption () const -{ - return FiniteElementBase::memory_consumption (); -} - - - - template const FiniteElement& FiniteElement::base_element(const unsigned index) const @@ -873,5 +851,4 @@ FiniteElement::base_element(const unsigned index) const /*------------------------------- Explicit Instantiations -------------*/ -template class FiniteElementBase; template class FiniteElement; diff --git a/deal.II/deal.II/source/fe/fe_dgp_monomial.cc b/deal.II/deal.II/source/fe/fe_dgp_monomial.cc index 32f0ea11c1..36746d0a8e 100644 --- a/deal.II/deal.II/source/fe/fe_dgp_monomial.cc +++ b/deal.II/deal.II/source/fe/fe_dgp_monomial.cc @@ -179,7 +179,7 @@ FE_DGPMonomial::clone() const template void FE_DGPMonomial:: -get_interpolation_matrix (const FiniteElementBase &source_fe, +get_interpolation_matrix (const FiniteElement &source_fe, FullMatrix &interpolation_matrix) const { const FE_DGPMonomial *source_dgp_monomial diff --git a/deal.II/deal.II/source/fe/fe_dgq.cc b/deal.II/deal.II/source/fe/fe_dgq.cc index 81dfb6658a..02bf00a1c2 100644 --- a/deal.II/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/deal.II/source/fe/fe_dgq.cc @@ -361,7 +361,7 @@ FE_DGQ::rotate_indices (std::vector &numbers, template void FE_DGQ:: -get_interpolation_matrix (const FiniteElementBase &x_source_fe, +get_interpolation_matrix (const FiniteElement &x_source_fe, FullMatrix &interpolation_matrix) const { // this is only implemented, if the @@ -370,7 +370,7 @@ get_interpolation_matrix (const FiniteElementBase &x_source_fe, AssertThrow ((x_source_fe.get_name().find ("FE_DGQ<") == 0) || (dynamic_cast*>(&x_source_fe) != 0), - typename FiniteElementBase:: + typename FiniteElement:: ExcInterpolationNotImplemented()); // ok, source is a Q element, so diff --git a/deal.II/deal.II/source/fe/fe_poly_tensor.cc b/deal.II/deal.II/source/fe/fe_poly_tensor.cc index 1849c87baa..7c92ab7d75 100644 --- a/deal.II/deal.II/source/fe/fe_poly_tensor.cc +++ b/deal.II/deal.II/source/fe/fe_poly_tensor.cc @@ -40,7 +40,7 @@ double FE_PolyTensor::shape_value ( const unsigned int, const Point &) const { - Assert(false, typename FiniteElementBase::ExcFENotPrimitive()); + Assert(false, typename FiniteElement::ExcFENotPrimitive()); return 0.; } @@ -78,7 +78,7 @@ Tensor<1,dim> FE_PolyTensor::shape_grad ( const unsigned int, const Point &) const { - Assert(false, typename FiniteElementBase::ExcFENotPrimitive()); + Assert(false, typename FiniteElement::ExcFENotPrimitive()); return Tensor<1,dim>(); } @@ -118,7 +118,7 @@ Tensor<2,dim> FE_PolyTensor::shape_grad_grad ( const unsigned int, const Point &) const { - Assert(false, typename FiniteElementBase::ExcFENotPrimitive()); + Assert(false, typename FiniteElement::ExcFENotPrimitive()); return Tensor<2,dim>(); } diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index d534e4e83f..bdddc26a3d 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -249,7 +249,7 @@ FE_Q::clone() const template void FE_Q:: -get_interpolation_matrix (const FiniteElementBase &x_source_fe, +get_interpolation_matrix (const FiniteElement &x_source_fe, FullMatrix &interpolation_matrix) const { // this is only implemented, if the @@ -258,7 +258,7 @@ get_interpolation_matrix (const FiniteElementBase &x_source_fe, AssertThrow ((x_source_fe.get_name().find ("FE_Q<") == 0) || (dynamic_cast*>(&x_source_fe) != 0), - typename FiniteElementBase:: + typename FiniteElement:: ExcInterpolationNotImplemented()); // ok, source is a Q element, so @@ -607,7 +607,7 @@ FE_Q<3>::initialize_constraints () // In the following the points x_i // are constructed in the order as // described in the documentation - // of the FiniteElementBase class + // of the FiniteElement class // (fe_base.h), i.e. // *--13--3--14--* // | | | @@ -765,7 +765,7 @@ FE_Q<3>::initialize_constraints () // // For a different explanation of // the problem, see the discussion - // in the FiniteElementBase class + // in the FiniteElement class // for constraint matrices in 3d. mirror[k] = (constraint_point(k) > 0.5); if (mirror[k]) diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc index e4b1df845e..6becd13f1a 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc @@ -287,7 +287,7 @@ FE_RaviartThomas::shape_grad_grad_component (const unsigned int i, template <> void FE_RaviartThomas<1>:: -get_interpolation_matrix (const FiniteElementBase<1> &, +get_interpolation_matrix (const FiniteElement<1> &, FullMatrix &) const { Assert (false, ExcImpossibleInDim(1)); @@ -299,7 +299,7 @@ get_interpolation_matrix (const FiniteElementBase<1> &, template void FE_RaviartThomas:: -get_interpolation_matrix (const FiniteElementBase &x_source_fe, +get_interpolation_matrix (const FiniteElement &x_source_fe, FullMatrix &interpolation_matrix) const { // this is only implemented, if the @@ -310,7 +310,7 @@ get_interpolation_matrix (const FiniteElementBase &x_source_fe, AssertThrow ((x_source_fe.get_name().find ("FE_RaviartThomas<") == 0) || (dynamic_cast*>(&x_source_fe) != 0), - typename FiniteElementBase:: + typename FiniteElement:: ExcInterpolationNotImplemented()); // ok, source is a RT element, so diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index f120aa76f1..9ed9b15973 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -62,7 +62,7 @@ FESystem::InternalData::~InternalData() template -typename FiniteElementBase::InternalDataBase & +typename FiniteElement::InternalDataBase & FESystem:: InternalData::get_fe_data (const unsigned int base_no) const { @@ -77,7 +77,7 @@ template void FESystem:: InternalData::set_fe_data (const unsigned int base_no, - typename FiniteElementBase::InternalDataBase *ptr) + typename FiniteElement::InternalDataBase *ptr) { Assert(base_no::InternalData::clear_first_cell () { // call respective function of base // class - FiniteElementBase::InternalDataBase::clear_first_cell (); + FiniteElement::InternalDataBase::clear_first_cell (); // then the functions of all the // sub-objects for (unsigned int i=0; i::shape_value (const unsigned int i, { Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); Assert (this->is_primitive(i), - typename FiniteElementBase::ExcShapeFunctionNotPrimitive(i)); + typename FiniteElement::ExcShapeFunctionNotPrimitive(i)); return (base_element(this->system_to_base_table[i].first.first) .shape_value(this->system_to_base_table[i].second, p)); @@ -351,7 +351,7 @@ FESystem::shape_grad (const unsigned int i, { Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); Assert (this->is_primitive(i), - typename FiniteElementBase::ExcShapeFunctionNotPrimitive(i)); + typename FiniteElement::ExcShapeFunctionNotPrimitive(i)); return (base_element(this->system_to_base_table[i].first.first) .shape_grad(this->system_to_base_table[i].second, p)); @@ -403,7 +403,7 @@ FESystem::shape_grad_grad (const unsigned int i, { Assert (idofs_per_cell, ExcIndexRange(i, 0, this->dofs_per_cell)); Assert (this->is_primitive(i), - typename FiniteElementBase::ExcShapeFunctionNotPrimitive(i)); + typename FiniteElement::ExcShapeFunctionNotPrimitive(i)); return (base_element(this->system_to_base_table[i].first.first) .shape_grad_grad(this->system_to_base_table[i].second, p)); @@ -451,7 +451,7 @@ FESystem::shape_grad_grad_component (const unsigned int i, template void FESystem:: -get_interpolation_matrix (const FiniteElementBase &x_source_fe, +get_interpolation_matrix (const FiniteElement &x_source_fe, FullMatrix &interpolation_matrix) const { Assert (interpolation_matrix.m() == this->dofs_per_cell, @@ -470,7 +470,7 @@ get_interpolation_matrix (const FiniteElementBase &x_source_fe, AssertThrow ((x_source_fe.get_name().find ("FESystem<") == 0) || (dynamic_cast*>(&x_source_fe) != 0), - typename FiniteElementBase:: + typename FiniteElement:: ExcInterpolationNotImplemented()); // ok, source is a system element, @@ -481,7 +481,7 @@ get_interpolation_matrix (const FiniteElementBase &x_source_fe, // condition 2: same number of // basis elements AssertThrow (n_base_elements() == source_fe.n_base_elements(), - typename FiniteElementBase:: + typename FiniteElement:: ExcInterpolationNotImplemented()); // condition 3: same number of @@ -489,7 +489,7 @@ get_interpolation_matrix (const FiniteElementBase &x_source_fe, for (unsigned int i=0; i:: + typename FiniteElement:: ExcInterpolationNotImplemented()); // ok, so let's try whether it @@ -617,8 +617,8 @@ FESystem::get_data (const UpdateFlags flags_, typename Mapping::InternalDataBase *base_fe_data_base = base_element(base_no).get_data(sub_flags, mapping, quadrature); - typename FiniteElementBase::InternalDataBase *base_fe_data = - dynamic_cast::InternalDataBase *> + typename FiniteElement::InternalDataBase *base_fe_data = + dynamic_cast::InternalDataBase *> (base_fe_data_base); data->set_fe_data(base_no, base_fe_data); @@ -823,7 +823,7 @@ compute_fill (const Mapping &mapping, { const FiniteElement & base_fe = base_element(base_no); - typename FiniteElementBase::InternalDataBase & + typename FiniteElement::InternalDataBase & base_fe_data = fe_data.get_fe_data(base_no); FEValuesData & base_data = fe_data.get_fe_values_data(base_no); @@ -2389,7 +2389,7 @@ FESystem::unit_support_point (const unsigned index) const ExcIndexRange (index, 0, this->dofs_per_cell)); Assert ((this->unit_support_points.size() == this->dofs_per_cell) || (this->unit_support_points.size() == 0), - typename FiniteElementBase::ExcFEHasNoSupportPoints ()); + typename FiniteElement::ExcFEHasNoSupportPoints ()); // let's see whether we have the // information pre-computed @@ -2413,7 +2413,7 @@ FESystem::unit_face_support_point (const unsigned index) const ExcIndexRange (index, 0, this->dofs_per_face)); Assert ((this->unit_face_support_points.size() == this->dofs_per_face) || (this->unit_face_support_points.size() == 0), - typename FiniteElementBase::ExcFEHasNoSupportPoints ()); + typename FiniteElement::ExcFEHasNoSupportPoints ()); // let's see whether we have the // information pre-computed diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index 92499e2ac0..961873332f 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -298,7 +298,7 @@ void FETools::get_interpolation_matrix (const FiniteElement &fe1, { gim_forwarder (fe1, fe2, interpolation_matrix); } - catch (typename FiniteElementBase::ExcInterpolationNotImplemented &) + catch (typename FiniteElement::ExcInterpolationNotImplemented &) { // too bad.... fe_implements_interpolation = false; @@ -321,7 +321,7 @@ void FETools::get_interpolation_matrix (const FiniteElement &fe1, fe2_support_points = fe2.get_unit_support_points (); Assert(fe2_support_points.size()==fe2.dofs_per_cell, - typename FiniteElementBase::ExcFEHasNoSupportPoints()); + typename FiniteElement::ExcFEHasNoSupportPoints()); for (unsigned int i=0; i +
  • Removed: The class FiniteElementBase has been removed and all its + functions are now in FiniteElement. +
    + (GK, 2005/08/25) +

    +
  • New: class DoFTools now has two functions compute_row_length_vector, one for equations and one for systems. These give a much fine estimate -- 2.39.5