I am implementing the solution to an elliptic PDE for which the diffusion coefficient has a different expression in different subdomains. In each subdomain, the expression of the coefficient is quite involved, so I need to define a class inheriting from UserExpression to define it. Maybe I should also mention that I have defined the subdomains using flags (PhysicalSurface) in gmsh. So far, I have implemented my subdomain-dependent coefficient as follows:
I define a user-defined expression using a class, for each expression that I need to define (so, if for example I have two subdomains, each with its own expression, I define two classes for the coefficient expressions)
I use the gmsh flags to define indicator functions which are 1 in a subdomain and 0 in the others; to define these I use functions in the finite element space ‘DG’ of degree 0
my PDE coefficient is the sum of the expressions from point 1. each multiplied by a coefficient from point 2.
For example, if I have 2 subdomains, I define two expressions expr1 and expr2 and two piecewise constant coefficients coeff1 and coeff2 as indicator functions, and my resulting PDE coefficient is
expr1*coeff1 + expr2*coeff2
Now, this works but the code is quite slow, and I expect it is because each of the expressions expr1 and expr2 is assembled on the whole domain.
Is there a more efficient way to define subdomain-dependent expressions?
I have seen the post https://fenicsproject.org/qa/5364/setting-different-expressions-different-subdomains-python/, but first I am not sure that going back and forth submeshes and extracting the mesh coordinates would be more efficient (not mentioning that maybe fenics uses a higher order quadrature rule), second, in principle I do not need access to the node coordinates and using an expression would be sufficient for my purposes.
The UserExpression class is slow in general, because it requires execution of Python code. Some ways to implement more complicated logic in JIT-compiled C++ Expressions can be found in this thread and another one (linked by my answer in the first).
not mentioning that maybe fenics uses a higher order quadrature rule
Note that you can manually set the quadrature degree used by FEniCS, to override its (sometimes excessive) default choice, e.g.,
dx = dx(metadata={"quadrature_degree":n})
(and likewise for ds and dS), where n is an integer giving the degree of polynomial up to which integration is exact.
Maybe I’m not fully understanding what you are trying to do, but if you are using the indicator functions to restrict the expressions, wouldn’t it be easier to just use ufl syntax for your expressions, e.g. like
x = SpatialCoordinate(mesh)
expr1 = sqrt(x[0]*x[1])
expr2 = x[0]*x[1]**2
I am not sure I fully get what you mean, but my expressions are quite complicated (I need to compute a small matrix and I use this matrix and its inverse to define matrix-valued coefficients), so this is why I am using a classes define the expressions.
My question was more if there is a way to tell fenics that one expression has to be computed on one subdomain only. But you hint is really helpful, thank you! I will try that and come back if I have issues. Just a short questions: at the moment, I am using another expression to the constructor of my expressions for the coefficients. It should be possible to do the same also in the C++ syntax, right?
Sorry, it’s me again… I installed pybind and I tried to recreate the example that you linked about the subdomain. I get the error
ImportError: generic_type: type “OmegaOut” referenced unknown base type “dolfin::SubDomain”
Which I realized is the same problem that people got here (I tried to recreate that example too, with the same problem). Do you know if there is some progress on this issue or if I should change something compared to the subdomain example that you linked? I have dolfin version 2019.1.0.
You shouldn’t need to manually install pybind11. (FEniCS wouldn’t work without it.) The code in the SubDomain example I linked is spread between several listings, which might be causing confusion. The following complete Python script runs without error for me, using version 2019.1.0:
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)
omega_out_code = """
#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
#include <math.h>
namespace py = pybind11;
#include <dolfin/function/Function.h>
#include <dolfin/mesh/SubDomain.h>
class OmegaOut : public dolfin::SubDomain
{
public:
// Create expression with 3 components
OmegaOut() : dolfin::SubDomain() {}
// Function for evaluating expression on each cell
bool inside(Eigen::Ref<const Eigen::VectorXd> x, bool on_boundary) const override
{
Eigen::Array<double, 1, 1> values;
phi.eval(values, x);
return values[0] > 0.5;
}
dolfin::Function phi;
};
PYBIND11_MODULE(SIGNATURE, m)
{
py::class_<OmegaOut, std::shared_ptr<OmegaOut>, dolfin::SubDomain>
(m, "OmegaOut")
.def(py::init<>())
.def_readwrite("phi", &OmegaOut::phi);
}
"""
c = compile_cpp_code(omega_out_code).OmegaOut()
c.phi = phi.cpp_object()
# (Change phi_expr definition to just "x[0]" for simpler testing.)
print(c.inside([0.75,0.5],False))
Thank you for resuming the code. I tried this version and still I get the same error
Traceback (most recent call last):
File “test.py”, line 91, in
c = compile_cpp_code(omega_out_code).OmegaOut()
File “/usr/lib/python3/dist-packages/dolfin/jit/pybind11jit.py”, line 96, in compile_cpp_code generate=jit_generate)
File “/usr/lib/python3/dist-packages/dolfin/jit/jit.py”, line 47, in mpi_jit
return local_jit(*args, **kwargs)
File “/usr/lib/python3/dist-packages/dolfin/jit/jit.py”, line 103, in dijitso_jit
return dijitso.jit(*args, **kwargs)
File “/usr/lib/python3/dist-packages/dijitso/jit.py”, line 212, in jit
lib = load_library(signature, cache_params)
File “/usr/lib/python3/dist-packages/dijitso/cache.py”, line 363, in load_library lib = import(signature)
ImportError: generic_type: type “OmegaOut” referenced unknown base type “dolfin::SubDomain”
I am now concerned that I messed up pybind by getting it from pip3, but then according to what you say also Fenics should not work, and I tried one of my codes and it runs ok. I don’t think it is relevant, but I am running Fenics on a virtual machine. It seems to be a problem of my system anyways.
We spotted a similar error involving pybind11 few months ago. See this topic.
I’m also working on a virtual machine (Ubuntu 18.04LTS) with python virtual environments (venv). I still see this error when I try to define a new expression in C++ with a new name or when I’m not working in the usual virtual environment.
As far as I know no fix has been proposed yet.
Hi, thank you for the useful info. For the moment I can at least compile the code on other machine where it works, but let’s hope they will fix the problem.
Hi, I managed to implement the cpp code as you suggested, using the cell indices to distinguish the subdomains. It worked and the code got 100 times faster! Thank you!!