added: python defined scalar|real|vec|tensor|stensor functions

cmake/Modules/FindIFEMDeps.cmake

@@ -254,6 +254,19 @@ if(IFEM_USE_ZOLTAN)
   endif()
 endif()
 
+# Python
+if(IFEM_USE_PYTHON)
+  find_package(Python3 COMPONENTS Development)
+  find_package(pybind11)
+  if(pybind11_FOUND)
+    include(CMakePrintHelpers)
+    list(APPEND IFEM_DEPLIBS pybind11::pybind11 Python3::Python)
+    list(APPEND IFEM_DEPINCLUDES ${pybind11_INCLUDE_DIRS})
+    list(APPEND IFEM_DEFINITIONS -DHAS_PYTHON=1)
+    message(STATUS "Python support enabled")
+  endif()
+endif()
+
 # Portability issues
 include(CheckFunctionExists)
 set(CMAKE_REQUIRED_DEFINITIONS)

cmake/Modules/IFEMOptions.cmake

@@ -13,6 +13,7 @@ OPTION(IFEM_USE_VTFWRITER      "Compile with VTFWriter support?"      ON)
 OPTION(IFEM_USE_UMFPACK        "Compile with UMFPACK support?"        ON)
 OPTION(IFEM_USE_CEREAL         "Compile with cereal support?"         ON)
 OPTION(IFEM_USE_ZOLTAN         "Compile with zoltan support?"         OFF)
+OPTION(IFEM_USE_PYTHON         "Compile with python support?"         OFF)
 OPTION(IFEM_AS_SUBMODULE       "Compile IFEM as a submodule of apps?" OFF)
 OPTION(IFEM_WHOLE_PROG_OPTIM   "Compile IFEM with link-time optimizations?" OFF)
 OPTION(IFEM_TEST_MEMCHECK      "Run tests through valgrind?"          OFF)

cmake/Scripts/IFEMTesting.cmake

@@ -71,6 +71,10 @@ macro(IFEM_add_unittests IFEM_PATH)
     list(REMOVE_ITEM TEST_SOURCES ${IFEM_PATH}/src/Utility/Test/TestFieldFunctionsLR.C)
   endif()
 
+  if(NOT pybind11_FOUND)
+    list(REMOVE_ITEM TEST_SOURCES ${IFEM_PATH}/src/Utility/Test/TestPythonFunctions.C)
+  endif()
+
   IFEM_add_test_app("${TEST_SOURCES}"
                     ${IFEM_PATH}
                     IFEM

src/Utility/Functions.C

@@ -15,6 +15,7 @@
 #include "Chebyshev.h"
 #include "ExprFunctions.h"
 #include "FieldFunctions.h"
+#include "PythonFunctions.h"
 #include "TensorFunction.h"
 #include "Vec3Oper.h"
 #include "IFEM.h"
@@ -681,6 +682,16 @@ const ScalarFunc* utl::parseTimeFunction (const char* type, char* cline, Real C)
     }
     return sf;
   }
+#ifdef HAS_PYTHON
+  else if (strncasecmp(type,"python",6) == 0 && cline != nullptr) {
+    std::string func(cline);
+    auto pos = func.find_first_of('{');
+    std::string module = func.substr(0,pos-1);
+    std::string params = func.substr(pos);
+    IFEM::cout << "\n\t\tModule = " << module << "\n\t\tParams = " << params << std::endl;
+    return new PythonFunc(module.c_str(), params.c_str());
+  }
+#endif
   else if (strncasecmp(type,"Ramp",4) == 0 || strcmp(type,"Tinit") == 0)
   {
     Real xmax = atof(strtok(cline," "));
@@ -798,6 +809,15 @@ RealFunc* utl::parseRealFunc (const std::string& func,
     }
     return rf;
   }
