]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
The FiniteElement class was extended by two methods, which create
authorkayser-herold <kayser-herold@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 27 Jul 2006 16:32:30 +0000 (16:32 +0000)
committerkayser-herold <kayser-herold@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 27 Jul 2006 16:32:30 +0000 (16:32 +0000)
interpolation matrices between element faces. These functions will
be needed for the hp FE method.

git-svn-id: https://svn.dealii.org/trunk@13450 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/fe/fe_base.h
deal.II/deal.II/include/fe/fe_q.h
deal.II/deal.II/source/dofs/dof_tools.cc
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_q.cc

index 7dd6811cdf7a9929e23eae1f27bf1fc02e8dea5a..37aaef223e4dffc0089aa50a876222f94e2aefd1 100644 (file)
@@ -896,7 +896,7 @@ class FiniteElement : public Subscriptor,
                                      */
     virtual void
     get_interpolation_matrix (const FiniteElement<dim> &source,
-                             FullMatrix<double>           &matrix) const;
+                             FullMatrix<double>       &matrix) const;
                                     //@}
 
                                     /**
@@ -904,6 +904,62 @@ class FiniteElement : public Subscriptor,
                                      * @{
                                      */
     
+
+                                    /**
+                                     * Return the matrix
+                                     * interpolating from a face of
+                                     * of one element to the face of
+                                     * the neighboring element. 
+                                     * The size of the matrix is
+                                     * then @p dofs_per_face times
+                                     * <tt>source.dofs_per_face</tt>.
+                                     *
+                                     * Derived elements will have to
+                                     * implement this function. They
+                                     * may only provide interpolation
+                                     * matrices for certain source
+                                     * finite elements, for example
+                                     * those from the same family. If
+                                     * they don't implement
+                                     * interpolation from a given
+                                     * element, then they must throw
+                                     * an exception of type
+                                     * FiniteElement<dim>::ExcInterpolationNotImplemented.
+                                     */
+    virtual void
+    get_face_interpolation_matrix (const FiniteElement<dim> &source,
+                                  FullMatrix<double>       &matrix) const;
+                                    //@}
+    
+
+                                    /**
+                                     * Return the matrix
+                                     * interpolating from a face of
+                                     * of one element to the subface of
+                                     * the neighboring element. 
+                                     * The size of the matrix is
+                                     * then @p dofs_per_face times
+                                     * <tt>source.dofs_per_face</tt>.
+                                     *
+                                     * Derived elements will have to
+                                     * implement this function. They
+                                     * may only provide interpolation
+                                     * matrices for certain source
+                                     * finite elements, for example
+                                     * those from the same family. If
+                                     * they don't implement
+                                     * interpolation from a given
+                                     * element, then they must throw
+                                     * an exception of type
+                                     * FiniteElement<dim>::ExcInterpolationNotImplemented.
+                                     */
+    virtual void
+    get_subface_interpolation_matrix (const FiniteElement<dim> &source,
+                                     const unsigned int        subface,
+                                     FullMatrix<double>       &matrix) const;
+                                    //@}
+    
+
                                     /**
                                      * If, on a vertex, several
                                      * finite elements are active,
index e48f7c0d2cb7f0d05d9e2f3c9638ea53bd301c8b..929528a24d262a4512342c3beaa84666c5eb4caa 100644 (file)
@@ -560,7 +560,9 @@ face_to_equivalent_cell_index (const unsigned int index) const
 }
 
 
-
+//TODO: This method produces results, that differ from those which
+// are returned bz "face_to_cell_index (index, 0)". This has to be
+// verified.
 template <>
 inline
 unsigned int
index 9ebe9dab1d8cd7694bb1907f168b2d7a39bd3b0e..a3e6915032cfd49ee43558e97fa6c3b2d93c07f0 100644 (file)
@@ -271,8 +271,64 @@ class FE_Q : public FE_Poly<TensorProductPolynomials<dim>,dim>
                                      */
     virtual void
     get_interpolation_matrix (const FiniteElement<dim> &source,
-                             FullMatrix<double>           &matrix) const;
+                             FullMatrix<double>       &matrix) const;
+
     
