Assume you have some mathematical function (such as cylindrical Bessel functions) from some C/C++ library (like the special function package from boost), and you want to use it within CoefficientFunctions. Now you can import it via a small, separate Python module. You can (but don't have to) export several versions of the function for real or complex evaluation, AVX acceleration, and you can also provide its derivative if you want to use it together with automatic differentiation. Note that Python is used only to build up the expression tree, the actual evaluation stays within the C++ code, possibly evaluated in parallel.
// ngscxx -c myfunctions.cpp ; ngsld -shared myfunctions.o -lngfem -lngstd -o myfunctions.so #include <fem.hpp> #include <python_ngstd.hpp> using namespace ngfem; // function class with overloaded ()-operator struct MyGenericSin { // real version double operator() (double x) const { return sin(x); } // complex version Complex operator() (Complex x) const { return sin(x); } // support automatic differntiation AutoDiff<1> operator() (AutoDiff<1> x) const { AutoDiff<1> result; result.Value() = sin(x.Value()); result.DValue(0) = cos(x.Value())*x.DValue(0); return result; } // speedup using vector functions (or dummy function instead) SIMD<double> operator() (SIMD<double> x) const { return SIMD<double>([&](int i)->double { return sin(x[i]); }); } static string Name() { return "mysin"; } template <typename T> T operator() (T x) const { throw Exception(ToString("MyGenericSin not implemented for type") + typeid(T).name()); } }; // function class with overloaded ()-operator struct MyGenericAtan2 { double operator() (double x, double y) const { return atan2(x,y); } template <typename T, typename Ty> T operator() (T x, Ty y) const { throw Exception(ToString("MyGenericAtan2 not implemented for type") + typeid(T).name()); } }; // new python module providing functions PYBIND11_PLUGIN(myfunctions) { cout << "loading myfunctions" << endl; py::module m("myfunctions", "pybind myfunctions"); ExportUnaryFunction<MyGenericSin> (m, "mysin"); ExportBinaryFunction<MyGenericAtan2> (m, "myatan2"); return m.ptr(); }
Its use from Python is shown here
from ngsolve import * from netgen.geom2d import unit_square from myfunctions import mysin mesh = Mesh(unit_square.GenerateMesh(maxh=0.1)) V = H1(mesh, order=3) u = GridFunction(V) u.Set (mysin(6.28*x)) Draw (u) print (mysin (0.2)) print (mysin(x))