From c53512359bc3b24c8daea108844dd3640dc96a26 Mon Sep 17 00:00:00 2001 From: David Wells Date: Thu, 14 Jan 2021 16:52:56 -0500 Subject: [PATCH] Implement Simplex::BarycentricPolynomials. This gives us a much more convenient way to set up polynomials on simplices. --- .../deal.II/simplex/barycentric_polynomials.h | 703 ++++++++++++++++++ source/simplex/CMakeLists.txt | 1 + source/simplex/barycentric_polynomials.cc | 341 +++++++++ source/simplex/fe_lib.cc | 1 + tests/simplex/barycentric_01.cc | 170 +++++ tests/simplex/barycentric_01.output | 51 ++ 6 files changed, 1267 insertions(+) create mode 100644 include/deal.II/simplex/barycentric_polynomials.h create mode 100644 source/simplex/barycentric_polynomials.cc create mode 100644 tests/simplex/barycentric_01.cc create mode 100644 tests/simplex/barycentric_01.output diff --git a/include/deal.II/simplex/barycentric_polynomials.h b/include/deal.II/simplex/barycentric_polynomials.h new file mode 100644 index 0000000000..6fe862635d --- /dev/null +++ b/include/deal.II/simplex/barycentric_polynomials.h @@ -0,0 +1,703 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 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_barycentric_polynomials_h +#define dealii_simplex_barycentric_polynomials_h + +#include + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Simplex +{ + /** + * Polynomial implemented in barycentric coordinates. + * + * Barycentric coordinates are a coordinate system defined on simplices that + * are particularly easy to work with since they express coordinates in the + * simplex as convex combinations of the vertices. For example, any point in a + * triangle can be written as + * + * @f[ + * (x, y) = c_0 (x_0, y_0) + c_1 (x_1, y_1) + c_2 (x_2, y_2). + * @f] + * + * where each value $c_i$ is the relative weight of each vertex (so the + * centroid is, in 2D, where each $c_i = 1/3$). Since we only consider convex + * combinations we can rewrite this equation as + * + * @f[ + * (x, y) = (1 - c_1 - c_2) (x_0, y_0) + c_1 (x_1, y_1) + c_2 (x_2, y_2). + * @f] + * + * This results in three polynomials that are equivalent to $P^1$ in 2D. More + * exactly, this class implements a polynomial space defined with the basis, + * in 2D, of + * @f{align*}{ + * t_0(x, y) &= 1 - x - y \\ + * t_1(x, y) &= x \\ + * t_2(x, y) &= y + * @f} + * and, in 3D, + * @f{align*}{ + * t_0(x, y) &= 1 - x - y - z \\ + * t_1(x, y) &= x \\ + * t_2(x, y) &= y \\ + * t_2(x, y) &= z + * @f} + * + * which is, in practice, a very convenient basis for defining simplex + * polynomials: for example, the fourth basis function of a TRI6 element is + * + * @f[ + * 4 * t_1(x, y) * t_2(x, y). + * @f] + * + * Barycentric polynomials in dim-dimensional space have + * dim + 1 variables in since t_0 can be written in + * terms of the other monomials. + * + * Monomials can be conveniently constructed with + * BarycentricPolynomial::monomial(). + * + * @ingroup Polynomials + */ + template + class BarycentricPolynomial + { + public: + /** + * Constructor for the zero polynomial. + */ + BarycentricPolynomial(); + + /** + * Constructor for a monomial. + */ + BarycentricPolynomial(const TableIndices &powers, + const Number coefficient); + + /** + * Return the specified monomial. + */ + static BarycentricPolynomial + monomial(const unsigned int d); + + /** + * Print the polynomial to the output stream with lowest-order terms first. + * For example, the first P6 basis function is printed as + * -1 * t0^1 + 2 * t0^2, where t0 is the first + * barycentric variable, t1 is the second, etc. + */ + void + print(std::ostream &out) const; + + /** + * Degree of each barycentric polynomial. + */ + TableIndices + degrees() const; + + /** + * Unary minus. + */ + BarycentricPolynomial + operator-() const; + + /** + * Add a scalar. + */ + template + BarycentricPolynomial + operator+(const Number2 &a) const; + + /** + * Subtract a scalar. + */ + template + BarycentricPolynomial + operator-(const Number2 &a) const; + + /** + * Multiply by a scalar. + */ + template + BarycentricPolynomial operator*(const Number2 &a) const; + + /** + * Divide by a scalar. + */ + template + BarycentricPolynomial + operator/(const Number2 &a) const; + + /** + * Add another barycentric polynomial. + */ + BarycentricPolynomial + operator+(const BarycentricPolynomial &augend) const; + + /** + * Subtract another barycentric polynomial. + */ + BarycentricPolynomial + operator-(const BarycentricPolynomial &augend) const; + + /** + * Multiply by another barycentric polynomial. + */ + BarycentricPolynomial + operator*(const BarycentricPolynomial &multiplicand) const; + + /** + * Differentiate in barycentric coordinates. + */ + BarycentricPolynomial + barycentric_derivative(const unsigned int coordinate) const; + + /** + * Differentiate in Cartesian coordinates. + */ + BarycentricPolynomial + derivative(const unsigned int coordinate) const; + + /** + * Evaluate the polynomial. + */ + Number + value(const Point &point) const; + + /** + * Return an estimate, in bytes, of the memory usage of the object. + */ + std::size_t + memory_consumption() const; + + protected: + /** + * Coefficients of the polynomial. The exponents are the integer indexes. + */ + Table coefficients; + + /** + * Utility function for barycentric polynomials - its convenient to loop + * over all the indices at once in a dimension-independent way, but we also + * need to access the actual indices of the underlying Table object. This + * utility function converts an integral index into the equivalent + * TableIndices array (which are also the implicitly stored polynomial + * exponents). + */ + static TableIndices + index_to_indices(const std::size_t & index, + const TableIndices &extent); + }; + + /** + * Scalar polynomial space based on barycentric polynomials. + */ + template + class BarycentricPolynomials : public ScalarPolynomialsBase + { + public: + /** + * Make the dimension available to the outside. + */ + static const unsigned int dimension = dim; + + /** + * Get the standard Lagrange basis for a specified degree. + */ + static BarycentricPolynomials + get_fe_p_basis(const unsigned int degree); + + /* + * Constructor taking the polynomial @p degree as input. + */ + BarycentricPolynomials( + const std::vector> &polynomials); + + /** + * @copydoc ScalarPolynomialsBase::evaluate() + */ + 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_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() + */ + Tensor<3, dim> + compute_3rd_derivative(const unsigned int i, + const Point & p) const override; + + /** + * @copydoc ScalarPolynomialsBase::compute_4th_derivative() + */ + Tensor<4, dim> + compute_4th_derivative(const unsigned int i, + const Point & p) const override; + + /** + * @copydoc ScalarPolynomialsBase::compute_grad() + */ + Tensor<1, dim> + compute_grad(const unsigned int i, const Point &p) const override; + + /** + * @copydoc ScalarPolynomialsBase::compute_grad_grad() + */ + Tensor<2, dim> + compute_grad_grad(const unsigned int i, const Point &p) const override; + + /** + * @copydoc ScalarPolynomialsBase::memory_consumption() + */ + virtual std::size_t + memory_consumption() const override; + + /** + * @copydoc ScalarPolynomialsBase::name() + */ + std::string + name() const override; + + /** + * @copydoc ScalarPolynomialsBase::clone() + */ + virtual std::unique_ptr> + clone() const override; + + protected: + std::vector> polys; + + Table<2, BarycentricPolynomial> poly_grads; + + Table<3, BarycentricPolynomial> poly_hessians; + + Table<4, BarycentricPolynomial> poly_third_derivatives; + + Table<5, BarycentricPolynomial> poly_fourth_derivatives; + }; + + // non-member template functions for algebra + + /** + * Multiply a Simplex::BarycentricPolynomial by a constant. + */ + template + BarycentricPolynomial + operator*(const Number2 &a, const BarycentricPolynomial &bp) + { + return bp * Number1(a); + } + + /** + * Add a constant to a Simplex::BarycentricPolynomial. + */ + template + BarycentricPolynomial + operator+(const Number2 &a, const BarycentricPolynomial &bp) + { + return bp + Number1(a); + } + + /** + * Subtract a Simplex::BarycentricPolynomial from a constant. + */ + template + BarycentricPolynomial + operator-(const Number2 &a, const BarycentricPolynomial &bp) + { + return bp - Number1(a); + } + + /** + * Write a Simplex::BarycentricPolynomial to the provided output stream. + */ + template + std::ostream & + operator<<(std::ostream &out, const BarycentricPolynomial &bp) + { + bp.print(out); + return out; + } +} // namespace Simplex + +// Template function definitions + +namespace Simplex +{ + // BarycentricPolynomial: + template + BarycentricPolynomial::BarycentricPolynomial() + { + TableIndices extents; + for (unsigned int d = 0; d < dim + 1; ++d) + extents[d] = 1; + coefficients.reinit(extents); + + coefficients(TableIndices{}) = Number(); + } + + + + template + BarycentricPolynomial::BarycentricPolynomial( + const TableIndices &powers, + const Number coefficient) + { + TableIndices extents; + for (unsigned int d = 0; d < dim + 1; ++d) + extents[d] = powers[d] + 1; + coefficients.reinit(extents); + + coefficients(powers) = coefficient; + } + + + + template + BarycentricPolynomial + BarycentricPolynomial::monomial(const unsigned int d) + { + AssertIndexRange(d, dim + 1); + TableIndices indices; + indices[d] = 1; + return BarycentricPolynomial(indices, Number(1)); + } + + + + template + void + BarycentricPolynomial::print(std::ostream &out) const + { + const auto &coeffs = this->coefficients; + auto first = index_to_indices(0, coeffs.size()); + bool print_plus = false; + if (coeffs(first) != Number()) + { + out << coeffs(first); + print_plus = true; + } + for (std::size_t i = 1; i < coeffs.n_elements(); ++i) + { + const auto indices = index_to_indices(i, coeffs.size()); + if (coeffs(indices) == Number()) + continue; + if (print_plus) + out << " + "; + out << coeffs(indices); + for (unsigned int d = 0; d < dim + 1; ++d) + { + if (indices[d] != 0) + out << " * t" << d << '^' << indices[d]; + } + print_plus = true; + } + + if (!print_plus) + out << Number(); + } + + + + template + TableIndices + BarycentricPolynomial::degrees() const + { + auto deg = coefficients.size(); + for (unsigned int d = 0; d < dim + 1; ++d) + deg[d] -= 1; + return deg; + } + + + + template + BarycentricPolynomial + BarycentricPolynomial::operator-() const + { + return *this * Number(-1); + } + + + + template + template + BarycentricPolynomial + BarycentricPolynomial::operator+(const Number2 &a) const + { + BarycentricPolynomial result(*this); + result.coefficients(index_to_indices(0, result.coefficients.size())) += a; + + return result; + } + + + + template + template + BarycentricPolynomial + BarycentricPolynomial::operator-(const Number2 &a) const + { + return *this + (-a); + } + + + + template + template + BarycentricPolynomial BarycentricPolynomial:: + operator*(const Number2 &a) const + { + if (a == Number2()) + { + return BarycentricPolynomial(); + } + + BarycentricPolynomial result(*this); + for (std::size_t i = 0; i < result.coefficients.n_elements(); ++i) + { + const auto index = index_to_indices(i, result.coefficients.size()); + result.coefficients(index) *= a; + } + + return result; + } + + + + template + template + BarycentricPolynomial + BarycentricPolynomial::operator/(const Number2 &a) const + { + Assert(a != Number2(), ExcDivideByZero()); + return *this * (Number(1) / Number(a)); + } + + + + template + BarycentricPolynomial + BarycentricPolynomial:: + operator+(const BarycentricPolynomial &augend) const + { + TableIndices deg; + for (unsigned int d = 0; d < dim + 1; ++d) + { + deg[d] = std::max(degrees()[d], augend.degrees()[d]); + } + + BarycentricPolynomial result(deg, Number()); + + auto add_coefficients = [&](const Table &in) { + for (std::size_t i = 0; i < in.n_elements(); ++i) + { + const auto index = index_to_indices(i, in.size()); + result.coefficients(index) += in(index); + } + }; + + add_coefficients(this->coefficients); + add_coefficients(augend.coefficients); + return result; + } + + + + template + BarycentricPolynomial + BarycentricPolynomial:: + operator-(const BarycentricPolynomial &augend) const + { + return *this + (-augend); + } + + + + template + BarycentricPolynomial BarycentricPolynomial:: + operator*(const BarycentricPolynomial &multiplicand) const + { + TableIndices deg; + for (unsigned int d = 0; d < dim + 1; ++d) + { + deg[d] = multiplicand.degrees()[d] + degrees()[d]; + } + + BarycentricPolynomial result(deg, Number()); + + const auto &coef_1 = this->coefficients; + const auto &coef_2 = multiplicand.coefficients; + auto & coef_out = result.coefficients; + + for (std::size_t i1 = 0; i1 < coef_1.n_elements(); ++i1) + { + const auto index_1 = index_to_indices(i1, coef_1.size()); + for (std::size_t i2 = 0; i2 < coef_2.n_elements(); ++i2) + { + const auto index_2 = index_to_indices(i2, coef_2.size()); + + TableIndices index_out; + for (unsigned int d = 0; d < dim + 1; ++d) + index_out[d] = index_1[d] + index_2[d]; + coef_out(index_out) += coef_1(index_1) * coef_2(index_2); + } + } + + return result; + } + + + + template + BarycentricPolynomial + BarycentricPolynomial::barycentric_derivative( + const unsigned int coordinate) const + { + AssertIndexRange(coordinate, dim + 1); + + if (degrees()[coordinate] == 0) + return BarycentricPolynomial(); + + auto deg = degrees(); + deg[coordinate] -= 1; + BarycentricPolynomial result( + deg, std::numeric_limits::max()); + const auto &coeffs_in = coefficients; + auto & coeffs_out = result.coefficients; + for (std::size_t i = 0; i < coeffs_out.n_elements(); ++i) + { + const auto out_index = index_to_indices(i, coeffs_out.size()); + auto input_index = out_index; + input_index[coordinate] += 1; + + coeffs_out(out_index) = + coeffs_in(input_index) * input_index[coordinate]; + } + + return result; + } + + + + template + BarycentricPolynomial + BarycentricPolynomial::derivative( + const unsigned int coordinate) const + { + AssertIndexRange(coordinate, dim); + return -barycentric_derivative(0) + barycentric_derivative(coordinate + 1); + } + + + + template + Number + BarycentricPolynomial::value(const Point &point) const + { + // TODO: this is probably not numerically stable for higher order. + // We really need some version of Horner's method. + Number result = {}; + + // Begin by converting point (which is in Cartesian coordinates) to + // barycentric coordinates: + std::array b_point; + b_point[0] = 1.0; + for (unsigned int d = 0; d < dim; ++d) + { + b_point[0] -= point[d]; + b_point[d + 1] = point[d]; + } + + // Now evaluate the polynomial at the computed barycentric point: + for (std::size_t i = 0; i < coefficients.n_elements(); ++i) + { + const auto indices = index_to_indices(i, coefficients.size()); + const auto coef = coefficients(indices); + if (coef == Number()) + continue; + + auto temp = Number(1); + for (unsigned int d = 0; d < dim + 1; ++d) + temp *= std::pow(b_point[d], indices[d]); + result += coef * temp; + } + + return result; + } + + template + std::size_t + BarycentricPolynomial::memory_consumption() const + { + return coefficients.memory_consumption(); + } + + template + TableIndices + BarycentricPolynomial::index_to_indices( + const std::size_t & index, + const TableIndices &extent) + { + TableIndices result; + auto temp = index; + + for (unsigned int n = 0; n < dim + 1; ++n) + { + std::size_t slice_size = 1; + for (unsigned int n2 = n + 1; n2 < dim + 1; ++n2) + slice_size *= extent[n2]; + result[n] = temp / slice_size; + temp %= slice_size; + } + return result; + } +} // namespace Simplex + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/simplex/CMakeLists.txt b/source/simplex/CMakeLists.txt index f547691c5a..398736607d 100644 --- a/source/simplex/CMakeLists.txt +++ b/source/simplex/CMakeLists.txt @@ -18,6 +18,7 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) SET(_unity_include_src fe_lib.cc grid_generator.cc + barycentric_polynomials.cc polynomials.cc quadrature_lib.cc ) diff --git a/source/simplex/barycentric_polynomials.cc b/source/simplex/barycentric_polynomials.cc new file mode 100644 index 0000000000..8e4e5416a5 --- /dev/null +++ b/source/simplex/barycentric_polynomials.cc @@ -0,0 +1,341 @@ +// --------------------------------------------------------------------- +// +// 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 internal + { + /** + * Get the highest degree of the barycentric polynomial (in Cartesian + * coordinates). + */ + template + unsigned int + get_degree(const std::vector> &polys) + { + // Since the first variable in a simplex polynomial is, e.g., in 2D, + // + // t0 = 1 - x - y + // + // (that is, it depends on the Cartesian variables), we have to compute + // its degree separately. An example: t0*t1*t2 has degree 1 in the affine + // polynomial basis but is degree 2 in the Cartesian polynomial basis. + std::size_t max_degree = 0; + for (const auto &poly : polys) + { + const TableIndices degrees = poly.degrees(); + + const auto degree_0 = degrees[0]; + std::size_t degree_d = 0; + for (unsigned int d = 1; d < dim + 1; ++d) + degree_d = std::max(degree_d, degrees[d]); + + max_degree = std::max(max_degree, degree_d + degree_0); + } + + return max_degree; + } + } // namespace internal + + + template + BarycentricPolynomials + BarycentricPolynomials::get_fe_p_basis(const unsigned int degree) + { + std::vector> polys; + + auto M = [](const unsigned int d) { + return BarycentricPolynomial::monomial(d); + }; + switch (degree) + { + case 0: + polys.push_back(0 * M(0) + 1); + break; + case 1: + { + for (unsigned int d = 0; d < dim + 1; ++d) + polys.push_back(M(d)); + break; + } + case 2: + { + for (unsigned int d = 0; d < dim + 1; ++d) + polys.push_back(M(d) * (2 * M(d) - 1)); + polys.push_back(4 * M(1) * M(0)); + if (dim >= 2) + { + polys.push_back(4 * M(1) * M(2)); + polys.push_back(4 * M(2) * M(0)); + } + if (dim == 3) + { + polys.push_back(4 * M(3) * M(0)); + polys.push_back(4 * M(1) * M(3)); + polys.push_back(4 * M(2) * M(3)); + } + break; + } + default: + Assert(false, ExcNotImplemented()); + } + + return BarycentricPolynomials(polys); + } + + + + template + BarycentricPolynomials::BarycentricPolynomials( + const std::vector> &polynomials) + : ScalarPolynomialsBase(internal::get_degree(polynomials), + polynomials.size()) + { + polys = polynomials; + + poly_grads.reinit({polynomials.size(), dim}); + poly_hessians.reinit({polynomials.size(), dim, dim}); + poly_third_derivatives.reinit({polynomials.size(), dim, dim, dim}); + poly_fourth_derivatives.reinit({polynomials.size(), dim, dim, dim, dim}); + + for (std::size_t i = 0; i < polynomials.size(); ++i) + { + // gradients + for (unsigned int d = 0; d < dim; ++d) + poly_grads[i][d] = polynomials[i].derivative(d); + + // hessians + for (unsigned int d0 = 0; d0 < dim; ++d0) + for (unsigned int d1 = 0; d1 < dim; ++d1) + poly_hessians[i][d0][d1] = poly_grads[i][d0].derivative(d1); + + // third derivatives + for (unsigned int d0 = 0; d0 < dim; ++d0) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + poly_third_derivatives[i][d0][d1][d2] = + poly_hessians[i][d0][d1].derivative(d2); + + // fourth derivatives + for (unsigned int d0 = 0; d0 < dim; ++d0) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + for (unsigned int d3 = 0; d3 < dim; ++d3) + poly_fourth_derivatives[i][d0][d1][d2][d3] = + poly_third_derivatives[i][d0][d1][d2].derivative(d3); + } + } + + + + template + void + BarycentricPolynomials::evaluate( + const Point & unit_point, + std::vector & values, + std::vector> &grads, + std::vector> &grad_grads, + std::vector> &third_derivatives, + std::vector> &fourth_derivatives) const + { + Assert(values.size() == this->n() || values.size() == 0, + ExcDimensionMismatch2(values.size(), this->n(), 0)); + Assert(grads.size() == this->n() || grads.size() == 0, + ExcDimensionMismatch2(grads.size(), this->n(), 0)); + Assert(grad_grads.size() == this->n() || grad_grads.size() == 0, + ExcDimensionMismatch2(grad_grads.size(), this->n(), 0)); + Assert(third_derivatives.size() == this->n() || + third_derivatives.size() == 0, + ExcDimensionMismatch2(third_derivatives.size(), this->n(), 0)); + Assert(fourth_derivatives.size() == this->n() || + fourth_derivatives.size() == 0, + ExcDimensionMismatch2(fourth_derivatives.size(), this->n(), 0)); + + for (std::size_t i = 0; i < polys.size(); ++i) + { + if (values.size() == this->n()) + values[i] = polys[i].value(unit_point); + + // gradients + if (grads.size() == this->n()) + for (unsigned int d = 0; d < dim; ++d) + grads[i][d] = poly_grads[i][d].value(unit_point); + + // hessians + if (grad_grads.size() == this->n()) + for (unsigned int d0 = 0; d0 < dim; ++d0) + for (unsigned int d1 = 0; d1 < dim; ++d1) + grad_grads[i][d0][d1] = + poly_hessians[i][d0][d1].value(unit_point); + + // third derivatives + if (third_derivatives.size() == this->n()) + for (unsigned int d0 = 0; d0 < dim; ++d0) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + third_derivatives[i][d0][d1][d2] = + poly_third_derivatives[i][d0][d1][d2].value(unit_point); + + // fourth derivatives + if (fourth_derivatives.size() == this->n()) + for (unsigned int d0 = 0; d0 < dim; ++d0) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + for (unsigned int d3 = 0; d3 < dim; ++d3) + fourth_derivatives[i][d0][d1][d2][d3] = + poly_fourth_derivatives[i][d0][d1][d2][d3].value( + unit_point); + } + } + + + + template + double + BarycentricPolynomials::compute_value(const unsigned int i, + const Point & p) const + { + AssertIndexRange(i, this->n()); + return polys[i].value(p); + } + + + + template + Tensor<1, dim> + BarycentricPolynomials::compute_1st_derivative(const unsigned int i, + const Point &p) const + { + Tensor<1, dim> result; + for (unsigned int d = 0; d < dim; ++d) + result[d] = poly_grads[i][d].value(p); + return result; + } + + + + template + Tensor<2, dim> + BarycentricPolynomials::compute_2nd_derivative(const unsigned int i, + const Point &p) const + { + Tensor<2, dim> result; + for (unsigned int d0 = 0; d0 < dim; ++d0) + for (unsigned int d1 = 0; d1 < dim; ++d1) + result[d0][d1] = poly_hessians[i][d0][d1].value(p); + + return result; + } + + + + template + Tensor<3, dim> + BarycentricPolynomials::compute_3rd_derivative(const unsigned int i, + const Point &p) const + { + Tensor<3, dim> result; + for (unsigned int d0 = 0; d0 < dim; ++d0) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + result[d0][d1][d2] = poly_third_derivatives[i][d0][d1][d2].value(p); + + return result; + } + + + + template + Tensor<4, dim> + BarycentricPolynomials::compute_4th_derivative(const unsigned int i, + const Point &p) const + { + Tensor<4, dim> result; + for (unsigned int d0 = 0; d0 < dim; ++d0) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + for (unsigned int d3 = 0; d3 < dim; ++d3) + result[d0][d1][d2][d3] = + poly_fourth_derivatives[i][d0][d1][d2][d3].value(p); + + return result; + } + + + + template + Tensor<1, dim> + BarycentricPolynomials::compute_grad(const unsigned int i, + const Point & p) const + { + return compute_1st_derivative(i, p); + } + + + + template + Tensor<2, dim> + BarycentricPolynomials::compute_grad_grad(const unsigned int i, + const Point & p) const + { + return compute_2nd_derivative(i, p); + } + + + + template + std::unique_ptr> + BarycentricPolynomials::clone() const + { + return std::make_unique>(*this); + } + + + + template + std::string + BarycentricPolynomials::name() const + { + return "BarycentricPolynomials<" + std::to_string(dim) + ">"; + } + + + + template + std::size_t + BarycentricPolynomials::memory_consumption() const + { + std::size_t poly_memory = 0; + for (const auto &poly : polys) + poly_memory += poly.memory_consumption(); + return ScalarPolynomialsBase::memory_consumption() + poly_memory + + poly_grads.memory_consumption() + + poly_hessians.memory_consumption() + + poly_third_derivatives.memory_consumption() + + poly_fourth_derivatives.memory_consumption(); + } + + template class BarycentricPolynomials<1>; + template class BarycentricPolynomials<2>; + template class BarycentricPolynomials<3>; + +} // namespace Simplex + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/simplex/fe_lib.cc b/source/simplex/fe_lib.cc index 83f965f66c..d06188f983 100644 --- a/source/simplex/fe_lib.cc +++ b/source/simplex/fe_lib.cc @@ -20,6 +20,7 @@ #include #include +#include #include DEAL_II_NAMESPACE_OPEN diff --git a/tests/simplex/barycentric_01.cc b/tests/simplex/barycentric_01.cc new file mode 100644 index 0000000000..6bd73fbff9 --- /dev/null +++ b/tests/simplex/barycentric_01.cc @@ -0,0 +1,170 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +// Test Simplex::BarycentricPolynomial and Simplex::BarycentricPolynomials. + +#include +#include + +#include +#include +#include + +#include "../tests.h" + +using namespace dealii; + +int +main() +{ + initlog(); + + Simplex::BarycentricPolynomial<2> bp2({1, 0, 0}, 1.0); + deallog << bp2 << std::endl; + + // test some basic algebra with barycentric polynomials + { + deallog << "1D:" << std::endl; + const auto bp1_0 = Simplex::BarycentricPolynomial<1>::monomial(0); + const auto bp1_1 = Simplex::BarycentricPolynomial<1>::monomial(1); + + deallog << "bp1_0 = " << bp1_0 << std::endl; + deallog << "bp1_1 = " << bp1_1 << std::endl; + deallog << "bp1_0 * 2 * bp1_1 / 2 = " << bp1_0 * 2 * bp1_1 / 2 << std::endl + << std::endl; + } + + { + deallog << std::endl << "2D:" << std::endl; + const auto bp2_0 = Simplex::BarycentricPolynomial<2>::monomial(0) * 2; + deallog << "bp2_0 = " << bp2_0 << std::endl; + + const auto bp2_1 = 3.0 * Simplex::BarycentricPolynomial<2>::monomial(1); + deallog << "bp2_1 = " << bp2_1 << std::endl; + + const auto bp2_2 = Simplex::BarycentricPolynomial<2>::monomial(2); + deallog << "bp2_2 = " << bp2_2 << std::endl; + + const auto prod1 = bp2_0 + bp2_1; + deallog << "bp2_0 + bp2_1 = " << prod1 << std::endl; + + const auto prod2 = prod1 * bp2_0; + deallog << "(bp2_0 + bp2_1) * bp2_0 = " << prod2 << std::endl; + deallog << "bp2_0 * bp2_0 + bp2_1 * bp2_0 = " + << bp2_0 * bp2_0 + bp2_1 * bp2_0 << std::endl; + deallog << "bp2_1 * bp2_0 + bp2_0 * bp2_0 = " + << bp2_1 * bp2_0 + bp2_0 * bp2_0 << std::endl; + + // test derivatives + deallog << "d/dx bp2_0 = " << bp2_0.derivative(0) << std::endl; + deallog << "d/dy bp2_0 = " << bp2_0.derivative(1) << std::endl; + + deallog << "d/dx bp2_2 = " << bp2_2.derivative(0) << std::endl; + deallog << "d/dy bp2_2 = " << bp2_2.derivative(1) << std::endl; + } + + // test various finite element spaces + { + deallog << std::endl << "Test with TRI6" << std::endl; + + const auto t1 = Simplex::BarycentricPolynomial<2>::monomial(0); + const auto t2 = Simplex::BarycentricPolynomial<2>::monomial(1); + const auto t3 = Simplex::BarycentricPolynomial<2>::monomial(2); + + std::vector> p2; + p2.push_back(t1 * (2 * t1 - 1)); + p2.push_back(t2 * (2 * t2 - 1)); + p2.push_back(t3 * (2 * t3 - 1)); + p2.push_back(4 * t2 * t1); + p2.push_back(4 * t2 * t3); + p2.push_back(4 * t3 * t1); + + Simplex::FE_P<2> fe(2); + for (unsigned int i = 0; i < 6; ++i) + { + deallog << "p = " << p2[i] << std::endl; + deallog << "p_x = " << p2[i].derivative(0) << std::endl; + deallog << "p_y = " << p2[i].derivative(1) << std::endl; + for (unsigned int j = 0; j < 6; ++j) + { + Assert(std::abs(p2[i].value(fe.get_unit_support_points()[j]) - + double(i == j)) < 1e-12, + ExcInternalError()); + } + deallog << std::endl; + } + } + + { + deallog << std::endl << "Test with TET4" << std::endl; + const auto tet4 = Simplex::BarycentricPolynomials<3>::get_fe_p_basis(1); + + Simplex::FE_P<3> fe(1); + const auto & points = fe.get_unit_support_points(); + for (unsigned int i = 0; i < 4; ++i) + { + Assert(points.size() == 4, ExcInternalError()); + for (unsigned int j = 0; j < 4; ++j) + { + Assert(std::abs(tet4.compute_value(i, points[j]) - double(i == j)) < + 1e-12, + ExcInternalError()); + + // first derivatives should be constant + Assert((tet4.compute_grad(i, points[0]) - + tet4.compute_grad(i, points[j])) + .norm() == 0.0, + ExcInternalError()); + Assert(tet4.compute_2nd_derivative(i, points[j]).norm() == 0.0, + ExcInternalError()); + Assert(tet4.compute_3rd_derivative(i, points[j]).norm() == 0.0, + ExcInternalError()); + Assert(tet4.compute_4th_derivative(i, points[j]).norm() == 0.0, + ExcInternalError()); + } + } + deallog << "Test with TET4 - Success" << std::endl; + } + + { + deallog << "Test with TET10" << std::endl; + const auto tet10 = Simplex::BarycentricPolynomials<3>::get_fe_p_basis(2); + + Simplex::FE_P<3> fe(2); + const auto & points = fe.get_unit_support_points(); + for (unsigned int i = 0; i < 10; ++i) + { + Assert(points.size() == 10, ExcInternalError()); + for (unsigned int j = 0; j < 10; ++j) + { + Assert(std::abs(tet10.compute_value(i, points[j]) - + double(i == j)) < 1e-12, + ExcInternalError()); + + // second derivatives should be constant + Assert((tet10.compute_2nd_derivative(i, points[0]) - + tet10.compute_2nd_derivative(i, points[j])) + .norm() == 0.0, + ExcInternalError()); + + Assert(tet10.compute_3rd_derivative(i, points[j]).norm() == 0.0, + ExcInternalError()); + Assert(tet10.compute_4th_derivative(i, points[j]).norm() == 0.0, + ExcInternalError()); + } + } + deallog << "Test with TET10 - Success" << std::endl; + } +} diff --git a/tests/simplex/barycentric_01.output b/tests/simplex/barycentric_01.output new file mode 100644 index 0000000000..f083f6c228 --- /dev/null +++ b/tests/simplex/barycentric_01.output @@ -0,0 +1,51 @@ + +DEAL::1.00000 * t0^1 +DEAL::1D: +DEAL::bp1_0 = 1.00000 * t0^1 +DEAL::bp1_1 = 1.00000 * t1^1 +DEAL::bp1_0 * 2 * bp1_1 / 2 = 1.00000 * t0^1 * t1^1 +DEAL:: +DEAL:: +DEAL::2D: +DEAL::bp2_0 = 2.00000 * t0^1 +DEAL::bp2_1 = 3.00000 * t1^1 +DEAL::bp2_2 = 1.00000 * t2^1 +DEAL::bp2_0 + bp2_1 = 3.00000 * t1^1 + 2.00000 * t0^1 +DEAL::(bp2_0 + bp2_1) * bp2_0 = 6.00000 * t0^1 * t1^1 + 4.00000 * t0^2 +DEAL::bp2_0 * bp2_0 + bp2_1 * bp2_0 = 6.00000 * t0^1 * t1^1 + 4.00000 * t0^2 +DEAL::bp2_1 * bp2_0 + bp2_0 * bp2_0 = 6.00000 * t0^1 * t1^1 + 4.00000 * t0^2 +DEAL::d/dx bp2_0 = -2.00000 +DEAL::d/dy bp2_0 = -2.00000 +DEAL::d/dx bp2_2 = 0.00000 +DEAL::d/dy bp2_2 = 1.00000 +DEAL:: +DEAL::Test with TRI6 +DEAL::p = -1.00000 * t0^1 + 2.00000 * t0^2 +DEAL::p_x = 1.00000 + -4.00000 * t0^1 +DEAL::p_y = 1.00000 + -4.00000 * t0^1 +DEAL:: +DEAL::p = -1.00000 * t1^1 + 2.00000 * t1^2 +DEAL::p_x = -1.00000 + 4.00000 * t1^1 +DEAL::p_y = 0.00000 +DEAL:: +DEAL::p = -1.00000 * t2^1 + 2.00000 * t2^2 +DEAL::p_x = 0.00000 +DEAL::p_y = -1.00000 + 4.00000 * t2^1 +DEAL:: +DEAL::p = 4.00000 * t0^1 * t1^1 +DEAL::p_x = -4.00000 * t1^1 + 4.00000 * t0^1 +DEAL::p_y = -4.00000 * t1^1 +DEAL:: +DEAL::p = 4.00000 * t1^1 * t2^1 +DEAL::p_x = 4.00000 * t2^1 +DEAL::p_y = 4.00000 * t1^1 +DEAL:: +DEAL::p = 4.00000 * t0^1 * t2^1 +DEAL::p_x = -4.00000 * t2^1 +DEAL::p_y = -4.00000 * t2^1 + 4.00000 * t0^1 +DEAL:: +DEAL:: +DEAL::Test with TET4 +DEAL::Test with TET4 - Success +DEAL::Test with TET10 +DEAL::Test with TET10 - Success -- 2.39.5