diff --git a/cmake/Modules/FindIFEMDeps.cmake b/cmake/Modules/FindIFEMDeps.cmake index 41ea63b4f..da804d7d5 100644 --- a/cmake/Modules/FindIFEMDeps.cmake +++ b/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) diff --git a/cmake/Modules/IFEMOptions.cmake b/cmake/Modules/IFEMOptions.cmake index 34c545e19..37a1d6b42 100644 --- a/cmake/Modules/IFEMOptions.cmake +++ b/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) diff --git a/cmake/Scripts/IFEMTesting.cmake b/cmake/Scripts/IFEMTesting.cmake index 67e6a3808..b1a0c8cb3 100644 --- a/cmake/Scripts/IFEMTesting.cmake +++ b/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 diff --git a/src/Utility/Functions.C b/src/Utility/Functions.C index c9e9a4fb6..e254baea6 100644 --- a/src/Utility/Functions.C +++ b/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; } diff --git a/src/Utility/PythonFunctions.C b/src/Utility/PythonFunctions.C new file mode 100644 index 000000000..14b517993 --- /dev/null +++ b/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 + +#include + +std::shared_ptr pyInterp; + + +static std::shared_ptr initInterpreter() +{ + if (!pyInterp) + pyInterp = std::make_shared(); + + 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(); + } + return d; +} + + +Real PythonFunction::evaluate (const Vec3& x) const +{ + const Vec4* Xt = dynamic_cast(&x); + std::vector 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 X(data.size(), data.data()); + auto res = eval(X); + d = res.cast(); + } + return d; +} + + +Vec3 PythonVecFunc::evaluate (const Vec3& x) const +{ + const Vec4* Xt = dynamic_cast(&x); + std::vector 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 X(data.size(), data.data()); + auto res = eval(X); + size_t i = 0; + for (auto item : res) { + d[i++] = item.cast(); + if (i > 2) + break; + } + } + return d; +} + + +Tensor PythonTensorFunc::evaluate (const Vec3& x) const +{ + const Vec4* Xt = dynamic_cast(&x); + std::vector 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 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 f = res.cast>(); + 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(); + ++i; + if (i >= cmps) + break; + } + } + return d; +} + + +SymmTensor PythonSTensorFunc::evaluate (const Vec3& x) const +{ + const Vec4* Xt = dynamic_cast(&x); + std::vector 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 X(data.size(), data.data()); + auto res = eval(X); + size_t i = 0; + std::vector& svec = d; + for (auto item : res) { + svec[i++] = item.cast(); + if (i >= svec.size()) + break; + } + } + return d; +} + + +#endif diff --git a/src/Utility/PythonFunctions.h b/src/Utility/PythonFunctions.h new file mode 100644 index 000000000..c1eb35dae --- /dev/null +++ b/src/Utility/PythonFunctions.h @@ -0,0 +1,212 @@ +// $Id$ +//============================================================================== +//! +//! \file PythonFunctions.h +//! +//! \date Sep 7 2021 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Python function implementations. +//! +//============================================================================== + +#ifndef _PYTHON_FUNCTIONS_H +#define _PYTHON_FUNCTIONS_H + +#ifdef HAS_PYTHON + +#include + +#include "Function.h" +#include "TensorFunction.h" + +/*! + * \brief RAII style class for managing the interpreter in use. +*/ + +class InterpreterRAII { +public: + //! \brief The constructor allocates the interpreter. + InterpreterRAII() + { + pybind11::initialize_interpreter(); + } + + //! \brief The destructor finalizes the interpreter. + ~InterpreterRAII() + { + pybind11::finalize_interpreter(); + } +}; + + +extern std::shared_ptr pyInterp; + + +/*! + \brief Base class for python module functions. +*/ + +class __attribute__((visibility("hidden"))) PythonBaseFunc { +protected: + //! \brief The constructor sets module to use and its parameters. + //! \param module Name of python module to use + //! \param params Parameters for python module in JSON format + PythonBaseFunc(const char* module, const char* params); + + //! \brief Empty destructor. + virtual ~PythonBaseFunc() = default; + + std::shared_ptr myInterp; //!< Ref-counted pointer to global interpreter + + pybind11::module myModule; //!< Imported python module + pybind11::object myInstance; //!< Instance of module + pybind11::object setParams; //!< Method to set parameters in instance + pybind11::object eval; //!< Method to evaluate instance +}; + + +/*! + \brief A scalar-valued function, python module. +*/ + +class __attribute__((visibility("hidden"))) PythonFunc : + public ScalarFunc, public PythonBaseFunc +{ +public: + //! \brief The constructor sets the module to use and its parameters. + //! \param module Name of python module to use + //! \param params Parameters for python module in JSON format + PythonFunc(const char* function, const char* params) : + PythonBaseFunc(function,params) + {} + + //! \brief The destructor frees the dynamically allocated objects. + virtual ~PythonFunc() = default; + + //! \brief Returns whether the function is time-independent or not. + bool isConstant() const override { return false; } + +protected: + //! \brief Evaluates the function expression. + Real evaluate(const Real& X) const override; +}; + + +/*! + \brief A scalar-valued spatial function, python module. +*/ + +class __attribute__((visibility("hidden"))) PythonFunction : + public RealFunc, public PythonBaseFunc +{ +public: + //! \brief The constructor sets the module to use and its parameters. + //! \param module Name of python module to use + //! \param params Parameters for python module in JSON format + PythonFunction(const char* function, const char* params) : + PythonBaseFunc(function,params) + {} + + //! \brief The destructor frees the dynamically allocated objects. + virtual ~PythonFunction() = default; + + //! \brief Returns whether the function is time-independent or not. + bool isConstant() const override { return false; } + +protected: + //! \brief Evaluates the function expression. + Real evaluate(const Vec3& X) const override; +}; + + +/*! + \brief A vector-valued spatial function, python module. +*/ + +class __attribute__((visibility("hidden"))) PythonVecFunc : + public VecFunc, public PythonBaseFunc +{ +public: + //! \brief The constructor sets the module to use and its parameters. + //! \param module Name of python module to use + //! \param params Parameters for python module in JSON format + PythonVecFunc(const char* function, const char* params) : + PythonBaseFunc(function,params) + {} + + //! \brief The destructor frees the dynamically allocated objects. + virtual ~PythonVecFunc() = default; + + //! \brief Returns whether the function is time-independent or not. + bool isConstant() const override { return false; } + +protected: + //! \brief Evaluates the function expression. + Vec3 evaluate(const Vec3& X) const override; +}; + + +/*! + \brief A tensor-valued spatial function, python module. +*/ + +class __attribute__((visibility("hidden"))) PythonTensorFunc : + public TensorFunc, public PythonBaseFunc +{ +public: + //! \brief The constructor sets the module to use and its parameters. + //! \param n Number of spatial dimensions + //! \param module Name of python module to use + //! \param params Parameters for python module in JSON format + PythonTensorFunc(const char* function, const char* params, size_t n = 0) : + TensorFunc(n), PythonBaseFunc(function,params) + {} + + //! \brief The destructor frees the dynamically allocated objects. + virtual ~PythonTensorFunc() = default; + + //! \brief Returns whether the function is time-independent or not. + bool isConstant() const override { return false; } + +protected: + //! \brief Evaluates the function expression. + Tensor evaluate(const Vec3& X) const override; +}; + + +/*! + \brief A symmetric tensor-valued spatial function, python module. +*/ + +class __attribute__((visibility("hidden"))) PythonSTensorFunc : + public STensorFunc, public PythonBaseFunc +{ +public: + //! \brief The constructor sets the module to use and its parameters. + //! \param n Number of spatial dimensions + //! \param module Name of python module to use + //! \param params Parameters for python module in JSON format + PythonSTensorFunc(size_t n, const char* function, const char* params, bool with33 = false) : + STensorFunc(n,with33), PythonBaseFunc(function,params), Result(n,with33) + {} + + //! \brief The destructor frees the dynamically allocated objects. + virtual ~PythonSTensorFunc() = default; + + //! \brief Returns whether the function is time-independent or not. + bool isConstant() const override { return false; } + +protected: + //! \brief Evaluates the function expression. + SymmTensor evaluate(const Vec3& X) const override; + + SymmTensor Result; //!< Cloned for results +}; + + + +#endif + +#endif diff --git a/src/Utility/Test/TestPythonFunctions.C b/src/Utility/Test/TestPythonFunctions.C new file mode 100644 index 000000000..e9a7421b6 --- /dev/null +++ b/src/Utility/Test/TestPythonFunctions.C @@ -0,0 +1,256 @@ +//============================================================================== +//! +//! \file TestPythonFunctions.C +//! +//! \date Sep 7 2021 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Tests for python functions. +//! +//============================================================================== + +#include "PythonFunctions.h" +#include + +#include "gtest/gtest.h" + + +TEST(TestPythonFunctions, ScalarFunc) +{ + if (!pyInterp) + pyInterp = std::make_shared(); + pybind11::module sys = pybind11::module::import("sys"); + pybind11::list path = sys.attr("path"); + path.append("./src/Utility/Test/refdata"); + double x = 2.0; + double scale = 1.0; + + PythonFunc f("scalar", R"({"scale" : 1.0})"); + EXPECT_DOUBLE_EQ(f(x), scale*x); + + PythonFunc f2("scalar", R"({"scale" : 10.0})"); + x = 3.0; + scale = 10.0; + EXPECT_DOUBLE_EQ(f2(x), scale*x); +} + + +TEST(TestPythonFunctions, RealFunc) +{ + if (!pyInterp) + pyInterp = std::make_shared(); + + pybind11::module sys = pybind11::module::import("sys"); + pybind11::list path = sys.attr("path"); + path.append("./src/Utility/Test/refdata"); + + PythonFunction f("real", R"({"scale" : 1.0})"); + Vec3 X; + const double& x = X[0] = 1.0; + const double& y = X[1] = 2.0; + double scale = 1.0; + + EXPECT_DOUBLE_EQ(f(X), scale*(x + 100.0*y)); + scale = 10.0; + + PythonFunction f2("real", R"({"scale" : 10.0})"); + EXPECT_DOUBLE_EQ(f2(X), scale*(x + 100.0*y)); +} + + +TEST(TestPythonFunctions, VecFunc) +{ + if (!pyInterp) + pyInterp = std::make_shared(); + + pybind11::module sys = pybind11::module::import("sys"); + pybind11::list path = sys.attr("path"); + path.append("./src/Utility/Test/refdata"); + + PythonVecFunc f("vec", R"({"scale" : 1.0})"); + Vec4 X; + const double& x = X[0] = 1.0; + const double& y = X[1] = 2.0; + const double& z = X[2] = 3.0; + const double& t = X.t = 4.0; + double scale = 1.0; + + Vec3 res = f(X); + + EXPECT_DOUBLE_EQ(res[0], scale*t*x); + EXPECT_DOUBLE_EQ(res[1], scale*t*y); + EXPECT_DOUBLE_EQ(res[2], scale*t*z); + + PythonVecFunc f2("vec", R"({"scale" : 10.0})"); + + res = f2(X); + scale = 10.0; + + EXPECT_DOUBLE_EQ(res[0], scale*t*x); + EXPECT_DOUBLE_EQ(res[1], scale*t*y); + EXPECT_DOUBLE_EQ(res[2], scale*t*z); +} + + +TEST(TestPythonFunctions, TensorFunc2D) +{ + if (!pyInterp) + pyInterp = std::make_shared(); + + pybind11::module sys = pybind11::module::import("sys"); + pybind11::list path = sys.attr("path"); + path.append("./src/Utility/Test/refdata"); + + PythonTensorFunc f("tensor", R"({"scale" : 1.0})"); + Vec4 X; + const double& x = X[0] = 1.0; + const double& y = X[1] = 2.0; + const double& t = X.t = 4.0; + double scale = 1.0; + + Tensor res = f(X); + + EXPECT_DOUBLE_EQ(res(1,1), scale*t*x*x); + EXPECT_DOUBLE_EQ(res(1,2), scale*t*x*y); + EXPECT_DOUBLE_EQ(res(2,1), scale*t*y*x); + EXPECT_DOUBLE_EQ(res(2,2), scale*t*y*y); + + PythonTensorFunc f2("tensor", R"({"scale" : 10.0})"); + + res = f2(X); + scale = 10.0; + + EXPECT_DOUBLE_EQ(res(1,1), scale*t*x*x); + EXPECT_DOUBLE_EQ(res(1,2), scale*t*x*y); + EXPECT_DOUBLE_EQ(res(2,1), scale*t*y*x); + EXPECT_DOUBLE_EQ(res(2,2), scale*t*y*y); +} + + +TEST(TestPythonFunctions, TensorFunc3D) +{ + if (!pyInterp) + pyInterp = std::make_shared(); + + pybind11::module sys = pybind11::module::import("sys"); + pybind11::list path = sys.attr("path"); + path.append("./src/Utility/Test/refdata"); + + PythonTensorFunc f("tensor", R"({"scale" : 1.0})"); + Vec4 X; + const double& x = X[0] = 1.0; + const double& y = X[1] = 2.0; + const double& z = X[2] = 3.0; + const double& t = X.t = 4.0; + double scale = 1.0; + + Tensor res = f(X); + + EXPECT_DOUBLE_EQ(res(1,1), scale*t*x*x); + EXPECT_DOUBLE_EQ(res(1,2), scale*t*x*y); + EXPECT_DOUBLE_EQ(res(1,3), scale*t*x*z); + EXPECT_DOUBLE_EQ(res(2,1), scale*t*y*x); + EXPECT_DOUBLE_EQ(res(2,2), scale*t*y*y); + EXPECT_DOUBLE_EQ(res(2,3), scale*t*y*z); + EXPECT_DOUBLE_EQ(res(3,1), scale*t*z*x); + EXPECT_DOUBLE_EQ(res(3,2), scale*t*z*y); + EXPECT_DOUBLE_EQ(res(3,3), scale*t*z*z); + + PythonTensorFunc f2("tensor", R"({"scale" : 10.0})"); + + res = f2(X); + scale = 10.0; + + EXPECT_DOUBLE_EQ(res(1,1), scale*t*x*x); + EXPECT_DOUBLE_EQ(res(1,2), scale*t*x*y); + EXPECT_DOUBLE_EQ(res(1,3), scale*t*x*z); + EXPECT_DOUBLE_EQ(res(2,1), scale*t*y*x); + EXPECT_DOUBLE_EQ(res(2,2), scale*t*y*y); + EXPECT_DOUBLE_EQ(res(2,3), scale*t*y*z); + EXPECT_DOUBLE_EQ(res(3,1), scale*t*z*x); + EXPECT_DOUBLE_EQ(res(3,2), scale*t*z*y); + EXPECT_DOUBLE_EQ(res(3,3), scale*t*z*z); +} + + +TEST(TestPythonFunctions, STensorFunc2D) +{ + if (!pyInterp) + pyInterp = std::make_shared(); + + pybind11::module sys = pybind11::module::import("sys"); + pybind11::list path = sys.attr("path"); + path.append("./src/Utility/Test/refdata"); + + PythonSTensorFunc f(2, "symmtensor", R"({"scale" : 1.0})"); + Vec4 X; + double& x = X[0] = 1.0; + double& y = X[1] = 2.0; + double& t = X.t = 4.0; + double scale = 1.0; + + SymmTensor res = f(X); + + EXPECT_DOUBLE_EQ(res(1,1), scale*t*x*x); + EXPECT_DOUBLE_EQ(res(1,2), scale*t*x*y); + EXPECT_DOUBLE_EQ(res(2,1), scale*t*y*x); + EXPECT_DOUBLE_EQ(res(2,2), scale*t*y*y); + + PythonSTensorFunc f2(2, "symmtensor", R"({"scale" : 10.0, "withtt" : true})", true); + + auto res2 = f2(X); + scale = 10.0; + + EXPECT_DOUBLE_EQ(res2(1,1), scale*t*x*x); + EXPECT_DOUBLE_EQ(res2(1,2), scale*t*x*y); + EXPECT_DOUBLE_EQ(res2(2,1), scale*t*y*x); + EXPECT_DOUBLE_EQ(res2(2,2), scale*t*y*y); + EXPECT_DOUBLE_EQ(res2(3,3), scale*t*3.0); +} + + +TEST(TestPythonFunctions, STensorFunc3D) +{ + if (!pyInterp) + pyInterp = std::make_shared(); + + pybind11::module sys = pybind11::module::import("sys"); + pybind11::list path = sys.attr("path"); + path.append("./src/Utility/Test/refdata"); + + PythonSTensorFunc f(3, "symmtensor", R"({"scale" : 1.0})"); + Vec4 X; + double& x = X[0] = 1.0; + double& y = X[1] = 2.0; + double& z = X[2] = 3.0; + double& t = X.t = 4.0; + double scale = 1.0; + + SymmTensor res = f(X); + + EXPECT_DOUBLE_EQ(res(1,1), scale*t*x*x); + EXPECT_DOUBLE_EQ(res(1,2), scale*t*x*y); + EXPECT_DOUBLE_EQ(res(1,3), scale*t*x*z); + EXPECT_DOUBLE_EQ(res(2,1), scale*t*y*x); + EXPECT_DOUBLE_EQ(res(2,2), scale*t*y*y); + EXPECT_DOUBLE_EQ(res(2,3), scale*t*y*z); + EXPECT_DOUBLE_EQ(res(3,1), scale*t*z*x); + EXPECT_DOUBLE_EQ(res(3,2), scale*t*z*y); + EXPECT_DOUBLE_EQ(res(3,3), scale*t*z*z); + + PythonSTensorFunc f2(3, "symmtensor", R"({"scale" : 10.0})"); + + auto res2 = f2(X); + scale = 10.0; + + EXPECT_DOUBLE_EQ(res2(1,1), scale*t*x*x); + EXPECT_DOUBLE_EQ(res2(1,2), scale*t*x*y); + EXPECT_DOUBLE_EQ(res2(1,3), scale*t*x*z); + EXPECT_DOUBLE_EQ(res2(2,1), scale*t*y*x); + EXPECT_DOUBLE_EQ(res2(2,2), scale*t*y*y); + EXPECT_DOUBLE_EQ(res2(2,3), scale*t*y*z); + EXPECT_DOUBLE_EQ(res2(3,1), scale*t*z*x); + EXPECT_DOUBLE_EQ(res2(3,2), scale*t*z*y); + EXPECT_DOUBLE_EQ(res2(3,3), scale*t*z*z); +} diff --git a/src/Utility/Test/refdata/real/__init__.py b/src/Utility/Test/refdata/real/__init__.py new file mode 100644 index 000000000..cfa21b819 --- /dev/null +++ b/src/Utility/Test/refdata/real/__init__.py @@ -0,0 +1,16 @@ +import numpy as np +from math import pi, cos +import json + +def initialize(params): + return SourceTerm(json.loads(params)) + +class SourceTerm: + def __init__(self, params): + self._scale = params['scale'] + + def setParams(self, mu): + _mu = mu + + def eval(self, X): + return self._scale * (X[0] + X[1]*100.0) diff --git a/src/Utility/Test/refdata/scalar/__init__.py b/src/Utility/Test/refdata/scalar/__init__.py new file mode 100644 index 000000000..7dfeb6d68 --- /dev/null +++ b/src/Utility/Test/refdata/scalar/__init__.py @@ -0,0 +1,16 @@ +import numpy as np +from math import pi, cos +import json + +def initialize(params): + return SourceTerm(json.loads(params)) + +class SourceTerm: + def __init__(self, params): + self._scale = params['scale'] + + def setParams(self, mu): + _mu = mu + + def eval(self, x): + return self._scale * x diff --git a/src/Utility/Test/refdata/symmtensor/__init__.py b/src/Utility/Test/refdata/symmtensor/__init__.py new file mode 100644 index 000000000..93d108180 --- /dev/null +++ b/src/Utility/Test/refdata/symmtensor/__init__.py @@ -0,0 +1,31 @@ +import numpy as np +from math import pi, cos +import json + +def initialize(params): + return SourceTerm(json.loads(params)) + +class SourceTerm: + def __init__(self, params): + self._scale = params['scale'] + try: + self._with33 = params["withtt"] + except: + self._with33 = False + + def setParams(self, mu): + _mu = mu + + def eval(self, X): + x = X[0] + y = X[1] + z = X[2] + t = X[3] + if z == 0.0: + if self._with33: + return self._scale * t * np.array([x*x, y*y, 3.0, x*y]) + else: + return self._scale * t * np.array([x*x, y*y, x*y]) + else: + return self._scale * t * np.array([x*x, y*y, z*z, x*y, y*z, x*z]) + diff --git a/src/Utility/Test/refdata/tensor/__init__.py b/src/Utility/Test/refdata/tensor/__init__.py new file mode 100644 index 000000000..c4037991f --- /dev/null +++ b/src/Utility/Test/refdata/tensor/__init__.py @@ -0,0 +1,23 @@ +import numpy as np +from math import pi, cos +import json + +def initialize(params): + return SourceTerm(json.loads(params)) + +class SourceTerm: + def __init__(self, params): + self._scale = params['scale'] + + def setParams(self, mu): + _mu = mu + + def eval(self, X): + x = X[0] + y = X[1] + z = X[2] + t = X[3] + if z == 0.0: + return self._scale * t * np.array([1.0*x, 2.0*x, 1.0*y, 2.0*y]) + else: + return self._scale * t * np.array([1.0*x, 2.0*x, 3.0*x, 1.0*y, 2.0*y, 3.0*y, 1.0*z, 2.0*z, 3.0*z]) diff --git a/src/Utility/Test/refdata/vec/__init__.py b/src/Utility/Test/refdata/vec/__init__.py new file mode 100644 index 000000000..40ad82d20 --- /dev/null +++ b/src/Utility/Test/refdata/vec/__init__.py @@ -0,0 +1,16 @@ +import numpy as np +from math import pi, cos +import json + +def initialize(params): + return SourceTerm(json.loads(params)) + +class SourceTerm: + def __init__(self, params): + self._scale = params['scale'] + + def setParams(self, mu): + _mu = mu + + def eval(self, X): + return self._scale * X[3] * np.array([X[0], X[1], X[2]])