From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 7 Apr 2003 15:43:52 +0000 (+0000)
Subject: Add initial versions of these files.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c417172df106750daeea12a79c8b6c92b51bee33;p=dealii-svn.git

Add initial versions of these files.


git-svn-id: https://svn.dealii.org/trunk@7370 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
new file mode 100644
index 0000000000..73ee24cc7a
--- /dev/null
+++ b/deal.II/deal.II/include/fe/fe_raviart_thomas.h
@@ -0,0 +1,581 @@
+//---------------------------------------------------------------
+//    $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
new file mode 100644
index 0000000000..64e44f5f5e
--- /dev/null
+++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc
@@ -0,0 +1,1237 @@
+//----------------------------------------------------------------
+//    $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
new file mode 100644
index 0000000000..a47a21a4a6
--- /dev/null
+++ b/deal.II/deal.II/source/fe/fe_raviart_thomas_1d.cc
@@ -0,0 +1,53 @@
+//----------------------------------------------------------------
+//    $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
new file mode 100644
index 0000000000..70e9484b54
--- /dev/null
+++ b/deal.II/deal.II/source/fe/fe_raviart_thomas_2d.cc
@@ -0,0 +1,137 @@
+//----------------------------------------------------------------
+//    $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
new file mode 100644
index 0000000000..b5d60f4c43
--- /dev/null
+++ b/deal.II/deal.II/source/fe/fe_raviart_thomas_3d.cc
@@ -0,0 +1,242 @@
+//----------------------------------------------------------------
+//    $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