CSpline: add thread mutex

diff --git a/include/deal.II/base/function_cspline.h b/include/deal.II/base/function_cspline.h
index 1653ba0fd29fcbea770bec038122d7ae2b8c7f60..0338d8f268b39d3104405362ccdc66609b873a0c 100644 (file)
--- a/include/deal.II/base/function_cspline.h
+++ b/include/deal.II/base/function_cspline.h
@@ -21,6 +21,7 @@
 #ifdef DEAL_II_WITH_GSL
 #include <deal.II/base/function.h>
 #include <deal.II/base/point.h>
+#include <deal.II/base/thread_management.h>
 #include <gsl/gsl_spline.h>
 
 DEAL_II_NAMESPACE_OPEN
@@ -115,6 +116,11 @@ namespace Functions
      * GSL cubic spline interpolator
      */
     gsl_spline *cspline;
+
+    /**
+     * A mutex for accelerator object.
+     */
+    mutable Threads::Mutex acc_mutex;
   };
 }
 

diff --git a/source/base/function_cspline.cc b/source/base/function_cspline.cc
index abb8d87afcebeb5c67db37e2fc4669b80c251e2a..ee8bb54fee5fdca8fbb144d22e9e416017b063a9 100644 (file)
--- a/source/base/function_cspline.cc
+++ b/source/base/function_cspline.cc
@@ -67,6 +67,11 @@ namespace Functions
   CSpline<dim>::value (const Point<dim>   &p,
                        const unsigned int) const
   {
+    // GSL functions may modify gsl_interp_accel *acc object (last argument).
+    // This can only work in multithreaded applications if we lock the data
+    // structures via a mutex.
+    Threads::Mutex::ScopedLock lock (acc_mutex);
+
     const double &x = p[0];
     Assert (x >= interpolation_points.front() && x <= interpolation_points.back(),
             ExcCSplineRange(x,interpolation_points.front(),interpolation_points.back()));
@@ -81,6 +86,11 @@ namespace Functions
   CSpline<dim>::gradient (const Point<dim>   &p,
                           const unsigned int) const
   {
+    // GSL functions may modify gsl_interp_accel *acc object (last argument).
+    // This can only work in multithreaded applications if we lock the data
+    // structures via a mutex.
+    Threads::Mutex::ScopedLock lock (acc_mutex);
+
     const double &x = p[0];
     Assert (x >= interpolation_points.front() && x <= interpolation_points.back(),
             ExcCSplineRange(x,interpolation_points.front(),interpolation_points.back()));
@@ -98,6 +108,11 @@ namespace Functions
   CSpline<dim>::laplacian (const Point<dim>   &p,
                            const unsigned int) const
   {
+    // GSL functions may modify gsl_interp_accel *acc object (last argument).
+    // This can only work in multithreaded applications if we lock the data
+    // structures via a mutex.
+    Threads::Mutex::ScopedLock lock (acc_mutex);
+
     const double &x = p[0];
     Assert (x >= interpolation_points.front() && x <= interpolation_points.back(),
             ExcCSplineRange(x,interpolation_points.front(),interpolation_points.back()));