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
a_2 = inner(grad(u), grad(v))*dx
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?