]> https://gitweb.dealii.org/ - dealii.git/commitdiff
implement Raviart-Thomas basis functions for standard RT interpolation
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 25 Oct 2005 08:19:16 +0000 (08:19 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 25 Oct 2005 08:19:16 +0000 (08:19 +0000)
git-svn-id: https://svn.dealii.org/trunk@11657 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_raviart_thomas.h
deal.II/deal.II/source/fe/fe_raviart_thomas.cc
deal.II/deal.II/source/fe/fe_raviart_thomas_nodal.cc

index 4b08664d9f260c540e02530e921bff1d8d4812bf..70f020661159d4ad0c80b75085640bfddda88369 100644 (file)
@@ -136,7 +136,9 @@ template <int dim> class MappingQ;
  * @author Wolfgang Bangerth, 2003
  */
 template <int dim>
-class FE_RaviartThomas : public FiniteElement<dim>
+class FE_RaviartThomas
+  :
+  public FE_PolyTensor<PolynomialsRaviartThomas<dim>, dim>
 {
   public:
                                     /**
@@ -156,92 +158,6 @@ class FE_RaviartThomas : public FiniteElement<dim>
                                      */
     virtual std::string get_name () const;
 
-                                    /**
-                                     * Return the value of the
-                                     * @p componentth vector
-                                     * component of the @p ith shape
-                                     * function at the point
-                                     * @p p. See the
-                                     * FiniteElement base
-                                     * class for more information
-                                     * about the semantics of this
-                                     * function.
-                                     */
-    virtual double shape_value_component (const unsigned int i,
-                                         const Point<dim> &p,
-                                         const unsigned int component) const;
-
-                                    /**
-                                     * Return the gradient of the
-                                     * @p componentth vector
-                                     * component of the @p ith shape
-                                     * function at the point
-                                     * @p p. See the
-                                     * FiniteElement base
-                                     * class for more information
-                                     * about the semantics of this
-                                     * function.
-                                     */
-    virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
-                                               const Point<dim> &p,
-                                               const unsigned int component) const;
-
-                                    /**
-                                     * Return the second derivative
-                                     * of the @p componentth vector
-                                     * component of the @p ith shape
-                                     * function at the point
-                                     * @p p. See the
-                                     * FiniteElement base
-                                     * class for more information
-                                     * about the semantics of this
-                                     * function.
-                                     */
-    virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
-                                                    const Point<dim> &p,
-                                                    const unsigned int component) const;
-
-                                    /**
-                                     * Return the order
-                                     * of this finite element,
-                                     * i.e. the value passed to the
-                                     * constructor.
-                                     *
-                                     * Note that for this element,
-                                     * the order is actually one
-                                     * lower than the maximal
-                                     * polynomial degree, unlike for
-                                     * most of the other
-                                     * elements. For example, the RT0
-                                     * element as piecewise linear
-                                     * shape functions, even though
-                                     * the normal component of them
-                                     * is piecewise constant on each
-                                     * face (the latter being the
-                                     * property that defines the
-                                     * order of the element).
-                                     */
-    unsigned int get_degree () const;
-    
-                                    /**
-                                     * Return the matrix
-                                     * interpolating from the given
-                                     * finite element to the present
-                                     * one. The size of the matrix is
-                                     * then @p dofs_per_cell times
-                                     * <tt>source.dofs_per_cell</tt>.
-                                     *
-                                     * These matrices are only
-                                     * available if the source
-                                     * element is also a Raviart
-                                     * Thomas element. Otherwise, an
-                                     * exception of type
-                                     * FiniteElement<dim>::ExcInterpolationNotImplemented
-                                     * is thrown.
-                                     */
-    virtual void
-    get_interpolation_matrix (const FiniteElement<dim> &source,
-                             FullMatrix<double>           &matrix) const;
 
                                     /**
                                      * Number of base elements in a
@@ -282,66 +198,18 @@ class FE_RaviartThomas : public FiniteElement<dim>
                                      */
     virtual bool has_support_on_face (const unsigned int shape_index,
                                      const unsigned int face_index) const;
-
+    
     virtual void interpolate(std::vector<double>&                local_dofs,
                             const std::vector<double>& values) const;
     virtual void interpolate(std::vector<double>&                local_dofs,
                             const std::vector<Vector<double> >& values,
                             unsigned int offset = 0) const;
-    
     virtual void interpolate(
       std::vector<double>& local_dofs,
       const VectorSlice<const std::vector<std::vector<double> > >& values) const;
-                                    /**
-                                     * Determine an estimate for the
-                                     * memory consumption (in bytes)
-                                     * of this object.
-                                     *
-                                     * This function is made virtual,
-                                     * since finite element objects
-                                     * are usually accessed through
-                                     * pointers to their base class,
-                                     * rather than the class itself.
-                                     */
     virtual unsigned int memory_consumption () const;
-
-  protected:    
-    
     virtual FiniteElement<dim> * clone() const;
-  
-    virtual
-    typename Mapping<dim>::InternalDataBase *
-    get_data (const UpdateFlags,
-             const Mapping<dim>& mapping,
-             const Quadrature<dim>& quadrature) const ;
-
-    virtual void
-    fill_fe_values (const Mapping<dim> &mapping,
-                   const typename Triangulation<dim>::cell_iterator &cell,
-                   const Quadrature<dim>                &quadrature,
-                   typename Mapping<dim>::InternalDataBase      &mapping_internal,
-                   typename Mapping<dim>::InternalDataBase      &fe_internal,
-                   FEValuesData<dim>& data) const;
-    
-    virtual void
-    fill_fe_face_values (const Mapping<dim> &mapping,
-                        const typename Triangulation<dim>::cell_iterator &cell,
-                        const unsigned int                    face_no,
-                        const Quadrature<dim-1>                &quadrature,
-                        typename Mapping<dim>::InternalDataBase      &mapping_internal,
-                        typename Mapping<dim>::InternalDataBase      &fe_internal,
-                        FEValuesData<dim>& data) const;
     
-    virtual void
-    fill_fe_subface_values (const Mapping<dim> &mapping,
-                           const typename Triangulation<dim>::cell_iterator &cell,
-                           const unsigned int                    face_no,
-                           const unsigned int                    sub_no,
-                           const Quadrature<dim-1>                &quadrature,
-                           typename Mapping<dim>::InternalDataBase      &mapping_internal,
-                           typename Mapping<dim>::InternalDataBase      &fe_internal,
-                           FEValuesData<dim>& data) const;
-
   private:
                                     /**
                                      * The order of the
@@ -354,45 +222,6 @@ class FE_RaviartThomas : public FiniteElement<dim>
                                      */  
     const unsigned int rt_order;
 
-                                     /**
-                                      * Spaces describing the
-                                      * anisotropic polynomial spaces
-                                      * for each vector component,
-                                      * i.e. there are @p dim
-                                      * elements of this field. The
-                                      * values for this member are
-                                      * created in
-                                      * create_polynomials().
-                                      */
-    const std::vector<AnisotropicPolynomials<dim> > polynomials;
-
-                                     /**
-                                      * For each shape function, store
-                                      * to which vector component (on
-                                      * the unit cell, they are mixed
-                                      * on the real cell by the
-                                      * transformation) they belong,
-                                      * and which index they have
-                                      * within the anisotropic tensor
-                                      * product polynomial space
-                                      * describing this vector
-                                      * component.
-                                      *
-                                      * These values are computed by
-                                      * the compute_renumber()
-                                      * function.
-                                      */
-    const std::vector<std::pair<unsigned int, unsigned int> > renumber;
-    
-    
-                                     /**
-                                      * Generate the polynomial spaces
-                                      * for the polynomials()
-                                      * member.
-                                      */
-    static std::vector<AnisotropicPolynomials<dim> >
-    create_polynomials (const unsigned int degree);
-    
                                     /**
                                      * Only for internal use. Its
                                      * full name is
@@ -416,34 +245,6 @@ class FE_RaviartThomas : public FiniteElement<dim>
     static std::vector<bool>
     get_ria_vector (const unsigned int degree);
 
-                                     /**
-                                      * Compute the values of the
-                                      * @p renumber field.
-                                      */
-    static std::vector<std::pair<unsigned int, unsigned int> >
-    compute_renumber (const unsigned int);
-
-                                    /**
-                                     * Initialize the hanging node
-                                     * constraints matrices. Called
-                                     * from the constructor.
-                                     */
-    void initialize_constraints ();
-
-                                    /**
-                                     * Initialize the embedding
-                                     * matrices. Called from the
-                                     * constructor.
-                                     */
-    void initialize_embedding ();
-
-                                    /**
-                                     * Initialize the restriction
-                                     * matrices. Called from the
-                                     * constructor.
-                                     */
-    void initialize_restriction ();
-    
                                     /**
                                      * Initialize the @p
                                      * generalized_support_points
@@ -456,14 +257,6 @@ class FE_RaviartThomas : public FiniteElement<dim>
                                      */
     void initialize_support_points (const unsigned int rt_degree);
 
