/*@{*/
/**
- * Implementation of a scalar Lagrange finite element @p Qp that
- * yields the finite element space of continuous, piecewise polynomials
- * of degree @p p in each coordinate direction. This class is realized
- * using tensor product polynomials based on equidistant support
- * points.
+ * Implementation of a scalar Lagrange finite element @p Qp that yields the
+ * finite element space of continuous, piecewise polynomials of degree @p p in
+ * each coordinate direction. This class is realized using tensor product
+ * polynomials based on equidistant or given support points.
*
- * The constructor of this class takes the degree @p p of this finite
- * element.
+ * The standard constructor of this class takes the degree @p p of this finite
+ * element. Alternatively, it can take a quadrature formula @p points defining
+ * the support points of the Lagrange interpolation in one coordinate direction.
*
* <h3>Implementation</h3>
*
- * The constructor creates a TensorProductPolynomials object
- * that includes the tensor product of @p LagrangeEquidistant
- * polynomials of degree @p p. This @p TensorProductPolynomials
- * object provides all values and derivatives of the shape functions.
+ * The constructor creates a TensorProductPolynomials object that includes the
+ * tensor product of @p LagrangeEquidistant polynomials of degree @p p. This
+ * @p TensorProductPolynomials object provides all values and derivatives of
+ * the shape functions. In case a quadrature rule is given, the constructure
+ * creates a TensorProductPolynomials object that includes the tensor product
+ * of @p Lagrange polynomials with the support points from @p points.
*
* Furthermore the constructor filles the @p interface_constraints,
* the @p prolongation (embedding) and the @p restriction
* @endverbatim
* </ul>
*
- * @author Wolfgang Bangerth, 1998, 2003; Guido Kanschat, 2001; Ralf Hartmann, 2001, 2004, 2005; Oliver Kayser-Herold, 2004
+ * @author Wolfgang Bangerth, 1998, 2003; Guido Kanschat, 2001; Ralf Hartmann, 2001, 2004, 2005; Oliver Kayser-Herold, 2004; Katharina Kormann, 2008; Martin Kronbichler, 2008
*/
template <int dim>
class FE_Q : public FE_Poly<TensorProductPolynomials<dim>,dim>
* polynomials of degree @p p.
*/
FE_Q (const unsigned int p);
-
+
+ /**
+ * Constructor for tensor product
+ * polynomials with support points @p
+ * points based on a one-dimensional
+ * quadrature formula. The degree of the
+ * finite element is
+ * <tt>points.size()-1</tt>. Note that
+ * the first point has to be 0 and the
+ * last one 1.
+ */
+
+ FE_Q (const Quadrature<1> &points);
/**
* Return a string that uniquely
* identifies a finite
* from the constructor.
*/
void initialize_constraints ();
+
+ /**
+ * Initialize the hanging node
+ * constraints matrices. Called from the
+ * constructor in case the finite element
+ * is based on quadrature points.
+ */
+ void initialize_constraints (const Quadrature<1> &points);
/**
* Initialize the embedding
* constructor.
*/
void initialize_unit_support_points ();
+
+ /**
+ * Initialize the @p unit_support_points
+ * field of the FiniteElement
+ * class. Called from the constructor in
+ * case the finite element is based on
+ * quadrature points.
+ */
+ void initialize_unit_support_points (const Quadrature<1> &points);
/**
* Initialize the
* constructor.
*/
void initialize_unit_face_support_points ();
+
+ /**
+ * Initialize the @p
+ * unit_face_support_points field of the
+ * FiniteElement class. Called from the
+ * constructor in case the finite element
+ * is based on quadrature points.
+ */
+ void initialize_unit_face_support_points (const Quadrature<1> &points);
/**
* Initialize the
FE_Q<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int);
template <>
-void FE_Q<1>::initialize_constraints ();
+void FE_Q<1>::initialize_constraints (const Quadrature<1>&);
template <>
-void FE_Q<2>::initialize_constraints ();
+void FE_Q<2>::initialize_constraints (const Quadrature<1>&);
template <>
-void FE_Q<3>::initialize_constraints ();
+void FE_Q<3>::initialize_constraints (const Quadrature<1>&);
DEAL_II_NAMESPACE_CLOSE
#include <base/template_constraints.h>
#include <fe/fe_q.h>
#include <fe/fe_tools.h>
+#include <base/quadrature_lib.h>
#include <vector>
-
+#include <iostream>
#include <sstream>
DEAL_II_NAMESPACE_OPEN
+template <int dim>
+FE_Q<dim>::FE_Q (const Quadrature<1> &points)
+ :
+ FE_Poly<TensorProductPolynomials<dim>, dim> (
+ TensorProductPolynomials<dim>(Polynomials::Lagrange::generate_complete_basis(points.get_points())),
+ FiniteElementData<dim>(get_dpo_vector(points.n_quadrature_points-1),
+ 1, points.n_quadrature_points-1,
+ FiniteElementData<dim>::H1),
+ std::vector<bool> (1, false),
+ std::vector<std::vector<bool> >(1, std::vector<bool>(1,true))),
+ face_index_map(FE_Q_Helper::invert_numbering(face_lexicographic_to_hierarchic_numbering (points.n_quadrature_points-1)))
+{
+ const unsigned int degree = points.n_quadrature_points-1;
+
+ Assert (degree > 0,
+ ExcMessage ("This element can only be used for polynomial degrees "
+ "at least zero"));
+ Assert (points.point(0)(0) == 0,
+ ExcMessage ("The first support point has to be zero."));
+ Assert (points.point(degree)(0) == 1,
+ ExcMessage ("The last support point has to be one."));
+
+ std::vector<unsigned int> renumber (this->dofs_per_cell);
+ FETools::hierarchic_to_lexicographic_numbering (*this, renumber);
+ this->poly_space.set_numbering(renumber);
+
+ // finally fill in support points
+ // on cell and face
+ initialize_unit_support_points (points);
+ initialize_unit_face_support_points (points);
+
+ // compute constraint, embedding
+ // and restriction matrices
+ initialize_constraints (points);
+
+ // Reinit the vectors of
+ // restriction and prolongation
+ // matrices to the right sizes
+ this->reinit_restriction_and_prolongation_matrices();
+
+ // Fill prolongation matrices with
+ // embedding operators
+ FETools::compute_embedding_matrices (*this, this->prolongation);
+
+ // Fill restriction matrices with
+ // L2-projection
+ FETools::compute_projection_matrices (*this, this->restriction);
+ initialize_quad_dof_index_permutation();
+}
+
+
+
template <int dim>
std::string
FE_Q<dim>::get_name () const
// have to be kept in synch
std::ostringstream namebuf;
- namebuf << "FE_Q<" << dim << ">(" << this->degree << ")";
+ bool type = true;
+ const unsigned int n_points = this->degree +1;
+ std::vector<double> points(n_points);
+ const unsigned int dofs_per_cell = this->dofs_per_cell;
+ const std::vector<Point<dim> > &unit_support_points = this->unit_support_points;
+ unsigned int index = 0;
+
+ // Decode the support points
+ // in one coordinate direction.
+ for (unsigned int j=0;j<dofs_per_cell;j++)
+ {
+ if ((dim>1) ? (unit_support_points[j](1)==0 &&
+ ((dim>2) ? unit_support_points[j](2)==0: true)) : true)
+ {
+ if (index == 0)
+ points[index] = unit_support_points[j](0);
+ else if (index == 1)
+ points[n_points-1] = unit_support_points[j](0);
+ else
+ points[index-1] = unit_support_points[j](0);
+ index++;
+ }
+ }
+ Assert (index == n_points,
+ ExcMessage ("Could not decode support points in one coordinate direction."));
+
+ // Check whether the support
+ // points are equidistant.
+ for(unsigned int j=0;j<n_points;j++)
+ if (std::fabs(points[j] - (double)j/this->degree) > 1e-15)
+ {
+ type = false;
+ break;
+ }
+
+ if (type == true)
+ namebuf << "FE_Q<" << dim << ">(" << this->degree << ")";
+ else
+ {
+
+ // Check whether the support
+ // points come from QGaussLobatto.
+ const QGaussLobatto<1> points_gl(n_points);
+ type = true;
+ for(unsigned int j=0;j<n_points;j++)
+ if (points[j] != points_gl.point(j)(0))
+ {
+ type = false;
+ break;
+ }
+ if(type == true)
+ namebuf << "FE_Q<" << dim << ">(QGaussLobatto(" << this->degree+1 << "))";
+ else
+ namebuf << "FE_Q<" << dim << ">(QUnknownNodes(" << this->degree << "))";
+ }
return namebuf.str();
}
void
FE_Q<dim>::
get_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
- FullMatrix<double> &interpolation_matrix) const
+ FullMatrix<double> &interpolation_matrix) const
{
// this is only implemented, if the
// source FE is also a
const FE_Q<dim> &source_fe
= dynamic_cast<const FE_Q<dim>&>(x_source_fe);
- const std::vector<unsigned int> &index_map=
- this->poly_space.get_numbering();
-
// compute the interpolation
// matrices in much the same way as
// we do for the embedding matrices
source_fe.dofs_per_cell);
for (unsigned int j=0; j<this->dofs_per_cell; ++j)
{
- // generate a point on this
+ // read in a point on this
// cell and evaluate the
// shape functions there
- const Point<dim>
- p = FE_Q_Helper::generate_unit_point (index_map[j], this->dofs_per_cell,
- dealii::internal::int2type<dim>());
+ const Point<dim> p = this->unit_support_points[j];
for (unsigned int i=0; i<this->dofs_per_cell; ++i)
cell_interpolation(j,i) = this->poly_space.compute_value (i, p);
}
// then compute the
- // interpolation matrix matrix
+ // interpolation matrix
// for this coordinate
// direction
cell_interpolation.gauss_jordan ();
}
+
template <>
void
FE_Q<1>::
const Point<dim> &p = face_quadrature.point (i);
for (unsigned int j=0; j<this->dofs_per_face; ++j)
- {
+ {
double matrix_entry = this->shape_value (this->face_to_cell_index(j, 0), p);
- // Correct the interpolated
- // value. I.e. if it is close
- // to 1 or 0, make it exactly
- // 1 or 0. Unfortunately, this
- // is required to avoid problems
- // with higher order elements.
+ // Correct the interpolated
+ // value. I.e. if it is close
+ // to 1 or 0, make it exactly
+ // 1 or 0. Unfortunately, this
+ // is required to avoid problems
+ // with higher order elements.
if (fabs (matrix_entry - 1.0) < eps)
matrix_entry = 1.0;
if (fabs (matrix_entry) < eps)
}
+
template <int dim>
void
FE_Q<dim>::
get_subface_interpolation_matrix (const FiniteElement<dim> &x_source_fe,
- const unsigned int subface,
- FullMatrix<double> &interpolation_matrix) const
+ const unsigned int subface,
+ FullMatrix<double> &interpolation_matrix) const
{
// this is only implemented, if the
// source FE is also a
}
+
+template <int dim>
+void FE_Q<dim>::initialize_unit_support_points (const Quadrature<1> &points)
+{
+ // number of points: (degree+1)^dim
+ unsigned int n = this->degree+1;
+ for (unsigned int i=1; i<dim; ++i)
+ n *= this->degree+1;
+
+ this->unit_support_points.resize(n);
+
+ const std::vector<unsigned int> &index_map_inverse=
+ this->poly_space.get_numbering_inverse();
+
+ Quadrature<dim> support_quadrature(points);
+
+ Point<dim> p;
+
+ for (unsigned int k=0;k<n ;k++)
+ {
+ this->unit_support_points[index_map_inverse[k]] = support_quadrature.point(k);
+ }
+}
+
+
#if deal_II_dimension == 1
template <>
// no faces in 1d, so nothing to do
}
+template <>
+void FE_Q<1>::initialize_unit_face_support_points (const Quadrature<1> &/*points*/)
+{
+ // no faces in 1d, so nothing to do
+}
+
#endif
+template <int dim>
+void FE_Q<dim>::initialize_unit_face_support_points (const Quadrature<1> &points)
+{
+ const unsigned int codim = dim-1;
+
+ // number of points: (degree+1)^codim
+ unsigned int n = this->degree+1;
+ for (unsigned int i=1; i<codim; ++i)
+ n *= this->degree+1;
+
+ this->unit_face_support_points.resize(n);
+
+ const std::vector< Point<1> > edge = points.get_points();
+
+ const std::vector<unsigned int> &face_index_map_inverse=
+ FE_Q_Helper::invert_numbering(face_index_map);
+
+ Point<codim> p;
+
+ unsigned int k=0;
+ for (unsigned int iz=0; iz <= ((codim>2) ? this->degree : 0) ; ++iz)
+ for (unsigned int iy=0; iy <= ((codim>1) ? this->degree : 0) ; ++iy)
+ for (unsigned int ix=0; ix<=this->degree; ++ix)
+ {
+ p(0) = edge[ix](0);
+ if (codim>1)
+ p(1) = edge[iy](0);
+ if (codim>2)
+ p(2) = edge[iz](0);
+
+ this->unit_face_support_points[face_index_map_inverse[k++]] = p;
+ }
+}
+
+
+
template <int dim>
void
FE_Q<dim>::initialize_quad_dof_index_permutation ()
}
-
#if deal_II_dimension == 3
template <>
}
+
template <int dim>
std::vector<unsigned int>
FE_Q<dim>::face_lexicographic_to_hierarchic_numbering (const unsigned int degree)
#endif
+
+
+template <int dim>
+void
+FE_Q<dim>::initialize_constraints ()
+{
+ QTrapez<1> trapez;
+ QIterated<1> points (trapez,this->degree);
+ initialize_constraints (points);
+}
+
+
#if deal_II_dimension == 1
template <>
void
-FE_Q<1>::initialize_constraints ()
+FE_Q<1>::initialize_constraints (const Quadrature<1> &/*points*/)
{
// no constraints in 1d
}
#endif
+
#if deal_II_dimension == 2
template <>
void
-FE_Q<2>::initialize_constraints ()
+FE_Q<2>::initialize_constraints (const Quadrature<1> &points)
{
const unsigned int dim = 2;
// destination (child) and source (mother)
// dofs.
const std::vector<Polynomials::Polynomial<double> > polynomials=
- Polynomials::LagrangeEquidistant::generate_complete_basis(this->degree);
+ Polynomials::Lagrange::generate_complete_basis(points.get_points());
this->interface_constraints
.TableBase<2,double>::reinit (this->interface_constraints_size());
this->interface_constraints(i,j) =
polynomials[face_index_map[j]].value (constraint_points[i](0));
- // if the value is small up
- // to round-off, then
- // simply set it to zero to
- // avoid unwanted fill-in
- // of the constraint
- // matrices (which would
- // then increase the number
- // of other DoFs a
- // constrained DoF would
- // couple to)
+ // if the value is small up
+ // to round-off, then
+ // simply set it to zero to
+ // avoid unwanted fill-in
+ // of the constraint
+ // matrices (which would
+ // then increase the number
+ // of other DoFs a
+ // constrained DoF would
+ // couple to)
if (std::fabs(this->interface_constraints(i,j)) < 1e-14)
this->interface_constraints(i,j) = 0;
}
#endif
#if deal_II_dimension == 3
+
template <>
void
-FE_Q<3>::initialize_constraints ()
+FE_Q<3>::initialize_constraints (const Quadrature<1> &points)
{
const unsigned int dim = 3;
// Now construct relation between
// destination (child) and source (mother)
// dofs.
- const unsigned int pnts=(this->degree+1)*(this->degree+1);
+ const unsigned int pnts=(this->degree+1)*(this->degree+1);
const std::vector<Polynomials::Polynomial<double> > polynomial_basis=
- Polynomials::LagrangeEquidistant::generate_complete_basis(this->degree);
+ Polynomials::Lagrange::generate_complete_basis(points.get_points());
+ //const std::vector<Polynomials::Polynomial<double> > polynomial_basis=
+ //Polynomials::LagrangeEquidistant::generate_complete_basis(this->degree);
const TensorProductPolynomials<dim-1> face_polynomials(polynomial_basis);