From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Thu, 28 Jan 2010 16:28:25 +0000 (+0000)
Subject: Fix a race condition in PolynomialsBDM.
X-Git-Tag: v8.0.0~6575
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=67c8d7824e99201e4a1e63a31dd449fc7fceaf35;p=dealii.git

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 89ed5a7124..71df841afb 100644
--- 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 782348ee7d..c4ae72a235 100644
--- 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 a2bfc6fb4e..f3842f5a72 100644
--- 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)