]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Remove again from main branch until element is actually tested a little.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 7 Apr 2003 15:47:33 +0000 (15:47 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 7 Apr 2003 15:47:33 +0000 (15:47 +0000)
git-svn-id: https://svn.dealii.org/trunk@7371 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_raviart_thomas.h [deleted file]
deal.II/deal.II/source/fe/fe_raviart_thomas.cc [deleted file]
deal.II/deal.II/source/fe/fe_raviart_thomas_1d.cc [deleted file]
deal.II/deal.II/source/fe/fe_raviart_thomas_2d.cc [deleted file]
deal.II/deal.II/source/fe/fe_raviart_thomas_3d.cc [deleted file]

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

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.