+                                    /**
+                                     * Return the matrix
+                                     * interpolating from a face of
+                                     * of one element to the face of
+                                     * the neighboring element. 
+                                     * The size of the matrix is
+                                     * then @p dofs_per_face times
+                                     * <tt>source.dofs_per_face</tt>.
+                                     *
+                                     * Derived elements will have to
+                                     * implement this function. They
+                                     * may only provide interpolation
+                                     * matrices for certain source
+                                     * finite elements, for example
+                                     * those from the same family. If
+                                     * they don't implement
+                                     * interpolation from a given
+                                     * element, then they must throw
+                                     * an exception of type
+                                     * FiniteElement<dim>::ExcInterpolationNotImplemented.
+                                     */
+    virtual void
+    get_face_interpolation_matrix (const FiniteElement<dim> &source,
+                                  FullMatrix<double>       &matrix) const;
+                                    //@}
+    
+
+                                    /**
+                                     * Return the matrix
+                                     * interpolating from a face of
+                                     * of one element to the face of
+                                     * the neighboring element. 
+                                     * The size of the matrix is
+                                     * then @p dofs_per_face times
+                                     * <tt>source.dofs_per_face</tt>.
+                                     *
+                                     * Derived elements will have to
+                                     * implement this function. They
+                                     * may only provide interpolation
+                                     * matrices for certain source
+                                     * finite elements, for example
+                                     * those from the same family. If
+                                     * they don't implement
+                                     * interpolation from a given
+                                     * element, then they must throw
+                                     * an exception of type
+                                     * FiniteElement<dim>::ExcInterpolationNotImplemented.
+                                     */
+    virtual void
+    get_subface_interpolation_matrix (const FiniteElement<dim> &source,
+                                     const unsigned int        subface,
+                                     FullMatrix<double>       &matrix) const;
+                                    //@}
+    
+
                                     /**
                                      * Check for non-zero values on a face.
                                      *
index 7f2f251e3910e19fac9871e936fef10ada719b92..59534efa5c4c26d05b3070c83b121b1c8cbc2ca3 100644 (file)
@@ -1704,6 +1704,8 @@ namespace internal
                    = this_face->child(child)->dof_index(dof, fe_index);
              Assert (next_index == dofs_on_children.size(),
                      ExcInternalError());
+
+// TODO[OKH]: FE::get_subface_interpolation_matrix ()
          
                                               // for each row in the constraint
                                               // matrix for this line:
@@ -1729,6 +1731,8 @@ namespace internal
              Assert (cell->face(face)
                      ->fe_index_is_active(cell->active_fe_index()) == true,
                      ExcInternalError());
+
+// TODO[OKH]: FE::get_face_interpolation_matrix ()
            } 
     }
 #endif
index 13b1c8ea87e052ebaaa20a94248ebd4120aeb0b8..dcc9cc824b4b316d7eac60d3cc85517abc0d0f82 100644 (file)
@@ -442,6 +442,39 @@ get_interpolation_matrix (const FiniteElement<dim> &,
 
 
 
+template <int dim>
+void
+FiniteElement<dim>::
+get_face_interpolation_matrix (const FiniteElement<dim> &,
+                              FullMatrix<double>           &) const
+{
+                                  // by default, no interpolation
+                                  // implemented. so throw exception,
+                                  // as documentation says
+  AssertThrow (false,
+               typename FiniteElement<dim>::
+               ExcInterpolationNotImplemented());
+}
+                                   
+
+                                   
+template <int dim>
+void
+FiniteElement<dim>::
+get_subface_interpolation_matrix (const FiniteElement<dim> &,
+                                 const unsigned int,
+                                 FullMatrix<double>           &) const
+{
+                                  // by default, no interpolation
+                                  // implemented. so throw exception,
+                                  // as documentation says
+  AssertThrow (false,
+               typename FiniteElement<dim>::
+               ExcInterpolationNotImplemented());
+}
+                                   
+
+
 template <int dim>
 std::vector<std::pair<unsigned int, unsigned int> >
 FiniteElement<dim>::
index 3b2fbf8f0a6989944996de34d7e127132b9cb2ce..d25befe7a4fa2491201435a8d7a8efdc8fbac510 100644 (file)
@@ -324,6 +324,167 @@ get_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
 }
 
 
+template <int dim>
+void
+FE_Q<dim>::
+get_face_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
+                              FullMatrix<double>       &interpolation_matrix) const
+{
+                                  // this is only implemented, if the
+                                  // source FE is also a
+                                  // Q element
+  AssertThrow ((x_source_fe.get_name().find ("FE_Q<") == 0)
+               ||
+               (dynamic_cast<const FE_Q<dim>*>(&x_source_fe) != 0),
+               typename FiniteElement<dim>::
+               ExcInterpolationNotImplemented());
+  
+                                  // ok, source is a Q element, so
+                                  // we will be able to do the work
+  const FE_Q<dim> &source_fe
+    = dynamic_cast<const FE_Q<dim>&>(x_source_fe);
+
+  Assert (interpolation_matrix.m() == this->dofs_per_face,
+         ExcDimensionMismatch (interpolation_matrix.m(),
+                               this->dofs_per_face));
+  Assert (interpolation_matrix.n() == source_fe.dofs_per_face,
+         ExcDimensionMismatch (interpolation_matrix.m(),
+                               source_fe.dofs_per_face));
+
+  const std::vector<unsigned int> &index_map=
+    this->poly_space.get_numbering();
+  
+
+                                       // generate a point on this
+                                       // cell and evaluate the
+                                       // shape functions there
+  Quadrature<dim-1> quad_face_support (source_fe.get_unit_face_support_points ());
+
+
+                                  // compute the interpolation
+                                  // matrix by simply taking the
+                                   // value at the support points.
+  for (unsigned int i=0; i<source_fe.dofs_per_face; ++i)
+    {
+      //TODO: Verify that all faces are the same with respect to
+      // these support points. Furthermore, check if something has to
+      // be done for the face orientation flag in 3D.
+      Point<dim> p = QProjector<dim>::project_to_face (quad_face_support, 0).point (i);
+
+      for (unsigned int j=0; j<this->dofs_per_face; ++j)
+       interpolation_matrix(j,i) = this->shape_value (this->face_to_cell_index(j, 0), p);
+    }
+
+                                   // cut off very small values
+  for (unsigned int i=0; i<this->dofs_per_face; ++i)
+    for (unsigned int j=0; j<source_fe.dofs_per_face; ++j)
+      if (std::fabs(interpolation_matrix(i,j)) < 1e-15)
+        interpolation_matrix(i,j) = 0.;
+
+                                  // make sure that the column sum of
+                                  // each of the matrices is 1 at
+                                  // this point. this must be so
+                                  // since the shape functions sum up
+                                  // to 1
+  for (unsigned int j=0; j<source_fe.dofs_per_face; ++j)
+    {
+      double sum = 0.;
+
+      for (unsigned int i=0; i<this->dofs_per_face; ++i)
+        sum += interpolation_matrix(i,j);
+
+      Assert (std::fabs(sum-1) < 2e-14*this->degree*(dim-1),
+              ExcInternalError());
+    }
+}
+
+
+template <int dim>
+void
+FE_Q<dim>::
+get_subface_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
+                                 const unsigned int subface,
+                                 FullMatrix<double>           &interpolation_matrix) const
+{
+                                  // this is only implemented, if the
+                                  // source FE is also a
+                                  // Q element
+  AssertThrow ((x_source_fe.get_name().find ("FE_Q<") == 0)
+               ||
+               (dynamic_cast<const FE_Q<dim>*>(&x_source_fe) != 0),
+               typename FiniteElement<dim>::
+               ExcInterpolationNotImplemented());
+  
+                                  // ok, source is a Q element, so
+                                  // we will be able to do the work
+  const FE_Q<dim> &source_fe
+    = dynamic_cast<const FE_Q<dim>&>(x_source_fe);
+
+  Assert (interpolation_matrix.m() == this->dofs_per_cell,
+         ExcDimensionMismatch (interpolation_matrix.m(),
+                               this->dofs_per_cell));
+  Assert (interpolation_matrix.n() == source_fe.dofs_per_cell,
+         ExcDimensionMismatch (interpolation_matrix.m(),
+                               source_fe.dofs_per_cell));
+
+  const std::vector<unsigned int> &index_map=
+    this->poly_space.get_numbering();
+  
+                                  // compute the interpolation
+                                  // matrices in much the same way as
+                                  // we do for the embedding matrices
+                                  // from mother to child.
+  FullMatrix<double> cell_interpolation (this->dofs_per_cell,
+                                        this->dofs_per_cell);
+  FullMatrix<double> source_interpolation (this->dofs_per_cell,
+                                          source_fe.dofs_per_cell);
+  FullMatrix<double> tmp (this->dofs_per_cell,
+                         source_fe.dofs_per_cell);
+  for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+    {
+                                       // generate a point on this
+                                       // cell and evaluate the
+                                       // shape functions there
+      const Point<dim>
+       p = FE_Q_Helper::generate_unit_point (index_map[j], this->dofs_per_cell,
+                                             FE_Q_Helper::int2type<dim>());
+      for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+        cell_interpolation(j,i) = this->poly_space.compute_value (i, p);
+
+      for (unsigned int i=0; i<source_fe.dofs_per_cell; ++i)
+        source_interpolation(j,i) = source_fe.poly_space.compute_value (i, p);
+    }
+
+                                   // then compute the
+                                   // interpolation matrix matrix
+                                   // for this coordinate
+                                   // direction
+  cell_interpolation.gauss_jordan ();
+  cell_interpolation.mmult (interpolation_matrix,
+                            source_interpolation);
+
+                                   // cut off very small values
+  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+    for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
+      if (std::fabs(interpolation_matrix(i,j)) < 1e-15)
+        interpolation_matrix(i,j) = 0.;
+
+                                  // make sure that the row sum of
+                                  // each of the matrices is 1 at
+                                  // this point. this must be so
+                                  // since the shape functions sum up
+                                  // to 1
+  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+    {
+      double sum = 0.;
+      for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
+        sum += interpolation_matrix(i,j);
+
+      Assert (std::fabs(sum-1) < 2e-14*this->degree*dim,
+              ExcInternalError());
+    }
+}
+
 
 template <int dim>
 std::vector<std::pair<unsigned int, unsigned int> >

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.