Fix a race condition in PolynomialsBDM.

git-svn-id: https://svn.dealii.org/trunk@20467 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/base/include/base/polynomials_bdm.h b/deal.II/base/include/base/polynomials_bdm.h
index 89ed5a71243887cbe39c0bfce6fff9b0887308e9..71df841afbff48834f837049b3f61ed803c9cb80 100644 (file)
--- a/deal.II/base/include/base/polynomials_bdm.h
+++ b/deal.II/base/include/base/polynomials_bdm.h
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004, 2005, 2006, 2009 by the deal.II authors
+//    Copyright (C) 2004, 2005, 2006, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -21,6 +21,7 @@
 #include <base/polynomial.h>
 #include <base/polynomial_space.h>
 #include <base/table.h>
+#include <base/thread_management.h>
 
 #include <vector>
 
@@ -46,7 +47,7 @@ DEAL_II_NAMESPACE_OPEN
  * <dt> In 2D:
  * <dd> The 2D-curl of the functions <i>x<sup>k+1</sup>y</i>
  * and <i>xy<sup>k+1</sup></i>.
- * <dt>In 3D: 
+ * <dt>In 3D:
  * <dd> For any <i>i=0,...,k</i> the curls of
  * <i>(0,0,xy<sup>i+1</sup>z<sup>k-i</sup>)</i>,
  * <i>(x<sup>k-i</sup>yz<sup>i+1</sup>,0,0)</i> and
@@ -102,20 +103,20 @@ class PolynomialsBDM
     void compute (const Point<dim>            &unit_point,
                   std::vector<Tensor<1,dim> > &values,
                   std::vector<Tensor<2,dim> > &grads,
-                  std::vector<Tensor<3,dim> > &grad_grads) const;    
-    
+                  std::vector<Tensor<3,dim> > &grad_grads) const;
+
                                     /**
                                      * Returns the number of BDM polynomials.
                                      */
     unsigned int n () const;
-    
+
                                     /**
                                      * Returns the degree of the BDM
                                      * space, which is one less than
                                      * the highest polynomial degree.
                                      */
     unsigned int degree () const;
-    
+
                                     /**
                                      * Return the number of
                                      * polynomials in the space
@@ -126,7 +127,7 @@ class PolynomialsBDM
                                      * classes.
                                      */
     static unsigned int compute_n_pols(unsigned int degree);
-    
+
   private:
                                     /**
                                      * An object representing the
@@ -144,23 +145,29 @@ class PolynomialsBDM
                                      * degree zero to <i>k</i>.
                                      */
     std::vector<Polynomials::Polynomial<double> > monomials;
-    
+
                                     /**
                                      * Number of BDM
                                      * polynomials.
                                      */
     unsigned int n_pols;
-    
+
+                                    /**
+                                     * A mutex that guards the
+                                     * following scratch arrays.
+                                     */
+    mutable Threads::Mutex mutex;
+
                                     /**
                                      * Auxiliary memory.
                                      */
     mutable std::vector<double> p_values;
-    
+
                                     /**
                                      * Auxiliary memory.
                                      */
     mutable std::vector<Tensor<1,dim> > p_grads;
-    
+
                                     /**
                                      * Auxiliary memory.
                                      */

diff --git a/deal.II/base/source/polynomials_bdm.cc b/deal.II/base/source/polynomials_bdm.cc
index 782348ee7d904892401b653a951604c8cb2aaba4..c4ae72a235fefbfa78fd2964b70edda90bb8f9f3 100644 (file)
--- a/deal.II/base/source/polynomials_bdm.cc
+++ b/deal.II/base/source/polynomials_bdm.cc
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
+//    Copyright (C) 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -60,40 +60,41 @@ PolynomialsBDM<dim>::compute (const Point<dim>            &unit_point,
         ExcDimensionMismatch(grad_grads.size(), n_pols));
 
   const unsigned int n_sub = polynomial_space.n();
