virtual std::string get_name () 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, <tt>base_element(0)</tt> 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,
- * <tt>element_multiplicity(0)</tt>
- * returns one, and all other
- * indices will throw an error.
- */
- virtual unsigned int element_multiplicity (const unsigned int index) const;
-
/**
* Check whether a shape function
* may be non-zero on a face.
virtual FiniteElement<dim> * clone() const;
private:
- /**
- * The order of the
- * Raviart-Thomas element. The
- * lowest order elements are
- * usually referred to as RT0,
- * even though their shape
- * functions are piecewise
- * linears.
- */
- const unsigned int rt_order;
-
/**
* Only for internal use. Its
* full name is
*/
void initialize_restriction ();
- /**
- * 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;
-
/**
* Fields of cell-independent data.
*
virtual FiniteElementDomination::Domination
compare_for_face_domination (const FiniteElement<dim> &fe_other) const;
- virtual UpdateFlags update_once (const UpdateFlags flags) const;
- virtual UpdateFlags update_each (const UpdateFlags flags) const;
private:
/**
* Only for internal use. Its
UpdateFlags
FE_PolyTensor<POLY,dim,spacedim>::update_each (const UpdateFlags flags) const
{
- const bool values_once = (mapping_type == mapping_none);
-
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_hessians)
- out |= update_hessians | update_covariant_transformation;
- if (flags & update_cell_normal_vectors)
- out |= update_cell_normal_vectors | update_JxW_values;
-
+ switch (mapping_type)
+ {
+ case mapping_piola:
+ {
+ if (flags & update_values)
+ out |= update_values | update_piola;
+
+ if (flags & update_gradients)
+ out |= update_gradients | update_piola | update_covariant_transformation;
+
+ if (flags & update_hessians)
+ out |= update_hessians | update_piola | update_covariant_transformation;
+
+ break;
+ }
+
+ default:
+ {
+ Assert (false, ExcNotImplemented());
+ }
+ }
+
return out;
}
DEAL_II_NAMESPACE_OPEN
+
template <int dim>
FE_RaviartThomas<dim>::FE_RaviartThomas (const unsigned int deg)
:
dim, deg+1, FiniteElementData<dim>::Hdiv, 1),
std::vector<bool>(PolynomialsRaviartThomas<dim>::compute_n_pols(deg), true),
std::vector<std::vector<bool> >(PolynomialsRaviartThomas<dim>::compute_n_pols(deg),
- std::vector<bool>(dim,true))),
- rt_order(deg)
+ std::vector<bool>(dim,true)))
{
Assert (dim >= 2, ExcImpossibleInDim(dim));
const unsigned int n_dofs = this->dofs_per_cell;
//matrix, generating the correct
//basis functions from the raw
//ones.
+
+ // 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> M(n_dofs, n_dofs);
FETools::compute_node_matrix(M, *this);
this->inverse_node_matrix.reinit(n_dofs, n_dofs);
FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face];
for (unsigned int i=0; i<GeometryInfo<dim>::max_children_per_face; ++i)
face_embeddings[i].reinit (this->dofs_per_face, this->dofs_per_face);
- FETools::compute_face_embedding_matrices<dim,double>(*this,
- face_embeddings,
- 0, 0);
+ FETools::compute_face_embedding_matrices<dim,double>(*this, face_embeddings, 0, 0);
this->interface_constraints.reinit((1<<(dim-1)) * this->dofs_per_face,
this->dofs_per_face);
unsigned int target_row=0;
// this function returns, so they
// have to be kept in synch
+ // note that this->degree is the maximal
+ // polynomial degree and is thus one higher
+ // than the argument given to the
+ // constructor
std::ostringstream namebuf;
- namebuf << "FE_RaviartThomas<" << dim << ">(" << rt_order << ")";
+ namebuf << "FE_RaviartThomas<" << dim << ">(" << this->degree-1 << ")";
return namebuf.str();
}
-
template <int dim>
FiniteElement<dim> *
FE_RaviartThomas<dim>::clone() const
{
const unsigned int iso=RefinementCase<dim>::isotropic_refinement-1;
- QGauss<dim-1> q_base (rt_order+1);
+ QGauss<dim-1> q_base (this->degree);
const unsigned int n_face_points = q_base.size();
// First, compute interpolation on
// subfaces
}
}
- if (rt_order==0) return;
+ if (this->degree == 1) return;
// Create Legendre basis for the
// space D_xi Q_k. Here, we cannot
{
std::vector<std::vector<Polynomials::Polynomial<double> > > poly(dim);
for (unsigned int d=0;d<dim;++d)
- poly[d] = Polynomials::Legendre::generate_complete_basis(rt_order);
- poly[dd] = Polynomials::Legendre::generate_complete_basis(rt_order-1);
+ poly[d] = Polynomials::Legendre::generate_complete_basis(this->degree-1);
+ poly[dd] = Polynomials::Legendre::generate_complete_basis(this->degree-2);
polynomials[dd] = new AnisotropicPolynomials<dim>(poly);
}
- QGauss<dim> q_cell(rt_order+1);
+ QGauss<dim> q_cell(this->degree);
const unsigned int start_cell_dofs
= GeometryInfo<dim>::faces_per_cell*this->dofs_per_face;
-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_piola;
-
- if (flags & update_gradients)
- out |= update_gradients | update_piola | update_covariant_transformation;
-
- if (flags & update_hessians)
- out |= update_hessians | update_piola | update_covariant_transformation;
-
- return out;
-}
-
//---------------------------------------------------------------------------
// Data field initialization
//---------------------------------------------------------------------------
-
-
-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,
// Return computed values if we
// know them easily. Otherwise, it
// is always safe to return true.
- switch (rt_order)
+ switch (this->degree)
{
- case 0:
+ case 1:
{
switch (dim)
{
#include <sstream>
+
DEAL_II_NAMESPACE_OPEN
const unsigned int n_dofs = this->dofs_per_cell;
this->mapping_type = mapping_raviart_thomas;
- // These must be done first, since
- // they change the evaluation of
- // basis functions
-
- // Set up the generalized support
- // points
- initialize_support_points (deg);
- //Now compute the inverse node
+ // First, initialize the
+ // generalized support points and
+ // quadrature weights, since they
+ // are required for interpolation.
+ initialize_support_points(deg);
+ // Now compute the inverse node
//matrix, generating the correct
//basis functions from the raw
//ones.
}
+//---------------------------------------------------------------------------
+// Auxiliary and internal functions
+//---------------------------------------------------------------------------
+
+
+
+template <int dim>
+void
+FE_RaviartThomasNodal<dim>::initialize_support_points (const unsigned int deg)
+{
+ this->generalized_support_points.resize (this->dofs_per_cell);
+ this->generalized_face_support_points.resize (this->dofs_per_face);
+
+ // Number of the point being entered
+ 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.size() == this->dofs_per_face,
+ ExcInternalError());
+ for (unsigned int k=0;k<this->dofs_per_face;++k)
+ this->generalized_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<this->dofs_per_face*GeometryInfo<dim>::faces_per_cell;
+ ++k)
+ this->generalized_support_points[k] = faces.point(k+QProjector<dim>
+ ::DataSetDescriptor::face(0,
+ true,
+ false,
+ false,
+ this->dofs_per_face));
+
+ current = this->dofs_per_face*GeometryInfo<dim>::faces_per_cell;
+ }
+
+ if (deg==0) return;
+ // 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));
+ Assert(dim<=3, ExcNotImplemented());
+
+ for (unsigned int k=0;k<quadrature->size();++k)
+ this->generalized_support_points[current++] = quadrature->point(k);
+ delete quadrature;
+ }
+ Assert (current == this->dofs_per_cell, ExcInternalError());
+}
+
+
+#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=1; d<dim; ++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)
+{
+ const unsigned int dofs_per_cell = PolynomialsRaviartThomas<dim>::compute_n_pols(deg);
+ unsigned int dofs_per_face = deg+1;
+ for (unsigned int d=2;d<dim;++d)
+ dofs_per_face *= deg+1;
+ // 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
-#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=1; d<dim; ++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)
-{
- const unsigned int dofs_per_cell = PolynomialsRaviartThomas<dim>::compute_n_pols(deg);
- unsigned int dofs_per_face = deg+1;
- for (unsigned int d=2;d<dim;++d)
- dofs_per_face *= deg+1;
- // 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_support_points (const unsigned int deg)
-{
- this->generalized_support_points.resize (this->dofs_per_cell);
- this->generalized_face_support_points.resize (this->dofs_per_face);
-
- // Number of the point being entered
- 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.size() == this->dofs_per_face,
- ExcInternalError());
- for (unsigned int k=0;k<this->dofs_per_face;++k)
- this->generalized_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<this->dofs_per_face*GeometryInfo<dim>::faces_per_cell;
- ++k)
- this->generalized_support_points[k] = faces.point(k+QProjector<dim>
- ::DataSetDescriptor::face(0,
- true,
- false,
- false,
- this->dofs_per_face));
-
- current = this->dofs_per_face*GeometryInfo<dim>::faces_per_cell;
- }
-
- if (deg==0) return;
- // 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));
- Assert(dim<=3, ExcNotImplemented());
-
- for (unsigned int k=0;k<quadrature->size();++k)
- this->generalized_support_points[current++] = quadrature->point(k);
- delete quadrature;
- }
- Assert (current == this->dofs_per_cell, ExcInternalError());
-}
-
-
-template <int dim>
-UpdateFlags
-FE_RaviartThomasNodal<dim>::update_once (const UpdateFlags) const
-{
- return update_default;
-}
-
-
-template <int dim>
-UpdateFlags
-FE_RaviartThomasNodal<dim>::update_each (const UpdateFlags flags) const
-{
- UpdateFlags out = update_default;
-
- if (flags & update_values)
- out |= update_values | update_piola;
-
- if (flags & update_gradients)
- out |= update_gradients | update_piola | update_covariant_transformation;
-
- if (flags & update_hessians)
- out |= update_hessians | update_piola | update_covariant_transformation;
-
- return out;
-}
-
-
template class FE_RaviartThomasNodal<deal_II_dimension>;
DEAL_II_NAMESPACE_CLOSE