From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 7 Apr 2003 15:47:33 +0000 (+0000)
Subject: Remove again from main branch until element is actually tested a little.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fea705afd4de89870c6e7e9abefec463693be1a6;p=dealii-svn.git

Remove again from main branch until element is actually tested a little.


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

diff --git a/deal.II/deal.II/include/fe/fe_raviart_thomas.h b/deal.II/deal.II/include/fe/fe_raviart_thomas.h
deleted file mode 100644
index 73ee24cc7a..0000000000
--- a/deal.II/deal.II/include/fe/fe_raviart_thomas.h
+++ /dev/null
@@ -1,581 +0,0 @@
-//---------------------------------------------------------------
-//    $Id$
-//    Version: $Name$
-//
-//    Copyright (C) 2002, 2003 by the deal.II authors
-//
-//    This file is subject to QPL and may not be  distributed
-//    without copyright and license information. Please refer
-//    to the file deal.II/doc/license.html for the  text  and
-//    further information on this license.
-//
-//---------------------------------------------------------------
-#ifndef __deal2__fe_raviart_thomas_h
-#define __deal2__fe_raviart_thomas_h
-
-#include <base/config.h>
-#include <base/polynomial.h>
-#include <base/tensor_product_polynomials.h>
-#include <grid/geometry_info.h>
-#include <fe/fe.h>
-
-template <int dim> class TensorProductPolynomials;
-template <int dim> class MappingQ;
-
-
-
-/**
- * Implementation of continuous Raviart-Thomas elements for the space
- * H_div. Note, however, that continuity only concerns the normal
- * component of the vector field.
- *
- * The constructor of this class takes the degree @p{p} of this finite
- * element. However, presently, only lowest order elements
- * (i.e. @p{p==1}) are implemented.
- * 
- * 
- * @sect3{Interpolation to finer and coarser meshes}
- *
- * Each finite element class in deal.II provides matrices that are
- * used to interpolate from coarser to finer meshes and the other way
- * round. Interpolation from a mother cell to its children is usually
- * trivial, since finite element spaces are normally nested and this
- * kind of interpolation is therefore exact. On the other hand, when
- * we interpolate from child cells to the mother cell, we usually have
- * to throw away some information.
- *
- * For continuous elements, this transfer usually happens by
- * interpolating the values on the child cells at the support points
- * of the shape functions of the mother cell. However, for
- * discontinuous elements, we often use a projection from the child
- * cells to the mother cell. The projection approach is only possible
- * for discontinuous elements, since it cannot be guaranteed that the
- * values of the projected functions on one cell and its neighbor
- * match. In this case, only an interpolation can be
- * used. (Internally, whether the values of a shape function are
- * interpolated or projected, or better: whether the matrices the
- * finite element provides are to be treated with the properties of a
- * projection or of an interpolation, is controlled by the
- * @p{restriction_is_additive} flag. See there for more information.)
- *
- * Here, things are not so simple: since the element has some
- * continuity requirements across faces, we can only resort to some
- * kind of interpolation. On the other hand, for the lowest order
- * elements, the values of generating functionals are the (constant)
- * tangential values of the shape functions. We would therefore really
- * like to take the mean value of the tangential values of the child
- * faces, and make this the value of the mother face. Then, however,
- * taking a mean value of two piecewise constant function is not an
- * interpolation, but a restriction. Since this is not possible, we
- * cannot use this.
- *
- * To make a long story somewhat shorter, when interpolating from
- * refined edges to a coarse one, we do not take the mean value, but
- * pick only one (the one from the first child edge). While this is
- * not optimal, it is certainly a valid choice (using an interpolation
- * point that is not in the middle of the cell, but shifted to one
- * side), and it also preserves the order of the interpolation.
- * 
- *
- * @sect3{Numbering of the degrees of freedom (DoFs)}
- *
- * Nedelec elements have their degrees of freedom on edges, with shape
- * functions being vector valued and pointing in tangential
- * direction. We use the standard enumeration and direction of edges
- * in deal.II, yielding the following shape functions in 2d:
- *
- *   @begin{verbatim}
- *          2
- *      *--->---*
- *      |       |
- *     3^       ^1
- *      |       |
- *      *--->---*
- *          0
- *   @end{verbatim}
- *
- * For the 3d case, the ordering follows the same scheme: the lines
- * are numbered as described in the documentation of the
- * @ref{Triangulation} class, i.e.
- *   @begin{verbatim}
- *         *---6---*        *---6---*
- *        /|       |       /       /|
- *      11 |       5      11     10 5
- *      /  7       |     /       /  |
- *     *   |       |    *---2---*   |
- *     |   *---4---*    |       |   *
- *     |  /       /     |       1  /
- *     3 8       9      3       | 9
- *     |/       /       |       |/
- *     *---0---*        *---0---*
- *   @end{verbatim}
- * and their directions are as follows:
- *   @begin{verbatim}
- *         *--->---*        *--->---*
- *        /|       |       /       /|
- *       ^ |       ^      ^       ^ ^
- *      /  ^       |     /       /  |
- *     *   |       |    *--->---*   |
- *     |   *--->---*    |       |   *
- *     |  /       /     |       ^  /
- *     ^ ^       ^      ^       | ^
- *     |/       /       |       |/
- *     *--->---*        *--->---*
- *   @end{verbatim}
- *
- * The element does not make much sense in 1d, so it is not
- * implemented there.
- *
- *
- * @author Wolfgang Bangerth, 2003
- */
-template <int dim>
-class FE_RaviartThomas : public FiniteElement<dim>
-{
-  public:
-				     /**
-				      * Constructor for the Nedelec
-				      * element of degree @p{p}.
-				      */
-    FE_RaviartThomas (const unsigned int p);
-    
-				     /**
-				      * Return the value of the
-				      * @p{component}th vector
-				      * component of the @p{i}th shape
-				      * function at the point
-				      * @p{p}. See the
-				      * @ref{FiniteElementBase} 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{component}th vector
-				      * component of the @p{i}th shape
-				      * function at the point
-				      * @p{p}. See the
-				      * @ref{FiniteElementBase} 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{component}th vector
-				      * component of the @p{i}th shape
-				      * function at the point
-				      * @p{p}. See the
-				      * @ref{FiniteElementBase} 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 polynomial degree
-				      * of this finite element,
-				      * i.e. the value passed to the
-				      * constructor.
-				      */
-    unsigned int get_degree () const;
-    
-				     /**
-				      * Number of base elements in a
-				      * mixed discretization. Here,
-				      * this is of course equal to
-				      * one.
-				      */
-    virtual unsigned int n_base_elements () const;
-    
-				     /**
-				      * Access to base element
-				      * objects. Since this element is
-				      * atomic, @p{base_element(0)} is
-				      * @p{this}, and all other
-				      * indices throw an error.
-				      */
-    virtual const FiniteElement<dim> &
-    base_element (const unsigned int index) const;
-
-                                     /**
-                                      * Multiplicity of base element
-                                      * @p{index}. Since this is an
-                                      * atomic element,
-                                      * @p{element_multiplicity(0)}
-                                      * returns one, and all other
-                                      * indices will throw an error.
-                                      */
-    virtual unsigned int element_multiplicity (const unsigned int index) const;
-    
-				     /**
-				      * This function returns
-				      * @p{true}, if the shape
-				      * function @p{shape_index} has
-				      * non-zero values on the face
-				      * @p{face_index}. For the lowest
-				      * order Nedelec elements, this
-				      * is actually the case for the
-				      * one on which the shape
-				      * function is defined and all
-				      * neighboring ones.
-				      *
-				      * Implementation of the
-				      * interface in
-				      * @ref{FiniteElement}
-				      */
-    virtual bool has_support_on_face (const unsigned int shape_index,
-				      const unsigned int face_index) 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;
-
-
-				     /**
-				      * Declare a nested class which
-				      * will hold static definitions
-				      * of various matrices such as
-				      * constraint and embedding
-				      * matrices. The definition of
-				      * the various static fields are
-				      * in the files
-				      * @p{fe_raviart_thomas_[23]d.cc}
-				      * in the source directory.
-				      */
-    struct Matrices
-    {
-					 /**
-					  * Embedding matrices. For
-					  * each element type (the
-					  * first index) there are as
-					  * many embedding matrices as
-					  * there are children per
-					  * cell. The first index
-					  * starts with linear
-					  * elements and goes up in
-					  * polynomial degree. The
-					  * array may grow in the
-					  * future with the number of
-					  * elements for which these
-					  * matrices have been
-					  * computed. If for some
-					  * element, the matrices have
-					  * not been computed then you
-					  * may use the element
-					  * nevertheless but can not
-					  * access the respective
-					  * fields.
-					  */
-	static const double * const
-	embedding[][GeometryInfo<dim>::children_per_cell];
-
-					 /**
-					  * Number of elements (first
-					  * index) the above field
-					  * has. Equals the highest
-					  * polynomial degree for
-					  * which the embedding
-					  * matrices have been
-					  * computed.
-					  */
-	static const unsigned int n_embedding_matrices;
-
-					 /**
-					  * As the
-					  * @p{embedding_matrices}
-					  * field, but for the
-					  * interface constraints. One
-					  * for each element for which
-					  * it has been computed.
-					  */
-	static const double * const constraint_matrices[];
-
-					 /**
-					  * Like
-					  * @p{n_embedding_matrices},
-					  * but for the number of
-					  * interface constraint
-					  * matrices.
-					  */
-	static const unsigned int n_constraint_matrices;
-    };
-				     /**
-				      * Exception
-				      */
-    DeclException0 (ExcNotUsefulInThisDimension);
-    
-  protected:    
-				     /**
-				      * @p{clone} function instead of
-				      * a copy constructor.
-				      *
-				      * This function is needed by the
-				      * constructors of @p{FESystem}.
-				      */
-    virtual FiniteElement<dim> * clone() const;
-  
-				     /**
-				      * Prepare internal data
-				      * structures and fill in values
-				      * independent of the cell.
-				      */
-    virtual
-    typename Mapping<dim>::InternalDataBase *
-    get_data (const UpdateFlags,
-	      const Mapping<dim>& mapping,
-	      const Quadrature<dim>& quadrature) const ;
-
-				     /**
-				      * Implementation of the same
-				      * function in
-				      * @ref{FiniteElement}.
-				      */
-    virtual void
-    fill_fe_values (const Mapping<dim> &mapping,
-		    const typename DoFHandler<dim>::cell_iterator &cell,
-		    const Quadrature<dim>                &quadrature,
-		    typename Mapping<dim>::InternalDataBase      &mapping_internal,
-		    typename Mapping<dim>::InternalDataBase      &fe_internal,
-		    FEValuesData<dim>& data) const;
-    
-				     /**
-				      * Implementation of the same
-				      * function in
-				      * @ref{FiniteElement}.
-				      */
-    virtual void
-    fill_fe_face_values (const Mapping<dim> &mapping,
-			 const typename DoFHandler<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 ;
-    
-				     /**
-				      * Implementation of the same
-				      * function in
-				      * @ref{FiniteElement}.
-				      */
-    virtual void
-    fill_fe_subface_values (const Mapping<dim> &mapping,
-			    const typename DoFHandler<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:
-    
-				     /**
-				      * Only for internal use. Its
-				      * full name is
-				      * @p{get_dofs_per_object_vector}
-				      * function and it creates the
-				      * @p{dofs_per_object} vector that is
-				      * needed within the constructor to
-				      * be passed to the constructor of
-				      * @p{FiniteElementData}.
-				      */
-    static std::vector<unsigned int> get_dpo_vector(const unsigned int degree);
-
-				     /**
-				      * Initialize the
-				      * @p{unit_support_points} field
-				      * of the @ref{FiniteElementBase}
-				      * class. Called from the
-				      * constructor.
-				      */
-    void initialize_unit_support_points ();
-
-				     /**
-				      * Initialize the
-				      * @p{unit_face_support_points} field
-				      * of the @ref{FiniteElementBase}
-				      * class. Called from the
-				      * constructor.
-				      */
-    void initialize_unit_face_support_points ();
-    
-				     /**
-				      * Given a set of flags indicating
-				      * what quantities are requested
-				      * from a @p{FEValues} object,
-				      * return which of these can be
-				      * precomputed once and for
-				      * all. Often, the values of
-				      * shape function at quadrature
-				      * points can be precomputed, for
-				      * example, in which case the
-				      * return value of this function
-				      * would be the logical and of
-				      * the input @p{flags} and
-				      * @p{update_values}.
-				      *
-				      * For the present kind of finite
-				      * element, this is exactly the
-				      * case.
-				      */
-    virtual UpdateFlags update_once (const UpdateFlags flags) const;
-  
-				     /**
-				      * This is the opposite to the
-				      * above function: given a set of
-				      * flags indicating what we want
-				      * to know, return which of these
-				      * need to be computed each time
-				      * we visit a new cell.
-				      *
-				      * If for the computation of one
-				      * quantity something else is
-				      * also required (for example, we
-				      * often need the covariant
-				      * transformation when gradients
-				      * need to be computed), include
-				      * this in the result as well.
-				      */
-    virtual UpdateFlags update_each (const UpdateFlags flags) const;
-    
-				     /**
-				      * Degree of the polynomials.
-				      */  
-    const unsigned int degree;
-
-				     /**
-				      * Fields of cell-independent data.
-				      *
-				      * For information about the
-				      * general purpose of this class,
-				      * see the documentation of the
-				      * base class.
-				      */
-    class InternalData : public FiniteElementBase<dim>::InternalDataBase
-    {
-      public:
-					 /**
-					  * Array with shape function
-					  * values in quadrature
-					  * points. There is one row
-					  * for each shape function,
-					  * containing values for each
-					  * quadrature point. Since
-					  * the shape functions are
-					  * vector-valued (with as
-					  * many components as there
-					  * are space dimensions), the
-					  * value is a tensor.
-					  *
-					  * In this array, we store
-					  * the values of the shape
-					  * function in the quadrature
-					  * points on the unit
-					  * cell. The transformation
-					  * to the real space cell is
-					  * then simply done by
-					  * multiplication with the
-					  * Jacobian of the mapping.
-					  */
-	Table<2,Tensor<1,dim> > shape_values;
-
-					 /**
-					  * Array with shape function
-					  * gradients in quadrature
-					  * points. There is one
-					  * row for each shape
-					  * function, containing
-					  * values for each quadrature
-					  * point.
-					  *
-					  * We store the gradients in
-					  * the quadrature points on
-					  * the unit cell. We then
-					  * only have to apply the
-					  * transformation (which is a
-					  * matrix-vector
-					  * multiplication) when
-					  * visiting an actual cell.
-					  */
-	Table<2,Tensor<2,dim> > shape_gradients;
-    };
-    
-				     /**
-				      * Allow access from other
-				      * dimensions.
-				      */
-    template <int dim1> friend class FE_RaviartThomas;
-};
-
-
-/* -------------- declaration of explicit specializations ------------- */
-
-template <> void FE_RaviartThomas<1>::initialize_unit_face_support_points ();
-
-// declaration of explicit specializations of member variables, if the
-// compiler allows us to do that (the standard says we must)
-#ifndef DEAL_II_MEMBER_VAR_SPECIALIZATION_BUG
-template <> 
-const double * const 
-FE_RaviartThomas<1>::Matrices::embedding[][GeometryInfo<1>::children_per_cell];
-
-template <>
-const unsigned int FE_RaviartThomas<1>::Matrices::n_embedding_matrices;
-
-template <>
-const double * const FE_RaviartThomas<1>::Matrices::constraint_matrices[];
-
-template <>
-const unsigned int FE_RaviartThomas<1>::Matrices::n_constraint_matrices;
-
-template <> 
-const double * const 
-FE_RaviartThomas<2>::Matrices::embedding[][GeometryInfo<2>::children_per_cell];
-
-template <>
-const unsigned int FE_RaviartThomas<2>::Matrices::n_embedding_matrices;
-
-template <>
-const double * const FE_RaviartThomas<2>::Matrices::constraint_matrices[];
-
-template <>
-const unsigned int FE_RaviartThomas<2>::Matrices::n_constraint_matrices;
-
-template <> 
-const double * const 
-FE_RaviartThomas<3>::Matrices::embedding[][GeometryInfo<3>::children_per_cell];
-
-template <>
-const unsigned int FE_RaviartThomas<3>::Matrices::n_embedding_matrices;
-
-template <>
-const double * const FE_RaviartThomas<3>::Matrices::constraint_matrices[];
-
-template <>
-const unsigned int FE_RaviartThomas<3>::Matrices::n_constraint_matrices;
-
-#endif
-
-#endif
diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc
deleted file mode 100644
index 64e44f5f5e..0000000000
--- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc
+++ /dev/null
@@ -1,1237 +0,0 @@
-//----------------------------------------------------------------
-//    $Id$
-//    Version: $Name$
-//
-//    Copyright (C) 2003 by the deal.II authors
-//
-//    This file is subject to QPL and may not be  distributed
-//    without copyright and license information. Please refer
-//    to the file deal.II/doc/license.html for the  text  and
-//    further information on this license.
-//
-//----------------------------------------------------------------
-
-#include <base/quadrature.h>
-#include <base/table.h>
-#include <grid/tria.h>
-#include <grid/tria_iterator.h>
-#include <dofs/dof_accessor.h>
-#include <fe/fe.h>
-#include <fe/mapping.h>
-#include <fe/fe_raviart_thomas.h>
-#include <fe/fe_values.h>
-
-
-template <int dim>
-FE_RaviartThomas<dim>::FE_RaviartThomas (const unsigned int degree)
-		:
-		FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree),
-							   dim),
-				    std::vector<bool> (FiniteElementData<dim>(get_dpo_vector(degree),dim).dofs_per_cell,false),
-				    std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),dim).dofs_per_cell,
-								    std::vector<bool>(dim,true))),
-		degree(degree)
-{
-  Assert (dim >= 2, ExcNotUsefulInThisDimension());
-  
-				   // copy constraint matrices if they
-				   // are defined. otherwise leave
-				   // them at zero size
-  if (degree<Matrices::n_constraint_matrices+1)
-    {
-      this->interface_constraints.
-        TableBase<2,double>::reinit (this->interface_constraints_size());
-      this->interface_constraints.fill (Matrices::constraint_matrices[degree-1]);
-    };
-
-				   // next copy over embedding
-				   // matrices if they are defined
-  if ((degree < Matrices::n_embedding_matrices+1) &&
-      (Matrices::embedding[degree-1][0] != 0))
-    for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
-      {
-                                         // copy
-        this->prolongation[c].reinit (this->dofs_per_cell,
-                                      this->dofs_per_cell);
-        this->prolongation[c].fill (Matrices::embedding[degree-1][c]);
-                                         // and make sure that the row
-                                         // sum is 0.5 (for usual
-                                         // elements, the row sum must
-                                         // be 1, but here the shape
-                                         // function is multiplied by
-                                         // the inverse of the
-                                         // Jacobian, which introduces
-                                         // a factor of 1/2 when going
-                                         // from mother to child)
-        for (unsigned int row=0; row<this->dofs_per_cell; ++row)
-          {
-            double sum = 0;
-            for (unsigned int col=0; col<this->dofs_per_cell; ++col)
-              sum += this->prolongation[c](row,col);
-            Assert (std::fabs(sum-.5) < 1e-14,
-                    ExcInternalError());
-          };
-      };
-
-				   // then fill restriction
-				   // matrices. they are hardcoded for
-				   // the first few elements
-  switch (dim)
-    {
-      case 2:   // 2d
-      {
-	switch (degree)
-	  {
-	    case 1:
-	    {
-                                               // 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;
-	    };
-	    
-	    default:
-	    {
-					       // 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
-              break;
-	    };
-	  };
-	
-	break;
-      };
-
-
-      case 3:   // 3d
-      {
-	switch (degree)
-	  {
-	    case 1:
-	    {
-					       // same principle as in
-					       // 2d, take one child
-					       // cell to get at the
-					       // values of each of
-					       // the 12 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[0](3,3) = 2.;
-	      this->restriction[1](1,1) = 2.;
-	      this->restriction[3](2,2) = 2.;
-              
-	      this->restriction[4](4,4) = 2.;
-	      this->restriction[4](7,7) = 2.;
-	      this->restriction[5](5,5) = 2.;
-	      this->restriction[7](6,6) = 2.;
-              
-	      this->restriction[0](8,8) = 2.;
-	      this->restriction[1](9,9) = 2.;
-	      this->restriction[2](10,10) = 2.;
-	      this->restriction[3](11,11) = 2.;
-              
-	      break;
-	    };
-	    
-	    default:
-	    {
-					       // 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
-              break;
-	    };
-	  };
-	
-	break;
-      };
-      
-      default:
-	    Assert (false,ExcNotImplemented());
-    }
-
-				   // finally fill in support points
-				   // on cell and face
-  initialize_unit_support_points ();
-  initialize_unit_face_support_points ();
-
-                                   // then make
-                                   // system_to_component_table
-                                   // invalid, since this has no
-                                   // meaning for the present element
-  std::vector<std::pair<unsigned,unsigned> > tmp1, tmp2;
-  this->system_to_component_table.swap (tmp1);
-  this->face_system_to_component_table.swap (tmp2);
-}
-
-
-
-template <int dim>
-FiniteElement<dim> *
-FE_RaviartThomas<dim>::clone() const
-{
-  return new FE_RaviartThomas<dim>(degree);
-}
-
-
-#if deal_II_dimension == 1
-
-template <>
-double
-FE_RaviartThomas<1>::shape_value_component (const unsigned int ,
-                                            const Point<1>    &,
-                                            const unsigned int ) const
-{
-  Assert (false, ExcNotImplemented());
-  return 0.;
-}
-
-#endif
-
-#if deal_II_dimension == 2
-
-template <>
-double
-FE_RaviartThomas<2>::shape_value_component (const unsigned int i,
-                                            const Point<2>    &p,
-                                            const unsigned int component) const
-{
-  const unsigned int dim = 2;
-  
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  Assert (component < dim, ExcIndexRange (component, 0, dim));
-  
-  switch (degree)
-    {
-				       // first order Raviart-Thomas elements
-      case 1:
-      {
-	switch (i)
-	  {
-					     // (0, 1-y)
-	    case 0: return (component == 0 ? 0: 1-p(1));
-						   // (x,0)
-	    case 1: return (component == 0 ? p(0) : 0);
-						   // (0, y)
-	    case 2: return (component == 0 ? 0: p(1));
-						   // (1-x, 0)
-	    case 3: return (component == 0 ? 1-p(0) : 0);
-                        
-						   // there are only
-						   // four shape
-						   // functions!?
-	    default:
-		  Assert (false, ExcInternalError());
-		  return 0;
-	  };
-      };
-
-                                       // no other degrees
-                                       // implemented
-      default:
-	    Assert (false, ExcNotImplemented());
-    };
-  
-  return 0;
-}
-
-#endif
-
-#if deal_II_dimension == 3
-
-template <>
-double
-FE_RaviartThomas<3>::shape_value_component (const unsigned int i,
-                                            const Point<3>    &/*p*/,
-                                            const unsigned int component) const
-{
-  const unsigned int dim = 3;
-  
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  Assert (component < dim, ExcIndexRange (component, 0, dim));
-  
-  switch (degree)
-    {
-                                       // no other degrees
-                                       // implemented
-      default:
-	    Assert (false, ExcNotImplemented());
-    };
-  
-  return 0;
-}
-
-#endif
-
-#if deal_II_dimension == 1
-
-template <>
-Tensor<1,1>
-FE_RaviartThomas<1>::shape_grad_component (const unsigned int ,
-                                           const Point<1>    &,
-                                           const unsigned int ) const
-{
-  Assert (false, ExcNotImplemented());
-  return Tensor<1,1>();
-}
-
-#endif
-
-#if deal_II_dimension == 2
-
-template <>
-Tensor<1,2>
-FE_RaviartThomas<2>::shape_grad_component (const unsigned int i,
-                                           const Point<2>    &,
-                                           const unsigned int component) const
-{
-  const unsigned int dim = 2;
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  Assert (component < dim, ExcIndexRange (component, 0, dim));
-
-  switch (degree)
-    {
-				       // first order Raviart-Thomas elements
-      case 1:
-      {
-					 // on the unit cell, the
-					 // gradients of these shape
-					 // functions are constant, so
-					 // we pack them into a table
-					 // for simpler lookup
-					 //
-					 // the format is: first
-					 // index=shape function
-					 // number; second
-					 // index=vector component,
-					 // third index=component
-					 // within gradient
-	static const double unit_gradients[4][2][2]
-	  = { { {0.,0.} , {0.,-1.} },
-	      { {1.,0.} , {0.,0.}  },
-	      { {0.,0.} , {0.,+1.} },
-	      { {-1.,0.}, {0.,0.}  } };
-	return Tensor<1,dim>(unit_gradients[i][component]);
-      };
-
-                                       // no other degrees
-                                       // implemented
-      default:
-	    Assert (false, ExcNotImplemented());
-    };
-  
-  return Tensor<1,dim>();
-}
-
-#endif
-
-#if deal_II_dimension == 3
-
-template <>
-Tensor<1,3>
-FE_RaviartThomas<3>::shape_grad_component (const unsigned int i,
-                                           const Point<3>    &/*p*/,
-                                           const unsigned int component) const
-{
-  const unsigned int dim = 3;
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  Assert (component < dim, ExcIndexRange (component, 0, dim));
-
-  switch (degree)
-    {
-                                       // no other degrees
-                                       // implemented
-      default:
-	    Assert (false, ExcNotImplemented());
-    };
-  
-  return Tensor<1,dim>();
-}
-
-#endif
-
-
-#if deal_II_dimension == 1
-
-template <>
-Tensor<2,1>
-FE_RaviartThomas<1>::shape_grad_grad_component (const unsigned int ,
-                                                const Point<1>    &,
-                                                const unsigned int ) const
-{
-  Assert (false, ExcNotImplemented());
-  return Tensor<2,1>();
-}
-
-#endif
-
-
-#if deal_II_dimension == 2
-
-template <>
-Tensor<2,2>
-FE_RaviartThomas<2>::shape_grad_grad_component (const unsigned int i,
-                                                const Point<2> &/*p*/,
-                                                const unsigned int component) const
-{
-  const unsigned int dim = 2;
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  Assert (component < dim, ExcIndexRange (component, 0, dim));
-
-  switch (degree)
-    {
-				       // first order Raviart-Thomas
-				       // elements. their second
-				       // derivatives on the unit cell
-				       // are zero
-      case 1:
-      {
-	return Tensor<2,dim>();
-      };
-
-                                       // no other degrees
-                                       // implemented
-      default:
-	    Assert (false, ExcNotImplemented());
-    };
-
-  return Tensor<2,dim>();
-}
-
-#endif
-
-#if deal_II_dimension == 3
-
-template <>
-Tensor<2,3>
-FE_RaviartThomas<3>::shape_grad_grad_component (const unsigned int i,
-                                                const Point<3>    &/*p*/,
-                                                const unsigned int component) const
-{
-  const unsigned int dim = 3;
-  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
-  Assert (component < dim, ExcIndexRange (component, 0, dim));
-
-  switch (degree)
-    {
-                                       // no other degrees
-                                       // implemented
-      default:
-	    Assert (false, ExcNotImplemented());
-    };
-
-  return Tensor<2,dim>();
-}
-
-#endif
-
-//----------------------------------------------------------------------
-// Auxiliary functions
-//----------------------------------------------------------------------
-
-
-
-template <int dim>
-void FE_RaviartThomas<dim>::initialize_unit_support_points ()
-{
-  switch (degree)
-    {
-      case 1:
-      {
-					 // all degrees of freedom are
-					 // on edges, and their order
-					 // is the same as the edges
-					 // themselves
-	this->unit_support_points.resize(GeometryInfo<dim>::lines_per_cell);
-	for (unsigned int line=0; line<GeometryInfo<dim>::lines_per_cell; ++line)
-	  {
-	    const unsigned int
-	      vertex_index_0 = GeometryInfo<dim>::vertices_adjacent_to_line(line,0),
-	      vertex_index_1 = GeometryInfo<dim>::vertices_adjacent_to_line(line,1);
-	    
-	    const Point<dim>
-	      vertex_0 = GeometryInfo<dim>::unit_cell_vertex(vertex_index_0),
-	      vertex_1 = GeometryInfo<dim>::unit_cell_vertex(vertex_index_1);
-	    
-					     // place dofs right
-					     // between the vertices
-					     // of each line
-	    this->unit_support_points[line] = (vertex_0 + vertex_1) / 2;
-	  };
-	    
-	break;
-      };
-
-      default:
-					     // no higher order
-					     // elements implemented
-					     // right now
-	    Assert (false, ExcNotImplemented());
-    };
-}
-
-
-#if deal_II_dimension == 1
-
-template <>
-void FE_RaviartThomas<1>::initialize_unit_face_support_points ()
-{
-				   // no faces in 1d, so nothing to do
-}
-
-#endif
-
-
-template <int dim>
-void FE_RaviartThomas<dim>::initialize_unit_face_support_points ()
-{
-  switch (degree)
-    {
-      case 1:
-      {
-					 // do this the same as above, but
-					 // for one dimension less
-	this->unit_face_support_points.resize(GeometryInfo<dim-1>::lines_per_cell);
-	for (unsigned int line=0; line<GeometryInfo<dim-1>::lines_per_cell; ++line)
-	  {
-	    const unsigned int
-	      vertex_index_0 = GeometryInfo<dim-1>::vertices_adjacent_to_line(line,0),
-	      vertex_index_1 = GeometryInfo<dim-1>::vertices_adjacent_to_line(line,1);
-      
-	    const Point<dim-1>
-	      vertex_0 = GeometryInfo<dim-1>::unit_cell_vertex(vertex_index_0),
-	      vertex_1 = GeometryInfo<dim-1>::unit_cell_vertex(vertex_index_1);
-
-					     // place dofs right
-					     // between the vertices of each
-					     // line
-            this->unit_face_support_points[line] = (vertex_0 + vertex_1) / 2;
-	  };
-	break;
-      };
-
-      default:
-					     // no higher order
-					     // elements implemented
-					     // right now
-	    Assert (false, ExcNotImplemented());
-    };	    
-}
-
-
-
-template <int dim>
-std::vector<unsigned int>
-FE_RaviartThomas<dim>::get_dpo_vector(const unsigned int degree)
-{
-  Assert (degree == 1, ExcNotImplemented());
-
-				   // for degree==1, put all degrees
-				   // of freedom on the lines, and in
-				   // particular @p{degree} DoFs per
-				   // line:
-  std::vector<unsigned int> dpo(dim+1, 0U);
-  dpo[1] = degree;
-
-  return dpo;
-}
-
-
-
-template <int dim>
-UpdateFlags
-FE_RaviartThomas<dim>::update_once (const UpdateFlags) const
-{
-				   // even the values have to be
-				   // computed on the real cell, so
-				   // nothing can be done in advance
-  return update_default;
-}
-
-
-
-template <int dim>
-UpdateFlags
-FE_RaviartThomas<dim>::update_each (const UpdateFlags flags) const
-{
-  UpdateFlags out = update_default;
-
-  if (flags & update_values)
-    out |= update_values             | update_covariant_transformation;
-  if (flags & update_gradients)
-    out |= update_gradients          | update_covariant_transformation;
-  if (flags & update_second_derivatives)
-    out |= update_second_derivatives | update_covariant_transformation;
-
-  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.reinit (this->dofs_per_cell, n_q_points);
-
-  if (flags & update_gradients)
-    data->shape_gradients.reinit (this->dofs_per_cell, 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 DoFHandler<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(&*shape_values.begin(),
-                                      &*shape_values.end(),
-                                      fe_data.shape_values[k].begin(),
-                                      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.n_rows() == this->dofs_per_cell * dim,
-	      ExcInternalError());
-      Assert (data.shape_gradients.n_cols() == 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(&*shape_grads1.begin(),
-                                      &*shape_grads1.end(),
-                                      fe_data.shape_gradients[k].begin(),
-                                      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_grads1.begin(),
-                                      &*shape_grads1.end(),
-                                      &*shape_grads2.begin(),
-                                      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, 0, mapping_data, fe_data, data);
-}
-
-
-
-template <int dim>
-void
-FE_RaviartThomas<dim>::fill_fe_face_values (const Mapping<dim>                   &mapping,
-                                            const typename DoFHandler<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 unsigned int offset = 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.n_cols() ==
-              GeometryInfo<dim>::faces_per_cell * 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(&*shape_values.begin(),
-                                      &*shape_values.end(),
-                                      fe_data.shape_values[k].begin()+offset,
-                                      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.n_cols() ==
-              GeometryInfo<dim>::faces_per_cell * 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.n_rows() == this->dofs_per_cell * dim,
-	      ExcInternalError());
-      Assert (data.shape_gradients.n_cols() == 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(&*shape_grads1.begin(),
-                                      &*shape_grads1.end(),
-                                      fe_data.shape_gradients[k].begin()+offset,
-                                      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_grads1.begin(),
-                                      &*shape_grads1.end(),
-                                      &*shape_grads2.begin(),
-                                      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 DoFHandler<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 unsigned int offset = ((face * GeometryInfo<dim>::subfaces_per_face + subface)
-                               * 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.n_cols() ==
-              GeometryInfo<dim>::faces_per_cell * 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(&*shape_values.begin(),
-                                      &*shape_values.end(),
-                                      fe_data.shape_values[k].begin()+offset,
-                                      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.n_cols() ==
-              GeometryInfo<dim>::faces_per_cell * 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.n_rows() == this->dofs_per_cell * dim,
-	      ExcInternalError());
-      Assert (data.shape_gradients.n_cols() == 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(&*shape_grads1.begin(),
-                                      &*shape_grads1.end(),
-                                      fe_data.shape_gradients[k].begin()+offset,
-                                      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_grads1.begin(),
-                                      &*shape_grads1.end(),
-                                      &*shape_grads2.begin(),
-                                      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>
-unsigned int
-FE_RaviartThomas<dim>::n_base_elements () const
-{
-  return 1;
-}
-
-
-
-template <int dim>
-const FiniteElement<dim> &
-FE_RaviartThomas<dim>::base_element (const unsigned int index) const
-{
-  Assert (index==0, ExcIndexRange(index, 0, 1));
-  return *this;
-}
-
-
-
-template <int dim>
-unsigned int
-FE_RaviartThomas<dim>::element_multiplicity (const unsigned int index) const
-{
-  Assert (index==0, ExcIndexRange(index, 0, 1));
-  return 1;
-}
-
-
-
-template <int dim>
-bool
-FE_RaviartThomas<dim>::has_support_on_face (const unsigned int shape_index,
-                                            const unsigned int face_index) const
-{
-  Assert (shape_index < this->dofs_per_cell,
-	  ExcIndexRange (shape_index, 0, this->dofs_per_cell));
-  Assert (face_index < GeometryInfo<dim>::faces_per_cell,
-	  ExcIndexRange (face_index, 0, GeometryInfo<dim>::faces_per_cell));
-
-  switch (degree)
-    {
-      case 1:
-      {
-        switch (dim)
-          {
-            case 2:
-            {
-                                               // only on the one
-                                               // non-adjacent face
-                                               // are the values
-                                               // actually zero. list
-                                               // these in a table
-              const unsigned int
-                opposite_faces[GeometryInfo<2>::faces_per_cell]
-                = { 2, 3, 0, 1};
-              
-              return (face_index != opposite_faces[shape_index]);
-            };
-            
-            default: Assert (false, ExcNotImplemented());
-          };
-      };
-      
-      default:  // other degree
-            Assert (false, ExcNotImplemented());
-    };
-  
-  return true;
-}
-
-
-
-template <int dim>
-unsigned int
-FE_RaviartThomas<dim>::memory_consumption () const
-{
-  Assert (false, ExcNotImplemented ());
-  return 0;
-}
-
-
-
-template <int dim>
-unsigned int
-FE_RaviartThomas<dim>::get_degree () const
-{
-  return degree;
-}
-
-
-
-template class FE_RaviartThomas<deal_II_dimension>;
diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas_1d.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas_1d.cc
deleted file mode 100644
index a47a21a4a6..0000000000
--- a/deal.II/deal.II/source/fe/fe_raviart_thomas_1d.cc
+++ /dev/null
@@ -1,53 +0,0 @@
-//----------------------------------------------------------------
-//    $Id$
-//    Version: $Name$
-//
-//    Copyright (C) 2003 by the deal.II authors
-//
-//    This file is subject to QPL and may not be  distributed
-//    without copyright and license information. Please refer
-//    to the file deal.II/doc/license.html for the  text  and
-//    further information on this license.
-//
-//----------------------------------------------------------------
-
-
-
-// only compile this file if in 1d. note that Raviart-Thomas elements
-// do not make much sense in 1d anyway, so this file only contains
-// dummy implementations to avoid linker errors due to missing symbols
-#if deal_II_dimension == 1
-
-
-#include <fe/fe_raviart_thomas.h>
-
-
-template <>
-const double * const
-FE_RaviartThomas<1>::Matrices::embedding[][GeometryInfo<1>::children_per_cell] =
-{};
-
-
-template <>
-const unsigned int
-FE_RaviartThomas<1>::Matrices::n_embedding_matrices = 0;
-
-
-
-// No constraints in 1d
-template <>
-const unsigned int 
-FE_RaviartThomas<1>::Matrices::n_constraint_matrices = 0;
-
-
-template <>
-const double * const
-FE_RaviartThomas<1>::Matrices::constraint_matrices[] = {};
-
-
-#else // #if deal_II_dimension
-// On gcc2.95 on Alpha OSF1, the native assembler does not like empty
-// files, so provide some dummy code
-namespace { void dummy () {} }
-#endif // #if deal_II_dimension == 1
-
diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas_2d.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas_2d.cc
deleted file mode 100644
index 70e9484b54..0000000000
--- a/deal.II/deal.II/source/fe/fe_raviart_thomas_2d.cc
+++ /dev/null
@@ -1,137 +0,0 @@
-//----------------------------------------------------------------
-//    $Id$
-//    Version: $Name$
-//
-//    Copyright (C) 2003 by the deal.II authors
-//
-//    This file is subject to QPL and may not be  distributed
-//    without copyright and license information. Please refer
-//    to the file deal.II/doc/license.html for the  text  and
-//    further information on this license.
-//
-//----------------------------------------------------------------
-
-
-// only compile this file if in 2d
-#if deal_II_dimension == 2
-
-
-#include <fe/fe_raviart_thomas.h>
-
-// Transfer matrices for finite elements: have one matrix for each of
-// the four child cells which tells us how the degrees of freedom on
-// the child cell are obtained from the degrees of freedom on the
-// mother cell
-//
-// note the following: since the shape functions themselves and not
-// only the gradients are transformed using the mapping object from
-// the unit cell to the real cell, the actual values of the function
-// on the real cell is degree of freedom times value of the shape
-// function on the unit cell times inverse Jacobian. Thus, what has
-// the DoF value 1 on the mother cell must have the DoF value 1/2 on
-// the child cell since the latter is smaller by a (linear scaling)
-// factor of two.
-namespace FE_RaviartThomas_2d
-{
-  static const double q1_into_q1_refined_0[] =
-  {
-	.5,   0,   0 ,  0,
-	0,    0.25,0,   0.25,
-	0.25, 0,   0.25,0,
-	0,    0,   0,   .5 
-  };
-
-  static const double q1_into_q1_refined_1[] =
-  {
-  	.5,   0.,   0.,   0.,
-  	0.,   .5,   0.,   0.,
-  	0.25, 0.,   0.25, 0.,
-	0.,   0.25, 0.,   0.25,
-  };
-
-  static const double q1_into_q1_refined_2[] =
-  {
-  	0.25, 0.,   0.25, 0.,
- 	0.,   .5,   0.,   0.,
-	0.,   0.,   .5,   0.,
-  	0.,   0.25, 0.,   0.25,
-  };
-
-  static const double q1_into_q1_refined_3[] =
-  {
-  	0.25, 0.,   0.25, 0.,
-  	0.,   0.25, 0.,   0.25,
-	0.,   0.,   .5,   0.,
-  	0.,   0.,   0.,   .5,
-  };
-}  // namespace FE_RaviartThomas_2d
-
-
-// embedding matrices
-
-template <>
-const double * const 
-FE_RaviartThomas<2>::Matrices::embedding[][GeometryInfo<2>::children_per_cell] =
-{
-      { FE_RaviartThomas_2d::q1_into_q1_refined_0, FE_RaviartThomas_2d::q1_into_q1_refined_1,
-	FE_RaviartThomas_2d::q1_into_q1_refined_2, FE_RaviartThomas_2d::q1_into_q1_refined_3 }
-};
-
-
-template <>
-const unsigned int
-FE_RaviartThomas<2>::Matrices::n_embedding_matrices
-= sizeof(FE_RaviartThomas<2>::Matrices::embedding) /
-sizeof(FE_RaviartThomas<2>::Matrices::embedding[0]);
-
-
-// Constraint matrices: how do the new value on child faces depend on
-// the values on the mother face if that face has a hanging node
-//
-// Here, the same applies as for the embedding matrices: since the DoF
-// values are not only multiplied by the values of the shape function
-// on the unit cell, but also by the transformation, we have to
-// multiply the value on the large face by 1/2 to get the same value
-// back on the small face.  in other words, if a DoF has weight 1 on
-// the big cell, then it has to have weight 1/2 on the small ones, in
-// order to give the same value of the shape function in real space
-namespace FE_RaviartThomas_2d 
-{
-  static const double constraint_q1[] =
-  {
-					 // the function is constant
-					 // along each edge, so each
-					 // degree of freedom on the
-					 // refined edge has the same
-					 // value as that on the
-					 // coarse edge, modulo the
-					 // issue with the
-					 // transformation described
-					 // above
-  	1./2., 1./2.
-  };
-
-}
-
-
-template <>
-const double * const 
-FE_RaviartThomas<2>::Matrices::constraint_matrices[] =
-{
-      FE_RaviartThomas_2d::constraint_q1
-};
-
-
-template <>
-const unsigned int 
-FE_RaviartThomas<2>::Matrices::n_constraint_matrices
-= sizeof(FE_RaviartThomas<2>::Matrices::constraint_matrices) /
-sizeof(FE_RaviartThomas<2>::Matrices::constraint_matrices[0]);
-
-
-
-#else // #if deal_II_dimension
-// On gcc2.95 on Alpha OSF1, the native assembler does not like empty
-// files, so provide some dummy code
-namespace { void dummy () {} }
-#endif // #if deal_II_dimension == 2
diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas_3d.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas_3d.cc
deleted file mode 100644
index b5d60f4c43..0000000000
--- a/deal.II/deal.II/source/fe/fe_raviart_thomas_3d.cc
+++ /dev/null
@@ -1,242 +0,0 @@
-//----------------------------------------------------------------
-//    $Id$
-//    Version: $Name$
-//
-//    Copyright (C) 2003 by the deal.II authors
-//
-//    This file is subject to QPL and may not be  distributed
-//    without copyright and license information. Please refer
-//    to the file deal.II/doc/license.html for the  text  and
-//    further information on this license.
-//
-//----------------------------------------------------------------
-
-// Transfer matrices for finite elements
-
-
-// only compile this file if in 3d
-#if deal_II_dimension == 3
-
-#include <fe/fe_raviart_thomas.h>
-
-// Transfer matrices for finite elements: have one matrix for each of
-// the four child cells which tells us how the degrees of freedom on
-// the child cell are obtained from the degrees of freedom on the
-// mother cell
-//
-// note the following: since the shape functions themselves and not
-// only the gradients are transformed using the mapping object from
-// the unit cell to the real cell, the actual values of the function
-// on the real cell is degree of freedom times value of the shape
-// function on the unit cell times Jacobian. Thus, what has the DoF
-// value 1 on the mother cell must have the DoF value 2 on the child
-// cell since the latter is smaller by a (linear scaling) factor of
-// two.
-namespace FE_RaviartThomas_3d
-{
-  static const double q1_into_q1_refined_0[] =
-  {
-	.5, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
-	0., 0.25, 0., 0.25, 0., 0., 0., 0.,0.,0.,0.,0.,
-	0.25, 0., 0.25, 0., 0., 0., 0., 0.,0.,0.,0.,0.,
-	0., 0., 0., .5, 0., 0., 0., 0.,0.,0.,0.,0.,
-	0.25, 0., 0., 0., 0.25, 0., 0., 0.,0.,0.,0.,0.,
-	0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0.,
-	0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0.,
-	0., 0., 0., 0.25, 0., 0., 0., 0.25,0.,0.,0.,0.,
-	0., 0., 0., 0., 0., 0., 0., 0., .5, 0., 0., 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125,
-	0., 0., 0., 0., 0., 0., 0., 0., 0.25, .0, 0., 0.25,
-  };
-
-  static const double q1_into_q1_refined_1[] =
-  {
-	.5, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
-	0., .5, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
-	0.25, 0., 0.25, 0., 0., 0., 0., 0.,0.,0.,0.,0.,
-	0., 0.25, 0., 0.25, 0., 0., 0., 0.,0.,0.,0.,0.,
-	0.25, 0., 0., 0., 0.25, 0., 0., 0.,0.,0.,0.,0.,
-	0., 0.25, 0., 0., 0., 0.25, 0., 0.,0.,0.,0.,0.,
-	0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0.,
-	0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0., .5, 0., 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125,
-  };
-
-  static const double q1_into_q1_refined_2[] =
-  {
-	0.25, 0., 0.25, 0., 0., 0., 0., 0.,0.,0.,0.,0.,
-	0., .5, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
-	0., 0., .5, 0., 0., 0., 0., 0., 0., 0., 0., 0.,
-	0., 0.25, 0., 0.25, 0., 0., 0., 0.,0.,0.,0.,0.,
-	0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0.,
-	0., 0.25, 0., 0., 0., 0.25, 0., 0.,0.,0.,0.,0.,
-	0., 0., 0.25, 0., 0., 0., 0.25, 0.,0.,0.,0.,0.,
-	0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125,
-	0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .5, 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25,
-  };
-
-  static const double q1_into_q1_refined_3[] =
-  {
-	0.25, 0., 0.25, 0., 0., 0., 0., 0.,0.,0.,0.,0.,
-	0., 0.25, 0., 0.25, 0., 0., 0., 0.,0.,0.,0.,0.,
-	0., 0., .5, 0., 0., 0., 0., 0., 0., 0., 0., 0.,
-	0., 0., 0., .5, 0., 0., 0., 0., 0., 0., 0., 0.,
-	0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0.,
-	0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0.,
-	0., 0., 0.25, 0., 0., 0., 0.25, 0.,0.,0.,0.,0.,
-	0., 0., 0., 0.25, 0., 0., 0., 0.25,0.,0.,0.,0.,
-	0., 0., 0., 0., 0., 0., 0., 0.,0.25,0.,0.,0.25,
-	0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125,
-	0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25,
-	0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .5,
-  };
-
-  static const double q1_into_q1_refined_4[] =
-  {
-	0.25, 0., 0., 0., 0.25, 0., 0., 0.,0.,0.,0.,0.,
-	0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0.,
-	0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0.,
-	0., 0., 0., 0.25, 0., 0., 0., 0.25, 0., 0., 0., 0.,
-	0., 0., 0., 0., .5, 0., 0., 0., 0., 0., 0., 0.,
-	0., 0., 0., 0., 0., 0.25, 0., 0.25,0.,0.,0.,0.,
-	0., 0., 0., 0., 0.25, 0., 0.25, 0.,0.,0.,0.,0.,
-	0., 0., 0., 0., 0., 0., 0., .5, 0., 0., 0., 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., .5, 0., 0., 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125,
-	0., 0., 0., 0., 0., 0., 0., 0.,0.25,0.,0.,0.25,
-  };
-
-  static const double q1_into_q1_refined_5[] =
-  { 
-	0.25, 0., 0., 0., 0.25, 0., 0., 0.,0.,0.,0.,0.,
-	0., 0.25, 0., 0., 0., 0.25, 0., 0.,0.,0.,0.,0.,
-	0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0.,
-	0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0.,
-	0., 0., 0., 0., .5, 0., 0., 0., 0., 0., 0., 0.,
-	0., 0., 0., 0., 0., .5, 0., 0., 0., 0., 0., 0.,
-	0., 0., 0., 0., 0.25, 0., 0.25, 0.,0.,0.,0.,0.,
-	0., 0., 0., 0., 0., 0.25, 0., 0.25,0.,0.,0.,0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0., 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0., .5, 0., 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125,
-
-  };
-
-  static const double q1_into_q1_refined_6[] =
-  {
-	0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0.,
-	0., 0.25, 0., 0., 0., 0.25, 0., 0.,0.,0.,0.,0.,
-	0., 0., 0.25, 0., 0., 0., 0.25, 0.,0.,0.,0.,0.,
-	0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0.,
-	0., 0., 0., 0., 0.25, 0., 0.25, 0.,0.,0.,0.,0.,
-	0., 0., 0., 0., 0., .5, 0., 0., 0., 0., 0., 0.,
-	0., 0., 0., 0., 0., 0., .5, 0., 0., 0., 0., 0.,
-	0., 0., 0., 0., 0., 0.25, 0., 0.25, 0., 0., 0., 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125,
-	0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25, 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .5, 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25,
-	
-
-  };
-
-  static const double q1_into_q1_refined_7[] =
-  {
-	0.125, 0., 0.125, 0., 0.125, 0., 0.125, 0.,0.,0.,0.,0.,
-	0., 0.125, 0., 0.125, 0., 0.125, 0., 0.125,0.,0.,0.,0.,
-	0., 0., 0.25, 0., 0., 0., 0.25, 0.,0.,0.,0.,0.,
-	0., 0., 0., 0.25, 0., 0., 0., 0.25,0.,0.,0.,0.,
-	0., 0., 0., 0., 0.25, 0., 0.25, 0.,0.,0.,0.,0.,
-	0., 0., 0., 0., 0., 0.25, 0., 0.25,0.,0.,0.,0.,
-	0., 0., 0., 0., 0., 0., .5, 0., 0., 0., 0., 0.,
-	0., 0., 0., 0., 0., 0., 0., .5, 0., 0., 0., 0.,
-	0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0., 0., 0.25,
-	0., 0., 0., 0., 0., 0., 0., 0., 0.125, 0.125, 0.125, 0.125,
-	0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.25, 0.25,
-	0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5,
-  };
-
-}  // namespace FE_RaviartThomas_3d
-
-
-// embedding matrices
-
-template <>
-const double * const 
-FE_RaviartThomas<3>::Matrices::embedding[][GeometryInfo<3>::children_per_cell] =
-{
-      { FE_RaviartThomas_3d::q1_into_q1_refined_0, FE_RaviartThomas_3d::q1_into_q1_refined_1,
-	FE_RaviartThomas_3d::q1_into_q1_refined_2, FE_RaviartThomas_3d::q1_into_q1_refined_3,
-	FE_RaviartThomas_3d::q1_into_q1_refined_4, FE_RaviartThomas_3d::q1_into_q1_refined_5,
-	FE_RaviartThomas_3d::q1_into_q1_refined_6, FE_RaviartThomas_3d::q1_into_q1_refined_7 }
-};
-
-
-template <>
-const unsigned int
-FE_RaviartThomas<3>::Matrices::n_embedding_matrices
-= sizeof(FE_RaviartThomas<3>::Matrices::embedding) /
-sizeof(FE_RaviartThomas<3>::Matrices::embedding[0]);
-
-
-
-// Constraint matrices: how do the new value on child faces depend on
-// the values on the mother face if that face has a hanging node
-//
-// Here, the same applies as for the embedding matrices: since the DoF
-// values are not only multiplied by the values of the shape function
-// on the unit cell, but also by the transformation, we have to
-// multiply the value on the large face by 1/2 to get the same value
-// back on the small face
-namespace FE_RaviartThomas_3d 
-{
-  static const double constraint_q1[] =
-  {
-	0, .25, 0, .25,  // first the four interior lines
-	.25, 0, .25, 0,
-	0, .25, 0, .25,
-	.25, 0, .25, 0,
-	.5, 0, 0, 0,  // then the two child lines of each of the four outer
-	.5, 0, 0, 0,  // ones. since the shape functions are constant on each
-	0, .5, 0, 0,  // line, the two child lines get the same weights, modulo
-	0, .5, 0, 0,  // the issue with the division by length scaling
-	0, 0, .5, 0,
-	0, 0, .5, 0,
-	0, 0, 0, .5,
-	0, 0, 0, .5
-  };
-}
-
-
-
-template <>
-const double * const 
-FE_RaviartThomas<3>::Matrices::constraint_matrices[] =
-{
-      FE_RaviartThomas_3d::constraint_q1
-};
-
-
-
-template <>
-const unsigned int 
-FE_RaviartThomas<3>::Matrices::n_constraint_matrices
-= sizeof(FE_RaviartThomas<3>::Matrices::constraint_matrices) /
-sizeof(FE_RaviartThomas<3>::Matrices::constraint_matrices[0]);
-
-
-
-#else // #if deal_II_dimension
-// On gcc2.95 on Alpha OSF1, the native assembler does not like empty
-// files, so provide some dummy code
-namespace { void dummy () {} }
-#endif // #if deal_II_dimension == 3