From 901c10cd02eac8c1db5adfdc5d30f2826ff36673 Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 20 Jun 2005 21:26:48 +0000 Subject: [PATCH] compute_projection_from_quadrature_points_matrix compute_interpolation_to_quadrature_points_matrix git-svn-id: https://svn.dealii.org/trunk@10899 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_tools.h | 117 +++++++++++++++++++++++++- deal.II/deal.II/source/fe/fe_tools.cc | 75 +++++++++++++++++ 2 files changed, 191 insertions(+), 1 deletion(-) diff --git a/deal.II/deal.II/include/fe/fe_tools.h b/deal.II/deal.II/include/fe/fe_tools.h index 4c76d2e469..f694d7b266 100644 --- a/deal.II/deal.II/include/fe/fe_tools.h +++ b/deal.II/deal.II/include/fe/fe_tools.h @@ -48,7 +48,7 @@ class ConstraintMatrix; * $id-I_h$ that is needed for evaluating $(id-I_h)z$ for e.g. the * dual solution $z$. * - * @author Ralf Hartmann, 2000; Wolfgang Bangerth, 2003, Guido Kanschat, 2000, 2004 + * @author Ralf Hartmann, 2000; Wolfgang Bangerth, 2003, 2005, Guido Kanschat, 2000, 2004 */ class FETools { @@ -666,6 +666,121 @@ class FETools FiniteElement * get_fe_from_name (const std::string &name); + + /** + * Projects scalar data defined in + * quadrature points to a finite element + * space on a single cell. + * + * What this function does is the + * following: assume that there is scalar + * data uq, 0 <= q < + * Q:=quadrature.n_quadrature_points + * defined at the quadrature points of a + * cell, with the points defined by the + * given rhs_quadrature + * object. We may then want to ask for + * that finite element function (on a + * single cell) vh in + * the finite-dimensional space defined + * by the given FE object that is the + * projection of u in the + * following sense: + * + * Usually, the projection + * vh is that + * function that satisfies + * (vh,w)=(u,w) for + * all discrete test functions + * w. In the present case, we + * can't evaluate the right hand side, + * since u is only defined in + * the quadrature points given by + * rhs_quadrature, so we replace + * it by a quadrature + * approximation. Likewise, the left hand + * side is approximated using the + * lhs_quadrature object; if + * this quadrature object is chosen + * appropriately, then the integration of + * the left hand side can be done + * exactly, without any + * approximation. The use of different + * quadrature objects is necessary if the + * quadrature object for the right hand + * side has too few quadrature points -- + * for example, if data q is + * only defined at the cell center, then + * the corresponding one-point quadrature + * formula is obviously insufficient to + * approximate the scalar product on the + * left hand side by a definite form. + * + * After these quadrature approximations, + * we end up with a nodal representation + * Vh of + * vh that satisfies + * the following system of linear + * equations: M Vh = Q + * U, where + * Mij=(phi_i,phi_j) + * is the mass matrix approximated by + * lhs_quadrature, and + * Q is the matrix + * Qiq=phii(xq) + * wq where + * wq are quadrature + * weights; U is the vector of + * quadrature point data + * uq. + * + * In order to then get the nodal + * representation Vh + * of the projection of U, one + * computes Vh = X U, + * X=M-1 Q. The purpose + * of this function is to compute the + * matrix X and return it + * through the last argument of this + * function. + * + * Note that this function presently only + * supports scalar data. An extension of + * the mass matrix is of course trivial, + * but one has to define the order of + * data in the vector U if it + * contains vector valued data in all + * quadrature points. + * + * A use for this function is described + * in the introduction to the step-18 + * example program. + * + * The opposite of this function, + * interpolation of a finite element + * function onto quadrature points is + * essentially what the + * FEValues::get_function_values + * functions do; to make things a little + * simpler, the + * FETools::compute_interpolation_to_quadrature_points_matrix + * provides the matrix form of this. + */ + template + static + void + compute_projection_from_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature, + FullMatrix &X); + + template + static + void + compute_interpolation_to_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &quadrature, + FullMatrix &I_q); + /** * Exception */ diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index ed4254a026..3bba9f9d9e 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -1416,6 +1416,66 @@ FETools::get_fe_from_name_aux (const std::string &name) +template +void +FETools:: +compute_projection_from_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature, + FullMatrix &X) +{ + Assert (fe.n_components() == 1, ExcNotImplemented()); + + // first build the matrices M and Q + // described in the documentation + FullMatrix M (fe.dofs_per_cell, fe.dofs_per_cell); + FullMatrix Q (fe.dofs_per_cell, rhs_quadrature.n_quadrature_points); + + for (unsigned int i=0; i M_inverse (fe.dofs_per_cell, fe.dofs_per_cell); + M_inverse.invert (M); + + // finally compute the result + X.reinit (fe.dofs_per_cell, rhs_quadrature.n_quadrature_points); + M_inverse.mmult (X, Q); + + Assert (X.m() == fe.dofs_per_cell, ExcInternalError()); + Assert (X.n() == rhs_quadrature.n_quadrature_points, ExcInternalError()); +} + + + +template +void +FETools:: +compute_interpolation_to_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &quadrature, + FullMatrix &I_q) +{ + Assert (fe.n_components() == 1, ExcNotImplemented()); + Assert (I_q.m() == quadrature.n_quadrature_points, + ExcMessage ("Wrong matrix size")); + Assert (I_q.n() == fe.dofs_per_cell, ExcMessage ("Wrong matrix size")); + + for (unsigned int q=0; q * FETools::get_fe_from_name (const std::string &); +template +void +FETools:: +compute_projection_from_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature, + FullMatrix &X); + +template +void +FETools:: +compute_interpolation_to_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &quadrature, + FullMatrix &I_q); + /*---------------------------- fe_tools.cc ---------------------------*/ -- 2.39.5