Apply DirichletBCs to vector space, with direction normal to on mesh boundary

[gnikit]

Hi, I am having trouble coming up with a programmatic way to compute Dirichlet BCs normal to a mesh region, for a vector function space. I am guessing that the input arguments to the such a Function should be the vector space V, the subdomain id i, the vector magnitude m and potentially the mesh (to fetch the normals n of i). You would then interpolate m*n onto V and pass it to the DirichletBC, e.g.

n = fd.FacetNormal(mesh)
my_normal_direction = fd.Function(V).interpolate(fd.Constant((m,)) * n)
DirichletBC(V, g=my_normal_direction, [i])

Clearly that does not work. I suspect the limitation is my knowledge of firedrake and ufl syntax, rather than the actual capabilities of the libraries.

NOTE: This is a dummy test case, but if the surface is curved e.g. a pipe's outer diameter, specifying BCs in this manner is quite useful.

Any help would be greatly appreciated.

MWE - Stokes flow

import firedrake as fd

N = 10
mesh = fd.UnitSquareMesh(N, N)

V = fd.VectorFunctionSpace(mesh, "CG", 2)
W = fd.FunctionSpace(mesh, "CG", 1)
Z = V * W

u, p = fd.TrialFunctions(Z)
v, q = fd.TestFunctions(Z)

a = (
    fd.inner(fd.grad(u), fd.grad(v)) - fd.inner(p, fd.div(v)) + fd.inner(fd.div(u), q)
) * fd.dx
L = fd.inner(fd.Constant((0, 0)), v) * fd.dx

# TODO: programmatically calc the normal vector to subdomain 1 with mag 1.2
# n = fd.FacetNormal(mesh)
# norm_to_mesh = - fd.Function(Z.sub(0)).interpolate(fd.Constant((1.2,)) * n)
# bc_inlet = fd.DirichletBC(Z.sub(0), norm_to_mesh, [1])
bc_inlet = fd.DirichletBC(Z.sub(0), fd.Constant((1.2, 0)), [1])
bc_walls = fd.DirichletBC(Z.sub(0), fd.Constant((0, 0)), [3, 4])
bcs = [bc_inlet, bc_walls]

up = fd.Function(Z)

# Use direct solver
fd.solve(
    a == L,
    up,
    bcs=bcs,
    solver_parameters={
        "ksp_type": "gmres",
        "mat_type": "aij",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
    },
)

[dham]

There isn't a straightforward way to impose boundary conditions only on the normal component of an H1 vector function space such as CG2. This is because it's not straightforward to decompose degrees of freedom into normal and tangent degrees of freedom. There are two ways of achieving what you are trying to do. The nicer one is to use an H(div)-L2 discretisation of your fluid. The H(div) function spaces only have normal degrees of freedom on the boundary, so normal Dirichlet conditions just work. The less nice one is to use a Nitsche style (penalty) boundary condition to penalise normal flow.