Featured

## Exporting mathematical functions

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)
{
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))

Attachments
use myfunctions.py [252B]
Uploaded Wednesday, 11 January 2017 by Joachim Schoeberl
myfunctions.cpp [1.61Kb]
Uploaded Tuesday, 21 March 2017 by Joachim Schoeberl