Finite element spaces on manifolds do not carry information about their value_shape in physical space

[jorgensd]

Summarize the issue

The value check for interpolating on manifolds are wrong. Introduced in: https://github.com/FEniCS/dolfinx/pull/3500.

How to reproduce the bug

The following example runs on v0.9.0, but not nightly

Minimal Example (Python)

from mpi4py import MPI
import numpy as np
import basix.ufl
import dolfinx
import ufl

points = np.array([[0.0, 0.0, 0.0], [1.0, 0, 1.0], [1.0, 1.0, 0.0]])
cells = np.array([[0, 1, 2]], dtype=np.int32)
ud = ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(3,)))
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, points, ud)


el = basix.ufl.element("RT", mesh.basix_cell(), 2)
V = dolfinx.fem.functionspace(mesh, el)
w = dolfinx.fem.Function(V)


el_o = basix.ufl.element("DG", mesh.basix_cell(), 1, shape=(mesh.geometry.dim,))
Q = dolfinx.fem.functionspace(mesh, el_o)
q = dolfinx.fem.Function(Q)
q.interpolate(w)

Output (Python)

File "/root/shared/examples/eikonal/debug_interpolate_rt_manifold.py", line 21, in <module>
    q.interpolate(w)
  File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/function.py", line 460, in interpolate
    _interpolate(u0)
  File "/usr/lib/python3.12/functools.py", line 909, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/function.py", line 446, in _
    self._cpp_object.interpolate(u0._cpp_object, cells0, cells1)  # type: ignore
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Interpolation: elements have different value dimensions

Version

main branch

DOLFINx git commit

d45b65f6fff295e575f774385b8172054f36d908

Installation

Running DOLFINx nightly container

Additional information

No response

[michalhabera]

Is it just the check that is wrong or is the case of elements with non-scalar basis on manifolds broken?

[jorgensd]

It is in general broken. even a simpler example, such as

from mpi4py import MPI
import numpy as np
import ufl
import basix.ufl
import dolfinx

points = np.array([[0.0, 0.0, 0.0], [1.0, 0, 1.0], [1.0, 1.0, 0.0]])
cells = np.array([[0, 1, 2]], dtype=np.int32)
ud = ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(3,)))
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, points, ud)


el = basix.ufl.element("RT", mesh.basix_cell(), 2)
V = dolfinx.fem.functionspace(mesh, el)
w = dolfinx.fem.Function(V)
w.interpolate(lambda x: np.asarray([x[0], x[1], x[2]]))

fails on nightly

[jorgensd]

This might also be because the error handling with the "new" value size is wrong.

[michalhabera]

This might also be because the error handling with the "new" value size is wrong.

My gut feeling is that the difference between reference_value_shape and value_shape got lost, i.e. it now only exists for symmetric elements and is encoded in the symmetry keyword and block size. But tdim vs. gdim is not taken into account. Could it be?

[jorgensd]

This might also be because the error handling with the "new" value size is wrong.

My gut feeling is that the difference between reference_value_shape and value_shape got lost, i.e. it now only exists for symmetric elements and is encoded in the symmetry keyword and block size. But tdim vs. gdim is not taken into account. Could it be?

Yes, that seems to be the issue. Currently, it seems like reference_value_shape and value_shape only differs for blocked , quadrature and mixed elements:

/// @brief Value shape of the base (non-blocked) finite element field. /// /// For non-blocked elements, this function returns the same as /// FiniteElement::value_shape. For blocked and quadrature elements /// the returned shape will be {}. /// /// Mixed elements do not have a reference value shape. /// /// @throws Exception is thrown for a mixed element as mixed elements /// do not have a value shape. /// @return The value shape. std::span reference_value_shape() const;

which is not correct on manifolds. @garth-wells what do you think is the best way to restore manifold behavior?

[garth-wells]

The difference between reference_value_shape and value_shape is not documented or explained anywhere in the code. This means bugs because it's not clear which one should be used.

If the difference can be properly explained, supported by an example, we can then look at a fix.

[jorgensd]

The difference between reference_value_shape and value_shape is not documented or explained anywhere in the code. This means bugs because it's not clear which one should be used.

If the difference can be properly explained, supported by an example, we can then look at a fix.

Reference value shape for RT on a triangular manifold is (2,), while a physical value shape would be (3,).

This is equivalent to the ufl_shape of a function from the space.

In Python, function space has a value_shape, which is what one would like.

I feel like all of the examples I have given above illustrates what value shape was/should be. Do you need further clarification?

[garth-wells]

Thanks for the examples. Can you give a math explanation too? My initial expectation would be that the value shape remains (2,), i.e. remain in the plane of the element. How does it get from 2 to 3?

[jorgensd]

Thanks for the examples. Can you give a math explanation too? My initial expectation would be that the value shape remains (2,), i.e. remain in the plane of the element. How does it get from 2 to 3?

The degrees of freedom of an RT space relates to the normal at the facet of the manifold. This normal will have three components. You can for instance see this when we want to interpolate a vector field onto a manifold. Then we need to supply a 3D field, which is then «projected» into the normal plane on each element.

This is at least how value shape is interpreted in ufl https://github.com/FEniCS/ufl/blob/1ab30c2de2a975b9fc957c6b356541bd15b35e15/ufl/pullback.py#L160

[jorgensd]

Thinking even more about this, it’s clear that then name value_shape links to the shape of the basis values when pushed to the physical space. Thus it is natural that a RT basis function on a curved manifold (triangular) has value_shape (3,)

[michalhabera]

Thats correct. Another way to see this is that the RT element transforms using contravariant Piola $P = (1/|\det J|) J$ (this will be the push-forward that maps reference basis to the physical configuration). Since Jacobian $J_{ij} = \frac{\partial x_i}{\partial X_j}$ where $x_i$ are physical and $X_i$ reference coordinates, it is a rectangular matrix $J \in \mathbb R^{\text{gdim} \times {tdim}}$. So it maps (tdim, ) vector basis to (gdim, ) vectors in the physical configuration.

How value_shape differs from reference_value_shape is thus encoded in the type of element transformation, gdim and tdim. Symmetric elements have also custom logic for (physical) value_shape, https://github.com/FEniCS/ufl/blob/main/ufl/pullback.py#L525

I think a proper fix is to pass reference_value_shape and value_shape when constructing elements in dolfinx, or get it from basix if available. Any logic to try to deduce these will go against the idea basix as a single source of information about elements.

[garth-wells]

The data shouldn't go into Basix. Each class/object should carry the minimum required information otherwise we get into a mess. It should be set in DOLFINx.

[michalhabera]

@garth-wells Do you see a plan for a fix? I just want to stress once again that deducing value_shape from reference_value_shape, pullback, tdim and gdim is a nontrivial assumption. AFAIK the logic is already buried in https://github.com/FEniCS/basix/blob/main/cpp/basix/maps.h, so basix has that information, it just does not store it explicitly.

I also do not understand why UFL still contains the definition of pullbacks, I thought one would like to provide them in the concrete implementation of the element in basix. @mscroggs Do you know? This is error prone, since one inserts UFL-defined pullbacks but the implementation uses basix definition of the pullback. Meaning wrong results of these two mismatch.

[garth-wells]

Pass the reference value shape to the fem::FiniteElement constructor. And document well. fem::FiniteElement has a private member:

std::optional<std::vector<std::size_t>> _reference_value_shape;

[jorgensd]

Do we have any resolution to this? It is quite frustrating that manifolds are kind of broken.