From 6d9e7f88a66edebf8aad14f4068520ac60a8b9ed Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Mon, 30 Aug 2021 19:13:35 -0600
Subject: [PATCH] Improve documentation of PolynomialsRaviartThomas a bit.

---
 source/base/polynomials_raviart_thomas.cc | 31 +++++++++++++++++++----
 1 file changed, 26 insertions(+), 5 deletions(-)

diff --git a/source/base/polynomials_raviart_thomas.cc b/source/base/polynomials_raviart_thomas.cc
index a728e7378a..ddf589a918 100644
--- a/source/base/polynomials_raviart_thomas.cc
+++ b/source/base/polynomials_raviart_thomas.cc
@@ -46,14 +46,35 @@ template <int dim>
 std::vector<std::vector<Polynomials::Polynomial<double>>>
 PolynomialsRaviartThomas<dim>::create_polynomials(const unsigned int k)
 {
+  // Create a vector of polynomial spaces where the first element
+  // has degree k+1 and the rest has degree k. This corresponds to
+  // the space of single-variable polynomials from which we will create the
+  // space for the first component of the RT element by way of tensor
+  // product.
+  //
+  // The other components of the RT space can be created by rotating
+  // this vector of single-variable polynomials.
+  //
   std::vector<std::vector<Polynomials::Polynomial<double>>> pols(dim);
-  pols[0] = Polynomials::LagrangeEquidistant::generate_complete_basis(k + 1);
   if (k == 0)
-    for (unsigned int d = 1; d < dim; ++d)
-      pols[d] = Polynomials::Legendre::generate_complete_basis(0);
+    {
+      // We need to treat the case k=0 differently because there,
+      // the polynomial in x has degree 1 and so has node points
+      // equal to the end points of the interval [0,1] (i.e., it
+      // is a "Lagrange" polynomial), whereas the y- and z-polynomials
+      // only have the interval midpoint as a node (i.e., they are
+      // a "Legendre" polynomial).
+      pols[0] = Polynomials::LagrangeEquidistant::generate_complete_basis(1);
+      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);
+    {
+      pols[0] =
+        Polynomials::LagrangeEquidistant::generate_complete_basis(k + 1);
+      for (unsigned int d = 1; d < dim; ++d)
+        pols[d] = Polynomials::LagrangeEquidistant::generate_complete_basis(k);
+    }
 
   return pols;
 }
-- 
2.39.5