From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Mon, 20 Jun 2005 21:26:48 +0000 (+0000)
Subject: compute_projection_from_quadrature_points_matrix
X-Git-Tag: v8.0.0~13657
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=94ea8af626b167842c6a98c21ef20549a8554801;p=dealii.git

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
---

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<dim> *
     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 <tt>u<sub>q</sub>, 0 <= q <
+                                      * Q:=quadrature.n_quadrature_points</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.
+                                      */
+    template <int dim>
+    static
+    void
+    compute_projection_from_quadrature_points_matrix (const FiniteElement<dim> &fe,
+                                                      const Quadrature<dim>    &lhs_quadrature,
+                                                      const Quadrature<dim>    &rhs_quadrature,
+                                                      FullMatrix<double>       &X);
+
+    template <int dim>
+    static
+    void
+    compute_interpolation_to_quadrature_points_matrix (const FiniteElement<dim> &fe,
+                                                       const Quadrature<dim>    &quadrature,
+                                                       FullMatrix<double>       &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 <int dim>
+void
+FETools::
+compute_projection_from_quadrature_points_matrix (const FiniteElement<dim> &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.n_quadrature_points);
+
+  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.n_quadrature_points; ++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.n_quadrature_points; ++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.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 <int dim>
+void
+FETools::
+compute_interpolation_to_quadrature_points_matrix (const FiniteElement<dim> &fe,
+                                                   const Quadrature<dim>    &quadrature,
+                                                   FullMatrix<double>       &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<quadrature.n_quadrature_points; ++q)
+    for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+      I_q(q,i) = fe.shape_value (i, quadrature.point(q));
+}
+
+
+
 
 /*-------------- Explicit Instantiations -------------------------------*/
 
@@ -1662,5 +1722,20 @@ template
 FiniteElement<deal_II_dimension> *
 FETools::get_fe_from_name<deal_II_dimension> (const std::string &);
 
+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
+FETools::
+compute_interpolation_to_quadrature_points_matrix (const FiniteElement<deal_II_dimension> &fe,
+                                                   const Quadrature<deal_II_dimension>    &quadrature,
+                                                   FullMatrix<double>       &I_q);
+
 
 /*----------------------------   fe_tools.cc     ---------------------------*/