/**
- * 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.
*
* @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 <int dim, int spacedim=dim>
- class FEFactoryBase
- {
- public:
- /**
- * Create a FiniteElement and
- * return a pointer to it.
- */
- virtual FiniteElement<dim,spacedim>*
- 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 FE>
- class FEFactory : public FEFactoryBase<FE::dimension,FE::dimension>
- {
- public:
- /**
- * Create a FiniteElement and
- * return a pointer to it.
- */
- virtual FiniteElement<FE::dimension,FE::dimension>*
- 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<int dim, int spacedim>
- static void compute_component_wise(
- const FiniteElement<dim,spacedim>& fe,
- std::vector<unsigned int>& renumbering,
- std::vector<std::vector<unsigned int> >& 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<int dim, int spacedim>
- static void compute_block_renumbering (
- const FiniteElement<dim,spacedim>& fe,
- std::vector<unsigned int>& renumbering,
- std::vector<unsigned int>& 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
- * <tt>(fe2.dofs_per_cell, fe1.dofs_per_cell)</tt>.
- *
- * 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 <int dim, typename number, int spacedim>
- static
- void
- get_interpolation_matrix(const FiniteElement<dim,spacedim> &fe1,
- const FiniteElement<dim,spacedim> &fe2,
- FullMatrix<number> &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
- * <tt>(fe1.dofs_per_cell, fe1.dofs_per_cell)</tt>.
- *
- * 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 <int dim, typename number, int spacedim>
- static
- void
- get_back_interpolation_matrix(const FiniteElement<dim,spacedim> &fe1,
- const FiniteElement<dim,spacedim> &fe2,
- FullMatrix<number> &interpolation_matrix);
-
- /**
- * Gives the unit matrix minus the
- * back interpolation matrix.
- * The @p difference_matrix
- * needs to be of size
- * <tt>(fe1.dofs_per_cell, fe1.dofs_per_cell)</tt>.
- *
- * 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 <int dim, typename number, int spacedim>
- static
- void
- get_interpolation_difference_matrix(const FiniteElement<dim,spacedim> &fe1,
- const FiniteElement<dim,spacedim> &fe2,
- FullMatrix<number> &difference_matrix);
-
- /**
- * Compute the local
- * $L^2$-projection matrix from
- * fe1 to fe2.
- */
- template <int dim, typename number, int spacedim>
- static void get_projection_matrix(const FiniteElement<dim,spacedim> &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 <int dim, int spacedim=dim>
+ class FEFactoryBase
+ {
+ public:
+ /**
+ * Create a FiniteElement and
+ * return a pointer to it.
+ */
+ virtual FiniteElement<dim,spacedim>*
+ 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 FE>
+ class FEFactory : public FEFactoryBase<FE::dimension,FE::dimension>
+ {
+ public:
+ /**
+ * Create a FiniteElement and
+ * return a pointer to it.
+ */
+ virtual FiniteElement<FE::dimension,FE::dimension>*
+ 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<int dim, int spacedim>
+ void compute_component_wise(
+ const FiniteElement<dim,spacedim>& fe,
+ std::vector<unsigned int>& renumbering,
+ std::vector<std::vector<unsigned int> >& 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<int dim, int spacedim>
+ void compute_block_renumbering (
+ const FiniteElement<dim,spacedim>& fe,
+ std::vector<unsigned int>& renumbering,
+ std::vector<unsigned int>& 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
+ * <tt>(fe2.dofs_per_cell, fe1.dofs_per_cell)</tt>.
+ *
+ * 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 <int dim, typename number, int spacedim>
+ void
+ get_interpolation_matrix(const FiniteElement<dim,spacedim> &fe1,
+ const FiniteElement<dim,spacedim> &fe2,
+ FullMatrix<number> &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
+ * <tt>(fe1.dofs_per_cell, fe1.dofs_per_cell)</tt>.
+ *
+ * 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 <int dim, typename number, int spacedim>
+ void
+ get_back_interpolation_matrix(const FiniteElement<dim,spacedim> &fe1,
+ const FiniteElement<dim,spacedim> &fe2,
+ FullMatrix<number> &interpolation_matrix);
+
+ /**
+ * Gives the unit matrix minus the
+ * back interpolation matrix.
+ * The @p difference_matrix
+ * needs to be of size
+ * <tt>(fe1.dofs_per_cell, fe1.dofs_per_cell)</tt>.
+ *
+ * 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 <int dim, typename number, int spacedim>
+ void
+ get_interpolation_difference_matrix(const FiniteElement<dim,spacedim> &fe1,
const FiniteElement<dim,spacedim> &fe2,
- FullMatrix<number> &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.
- * <ol>
- *
- * <li>Define the space of shape
- * functions using an arbitrary
- * basis <i>w<sub>j</sub></i> and
- * compute the matrix <i>M</i> of
- * node functionals
- * <i>N<sub>i</sub></i> applied
- * to these basis functions.
- *
- * <li>Compute the basis
- * <i>v<sub>j</sub></i> of the
- * finite element shape function
- * space by applying
- * <i>M<sup>-1</sup></i> to the
- * basis <i>w<sub>j</sub></i>.
- * </ol>
- *
- * @note The FiniteElement must
- * provide generalized support
- * points and and interpolation
- * functions.
- */
- template <int dim, int spacedim>
- static void compute_node_matrix(FullMatrix<double>& M,
- const FiniteElement<dim,spacedim>& 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
- * <code>matrices[refinement_case-1][child]</code>
- * includes the embedding (or prolongation)
- * matrix of child
- * <code>child</code> for the
- * RefinementCase
- * <code>refinement_case</code>. Here,
- * we use
- * <code>refinement_case-1</code>
- * instead of
- * <code>refinement_case</code>
- * 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<dim>::isotropic_refinement
- * vectors of FullMatrix
- * objects. Each vector
- * corresponds to one
- * RefinementCase @p
- * refinement_case and is of the
- * vector size
- * GeometryInfo<dim>::n_children(refinement_case). This
- * is the format used in
- * FiniteElement, where we want
- * to use this function mostly.
- *
- * @param isotropic_only Set
- * to <code>true</code> if you only
- * want to compute matrices for
- * isotropic refinement.
- */
- template <int dim, typename number, int spacedim>
- static void compute_embedding_matrices(const FiniteElement<dim,spacedim> &fe,
- std::vector<std::vector<FullMatrix<number> > >& 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
- * <i>GeometryInfo<dim>::subfaces_per_face
- * = 2<sup>dim-1</sup></i>
- * 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 <int dim, typename number, int spacedim>
- static void
- compute_face_embedding_matrices(const FiniteElement<dim,spacedim>& fe,
- FullMatrix<number> (&matrices)[GeometryInfo<dim>::max_children_per_face],
- const unsigned int face_coarse,
- const unsigned int face_fine);
-
- /**
- * For all possible (isotropic
- * and anisotropic) refinement
- * cases compute the
- * <i>L<sup>2</sup></i>-projection
- * matrices from the children to
- * a coarse cell.
- *
- * Note, that
- * <code>matrices[refinement_case-1][child]</code>
- * includes the projection (or restriction)
- * matrix of child
- * <code>child</code> for the
- * RefinementCase
- * <code>refinement_case</code>. Here,
- * we use
- * <code>refinement_case-1</code>
- * instead of
- * <code>refinement_case</code>
- * 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
- * <tt>RefinementCase<dim>::isotropic_refinement</tt>
- * vectors of FullMatrix
- * objects. Each vector
- * corresponds to one
- * RefinementCase @p
- * refinement_case and is of the
- * vector size
- * <tt>GeometryInfo<dim>::n_children(refinement_case)</tt>. This
- * is the format used in
- * FiniteElement, where we want
- * to use this function mostly.
- *
- * @arg isotropic_only Set
- * to <code>true</code> if you only
- * want to compute matrices for
- * isotropic refinement.
- */
- template <int dim, typename number, int spacedim>
- static void compute_projection_matrices(
- const FiniteElement<dim,spacedim> &fe,
- std::vector<std::vector<FullMatrix<number> > >& 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 <tt>u<sub>q</sub>, 0 <= q <
- * Q:=quadrature.size()</tt>
- * defined at the quadrature points of a
- * cell, with the points defined by the
- * given <tt>rhs_quadrature</tt>
- * object. We may then want to ask for
- * that finite element function (on a
- * single cell) <tt>v<sub>h</sub></tt> in
- * the finite-dimensional space defined
- * by the given FE object that is the
- * projection of <tt>u</tt> in the
- * following sense:
- *
- * Usually, the projection
- * <tt>v<sub>h</sub></tt> is that
- * function that satisfies
- * <tt>(v<sub>h</sub>,w)=(u,w)</tt> for
- * all discrete test functions
- * <tt>w</tt>. In the present case, we
- * can't evaluate the right hand side,
- * since <tt>u</tt> is only defined in
- * the quadrature points given by
- * <tt>rhs_quadrature</tt>, so we replace
- * it by a quadrature
- * approximation. Likewise, the left hand
- * side is approximated using the
- * <tt>lhs_quadrature</tt> 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 <tt>q</tt> 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
- * <tt>V<sub>h</sub></tt> of
- * <tt>v<sub>h</sub></tt> that satisfies
- * the following system of linear
- * equations: <tt>M V<sub>h</sub> = Q
- * U</tt>, where
- * <tt>M<sub>ij</sub>=(phi_i,phi_j)</tt>
- * is the mass matrix approximated by
- * <tt>lhs_quadrature</tt>, and
- * <tt>Q</tt> is the matrix
- * <tt>Q<sub>iq</sub>=phi<sub>i</sub>(x<sub>q</sub>)
- * w<sub>q</sub></tt> where
- * <tt>w<sub>q</sub></tt> are quadrature
- * weights; <tt>U</tt> is the vector of
- * quadrature point data
- * <tt>u<sub>q</sub></tt>.
- *
- * In order to then get the nodal
- * representation <tt>V<sub>h</sub></tt>
- * of the projection of <tt>U</tt>, one
- * computes <tt>V<sub>h</sub> = X U,
- * X=M<sup>-1</sup> Q</tt>. The purpose
- * of this function is to compute the
- * matrix <tt>X</tt> 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 <tt>U</tt> 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
- * <tt>FEValues::get_function_values</tt>
- * functions do; to make things a little
- * simpler, the
- * <tt>FETools::compute_interpolation_to_quadrature_points_matrix</tt>
- * 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
- * <tt>ths_quadrature</tt> equal
- * <tt>fe.get_unit_support_points()</tt>,
- * 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 <tt>Q</tt> in
- * that case is the identity
- * matrix, the mass matrix
- * <tt>M</tt> is not equal to the
- * identity matrix, except for
- * the special case that the
- * quadrature formula
- * <tt>lhs_quadrature</tt> 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 <int dim, int spacedim>
- static
- void
- compute_projection_from_quadrature_points_matrix (const FiniteElement<dim,spacedim> &fe,
- const Quadrature<dim> &lhs_quadrature,
- const Quadrature<dim> &rhs_quadrature,
- FullMatrix<double> &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 <int dim, int spacedim>
- static
- void
- compute_interpolation_to_quadrature_points_matrix (const FiniteElement<dim,spacedim> &fe,
- const Quadrature<dim> &quadrature,
- FullMatrix<double> &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 <int dim>
- static
- void
- compute_projection_from_quadrature_points(
- const FullMatrix<double> &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 <int dim>
- static
- void
- compute_projection_from_quadrature_points(
- const FullMatrix<double> &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<dim-1> objects are required. See the
- * documentation on the Quadrature class for more information.
- */
- template <int dim, int spacedim>
- static
- void
- compute_projection_from_face_quadrature_points_matrix (const FiniteElement<dim, spacedim> &fe,
- const Quadrature<dim-1> &lhs_quadrature,
- const Quadrature<dim-1> &rhs_quadrature,
- const typename DoFHandler<dim, spacedim>::active_cell_iterator & cell,
- unsigned int face,
- FullMatrix<double> &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 <int dim, int spacedim,
- template <int,int> class DH1,
- template <int,int> class DH2,
- class InVector, class OutVector>
- static
- void
- interpolate (const DH1<dim,spacedim> &dof1,
- const InVector &u1,
- const DH2<dim,spacedim> &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 <int dim, int spacedim,
- template <int, int> class DH1,
- template <int, int> class DH2,
- class InVector, class OutVector>
- static void interpolate (const DH1<dim,spacedim> &dof1,
- const InVector &u1,
- const DH2<dim,spacedim> &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 <int dim, class InVector, class OutVector, int spacedim>
- static void back_interpolate (const DoFHandler<dim,spacedim> &dof1,
- const InVector &u1,
- const FiniteElement<dim,spacedim> &fe2,
- OutVector &u1_interpolated);
-
- /**
- * Same as last function, except
- * that the dof handler objects
- * might be of type
- * @p hp::DoFHandler.
- */
- template <int dim,
- template <int> class DH,
- class InVector, class OutVector, int spacedim>
- static void back_interpolate (const DH<dim> &dof1,
- const InVector &u1,
- const FiniteElement<dim,spacedim> &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 <int dim, class InVector, class OutVector, int spacedim>
- static void back_interpolate (const DoFHandler<dim,spacedim>& dof1,
- const ConstraintMatrix& constraints1,
- const InVector& u1,
- const DoFHandler<dim,spacedim>& 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 <int dim, class InVector, class OutVector, int spacedim>
- static void interpolation_difference(const DoFHandler<dim,spacedim> &dof1,
- const InVector &z1,
- const FiniteElement<dim,spacedim> &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 <int dim, class InVector, class OutVector, int spacedim>
- static void interpolation_difference(const DoFHandler<dim,spacedim>& dof1,
- const ConstraintMatrix& constraints1,
- const InVector& z1,
- const DoFHandler<dim,spacedim>& 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 <int dim, class InVector, class OutVector, int spacedim>
- static void project_dg (const DoFHandler<dim,spacedim>& dof1,
- const InVector& u1,
- const DoFHandler<dim,spacedim>& 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 <int dim, class InVector, class OutVector, int spacedim>
- static void extrapolate (const DoFHandler<dim,spacedim>& dof1,
- const InVector& z1,
- const DoFHandler<dim,spacedim>& 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 <int dim, class InVector, class OutVector, int spacedim>
- static void extrapolate (const DoFHandler<dim,spacedim>& dof1,
- const InVector& z1,
- const DoFHandler<dim,spacedim>& 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<dim> 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 <int dim>
- static void
- hierarchic_to_lexicographic_numbering (const FiniteElementData<dim> &fe_data,
- std::vector<unsigned int> &h2l);
-
- /**
- * Like the previous function but
- * instead of returning its
- * result through the last
- * argument return it as a value.
- */
- template <int dim>
- static
- std::vector<unsigned int>
- hierarchic_to_lexicographic_numbering (const FiniteElementData<dim> &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 <int dim>
- static void
- lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe_data,
- std::vector<unsigned int> &l2h);
-
- /**
- * Like the previous function but
- * instead of returning its
- * result through the last
- * argument return it as a value.
- */
- template <int dim>
- static
- std::vector<unsigned int>
- lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &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:
- *
- * <ul><li> Dimension template
- * parameters <2> etc. can
- * be omitted. Alternatively, the
- * explicit number can be
- * replaced by <tt>dim</tt> or
- * <tt>d</tt>. If a number is
- * given, it <b>must</b> match
- * the template parameter of this
- * function.
- *
- * <li> The powers used for
- * FESystem may either be numbers
- * or can be
- * replaced by <tt>dim</tt> or
- * <tt>d</tt>.
- * </ul>
- *
- * 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 <int dim>
- static
- FiniteElement<dim, dim> *
- 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
- * <tt>new</tt>. 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
- * <code>MyFE</code> be created
- * whenever the name
- * <code>my_fe</code> is given to
- * get_fe_from_name, you will
- * want the second argument to
- * this function be of type
- * FEFactory@<MyFE@>, 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 <int dim, int spacedim>
- static void add_fe_name (const std::string& name,
- const FEFactoryBase<dim,spacedim>* 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
- * <tt>dim</tt> or
- * <tt>d</tt>. 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<class FE>
-FiniteElement<FE::dimension, FE::dimension>*
-FETools::FEFactory<FE>::get (const unsigned int degree) const
-{
- return new FE(degree);
+ FullMatrix<number> &difference_matrix);
+
+ /**
+ * Compute the local
+ * $L^2$-projection matrix from
+ * fe1 to fe2.
+ */
+ template <int dim, typename number, int spacedim>
+ void get_projection_matrix(const FiniteElement<dim,spacedim> &fe1,
+ const FiniteElement<dim,spacedim> &fe2,
+ FullMatrix<number> &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.
+ * <ol>
+ *
+ * <li>Define the space of shape
+ * functions using an arbitrary
+ * basis <i>w<sub>j</sub></i> and
+ * compute the matrix <i>M</i> of
+ * node functionals
+ * <i>N<sub>i</sub></i> applied
+ * to these basis functions.
+ *
+ * <li>Compute the basis
+ * <i>v<sub>j</sub></i> of the
+ * finite element shape function
+ * space by applying
+ * <i>M<sup>-1</sup></i> to the
+ * basis <i>w<sub>j</sub></i>.
+ * </ol>
+ *
+ * @note The FiniteElement must
+ * provide generalized support
+ * points and and interpolation
+ * functions.
+ */
+ template <int dim, int spacedim>
+ void compute_node_matrix(FullMatrix<double>& M,
+ const FiniteElement<dim,spacedim>& 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
+ * <code>matrices[refinement_case-1][child]</code>
+ * includes the embedding (or prolongation)
+ * matrix of child
+ * <code>child</code> for the
+ * RefinementCase
+ * <code>refinement_case</code>. Here,
+ * we use
+ * <code>refinement_case-1</code>
+ * instead of
+ * <code>refinement_case</code>
+ * 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<dim>::isotropic_refinement
+ * vectors of FullMatrix
+ * objects. Each vector
+ * corresponds to one
+ * RefinementCase @p
+ * refinement_case and is of the
+ * vector size
+ * GeometryInfo<dim>::n_children(refinement_case). This
+ * is the format used in
+ * FiniteElement, where we want
+ * to use this function mostly.
+ *
+ * @param isotropic_only Set
+ * to <code>true</code> if you only
+ * want to compute matrices for
+ * isotropic refinement.
+ */
+ template <int dim, typename number, int spacedim>
+ void compute_embedding_matrices(const FiniteElement<dim,spacedim> &fe,
+ std::vector<std::vector<FullMatrix<number> > >& 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
+ * <i>GeometryInfo<dim>::subfaces_per_face
+ * = 2<sup>dim-1</sup></i>
+ * 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 <int dim, typename number, int spacedim>
+ void
+ compute_face_embedding_matrices(const FiniteElement<dim,spacedim>& fe,
+ FullMatrix<number> (&matrices)[GeometryInfo<dim>::max_children_per_face],
+ const unsigned int face_coarse,
+ const unsigned int face_fine);
+
+ /**
+ * For all possible (isotropic
+ * and anisotropic) refinement
+ * cases compute the
+ * <i>L<sup>2</sup></i>-projection
+ * matrices from the children to
+ * a coarse cell.
+ *
+ * Note, that
+ * <code>matrices[refinement_case-1][child]</code>
+ * includes the projection (or restriction)
+ * matrix of child
+ * <code>child</code> for the
+ * RefinementCase
+ * <code>refinement_case</code>. Here,
+ * we use
+ * <code>refinement_case-1</code>
+ * instead of
+ * <code>refinement_case</code>
+ * 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
+ * <tt>RefinementCase<dim>::isotropic_refinement</tt>
+ * vectors of FullMatrix
+ * objects. Each vector
+ * corresponds to one
+ * RefinementCase @p
+ * refinement_case and is of the
+ * vector size
+ * <tt>GeometryInfo<dim>::n_children(refinement_case)</tt>. This
+ * is the format used in
+ * FiniteElement, where we want
+ * to use this function mostly.
+ *
+ * @arg isotropic_only Set
+ * to <code>true</code> if you only
+ * want to compute matrices for
+ * isotropic refinement.
+ */
+ template <int dim, typename number, int spacedim>
+ void compute_projection_matrices(
+ const FiniteElement<dim,spacedim> &fe,
+ std::vector<std::vector<FullMatrix<number> > >& 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 <tt>u<sub>q</sub>, 0 <= q <
+ * Q:=quadrature.size()</tt>
+ * defined at the quadrature points of a
+ * cell, with the points defined by the
+ * given <tt>rhs_quadrature</tt>
+ * object. We may then want to ask for
+ * that finite element function (on a
+ * single cell) <tt>v<sub>h</sub></tt> in
+ * the finite-dimensional space defined
+ * by the given FE object that is the
+ * projection of <tt>u</tt> in the
+ * following sense:
+ *
+ * Usually, the projection
+ * <tt>v<sub>h</sub></tt> is that
+ * function that satisfies
+ * <tt>(v<sub>h</sub>,w)=(u,w)</tt> for
+ * all discrete test functions
+ * <tt>w</tt>. In the present case, we
+ * can't evaluate the right hand side,
+ * since <tt>u</tt> is only defined in
+ * the quadrature points given by
+ * <tt>rhs_quadrature</tt>, so we replace
+ * it by a quadrature
+ * approximation. Likewise, the left hand
+ * side is approximated using the
+ * <tt>lhs_quadrature</tt> 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 <tt>q</tt> 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
+ * <tt>V<sub>h</sub></tt> of
+ * <tt>v<sub>h</sub></tt> that satisfies
+ * the following system of linear
+ * equations: <tt>M V<sub>h</sub> = Q
+ * U</tt>, where
+ * <tt>M<sub>ij</sub>=(phi_i,phi_j)</tt>
+ * is the mass matrix approximated by
+ * <tt>lhs_quadrature</tt>, and
+ * <tt>Q</tt> is the matrix
+ * <tt>Q<sub>iq</sub>=phi<sub>i</sub>(x<sub>q</sub>)
+ * w<sub>q</sub></tt> where
+ * <tt>w<sub>q</sub></tt> are quadrature
+ * weights; <tt>U</tt> is the vector of
+ * quadrature point data
+ * <tt>u<sub>q</sub></tt>.
+ *
+ * In order to then get the nodal
+ * representation <tt>V<sub>h</sub></tt>
+ * of the projection of <tt>U</tt>, one
+ * computes <tt>V<sub>h</sub> = X U,
+ * X=M<sup>-1</sup> Q</tt>. The purpose
+ * of this function is to compute the
+ * matrix <tt>X</tt> 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 <tt>U</tt> 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
+ * <tt>FEValues::get_function_values</tt>
+ * functions do; to make things a little
+ * simpler, the
+ * <tt>FETools::compute_interpolation_to_quadrature_points_matrix</tt>
+ * 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
+ * <tt>ths_quadrature</tt> equal
+ * <tt>fe.get_unit_support_points()</tt>,
+ * 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 <tt>Q</tt> in
+ * that case is the identity
+ * matrix, the mass matrix
+ * <tt>M</tt> is not equal to the
+ * identity matrix, except for
+ * the special case that the
+ * quadrature formula
+ * <tt>lhs_quadrature</tt> 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 <int dim, int spacedim>
+ void
+ compute_projection_from_quadrature_points_matrix (const FiniteElement<dim,spacedim> &fe,
+ const Quadrature<dim> &lhs_quadrature,
+ const Quadrature<dim> &rhs_quadrature,
+ FullMatrix<double> &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 <int dim, int spacedim>
+ void
+ compute_interpolation_to_quadrature_points_matrix (const FiniteElement<dim,spacedim> &fe,
+ const Quadrature<dim> &quadrature,
+ FullMatrix<double> &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 <int dim>
+ void
+ compute_projection_from_quadrature_points(
+ const FullMatrix<double> &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 <int dim>
+ void
+ compute_projection_from_quadrature_points(
+ const FullMatrix<double> &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<dim-1> objects are required. See the
+ * documentation on the Quadrature class for more information.
+ */
+ template <int dim, int spacedim>
+ void
+ compute_projection_from_face_quadrature_points_matrix (const FiniteElement<dim, spacedim> &fe,
+ const Quadrature<dim-1> &lhs_quadrature,
+ const Quadrature<dim-1> &rhs_quadrature,
+ const typename DoFHandler<dim, spacedim>::active_cell_iterator & cell,
+ unsigned int face,
+ FullMatrix<double> &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 <int dim, int spacedim,
+ template <int,int> class DH1,
+ template <int,int> class DH2,
+ class InVector, class OutVector>
+ void
+ interpolate (const DH1<dim,spacedim> &dof1,
+ const InVector &u1,
+ const DH2<dim,spacedim> &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 <int dim, int spacedim,
+ template <int, int> class DH1,
+ template <int, int> class DH2,
+ class InVector, class OutVector>
+ void interpolate (const DH1<dim,spacedim> &dof1,
+ const InVector &u1,
+ const DH2<dim,spacedim> &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 <int dim, class InVector, class OutVector, int spacedim>
+ void back_interpolate (const DoFHandler<dim,spacedim> &dof1,
+ const InVector &u1,
+ const FiniteElement<dim,spacedim> &fe2,
+ OutVector &u1_interpolated);
+
+ /**
+ * Same as last function, except
+ * that the dof handler objects
+ * might be of type
+ * @p hp::DoFHandler.
+ */
+ template <int dim,
+ template <int> class DH,
+ class InVector, class OutVector, int spacedim>
+ void back_interpolate (const DH<dim> &dof1,
+ const InVector &u1,
+ const FiniteElement<dim,spacedim> &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 <int dim, class InVector, class OutVector, int spacedim>
+ void back_interpolate (const DoFHandler<dim,spacedim>& dof1,
+ const ConstraintMatrix& constraints1,
+ const InVector& u1,
+ const DoFHandler<dim,spacedim>& 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 <int dim, class InVector, class OutVector, int spacedim>
+ void interpolation_difference(const DoFHandler<dim,spacedim> &dof1,
+ const InVector &z1,
+ const FiniteElement<dim,spacedim> &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 <int dim, class InVector, class OutVector, int spacedim>
+ void interpolation_difference(const DoFHandler<dim,spacedim>& dof1,
+ const ConstraintMatrix& constraints1,
+ const InVector& z1,
+ const DoFHandler<dim,spacedim>& 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 <int dim, class InVector, class OutVector, int spacedim>
+ void project_dg (const DoFHandler<dim,spacedim>& dof1,
+ const InVector& u1,
+ const DoFHandler<dim,spacedim>& 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 <int dim, class InVector, class OutVector, int spacedim>
+ void extrapolate (const DoFHandler<dim,spacedim>& dof1,
+ const InVector& z1,
+ const DoFHandler<dim,spacedim>& 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 <int dim, class InVector, class OutVector, int spacedim>
+ void extrapolate (const DoFHandler<dim,spacedim>& dof1,
+ const InVector& z1,
+ const DoFHandler<dim,spacedim>& 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<dim> 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 <int dim>
+ void
+ hierarchic_to_lexicographic_numbering (const FiniteElementData<dim> &fe_data,
+ std::vector<unsigned int> &h2l);
+
+ /**
+ * Like the previous function but
+ * instead of returning its
+ * result through the last
+ * argument return it as a value.
+ */
+ template <int dim>
+ std::vector<unsigned int>
+ hierarchic_to_lexicographic_numbering (const FiniteElementData<dim> &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 <int dim>
+ void
+ lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe_data,
+ std::vector<unsigned int> &l2h);
+
+ /**
+ * Like the previous function but
+ * instead of returning its
+ * result through the last
+ * argument return it as a value.
+ */
+ template <int dim>
+ std::vector<unsigned int>
+ lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &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:
+ *
+ * <ul><li> Dimension template
+ * parameters <2> etc. can
+ * be omitted. Alternatively, the
+ * explicit number can be
+ * replaced by <tt>dim</tt> or
+ * <tt>d</tt>. If a number is
+ * given, it <b>must</b> match
+ * the template parameter of this
+ * function.
+ *
+ * <li> The powers used for
+ * FESystem may either be numbers
+ * or can be
+ * replaced by <tt>dim</tt> or
+ * <tt>d</tt>.
+ * </ul>
+ *
+ * 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 <int dim>
+ FiniteElement<dim, dim> *
+ 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
+ * <tt>new</tt>. 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
+ * <code>MyFE</code> be created
+ * whenever the name
+ * <code>my_fe</code> is given to
+ * get_fe_from_name, you will
+ * want the second argument to
+ * this function be of type
+ * FEFactory@<MyFE@>, 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 <int dim, int spacedim>
+ void add_fe_name (const std::string& name,
+ const FEFactoryBase<dim,spacedim>* 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
+ * <tt>dim</tt> or
+ * <tt>d</tt>. 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 <class FE>
+ FiniteElement<FE::dimension, FE::dimension>*
+ FEFactory<FE>::get (const unsigned int degree) const
+ {
+ return new FE(degree);
+ }
+}
+
+#endif
/*@}*/
}
-template <int dim, int spacedim>
-FETools::FEFactoryBase<dim,spacedim>::~FEFactoryBase()
-{}
+namespace FETools
+{
+ template <int dim, int spacedim>
+ FEFactoryBase<dim,spacedim>::~FEFactoryBase()
+ {}
-template<int dim, int spacedim>
-void FETools::compute_component_wise(
- const FiniteElement<dim,spacedim>& element,
- std::vector<unsigned int>& renumbering,
- std::vector<std::vector<unsigned int> >& comp_start)
-{
- Assert(renumbering.size() == element.dofs_per_cell,
- ExcDimensionMismatch(renumbering.size(),
- element.dofs_per_cell));
+ template<int dim, int spacedim>
+ void compute_component_wise(
+ const FiniteElement<dim,spacedim>& element,
+ std::vector<unsigned int>& renumbering,
+ std::vector<std::vector<unsigned int> >& 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<comp_start.size();++i)
- {
- comp_start[i].resize(element.element_multiplicity(i));
- const unsigned int increment
- = element.base_element(i).dofs_per_cell;
+ unsigned int k=0;
+ for (unsigned int i=0;i<comp_start.size();++i)
+ {
+ comp_start[i].resize(element.element_multiplicity(i));
+ const unsigned int increment
+ = element.base_element(i).dofs_per_cell;
+
+ for (unsigned int j=0;j<comp_start[i].size();++j)
+ {
+ comp_start[i][j] = k;
+ k += increment;
+ }
+ }
+
+ // For each index i of the
+ // unstructured cellwise
+ // numbering, renumbering
+ // contains the index of the
+ // cell-block numbering
+ for (unsigned int i=0;i<element.dofs_per_cell;++i)
+ {
+ std::pair<std::pair<unsigned int, unsigned int>, 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<comp_start[i].size();++j)
+
+
+ template<int dim, int spacedim>
+ void compute_block_renumbering (
+ const FiniteElement<dim,spacedim>& element,
+ std::vector<unsigned int>& renumbering,
+ std::vector<unsigned int>& 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<element.n_base_elements();++b)
+ for (unsigned int m=0;m<element.element_multiplicity(b);++m)
{
- comp_start[i][j] = k;
- k += increment;
+ block_data[i++] = (return_start_indices)
+ ? k
+ : (element.base_element(b).n_dofs_per_cell());
+ k += element.base_element(b).n_dofs_per_cell();
}
- }
+ Assert (i == element.n_blocks(), ExcInternalError());
- // For each index i of the
- // unstructured cellwise
- // numbering, renumbering
- // contains the index of the
- // cell-block numbering
- for (unsigned int i=0;i<element.dofs_per_cell;++i)
- {
- std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
- indices = element.system_to_base_index(i);
- renumbering[i] = comp_start[indices.first.first][indices.first.second]
- +indices.second;
- }
-}
+ std::vector<unsigned int> start_indices(block_data.size());
+ k = 0;
+ for (unsigned int i=0;i<block_data.size();++i)
+ if (return_start_indices)
+ start_indices[i] = block_data[i];
+ else
+ {
+ start_indices[i] = k;
+ k += block_data[i];
+ }
+//TODO:[GK] This does not work for a single RT
+ for (unsigned int i=0;i<element.dofs_per_cell;++i)
+ {
+ std::pair<unsigned int, unsigned int>
+ indices = element.system_to_block_index(i);
+ renumbering[i] = start_indices[indices.first]
+ +indices.second;
+ }
+ }
-template<int dim, int spacedim>
-void FETools::compute_block_renumbering (
- const FiniteElement<dim,spacedim>& element,
- std::vector<unsigned int>& renumbering,
- std::vector<unsigned int>& 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<element.n_base_elements();++b)
- for (unsigned int m=0;m<element.element_multiplicity(b);++m)
+
+ template <int dim, typename number, int spacedim>
+ void get_interpolation_matrix (const FiniteElement<dim,spacedim> &fe1,
+ const FiniteElement<dim,spacedim> &fe2,
+ FullMatrix<number> &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<unsigned int> start_indices(block_data.size());
- k = 0;
- for (unsigned int i=0;i<block_data.size();++i)
- if (return_start_indices)
- start_indices[i] = block_data[i];
- else
+ catch (typename FiniteElement<dim,spacedim>::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<element.dofs_per_cell;++i)
- {
- std::pair<unsigned int, unsigned int>
- 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<Point<dim> > &
+ fe2_support_points = fe2.get_unit_support_points ();
+
+ typedef FiniteElement<dim,spacedim> FEL;
+ Assert(fe2_support_points.size()==fe2.dofs_per_cell,
+ typename FEL::ExcFEHasNoSupportPoints());
+
+ for (unsigned int i=0; i<fe2.dofs_per_cell; ++i)
+ {
+ const unsigned int i1 = fe2.system_to_component_index(i).first;
+ for (unsigned int j=0; j<fe1.dofs_per_cell; ++j)
+ {
+ const unsigned int j1 = fe1.system_to_component_index(j).first;
+ if (i1==j1)
+ interpolation_matrix(i,j) = fe1.shape_value (j,fe2_support_points[i]);
+ else
+ interpolation_matrix(i,j)=0.;
+ }
+ }
+ }
-template <int dim, typename number, int spacedim>
-void FETools::get_interpolation_matrix (const FiniteElement<dim,spacedim> &fe1,
- const FiniteElement<dim,spacedim> &fe2,
- FullMatrix<number> &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<dim,spacedim>::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<Point<dim> > &
- fe2_support_points = fe2.get_unit_support_points ();
-
- typedef FiniteElement<dim,spacedim> FEL;
- Assert(fe2_support_points.size()==fe2.dofs_per_cell,
- typename FEL::ExcFEHasNoSupportPoints());
-
- for (unsigned int i=0; i<fe2.dofs_per_cell; ++i)
- {
- const unsigned int i1 = fe2.system_to_component_index(i).first;
- for (unsigned int j=0; j<fe1.dofs_per_cell; ++j)
- {
- const unsigned int j1 = fe1.system_to_component_index(j).first;
- if (i1==j1)
- interpolation_matrix(i,j) = fe1.shape_value (j,fe2_support_points[i]);
- else
- interpolation_matrix(i,j)=0.;
- }
- }
-}
+ template <int dim, typename number, int spacedim>
+ void get_back_interpolation_matrix(const FiniteElement<dim,spacedim> &fe1,
+ const FiniteElement<dim,spacedim> &fe2,
+ FullMatrix<number> &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<number> first_matrix (fe2.dofs_per_cell, fe1.dofs_per_cell);
+ FullMatrix<number> 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 <int dim, typename number, int spacedim>
-void FETools::get_back_interpolation_matrix(const FiniteElement<dim,spacedim> &fe1,
+ template <int dim, typename number, int spacedim>
+ void get_interpolation_difference_matrix (const FiniteElement<dim,spacedim> &fe1,
const FiniteElement<dim,spacedim> &fe2,
- FullMatrix<number> &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<number> first_matrix (fe2.dofs_per_cell, fe1.dofs_per_cell);
- FullMatrix<number> 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<number> &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<number> interpolation_matrix(fe1.dofs_per_cell);
+ get_back_interpolation_matrix(fe1, fe2, interpolation_matrix);
+
+ for (unsigned int i=0; i<fe1.dofs_per_cell; ++i)
+ difference_matrix(i,i) = 1.;
+
+ // compute difference
+ difference_matrix.add (-1, interpolation_matrix);
+ }
-template <int dim, typename number, int spacedim>
-void FETools::get_interpolation_difference_matrix (const FiniteElement<dim,spacedim> &fe1,
- const FiniteElement<dim,spacedim> &fe2,
- FullMatrix<number> &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<number> interpolation_matrix(fe1.dofs_per_cell);
- get_back_interpolation_matrix(fe1, fe2, interpolation_matrix);
-
- for (unsigned int i=0; i<fe1.dofs_per_cell; ++i)
- difference_matrix(i,i) = 1.;
-
- // compute difference
- difference_matrix.add (-1, interpolation_matrix);
-}
+ template <int dim, typename number, int spacedim>
+ void get_projection_matrix (const FiniteElement<dim,spacedim> &fe1,
+ const FiniteElement<dim,spacedim> &fe2,
+ FullMatrix<number> &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<dim,spacedim> 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<dim> quadrature(degree+1);
+ // Set up FEValues.
+ const UpdateFlags flags = update_values | update_quadrature_points | update_JxW_values;
+ FEValues<dim> val1 (fe1, quadrature, update_values);
+ val1.reinit (tr.begin_active());
+ FEValues<dim> val2 (fe2, quadrature, flags);
+ val2.reinit (tr.begin_active());
+
+ // Integrate and invert mass matrix
+ // This happens in the target space
+ FullMatrix<double> mass (n2, n2);
+
+ for (unsigned int k=0;k<quadrature.size();++k)
+ {
+ const double w = val2.JxW(k);
+ for (unsigned int i=0;i<n2;++i)
+ {
+ const double v = val2.shape_value(i,k);
+ for (unsigned int j=0;j<n2;++j)
+ mass(i,j) += w*v * val2.shape_value(j,k);
+ }
+ }
+ // Gauss-Jordan should be
+ // sufficient since we expect the
+ // mass matrix to be
+ // well-conditioned
+ mass.gauss_jordan();
+ // Now, test every function of fe1
+ // with test functions of fe2 and
+ // compute the projection of each
+ // unit vector.
+ Vector<double> b(n2);
+ Vector<double> x(n2);
+ for (unsigned int j=0;j<n1;++j)
+ {
+ b = 0.;
+ for (unsigned int i=0;i<n2;++i)
+ for (unsigned int k=0;k<quadrature.size();++k)
+ {
+ const double w = val2.JxW(k);
+ const double u = val1.shape_value(j,k);
+ const double v = val2.shape_value(i,k);
+ b(i) += u*v*w;
+ }
-template <int dim, typename number, int spacedim>
-void FETools::get_projection_matrix (const FiniteElement<dim,spacedim> &fe1,
- const FiniteElement<dim,spacedim> &fe2,
- FullMatrix<number> &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<dim,spacedim> 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<dim> quadrature(degree+1);
- // Set up FEValues.
- const UpdateFlags flags = update_values | update_quadrature_points | update_JxW_values;
- FEValues<dim> val1 (fe1, quadrature, update_values);
- val1.reinit (tr.begin_active());
- FEValues<dim> val2 (fe2, quadrature, flags);
- val2.reinit (tr.begin_active());
-
- // Integrate and invert mass matrix
- // This happens in the target space
- FullMatrix<double> mass (n2, n2);
-
- for (unsigned int k=0;k<quadrature.size();++k)
- {
- const double w = val2.JxW(k);
- for (unsigned int i=0;i<n2;++i)
- {
- const double v = val2.shape_value(i,k);
- for (unsigned int j=0;j<n2;++j)
- mass(i,j) += w*v * val2.shape_value(j,k);
- }
- }
- // Gauss-Jordan should be
- // sufficient since we expect the
- // mass matrix to be
- // well-conditioned
- mass.gauss_jordan();
-
- // Now, test every function of fe1
- // with test functions of fe2 and
- // compute the projection of each
- // unit vector.
- Vector<double> b(n2);
- Vector<double> x(n2);
-
- for (unsigned int j=0;j<n1;++j)
- {
- b = 0.;
- for (unsigned int i=0;i<n2;++i)
- for (unsigned int k=0;k<quadrature.size();++k)
- {
- const double w = val2.JxW(k);
- const double u = val1.shape_value(j,k);
- const double v = val2.shape_value(i,k);
- b(i) += u*v*w;
- }
-
- // Multiply by the inverse
- mass.vmult(x,b);
- for (unsigned int i=0;i<n2;++i)
- matrix(i,j) = x(i);
- }
-}
+ // Multiply by the inverse
+ mass.vmult(x,b);
+ for (unsigned int i=0;i<n2;++i)
+ matrix(i,j) = x(i);
+ }
+ }
-template<int dim, int spacedim>
-void
-FETools::compute_node_matrix(
- FullMatrix<double>& N,
- const FiniteElement<dim,spacedim>& 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<Point<dim> >& points = fe.get_generalized_support_points();
-
- // We need the values of the
- // polynomials in all generalized
- // support points.
- std::vector<std::vector<double> >
- values (dim, std::vector<double>(points.size()));
-
- // In this vector, we store the
- // result of the interpolation
- std::vector<double> 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<n_dofs;++i)
- {
- for (unsigned int k=0;k<values[0].size();++k)
- for (unsigned int d=0;d<dim;++d)
- values[d][k] = fe.shape_value_component(i,points[k],d);
- fe.interpolate(local_dofs, values);
- // Enter the interpolated dofs
- // into the matrix
- for (unsigned int j=0;j<n_dofs;++j)
- N(j,i) = local_dofs[j];
- }
-}
+ template<int dim, int spacedim>
+ void
+ compute_node_matrix(
+ FullMatrix<double>& N,
+ const FiniteElement<dim,spacedim>& 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<Point<dim> >& points = fe.get_generalized_support_points();
+
+ // We need the values of the
+ // polynomials in all generalized
+ // support points.
+ std::vector<std::vector<double> >
+ values (dim, std::vector<double>(points.size()));
+
+ // In this vector, we store the
+ // result of the interpolation
+ std::vector<double> 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<n_dofs;++i)
+ {
+ for (unsigned int k=0;k<values[0].size();++k)
+ for (unsigned int d=0;d<dim;++d)
+ values[d][k] = fe.shape_value_component(i,points[k],d);
+ fe.interpolate(local_dofs, values);
+ // Enter the interpolated dofs
+ // into the matrix
+ for (unsigned int j=0;j<n_dofs;++j)
+ N(j,i) = local_dofs[j];
+ }
+ }
#if deal_II_dimension == 1
-template<>
-void
-FETools::compute_embedding_matrices(const FiniteElement<1,2> &,
- std::vector<std::vector<FullMatrix<double> > > &,
- const bool)
-{
- Assert(false, ExcNotImplemented());
-}
+ template<>
+ void
+ compute_embedding_matrices(const FiniteElement<1,2> &,
+ std::vector<std::vector<FullMatrix<double> > > &,
+ const bool)
+ {
+ Assert(false, ExcNotImplemented());
+ }
#elif deal_II_dimension == 2
-template<>
-void
-FETools::compute_embedding_matrices(const FiniteElement<2,3>&,
- std::vector<std::vector<FullMatrix<double> > >&,
- const bool)
-{
- Assert(false, ExcNotImplemented());
-}
+ template<>
+ void
+ compute_embedding_matrices(const FiniteElement<2,3>&,
+ std::vector<std::vector<FullMatrix<double> > >&,
+ const bool)
+ {
+ Assert(false, ExcNotImplemented());
+ }
#endif
-namespace {
- template<int dim, typename number, int spacedim>
- void
- compute_embedding_matrices_for_refinement_case (const FiniteElement<dim, spacedim>& fe,
- std::vector<FullMatrix<number> >& matrices,
- const unsigned int ref_case)
- {
- const unsigned int n = fe.dofs_per_cell;
- const unsigned int nc = GeometryInfo<dim>::n_children(RefinementCase<dim>(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<int dim, typename number, int spacedim>
+ void
+ compute_embedding_matrices_for_refinement_case (const FiniteElement<dim, spacedim>& fe,
+ std::vector<FullMatrix<number> >& matrices,
+ const unsigned int ref_case)
+ {
+ const unsigned int n = fe.dofs_per_cell;
+ const unsigned int nc = GeometryInfo<dim>::n_children(RefinementCase<dim>(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<dim,spacedim> tria;
- GridGenerator::hyper_cube (tria, 0, 1);
- tria.begin_active()->set_refine_flag (RefinementCase<dim>(ref_case));
- tria.execute_coarsening_and_refinement ();
+ // Set up meshes, one with a single
+ // reference cell and refine it once
+ Triangulation<dim,spacedim> tria;
+ GridGenerator::hyper_cube (tria, 0, 1);
+ tria.begin_active()->set_refine_flag (RefinementCase<dim>(ref_case));
+ tria.execute_coarsening_and_refinement ();
- MappingCartesian<dim> mapping;
- const unsigned int degree = fe.degree;
- QGauss<dim> q_fine (degree+1);
- const unsigned int nq = q_fine.size();
+ MappingCartesian<dim> mapping;
+ const unsigned int degree = fe.degree;
+ QGauss<dim> q_fine (degree+1);
+ const unsigned int nq = q_fine.size();
- FEValues<dim> fine (mapping, fe, q_fine,
- update_quadrature_points |
- update_JxW_values |
- update_values);
+ FEValues<dim> 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
// This matrix is the same for all
// children.
- fine.reinit (tria.begin_active ());
- const unsigned int nd = fe.n_components ();
- FullMatrix<number> 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<double> H (A);
- static Threads::Mutex mutex;
- Vector<number> v_coarse (nq * nd);
- Vector<number> v_fine (n);
- unsigned int cell_number = 0;
-
- for (typename Triangulation<dim>::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<number> 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<double> H (A);
+ static Threads::Mutex mutex;
+ Vector<number> v_coarse (nq * nd);
+ Vector<number> v_fine (n);
+ unsigned int cell_number = 0;
+
+ for (typename Triangulation<dim>::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<dim> q_coarse (fine.get_quadrature_points (),
- fine.get_JxW_values ());
- FEValues<dim> coarse (mapping, fe, q_coarse, update_values);
-
- coarse.reinit (tria.begin (0));
+ const Quadrature<dim> q_coarse (fine.get_quadrature_points (),
+ fine.get_JxW_values ());
+ FEValues<dim> coarse (mapping, fe, q_coarse, update_values);
- FullMatrix<double> &this_matrix = matrices[cell_number];
-
- v_coarse = 0;
+ coarse.reinit (tria.begin (0));
+
+ FullMatrix<double> &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
// 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<dim>::n_children (RefinementCase<dim> (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 <int dim, typename number, int spacedim>
-void
-FETools::compute_embedding_matrices(const FiniteElement<dim,spacedim>& fe,
- std::vector<std::vector<FullMatrix<number> > >& matrices,
- const bool isotropic_only)
-{
- Threads::TaskGroup<void> task_group;
-
- // loop over all possible refinement cases
- unsigned int ref_case = (isotropic_only)
- ? RefinementCase<dim>::isotropic_refinement
- : RefinementCase<dim>::cut_x;
-
- for (;ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
- task_group += Threads::new_task (&compute_embedding_matrices_for_refinement_case<dim, number, spacedim>,
- fe, matrices[ref_case-1], ref_case);
-
- task_group.join_all ();
-}
+ Assert (cell_number == GeometryInfo<dim>::n_children (RefinementCase<dim> (ref_case)),
+ ExcInternalError ());
+ }
+ }
// This function is tested by tests/fe/internals, since it produces the matrices printed there
+ template <int dim, typename number, int spacedim>
+ void
+ compute_embedding_matrices(const FiniteElement<dim,spacedim>& fe,
+ std::vector<std::vector<FullMatrix<number> > >& matrices,
+ const bool isotropic_only)
+ {
+ Threads::TaskGroup<void> task_group;
-//TODO:[GK] Is this correct for vector valued?
-template <int dim, typename number, int spacedim>
-void
-FETools::compute_face_embedding_matrices(const FiniteElement<dim,spacedim>& fe,
- FullMatrix<number> (&matrices)[GeometryInfo<dim>::max_children_per_face],
- const unsigned int face_coarse,
- const unsigned int face_fine)
-{
- const unsigned int nc = GeometryInfo<dim>::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<dim>::isotropic_refinement
+ : RefinementCase<dim>::cut_x;
- 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));
- }
+ for (;ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
+ task_group += Threads::new_task (&compute_embedding_matrices_for_refinement_case<dim, number, spacedim>,
+ fe, matrices[ref_case-1], ref_case);
- // Set up meshes, one with a single
- // reference cell and refine it once
- Triangulation<dim,spacedim> tria;
- GridGenerator::hyper_cube (tria, 0, 1);
- tria.refine_global(1);
-
- MappingCartesian<dim> mapping;
- QGauss<dim-1> q_gauss(degree+1);
- const Quadrature<dim> q_fine = QProjector<dim>::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<unsigned int> face_c_dofs(n);
- std::vector<unsigned int> face_f_dofs(n);
- unsigned int k=0;
- for (unsigned int i=0;i<GeometryInfo<dim>::vertices_per_face;++i)
- {
- const unsigned int offset_c = GeometryInfo<dim>::face_to_cell_vertices(face_coarse, i)
- *fe.dofs_per_vertex;
- const unsigned int offset_f = GeometryInfo<dim>::face_to_cell_vertices(face_fine, i)
- *fe.dofs_per_vertex;
- for (unsigned int j=0;j<fe.dofs_per_vertex;++j)
- {
- face_c_dofs[k] = offset_c + j;
- face_f_dofs[k] = offset_f + j;
- ++k;
- }
- }
- for (unsigned int i=1;i<=GeometryInfo<dim>::lines_per_face;++i)
- {
- const unsigned int offset_c = fe.first_line_index
- + GeometryInfo<dim>::face_to_cell_lines(face_coarse, i-1)
- *fe.dofs_per_line;
- const unsigned int offset_f = fe.first_line_index
- + GeometryInfo<dim>::face_to_cell_lines(face_fine, i-1)
- *fe.dofs_per_line;
- for (unsigned int j=0;j<fe.dofs_per_line;++j)
- {
- face_c_dofs[k] = offset_c + j;
- face_f_dofs[k] = offset_f + j;
- ++k;
- }
- }
- for (unsigned int i=1;i<=GeometryInfo<dim>::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<fe.dofs_per_quad;++j)
- {
- face_c_dofs[k] = offset_c + j;
- face_f_dofs[k] = offset_f + j;
- ++k;
- }
- }
- Assert (k == fe.dofs_per_face, ExcInternalError());
+ task_group.join_all ();
+ }
- FEValues<dim> 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 <int dim, typename number, int spacedim>
+ void
+ compute_face_embedding_matrices(const FiniteElement<dim,spacedim>& fe,
+ FullMatrix<number> (&matrices)[GeometryInfo<dim>::max_children_per_face],
+ const unsigned int face_coarse,
+ const unsigned int face_fine)
+ {
+ const unsigned int nc = GeometryInfo<dim>::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<number> 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(face_f_dofs[j],k,d);
+ 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));
+ }
- Householder<double> H(A);
+ // Set up meshes, one with a single
+ // reference cell and refine it once
+ Triangulation<dim,spacedim> tria;
+ GridGenerator::hyper_cube (tria, 0, 1);
+ tria.refine_global(1);
- Vector<number> v_coarse(nq*nd);
- Vector<number> v_fine(n);
+ MappingCartesian<dim> mapping;
+ QGauss<dim-1> q_gauss(degree+1);
+ const Quadrature<dim> q_fine = QProjector<dim>::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<unsigned int> face_c_dofs(n);
+ std::vector<unsigned int> face_f_dofs(n);
+ unsigned int k=0;
+ for (unsigned int i=0;i<GeometryInfo<dim>::vertices_per_face;++i)
+ {
+ const unsigned int offset_c = GeometryInfo<dim>::face_to_cell_vertices(face_coarse, i)
+ *fe.dofs_per_vertex;
+ const unsigned int offset_f = GeometryInfo<dim>::face_to_cell_vertices(face_fine, i)
+ *fe.dofs_per_vertex;
+ for (unsigned int j=0;j<fe.dofs_per_vertex;++j)
+ {
+ face_c_dofs[k] = offset_c + j;
+ face_f_dofs[k] = offset_f + j;
+ ++k;
+ }
+ }
+ for (unsigned int i=1;i<=GeometryInfo<dim>::lines_per_face;++i)
+ {
+ const unsigned int offset_c = fe.first_line_index
+ + GeometryInfo<dim>::face_to_cell_lines(face_coarse, i-1)
+ *fe.dofs_per_line;
+ const unsigned int offset_f = fe.first_line_index
+ + GeometryInfo<dim>::face_to_cell_lines(face_fine, i-1)
+ *fe.dofs_per_line;
+ for (unsigned int j=0;j<fe.dofs_per_line;++j)
+ {
+ face_c_dofs[k] = offset_c + j;
+ face_f_dofs[k] = offset_f + j;
+ ++k;
+ }
+ }
+ for (unsigned int i=1;i<=GeometryInfo<dim>::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<fe.dofs_per_quad;++j)
+ {
+ face_c_dofs[k] = offset_c + j;
+ face_f_dofs[k] = offset_f + j;
+ ++k;
+ }
+ }
+ Assert (k == fe.dofs_per_face, ExcInternalError());
- for (unsigned int cell_number = 0; cell_number < GeometryInfo<dim>::max_children_per_face;
- ++cell_number)
- {
- const Quadrature<dim> q_coarse
- = QProjector<dim>::project_to_subface(q_gauss, face_coarse, cell_number);
- FEValues<dim> coarse (mapping, fe, q_coarse, update_values);
+ FEValues<dim> fine (mapping, fe, q_fine,
+ update_quadrature_points | update_JxW_values | update_values);
- typename Triangulation<dim,spacedim>::active_cell_iterator fine_cell
- = tria.begin(0)->child(GeometryInfo<dim>::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<double> &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<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.
- 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 (face_c_dofs[i],k,d);
-
- // solve the least squares
- // problem.
- 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
- // maps a coarse grid
- // function to a fine grid
- // function, the columns
- // are fine grid.
- for (unsigned int j=0;j<n;++j)
- this_matrix(j,i) = v_fine(j);
- }
- // 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.;
- }
-}
+ // This matrix is the same for all
+ // children.
+ fine.reinit(tria.begin_active());
+ FullMatrix<number> 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(face_f_dofs[j],k,d);
-#if deal_II_dimension == 1
-template <>
-void
-FETools::compute_projection_matrices(const FiniteElement<1,2>&,
- std::vector<std::vector<FullMatrix<double> > >&, bool)
-{
- Assert(false, ExcNotImplemented());
-}
+ Householder<double> H(A);
-#elif deal_II_dimension == 2
-template <>
-void
-FETools::compute_projection_matrices(const FiniteElement<2,3>&,
- std::vector<std::vector<FullMatrix<double> > >&, bool)
-{
- Assert(false, ExcNotImplemented());
-}
+ Vector<number> v_coarse(nq*nd);
+ Vector<number> v_fine(n);
-#endif
-template <int dim, typename number, int spacedim>
-void
-FETools::compute_projection_matrices(const FiniteElement<dim,spacedim>& fe,
- std::vector<std::vector<FullMatrix<number> > >& 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<dim> mapping;
- QGauss<dim> q_fine(degree+1);
- const unsigned int nq = q_fine.size();
-
- // create mass matrix on coarse cell.
- FullMatrix<number> mass(n, n);
- {
- // set up a triangulation for coarse cell
- Triangulation<dim,spacedim> tr;
- GridGenerator::hyper_cube (tr, 0, 1);
+ for (unsigned int cell_number = 0; cell_number < GeometryInfo<dim>::max_children_per_face;
+ ++cell_number)
+ {
+ const Quadrature<dim> q_coarse
+ = QProjector<dim>::project_to_subface(q_gauss, face_coarse, cell_number);
+ FEValues<dim> coarse (mapping, fe, q_coarse, update_values);
- FEValues<dim> coarse (mapping, fe, q_fine,
- update_JxW_values | update_values);
+ typename Triangulation<dim,spacedim>::active_cell_iterator fine_cell
+ = tria.begin(0)->child(GeometryInfo<dim>::child_cell_on_face(
+ tria.begin(0)->refinement_case(), face_coarse, cell_number));
+ fine.reinit(fine_cell);
+ coarse.reinit(tria.begin(0));
- typename Triangulation<dim,spacedim>::cell_iterator coarse_cell
- = tr.begin(0);
- coarse.reinit (coarse_cell);
+ FullMatrix<double> &this_matrix = matrices[cell_number];
- const std::vector<double> & JxW = coarse.get_JxW_values();
- for (unsigned int i=0;i<n;++i)
- for (unsigned int j=0;j<n;++j)
- if (fe.is_primitive())
+ // Compute this once for each
+ // coarse grid basis function
+ for (unsigned int i=0;i<n;++i)
{
- const double * coarse_i = &coarse.shape_value(i,0);
- const double * coarse_j = &coarse.shape_value(j,0);
- double mass_ij = 0;
- for (unsigned int k=0;k<nq;++k)
- mass_ij += JxW[k] * coarse_i[k] * coarse_j[k];
- mass(i,j) = mass_ij;
- }
- else
- {
- double mass_ij = 0;
+ // The right hand side of
+ // the least squares
+ // problem consists of the
+ // function values of the
+ // coarse grid function in
+ // each quadrature point.
for (unsigned int d=0;d<nd;++d)
for (unsigned int k=0;k<nq;++k)
- mass_ij += JxW[k] * coarse.shape_value_component(i,k,d)
- * coarse.shape_value_component(j,k,d);
- mass(i,j) = mass_ij;
+ v_coarse(k*nd+d) = coarse.shape_value_component (face_c_dofs[i],k,d);
+
+ // solve the least squares
+ // problem.
+ 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
+ // maps a coarse grid
+ // function to a fine grid
+ // function, the columns
+ // are fine grid.
+ for (unsigned int j=0;j<n;++j)
+ this_matrix(j,i) = v_fine(j);
}
+ // 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.;
+ }
+ }
- // invert mass matrix
- mass.gauss_jordan();
+#if deal_II_dimension == 1
+ template <>
+ void
+ compute_projection_matrices(const FiniteElement<1,2>&,
+ std::vector<std::vector<FullMatrix<double> > >&, bool)
+ {
+ Assert(false, ExcNotImplemented());
}
- // loop over all possible
- // refinement cases
- unsigned int ref_case = (isotropic_only)
- ? RefinementCase<dim>::isotropic_refinement
- : RefinementCase<dim>::cut_x;
- for (;ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
- {
- const unsigned int
- nc = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
+#elif deal_II_dimension == 2
+ template <>
+ void
+ compute_projection_matrices(const FiniteElement<2,3>&,
+ std::vector<std::vector<FullMatrix<double> > >&, bool)
+ {
+ Assert(false, ExcNotImplemented());
+ }
- for (unsigned int i=0;i<nc;++i)
- {
- Assert(matrices[ref_case-1][i].n() == n,
- ExcDimensionMismatch(matrices[ref_case-1][i].n(),n));
- Assert(matrices[ref_case-1][i].m() == n,
- ExcDimensionMismatch(matrices[ref_case-1][i].m(),n));
- }
+#endif
- // create a respective refinement on the
- // triangulation
+
+ template <int dim, typename number, int spacedim>
+ void
+ compute_projection_matrices(const FiniteElement<dim,spacedim>& fe,
+ std::vector<std::vector<FullMatrix<number> > >& 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<dim> mapping;
+ QGauss<dim> q_fine(degree+1);
+ const unsigned int nq = q_fine.size();
+
+ // create mass matrix on coarse cell.
+ FullMatrix<number> mass(n, n);
+ {
+ // set up a triangulation for coarse cell
Triangulation<dim,spacedim> tr;
GridGenerator::hyper_cube (tr, 0, 1);
- tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
- tr.execute_coarsening_and_refinement();
- FEValues<dim> fine (mapping, fe, q_fine,
- update_quadrature_points | update_JxW_values |
- update_values);
+ FEValues<dim> coarse (mapping, fe, q_fine,
+ update_JxW_values | update_values);
typename Triangulation<dim,spacedim>::cell_iterator coarse_cell
= tr.begin(0);
+ coarse.reinit (coarse_cell);
- Vector<number> v_coarse(n);
- Vector<number> v_fine(n);
+ const std::vector<double> & JxW = coarse.get_JxW_values();
+ for (unsigned int i=0;i<n;++i)
+ for (unsigned int j=0;j<n;++j)
+ if (fe.is_primitive())
+ {
+ const double * coarse_i = &coarse.shape_value(i,0);
+ const double * coarse_j = &coarse.shape_value(j,0);
+ double mass_ij = 0;
+ for (unsigned int k=0;k<nq;++k)
+ mass_ij += JxW[k] * coarse_i[k] * coarse_j[k];
+ mass(i,j) = mass_ij;
+ }
+ else
+ {
+ double mass_ij = 0;
+ for (unsigned int d=0;d<nd;++d)
+ for (unsigned int k=0;k<nq;++k)
+ mass_ij += JxW[k] * coarse.shape_value_component(i,k,d)
+ * coarse.shape_value_component(j,k,d);
+ mass(i,j) = mass_ij;
+ }
- for (unsigned int cell_number=0;cell_number<nc;++cell_number)
- {
- FullMatrix<double> &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<dim> q_coarse (fine.get_quadrature_points(),
- fine.get_JxW_values());
- FEValues<dim> 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<dim>::isotropic_refinement
+ : RefinementCase<dim>::cut_x;
+ for (;ref_case <= RefinementCase<dim>::isotropic_refinement; ++ref_case)
+ {
+ const unsigned int
+ nc = GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
- const std::vector<double> & JxW = fine.get_JxW_values();
+ for (unsigned int i=0;i<nc;++i)
+ {
+ Assert(matrices[ref_case-1][i].n() == n,
+ ExcDimensionMismatch(matrices[ref_case-1][i].n(),n));
+ Assert(matrices[ref_case-1][i].m() == n,
+ ExcDimensionMismatch(matrices[ref_case-1][i].m(),n));
+ }
- // Outer loop over all fine
- // grid shape functions phi_j
- for (unsigned int j=0;j<fe.dofs_per_cell;++j)
- {
- for (unsigned int i=0; i<fe.dofs_per_cell;++i)
- {
- if (fe.is_primitive())
- {
- const double * coarse_i = &coarse.shape_value(i,0);
- const double * fine_j = &fine.shape_value(j,0);
+ // create a respective refinement on the
+ // triangulation
+ Triangulation<dim,spacedim> tr;
+ GridGenerator::hyper_cube (tr, 0, 1);
+ tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
+ tr.execute_coarsening_and_refinement();
- double update = 0;
- for (unsigned int k=0; k<nq; ++k)
- update += JxW[k] * coarse_i[k] * fine_j[k];
- v_fine(i) = update;
- }
- else
- {
- double update = 0;
- for (unsigned int d=0; d<nd; ++d)
+ FEValues<dim> fine (mapping, fe, q_fine,
+ update_quadrature_points | update_JxW_values |
+ update_values);
+
+ typename Triangulation<dim,spacedim>::cell_iterator coarse_cell
+ = tr.begin(0);
+
+ Vector<number> v_coarse(n);
+ Vector<number> v_fine(n);
+
+ for (unsigned int cell_number=0;cell_number<nc;++cell_number)
+ {
+ FullMatrix<double> &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<dim> q_coarse (fine.get_quadrature_points(),
+ fine.get_JxW_values());
+ FEValues<dim> coarse (mapping, fe, q_coarse, update_values);
+ coarse.reinit(coarse_cell);
+
+ // Build RHS
+
+ const std::vector<double> & JxW = fine.get_JxW_values();
+
+ // Outer loop over all fine
+ // grid shape functions phi_j
+ for (unsigned int j=0;j<fe.dofs_per_cell;++j)
+ {
+ for (unsigned int i=0; i<fe.dofs_per_cell;++i)
+ {
+ if (fe.is_primitive())
+ {
+ const double * coarse_i = &coarse.shape_value(i,0);
+ const double * fine_j = &fine.shape_value(j,0);
+
+ double update = 0;
for (unsigned int k=0; k<nq; ++k)
- update += JxW[k] * coarse.shape_value_component(i,k,d)
- * fine.shape_value_component(j,k,d);
- v_fine(i) = update;
- }
- }
+ update += JxW[k] * coarse_i[k] * fine_j[k];
+ v_fine(i) = update;
+ }
+ else
+ {
+ double update = 0;
+ for (unsigned int d=0; d<nd; ++d)
+ for (unsigned int k=0; k<nq; ++k)
+ update += JxW[k] * coarse.shape_value_component(i,k,d)
+ * fine.shape_value_component(j,k,d);
+ v_fine(i) = update;
+ }
+ }
- // RHS ready. Solve system
- // and enter row into
- // matrix
- mass.vmult (v_coarse, v_fine);
- for (unsigned int i=0;i<fe.dofs_per_cell;++i)
- this_matrix(i,j) = v_coarse(i);
- }
+ // RHS ready. Solve system
+ // and enter row into
+ // matrix
+ mass.vmult (v_coarse, v_fine);
+ for (unsigned int i=0;i<fe.dofs_per_cell;++i)
+ this_matrix(i,j) = v_coarse(i);
+ }
- // 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.;
- }
- }
-}
+ // 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.;
+ }
+ }
+ }
-template <int dim, int spacedim,
- template <int, int> class DH1,
- template <int, int> class DH2,
- class InVector, class OutVector>
-void
-FETools::interpolate(const DH1<dim, spacedim> &dof1,
- const InVector &u1,
- const DH2<dim, spacedim> &dof2,
- OutVector &u2)
-{
- ConstraintMatrix dummy;
- dummy.close();
- interpolate(dof1, u1, dof2, dummy, u2);
-}
+ template <int dim, int spacedim,
+ template <int, int> class DH1,
+ template <int, int> class DH2,
+ class InVector, class OutVector>
+ void
+ interpolate(const DH1<dim, spacedim> &dof1,
+ const InVector &u1,
+ const DH2<dim, spacedim> &dof2,
+ OutVector &u2)
+ {
+ ConstraintMatrix dummy;
+ dummy.close();
+ interpolate(dof1, u1, dof2, dummy, u2);
+ }
template <int, int> class DH2,
class InVector, class OutVector>
void
-FETools::interpolate (const DH1<dim, spacedim> &dof1,
- const InVector &u1,
- const DH2<dim, spacedim> &dof2,
- const ConstraintMatrix &constraints,
- OutVector &u2)
+interpolate (const DH1<dim, spacedim> &dof1,
+ const InVector &u1,
+ const DH2<dim, spacedim> &dof2,
+ const ConstraintMatrix &constraints,
+ OutVector &u2)
{
Assert(&dof1.get_tria()==&dof2.get_tria(), ExcTriangulationMismatch());
// matrices. shared_ptr make sure
// that memory is released again
std::map<const FiniteElement<dim,spacedim> *,
- std::map<const FiniteElement<dim,spacedim> *,
- std_cxx1x::shared_ptr<FullMatrix<double> > > >
+ std::map<const FiniteElement<dim,spacedim> *,
+ std_cxx1x::shared_ptr<FullMatrix<double> > > >
interpolation_matrices;
typename DH1<dim,spacedim>::active_cell_iterator cell1 = dof1.begin_active(),
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);
-template <int dim, class InVector, class OutVector, int spacedim>
-void
-FETools::back_interpolate(const DoFHandler<dim,spacedim> &dof1,
- const InVector &u1,
- const FiniteElement<dim,spacedim> &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<typename OutVector::value_type> u1_local(dofs_per_cell1);
- Vector<typename OutVector::value_type> u1_int_local(dofs_per_cell1);
-
- typename DoFHandler<dim,spacedim>::active_cell_iterator cell = dof1.begin_active(),
- endc = dof1.end();
-
- FullMatrix<double> 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<GeometryInfo<dim>::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 <int dim,
- template <int> class DH,
- class InVector, class OutVector, int spacedim>
-void
-FETools::back_interpolate(const DH<dim> &dof1,
- const InVector &u1,
- const FiniteElement<dim,spacedim> &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<typename OutVector::value_type> u1_local(DoFTools::max_dofs_per_cell(dof1));
- Vector<typename OutVector::value_type> u1_int_local(DoFTools::max_dofs_per_cell(dof1));
-
- typename DH<dim>::active_cell_iterator cell = dof1.begin_active(),
- endc = dof1.end();
-
- // map from possible fe objects in
- // dof1 to the back_interpolation
- // matrices
- std::map<const FiniteElement<dim> *,
- std_cxx1x::shared_ptr<FullMatrix<double> > > 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<GeometryInfo<dim>::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<FullMatrix<double> >
- (new FullMatrix<double>(dofs_per_cell1, dofs_per_cell1));
- get_back_interpolation_matrix(dof1.get_fe(), fe2,
- *interpolation_matrices[&cell->get_fe()]);
- }
+ template <int dim, class InVector, class OutVector, int spacedim>
+ void
+ back_interpolate(const DoFHandler<dim,spacedim> &dof1,
+ const InVector &u1,
+ const FiniteElement<dim,spacedim> &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<typename OutVector::value_type> u1_local(dofs_per_cell1);
+ Vector<typename OutVector::value_type> u1_int_local(dofs_per_cell1);
+
+ typename DoFHandler<dim,spacedim>::active_cell_iterator cell = dof1.begin_active(),
+ endc = dof1.end();
+
+ FullMatrix<double> 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<GeometryInfo<dim>::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 <int dim,
+ template <int> class DH,
+ class InVector, class OutVector, int spacedim>
+ void
+ back_interpolate(const DH<dim> &dof1,
+ const InVector &u1,
+ const FiniteElement<dim,spacedim> &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<typename OutVector::value_type> u1_local(DoFTools::max_dofs_per_cell(dof1));
+ Vector<typename OutVector::value_type> u1_int_local(DoFTools::max_dofs_per_cell(dof1));
-template <int dim, class InVector, class OutVector, int spacedim>
-void FETools::back_interpolate(const DoFHandler<dim,spacedim> &dof1,
- const ConstraintMatrix &constraints1,
- const InVector &u1,
- const DoFHandler<dim,spacedim> &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<typename OutVector::value_type> u2(dof2.n_dofs());
- interpolate(dof1, u1, dof2, constraints2, u2);
- interpolate(dof2, u2, dof1, constraints1, u1_interpolated);
- }
-}
+ typename DH<dim>::active_cell_iterator cell = dof1.begin_active(),
+ endc = dof1.end();
+ // map from possible fe objects in
+ // dof1 to the back_interpolation
+ // matrices
+ std::map<const FiniteElement<dim> *,
+ std_cxx1x::shared_ptr<FullMatrix<double> > > 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<GeometryInfo<dim>::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<FullMatrix<double> >
+ (new FullMatrix<double>(dofs_per_cell1, dofs_per_cell1));
+ get_back_interpolation_matrix(dof1.get_fe(), fe2,
+ *interpolation_matrices[&cell->get_fe()]);
+ }
-template <int dim, class InVector, class OutVector, int spacedim>
-void FETools::interpolation_difference (const DoFHandler<dim,spacedim> &dof1,
- const InVector &u1,
- const FiniteElement<dim,spacedim> &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<typename OutVector::value_type> u1_local(dofs_per_cell);
- Vector<typename OutVector::value_type> u1_diff_local(dofs_per_cell);
-
- FullMatrix<double> difference_matrix(dofs_per_cell, dofs_per_cell);
- FETools::get_interpolation_difference_matrix(dof1.get_fe(), fe2,
- difference_matrix);
-
- typename DoFHandler<dim,spacedim>::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<GeometryInfo<dim>::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 <int dim, class InVector, class OutVector, int spacedim>
-void FETools::interpolation_difference(const DoFHandler<dim,spacedim> &dof1,
- const ConstraintMatrix &constraints1,
- const InVector &u1,
- const DoFHandler<dim,spacedim> &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 <int dim, class InVector, class OutVector, int spacedim>
+ void back_interpolate(const DoFHandler<dim,spacedim> &dof1,
+ const ConstraintMatrix &constraints1,
+ const InVector &u1,
+ const DoFHandler<dim,spacedim> &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<typename OutVector::value_type> u2(dof2.n_dofs());
+ interpolate(dof1, u1, dof2, constraints2, u2);
+ interpolate(dof2, u2, dof1, constraints1, u1_interpolated);
+ }
+ }
-template <int dim, class InVector, class OutVector, int spacedim>
-void FETools::project_dg(const DoFHandler<dim,spacedim> &dof1,
- const InVector &u1,
- const DoFHandler<dim,spacedim> &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<dim,spacedim>::active_cell_iterator cell1 = dof1.begin_active();
- typename DoFHandler<dim,spacedim>::active_cell_iterator cell2 = dof2.begin_active();
- typename DoFHandler<dim,spacedim>::active_cell_iterator end = dof2.end();
+ template <int dim, class InVector, class OutVector, int spacedim>
+ void interpolation_difference (const DoFHandler<dim,spacedim> &dof1,
+ const InVector &u1,
+ const FiniteElement<dim,spacedim> &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<typename OutVector::value_type> u1_local(dofs_per_cell);
+ Vector<typename OutVector::value_type> u1_diff_local(dofs_per_cell);
+
+ FullMatrix<double> difference_matrix(dofs_per_cell, dofs_per_cell);
+ get_interpolation_difference_matrix(dof1.get_fe(), fe2,
+ difference_matrix);
+
+ typename DoFHandler<dim,spacedim>::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<GeometryInfo<dim>::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<double> u1_local(n1);
- Vector<double> u2_local(n2);
- std::vector<unsigned int> dofs(n2);
- FullMatrix<double> matrix(n2,n1);
- get_projection_matrix(dof1.get_fe(), dof2.get_fe(), matrix);
+ template <int dim, class InVector, class OutVector, int spacedim>
+ void interpolation_difference(const DoFHandler<dim,spacedim> &dof1,
+ const ConstraintMatrix &constraints1,
+ const InVector &u1,
+ const DoFHandler<dim,spacedim> &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<n2; ++i)
- {
- u2(dofs[i])+=u2_local(i);
- }
- ++cell1;
- ++cell2;
- }
-}
+ template <int dim, class InVector, class OutVector, int spacedim>
+ void project_dg(const DoFHandler<dim,spacedim> &dof1,
+ const InVector &u1,
+ const DoFHandler<dim,spacedim> &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<dim,spacedim>::active_cell_iterator cell1 = dof1.begin_active();
+ typename DoFHandler<dim,spacedim>::active_cell_iterator cell2 = dof2.begin_active();
+ typename DoFHandler<dim,spacedim>::active_cell_iterator end = dof2.end();
-template <int dim, class InVector, class OutVector, int spacedim>
-void FETools::extrapolate(const DoFHandler<dim,spacedim> &dof1,
- const InVector &u1,
- const DoFHandler<dim,spacedim> &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<double> u1_local(n1);
+ Vector<double> u2_local(n2);
+ std::vector<unsigned int> dofs(n2);
+ FullMatrix<double> matrix(n2,n1);
+ get_projection_matrix(dof1.get_fe(), dof2.get_fe(), matrix);
-template <int dim, class InVector, class OutVector, int spacedim>
-void FETools::extrapolate(const DoFHandler<dim,spacedim> &dof1,
- const InVector &u1,
- const DoFHandler<dim,spacedim> &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<n2; ++i)
+ {
+ u2(dofs[i])+=u2_local(i);
+ }
- OutVector u3;
- u3.reinit(u2);
- interpolate(dof1, u1, dof2, constraints, u3);
+ ++cell1;
+ ++cell2;
+ }
+ }
- const unsigned int dofs_per_cell = dof2.get_fe().dofs_per_cell;
- Vector<typename OutVector::value_type> 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 <int dim, class InVector, class OutVector, int spacedim>
+ void extrapolate(const DoFHandler<dim,spacedim> &dof1,
+ const InVector &u1,
+ const DoFHandler<dim,spacedim> &dof2,
+ OutVector &u2)
{
- typename DoFHandler<dim,spacedim>::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<dof1.get_tria().n_levels()-1; ++level)
- {
- typename DoFHandler<dim,spacedim>::cell_iterator cell=dof2.begin(level),
- endc=dof2.end(level);
+
+ template <int dim, class InVector, class OutVector, int spacedim>
+ void extrapolate(const DoFHandler<dim,spacedim> &dof1,
+ const InVector &u1,
+ const DoFHandler<dim,spacedim> &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<typename OutVector::value_type> 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<dim,spacedim>::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_n<cell->n_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<dof1.get_tria().n_levels()-1; ++level)
+ {
+ typename DoFHandler<dim,spacedim>::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_n<cell->n_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 <int dim, int spacedim>
-void
-FETools::add_fe_name(const std::string& parameter_name,
- const FEFactoryBase<dim,spacedim>* 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 <int dim, int spacedim>
+ void
+ add_fe_name(const std::string& parameter_name,
+ const FEFactoryBase<dim,spacedim>* 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 <int dim, int spacedim>
- FiniteElement<dim,spacedim>*
- 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<FiniteElement<dim,spacedim>*> base_fes;
- std::vector<unsigned int> 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<dim,spacedim> (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<int,unsigned int> 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<dim,spacedim> *system_element = 0;
- switch (base_fes.size())
- {
- case 1:
+ template <int dim, int spacedim>
+ FiniteElement<dim,spacedim>*
+ 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<FiniteElement<dim,spacedim>*> base_fes;
+ std::vector<unsigned int> 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<dim>(*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<dim,spacedim> (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<int,unsigned int> 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<dim,spacedim> *system_element = 0;
+ switch (base_fes.size())
{
- system_element = new FESystem<dim>(*base_fes[0],
- base_multiplicities[0],
- *base_fes[1],
- base_multiplicities[1]);
- break;
- }
+ case 1:
+ {
+ system_element = new FESystem<dim>(*base_fes[0],
+ base_multiplicities[0]);
+ break;
+ }
- case 3:
- {
- system_element = new FESystem<dim>(*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<dim>(*base_fes[0],
+ base_multiplicities[0],
+ *base_fes[1],
+ base_multiplicities[1]);
+ break;
+ }
- default:
- AssertThrow (false, ExcNotImplemented());
- }
+ case 3:
+ {
+ system_element = new FESystem<dim>(*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<base_fes.size(); ++i)
- delete base_fes[i];
+ default:
+ AssertThrow (false, ExcNotImplemented());
+ }
- // finally return our
- // findings
- // Add the closing ']' to
- // the length
- return system_element;
+ // now we don't need the
+ // list of base elements
+ // any more
+ for (unsigned int i=0; i<base_fes.size(); ++i)
+ delete base_fes[i];
- }
- 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<base_fes.size(); ++i)
- delete base_fes[i];
- throw;
- }
+ // finally return our
+ // findings
+ // Add the closing ']' to
+ // the length
+ return system_element;
- // this is a place where we
- // should really never get,
- // since above we have either
- // returned from the
- // try-clause, or have
- // re-thrown in the catch
- // clause. check that we
- // never get here
- Assert (false, ExcInternalError());
- }
- else
- {
- // Make sure no other thread
- // is just adding an element
- Threads::ThreadMutex::ScopedLock lock (fe_name_map_lock);
-
- AssertThrow (fe_name_map.find(name_part) != fe_name_map.end(),
- FETools::ExcInvalidFEName(name));
- // Now, just the (degree)
- // or (Quadrature<1>(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<int,unsigned int> 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<int,unsigned int> 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<base_fes.size(); ++i)
+ delete base_fes[i];
+ throw;
+ }
+
+ // this is a place where we
+ // should really never get,
+ // since above we have either
+ // returned from the
+ // try-clause, or have
+ // re-thrown in the catch
+ // clause. check that we
+ // never get here
+ Assert (false, ExcInternalError());
+ }
+ else
+ {
+ // Make sure no other thread
+ // is just adding an element
+ Threads::ThreadMutex::ScopedLock lock (fe_name_map_lock);
+
+ AssertThrow (fe_name_map.find(name_part) != fe_name_map.end(),
+ ExcInvalidFEName(name));
+ // Now, just the (degree)
+ // or (Quadrature<1>(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<int,unsigned int> 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<int,unsigned int> 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 <int dim>
-FiniteElement<dim, dim> *
-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 <int dim>
+ FiniteElement<dim, dim> *
+ 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<dim,dim> *fe = internal::get_fe_from_name<dim,dim> (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<dim,dim> *fe = internal::get_fe_from_name<dim,dim> (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 <int dim>
// FiniteElement<dim> *
-// 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<dim,dim>(parameter_name);
// }
-template <int dim, int spacedim>
-void
-FETools::
-compute_projection_from_quadrature_points_matrix (const FiniteElement<dim,spacedim> &fe,
- const Quadrature<dim> &lhs_quadrature,
- const Quadrature<dim> &rhs_quadrature,
- FullMatrix<double> &X)
-{
- Assert (fe.n_components() == 1, ExcNotImplemented());
-
- // first build the matrices M and Q
- // described in the documentation
- FullMatrix<double> M (fe.dofs_per_cell, fe.dofs_per_cell);
- FullMatrix<double> Q (fe.dofs_per_cell, rhs_quadrature.size());
-
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
- for (unsigned int q=0; q<lhs_quadrature.size(); ++q)
- M(i,j) += fe.shape_value (i, lhs_quadrature.point(q)) *
- fe.shape_value (j, lhs_quadrature.point(q)) *
- lhs_quadrature.weight(q);
-
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- for (unsigned int q=0; q<rhs_quadrature.size(); ++q)
- Q(i,q) += fe.shape_value (i, rhs_quadrature.point(q)) *
- rhs_quadrature.weight(q);
-
- // then invert M
- FullMatrix<double> 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 <int dim, int spacedim>
+ void
+ compute_projection_from_quadrature_points_matrix (const FiniteElement<dim,spacedim> &fe,
+ const Quadrature<dim> &lhs_quadrature,
+ const Quadrature<dim> &rhs_quadrature,
+ FullMatrix<double> &X)
+ {
+ Assert (fe.n_components() == 1, ExcNotImplemented());
+ // first build the matrices M and Q
+ // described in the documentation
+ FullMatrix<double> M (fe.dofs_per_cell, fe.dofs_per_cell);
+ FullMatrix<double> Q (fe.dofs_per_cell, rhs_quadrature.size());
-template <int dim, int spacedim>
-void
-FETools::
-compute_interpolation_to_quadrature_points_matrix (const FiniteElement<dim,spacedim> &fe,
- const Quadrature<dim> &quadrature,
- FullMatrix<double> &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<fe.dofs_per_cell; ++i)
+ for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ for (unsigned int q=0; q<lhs_quadrature.size(); ++q)
+ M(i,j) += fe.shape_value (i, lhs_quadrature.point(q)) *
+ fe.shape_value (j, lhs_quadrature.point(q)) *
+ lhs_quadrature.weight(q);
- for (unsigned int q=0; q<quadrature.size(); ++q)
for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- I_q(q,i) = fe.shape_value (i, quadrature.point(q));
-}
+ for (unsigned int q=0; q<rhs_quadrature.size(); ++q)
+ Q(i,q) += fe.shape_value (i, rhs_quadrature.point(q)) *
+ rhs_quadrature.weight(q);
+ // then invert M
+ FullMatrix<double> 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 <int dim>
-void
-FETools::compute_projection_from_quadrature_points(
- const FullMatrix<double> &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 <int dim, int spacedim>
+ void
+ compute_interpolation_to_quadrature_points_matrix (const FiniteElement<dim,spacedim> &fe,
+ const Quadrature<dim> &quadrature,
+ FullMatrix<double> &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<quadrature.size(); ++q)
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ I_q(q,i) = fe.shape_value (i, quadrature.point(q));
+ }
+
+
+
+ template <int dim>
+ void
+ compute_projection_from_quadrature_points(
+ const FullMatrix<double> &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<double> component_at_node(n_support_points);
- // component at the quadrature point
+ // component at the quadrature point
Vector<double> 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 <n_support_points; ++nn) {
- vector_of_tensors_at_nodes[nn][ii] = component_at_node(nn);
- }
+ // rewrite the projection of the components
+ // back into the vector of tensors
+ for (unsigned int nn =0; nn <n_support_points; ++nn) {
+ vector_of_tensors_at_nodes[nn][ii] = component_at_node(nn);
+ }
}
-}
+ }
-template <int dim>
-void
-FETools::compute_projection_from_quadrature_points(
- const FullMatrix<double> &projection_matrix,
- const std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_qp,
- std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_nodes)
-{
+ template <int dim>
+ void
+ compute_projection_from_quadrature_points(
+ const FullMatrix<double> &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<double> component_at_node(n_support_points);
- // component at the quadrature point
+ // component at the quadrature point
Vector<double> 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 <n_support_points; ++nn) {
- (vector_of_tensors_at_nodes[nn])[row][column] = component_at_node(nn);
- }
+ // rewrite the projection of the components back into the vector of tensors
+ for (unsigned int nn =0; nn <n_support_points; ++nn) {
+ (vector_of_tensors_at_nodes[nn])[row][column] = component_at_node(nn);
+ }
}
-}
+ }
-template <int dim, int spacedim>
-void
-FETools::
-compute_projection_from_face_quadrature_points_matrix (const FiniteElement<dim, spacedim> &fe,
- const Quadrature<dim-1> &lhs_quadrature,
- const Quadrature<dim-1> &rhs_quadrature,
- const typename DoFHandler<dim, spacedim>::active_cell_iterator & cell,
- unsigned int face,
- FullMatrix<double> &X)
-{
- Assert (fe.n_components() == 1, ExcNotImplemented());
- Assert (lhs_quadrature.size () > fe.degree, ExcNotGreaterThan (lhs_quadrature.size (), fe.degree));
+ template <int dim, int spacedim>
+ void
+ compute_projection_from_face_quadrature_points_matrix (const FiniteElement<dim, spacedim> &fe,
+ const Quadrature<dim-1> &lhs_quadrature,
+ const Quadrature<dim-1> &rhs_quadrature,
+ const typename DoFHandler<dim, spacedim>::active_cell_iterator & cell,
+ unsigned int face,
+ FullMatrix<double> &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<double> M (fe.dofs_per_cell, fe.dofs_per_cell);
- FullMatrix<double> Q (fe.dofs_per_cell, rhs_quadrature.size());
+ // build the matrices M and Q
+ // described in the documentation
+ FullMatrix<double> M (fe.dofs_per_cell, fe.dofs_per_cell);
+ FullMatrix<double> Q (fe.dofs_per_cell, rhs_quadrature.size());
- {
- // need an FEFaceValues object to evaluate shape function
- // values on the specified face.
- FEFaceValues <dim> 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 <dim> 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.dofs_per_cell; ++i)
+ for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+ for (unsigned int q=0; q<lhs_quadrature.size(); ++q)
+ M(i,j) += fe_face_values.shape_value (i, q) *
+ fe_face_values.shape_value (j, q) *
+ lhs_quadrature.weight(q);
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ {
+ M(i,i) = (M(i,i) == 0 ? 1 : M(i,i));
+ }
+ }
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
- for (unsigned int q=0; q<lhs_quadrature.size(); ++q)
- M(i,j) += fe_face_values.shape_value (i, q) *
- fe_face_values.shape_value (j, q) *
- lhs_quadrature.weight(q);
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
{
- M(i,i) = (M(i,i) == 0 ? 1 : M(i,i));
+ FEFaceValues <dim> 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<fe.dofs_per_cell; ++i)
+ for (unsigned int q=0; q<rhs_quadrature.size(); ++q)
+ Q(i,q) += fe_face_values.shape_value (i, q) *
+ rhs_quadrature.weight(q);
}
- }
+ // then invert M
+ FullMatrix<double> M_inverse (fe.dofs_per_cell, fe.dofs_per_cell);
+ M_inverse.invert (M);
- {
- FEFaceValues <dim> 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<fe.dofs_per_cell; ++i)
- for (unsigned int q=0; q<rhs_quadrature.size(); ++q)
- Q(i,q) += fe_face_values.shape_value (i, q) *
- rhs_quadrature.weight(q);
+ Assert (X.m() == fe.dofs_per_cell, ExcInternalError());
+ Assert (X.n() == rhs_quadrature.size(), ExcInternalError());
}
- // then invert M
- FullMatrix<double> 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<deal_II_dimension>;
-
-template
-void FETools::compute_node_matrix(
- FullMatrix<double>&,
- const FiniteElement<deal_II_dimension>&);
-
-template
-void FETools::compute_component_wise(
- const FiniteElement<deal_II_dimension>& element,
- std::vector<unsigned int>&, std::vector<std::vector<unsigned int> >&);
-template
-void FETools::compute_block_renumbering (
- const FiniteElement<deal_II_dimension>& element,
- std::vector<unsigned int>&, std::vector<unsigned int>&_indices, bool);
-template
-void FETools::get_interpolation_matrix<deal_II_dimension>
-(const FiniteElement<deal_II_dimension> &,
- const FiniteElement<deal_II_dimension> &,
- FullMatrix<double> &);
-template
-void FETools::get_back_interpolation_matrix<deal_II_dimension>
-(const FiniteElement<deal_II_dimension> &,
- const FiniteElement<deal_II_dimension> &,
- FullMatrix<double> &);
-template
-void FETools::get_interpolation_difference_matrix<deal_II_dimension>
-(const FiniteElement<deal_II_dimension> &,
- const FiniteElement<deal_II_dimension> &,
- FullMatrix<double> &);
-template
-void FETools::get_interpolation_matrix<deal_II_dimension>
-(const FiniteElement<deal_II_dimension> &,
- const FiniteElement<deal_II_dimension> &,
- FullMatrix<float> &);
-template
-void FETools::get_back_interpolation_matrix<deal_II_dimension>
-(const FiniteElement<deal_II_dimension> &,
- const FiniteElement<deal_II_dimension> &,
- FullMatrix<float> &);
-template
-void FETools::get_interpolation_difference_matrix<deal_II_dimension>
-(const FiniteElement<deal_II_dimension> &,
- const FiniteElement<deal_II_dimension> &,
- FullMatrix<float> &);
-
-template
-void FETools::get_projection_matrix<deal_II_dimension>
-(const FiniteElement<deal_II_dimension> &,
- const FiniteElement<deal_II_dimension> &,
- FullMatrix<double> &);
-
-template
-void FETools::compute_embedding_matrices<deal_II_dimension>
-(const FiniteElement<deal_II_dimension> &, std::vector<std::vector<FullMatrix<double> > >&,bool);
-
-template
-void FETools::compute_face_embedding_matrices<deal_II_dimension,double>
-(const FiniteElement<deal_II_dimension> &, FullMatrix<double> (&matrices)[GeometryInfo<deal_II_dimension>::max_children_per_face],
- unsigned int, unsigned int);
-
-template
-void FETools::compute_projection_matrices<deal_II_dimension>
-(const FiniteElement<deal_II_dimension> &, std::vector<std::vector<FullMatrix<double> > >&, bool);
-
-template
-void FETools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const Vector<double> &,
- const DoFHandler<deal_II_dimension> &, Vector<double> &);
-template
-void FETools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const Vector<double> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- Vector<double> &);
+namespace FETools
+{
+ template class FEFactoryBase<deal_II_dimension>;
+
+ template
+ void compute_node_matrix(
+ FullMatrix<double>&,
+ const FiniteElement<deal_II_dimension>&);
+
+ template
+ void compute_component_wise(
+ const FiniteElement<deal_II_dimension>& element,
+ std::vector<unsigned int>&, std::vector<std::vector<unsigned int> >&);
+ template
+ void compute_block_renumbering (
+ const FiniteElement<deal_II_dimension>& element,
+ std::vector<unsigned int>&, std::vector<unsigned int>&_indices, bool);
+ template
+ void get_interpolation_matrix<deal_II_dimension>
+ (const FiniteElement<deal_II_dimension> &,
+ const FiniteElement<deal_II_dimension> &,
+ FullMatrix<double> &);
+ template
+ void get_back_interpolation_matrix<deal_II_dimension>
+ (const FiniteElement<deal_II_dimension> &,
+ const FiniteElement<deal_II_dimension> &,
+ FullMatrix<double> &);
+ template
+ void get_interpolation_difference_matrix<deal_II_dimension>
+ (const FiniteElement<deal_II_dimension> &,
+ const FiniteElement<deal_II_dimension> &,
+ FullMatrix<double> &);
+ template
+ void get_interpolation_matrix<deal_II_dimension>
+ (const FiniteElement<deal_II_dimension> &,
+ const FiniteElement<deal_II_dimension> &,
+ FullMatrix<float> &);
+ template
+ void get_back_interpolation_matrix<deal_II_dimension>
+ (const FiniteElement<deal_II_dimension> &,
+ const FiniteElement<deal_II_dimension> &,
+ FullMatrix<float> &);
+ template
+ void get_interpolation_difference_matrix<deal_II_dimension>
+ (const FiniteElement<deal_II_dimension> &,
+ const FiniteElement<deal_II_dimension> &,
+ FullMatrix<float> &);
+
+ template
+ void get_projection_matrix<deal_II_dimension>
+ (const FiniteElement<deal_II_dimension> &,
+ const FiniteElement<deal_II_dimension> &,
+ FullMatrix<double> &);
+
+ template
+ void compute_embedding_matrices<deal_II_dimension>
+ (const FiniteElement<deal_II_dimension> &, std::vector<std::vector<FullMatrix<double> > >&,bool);
+
+ template
+ void compute_face_embedding_matrices<deal_II_dimension,double>
+ (const FiniteElement<deal_II_dimension> &, FullMatrix<double> (&matrices)[GeometryInfo<deal_II_dimension>::max_children_per_face],
+ unsigned int, unsigned int);
+
+ template
+ void compute_projection_matrices<deal_II_dimension>
+ (const FiniteElement<deal_II_dimension> &, std::vector<std::vector<FullMatrix<double> > >&, bool);
+
+ template
+ void interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const Vector<double> &,
+ const DoFHandler<deal_II_dimension> &, Vector<double> &);
+ template
+ void interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const Vector<double> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ Vector<double> &);
#if deal_II_dimension != 3
-template
-void FETools::compute_block_renumbering (
- const FiniteElement<deal_II_dimension,deal_II_dimension+1>& element,
- std::vector<unsigned int>&, std::vector<unsigned int>&_indices, bool);
-template
-void FETools::interpolate<deal_II_dimension,deal_II_dimension+1>
-(const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const Vector<double> &,
- const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, Vector<double> &);
-template
-void FETools::interpolate<deal_II_dimension,deal_II_dimension+1>
-(const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const Vector<double> &,
- const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const ConstraintMatrix &,
- Vector<double> &);
+ template
+ void compute_block_renumbering (
+ const FiniteElement<deal_II_dimension,deal_II_dimension+1>& element,
+ std::vector<unsigned int>&, std::vector<unsigned int>&_indices, bool);
+ template
+ void interpolate<deal_II_dimension,deal_II_dimension+1>
+ (const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const Vector<double> &,
+ const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, Vector<double> &);
+ template
+ void interpolate<deal_II_dimension,deal_II_dimension+1>
+ (const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const Vector<double> &,
+ const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const ConstraintMatrix &,
+ Vector<double> &);
#endif
-template
-void FETools::back_interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const Vector<double> &,
- const FiniteElement<deal_II_dimension> &, Vector<double> &);
-template
-void FETools::back_interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- const Vector<double> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- Vector<double> &);
-template
-void FETools::interpolation_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const Vector<double> &,
- const FiniteElement<deal_II_dimension> &, Vector<double> &);
-template
-void FETools::interpolation_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- const Vector<double> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- Vector<double> &);
-template
-void FETools::project_dg<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const Vector<double> &,
- const DoFHandler<deal_II_dimension> &, Vector<double> &);
-template
-void FETools::extrapolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const Vector<double> &,
- const DoFHandler<deal_II_dimension> &, Vector<double> &);
-template
-void FETools::extrapolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const Vector<double> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- Vector<double> &);
-
-
-template
-void FETools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const Vector<float> &,
- const DoFHandler<deal_II_dimension> &, Vector<float> &);
-template
-void FETools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const Vector<float> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- Vector<float> &);
+ template
+ void back_interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const Vector<double> &,
+ const FiniteElement<deal_II_dimension> &, Vector<double> &);
+ template
+ void back_interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ const Vector<double> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ Vector<double> &);
+ template
+ void interpolation_difference<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const Vector<double> &,
+ const FiniteElement<deal_II_dimension> &, Vector<double> &);
+ template
+ void interpolation_difference<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ const Vector<double> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ Vector<double> &);
+ template
+ void project_dg<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const Vector<double> &,
+ const DoFHandler<deal_II_dimension> &, Vector<double> &);
+ template
+ void extrapolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const Vector<double> &,
+ const DoFHandler<deal_II_dimension> &, Vector<double> &);
+ template
+ void extrapolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const Vector<double> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ Vector<double> &);
+
+
+ template
+ void interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const Vector<float> &,
+ const DoFHandler<deal_II_dimension> &, Vector<float> &);
+ template
+ void interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const Vector<float> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ Vector<float> &);
#if deal_II_dimension != 3
-template
-void FETools::interpolate<deal_II_dimension,deal_II_dimension+1>
-(const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const Vector<float> &,
- const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, Vector<float> &);
-template
-void FETools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const Vector<float> &,
- const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const ConstraintMatrix &,
- Vector<float> &);
+ template
+ void interpolate<deal_II_dimension,deal_II_dimension+1>
+ (const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const Vector<float> &,
+ const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, Vector<float> &);
+ template
+ void interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const Vector<float> &,
+ const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const ConstraintMatrix &,
+ Vector<float> &);
#endif
-template
-void FETools::back_interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const Vector<float> &,
- const FiniteElement<deal_II_dimension> &, Vector<float> &);
-template
-void FETools::back_interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- const Vector<float> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- Vector<float> &);
-template
-void FETools::interpolation_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const Vector<float> &,
- const FiniteElement<deal_II_dimension> &, Vector<float> &);
-template
-void FETools::interpolation_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- const Vector<float> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- Vector<float> &);
-template
-void FETools::project_dg<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const Vector<float> &,
- const DoFHandler<deal_II_dimension> &, Vector<float> &);
-template
-void FETools::extrapolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const Vector<float> &,
- const DoFHandler<deal_II_dimension> &, Vector<float> &);
-template
-void FETools::extrapolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const Vector<float> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- Vector<float> &);
+ template
+ void back_interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const Vector<float> &,
+ const FiniteElement<deal_II_dimension> &, Vector<float> &);
+ template
+ void back_interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ const Vector<float> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ Vector<float> &);
+ template
+ void interpolation_difference<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const Vector<float> &,
+ const FiniteElement<deal_II_dimension> &, Vector<float> &);
+ template
+ void interpolation_difference<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ const Vector<float> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ Vector<float> &);
+ template
+ void project_dg<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const Vector<float> &,
+ const DoFHandler<deal_II_dimension> &, Vector<float> &);
+ template
+ void extrapolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const Vector<float> &,
+ const DoFHandler<deal_II_dimension> &, Vector<float> &);
+ template
+ void extrapolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const Vector<float> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ Vector<float> &);
#ifdef DEAL_II_USE_TRILINOS
-template
-void FETools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::Vector &,
- const DoFHandler<deal_II_dimension> &, TrilinosWrappers::Vector &);
-template
-void FETools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::Vector &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- TrilinosWrappers::Vector &);
-template
-void FETools::back_interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::Vector &,
- const FiniteElement<deal_II_dimension> &, TrilinosWrappers::Vector &);
-template
-void FETools::back_interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- const TrilinosWrappers::Vector &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- TrilinosWrappers::Vector &);
-template
-void FETools::interpolation_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::Vector &,
- const FiniteElement<deal_II_dimension> &, TrilinosWrappers::Vector &);
-template
-void FETools::interpolation_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- const TrilinosWrappers::Vector &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- TrilinosWrappers::Vector &);
-template
-void FETools::project_dg<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::Vector &,
- const DoFHandler<deal_II_dimension> &, TrilinosWrappers::Vector &);
-template
-void FETools::extrapolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::Vector &,
- const DoFHandler<deal_II_dimension> &, TrilinosWrappers::Vector &);
-template
-void FETools::extrapolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::Vector &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- TrilinosWrappers::Vector &);
+ template
+ void interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::Vector &,
+ const DoFHandler<deal_II_dimension> &, TrilinosWrappers::Vector &);
+ template
+ void interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::Vector &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ TrilinosWrappers::Vector &);
+ template
+ void back_interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::Vector &,
+ const FiniteElement<deal_II_dimension> &, TrilinosWrappers::Vector &);
+ template
+ void back_interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ const TrilinosWrappers::Vector &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ TrilinosWrappers::Vector &);
+ template
+ void interpolation_difference<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::Vector &,
+ const FiniteElement<deal_II_dimension> &, TrilinosWrappers::Vector &);
+ template
+ void interpolation_difference<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ const TrilinosWrappers::Vector &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ TrilinosWrappers::Vector &);
+ template
+ void project_dg<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::Vector &,
+ const DoFHandler<deal_II_dimension> &, TrilinosWrappers::Vector &);
+ template
+ void extrapolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::Vector &,
+ const DoFHandler<deal_II_dimension> &, TrilinosWrappers::Vector &);
+ template
+ void extrapolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::Vector &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ TrilinosWrappers::Vector &);
#endif
-template
-void FETools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
- const DoFHandler<deal_II_dimension> &, BlockVector<double> &);
-template
-void FETools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- BlockVector<double> &);
+ template
+ void interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
+ const DoFHandler<deal_II_dimension> &, BlockVector<double> &);
+ template
+ void interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ BlockVector<double> &);
#if deal_II_dimension != 3
-template
-void FETools::interpolate<deal_II_dimension,deal_II_dimension+1>
-(const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const BlockVector<double> &,
- const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, BlockVector<double> &);
-template
-void FETools::interpolate<deal_II_dimension,deal_II_dimension+1>
-(const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const BlockVector<double> &,
- const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const ConstraintMatrix &,
- BlockVector<double> &);
+ template
+ void interpolate<deal_II_dimension,deal_II_dimension+1>
+ (const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const BlockVector<double> &,
+ const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, BlockVector<double> &);
+ template
+ void interpolate<deal_II_dimension,deal_II_dimension+1>
+ (const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const BlockVector<double> &,
+ const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const ConstraintMatrix &,
+ BlockVector<double> &);
#endif
-template
-void FETools::back_interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
- const FiniteElement<deal_II_dimension> &, BlockVector<double> &);
-template
-void FETools::back_interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- const BlockVector<double> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- BlockVector<double> &);
-template
-void FETools::interpolation_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
- const FiniteElement<deal_II_dimension> &, BlockVector<double> &);
-template
-void FETools::interpolation_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- const BlockVector<double> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- BlockVector<double> &);
-template
-void FETools::project_dg<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
- const DoFHandler<deal_II_dimension> &, BlockVector<double> &);
-template
-void FETools::extrapolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
- const DoFHandler<deal_II_dimension> &, BlockVector<double> &);
-template
-void FETools::extrapolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- BlockVector<double> &);
-template
-void FETools::extrapolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
- const DoFHandler<deal_II_dimension> &, Vector<double> &);
-template
-void FETools::extrapolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- Vector<double> &);
-
-
-template
-void FETools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
- const DoFHandler<deal_II_dimension> &, BlockVector<float> &);
-template
-void FETools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- BlockVector<float> &);
+ template
+ void back_interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
+ const FiniteElement<deal_II_dimension> &, BlockVector<double> &);
+ template
+ void back_interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ const BlockVector<double> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ BlockVector<double> &);
+ template
+ void interpolation_difference<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
+ const FiniteElement<deal_II_dimension> &, BlockVector<double> &);
+ template
+ void interpolation_difference<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ const BlockVector<double> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ BlockVector<double> &);
+ template
+ void project_dg<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
+ const DoFHandler<deal_II_dimension> &, BlockVector<double> &);
+ template
+ void extrapolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
+ const DoFHandler<deal_II_dimension> &, BlockVector<double> &);
+ template
+ void extrapolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ BlockVector<double> &);
+ template
+ void extrapolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
+ const DoFHandler<deal_II_dimension> &, Vector<double> &);
+ template
+ void extrapolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ Vector<double> &);
+
+
+ template
+ void interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
+ const DoFHandler<deal_II_dimension> &, BlockVector<float> &);
+ template
+ void interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ BlockVector<float> &);
#if deal_II_dimension != 3
-template
-void FETools::interpolate<deal_II_dimension,deal_II_dimension+1>
-(const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const BlockVector<float> &,
- const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, BlockVector<float> &);
-template
-void FETools::interpolate<deal_II_dimension,deal_II_dimension+1>
-(const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const BlockVector<float> &,
- const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const ConstraintMatrix &,
- BlockVector<float> &);
+ template
+ void interpolate<deal_II_dimension,deal_II_dimension+1>
+ (const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const BlockVector<float> &,
+ const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, BlockVector<float> &);
+ template
+ void interpolate<deal_II_dimension,deal_II_dimension+1>
+ (const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const BlockVector<float> &,
+ const DoFHandler<deal_II_dimension,deal_II_dimension+1> &, const ConstraintMatrix &,
+ BlockVector<float> &);
#endif
-template
-void FETools::back_interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
- const FiniteElement<deal_II_dimension> &, BlockVector<float> &);
-template
-void FETools::back_interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- const BlockVector<float> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- BlockVector<float> &);
-template
-void FETools::interpolation_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
- const FiniteElement<deal_II_dimension> &, BlockVector<float> &);
-template
-void FETools::interpolation_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- const BlockVector<float> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- BlockVector<float> &);
-template
-void FETools::project_dg<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
- const DoFHandler<deal_II_dimension> &, BlockVector<float> &);
-template
-void FETools::extrapolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
- const DoFHandler<deal_II_dimension> &, BlockVector<float> &);
-template
-void FETools::extrapolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- BlockVector<float> &);
-template
-void FETools::extrapolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
- const DoFHandler<deal_II_dimension> &, Vector<float> &);
-template
-void FETools::extrapolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- Vector<float> &);
+ template
+ void back_interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
+ const FiniteElement<deal_II_dimension> &, BlockVector<float> &);
+ template
+ void back_interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ const BlockVector<float> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ BlockVector<float> &);
+ template
+ void interpolation_difference<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
+ const FiniteElement<deal_II_dimension> &, BlockVector<float> &);
+ template
+ void interpolation_difference<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ const BlockVector<float> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ BlockVector<float> &);
+ template
+ void project_dg<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
+ const DoFHandler<deal_II_dimension> &, BlockVector<float> &);
+ template
+ void extrapolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
+ const DoFHandler<deal_II_dimension> &, BlockVector<float> &);
+ template
+ void extrapolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ BlockVector<float> &);
+ template
+ void extrapolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
+ const DoFHandler<deal_II_dimension> &, Vector<float> &);
+ template
+ void extrapolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ Vector<float> &);
#ifdef DEAL_II_USE_TRILINOS
-template
-void FETools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::BlockVector &,
- const DoFHandler<deal_II_dimension> &, TrilinosWrappers::BlockVector &);
-template
-void FETools::interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::BlockVector &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- TrilinosWrappers::BlockVector &);
-template
-void FETools::back_interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::BlockVector &,
- const FiniteElement<deal_II_dimension> &, TrilinosWrappers::BlockVector &);
-template
-void FETools::back_interpolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- const TrilinosWrappers::BlockVector &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- TrilinosWrappers::BlockVector &);
-template
-void FETools::interpolation_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::BlockVector &,
- const FiniteElement<deal_II_dimension> &, TrilinosWrappers::BlockVector &);
-template
-void FETools::interpolation_difference<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- const TrilinosWrappers::BlockVector &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- TrilinosWrappers::BlockVector &);
-template
-void FETools::project_dg<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::BlockVector &,
- const DoFHandler<deal_II_dimension> &, TrilinosWrappers::BlockVector &);
-template
-void FETools::extrapolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::BlockVector &,
- const DoFHandler<deal_II_dimension> &, TrilinosWrappers::BlockVector &);
-template
-void FETools::extrapolate<deal_II_dimension>
-(const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::BlockVector &,
- const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- TrilinosWrappers::BlockVector &);
+ template
+ void interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::BlockVector &,
+ const DoFHandler<deal_II_dimension> &, TrilinosWrappers::BlockVector &);
+ template
+ void interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::BlockVector &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ TrilinosWrappers::BlockVector &);
+ template
+ void back_interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::BlockVector &,
+ const FiniteElement<deal_II_dimension> &, TrilinosWrappers::BlockVector &);
+ template
+ void back_interpolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ const TrilinosWrappers::BlockVector &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ TrilinosWrappers::BlockVector &);
+ template
+ void interpolation_difference<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::BlockVector &,
+ const FiniteElement<deal_II_dimension> &, TrilinosWrappers::BlockVector &);
+ template
+ void interpolation_difference<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ const TrilinosWrappers::BlockVector &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ TrilinosWrappers::BlockVector &);
+ template
+ void project_dg<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::BlockVector &,
+ const DoFHandler<deal_II_dimension> &, TrilinosWrappers::BlockVector &);
+ template
+ void extrapolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::BlockVector &,
+ const DoFHandler<deal_II_dimension> &, TrilinosWrappers::BlockVector &);
+ template
+ void extrapolate<deal_II_dimension>
+ (const DoFHandler<deal_II_dimension> &, const TrilinosWrappers::BlockVector &,
+ const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ TrilinosWrappers::BlockVector &);
#endif
-template
-void FETools::interpolate<deal_II_dimension>
-(const hp::DoFHandler<deal_II_dimension> &, const Vector<double> &,
- const hp::DoFHandler<deal_II_dimension> &, Vector<double> &);
-template
-void FETools::interpolate<deal_II_dimension>
-(const hp::DoFHandler<deal_II_dimension> &, const Vector<double> &,
- const hp::DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- Vector<double> &);
-template
-void FETools::interpolate<deal_II_dimension>
-(const hp::DoFHandler<deal_II_dimension> &, const Vector<float> &,
- const hp::DoFHandler<deal_II_dimension> &, Vector<float> &);
-template
-void FETools::interpolate<deal_II_dimension>
-(const hp::DoFHandler<deal_II_dimension> &, const Vector<float> &,
- const hp::DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
- Vector<float> &);
-
-
-template FiniteElement<deal_II_dimension,deal_II_dimension> *
-FETools::get_fe_from_name<deal_II_dimension> (const std::string &);
-
-
-template
-void FETools::add_fe_name<deal_II_dimension>(
- const std::string& name,
- const FEFactoryBase<deal_II_dimension>* factory);
-
-template
-void
-FETools::
-compute_projection_from_quadrature_points_matrix (const FiniteElement<deal_II_dimension> &fe,
- const Quadrature<deal_II_dimension> &lhs_quadrature,
- const Quadrature<deal_II_dimension> &rhs_quadrature,
- FullMatrix<double> &X);
+ template
+ void interpolate<deal_II_dimension>
+ (const hp::DoFHandler<deal_II_dimension> &, const Vector<double> &,
+ const hp::DoFHandler<deal_II_dimension> &, Vector<double> &);
+ template
+ void interpolate<deal_II_dimension>
+ (const hp::DoFHandler<deal_II_dimension> &, const Vector<double> &,
+ const hp::DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ Vector<double> &);
+ template
+ void interpolate<deal_II_dimension>
+ (const hp::DoFHandler<deal_II_dimension> &, const Vector<float> &,
+ const hp::DoFHandler<deal_II_dimension> &, Vector<float> &);
+ template
+ void interpolate<deal_II_dimension>
+ (const hp::DoFHandler<deal_II_dimension> &, const Vector<float> &,
+ const hp::DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
+ Vector<float> &);
+
+
+ template FiniteElement<deal_II_dimension,deal_II_dimension> *
+ get_fe_from_name<deal_II_dimension> (const std::string &);
+
+
+ template
+ void add_fe_name<deal_II_dimension>(
+ const std::string& name,
+ const FEFactoryBase<deal_II_dimension>* factory);
+
+ template
+ void
+ compute_projection_from_quadrature_points_matrix (const FiniteElement<deal_II_dimension> &fe,
+ const Quadrature<deal_II_dimension> &lhs_quadrature,
+ const Quadrature<deal_II_dimension> &rhs_quadrature,
+ FullMatrix<double> &X);
-template
-void
-FETools::
-compute_projection_from_quadrature_points(
- const FullMatrix<double> &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<double> &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<double> &projection_matrix,
- const std::vector<SymmetricTensor<2, deal_II_dimension> > &vector_of_tensors_at_qp,
- std::vector<SymmetricTensor<2, deal_II_dimension> > &vector_of_tensors_at_nodes);
+ template
+ void
+ compute_projection_from_quadrature_points(
+ const FullMatrix<double> &projection_matrix,
+ const std::vector<SymmetricTensor<2, deal_II_dimension> > &vector_of_tensors_at_qp,
+ std::vector<SymmetricTensor<2, deal_II_dimension> > &vector_of_tensors_at_nodes);
-template
-void
-FETools::
-compute_interpolation_to_quadrature_points_matrix (const FiniteElement<deal_II_dimension> &fe,
- const Quadrature<deal_II_dimension> &quadrature,
- FullMatrix<double> &I_q);
+ template
+ void
+ compute_interpolation_to_quadrature_points_matrix (const FiniteElement<deal_II_dimension> &fe,
+ const Quadrature<deal_II_dimension> &quadrature,
+ FullMatrix<double> &I_q);
#if deal_II_dimension != 1
-template
-void
-FETools::
-compute_projection_from_face_quadrature_points_matrix (const FiniteElement<deal_II_dimension> &fe,
- const Quadrature<deal_II_dimension-1> &lhs_quadrature,
- const Quadrature<deal_II_dimension-1> &rhs_quadrature,
- const DoFHandler<deal_II_dimension>::active_cell_iterator & cell,
- unsigned int face,
- FullMatrix<double> &X);
+ template
+ void
+ compute_projection_from_face_quadrature_points_matrix (const FiniteElement<deal_II_dimension> &fe,
+ const Quadrature<deal_II_dimension-1> &lhs_quadrature,
+ const Quadrature<deal_II_dimension-1> &rhs_quadrature,
+ const DoFHandler<deal_II_dimension>::active_cell_iterator & cell,
+ unsigned int face,
+ FullMatrix<double> &X);
#endif
+}
+
/*---------------------------- fe_tools.cc ---------------------------*/