]> https://gitweb.dealii.org/ - dealii.git/commitdiff
new computation routine for embedding matrices
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 8 Feb 2005 11:10:39 +0000 (11:10 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 8 Feb 2005 11:10:39 +0000 (11:10 +0000)
git-svn-id: https://svn.dealii.org/trunk@9901 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_tools.h
deal.II/deal.II/source/fe/fe_tools.cc
deal.II/doc/news/changes.html

index da4b702484ae5e403ec1e80169f8197b70eb007a..f1281424b2fa7023f8addb715b91dd0e79372316 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2000, 2001, 2002, 2003, 2004 by the deal.II authors
+//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -160,6 +160,36 @@ class FETools
                                      const FiniteElement<dim> &fe2,
                                      FullMatrix<number> &matrix);
 
+                                    /**
+                                     * Compute the embedding matrices
+                                     * from a coarse cell to
+                                     * 2<sup>dim</sup> child cells.
+                                     *
+                                     * 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.
+                                     *
+                                     * @arg fe The finite element
+                                     * class for which we compute the
+                                     * embedding matrices.
+                                     * @arg matrices A pointer to
+                                     * 2<sup>dim</sup> FullMatrix
+                                     * objects. This is the format
+                                     * used in FiniteElementBase,
+                                     * where we want to use ths
+                                     * function mostly.
+                                     */
+    template <int dim, typename number>
+    static void compute_embedding_matrices(const FiniteElement<dim> &fe,
+                                          FullMatrix<number>* matrices);
     
                                     /**
                                      * Gives the interpolation of a the
@@ -666,6 +696,14 @@ class FETools
     static
     std::pair<FiniteElement<dim> *, unsigned int>
     get_fe_from_name_aux (const std::string &name);    
+
+                                    /**
+                                     * Exception thrown if an
+                                     * embedding matrix was computed
+                                     * inaccurately.
+                                     */
+    DeclException1(ExcLeastSquaresError, double,
+                  << "Least squares fit leaves a gap of " << arg1);
 };
 
 /*@}*/
index 12ed9265aa591b1862aad7c1f6c4f9c64c83abab..0d2bf7d48d9669dd511c5d7827547a4285d50438 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2000, 2001, 2002, 2003, 2004 by the deal.II authors
+//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -13,6 +13,7 @@
 
 
 #include <base/quadrature_lib.h>
+#include <base/logstream.h>
 #include <lac/full_matrix.h>
 #include <lac/vector.h>
 #include <lac/block_vector.h>
@@ -30,6 +31,7 @@
 #include <fe/fe_raviart_thomas.h>
 #include <fe/fe_system.h>
 #include <fe/fe_values.h>
+#include <fe/mapping_cartesian.h>
 #include <dofs/dof_constraints.h>
 #include <dofs/dof_handler.h>
 #include <dofs/dof_accessor.h>
@@ -469,6 +471,117 @@ void FETools::get_projection_matrix (const FiniteElement<dim> &fe1,
 }
 
 
+template<int dim, typename number>
+void
+FETools::compute_embedding_matrices(const FiniteElement<dim>& fe,
+                                   FullMatrix<number>* matrices)
+{
+  const unsigned int nc = GeometryInfo<dim>::children_per_cell;
+  const unsigned int n  = fe.dofs_per_cell;
+  const unsigned int degree = fe.degree;
+  
+  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 two meshes, one with a
+                                   * single reference cell and the
+                                   * other refined, together with
+                                   * DoFHandler and finite elements.
+                                   */
+  Triangulation<dim> tr_coarse;
+  Triangulation<dim> tr_fine;
+  GridGenerator::hyper_cube (tr_coarse, 0, 1);
+  GridGenerator::hyper_cube (tr_fine, 0, 1);
+  tr_fine.refine_global(1);
+  DoFHandler<dim> dof_coarse(tr_coarse);
+  dof_coarse.distribute_dofs(fe);
+  DoFHandler<dim> dof_fine(tr_fine);
+  dof_fine.distribute_dofs(fe);
+
+  MappingCartesian<dim> mapping;
+  QGauss<dim> q_fine(degree+1);
+  const unsigned int nq = q_fine.n_quadrature_points;
+  
+  FEValues<dim> fine (mapping, fe, q_fine,
+                     update_q_points | update_JxW_values | update_values);
+  
+  typename DoFHandler<dim>::active_cell_iterator coarse_cell
+    = dof_coarse.begin_active();
+  typename DoFHandler<dim>::active_cell_iterator fine_cell;
+
+                                  /**
+                                   * We search for the polynomial on
+                                   * the small cell, being equal to
+                                   * the coarse polynomial in all
+                                   * quadrature points.
+                                   *
+                                   * First build the matrix for this
+                                   * least squares problem. This
+                                   * contains the values of the fine
+                                   * cell polynomials in the fine
+                                   * cell grid points.
+                                   *
+                                   * This matrix is the same for all
+                                   * children.
+                                   */
+  fine.reinit(dof_fine.begin_active());
+  FullMatrix<number> A(nq, n);
+  for (unsigned int i=0;i<nq;++i)
+    for (unsigned int j=0;j<n;++j)
+      A(i,j) = fine.shape_value(j,i);
+  
+  Vector<number> v_coarse(nq);
+  Vector<number> v_fine(n);
+  
+  unsigned int cell_number = 0;
+  for (fine_cell = dof_fine.begin_active();
+       fine_cell != dof_fine.end();
+       ++fine_cell, ++cell_number)
+    {
+      FullMatrix<number>& matrix = matrices[cell_number];
+      fine.reinit(fine_cell);
+      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);
+
+                                      // 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 k=0;k<nq;++k)
+           v_coarse(k) = coarse.shape_value(i,k);
+
+                                          // Copy the matrix and
+                                          // solve the least squares
+                                          // problem.
+         FullMatrix<number> B = A;
+         const double result = B.least_squares(v_fine, v_coarse);
+         Assert (result < 1.e-10, ExcLeastSquaresError(result));
+           
+                                          // Finally 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)
+           matrix(j,i) = v_fine(j);
+       }
+    }
+  
+}
+
 
 template <int dim,
           class InVector, class OutVector>
@@ -1659,6 +1772,10 @@ void FETools::get_projection_matrix<deal_II_dimension>
  const FiniteElement<deal_II_dimension> &,
  FullMatrix<double> &);
 
+template
+void FETools::compute_embedding_matrices<deal_II_dimension>
+(const FiniteElement<deal_II_dimension> &, FullMatrix<double>*);
+
 template
 void FETools::interpolate<deal_II_dimension>
 (const DoFHandler<deal_II_dimension> &, const Vector<double> &,
index 6d2014b9ee01cdadaa87b26cbd7c6bb5fc657685..39c7a298415584ba569e19ce7a8d142e04f03c4a 100644 (file)
@@ -56,6 +56,16 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p>
+       New: Function <code class="class">FETools</code>::<code
+       class="member">compute_embedding_matrices</code> allows for
+       automatic computation of embedding matrices under the sole
+       assumption that the coarse grid space is embedded into the fine
+       grid space. In particular, no interpolation points are required.
+       <br> 
+       (GK, 2005/02/08)
+       </p>
+
   <li> <p>
        Fixed: Several wrong assertions in the Raviart-Thomas finite element
        class have been fixed, that were triggered when integrating face terms.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.