//
// ---------------------------------------------------------------------
-
#include <deal.II/base/polynomials_barycentric.h>
+#include <deal.II/grid/reference_cell.h>
+
DEAL_II_NAMESPACE_OPEN
namespace internal
{
std::vector<PolyType> polys;
- auto M = [](const unsigned int d) {
- return BarycentricPolynomial<dim, double>::monomial(d);
- };
+ const auto reference_cell = ReferenceCells::get_simplex<dim>();
+
switch (degree)
{
case 0:
- polys.push_back(0 * M(0) + 1);
+ polys.push_back(0 * BarycentricPolynomial<dim, double>::monomial(0) +
+ 1);
break;
case 1:
{
- for (unsigned int d = 0; d < dim + 1; ++d)
- polys.push_back(M(d));
+ for (unsigned int v : reference_cell.vertex_indices())
+ polys.push_back(BarycentricPolynomial<dim, double>::monomial(v));
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)
+ // vertices, then lines:
+ for (unsigned int v : reference_cell.vertex_indices())
+ polys.push_back(
+ BarycentricPolynomial<dim, double>::monomial(v) *
+ (2 * BarycentricPolynomial<dim, double>::monomial(v) - 1));
+ for (unsigned int l : reference_cell.line_indices())
{
- polys.push_back(4 * M(3) * M(0));
- polys.push_back(4 * M(1) * M(3));
- polys.push_back(4 * M(2) * M(3));
+ const auto v0 = reference_cell.line_to_cell_vertices(l, 0);
+ const auto v1 = reference_cell.line_to_cell_vertices(l, 1);
+ polys.push_back(4 *
+ BarycentricPolynomial<dim, double>::monomial(v0) *
+ BarycentricPolynomial<dim, double>::monomial(v1));
}
break;
}