From de55f6fe5fd7bf8c5a9bee8f8b43c41677a000f3 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 8 Jan 2017 21:28:25 -0700 Subject: [PATCH] Fix the computation of the number of polynomials the ABF element has. --- source/base/polynomials_abf.cc | 36 +++++++++++++++++++++++++++------- 1 file changed, 29 insertions(+), 7 deletions(-) diff --git a/source/base/polynomials_abf.cc b/source/base/polynomials_abf.cc index f2718a00a2..ae13827fa4 100644 --- a/source/base/polynomials_abf.cc +++ b/source/base/polynomials_abf.cc @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2004 - 2015 by the deal.II authors +// Copyright (C) 2004 - 2017 by the deal.II authors // // This file is part of the deal.II library. // @@ -31,13 +31,21 @@ PolynomialsABF::PolynomialsABF (const unsigned int k) { std::vector > > pols(dim); pols[0] = Polynomials::LagrangeEquidistant::generate_complete_basis(k+2); + if (k == 0) for (unsigned int d=1; d(pols); + + // check that the dimensions match. we only store one of the 'dim' + // anisotropic polynomials that make up the vector-valued space, so + // multiply by 'dim' + Assert (dim * polynomial_space->n() == compute_n_pols(k), + ExcInternalError()); } @@ -135,14 +143,28 @@ PolynomialsABF::compute (const Point &unit_point, template unsigned int -PolynomialsABF::compute_n_pols(unsigned int k) +PolynomialsABF::compute_n_pols(const unsigned int k) { - if (dim == 1) return k+1; - if (dim == 2) return 2*(k+1)*(k+3); - //TODO:Check what are the correct numbers ... - if (dim == 3) return 3*(k+1)*(k+1)*(k+2); + switch (dim) + { + case 1: + // in 1d, we simply have Q_{k+2}, which has dimension k+3 + return k+3; + + case 2: + // the polynomial space is Q_{k+2,k} \times Q_{k,k+2}, which has + // 2(k+3)(k+1) DoFs + return 2*(k+3)*(k+1); + + case 3: + // the polynomial space is Q_{k+2,k,k} \times Q_{k,k+2,k} \times Q_{k,k,k+2}, + // which has 3(k+3)(k+1)(k+1) DoFs + return 3*(k+3)*(k+1)*(k+1); + + default: + Assert(false, ExcNotImplemented()); + } - Assert(false, ExcNotImplemented()); return 0; } -- 2.39.5