Boundary conditions for two different variables

1

Hi all,

I am trying to implement a time-dependent boundary condition for velocity and pressure.
The element space for the velocity is order 2, whereas the element space for pressure is order 1.

Please, find below a minimal example:

import dolfinx
from dolfinx.io import XDMFFile
from dolfinx.fem import (Constant, DirichletBC, Function, FunctionSpace, apply_lifting, assemble_matrix, 
                         assemble_scalar, assemble_vector, create_vector, locate_dofs_topological, set_bc)
import gmsh
import ufl
from ufl import inner, div, dx
from mpi4py import MPI
import numpy as np
from petsc4py import PETSc

# Domain
gdim = 3
mesh = dolfinx.UnitCubeMesh(MPI.COMM_WORLD, 10, 10, 10)

# Velocity element space
v_cg2 = ufl.VectorElement('CG', mesh.ufl_cell(), 2)
V = dolfinx.FunctionSpace(mesh, v_cg2)
u = Function(V)
# Pressure element space
p_cg1 = ufl.FiniteElement('CG', mesh.ufl_cell(), 1)
Q = dolfinx.FunctionSpace(mesh, p_cg1)
p = Function(Q)
# Mixed element space
mel = ufl.MixedElement([v_cg2, p_cg1])
Y = dolfinx.FunctionSpace(mesh,mel)
q = Function(Y)

# Define boundary function (pulse function)
t = 0
class BoundaryCondition():
    def __init__(self, t):
        self.t = t

    def __call__(self, x):
        values = np.zeros((gdim, x.shape[1]),dtype=PETSc.ScalarType)
        if t == 0:
            values[0] = 1
        else:
            values[0] = 0
        return values

# Get dofs of the left plane
def left(x):
    return np.isclose(x[0], 0)
left_dofsu = dolfinx.fem.locate_dofs_geometrical(V, left) # Velocity
left_dofsp = dolfinx.fem.locate_dofs_geometrical(Q, left) # Pressure

# Define boundary condition for velocity
u_left = Function(V)
left_velocity = BoundaryCondition(t)
u_left.interpolate(left_velocity)
bcu = dolfinx.DirichletBC(u_left, left_dofsu)
# Define boundary condition for pressure
p_left = Function(Q)
left_velocity = BoundaryCondition(t)
p_left.interpolate(left_velocity)
bcp = dolfinx.DirichletBC(p_left, left_dofsp)
# Boundary condition
bc = [bcu,bcp]
dolfinx.fem.set_bc(q.vector, [bc])

If I just consider the velocity, the code works perfectly and I can visualise the correct value of the velocity in the left boundary. However, when I try to implement the same condition for the velocity and pressure using the previous code, I get the following error:

Traceback (most recent call last):
File “/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/function.py”, line 301, in _interpolate
self._cpp_object.interpolate(u._cpp_object)
AttributeError: ‘BoundaryCondition’ object has no attribute ‘_cpp_object’
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File “/root/dolfinX/Test2.py”, line 57, in
p_left.interpolate(left_pressure)
File “/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/function.py”, line 309, in interpolate
_interpolate(u)
File “/usr/lib/python3.9/functools.py”, line 877, in wrapper
return dispatch(args[0].__class__)(*args, **kw)
File “/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/function.py”, line 303, in _interpolate
self._cpp_object.interpolate(u)
RuntimeError: Interpolation data has the wrong shape.

Does anybody know what I am doing wrong?
I guess that the problem is due to the different function spaces that I am using for the velocity and pressure.

Thanks in advance!

2

Doesn’t

values = np.zeros((gdim, x.shape[1]),dtype=PETSc.ScalarType)

return a 3x1 vector, which is fine for the velocity, however for the pressure you need a scalar?

Therefore adding a second boundary condition specifically for the pressure could solve your issue:

class BoundaryCondition_Pressure():
    def __init__(self, t):
        self.t = t

    def __call__(self, x):
        values = np.zeros((1, x.shape[1]),dtype=PETSc.ScalarType)
        if t == 0:
            values[0] = 1
        else:
            values[0] = 0
        return values
3

Yes, you are absolutely right!

Many thanks