From 6beb171a918da3707ab839cb6397925fa2bdcc64 Mon Sep 17 00:00:00 2001 From: kormann Date: Sun, 7 Sep 2008 19:44:38 +0000 Subject: [PATCH] Added a second constructor to FE_Q that allows to generate it based on support points from a one-dimensional quadrature formula. Added get_name info for Gauss-Lobatto based finite elements even to FE_DGQ. Gauss-Lobatto based FE_Qs and FE_DGQs can also be recognized by FE_Tools::get_fe_from_name but not constructed. git-svn-id: https://svn.dealii.org/trunk@16756 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/source/quadrature_lib.cc | 13 ++ deal.II/deal.II/include/fe/fe_q.h | 72 +++++-- deal.II/deal.II/source/fe/fe_dgq.cc | 55 +++++- deal.II/deal.II/source/fe/fe_q.cc | 261 ++++++++++++++++++++++---- deal.II/deal.II/source/fe/fe_tools.cc | 33 +++- deal.II/doc/authors.html | 3 + deal.II/doc/news/changes.h | 8 + 7 files changed, 385 insertions(+), 60 deletions(-) diff --git a/deal.II/base/source/quadrature_lib.cc b/deal.II/base/source/quadrature_lib.cc index 05d3d2ec36..72fd248008 100644 --- a/deal.II/base/source/quadrature_lib.cc +++ b/deal.II/base/source/quadrature_lib.cc @@ -50,6 +50,19 @@ QGauss<0>::QGauss (const unsigned int) +template <> +QGaussLobatto<0>::QGaussLobatto (const unsigned int) + : + Quadrature<0> (0) +{ + // this function has to be provided to + // avoid certain linker failures, but it + // should never be called + Assert (false, ExcInternalError()); +} + + + template <> QGauss<1>::QGauss (const unsigned int n) : diff --git a/deal.II/deal.II/include/fe/fe_q.h b/deal.II/deal.II/include/fe/fe_q.h index 0975f0dbed..35c147ff66 100644 --- a/deal.II/deal.II/include/fe/fe_q.h +++ b/deal.II/deal.II/include/fe/fe_q.h @@ -24,21 +24,23 @@ DEAL_II_NAMESPACE_OPEN /*@{*/ /** - * 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. * *

Implementation

