--- /dev/null
+New: The new class Simplex::ScalarPolynomial provides polynomials defined on
+simplices.
+<br>
+(Peter Munch, 2020/07/02)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+
+#include <deal.II/base/scalar_polynomials_base.h>
+
+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 <int dim>
+ class ScalarPolynomial : public ScalarPolynomialsBase<dim>
+ {
+ 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<dim> & unit_point,
+ std::vector<double> & values,
+ std::vector<Tensor<1, dim>> &grads,
+ std::vector<Tensor<2, dim>> &grad_grads,
+ std::vector<Tensor<3, dim>> &third_derivatives,
+ std::vector<Tensor<4, dim>> &fourth_derivatives) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_value()
+ */
+ double
+ compute_value(const unsigned int i, const Point<dim> &p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_derivative()
+ *
+ * @note Currently, only implemented for first derivative.
+ */
+ template <int order>
+ Tensor<order, dim>
+ compute_derivative(const unsigned int i, const Point<dim> &p) const;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_1st_derivative()
+ */
+ Tensor<1, dim>
+ compute_1st_derivative(const unsigned int i,
+ const Point<dim> & p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_2nd_derivative()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<2, dim>
+ compute_2nd_derivative(const unsigned int i,
+ const Point<dim> & p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_3rd_derivative()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<3, dim>
+ compute_3rd_derivative(const unsigned int i,
+ const Point<dim> & p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_4th_derivative()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<4, dim>
+ compute_4th_derivative(const unsigned int i,
+ const Point<dim> & p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_grad()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<1, dim>
+ compute_grad(const unsigned int i, const Point<dim> &p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::compute_grad_grad()
+ *
+ * @note Not implemented yet.
+ */
+ Tensor<2, dim>
+ compute_grad_grad(const unsigned int i, const Point<dim> &p) const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::name()
+ */
+ std::string
+ name() const override;
+
+ /**
+ * @copydoc ScalarPolynomialsBase::clone()
+ */
+ virtual std::unique_ptr<ScalarPolynomialsBase<dim>>
+ clone() const override;
+ };
+
+ // template functions
+ template <int dim>
+ template <int order>
+ Tensor<order, dim>
+ ScalarPolynomial<dim>::compute_derivative(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ Tensor<order, dim> 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
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
SET(_unity_include_src
+ polynomials.cc
quadrature_lib.cc
)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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/simplex/polynomials.h>
+
+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 <int dim>
+ ScalarPolynomial<dim>::ScalarPolynomial(const unsigned int degree)
+ : ScalarPolynomialsBase<dim>(degree, compute_n_polynomials(dim, degree))
+ {}
+
+
+
+ template <int dim>
+ double
+ ScalarPolynomial<dim>::compute_value(const unsigned int i,
+ const Point<dim> & 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 <int dim>
+ Tensor<1, dim>
+ ScalarPolynomial<dim>::compute_grad(const unsigned int i,
+ const Point<dim> & 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 <int dim>
+ Tensor<2, dim>
+ ScalarPolynomial<dim>::compute_grad_grad(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ (void)i;
+ (void)p;
+
+ Assert(false, ExcNotImplemented());
+ return Tensor<2, dim>();
+ }
+
+
+
+ template <int dim>
+ void
+ ScalarPolynomial<dim>::evaluate(
+ const Point<dim> & unit_point,
+ std::vector<double> & values,
+ std::vector<Tensor<1, dim>> &grads,
+ std::vector<Tensor<2, dim>> &grad_grads,
+ std::vector<Tensor<3, dim>> &third_derivatives,
+ std::vector<Tensor<4, dim>> &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 <int dim>
+ Tensor<1, dim>
+ ScalarPolynomial<dim>::compute_1st_derivative(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ return compute_grad(i, p);
+ }
+
+
+
+ template <int dim>
+ Tensor<2, dim>
+ ScalarPolynomial<dim>::compute_2nd_derivative(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ (void)i;
+ (void)p;
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+
+
+ template <int dim>
+ Tensor<3, dim>
+ ScalarPolynomial<dim>::compute_3rd_derivative(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ (void)i;
+ (void)p;
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+
+
+ template <int dim>
+ Tensor<4, dim>
+ ScalarPolynomial<dim>::compute_4th_derivative(const unsigned int i,
+ const Point<dim> & p) const
+ {
+ (void)i;
+ (void)p;
+
+ Assert(false, ExcNotImplemented());
+
+ return {};
+ }
+
+
+
+ template <int dim>
+ std::string
+ ScalarPolynomial<dim>::name() const
+ {
+ return "Simplex";
+ }
+
+
+
+ template <int dim>
+ std::unique_ptr<ScalarPolynomialsBase<dim>>
+ ScalarPolynomial<dim>::clone() const
+ {
+ return std::make_unique<ScalarPolynomial<dim>>(*this);
+ }
+
+ template class ScalarPolynomial<2>;
+ template class ScalarPolynomial<3>;
+
+} // namespace Simplex
+
+DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/simplex/polynomials.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+void
+test(const unsigned int degree)
+{
+ Simplex::ScalarPolynomial<dim> poly(degree);
+ QSimplex<dim> quad(QGauss<dim>(degree + 1));
+
+ std::vector<double> values(poly.n());
+ std::vector<Tensor<1, dim>> grads(poly.n());
+ std::vector<Tensor<2, dim>> grad_grads;
+ std::vector<Tensor<3, dim>> third_derivatives;
+ std::vector<Tensor<4, dim>> 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();
+ }
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff