From ec43010240f66c7c7c6d96288801f8ca3a1c9f7c Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Mon, 9 Jan 2017 08:52:15 -0700
Subject: [PATCH] Do not use a plain pointer in PolynomialsABF.

These objects are being copied by the FE_Poly* classes, so plain pointers without
dedicated copy constructors and operators are likely going to lead to memory
corruption. It's not clear to me how this ever worked, but it's easy to fix.
---
 include/deal.II/base/polynomials_abf.h | 13 +++----
 source/base/polynomials_abf.cc         | 48 ++++++++++++++------------
 2 files changed, 31 insertions(+), 30 deletions(-)

diff --git a/include/deal.II/base/polynomials_abf.h b/include/deal.II/base/polynomials_abf.h
index 3ce995f56e..f183c406ee 100644
--- a/include/deal.II/base/polynomials_abf.h
+++ b/include/deal.II/base/polynomials_abf.h
@@ -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.
 //
@@ -27,6 +27,7 @@
 #include <deal.II/base/table.h>
 #include <deal.II/base/thread_management.h>
 
+#include <deal.II/base/std_cxx11/shared_ptr.h>
 #include <vector>
 
 DEAL_II_NAMESPACE_OPEN
@@ -63,11 +64,6 @@ public:
    */
   PolynomialsABF (const unsigned int k);
 
-  /**
-   * Destructor deleting the polynomials.
-   */
-  ~PolynomialsABF ();
-
   /**
    * Compute the value and the first and second derivatives of each Raviart-
    * Thomas polynomial at @p unit_point.
@@ -118,9 +114,10 @@ private:
 
   /**
    * An object representing the polynomial space for a single component. We
-   * can re-use it by rotating the coordinates of the evaluation point.
+   * can re-use it for the other vector components by rotating the
+   * coordinates of the evaluation point.
    */
-  AnisotropicPolynomials<dim> *polynomial_space;
+  const AnisotropicPolynomials<dim> polynomial_space;
 
   /**
    * Number of Raviart-Thomas polynomials.
diff --git a/source/base/polynomials_abf.cc b/source/base/polynomials_abf.cc
index ae13827fa4..c420338a97 100644
--- a/source/base/polynomials_abf.cc
+++ b/source/base/polynomials_abf.cc
@@ -23,38 +23,42 @@
 DEAL_II_NAMESPACE_OPEN
 
 
+
+namespace
+{
+  template <int dim>
+  std::vector<std::vector< Polynomials::Polynomial< double > > >
+  get_abf_polynomials (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);
+
+    return pols;
+  }
+}
+
 template <int dim>
 PolynomialsABF<dim>::PolynomialsABF (const unsigned int k)
   :
   my_degree(k),
+  polynomial_space(get_abf_polynomials<dim>(k)),
   n_pols(compute_n_pols(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),
+  Assert (dim * polynomial_space.n() == compute_n_pols(k),
           ExcInternalError());
 }
 
 
-template <int dim>
-PolynomialsABF<dim>::~PolynomialsABF ()
-{
-  delete polynomial_space;
-}
-
 
 template <int dim>
 void
@@ -76,7 +80,7 @@ PolynomialsABF<dim>::compute (const Point<dim>            &unit_point,
   Assert(fourth_derivatives.size()==n_pols|| fourth_derivatives.size()==0,
          ExcDimensionMismatch(fourth_derivatives.size(), n_pols));
 
-  const unsigned int n_sub = polynomial_space->n();
+  const unsigned int n_sub = polynomial_space.n();
   // guard access to the scratch
   // arrays in the following block
   // using a mutex to make sure they
@@ -107,8 +111,8 @@ PolynomialsABF<dim>::compute (const Point<dim>            &unit_point,
       for (unsigned int c=0; c<dim; ++c)
         p(c) = unit_point((c+d)%dim);
 
-      polynomial_space->compute (p, p_values, p_grads, p_grad_grads,
-                                 p_third_derivatives, p_fourth_derivatives);
+      polynomial_space.compute (p, p_values, p_grads, p_grad_grads,
+                                p_third_derivatives, p_fourth_derivatives);
 
       for (unsigned int i=0; i<p_values.size(); ++i)
         values[i+d*n_sub][d] = p_values[i];
-- 
2.39.5