public:
/**
* Constructor.
+ *
+ * @arg @c degree: constructor
+ * argument for poly. May be
+ * different from @p
+ * fe_data.degree.
*/
FE_PolyTensor (unsigned int degree,
const FiniteElementData<dim> &fe_data,
*/
virtual unsigned int element_multiplicity (const unsigned int index) const;
+ /**
+ * Given <tt>flags</tt>,
+ * determines the values which
+ * must be computed only for the
+ * reference cell. Make sure,
+ * that #mapping_type is set by
+ * the derived class, such that
+ * this function can operate
+ * correctly.
+ */
+ virtual UpdateFlags update_once (const UpdateFlags flags) const;
+ /**
+ * Given <tt>flags</tt>,
+ * determines the values which
+ * must be computed in each cell
+ * cell. Make sure, that
+ * #mapping_type is set by the
+ * derived class, such that this
+ * function can operate
+ * correctly.
+ */
+ virtual UpdateFlags update_each (const UpdateFlags flags) const;
protected:
/**
* possible to avoid the mapping.
*/
enum MappingType {
- /// Shape functions do not depend on actual mesh cell
+ /**
+ * No mapping has been
+ * selected, throw an error
+ * if needed.
+ */
+ no_mapping,
+ /**
+ * Shape functions do not
+ * depend on actual mesh
+ * cell
+ */
independent,
- /// Shape functions are transformed covariant.
+ /**
+ * Shape functions do not
+ * depend on actual mesh
+ * cell. The mapping class
+ * must be
+ * MappingCartesian.
+ */
+ independent_on_cartesian,
+ /**
+ * Shape functions are
+ * transformed covariant.
+ */
covariant,
- /// Shape functions are transformed contravariant.
+ /**
+ * Shape functions are
+ * transformed
+ * contravariant.
+ */
contravariant
};
std::vector<std::vector<Tensor<2,dim> > > shape_grads;
};
- /**
- * Degree of the polynomials.
- */
- unsigned int degree;
-
/**
* The polynomial space. Its type
* is given by the template
* parameter POLY.
*/
POLY poly_space;
+
/**
* The inverse of the matrix
* <i>a<sub>ij</sub></i> of node
#define __deal2__fe_raviart_thomas_h
#include <base/config.h>
+#include <base/polynomials_raviart_thomas.h>
#include <base/tensor_product_polynomials.h>
#include <grid/geometry_info.h>
#include <fe/fe.h>
+#include <fe/fe_poly_tensor.h>
template <int dim> class MappingQ;
/**
* Check whether a shape function
- * is non-zero on a face.
+ * may be non-zero on a face.
*
* Right now, this is only
* implemented for RT0 in
* 1D. Otherwise, returns always
* @p true.
- *
- * Implementation of the
- * interface in
- * FiniteElement
*/
virtual bool has_support_on_face (const unsigned int shape_index,
const unsigned int face_index) const;
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
- * FiniteElement.
- */
virtual void
fill_fe_values (const Mapping<dim> &mapping,
const typename Triangulation<dim>::cell_iterator &cell,
typename Mapping<dim>::InternalDataBase &fe_internal,
FEValuesData<dim>& data) const;
- /**
- * Implementation of the same
- * function in
- * FiniteElement.
- */
virtual void
fill_fe_face_values (const Mapping<dim> &mapping,
const typename Triangulation<dim>::cell_iterator &cell,
typename Mapping<dim>::InternalDataBase &fe_internal,
FEValuesData<dim>& data) const;
- /**
- * Implementation of the same
- * function in
- * FiniteElement.
- */
virtual void
fill_fe_subface_values (const Mapping<dim> &mapping,
const typename Triangulation<dim>::cell_iterator &cell,
template <int dim1> friend class FE_RaviartThomas;
};
+
+
+/**
+ * The Raviart-Thomas elements with node functionals defined as point
+ * values in Gauss points.
+ *
+ * <h3>Description of node values</h3>
+ *
+ * For this Raviart-Thomas element, the node values are not cell and
+ * face moments with respect to certain polynomials, but the values in
+ * quadrature points.
+ *
+ * For an RT-element of degree <i>k</i>, we choose
+ * <i>k+1<sup>d-1</sup></i> Gauss points on each face. This way, the
+ * normal component which is in <i>Q<sub>k</sub></i> is uniquely
+ * determined. Furthermore, since this Gauss-formula is exact on
+ * <i>Q<sub>2k+1</sub></i>, these node values correspond to the exact
+ * integration of the moments of the RT-space.
+ *
+ * In the interior of the cells, the moments are with respect to an
+ * anisotropic <i>Q<sub>k</sub></i> space, where the test functions
+ * are one degree lower in the direction corresponding to the vector
+ * component under consideration. This can be emulated by using an
+ * anisotropic Gauss formula for integration.
+ *
+ * @warning The degree stored in the member variable
+ * FiniteElementData<dim>::degree is higher by one than the
+ * constructor argument!
+ *
+ * @author Guido Kanschat, 2005
+ */
+template <int dim>
+class FE_RaviartThomasNodal
+ :
+ public FE_PolyTensor<PolynomialsRaviartThomas<dim>, dim>
+{
+ /**
+ * Constructor for the Raviart-Thomas
+ * element of degree @p p.
+ */
+ FE_RaviartThomasNodal (const unsigned int p);
+
+ /**
+ * Return a string that uniquely
+ * identifies a finite
+ * element. This class returns
+ * <tt>FE_RaviartThomasNodal<dim>(degree)</tt>, with
+ * @p dim and @p degree
+ * replaced by appropriate
+ * values.
+ */
+ virtual std::string get_name () const;
+
+ virtual FiniteElement<dim>* clone () const;
+
+ /**
+ * Check whether a shape function
+ * may be non-zero on a face.
+ *
+ * Right now, always returns
+ * @p true.
+ */
+ virtual bool has_support_on_face (const unsigned int shape_index,
+ const unsigned int face_index) 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);
+
+ /**
+ * Compute the vector used for
+ * the
+ * @p restriction_is_additive
+ * field passed to the base
+ * class's constructor.
+ */
+ static std::vector<bool>
+ get_ria_vector (const unsigned int degree);
+ /**
+ * Initialize the
+ * FiniteElementBase<dim>::unit_support_points
+ * and FiniteElementBase<dim>::unit_face_support_points
+ * fields. Called from the
+ * constructor.
+ */
+ void initialize_unit_support_points (const unsigned int degree);
+
+ /**
+ * Initialize the
+ * #inverse_node_matrix
+ * field. Called from the
+ * constructor.
+ */
+ void initialize_node_matrix ();
+};
+
+
/*@}*/
/* -------------- declaration of explicit specializations ------------- */
FiniteElement<dim> (fe_data,
restriction_is_additive_flags,
nonzero_components),
- degree(degree),
poly_space(POLY(degree))
{
cached_point(0) = -1;
template <class POLY, int dim>
void
-FE_PolyTensor<POLY,dim>::fill_fe_values (const Mapping<dim> &mapping,
- const typename Triangulation<dim>::cell_iterator &cell,
- const Quadrature<dim> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename Mapping<dim>::InternalDataBase &fedata,
- FEValuesData<dim> &data) const
+FE_PolyTensor<POLY,dim>::fill_fe_values (
+ const Mapping<dim> &mapping,
+ const typename Triangulation<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
const unsigned int n_quad = quadrature.n_quadrature_points;
const UpdateFlags flags(fe_data.current_update_flags());
- Assert(mapping_type == independent, ExcNotImplemented());
+ Assert(mapping_type == independent || mapping_type == independent_on_cartesian,
+ ExcNotImplemented());
Assert(!(flags & update_values) || fe_data.shape_values.size() == n_dofs,
ExcDimensionMismatch(fe_data.shape_values.size(), n_dofs));
template <class POLY, int dim>
void
-FE_PolyTensor<POLY,dim>::fill_fe_face_values (const Mapping<dim> &mapping,
- const typename Triangulation<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
+FE_PolyTensor<POLY,dim>::fill_fe_face_values (
+ const Mapping<dim> &mapping,
+ const typename Triangulation<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
const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
- Assert(mapping_type == independent, ExcNotImplemented());
+ Assert(mapping_type == independent || mapping_type == independent_on_cartesian,
+ ExcNotImplemented());
//TODO: Size assertions
for (unsigned int i=0; i<n_dofs; ++i)
template <class POLY, int dim>
void
-FE_PolyTensor<POLY,dim>::fill_fe_subface_values (const Mapping<dim> &mapping,
- const typename Triangulation<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
+FE_PolyTensor<POLY,dim>::fill_fe_subface_values (
+ const Mapping<dim> &mapping,
+ const typename Triangulation<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
const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
- Assert(mapping_type == independent, ExcNotImplemented());
+ Assert(mapping_type == independent || mapping_type == independent_on_cartesian,
+ ExcNotImplemented());
//TODO: Size assertions
for (unsigned int i=0; i<n_dofs; ++i)
}
+template <class POLY, int dim>
+UpdateFlags
+FE_PolyTensor<POLY,dim>::update_once (const UpdateFlags flags) const
+{
+ Assert (mapping_type != no_mapping, ExcNotInitialized());
+ const bool values_once = (mapping_type == independent);
+
+ UpdateFlags out = update_default;
+
+ if (values_once && (flags & update_values))
+ out |= update_values;
+
+ return out;
+}
+
+
+template <class POLY, int dim>
+UpdateFlags
+FE_PolyTensor<POLY,dim>::update_each (const UpdateFlags flags) const
+{
+ Assert (mapping_type != no_mapping, ExcNotInitialized());
+ const bool values_once = (mapping_type == independent);
+
+ UpdateFlags out = update_default;
+
+ if (!values_once && (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;
+}
+
template class FE_PolyTensor<PolynomialsRaviartThomas<deal_II_dimension>,deal_II_dimension>;
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2003, 2004, 2005 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_lib.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>
+#include <fe/fe_tools.h>
+
+#ifdef HAVE_STD_STRINGSTREAM
+# include <sstream>
+#else
+# include <strstream>
+#endif
+
+
+template <int dim>
+FE_RaviartThomasNodal<dim>::FE_RaviartThomasNodal (const unsigned int deg)
+ :
+ FE_PolyTensor<PolynomialsRaviartThomas<dim>, dim> (
+ deg,
+ FiniteElementData<dim>(get_dpo_vector(deg),
+ dim, deg+1),
+ get_ria_vector (deg),
+ std::vector<std::vector<bool> >(
+ FiniteElementData<dim>(get_dpo_vector(deg),
+ dim,deg+1).dofs_per_cell,
+ std::vector<bool>(dim,true)))
+{
+ Assert (dim >= 2, ExcImpossibleInDim(dim));
+
+ this->mapping_type = independent_on_cartesian;
+
+ for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i)
+ this->prolongation[i].reinit (this->dofs_per_cell,
+ this->dofs_per_cell);
+ FETools::compute_embedding_matrices (*this, &this->prolongation[0]);
+ // finally fill in support points
+ // on cell and face
+ initialize_unit_support_points (deg);
+ initialize_node_matrix();
+}
+
+
+
+template <int dim>
+std::string
+FE_RaviartThomasNodal<dim>::get_name () const
+{
+ // note that the
+ // FETools::get_fe_from_name
+ // function depends on the
+ // particular format of the string
+ // this function returns, so they
+ // have to be kept in synch
+
+#ifdef HAVE_STD_STRINGSTREAM
+ std::ostringstream namebuf;
+#else
+ std::ostrstream namebuf;
+#endif
+
+ namebuf << "FE_RaviartThomasNodal<" << dim << ">(" << degree-1 << ")";
+
+#ifndef HAVE_STD_STRINGSTREAM
+ namebuf << std::ends;
+#endif
+ return namebuf.str();
+}
+
+
+
+template <int dim>
+bool
+FE_RaviartThomasNodal<dim>::has_support_on_face (unsigned int, unsigned int) const
+{
+ return true;
+}
+
+
+
+template <int dim>
+FiniteElement<dim> *
+FE_RaviartThomasNodal<dim>::clone() const
+{
+ return new FE_RaviartThomasNodal<dim>(degree-1);
+}
+
+
+#if deal_II_dimension == 1
+
+template <>
+std::vector<unsigned int>
+FE_RaviartThomasNodal<1>::get_dpo_vector (const unsigned int deg)
+{
+ std::vector<unsigned int> dpo(2);
+ dpo[0] = 1;
+ dpo[1] = deg;
+ return dpo;
+}
+
+#endif
+
+
+template <int dim>
+std::vector<unsigned int>
+FE_RaviartThomasNodal<dim>::get_dpo_vector (const unsigned int deg)
+{
+ // the element is face-based and we have
+ // (deg+1)^(dim-1) DoFs per face
+ unsigned int dofs_per_face = 1;
+ for (unsigned int d=0; d<dim-1; ++d)
+ dofs_per_face *= deg+1;
+
+ // and then there are interior dofs
+ const unsigned int
+ interior_dofs = dim*deg*dofs_per_face;
+
+ std::vector<unsigned int> dpo(dim+1);
+ dpo[dim-1] = dofs_per_face;
+ dpo[dim] = interior_dofs;
+
+ return dpo;
+}
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+std::vector<bool>
+FE_RaviartThomasNodal<1>::get_ria_vector (const unsigned int)
+{
+ Assert (false, ExcImpossibleInDim(1));
+ return std::vector<bool>();
+}
+
+#endif
+
+
+template <int dim>
+std::vector<bool>
+FE_RaviartThomasNodal<dim>::get_ria_vector (const unsigned int deg)
+{
+ unsigned int dofs_per_cell, dofs_per_face;
+ switch (dim)
+ {
+ case 2:
+ dofs_per_face = deg+1;
+ dofs_per_cell = 2*(deg+1)*(deg+2);
+ break;
+ case 3:
+ dofs_per_face = (deg+1)*(deg+1);
+ dofs_per_cell = 3*(deg+1)*(deg+1)*(deg+2);
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ Assert (FiniteElementData<dim>(get_dpo_vector(deg),dim).dofs_per_cell ==
+ dofs_per_cell,
+ ExcInternalError());
+ Assert (FiniteElementData<dim>(get_dpo_vector(deg),dim).dofs_per_face ==
+ dofs_per_face,
+ ExcInternalError());
+
+ // all face dofs need to be
+ // non-additive, since they have
+ // continuity requirements.
+ // however, the interior dofs are
+ // made additive
+ std::vector<bool> ret_val(dofs_per_cell,false);
+ for (unsigned int i=GeometryInfo<dim>::faces_per_cell*dofs_per_face;
+ i < dofs_per_cell; ++i)
+ ret_val[i] = true;
+
+ return ret_val;
+}
+
+
+template <int dim>
+void
+FE_RaviartThomasNodal<dim>::initialize_unit_support_points (const unsigned int deg)
+{
+ this->unit_support_points.resize (this->dofs_per_cell);
+ this->unit_face_support_points.resize (this->dofs_per_face);
+
+ unsigned int current = 0;
+
+ // On the faces, we choose as many
+ // Gauss points as necessary to
+ // determine the normal component
+ // uniquely. This is the deg of
+ // the Raviart-Thomas element plus
+ // one.
+ if (dim>1)
+ {
+ QGauss<dim-1> face_points (deg+1);
+ Assert (face_points.n_quadrature_points == this->dofs_per_face,
+ ExcInternalError());
+ for (unsigned int k=0;k<this->dofs_per_face;++k)
+ this->unit_face_support_points[k] = face_points.point(k);
+ Quadrature<dim> faces = QProjector<dim>::project_to_all_faces(face_points);
+ for (unsigned int k=0;k<faces.n_quadrature_points;++k)
+ this->unit_support_points[k] = faces.point(k);
+
+ current = faces.n_quadrature_points;
+ }
+ // In the interior, we need
+ // anisotropic Gauss quadratures,
+ // different for each direction.
+ QGauss<1> high(deg+1);
+ QGauss<1> low(deg);
+
+ for (unsigned int d=0;d<dim;++d)
+ {
+ QAnisotropic<dim>* quadrature;
+ if (dim == 1) quadrature = new QAnisotropic<dim>(high);
+ if (dim == 2) quadrature = new QAnisotropic<dim>(((d==0) ? low : high),
+ ((d==1) ? low : high));
+ if (dim == 3) quadrature = new QAnisotropic<dim>(((d==0) ? low : high),
+ ((d==1) ? low : high),
+ ((d==2) ? low : high));
+ for (unsigned int k=0;k<quadrature->n_quadrature_points;++k)
+ this->unit_support_points[current++] = quadrature->point(k);
+ delete quadrature;
+ }
+ Assert (current == this->dofs_per_cell, ExcInternalError());
+}
+
+
+template <int dim>
+void
+FE_RaviartThomasNodal<dim>::initialize_node_matrix ()
+{
+ const unsigned int n_dofs = this->dofs_per_cell;
+
+ // We use an auxiliary matrix in
+ // this function. Therefore,
+ // inverse_node_matrix is still
+ // empty and shape_value_component
+ // returns the 'raw' shape values.
+ FullMatrix<double> N(n_dofs, n_dofs);
+
+ // The curent node functional index
+ unsigned int current = 0;
+
+ // For each face and all quadrature
+ // points on it, the node value is
+ // the normal component of the
+ // shape function, possibly
+ // pointing in negative direction.
+ for (unsigned int face = 0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ for (unsigned int k=0;k<this->dofs_per_face;++k)
+ {
+ for (unsigned int i=0;i<n_dofs;++i)
+ N(current,i) = this->shape_value_component(
+ i, unit_support_point(current),
+ GeometryInfo< dim >::unit_normal_direction[face])
+ * GeometryInfo< dim >::unit_normal_orientation[face];
+ ++current;
+ }
+ // Interior degrees of freedom in each direction
+ const unsigned int n_cell = (n_dofs - current) / dim;
+
+ for (unsigned int d=0;d<dim;++d)
+ for (unsigned int k=0;k<n_cell;++k)
+ {
+ for (unsigned int i=0;i<n_dofs;++i)
+ N(current,i) = this->shape_value_component(i, unit_support_point(current), d);
+ ++current;
+ }
+ Assert (current == n_dofs, ExcInternalError());
+
+ inverse_node_matrix.invert(N);
+}
+
+
+
+template FE_RaviartThomasNodal<deal_II_dimension>;