-  p_values.resize((values.size() == 0) ? 0 : n_sub);
-  p_grads.resize((grads.size() == 0) ? 0 : n_sub);
-  p_grad_grads.resize((grad_grads.size() == 0) ? 0 : n_sub);
-
-                                  // Compute values of complete space
-                                  // and insert into tensors.  Result
-                                  // will have first all polynomials
-                                  // in the x-component, then y and
-                                  // z.
-  polynomial_space.compute (unit_point, p_values, p_grads, p_grad_grads);
-  
-  std::fill(values.begin(), values.end(), Tensor<1,dim>());
-  for (unsigned int i=0;i<p_values.size();++i)
-    {
+
+                                  // guard access to the scratch
+                                  // arrays in the following block
+                                  // using a mutex to make sure they
+                                  // are not used by multiple threads
+                                  // at once
+  {
+    Threads::Mutex::ScopedLock lock(mutex);
+
+    p_values.resize((values.size() == 0) ? 0 : n_sub);
+    p_grads.resize((grads.size() == 0) ? 0 : n_sub);
+    p_grad_grads.resize((grad_grads.size() == 0) ? 0 : n_sub);
+
+                                    // Compute values of complete space
+                                    // and insert into tensors.  Result
+                                    // will have first all polynomials
+                                    // in the x-component, then y and
+                                    // z.
+    polynomial_space.compute (unit_point, p_values, p_grads, p_grad_grads);
+
+    std::fill(values.begin(), values.end(), Tensor<1,dim>());
+    for (unsigned int i=0;i<p_values.size();++i)
       for (unsigned int j=0;j<dim;++j)
-       {
-         values[i+j*n_sub][j] = p_values[i];
-       }
-      
-    }
-  
-  std::fill(grads.begin(), grads.end(), Tensor<2,dim>());
-  for (unsigned int i=0;i<p_grads.size();++i)
-    {
+       values[i+j*n_sub][j] = p_values[i];
+
+    std::fill(grads.begin(), grads.end(), Tensor<2,dim>());
+    for (unsigned int i=0;i<p_grads.size();++i)
       for (unsigned int j=0;j<dim;++j)
        grads[i+j*n_sub][j] = p_grads[i];
-    }
-  
-  std::fill(grad_grads.begin(), grad_grads.end(), Tensor<3,dim>());
-  for (unsigned int i=0;i<p_grad_grads.size();++i)
-    {
+
+    std::fill(grad_grads.begin(), grad_grads.end(), Tensor<3,dim>());
+    for (unsigned int i=0;i<p_grad_grads.size();++i)
       for (unsigned int j=0;j<dim;++j)
        grad_grads[i+j*n_sub][j] = p_grad_grads[i];
-    }
+  }
 
                                   // This is the first polynomial not
                                   // covered by the P_k subspace
@@ -104,7 +105,7 @@ PolynomialsBDM<dim>::compute (const Point<dim>            &unit_point,
                                   // derivatives
   std::vector<std::vector<double> > monovali(dim, std::vector<double>(4));
   std::vector<std::vector<double> > monovalk(dim, std::vector<double>(4));
-  
+
   if (dim == 2)
     {
       for (unsigned int d=0;d<dim;++d)
@@ -252,7 +253,7 @@ PolynomialsBDM<dim>::compute_node_matrix (Table<2,double>& A) const
       double orientation = 1.;
       if ((face==0) || (face==3))
        orientation = -1.;
-      
+
       for (unsigned int k=0;k<qface.size();++k)
        {
          const double w = qface.weight(k) * orientation;
@@ -269,15 +270,15 @@ PolynomialsBDM<dim>::compute_node_matrix (Table<2,double>& A) const
                p(0) = 1.;
              case 3:
                p(1) = x;
-               break;        
+               break;
            }
 //     std::cerr << p
 //               << '\t' << moment_weight[0].value(x)
 //               << '\t' << moment_weight[1].value(x)
 //       ;
-         
+
          compute (p, values, grads, grad_grads);
-         
+
          for (unsigned int i=0;i<n();++i)
            {
 //         std::cerr << '\t' << std::setw(6) << values[i][1-face%2];
@@ -290,7 +291,7 @@ PolynomialsBDM<dim>::compute_node_matrix (Table<2,double>& A) const
 //     std::cerr << std::endl;
        }
     }
-  
+
                                   // Volume integrals are missing
                                   //
                                   // This degree is one larger

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index a2bfc6fb4e0782add60c82210786982ed340e24b..f3842f5a72edaf4c65783707cb57380f65bb43f8 100644 (file)
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -246,6 +246,12 @@ inconvenience this causes.
 <h3>base</h3>
 
 <ol>
+  <li><p>Fixed: The PolynomialsBDM had a race condition in multithreaded
+  programs. This now fixed.
+  <br>
+  (WB 2010/01/27)
+  </p></li>
+
   <li><p>New: The new class IndexSet can represent sets and ranges of indices.
   <br>
   (WB 2009/10/09)