-                                    /**
-                                     * Initialize the
-                                     * @p unit_face_support_points field
-                                     * of the FiniteElement
-                                     * class. Called from the
-                                     * constructor.
-                                     */
-    void initialize_face_support_points ();
                                     /**
                                      * Given a set of flags indicating
                                      * what quantities are requested
@@ -658,15 +451,6 @@ class FE_RaviartThomasNodal
 
     virtual FiniteElement<dim>* clone () const;
 
-                                    /**
-                                     * Check whether a shape function
-                                     * may be non-zero on a face.
-                                     *
-                                     * Right now, always returns
-                                     * @p true.
-                                     */
-    virtual bool has_support_on_face (const unsigned int shape_index,
-                                     const unsigned int face_index) const;    
     virtual void interpolate(std::vector<double>&                local_dofs,
                             const std::vector<double>& values) const;
     virtual void interpolate(std::vector<double>&                local_dofs,
@@ -715,70 +499,9 @@ class FE_RaviartThomasNodal
 
 #ifndef DOXYGEN
 
-template <>
-void FE_RaviartThomas<1>::initialize_face_support_points ();
-
 template <>
 std::vector<unsigned int> FE_RaviartThomas<1>::get_dpo_vector (const unsigned int);
 
-template <>
-std::vector<AnisotropicPolynomials<1> >
-FE_RaviartThomas<1>::create_polynomials (const unsigned int);
-
-template <>
-std::vector<AnisotropicPolynomials<2> >
-FE_RaviartThomas<2>::create_polynomials (const unsigned int);
-
-template <>
-std::vector<AnisotropicPolynomials<3> >
-FE_RaviartThomas<3>::create_polynomials (const unsigned int);
-
-template <>
-std::vector<std::pair<unsigned int, unsigned int> >
-FE_RaviartThomas<1>::compute_renumber (const unsigned int);
-
-template <>
-std::vector<std::pair<unsigned int, unsigned int> >
-FE_RaviartThomas<2>::compute_renumber (const unsigned int);
-
-template <>
-std::vector<std::pair<unsigned int, unsigned int> >
-FE_RaviartThomas<3>::compute_renumber (const unsigned int);
-
-template <>
-void
-FE_RaviartThomas<1>::initialize_constraints ();
-
-template <>
-void
-FE_RaviartThomas<2>::initialize_constraints ();
-
-template <>
-void
-FE_RaviartThomas<3>::initialize_constraints ();
-
-template <>
-void
-FE_RaviartThomas<1>::initialize_embedding ();
-
-template <>
-void
-FE_RaviartThomas<1>::initialize_restriction ();
-
-template <>
-void
-FE_RaviartThomas<2>::initialize_restriction ();
-
-template <>
-void
-FE_RaviartThomas<3>::initialize_restriction ();
-
-template <>
-void
-FE_RaviartThomas<1>::
-get_interpolation_matrix (const FiniteElement<1> &,
-                         FullMatrix<double>         &) const;
-
 #endif // DOXYGEN
 
 #endif
index f69681b45a542005d6008339252fcea36afe6cbf..e8c05130429746fcab39f89d4655dd778f8bcb3f 100644 (file)
 #  include <strstream>
 #endif
 
-
-// namespace for some functions that are used in this file. they are
-// specific to numbering conventions used for the FE_RT element, and
-// are thus not very interesting to the outside world
-namespace 
-{
-                                  // auxiliary type to allow for some
-                                  // kind of explicit template
-                                  // specialization of the following
-                                  // functions
-  template <int dim> struct int2type {};
-
-  
-                                  // generate the j-th out of a total
-                                  // of N points on the unit square
-                                  // in 2d. N needs not be a square
-                                  // number, but must be the product
-                                  // of two integers
-                                  //
-                                  // there is one complication: we
-                                  // want to generate interpolation
-                                  // points on the unit square for
-                                  // the shape functions for this
-                                  // element, but for that we need to
-                                  // make sure that these
-                                  // interpolation points make the
-                                  // resulting matrix rows linearly
-                                  // independent. this is a problem
-                                  // since we have anisotropic
-                                  // polynomials, so for example for
-                                  // the lowest order elements, we
-                                  // have as polynomials in for the
-                                  // x-component of the shape
-                                  // functions only "x" and "1-x",
-                                  // i.e. no y-dependence. if we
-                                  // select as interpolation points
-                                  // the points (.5,0) and (.5,1),
-                                  // we're hosed!
-                                  //
-                                  // thus, the third parameter gives
-                                  // the coordinate direction in
-                                  // which the polynomial degree is
-                                  // highest. we use this to select
-                                  // interpolation points primarily
-                                  // in this direction then
-  inline
-  Point<2> generate_unit_point (const unsigned int j,
-                               const unsigned int N,
-                               const unsigned int d,
-                               const int2type<2>  &)
-  {
-    Assert (d<2, ExcInternalError());
-    
-                                    // factorize N int N1*N2. note
-                                    // that we always have N1<=N2,
-                                    // since the square root is
-                                    // rounded down
-    const unsigned int N1 = static_cast<unsigned int>(std::sqrt(1.*N));
-    const unsigned int N2 = N/N1;
-    Assert (N1*N2 == N, ExcInternalError());
-
-    const unsigned int Nx = (d==0 ? N2 : N1),
-                      Ny = (d==1 ? N2 : N1);
-    
-    return Point<2> (Nx == 1 ? .5 : 1.*(j%Nx)/(Nx-1),
-                    Ny == 1 ? .5 : 1.*(j/Nx)/(Ny-1));
-  }
-  
-
-                                  // generate the j-th out of a total
-                                  // of N points on the unit cube
-                                  // in 3d. N needs not be a cube
-                                  // number, but must be the product
-                                  // of three integers
-                                  //
-                                  // the same applies as above for
-                                  // the meaning of the parameter "d"
-  inline
-  Point<3> generate_unit_point (const unsigned int /*j*/,
-                               const unsigned int N,
-                               const unsigned int d,
-                               const int2type<3>  &)
-  {
-    Assert (d<3, ExcInternalError());
-
-    const unsigned int N1 = static_cast<unsigned int>(std::pow(1.*N, 1./3.));
-    const unsigned int N2 = static_cast<unsigned int>(std::sqrt(1.*N/N1));
-    const unsigned int N3 = N/(N1*N2);
-    Assert (N1*N2*N3 == N, ExcInternalError());
-
-    Assert (false, ExcNotImplemented());
-
-    return Point<3> ();
-  }
-  
-}
-
+#include <iostream>
+using namespace std;
 
 
 template <int dim>
-FE_RaviartThomas<dim>::FE_RaviartThomas (const unsigned int rt_order)
+FE_RaviartThomas<dim>::FE_RaviartThomas (const unsigned int deg)
                :
-               FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(rt_order),
-                                                          dim, rt_order+1, FiniteElementData<dim>::Hdiv),
-                                   get_ria_vector (rt_order),
-                                   std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(rt_order),dim,rt_order+1).dofs_per_cell,
-                                                                   std::vector<bool>(dim,true))),
-               rt_order(rt_order),
-                polynomials (create_polynomials(rt_order)),
-                renumber (compute_renumber(rt_order))
+               FE_PolyTensor<PolynomialsRaviartThomas<dim>, dim> (
+                 deg,
+                 FiniteElementData<dim>(get_dpo_vector(deg),
+                                        dim, deg+1, FiniteElementData<dim>::Hdiv),
+                 get_ria_vector (deg),
+                 std::vector<std::vector<bool> >(
+                   FiniteElementData<dim>(get_dpo_vector(deg),
+                                          dim,deg+1).dofs_per_cell,
+                   std::vector<bool>(dim,true))),
+               rt_order(deg)
 {
   Assert (dim >= 2, ExcImpossibleInDim(dim));
-
-                                   // check formula (III.3.22) in the
-                                   // book by Brezzi & Fortin about
-                                   // the number of degrees of freedom
-                                   // per cell
-  Assert (((dim==2) &&
-           (this->dofs_per_cell == 2*(rt_order+1)*(rt_order+2)))
-          ||
-          ((dim==3) &&
-           (this->dofs_per_cell == 3*(rt_order+1)*(rt_order+1)*(rt_order+2))),
-          ExcInternalError());
-  Assert (renumber.size() == this->dofs_per_cell,
-          ExcInternalError());
+  const unsigned int n_dofs = this->dofs_per_cell;
+  
+                                  // First, initialize the
+                                  // generalized support points and
+                                  // quadrature weights, since they
+                                  // are required for interpolation.
+  initialize_support_points(deg);
+                                  // Now compute the inverse node
+                                  //matrix, generating the correct
+                                  //basis functions from the raw
+                                  //ones.
+  FullMatrix<double> M(n_dofs, n_dofs);
+  FETools::compute_node_matrix(M, *this);
+  this->inverse_node_matrix.reinit(n_dofs, n_dofs);
+  this->inverse_node_matrix.invert(M);
+                                  // From now on, the shape functions
+                                  // will be the correct ones, not
+                                  // the raw shape functions anymore.
+  
 
                                   // initialize the various matrices
-  initialize_constraints ();
-
   for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i)
