From 7acd5ab7ab9db3f9020f2542769527757f9505e0 Mon Sep 17 00:00:00 2001 From: David Wells Date: Thu, 14 Jan 2021 19:31:13 -0500 Subject: [PATCH] Remove Simplex::ScalarPolynomial. --- include/deal.II/simplex/polynomials.h | 160 +------ source/simplex/polynomials.cc | 578 -------------------------- 2 files changed, 3 insertions(+), 735 deletions(-) diff --git a/include/deal.II/simplex/polynomials.h b/include/deal.II/simplex/polynomials.h index 31acf50362..558e1f97a4 100644 --- a/include/deal.II/simplex/polynomials.h +++ b/include/deal.II/simplex/polynomials.h @@ -33,129 +33,14 @@ DEAL_II_NAMESPACE_OPEN */ namespace Simplex { - /** - * Polynomials defined on dim-dimensional simplex entities. This class is - * basis of Simplex::FE_P. - * - * @ingroup simplex - */ - template - class ScalarPolynomial : public ScalarPolynomialsBase - { - public: - /** - * Make the dimension available to the outside. - */ - static const unsigned int dimension = dim; - - /* - * Constructor taking the polynomial @p degree as input. - * - * @note Currently, only linear (degree=1) and quadratic polynomials - * (degree=2) are implemented. - */ - ScalarPolynomial(const unsigned int degree); - - /** - * @copydoc ScalarPolynomialsBase::evaluate() - * - * @note Currently, only the vectors @p values, @p grads, and @p grad_grads - * are filled. - */ - void - evaluate(const Point & unit_point, - std::vector & values, - std::vector> &grads, - std::vector> &grad_grads, - std::vector> &third_derivatives, - std::vector> &fourth_derivatives) const override; - - /** - * @copydoc ScalarPolynomialsBase::compute_value() - */ - double - compute_value(const unsigned int i, const Point &p) const override; - - /** - * @copydoc ScalarPolynomialsBase::compute_derivative() - * - * @note Currently, only implemented for first and second derivative. - */ - template - Tensor - compute_derivative(const unsigned int i, const Point &p) const; - - /** - * @copydoc ScalarPolynomialsBase::compute_1st_derivative() - */ - Tensor<1, dim> - compute_1st_derivative(const unsigned int i, - const Point & p) const override; - - /** - * @copydoc ScalarPolynomialsBase::compute_2nd_derivative() - */ - Tensor<2, dim> - compute_2nd_derivative(const unsigned int i, - const Point & p) const override; - - /** - * @copydoc ScalarPolynomialsBase::compute_3rd_derivative() - * - * @note Not implemented yet. - */ - Tensor<3, dim> - compute_3rd_derivative(const unsigned int i, - const Point & p) const override; - - /** - * @copydoc ScalarPolynomialsBase::compute_4th_derivative() - * - * @note Not implemented yet. - */ - Tensor<4, dim> - compute_4th_derivative(const unsigned int i, - const Point & p) const override; - - /** - * @copydoc ScalarPolynomialsBase::compute_grad() - * - * @note Not implemented yet. - */ - Tensor<1, dim> - compute_grad(const unsigned int i, const Point &p) const override; - - /** - * @copydoc ScalarPolynomialsBase::compute_grad_grad() - * - * @note Not implemented yet. - */ - Tensor<2, dim> - compute_grad_grad(const unsigned int i, const Point &p) const override; - - /** - * @copydoc ScalarPolynomialsBase::name() - */ - std::string - name() const override; - - /** - * @copydoc ScalarPolynomialsBase::clone() - */ - virtual std::unique_ptr> - clone() const override; - }; - - - /** * Polynomials defined on wedge entities. This class is basis of * Simplex::FE_WedgeP. * * The polynomials are created via a tensor product of a - * Simplex::ScalarPolynomial<2>(degree) and a - * Simplex::ScalarPolynomial<1>(degree), however, are re-numerated to better - * match the definition of FiniteElement. + * Simplex::BarycentricPolynomials<2>::get_fe_p_basis(degree) and a + * Simplex::BarycentricPolynomials<1>::get_fe_p_basis(degree), however, are + * re-numerated to better match the definition of FiniteElement. */ template class ScalarWedgePolynomial : public ScalarPolynomialsBase @@ -389,45 +274,6 @@ namespace Simplex - // template functions - template - template - Tensor - ScalarPolynomial::compute_derivative(const unsigned int i, - const Point & p) const - { - Tensor derivative; - - if (order == 1) - { - Tensor<1, dim> &derivative_1 = - *reinterpret_cast *>(&derivative); - - const auto grad = compute_grad(i, p); - for (unsigned int i = 0; i < dim; ++i) - derivative_1[i] = grad[i]; - } - else if (order == 2) - { - Tensor<2, dim> &derivative_2 = - *reinterpret_cast *>(&derivative); - - const auto grad_grad = compute_grad_grad(i, p); - - for (unsigned int i = 0; i < dim; ++i) - for (unsigned int j = 0; j < dim; ++j) - derivative_2[i][j] = grad_grad[i][j]; - } - else - { - Assert(false, ExcNotImplemented()); - } - - return derivative; - } - - - template template Tensor diff --git a/source/simplex/polynomials.cc b/source/simplex/polynomials.cc index 58910a9a87..e423e5e7ff 100644 --- a/source/simplex/polynomials.cc +++ b/source/simplex/polynomials.cc @@ -23,36 +23,6 @@ namespace Simplex { namespace { - unsigned int - compute_n_polynomials(const unsigned int dim, const unsigned int degree) - { - if (dim == 1) - { - if (degree == 1) - return 2; - if (degree == 2) - return 3; - } - else if (dim == 2) - { - if (degree == 1) - return 3; - if (degree == 2) - return 6; - } - else if (dim == 3) - { - if (degree == 1) - return 4; - if (degree == 2) - return 10; - } - - Assert(false, ExcNotImplemented()); - - return 0; - } - unsigned int compute_n_polynomials_pyramid(const unsigned int dim, const unsigned int degree) @@ -88,551 +58,6 @@ namespace Simplex - template - ScalarPolynomial::ScalarPolynomial(const unsigned int degree) - : ScalarPolynomialsBase(degree, compute_n_polynomials(dim, degree)) - {} - - - - template - double - ScalarPolynomial::compute_value(const unsigned int i, - const Point & p) const - { - if (dim == 1) - { - if (this->degree() == 1) - { - if (i == 0) - return 1.0 - p[0]; - else if (i == 1) - return p[0]; - } - else if (this->degree() == 2) - { - if (i == 0) - return 2.0 * p[0] * p[0] - 3.0 * p[0] + 1; - else if (i == 1) - return 2.0 * p[0] * p[0] - p[0]; - else if (i == 2) - return -4.0 * p[0] * p[0] + 4.0 * p[0]; - } - } - else if (dim == 2) - { - if (this->degree() == 1) - { - if (i == 0) - return 1.0 - p[0] - p[1]; - else if (i == 1) - return p[0]; - else if (i == 2) - return p[1]; - } - else if (this->degree() == 2) - { - const double t1 = 1.0 - p[0] - p[1]; - const double t2 = p[0]; - const double t3 = p[1]; - - if (i == 0) - return t1 * (2.0 * t1 - 1.0); - else if (i == 1) - return t2 * (2.0 * t2 - 1.0); - else if (i == 2) - return t3 * (2.0 * t3 - 1.0); - else if (i == 3) - return 4.0 * t2 * t1; - else if (i == 4) - return 4.0 * t2 * t3; - else if (i == 5) - return 4.0 * t3 * t1; - } - } - else if (dim == 3) - { - if (this->degree() == 1) - { - if (i == 0) - return 1.0 - p[0] - p[1] - p[2]; - else if (i == 1) - return p[0]; - else if (i == 2) - return p[1]; - else if (i == 3) - return p[2]; - } - else if (this->degree() == 2) - { - const double r = p[0]; - const double s = p[1]; - const double t = p[2]; - const double u = 1.0 - p[0] - p[1] - p[2]; - if (i == 0) - return u * (2.0 * u - 1.0); - else if (i == 1) - return r * (2.0 * r - 1.0); - else if (i == 2) - return s * (2.0 * s - 1.0); - else if (i == 3) - return t * (2.0 * t - 1.0); - else if (i == 4) - return 4.0 * r * u; - else if (i == 5) - return 4.0 * r * s; - else if (i == 6) - return 4.0 * s * u; - else if (i == 7) - return 4.0 * t * u; - else if (i == 8) - return 4.0 * r * t; - else if (i == 9) - return 4.0 * s * t; - } - } - - Assert(false, ExcNotImplemented()); - - return 0; - } - - - - template - Tensor<1, dim> - ScalarPolynomial::compute_grad(const unsigned int i, - const Point & p) const - { - Tensor<1, dim> grad; - - if (dim == 1) - { - if (this->degree() == 1) - { - if (i == 0) - grad[0] = -1.0; - else if (i == 1) - grad[0] = 1.0; - } - else if (this->degree() == 2) - { - if (i == 0) - grad[0] = 4.0 * p[0] - 3.0; - else if (i == 1) - grad[0] = 4.0 * p[0] - 1.0; - else if (i == 2) - grad[0] = -8.0 * p[0] + 4.0; - } - else - { - Assert(false, ExcNotImplemented()); - } - } - else if (dim == 2) - { - if (this->degree() == 1) - { - if (i == 0) - { - grad[0] = -1.0; - grad[1] = -1.0; - } - else if (i == 1) - { - grad[0] = +1.0; - grad[1] = +0.0; - } - else if (i == 2) - { - grad[0] = +0.0; - grad[1] = +1.0; - } - else - { - Assert(false, ExcNotImplemented()); - } - } - else if (this->degree() == 2) - { - if (i == 0) - { - grad[0] = -3.0 + 4.0 * (p[0] + p[1]); - grad[1] = -3.0 + 4.0 * (p[0] + p[1]); - } - else if (i == 1) - { - grad[0] = 4.0 * p[0] - 1.0; - grad[1] = 0.0; - } - else if (i == 2) - { - grad[0] = 0.0; - grad[1] = 4.0 * p[1] - 1.0; - } - else if (i == 3) - { - grad[0] = 4.0 * (1.0 - 2.0 * p[0] - p[1]); - grad[1] = -4.0 * p[0]; - } - else if (i == 4) - { - grad[0] = 4.0 * p[1]; - grad[1] = 4.0 * p[0]; - } - else if (i == 5) - { - grad[0] = -4.0 * p[1]; - grad[1] = 4.0 * (1.0 - p[0] - 2.0 * p[1]); - } - else - { - Assert(false, ExcNotImplemented()); - } - } - else - { - Assert(false, ExcNotImplemented()); - } - } - else if (dim == 3) - { - if (this->degree() == 1) - { - if (i == 0) - { - grad[0] = -1.0; - grad[1] = -1.0; - grad[2] = -1.0; - } - else if (i == 1) - { - grad[0] = +1.0; - grad[1] = +0.0; - grad[2] = +0.0; - } - else if (i == 2) - { - grad[0] = +0.0; - grad[1] = +1.0; - grad[2] = +0.0; - } - else if (i == 3) - { - grad[0] = +0.0; - grad[1] = +0.0; - grad[2] = +1.0; - } - } - else if (this->degree() == 2) - { - const double r = p[0]; - const double s = p[1]; - const double t = p[2]; - const double u = 1.0 - p[0] - p[1] - p[2]; - - if (i == 0) - { - grad[0] = -4.0 * u + 1.; - grad[1] = grad[0]; - grad[2] = grad[0]; - } - else if (i == 1) - { - grad[0] = +4.0 * r - 1.; - grad[1] = +0.0; - grad[2] = +0.0; - } - else if (i == 2) - { - grad[0] = +0.0; - grad[1] = +4.0 * s - 1.; - grad[2] = +0.0; - } - else if (i == 3) - { - grad[0] = +0.0; - grad[1] = +0.0; - grad[2] = +4.0 * t - 1.; - } - else if (i == 4) - { - grad[0] = +4.0 * (u - r); - grad[1] = -4.0 * r; - grad[2] = -4.0 * r; - } - else if (i == 5) - { - grad[0] = +4.0 * s; - grad[1] = +4.0 * r; - grad[2] = +0.0; - } - else if (i == 6) - { - grad[0] = -4.0 * s; - grad[1] = +4.0 * (u - s); - grad[2] = -4.0 * s; - } - else if (i == 7) - { - grad[0] = -4.0 * t; - grad[1] = -4.0 * t; - grad[2] = +4.0 * (u - t); - } - else if (i == 8) - { - grad[0] = +4.0 * t; - grad[1] = +0.0; - grad[2] = +4.0 * r; - } - else if (i == 9) - { - grad[0] = +0.0; - grad[1] = +4.0 * t; - grad[2] = +4.0 * s; - } - } - else - { - Assert(false, ExcNotImplemented()); - } - } - else - { - Assert(false, ExcNotImplemented()); - } - - return grad; - } - - - - template - Tensor<2, dim> - ScalarPolynomial::compute_grad_grad(const unsigned int i, - const Point & p) const - { - (void)p; - - Tensor<2, dim> result; - - if (this->degree() < 2) - return result; - - if (dim == 1) - { - if (i == 0) - { - result[0][0] = 4.0; - } - else if (i == 1) - { - result[0][0] = 4.0; - } - else if (i == 2) - { - result[0][0] = -8.0; - } - else - { - Assert(false, ExcNotImplemented()); - } - } - else if (dim == 2) - { - if (i == 0) - { - for (unsigned j = 0; j < dim; ++j) - for (unsigned k = 0; k < dim; ++k) - result[j][k] = 4.0; - } - else if (i == 1) - { - result[0][0] = 4.0; - } - else if (i == 2) - { - result[1][1] = 4.0; - } - else if (i == 3) - { - result[0][0] = -8.0; - result[0][1] = result[1][0] = -4.0; - } - else if (i == 4) - { - result[0][1] = result[1][0] = 4.0; - } - else if (i == 5) - { - result[1][1] = -8.0; - result[0][1] = result[1][0] = -4.0; - } - else - { - Assert(false, ExcNotImplemented()); - } - } - else if (dim == 3) - { - if (i == 0) - { - for (unsigned j = 0; j < dim; ++j) - for (unsigned k = 0; k < dim; ++k) - result[j][k] = 4.0; - } - else if (i == 1) - { - result[0][0] = 4.0; - } - else if (i == 2) - { - result[1][1] = 4.0; - } - else if (i == 3) - { - result[2][2] = 4.0; - } - else if (i == 4) - { - result[0][0] = -8.0; - result[0][1] = result[0][2] = result[1][0] = result[2][0] = -4.0; - } - else if (i == 5) - { - result[0][1] = result[1][0] = 4.0; - } - else if (i == 6) - { - result[1][1] = -8.0; - result[0][1] = result[1][0] = result[1][2] = result[2][1] = -4.0; - } - else if (i == 7) - { - result[2][2] = -8.0; - result[2][0] = result[2][1] = result[0][2] = result[1][2] = -4.0; - } - else if (i == 8) - { - result[0][2] = result[2][0] = 4.0; - } - else if (i == 9) - { - result[1][2] = result[2][1] = 4.0; - } - else - { - Assert(false, ExcNotImplemented()); - } - } - else - { - Assert(false, ExcNotImplemented()); - } - - return result; - } - - - - template - void - ScalarPolynomial::evaluate( - const Point & unit_point, - std::vector & values, - std::vector> &grads, - std::vector> &grad_grads, - std::vector> &third_derivatives, - std::vector> &fourth_derivatives) const - { - (void)grads; - (void)grad_grads; - (void)third_derivatives; - (void)fourth_derivatives; - - if (values.size() == this->n()) - for (unsigned int i = 0; i < this->n(); i++) - values[i] = compute_value(i, unit_point); - - if (grads.size() == this->n()) - for (unsigned int i = 0; i < this->n(); i++) - grads[i] = compute_grad(i, unit_point); - - if (grad_grads.size() == this->n()) - for (unsigned int i = 0; i < this->n(); i++) - grad_grads[i] = compute_grad_grad(i, unit_point); - } - - - - template - Tensor<1, dim> - ScalarPolynomial::compute_1st_derivative(const unsigned int i, - const Point & p) const - { - return compute_grad(i, p); - } - - - - template - Tensor<2, dim> - ScalarPolynomial::compute_2nd_derivative(const unsigned int i, - const Point & p) const - { - return compute_grad_grad(i, p); - } - - - - template - Tensor<3, dim> - ScalarPolynomial::compute_3rd_derivative(const unsigned int i, - const Point & p) const - { - (void)i; - (void)p; - - Assert(false, ExcNotImplemented()); - - return {}; - } - - - - template - Tensor<4, dim> - ScalarPolynomial::compute_4th_derivative(const unsigned int i, - const Point & p) const - { - (void)i; - (void)p; - - Assert(false, ExcNotImplemented()); - - return {}; - } - - - - template - std::string - ScalarPolynomial::name() const - { - return "Simplex"; - } - - - - template - std::unique_ptr> - ScalarPolynomial::clone() const - { - return std::make_unique>(*this); - } - - - template ScalarWedgePolynomial::ScalarWedgePolynomial(const unsigned int degree) : ScalarPolynomialsBase(degree, @@ -1076,9 +501,6 @@ namespace Simplex - template class ScalarPolynomial<1>; - template class ScalarPolynomial<2>; - template class ScalarPolynomial<3>; template class ScalarWedgePolynomial<1>; template class ScalarWedgePolynomial<2>; template class ScalarWedgePolynomial<3>; -- 2.39.5