From d3b6835acb694995a66159f7c39c0946b8177b9b Mon Sep 17 00:00:00 2001 From: bangerth Date: Sun, 26 Sep 2010 02:57:21 +0000 Subject: [PATCH] Make FETools a namespace, rather than a class with all static member functions. git-svn-id: https://svn.dealii.org/trunk@22162 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_tools.h | 2764 ++++++------ .../source/fe/fe_tools.all_dimensions.cc | 525 +-- deal.II/deal.II/source/fe/fe_tools.cc | 4020 ++++++++--------- deal.II/doc/news/changes.h | 6 + 4 files changed, 3659 insertions(+), 3656 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_tools.h b/deal.II/deal.II/include/fe/fe_tools.h index a2e2da8ed5..20dd758f09 100644 --- a/deal.II/deal.II/include/fe/fe_tools.h +++ b/deal.II/deal.II/include/fe/fe_tools.h @@ -47,7 +47,7 @@ class ConstraintMatrix; /** - * This class performs interpolations and extrapolations of discrete + * This namespace offers interpolations and extrapolations of discrete * functions of one @p FiniteElement @p fe1 to another @p FiniteElement * @p fe2. * @@ -63,1396 +63,1388 @@ class ConstraintMatrix; * @author Wolfgang Bangerth, Ralf Hartmann, Guido Kanschat; * 2000, 2003, 2004, 2005, 2006 */ -class FETools +namespace FETools { - public: - /** - * A base class for factory - * objects creating finite - * elements of a given - * degree. Derived classes are - * called whenever one wants to - * have a transparent way to - * create a finite element - * object. - * - * This class is used in the - * FETools::get_fe_from_name() - * and FETools::add_fe_name() - * functions. - * - * @author Guido Kanschat, 2006 - */ - template - class FEFactoryBase - { - public: - /** - * Create a FiniteElement and - * return a pointer to it. - */ - virtual FiniteElement* - get (const unsigned int degree) const = 0; - - /** - * Virtual destructor doing - * nothing but making the - * compiler happy. - */ - virtual ~FEFactoryBase(); - }; - - /** - * A concrete class for factory - * objects creating finite - * elements of a given degree. - * - * The class's get() function - * generates a finite element - * object of the type given as - * template argument, and with - * the degree (however the finite - * element class wishes to - * interpret this number) given - * as argument to get(). - * - * @author Guido Kanschat, 2006 - */ - template - class FEFactory : public FEFactoryBase - { - public: - /** - * Create a FiniteElement and - * return a pointer to it. - */ - virtual FiniteElement* - get (const unsigned int degree) const; - }; - - /** - * @warning In most cases, you - * will probably want to use - * compute_base_renumbering(). - * - * Compute the vector required to - * renumber the dofs of a cell by - * component. Furthermore, - * compute the vector storing the - * start indices of each - * component in the local block - * vector. - * - * The second vector is organized - * such that there is a vector - * for each base element - * containing the start index for - * each component served by this - * base element. - * - * While the first vector is - * checked to have the correct - * size, the second one is - * reinitialized for convenience. - */ - template - static void compute_component_wise( - const FiniteElement& fe, - std::vector& renumbering, - std::vector >& start_indices); - - /** - * Compute the vector required to - * renumber the dofs of a cell by - * block. Furthermore, compute - * the vector storing either the - * start indices or the size of - * each local block vector. - * - * If the @p bool parameter is - * true, @p block_data is filled - * with the start indices of each - * local block. If it is false, - * then the block sizes are - * returned. - * - * @todo Which way does this - * vector map the numbers? - */ - template - static void compute_block_renumbering ( - const FiniteElement& fe, - std::vector& renumbering, - std::vector& block_data, - bool return_start_indices = true); - - /** - * @name Generation of local matrices - * @{ - */ - /** - * Gives the interpolation matrix - * that interpolates a @p fe1- - * function to a @p fe2-function on - * each cell. The interpolation_matrix - * needs to be of size - * (fe2.dofs_per_cell, fe1.dofs_per_cell). - * - * Note, that if the finite element - * space @p fe1 is a subset of - * the finite element space - * @p fe2 then the @p interpolation_matrix - * is an embedding matrix. - */ - template - static - void - get_interpolation_matrix(const FiniteElement &fe1, - const FiniteElement &fe2, - FullMatrix &interpolation_matrix); - - /** - * Gives the interpolation matrix - * that interpolates a @p fe1- - * function to a @p fe2-function, and - * interpolates this to a second - * @p fe1-function on - * each cell. The interpolation_matrix - * needs to be of size - * (fe1.dofs_per_cell, fe1.dofs_per_cell). - * - * Note, that this function only - * makes sense if the finite element - * space due to @p fe1 is not a subset of - * the finite element space due to - * @p fe2, as if it were a subset then - * the @p interpolation_matrix would be - * only the unit matrix. - */ - template - static - void - get_back_interpolation_matrix(const FiniteElement &fe1, - const FiniteElement &fe2, - FullMatrix &interpolation_matrix); - - /** - * Gives the unit matrix minus the - * back interpolation matrix. - * The @p difference_matrix - * needs to be of size - * (fe1.dofs_per_cell, fe1.dofs_per_cell). - * - * This function gives - * the matrix that transforms a - * @p fe1 function $z$ to $z-I_hz$ - * where $I_h$ denotes the interpolation - * operator from the @p fe1 space to - * the @p fe2 space. This matrix hence - * is useful to evaluate - * error-representations where $z$ - * denotes the dual solution. - */ - template - static - void - get_interpolation_difference_matrix(const FiniteElement &fe1, - const FiniteElement &fe2, - FullMatrix &difference_matrix); - - /** - * Compute the local - * $L^2$-projection matrix from - * fe1 to fe2. - */ - template - static void get_projection_matrix(const FiniteElement &fe1, + /** + * A base class for factory + * objects creating finite + * elements of a given + * degree. Derived classes are + * called whenever one wants to + * have a transparent way to + * create a finite element + * object. + * + * This class is used in the + * FETools::get_fe_from_name() + * and FETools::add_fe_name() + * functions. + * + * @author Guido Kanschat, 2006 + */ + template + class FEFactoryBase + { + public: + /** + * Create a FiniteElement and + * return a pointer to it. + */ + virtual FiniteElement* + get (const unsigned int degree) const = 0; + + /** + * Virtual destructor doing + * nothing but making the + * compiler happy. + */ + virtual ~FEFactoryBase(); + }; + + /** + * A concrete class for factory + * objects creating finite + * elements of a given degree. + * + * The class's get() function + * generates a finite element + * object of the type given as + * template argument, and with + * the degree (however the finite + * element class wishes to + * interpret this number) given + * as argument to get(). + * + * @author Guido Kanschat, 2006 + */ + template + class FEFactory : public FEFactoryBase + { + public: + /** + * Create a FiniteElement and + * return a pointer to it. + */ + virtual FiniteElement* + get (const unsigned int degree) const; + }; + + /** + * @warning In most cases, you + * will probably want to use + * compute_base_renumbering(). + * + * Compute the vector required to + * renumber the dofs of a cell by + * component. Furthermore, + * compute the vector storing the + * start indices of each + * component in the local block + * vector. + * + * The second vector is organized + * such that there is a vector + * for each base element + * containing the start index for + * each component served by this + * base element. + * + * While the first vector is + * checked to have the correct + * size, the second one is + * reinitialized for convenience. + */ + template + void compute_component_wise( + const FiniteElement& fe, + std::vector& renumbering, + std::vector >& start_indices); + + /** + * Compute the vector required to + * renumber the dofs of a cell by + * block. Furthermore, compute + * the vector storing either the + * start indices or the size of + * each local block vector. + * + * If the @p bool parameter is + * true, @p block_data is filled + * with the start indices of each + * local block. If it is false, + * then the block sizes are + * returned. + * + * @todo Which way does this + * vector map the numbers? + */ + template + void compute_block_renumbering ( + const FiniteElement& fe, + std::vector& renumbering, + std::vector& block_data, + bool return_start_indices = true); + + /** + * @name Generation of local matrices + * @{ + */ + /** + * Gives the interpolation matrix + * that interpolates a @p fe1- + * function to a @p fe2-function on + * each cell. The interpolation_matrix + * needs to be of size + * (fe2.dofs_per_cell, fe1.dofs_per_cell). + * + * Note, that if the finite element + * space @p fe1 is a subset of + * the finite element space + * @p fe2 then the @p interpolation_matrix + * is an embedding matrix. + */ + template + void + get_interpolation_matrix(const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix); + + /** + * Gives the interpolation matrix + * that interpolates a @p fe1- + * function to a @p fe2-function, and + * interpolates this to a second + * @p fe1-function on + * each cell. The interpolation_matrix + * needs to be of size + * (fe1.dofs_per_cell, fe1.dofs_per_cell). + * + * Note, that this function only + * makes sense if the finite element + * space due to @p fe1 is not a subset of + * the finite element space due to + * @p fe2, as if it were a subset then + * the @p interpolation_matrix would be + * only the unit matrix. + */ + template + void + get_back_interpolation_matrix(const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix); + + /** + * Gives the unit matrix minus the + * back interpolation matrix. + * The @p difference_matrix + * needs to be of size + * (fe1.dofs_per_cell, fe1.dofs_per_cell). + * + * This function gives + * the matrix that transforms a + * @p fe1 function $z$ to $z-I_hz$ + * where $I_h$ denotes the interpolation + * operator from the @p fe1 space to + * the @p fe2 space. This matrix hence + * is useful to evaluate + * error-representations where $z$ + * denotes the dual solution. + */ + template + void + get_interpolation_difference_matrix(const FiniteElement &fe1, const FiniteElement &fe2, - FullMatrix &matrix); - - /** - * Compute the matrix of nodal - * values of a finite element - * applied to all its shape - * functions. - * - * This function is supposed to - * help building finite elements - * from polynomial spaces and - * should be called inside the - * constructor of an - * element. Applied to a - * completely initialized finite - * element, the result should be - * the unit matrix by definition - * of the node values. - * - * Using this matrix allows the - * construction of the basis of - * shape functions in two steps. - *
    - * - *
  1. Define the space of shape - * functions using an arbitrary - * basis wj and - * compute the matrix M of - * node functionals - * Ni applied - * to these basis functions. - * - *
  2. Compute the basis - * vj of the - * finite element shape function - * space by applying - * M-1 to the - * basis wj. - *
- * - * @note The FiniteElement must - * provide generalized support - * points and and interpolation - * functions. - */ - template - static void compute_node_matrix(FullMatrix& M, - const FiniteElement& fe); - - /** - * For all possible (isotropic - * and anisotropic) refinement - * cases compute the embedding - * matrices from a coarse cell to - * the child cells. Each column - * of the resulting matrices - * contains the representation of - * a coarse grid basis functon by - * the fine grid basis; the - * matrices are split such that - * there is one matrix for every - * child. - * - * This function computes the - * coarse grid function in a - * sufficiently large number of - * quadrature points and fits the - * fine grid functions using - * least squares - * approximation. Therefore, the - * use of this function is - * restricted to the case that - * the finite element spaces are - * actually nested. - * - * Note, that - * matrices[refinement_case-1][child] - * includes the embedding (or prolongation) - * matrix of child - * child for the - * RefinementCase - * refinement_case. Here, - * we use - * refinement_case-1 - * instead of - * refinement_case - * as for - * RefinementCase::no_refinement(=0) - * there are no prolongation - * matrices available. - * - * Typically this function is - * called by the various - * implementations of - * FiniteElement classes in order - * to fill the respective - * FiniteElement::prolongation - * matrices. - * - * @param fe The finite element - * class for which we compute the - * embedding matrices. - * - * @param matrices A reference to - * RefinementCase::isotropic_refinement - * vectors of FullMatrix - * objects. Each vector - * corresponds to one - * RefinementCase @p - * refinement_case and is of the - * vector size - * GeometryInfo::n_children(refinement_case). This - * is the format used in - * FiniteElement, where we want - * to use this function mostly. - * - * @param isotropic_only Set - * to true if you only - * want to compute matrices for - * isotropic refinement. - */ - template - static void compute_embedding_matrices(const FiniteElement &fe, - std::vector > >& matrices, - const bool isotropic_only = false); - - /** - * Compute the embedding matrices - * on faces needed for constraint - * matrices. - * - * @param fe The finite element - * for which to compute these - * matrices. @param matrices An - * array of - * GeometryInfo::subfaces_per_face - * = 2dim-1 - * FullMatrix objects,holding the - * embedding matrix for each - * subface. @param face_coarse - * The number of the face on the - * coarse side of the face for - * which this is computed. - * @param face_fine The number of - * the face on the refined side - * of the face for which this is - * computed. - * - * @warning This function will be - * used in computing constraint - * matrices. It is not - * sufficiently tested yet. - */ - template - static void - compute_face_embedding_matrices(const FiniteElement& fe, - FullMatrix (&matrices)[GeometryInfo::max_children_per_face], - const unsigned int face_coarse, - const unsigned int face_fine); - - /** - * For all possible (isotropic - * and anisotropic) refinement - * cases compute the - * L2-projection - * matrices from the children to - * a coarse cell. - * - * Note, that - * matrices[refinement_case-1][child] - * includes the projection (or restriction) - * matrix of child - * child for the - * RefinementCase - * refinement_case. Here, - * we use - * refinement_case-1 - * instead of - * refinement_case - * as for - * RefinementCase::no_refinement(=0) - * there are no projection - * matrices available. - * - * Typically this function is - * called by the various - * implementations of - * FiniteElement classes in order - * to fill the respective - * FiniteElement::restriction - * matrices. - * - * @arg fe The finite element - * class for which we compute the - * projection matrices. @arg - * matrices A reference to - * RefinementCase::isotropic_refinement - * vectors of FullMatrix - * objects. Each vector - * corresponds to one - * RefinementCase @p - * refinement_case and is of the - * vector size - * GeometryInfo::n_children(refinement_case). This - * is the format used in - * FiniteElement, where we want - * to use this function mostly. - * - * @arg isotropic_only Set - * to true if you only - * want to compute matrices for - * isotropic refinement. - */ - template - static void compute_projection_matrices( - const FiniteElement &fe, - std::vector > >& matrices, - const bool isotropic_only = false); - - /** - * Projects scalar data defined in - * quadrature points to a finite element - * space on a single cell. - * - * What this function does is the - * following: assume that there is scalar - * data uq, 0 <= q < - * Q:=quadrature.size() - * defined at the quadrature points of a - * cell, with the points defined by the - * given rhs_quadrature - * object. We may then want to ask for - * that finite element function (on a - * single cell) vh in - * the finite-dimensional space defined - * by the given FE object that is the - * projection of u in the - * following sense: - * - * Usually, the projection - * vh is that - * function that satisfies - * (vh,w)=(u,w) for - * all discrete test functions - * w. In the present case, we - * can't evaluate the right hand side, - * since u is only defined in - * the quadrature points given by - * rhs_quadrature, so we replace - * it by a quadrature - * approximation. Likewise, the left hand - * side is approximated using the - * lhs_quadrature object; if - * this quadrature object is chosen - * appropriately, then the integration of - * the left hand side can be done - * exactly, without any - * approximation. The use of different - * quadrature objects is necessary if the - * quadrature object for the right hand - * side has too few quadrature points -- - * for example, if data q is - * only defined at the cell center, then - * the corresponding one-point quadrature - * formula is obviously insufficient to - * approximate the scalar product on the - * left hand side by a definite form. - * - * After these quadrature approximations, - * we end up with a nodal representation - * Vh of - * vh that satisfies - * the following system of linear - * equations: M Vh = Q - * U, where - * Mij=(phi_i,phi_j) - * is the mass matrix approximated by - * lhs_quadrature, and - * Q is the matrix - * Qiq=phii(xq) - * wq where - * wq are quadrature - * weights; U is the vector of - * quadrature point data - * uq. - * - * In order to then get the nodal - * representation Vh - * of the projection of U, one - * computes Vh = X U, - * X=M-1 Q. The purpose - * of this function is to compute the - * matrix X and return it - * through the last argument of this - * function. - * - * Note that this function presently only - * supports scalar data. An extension of - * the mass matrix is of course trivial, - * but one has to define the order of - * data in the vector U if it - * contains vector valued data in all - * quadrature points. - * - * A use for this function is described - * in the introduction to the step-18 - * example program. - * - * The opposite of this function, - * interpolation of a finite element - * function onto quadrature points is - * essentially what the - * FEValues::get_function_values - * functions do; to make things a little - * simpler, the - * FETools::compute_interpolation_to_quadrature_points_matrix - * provides the matrix form of this. - * - * Note that this function works - * on a single cell, rather than - * an entire triangulation. In - * effect, it therefore doesn't - * matter if you use a continuous - * or discontinuous version of - * the finite element. - * - * It is worth noting that there - * are a few confusing cases of - * this function. The first one - * is that it really only makes - * sense to project onto a finite - * element that has at most as - * many degrees of freedom per - * cell as there are quadrature - * points; the projection of N - * quadrature point data into a - * space with M>N unknowns is - * well-defined, but often yields - * funny and non-intuitive - * results. Secondly, one would - * think that if the quadrature - * point data is defined in the - * support points of the finite - * element, i.e. the quadrature - * points of - * ths_quadrature equal - * fe.get_unit_support_points(), - * then the projection should be - * the identity, i.e. each degree - * of freedom of the finite - * element equals the value of - * the given data in the support - * point of the corresponding - * shape function. However, this - * is not generally the case: - * while the matrix Q in - * that case is the identity - * matrix, the mass matrix - * M is not equal to the - * identity matrix, except for - * the special case that the - * quadrature formula - * lhs_quadrature also - * has its quadrature points in - * the support points of the - * finite element. - * - * Finally, this function only defines a - * cell wise projection, while one - * frequently wants to apply it to all - * cells in a triangulation. However, if - * it is applied to one cell after the - * other, the results from later cells - * may overwrite nodal values computed - * already from previous cells if degrees - * of freedom live on the interfaces - * between cells. The function is - * therefore most useful for - * discontinuous elements. - */ - template - static - void - compute_projection_from_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &lhs_quadrature, - const Quadrature &rhs_quadrature, - FullMatrix &X); - - /** - * Given a (scalar) local finite element - * function, compute the matrix that maps - * the vector of nodal values onto the - * vector of values of this function at - * quadrature points as given by the - * second argument. In a sense, this - * function does the opposite of the @p - * compute_projection_from_quadrature_points_matrix - * function. - */ - template - static - void - compute_interpolation_to_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &quadrature, - FullMatrix &I_q); - - /** - * Computes the projection of tensorial - * (first-order tensor) - * data stored at the quadrature points - * @p vector_of_tensors_at_qp - * to data @p vector_of_tensors_at_nodes - * at the support points of the cell. - * The data in - * @p vector_of_tensors_at_qp - * is ordered sequentially following the - * quadrature point numbering. - * The size of - * @p vector_of_tensors_at_qp - * must correspond to the number of columns - * of @p projection_matrix. - * The size of @p vector_of_tensors_at_nodes - * must correspond to the number of rows of - * @p vector_of_tensors_at_nodes . - * The projection matrix - * @p projection_matrix desribes the - * projection of scalar data from the - * quadrature points and can be obtained - * from the - * FETools::compute_projection_from_quadrature_points_matrix - * function. - */ - template - static - void - compute_projection_from_quadrature_points( - const FullMatrix &projection_matrix, - const std::vector< Tensor<1, dim > > &vector_of_tensors_at_qp, - std::vector< Tensor<1, dim > > &vector_of_tensors_at_nodes); - - - - /** - * same as last function but for a - * @p SymmetricTensor . - */ - template - static - void - compute_projection_from_quadrature_points( - const FullMatrix &projection_matrix, - const std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_qp, - std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_nodes); - - - - - /** - * This method implements the - * FETools::compute_projection_from_quadrature_points_matrix - * method for faces of a mesh. - * The matrix that it returns, X, is face specific - * and its size is fe.dofs_per_cell by - * rhs_quadrature.size(). - * The dimension, dim must be larger than 1 for this class, - * since Quadrature objects are required. See the - * documentation on the Quadrature class for more information. - */ - template - static - void - compute_projection_from_face_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &lhs_quadrature, - const Quadrature &rhs_quadrature, - const typename DoFHandler::active_cell_iterator & cell, - unsigned int face, - FullMatrix &X); - - - - //@} - /** - * @name Functions which should be in DoFTools - */ - //@{ - /** - * Gives the interpolation of a the - * @p dof1-function @p u1 to a - * @p dof2-function @p u2. @p dof1 and - * @p dof2 need to be DoFHandlers - * based on the same triangulation. - * - * If the elements @p fe1 and @p fe2 - * are either both continuous or - * both discontinuous then this - * interpolation is the usual point - * interpolation. The same is true - * if @p fe1 is a continuous and - * @p fe2 is a discontinuous finite - * element. For the case that @p fe1 - * is a discontinuous and @p fe2 is - * a continuous finite element - * there is no point interpolation - * defined at the discontinuities. - * Therefore the meanvalue is taken - * at the DoF values on the - * discontinuities. - * - * Note that for continuous - * elements on grids with hanging - * nodes (i.e. locally refined - * grids) this function does not - * give the expected output. - * Indeed, the resulting output - * vector does not necessarily - * respect continuity - * requirements at hanging nodes: - * if, for example, you are - * interpolating a Q2 field to a - * Q1 field, then at hanging - * nodes the output field will - * have the function value of the - * input field, which however is - * not usually the mean value of - * the two adjacent nodes. It is - * thus not part of the Q1 - * function space on the whole - * triangulation, although it is - * of course Q1 on each cell. - * - * For this case (continuous - * elements on grids with hanging - * nodes), please use the - * @p interpolate function with - * an additional - * @p ConstraintMatrix argument, - * see below, or make the field - * conforming yourself by calling - * the @p distribute function of - * your hanging node constraints - * object. - */ - template class DH1, - template class DH2, - class InVector, class OutVector> - static - void - interpolate (const DH1 &dof1, - const InVector &u1, - const DH2 &dof2, - OutVector &u2); - - /** - * Gives the interpolation of a - * the @p dof1-function @p u1 to - * a @p dof2-function @p u2. @p - * dof1 and @p dof2 need to be - * DoFHandlers (or - * hp::DoFHandlers) based on the - * same triangulation. @p - * constraints is a hanging node - * constraints object - * corresponding to @p dof2. This - * object is particular important - * when interpolating onto - * continuous elements on grids - * with hanging nodes (locally - * refined grids). - * - * If the elements @p fe1 and @p fe2 - * are either both continuous or - * both discontinuous then this - * interpolation is the usual point - * interpolation. The same is true - * if @p fe1 is a continuous and - * @p fe2 is a discontinuous finite - * element. For the case that @p fe1 - * is a discontinuous and @p fe2 is - * a continuous finite element - * there is no point interpolation - * defined at the discontinuities. - * Therefore the meanvalue is taken - * at the DoF values on the - * discontinuities. - */ - template class DH1, - template class DH2, - class InVector, class OutVector> - static void interpolate (const DH1 &dof1, - const InVector &u1, - const DH2 &dof2, - const ConstraintMatrix &constraints, - OutVector& u2); - - /** - * Gives the interpolation of the - * @p fe1-function @p u1 to a - * @p fe2-function, and - * interpolates this to a second - * @p fe1-function named - * @p u1_interpolated. - * - * Note, that this function does - * not work on continuous - * elements at hanging nodes. For - * that case use the - * @p back_interpolate function, - * below, that takes an - * additional - * @p ConstraintMatrix object. - * - * Furthermore note, that for the - * specific case when the finite - * element space corresponding to - * @p fe1 is a subset of the - * finite element space - * corresponding to @p fe2, this - * function is simply an identity - * mapping. - */ - template - static void back_interpolate (const DoFHandler &dof1, - const InVector &u1, - const FiniteElement &fe2, - OutVector &u1_interpolated); - - /** - * Same as last function, except - * that the dof handler objects - * might be of type - * @p hp::DoFHandler. - */ - template class DH, - class InVector, class OutVector, int spacedim> - static void back_interpolate (const DH &dof1, - const InVector &u1, - const FiniteElement &fe2, - OutVector &u1_interpolated); - - /** - * Gives the interpolation of the - * @p dof1-function @p u1 to a - * @p dof2-function, and - * interpolates this to a second - * @p dof1-function named - * @p u1_interpolated. - * @p constraints1 and - * @p constraints2 are the - * hanging node constraints - * corresponding to @p dof1 and - * @p dof2, respectively. These - * objects are particular - * important when continuous - * elements on grids with hanging - * nodes (locally refined grids) - * are involved. - * - * Furthermore note, that for the - * specific case when the finite - * element space corresponding to - * @p dof1 is a subset of the - * finite element space - * corresponding to @p dof2, this - * function is simply an identity - * mapping. - */ - template - static void back_interpolate (const DoFHandler& dof1, - const ConstraintMatrix& constraints1, - const InVector& u1, - const DoFHandler& dof2, - const ConstraintMatrix& constraints2, - OutVector& u1_interpolated); - - /** - * Gives $(Id-I_h)z_1$ for a given - * @p dof1-function $z_1$, where $I_h$ - * is the interpolation from @p fe1 - * to @p fe2. The result $(Id-I_h)z_1$ is - * written into @p z1_difference. - * - * Note, that this function does - * not work for continuous - * elements at hanging nodes. For - * that case use the - * @p interpolation_difference - * function, below, that takes an - * additional - * @p ConstraintMatrix object. - */ - template - static void interpolation_difference(const DoFHandler &dof1, - const InVector &z1, - const FiniteElement &fe2, - OutVector &z1_difference); - - /** - * Gives $(Id-I_h)z_1$ for a given - * @p dof1-function $z_1$, where $I_h$ - * is the interpolation from @p fe1 - * to @p fe2. The result $(Id-I_h)z_1$ is - * written into @p z1_difference. - * @p constraints1 and - * @p constraints2 are the - * hanging node constraints - * corresponding to @p dof1 and - * @p dof2, respectively. These - * objects are particular - * important when continuous - * elements on grids with hanging - * nodes (locally refined grids) - * are involved. - */ - template - static void interpolation_difference(const DoFHandler& dof1, - const ConstraintMatrix& constraints1, - const InVector& z1, - const DoFHandler& dof2, - const ConstraintMatrix& constraints2, - OutVector& z1_difference); - - /** - * $L^2$ projection for - * discontinuous - * elements. Operates the same - * direction as interpolate. - * - * The global projection can be - * computed by local matrices if - * the finite element spaces are - * discontinuous. With continuous - * elements, this is impossible, - * since a global mass matrix - * must be inverted. - */ - template - static void project_dg (const DoFHandler& dof1, - const InVector& u1, - const DoFHandler& dof2, - OutVector& u2); - - /** - * Gives the patchwise - * extrapolation of a @p dof1 - * function @p z1 to a @p dof2 - * function @p z2. @p dof1 and - * @p dof2 need to be DoFHandler - * based on the same triangulation. - * - * This function is interesting - * for e.g. extrapolating - * patchwise a piecewise linear - * solution to a piecewise - * quadratic solution. - * - * Note that the resulting field - * does not satisfy continuity - * requirements of the given - * finite elements. - * - * When you use continuous - * elements on grids with hanging - * nodes, please use the - * @p extrapolate function with - * an additional - * ConstraintMatrix argument, - * see below. - * - * Since this function operates - * on patches of cells, it is - * required that the underlying - * grid is refined at least once - * for every coarse grid cell. If - * this is not the case, an - * exception will be raised. - */ - template - static void extrapolate (const DoFHandler& dof1, - const InVector& z1, - const DoFHandler& dof2, - OutVector& z2); - - /** - * Gives the patchwise - * extrapolation of a @p dof1 - * function @p z1 to a @p dof2 - * function @p z2. @p dof1 and - * @p dof2 need to be DoFHandler - * based on the same triangulation. - * @p constraints is a hanging - * node constraints object - * corresponding to - * @p dof2. This object is - * particular important when - * interpolating onto continuous - * elements on grids with hanging - * nodes (locally refined grids). - * - * Otherwise, the same holds as - * for the other @p extrapolate - * function. - */ - template - static void extrapolate (const DoFHandler& dof1, - const InVector& z1, - const DoFHandler& dof2, - const ConstraintMatrix& constraints, - OutVector& z2); - //@} - /** - * The numbering of the degrees - * of freedom in continous finite - * elements is hierarchic, - * i.e. in such a way that we - * first number the vertex dofs, - * in the order of the vertices - * as defined by the - * triangulation, then the line - * dofs in the order and - * respecting the direction of - * the lines, then the dofs on - * quads, etc. However, we could - * have, as well, numbered them - * in a lexicographic way, - * i.e. with indices first - * running in x-direction, then - * in y-direction and finally in - * z-direction. Discontinuous - * elements of class FE_DGQ() - * are numbered in this way, for - * example. - * - * This function constructs a - * table which lexicographic - * index each degree of freedom - * in the hierarchic numbering - * would have. It operates on the - * continuous finite element - * given as first argument, and - * outputs the lexicographic - * indices in the second. - * - * Note that since this function - * uses specifics of the - * continuous finite elements, it - * can only operate on - * FiniteElementData objects - * inherent in FE_Q(). However, - * this function does not take a - * FE_Q object as it is also - * invoked by the FE_Q() - * constructor. - * - * It is assumed that the size of - * the output argument already - * matches the correct size, - * which is equal to the number - * of degrees of freedom in the - * finite element. - */ - template - static void - hierarchic_to_lexicographic_numbering (const FiniteElementData &fe_data, - std::vector &h2l); - - /** - * Like the previous function but - * instead of returning its - * result through the last - * argument return it as a value. - */ - template - static - std::vector - hierarchic_to_lexicographic_numbering (const FiniteElementData &fe_data); - - /** - * This is the reverse function - * to the above one, generating - * the map from the lexicographic - * to the hierarchical - * numbering. All the remarks - * made about the above function - * are also valid here. - */ - template - static void - lexicographic_to_hierarchic_numbering (const FiniteElementData &fe_data, - std::vector &l2h); - - /** - * Like the previous function but - * instead of returning its - * result through the last - * argument return it as a value. - */ - template - static - std::vector - lexicographic_to_hierarchic_numbering (const FiniteElementData &fe_data); - - /** - * Parse the name of a finite - * element and generate a finite - * element object accordingly. - * - * The name must be in the form which - * is returned by the - * FiniteElement::get_name - * function, where a few - * modifications are allowed: - * - *
  • Dimension template - * parameters <2> etc. can - * be omitted. Alternatively, the - * explicit number can be - * replaced by dim or - * d. If a number is - * given, it must match - * the template parameter of this - * function. - * - *
  • The powers used for - * FESystem may either be numbers - * or can be - * replaced by dim or - * d. - *
