From: grahambenharper Date: Sat, 7 Sep 2019 23:49:59 +0000 (-0600) Subject: Remove PolynomialType template argument from FE_Poly X-Git-Tag: v9.3.0-rc1~1594^2~12 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3b3e9b2b4d1f609718b99d449ce8793df087ed97;p=dealii.git Remove PolynomialType template argument from FE_Poly --- diff --git a/include/deal.II/base/polynomial_space.h b/include/deal.II/base/polynomial_space.h index ed1a0a5cfd..4e050b6225 100644 --- a/include/deal.II/base/polynomial_space.h +++ b/include/deal.II/base/polynomial_space.h @@ -157,7 +157,7 @@ public: * Consider using evaluate() instead. */ double - compute_value(const unsigned int i, const Point &p) const; + compute_value(const unsigned int i, const Point &p) const override; /** * Compute the orderth derivative of the ith polynomial @@ -178,7 +178,7 @@ public: * Consider using evaluate() instead. */ Tensor<1, dim> - compute_grad(const unsigned int i, const Point &p) const; + compute_grad(const unsigned int i, const Point &p) const override; /** * Compute the second derivative (grad_grad) of the ith polynomial @@ -187,7 +187,7 @@ public: * Consider using evaluate() instead. */ Tensor<2, dim> - compute_grad_grad(const unsigned int i, const Point &p) const; + compute_grad_grad(const unsigned int i, const Point &p) const override; /** * Return the number of polynomials spanning the space represented by this diff --git a/include/deal.II/base/polynomials_adini.h b/include/deal.II/base/polynomials_adini.h index 4b9c0b6a06..64daffad89 100644 --- a/include/deal.II/base/polynomials_adini.h +++ b/include/deal.II/base/polynomials_adini.h @@ -78,7 +78,7 @@ public: * Consider using evaluate() instead. */ double - compute_value(const unsigned int i, const Point &p) const; + compute_value(const unsigned int i, const Point &p) const override; /** * Compute the gradient of the ith polynomial at @@ -87,7 +87,7 @@ public: * Consider using evaluate() instead. */ Tensor<1, dim> - compute_grad(const unsigned int i, const Point &p) const; + compute_grad(const unsigned int i, const Point &p) const override; /** * Compute the second derivative (grad_grad) of the ith polynomial @@ -96,7 +96,7 @@ public: * Consider using evaluate() instead. */ Tensor<2, dim> - compute_grad_grad(const unsigned int i, const Point &p) const; + compute_grad_grad(const unsigned int i, const Point &p) const override; /** * Return the name of the space, which is PolynomialsAdini. diff --git a/include/deal.II/base/polynomials_rannacher_turek.h b/include/deal.II/base/polynomials_rannacher_turek.h index 0b56598ba1..4fb378458b 100644 --- a/include/deal.II/base/polynomials_rannacher_turek.h +++ b/include/deal.II/base/polynomials_rannacher_turek.h @@ -58,7 +58,7 @@ public: * Value of basis function @p i at @p p. */ double - compute_value(const unsigned int i, const Point &p) const; + compute_value(const unsigned int i, const Point &p) const override; /** * order-th of basis function @p i at @p p. @@ -73,13 +73,13 @@ public: * Gradient of basis function @p i at @p p. */ Tensor<1, dim> - compute_grad(const unsigned int i, const Point &p) const; + compute_grad(const unsigned int i, const Point &p) const override; /** * Gradient of gradient of basis function @p i at @p p. */ Tensor<2, dim> - compute_grad_grad(const unsigned int i, const Point &p) const; + compute_grad_grad(const unsigned int i, const Point &p) const override; /** * Compute values and derivatives of all basis functions at @p unit_point. diff --git a/include/deal.II/base/scalar_polynomials_base.h b/include/deal.II/base/scalar_polynomials_base.h index 5b192a885c..dc702bbc06 100644 --- a/include/deal.II/base/scalar_polynomials_base.h +++ b/include/deal.II/base/scalar_polynomials_base.h @@ -106,6 +106,64 @@ public: std::vector> &third_derivatives, std::vector> &fourth_derivatives) const = 0; + /** + * Compute the value of the ith polynomial at unit point + * p. + * + * Consider using compute() instead. + */ + virtual double + compute_value(const unsigned int /*i*/, const Point & /*p*/) const + { + Assert(false, ExcNotImplemented()); + return 0; + } + + /** + * Compute the orderth derivative of the ith polynomial + * at unit point p. + * + * Consider using compute() instead. + * + * @tparam order The order of the derivative. + */ + template + Tensor + compute_derivative(const unsigned int /*i*/, const Point & /*p*/) const + { + Assert(false, ExcNotImplemented()); + Tensor empty; + return empty; + } + + /** + * Compute the gradient of the ith polynomial at unit point + * p. + * + * Consider using compute() instead. + */ + virtual Tensor<1, dim> + compute_grad(const unsigned int /*i*/, const Point & /*p*/) const + { + Assert(false, ExcNotImplemented()); + Tensor<1, dim> empty; + return empty; + } + + /** + * Compute the second derivative (grad_grad) of the ith polynomial + * at unit point p. + * + * Consider using compute() instead. + */ + virtual Tensor<2, dim> + compute_grad_grad(const unsigned int /*i*/, const Point & /*p*/) const + { + Assert(false, ExcNotImplemented()); + Tensor<2, dim> empty; + return empty; + } + /** * Return the number of polynomials. */ diff --git a/include/deal.II/base/tensor_product_polynomials.h b/include/deal.II/base/tensor_product_polynomials.h index c6bec0062d..e024553d69 100644 --- a/include/deal.II/base/tensor_product_polynomials.h +++ b/include/deal.II/base/tensor_product_polynomials.h @@ -147,7 +147,7 @@ public: * polynomials all at once and in a much more efficient way. */ double - compute_value(const unsigned int i, const Point &p) const; + compute_value(const unsigned int i, const Point &p) const override; /** * Compute the orderth derivative of the ith tensor @@ -180,7 +180,7 @@ public: * polynomials all at once and in a much more efficient way. */ Tensor<1, dim> - compute_grad(const unsigned int i, const Point &p) const; + compute_grad(const unsigned int i, const Point &p) const override; /** * Compute the second derivative (grad_grad) of the ith tensor @@ -195,7 +195,7 @@ public: * polynomials all at once and in a much more efficient way. */ Tensor<2, dim> - compute_grad_grad(const unsigned int i, const Point &p) const; + compute_grad_grad(const unsigned int i, const Point &p) const override; /** * Return the name of the space, which is TensorProductPolynomials. @@ -340,7 +340,7 @@ public: * tensor polynomials all at once and in a much more efficient way. */ double - compute_value(const unsigned int i, const Point &p) const; + compute_value(const unsigned int i, const Point &p) const override; /** * Compute the orderth derivative of the ith tensor @@ -373,7 +373,7 @@ public: * tensor polynomials all at once and in a much more efficient way. */ Tensor<1, dim> - compute_grad(const unsigned int i, const Point &p) const; + compute_grad(const unsigned int i, const Point &p) const override; /** * Compute the second derivative (grad_grad) of the ith tensor @@ -388,7 +388,7 @@ public: * all tensor polynomials all at once and in a much more efficient way. */ Tensor<2, dim> - compute_grad_grad(const unsigned int i, const Point &p) const; + compute_grad_grad(const unsigned int i, const Point &p) const override; /** * Return the name of the space, which is AnisotropicPolynomials. diff --git a/include/deal.II/base/tensor_product_polynomials_bubbles.h b/include/deal.II/base/tensor_product_polynomials_bubbles.h index 1bfcaaf4d4..0e8d28823f 100644 --- a/include/deal.II/base/tensor_product_polynomials_bubbles.h +++ b/include/deal.II/base/tensor_product_polynomials_bubbles.h @@ -120,7 +120,7 @@ public: * polynomials all at once and in a much more efficient way. */ double - compute_value(const unsigned int i, const Point &p) const; + compute_value(const unsigned int i, const Point &p) const override; /** * Compute the order @p order derivative of the ith tensor product @@ -151,7 +151,7 @@ public: * polynomials all at once and in a much more efficient way. */ Tensor<1, dim> - compute_grad(const unsigned int i, const Point &p) const; + compute_grad(const unsigned int i, const Point &p) const override; /** * Compute the second derivative (grad_grad) of the ith tensor @@ -166,7 +166,7 @@ public: * polynomials all at once and in a much more efficient way. */ Tensor<2, dim> - compute_grad_grad(const unsigned int i, const Point &p) const; + compute_grad_grad(const unsigned int i, const Point &p) const override; /** * Return the number of tensor product polynomials plus the bubble diff --git a/include/deal.II/base/tensor_product_polynomials_const.h b/include/deal.II/base/tensor_product_polynomials_const.h index dc069adef4..0c75e2f59c 100644 --- a/include/deal.II/base/tensor_product_polynomials_const.h +++ b/include/deal.II/base/tensor_product_polynomials_const.h @@ -121,7 +121,7 @@ public: * polynomials all at once and in a much more efficient way. */ double - compute_value(const unsigned int i, const Point &p) const; + compute_value(const unsigned int i, const Point &p) const override; /** * Compute the orderth derivative of the ith tensor @@ -154,7 +154,7 @@ public: * polynomials all at once and in a much more efficient way. */ Tensor<1, dim> - compute_grad(const unsigned int i, const Point &p) const; + compute_grad(const unsigned int i, const Point &p) const override; /** * Compute the second derivative (grad_grad) of the ith tensor @@ -169,7 +169,7 @@ public: * polynomials all at once and in a much more efficient way. */ Tensor<2, dim> - compute_grad_grad(const unsigned int i, const Point &p) const; + compute_grad_grad(const unsigned int i, const Point &p) const override; /** * Return the number of tensor product polynomials plus the constant diff --git a/include/deal.II/fe/fe_dgp.h b/include/deal.II/fe/fe_dgp.h index 355252111a..dcfbac266c 100644 --- a/include/deal.II/fe/fe_dgp.h +++ b/include/deal.II/fe/fe_dgp.h @@ -308,7 +308,7 @@ DEAL_II_NAMESPACE_OPEN * @author Guido Kanschat, 2001, 2002, Ralf Hartmann 2004 */ template -class FE_DGP : public FE_Poly, dim, spacedim> +class FE_DGP : public FE_Poly { public: /** diff --git a/include/deal.II/fe/fe_dgp_monomial.h b/include/deal.II/fe/fe_dgp_monomial.h index 16dffa8ed8..fca49e7da9 100644 --- a/include/deal.II/fe/fe_dgp_monomial.h +++ b/include/deal.II/fe/fe_dgp_monomial.h @@ -283,7 +283,7 @@ DEAL_II_NAMESPACE_OPEN * @author Ralf Hartmann, 2004 */ template -class FE_DGPMonomial : public FE_Poly, dim> +class FE_DGPMonomial : public FE_Poly { public: /** diff --git a/include/deal.II/fe/fe_dgq.h b/include/deal.II/fe/fe_dgq.h index 26705d2c3c..2c9a697e44 100644 --- a/include/deal.II/fe/fe_dgq.h +++ b/include/deal.II/fe/fe_dgq.h @@ -109,7 +109,7 @@ class Quadrature; * @author Ralf Hartmann, Guido Kanschat 2001, 2004 */ template -class FE_DGQ : public FE_Poly, dim, spacedim> +class FE_DGQ : public FE_Poly { public: /** diff --git a/include/deal.II/fe/fe_poly.h b/include/deal.II/fe/fe_poly.h index 77efe6a4a5..cfa174d72d 100644 --- a/include/deal.II/fe/fe_poly.h +++ b/include/deal.II/fe/fe_poly.h @@ -20,6 +20,7 @@ #include #include +#include #include #include @@ -70,20 +71,23 @@ DEAL_II_NAMESPACE_OPEN * @author Ralf Hartmann 2004, Guido Kanschat, 2009 */ -template +template class FE_Poly : public FiniteElement { public: /** * Constructor. */ - FE_Poly(const PolynomialType & poly_space, + FE_Poly(const ScalarPolynomialsBase &poly_space, const FiniteElementData & fe_data, const std::vector & restriction_is_additive_flags, const std::vector &nonzero_components); + /** + * Copy constructor. + */ + FE_Poly(const FE_Poly &fe); + /** * Return the polynomial degree of this finite element, i.e. the value * passed to the constructor. @@ -311,12 +315,12 @@ protected: update_3rd_derivatives)) for (unsigned int i = 0; i < n_q_points; ++i) { - poly_space.evaluate(quadrature.point(i), - values, - grads, - grad_grads, - third_derivatives, - fourth_derivatives); + poly_space->evaluate(quadrature.point(i), + values, + grads, + grad_grads, + third_derivatives, + fourth_derivatives); // the values of shape functions at quadrature points don't change. // consequently, write these values right into the output array if @@ -510,10 +514,9 @@ protected: /** - * The polynomial space. Its type is given by the template parameter - * PolynomialType. + * The polynomial space. */ - PolynomialType poly_space; + const std::unique_ptr> poly_space; }; /*@}*/ diff --git a/include/deal.II/fe/fe_poly.templates.h b/include/deal.II/fe/fe_poly.templates.h index 3cea1551b1..4285aeae0a 100644 --- a/include/deal.II/fe/fe_poly.templates.h +++ b/include/deal.II/fe/fe_poly.templates.h @@ -32,42 +32,40 @@ DEAL_II_NAMESPACE_OPEN -template -FE_Poly::FE_Poly( - const PolynomialType & poly_space, +template +FE_Poly::FE_Poly( + const ScalarPolynomialsBase &poly_space, const FiniteElementData & fe_data, const std::vector & restriction_is_additive_flags, const std::vector &nonzero_components) : FiniteElement(fe_data, restriction_is_additive_flags, nonzero_components) - , poly_space(poly_space) -{ - AssertDimension(dim, PolynomialType::dimension); -} + , poly_space(poly_space.clone()) +{} -template +template unsigned int -FE_Poly::get_degree() const +FE_Poly::get_degree() const { return this->degree; } -template +template double -FE_Poly::shape_value(const unsigned int i, - const Point & p) const +FE_Poly::shape_value(const unsigned int i, + const Point & p) const { AssertIndexRange(i, this->dofs_per_cell); - return poly_space.compute_value(i, p); + return poly_space->compute_value(i, p); } -template +template double -FE_Poly::shape_value_component( +FE_Poly::shape_value_component( const unsigned int i, const Point & p, const unsigned int component) const @@ -75,52 +73,50 @@ FE_Poly::shape_value_component( (void)component; AssertIndexRange(i, this->dofs_per_cell); AssertIndexRange(component, 1); - return poly_space.compute_value(i, p); + return poly_space->compute_value(i, p); } -template +template Tensor<1, dim> -FE_Poly::shape_grad(const unsigned int i, - const Point & p) const +FE_Poly::shape_grad(const unsigned int i, + const Point & p) const { AssertIndexRange(i, this->dofs_per_cell); - return poly_space.template compute_derivative<1>(i, p); + return poly_space->template compute_derivative<1>(i, p); } -template +template Tensor<1, dim> -FE_Poly::shape_grad_component( - const unsigned int i, - const Point & p, - const unsigned int component) const +FE_Poly::shape_grad_component(const unsigned int i, + const Point & p, + const unsigned int component) const { (void)component; AssertIndexRange(i, this->dofs_per_cell); AssertIndexRange(component, 1); - return poly_space.template compute_derivative<1>(i, p); + return poly_space->template compute_derivative<1>(i, p); } -template +template Tensor<2, dim> -FE_Poly::shape_grad_grad( - const unsigned int i, - const Point & p) const +FE_Poly::shape_grad_grad(const unsigned int i, + const Point & p) const { AssertIndexRange(i, this->dofs_per_cell); - return poly_space.template compute_derivative<2>(i, p); + return poly_space->template compute_derivative<2>(i, p); } -template +template Tensor<2, dim> -FE_Poly::shape_grad_grad_component( +FE_Poly::shape_grad_grad_component( const unsigned int i, const Point & p, const unsigned int component) const @@ -128,26 +124,25 @@ FE_Poly::shape_grad_grad_component( (void)component; AssertIndexRange(i, this->dofs_per_cell); AssertIndexRange(component, 1); - return poly_space.template compute_derivative<2>(i, p); + return poly_space->template compute_derivative<2>(i, p); } -template +template Tensor<3, dim> -FE_Poly::shape_3rd_derivative( - const unsigned int i, - const Point & p) const +FE_Poly::shape_3rd_derivative(const unsigned int i, + const Point & p) const { AssertIndexRange(i, this->dofs_per_cell); - return poly_space.template compute_derivative<3>(i, p); + return poly_space->template compute_derivative<3>(i, p); } -template +template Tensor<3, dim> -FE_Poly::shape_3rd_derivative_component( +FE_Poly::shape_3rd_derivative_component( const unsigned int i, const Point & p, const unsigned int component) const @@ -155,26 +150,25 @@ FE_Poly::shape_3rd_derivative_component( (void)component; AssertIndexRange(i, this->dofs_per_cell); AssertIndexRange(component, 1); - return poly_space.template compute_derivative<3>(i, p); + return poly_space->template compute_derivative<3>(i, p); } -template +template Tensor<4, dim> -FE_Poly::shape_4th_derivative( - const unsigned int i, - const Point & p) const +FE_Poly::shape_4th_derivative(const unsigned int i, + const Point & p) const { AssertIndexRange(i, this->dofs_per_cell); - return poly_space.template compute_derivative<4>(i, p); + return poly_space->template compute_derivative<4>(i, p); } -template +template Tensor<4, dim> -FE_Poly::shape_4th_derivative_component( +FE_Poly::shape_4th_derivative_component( const unsigned int i, const Point & p, const unsigned int component) const @@ -182,7 +176,7 @@ FE_Poly::shape_4th_derivative_component( (void)component; AssertIndexRange(i, this->dofs_per_cell); AssertIndexRange(component, 1); - return poly_space.template compute_derivative<4>(i, p); + return poly_space->template compute_derivative<4>(i, p); } @@ -192,10 +186,9 @@ FE_Poly::shape_4th_derivative_component( //--------------------------------------------------------------------------- -template +template UpdateFlags -FE_Poly::requires_update_flags( - const UpdateFlags flags) const +FE_Poly::requires_update_flags(const UpdateFlags flags) const { UpdateFlags out = update_default; @@ -264,9 +257,9 @@ higher_derivatives_need_correcting( -template +template void -FE_Poly::fill_fe_values( +FE_Poly::fill_fe_values( const typename Triangulation::cell_iterator &, const CellSimilarity::Similarity cell_similarity, const Quadrature & quadrature, @@ -336,9 +329,9 @@ FE_Poly::fill_fe_values( -template +template void -FE_Poly::fill_fe_face_values( +FE_Poly::fill_fe_face_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const Quadrature & quadrature, @@ -427,9 +420,9 @@ FE_Poly::fill_fe_face_values( -template +template void -FE_Poly::fill_fe_subface_values( +FE_Poly::fill_fe_subface_values( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, @@ -521,9 +514,9 @@ FE_Poly::fill_fe_subface_values( -template +template inline void -FE_Poly::correct_hessians( +FE_Poly::correct_hessians( internal::FEValuesImplementation::FiniteElementRelatedData &output_data, const internal::FEValuesImplementation::MappingRelatedData @@ -540,9 +533,9 @@ FE_Poly::correct_hessians( -template +template inline void -FE_Poly::correct_third_derivatives( +FE_Poly::correct_third_derivatives( internal::FEValuesImplementation::FiniteElementRelatedData &output_data, const internal::FEValuesImplementation::MappingRelatedData @@ -569,82 +562,62 @@ FE_Poly::correct_third_derivatives( } } -namespace internal -{ - template - inline std::vector - get_poly_space_numbering(const PolynomialType &) - { - Assert(false, ExcNotImplemented()); - return std::vector(); - } - - template - inline std::vector - get_poly_space_numbering_inverse(const PolynomialType &) - { - Assert(false, ExcNotImplemented()); - return std::vector(); - } - - template - inline std::vector - get_poly_space_numbering( - const TensorProductPolynomials &poly) - { - return poly.get_numbering(); - } - - template - inline std::vector - get_poly_space_numbering_inverse( - const TensorProductPolynomials &poly) - { - return poly.get_numbering_inverse(); - } - - template - inline std::vector - get_poly_space_numbering(const TensorProductPolynomialsConst &poly) - { - return poly.get_numbering(); - } - - template - inline std::vector - get_poly_space_numbering_inverse( - const TensorProductPolynomialsConst &poly) - { - return poly.get_numbering_inverse(); - } -} // namespace internal - - - -template + + +template std::vector -FE_Poly::get_poly_space_numbering() const +FE_Poly::get_poly_space_numbering() const { - return internal::get_poly_space_numbering(poly_space); + auto *const space_tensor_prod = + dynamic_cast> *>( + this->poly_space.get()); + + if (space_tensor_prod != nullptr) + return space_tensor_prod->get_numbering(); + + auto *const space_tensor_prod_const = + dynamic_cast> *>( + this->poly_space.get()); + + if (space_tensor_prod_const != nullptr) + return space_tensor_prod_const->get_numbering(); + + Assert(false, ExcNotImplemented()); + return std::vector(); } -template +template std::vector -FE_Poly::get_poly_space_numbering_inverse() const +FE_Poly::get_poly_space_numbering_inverse() const { - return internal::get_poly_space_numbering_inverse(poly_space); + auto *const space_tensor_prod = + dynamic_cast> *>( + this->poly_space.get()); + + if (space_tensor_prod != nullptr) + return space_tensor_prod->get_numbering_inverse(); + + auto *const space_tensor_prod_const = + dynamic_cast> *>( + this->poly_space.get()); + + if (space_tensor_prod_const != nullptr) + return space_tensor_prod_const->get_numbering_inverse(); + + Assert(false, ExcNotImplemented()); + return std::vector(); } -template +template std::size_t -FE_Poly::memory_consumption() const +FE_Poly::memory_consumption() const { return FiniteElement::memory_consumption() + - poly_space.memory_consumption(); + poly_space->memory_consumption(); } diff --git a/include/deal.II/fe/fe_q_base.h b/include/deal.II/fe/fe_q_base.h index 5828c6e970..998c928f56 100644 --- a/include/deal.II/fe/fe_q_base.h +++ b/include/deal.II/fe/fe_q_base.h @@ -41,7 +41,7 @@ DEAL_II_NAMESPACE_OPEN template -class FE_Q_Base : public FE_Poly +class FE_Q_Base : public FE_Poly { public: /** diff --git a/include/deal.II/fe/fe_q_hierarchical.h b/include/deal.II/fe/fe_q_hierarchical.h index 250e98bb59..738ce3aa5d 100644 --- a/include/deal.II/fe/fe_q_hierarchical.h +++ b/include/deal.II/fe/fe_q_hierarchical.h @@ -540,7 +540,7 @@ DEAL_II_NAMESPACE_OPEN * @author Brian Carnes, 2002, Ralf Hartmann 2004, 2005, Denis Davydov, 2015 */ template -class FE_Q_Hierarchical : public FE_Poly, dim> +class FE_Q_Hierarchical : public FE_Poly { public: /** diff --git a/include/deal.II/fe/fe_rannacher_turek.h b/include/deal.II/fe/fe_rannacher_turek.h index fedaa1e8a0..0ca8a0204f 100644 --- a/include/deal.II/fe/fe_rannacher_turek.h +++ b/include/deal.II/fe/fe_rannacher_turek.h @@ -81,7 +81,7 @@ DEAL_II_NAMESPACE_OPEN * @date 2015 */ template -class FE_RannacherTurek : public FE_Poly, dim> +class FE_RannacherTurek : public FE_Poly { public: /** diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index ca57050fc0..5827480fba 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -555,13 +555,9 @@ MatrixFree::is_supported( return false; // then check of the base element is supported - if (dynamic_cast, dim, spacedim> - *>(fe_ptr) != nullptr) + if (dynamic_cast *>(fe_ptr) != nullptr) return true; - if (dynamic_cast>, - dim, - spacedim> *>(fe_ptr) != nullptr) + if (dynamic_cast *>(fe_ptr) != nullptr) return true; if (dynamic_cast *>(fe_ptr) != nullptr) return true; diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 542bf4d6b2..4749787ebf 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -136,20 +136,11 @@ namespace internal Assert(fe->n_components() == 1, ExcMessage("Expected a scalar element")); - const FE_Poly, dim, dim> *fe_poly = - dynamic_cast< - const FE_Poly, dim, dim> *>(fe); - - const FE_Poly< - TensorProductPolynomials>, - dim, - dim> *fe_poly_piece = - dynamic_cast>, - dim, - dim> *>(fe); + const FE_Poly *fe_poly = + dynamic_cast *>(fe); + + const FE_Poly *fe_poly_piece = + dynamic_cast *>(fe); const FE_DGP *fe_dgp = dynamic_cast *>(fe); diff --git a/source/fe/fe_dgp.cc b/source/fe/fe_dgp.cc index 208f8705cb..04fe320a73 100644 --- a/source/fe/fe_dgp.cc +++ b/source/fe/fe_dgp.cc @@ -27,7 +27,7 @@ DEAL_II_NAMESPACE_OPEN template FE_DGP::FE_DGP(const unsigned int degree) - : FE_Poly, dim, spacedim>( + : FE_Poly( PolynomialSpace( Polynomials::Legendre::generate_complete_basis(degree)), FiniteElementData(get_dpo_vector(degree), diff --git a/source/fe/fe_dgp_monomial.cc b/source/fe/fe_dgp_monomial.cc index 314950e33d..f505bb986f 100644 --- a/source/fe/fe_dgp_monomial.cc +++ b/source/fe/fe_dgp_monomial.cc @@ -130,7 +130,7 @@ namespace internal template FE_DGPMonomial::FE_DGPMonomial(const unsigned int degree) - : FE_Poly, dim>( + : FE_Poly( PolynomialsP(degree), FiniteElementData(get_dpo_vector(degree), 1, @@ -143,8 +143,8 @@ FE_DGPMonomial::FE_DGPMonomial(const unsigned int degree) FiniteElementData(get_dpo_vector(degree), 1, degree).dofs_per_cell, std::vector(1, true))) { - Assert(this->poly_space.n() == this->dofs_per_cell, ExcInternalError()); - Assert(this->poly_space.degree() == this->degree, ExcInternalError()); + Assert(this->poly_space->n() == this->dofs_per_cell, ExcInternalError()); + Assert(this->poly_space->degree() == this->degree, ExcInternalError()); // DG doesn't have constraints, so // leave them empty @@ -233,9 +233,13 @@ FE_DGPMonomial::get_interpolation_matrix( source_fe_matrix(k, j) = source_fe.shape_value(j, unit_points[k]); FullMatrix this_matrix(this->dofs_per_cell, this->dofs_per_cell); + auto *const polynomial_space_p = + dynamic_cast *>(this->poly_space.get()); + Assert(polynomial_space_p != nullptr, ExcInternalError()); for (unsigned int j = 0; j < this->dofs_per_cell; ++j) for (unsigned int k = 0; k < unit_points.size(); ++k) - this_matrix(k, j) = this->poly_space.compute_value(j, unit_points[k]); + this_matrix(k, j) = + polynomial_space_p->compute_value(j, unit_points[k]); this_matrix.gauss_jordan(); @@ -454,8 +458,11 @@ FE_DGPMonomial<2>::has_support_on_face(const unsigned int shape_index, support_on_face = true; else { + auto *const polynomial_space_p = + dynamic_cast *>(this->poly_space.get()); + Assert(polynomial_space_p != nullptr, ExcInternalError()); const std::array degrees = - this->poly_space.directional_degrees(shape_index); + polynomial_space_p->directional_degrees(shape_index); if ((face_index == 0 && degrees[1] == 0) || (face_index == 3 && degrees[0] == 0)) @@ -476,8 +483,11 @@ FE_DGPMonomial<3>::has_support_on_face(const unsigned int shape_index, support_on_face = true; else { + auto *const polynomial_space_p = + dynamic_cast *>(this->poly_space.get()); + Assert(polynomial_space_p != nullptr, ExcInternalError()); const std::array degrees = - this->poly_space.directional_degrees(shape_index); + polynomial_space_p->directional_degrees(shape_index); if ((face_index == 0 && degrees[1] == 0) || (face_index == 2 && degrees[2] == 0) || diff --git a/source/fe/fe_dgq.cc b/source/fe/fe_dgq.cc index 818ddb0b6f..f54154595d 100644 --- a/source/fe/fe_dgq.cc +++ b/source/fe/fe_dgq.cc @@ -60,7 +60,7 @@ namespace internal template FE_DGQ::FE_DGQ(const unsigned int degree) - : FE_Poly, dim, spacedim>( + : FE_Poly( TensorProductPolynomials( Polynomials::generate_complete_Lagrange_basis( internal::FE_DGQ::get_QGaussLobatto_points(degree))), @@ -93,7 +93,7 @@ FE_DGQ::FE_DGQ(const unsigned int degree) template FE_DGQ::FE_DGQ( const std::vector> &polynomials) - : FE_Poly, dim, spacedim>( + : FE_Poly( TensorProductPolynomials(polynomials), FiniteElementData(get_dpo_vector(polynomials.size() - 1), 1, @@ -303,11 +303,14 @@ FE_DGQ::get_interpolation_matrix( // cell and evaluate the // shape functions there const Point p = this->unit_support_points[j]; + auto *const polynomial_space_p = + dynamic_cast *>(this->poly_space.get()); + Assert(polynomial_space_p != nullptr, ExcInternalError()); for (unsigned int i = 0; i < this->dofs_per_cell; ++i) - cell_interpolation(j, i) = this->poly_space.compute_value(i, p); + cell_interpolation(j, i) = polynomial_space_p->compute_value(i, p); for (unsigned int i = 0; i < source_fe.dofs_per_cell; ++i) - source_interpolation(j, i) = source_fe.poly_space.compute_value(i, p); + source_interpolation(j, i) = polynomial_space_p->compute_value(i, p); } // then compute the @@ -832,8 +835,11 @@ FE_DGQArbitraryNodes::get_name() const bool equidistant = true; std::vector points(this->degree + 1); + auto *const polynomial_space_p = + dynamic_cast *>(this->poly_space.get()); + Assert(polynomial_space_p != nullptr, ExcInternalError()); std::vector lexicographic = - this->poly_space.get_numbering_inverse(); + polynomial_space_p->get_numbering_inverse(); for (unsigned int j = 0; j <= this->degree; j++) points[j] = this->unit_support_points[lexicographic[j]][0]; @@ -945,9 +951,12 @@ std::unique_ptr> FE_DGQArbitraryNodes::clone() const { // Construct a dummy quadrature formula containing the FE's nodes: - std::vector> qpoints(this->degree + 1); + std::vector> qpoints(this->degree + 1); + auto *const polynomial_space_p = + dynamic_cast *>(this->poly_space.get()); + Assert(polynomial_space_p != nullptr, ExcInternalError()); std::vector lexicographic = - this->poly_space.get_numbering_inverse(); + polynomial_space_p->get_numbering_inverse(); for (unsigned int i = 0; i <= this->degree; ++i) qpoints[i] = Point<1>(this->unit_support_points[lexicographic[i]][0]); Quadrature<1> pquadrature(qpoints); diff --git a/source/fe/fe_poly.cc b/source/fe/fe_poly.cc index c17fd45f7d..bee64acba8 100644 --- a/source/fe/fe_poly.cc +++ b/source/fe/fe_poly.cc @@ -29,139 +29,12 @@ DEAL_II_NAMESPACE_OPEN -template <> -void -FE_Poly, 1, 2>::fill_fe_values( - const Triangulation<1, 2>::cell_iterator &, - const CellSimilarity::Similarity cell_similarity, - const Quadrature<1> & quadrature, - const Mapping<1, 2> & mapping, - const Mapping<1, 2>::InternalDataBase &mapping_internal, - const dealii::internal::FEValuesImplementation::MappingRelatedData<1, 2> - & mapping_data, - const FiniteElement<1, 2>::InternalDataBase &fe_internal, - dealii::internal::FEValuesImplementation::FiniteElementRelatedData<1, 2> - &output_data) const -{ - // convert data object to internal data for this class. fails with an - // exception if that is not possible - Assert(dynamic_cast(&fe_internal) != nullptr, - ExcInternalError()); - const InternalData &fe_data = - static_cast(fe_internal); // NOLINT +template +FE_Poly::FE_Poly(const FE_Poly &fe) + : FiniteElement(fe) + , poly_space(fe.poly_space->clone()) +{} - const bool need_to_correct_higher_derivatives = - higher_derivatives_need_correcting(mapping, - mapping_data, - quadrature.size(), - fe_data.update_each); - - // transform gradients and higher derivatives. there is nothing to do - // for values since we already emplaced them into output_data when - // we were in get_data() - if (fe_data.update_each & update_gradients && - cell_similarity != CellSimilarity::translation) - for (unsigned int k = 0; k < this->dofs_per_cell; ++k) - mapping.transform(make_array_view(fe_data.shape_gradients, k), - mapping_covariant, - mapping_internal, - make_array_view(output_data.shape_gradients, k)); - - if (fe_data.update_each & update_hessians && - cell_similarity != CellSimilarity::translation) - { - for (unsigned int k = 0; k < this->dofs_per_cell; ++k) - mapping.transform(make_array_view(fe_data.shape_hessians, k), - mapping_covariant_gradient, - mapping_internal, - make_array_view(output_data.shape_hessians, k)); - - if (need_to_correct_higher_derivatives) - correct_hessians(output_data, mapping_data, quadrature.size()); - } - - if (fe_data.update_each & update_3rd_derivatives && - cell_similarity != CellSimilarity::translation) - { - for (unsigned int k = 0; k < this->dofs_per_cell; ++k) - mapping.transform(make_array_view(fe_data.shape_3rd_derivatives, k), - mapping_covariant_hessian, - mapping_internal, - make_array_view(output_data.shape_3rd_derivatives, - k)); - - if (need_to_correct_higher_derivatives) - correct_third_derivatives(output_data, mapping_data, quadrature.size()); - } -} - - - -template <> -void -FE_Poly, 2, 3>::fill_fe_values( - const Triangulation<2, 3>::cell_iterator &, - const CellSimilarity::Similarity cell_similarity, - const Quadrature<2> & quadrature, - const Mapping<2, 3> & mapping, - const Mapping<2, 3>::InternalDataBase &mapping_internal, - const dealii::internal::FEValuesImplementation::MappingRelatedData<2, 3> - & mapping_data, - const FiniteElement<2, 3>::InternalDataBase &fe_internal, - dealii::internal::FEValuesImplementation::FiniteElementRelatedData<2, 3> - &output_data) const -{ - // convert data object to internal data for this class. fails with an - // exception if that is not possible - Assert(dynamic_cast(&fe_internal) != nullptr, - ExcInternalError()); - const InternalData &fe_data = - static_cast(fe_internal); // NOLINT - - const bool need_to_correct_higher_derivatives = - higher_derivatives_need_correcting(mapping, - mapping_data, - quadrature.size(), - fe_data.update_each); - - // transform gradients and higher derivatives. there is nothing to do - // for values since we already emplaced them into output_data when - // we were in get_data() - if (fe_data.update_each & update_gradients && - cell_similarity != CellSimilarity::translation) - for (unsigned int k = 0; k < this->dofs_per_cell; ++k) - mapping.transform(make_array_view(fe_data.shape_gradients, k), - mapping_covariant, - mapping_internal, - make_array_view(output_data.shape_gradients, k)); - - if (fe_data.update_each & update_hessians && - cell_similarity != CellSimilarity::translation) - { - for (unsigned int k = 0; k < this->dofs_per_cell; ++k) - mapping.transform(make_array_view(fe_data.shape_hessians, k), - mapping_covariant_gradient, - mapping_internal, - make_array_view(output_data.shape_hessians, k)); - - if (need_to_correct_higher_derivatives) - correct_hessians(output_data, mapping_data, quadrature.size()); - } - - if (fe_data.update_each & update_3rd_derivatives && - cell_similarity != CellSimilarity::translation) - { - for (unsigned int k = 0; k < this->dofs_per_cell; ++k) - mapping.transform(make_array_view(fe_data.shape_3rd_derivatives, k), - mapping_covariant_hessian, - mapping_internal, - make_array_view(output_data.shape_3rd_derivatives, - k)); - - if (need_to_correct_higher_derivatives) - correct_third_derivatives(output_data, mapping_data, quadrature.size()); - } -} #include "fe_poly.inst" diff --git a/source/fe/fe_poly.inst.in b/source/fe/fe_poly.inst.in index 486fbd92dd..dcaa87fc0d 100644 --- a/source/fe/fe_poly.inst.in +++ b/source/fe/fe_poly.inst.in @@ -18,28 +18,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) { #if deal_II_dimension <= deal_II_space_dimension - template class FE_Poly, - deal_II_dimension, - deal_II_space_dimension>; - template class FE_Poly, - deal_II_dimension, - deal_II_space_dimension>; - template class FE_Poly, - deal_II_dimension, - deal_II_space_dimension>; - template class FE_Poly< - TensorProductPolynomials>, - deal_II_dimension, - deal_II_space_dimension>; - template class FE_Poly, - deal_II_dimension, - deal_II_space_dimension>; - template class FE_Poly, - deal_II_dimension, - deal_II_space_dimension>; - template class FE_Poly, - deal_II_dimension, - deal_II_space_dimension>; + template class FE_Poly; #endif } diff --git a/source/fe/fe_q.cc b/source/fe/fe_q.cc index 006f9e1c4e..db825ba22d 100644 --- a/source/fe/fe_q.cc +++ b/source/fe/fe_q.cc @@ -99,8 +99,10 @@ FE_Q::get_name() const std::vector points(this->degree + 1); // Decode the support points in one coordinate direction. + TensorProductPolynomials *poly_space_derived_ptr = + dynamic_cast *>(this->poly_space.get()); std::vector lexicographic = - this->poly_space.get_numbering_inverse(); + poly_space_derived_ptr->get_numbering_inverse(); for (unsigned int j = 0; j <= this->degree; j++) points[j] = this->unit_support_points[lexicographic[j]][0]; diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index 4ab6de7034..dc183f5151 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -180,12 +180,14 @@ struct FE_Q_Base::Implementation // use that the element evaluates to 1 at index 0 and along the line at // zero + TensorProductPolynomials *poly_space_derived_ptr = + dynamic_cast *>(fe.poly_space.get()); const std::vector &index_map_inverse = - fe.poly_space.get_numbering_inverse(); + poly_space_derived_ptr->get_numbering_inverse(); const std::vector face_index_map = FETools::lexicographic_to_hierarchic_numbering(q_deg); Assert(std::abs( - fe.poly_space.compute_value(index_map_inverse[0], Point()) - + fe.poly_space->compute_value(index_map_inverse[0], Point()) - 1.) < 1e-14, ExcInternalError()); @@ -195,7 +197,7 @@ struct FE_Q_Base::Implementation Point p; p[0] = constraint_points[i](0); fe.interface_constraints(i, face_index_map[j]) = - fe.poly_space.compute_value(index_map_inverse[j], p); + fe.poly_space->compute_value(index_map_inverse[j], p); // 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 @@ -315,12 +317,14 @@ struct FE_Q_Base::Implementation // use that the element evaluates to 1 at index 0 and along the line at // zero + TensorProductPolynomials *poly_space_derived_ptr = + dynamic_cast *>(fe.poly_space.get()); const std::vector &index_map_inverse = - fe.poly_space.get_numbering_inverse(); + poly_space_derived_ptr->get_numbering_inverse(); const std::vector face_index_map = FETools::lexicographic_to_hierarchic_numbering(q_deg); Assert(std::abs( - fe.poly_space.compute_value(index_map_inverse[0], Point()) - + fe.poly_space->compute_value(index_map_inverse[0], Point()) - 1.) < 1e-14, ExcInternalError()); @@ -385,8 +389,8 @@ struct FE_Q_Base::Implementation indices[1] * (q_deg + 1) + indices[0]; fe.interface_constraints(i, face_index_map[j]) = - fe.poly_space.compute_value(index_map_inverse[new_index], - constraint_point); + fe.poly_space->compute_value(index_map_inverse[new_index], + constraint_point); // if the value is small up to round-off, then simply set it to // zero to avoid unwanted fill-in of the constraint matrices @@ -407,7 +411,7 @@ FE_Q_Base::FE_Q_Base( const PolynomialType & poly_space, const FiniteElementData &fe_data, const std::vector & restriction_is_additive_flags) - : FE_Poly( + : FE_Poly( poly_space, fe_data, restriction_is_additive_flags, @@ -445,7 +449,9 @@ FE_Q_Base::initialize( FETools::hierarchic_to_lexicographic_numbering(q_degree); for (unsigned int i = q_dofs_per_cell; i < this->dofs_per_cell; ++i) renumber.push_back(i); - this->poly_space.set_numbering(renumber); + TensorProductPolynomials *poly_space_derived_ptr = + dynamic_cast *>(this->poly_space.get()); + poly_space_derived_ptr->set_numbering(renumber); } // Finally fill in support points on cell and face and initialize @@ -502,12 +508,12 @@ FE_Q_Base::get_interpolation_matrix( // FE_Q element evaluates to 1 in unit support point and to zero in // all other points by construction - Assert(std::abs(this->poly_space.compute_value(j, p) - 1.) < 1e-13, + Assert(std::abs(this->poly_space->compute_value(j, p) - 1.) < 1e-13, ExcInternalError()); for (unsigned int i = 0; i < source_q_dofs_per_cell; ++i) interpolation_matrix(j, i) = - source_fe->poly_space.compute_value(i, p); + source_fe->poly_space->compute_value(i, p); } // for FE_Q_DG0, add one last row of identity @@ -748,10 +754,15 @@ FE_Q_Base::hp_line_dof_identities( std::vector> identities; + TensorProductPolynomials *poly_space_derived_ptr = + dynamic_cast *>(this->poly_space.get()); const std::vector &index_map_inverse = - this->poly_space.get_numbering_inverse(); + poly_space_derived_ptr->get_numbering_inverse(); + TensorProductPolynomials *poly_space_derived_ptr_other = + dynamic_cast *>( + fe_q_other->poly_space.get()); const std::vector &index_map_inverse_other = - fe_q_other->poly_space.get_numbering_inverse(); + poly_space_derived_ptr_other->get_numbering_inverse(); for (unsigned int i = 0; i < p - 1; ++i) for (unsigned int j = 0; j < q - 1; ++j) @@ -810,10 +821,15 @@ FE_Q_Base::hp_quad_dof_identities( std::vector> identities; + TensorProductPolynomials *poly_space_derived_ptr = + dynamic_cast *>(this->poly_space.get()); const std::vector &index_map_inverse = - this->poly_space.get_numbering_inverse(); + poly_space_derived_ptr->get_numbering_inverse(); + TensorProductPolynomials *poly_space_derived_ptr_other = + dynamic_cast *>( + fe_q_other->poly_space.get()); const std::vector &index_map_inverse_other = - fe_q_other->poly_space.get_numbering_inverse(); + poly_space_derived_ptr_other->get_numbering_inverse(); for (unsigned int i1 = 0; i1 < p - 1; ++i1) for (unsigned int i2 = 0; i2 < p - 1; ++i2) @@ -870,8 +886,10 @@ void FE_Q_Base::initialize_unit_support_points( const std::vector> &points) { + TensorProductPolynomials *poly_space_derived_ptr = + dynamic_cast *>(this->poly_space.get()); const std::vector &index_map_inverse = - this->poly_space.get_numbering_inverse(); + poly_space_derived_ptr->get_numbering_inverse(); // We can compute the support points by computing the tensor // product of the 1d set of points. We could do this by hand, but it's @@ -1175,7 +1193,7 @@ FE_Q_Base::get_prolongation_matrix( // current numbering gives the identity operation for (unsigned int i = 0; i < q_dofs_per_cell; ++i) { - Assert(std::fabs(1. - this->poly_space.compute_value( + Assert(std::fabs(1. - this->poly_space->compute_value( i, this->unit_support_points[i])) < eps, ExcInternalError("The Lagrange polynomial does not evaluate " "to one or zero in a nodal point. " @@ -1185,7 +1203,7 @@ FE_Q_Base::get_prolongation_matrix( "prevents the sum to be one.")); for (unsigned int j = 0; j < q_dofs_per_cell; ++j) if (j != i) - Assert(std::fabs(this->poly_space.compute_value( + Assert(std::fabs(this->poly_space->compute_value( i, this->unit_support_points[j])) < eps, ExcInternalError( "The Lagrange polynomial does not evaluate " @@ -1204,8 +1222,11 @@ FE_Q_Base::get_prolongation_matrix( const unsigned int dofs1d = q_degree + 1; std::vector> subcell_evaluations( dim, Table<2, double>(dofs1d, dofs1d)); + + TensorProductPolynomials *poly_space_derived_ptr = + dynamic_cast *>(this->poly_space.get()); const std::vector &index_map_inverse = - this->poly_space.get_numbering_inverse(); + poly_space_derived_ptr->get_numbering_inverse(); // helper value: step size how to walk through diagonal and how many // points we have left apart from the first dimension @@ -1238,7 +1259,7 @@ FE_Q_Base::get_prolongation_matrix( Point point; point[0] = p_cell[d]; const double cell_value = - this->poly_space.compute_value(index_map_inverse[i], point); + this->poly_space->compute_value(index_map_inverse[i], point); // cut off values that are too small. note that we have here // Lagrange interpolation functions, so they should be zero at @@ -1280,7 +1301,7 @@ FE_Q_Base::get_prolongation_matrix( subcell_evaluations[d](j_indices[d - 1], i_indices[d - 1]); // innermost sum where we actually compute. the same as - // prolongate(j,i) = this->poly_space.compute_value (i, p_cell) + // prolongate(j,i) = this->poly_space->compute_value (i, p_cell) for (unsigned int jj = 0; jj < dofs1d; ++jj) { const unsigned int j_ind = index_map_inverse[j + jj]; @@ -1379,9 +1400,12 @@ FE_Q_Base::get_restriction_matrix( // assumption that whenever a row makes a non-zero contribution to the // mother's residual, the correct value is interpolated. - const double eps = 1e-15 * q_degree * dim; + const double eps = 1e-15 * q_degree * dim; + + TensorProductPolynomials *poly_space_derived_ptr = + dynamic_cast *>(this->poly_space.get()); const std::vector &index_map_inverse = - this->poly_space.get_numbering_inverse(); + poly_space_derived_ptr->get_numbering_inverse(); const unsigned int dofs1d = q_degree + 1; std::vector> evaluations1d(dofs1d); @@ -1411,8 +1435,8 @@ FE_Q_Base::get_restriction_matrix( Point point; point[0] = p_subcell[d]; evaluations1d[j][d] = - this->poly_space.compute_value(index_map_inverse[j], - point); + this->poly_space->compute_value(index_map_inverse[j], + point); } unsigned int j_indices[dim]; internal::FE_Q_Base::zero_indices(j_indices); diff --git a/source/fe/fe_q_hierarchical.cc b/source/fe/fe_q_hierarchical.cc index d77f389ad5..41e7021282 100644 --- a/source/fe/fe_q_hierarchical.cc +++ b/source/fe/fe_q_hierarchical.cc @@ -58,8 +58,9 @@ namespace internal template FE_Q_Hierarchical::FE_Q_Hierarchical(const unsigned int degree) - : FE_Poly, dim>( - Polynomials::Hierarchical::generate_complete_basis(degree), + : FE_Poly( + TensorProductPolynomials( + Polynomials::Hierarchical::generate_complete_basis(degree)), FiniteElementData(get_dpo_vector(degree), 1, degree, @@ -72,7 +73,9 @@ FE_Q_Hierarchical::FE_Q_Hierarchical(const unsigned int degree) std::vector(1, true))) , face_renumber(face_fe_q_hierarchical_to_hierarchic_numbering(degree)) { - this->poly_space.set_numbering( + TensorProductPolynomials *poly_space_derived_ptr = + dynamic_cast *>(this->poly_space.get()); + poly_space_derived_ptr->set_numbering( hierarchic_to_fe_q_hierarchical_numbering(*this)); // The matrix @p{dofs_cell} contains the @@ -673,7 +676,10 @@ FE_Q_Hierarchical::initialize_embedding_and_restriction( unsigned int iso = RefinementCase::isotropic_refinement - 1; const unsigned int dofs_1d = 2 * this->dofs_per_vertex + this->dofs_per_line; - const std::vector &renumber = this->poly_space.get_numbering(); + TensorProductPolynomials *poly_space_derived_ptr = + dynamic_cast *>(this->poly_space.get()); + const std::vector &renumber = + poly_space_derived_ptr->get_numbering(); for (unsigned int c = 0; c < GeometryInfo::max_children_per_cell; ++c) { @@ -810,8 +816,10 @@ FE_Q_Hierarchical::initialize_generalized_support_points() this->generalized_support_points.resize(n); + TensorProductPolynomials *poly_space_derived_ptr = + dynamic_cast *>(this->poly_space.get()); const std::vector &index_map_inverse = - this->poly_space.get_numbering_inverse(); + poly_space_derived_ptr->get_numbering_inverse(); Point p; // the method of numbering allows diff --git a/source/fe/fe_rannacher_turek.cc b/source/fe/fe_rannacher_turek.cc index 24d4a5714b..dded3b04f5 100644 --- a/source/fe/fe_rannacher_turek.cc +++ b/source/fe/fe_rannacher_turek.cc @@ -32,14 +32,13 @@ template FE_RannacherTurek::FE_RannacherTurek( const unsigned int order, const unsigned int n_face_support_points) - : FE_Poly, dim>( - PolynomialsRannacherTurek(), - FiniteElementData(this->get_dpo_vector(), - 1, - 2, - FiniteElementData::L2), - std::vector(4, false), // restriction not implemented - std::vector(4, std::vector(1, true))) + : FE_Poly(PolynomialsRannacherTurek(), + FiniteElementData(this->get_dpo_vector(), + 1, + 2, + FiniteElementData::L2), + std::vector(4, false), // restriction not implemented + std::vector(4, std::vector(1, true))) , order(order) , n_face_support_points(n_face_support_points) {