* - * 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 @@ -232,7 +234,7 @@ DEAL_II_NAMESPACE_OPEN * @endverbatim * * - * @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 class FE_Q : public FE_Poly,dim> @@ -243,7 +245,19 @@ class FE_Q : public FE_Poly,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 + * points.size()-1. 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 @@ -494,6 +508,14 @@ class FE_Q : public FE_Poly,dim> * 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 @@ -517,6 +539,15 @@ class FE_Q : public FE_Poly,dim> * 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 @@ -526,6 +557,15 @@ class FE_Q : public FE_Poly,dim> * 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 @@ -571,13 +611,13 @@ std::vector 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 diff --git a/deal.II/deal.II/source/fe/fe_dgq.cc b/deal.II/deal.II/source/fe/fe_dgq.cc index 4b06ab4786..e0f9bfc951 100644 --- a/deal.II/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/deal.II/source/fe/fe_dgq.cc @@ -12,6 +12,7 @@ //--------------------------------------------------------------------------- #include +#include #include #include #include @@ -231,7 +232,6 @@ FE_DGQ::get_name () const std::ostringstream namebuf; namebuf << "FE_DGQ<" << dim << ">(" << this->degree << ")"; - return namebuf.str(); } @@ -666,8 +666,57 @@ FE_DGQArbitraryNodes::get_name () const // FE_DGQArbitraryNodes since // there is no initialization by // a degree value. - std::ostringstream namebuf; - namebuf << "FE_DGQArbitraryNodes<" << dim << ">(" << this->degree << ")"; + std::ostringstream namebuf; + + bool type = true; + const unsigned int n_points = this->degree +1; + std::vector points(n_points); + const unsigned int dofs_per_cell = this->dofs_per_cell; + const std::vector > &unit_support_points = this->unit_support_points; + unsigned int index = 0; + + // Decode the support points + // in one coordinate direction. + for (unsigned int j=0;j1) ? (unit_support_points[j](1)==0 && + ((dim>2) ? unit_support_points[j](2)==0: true)) : true) + { + points[index++] = unit_support_points[j](0); + } + } + 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;jdegree) > 1e-15) + { + type = false; + break; + } + + if (type == true) + namebuf << "FE_DGQ<" << 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(QGaussLobatto(" << this->degree+1 << "))"; + else + namebuf << "FE_DGQArbitraryNodes<" << dim << ">(QUnknownNodes(" << this->degree << "))"; + } return namebuf.str(); } diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 29aa0d8eca..5377a8fb47 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -16,9 +16,10 @@ #include #include #include +#include #include - +#include #include DEAL_II_NAMESPACE_OPEN @@ -209,6 +210,58 @@ FE_Q::FE_Q (const unsigned int degree) +template +FE_Q::FE_Q (const Quadrature<1> &points) + : + FE_Poly, dim> ( + TensorProductPolynomials(Polynomials::Lagrange::generate_complete_basis(points.get_points())), + FiniteElementData(get_dpo_vector(points.n_quadrature_points-1), + 1, points.n_quadrature_points-1, + FiniteElementData::H1), + std::vector (1, false), + std::vector >(1, std::vector(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 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 std::string FE_Q::get_name () const @@ -221,8 +274,62 @@ FE_Q::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 points(n_points); + const unsigned int dofs_per_cell = this->dofs_per_cell; + const std::vector > &unit_support_points = this->unit_support_points; + unsigned int index = 0; + + // Decode the support points + // in one coordinate direction. + for (unsigned int j=0;j1) ? (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;jdegree) > 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(QGaussLobatto(" << this->degree+1 << "))"; + else + namebuf << "FE_Q<" << dim << ">(QUnknownNodes(" << this->degree << "))"; + } return namebuf.str(); } @@ -241,7 +348,7 @@ template void FE_Q:: get_interpolation_matrix (const FiniteElement &x_source_fe, - FullMatrix &interpolation_matrix) const + FullMatrix &interpolation_matrix) const { // this is only implemented, if the // source FE is also a @@ -264,9 +371,6 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, const FE_Q &source_fe = dynamic_cast&>(x_source_fe); - const std::vector &index_map= - this->poly_space.get_numbering(); - // compute the interpolation // matrices in much the same way as // we do for the embedding matrices @@ -279,12 +383,10 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, source_fe.dofs_per_cell); for (unsigned int j=0; jdofs_per_cell; ++j) { - // generate a point on this + // read in a point on this // cell and evaluate the // shape functions there - const Point - p = FE_Q_Helper::generate_unit_point (index_map[j], this->dofs_per_cell, - dealii::internal::int2type()); + const Point p = this->unit_support_points[j]; for (unsigned int i=0; idofs_per_cell; ++i) cell_interpolation(j,i) = this->poly_space.compute_value (i, p); @@ -293,7 +395,7 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, } // then compute the - // interpolation matrix matrix + // interpolation matrix // for this coordinate // direction cell_interpolation.gauss_jordan (); @@ -337,6 +439,7 @@ get_face_interpolation_matrix (const FiniteElement<1> &/*x_source_fe*/, } + template <> void FE_Q<1>:: @@ -424,15 +527,15 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, const Point &p = face_quadrature.point (i); for (unsigned int j=0; jdofs_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) @@ -460,12 +563,13 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, } + template void FE_Q:: get_subface_interpolation_matrix (const FiniteElement &x_source_fe, - const unsigned int subface, - FullMatrix &interpolation_matrix) const + const unsigned int subface, + FullMatrix &interpolation_matrix) const { // this is only implemented, if the // source FE is also a @@ -747,6 +851,31 @@ void FE_Q::initialize_unit_support_points () } + +template +void FE_Q::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; idegree+1; + + this->unit_support_points.resize(n); + + const std::vector &index_map_inverse= + this->poly_space.get_numbering_inverse(); + + Quadrature support_quadrature(points); + + Point p; + + for (unsigned int k=0;kunit_support_points[index_map_inverse[k]] = support_quadrature.point(k); + } +} + + #if deal_II_dimension == 1 template <> @@ -755,6 +884,12 @@ void FE_Q<1>::initialize_unit_face_support_points () // 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 @@ -793,6 +928,42 @@ void FE_Q::initialize_unit_face_support_points () +template +void FE_Q::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; idegree+1; + + this->unit_face_support_points.resize(n); + + const std::vector< Point<1> > edge = points.get_points(); + + const std::vector &face_index_map_inverse= + FE_Q_Helper::invert_numbering(face_index_map); + + Point 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 void FE_Q::initialize_quad_dof_index_permutation () @@ -801,7 +972,6 @@ FE_Q::initialize_quad_dof_index_permutation () } - #if deal_II_dimension == 3 template <> @@ -885,6 +1055,7 @@ FE_Q::get_dpo_vector(const unsigned int deg) } + template std::vector FE_Q::face_lexicographic_to_hierarchic_numbering (const unsigned int degree) @@ -907,22 +1078,35 @@ FE_Q<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int) #endif + + +template +void +FE_Q::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; @@ -1022,7 +1206,7 @@ FE_Q<2>::initialize_constraints () // destination (child) and source (mother) // dofs. const std::vector > 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()); @@ -1033,16 +1217,16 @@ FE_Q<2>::initialize_constraints () 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; } @@ -1051,9 +1235,10 @@ FE_Q<2>::initialize_constraints () #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; @@ -1148,9 +1333,11 @@ FE_Q<3>::initialize_constraints () // 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 > polynomial_basis= - Polynomials::LagrangeEquidistant::generate_complete_basis(this->degree); + Polynomials::Lagrange::generate_complete_basis(points.get_points()); + //const std::vector > polynomial_basis= + //Polynomials::LagrangeEquidistant::generate_complete_basis(this->degree); const TensorProductPolynomials face_polynomials(polynomial_basis); diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index 83f3e21908..42dd5e1872 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -92,6 +92,8 @@ namespace = FEFactoryPointer(new FETools::FEFactory >); default_map["FE_DGQ"] = FEFactoryPointer(new FETools::FEFactory >); + default_map["FE_DGQArbitraryNodes"] + = FEFactoryPointer(new FETools::FEFactory >); default_map["FE_Q"] = FEFactoryPointer(new FETools::FEFactory >); @@ -1751,14 +1753,37 @@ FETools::get_fe_from_name_aux (std::string &name) AssertThrow (fe_name_map.find(name_part) != fe_name_map.end(), ExcInvalidFEName(name)); // Now, just the (degree) + // or (Quadrature<1>(degree+1)) // part should be left. if (name.size() == 0 || name[0] != '(') throw (std::string("Invalid first character in ") + name); name.erase(0,1); - const std::pair tmp - = Utilities::get_integer_at_position (name, 0); - name.erase(0, tmp.second+1); - return fe_name_map.find(name_part)->second->get(tmp.first); + if (name[0] != 'Q') + { + const std::pair tmp + = Utilities::get_integer_at_position (name, 0); + name.erase(0, tmp.second+1); + return fe_name_map.find(name_part)->second->get(tmp.first); + } + else + { + unsigned int position = name.find('('); + const std::string quadrature_name(name, 0, position-1); + name.erase(0,position); + if (quadrature_name.compare("QGaussLobatto") == 0) + { + const std::pair tmp + = Utilities::get_integer_at_position (name, 0); + name.erase(0, tmp.second+1); +//TODO: Implement a get function taking Quadrature<1> in fe_tools.h. + //return fe_name_map.find(name_part)->second->get(QGaussLobatto<1>(tmp.first)); + AssertThrow (false, ExcNotImplemented()); + } + else + { + AssertThrow (false,ExcNotImplemented()); + } + } } diff --git a/deal.II/doc/authors.html b/deal.II/doc/authors.html index 72d326ba7c..6d2d659e29 100644 --- a/deal.II/doc/authors.html +++ b/deal.II/doc/authors.html @@ -95,6 +95,9 @@ alphabetical order):
  • Benjamin Shelton Kirk: Tecplot output. +
  • Katharina Kormann: + Support for arbitrary nodes in FE_Q. +
  • Martin Kronbichler: Parts of step-22 and step-31 tutorial programs, interfaces to Trilinos (partly), and a few enhancements. diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 18222d8231..3936007689 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -278,6 +278,14 @@ inconvenience this causes.

    deal.II

      +
    1. +

      + New: The class FE_Q can now alternatively be constructed based on + support points from a given one-dimensional quadrature rule. +
      + (Katharina Kormann, Martin Kronbichler, 2008/09/07) +

      +
    2. Fixed: Using the ConstraintMatrix class, when a degree of freedom was -- 2.39.5