It seems to me that the interpolation of UserExpression is very costly as soon as if -statements come into play. For instance, I am using the following implementation for a characteristic function of some set of class dolfin.SubDomain :
class G(dolfin.UserExpression):
def __init__(self, SubDomain, degree = 2, **kwargs):
try: super(G, self).__init__(**kwargs)
except: pass
self.degree = degree
self.SubDomain = SubDomain
def eval(self, value, x):
if self.SubDomain.inside(x, True):
value[0] = 1.0
else:
value[0] = 0.0
def value_shape(self):
return ()
Is there a better way to implement a characteristic function such that the interpolation into a dolfin.FunctionSpace is not so costly? It is important for me that I obtain a dolfin.Function . Please find attached a small MWE which compares the interpolation with and without the if-statement. Any help is very much appreciated.
import dolfin
import time
class G(dolfin.UserExpression):
def __init__(self, SubDomain, degree = 2, **kwargs):
try: super(G, self).__init__(**kwargs)
except: pass
self.degree = degree
self.SubDomain = SubDomain
def eval(self, value, x):
if self.SubDomain.inside(x, True):
value[0] = 1.0
else:
value[0] = 0.0
def value_shape(self):
return ()
class G_fast(dolfin.UserExpression):
def __init__(self, SubDomain, degree = 2, **kwargs):
try: super(G_fast, self).__init__(**kwargs)
except: pass
self.degree = degree
self.SubDomain = SubDomain
def eval(self, value, x):
value[0] = 1.0
def value_shape(self):
return ()
class Omega(dolfin.SubDomain):
def inside(self, x, on_boundary):
radius = dolfin.sqrt(x[0]**2 + x[1]**2)
return radius <= 0.2
if __name__ == "__main__":
nx = ny = 100
mesh = dolfin.UnitSquareMesh(nx,ny)
V = dolfin.FunctionSpace(mesh, "DG", 0)
omega = Omega()
g = G(omega)
g_fast = G_fast(omega)
start_1 = time.process_time()
g_inter = dolfin.interpolate(g,V)
end_1 = time.process_time()
start_2 = time.process_time()
g_fast_inter = dolfin.interpolate(g_fast, V)
end_2 = time.process_time()
print('slow: {:5.10f}s'.format(end_1 - start_1))
print('fast: {:5.10f}s'.format(end_2 - start_2))
plugged
November 26, 2019, 12:49pm
2
Maybe take a look at these, and try to implement your code in C++, this will make it much faster than using a UserExpression, that uses python to do the computations
Thanks Kamensky!
Just to add: Using the Eigen::Ref<> template instead of py::array_t<> , pybind11 casts numpy arrays directly and also gives access to the Eigen library methods.
i.e Change npArray definition to
typedef Eigen::Ref<Eigen::VectorXd> npArray;
arr declaration to npArray arr;
and the constructor method to
test_exp(npArray a) : dolfin::Expression(), arr(a) {}
This can help in using Eigen Methods e.g arr.sum() in the eval method.
My goal is to speed up the evaluation of the SubDomain Gamma
from dolfin import *
mesh = UnitSquareMesh(50,50, 'left/right')
PHI = FunctionSpace(mesh, 'CG', 1)
phi_expr = Expression('2.0*(1.0 - x[0])*x[0] + 0.2 - x[1]', degree=2)
phi = interpolate(phi_expr , PHI)
File("phi_initial.pvd") << phi
class Gamma(SubDomain):
def inside(self, x, on_boundary):
return phi(x) > 0.5
with a C++ expression as follows
omega_out_code = """
#include <pybin…
1 Like
Thank you for the links. It really helped me a lot! However, the code I put together still gives me an Error that I cannot resolve. I guess, the problem is the pybind for cpp_code_charFunc . Since I get the following error:
ImportError: generic_type: type "CharFunc" is already registered!
Please find below the code. Any further help is very much appreciated.
import dolfin
import numpy as np
cpp_code_charFunc = """
#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
namespace py = pybind11;
#include <dolfin/function/Expression.h>
#include <dolfin/function/Function.h>
#include <dolfin/function/Constant.h>
#include <dolfin/mesh/SubDomain.h>
class CharFunc : public dolfin::Expression
{
public:
std::shared_ptr<dolfin::Constant> jc;
std::shared_ptr<dolfin::SubDomain> omega;
CharFunc(std::shared_ptr<dolfin::Constant> jc_, std::shared_ptr<dolfin::SubDomain> omega_) : dolfin::Expression(){
jc = jc_;
omega = omega_;
}
void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x) const override
{
if (omega->inside(x, 1>0)){
values[0] = jc->values()[0];
}
else{
values[0] = 0.0;
}
}
};
PYBIND11_MODULE(SIGNATURE, m)
{
py::class_<CharFunc, std::shared_ptr<CharFunc>, dolfin::Expression>(m, "CharFunc")
.def(py::init<std::shared_ptr<dolfin::Constant>, std::shared_ptr<dolfin::SubDomain> >())
.def_readwrite("jc", &CharFunc::jc)
.def_readwrite("omega", &CharFunc::omega);
}
"""
cpp_code_Omega = """
#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
namespace py = pybind11;
#include <dolfin/function/Function.h>
#include <dolfin/mesh/SubDomain.h>
class Omega : public dolfin::SubDomain
{
public:
// Create expression with 3 components
Omega() : dolfin::SubDomain() {}
// Function for evaluating expression on each cell
bool inside(Eigen::Ref<const Eigen::VectorXd> x, bool on_boundary) const override
{
return sqrt( pow(x[0], 2) + pow(x[1], 2)) < 0.2;
}
};
PYBIND11_MODULE(SIGNATURE, m)
{
py::class_<Omega, std::shared_ptr<Omega>, dolfin::SubDomain>
(m, "Omega")
.def(py::init<>());
}
"""
mesh = dolfin.UnitSquareMesh(10,10)
V = dolfin.FunctionSpace(mesh,"DG",0)
omega = dolfin.compile_cpp_code(cpp_code_Omega).Omega()
#omega = Omega()
jc = dolfin.Constant(2.0)
g = dolfin.compile_cpp_code(cpp_code_charFunc)
expr = dolfin.CompiledExpression(g.CharFunc(jc, omega), degree=1)
expr_inter = dolfin.interpolate(expr, V)
dolfin.File("expr.pvd") << expr_inter
I resolved the problem myself. Replace the line
with
expr = dolfin.CompiledExpression(g.CharFunc(jc.cpp_object(), omega), degree=1)