From f268d17864db0e9bb31d4f6bb8e907e5935a1be4 Mon Sep 17 00:00:00 2001
From: Denis Davydov <davydden@gmail.com>
Date: Tue, 31 Jan 2017 09:17:52 +0100
Subject: [PATCH] CSpline: add thread mutex

---
 include/deal.II/base/function_cspline.h |  6 ++++++
 source/base/function_cspline.cc         | 15 +++++++++++++++
 2 files changed, 21 insertions(+)

diff --git a/include/deal.II/base/function_cspline.h b/include/deal.II/base/function_cspline.h
index 1653ba0fd2..0338d8f268 100644
--- 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 abb8d87afc..ee8bb54fee 100644
--- 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()));
-- 
2.39.5