-    this->prolongation[i].reinit (this->dofs_per_cell,
-                                 this->dofs_per_cell);
+    this->prolongation[i].reinit (n_dofs,
+                                 n_dofs);
   FETools::compute_embedding_matrices (*this, &this->prolongation[0]);
   
-  initialize_restriction ();
-
-                                  // finally fill in support points
-                                  // on cell and face
-  initialize_support_points (rt_order);
-  initialize_face_support_points ();
-
+  std::vector<FullMatrix<double> >
+    face_embeddings(1<<(dim-1), FullMatrix<double>(this->dofs_per_face,
+                                                  this->dofs_per_face));
+  FETools::compute_face_embedding_matrices(*this, &face_embeddings[0], 0, 0);
+  this->interface_constraints.reinit((1<<(dim-1)) * this->dofs_per_face,
+                                    this->dofs_per_face);
+  unsigned int target_row=0;
+  for (unsigned int d=0;d<face_embeddings.size();++d)
+    for (unsigned int i=0;i<face_embeddings[d].m();++i)
+      {
+       for (unsigned int j=0;j<face_embeddings[d].n();++j)
+         this->interface_constraints(target_row,j) = face_embeddings[d](i,j);
+       ++target_row;
+      }
+//TODO:[WB] What is this?
                                    // then make
                                    // system_to_component_table
                                    // invalid, since this has no
@@ -217,599 +136,11 @@ FE_RaviartThomas<dim>::clone() const
 }
 
 
-template <int dim>
-double
-FE_RaviartThomas<dim>::shape_value_component (const unsigned int i,
-                                              const Point<dim>    &p,
-                                              const unsigned int component) const
-{
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  Assert (component < dim, ExcIndexRange (component, 0, dim));
-
-                                   // check whether this shape
-                                   // function has a contribution in
-                                   // this component at all, and if so
-                                   // delegate to the respective
-                                   // polynomial
-  if (component == renumber[i].first)
-    return polynomials[component].compute_value(renumber[i].second, p);
-  else
-    return 0;
-}
-
-
-
-template <int dim>
-Tensor<1,dim>
-FE_RaviartThomas<dim>::shape_grad_component (const unsigned int  i,
-                                             const Point<dim>   &p,
-                                             const unsigned int  component) const
-{
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  Assert (component < dim, ExcIndexRange (component, 0, dim));
-
-                                   // check whether this shape
-                                   // function has a contribution in
-                                   // this component at all, and if so
-                                   // delegate to the respective
-                                   // polynomial
-  if (component == renumber[i].first)
-    return polynomials[component].compute_grad(renumber[i].second, p);
-  else
-    return Tensor<1,dim>();
-}
-
-
-
-template <int dim>
-Tensor<2,dim>
-FE_RaviartThomas<dim>::shape_grad_grad_component (const unsigned int i,
-                                                  const Point<dim>  &p,
-                                                  const unsigned int component) const
-{
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  Assert (component < dim, ExcIndexRange (component, 0, dim));
-
-                                   // check whether this shape
-                                   // function has a contribution in
-                                   // this component at all, and if so
-                                   // delegate to the respective
-                                   // polynomial
-  if (component == renumber[i].first)
-    return polynomials[component].compute_grad_grad(renumber[i].second, p);
-  else
-    return Tensor<2,dim>();
-}
-
-
-
-#if deal_II_dimension == 1
-
-template <>
-void
-FE_RaviartThomas<1>::
-get_interpolation_matrix (const FiniteElement<1> &,
-                         FullMatrix<double>         &) const
-{
-  Assert (false, ExcImpossibleInDim(1));
-}
-
-#endif
-
-
-template <int dim>
-void
-FE_RaviartThomas<dim>::
-get_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
-                         FullMatrix<double>           &interpolation_matrix) const
-{
-                                  // this is only implemented, if the
-                                  // source FE is also a
-                                  // Raviart-Thomas element,
-                                  // otherwise throw an exception, as
-                                  // the documentation says
-  AssertThrow ((x_source_fe.get_name().find ("FE_RaviartThomas<") == 0)
-               ||
-               (dynamic_cast<const FE_RaviartThomas<dim>*>(&x_source_fe) != 0),
-               typename FiniteElement<dim>::
-               ExcInterpolationNotImplemented());
-  
-                                  // ok, source is a RT element, so
-                                  // we will be able to do the work
-  const FE_RaviartThomas<dim> &source_fe
-    = dynamic_cast<const FE_RaviartThomas<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));
-  
-  
-                                  // compute the interpolation
-                                  // matrices in much the same way as
-                                  // we do for the embedding matrices
-                                  // from mother to child.
-  const unsigned int dofs_per_coordinate = this->dofs_per_cell/dim;
-  Assert (dofs_per_coordinate*dim == this->dofs_per_cell,
-         ExcInternalError());
-  for (unsigned int d=0; d<dim; ++d)
-    Assert (polynomials[d].n() == dofs_per_coordinate, ExcInternalError());
-
-  const unsigned int source_dofs_per_coordinate = source_fe.dofs_per_cell/dim;
-  Assert (source_dofs_per_coordinate*dim == source_fe.dofs_per_cell,
-         ExcInternalError());
-  for (unsigned int d=0; d<dim; ++d)
-    Assert (source_fe.polynomials[d].n() == source_dofs_per_coordinate, ExcInternalError());
-  
-  FullMatrix<double> cell_interpolation (dofs_per_coordinate,
-                                        dofs_per_coordinate);
-  FullMatrix<double> source_interpolation (dofs_per_coordinate,
-                                          source_dofs_per_coordinate);
-  FullMatrix<double> tmp (dofs_per_coordinate,
-                         source_dofs_per_coordinate);
-  for (unsigned int d=0; d<dim; ++d)
-    {
-      for (unsigned int j=0; j<dofs_per_coordinate; ++j)
-       {
-                                          // generate a point on this
-                                          // cell and evaluate the
-                                          // shape functions there
-                                          //
-                                          // see the comment for that
-                                          // function to see why the
-                                          // third parameter is
-                                          // necessary
-         const Point<dim> p = generate_unit_point (j, dofs_per_coordinate,
-                                                   d, int2type<dim>());
-         for (unsigned int i=0; i<dofs_per_coordinate; ++i)
-           cell_interpolation(j,i) = polynomials[d].compute_value (i, p);
-
-         for (unsigned int i=0; i<source_dofs_per_coordinate; ++i)
-           source_interpolation(j,i) = source_fe.polynomials[d].compute_value (i, p);
-       }
-
-                                      // then compute the
-                                      // interpolation matrix matrix
-                                      // for this coordinate
-                                      // direction
-      cell_interpolation.gauss_jordan ();
-      cell_interpolation.mmult (tmp, source_interpolation);
-
-                                      // finally transfer the
-                                      // results for this
-                                      // coordinate into the matrix
-                                      // corresponding to the
-                                      // entire space on this
-                                      // cell. cut off very small
-                                      // values here
-      for (unsigned int i=0; i<this->dofs_per_cell; ++i)
-       if (renumber[i].first == d)
-         for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
-           if (source_fe.renumber[j].first == d)
-             if (std::fabs(tmp(renumber[i].second,
-                               source_fe.renumber[j].second)) > 1e-15)
-               interpolation_matrix(i,j) = tmp(renumber[i].second,
-                                               source_fe.renumber[j].second);
-    }
-
-                                  // if this were a Lagrange
-                                  // interpolation element, we could
-                                  // make sure that the row sum of
-                                  // each of the matrices is 1 at
-                                  // this point. note that this won't
-                                  // work here, since we are working
-                                  // with hierarchical elements for
-                                  // which the shape functions don't
-                                  // sum up to 1
-                                  //
-                                  // however, we can make sure that
-                                  // only components couple that have
-                                  // the same vector component
-  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
-    for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
-      Assert ((interpolation_matrix(i,j) == 0.) ||
-             (renumber[i].first == source_fe.renumber[j].first),
-             ExcInternalError());
-}
-
-
-
 //---------------------------------------------------------------------------
 // Auxiliary and internal functions
 //---------------------------------------------------------------------------
 
 
