From 0cc30f6a1d072099aef45e32593fbc2af30bd1ca Mon Sep 17 00:00:00 2001
From: Norihiro Watanabe <norihiro.watanabe@aist.go.jp>
Date: Thu, 31 Jan 2019 00:18:06 +0900
Subject: [PATCH] [PL] add FunctionParameter

---
 ProcessLib/Parameter/FunctionParameter.cpp |  37 +++++++
 ProcessLib/Parameter/FunctionParameter.h   | 116 +++++++++++++++++++++
 ProcessLib/Parameter/Parameter.cpp         |   7 ++
 3 files changed, 160 insertions(+)
 create mode 100644 ProcessLib/Parameter/FunctionParameter.cpp
 create mode 100644 ProcessLib/Parameter/FunctionParameter.h

diff --git a/ProcessLib/Parameter/FunctionParameter.cpp b/ProcessLib/Parameter/FunctionParameter.cpp
new file mode 100644
index 00000000000..518867f0bba
--- /dev/null
+++ b/ProcessLib/Parameter/FunctionParameter.cpp
@@ -0,0 +1,37 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "FunctionParameter.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "MeshLib/Mesh.h"
+
+namespace ProcessLib
+{
+std::unique_ptr<ParameterBase> createFunctionParameter(
+    std::string const& name, BaseLib::ConfigTree const& config,
+    MeshLib::Mesh const& mesh)
+{
+    //! \ogs_file_param{prj__parameters__parameter__type}
+    config.checkConfigParameter("type", "Function");
+
+    std::vector<std::string> vec_expressions;
+
+    //! \ogs_file_param{prj__parameters__parameter__Function__expression}
+    for (auto p : config.getConfigSubtreeList("expression"))
+    {
+        std::string const expression_str = p.getValue<std::string>();
+        vec_expressions.emplace_back(expression_str);
+    }
+
+    return std::make_unique<FunctionParameter<double>>(
+        name, mesh, vec_expressions);
+}
+
+}  // ProcessLib
diff --git a/ProcessLib/Parameter/FunctionParameter.h b/ProcessLib/Parameter/FunctionParameter.h
new file mode 100644
index 00000000000..39340a5d51c
--- /dev/null
+++ b/ProcessLib/Parameter/FunctionParameter.h
@@ -0,0 +1,116 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <utility>
+#include <vector>
+
+#include <exprtk.hpp>
+
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Node.h"
+
+#include "Parameter.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+namespace ProcessLib
+{
+
+/// A parameter class evaluating functons defined by
+/// user-provided mathematical expressions.
+///
+/// Currently, x, y, and z are supported as variables
+/// of the functions.
+template <typename T>
+struct FunctionParameter final : public Parameter<T>
+{
+    typedef exprtk::symbol_table<T> symbol_table_t;
+    typedef exprtk::expression<T> expression_t;
+    typedef exprtk::parser<T> parser_t;
+    typedef exprtk::parser_error::type error_t;
+
+    /**
+     * Constructing from a vector of expressions
+     *
+     * @param name_       the parameter's name
+     * @param mesh_       a mesh object
+     * @param vec_expression_str_  a vector of mathematical expressions
+     * The vector size specifies the number of components of the parameter.
+     */
+    FunctionParameter(std::string const& name_,
+                      MeshLib::Mesh const& mesh_,
+                      std::vector<std::string> const& vec_expression_str_)
+        : Parameter<T>(name_), _mesh(mesh_), _vec_expression_str(vec_expression_str_)
+    {
+        _symbol_table.add_constants();
+        _symbol_table.create_variable("x");
+        _symbol_table.create_variable("y");
+        _symbol_table.create_variable("z");
+
+        _vec_expression.resize(_vec_expression_str.size());
+        for (unsigned i=0; i<_vec_expression_str.size(); i++)
+        {
+            _vec_expression[i].register_symbol_table(_symbol_table);
+            parser_t parser;
+            if (!parser.compile(_vec_expression_str[i], _vec_expression[i]))
+            {
+                OGS_FATAL("Error: %s\tExpression: %s\n", parser.error().c_str(),
+                    _vec_expression_str[i].c_str());
+            }
+        }
+    }
+
+    bool isTimeDependent() const override { return false; }
+
+    int getNumberOfComponents() const override
+    {
+        return _vec_expression.size();
+    }
+
+    std::vector<T> operator()(double const /*t*/,
+                                     SpatialPosition const& pos) const override
+    {
+        std::vector<T> cache(getNumberOfComponents());
+        auto& x = _symbol_table.get_variable("x")->ref();
+        auto& y = _symbol_table.get_variable("y")->ref();
+        auto& z = _symbol_table.get_variable("z")->ref();
+        if (pos.getCoordinates())
+        {
+            auto const coords = pos.getCoordinates().get();
+            x = coords[0];
+            y = coords[1];
+            z = coords[2];
+        }
+        else if (pos.getNodeID())
+        {
+            auto const& node = *_mesh.getNode(pos.getNodeID().get());
+            x = node[0];
+            y = node[1];
+            z = node[2];
+        }
+
+        for (unsigned i=0; i<_vec_expression.size(); i++)
+            cache[i] = _vec_expression[i].value();
+
+        return cache;
+    }
+
+private:
+    MeshLib::Mesh const& _mesh;
+    std::vector<std::string> _vec_expression_str;
+    symbol_table_t _symbol_table;
+    std::vector<expression_t> _vec_expression;
+};
+
+std::unique_ptr<ParameterBase> createFunctionParameter(
+    std::string const& name, BaseLib::ConfigTree const& config,
+    MeshLib::Mesh const& mesh);
+
+}  // ProcessLib
diff --git a/ProcessLib/Parameter/Parameter.cpp b/ProcessLib/Parameter/Parameter.cpp
index 03013ead1c5..0061e4ab436 100644
--- a/ProcessLib/Parameter/Parameter.cpp
+++ b/ProcessLib/Parameter/Parameter.cpp
@@ -13,6 +13,7 @@
 
 #include "ConstantParameter.h"
 #include "CurveScaledParameter.h"
+#include "FunctionParameter.h"
 #include "GroupBasedParameter.h"
 #include "MeshElementParameter.h"
 #include "MeshNodeParameter.h"
@@ -45,6 +46,12 @@ std::unique_ptr<ParameterBase> createParameter(
         auto param = createCurveScaledParameter(name, config, curves);
         return param;
     }
+    if (type == "Function")
+    {
+        INFO("FunctionParameter: %s", name.c_str());
+        auto param = createFunctionParameter(name, config, *meshes.front());
+        return param;
+    }
     if (type == "Group")
     {
         INFO("GroupBasedParameter: %s", name.c_str());
-- 
GitLab