Dirichlet Boundary Conditions in Mini-Elements

Question

Hello everyone,

I am trying to use mini-elements to solve the Stokes flow problem (here is an example in the FEniCS documentation: link), and I would like to apply a free-slip boundary condition in my domain. Basically, it is a Dirichlet boundary condition applied in only one dimension in the vector function space.

Generally, you can apply the free-slip boundary condition like this:

scalar = FunctionSpace(mesh, "CG", 1)
vector = VectorFunctionSpace(mesh, "CG", 1)
system = vector*scalar
...
BC = DirichletBC(system.sub(0).sub(1), Constant(0.0), boundary_parts, 1)

Unfortunately, this does not give a stable solution to the Stokes flow problem and that’s why I’m trying the mini-element. My setup now looks like this:

scalar = FunctionSpace(mesh, "CG", 1)
vector = VectorFunctionSpace(mesh, "CG", 1)
bubble = VectorFunctionSpace(mesh, "Bubble", 3)
system = (vector+bubble)*scalar
...
BC = DirichletBC(system.sub(0).sub(1), Constant(0.0), boundary_parts, 1)

This gives me the following error: ValueError: Can only extract SubSpaces with i = 0 ... -1. It seems that the vector+bubble function space cannot be separated into subspaces for each spatial dimension, and thus I cannot apply my free-slip boundary condition. Does anybody know of a solution or a workaround?

Thanks in advance.

Answer 1

You could try to impose the boundary condition weakly with a Nitsche-type method.

Something similar to the following formulation may work. $u\cdot n = 0$ is weakly enforced on the boundary, and implicitly the normal derivative of the tangential component also vanishes.

P1 = FiniteElement("Lagrange", triangle, 1)
B3 = FiniteElement("Bubble", triangle, 3)

V_element = (P1+B3)*(P1+B3)
Q_element = P1

W = FunctionSpace(mesh, V_element * Q_element)

u, p = TrialFunctions(W)
v, q = TestFunctions(W)

h = CellSize(mesh)
n = FacetNormal(mesh)
g = Constant(0.)

gamma = Constant(10.)

# standard stokes forms
a = (inner(grad(u), grad(v)) + p * div(v) + q * div(u)) * dx
L = inner(f, v) * dx

# integration by parts terms
a -= (inner(grad(u)*n, n) * inner(v, n) + inner(v, n*p)) * ds

# symmetrize
a -= (inner(grad(v)*n, n) * inner(u, n) + inner(u, n*q)) * ds
L -= (inner(grad(v)*n, n) + q)* g * ds

# stabilize
a += gamma * h**-1 * inner(u, n) * inner(v, n) * ds
L += gamma * h**-1 * g * inner(v, n) * ds

edit: fixed formulation, it was imposing noslip instead free-slip.