From 10512d80e3505eb8b2fc12617fd4f4acaafb001f Mon Sep 17 00:00:00 2001 From: grahambenharper Date: Mon, 16 Sep 2019 10:13:57 -0600 Subject: [PATCH] Change TensorProductPolynomialsConst to no longer derive from TensorProductPolynomials --- .../base/tensor_product_polynomials_const.h | 78 ++++++++++++++--- .../base/tensor_product_polynomials_const.cc | 85 +++++++++++++------ 2 files changed, 124 insertions(+), 39 deletions(-) diff --git a/include/deal.II/base/tensor_product_polynomials_const.h b/include/deal.II/base/tensor_product_polynomials_const.h index a9b1e131f8..0f8cb93d0c 100644 --- a/include/deal.II/base/tensor_product_polynomials_const.h +++ b/include/deal.II/base/tensor_product_polynomials_const.h @@ -45,7 +45,7 @@ DEAL_II_NAMESPACE_OPEN * @author Timo Heister, 2012 */ template -class TensorProductPolynomialsConst : public TensorProductPolynomials +class TensorProductPolynomialsConst { public: /** @@ -62,6 +62,32 @@ public: template TensorProductPolynomialsConst(const std::vector &pols); + /** + * Print the list of tensor_polys indices to out. + */ + void + output_indices(std::ostream &out) const; + + /** + * Set the ordering of the polynomials. Requires + * renumber.size()==tensor_polys.n(). Stores a copy of + * renumber. + */ + void + set_numbering(const std::vector &renumber); + + /** + * Give read access to the renumber vector. + */ + const std::vector & + get_numbering() const; + + /** + * Give read access to the inverse renumber vector. + */ + const std::vector & + get_numbering_inverse() const; + /** * Compute the value and the first and second derivatives of each tensor * product polynomial at unit_point. @@ -151,6 +177,22 @@ public: */ unsigned int n() const; + +private: + /** + * The TensorProductPolynomials object + */ + TensorProductPolynomials tensor_polys; + + /** + * Index map for reordering the polynomials. + */ + std::vector index_map; + + /** + * Index map for reordering the polynomials. + */ + std::vector index_map_inverse; }; /** @} */ @@ -164,12 +206,10 @@ template template inline TensorProductPolynomialsConst::TensorProductPolynomialsConst( const std::vector &pols) - : TensorProductPolynomials(pols) -{ - // append index for renumbering - this->index_map.push_back(TensorProductPolynomials::n()); - this->index_map_inverse.push_back(TensorProductPolynomials::n()); -} + : tensor_polys(pols) + , index_map(tensor_polys.n() + 1) + , index_map_inverse(tensor_polys.n() + 1) +{} @@ -177,7 +217,24 @@ template inline unsigned int TensorProductPolynomialsConst::n() const { - return TensorProductPolynomials::n() + 1; + return tensor_polys.n() + 1; +} + + + +template +inline const std::vector & +TensorProductPolynomialsConst::get_numbering() const +{ + return index_map; +} + + +template +inline const std::vector & +TensorProductPolynomialsConst::get_numbering_inverse() const +{ + return index_map_inverse; } @@ -196,13 +253,12 @@ TensorProductPolynomialsConst::compute_derivative( const unsigned int i, const Point & p) const { - const unsigned int max_indices = TensorProductPolynomials::n(); + const unsigned int max_indices = tensor_polys.n(); Assert(i <= max_indices, ExcInternalError()); // treat the regular basis functions if (i < max_indices) - return this - ->TensorProductPolynomials::template compute_derivative(i, p); + return tensor_polys.template compute_derivative(i, p); else // this is for the constant function return Tensor(); diff --git a/source/base/tensor_product_polynomials_const.cc b/source/base/tensor_product_polynomials_const.cc index 54286ba84c..c9b988073d 100644 --- a/source/base/tensor_product_polynomials_const.cc +++ b/source/base/tensor_product_polynomials_const.cc @@ -24,17 +24,55 @@ DEAL_II_NAMESPACE_OPEN /* ------------------- TensorProductPolynomialsConst -------------- */ + +template +void +TensorProductPolynomialsConst::output_indices(std::ostream &out) const +{ + unsigned int ix[dim]; + for (unsigned int i = 0; i < tensor_polys.n(); ++i) + { + tensor_polys.compute_index(i, ix); + out << i << "\t"; + for (unsigned int d = 0; d < dim; ++d) + out << ix[d] << " "; + out << std::endl; + } +} + + + +template +void +TensorProductPolynomialsConst::set_numbering( + const std::vector &renumber) +{ + Assert(renumber.size() == index_map.size(), + ExcDimensionMismatch(renumber.size(), index_map.size())); + + index_map = renumber; + for (unsigned int i = 0; i < index_map.size(); ++i) + index_map_inverse[index_map[i]] = i; + + std::vector renumber_base; + for (unsigned int i = 0; i < tensor_polys.n(); ++i) + renumber_base.push_back(renumber[i]); + + tensor_polys.set_numbering(renumber_base); +} + + template double TensorProductPolynomialsConst::compute_value(const unsigned int i, const Point & p) const { - const unsigned int max_indices = TensorProductPolynomials::n(); + const unsigned int max_indices = tensor_polys.n(); Assert(i <= max_indices, ExcInternalError()); // treat the regular basis functions if (i < max_indices) - return this->TensorProductPolynomials::compute_value(i, p); + return tensor_polys.compute_value(i, p); else // this is for the constant function return 1.; @@ -57,12 +95,12 @@ Tensor<1, dim> TensorProductPolynomialsConst::compute_grad(const unsigned int i, const Point & p) const { - const unsigned int max_indices = TensorProductPolynomials::n(); + const unsigned int max_indices = tensor_polys.n(); Assert(i <= max_indices, ExcInternalError()); // treat the regular basis functions if (i < max_indices) - return this->TensorProductPolynomials::compute_grad(i, p); + return tensor_polys.compute_grad(i, p); else // this is for the constant function return Tensor<1, dim>(); @@ -73,12 +111,12 @@ Tensor<2, dim> TensorProductPolynomialsConst::compute_grad_grad(const unsigned int i, const Point &p) const { - const unsigned int max_indices = TensorProductPolynomials::n(); + const unsigned int max_indices = tensor_polys.n(); Assert(i <= max_indices, ExcInternalError()); // treat the regular basis functions if (i < max_indices) - return this->TensorProductPolynomials::compute_grad_grad(i, p); + return tensor_polys.compute_grad_grad(i, p); else // this is for the constant function return Tensor<2, dim>(); @@ -94,30 +132,21 @@ TensorProductPolynomialsConst::evaluate( std::vector> &third_derivatives, std::vector> &fourth_derivatives) const { - Assert(values.size() == TensorProductPolynomials::n() + 1 || - values.size() == 0, - ExcDimensionMismatch2(values.size(), - TensorProductPolynomials::n() + 1, - 0)); - Assert(grads.size() == TensorProductPolynomials::n() + 1 || - grads.size() == 0, - ExcDimensionMismatch2(grads.size(), - TensorProductPolynomials::n() + 1, - 0)); - Assert(grad_grads.size() == TensorProductPolynomials::n() + 1 || - grad_grads.size() == 0, - ExcDimensionMismatch2(grad_grads.size(), - TensorProductPolynomials::n() + 1, - 0)); - Assert(third_derivatives.size() == TensorProductPolynomials::n() + 1 || + Assert(values.size() == tensor_polys.n() + 1 || values.size() == 0, + ExcDimensionMismatch2(values.size(), tensor_polys.n() + 1, 0)); + Assert(grads.size() == tensor_polys.n() + 1 || grads.size() == 0, + ExcDimensionMismatch2(grads.size(), tensor_polys.n() + 1, 0)); + Assert(grad_grads.size() == tensor_polys.n() + 1 || grad_grads.size() == 0, + ExcDimensionMismatch2(grad_grads.size(), tensor_polys.n() + 1, 0)); + Assert(third_derivatives.size() == tensor_polys.n() + 1 || third_derivatives.size() == 0, ExcDimensionMismatch2(third_derivatives.size(), - TensorProductPolynomials::n() + 1, + tensor_polys.n() + 1, 0)); - Assert(fourth_derivatives.size() == TensorProductPolynomials::n() + 1 || + Assert(fourth_derivatives.size() == tensor_polys.n() + 1 || fourth_derivatives.size() == 0, ExcDimensionMismatch2(fourth_derivatives.size(), - TensorProductPolynomials::n() + 1, + tensor_polys.n() + 1, 0)); // remove slot for const value, go into the base class compute method and @@ -141,16 +170,16 @@ TensorProductPolynomialsConst::evaluate( } if (third_derivatives.empty() == false) { - third_derivatives.resize(TensorProductPolynomials::n()); + third_derivatives.resize(tensor_polys.n()); do_3rd_derivatives = true; } if (fourth_derivatives.empty() == false) { - fourth_derivatives.resize(TensorProductPolynomials::n()); + fourth_derivatives.resize(tensor_polys.n()); do_4th_derivatives = true; } - this->TensorProductPolynomials::evaluate( + tensor_polys.evaluate( p, values, grads, grad_grads, third_derivatives, fourth_derivatives); // for dgq node: values =1, grads=0, grads_grads=0, third_derivatives=0, -- 2.39.5