]> https://gitweb.dealii.org/ - dealii.git/commitdiff
new function FETools::compute_face_embedding_matrices
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 12 Jul 2005 22:11:29 +0000 (22:11 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 12 Jul 2005 22:11:29 +0000 (22:11 +0000)
git-svn-id: https://svn.dealii.org/trunk@11135 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_tools.h
deal.II/deal.II/source/fe/fe_tools.cc

index d504f7fdd96ae8a2a0069d78170d4af81e3b4cfc..e79f2194dfd0c57f950c314eca0eae80ab61b6e5 100644 (file)
@@ -185,11 +185,11 @@ class FETools
                                      * the finite element spaces are
                                      * actually nested.
                                      *
-                                     * @arg fe The finite element
+                                     * @param fe The finite element
                                      * class for which we compute the
                                      * embedding matrices.
-                                     * @arg matrices A pointer to
-                                     * 2<sup>dim</sup> FullMatrix
+                                     * @param matrices A pointer to
+                                     * <i>2<sup>dim</sup></i> FullMatrix
                                      * objects. This is the format
                                      * used in FiniteElementBase,
                                      * where we want to use ths
@@ -198,7 +198,34 @@ class FETools
     template <int dim, typename number>
     static void compute_embedding_matrices(const FiniteElement<dim> &fe,
                                           FullMatrix<number>* matrices);
-    
+
+                                    /**
+                                     * 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>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.
+                                     */
+    template<int dim, typename number>
+    static void compute_face_embedding_matrices(const FiniteElement<dim>& fe,
+                                               FullMatrix<number>* matrices,
+                                               unsigned int face_coarse,
+                                               unsigned int face_fine);
+
                                     /**
                                      * Compute the
                                      * <i>L<sup>2</sup></i>-projection
index 41e682a1cd11c18b86cc755e349aa8a38cb3d86d..e6285240b3a1c9c908a6de224db4403d19e028ee 100644 (file)
@@ -13,6 +13,7 @@
 
 
 #include <base/quadrature_lib.h>
+#include <base/qprojector.h>
 #include <base/logstream.h>
 #include <lac/full_matrix.h>
 #include <lac/householder.h>
@@ -586,6 +587,175 @@ FETools::compute_embedding_matrices(const FiniteElement<dim>& fe,
 }
 
 
+// This function is tested by tests/fe/internals, since it produces the matrices printed there
+
+//TODO:[GK] Is this correct for vector valued?
+template<int dim, typename number>
+void
+FETools::compute_face_embedding_matrices(const FiniteElement<dim>& fe,
+                                        FullMatrix<number>* matrices,
+                                        unsigned int face_coarse,
+                                        unsigned int face_fine)
+{
+  const unsigned int nc = GeometryInfo<dim>::subfaces_per_face;
+  const unsigned int n  = fe.dofs_per_face;
+  const unsigned int nd = fe.n_components();
+  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 a meshes, one with a single
+                                   // reference cell and refine it once
+  Triangulation<dim> 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.n_quadrature_points;
+
+                                  // 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());
+  
+  FEValues<dim> fine (mapping, fe, q_fine,
+                     update_q_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.
+                                   
+                                  // 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(tria.begin_active());
+  FullMatrix<number> A(nq*nd, n);
+  for (unsigned int d=0;d<nd;++d)
+    for (unsigned int k=0;k<nq;++k)
+      for (unsigned int j=0;j<n;++j)
+       A(k*nd+d,j) = fine.shape_value_component(face_f_dofs[j],k,d);
+
+  Householder<double> H(A);
+  
+  Vector<number> v_coarse(nq*nd);
+  Vector<number> v_fine(n);
+  
+  
+  
+  for (unsigned int cell_number = 0; cell_number < GeometryInfo<dim>::subfaces_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);
+      
+      typename Triangulation<dim>::active_cell_iterator fine_cell
+       = tria.begin(0)->child(GeometryInfo<dim>::child_cell_on_face(face_coarse,
+                                                                   cell_number));
+      fine.reinit(fine_cell);
+      coarse.reinit(tria.begin(0));
+      
+      FullMatrix<double> &this_matrix = matrices[cell_number];
+      
+                                      // 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 function is tested by tests/fe/internals, since it produces the matrices printed there
 template<int dim, typename number>
 void
@@ -1553,6 +1723,11 @@ template
 void FETools::compute_embedding_matrices<deal_II_dimension>
 (const FiniteElement<deal_II_dimension> &, FullMatrix<double>*);
 
+template
+void FETools::compute_face_embedding_matrices<deal_II_dimension>
+(const FiniteElement<deal_II_dimension> &, FullMatrix<double>*,
+ unsigned int, unsigned int);
+
 template
 void FETools::compute_projection_matrices<deal_II_dimension>
 (const FiniteElement<deal_II_dimension> &, FullMatrix<double>*);

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.