Mixed element with a single element

[jorgensd]

Summarize the issue

Since #3500 one can no longer make a mixed function space with a single element.
I see no sensible justification for this, as it massively complicates packages that uses an arbitrary number of mixed elements (which can be >=1).
First reported at festim-dev/FESTIM#922 (comment)

How to reproduce the bug

Run the below on nightly

Minimal Example (Python)

from mpi4py import MPI
import dolfinx

import basix.ufl


mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10)
el = basix.ufl.element("Lagrange", mesh.basix_cell(), 1)
me = basix.ufl.mixed_element([el])

V = dolfinx.fem.functionspace(mesh, me)

Output (Python)

Traceback (most recent call last):
  File "/root/shared/mwe_mixed.py", line 11, in <module>
    V = dolfinx.fem.functionspace(mesh, me)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/function.py", line 630, in functionspace
    cpp_element = _create_dolfinx_element(mesh.topology.cell_type, ufl_e, dtype)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/function.py", line 578, in _create_dolfinx_element
    return CppElement(elements)
           ^^^^^^^^^^^^^^^^^^^^
RuntimeError: FiniteElement constructor for mixed elements called with a single element.

Version

main branch

DOLFINx git commit

No response

Installation

No response

Additional information

No response

[nate-sime]

This change has been quite frustrating for some of my applications also. Particularly for mixed domain problems, e.g., HDG-like formulations.

[RemDelaporteMathurin]

We would like this too in our applications (see festim-dev/FESTIM#922 (comment))

[garth-wells]

The title is misleading; the underlying wish seems to be the ability to create an element from a list of one element. The 'mixed' part is not relevant.

See resolution in FEniCS/basix#882. Not yet tested with DOLFINx.

[jorgensd]

The title is misleading; the underlying wish seems to be the ability to create an element from a list of one element. The 'mixed' part is not relevant.

See resolution in FEniCS/basix#882. Not yet tested with DOLFINx.

The mixed part is relevant.
From a high level perspective image you want to couple N>=1 equations in a domain. You can set up relations between each of these equations in an algorithmic fashion. One would then rely on .sub functionality for boundary conditions and interpolation.

Removing the capability of a mixed element with one component means that one needs to branch code, as a normal element cannot do sub(0) and return itself (it might even return a component of the element if it is blocked).

[RemDelaporteMathurin]

@garth-wells I hope this explains the problem a bit better:

from mpi4py import MPI
import dolfinx
from dolfinx import fem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.fem.petsc import NonlinearProblem
import basix.ufl
import ufl

nb_elements = 2  # assume that this can be an arbitrary number >=1


mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_connectivity(fdim, tdim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)

el = basix.ufl.element("Lagrange", mesh.basix_cell(), 1)

# first check
if nb_elements == 1:
    me = el
else:
    me = basix.ufl.mixed_element([el] * nb_elements)

V = dolfinx.fem.functionspace(mesh, me)

u = dolfinx.fem.Function(V)

# second check
if nb_elements == 1:
    v = ufl.TestFunction(V)
    sub_test_functions = [v]
    sub_functions = [u]
else:
    v = ufl.TestFunctions(V)
    sub_test_functions = list(v)
    sub_functions = list(ufl.split(u))

F = 0
f = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(-1))
for i in range(nb_elements):
    u_i = sub_functions[i]
    v_i = sub_test_functions[i]
    F += ufl.inner(ufl.grad(u_i), ufl.grad(v_i)) * ufl.dx
    F += f * v_i * ufl.dx

bcs = []

for i in range(nb_elements):
    # make a function
    # third check
    if nb_elements > 1:
        function_space_bc_function, _ = V.sub(i).collapse()
    else:
        function_space_bc_function = V

    u_bc = dolfinx.fem.Function(function_space_bc_function)
    u_bc.interpolate(lambda x: x[0])

    # find the dofs on the boundary
    if nb_elements > 1:  # fourth check
        V_collapsed, _ = V.sub(i).collapse()
        function_space_dofs = (V.sub(i), V_collapsed)
    else:
        function_space_dofs = V

    bc_dofs = fem.locate_dofs_topological(function_space_dofs, fdim, boundary_facets)

    # make the BC form
    if nb_elements == 1:  # fourth check
        function_space_form = None
    else:
        function_space_form = V.sub(i)
    form = fem.dirichletbc(
        value=u_bc,
        dofs=bc_dofs,
        V=function_space_form,
    )
    bcs.append(form)

# solve
dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
problem = NonlinearProblem(F, u, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)

solver.solve(u)

Right now in dolfinx 0.9.0 we can do the following instead, without all these checks:

from mpi4py import MPI
import dolfinx
from dolfinx import fem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.fem.petsc import NonlinearProblem
import basix.ufl
import ufl

nb_elements = 1  # assume that this can be an arbitrary number >=1


mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_connectivity(fdim, tdim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)

el = basix.ufl.element("Lagrange", mesh.basix_cell(), 1)

me = basix.ufl.mixed_element([el] * nb_elements)

V = dolfinx.fem.functionspace(mesh, me)

u = dolfinx.fem.Function(V)


v = ufl.TestFunctions(V)
sub_test_functions = list(v)
sub_functions = list(ufl.split(u))

F = 0
f = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(-1))
for i in range(nb_elements):
    u_i = sub_functions[i]
    v_i = sub_test_functions[i]
    F += ufl.inner(ufl.grad(u_i), ufl.grad(v_i)) * ufl.dx
    F += f * v_i * ufl.dx

bcs = []

for i in range(nb_elements):
    u_bc = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0))
    sub_V = V.sub(i)

    bc_dofs = fem.locate_dofs_topological(
        sub_V,
        fdim,
        boundary_facets,
    )

    form = fem.dirichletbc(
        value=u_bc,
        dofs=bc_dofs,
        V=sub_V,
    )
    bcs.append(form)

# solve
dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
problem = NonlinearProblem(F, u, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)

solver.solve(u)

Does this help understanding the need? Let me know if you need more info

[jorgensd]

Just to comment on the above example by Rémi.
That particular example can be done with a blocked element as el = basix.ufl.element("Lagrange", mesh.basix_cell(), 1) is identical for all species. However, this will in general not be the case, as different concentrations can live in different spaces.

[garth-wells]

Thanks @RemDelaporteMathurin for the clear use case. I'm happy for dolfinx::fem::FiniteElement to allow elements with one sub-element to support this use case.

Over time we can work towards a better design. The underlying need in the use case is for Functions to be iterable. It would be nicer to to able to something like:

e = basix.ufl.element([el] * nb_elements)
V = dolfinx.fem.functionspace(mesh, e)
u = dolfinx.fem.Function(V)
v = ufl.TestFunctions(V)
for ui, vi in zip(u, v):
    F += ufl.inner(ufl.grad(ui), ufl.grad(vi)) * ufl.dx
    F += f * vi * ufl.dx

for i in range(len(e)):
    _e = e[i]

and to remove the concept of a 'mixed element'.

[RemDelaporteMathurin]

Being able to iterate through functions and test functions like this would be useful, but we need to be careful and treat split(u) and u.split() differently, no? (One for creating the form the other for post processing)

[RemDelaporteMathurin]

treat split(u) and u.split() differently

Or rather split(u) and u.sub(i)