- * - * If no finite element can be - * reconstructed from this - * string, an exception of type - * @p FETools::ExcInvalidFEName - * is thrown. - * - * The function returns a pointer - * to a newly create finite - * element. It is in the caller's - * responsibility to destroy the - * object pointed to at an - * appropriate later time. - * - * Since the value of the template - * argument can't be deduced from the - * (string) argument given to this - * function, you have to explicitly - * specify it when you call this - * function. - * - * This function knows about all - * the standard elements defined - * in the library. However, it - * doesn't by default know about - * elements that you may have - * defined in your program. To - * make your own elements known - * to this function, use the - * add_fe_name() function. - * This function does not work - * if one wants to get a codimension - * 1 finite element. - */ - template - static - FiniteElement * - get_fe_from_name (const std::string &name); - - - /** - * Extend the list of finite - * elements that can be generated - * by get_fe_from_name() by the - * one given as @p name. If - * get_fe_from_name() is later - * called with this name, it will - * use the object given as second - * argument to create a finite - * element object. - * - * The format of the @p name - * parameter should include the - * name of a finite - * element. However, it is safe - * to use either the class name - * alone or to use the result of - * FiniteElement::get_name (which - * includes the space dimension - * as well as the polynomial - * degree), since everything - * after the first non-name - * character will be ignored. - * - * The FEFactory object should be - * an object newly created with - * new. FETools will - * take ownership of this object - * and delete it once it is not - * used anymore. - * - * In most cases, if you want - * objects of type - * MyFE be created - * whenever the name - * my_fe is given to - * get_fe_from_name, you will - * want the second argument to - * this function be of type - * FEFactory@, but you can - * of course create your custom - * finite element factory class. - * - * This function takes over - * ownership of the object given - * as second argument, i.e. you - * should never attempt to - * destroy it later on. The - * object will be deleted at the - * end of the program's lifetime. - * - * If the name of the element - * is already in use, an exception - * is thrown. Thus, functionality - * of get_fe_from_name() can only - * be added, not changed. - * - * @note This function - * manipulates a global table - * (one table for each space - * dimension). It is thread safe - * in the sense that every access - * to this table is secured by a - * lock. Nevertheless, since each - * name can be added only once, - * user code has to make sure - * that only one thread adds a - * new element. - * - * Note also that this table - * exists once for each space - * dimension. If you have a - * program that works with finite - * elements in different space - * dimensions (for example, @ref - * step_4 "step-4" does something - * like this), then you should - * call this function for each - * space dimension for which you - * want your finite element added - * to the map. - */ - template - static void add_fe_name (const std::string& name, - const FEFactoryBase* factory); - - /** - * The string used for - * get_fe_from_name() cannot be - * translated to a finite - * element. - * - * Either the string is badly - * formatted or you are using a - * custom element that must be - * added using add_fe_name() - * first. - * - * @ingroup Exceptions - */ - DeclException1 (ExcInvalidFEName, - std::string, - << "Can't re-generate a finite element from the string '" - << arg1 << "'."); - - /** - * The string used for - * get_fe_from_name() cannot be - * translated to a finite - * element. - * - * Dimension arguments in finite - * element names should be - * avoided. If they are there, - * the dimension should be - * dim or - * d. Here, you gave a - * numeric dimension argument, - * which does not match the - * template dimension of the - * finite element class. - * - * @ingroup Exceptions - */ - DeclException2 (ExcInvalidFEDimension, - char, int, - << "The dimension " << arg1 - << " in the finite element string must match " - << "the space dimension " - << arg2 << "."); - - /** - * Exception - * - * @ingroup Exceptions - */ - DeclException0 (ExcInvalidFE); - - /** - * The finite element must be - * @ref GlossPrimitive "primitive". - * - * @ingroup Exceptions - */ - DeclException0 (ExcFENotPrimitive); - /** - * Exception - * - * @ingroup Exceptions - */ - DeclException0 (ExcTriangulationMismatch); - - /** - * A continuous element is used - * on a mesh with hanging nodes, - * but the constraint matrices - * are missing. - * - * @ingroup Exceptions - */ - DeclException1 (ExcHangingNodesNotAllowed, - int, - << "You are using continuous elements on a grid with " - << "hanging nodes but without providing hanging node " - << "constraints. Use the respective function with " - << "additional ConstraintMatrix argument(s), instead."); - /** - * You need at least two grid levels. - * - * @ingroup Exceptions - */ - DeclException0 (ExcGridNotRefinedAtLeastOnce); - /** - * The dimensions of the matrix - * used did not match the - * expected dimensions. - * - * @ingroup Exceptions - */ - DeclException4 (ExcMatrixDimensionMismatch, - int, int, int, int, - << "This is a " << arg1 << "x" << arg2 << " matrix, " - << "but should be a " << arg3 << "x" << arg4 << " matrix."); - - /** - * Exception thrown if an - * embedding matrix was computed - * inaccurately. - * - * @ingroup Exceptions - */ - DeclException1(ExcLeastSquaresError, double, - << "Least squares fit leaves a gap of " << arg1); - - /** - * Exception thrown if one variable - * may not be greater than another. - * - * @ingroup Exceptions - */ - DeclException2 (ExcNotGreaterThan, - int, int, - << arg1 << " must be greater than " << arg2); -}; - - -template -FiniteElement* -FETools::FEFactory::get (const unsigned int degree) const -{ - return new FE(degree); + FullMatrix &difference_matrix); + + /** + * Compute the local + * $L^2$-projection matrix from + * fe1 to fe2. + */ + template + void get_projection_matrix(const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &matrix); + + /** + * Compute the matrix of nodal + * values of a finite element + * applied to all its shape + * functions. + * + * This function is supposed to + * help building finite elements + * from polynomial spaces and + * should be called inside the + * constructor of an + * element. Applied to a + * completely initialized finite + * element, the result should be + * the unit matrix by definition + * of the node values. + * + * Using this matrix allows the + * construction of the basis of + * shape functions in two steps. + *
    + * + *
  1. Define the space of shape + * functions using an arbitrary + * basis wj and + * compute the matrix M of + * node functionals + * Ni applied + * to these basis functions. + * + *
  2. Compute the basis + * vj of the + * finite element shape function + * space by applying + * M-1 to the + * basis wj. + *
+ * + * @note The FiniteElement must + * provide generalized support + * points and and interpolation + * functions. + */ + template + void compute_node_matrix(FullMatrix& M, + const FiniteElement& fe); + + /** + * For all possible (isotropic + * and anisotropic) refinement + * cases compute the embedding + * matrices from a coarse cell to + * the child cells. Each column + * of the resulting matrices + * contains the representation of + * a coarse grid basis functon by + * the fine grid basis; the + * matrices are split such that + * there is one matrix for every + * child. + * + * This function computes the + * coarse grid function in a + * sufficiently large number of + * quadrature points and fits the + * fine grid functions using + * least squares + * approximation. Therefore, the + * use of this function is + * restricted to the case that + * the finite element spaces are + * actually nested. + * + * Note, that + * matrices[refinement_case-1][child] + * includes the embedding (or prolongation) + * matrix of child + * child for the + * RefinementCase + * refinement_case. Here, + * we use + * refinement_case-1 + * instead of + * refinement_case + * as for + * RefinementCase::no_refinement(=0) + * there are no prolongation + * matrices available. + * + * Typically this function is + * called by the various + * implementations of + * FiniteElement classes in order + * to fill the respective + * FiniteElement::prolongation + * matrices. + * + * @param fe The finite element + * class for which we compute the + * embedding matrices. + * + * @param matrices A reference to + * RefinementCase::isotropic_refinement + * vectors of FullMatrix + * objects. Each vector + * corresponds to one + * RefinementCase @p + * refinement_case and is of the + * vector size + * GeometryInfo::n_children(refinement_case). This + * is the format used in + * FiniteElement, where we want + * to use this function mostly. + * + * @param isotropic_only Set + * to true if you only + * want to compute matrices for + * isotropic refinement. + */ + template + void compute_embedding_matrices(const FiniteElement &fe, + std::vector > >& matrices, + const bool isotropic_only = false); + + /** + * Compute the embedding matrices + * on faces needed for constraint + * matrices. + * + * @param fe The finite element + * for which to compute these + * matrices. @param matrices An + * array of + * GeometryInfo::subfaces_per_face + * = 2dim-1 + * FullMatrix objects,holding the + * embedding matrix for each + * subface. @param face_coarse + * The number of the face on the + * coarse side of the face for + * which this is computed. + * @param face_fine The number of + * the face on the refined side + * of the face for which this is + * computed. + * + * @warning This function will be + * used in computing constraint + * matrices. It is not + * sufficiently tested yet. + */ + template + void + compute_face_embedding_matrices(const FiniteElement& fe, + FullMatrix (&matrices)[GeometryInfo::max_children_per_face], + const unsigned int face_coarse, + const unsigned int face_fine); + + /** + * For all possible (isotropic + * and anisotropic) refinement + * cases compute the + * L2-projection + * matrices from the children to + * a coarse cell. + * + * Note, that + * matrices[refinement_case-1][child] + * includes the projection (or restriction) + * matrix of child + * child for the + * RefinementCase + * refinement_case. Here, + * we use + * refinement_case-1 + * instead of + * refinement_case + * as for + * RefinementCase::no_refinement(=0) + * there are no projection + * matrices available. + * + * Typically this function is + * called by the various + * implementations of + * FiniteElement classes in order + * to fill the respective + * FiniteElement::restriction + * matrices. + * + * @arg fe The finite element + * class for which we compute the + * projection matrices. @arg + * matrices A reference to + * RefinementCase::isotropic_refinement + * vectors of FullMatrix + * objects. Each vector + * corresponds to one + * RefinementCase @p + * refinement_case and is of the + * vector size + * GeometryInfo::n_children(refinement_case). This + * is the format used in + * FiniteElement, where we want + * to use this function mostly. + * + * @arg isotropic_only Set + * to true if you only + * want to compute matrices for + * isotropic refinement. + */ + template + void compute_projection_matrices( + const FiniteElement &fe, + std::vector > >& matrices, + const bool isotropic_only = false); + + /** + * Projects scalar data defined in + * quadrature points to a finite element + * space on a single cell. + * + * What this function does is the + * following: assume that there is scalar + * data uq, 0 <= q < + * Q:=quadrature.size() + * defined at the quadrature points of a + * cell, with the points defined by the + * given rhs_quadrature + * object. We may then want to ask for + * that finite element function (on a + * single cell) vh in + * the finite-dimensional space defined + * by the given FE object that is the + * projection of u in the + * following sense: + * + * Usually, the projection + * vh is that + * function that satisfies + * (vh,w)=(u,w) for + * all discrete test functions + * w. In the present case, we + * can't evaluate the right hand side, + * since u is only defined in + * the quadrature points given by + * rhs_quadrature, so we replace + * it by a quadrature + * approximation. Likewise, the left hand + * side is approximated using the + * lhs_quadrature object; if + * this quadrature object is chosen + * appropriately, then the integration of + * the left hand side can be done + * exactly, without any + * approximation. The use of different + * quadrature objects is necessary if the + * quadrature object for the right hand + * side has too few quadrature points -- + * for example, if data q is + * only defined at the cell center, then + * the corresponding one-point quadrature + * formula is obviously insufficient to + * approximate the scalar product on the + * left hand side by a definite form. + * + * After these quadrature approximations, + * we end up with a nodal representation + * Vh of + * vh that satisfies + * the following system of linear + * equations: M Vh = Q + * U, where + * Mij=(phi_i,phi_j) + * is the mass matrix approximated by + * lhs_quadrature, and + * Q is the matrix + * Qiq=phii(xq) + * wq where + * wq are quadrature + * weights; U is the vector of + * quadrature point data + * uq. + * + * In order to then get the nodal + * representation Vh + * of the projection of U, one + * computes Vh = X U, + * X=M-1 Q. The purpose + * of this function is to compute the + * matrix X and return it + * through the last argument of this + * function. + * + * Note that this function presently only + * supports scalar data. An extension of + * the mass matrix is of course trivial, + * but one has to define the order of + * data in the vector U if it + * contains vector valued data in all + * quadrature points. + * + * A use for this function is described + * in the introduction to the step-18 + * example program. + * + * The opposite of this function, + * interpolation of a finite element + * function onto quadrature points is + * essentially what the + * FEValues::get_function_values + * functions do; to make things a little + * simpler, the + * FETools::compute_interpolation_to_quadrature_points_matrix + * provides the matrix form of this. + * + * Note that this function works + * on a single cell, rather than + * an entire triangulation. In + * effect, it therefore doesn't + * matter if you use a continuous + * or discontinuous version of + * the finite element. + * + * It is worth noting that there + * are a few confusing cases of + * this function. The first one + * is that it really only makes + * sense to project onto a finite + * element that has at most as + * many degrees of freedom per + * cell as there are quadrature + * points; the projection of N + * quadrature point data into a + * space with M>N unknowns is + * well-defined, but often yields + * funny and non-intuitive + * results. Secondly, one would + * think that if the quadrature + * point data is defined in the + * support points of the finite + * element, i.e. the quadrature + * points of + * ths_quadrature equal + * fe.get_unit_support_points(), + * then the projection should be + * the identity, i.e. each degree + * of freedom of the finite + * element equals the value of + * the given data in the support + * point of the corresponding + * shape function. However, this + * is not generally the case: + * while the matrix Q in + * that case is the identity + * matrix, the mass matrix + * M is not equal to the + * identity matrix, except for + * the special case that the + * quadrature formula + * lhs_quadrature also + * has its quadrature points in + * the support points of the + * finite element. + * + * Finally, this function only defines a + * cell wise projection, while one + * frequently wants to apply it to all + * cells in a triangulation. However, if + * it is applied to one cell after the + * other, the results from later cells + * may overwrite nodal values computed + * already from previous cells if degrees + * of freedom live on the interfaces + * between cells. The function is + * therefore most useful for + * discontinuous elements. + */ + template + void + compute_projection_from_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature, + FullMatrix &X); + + /** + * Given a (scalar) local finite element + * function, compute the matrix that maps + * the vector of nodal values onto the + * vector of values of this function at + * quadrature points as given by the + * second argument. In a sense, this + * function does the opposite of the @p + * compute_projection_from_quadrature_points_matrix + * function. + */ + template + void + compute_interpolation_to_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &quadrature, + FullMatrix &I_q); + + /** + * Computes the projection of tensorial + * (first-order tensor) + * data stored at the quadrature points + * @p vector_of_tensors_at_qp + * to data @p vector_of_tensors_at_nodes + * at the support points of the cell. + * The data in + * @p vector_of_tensors_at_qp + * is ordered sequentially following the + * quadrature point numbering. + * The size of + * @p vector_of_tensors_at_qp + * must correspond to the number of columns + * of @p projection_matrix. + * The size of @p vector_of_tensors_at_nodes + * must correspond to the number of rows of + * @p vector_of_tensors_at_nodes . + * The projection matrix + * @p projection_matrix desribes the + * projection of scalar data from the + * quadrature points and can be obtained + * from the + * FETools::compute_projection_from_quadrature_points_matrix + * function. + */ + template + void + compute_projection_from_quadrature_points( + const FullMatrix &projection_matrix, + const std::vector< Tensor<1, dim > > &vector_of_tensors_at_qp, + std::vector< Tensor<1, dim > > &vector_of_tensors_at_nodes); + + + + /** + * same as last function but for a + * @p SymmetricTensor . + */ + template + void + compute_projection_from_quadrature_points( + const FullMatrix &projection_matrix, + const std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_qp, + std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_nodes); + + + + + /** + * This method implements the + * FETools::compute_projection_from_quadrature_points_matrix + * method for faces of a mesh. + * The matrix that it returns, X, is face specific + * and its size is fe.dofs_per_cell by + * rhs_quadrature.size(). + * The dimension, dim must be larger than 1 for this class, + * since Quadrature objects are required. See the + * documentation on the Quadrature class for more information. + */ + template + void + compute_projection_from_face_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature, + const typename DoFHandler::active_cell_iterator & cell, + unsigned int face, + FullMatrix &X); + + + + //@} + /** + * @name Functions which should be in DoFTools + */ + //@{ + /** + * Gives the interpolation of a the + * @p dof1-function @p u1 to a + * @p dof2-function @p u2. @p dof1 and + * @p dof2 need to be DoFHandlers + * based on the same triangulation. + * + * If the elements @p fe1 and @p fe2 + * are either both continuous or + * both discontinuous then this + * interpolation is the usual point + * interpolation. The same is true + * if @p fe1 is a continuous and + * @p fe2 is a discontinuous finite + * element. For the case that @p fe1 + * is a discontinuous and @p fe2 is + * a continuous finite element + * there is no point interpolation + * defined at the discontinuities. + * Therefore the meanvalue is taken + * at the DoF values on the + * discontinuities. + * + * Note that for continuous + * elements on grids with hanging + * nodes (i.e. locally refined + * grids) this function does not + * give the expected output. + * Indeed, the resulting output + * vector does not necessarily + * respect continuity + * requirements at hanging nodes: + * if, for example, you are + * interpolating a Q2 field to a + * Q1 field, then at hanging + * nodes the output field will + * have the function value of the + * input field, which however is + * not usually the mean value of + * the two adjacent nodes. It is + * thus not part of the Q1 + * function space on the whole + * triangulation, although it is + * of course Q1 on each cell. + * + * For this case (continuous + * elements on grids with hanging + * nodes), please use the + * @p interpolate function with + * an additional + * @p ConstraintMatrix argument, + * see below, or make the field + * conforming yourself by calling + * the @p distribute function of + * your hanging node constraints + * object. + */ + template class DH1, + template class DH2, + class InVector, class OutVector> + void + interpolate (const DH1 &dof1, + const InVector &u1, + const DH2 &dof2, + OutVector &u2); + + /** + * Gives the interpolation of a + * the @p dof1-function @p u1 to + * a @p dof2-function @p u2. @p + * dof1 and @p dof2 need to be + * DoFHandlers (or + * hp::DoFHandlers) based on the + * same triangulation. @p + * constraints is a hanging node + * constraints object + * corresponding to @p dof2. This + * object is particular important + * when interpolating onto + * continuous elements on grids + * with hanging nodes (locally + * refined grids). + * + * If the elements @p fe1 and @p fe2 + * are either both continuous or + * both discontinuous then this + * interpolation is the usual point + * interpolation. The same is true + * if @p fe1 is a continuous and + * @p fe2 is a discontinuous finite + * element. For the case that @p fe1 + * is a discontinuous and @p fe2 is + * a continuous finite element + * there is no point interpolation + * defined at the discontinuities. + * Therefore the meanvalue is taken + * at the DoF values on the + * discontinuities. + */ + template class DH1, + template class DH2, + class InVector, class OutVector> + void interpolate (const DH1 &dof1, + const InVector &u1, + const DH2 &dof2, + const ConstraintMatrix &constraints, + OutVector& u2); + + /** + * Gives the interpolation of the + * @p fe1-function @p u1 to a + * @p fe2-function, and + * interpolates this to a second + * @p fe1-function named + * @p u1_interpolated. + * + * Note, that this function does + * not work on continuous + * elements at hanging nodes. For + * that case use the + * @p back_interpolate function, + * below, that takes an + * additional + * @p ConstraintMatrix object. + * + * Furthermore note, that for the + * specific case when the finite + * element space corresponding to + * @p fe1 is a subset of the + * finite element space + * corresponding to @p fe2, this + * function is simply an identity + * mapping. + */ + template + void back_interpolate (const DoFHandler &dof1, + const InVector &u1, + const FiniteElement &fe2, + OutVector &u1_interpolated); + + /** + * Same as last function, except + * that the dof handler objects + * might be of type + * @p hp::DoFHandler. + */ + template class DH, + class InVector, class OutVector, int spacedim> + void back_interpolate (const DH &dof1, + const InVector &u1, + const FiniteElement &fe2, + OutVector &u1_interpolated); + + /** + * Gives the interpolation of the + * @p dof1-function @p u1 to a + * @p dof2-function, and + * interpolates this to a second + * @p dof1-function named + * @p u1_interpolated. + * @p constraints1 and + * @p constraints2 are the + * hanging node constraints + * corresponding to @p dof1 and + * @p dof2, respectively. These + * objects are particular + * important when continuous + * elements on grids with hanging + * nodes (locally refined grids) + * are involved. + * + * Furthermore note, that for the + * specific case when the finite + * element space corresponding to + * @p dof1 is a subset of the + * finite element space + * corresponding to @p dof2, this + * function is simply an identity + * mapping. + */ + template + void back_interpolate (const DoFHandler& dof1, + const ConstraintMatrix& constraints1, + const InVector& u1, + const DoFHandler& dof2, + const ConstraintMatrix& constraints2, + OutVector& u1_interpolated); + + /** + * Gives $(Id-I_h)z_1$ for a given + * @p dof1-function $z_1$, where $I_h$ + * is the interpolation from @p fe1 + * to @p fe2. The result $(Id-I_h)z_1$ is + * written into @p z1_difference. + * + * Note, that this function does + * not work for continuous + * elements at hanging nodes. For + * that case use the + * @p interpolation_difference + * function, below, that takes an + * additional + * @p ConstraintMatrix object. + */ + template + void interpolation_difference(const DoFHandler &dof1, + const InVector &z1, + const FiniteElement &fe2, + OutVector &z1_difference); + + /** + * Gives $(Id-I_h)z_1$ for a given + * @p dof1-function $z_1$, where $I_h$ + * is the interpolation from @p fe1 + * to @p fe2. The result $(Id-I_h)z_1$ is + * written into @p z1_difference. + * @p constraints1 and + * @p constraints2 are the + * hanging node constraints + * corresponding to @p dof1 and + * @p dof2, respectively. These + * objects are particular + * important when continuous + * elements on grids with hanging + * nodes (locally refined grids) + * are involved. + */ + template + void interpolation_difference(const DoFHandler& dof1, + const ConstraintMatrix& constraints1, + const InVector& z1, + const DoFHandler& dof2, + const ConstraintMatrix& constraints2, + OutVector& z1_difference); + + /** + * $L^2$ projection for + * discontinuous + * elements. Operates the same + * direction as interpolate. + * + * The global projection can be + * computed by local matrices if + * the finite element spaces are + * discontinuous. With continuous + * elements, this is impossible, + * since a global mass matrix + * must be inverted. + */ + template + void project_dg (const DoFHandler& dof1, + const InVector& u1, + const DoFHandler& dof2, + OutVector& u2); + + /** + * Gives the patchwise + * extrapolation of a @p dof1 + * function @p z1 to a @p dof2 + * function @p z2. @p dof1 and + * @p dof2 need to be DoFHandler + * based on the same triangulation. + * + * This function is interesting + * for e.g. extrapolating + * patchwise a piecewise linear + * solution to a piecewise + * quadratic solution. + * + * Note that the resulting field + * does not satisfy continuity + * requirements of the given + * finite elements. + * + * When you use continuous + * elements on grids with hanging + * nodes, please use the + * @p extrapolate function with + * an additional + * ConstraintMatrix argument, + * see below. + * + * Since this function operates + * on patches of cells, it is + * required that the underlying + * grid is refined at least once + * for every coarse grid cell. If + * this is not the case, an + * exception will be raised. + */ + template + void extrapolate (const DoFHandler& dof1, + const InVector& z1, + const DoFHandler& dof2, + OutVector& z2); + + /** + * Gives the patchwise + * extrapolation of a @p dof1 + * function @p z1 to a @p dof2 + * function @p z2. @p dof1 and + * @p dof2 need to be DoFHandler + * based on the same triangulation. + * @p constraints is a hanging + * node constraints object + * corresponding to + * @p dof2. This object is + * particular important when + * interpolating onto continuous + * elements on grids with hanging + * nodes (locally refined grids). + * + * Otherwise, the same holds as + * for the other @p extrapolate + * function. + */ + template + void extrapolate (const DoFHandler& dof1, + const InVector& z1, + const DoFHandler& dof2, + const ConstraintMatrix& constraints, + OutVector& z2); + //@} + /** + * The numbering of the degrees + * of freedom in continous finite + * elements is hierarchic, + * i.e. in such a way that we + * first number the vertex dofs, + * in the order of the vertices + * as defined by the + * triangulation, then the line + * dofs in the order and + * respecting the direction of + * the lines, then the dofs on + * quads, etc. However, we could + * have, as well, numbered them + * in a lexicographic way, + * i.e. with indices first + * running in x-direction, then + * in y-direction and finally in + * z-direction. Discontinuous + * elements of class FE_DGQ() + * are numbered in this way, for + * example. + * + * This function constructs a + * table which lexicographic + * index each degree of freedom + * in the hierarchic numbering + * would have. It operates on the + * continuous finite element + * given as first argument, and + * outputs the lexicographic + * indices in the second. + * + * Note that since this function + * uses specifics of the + * continuous finite elements, it + * can only operate on + * FiniteElementData objects + * inherent in FE_Q(). However, + * this function does not take a + * FE_Q object as it is also + * invoked by the FE_Q() + * constructor. + * + * It is assumed that the size of + * the output argument already + * matches the correct size, + * which is equal to the number + * of degrees of freedom in the + * finite element. + */ + template + void + hierarchic_to_lexicographic_numbering (const FiniteElementData &fe_data, + std::vector &h2l); + + /** + * Like the previous function but + * instead of returning its + * result through the last + * argument return it as a value. + */ + template + std::vector + hierarchic_to_lexicographic_numbering (const FiniteElementData &fe_data); + + /** + * This is the reverse function + * to the above one, generating + * the map from the lexicographic + * to the hierarchical + * numbering. All the remarks + * made about the above function + * are also valid here. + */ + template + void + lexicographic_to_hierarchic_numbering (const FiniteElementData &fe_data, + std::vector &l2h); + + /** + * Like the previous function but + * instead of returning its + * result through the last + * argument return it as a value. + */ + template + std::vector + lexicographic_to_hierarchic_numbering (const FiniteElementData &fe_data); + + /** + * Parse the name of a finite + * element and generate a finite + * element object accordingly. + * + * The name must be in the form which + * is returned by the + * FiniteElement::get_name + * function, where a few + * modifications are allowed: + * + *
  • Dimension template + * parameters <2> etc. can + * be omitted. Alternatively, the + * explicit number can be + * replaced by dim or + * d. If a number is + * given, it must match + * the template parameter of this + * function. + * + *
  • The powers used for + * FESystem may either be numbers + * or can be + * replaced by dim or + * d. + *
