From 24cdb5bfd56b22d141f80c7a8c7a87352f02b145 Mon Sep 17 00:00:00 2001
From: alberto sartori <alberto.sartori86@gmail.com>
Date: Tue, 9 Feb 2016 12:58:47 +0100
Subject: [PATCH] added rand and rand_seed functions to FunctionParser

---
 doc/news/changes.h                   |  9 ++++-
 source/base/function_parser.cc       | 42 ++++++++++++++++++++-
 tests/base/function_parser_09.cc     | 56 ++++++++++++++++++++++++++++
 tests/base/function_parser_09.output |  3 ++
 4 files changed, 107 insertions(+), 3 deletions(-)
 create mode 100644 tests/base/function_parser_09.cc
 create mode 100644 tests/base/function_parser_09.output

diff --git a/doc/news/changes.h b/doc/news/changes.h
index eb4562206c..4f2b68b021 100644
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -512,7 +512,14 @@ inconvenience this causes.
 
 
 
-<ol>
+<ol> 
+  <li> New: FunctionParser now supports <code>rand()</code> and
+  <code>rand_seed(number)</code>, which return a random value in the
+  range [0,1].
+  <br>
+  (Wolfgang Bangerth, Luca Heltai, Alberto Sartori, 2016/02/09)
+  </li>
+
   <li> Fixed: FullMatrix::TmTmult for matrix multiplication used to compute
   wrong results for larger matrix sizes where external BLAS is called. This
   has been fixed.
diff --git a/source/base/function_parser.cc b/source/base/function_parser.cc
index badd034b01..d5ff52133e 100644
--- a/source/base/function_parser.cc
+++ b/source/base/function_parser.cc
@@ -16,9 +16,11 @@
 
 #include <deal.II/base/function_parser.h>
 #include <deal.II/base/utilities.h>
+#include <deal.II/base/thread_management.h>
 #include <deal.II/lac/vector.h>
 #include <cmath>
-
+#include <boost/random.hpp>
+#include <map>
 
 #ifdef DEAL_II_WITH_MUPARSER
 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
@@ -26,6 +28,8 @@ DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
 #else
 
+
+
 namespace fparser
 {
   class FunctionParser
@@ -176,6 +180,36 @@ namespace internal
   {
     return erfc(value);
   }
+
+  // returns a random value in the range [0,1] initializing the generator
+  // with the given seed
+  double mu_rand_seed(double seed)
+  {
+    static Threads::Mutex rand_mutex;
+    Threads::Mutex::ScopedLock lock(rand_mutex);
+
+    static boost::random::uniform_real_distribution<> uniform_distribution(0,1);
+
+    // for each seed an unique random number generator is created,
+    // which is initialized with the seed itself
+    static std::map<double, boost::random::mt19937> rng_map;
+
+    if (rng_map.find(seed) == rng_map.end())
+      rng_map[seed] = boost::random::mt19937(static_cast<unsigned int>(seed));
+
+    return uniform_distribution(rng_map[seed]);
+  }
+
+  // returns a random value in the range [0,1]
+  double mu_rand()
+  {
+    static Threads::Mutex rand_mutex;
+    Threads::Mutex::ScopedLock lock(rand_mutex);
+    static boost::random::uniform_real_distribution<> uniform_distribution(0,1);
+    static boost::random::mt19937 rng(static_cast<unsigned long>(std::time(0)));
+    return uniform_distribution(rng);
+  }
+
 }
 
 
@@ -215,6 +249,8 @@ void FunctionParser<dim>:: init_muparser() const
       fp.get()[component].DefineFun("log", internal::mu_log, true);
       fp.get()[component].DefineFun("pow", internal::mu_pow, true);
       fp.get()[component].DefineFun("erfc", internal::mu_erfc, true);
+      fp.get()[component].DefineFun("rand_seed", internal::mu_rand_seed, true);
+      fp.get()[component].DefineFun("rand", internal::mu_rand, true);
 
       try
         {
@@ -265,7 +301,9 @@ void FunctionParser<dim>:: init_muparser() const
             "floor",
             "sec",
             "pow",
-            "erfc"
+            "erfc",
+            "rand",
+            "rand_seed"
           };
           for (unsigned int f=0; f<sizeof(function_names)/sizeof(function_names[0]); ++f)
             {
diff --git a/tests/base/function_parser_09.cc b/tests/base/function_parser_09.cc
new file mode 100644
index 0000000000..03d294e519
--- /dev/null
+++ b/tests/base/function_parser_09.cc
@@ -0,0 +1,56 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2005 - 2016 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// test rand function
+
+#include "../tests.h"
+#include <fstream>
+#include <iomanip>
+#include <map>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/point.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/base/function_parser.h>
+
+
+double eval(const std::string &exp)
+{
+  std::string variables = "x,y";
+  std::map<std::string,double> constants;
+
+  FunctionParser<2> fp(1);
+  fp.initialize(variables,
+                exp,
+                constants);
+
+  Point<2> p;
+
+  return fp.value(p);
+}
+
+
+int main()
+{
+  initlog();
+
+  double random = eval("rand()"); //random seed
+
+  if (0.0 <= random && random <= 1.0)
+    deallog << "OK" << std::endl;
+
+  deallog << eval("rand_seed(10)") << std::endl;
+
+}
diff --git a/tests/base/function_parser_09.output b/tests/base/function_parser_09.output
new file mode 100644
index 0000000000..59f845a429
--- /dev/null
+++ b/tests/base/function_parser_09.output
@@ -0,0 +1,3 @@
+
+DEAL::OK
+DEAL::0.771321
-- 
2.39.5