From de55f6fe5fd7bf8c5a9bee8f8b43c41677a000f3 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
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<dim>::PolynomialsABF (const unsigned int k)
 {
   std::vector<std::vector< Polynomials::Polynomial< double > > > pols(dim);
   pols[0] = Polynomials::LagrangeEquidistant::generate_complete_basis(k+2);
+
   if (k == 0)
     for (unsigned int d=1; d<dim; ++d)
       pols[d] = Polynomials::Legendre::generate_complete_basis(0);
   else
     for (unsigned int d=1; d<dim; ++d)
       pols[d] = Polynomials::LagrangeEquidistant::generate_complete_basis(k);
+
   polynomial_space = new AnisotropicPolynomials<dim>(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<dim>::compute (const Point<dim>            &unit_point,
 
 template <int dim>
 unsigned int
-PolynomialsABF<dim>::compute_n_pols(unsigned int k)
+PolynomialsABF<dim>::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