-
-
-#if deal_II_dimension == 1
-
-template <>
-void
-FE_RaviartThomas<1>::initialize_constraints ()
-{
-  Assert (false, ExcImpossibleInDim(1));
-}
-
-#endif
-
-#if deal_II_dimension == 2
-
-template <>
-void
-FE_RaviartThomas<2>::initialize_constraints ()
-{
-  const unsigned int dim = 2;
-  
-  this->interface_constraints.
-    TableBase<2,double>::reinit (this->interface_constraints_size());
-
-                                  // this case is too easy, so
-                                  // special case it
-  if (rt_order == 0)
-    {
-      this->interface_constraints(0,0) = this->interface_constraints(1,0) = .5;
-      return;
-    }
-
-                                  // for higher orders of the
-                                  // Raviart-Thomas element:
-    
-                                  // restricted to each face, the
-                                  // normal component of the shape
-                                  // functions is an element of P_{k}
-                                  // (in 2d), or Q_{k} (in 3d), where
-                                  // k is the degree of the element
-                                  //
-                                  // from this, we interpolate
-                                  // between mother and cell
-                                  // face. this is slightly
-                                  // complicated by the fact that we
-                                  // don't use Lagrange interpolation
-                                  // polynomials, but rather
-                                  // hierarchical polynomials, so we
-                                  // can't just use point
-                                  // interpolation. what we do
-                                  // instead is to evaluate at a
-                                  // number of points and then invert
-                                  // the interpolation matrix
-
-                                  // mathematically speaking, this
-                                  // works in the following way: on
-                                  // each subface, we want that
-                                  // finite element solututions from
-                                  // both sides coincide. i.e. if a
-                                  // and b are expansion coefficients
-                                  // for the shape functions from
-                                  // both sides, we seek a relation
-                                  // between x and y such that
-                                  //   sum_i a_i phi^c_i(x)
-                                  //   == sum_j b_j phi_j(x)
-                                  // for all points x on the
-                                  // interface. here, phi^c_i are the
-                                  // shape functions on the small
-                                  // cell on one side of the face,
-                                  // and phi_j those on the big cell
-                                  // on the other side. To get this
-                                  // relation, it suffices to look at
-                                  // a sufficient number of points
-                                  // for which this has to hold. if
-                                  // there are n functions, then we
-                                  // need n evaluation points, and we
-                                  // choose them equidistantly.
-                                  //
-                                  // what one then gets is a matrix
-                                  // system
-                                  //    a A  ==  b B
-                                  // where
-                                  //    A_ij = phi^c_i(x_j)
-                                  //    B_ij = phi_i(x_j)
-                                  // and the relation we are looking for
-                                  // is
-                                  //    a = (A^T)^-1 B^T b
-                                  //
-                                  // below, we build up these
-                                  // matrices, but rather than
-                                  // transposing them after the
-                                  // fact, we do so while building
-                                  // them. A will be
-                                  // subface_interpolation, B will be
-                                  // face_interpolation. note that we
-                                  // build up these matrices for all
-                                  // faces at once, rather than
-                                  // considering them separately. the
-                                  // reason is that we finally will
-                                  // want to have them in this order
-                                  // anyway, as this is the format we
-                                  // need inside deal.II
-  const std::vector<Polynomials::Polynomial<double> >
-    face_polynomials (Polynomials::Hierarchical::
-                     generate_complete_basis (rt_order));
-  Assert (face_polynomials.size() == this->dofs_per_face, ExcInternalError());
-  
-  FullMatrix<double> face_interpolation (2*this->dofs_per_face, this->dofs_per_face);
-  FullMatrix<double> subface_interpolation (2*this->dofs_per_face, 2*this->dofs_per_face);
-  
-                                  // generate the matrix for the
-                                  // evaluation points on the big
-                                  // face, and the corresponding
-                                  // points in the coordinate system
-                                  // of the small face. order the
-                                  // shape functions in the same way
-                                  // we want to have them in the
-                                  // final matrix. extend shape
-                                  // functions on the small faces by
-                                  // zero to the other face on which
-                                  // they are not defined (we do this
-                                  // by simply not considering these
-                                  // entries in the matrix)
-                                  //
-                                  // note the agreeable fact that for
-                                  // this element, all the shape
-                                  // functions we presently care for
-                                  // are face-based (i.e. not vertex
-                                  // shape functions); thus, for this
-                                  // element, we can skip the
-                                  // annoying index shifting for the
-                                  // constraints matrix due to its
-                                  // weird format
-  for (unsigned int subface=0; subface<GeometryInfo<dim>::subfaces_per_face; ++subface)
-    for (unsigned int i=0; i<this->dofs_per_face; ++i)
-      {
-       const double p_face (1.*i/rt_order/2 + (subface == 0 ? 0. : .5));
-       const double p_subface (1.*i/rt_order);
-
-       for (unsigned int j=0; j<this->dofs_per_face; ++j)
-         {
-           face_interpolation(subface*this->dofs_per_face+i,
-                              j)
-             = face_polynomials[j].value(p_face);
-           subface_interpolation(subface*this->dofs_per_face+i,
-                                 subface*this->dofs_per_face+j)
-             = face_polynomials[j].value(p_subface);
-         }
-      }
-      
-  subface_interpolation.gauss_jordan ();
-  subface_interpolation.mmult (this->interface_constraints,
-                              face_interpolation);
-    
-                                  // there is one additional thing to
-                                  // be considered: since the shape
-                                  // functions on the real cell
-                                  // contain the Jacobian (actually,
-                                  // the determinant of the inverse),
-                                  // there is an additional factor of
-                                  // 2 when going from the big to the
-                                  // small cell:
-  this->interface_constraints *= 1./2;
-
-                                  // finally: constraints become
-                                  // really messy if the matrix in
-                                  // question has some entries that
-                                  // are almost zero, but not
-                                  // quite. this will happen in the
-                                  // above procedure due to
-                                  // round-off. let us simply delete
-                                  // these entries
-  for (unsigned int i=0; i<this->interface_constraints.m(); ++i)
-    for (unsigned int j=0; j<this->interface_constraints.n(); ++j)
-      if (std::fabs(this->interface_constraints(i,j)) < 1e-14)
-       this->interface_constraints(i,j) = 0.;
-}
-
-#endif
-
-#if deal_II_dimension == 3
-
-template <>
-void
-FE_RaviartThomas<3>::initialize_constraints ()
-{
-  Assert (false, ExcNotImplemented());
-}
-
-#endif
-
-
-#if deal_II_dimension == 1
-
-template <>
-void
-FE_RaviartThomas<1>::initialize_restriction ()
-{}
-
-#endif
-
-
-#if deal_II_dimension == 2
-
-template <>
-void
-FE_RaviartThomas<2>::initialize_restriction ()
-{
-  const unsigned int dim = 2;
-  switch (rt_order)
-    {
-      case 0:
-      {
-                                        // this is a strange element,
-                                        // since it is both additive
-                                        // and then it is also
-                                        // not. ideally, we would
-                                        // like to have the value of
-                                        // the shape function on the
-                                        // coarse line to be the mean
-                                        // value of that on the two
-                                        // child ones. thus, one
-                                        // should make it
-                                        // additive. however,
-                                        // additivity only works if
-                                        // an element does not have
-                                        // any continuity
-                                        // requirements, since
-                                        // otherwise degrees of
-                                        // freedom are shared between
-                                        // adjacent elements, and
-                                        // when we make the element
-                                        // additive, that would mean
-                                        // that we end up adding up
-                                        // contributions not only
-                                        // from the child cells of
-                                        // this cell, but also from
-                                        // the child cells of the
-                                        // neighbor, and since we
-                                        // cannot know whether there
-                                        // even exists a neighbor we
-                                        // cannot simply make the
-                                        // element additive.
-                                        //
-                                        // so, until someone comes
-                                        // along with a better
-                                        // alternative, we do the
-                                        // following: make the
-                                        // element non-additive, and
-                                        // simply pick the value of
-                                        // one of the child lines for
-                                        // the value of the mother
-                                        // line (note that we have to
-                                        // multiply by two, since the
-                                        // shape functions scale with
-                                        // the inverse Jacobian). we
-                                        // thus throw away the
-                                        // information of one of the
-                                        // child lines, but there
-                                        // seems to be no other way
-                                        // than that...
-                                        //
-                                        // note: to make things
-                                        // consistent, and
-                                        // restriction independent of
-                                        // the order in which we
-                                        // travel across the cells of
-                                        // the coarse grid, we have
-                                        // to make sure that we take
-                                        // the same small line when
-                                        // visiting its two
-                                        // neighbors, to get the
-                                        // value for the mother
-                                        // line. we take the first
-                                        // line always, in the
-                                        // canonical direction of
-                                        // lines
-       for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
-         this->restriction[c].reinit (this->dofs_per_cell,
-                                      this->dofs_per_cell);
-              
-       this->restriction[0](0,0) = 2.;
-       this->restriction[1](1,1) = 2.;
-       this->restriction[3](2,2) = 2.;
-       this->restriction[0](3,3) = 2.;
-
-       break;
-      };
-
-
-      case 1:
-      {
-       for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
-         this->restriction[c].reinit (this->dofs_per_cell,
-                                      this->dofs_per_cell);
-
-                                        // first set the corner
-                                        // nodes. note that they are
-                                        // non-additive
-       this->restriction[0](0,0) = 2.;
-       this->restriction[0](6,6) = 2.;
-
-       this->restriction[1](1,1) = 2.;
-       this->restriction[1](2,2) = 2.;
-
-       this->restriction[2](3,3) = 2.;
-       this->restriction[2](5,5) = 2.;
-
-       this->restriction[3](4,4) = 2.;
-       this->restriction[3](7,7) = 2.;
-
-                                        // then also set the bubble
-                                        // nodes. they _are_
-                                        // additive. to understand
-                                        // what's going on, recall
-                                        // that the bubble shape
-                                        // functions have value -1
-                                        // (!) at the center point,
-                                        // by construction of the
-                                        // polynomials, and that the
-                                        // corner nodes have values
-                                        // 1/2 there since they are
-                                        // just the linears, and not
-                                        // some interpolating
-                                        // polynomial
-                                        //
-                                        // (actually, the
-                                        // additive/non-additive
-                                        // business shouldn't make
-                                        // that much of a difference:
-                                        // node 4 on cell 0 and node
-                                        // 0 on cell 3 must have the
-                                        // same value, since normal
-                                        // components are
-                                        // continuous. so we could
-                                        // pick either and make these
-                                        // shape functions
-                                        // non-additive as well. we
-                                        // choose to take the mean
-                                        // value, which should be the
-                                        // same as either value, and
-                                        // make the shape function
-                                        // additive)
-       this->restriction[0](10,0) = 1.;
-       this->restriction[0](10,4) = -1.;
-       this->restriction[3](10,0) = -1.;
-       this->restriction[3](10,4) = 1.;
-
-       this->restriction[1](11,1) = 1.;
-       this->restriction[1](11,5) = -1.;
-       this->restriction[2](11,1) = -1.;
-       this->restriction[2](11,5) = 1.;
-       
-       this->restriction[0](8,6) = 1.;
-       this->restriction[0](8,2) = -1.;
-       this->restriction[1](8,6) = -1.;
-       this->restriction[1](8,2) = 1.;
-       
-       this->restriction[3](9,7) = 1.;
-       this->restriction[3](9,3) = -1.;
-       this->restriction[2](9,7) = -1.;
-       this->restriction[2](9,3) = 1.;
-       
-       break;
-      };
-       
-                                       // in case we don't have the
-                                       // matrices (yet), leave them
-                                       // empty. this does not
-                                       // prevent the use of this FE,
-                                       // but will prevent the use of
-                                       // these matrices
-    };
-}
-
-#endif
-
-#if deal_II_dimension == 3
-
-template <>
-void
-FE_RaviartThomas<3>::initialize_restriction ()
-{
-  Assert (false, ExcNotImplemented());
-}
-
-#endif
-
-
 #if deal_II_dimension == 1
 
 template <int dim>
