Hello all! Hope you’re having a good one.

To illustrate my question I will use the 2d Poisson equation as a toy example, and will use the example code, which I also include here:

```
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)
# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx
# Compute solution
u = Function(V)
solve(a == L, u, bc)
# Plot solution and mesh
plot(u)
plt.show()
```

Suppose now, however, that we do not know the closed form of the Dirichlet boundary conditions, but that we have a python function

```
def py_fun(x,y):
return 1 + x**2 + 2*y**2
```

(this is for illustration purposes only, here we assume that in general we do **not** know the closed form of py_fun(x,y). This function may give some interpolation from numerically calculated data, e.g.)

How would we change the

```
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
```

part in order to implement the boundary conditions to be taken from py_fun(x,y)?

This topic might be related, but it only imposes constant boundary conditions.

Any input will be very helpful!

Thanks!

EDIT: So, as per dokken’s suggestion, writing

```
class MyExpression0(UserExpression):
def eval(self, value, x):
value[0] = py_fun(x[0],x[1])
u_D = MyExpression0()
```

Does the job.