Extend the `define_dirichlet_boundary_condition` function to work with more general function spaces

[matt-graham]

@phillu99 requested for a way to define Dirichlet boundary conditions on more general function spaces - for example in cases using mixed finite element formulations. @phillu99 - would it be possible for you to give an example of what it would be useful for us to support here? I might have not captured what you meant fully - apologies if so!

[Phillu99]

Hi Matt, the problem example is like this:

First I defined a mixed function space VW and extracting subspaces V and W as the following code

mesh = create_unit_square(MPI.COMM_WORLD, nele, nele, CellType.triangle)
V_el = element('CG', mesh.basix_cell(), p_v)
W_el = element('CG', mesh.basix_cell(), p_w)
VW = mixed_element([V_el,W_el])
VW = fem.functionspace(mesh, VW)
V, W = VW.sub(0), VW.sub(1)

Then I want to define different boundary conditions for functions in V and W, but it turns out I cannot do that directly with the subspaces I defined above. They are 'subspace' object rather than usual 'functionspace'. The way I get around this is as following: (in the example I want to define Dirichlet boundaries for V and W on different sides of the geometry)

fdim = mesh.topology.dim -1
face_V = locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[0], 0.0)|np.isclose(x[1],0.0))
face_W = locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[0], 1.0)|np.isclose(x[1],1.0))
P,_ = V.collapse()
Q,_ = W.collapse()
bcf_V,bcf_W = fem.Function(P),fem.Function(Q)
bcf_V.x.array[:] = 0.0
bcf_W.x.array[:] = 0.0
dofs_V = fem.locate_dofs_topological((V, P), fdim, face_V)
dofs_W = fem.locate_dofs_topological((W, Q), fdim, face_W)
bc_V = fem.dirichletbc(value = bcf_V, dofs = dofs_V, V= V)
bc_W = fem.dirichletbc(value = bcf_W, dofs = dofs_W, V= W)
bcs = [bc_V,bc_W]

[Phillu99]

So I'm wondering if it is possible that this dxh function can deal with VW directly with argument that we can indicate what boundary conditions we would like to pose on which subspace

[matt-graham]

Thanks @Phillu99 for the detailed explanation and example - I think what you suggest should be possible. I'll have a go at implementing before we meet next week!

[matt-graham]

Hi @Phillu99. The following adaptation of your code snippet above seems to work already with DOLFINx v0.8

import basix
import dolfinx
import dxh
import numpy as np
from mpi4py import MPI

n_element_per_dimension = 10
element_degrees = [2, 3]

mesh = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD,
    n_element_per_dimension,
    n_element_per_dimension ,
    dolfinx.mesh.CellType.triangle
)

elements = [basix.ufl.element("CG", mesh.basix_cell(), d) for d in element_degrees]
mixed_element = basix.ufl.mixed_element(elements)
mixed_function_space = dolfinx.fem.functionspace(mesh, mixed_element)
sub_function_spaces = [mixed_function_space.sub(i) for i in range(len(elements))]

boundary_conditions = [
    dxh.define_dirichlet_boundary_condition(
        boundary_value=0.,
        function_space=function_space,
        boundary_indicator_function=lambda x: np.isclose(x, 1).any(axis=0)
    )
    for function_space in sub_function_spaces 
]

I'm not sure if the behaviour of sub method has changed in a recent version as mixed_function_space.sub seems to return a FunctionSpace object.

We could still potentially allow for directly passing in mixed_function_space to dxh.define_dirichlet_boundary_condition so something like the following would work

...

mixed_function_space = dolfinx.fem.functionspace(mesh, mixed_element)

boundary_conditions = dxh.define_dirichlet_boundary_condition(
    boundary_value=[0., 0.],
    function_space=mixed_function_space,
    boundary_indicator_function=[
        lambda x: np.isclose(x, 1).any(axis=0),
        lambda x: np.isclose(x, 1).any(axis=0)
    ]
)

with boundary_value accepting an iterable of boundary values for each subspace of mixed_function_space and boundary_indicator_function accepting an iterable indicator functions for each subspace of mixed_function_space. We could potentially also allow a single value to be specified for either argument when the same boundary value / indicator function should be shared across all subspaces.

Would this fit with your requirements?

[Phillu99]

@matt-graham It seems the code with DOLFINX 0.8.0 is not working on my computer. It shows an error "RuntimeError: Cannot tabulate coordinates for a FunctionSpace that is a subspace."

Although when I check the object of mixed_function_space.sub, it also returns a FunctionSpace.

Yes I would think the wrapped up version be very helpful. It seems to me that when dolfinx deals with subspaces, it may not be possible to define boundary conditions for "v" or "w" only when dealing with functions in subspace looking like (v,0) or (0,w), so the collapse method like the following may need to be used.

P,_ = V.collapse()
Q,_ = W.collapse()