@@ -873,8 +204,8 @@ FE_RaviartThomas<dim>::initialize_support_points (const unsigned int deg)
 
       boundary_weights.reinit(n_face_points, legendre.n());
       
-      Assert (face_points.n_quadrature_points == this->dofs_per_face,
-             ExcInternalError());
+//       Assert (face_points.n_quadrature_points == this->dofs_per_face,
+//           ExcInternalError());
       
       for (unsigned int k=0;k<n_face_points;++k)
        {
@@ -889,7 +220,7 @@ FE_RaviartThomas<dim>::initialize_support_points (const unsigned int deg)
                * legendre.compute_value(i, face_points.point(k));
            }
        }
-      
+
       Quadrature<dim> faces = QProjector<dim>::project_to_all_faces(face_points);
       for (;current<GeometryInfo<dim>::faces_per_cell*n_face_points;
           ++current)
@@ -935,34 +266,6 @@ FE_RaviartThomas<dim>::initialize_support_points (const unsigned int deg)
 
 #endif
 
-#if deal_II_dimension == 1
-
-template <>
-void FE_RaviartThomas<1>::initialize_face_support_points ()
-{
-                                  // no faces in 1d, so nothing to do
-}
-
-#endif
-
-
-template <int dim>
-void FE_RaviartThomas<dim>::initialize_face_support_points ()
-{
-  this->unit_face_support_points.resize (this->dofs_per_face);
-
-                                         // like with cell
-                                         // unit_support_points:
-                                         // associate all of the in
-                                         // the face mid-point, since
-                                         // there is no other useful
-                                         // way
-  for (unsigned int i=0; i<this->dofs_per_face; ++i)
-    this->unit_face_support_points[i] = (dim == 2 ?
-                                         Point<dim-1>(.5) :
-                                         Point<dim-1>(.5,.5));
-}
-
 
 #if deal_II_dimension == 1
 
@@ -1055,263 +358,6 @@ FE_RaviartThomas<dim>::get_ria_vector (const unsigned int rt_order)
 }
 
 
