How do I use a Python function to define Dirichlet boundary conditions which are not known in closed form?

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! :slight_smile:

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.

Have a look at UserExpressions, for instance:

Here, you can substitute self.T with your function, which has to be evaluated at x[0],x[1],x[2]

3 Likes