What is the correct mixed element formulation?

General
1

I have to echo some frustration that is shown in this post → FEniCS version and compatibility problems (chaos?)

What is the correct way of dealing with mixed elements?

Is is like in @dokken 's Demo for Mixed formulation for the Poisson equation?

Q_el = element("BDMCF", msh.basix_cell(), k, dtype=default_real_type)
P_el = element("DG", msh.basix_cell(), k - 1, dtype=default_real_type)
V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(msh, V_el)

This method is also shown in Mixed finite element problems — FEniCS Workshop …

el_u = basix.ufl.element("Lagrange", mesh.basix_cell(), 3, shape=(mesh.geometry.dim,))
el_p = basix.ufl.element("Lagrange", mesh.basix_cell(), 2)
el_mixed = basix.ufl.mixed_element([el_u, el_p])

Or is it as @jpdean states in his GitHub - jpdean/mixed_domain_demos at 65975c7e073fff490bcd9a3809f9234c6c66b771 " Mixed-domain functionality in FEniCSx is under active development and is not ready for production use."

Or should it be done as shown here (FEniCS release notes | FEniCS Project)

V = dolfinx.fem.functionspace(domain, ("Lagrange", 2))
Q = dolfinx.fem.functionspace(submesh, ("Lagrange", 1))

# Create mixed problem residual F
W = ufl.MixedFunctionSpace(V, Q)

I am constantly running into errors when I try to use mixed elements. For example (MWE):

from dolfinx import fem, mesh
from mpi4py import MPI
from basix.ufl import element, mixed_element

nx, ny = 5, 5
domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.triangle)

U   = element("DG", domain.basix_cell(), 2, shape=(2, )) # velocity
V   = element("DG", domain.basix_cell(), 2, shape=(3, )) # deformation
X = fem.functionspace(domain, mixed_element([U,V]))


def interpolation_function(x):
        return 1 + x[0]**2 + x[1]**2

source = fem.Function(X)

# LINE BLOW THROWS ERROR:
# RuntimeError: Cannot get interpolation points - no Basix element available. Maybe this is a mixed element?
source.interpolate(interpolation_function)

I am just trying to port @garth and Steven Vandekerckhove’s Automatic calibration of damping layers in finite element time domain simulations from FEniCS to FEniCSx, but it feels like trying to hit a moving target.

If the answer is that FEniCSx isn’t ready for mixed element problems, that would be very helpful to know.

Thank you for your help:

VERSION INFORMATION

either tutorial compatable docker images

2

Hello,

The best place to find the most up-to-date and recommended ways of doing things are the demos in the official DOLFINx repository here. FEniCSx is capable of solving mixed-element problems (see, for example, the mixed Poisson demo or the Stokes demo).

The comment in my repository about mixed-domain functionality refers to a separate feature that is unrelated to mixed elements.

Hope this clears things up!

Thanks,

Joe

1 Like
3

Thanks @jpdean that sort of clears things up. Although I find it strange that the FEniCS release notes suggest ufl.MixedfunctionSpace() while the mixed Poisson demo we both linked to uses

V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(msh, V_el)

I am correct in saying the release notes are less up to date then the dolfinx demos?

Also, when I run that particular Poisson demo (in a Jupyter notebook with the version information in the original post) I get the following output when viewed in Paraview:

Screenshot-2024-12-10_10:22:48
Screenshot-2024-12-10_10:22:483746×1898 213 KB

Just a red square with u = 1 everywhere. Is that expected?

Finally, is it possible to interpolate a mixed function space in FEniCSx?

4

You can interpolate into each sub component of a mixed space.

Here you are trying to interpolate a scalar value into a mixed space with in total 6 components.

Say you want to interpolate into the U space. Then you can write something like:

source.sub(0).interpolate(lambda x: (x[0], x[1]))

to interpolate (x, y) into the velocity space.
For the deformation space, that has three components, you can call

source.sub(1).interpolate(lambda x: (np.sin(x[1]), np.cos(x[0]), x[1])))

to interpolate (sin(y), cos(x), y) into the deformation space.

1 Like
5

Thank you @dokken . That will hopefully get me a bit further along.

I’ve noticed that the mixed Poisson demo discussed above hasn’t been updated since the v0.9.0 release notes. Could you speak to if I should be using

W = ufl.MixedFunctionSpace(V, Q)

or

V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(msh, V_el)

thanks again

6

It is not 1 everywhere. This is an issue introduced in later versions of “superlu_dist”.

Changing it to:

problem = LinearProblem(
    a,
    L,
    bcs=bcs,
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
        "ksp_error_if_not_converged": True,
    },
)

resolves the issue and you get

image
image1736×1391 21.1 KB

I’ve made a pull request adding these changes to the demo to make it clearer in the future: Change solver in mixed Poisson by jorgensd · Pull Request #3558 · FEniCS/dolfinx · GitHub

1 Like
7

Thanks again @dokken. I made the change and that fixed the mixed Poisson demo, just as you showed.

I do still have the lingering doubt as to whether I should use the mixed element formulation shown in this demo, or the one shown in the v0.9.0 release notes.

Since the release notes are more up to date, I would assume I should use ufl.MixedFunctionSpace() but this goes against @jpdean 's advice. (the demo was updated 3 months ago and the release notes are 2 months old, perhaps that is causing the confusion?)

8

This is no longer a true statement,
Mixed domain is stabilized and well tested.

For mixed element formulations on a single grid, I would probably stick with basix.ufl.mixed_element over ufl.MixedFunctionSpace.

MixedFunctionSpace is especially neat when integrating over a mesh and a submesh.

1 Like
9

Thank you for the very helpful information.

I realize it must be difficult to maintain an open source project like this. After spending a little bit of time here, I realize that you and @jpdean are really working hard to try to get people up to speed, and I very much appreciate it. I just wish Simula would hire 10 more people like you to help maintain this great project and make it more viable for a user like me.