-#if deal_II_dimension == 1
-
-template <>
-std::vector<AnisotropicPolynomials<1> >
-FE_RaviartThomas<1>::create_polynomials (const unsigned int)
-{
-  Assert (false, ExcImpossibleInDim(1));
-  return std::vector<AnisotropicPolynomials<1> > ();
-}
-
-#endif
-
-
-#if deal_II_dimension == 2
-
-template <>
-std::vector<AnisotropicPolynomials<2> >
-FE_RaviartThomas<2>::create_polynomials (const unsigned int rt_order)
-{
-  const unsigned int dim = 2;
-  
-                                   // use the fact that the RT(k)
-                                   // spaces are spanned by the
-                                   // functions
-                                   // P_{k+1,k} \times P_{k,k+1},
-                                   // see the book by Brezzi and
-                                   // Fortin
-  const std::vector<Polynomials::Polynomial<double> > pols[2]
-    = { Polynomials::Hierarchical::generate_complete_basis (rt_order+1),
-        Polynomials::Hierarchical::generate_complete_basis (rt_order)};
-
-                                   // create spaces (k+1,k) and (k,k+1)
-  std::vector<std::vector<Polynomials::Polynomial<double> > >
-    pols_vector_1(dim), pols_vector_2(dim);
-  pols_vector_1[0] = pols[0];
-  pols_vector_1[1] = pols[1];
-
-  pols_vector_2[0] = pols[1];
-  pols_vector_2[1] = pols[0];
-  
-  const AnisotropicPolynomials<dim> anisotropic[dim]
-    = { AnisotropicPolynomials<dim> (pols_vector_1),
-        AnisotropicPolynomials<dim> (pols_vector_2) };
-
-                                   // work around a stupid bug in
-                                   // gcc2.95 where the compiler
-                                   // complains about reaching the end
-                                   // of a non-void function when we
-                                   // simply return the following
-                                   // object unnamed, rather than
-                                   // first creating a named object
-                                   // and then returning it...
-  const std::vector<AnisotropicPolynomials<dim> >
-    ret_val (&anisotropic[0], &anisotropic[dim]);
-  return ret_val;
-}
-
-#endif
-
-
-#if deal_II_dimension == 3
-
-template <>
-std::vector<AnisotropicPolynomials<3> >
-FE_RaviartThomas<3>::create_polynomials (const unsigned int rt_order)
-{
-  const unsigned int dim = 3;
-  
-                                   // use the fact that the RT(k)
-                                   // spaces are spanned by the
-                                   // functions
-                                   // P_{k+1,k,k} \times P_{k,k+1,k}
-                                   // \times P_{k,k,k+1},
-                                   // see the book by Brezzi and
-                                   // Fortin
-  const std::vector<Polynomials::Polynomial<double> > pols[2]
-    = { Polynomials::Hierarchical::generate_complete_basis (rt_order+1),
-        Polynomials::Hierarchical::generate_complete_basis (rt_order)};
-
-                                   // create spaces (k+1,k,k),
-                                   // (k,k+1,k) and (k,k,k+1)
-  std::vector<std::vector<Polynomials::Polynomial<double> > >
-    pols_vector_1(dim), pols_vector_2(dim), pols_vector_3(dim);
-  pols_vector_1[0] = pols[0];
-  pols_vector_1[1] = pols[1];
-  pols_vector_1[2] = pols[1];
-
-  pols_vector_2[0] = pols[1];
-  pols_vector_2[1] = pols[0];
-  pols_vector_2[2] = pols[1];
-  
-  pols_vector_3[0] = pols[1];
-  pols_vector_3[1] = pols[1];
-  pols_vector_3[2] = pols[0];
-  
-  const AnisotropicPolynomials<dim> anisotropic[dim]
-    = { AnisotropicPolynomials<dim> (pols_vector_1),
-        AnisotropicPolynomials<dim> (pols_vector_2),
-        AnisotropicPolynomials<dim> (pols_vector_3) };
-
-                                   // work around a stupid bug in
-                                   // gcc2.95 where the compiler
-                                   // complains about reaching the end
-                                   // of a non-void function when we
-                                   // simply return the following
-                                   // object unnamed, rather than
-                                   // first creating a named object
-                                   // and then returning it...
-  const std::vector<AnisotropicPolynomials<dim> >
-    ret_val (&anisotropic[0], &anisotropic[dim]);
-  return ret_val;
-}
-
-#endif
-
-
-
-#if deal_II_dimension == 1
-
-template <>
-std::vector<std::pair<unsigned int, unsigned int> >
-FE_RaviartThomas<1>::compute_renumber (const unsigned int)
-{
-  Assert (false, ExcImpossibleInDim(1));
-  return std::vector<std::pair<unsigned int, unsigned int> > ();
-}
-
-#endif
-
-
-#if deal_II_dimension == 2
-
-template <>
-std::vector<std::pair<unsigned int, unsigned int> >
-FE_RaviartThomas<2>::compute_renumber (const unsigned int rt_order)
-{
-  const unsigned int dim = 2;
-  
-  std::vector<std::pair<unsigned int, unsigned int> > ret_val;
-  
-                                   // to explain the following: the
-                                   // first (rt_order+1) shape functions
-                                   // are on face 0, and point in
-                                   // y-direction, so are for the
-                                   // second vector component. then
-                                   // there are (rt_order+1) shape
-                                   // functions on face 1, which is
-                                   // for the x vector component, and
-                                   // so on. since the order of face
-                                   // rt_orders of freedom is arbitrary,
-                                   // we simply use the same order as
-                                   // that provided by the 1d
-                                   // polynomial class on which this
-                                   // element is based. after
-                                   // 4*(rt_order+1), the remaining
-                                   // shape functions are all bubbles,
-                                   // so we can number them in any way
-                                   // we want. we do so by first
-                                   // numbering the x-vectors, then
-                                   // the y-vectors
-                                   //
-                                   // now, we have to find a mapping
-                                   // from the above ordering to:
-                                   // first which vector component
-                                   // they belong to (easy), and
-                                   // second the index within this
-                                   // component as provided by the
-                                   // AnisotropicPolynomials class
-                                   //
-                                   // this is mostly a counting
-                                   // argument, tedious and error
-                                   // prone, and so boring to explain
-                                   // that we rather not try to do so
-                                   // here (it's simple, but boring,
-                                   // as said), aside from a few
-                                   // comments below
-
-                                   // face 0
-  for (unsigned int i=0; i<rt_order+1; ++i)
-    ret_val.push_back (std::make_pair (1U, i));
-  
-                                   // face 1
-  for (unsigned int i=0; i<rt_order+1; ++i)
-    ret_val.push_back (std::make_pair (0U, (rt_order+2)*i+1));
-  
-                                   // face 2
-  for (unsigned int i=0; i<rt_order+1; ++i)
-    ret_val.push_back (std::make_pair (1U, (rt_order+1)+i));
-  
-                                   // face 3
-  for (unsigned int i=0; i<rt_order+1; ++i)
-    ret_val.push_back (std::make_pair (0U, (rt_order+2)*i));
-
-                                   // then go on with interior bubble
-                                   // functions, first for the
-                                   // x-direction, then for the
-                                   // y-direction
-  for (unsigned int x=0; x<rt_order; ++x)
-    for (unsigned int y=0; y<rt_order+1; ++y)
-      {
-        const unsigned int index_in_component = (x+2) + y*(rt_order+2);
-        Assert (index_in_component < (rt_order+1)*(rt_order+2),
-                ExcInternalError());
-        ret_val.push_back (std::make_pair(0U, index_in_component));
-      }
-  for (unsigned int x=0; x<rt_order+1; ++x)
-    for (unsigned int y=0; y<rt_order; ++y)
-      {
-        const unsigned int index_in_component = 2*(rt_order+1) + y + x*rt_order;
-        Assert (index_in_component < (rt_order+1)*(rt_order+2),
-                ExcInternalError());
-        ret_val.push_back (std::make_pair(1U, index_in_component));
-      }
-
-#ifdef DEBUG  
-                                   // make sure we have actually used
-                                   // up all elements of the tensor
-                                   // product polynomial
-  Assert (ret_val.size() == 2*(rt_order+1)*(rt_order+2),
-          ExcInternalError());
-  std::vector<bool> test[dim] = { std::vector<bool>(ret_val.size()/dim, false),
-                                  std::vector<bool>(ret_val.size()/dim, false) };
-  for (unsigned int i=0; i<ret_val.size(); ++i)
-    {
-      Assert (ret_val[i].first < dim, ExcInternalError());
-      Assert (ret_val[i].second < test[ret_val[i].first].size(),
-             ExcInternalError());
-      Assert (test[ret_val[i].first][ret_val[i].second] == false,
-              ExcInternalError());
-      
-      test[ret_val[i].first][ret_val[i].second] = true;
-    }
-  for (unsigned int d=0; d<dim; ++d)
-    Assert (std::find (test[d].begin(), test[d].end(), false) == test[d].end(),
-            ExcInternalError());
-#endif
-  
-  return ret_val;
-}
-
-#endif
-
-
-#if deal_II_dimension == 3
-
-template <>
-std::vector<std::pair<unsigned int, unsigned int> >
-FE_RaviartThomas<3>::compute_renumber (const unsigned int /*rt_order*/)
-{
-  Assert (false, ExcNotImplemented());
-  return std::vector<std::pair<unsigned int, unsigned int> > ();  
-}
-
-#endif
-
-
-
 
 template <int dim>
 UpdateFlags
@@ -1341,489 +387,10 @@ FE_RaviartThomas<dim>::update_each (const UpdateFlags flags) const
   return out;
 }
 
-
-
 //---------------------------------------------------------------------------
 // Data field initialization
 //---------------------------------------------------------------------------
 
