--- /dev/null
+//---------------------------------------------------------------
+// $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
--- /dev/null
+//----------------------------------------------------------------
+// $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>;
--- /dev/null
+//----------------------------------------------------------------
+// $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
+
--- /dev/null
+//----------------------------------------------------------------
+// $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
--- /dev/null
+//----------------------------------------------------------------
+// $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