# Inhomogeneous Poisson equation with user defined expression

Hi guys, I am writing code for the inhomogeneous Poisson equation - \nabla \cdot K(x) \nabla u = f, where K(x) is a piecewise constant function (defined in the code below).

The issue I encountered is that the assembling of the weak form is extremely time-consuming. I believe it is due to the way I defined the expression for K(x). It’s probably like iterating over a NumPy array using a for loop, which is very slow in python.

Here is the minimal working code. I made a time comparison between running assembling (\nabla u, \nabla v) and (K(x) \nabla u, \nabla v ) in the code.

from fenics import *
import time
import numpy as np
import math

# Define Layer structure dielectric constant
dielectric_consts = [1,50,1,50,1,50]
class MyExpression(UserExpression):
def eval(self, values, x):
num_layer = len(dielectric_consts)
width = 1.0/num_layer

if x[1] <= width:
values[0] = dielectric_consts[0]
elif x[1] <= 2*width:
values[0] = dielectric_consts[1]
elif x[1] <= 3*width:
values[0] = dielectric_consts[2]
elif x[1] <= 4*width:
values[0] = dielectric_consts[3]
elif x[1] <= 5*width:
values[0] = dielectric_consts[4]
else:
values[0] = dielectric_consts[5]
def value_shape(self):
return() #scalar function

mesh = UnitSquareMesh(400,400)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("100*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.01)", degree=2)

# Define weak form 1
K = MyExpression()
a_1 = inner(K*grad(u), grad(v))*dx # Assmbling is super slow with K, not optimized.
L_1 = f*v*dx

# Define weak form 2

begin = time.time()
A, b = assemble_system(a_1, L_1,None)
end = time.time() - begin
print("With K ... time taken: ",end)

begin = time.time()
A, b = assemble_system(a_2, L_1,None)
end = time.time() - begin
print("Without K ... time taken: ",end)


Any ideas on how to optimize this in FEniCS?

1 Like

This issue was addressed in the book Solving PDEs in Python
– The FEniCS Tutorial Volume 1, 4.3.3 Chaper 4.

The solution is to use C++ code and pass it to python expression.

Update:
However, I am not able to compile to code after using C++ code in python. Still need help

To reproduce the error, please run the following code. Thank you!

from fenics import *

cppcode = """
class DielectricConstant : public Expression
{
void eval(Array<double>& values, const Array<double>& x) const
{
int num_layer =  dielectric.size();
double width = 1.0/num_layer;
for (int i = 0; i < num_layer; i++){
if( (x[1] >= i*width) && (x[1] <(i+1)*width) ){
values[0] = dielectric[i];
break;
}
}
}
public:
std::vector<double> dielectric {1,50,1,50,1,50};
};
"""
mesh = UnitSquareMesh(200,200)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("100*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.01)", degree=1)

K_cpp = Expression(cpp_code = cppcode,degree = 1)

# Define weak form