Enriched element issue

variational formulation
1

Hi all,

Now that I understand my previous issue Lid-driven cavity simulation issue (thanks to Jack Hale) I’m trying to look at my problem with an enriched element.

One of my options is to use the (P2 + Bubble3) * P1 element :

P2 = FiniteElement("CG", mesh.ufl_cell(), 2)
B = FiniteElement("B", mesh.ufl_cell(), 3)
P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
P2B = P2 + B

element = MixedElement([VectorElement(P2B, dim = 2), P1])

W = FunctionSpace(mesh, element)

When I use this element it raises a new warning :

*** FFC warning: evaluate_dof(s) for enriched element not implemented.

Is there any way to bypass this problem ?

Any help will be appreciated !

Here is the full code using the latest Fenics version

from dolfin import *

n = 50
Re= 100
nu = Constant(1.0 / Re)
T = 5
steps = 500
dt = T / steps
t = 0.0

ufile = File("velocity.pvd")
pfile = File("pressure.pvd")

mesh = UnitSquareMesh(n, n)

def MarkDomain(mesh):
    boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
    boundary_parts.set_all(0)

    class Inflow(SubDomain):
        def inside(self, x, on_boundary):
            return near(x[1], 1.0) and on_boundary
    inflow = Inflow()
    inflow.mark(boundary_parts, 1)

    class Walls(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and (near(x[0], 0.0) or near(x[1], 0.0) or near(x[0], 1.0))
    walls = Walls()
    walls.mark(boundary_parts, 2)

    return boundary_parts

P2 = FiniteElement("CG", mesh.ufl_cell(), 2)
B = FiniteElement("B", mesh.ufl_cell(), 3)
P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
P2B = P2 + B

element = MixedElement([VectorElement(P2B, dim = 2), P1])

W = FunctionSpace(mesh, element)

v, q = TestFunctions(W)
w = Function(W)
u, p = split(w)
w_prec = Function(W)
u_prec, v_prec = split(w_prec)

boundary_parts = MarkDomain(mesh)

lid = DirichletBC(W.sub(0), Constant((1.0,0.0)), boundary_parts, 1)
noslip = DirichletBC(W.sub(0), Constant((0.0,0.0)), boundary_parts, 2) 
bc = [noslip, lid]

class InitialConditions(UserExpression):
    def eval(self, values, x):
        values[0] = 0.0
        values[1] = 0.0
        values[2] = 0.0
    def value_shape(self):
        return (3,)
w.interpolate(InitialConditions())

F = (inner(u - u_prec, v) / dt + inner(grad(u) * u_prec, v) + nu * inner(grad(u), grad(v)) - inner(p, div(v)) - inner(q, div(u))) * dx

for step in range(1, steps):
    print("Step = ", step)
    assign(w_prec, w)
    t += dt
    solve(F == 0, w, bc)
    u, p = w.split(deepcopy = True)
    ufile << u
    pfile << p
2

Hi, maybe you can try ‘NodalEnrichedElement’. More information can be found: