Hi all, I am attempting to simulate a chaotic advection equation where the velocity vector field is based on the Arnold–Beltrami–Childress (ABC) flow equation: v_x = A*sin(z) + C*cos(y) v_x = B*sin(x) + A*cos(z) v_x = C*sin(y) + B*cos(x) Right now I have the following: from firedrake import * mesh = UnitCubeMesh(20,20,20) A = 0.08 B = 0.1 C = 0.15 V = VectorFunctionSpace(mesh, "CG", 1) velocity = Function(V) velocity.interpolate(Expression(( "A*sin(2*pi*x[2])+C*cos(2*pi*x[1])", "B*sin(2*pi*x[0])+A*cos(2*pi*x[2])", "C*sin(8*pi*x[1])+B*cos(8*pi*x[0])"), A=A,B=B,C=C)) Where A,B,C are user-defined constants (will add such constants inside all sin and cos functions later). I want to have DirichletBC's only applied to faces and regions where the velocity points into the cube (i.e., the "upstream" location). The region locations and sizes will vary depending on the constants chosen. Yes I could create separate GMSH files and manually determine where these "inflow" regions will be for every set of A/B/C/etc constants, but I was wondering if there's a way to automate this the Firedrake way, particularly if these constants will vary as a function of time. Thanks, Justin