+ * + * If no finite element can be + * reconstructed from this + * string, an exception of type + * @p FETools::ExcInvalidFEName + * is thrown. + * + * The function returns a pointer + * to a newly create finite + * element. It is in the caller's + * responsibility to destroy the + * object pointed to at an + * appropriate later time. + * + * Since the value of the template + * argument can't be deduced from the + * (string) argument given to this + * function, you have to explicitly + * specify it when you call this + * function. + * + * This function knows about all + * the standard elements defined + * in the library. However, it + * doesn't by default know about + * elements that you may have + * defined in your program. To + * make your own elements known + * to this function, use the + * add_fe_name() function. + * This function does not work + * if one wants to get a codimension + * 1 finite element. + */ + template + FiniteElement * + get_fe_from_name (const std::string &name); + + + /** + * Extend the list of finite + * elements that can be generated + * by get_fe_from_name() by the + * one given as @p name. If + * get_fe_from_name() is later + * called with this name, it will + * use the object given as second + * argument to create a finite + * element object. + * + * The format of the @p name + * parameter should include the + * name of a finite + * element. However, it is safe + * to use either the class name + * alone or to use the result of + * FiniteElement::get_name (which + * includes the space dimension + * as well as the polynomial + * degree), since everything + * after the first non-name + * character will be ignored. + * + * The FEFactory object should be + * an object newly created with + * new. FETools will + * take ownership of this object + * and delete it once it is not + * used anymore. + * + * In most cases, if you want + * objects of type + * MyFE be created + * whenever the name + * my_fe is given to + * get_fe_from_name, you will + * want the second argument to + * this function be of type + * FEFactory@, but you can + * of course create your custom + * finite element factory class. + * + * This function takes over + * ownership of the object given + * as second argument, i.e. you + * should never attempt to + * destroy it later on. The + * object will be deleted at the + * end of the program's lifetime. + * + * If the name of the element + * is already in use, an exception + * is thrown. Thus, functionality + * of get_fe_from_name() can only + * be added, not changed. + * + * @note This function + * manipulates a global table + * (one table for each space + * dimension). It is thread safe + * in the sense that every access + * to this table is secured by a + * lock. Nevertheless, since each + * name can be added only once, + * user code has to make sure + * that only one thread adds a + * new element. + * + * Note also that this table + * exists once for each space + * dimension. If you have a + * program that works with finite + * elements in different space + * dimensions (for example, @ref + * step_4 "step-4" does something + * like this), then you should + * call this function for each + * space dimension for which you + * want your finite element added + * to the map. + */ + template + void add_fe_name (const std::string& name, + const FEFactoryBase* factory); + + /** + * The string used for + * get_fe_from_name() cannot be + * translated to a finite + * element. + * + * Either the string is badly + * formatted or you are using a + * custom element that must be + * added using add_fe_name() + * first. + * + * @ingroup Exceptions + */ + DeclException1 (ExcInvalidFEName, + std::string, + << "Can't re-generate a finite element from the string '" + << arg1 << "'."); + + /** + * The string used for + * get_fe_from_name() cannot be + * translated to a finite + * element. + * + * Dimension arguments in finite + * element names should be + * avoided. If they are there, + * the dimension should be + * dim or + * d. Here, you gave a + * numeric dimension argument, + * which does not match the + * template dimension of the + * finite element class. + * + * @ingroup Exceptions + */ + DeclException2 (ExcInvalidFEDimension, + char, int, + << "The dimension " << arg1 + << " in the finite element string must match " + << "the space dimension " + << arg2 << "."); + + /** + * Exception + * + * @ingroup Exceptions + */ + DeclException0 (ExcInvalidFE); + + /** + * The finite element must be + * @ref GlossPrimitive "primitive". + * + * @ingroup Exceptions + */ + DeclException0 (ExcFENotPrimitive); + /** + * Exception + * + * @ingroup Exceptions + */ + DeclException0 (ExcTriangulationMismatch); + + /** + * A continuous element is used + * on a mesh with hanging nodes, + * but the constraint matrices + * are missing. + * + * @ingroup Exceptions + */ + DeclException1 (ExcHangingNodesNotAllowed, + int, + << "You are using continuous elements on a grid with " + << "hanging nodes but without providing hanging node " + << "constraints. Use the respective function with " + << "additional ConstraintMatrix argument(s), instead."); + /** + * You need at least two grid levels. + * + * @ingroup Exceptions + */ + DeclException0 (ExcGridNotRefinedAtLeastOnce); + /** + * The dimensions of the matrix + * used did not match the + * expected dimensions. + * + * @ingroup Exceptions + */ + DeclException4 (ExcMatrixDimensionMismatch, + int, int, int, int, + << "This is a " << arg1 << "x" << arg2 << " matrix, " + << "but should be a " << arg3 << "x" << arg4 << " matrix."); + + /** + * Exception thrown if an + * embedding matrix was computed + * inaccurately. + * + * @ingroup Exceptions + */ + DeclException1(ExcLeastSquaresError, double, + << "Least squares fit leaves a gap of " << arg1); + + /** + * Exception thrown if one variable + * may not be greater than another. + * + * @ingroup Exceptions + */ + DeclException2 (ExcNotGreaterThan, + int, int, + << arg1 << " must be greater than " << arg2); } +#ifndef DOXYGEN + +namespace FETools +{ + template + FiniteElement* + FEFactory::get (const unsigned int degree) const + { + return new FE(degree); + } +} + +#endif /*@}*/ diff --git a/deal.II/deal.II/source/fe/fe_tools.all_dimensions.cc b/deal.II/deal.II/source/fe/fe_tools.all_dimensions.cc index 18cceafaa9..4c40c955fb 100644 --- a/deal.II/deal.II/source/fe/fe_tools.all_dimensions.cc +++ b/deal.II/deal.II/source/fe/fe_tools.all_dimensions.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2005, 2006, 2008 by the deal.II authors +// Copyright (C) 2005, 2006, 2008, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -18,280 +18,285 @@ DEAL_II_NAMESPACE_OPEN - -template -void -FETools::hierarchic_to_lexicographic_numbering (const FiniteElementData &fe, - std::vector &h2l) +namespace FETools { - Assert (h2l.size() == fe.dofs_per_cell, - ExcDimensionMismatch (h2l.size(), fe.dofs_per_cell)); - h2l = hierarchic_to_lexicographic_numbering (fe); -} - - -template -std::vector -FETools::hierarchic_to_lexicographic_numbering (const FiniteElementData &fe) -{ - Assert (fe.n_components() == 1, ExcInvalidFE()); - - std::vector h2l (fe.dofs_per_cell); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - // polynomial degree - const unsigned int degree = fe.dofs_per_line+1; - // number of grid points in each - // direction - const unsigned int n = degree+1; - - // the following lines of code are - // somewhat odd, due to the way the - // hierarchic numbering is - // organized. if someone would - // really want to understand these - // lines, you better draw some - // pictures where you indicate the - // indices and orders of vertices, - // lines, etc, along with the - // numbers of the degrees of - // freedom in hierarchical and - // lexicographical order - switch (dim) - { - case 1: + template + void + hierarchic_to_lexicographic_numbering (const FiniteElementData &fe, + std::vector &h2l) + { + Assert (h2l.size() == fe.dofs_per_cell, + ExcDimensionMismatch (h2l.size(), fe.dofs_per_cell)); + h2l = hierarchic_to_lexicographic_numbering (fe); + } + + + + template + std::vector + hierarchic_to_lexicographic_numbering (const FiniteElementData &fe) + { + Assert (fe.n_components() == 1, ExcInvalidFE()); + + std::vector h2l (fe.dofs_per_cell); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + // polynomial degree + const unsigned int degree = fe.dofs_per_line+1; + // number of grid points in each + // direction + const unsigned int n = degree+1; + + // the following lines of code are + // somewhat odd, due to the way the + // hierarchic numbering is + // organized. if someone would + // really want to understand these + // lines, you better draw some + // pictures where you indicate the + // indices and orders of vertices, + // lines, etc, along with the + // numbers of the degrees of + // freedom in hierarchical and + // lexicographical order + switch (dim) { - h2l[0] = 0; - h2l[1] = dofs_per_cell-1; - for (unsigned int i=2; i + void + lexicographic_to_hierarchic_numbering (const FiniteElementData &fe, + std::vector &l2h) + { + l2h = lexicographic_to_hierarchic_numbering (fe); + } -template -void -FETools::lexicographic_to_hierarchic_numbering (const FiniteElementData &fe, - std::vector &l2h) -{ - l2h = lexicographic_to_hierarchic_numbering (fe); -} + template + std::vector + lexicographic_to_hierarchic_numbering (const FiniteElementData &fe) + { + return Utilities::invert_permutation(hierarchic_to_lexicographic_numbering (fe)); + } -template -std::vector -FETools::lexicographic_to_hierarchic_numbering (const FiniteElementData &fe) -{ - return Utilities::invert_permutation(hierarchic_to_lexicographic_numbering (fe)); } - - -template -void -FETools::hierarchic_to_lexicographic_numbering<1> -(const FiniteElementData<1> &fe, - std::vector &h2l); -template -void -FETools::hierarchic_to_lexicographic_numbering<2> -(const FiniteElementData<2> &fe, - std::vector &h2l); -template -void -FETools::hierarchic_to_lexicographic_numbering<3> -(const FiniteElementData<3> &fe, - std::vector &h2l); - - -template -void -FETools::lexicographic_to_hierarchic_numbering<1> -(const FiniteElementData<1> &fe, - std::vector &l2h); -template -void -FETools::lexicographic_to_hierarchic_numbering<2> -(const FiniteElementData<2> &fe, - std::vector &l2h); -template -void -FETools::lexicographic_to_hierarchic_numbering<3> -(const FiniteElementData<3> &fe, - std::vector &l2h); - - - -template -std::vector -FETools::hierarchic_to_lexicographic_numbering<1> -(const FiniteElementData<1> &fe); -template -std::vector -FETools::hierarchic_to_lexicographic_numbering<2> -(const FiniteElementData<2> &fe); -template -std::vector -FETools::hierarchic_to_lexicographic_numbering<3> -(const FiniteElementData<3> &fe); - - -template -std::vector -FETools::lexicographic_to_hierarchic_numbering<1> -(const FiniteElementData<1> &fe); -template -std::vector -FETools::lexicographic_to_hierarchic_numbering<2> -(const FiniteElementData<2> &fe); -template -std::vector -FETools::lexicographic_to_hierarchic_numbering<3> -(const FiniteElementData<3> &fe); - +// explicit instantiations +namespace FETools +{ + template + void + hierarchic_to_lexicographic_numbering<1> + (const FiniteElementData<1> &fe, + std::vector &h2l); + template + void + hierarchic_to_lexicographic_numbering<2> + (const FiniteElementData<2> &fe, + std::vector &h2l); + template + void + hierarchic_to_lexicographic_numbering<3> + (const FiniteElementData<3> &fe, + std::vector &h2l); + + + template + void + lexicographic_to_hierarchic_numbering<1> + (const FiniteElementData<1> &fe, + std::vector &l2h); + template + void + lexicographic_to_hierarchic_numbering<2> + (const FiniteElementData<2> &fe, + std::vector &l2h); + template + void + lexicographic_to_hierarchic_numbering<3> + (const FiniteElementData<3> &fe, + std::vector &l2h); + + + + template + std::vector + hierarchic_to_lexicographic_numbering<1> + (const FiniteElementData<1> &fe); + template + std::vector + hierarchic_to_lexicographic_numbering<2> + (const FiniteElementData<2> &fe); + template + std::vector + hierarchic_to_lexicographic_numbering<3> + (const FiniteElementData<3> &fe); + + + template + std::vector + lexicographic_to_hierarchic_numbering<1> + (const FiniteElementData<1> &fe); + template + std::vector + lexicographic_to_hierarchic_numbering<2> + (const FiniteElementData<2> &fe); + template + std::vector + lexicographic_to_hierarchic_numbering<3> + (const FiniteElementData<3> &fe); +} DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index 76082357d0..fb7242f7fb 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -238,397 +238,399 @@ namespace } -template -FETools::FEFactoryBase::~FEFactoryBase() -{} +namespace FETools +{ + template + FEFactoryBase::~FEFactoryBase() + {} -template -void FETools::compute_component_wise( - const FiniteElement& element, - std::vector& renumbering, - std::vector >& comp_start) -{ - Assert(renumbering.size() == element.dofs_per_cell, - ExcDimensionMismatch(renumbering.size(), - element.dofs_per_cell)); + template + void compute_component_wise( + const FiniteElement& element, + std::vector& renumbering, + std::vector >& comp_start) + { + Assert(renumbering.size() == element.dofs_per_cell, + ExcDimensionMismatch(renumbering.size(), + element.dofs_per_cell)); - comp_start.resize(element.n_base_elements()); + comp_start.resize(element.n_base_elements()); - unsigned int k=0; - for (unsigned int i=0;i, unsigned int> + indices = element.system_to_base_index(i); + renumbering[i] = comp_start[indices.first.first][indices.first.second] + +indices.second; + } + } - for (unsigned int j=0;j + void compute_block_renumbering ( + const FiniteElement& element, + std::vector& renumbering, + std::vector& block_data, + bool return_start_indices) + { + Assert(renumbering.size() == element.dofs_per_cell, + ExcDimensionMismatch(renumbering.size(), + element.dofs_per_cell)); + Assert(block_data.size() == element.n_blocks(), + ExcDimensionMismatch(block_data.size(), + element.n_blocks())); + + unsigned int k=0; + unsigned int i=0; + for (unsigned int b=0;b, unsigned int> - indices = element.system_to_base_index(i); - renumbering[i] = comp_start[indices.first.first][indices.first.second] - +indices.second; - } -} + std::vector start_indices(block_data.size()); + k = 0; + for (unsigned int i=0;i + indices = element.system_to_block_index(i); + renumbering[i] = start_indices[indices.first] + +indices.second; + } + } -template -void FETools::compute_block_renumbering ( - const FiniteElement& element, - std::vector& renumbering, - std::vector& block_data, - bool return_start_indices) -{ - Assert(renumbering.size() == element.dofs_per_cell, - ExcDimensionMismatch(renumbering.size(), - element.dofs_per_cell)); - Assert(block_data.size() == element.n_blocks(), - ExcDimensionMismatch(block_data.size(), - element.n_blocks())); - - unsigned int k=0; - unsigned int i=0; - for (unsigned int b=0;b + void get_interpolation_matrix (const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix) + { + Assert (fe1.n_components() == fe2.n_components(), + ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); + Assert(interpolation_matrix.m()==fe2.dofs_per_cell && + interpolation_matrix.n()==fe1.dofs_per_cell, + ExcMatrixDimensionMismatch(interpolation_matrix.m(), + interpolation_matrix.n(), + fe2.dofs_per_cell, + fe1.dofs_per_cell)); + + // first try the easy way: maybe + // the FE wants to implement things + // itself: + bool fe_implements_interpolation = true; + try { - block_data[i++] = (return_start_indices) - ? k - : (element.base_element(b).n_dofs_per_cell()); - k += element.base_element(b).n_dofs_per_cell(); + gim_forwarder (fe1, fe2, interpolation_matrix); } - Assert (i == element.n_blocks(), ExcInternalError()); - - std::vector start_indices(block_data.size()); - k = 0; - for (unsigned int i=0;i::ExcInterpolationNotImplemented &) { - start_indices[i] = k; - k += block_data[i]; + // too bad.... + fe_implements_interpolation = false; } - -//TODO:[GK] This does not work for a single RT - for (unsigned int i=0;i - indices = element.system_to_block_index(i); - renumbering[i] = start_indices[indices.first] - +indices.second; - } -} + if (fe_implements_interpolation == true) + return; + + // uh, so this was not the + // case. hm. then do it the hard + // way. note that this will only + // work if the element is + // primitive, so check this first + Assert (fe1.is_primitive() == true, ExcFENotPrimitive()); + Assert (fe2.is_primitive() == true, ExcFENotPrimitive()); + + // Initialize FEValues for fe1 at + // the unit support points of the + // fe2 element. + const std::vector > & + fe2_support_points = fe2.get_unit_support_points (); + + typedef FiniteElement FEL; + Assert(fe2_support_points.size()==fe2.dofs_per_cell, + typename FEL::ExcFEHasNoSupportPoints()); + + for (unsigned int i=0; i -void FETools::get_interpolation_matrix (const FiniteElement &fe1, - const FiniteElement &fe2, - FullMatrix &interpolation_matrix) -{ - Assert (fe1.n_components() == fe2.n_components(), - ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); - Assert(interpolation_matrix.m()==fe2.dofs_per_cell && - interpolation_matrix.n()==fe1.dofs_per_cell, - ExcMatrixDimensionMismatch(interpolation_matrix.m(), - interpolation_matrix.n(), - fe2.dofs_per_cell, - fe1.dofs_per_cell)); - - // first try the easy way: maybe - // the FE wants to implement things - // itself: - bool fe_implements_interpolation = true; - try - { - gim_forwarder (fe1, fe2, interpolation_matrix); - } - catch (typename FiniteElement::ExcInterpolationNotImplemented &) - { - // too bad.... - fe_implements_interpolation = false; - } - if (fe_implements_interpolation == true) - return; - - // uh, so this was not the - // case. hm. then do it the hard - // way. note that this will only - // work if the element is - // primitive, so check this first - Assert (fe1.is_primitive() == true, ExcFENotPrimitive()); - Assert (fe2.is_primitive() == true, ExcFENotPrimitive()); - - // Initialize FEValues for fe1 at - // the unit support points of the - // fe2 element. - const std::vector > & - fe2_support_points = fe2.get_unit_support_points (); - - typedef FiniteElement FEL; - Assert(fe2_support_points.size()==fe2.dofs_per_cell, - typename FEL::ExcFEHasNoSupportPoints()); - - for (unsigned int i=0; i + void get_back_interpolation_matrix(const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix) + { + Assert (fe1.n_components() == fe2.n_components(), + ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); + Assert(interpolation_matrix.m()==fe1.dofs_per_cell && + interpolation_matrix.n()==fe1.dofs_per_cell, + ExcMatrixDimensionMismatch(interpolation_matrix.m(), + interpolation_matrix.n(), + fe1.dofs_per_cell, + fe1.dofs_per_cell)); + + FullMatrix first_matrix (fe2.dofs_per_cell, fe1.dofs_per_cell); + FullMatrix second_matrix(fe1.dofs_per_cell, fe2.dofs_per_cell); + + get_interpolation_matrix(fe1, fe2, first_matrix); + get_interpolation_matrix(fe2, fe1, second_matrix); + + // int_matrix=second_matrix*first_matrix + second_matrix.mmult(interpolation_matrix, first_matrix); + } -template -void FETools::get_back_interpolation_matrix(const FiniteElement &fe1, + template + void get_interpolation_difference_matrix (const FiniteElement &fe1, const FiniteElement &fe2, - FullMatrix &interpolation_matrix) -{ - Assert (fe1.n_components() == fe2.n_components(), - ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); - Assert(interpolation_matrix.m()==fe1.dofs_per_cell && - interpolation_matrix.n()==fe1.dofs_per_cell, - ExcMatrixDimensionMismatch(interpolation_matrix.m(), - interpolation_matrix.n(), - fe1.dofs_per_cell, - fe1.dofs_per_cell)); - - FullMatrix first_matrix (fe2.dofs_per_cell, fe1.dofs_per_cell); - FullMatrix second_matrix(fe1.dofs_per_cell, fe2.dofs_per_cell); - - get_interpolation_matrix(fe1, fe2, first_matrix); - get_interpolation_matrix(fe2, fe1, second_matrix); - - // int_matrix=second_matrix*first_matrix - second_matrix.mmult(interpolation_matrix, first_matrix); -} + FullMatrix &difference_matrix) + { + Assert (fe1.n_components() == fe2.n_components(), + ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); + Assert(difference_matrix.m()==fe1.dofs_per_cell && + difference_matrix.n()==fe1.dofs_per_cell, + ExcMatrixDimensionMismatch(difference_matrix.m(), + difference_matrix.n(), + fe1.dofs_per_cell, + fe1.dofs_per_cell)); + + FullMatrix interpolation_matrix(fe1.dofs_per_cell); + get_back_interpolation_matrix(fe1, fe2, interpolation_matrix); + + for (unsigned int i=0; i -void FETools::get_interpolation_difference_matrix (const FiniteElement &fe1, - const FiniteElement &fe2, - FullMatrix &difference_matrix) -{ - Assert (fe1.n_components() == fe2.n_components(), - ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); - Assert(difference_matrix.m()==fe1.dofs_per_cell && - difference_matrix.n()==fe1.dofs_per_cell, - ExcMatrixDimensionMismatch(difference_matrix.m(), - difference_matrix.n(), - fe1.dofs_per_cell, - fe1.dofs_per_cell)); - - FullMatrix interpolation_matrix(fe1.dofs_per_cell); - get_back_interpolation_matrix(fe1, fe2, interpolation_matrix); - - for (unsigned int i=0; i + void get_projection_matrix (const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &matrix) + { + Assert (fe1.n_components() == 1, ExcNotImplemented()); + Assert (fe1.n_components() == fe2.n_components(), + ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); + Assert(matrix.m()==fe2.dofs_per_cell && matrix.n()==fe1.dofs_per_cell, + ExcMatrixDimensionMismatch(matrix.m(), matrix.n(), + fe2.dofs_per_cell, + fe1.dofs_per_cell)); + matrix = 0; + + unsigned int n1 = fe1.dofs_per_cell; + unsigned int n2 = fe2.dofs_per_cell; + + // First, create a local mass matrix for + // the unit cell + Triangulation tr; + GridGenerator::hyper_cube(tr); + + // Choose a quadrature rule + // Gauss is exact up to degree 2n-1 + const unsigned int degree = std::max(fe1.tensor_degree(), fe2.tensor_degree()); + Assert (degree != numbers::invalid_unsigned_int, + ExcNotImplemented()); + + QGauss quadrature(degree+1); + // Set up FEValues. + const UpdateFlags flags = update_values | update_quadrature_points | update_JxW_values; + FEValues val1 (fe1, quadrature, update_values); + val1.reinit (tr.begin_active()); + FEValues val2 (fe2, quadrature, flags); + val2.reinit (tr.begin_active()); + + // Integrate and invert mass matrix + // This happens in the target space + FullMatrix mass (n2, n2); + + for (unsigned int k=0;k b(n2); + Vector x(n2); + for (unsigned int j=0;j -void FETools::get_projection_matrix (const FiniteElement &fe1, - const FiniteElement &fe2, - FullMatrix &matrix) -{ - Assert (fe1.n_components() == 1, ExcNotImplemented()); - Assert (fe1.n_components() == fe2.n_components(), - ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); - Assert(matrix.m()==fe2.dofs_per_cell && matrix.n()==fe1.dofs_per_cell, - ExcMatrixDimensionMismatch(matrix.m(), matrix.n(), - fe2.dofs_per_cell, - fe1.dofs_per_cell)); - matrix = 0; - - unsigned int n1 = fe1.dofs_per_cell; - unsigned int n2 = fe2.dofs_per_cell; - - // First, create a local mass matrix for - // the unit cell - Triangulation tr; - GridGenerator::hyper_cube(tr); - - // Choose a quadrature rule - // Gauss is exact up to degree 2n-1 - const unsigned int degree = std::max(fe1.tensor_degree(), fe2.tensor_degree()); - Assert (degree != numbers::invalid_unsigned_int, - ExcNotImplemented()); - - QGauss quadrature(degree+1); - // Set up FEValues. - const UpdateFlags flags = update_values | update_quadrature_points | update_JxW_values; - FEValues val1 (fe1, quadrature, update_values); - val1.reinit (tr.begin_active()); - FEValues val2 (fe2, quadrature, flags); - val2.reinit (tr.begin_active()); - - // Integrate and invert mass matrix - // This happens in the target space - FullMatrix mass (n2, n2); - - for (unsigned int k=0;k b(n2); - Vector x(n2); - - for (unsigned int j=0;j -void -FETools::compute_node_matrix( - FullMatrix& N, - const FiniteElement& fe) -{ - const unsigned int n_dofs = fe.dofs_per_cell; - Assert (fe.has_generalized_support_points(), ExcNotInitialized()); - Assert (N.n()==n_dofs, ExcDimensionMismatch(N.n(), n_dofs)); - Assert (N.m()==n_dofs, ExcDimensionMismatch(N.m(), n_dofs)); - - const std::vector >& points = fe.get_generalized_support_points(); - - // We need the values of the - // polynomials in all generalized - // support points. - std::vector > - values (dim, std::vector(points.size())); - - // In this vector, we store the - // result of the interpolation - std::vector local_dofs(n_dofs); - - // One row per shape - // function. Remember that these - // are the 'raw' shape functions - // where the inverse node matrix is - // empty. Otherwise, this would - // yield identity. - for (unsigned int i=0;i + void + compute_node_matrix( + FullMatrix& N, + const FiniteElement& fe) + { + const unsigned int n_dofs = fe.dofs_per_cell; + Assert (fe.has_generalized_support_points(), ExcNotInitialized()); + Assert (N.n()==n_dofs, ExcDimensionMismatch(N.n(), n_dofs)); + Assert (N.m()==n_dofs, ExcDimensionMismatch(N.m(), n_dofs)); + + const std::vector >& points = fe.get_generalized_support_points(); + + // We need the values of the + // polynomials in all generalized + // support points. + std::vector > + values (dim, std::vector(points.size())); + + // In this vector, we store the + // result of the interpolation + std::vector local_dofs(n_dofs); + + // One row per shape + // function. Remember that these + // are the 'raw' shape functions + // where the inverse node matrix is + // empty. Otherwise, this would + // yield identity. + for (unsigned int i=0;i -void -FETools::compute_embedding_matrices(const FiniteElement<1,2> &, - std::vector > > &, - const bool) -{ - Assert(false, ExcNotImplemented()); -} + template<> + void + compute_embedding_matrices(const FiniteElement<1,2> &, + std::vector > > &, + const bool) + { + Assert(false, ExcNotImplemented()); + } #elif deal_II_dimension == 2 -template<> -void -FETools::compute_embedding_matrices(const FiniteElement<2,3>&, - std::vector > >&, - const bool) -{ - Assert(false, ExcNotImplemented()); -} + template<> + void + compute_embedding_matrices(const FiniteElement<2,3>&, + std::vector > >&, + const bool) + { + Assert(false, ExcNotImplemented()); + } #endif -namespace { - template - void - compute_embedding_matrices_for_refinement_case (const FiniteElement& fe, - std::vector >& matrices, - const unsigned int ref_case) - { - const unsigned int n = fe.dofs_per_cell; - const unsigned int nc = GeometryInfo::n_children(RefinementCase(ref_case)); - for (unsigned int i = 0; i < nc; ++i) - { - Assert(matrices[i].n() == n, ExcDimensionMismatch(matrices[i].n (), n)); - Assert(matrices[i].m() == n, ExcDimensionMismatch(matrices[i].m (), n)); - } + namespace { + template + void + compute_embedding_matrices_for_refinement_case (const FiniteElement& fe, + std::vector >& matrices, + const unsigned int ref_case) + { + const unsigned int n = fe.dofs_per_cell; + const unsigned int nc = GeometryInfo::n_children(RefinementCase(ref_case)); + for (unsigned int i = 0; i < nc; ++i) + { + Assert(matrices[i].n() == n, ExcDimensionMismatch(matrices[i].n (), n)); + Assert(matrices[i].m() == n, ExcDimensionMismatch(matrices[i].m (), n)); + } - // Set up meshes, one with a single - // reference cell and refine it once - Triangulation tria; - GridGenerator::hyper_cube (tria, 0, 1); - tria.begin_active()->set_refine_flag (RefinementCase(ref_case)); - tria.execute_coarsening_and_refinement (); + // Set up meshes, one with a single + // reference cell and refine it once + Triangulation tria; + GridGenerator::hyper_cube (tria, 0, 1); + tria.begin_active()->set_refine_flag (RefinementCase(ref_case)); + tria.execute_coarsening_and_refinement (); - MappingCartesian mapping; - const unsigned int degree = fe.degree; - QGauss q_fine (degree+1); - const unsigned int nq = q_fine.size(); + MappingCartesian mapping; + const unsigned int degree = fe.degree; + QGauss q_fine (degree+1); + const unsigned int nq = q_fine.size(); - FEValues fine (mapping, fe, q_fine, - update_quadrature_points | - update_JxW_values | - update_values); + FEValues fine (mapping, fe, q_fine, + update_quadrature_points | + update_JxW_values | + update_values); // We search for the polynomial on // the small cell, being equal to @@ -643,70 +645,70 @@ namespace { // This matrix is the same for all // children. - fine.reinit (tria.begin_active ()); - const unsigned int nd = fe.n_components (); - FullMatrix A (nq*nd, n); - - for (unsigned int j = 0; j < n; ++j) - for (unsigned int d = 0; d < nd; ++d) - for (unsigned int k = 0; k < nq; ++k) - A (k * nd + d, j) = fine.shape_value_component (j, k, d); - - Householder H (A); - static Threads::Mutex mutex; - Vector v_coarse (nq * nd); - Vector v_fine (n); - unsigned int cell_number = 0; - - for (typename Triangulation::active_cell_iterator - fine_cell = tria.begin_active (); fine_cell != tria.end (); - ++fine_cell, ++cell_number) - { - fine.reinit (fine_cell); + fine.reinit (tria.begin_active ()); + const unsigned int nd = fe.n_components (); + FullMatrix A (nq*nd, n); + + for (unsigned int j = 0; j < n; ++j) + for (unsigned int d = 0; d < nd; ++d) + for (unsigned int k = 0; k < nq; ++k) + A (k * nd + d, j) = fine.shape_value_component (j, k, d); + + Householder H (A); + static Threads::Mutex mutex; + Vector v_coarse (nq * nd); + Vector v_fine (n); + unsigned int cell_number = 0; + + for (typename Triangulation::active_cell_iterator + fine_cell = tria.begin_active (); fine_cell != tria.end (); + ++fine_cell, ++cell_number) + { + fine.reinit (fine_cell); // evaluate on the coarse cell (which // is the first -- inactive -- cell on // the lowest level of the // triangulation we have created) - const Quadrature q_coarse (fine.get_quadrature_points (), - fine.get_JxW_values ()); - FEValues coarse (mapping, fe, q_coarse, update_values); - - coarse.reinit (tria.begin (0)); + const Quadrature q_coarse (fine.get_quadrature_points (), + fine.get_JxW_values ()); + FEValues coarse (mapping, fe, q_coarse, update_values); - FullMatrix &this_matrix = matrices[cell_number]; - - v_coarse = 0; + coarse.reinit (tria.begin (0)); + + FullMatrix &this_matrix = matrices[cell_number]; + + v_coarse = 0; // Compute this once for each // coarse grid basis function - for (unsigned int i = 0;i < n; ++i) - { + for (unsigned int i = 0;i < n; ++i) + { // The right hand side of // the least squares // problem consists of the // function values of the // coarse grid function in // each quadrature point. - if (fe.is_primitive ()) - { - const unsigned int - d = fe.system_to_component_index (i).first; - const double* phi_i = &coarse.shape_value (i, 0); - - for (unsigned int k = 0; k < nq; ++k) - v_coarse (k * nd + d) = phi_i[k]; - } - - else - for (unsigned int d = 0; d < nd; ++d) - for (unsigned int k = 0; k < nq; ++k) - v_coarse (k * nd + d) = coarse.shape_value_component (i, k, d); + if (fe.is_primitive ()) + { + const unsigned int + d = fe.system_to_component_index (i).first; + const double* phi_i = &coarse.shape_value (i, 0); + + for (unsigned int k = 0; k < nq; ++k) + v_coarse (k * nd + d) = phi_i[k]; + } + + else + for (unsigned int d = 0; d < nd; ++d) + for (unsigned int k = 0; k < nq; ++k) + v_coarse (k * nd + d) = coarse.shape_value_component (i, k, d); // solve the least squares // problem. - const double result = H.least_squares (v_fine, v_coarse); - Assert (result < 1.e-12, FETools::ExcLeastSquaresError (result)); + const double result = H.least_squares (v_fine, v_coarse); + Assert (result < 1.e-12, ExcLeastSquaresError (result)); // Copy into the result // matrix. Since the matrix @@ -714,411 +716,411 @@ namespace { // function to a fine grid // function, the columns // are fine grid. - mutex.acquire (); - - for (unsigned int j = 0; j < n; ++j) - this_matrix(j, i) = v_fine(j); - - mutex.release (); - } - - mutex.acquire (); + mutex.acquire (); + + for (unsigned int j = 0; j < n; ++j) + this_matrix(j, i) = v_fine(j); + + mutex.release (); + } + + mutex.acquire (); // Remove small entries from // the matrix - for (unsigned int i = 0; i < this_matrix.m (); ++i) - for (unsigned int j = 0; j < this_matrix.n (); ++j) - if (std::fabs (this_matrix (i, j)) < 1e-12) - this_matrix (i, j) = 0.; - - mutex.release (); - } - - Assert (cell_number == GeometryInfo::n_children (RefinementCase (ref_case)), - ExcInternalError ()); - } -} + for (unsigned int i = 0; i < this_matrix.m (); ++i) + for (unsigned int j = 0; j < this_matrix.n (); ++j) + if (std::fabs (this_matrix (i, j)) < 1e-12) + this_matrix (i, j) = 0.; + mutex.release (); + } -// This function is tested by tests/fe/internals, since it produces the matrices printed there -template -void -FETools::compute_embedding_matrices(const FiniteElement& fe, - std::vector > >& matrices, - const bool isotropic_only) -{ - Threads::TaskGroup task_group; - - // loop over all possible refinement cases - unsigned int ref_case = (isotropic_only) - ? RefinementCase::isotropic_refinement - : RefinementCase::cut_x; - - for (;ref_case <= RefinementCase::isotropic_refinement; ++ref_case) - task_group += Threads::new_task (&compute_embedding_matrices_for_refinement_case, - fe, matrices[ref_case-1], ref_case); - - task_group.join_all (); -} + Assert (cell_number == GeometryInfo::n_children (RefinementCase (ref_case)), + ExcInternalError ()); + } + } // This function is tested by tests/fe/internals, since it produces the matrices printed there + template + void + compute_embedding_matrices(const FiniteElement& fe, + std::vector > >& matrices, + const bool isotropic_only) + { + Threads::TaskGroup task_group; -//TODO:[GK] Is this correct for vector valued? -template -void -FETools::compute_face_embedding_matrices(const FiniteElement& fe, - FullMatrix (&matrices)[GeometryInfo::max_children_per_face], - const unsigned int face_coarse, - const unsigned int face_fine) -{ - const unsigned int nc = GeometryInfo::max_children_per_face; - const unsigned int n = fe.dofs_per_face; - const unsigned int nd = fe.n_components(); - const unsigned int degree = fe.degree; + // loop over all possible refinement cases + unsigned int ref_case = (isotropic_only) + ? RefinementCase::isotropic_refinement + : RefinementCase::cut_x; - for (unsigned int i=0;i::isotropic_refinement; ++ref_case) + task_group += Threads::new_task (&compute_embedding_matrices_for_refinement_case, + fe, matrices[ref_case-1], ref_case); - // Set up meshes, one with a single - // reference cell and refine it once - Triangulation tria; - GridGenerator::hyper_cube (tria, 0, 1); - tria.refine_global(1); - - MappingCartesian mapping; - QGauss q_gauss(degree+1); - const Quadrature q_fine = QProjector::project_to_face(q_gauss, face_fine); - - const unsigned int nq = q_fine.size(); - - // In order to make the loops below - // simpler, we introduce vectors - // containing for indices 0-n the - // number of the corresponding - // shape value on the cell. - std::vector face_c_dofs(n); - std::vector face_f_dofs(n); - unsigned int k=0; - for (unsigned int i=0;i::vertices_per_face;++i) - { - const unsigned int offset_c = GeometryInfo::face_to_cell_vertices(face_coarse, i) - *fe.dofs_per_vertex; - const unsigned int offset_f = GeometryInfo::face_to_cell_vertices(face_fine, i) - *fe.dofs_per_vertex; - for (unsigned int j=0;j::lines_per_face;++i) - { - const unsigned int offset_c = fe.first_line_index - + GeometryInfo::face_to_cell_lines(face_coarse, i-1) - *fe.dofs_per_line; - const unsigned int offset_f = fe.first_line_index - + GeometryInfo::face_to_cell_lines(face_fine, i-1) - *fe.dofs_per_line; - for (unsigned int j=0;j::quads_per_face;++i) - { - const unsigned int offset_c = fe.first_quad_index - + face_coarse - *fe.dofs_per_quad; - const unsigned int offset_f = fe.first_quad_index - + face_fine - *fe.dofs_per_quad; - for (unsigned int j=0;j fine (mapping, fe, q_fine, - update_quadrature_points | update_JxW_values | update_values); - // We search for the polynomial on - // the small cell, being equal to - // the coarse polynomial in all - // quadrature points. +// This function is tested by tests/fe/internals, since it produces the matrices printed there - // First build the matrix for this - // least squares problem. This - // contains the values of the fine - // cell polynomials in the fine - // cell grid points. +//TODO:[GK] Is this correct for vector valued? + template + void + compute_face_embedding_matrices(const FiniteElement& fe, + FullMatrix (&matrices)[GeometryInfo::max_children_per_face], + const unsigned int face_coarse, + const unsigned int face_fine) + { + const unsigned int nc = GeometryInfo::max_children_per_face; + const unsigned int n = fe.dofs_per_face; + const unsigned int nd = fe.n_components(); + const unsigned int degree = fe.degree; - // This matrix is the same for all - // children. - fine.reinit(tria.begin_active()); - FullMatrix A(nq*nd, n); - for (unsigned int j=0;j H(A); + // Set up meshes, one with a single + // reference cell and refine it once + Triangulation tria; + GridGenerator::hyper_cube (tria, 0, 1); + tria.refine_global(1); - Vector v_coarse(nq*nd); - Vector v_fine(n); + MappingCartesian mapping; + QGauss q_gauss(degree+1); + const Quadrature q_fine = QProjector::project_to_face(q_gauss, face_fine); + const unsigned int nq = q_fine.size(); + // In order to make the loops below + // simpler, we introduce vectors + // containing for indices 0-n the + // number of the corresponding + // shape value on the cell. + std::vector face_c_dofs(n); + std::vector face_f_dofs(n); + unsigned int k=0; + for (unsigned int i=0;i::vertices_per_face;++i) + { + const unsigned int offset_c = GeometryInfo::face_to_cell_vertices(face_coarse, i) + *fe.dofs_per_vertex; + const unsigned int offset_f = GeometryInfo::face_to_cell_vertices(face_fine, i) + *fe.dofs_per_vertex; + for (unsigned int j=0;j::lines_per_face;++i) + { + const unsigned int offset_c = fe.first_line_index + + GeometryInfo::face_to_cell_lines(face_coarse, i-1) + *fe.dofs_per_line; + const unsigned int offset_f = fe.first_line_index + + GeometryInfo::face_to_cell_lines(face_fine, i-1) + *fe.dofs_per_line; + for (unsigned int j=0;j::quads_per_face;++i) + { + const unsigned int offset_c = fe.first_quad_index + + face_coarse + *fe.dofs_per_quad; + const unsigned int offset_f = fe.first_quad_index + + face_fine + *fe.dofs_per_quad; + for (unsigned int j=0;j::max_children_per_face; - ++cell_number) - { - const Quadrature q_coarse - = QProjector::project_to_subface(q_gauss, face_coarse, cell_number); - FEValues coarse (mapping, fe, q_coarse, update_values); + FEValues fine (mapping, fe, q_fine, + update_quadrature_points | update_JxW_values | update_values); - typename Triangulation::active_cell_iterator fine_cell - = tria.begin(0)->child(GeometryInfo::child_cell_on_face( - tria.begin(0)->refinement_case(), face_coarse, cell_number)); - fine.reinit(fine_cell); - coarse.reinit(tria.begin(0)); + // We search for the polynomial on + // the small cell, being equal to + // the coarse polynomial in all + // quadrature points. - FullMatrix &this_matrix = matrices[cell_number]; + // First build the matrix for this + // least squares problem. This + // contains the values of the fine + // cell polynomials in the fine + // cell grid points. - // Compute this once for each - // coarse grid basis function - for (unsigned int i=0;i A(nq*nd, n); + for (unsigned int j=0;j -void -FETools::compute_projection_matrices(const FiniteElement<1,2>&, - std::vector > >&, bool) -{ - Assert(false, ExcNotImplemented()); -} + Householder H(A); -#elif deal_II_dimension == 2 -template <> -void -FETools::compute_projection_matrices(const FiniteElement<2,3>&, - std::vector > >&, bool) -{ - Assert(false, ExcNotImplemented()); -} + Vector v_coarse(nq*nd); + Vector v_fine(n); -#endif -template -void -FETools::compute_projection_matrices(const FiniteElement& fe, - std::vector > >& matrices, - const bool isotropic_only) -{ - const unsigned int n = fe.dofs_per_cell; - const unsigned int nd = fe.n_components(); - const unsigned int degree = fe.degree; - - // prepare FEValues, quadrature etc on - // coarse cell - MappingCartesian mapping; - QGauss q_fine(degree+1); - const unsigned int nq = q_fine.size(); - - // create mass matrix on coarse cell. - FullMatrix mass(n, n); - { - // set up a triangulation for coarse cell - Triangulation tr; - GridGenerator::hyper_cube (tr, 0, 1); + for (unsigned int cell_number = 0; cell_number < GeometryInfo::max_children_per_face; + ++cell_number) + { + const Quadrature q_coarse + = QProjector::project_to_subface(q_gauss, face_coarse, cell_number); + FEValues coarse (mapping, fe, q_coarse, update_values); - FEValues coarse (mapping, fe, q_fine, - update_JxW_values | update_values); + typename Triangulation::active_cell_iterator fine_cell + = tria.begin(0)->child(GeometryInfo::child_cell_on_face( + tria.begin(0)->refinement_case(), face_coarse, cell_number)); + fine.reinit(fine_cell); + coarse.reinit(tria.begin(0)); - typename Triangulation::cell_iterator coarse_cell - = tr.begin(0); - coarse.reinit (coarse_cell); + FullMatrix &this_matrix = matrices[cell_number]; - const std::vector & JxW = coarse.get_JxW_values(); - for (unsigned int i=0;i + void + compute_projection_matrices(const FiniteElement<1,2>&, + std::vector > >&, bool) + { + Assert(false, ExcNotImplemented()); } - // loop over all possible - // refinement cases - unsigned int ref_case = (isotropic_only) - ? RefinementCase::isotropic_refinement - : RefinementCase::cut_x; - for (;ref_case <= RefinementCase::isotropic_refinement; ++ref_case) - { - const unsigned int - nc = GeometryInfo::n_children(RefinementCase(ref_case)); +#elif deal_II_dimension == 2 + template <> + void + compute_projection_matrices(const FiniteElement<2,3>&, + std::vector > >&, bool) + { + Assert(false, ExcNotImplemented()); + } - for (unsigned int i=0;i + void + compute_projection_matrices(const FiniteElement& fe, + std::vector > >& matrices, + const bool isotropic_only) + { + const unsigned int n = fe.dofs_per_cell; + const unsigned int nd = fe.n_components(); + const unsigned int degree = fe.degree; + + // prepare FEValues, quadrature etc on + // coarse cell + MappingCartesian mapping; + QGauss q_fine(degree+1); + const unsigned int nq = q_fine.size(); + + // create mass matrix on coarse cell. + FullMatrix mass(n, n); + { + // set up a triangulation for coarse cell Triangulation tr; GridGenerator::hyper_cube (tr, 0, 1); - tr.begin_active()->set_refine_flag(RefinementCase(ref_case)); - tr.execute_coarsening_and_refinement(); - FEValues fine (mapping, fe, q_fine, - update_quadrature_points | update_JxW_values | - update_values); + FEValues coarse (mapping, fe, q_fine, + update_JxW_values | update_values); typename Triangulation::cell_iterator coarse_cell = tr.begin(0); + coarse.reinit (coarse_cell); - Vector v_coarse(n); - Vector v_fine(n); + const std::vector & JxW = coarse.get_JxW_values(); + for (unsigned int i=0;i &this_matrix = matrices[ref_case-1][cell_number]; - - // Compute right hand side, - // which is a fine level basis - // function tested with the - // coarse level functions. - fine.reinit(coarse_cell->child(cell_number)); - Quadrature q_coarse (fine.get_quadrature_points(), - fine.get_JxW_values()); - FEValues coarse (mapping, fe, q_coarse, update_values); - coarse.reinit(coarse_cell); + // invert mass matrix + mass.gauss_jordan(); + } - // Build RHS + // loop over all possible + // refinement cases + unsigned int ref_case = (isotropic_only) + ? RefinementCase::isotropic_refinement + : RefinementCase::cut_x; + for (;ref_case <= RefinementCase::isotropic_refinement; ++ref_case) + { + const unsigned int + nc = GeometryInfo::n_children(RefinementCase(ref_case)); - const std::vector & JxW = fine.get_JxW_values(); + for (unsigned int i=0;i tr; + GridGenerator::hyper_cube (tr, 0, 1); + tr.begin_active()->set_refine_flag(RefinementCase(ref_case)); + tr.execute_coarsening_and_refinement(); - double update = 0; - for (unsigned int k=0; k fine (mapping, fe, q_fine, + update_quadrature_points | update_JxW_values | + update_values); + + typename Triangulation::cell_iterator coarse_cell + = tr.begin(0); + + Vector v_coarse(n); + Vector v_fine(n); + + for (unsigned int cell_number=0;cell_number &this_matrix = matrices[ref_case-1][cell_number]; + + // Compute right hand side, + // which is a fine level basis + // function tested with the + // coarse level functions. + fine.reinit(coarse_cell->child(cell_number)); + Quadrature q_coarse (fine.get_quadrature_points(), + fine.get_JxW_values()); + FEValues coarse (mapping, fe, q_coarse, update_values); + coarse.reinit(coarse_cell); + + // Build RHS + + const std::vector & JxW = fine.get_JxW_values(); + + // Outer loop over all fine + // grid shape functions phi_j + for (unsigned int j=0;j class DH1, - template class DH2, - class InVector, class OutVector> -void -FETools::interpolate(const DH1 &dof1, - const InVector &u1, - const DH2 &dof2, - OutVector &u2) -{ - ConstraintMatrix dummy; - dummy.close(); - interpolate(dof1, u1, dof2, dummy, u2); -} + template class DH1, + template class DH2, + class InVector, class OutVector> + void + interpolate(const DH1 &dof1, + const InVector &u1, + const DH2 &dof2, + OutVector &u2) + { + ConstraintMatrix dummy; + dummy.close(); + interpolate(dof1, u1, dof2, dummy, u2); + } @@ -1127,11 +1129,11 @@ template class DH2, class InVector, class OutVector> void -FETools::interpolate (const DH1 &dof1, - const InVector &u1, - const DH2 &dof2, - const ConstraintMatrix &constraints, - OutVector &u2) +interpolate (const DH1 &dof1, + const InVector &u1, + const DH2 &dof2, + const ConstraintMatrix &constraints, + OutVector &u2) { Assert(&dof1.get_tria()==&dof2.get_tria(), ExcTriangulationMismatch()); @@ -1152,8 +1154,8 @@ FETools::interpolate (const DH1 &dof1, // matrices. shared_ptr make sure // that memory is released again std::map *, - std::map *, - std_cxx1x::shared_ptr > > > + std::map *, + std_cxx1x::shared_ptr > > > interpolation_matrices; typename DH1::active_cell_iterator cell1 = dof1.begin_active(), @@ -1208,9 +1210,9 @@ FETools::interpolate (const DH1 &dof1, interpolation_matrices[&cell1->get_fe()][&cell2->get_fe()] = interpolation_matrix; - FETools::get_interpolation_matrix(cell1->get_fe(), - cell2->get_fe(), - *interpolation_matrix); + get_interpolation_matrix(cell1->get_fe(), + cell2->get_fe(), + *interpolation_matrix); } cell1->get_dof_values(u1, u1_local); @@ -1245,1442 +1247,1440 @@ FETools::interpolate (const DH1 &dof1, -template -void -FETools::back_interpolate(const DoFHandler &dof1, - const InVector &u1, - const FiniteElement &fe2, - OutVector &u1_interpolated) -{ - Assert(dof1.get_fe().n_components() == fe2.n_components(), - ExcDimensionMismatch(dof1.get_fe().n_components(), fe2.n_components())); - Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); - Assert(u1_interpolated.size()==dof1.n_dofs(), - ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs())); - - // For continuous elements on grids - // with hanging nodes we need - // hanging node - // constraints. Consequently, when - // the elements are continuous no - // hanging node constraints are - // allowed. - const bool hanging_nodes_not_allowed= - (dof1.get_fe().dofs_per_vertex != 0) || (fe2.dofs_per_vertex != 0); - - const unsigned int dofs_per_cell1=dof1.get_fe().dofs_per_cell; - - Vector u1_local(dofs_per_cell1); - Vector u1_int_local(dofs_per_cell1); - - typename DoFHandler::active_cell_iterator cell = dof1.begin_active(), - endc = dof1.end(); - - FullMatrix interpolation_matrix(dofs_per_cell1, dofs_per_cell1); - FETools::get_back_interpolation_matrix(dof1.get_fe(), fe2, - interpolation_matrix); - for (; cell!=endc; ++cell) - { - if (hanging_nodes_not_allowed) - for (unsigned int face=0; face::faces_per_cell; ++face) - Assert (cell->at_boundary(face) || - cell->neighbor(face)->level() == cell->level(), - ExcHangingNodesNotAllowed(0)); - - cell->get_dof_values(u1, u1_local); - interpolation_matrix.vmult(u1_int_local, u1_local); - cell->set_dof_values(u1_int_local, u1_interpolated); - } -} - - - -template class DH, - class InVector, class OutVector, int spacedim> -void -FETools::back_interpolate(const DH &dof1, - const InVector &u1, - const FiniteElement &fe2, - OutVector &u1_interpolated) -{ - Assert(u1.size() == dof1.n_dofs(), - ExcDimensionMismatch(u1.size(), dof1.n_dofs())); - Assert(u1_interpolated.size() == dof1.n_dofs(), - ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs())); - - Vector u1_local(DoFTools::max_dofs_per_cell(dof1)); - Vector u1_int_local(DoFTools::max_dofs_per_cell(dof1)); - - typename DH::active_cell_iterator cell = dof1.begin_active(), - endc = dof1.end(); - - // map from possible fe objects in - // dof1 to the back_interpolation - // matrices - std::map *, - std_cxx1x::shared_ptr > > interpolation_matrices; - - for (; cell!=endc; ++cell) - { - Assert(cell->get_fe().n_components() == fe2.n_components(), - ExcDimensionMismatch(cell->get_fe().n_components(), - fe2.n_components())); - - // For continuous elements on - // grids with hanging nodes we - // need hanging node - // constraints. Consequently, - // when the elements are - // continuous no hanging node - // constraints are allowed. - const bool hanging_nodes_not_allowed= - (cell->get_fe().dofs_per_vertex != 0) || (fe2.dofs_per_vertex != 0); - - if (hanging_nodes_not_allowed) - for (unsigned int face=0; face::faces_per_cell; ++face) - Assert (cell->at_boundary(face) || - cell->neighbor(face)->level() == cell->level(), - ExcHangingNodesNotAllowed(0)); - - const unsigned int dofs_per_cell1 = cell->get_fe().dofs_per_cell; - - // make sure back_interpolation - // matrix is available - if (interpolation_matrices[&cell->get_fe()] != 0) - { - interpolation_matrices[&cell->get_fe()] = - std_cxx1x::shared_ptr > - (new FullMatrix(dofs_per_cell1, dofs_per_cell1)); - get_back_interpolation_matrix(dof1.get_fe(), fe2, - *interpolation_matrices[&cell->get_fe()]); - } + template + void + back_interpolate(const DoFHandler &dof1, + const InVector &u1, + const FiniteElement &fe2, + OutVector &u1_interpolated) + { + Assert(dof1.get_fe().n_components() == fe2.n_components(), + ExcDimensionMismatch(dof1.get_fe().n_components(), fe2.n_components())); + Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); + Assert(u1_interpolated.size()==dof1.n_dofs(), + ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs())); + + // For continuous elements on grids + // with hanging nodes we need + // hanging node + // constraints. Consequently, when + // the elements are continuous no + // hanging node constraints are + // allowed. + const bool hanging_nodes_not_allowed= + (dof1.get_fe().dofs_per_vertex != 0) || (fe2.dofs_per_vertex != 0); + + const unsigned int dofs_per_cell1=dof1.get_fe().dofs_per_cell; + + Vector u1_local(dofs_per_cell1); + Vector u1_int_local(dofs_per_cell1); + + typename DoFHandler::active_cell_iterator cell = dof1.begin_active(), + endc = dof1.end(); + + FullMatrix interpolation_matrix(dofs_per_cell1, dofs_per_cell1); + get_back_interpolation_matrix(dof1.get_fe(), fe2, + interpolation_matrix); + for (; cell!=endc; ++cell) + { + if (hanging_nodes_not_allowed) + for (unsigned int face=0; face::faces_per_cell; ++face) + Assert (cell->at_boundary(face) || + cell->neighbor(face)->level() == cell->level(), + ExcHangingNodesNotAllowed(0)); + + cell->get_dof_values(u1, u1_local); + interpolation_matrix.vmult(u1_int_local, u1_local); + cell->set_dof_values(u1_int_local, u1_interpolated); + } + } - u1_local.reinit (dofs_per_cell1); - u1_int_local.reinit (dofs_per_cell1); - cell->get_dof_values(u1, u1_local); - interpolation_matrices[&cell->get_fe()]->vmult(u1_int_local, u1_local); - cell->set_dof_values(u1_int_local, u1_interpolated); - } -} + template class DH, + class InVector, class OutVector, int spacedim> + void + back_interpolate(const DH &dof1, + const InVector &u1, + const FiniteElement &fe2, + OutVector &u1_interpolated) + { + Assert(u1.size() == dof1.n_dofs(), + ExcDimensionMismatch(u1.size(), dof1.n_dofs())); + Assert(u1_interpolated.size() == dof1.n_dofs(), + ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs())); + Vector u1_local(DoFTools::max_dofs_per_cell(dof1)); + Vector u1_int_local(DoFTools::max_dofs_per_cell(dof1)); -template -void FETools::back_interpolate(const DoFHandler &dof1, - const ConstraintMatrix &constraints1, - const InVector &u1, - const DoFHandler &dof2, - const ConstraintMatrix &constraints2, - OutVector &u1_interpolated) -{ - // For discontinuous elements - // without constraints take the - // simpler version of the - // back_interpolate function. - if (dof1.get_fe().dofs_per_vertex==0 && dof2.get_fe().dofs_per_vertex==0 - && constraints1.n_constraints()==0 && constraints2.n_constraints()==0) - back_interpolate(dof1, u1, dof2.get_fe(), u1_interpolated); - else - { - Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(), - ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components())); - Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); - Assert(u1_interpolated.size()==dof1.n_dofs(), - ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs())); - - // For continuous elements - // first interpolate to dof2, - // taking into account - // constraints2, and then - // interpolate back to dof1 - // taking into account - // constraints1 - Vector u2(dof2.n_dofs()); - interpolate(dof1, u1, dof2, constraints2, u2); - interpolate(dof2, u2, dof1, constraints1, u1_interpolated); - } -} + typename DH::active_cell_iterator cell = dof1.begin_active(), + endc = dof1.end(); + // map from possible fe objects in + // dof1 to the back_interpolation + // matrices + std::map *, + std_cxx1x::shared_ptr > > interpolation_matrices; +for (; cell!=endc; ++cell) + { + Assert(cell->get_fe().n_components() == fe2.n_components(), + ExcDimensionMismatch(cell->get_fe().n_components(), + fe2.n_components())); + + // For continuous elements on + // grids with hanging nodes we + // need hanging node + // constraints. Consequently, + // when the elements are + // continuous no hanging node + // constraints are allowed. + const bool hanging_nodes_not_allowed= + (cell->get_fe().dofs_per_vertex != 0) || (fe2.dofs_per_vertex != 0); + + if (hanging_nodes_not_allowed) + for (unsigned int face=0; face::faces_per_cell; ++face) + Assert (cell->at_boundary(face) || + cell->neighbor(face)->level() == cell->level(), + ExcHangingNodesNotAllowed(0)); + + const unsigned int dofs_per_cell1 = cell->get_fe().dofs_per_cell; + + // make sure back_interpolation + // matrix is available + if (interpolation_matrices[&cell->get_fe()] != 0) + { + interpolation_matrices[&cell->get_fe()] = + std_cxx1x::shared_ptr > + (new FullMatrix(dofs_per_cell1, dofs_per_cell1)); + get_back_interpolation_matrix(dof1.get_fe(), fe2, + *interpolation_matrices[&cell->get_fe()]); + } -template -void FETools::interpolation_difference (const DoFHandler &dof1, - const InVector &u1, - const FiniteElement &fe2, - OutVector &u1_difference) -{ - Assert(dof1.get_fe().n_components() == fe2.n_components(), - ExcDimensionMismatch(dof1.get_fe().n_components(), fe2.n_components())); - Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); - Assert(u1_difference.size()==dof1.n_dofs(), - ExcDimensionMismatch(u1_difference.size(), dof1.n_dofs())); - - // For continuous elements on grids - // with hanging nodes we need - // hnaging node - // constraints. Consequently, when - // the elements are continuous no - // hanging node constraints are - // allowed. - const bool hanging_nodes_not_allowed= - (dof1.get_fe().dofs_per_vertex != 0) || (fe2.dofs_per_vertex != 0); - - const unsigned int dofs_per_cell=dof1.get_fe().dofs_per_cell; - - Vector u1_local(dofs_per_cell); - Vector u1_diff_local(dofs_per_cell); - - FullMatrix difference_matrix(dofs_per_cell, dofs_per_cell); - FETools::get_interpolation_difference_matrix(dof1.get_fe(), fe2, - difference_matrix); - - typename DoFHandler::active_cell_iterator cell = dof1.begin_active(), - endc = dof1.end(); - - for (; cell!=endc; ++cell) - { - if (hanging_nodes_not_allowed) - for (unsigned int face=0; face::faces_per_cell; ++face) - Assert (cell->at_boundary(face) || - cell->neighbor(face)->level() == cell->level(), - ExcHangingNodesNotAllowed(0)); + u1_local.reinit (dofs_per_cell1); + u1_int_local.reinit (dofs_per_cell1); - cell->get_dof_values(u1, u1_local); - difference_matrix.vmult(u1_diff_local, u1_local); - cell->set_dof_values(u1_diff_local, u1_difference); - } + cell->get_dof_values(u1, u1_local); + interpolation_matrices[&cell->get_fe()]->vmult(u1_int_local, u1_local); + cell->set_dof_values(u1_int_local, u1_interpolated); + } } -template -void FETools::interpolation_difference(const DoFHandler &dof1, - const ConstraintMatrix &constraints1, - const InVector &u1, - const DoFHandler &dof2, - const ConstraintMatrix &constraints2, - OutVector &u1_difference) -{ - // For discontinuous elements - // without constraints take the - // cheaper version of the - // interpolation_difference function. - if (dof1.get_fe().dofs_per_vertex==0 && dof2.get_fe().dofs_per_vertex==0 - && constraints1.n_constraints()==0 && constraints2.n_constraints()==0) - interpolation_difference(dof1, u1, dof2.get_fe(), u1_difference); - else - { - back_interpolate(dof1, constraints1, u1, dof2, constraints2, u1_difference); - u1_difference.sadd(-1, u1); - } -} + template + void back_interpolate(const DoFHandler &dof1, + const ConstraintMatrix &constraints1, + const InVector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + OutVector &u1_interpolated) + { + // For discontinuous elements + // without constraints take the + // simpler version of the + // back_interpolate function. + if (dof1.get_fe().dofs_per_vertex==0 && dof2.get_fe().dofs_per_vertex==0 + && constraints1.n_constraints()==0 && constraints2.n_constraints()==0) + back_interpolate(dof1, u1, dof2.get_fe(), u1_interpolated); + else + { + Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(), + ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components())); + Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); + Assert(u1_interpolated.size()==dof1.n_dofs(), + ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs())); + + // For continuous elements + // first interpolate to dof2, + // taking into account + // constraints2, and then + // interpolate back to dof1 + // taking into account + // constraints1 + Vector u2(dof2.n_dofs()); + interpolate(dof1, u1, dof2, constraints2, u2); + interpolate(dof2, u2, dof1, constraints1, u1_interpolated); + } + } -template -void FETools::project_dg(const DoFHandler &dof1, - const InVector &u1, - const DoFHandler &dof2, - OutVector &u2) -{ - Assert(&dof1.get_tria()==&dof2.get_tria(), ExcTriangulationMismatch()); - Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(), - ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components())); - Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); - Assert(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs())); - typename DoFHandler::active_cell_iterator cell1 = dof1.begin_active(); - typename DoFHandler::active_cell_iterator cell2 = dof2.begin_active(); - typename DoFHandler::active_cell_iterator end = dof2.end(); + template + void interpolation_difference (const DoFHandler &dof1, + const InVector &u1, + const FiniteElement &fe2, + OutVector &u1_difference) + { + Assert(dof1.get_fe().n_components() == fe2.n_components(), + ExcDimensionMismatch(dof1.get_fe().n_components(), fe2.n_components())); + Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); + Assert(u1_difference.size()==dof1.n_dofs(), + ExcDimensionMismatch(u1_difference.size(), dof1.n_dofs())); + + // For continuous elements on grids + // with hanging nodes we need + // hnaging node + // constraints. Consequently, when + // the elements are continuous no + // hanging node constraints are + // allowed. + const bool hanging_nodes_not_allowed= + (dof1.get_fe().dofs_per_vertex != 0) || (fe2.dofs_per_vertex != 0); + + const unsigned int dofs_per_cell=dof1.get_fe().dofs_per_cell; + + Vector u1_local(dofs_per_cell); + Vector u1_diff_local(dofs_per_cell); + + FullMatrix difference_matrix(dofs_per_cell, dofs_per_cell); + get_interpolation_difference_matrix(dof1.get_fe(), fe2, + difference_matrix); + + typename DoFHandler::active_cell_iterator cell = dof1.begin_active(), + endc = dof1.end(); + + for (; cell!=endc; ++cell) + { + if (hanging_nodes_not_allowed) + for (unsigned int face=0; face::faces_per_cell; ++face) + Assert (cell->at_boundary(face) || + cell->neighbor(face)->level() == cell->level(), + ExcHangingNodesNotAllowed(0)); + + cell->get_dof_values(u1, u1_local); + difference_matrix.vmult(u1_diff_local, u1_local); + cell->set_dof_values(u1_diff_local, u1_difference); + } + } - const unsigned int n1 = dof1.get_fe().dofs_per_cell; - const unsigned int n2 = dof2.get_fe().dofs_per_cell; - Vector u1_local(n1); - Vector u2_local(n2); - std::vector dofs(n2); - FullMatrix matrix(n2,n1); - get_projection_matrix(dof1.get_fe(), dof2.get_fe(), matrix); + template + void interpolation_difference(const DoFHandler &dof1, + const ConstraintMatrix &constraints1, + const InVector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints2, + OutVector &u1_difference) + { + // For discontinuous elements + // without constraints take the + // cheaper version of the + // interpolation_difference function. + if (dof1.get_fe().dofs_per_vertex==0 && dof2.get_fe().dofs_per_vertex==0 + && constraints1.n_constraints()==0 && constraints2.n_constraints()==0) + interpolation_difference(dof1, u1, dof2.get_fe(), u1_difference); + else + { + back_interpolate(dof1, constraints1, u1, dof2, constraints2, u1_difference); + u1_difference.sadd(-1, u1); + } + } - while (cell2 != end) - { - cell1->get_dof_values(u1, u1_local); - matrix.vmult(u2_local, u1_local); - cell2->get_dof_indices(dofs); - for (unsigned int i=0; i + void project_dg(const DoFHandler &dof1, + const InVector &u1, + const DoFHandler &dof2, + OutVector &u2) + { + Assert(&dof1.get_tria()==&dof2.get_tria(), ExcTriangulationMismatch()); + Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(), + ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components())); + Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); + Assert(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs())); + typename DoFHandler::active_cell_iterator cell1 = dof1.begin_active(); + typename DoFHandler::active_cell_iterator cell2 = dof2.begin_active(); + typename DoFHandler::active_cell_iterator end = dof2.end(); -template -void FETools::extrapolate(const DoFHandler &dof1, - const InVector &u1, - const DoFHandler &dof2, - OutVector &u2) -{ - ConstraintMatrix dummy; - dummy.close(); - extrapolate(dof1, u1, dof2, dummy, u2); -} + const unsigned int n1 = dof1.get_fe().dofs_per_cell; + const unsigned int n2 = dof2.get_fe().dofs_per_cell; + Vector u1_local(n1); + Vector u2_local(n2); + std::vector dofs(n2); + FullMatrix matrix(n2,n1); + get_projection_matrix(dof1.get_fe(), dof2.get_fe(), matrix); -template -void FETools::extrapolate(const DoFHandler &dof1, - const InVector &u1, - const DoFHandler &dof2, - const ConstraintMatrix &constraints, - OutVector &u2) -{ - Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(), - ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components())); - Assert(&dof1.get_tria()==&dof2.get_tria(), ExcTriangulationMismatch()); - Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); - Assert(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs())); + while (cell2 != end) + { + cell1->get_dof_values(u1, u1_local); + matrix.vmult(u2_local, u1_local); + cell2->get_dof_indices(dofs); + for (unsigned int i=0; i dof_values(dofs_per_cell); - // make sure that each cell on the - // coarsest level is at least once - // refined. otherwise, we can't - // treat these cells and would - // generate a bogus result + template + void extrapolate(const DoFHandler &dof1, + const InVector &u1, + const DoFHandler &dof2, + OutVector &u2) { - typename DoFHandler::cell_iterator cell = dof2.begin(0), - endc = dof2.end(0); - for (; cell!=endc; ++cell) - Assert (cell->has_children(), ExcGridNotRefinedAtLeastOnce()); + ConstraintMatrix dummy; + dummy.close(); + extrapolate(dof1, u1, dof2, dummy, u2); } - // then traverse grid bottom up - for (unsigned int level=0; level::cell_iterator cell=dof2.begin(level), - endc=dof2.end(level); + + template + void extrapolate(const DoFHandler &dof1, + const InVector &u1, + const DoFHandler &dof2, + const ConstraintMatrix &constraints, + OutVector &u2) + { + Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(), + ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components())); + Assert(&dof1.get_tria()==&dof2.get_tria(), ExcTriangulationMismatch()); + Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs())); + Assert(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs())); + + OutVector u3; + u3.reinit(u2); + interpolate(dof1, u1, dof2, constraints, u3); + + const unsigned int dofs_per_cell = dof2.get_fe().dofs_per_cell; + Vector dof_values(dofs_per_cell); + + // make sure that each cell on the + // coarsest level is at least once + // refined. otherwise, we can't + // treat these cells and would + // generate a bogus result + { + typename DoFHandler::cell_iterator cell = dof2.begin(0), + endc = dof2.end(0); for (; cell!=endc; ++cell) - if (!cell->active()) - { - // check whether this - // cell has active - // children - bool active_children=false; - for (unsigned int child_n=0; child_nn_children(); ++child_n) - if (cell->child(child_n)->active()) + Assert (cell->has_children(), ExcGridNotRefinedAtLeastOnce()); + } + + // then traverse grid bottom up + for (unsigned int level=0; level::cell_iterator cell=dof2.begin(level), + endc=dof2.end(level); + + for (; cell!=endc; ++cell) + if (!cell->active()) + { + // check whether this + // cell has active + // children + bool active_children=false; + for (unsigned int child_n=0; child_nn_children(); ++child_n) + if (cell->child(child_n)->active()) + { + active_children=true; + break; + } + + // if there are active + // children, the we have + // to work on this + // cell. get the data + // from the one vector + // and set it on the + // other + if (active_children) { - active_children=true; - break; + cell->get_interpolated_dof_values(u3, dof_values); + cell->set_dof_values_by_interpolation(dof_values, u2); } + } + } - // if there are active - // children, the we have - // to work on this - // cell. get the data - // from the one vector - // and set it on the - // other - if (active_children) - { - cell->get_interpolated_dof_values(u3, dof_values); - cell->set_dof_values_by_interpolation(dof_values, u2); - } - } - } - - // Apply hanging node constraints. - constraints.distribute(u2); -} + // Apply hanging node constraints. + constraints.distribute(u2); + } -template -void -FETools::add_fe_name(const std::string& parameter_name, - const FEFactoryBase* factory) -{ - // Erase everything after the - // actual class name - std::string name = parameter_name; - unsigned int name_end = - name.find_first_not_of(std::string("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_")); - if (name_end < name.size()) - name.erase(name_end); - // first make sure that no other - // thread intercepts the - // operation of this function; - // for this, acquire the lock - // until we quit this function - Threads::ThreadMutex::ScopedLock lock(fe_name_map_lock); - - Assert(fe_name_map.find(name) == fe_name_map.end(), - ExcMessage("Cannot change existing element in finite element name list")); - - // Insert the normalized name into - // the map - fe_name_map[name] = FEFactoryPointer(factory); -} + template + void + add_fe_name(const std::string& parameter_name, + const FEFactoryBase* factory) + { + // Erase everything after the + // actual class name + std::string name = parameter_name; + unsigned int name_end = + name.find_first_not_of(std::string("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_")); + if (name_end < name.size()) + name.erase(name_end); + // first make sure that no other + // thread intercepts the + // operation of this function; + // for this, acquire the lock + // until we quit this function + Threads::ThreadMutex::ScopedLock lock(fe_name_map_lock); + + Assert(fe_name_map.find(name) == fe_name_map.end(), + ExcMessage("Cannot change existing element in finite element name list")); + + // Insert the normalized name into + // the map + fe_name_map[name] = FEFactoryPointer(factory); + } -namespace internal -{ - namespace + namespace internal { - template - FiniteElement* - get_fe_from_name (std::string &name) + namespace { - // Extract the name of the - // finite element class, which only - // contains characters, numbers and - // underscores. - unsigned int name_end = - name.find_first_not_of(std::string("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_")); - const std::string name_part(name, 0, name_end); - name.erase(0, name_part.size()); - - // now things get a little more - // complicated: FESystem. it's - // more complicated, since we - // have to figure out what the - // base elements are. this can - // only be done recursively - if (name_part == "FESystem") - { - // next we have to get at the - // base elements. start with - // the first. wrap the whole - // block into try-catch to - // make sure we destroy the - // pointers we got from - // recursive calls if one of - // these calls should throw - // an exception - std::vector*> base_fes; - std::vector base_multiplicities; - try - { - // Now, just the [...] - // part should be left. - if (name.size() == 0 || name[0] != '[') - throw (std::string("Invalid first character in ") + name); - do - { - // Erase the - // leading '[' or '-' - name.erase(0,1); - // Now, the name of the - // first base element is - // first... Let's get it - base_fes.push_back (get_fe_from_name (name)); - // next check whether - // FESystem placed a - // multiplicity after - // the element name - if (name[0] == '^') - { - // yes. Delete the '^' - // and read this - // multiplicity - name.erase(0,1); - - const std::pair tmp - = Utilities::get_integer_at_position (name, 0); - name.erase(0, tmp.second); - // add to length, - // including the '^' - base_multiplicities.push_back (tmp.first); - } - else - // no, so - // multiplicity is - // 1 - base_multiplicities.push_back (1); - - // so that's it for - // this base - // element. base - // elements are - // separated by '-', - // and the list is - // terminated by ']', - // so loop while the - // next character is - // '-' - } - while (name[0] == '-'); - - // so we got to the end - // of the '-' separated - // list. make sure that - // we actually had a ']' - // there - if (name.size() == 0 || name[0] != ']') - throw (std::string("Invalid first character in ") + name); - name.erase(0,1); - // just one more sanity check - Assert ((base_fes.size() == base_multiplicities.size()) - && - (base_fes.size() > 0), - ExcInternalError()); - - // ok, apparently - // everything went ok. so - // generate the composed - // element - FiniteElement *system_element = 0; - switch (base_fes.size()) - { - case 1: + template + FiniteElement* + get_fe_from_name (std::string &name) + { + // Extract the name of the + // finite element class, which only + // contains characters, numbers and + // underscores. + unsigned int name_end = + name.find_first_not_of(std::string("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_")); + const std::string name_part(name, 0, name_end); + name.erase(0, name_part.size()); + + // now things get a little more + // complicated: FESystem. it's + // more complicated, since we + // have to figure out what the + // base elements are. this can + // only be done recursively + if (name_part == "FESystem") + { + // next we have to get at the + // base elements. start with + // the first. wrap the whole + // block into try-catch to + // make sure we destroy the + // pointers we got from + // recursive calls if one of + // these calls should throw + // an exception + std::vector*> base_fes; + std::vector base_multiplicities; + try + { + // Now, just the [...] + // part should be left. + if (name.size() == 0 || name[0] != '[') + throw (std::string("Invalid first character in ") + name); + do { - system_element = new FESystem(*base_fes[0], - base_multiplicities[0]); - break; + // Erase the + // leading '[' or '-' + name.erase(0,1); + // Now, the name of the + // first base element is + // first... Let's get it + base_fes.push_back (get_fe_from_name (name)); + // next check whether + // FESystem placed a + // multiplicity after + // the element name + if (name[0] == '^') + { + // yes. Delete the '^' + // and read this + // multiplicity + name.erase(0,1); + + const std::pair tmp + = Utilities::get_integer_at_position (name, 0); + name.erase(0, tmp.second); + // add to length, + // including the '^' + base_multiplicities.push_back (tmp.first); + } + else + // no, so + // multiplicity is + // 1 + base_multiplicities.push_back (1); + + // so that's it for + // this base + // element. base + // elements are + // separated by '-', + // and the list is + // terminated by ']', + // so loop while the + // next character is + // '-' } - - case 2: + while (name[0] == '-'); + + // so we got to the end + // of the '-' separated + // list. make sure that + // we actually had a ']' + // there + if (name.size() == 0 || name[0] != ']') + throw (std::string("Invalid first character in ") + name); + name.erase(0,1); + // just one more sanity check + Assert ((base_fes.size() == base_multiplicities.size()) + && + (base_fes.size() > 0), + ExcInternalError()); + + // ok, apparently + // everything went ok. so + // generate the composed + // element + FiniteElement *system_element = 0; + switch (base_fes.size()) { - system_element = new FESystem(*base_fes[0], - base_multiplicities[0], - *base_fes[1], - base_multiplicities[1]); - break; - } + case 1: + { + system_element = new FESystem(*base_fes[0], + base_multiplicities[0]); + break; + } - case 3: - { - system_element = new FESystem(*base_fes[0], - base_multiplicities[0], - *base_fes[1], - base_multiplicities[1], - *base_fes[2], - base_multiplicities[2]); - break; - } + case 2: + { + system_element = new FESystem(*base_fes[0], + base_multiplicities[0], + *base_fes[1], + base_multiplicities[1]); + break; + } - default: - AssertThrow (false, ExcNotImplemented()); - } + case 3: + { + system_element = new FESystem(*base_fes[0], + base_multiplicities[0], + *base_fes[1], + base_multiplicities[1], + *base_fes[2], + base_multiplicities[2]); + break; + } - // now we don't need the - // list of base elements - // any more - for (unsigned int i=0; i(degree+1)) - // part should be left. - if (name.size() == 0 || name[0] != '(') - throw (std::string("Invalid first character in ") + name); - name.erase(0,1); - if (name[0] != 'Q') - { - const std::pair tmp - = Utilities::get_integer_at_position (name, 0); - name.erase(0, tmp.second+1); - return fe_name_map.find(name_part)->second->get(tmp.first); - } - else - { - unsigned int position = name.find('('); - const std::string quadrature_name(name, 0, position-1); - name.erase(0,position); - if (quadrature_name.compare("QGaussLobatto") == 0) - { - const std::pair tmp - = Utilities::get_integer_at_position (name, 0); - name.erase(0, tmp.second+1); + } + catch (...) + { + // ups, some exception + // was thrown. prevent a + // memory leak, and then + // pass on the exception + // to the caller + for (unsigned int i=0; i(degree+1)) + // part should be left. + if (name.size() == 0 || name[0] != '(') + throw (std::string("Invalid first character in ") + name); + name.erase(0,1); + if (name[0] != 'Q') + { + const std::pair tmp + = Utilities::get_integer_at_position (name, 0); + name.erase(0, tmp.second+1); + return fe_name_map.find(name_part)->second->get(tmp.first); + } + else + { + unsigned int position = name.find('('); + const std::string quadrature_name(name, 0, position-1); + name.erase(0,position); + if (quadrature_name.compare("QGaussLobatto") == 0) + { + const std::pair tmp + = Utilities::get_integer_at_position (name, 0); + name.erase(0, tmp.second+1); //TODO: Implement a get function taking Quadrature<1> in fe_tools.h. //return fe_name_map.find(name_part)->second->get(QGaussLobatto<1>(tmp.first)); - AssertThrow (false, ExcNotImplemented()); - } - else - { - AssertThrow (false,ExcNotImplemented()); - } - } - } + AssertThrow (false, ExcNotImplemented()); + } + else + { + AssertThrow (false,ExcNotImplemented()); + } + } + } - // hm, if we have come thus far, we - // didn't know what to do with the - // string we got. so do as the docs - // say: raise an exception - AssertThrow (false, FETools::ExcInvalidFEName(name)); + // hm, if we have come thus far, we + // didn't know what to do with the + // string we got. so do as the docs + // say: raise an exception + AssertThrow (false, ExcInvalidFEName(name)); - // make some compilers happy that - // do not realize that we can't get - // here after throwing - return 0; + // make some compilers happy that + // do not realize that we can't get + // here after throwing + return 0; + } } } -} -template -FiniteElement * -FETools::get_fe_from_name (const std::string ¶meter_name) -{ - // Create a version of the name - // string where all template - // parameters are eliminated. - std::string name = parameter_name; - for (unsigned int pos1 = name.find('<'); - pos1 < name.size(); - pos1 = name.find('<')) - { + template + FiniteElement * + get_fe_from_name (const std::string ¶meter_name) + { + // Create a version of the name + // string where all template + // parameters are eliminated. + std::string name = parameter_name; + for (unsigned int pos1 = name.find('<'); + pos1 < name.size(); + pos1 = name.find('<')) + { - const unsigned int pos2 = name.find('>'); - // If there is only a single - // character between those two, - // it should be 'd' or the number - // representing the dimension. - if (pos2-pos1 == 2) - { - const char dimchar = '0' + dim; - if (name.at(pos1+1) != 'd') - Assert (name.at(pos1+1) == dimchar, - ExcInvalidFEDimension(name.at(pos1+1), dim)); - } - else - Assert(pos2-pos1 == 4, ExcInvalidFEName(name)); + const unsigned int pos2 = name.find('>'); + // If there is only a single + // character between those two, + // it should be 'd' or the number + // representing the dimension. + if (pos2-pos1 == 2) + { + const char dimchar = '0' + dim; + if (name.at(pos1+1) != 'd') + Assert (name.at(pos1+1) == dimchar, + ExcInvalidFEDimension(name.at(pos1+1), dim)); + } + else + Assert(pos2-pos1 == 4, ExcInvalidFEName(name)); - // If pos1==pos2, then we are - // probably at the end of the - // string - if (pos2 != pos1) - name.erase(pos1, pos2-pos1+1); - } - // Replace all occurences of "^dim" - // by "^d" to be handled by the - // next loop - for (unsigned int pos = name.find("^dim"); - pos < name.size(); - pos = name.find("^dim")) - name.erase(pos+2, 2); - - // Replace all occurences of "^d" - // by using the actual dimension - for (unsigned int pos = name.find("^d"); - pos < name.size(); - pos = name.find("^d")) - name.at(pos+1) = '0' + dim; - - try - { - FiniteElement *fe = internal::get_fe_from_name (name); - - // Make sure the auxiliary function - // ate up all characters of the name. - AssertThrow (name.size() == 0, - ExcInvalidFEName(parameter_name - + std::string(" extra characters after " - "end of name"))); - return fe; - } - catch (const std::string &errline) - { - AssertThrow(false, ExcInvalidFEName(parameter_name - + std::string(" at ") - + errline)); - return 0; - } -} + // If pos1==pos2, then we are + // probably at the end of the + // string + if (pos2 != pos1) + name.erase(pos1, pos2-pos1+1); + } + // Replace all occurences of "^dim" + // by "^d" to be handled by the + // next loop + for (unsigned int pos = name.find("^dim"); + pos < name.size(); + pos = name.find("^dim")) + name.erase(pos+2, 2); + + // Replace all occurences of "^d" + // by using the actual dimension + for (unsigned int pos = name.find("^d"); + pos < name.size(); + pos = name.find("^d")) + name.at(pos+1) = '0' + dim; + + try + { + FiniteElement *fe = internal::get_fe_from_name (name); + + // Make sure the auxiliary function + // ate up all characters of the name. + AssertThrow (name.size() == 0, + ExcInvalidFEName(parameter_name + + std::string(" extra characters after " + "end of name"))); + return fe; + } + catch (const std::string &errline) + { + AssertThrow(false, ExcInvalidFEName(parameter_name + + std::string(" at ") + + errline)); + return 0; + } + } // template // FiniteElement * -// FETools::get_fe_from_name (const std::string ¶meter_name) +// get_fe_from_name (const std::string ¶meter_name) // { // return internal::get_fe_from_name(parameter_name); // } -template -void -FETools:: -compute_projection_from_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &lhs_quadrature, - const Quadrature &rhs_quadrature, - FullMatrix &X) -{ - Assert (fe.n_components() == 1, ExcNotImplemented()); - - // first build the matrices M and Q - // described in the documentation - FullMatrix M (fe.dofs_per_cell, fe.dofs_per_cell); - FullMatrix Q (fe.dofs_per_cell, rhs_quadrature.size()); - - for (unsigned int i=0; i M_inverse (fe.dofs_per_cell, fe.dofs_per_cell); - M_inverse.invert (M); - - // finally compute the result - X.reinit (fe.dofs_per_cell, rhs_quadrature.size()); - M_inverse.mmult (X, Q); - - Assert (X.m() == fe.dofs_per_cell, ExcInternalError()); - Assert (X.n() == rhs_quadrature.size(), ExcInternalError()); -} + template + void + compute_projection_from_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature, + FullMatrix &X) + { + Assert (fe.n_components() == 1, ExcNotImplemented()); + // first build the matrices M and Q + // described in the documentation + FullMatrix M (fe.dofs_per_cell, fe.dofs_per_cell); + FullMatrix Q (fe.dofs_per_cell, rhs_quadrature.size()); -template -void -FETools:: -compute_interpolation_to_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &quadrature, - FullMatrix &I_q) -{ - Assert (fe.n_components() == 1, ExcNotImplemented()); - Assert (I_q.m() == quadrature.size(), - ExcMessage ("Wrong matrix size")); - Assert (I_q.n() == fe.dofs_per_cell, ExcMessage ("Wrong matrix size")); + for (unsigned int i=0; i M_inverse (fe.dofs_per_cell, fe.dofs_per_cell); + M_inverse.invert (M); + // finally compute the result + X.reinit (fe.dofs_per_cell, rhs_quadrature.size()); + M_inverse.mmult (X, Q); -template -void -FETools::compute_projection_from_quadrature_points( - const FullMatrix &projection_matrix, - const std::vector< Tensor<1, dim > > &vector_of_tensors_at_qp, - std::vector< Tensor<1, dim > > &vector_of_tensors_at_nodes) -{ + Assert (X.m() == fe.dofs_per_cell, ExcInternalError()); + Assert (X.n() == rhs_quadrature.size(), ExcInternalError()); + } + + + + template + void + compute_interpolation_to_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &quadrature, + FullMatrix &I_q) + { + Assert (fe.n_components() == 1, ExcNotImplemented()); + Assert (I_q.m() == quadrature.size(), + ExcMessage ("Wrong matrix size")); + Assert (I_q.n() == fe.dofs_per_cell, ExcMessage ("Wrong matrix size")); + + for (unsigned int q=0; q + void + compute_projection_from_quadrature_points( + const FullMatrix &projection_matrix, + const std::vector< Tensor<1, dim > > &vector_of_tensors_at_qp, + std::vector< Tensor<1, dim > > &vector_of_tensors_at_nodes) + { - // check that the number columns of the projection_matrix - // matches the size of the vector_of_tensors_at_qp + // check that the number columns of the projection_matrix + // matches the size of the vector_of_tensors_at_qp Assert(projection_matrix.n_cols() == vector_of_tensors_at_qp.size(), - ExcDimensionMismatch(projection_matrix.n_cols(), - vector_of_tensors_at_qp.size())); + ExcDimensionMismatch(projection_matrix.n_cols(), + vector_of_tensors_at_qp.size())); - // check that the number rows of the projection_matrix - // matches the size of the vector_of_tensors_at_nodes + // check that the number rows of the projection_matrix + // matches the size of the vector_of_tensors_at_nodes Assert(projection_matrix.n_rows() == vector_of_tensors_at_nodes.size(), - ExcDimensionMismatch(projection_matrix.n_rows(), - vector_of_tensors_at_nodes.size())); + ExcDimensionMismatch(projection_matrix.n_rows(), + vector_of_tensors_at_nodes.size())); - // number of support points (nodes) to project to + // number of support points (nodes) to project to const unsigned int n_support_points = projection_matrix.n_rows(); - // number of quadrature points to project from + // number of quadrature points to project from const unsigned int n_quad_points = projection_matrix.n_cols(); - // component projected to the nodes + // component projected to the nodes Vector component_at_node(n_support_points); - // component at the quadrature point + // component at the quadrature point Vector component_at_qp(n_quad_points); for (unsigned int ii = 0; ii < dim; ++ii) { - component_at_qp = 0; + component_at_qp = 0; - // populate the vector of components at the qps - // from vector_of_tensors_at_qp - // vector_of_tensors_at_qp data is in form: - // columns: 0, 1, ..., dim - // rows: 0,1,...., n_quad_points - // so extract the ii'th column of vector_of_tensors_at_qp - for (unsigned int q = 0; q < n_quad_points; ++q) { - component_at_qp(q) = vector_of_tensors_at_qp[q][ii]; - } + // populate the vector of components at the qps + // from vector_of_tensors_at_qp + // vector_of_tensors_at_qp data is in form: + // columns: 0, 1, ..., dim + // rows: 0,1,...., n_quad_points + // so extract the ii'th column of vector_of_tensors_at_qp + for (unsigned int q = 0; q < n_quad_points; ++q) { + component_at_qp(q) = vector_of_tensors_at_qp[q][ii]; + } - // project from the qps -> nodes - // component_at_node = projection_matrix_u * component_at_qp - projection_matrix.vmult(component_at_node, component_at_qp); + // project from the qps -> nodes + // component_at_node = projection_matrix_u * component_at_qp + projection_matrix.vmult(component_at_node, component_at_qp); - // rewrite the projection of the components - // back into the vector of tensors - for (unsigned int nn =0; nn -void -FETools::compute_projection_from_quadrature_points( - const FullMatrix &projection_matrix, - const std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_qp, - std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_nodes) -{ + template + void + compute_projection_from_quadrature_points( + const FullMatrix &projection_matrix, + const std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_qp, + std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_nodes) + { - // check that the number columns of the projection_matrix - // matches the size of the vector_of_tensors_at_qp + // check that the number columns of the projection_matrix + // matches the size of the vector_of_tensors_at_qp Assert(projection_matrix.n_cols() == vector_of_tensors_at_qp.size(), - ExcDimensionMismatch(projection_matrix.n_cols(), - vector_of_tensors_at_qp.size())); + ExcDimensionMismatch(projection_matrix.n_cols(), + vector_of_tensors_at_qp.size())); - // check that the number rows of the projection_matrix - // matches the size of the vector_of_tensors_at_nodes + // check that the number rows of the projection_matrix + // matches the size of the vector_of_tensors_at_nodes Assert(projection_matrix.n_rows() == vector_of_tensors_at_nodes.size(), - ExcDimensionMismatch(projection_matrix.n_rows(), - vector_of_tensors_at_nodes.size())); + ExcDimensionMismatch(projection_matrix.n_rows(), + vector_of_tensors_at_nodes.size())); - // number of support points (nodes) + // number of support points (nodes) const unsigned int n_support_points = projection_matrix.n_rows(); - // number of quadrature points to project from + // number of quadrature points to project from const unsigned int n_quad_points = projection_matrix.n_cols(); - // number of unique entries in a symmetric second-order tensor + // number of unique entries in a symmetric second-order tensor const unsigned int n_independent_components = - SymmetricTensor<2, dim >::n_independent_components; + SymmetricTensor<2, dim >::n_independent_components; - // component projected to the nodes + // component projected to the nodes Vector component_at_node(n_support_points); - // component at the quadrature point + // component at the quadrature point Vector component_at_qp(n_quad_points); - // loop over the number of unique dimensions of the tensor + // loop over the number of unique dimensions of the tensor for (unsigned int ii = 0; ii < n_independent_components; ++ii) { - component_at_qp = 0; - - // row-column entry of tensor corresponding the unrolled index - TableIndices<2> row_column_index = SymmetricTensor< 2, dim >::unrolled_to_component_indices(ii); - const unsigned int row = row_column_index[0]; - const unsigned int column = row_column_index[1]; - - // populate the vector of components at the qps - // from vector_of_tensors_at_qp - // vector_of_tensors_at_qp is in form: - // columns: 0, 1, ..., n_independent_components - // rows: 0,1,...., n_quad_points - // so extract the ii'th column of vector_of_tensors_at_qp - for (unsigned int q = 0; q < n_quad_points; ++q) { - component_at_qp(q) = (vector_of_tensors_at_qp[q])[row][column]; - } + component_at_qp = 0; + + // row-column entry of tensor corresponding the unrolled index + TableIndices<2> row_column_index = SymmetricTensor< 2, dim >::unrolled_to_component_indices(ii); + const unsigned int row = row_column_index[0]; + const unsigned int column = row_column_index[1]; + + // populate the vector of components at the qps + // from vector_of_tensors_at_qp + // vector_of_tensors_at_qp is in form: + // columns: 0, 1, ..., n_independent_components + // rows: 0,1,...., n_quad_points + // so extract the ii'th column of vector_of_tensors_at_qp + for (unsigned int q = 0; q < n_quad_points; ++q) { + component_at_qp(q) = (vector_of_tensors_at_qp[q])[row][column]; + } - // project from the qps -> nodes - // component_at_node = projection_matrix_u * component_at_qp - projection_matrix.vmult(component_at_node, component_at_qp); + // project from the qps -> nodes + // component_at_node = projection_matrix_u * component_at_qp + projection_matrix.vmult(component_at_node, component_at_qp); - // rewrite the projection of the components back into the vector of tensors - for (unsigned int nn =0; nn -void -FETools:: -compute_projection_from_face_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &lhs_quadrature, - const Quadrature &rhs_quadrature, - const typename DoFHandler::active_cell_iterator & cell, - unsigned int face, - FullMatrix &X) -{ - Assert (fe.n_components() == 1, ExcNotImplemented()); - Assert (lhs_quadrature.size () > fe.degree, ExcNotGreaterThan (lhs_quadrature.size (), fe.degree)); + template + void + compute_projection_from_face_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature, + const typename DoFHandler::active_cell_iterator & cell, + unsigned int face, + FullMatrix &X) + { + Assert (fe.n_components() == 1, ExcNotImplemented()); + Assert (lhs_quadrature.size () > fe.degree, ExcNotGreaterThan (lhs_quadrature.size (), fe.degree)); - // build the matrices M and Q - // described in the documentation - FullMatrix M (fe.dofs_per_cell, fe.dofs_per_cell); - FullMatrix Q (fe.dofs_per_cell, rhs_quadrature.size()); + // build the matrices M and Q + // described in the documentation + FullMatrix M (fe.dofs_per_cell, fe.dofs_per_cell); + FullMatrix Q (fe.dofs_per_cell, rhs_quadrature.size()); - { - // need an FEFaceValues object to evaluate shape function - // values on the specified face. - FEFaceValues fe_face_values (fe, lhs_quadrature, update_values); - fe_face_values.reinit (cell, face); // setup shape_value on this face. + { + // need an FEFaceValues object to evaluate shape function + // values on the specified face. + FEFaceValues fe_face_values (fe, lhs_quadrature, update_values); + fe_face_values.reinit (cell, face); // setup shape_value on this face. + + for (unsigned int i=0; i fe_face_values (fe, rhs_quadrature, update_values); + fe_face_values.reinit (cell, face); // setup shape_value on this face. + + for (unsigned int i=0; i M_inverse (fe.dofs_per_cell, fe.dofs_per_cell); + M_inverse.invert (M); - { - FEFaceValues fe_face_values (fe, rhs_quadrature, update_values); - fe_face_values.reinit (cell, face); // setup shape_value on this face. + // finally compute the result + X.reinit (fe.dofs_per_cell, rhs_quadrature.size()); + M_inverse.mmult (X, Q); - for (unsigned int i=0; i M_inverse (fe.dofs_per_cell, fe.dofs_per_cell); - M_inverse.invert (M); - // finally compute the result - X.reinit (fe.dofs_per_cell, rhs_quadrature.size()); - M_inverse.mmult (X, Q); - - Assert (X.m() == fe.dofs_per_cell, ExcInternalError()); - Assert (X.n() == rhs_quadrature.size(), ExcInternalError()); } - /*-------------- Explicit Instantiations -------------------------------*/ - -template class FETools::FEFactoryBase; - -template -void FETools::compute_node_matrix( - FullMatrix&, - const FiniteElement&); - -template -void FETools::compute_component_wise( - const FiniteElement& element, - std::vector&, std::vector >&); -template -void FETools::compute_block_renumbering ( - const FiniteElement& element, - std::vector&, std::vector&_indices, bool); -template -void FETools::get_interpolation_matrix -(const FiniteElement &, - const FiniteElement &, - FullMatrix &); -template -void FETools::get_back_interpolation_matrix -(const FiniteElement &, - const FiniteElement &, - FullMatrix &); -template -void FETools::get_interpolation_difference_matrix -(const FiniteElement &, - const FiniteElement &, - FullMatrix &); -template -void FETools::get_interpolation_matrix -(const FiniteElement &, - const FiniteElement &, - FullMatrix &); -template -void FETools::get_back_interpolation_matrix -(const FiniteElement &, - const FiniteElement &, - FullMatrix &); -template -void FETools::get_interpolation_difference_matrix -(const FiniteElement &, - const FiniteElement &, - FullMatrix &); - -template -void FETools::get_projection_matrix -(const FiniteElement &, - const FiniteElement &, - FullMatrix &); - -template -void FETools::compute_embedding_matrices -(const FiniteElement &, std::vector > >&,bool); - -template -void FETools::compute_face_embedding_matrices -(const FiniteElement &, FullMatrix (&matrices)[GeometryInfo::max_children_per_face], - unsigned int, unsigned int); - -template -void FETools::compute_projection_matrices -(const FiniteElement &, std::vector > >&, bool); - -template -void FETools::interpolate -(const DoFHandler &, const Vector &, - const DoFHandler &, Vector &); -template -void FETools::interpolate -(const DoFHandler &, const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); +namespace FETools +{ + template class FEFactoryBase; + + template + void compute_node_matrix( + FullMatrix&, + const FiniteElement&); + + template + void compute_component_wise( + const FiniteElement& element, + std::vector&, std::vector >&); + template + void compute_block_renumbering ( + const FiniteElement& element, + std::vector&, std::vector&_indices, bool); + template + void get_interpolation_matrix + (const FiniteElement &, + const FiniteElement &, + FullMatrix &); + template + void get_back_interpolation_matrix + (const FiniteElement &, + const FiniteElement &, + FullMatrix &); + template + void get_interpolation_difference_matrix + (const FiniteElement &, + const FiniteElement &, + FullMatrix &); + template + void get_interpolation_matrix + (const FiniteElement &, + const FiniteElement &, + FullMatrix &); + template + void get_back_interpolation_matrix + (const FiniteElement &, + const FiniteElement &, + FullMatrix &); + template + void get_interpolation_difference_matrix + (const FiniteElement &, + const FiniteElement &, + FullMatrix &); + + template + void get_projection_matrix + (const FiniteElement &, + const FiniteElement &, + FullMatrix &); + + template + void compute_embedding_matrices + (const FiniteElement &, std::vector > >&,bool); + + template + void compute_face_embedding_matrices + (const FiniteElement &, FullMatrix (&matrices)[GeometryInfo::max_children_per_face], + unsigned int, unsigned int); + + template + void compute_projection_matrices + (const FiniteElement &, std::vector > >&, bool); + + template + void interpolate + (const DoFHandler &, const Vector &, + const DoFHandler &, Vector &); + template + void interpolate + (const DoFHandler &, const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); #if deal_II_dimension != 3 -template -void FETools::compute_block_renumbering ( - const FiniteElement& element, - std::vector&, std::vector&_indices, bool); -template -void FETools::interpolate -(const DoFHandler &, const Vector &, - const DoFHandler &, Vector &); -template -void FETools::interpolate -(const DoFHandler &, const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); + template + void compute_block_renumbering ( + const FiniteElement& element, + std::vector&, std::vector&_indices, bool); + template + void interpolate + (const DoFHandler &, const Vector &, + const DoFHandler &, Vector &); + template + void interpolate + (const DoFHandler &, const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); #endif -template -void FETools::back_interpolate -(const DoFHandler &, const Vector &, - const FiniteElement &, Vector &); -template -void FETools::back_interpolate -(const DoFHandler &, const ConstraintMatrix &, - const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); -template -void FETools::interpolation_difference -(const DoFHandler &, const Vector &, - const FiniteElement &, Vector &); -template -void FETools::interpolation_difference -(const DoFHandler &, const ConstraintMatrix &, - const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); -template -void FETools::project_dg -(const DoFHandler &, const Vector &, - const DoFHandler &, Vector &); -template -void FETools::extrapolate -(const DoFHandler &, const Vector &, - const DoFHandler &, Vector &); -template -void FETools::extrapolate -(const DoFHandler &, const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); - - -template -void FETools::interpolate -(const DoFHandler &, const Vector &, - const DoFHandler &, Vector &); -template -void FETools::interpolate -(const DoFHandler &, const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); + template + void back_interpolate + (const DoFHandler &, const Vector &, + const FiniteElement &, Vector &); + template + void back_interpolate + (const DoFHandler &, const ConstraintMatrix &, + const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); + template + void interpolation_difference + (const DoFHandler &, const Vector &, + const FiniteElement &, Vector &); + template + void interpolation_difference + (const DoFHandler &, const ConstraintMatrix &, + const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); + template + void project_dg + (const DoFHandler &, const Vector &, + const DoFHandler &, Vector &); + template + void extrapolate + (const DoFHandler &, const Vector &, + const DoFHandler &, Vector &); + template + void extrapolate + (const DoFHandler &, const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); + + + template + void interpolate + (const DoFHandler &, const Vector &, + const DoFHandler &, Vector &); + template + void interpolate + (const DoFHandler &, const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); #if deal_II_dimension != 3 -template -void FETools::interpolate -(const DoFHandler &, const Vector &, - const DoFHandler &, Vector &); -template -void FETools::interpolate -(const DoFHandler &, const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); + template + void interpolate + (const DoFHandler &, const Vector &, + const DoFHandler &, Vector &); + template + void interpolate + (const DoFHandler &, const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); #endif -template -void FETools::back_interpolate -(const DoFHandler &, const Vector &, - const FiniteElement &, Vector &); -template -void FETools::back_interpolate -(const DoFHandler &, const ConstraintMatrix &, - const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); -template -void FETools::interpolation_difference -(const DoFHandler &, const Vector &, - const FiniteElement &, Vector &); -template -void FETools::interpolation_difference -(const DoFHandler &, const ConstraintMatrix &, - const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); -template -void FETools::project_dg -(const DoFHandler &, const Vector &, - const DoFHandler &, Vector &); -template -void FETools::extrapolate -(const DoFHandler &, const Vector &, - const DoFHandler &, Vector &); -template -void FETools::extrapolate -(const DoFHandler &, const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); + template + void back_interpolate + (const DoFHandler &, const Vector &, + const FiniteElement &, Vector &); + template + void back_interpolate + (const DoFHandler &, const ConstraintMatrix &, + const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); + template + void interpolation_difference + (const DoFHandler &, const Vector &, + const FiniteElement &, Vector &); + template + void interpolation_difference + (const DoFHandler &, const ConstraintMatrix &, + const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); + template + void project_dg + (const DoFHandler &, const Vector &, + const DoFHandler &, Vector &); + template + void extrapolate + (const DoFHandler &, const Vector &, + const DoFHandler &, Vector &); + template + void extrapolate + (const DoFHandler &, const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); #ifdef DEAL_II_USE_TRILINOS -template -void FETools::interpolate -(const DoFHandler &, const TrilinosWrappers::Vector &, - const DoFHandler &, TrilinosWrappers::Vector &); -template -void FETools::interpolate -(const DoFHandler &, const TrilinosWrappers::Vector &, - const DoFHandler &, const ConstraintMatrix &, - TrilinosWrappers::Vector &); -template -void FETools::back_interpolate -(const DoFHandler &, const TrilinosWrappers::Vector &, - const FiniteElement &, TrilinosWrappers::Vector &); -template -void FETools::back_interpolate -(const DoFHandler &, const ConstraintMatrix &, - const TrilinosWrappers::Vector &, - const DoFHandler &, const ConstraintMatrix &, - TrilinosWrappers::Vector &); -template -void FETools::interpolation_difference -(const DoFHandler &, const TrilinosWrappers::Vector &, - const FiniteElement &, TrilinosWrappers::Vector &); -template -void FETools::interpolation_difference -(const DoFHandler &, const ConstraintMatrix &, - const TrilinosWrappers::Vector &, - const DoFHandler &, const ConstraintMatrix &, - TrilinosWrappers::Vector &); -template -void FETools::project_dg -(const DoFHandler &, const TrilinosWrappers::Vector &, - const DoFHandler &, TrilinosWrappers::Vector &); -template -void FETools::extrapolate -(const DoFHandler &, const TrilinosWrappers::Vector &, - const DoFHandler &, TrilinosWrappers::Vector &); -template -void FETools::extrapolate -(const DoFHandler &, const TrilinosWrappers::Vector &, - const DoFHandler &, const ConstraintMatrix &, - TrilinosWrappers::Vector &); + template + void interpolate + (const DoFHandler &, const TrilinosWrappers::Vector &, + const DoFHandler &, TrilinosWrappers::Vector &); + template + void interpolate + (const DoFHandler &, const TrilinosWrappers::Vector &, + const DoFHandler &, const ConstraintMatrix &, + TrilinosWrappers::Vector &); + template + void back_interpolate + (const DoFHandler &, const TrilinosWrappers::Vector &, + const FiniteElement &, TrilinosWrappers::Vector &); + template + void back_interpolate + (const DoFHandler &, const ConstraintMatrix &, + const TrilinosWrappers::Vector &, + const DoFHandler &, const ConstraintMatrix &, + TrilinosWrappers::Vector &); + template + void interpolation_difference + (const DoFHandler &, const TrilinosWrappers::Vector &, + const FiniteElement &, TrilinosWrappers::Vector &); + template + void interpolation_difference + (const DoFHandler &, const ConstraintMatrix &, + const TrilinosWrappers::Vector &, + const DoFHandler &, const ConstraintMatrix &, + TrilinosWrappers::Vector &); + template + void project_dg + (const DoFHandler &, const TrilinosWrappers::Vector &, + const DoFHandler &, TrilinosWrappers::Vector &); + template + void extrapolate + (const DoFHandler &, const TrilinosWrappers::Vector &, + const DoFHandler &, TrilinosWrappers::Vector &); + template + void extrapolate + (const DoFHandler &, const TrilinosWrappers::Vector &, + const DoFHandler &, const ConstraintMatrix &, + TrilinosWrappers::Vector &); #endif -template -void FETools::interpolate -(const DoFHandler &, const BlockVector &, - const DoFHandler &, BlockVector &); -template -void FETools::interpolate -(const DoFHandler &, const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); + template + void interpolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, BlockVector &); + template + void interpolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); #if deal_II_dimension != 3 -template -void FETools::interpolate -(const DoFHandler &, const BlockVector &, - const DoFHandler &, BlockVector &); -template -void FETools::interpolate -(const DoFHandler &, const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); + template + void interpolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, BlockVector &); + template + void interpolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); #endif -template -void FETools::back_interpolate -(const DoFHandler &, const BlockVector &, - const FiniteElement &, BlockVector &); -template -void FETools::back_interpolate -(const DoFHandler &, const ConstraintMatrix &, - const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); -template -void FETools::interpolation_difference -(const DoFHandler &, const BlockVector &, - const FiniteElement &, BlockVector &); -template -void FETools::interpolation_difference -(const DoFHandler &, const ConstraintMatrix &, - const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); -template -void FETools::project_dg -(const DoFHandler &, const BlockVector &, - const DoFHandler &, BlockVector &); -template -void FETools::extrapolate -(const DoFHandler &, const BlockVector &, - const DoFHandler &, BlockVector &); -template -void FETools::extrapolate -(const DoFHandler &, const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); -template -void FETools::extrapolate -(const DoFHandler &, const BlockVector &, - const DoFHandler &, Vector &); -template -void FETools::extrapolate -(const DoFHandler &, const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); - - -template -void FETools::interpolate -(const DoFHandler &, const BlockVector &, - const DoFHandler &, BlockVector &); -template -void FETools::interpolate -(const DoFHandler &, const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); + template + void back_interpolate + (const DoFHandler &, const BlockVector &, + const FiniteElement &, BlockVector &); + template + void back_interpolate + (const DoFHandler &, const ConstraintMatrix &, + const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); + template + void interpolation_difference + (const DoFHandler &, const BlockVector &, + const FiniteElement &, BlockVector &); + template + void interpolation_difference + (const DoFHandler &, const ConstraintMatrix &, + const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); + template + void project_dg + (const DoFHandler &, const BlockVector &, + const DoFHandler &, BlockVector &); + template + void extrapolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, BlockVector &); + template + void extrapolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); + template + void extrapolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, Vector &); + template + void extrapolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); + + + template + void interpolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, BlockVector &); + template + void interpolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); #if deal_II_dimension != 3 -template -void FETools::interpolate -(const DoFHandler &, const BlockVector &, - const DoFHandler &, BlockVector &); -template -void FETools::interpolate -(const DoFHandler &, const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); + template + void interpolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, BlockVector &); + template + void interpolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); #endif -template -void FETools::back_interpolate -(const DoFHandler &, const BlockVector &, - const FiniteElement &, BlockVector &); -template -void FETools::back_interpolate -(const DoFHandler &, const ConstraintMatrix &, - const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); -template -void FETools::interpolation_difference -(const DoFHandler &, const BlockVector &, - const FiniteElement &, BlockVector &); -template -void FETools::interpolation_difference -(const DoFHandler &, const ConstraintMatrix &, - const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); -template -void FETools::project_dg -(const DoFHandler &, const BlockVector &, - const DoFHandler &, BlockVector &); -template -void FETools::extrapolate -(const DoFHandler &, const BlockVector &, - const DoFHandler &, BlockVector &); -template -void FETools::extrapolate -(const DoFHandler &, const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); -template -void FETools::extrapolate -(const DoFHandler &, const BlockVector &, - const DoFHandler &, Vector &); -template -void FETools::extrapolate -(const DoFHandler &, const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); + template + void back_interpolate + (const DoFHandler &, const BlockVector &, + const FiniteElement &, BlockVector &); + template + void back_interpolate + (const DoFHandler &, const ConstraintMatrix &, + const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); + template + void interpolation_difference + (const DoFHandler &, const BlockVector &, + const FiniteElement &, BlockVector &); + template + void interpolation_difference + (const DoFHandler &, const ConstraintMatrix &, + const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); + template + void project_dg + (const DoFHandler &, const BlockVector &, + const DoFHandler &, BlockVector &); + template + void extrapolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, BlockVector &); + template + void extrapolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); + template + void extrapolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, Vector &); + template + void extrapolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); #ifdef DEAL_II_USE_TRILINOS -template -void FETools::interpolate -(const DoFHandler &, const TrilinosWrappers::BlockVector &, - const DoFHandler &, TrilinosWrappers::BlockVector &); -template -void FETools::interpolate -(const DoFHandler &, const TrilinosWrappers::BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - TrilinosWrappers::BlockVector &); -template -void FETools::back_interpolate -(const DoFHandler &, const TrilinosWrappers::BlockVector &, - const FiniteElement &, TrilinosWrappers::BlockVector &); -template -void FETools::back_interpolate -(const DoFHandler &, const ConstraintMatrix &, - const TrilinosWrappers::BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - TrilinosWrappers::BlockVector &); -template -void FETools::interpolation_difference -(const DoFHandler &, const TrilinosWrappers::BlockVector &, - const FiniteElement &, TrilinosWrappers::BlockVector &); -template -void FETools::interpolation_difference -(const DoFHandler &, const ConstraintMatrix &, - const TrilinosWrappers::BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - TrilinosWrappers::BlockVector &); -template -void FETools::project_dg -(const DoFHandler &, const TrilinosWrappers::BlockVector &, - const DoFHandler &, TrilinosWrappers::BlockVector &); -template -void FETools::extrapolate -(const DoFHandler &, const TrilinosWrappers::BlockVector &, - const DoFHandler &, TrilinosWrappers::BlockVector &); -template -void FETools::extrapolate -(const DoFHandler &, const TrilinosWrappers::BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - TrilinosWrappers::BlockVector &); + template + void interpolate + (const DoFHandler &, const TrilinosWrappers::BlockVector &, + const DoFHandler &, TrilinosWrappers::BlockVector &); + template + void interpolate + (const DoFHandler &, const TrilinosWrappers::BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + TrilinosWrappers::BlockVector &); + template + void back_interpolate + (const DoFHandler &, const TrilinosWrappers::BlockVector &, + const FiniteElement &, TrilinosWrappers::BlockVector &); + template + void back_interpolate + (const DoFHandler &, const ConstraintMatrix &, + const TrilinosWrappers::BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + TrilinosWrappers::BlockVector &); + template + void interpolation_difference + (const DoFHandler &, const TrilinosWrappers::BlockVector &, + const FiniteElement &, TrilinosWrappers::BlockVector &); + template + void interpolation_difference + (const DoFHandler &, const ConstraintMatrix &, + const TrilinosWrappers::BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + TrilinosWrappers::BlockVector &); + template + void project_dg + (const DoFHandler &, const TrilinosWrappers::BlockVector &, + const DoFHandler &, TrilinosWrappers::BlockVector &); + template + void extrapolate + (const DoFHandler &, const TrilinosWrappers::BlockVector &, + const DoFHandler &, TrilinosWrappers::BlockVector &); + template + void extrapolate + (const DoFHandler &, const TrilinosWrappers::BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + TrilinosWrappers::BlockVector &); #endif -template -void FETools::interpolate -(const hp::DoFHandler &, const Vector &, - const hp::DoFHandler &, Vector &); -template -void FETools::interpolate -(const hp::DoFHandler &, const Vector &, - const hp::DoFHandler &, const ConstraintMatrix &, - Vector &); -template -void FETools::interpolate -(const hp::DoFHandler &, const Vector &, - const hp::DoFHandler &, Vector &); -template -void FETools::interpolate -(const hp::DoFHandler &, const Vector &, - const hp::DoFHandler &, const ConstraintMatrix &, - Vector &); - - -template FiniteElement * -FETools::get_fe_from_name (const std::string &); - - -template -void FETools::add_fe_name( - const std::string& name, - const FEFactoryBase* factory); - -template -void -FETools:: -compute_projection_from_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &lhs_quadrature, - const Quadrature &rhs_quadrature, - FullMatrix &X); + template + void interpolate + (const hp::DoFHandler &, const Vector &, + const hp::DoFHandler &, Vector &); + template + void interpolate + (const hp::DoFHandler &, const Vector &, + const hp::DoFHandler &, const ConstraintMatrix &, + Vector &); + template + void interpolate + (const hp::DoFHandler &, const Vector &, + const hp::DoFHandler &, Vector &); + template + void interpolate + (const hp::DoFHandler &, const Vector &, + const hp::DoFHandler &, const ConstraintMatrix &, + Vector &); + + + template FiniteElement * + get_fe_from_name (const std::string &); + + + template + void add_fe_name( + const std::string& name, + const FEFactoryBase* factory); + + template + void + compute_projection_from_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature, + FullMatrix &X); -template -void -FETools:: -compute_projection_from_quadrature_points( - const FullMatrix &projection_matrix, - const std::vector< Tensor<1, deal_II_dimension > > &vector_of_tensors_at_qp, - std::vector< Tensor<1, deal_II_dimension > > &vector_of_tensors_at_nodes); + template + void + compute_projection_from_quadrature_points( + const FullMatrix &projection_matrix, + const std::vector< Tensor<1, deal_II_dimension > > &vector_of_tensors_at_qp, + std::vector< Tensor<1, deal_II_dimension > > &vector_of_tensors_at_nodes); -template -void -FETools::compute_projection_from_quadrature_points( - const FullMatrix &projection_matrix, - const std::vector > &vector_of_tensors_at_qp, - std::vector > &vector_of_tensors_at_nodes); + template + void + compute_projection_from_quadrature_points( + const FullMatrix &projection_matrix, + const std::vector > &vector_of_tensors_at_qp, + std::vector > &vector_of_tensors_at_nodes); -template -void -FETools:: -compute_interpolation_to_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &quadrature, - FullMatrix &I_q); + template + void + compute_interpolation_to_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &quadrature, + FullMatrix &I_q); #if deal_II_dimension != 1 -template -void -FETools:: -compute_projection_from_face_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &lhs_quadrature, - const Quadrature &rhs_quadrature, - const DoFHandler::active_cell_iterator & cell, - unsigned int face, - FullMatrix &X); + template + void + compute_projection_from_face_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature, + const DoFHandler::active_cell_iterator & cell, + unsigned int face, + FullMatrix &X); #endif +} + /*---------------------------- fe_tools.cc ---------------------------*/ diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 077d799027..fdfca2971f 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -281,6 +281,12 @@ through DoFHandler::get_tria() and DoFHandler::get_fe(), respectively.
    +
  1. Changed: FETools is now a namespace rather than a class with only + static member functions. +
    + (WB 2010/09/25) +

  2. +
  3. New: FE_DGVector implements discontinuous elements based on vector valued polynomial spaces.
    -- 2.39.5