#define __deal2__fe_raviart_thomas_h
#include <base/config.h>
+#include <base/table.h>
#include <base/polynomials_raviart_thomas.h>
#include <base/polynomial.h>
#include <base/tensor_product_polynomials.h>
* space H<sup>div</sup>. These elements generate vector fields with
* normel components continuous between mesh cells.
*
+ * @todo Right now, the description of node values and interpolation
+ * does not match the actual state of this class. This will be
+ * improved after some discussion.
+ *
* We follow the usual definition of the degree of RT elements, which
* denotes the polynomial degree of the largest complete polynomial
* subspace contained in the RT space. Then, approciamtion order of
virtual bool has_support_on_face (const unsigned int shape_index,
const unsigned int face_index) const;
+ virtual void interpolate(std::vector<double>& local_dofs,
+ const std::vector<double>& values) const;
+ virtual void interpolate(std::vector<double>& local_dofs,
+ const std::vector<Vector<double> >& values,
+ unsigned int offset = 0) const;
+
+ virtual void interpolate(
+ std::vector<double>& local_dofs,
+ const VectorSlice<const std::vector<std::vector<double> > >& values) const;
/**
* Determine an estimate for the
* memory consumption (in bytes)
typename Mapping<dim>::InternalDataBase &fe_internal,
FEValuesData<dim>& data) const;
- virtual void interpolate(std::vector<double>& local_dofs,
- const std::vector<double>& values) const;
- virtual void interpolate(std::vector<double>& local_dofs,
- const std::vector<Vector<double> >& values,
- unsigned int offset = 0) const;
-
- virtual void interpolate(
- std::vector<double>& local_dofs,
- const VectorSlice<const std::vector<std::vector<double> > >& values) const;
private:
/**
* The order of the
void initialize_restriction ();
/**
- * Initialize the
- * @p unit_support_points field
- * of the FiniteElement
- * class. Called from the
- * constructor.
+ * Initialize the @p
+ * generalized_support_points
+ * field of the FiniteElement
+ * class and fill the tables with
+ * interpolation weights
+ * (#boundary_weights and
+ * #interior_weights). Called
+ * from the constructor.
*/
void initialize_support_points (const unsigned int rt_degree);
};
/**
- * The quadrature formula used to
- * generate support points on
- * faces and computing the
- * moments on faces. Its number
- * of points is one order higher
- * than the degree of the RT
- * space.
+ * These are the factors
+ * multiplied to a function in
+ * the
+ * #generalized_face_support_points
+ * when computing the
+ * integration. They are
+ * organized such that there is
+ * one row for each generalized
+ * face support point and one
+ * column for each degree of
+ * freedom on the face.
*/
- QGauss<dim-1> face_quadrature;
-
+ Table<2, double> boundary_weights;
/**
- * Legendre polynomials are used
- * for computing the moments.
+ * Precomputed factors for
+ * interpolation of interior
+ * degrees of freedom. The
+ * rationale for this Table is
+ * the same as for
+ * #boundary_weights. Only, this
+ * table has a third coordinate
+ * for the space direction of the
+ * component evaluated.
*/
- std::vector<Polynomials::Polynomial<double> > legendre_polynomials;
+ Table<3, double> interior_weights;
/**
* Allow access from other
std::vector<bool>(dim,true))),
rt_order(rt_order),
polynomials (create_polynomials(rt_order)),
- renumber (compute_renumber(rt_order)),
- face_quadrature((dim>1) ? rt_order+1 : 0),
- legendre_polynomials(rt_order)
+ renumber (compute_renumber(rt_order))
{
Assert (dim >= 2, ExcImpossibleInDim(dim));
#endif
+#if deal_II_dimension == 1
template <int dim>
void
FE_RaviartThomas<dim>::initialize_support_points (const unsigned int deg)
{
+ return;
+
+ Assert (false, ExcNotImplemented());
+
QGauss<dim> cell_quadrature(deg+1);
const unsigned int n_interior_points
- = (deg>1) ? cell_quadrature.n_quadrature_points : 0;
+ = (deg>0) ? cell_quadrature.n_quadrature_points : 0;
+
+ this->generalized_support_points.resize (2 + n_interior_points);
+
+ // Number of the point being entered
+ unsigned int current = 0;
+
+
+ if (deg==0) return;
+
+ interior_weights.reinit(TableIndices<3>(2+n_interior_points, 0, dim));
+
+ for (unsigned int k=0;k<cell_quadrature.n_quadrature_points;++k)
+ this->generalized_support_points[current++] = cell_quadrature.point(k);
+
+ Assert (current == this->generalized_support_points.size(),
+ ExcInternalError());
+}
+
+#else
+
+// Version for 2d and higher. See above for 1d version
+template <int dim>
+void
+FE_RaviartThomas<dim>::initialize_support_points (const unsigned int deg)
+{
+ QGauss<dim> cell_quadrature(deg+1);
+ const unsigned int n_interior_points
+ = (deg>0) ? cell_quadrature.n_quadrature_points : 0;
unsigned int n_face_points = (dim>1) ? 1 : 0;
// compute (deg+1)^(dim-1)
this->generalized_support_points.resize (GeometryInfo<dim>::faces_per_cell*n_face_points
- +cell_quadrature.n_quadrature_points);
+ + n_interior_points);
this->generalized_face_support_points.resize (n_face_points);
// Number of the point being entered
if (dim>1)
{
QGauss<dim-1> face_points (deg+1);
+ TensorProductPolynomials<dim-1> legendre
+ = Polynomials::Legendre::generate_complete_basis(deg);
+
+ boundary_weights.reinit(n_face_points, legendre.n());
+
Assert (face_points.n_quadrature_points == 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);
+ for (unsigned int k=0;k<n_face_points;++k)
+ {
+ this->generalized_face_support_points[k] = face_points.point(k);
+ // Compute its quadrature
+ // contribution for each
+ // moment.
+ for (unsigned int i=0;i<legendre.n();++i)
+ {
+ boundary_weights(k, i)
+ = face_points.weight(k)
+ * legendre.compute_value(i, 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);
-
- current = this->dofs_per_face*GeometryInfo<dim>::faces_per_cell;
+ for (;current<GeometryInfo<dim>::faces_per_cell*n_face_points;
+ ++current)
+ {
+ // Enter the support point
+ // into the vector
+ this->generalized_support_points[current] = faces.point(current);
+ }
}
if (deg==0) return;
- for (unsigned int k=0;k<cell_quadrature.n_quadrature_points;++k)
- this->generalized_support_points[current++] = cell_quadrature.point(k);
-
- Assert (current == this->generalized_support_points.size(),
- ExcInternalError());
+ // Create Legendre basis for the
+ // space D_xi Q_k
+ std::vector<AnisotropicPolynomials<dim>* > polynomials(dim);
+ for (unsigned int dd=0;dd<dim;++dd)
+ {
+ std::vector<std::vector<Polynomials::Polynomial<double> > > poly(dim);
+ for (unsigned int d=0;d<dim;++d)
+ poly[d] = Polynomials::Legendre::generate_complete_basis(deg);
+ poly[dd] = Polynomials::Legendre::generate_complete_basis(deg-1);
-//TODO: Unit support points should go away according to new definition.
+ polynomials[dd] = new AnisotropicPolynomials<dim>(poly);
+ }
+
+ interior_weights.reinit(TableIndices<3>(n_interior_points, polynomials[0]->n(), dim));
- this->unit_support_points.resize (this->dofs_per_cell);
- switch (dim)
+ for (unsigned int k=0;k<cell_quadrature.n_quadrature_points;++k)
{
- case 2:
- {
- Assert (rt_order+1 == this->dofs_per_face, ExcInternalError());
-
- // associate support points
- // with mid-face points if a
- // shape function has a
- // non-zero normal component
- // there, otherwise with the
- // cell center. the reason
- // for this non-unique
- // support point is that we
- // use hierarchical shape
- // functions, rather than
- // Lagrange functions, for
- // which we get into the same
- // trouble as in the
- // FE_Q_Hierarchical element;
- // see the respective
- // function there
-
- // start with the face shape
- // functions
- for (unsigned int i=0; i<this->dofs_per_face; ++i)
- this->unit_support_points[0*this->dofs_per_face+i] = Point<dim>(.5, .0);
- for (unsigned int i=0; i<this->dofs_per_face; ++i)
- this->unit_support_points[1*this->dofs_per_face+i] = Point<dim>(1., .5);
- for (unsigned int i=0; i<this->dofs_per_face; ++i)
- this->unit_support_points[2*this->dofs_per_face+i] = Point<dim>(.5, 1.);
- for (unsigned int i=0; i<this->dofs_per_face; ++i)
- this->unit_support_points[3*this->dofs_per_face+i] = Point<dim>(.0, .5);
-
- // associate the rest with
- // the cell center
- for (unsigned int i=4*this->dofs_per_face; i<this->dofs_per_cell; ++i)
- this->unit_support_points[i] = Point<dim>(.5, .5);
-
- break;
- }
-
- case 3:
- {
- // same as in 2d
- Assert ((rt_order+1)*(rt_order+1) == this->dofs_per_face, ExcInternalError());
-
- // start with the face shape
- // functions
- for (unsigned int i=0; i<this->dofs_per_face; ++i)
- this->unit_support_points[0*this->dofs_per_face+i] = Point<dim>(.5, .0, .5);
- for (unsigned int i=0; i<this->dofs_per_face; ++i)
- this->unit_support_points[1*this->dofs_per_face+i] = Point<dim>(.5, 1., .5);
- for (unsigned int i=0; i<this->dofs_per_face; ++i)
- this->unit_support_points[2*this->dofs_per_face+i] = Point<dim>(.5, .5, 0.);
- for (unsigned int i=0; i<this->dofs_per_face; ++i)
- this->unit_support_points[3*this->dofs_per_face+i] = Point<dim>(1., .5, .5);
- for (unsigned int i=0; i<this->dofs_per_face; ++i)
- this->unit_support_points[4*this->dofs_per_face+i] = Point<dim>(.5, .5, 1.);
- for (unsigned int i=0; i<this->dofs_per_face; ++i)
- this->unit_support_points[5*this->dofs_per_face+i] = Point<dim>(.0, .5, .5);
-
- // associate the rest with
- // the cell center
- for (unsigned int i=6*this->dofs_per_face; i<this->dofs_per_cell; ++i)
- this->unit_support_points[i] = Point<dim>(.5, .5, .5);
-
- break;
- }
+ this->generalized_support_points[current++] = cell_quadrature.point(k);
+ for (unsigned int i=0;i<polynomials[0]->n();++i)
+ for (unsigned int d=0;d<dim;++d)
+ interior_weights(k,i,d) = cell_quadrature.weight(k)
+ * polynomials[d]->compute_value(i,cell_quadrature.point(k));
+ }
- default:
- Assert (false, ExcNotImplemented());
- };
+ for (unsigned int d=0;d<dim;++d)
+ delete polynomials[d];
+
+ Assert (current == this->generalized_support_points.size(),
+ ExcInternalError());
}
+#endif
#if deal_II_dimension == 1
ExcDimensionMismatch(values[0].size(),offset+this->n_components()));
std::fill(local_dofs.begin(), local_dofs.end(), 0.);
+
+ const unsigned int n_face_points = boundary_weights.size(0);
+ for (unsigned int face=0;face<GeometryInfo<dim>::faces_per_cell;++face)
+ for (unsigned int k=0;k<n_face_points;++k)
+ for (unsigned int i=0;i<boundary_weights.size(1);++i)
+ {
+ local_dofs[i] += boundary_weights(k,i)
+ * values[face*n_face_points+k](offset+GeometryInfo<dim>::unit_normal_direction[face]);
+ }
}
template <int dim>
void
FE_RaviartThomas<dim>::interpolate(
- std::vector<double>& /*local_dofs*/,
- const VectorSlice<const std::vector<std::vector<double> > >& /*values*/) const
-{}
+ std::vector<double>& local_dofs,
+ const VectorSlice<const std::vector<std::vector<double> > >& values) const
+{
+ Assert (values.size() == this->n_components(),
+ ExcDimensionMismatch(values.size(), this->n_components()));
+ Assert (values[0].size() == this->generalized_support_points.size(),
+ ExcDimensionMismatch(values[0].size(), this->generalized_support_points.size()));
+ Assert (local_dofs.size() == this->dofs_per_cell,
+ ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
+
+ std::fill(local_dofs.begin(), local_dofs.end(), 0.);
+
+ const unsigned int n_face_points = boundary_weights.size(0);
+ for (unsigned int face=0;face<GeometryInfo<dim>::faces_per_cell;++face)
+ for (unsigned int k=0;k<n_face_points;++k)
+ for (unsigned int i=0;i<boundary_weights.size(1);++i)
+ {
+ local_dofs[i] += boundary_weights(k,i)
+ * values[GeometryInfo<dim>::unit_normal_direction[face]][face*n_face_points+k];
+ }
+}
template <int dim>