-template <int dim>
-typename Mapping<dim>::InternalDataBase *
-FE_RaviartThomas<dim>::get_data (const UpdateFlags      update_flags,
-                                 const Mapping<dim>    &mapping,
-                                 const Quadrature<dim> &quadrature) const
-{
-                                  // generate a new data object and
-                                  // initialize some fields
-  InternalData* data = new InternalData;
-
-                                  // check what needs to be
-                                  // initialized only once and what
-                                  // on every cell/face/subface we
-                                  // visit
-  data->update_once = update_once(update_flags);
-  data->update_each = update_each(update_flags);
-  data->update_flags = data->update_once | data->update_each;
-
-  const UpdateFlags flags(data->update_flags);
-  const unsigned int n_q_points = quadrature.n_quadrature_points;
-
-                                  // initialize fields only if really
-                                  // necessary. otherwise, don't
-                                  // allocate memory
-  if (flags & update_values)
-    data->shape_values.resize (this->dofs_per_cell,
-                               std::vector<Tensor<1,dim> > (n_q_points));
-
-  if (flags & update_gradients)
-    data->shape_gradients.resize (this->dofs_per_cell,
-                                  std::vector<Tensor<2,dim> > (n_q_points));
-
-                                  // if second derivatives through
-                                  // finite differencing is required,
-                                  // then initialize some objects for
-                                  // that
-  if (flags & update_second_derivatives)
-    data->initialize_2nd (this, mapping, quadrature);
-
-                                  // next already fill those fields
-                                  // of which we have information by
-                                  // now. note that the shape values
-                                  // and gradients are only those on
-                                  // the unit cell, and need to be
-                                  // transformed when visiting an
-                                  // actual cell
-  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
-    for (unsigned int q=0; q<n_q_points; ++q)
-      {
-        if (flags & update_values)
-          for (unsigned int c=0; c<dim; ++c)
-            data->shape_values[i][q][c]
-              = shape_value_component(i,quadrature.point(q),c);
-       
-        if (flags & update_gradients)
-          for (unsigned int c=0; c<dim; ++c)
-            data->shape_gradients[i][q][c]
-              = shape_grad_component(i,quadrature.point(q),c);
-      }
-   
-  return data;
-}
-
-
-
-
-//---------------------------------------------------------------------------
-// Fill data of FEValues
-//---------------------------------------------------------------------------
-
-template <int dim>
-void
-FE_RaviartThomas<dim>::fill_fe_values (const Mapping<dim>                   &mapping,
-                                       const typename Triangulation<dim>::cell_iterator &cell,
-                                       const Quadrature<dim>                &quadrature,
-                                       typename Mapping<dim>::InternalDataBase &mapping_data,
-                                       typename Mapping<dim>::InternalDataBase &fedata,
-                                       FEValuesData<dim>                    &data) const
-{
-                                  // convert data object to internal
-                                  // data for this class. fails with
-                                  // an exception if that is not
-                                  // possible
-  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-
-                                  // get the flags indicating the
-                                  // fields that have to be filled
-  const UpdateFlags flags(fe_data.current_update_flags());
-
-  const unsigned int n_q_points = quadrature.n_quadrature_points;
-                                 
-                                  // fill shape function
-                                  // values. these are vector-valued,
-                                  // so we have to transform
-                                  // them. since the output format
-                                  // (in data.shape_values) is a
-                                  // sequence of doubles (one for
-                                  // each non-zero shape function
-                                  // value, and for each quadrature
-                                  // point, rather than a sequence of
-                                  // small vectors, we have to use a
-                                  // number of conversions
-  if (flags & update_values)
-    {
-      std::vector<Tensor<1,dim> > shape_values (n_q_points);
-
-      Assert (data.shape_values.n_rows() == this->dofs_per_cell * dim,
-             ExcInternalError());
-      Assert (data.shape_values.n_cols() == n_q_points,
-             ExcInternalError());
-      
-      for (unsigned int k=0; k<this->dofs_per_cell; ++k)
-       {
-                                          // first transform shape
-                                          // values...
-         Assert (fe_data.shape_values[k].size() == n_q_points,
-                 ExcInternalError());
-         mapping.transform_covariant(fe_data.shape_values[k], 0,
-                                      shape_values,
-                                      mapping_data);
-
-                                          // then copy over to target:
-         for (unsigned int q=0; q<n_q_points; ++q)
-           for (unsigned int d=0; d<dim; ++d)
-             data.shape_values[k*dim+d][q] = shape_values[q][d];
-       };
-    };
-  
-      
-  if (flags & update_gradients)
-    {
-      std::vector<Tensor<2,dim> > shape_grads1 (n_q_points);
-      std::vector<Tensor<2,dim> > shape_grads2 (n_q_points);
-
-      Assert (data.shape_gradients.size() == this->dofs_per_cell * dim,
-             ExcInternalError());
-      Assert (data.shape_gradients[0].size() == n_q_points,
-             ExcInternalError());
-
-                                       // loop over all shape
-                                       // functions, and treat the
-                                       // gradients of each shape
-                                       // function at all quadrature
-                                       // points
-      for (unsigned int k=0; k<this->dofs_per_cell; ++k)
-       {
-                                           // treat the gradients of
-                                           // this particular shape
-                                           // function at all
-                                           // q-points. if Dv is the
-                                           // gradient of the shape
-                                           // function on the unit
-                                           // cell, then
-                                           // (J^-T)Dv(J^-1) is the
-                                           // value we want to have on
-                                           // the real cell. so, we
-                                           // will have to apply a
-                                           // covariant transformation
-                                           // to Dv twice. since the
-                                           // interface only allows
-                                           // multiplication with
-                                           // (J^-1) from the right,
-                                           // we have to trick a
-                                           // little in between
-         Assert (fe_data.shape_gradients[k].size() == n_q_points,
-                 ExcInternalError());
-                                           // do first transformation
-         mapping.transform_covariant(fe_data.shape_gradients[k], 0,
-                                      shape_grads1,
-                                      mapping_data);
-                                           // transpose matrix
-          for (unsigned int q=0; q<n_q_points; ++q)
-            shape_grads2[q] = transpose(shape_grads1[q]);
-                                           // do second transformation
-         mapping.transform_covariant(shape_grads2, 0, shape_grads1,
-                                      mapping_data);
-                                           // transpose back
-          for (unsigned int q=0; q<n_q_points; ++q)
-            shape_grads2[q] = transpose(shape_grads1[q]);
-          
-                                          // then copy over to target:
-         for (unsigned int q=0; q<n_q_points; ++q)
-           for (unsigned int d=0; d<dim; ++d)
-             data.shape_gradients[k*dim+d][q] = shape_grads2[q][d];
-       };
-    }
-
-  if (flags & update_second_derivatives)
-    this->compute_2nd (mapping, cell,
-                       QProjector<dim>::DataSetDescriptor::cell(),
-                       mapping_data, fe_data, data);
-}
-
-
-
-template <int dim>
-void
-FE_RaviartThomas<dim>::fill_fe_face_values (const Mapping<dim>                   &mapping,
-                                            const typename Triangulation<dim>::cell_iterator &cell,
-                                            const unsigned int                    face,
-                                            const Quadrature<dim-1>              &quadrature,
-                                            typename Mapping<dim>::InternalDataBase       &mapping_data,
-                                            typename Mapping<dim>::InternalDataBase       &fedata,
-                                            FEValuesData<dim>                    &data) const
-{
-                                  // convert data object to internal
-                                  // data for this class. fails with
-                                  // an exception if that is not
-                                  // possible
-  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-
-                                   // offset determines which data set
-                                  // to take (all data sets for all
-                                  // faces are stored contiguously)
-  const typename QProjector<dim>::DataSetDescriptor offset
-    = (QProjector<dim>::DataSetDescriptor::
-       face (face, cell->face_orientation(face),
-             quadrature.n_quadrature_points));
-
-                                  // get the flags indicating the
-                                  // fields that have to be filled
-  const UpdateFlags flags(fe_data.current_update_flags());
-
-  const unsigned int n_q_points = quadrature.n_quadrature_points;
-                                 
-                                  // fill shape function
-                                  // values. these are vector-valued,
-                                  // so we have to transform
-                                  // them. since the output format
-                                  // (in data.shape_values) is a
-                                  // sequence of doubles (one for
-                                  // each non-zero shape function
-                                  // value, and for each quadrature
-                                  // point, rather than a sequence of
-                                  // small vectors, we have to use a
-                                  // number of conversions
-  if (flags & update_values)
-    {
-      Assert (fe_data.shape_values.size() == this->dofs_per_cell,
-              ExcInternalError());
-      Assert (fe_data.shape_values[0].size() ==
-              GeometryInfo<dim>::faces_per_cell * n_q_points *
-             (dim == 3 ? 2 : 1),
-              ExcInternalError());
-      
-      std::vector<Tensor<1,dim> > shape_values (n_q_points);
-
-      Assert (data.shape_values.n_rows() == this->dofs_per_cell * dim,
-             ExcInternalError());
-      Assert (data.shape_values.n_cols() == n_q_points,
-             ExcInternalError());
-      
-      for (unsigned int k=0; k<this->dofs_per_cell; ++k)
-       {
-                                          // first transform shape
-                                          // values...
-         mapping.transform_covariant(fe_data.shape_values[k], offset,
-                                      shape_values,
-                                      mapping_data);
-
-                                          // then copy over to target:
-         for (unsigned int q=0; q<n_q_points; ++q)
-           for (unsigned int d=0; d<dim; ++d)
-             data.shape_values[k*dim+d][q] = shape_values[q][d];
-       };
-    };
-  
-      
-  if (flags & update_gradients)
-    {
-      Assert (fe_data.shape_values.size() == this->dofs_per_cell,
-              ExcInternalError());
-      Assert (fe_data.shape_gradients[0].size() ==
-              GeometryInfo<dim>::faces_per_cell * n_q_points *
-             (dim == 3 ? 2 : 1),
-              ExcInternalError());
-
-      std::vector<Tensor<2,dim> > shape_grads1 (n_q_points);
-      std::vector<Tensor<2,dim> > shape_grads2 (n_q_points);
-
-      Assert (data.shape_gradients.size() == this->dofs_per_cell * dim,
-             ExcInternalError());
-      Assert (data.shape_gradients[0].size() == n_q_points,
-             ExcInternalError());
-
-                                       // loop over all shape
-                                       // functions, and treat the
-                                       // gradients of each shape
-                                       // function at all quadrature
-                                       // points
-      for (unsigned int k=0; k<this->dofs_per_cell; ++k)
-       {
-                                           // treat the gradients of
-                                           // this particular shape
-                                           // function at all
-                                           // q-points. if Dv is the
-                                           // gradient of the shape
-                                           // function on the unit
-                                           // cell, then
-                                           // (J^-T)Dv(J^-1) is the
-                                           // value we want to have on
-                                           // the real cell. so, we
-                                           // will have to apply a
-                                           // covariant transformation
-                                           // to Dv twice. since the
-                                           // interface only allows
-                                           // multiplication with
-                                           // (J^-1) from the right,
-                                           // we have to trick a
-                                           // little in between
-                                           // 
-                                           // do first transformation
-         mapping.transform_covariant(fe_data.shape_gradients[k], offset,
-                                      shape_grads1,
-                                      mapping_data);
-                                           // transpose matrix
-          for (unsigned int q=0; q<n_q_points; ++q)
-            shape_grads2[q] = transpose(shape_grads1[q]);
-                                           // do second transformation
-         mapping.transform_covariant(shape_grads2, 0, shape_grads1,
-                                      mapping_data);
-                                           // transpose back
-          for (unsigned int q=0; q<n_q_points; ++q)
-            shape_grads2[q] = transpose(shape_grads1[q]);
-          
-                                          // then copy over to target:
-         for (unsigned int q=0; q<n_q_points; ++q)
-           for (unsigned int d=0; d<dim; ++d)
-             data.shape_gradients[k*dim+d][q] = shape_grads2[q][d];
-       };
-    }
-
-  if (flags & update_second_derivatives)
-    this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
-}
-
-
-
-template <int dim>
-void
-FE_RaviartThomas<dim>::fill_fe_subface_values (const Mapping<dim>                   &mapping,
-                                               const typename Triangulation<dim>::cell_iterator &cell,
-                                               const unsigned int                    face,
-                                               const unsigned int                    subface,
-                                               const Quadrature<dim-1>              &quadrature,
-                                               typename Mapping<dim>::InternalDataBase       &mapping_data,
-                                               typename Mapping<dim>::InternalDataBase       &fedata,
-                                               FEValuesData<dim>                    &data) const
-{
-                                  // convert data object to internal
-                                  // data for this class. fails with
-                                  // an exception if that is not
-                                  // possible
-  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-
-                                   // offset determines which data set
-                                  // to take (all data sets for all
-                                  // faces are stored contiguously)
-  const typename QProjector<dim>::DataSetDescriptor offset
-    = (QProjector<dim>::DataSetDescriptor::
-       sub_face (face, subface, cell->face_orientation(face),
-                 quadrature.n_quadrature_points));
-
-                                  // get the flags indicating the
-                                  // fields that have to be filled
-  const UpdateFlags flags(fe_data.current_update_flags());
-
-  const unsigned int n_q_points = quadrature.n_quadrature_points;
-                                 
-                                  // fill shape function
-                                  // values. these are vector-valued,
-                                  // so we have to transform
-                                  // them. since the output format
-                                  // (in data.shape_values) is a
-                                  // sequence of doubles (one for
-                                  // each non-zero shape function
-                                  // value, and for each quadrature
-                                  // point, rather than a sequence of
-                                  // small vectors, we have to use a
-                                  // number of conversions
-  if (flags & update_values)
-    {
-      Assert (fe_data.shape_values[0].size() ==
-              GeometryInfo<dim>::faces_per_cell *
-             GeometryInfo<dim>::subfaces_per_face *
-             n_q_points,
-              ExcInternalError());
-      
-      std::vector<Tensor<1,dim> > shape_values (n_q_points);
-
-      Assert (data.shape_values.n_rows() == this->dofs_per_cell * dim,
-             ExcInternalError());
-      Assert (data.shape_values.n_cols() == n_q_points,
-             ExcInternalError());
-      
-      for (unsigned int k=0; k<this->dofs_per_cell; ++k)
-       {
-                                          // first transform shape
-                                          // values...
-         mapping.transform_covariant(fe_data.shape_values[k], offset,
-                                      shape_values,
-                                      mapping_data);
-
-                                          // then copy over to target:
-         for (unsigned int q=0; q<n_q_points; ++q)
-           for (unsigned int d=0; d<dim; ++d)
-             data.shape_values[k*dim+d][q] = shape_values[q][d];
-       };
-    };
-  
-      
-  if (flags & update_gradients)
-    {
-      Assert (fe_data.shape_gradients.size() ==
-              GeometryInfo<dim>::faces_per_cell *
-             GeometryInfo<dim>::subfaces_per_face *
-             n_q_points,
-              ExcInternalError());
-
-      std::vector<Tensor<2,dim> > shape_grads1 (n_q_points);
-      std::vector<Tensor<2,dim> > shape_grads2 (n_q_points);
-
-      Assert (data.shape_gradients.size() == this->dofs_per_cell * dim,
-             ExcInternalError());
-      Assert (data.shape_gradients[0].size() == n_q_points,
-             ExcInternalError());
-
-                                       // loop over all shape
-                                       // functions, and treat the
-                                       // gradients of each shape
-                                       // function at all quadrature
-                                       // points
-      for (unsigned int k=0; k<this->dofs_per_cell; ++k)
-       {
-                                           // treat the gradients of
-                                           // this particular shape
-                                           // function at all
-                                           // q-points. if Dv is the
-                                           // gradient of the shape
-                                           // function on the unit
-                                           // cell, then
-                                           // (J^-T)Dv(J^-1) is the
-                                           // value we want to have on
-                                           // the real cell. so, we
-                                           // will have to apply a
-                                           // covariant transformation
-                                           // to Dv twice. since the
-                                           // interface only allows
-                                           // multiplication with
-                                           // (J^-1) from the right,
-                                           // we have to trick a
-                                           // little in between
-                                           // 
-                                           // do first transformation
-         mapping.transform_covariant(fe_data.shape_gradients[k], offset,
-                                      shape_grads1,
-                                      mapping_data);
-                                           // transpose matrix
-          for (unsigned int q=0; q<n_q_points; ++q)
-            shape_grads2[q] = transpose(shape_grads1[q]);
-                                           // do second transformation
-         mapping.transform_covariant(shape_grads2, 0, shape_grads1,
-                                      mapping_data);
-                                           // transpose back
-          for (unsigned int q=0; q<n_q_points; ++q)
-            shape_grads2[q] = transpose(shape_grads1[q]);
-          
-                                          // then copy over to target:
-         for (unsigned int q=0; q<n_q_points; ++q)
-           for (unsigned int d=0; d<dim; ++d)
-             data.shape_gradients[k*dim+d][q] = shape_grads2[q][d];
-       };
-    }
-
-  if (flags & update_second_derivatives)
-    this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
-}
 
 
 
