]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add initial versions of these files.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 7 Apr 2003 15:43:52 +0000 (15:43 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 7 Apr 2003 15:43:52 +0000 (15:43 +0000)
git-svn-id: https://svn.dealii.org/trunk@7370 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_raviart_thomas.h [new file with mode: 0644]
deal.II/deal.II/source/fe/fe_raviart_thomas.cc [new file with mode: 0644]
deal.II/deal.II/source/fe/fe_raviart_thomas_1d.cc [new file with mode: 0644]
deal.II/deal.II/source/fe/fe_raviart_thomas_2d.cc [new file with mode: 0644]
deal.II/deal.II/source/fe/fe_raviart_thomas_3d.cc [new file with mode: 0644]

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 (file)
index 0000000..73ee24c
--- /dev/null
@@ -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 (file)
index 0000000..64e44f5
--- /dev/null
@@ -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 (file)
index 0000000..a47a21a
--- /dev/null
@@ -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 (file)
index 0000000..70e9484
--- /dev/null
@@ -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 (file)
index 0000000..b5d60f4
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.