From db0ba4dd6c5f8c093e2548aa85d24f45e306aadd Mon Sep 17 00:00:00 2001 From: grahambenharper Date: Wed, 11 Sep 2019 22:22:31 -0600 Subject: [PATCH] Change TensorProductPolynomialsBubbles to no longer derive from TensorProductPolynomials --- .../deal.II/base/tensor_product_polynomials.h | 9 ++ .../base/tensor_product_polynomials_bubbles.h | 95 ++++++++++++++++--- .../tensor_product_polynomials_bubbles.cc | 85 ++++++++++++----- 3 files changed, 151 insertions(+), 38 deletions(-) diff --git a/include/deal.II/base/tensor_product_polynomials.h b/include/deal.II/base/tensor_product_polynomials.h index 9c39dec516..3fd35d075c 100644 --- a/include/deal.II/base/tensor_product_polynomials.h +++ b/include/deal.II/base/tensor_product_polynomials.h @@ -29,6 +29,9 @@ DEAL_II_NAMESPACE_OPEN +template +class TensorProductPolynomialsBubbles; + /** * @addtogroup Polynomials * @{ @@ -227,6 +230,12 @@ protected: void compute_index(const unsigned int i, unsigned int (&indices)[(dim > 0 ? dim : 1)]) const; + + /** + * TensorProductPolynomialsBubbles has a TensorProductPolynomials class + * so we declare it as a friend class. + */ + friend class TensorProductPolynomialsBubbles; }; diff --git a/include/deal.II/base/tensor_product_polynomials_bubbles.h b/include/deal.II/base/tensor_product_polynomials_bubbles.h index 7a0f3dd018..af66095a6d 100644 --- a/include/deal.II/base/tensor_product_polynomials_bubbles.h +++ b/include/deal.II/base/tensor_product_polynomials_bubbles.h @@ -30,7 +30,6 @@ DEAL_II_NAMESPACE_OPEN - /** * @addtogroup Polynomials * @{ @@ -45,9 +44,15 @@ DEAL_II_NAMESPACE_OPEN * @author Daniel Arndt, 2015 */ template -class TensorProductPolynomialsBubbles : public TensorProductPolynomials +class TensorProductPolynomialsBubbles { public: + /** + * Access to the dimension of this object, for checking and automatic + * setting of dimension in other classes. + */ + static const unsigned int dimension = dim; + /** * Constructor. pols is a vector of objects that should be derived * or otherwise convertible to one-dimensional polynomial objects. It will @@ -56,6 +61,32 @@ public: template TensorProductPolynomialsBubbles(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. @@ -145,6 +176,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; }; /** @} */ @@ -158,29 +205,31 @@ template template inline TensorProductPolynomialsBubbles::TensorProductPolynomialsBubbles( const std::vector &pols) - : TensorProductPolynomials(pols) + : tensor_polys(pols) + , index_map(tensor_polys.n() + + ((tensor_polys.polynomials.size() <= 2) ? 1 : dim)) + , index_map_inverse(tensor_polys.n() + + ((tensor_polys.polynomials.size() <= 2) ? 1 : dim)) { - const unsigned int q_degree = this->polynomials.size() - 1; + const unsigned int q_degree = tensor_polys.polynomials.size() - 1; const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim); // append index for renumbering - for (unsigned int i = 0; i < n_bubbles; ++i) + for (unsigned int i = 0; i < tensor_polys.n() + n_bubbles; ++i) { - this->index_map.push_back(i + this->n_tensor_pols); - this->index_map_inverse.push_back(i + this->n_tensor_pols); + index_map[i] = i; + index_map_inverse[i] = i; } } - template inline unsigned int TensorProductPolynomialsBubbles::n() const { - return this->n_tensor_pols + dim; + return tensor_polys.n() + dim; } - template <> inline unsigned int TensorProductPolynomialsBubbles<0>::n() const @@ -188,6 +237,23 @@ TensorProductPolynomialsBubbles<0>::n() const return numbers::invalid_unsigned_int; } + +template +inline const std::vector & +TensorProductPolynomialsBubbles::get_numbering() const +{ + return index_map; +} + + +template +inline const std::vector & +TensorProductPolynomialsBubbles::get_numbering_inverse() const +{ + return index_map_inverse; +} + + template template Tensor @@ -195,18 +261,17 @@ TensorProductPolynomialsBubbles::compute_derivative( const unsigned int i, const Point & p) const { - const unsigned int q_degree = this->polynomials.size() - 1; - const unsigned int max_q_indices = this->n_tensor_pols; + const unsigned int q_degree = tensor_polys.polynomials.size() - 1; + const unsigned int max_q_indices = tensor_polys.n(); const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim); (void)n_bubbles; Assert(i < max_q_indices + n_bubbles, ExcInternalError()); // treat the regular basis functions if (i < max_q_indices) - return this - ->TensorProductPolynomials::template compute_derivative(i, p); + return tensor_polys.template compute_derivative(i, p); - const unsigned int comp = i - this->n_tensor_pols; + const unsigned int comp = i - tensor_polys.n(); Tensor derivative; switch (order) diff --git a/source/base/tensor_product_polynomials_bubbles.cc b/source/base/tensor_product_polynomials_bubbles.cc index afd88bfb82..dd231f8304 100644 --- a/source/base/tensor_product_polynomials_bubbles.cc +++ b/source/base/tensor_product_polynomials_bubbles.cc @@ -24,22 +24,60 @@ DEAL_II_NAMESPACE_OPEN /* ------------------- TensorProductPolynomialsBubbles -------------- */ + +template +void +TensorProductPolynomialsBubbles::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 +TensorProductPolynomialsBubbles::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 TensorProductPolynomialsBubbles::compute_value(const unsigned int i, const Point & p) const { - const unsigned int q_degree = this->polynomials.size() - 1; - const unsigned int max_q_indices = this->n_tensor_pols; + const unsigned int q_degree = tensor_polys.polynomials.size() - 1; + const unsigned int max_q_indices = tensor_polys.n(); const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim); (void)n_bubbles; Assert(i < max_q_indices + n_bubbles, ExcInternalError()); // treat the regular basis functions if (i < max_q_indices) - return this->TensorProductPolynomials::compute_value(i, p); + return tensor_polys.compute_value(i, p); - const unsigned int comp = i - this->n_tensor_pols; + const unsigned int comp = i - tensor_polys.n(); // compute \prod_{i=1}^d 4*(1-x_i^2)(p) double value = 1.; @@ -63,22 +101,23 @@ TensorProductPolynomialsBubbles<0>::compute_value(const unsigned int, } + template Tensor<1, dim> TensorProductPolynomialsBubbles::compute_grad(const unsigned int i, const Point & p) const { - const unsigned int q_degree = this->polynomials.size() - 1; - const unsigned int max_q_indices = this->n_tensor_pols; + const unsigned int q_degree = tensor_polys.polynomials.size() - 1; + const unsigned int max_q_indices = tensor_polys.n(); const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim); (void)n_bubbles; Assert(i < max_q_indices + n_bubbles, ExcInternalError()); // treat the regular basis functions if (i < max_q_indices) - return this->TensorProductPolynomials::compute_grad(i, p); + return tensor_polys.compute_grad(i, p); - const unsigned int comp = i - this->n_tensor_pols; + const unsigned int comp = i - tensor_polys.n(); Tensor<1, dim> grad; for (unsigned int d = 0; d < dim; ++d) @@ -116,17 +155,17 @@ TensorProductPolynomialsBubbles::compute_grad_grad( const unsigned int i, const Point & p) const { - const unsigned int q_degree = this->polynomials.size() - 1; - const unsigned int max_q_indices = this->n_tensor_pols; + const unsigned int q_degree = tensor_polys.polynomials.size() - 1; + const unsigned int max_q_indices = tensor_polys.n(); const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim); (void)n_bubbles; Assert(i < max_q_indices + n_bubbles, ExcInternalError()); // treat the regular basis functions if (i < max_q_indices) - return this->TensorProductPolynomials::compute_grad_grad(i, p); + return tensor_polys.compute_grad_grad(i, p); - const unsigned int comp = i - this->n_tensor_pols; + const unsigned int comp = i - tensor_polys.n(); double v[dim + 1][3]; { @@ -213,6 +252,8 @@ TensorProductPolynomialsBubbles::compute_grad_grad( return grad_grad; } + + template void TensorProductPolynomialsBubbles::evaluate( @@ -223,8 +264,8 @@ TensorProductPolynomialsBubbles::evaluate( std::vector> &third_derivatives, std::vector> &fourth_derivatives) const { - const unsigned int q_degree = this->polynomials.size() - 1; - const unsigned int max_q_indices = this->n_tensor_pols; + const unsigned int q_degree = tensor_polys.polynomials.size() - 1; + const unsigned int max_q_indices = tensor_polys.n(); (void)max_q_indices; const unsigned int n_bubbles = ((q_degree <= 1) ? 1 : dim); Assert(values.size() == max_q_indices + n_bubbles || values.size() == 0, @@ -249,36 +290,34 @@ TensorProductPolynomialsBubbles::evaluate( bool do_3rd_derivatives = false, do_4th_derivatives = false; if (values.empty() == false) { - values.resize(this->n_tensor_pols); + values.resize(tensor_polys.n()); do_values = true; } if (grads.empty() == false) { - grads.resize(this->n_tensor_pols); + grads.resize(tensor_polys.n()); do_grads = true; } if (grad_grads.empty() == false) { - grad_grads.resize(this->n_tensor_pols); + grad_grads.resize(tensor_polys.n()); do_grad_grads = true; } if (third_derivatives.empty() == false) { - third_derivatives.resize(this->n_tensor_pols); + third_derivatives.resize(tensor_polys.n()); do_3rd_derivatives = true; } if (fourth_derivatives.empty() == false) { - fourth_derivatives.resize(this->n_tensor_pols); + 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 (unsigned int i = this->n_tensor_pols; - i < this->n_tensor_pols + n_bubbles; - ++i) + for (unsigned int i = tensor_polys.n(); i < tensor_polys.n() + n_bubbles; ++i) { if (do_values) values.push_back(compute_value(i, p)); -- 2.39.5