+++ /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