@@ -1935,9 +502,17 @@ FE_RaviartThomas<dim>::interpolate(
     for (unsigned int k=0;k<n_face_points;++k)
       for (unsigned int i=0;i<boundary_weights.size(1);++i)
       {
-       local_dofs[i] += boundary_weights(k,i)
-                        * values[face*n_face_points+k](offset+GeometryInfo<dim>::unit_normal_direction[face]);
+       local_dofs[i+face*this->dofs_per_face] += boundary_weights(k,i)
+                        * values[face*n_face_points+k](GeometryInfo<dim>::unit_normal_direction[face]+offset);
       }
+  
+  const unsigned start_cell_dofs = GeometryInfo<dim>::faces_per_cell*this->dofs_per_face;
+  const unsigned start_cell_points = GeometryInfo<dim>::faces_per_cell*n_face_points;
+  
+  for (unsigned int k=0;k<interior_weights.size(0);++k)
+    for (unsigned int i=0;i<interior_weights.size(1);++i)
+      for (unsigned int d=0;d<dim;++d)
+       local_dofs[start_cell_dofs+i*dim+d] += interior_weights(k,i,d) * values[k+start_cell_points](d+offset);
 }
 
 
@@ -1961,9 +536,17 @@ FE_RaviartThomas<dim>::interpolate(
     for (unsigned int k=0;k<n_face_points;++k)
       for (unsigned int i=0;i<boundary_weights.size(1);++i)
       {
-       local_dofs[i] += boundary_weights(k,i)
+       local_dofs[i+face*this->dofs_per_face] += boundary_weights(k,i)
                         * values[GeometryInfo<dim>::unit_normal_direction[face]][face*n_face_points+k];
       }
+  
+  const unsigned start_cell_dofs = GeometryInfo<dim>::faces_per_cell*this->dofs_per_face;
+  const unsigned start_cell_points = GeometryInfo<dim>::faces_per_cell*n_face_points;
+  
+  for (unsigned int k=0;k<interior_weights.size(0);++k)
+    for (unsigned int i=0;i<interior_weights.size(1);++i)
+      for (unsigned int d=0;d<dim;++d)
+       local_dofs[start_cell_dofs+i*dim+d] += interior_weights(k,i,d) * values[d][k+start_cell_points];
 }
 
 
@@ -1976,14 +559,4 @@ FE_RaviartThomas<dim>::memory_consumption () const
 }
 
 
-
-template <int dim>
-unsigned int
-FE_RaviartThomas<dim>::get_degree () const
-{
-  return rt_order;
-}
-
-
-
 template class FE_RaviartThomas<deal_II_dimension>;
index 3c978fa9d21f48f158e0651dc7cbcffdcee36367..a78698922286d88718b6ade443065fae4eb414b3 100644 (file)
@@ -29,8 +29,6 @@
 #  include <strstream>
 #endif
 
-#include <iostream>
-
 template <int dim>
 FE_RaviartThomasNodal<dim>::FE_RaviartThomasNodal (const unsigned int deg)
                :
@@ -55,7 +53,6 @@ FE_RaviartThomasNodal<dim>::FE_RaviartThomasNodal (const unsigned int deg)
                                   // Set up the generalized support
                                   // points
   initialize_support_points (deg);
-
                                   //Now compute the inverse node
                                   //matrix, generating the correct
                                   //basis functions from the raw
@@ -123,15 +120,15 @@ FE_RaviartThomasNodal<dim>::get_name () const
 }
 
 
-
 template <int dim>
-bool
-FE_RaviartThomasNodal<dim>::has_support_on_face (unsigned int, unsigned int) const
+FiniteElement<dim> *
+FE_RaviartThomasNodal<dim>::clone() const
 {
-  return true;
+  return new FE_RaviartThomasNodal<dim>(this->degree-1);
 }
 
 
+
 template <int dim>
 void
 FE_RaviartThomasNodal<dim>::interpolate(
@@ -238,14 +235,6 @@ FE_RaviartThomasNodal<dim>::interpolate(
 
 
 
-template <int dim>
-FiniteElement<dim> *
-FE_RaviartThomasNodal<dim>::clone() const
-{
-  return new FE_RaviartThomasNodal<dim>(this->degree-1);
-}
-
-
 #if deal_II_dimension == 1
 
 template <>

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.