From 3fe5311e2274765e382ce961239ae332a808db93 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Thu, 2 Jul 2020 10:19:03 +0200 Subject: [PATCH] Add Simplex::ScalarPolynomial --- doc/news/changes/minor/20200704Munch | 4 + include/deal.II/simplex/polynomials.h | 167 ++++++++++ source/simplex/CMakeLists.txt | 1 + source/simplex/polynomials.cc | 444 ++++++++++++++++++++++++++ tests/simplex/polynomials_01.cc | 86 +++++ tests/simplex/polynomials_01.output | 29 ++ 6 files changed, 731 insertions(+) create mode 100644 doc/news/changes/minor/20200704Munch create mode 100644 include/deal.II/simplex/polynomials.h create mode 100644 source/simplex/polynomials.cc create mode 100644 tests/simplex/polynomials_01.cc create mode 100644 tests/simplex/polynomials_01.output diff --git a/doc/news/changes/minor/20200704Munch b/doc/news/changes/minor/20200704Munch new file mode 100644 index 0000000000..858dbcb751 --- /dev/null +++ b/doc/news/changes/minor/20200704Munch @@ -0,0 +1,4 @@ +New: The new class Simplex::ScalarPolynomial provides polynomials defined on +simplices. +
+(Peter Munch, 2020/07/02) diff --git a/include/deal.II/simplex/polynomials.h b/include/deal.II/simplex/polynomials.h new file mode 100644 index 0000000000..a99ceb2c31 --- /dev/null +++ b/include/deal.II/simplex/polynomials.h @@ -0,0 +1,167 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_simplex_polynomials_h +#define dealii_simplex_polynomials_h + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +/** + * A namespace for functions and classes that provide support for simplex + * reference cell entities, i.e., triangles and tetrahedrons. + */ +namespace Simplex +{ + /** + * Polynomials defined on dim-dimensional simplex entities. This class is + * basis of Simplex::FE_P. + */ + template + class ScalarPolynomial : public ScalarPolynomialsBase + { + static_assert(dim == 2 || dim == 3, "Dimension not supported!"); + + 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 and @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 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() + * + * @note Not implemented yet. + */ + 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; + }; + + // template functions + template + template + Tensor + ScalarPolynomial::compute_derivative(const unsigned int i, + const Point & p) const + { + Tensor der; + + AssertDimension(order, 1); + const auto grad = compute_grad(i, p); + + for (unsigned int i = 0; i < dim; i++) + der[i] = grad[i]; + + return der; + } +} // namespace Simplex + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/simplex/CMakeLists.txt b/source/simplex/CMakeLists.txt index cfc0183c55..f0c06e1d0d 100644 --- a/source/simplex/CMakeLists.txt +++ b/source/simplex/CMakeLists.txt @@ -16,6 +16,7 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) SET(_unity_include_src + polynomials.cc quadrature_lib.cc ) diff --git a/source/simplex/polynomials.cc b/source/simplex/polynomials.cc new file mode 100644 index 0000000000..c1823cfdff --- /dev/null +++ b/source/simplex/polynomials.cc @@ -0,0 +1,444 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Simplex +{ + namespace + { + unsigned int + compute_n_polynomials(const unsigned int dim, const unsigned int degree) + { + 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; + } + } // namespace + + + + 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 == 2) + { + if (this->degree() == 1) + { + if (i == 0) + return p[0]; + else if (i == 1) + return p[1]; + else if (i == 2) + return 1.0 - p[0] - p[1]; + } + else if (this->degree() == 2) + { + const double t1 = p[0]; + const double t2 = p[1]; + const double t3 = 1.0 - p[0] - 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 == 2) + { + if (this->degree() == 1) + { + if (i == 0) + { + grad[0] = +1.0; + grad[1] = +0.0; + } + else if (i == 1) + { + grad[0] = +0.0; + grad[1] = +1.0; + } + else if (i == 2) + { + grad[0] = -1.0; + grad[1] = -1.0; + } + else + { + Assert(false, ExcNotImplemented()); + } + } + else if (this->degree() == 2) + { + if (i == 2) + { + grad[0] = -3.0 + 4.0 * (p[0] + p[1]); + grad[1] = -3.0 + 4.0 * (p[0] + p[1]); + } + else if (i == 0) + { + grad[0] = 4.0 * p[0] - 1.0; + grad[1] = 0.0; + } + else if (i == 1) + { + grad[0] = 0.0; + grad[1] = 4.0 * p[1] - 1.0; + } + else if (i == 5) + { + grad[0] = 4.0 * (1.0 - 2.0 * p[0] - p[1]); + grad[1] = -4.0 * p[0]; + } + else if (i == 3) + { + grad[0] = 4.0 * p[1]; + grad[1] = 4.0 * p[0]; + } + else if (i == 4) + { + 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)i; + (void)p; + + Assert(false, ExcNotImplemented()); + return Tensor<2, dim>(); + } + + + + 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); + } + + + + 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 + { + (void)i; + (void)p; + + Assert(false, ExcNotImplemented()); + + return {}; + } + + + + 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 class ScalarPolynomial<2>; + template class ScalarPolynomial<3>; + +} // namespace Simplex + +DEAL_II_NAMESPACE_CLOSE diff --git a/tests/simplex/polynomials_01.cc b/tests/simplex/polynomials_01.cc new file mode 100644 index 0000000000..86c276ef9e --- /dev/null +++ b/tests/simplex/polynomials_01.cc @@ -0,0 +1,86 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test Simplex::ScalarPolynomial on an the points of an arbitrary quadrature +// rule. + + +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template +void +test(const unsigned int degree) +{ + Simplex::ScalarPolynomial poly(degree); + QSimplex quad(QGauss(degree + 1)); + + std::vector values(poly.n()); + std::vector> grads(poly.n()); + std::vector> grad_grads; + std::vector> third_derivatives; + std::vector> fourth_derivatives; + + for (unsigned int i = 0; i < quad.size(); i++) + { + poly.evaluate(quad.point(i), + values, + grads, + grad_grads, + third_derivatives, + fourth_derivatives); + + for (auto v : values) + deallog << v << " "; + deallog << std::endl; + + for (auto v : grads) + deallog << v << " "; + deallog << std::endl; + } +} + +int +main() +{ + initlog(); + + { + deallog.push("2d-1"); + test<2>(1); + deallog.pop(); + } + { + deallog.push("2d-2"); + test<2>(2); + deallog.pop(); + } + { + deallog.push("3d-1"); + test<3>(1); + deallog.pop(); + } + { + deallog.push("3d-2"); + test<3>(2); + deallog.pop(); + } +} diff --git a/tests/simplex/polynomials_01.output b/tests/simplex/polynomials_01.output new file mode 100644 index 0000000000..fe28ac73ea --- /dev/null +++ b/tests/simplex/polynomials_01.output @@ -0,0 +1,29 @@ + +DEAL:2d-1::0.211325 0.211325 0.577350 +DEAL:2d-1::1.00000 0.00000 0.00000 1.00000 -1.00000 -1.00000 +DEAL:2d-1::0.788675 0.211325 2.77556e-17 +DEAL:2d-1::1.00000 0.00000 0.00000 1.00000 -1.00000 -1.00000 +DEAL:2d-1::0.211325 0.788675 0.00000 +DEAL:2d-1::1.00000 0.00000 0.00000 1.00000 -1.00000 -1.00000 +DEAL:2d-2::-0.0872983 -0.0872983 0.425403 0.0508067 0.349193 0.349193 +DEAL:2d-2::-0.549193 0.00000 0.00000 -0.549193 -2.09839 -2.09839 0.450807 0.450807 -0.450807 2.64758 2.64758 -0.450807 +DEAL:2d-2::0.00000 -0.0872983 -0.0872983 0.225403 0.174597 0.774597 +DEAL:2d-2::1.00000 0.00000 0.00000 -0.549193 -0.549193 -0.549193 0.450807 2.00000 -0.450807 1.09839 -0.450807 -2.00000 +DEAL:2d-2::0.687298 -0.0872983 1.38778e-17 0.400000 -6.25620e-18 -4.92550e-17 +DEAL:2d-2::2.54919 0.00000 0.00000 -0.549193 1.00000 1.00000 0.450807 3.54919 -0.450807 -0.450807 -3.54919 -3.54919 +DEAL:2d-2::-0.0872983 0.00000 -0.0872983 0.225403 0.774597 0.174597 +DEAL:2d-2::-0.549193 0.00000 0.00000 1.00000 -0.549193 -0.549193 2.00000 0.450807 -2.00000 -0.450807 1.09839 -0.450807 +DEAL:2d-2::0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 +DEAL:2d-2::1.00000 0.00000 0.00000 1.00000 1.00000 1.00000 2.00000 2.00000 -2.00000 -2.00000 -2.00000 -2.00000 +DEAL:2d-2::-0.0872983 0.687298 0.00000 0.400000 0.00000 0.00000 +DEAL:2d-2::-0.549193 0.00000 0.00000 2.54919 1.00000 1.00000 3.54919 0.450807 -3.54919 -3.54919 -0.450807 -0.450807 +DEAL:3d-1::0.366025 0.211325 0.211325 0.211325 +DEAL:3d-1::-1.00000 -1.00000 -1.00000 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000 +DEAL:3d-2::0.214315 -0.0872983 -0.0872983 -0.0872983 0.298387 0.0508067 0.298387 0.298387 0.0508067 0.0508067 +DEAL:3d-2::-1.64758 -1.64758 -1.64758 -0.549193 0.00000 0.00000 0.00000 -0.549193 0.00000 0.00000 0.00000 -0.549193 2.19677 -0.450807 -0.450807 0.450807 0.450807 0.00000 -0.450807 2.19677 -0.450807 -0.450807 -0.450807 2.19677 0.450807 0.00000 0.450807 0.00000 0.450807 0.450807 +DEAL:3d-2::-0.123790 0.00000 -0.0872983 -0.0872983 0.549193 0.225403 0.123790 0.123790 0.225403 0.0508067 +DEAL:3d-2::-0.0983867 -0.0983867 -0.0983867 1.00000 0.00000 0.00000 0.00000 -0.549193 0.00000 0.00000 0.00000 -0.549193 -0.901613 -2.00000 -2.00000 0.450807 2.00000 0.00000 -0.450807 0.647580 -0.450807 -0.450807 -0.450807 0.647580 0.450807 0.00000 2.00000 0.00000 0.450807 0.450807 +DEAL:3d-2::-0.123790 -0.0872983 0.00000 -0.0872983 0.123790 0.225403 0.549193 0.123790 0.0508067 0.225403 +DEAL:3d-2::-0.0983867 -0.0983867 -0.0983867 -0.549193 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 -0.549193 0.647580 -0.450807 -0.450807 2.00000 0.450807 0.00000 -2.00000 -0.901613 -2.00000 -0.450807 -0.450807 0.647580 0.450807 0.00000 0.450807 0.00000 0.450807 2.00000 +DEAL:3d-2::-0.123790 -0.0872983 -0.0872983 0.00000 0.123790 0.0508067 0.123790 0.549193 0.225403 0.225403 +DEAL:3d-2::-0.0983867 -0.0983867 -0.0983867 -0.549193 0.00000 0.00000 0.00000 -0.549193 0.00000 0.00000 0.00000 1.00000 0.647580 -0.450807 -0.450807 0.450807 0.450807 0.00000 -0.450807 0.647580 -0.450807 -2.00000 -2.00000 -0.901613 2.00000 0.00000 0.450807 0.00000 2.00000 0.450807 -- 2.39.5