+#ifdef HAS_PYTHON
+  else if (type == "python") {
+    auto pos = func.find_first_of('{');
+    std::string module = func.substr(0,pos-1);
+    std::string params = func.substr(pos);
+    IFEM::cout << "\n\t\tModule = " << module << "\n\t\tParams = " << params << std::endl;
+    return new PythonFunction(module.c_str(), params.c_str());
+  }
+#endif
   else if (type == "linear")
   {
     p = atof(func.c_str());
@@ -867,6 +887,15 @@ VecFunc* utl::parseVecFunc (const std::string& func, const std::string& type,
     return new ChebyshevVecFunc(splitValue(func),false);
   else if (type == "chebyshev2")
     return new ChebyshevVecFunc(splitValue(func),true);
+#ifdef HAS_PYTHON
+  else if (type == "python") {
+    auto pos = func.find_first_of('{');
+    std::string module = func.substr(0,pos-1);
+    std::string params = func.substr(pos);
+    IFEM::cout << "\n\t\tModule = " << module << "\n\t\tParams = " << params << std::endl;
+    return new PythonVecFunc(module.c_str(), params.c_str());
+  }
+#endif
 
   return nullptr;
 }
@@ -879,6 +908,15 @@ TensorFunc* utl::parseTensorFunc (const std::string& func,
     return new ChebyshevTensorFunc(splitValue(func),false);
   else if (type == "chebyshev2")
     return new ChebyshevTensorFunc(splitValue(func),true);
+#ifdef HAS_PYTHON
+  else if (type == "python") {
+    auto pos = func.find_first_of('{');
+    std::string module = func.substr(0,pos-1);
+    std::string params = func.substr(pos);
+    IFEM::cout << "\n\t\tModule = " << module << "\n\t\tParams = " << params << std::endl;
+    return new PythonTensorFunc(module.c_str(), params.c_str());
+  }
+#endif
 
   return nullptr;
 }

src/Utility/PythonFunctions.C

@@ -0,0 +1,168 @@
+// $Id$
+//==============================================================================
+//!
+//! \file PythonFunctions.C
+//!
+//! \date Sep 7 2021
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Python function implementations.
+//!
+//==============================================================================
+
+#include "PythonFunctions.h"
+
+#ifdef HAS_PYTHON
+
+#include <Vec3.h>
+
+#include <pybind11/numpy.h>
+
+std::shared_ptr<InterpreterRAII> pyInterp;
+
+
+static std::shared_ptr<InterpreterRAII> initInterpreter()
+{
+  if (!pyInterp)
+    pyInterp = std::make_shared<InterpreterRAII>();
+
+  return pyInterp;
+}
+
+
+PythonBaseFunc::PythonBaseFunc(const char* module, const char* params) :
+  myInterp(initInterpreter())
+{
+  myModule = pybind11::module::import(module);
+  pybind11::object init = myModule.attr("initialize");
+  myInstance = init(params);
+  eval = myInstance.attr("eval");
+}
+
+
+Real PythonFunc::evaluate (const Real& x) const
+{
+  Real d;
+
+  // python code cannot be called on multiple threads - at least not easily.
+  // for now serialize it
+#pragma omp critical
+  {
+    auto res = eval(x);
+    d = res.cast<double>();
+  }
+  return d;
+}
+
+
+Real PythonFunction::evaluate (const Vec3& x) const
+{
+  const Vec4* Xt = dynamic_cast<const Vec4*>(&x);
+  std::vector<double> data(x.ptr(), x.ptr()+3);
+  if (Xt)
+    data.push_back(Xt->t);
+
+  double d;
+  // python code cannot be called on multiple threads - at least not easily.
+  // for now serialize it
+#pragma omp critical
+  {
+    pybind11::array_t<double> X(data.size(), data.data());
+    auto res = eval(X);
+    d = res.cast<double>();
+  }
+  return d;
+}
+
+
+Vec3 PythonVecFunc::evaluate (const Vec3& x) const
+{
+  const Vec4* Xt = dynamic_cast<const Vec4*>(&x);
+  std::vector<double> data(x.ptr(), x.ptr()+3);
+  if (Xt)
+    data.push_back(Xt->t);
+
+  Vec3 d;
+
+  // python code cannot be called on multiple threads - at least not easily.
+  // for now serialize it
+#pragma omp critical
+  {
+    pybind11::array_t<double> X(data.size(), data.data());
+    auto res = eval(X);
+    size_t i = 0;
+    for (auto item : res) {
+      d[i++] = item.cast<double>();
+      if (i > 2)
+        break;
+    }
+  }
+  return d;
+}
+
+
+Tensor PythonTensorFunc::evaluate (const Vec3& x) const
+{
+  const Vec4* Xt = dynamic_cast<const Vec4*>(&x);
+  std::vector<double> data(x.ptr(), x.ptr()+3);
+  if (Xt)
+    data.push_back(Xt->t);
+
+  Tensor d(3);
+
+  // python code cannot be called on multiple threads - at least not easily.
+  // for now serialize it
+#pragma omp critical
+  {
+    pybind11::array_t<double> X(data.size(), data.data());
+    auto res = eval(X);
+    size_t i = 0;
+    size_t nsd;
+    size_t cmps = ncmp;
+    if (ncmp == 0) {
+      pybind11::array_t<double> f = res.cast<pybind11::array_t<double>>();
+      cmps = f.size();
+      nsd = sqrt(cmps);
+    } else
+      nsd = sqrt(ncmp);
+    d = Tensor(nsd); // funkyness needed to reset dimensionality
+    for (auto item : res) {
+      d(1 + i / nsd, 1 + (i % nsd)) = item.cast<double>();
+      ++i;
+      if (i >= cmps)
+        break;
+    }
+  }
+  return d;
+}
+
+
+SymmTensor PythonSTensorFunc::evaluate (const Vec3& x) const
+{
+  const Vec4* Xt = dynamic_cast<const Vec4*>(&x);
+  std::vector<double> data(x.ptr(), x.ptr()+3);
+  if (Xt)
+    data.push_back(Xt->t);
+
+  auto d = this->Result;
+
+  // python code cannot be called on multiple threads - at least not easily.
+  // for now serialize it
+#pragma omp critical
+  {
+    pybind11::array_t<double> X(data.size(), data.data());
+    auto res = eval(X);
+    size_t i = 0;
+    std::vector<Real>& svec = d;
+    for (auto item : res) {
+      svec[i++] = item.cast<double>();
+      if (i >= svec.size())
+        break;
+    }
+  }
+  return d;
+}
+
+
+#endif