Ulf issue for higher order derivatives on quadratic elements

Dear community,

I’m giving a following question to the following post : link.
The problem was solved for a second order derivative.
Here is my mwe that provides the same error

import gmsh
from mpi4py import MPI
from dolfinx.io import gmshio

from ufl import (TestFunction, TrialFunction,
                 Measure,
                 FacetNormal,
                 inner, grad)
from basix.ufl import element
from dolfinx.fem import functionspace, form, petsc

def generate_cube_gmshio(side_length=1.0, filename="cube.msh"):
    """
    Generate a cube with side length `side_length`,
    export the mesh to `filename`, and
    return (final_mesh, cell_tags, facet_tags).
    """
    # Initialize the Gmsh API
    gmsh.initialize()

    # Retrieve the MPI communicator and rank
    comm = MPI.COMM_WORLD
    model_rank = 0

    # Reference the Gmsh model
    model = gmsh.model
    model.add("cube")

    # Create the cube with the OpenCASCADE kernel
    # (x, y, z, dx, dy, dz)
    box_tag = model.occ.addBox(0, 0, 0, side_length, side_length, side_length)

    # Synchronize the geometry
    model.occ.synchronize()

    # (Optional) add a physical group for the volume (tag=1)
    model.addPhysicalGroup(3, [box_tag], tag=1)

    # Generate the 3D mesh with second-order elements
    gmsh.option.setNumber("Mesh.ElementOrder", 2)
    gmsh.option.setNumber("Mesh.HighOrderOptimize", 2)
    model.mesh.generate(3)

    # Save the mesh to a .msh file
    gmsh.write(filename)

    # Convert the Gmsh mesh into a FEniCSx mesh
    final_mesh, cell_tags, facet_tags = gmshio.model_to_mesh(
        model, comm, model_rank
    )

    # Finalize Gmsh
    gmsh.finalize()

    # Return the objects describing the mesh in FEniCSx
    return final_mesh, cell_tags, facet_tags

# Example usage:
if __name__ == "__main__":
    # Generate a cube mesh and retrieve mesh data
    final_mesh, cell_tags, facet_tags = generate_cube_gmshio(1.0, "cube.msh")

    # Define volume (dx) and surface (ds) measures
    dx = Measure("dx", domain=final_mesh)
    ds = Measure("ds", domain=final_mesh)

    # Define a Lagrange element of degree 3
    P1 = element("Lagrange", final_mesh.basix_cell(), 3)
    P = functionspace(final_mesh, P1)

    # Define trial and test functions
    p, v = TrialFunction(P), TestFunction(P)

    # Define the volumetric bilinear form
    k = inner(grad(p), grad(v)) * dx

    # Define the facet normal
    n = FacetNormal(final_mesh)

    # dp/dn = grad(p) · n
    dp = inner(grad(p), n)

    # d^2p/dn^2 = grad(dp/dn) · n = grad(grad(p) · n) · n
    ddp = inner(grad(dp), n)

    # Third derivative along the normal
    dddp = inner(grad(ddp), n)

    # Define the boundary integral
    c = inner(dddp, v) * ds

    # Assemble the total form
    A = form(k + c)

Providing the following error :

File "/opt/anaconda3/envs/env3-10-complex/lib/python3.10/site-packages/ufl/algorithms/apply_derivatives.py", line 680, in reference_grad
    raise ValueError("ReferenceGrad can only wrap a reference frame type!")
ValueError: ReferenceGrad can only wrap a reference frame type!

This error happens only when I use high order elements, and for the third normal derivatives (not before).

i’ve explored many ways to express this third order derivatives to go around this problem, but I need to use this and I don’t know what to do now…

Thanks,
Pierre.

So we can reduce this to a “pure” ufl issue by the following:

from ufl import (
    FunctionSpace,
    Mesh,
    TestFunction,
    TrialFunction,
    tetrahedron,
    dx,
    grad,
    inner,
    FacetNormal,
    ds,
)
from ufl.algorithms import compute_form_data
from ufl.finiteelement import FiniteElement
from ufl.pullback import identity_pullback
from ufl.sobolevspace import H1

el = FiniteElement("Lagrange", tetrahedron, 3, (), identity_pullback, H1)
domain = Mesh(FiniteElement("Lagrange", tetrahedron, 2, (3,), identity_pullback, H1))
P = FunctionSpace(domain, el)


# Define trial and test functions
p, v = TrialFunction(P), TestFunction(P)

# Define the volumetric bilinear form
k = inner(grad(p), grad(v)) * dx

# Define the facet normal
n = FacetNormal(domain)

# dp/dn = grad(p) · n
dp = inner(grad(p), n)

# d^2p/dn^2 = grad(dp/dn) · n = grad(grad(p) · n) · n
ddp = inner(grad(dp), n)

# Third derivative along the normal
dddp = inner(grad(ddp), n)

# Define the boundary integral
c = inner(dddp, v) * ds


A = compute_form_data(k + c)

returning

  File "/root/shared/mwe.py", line 45, in <module>
    A = compute_form_data(k + c)
        ^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/compute_form_data.py", line 278, in compute_form_data
    form = preprocess_form(form, complex_mode)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/compute_form_data.py", line 244, in preprocess_form
    form = apply_derivatives(form)
           ^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/apply_derivatives.py", line 1568, in apply_derivatives
    dexpression_dvar = map_integrand_dags(rules, expression)
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py", line 86, in map_integrand_dags
    return map_integrands(
           ^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py", line 30, in map_integrands
    map_integrands(function, itg, only_integral_type) for itg in form.integrals()
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py", line 39, in map_integrands
    return itg.reconstruct(function(itg.integrand()))
                           ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/map_integrands.py", line 87, in <lambda>
    lambda expr: map_expr_dag(function, expr, compress), form, only_integral_type
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py", line 35, in map_expr_dag
    (result,) = map_expr_dags(
                ^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py", line 103, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/apply_derivatives.py", line 1371, in grad
    return map_expr_dag(rules, f, vcache=self.vcaches[key], rcache=self.rcaches[key])
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py", line 35, in map_expr_dag
    (result,) = map_expr_dags(
                ^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/corealg/map_dag.py", line 101, in map_expr_dags
    r = handlers[v._ufl_typecode_](v)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/dolfinx-env/lib/python3.12/site-packages/ufl/algorithms/apply_derivatives.py", line 680, in reference_grad
    raise ValueError("ReferenceGrad can only wrap a reference frame type!")
ValueError: ReferenceGrad can only wrap a reference frame type!

Indeed, changing the domain to a linear one removes the error.

Reducing it further, you will see that it is an issue with the double derivative of the facet normal.

from ufl import (
    Mesh,
    tetrahedron,
    FacetNormal,
    ds,
)
from ufl.algorithms import compute_form_data
from ufl.finiteelement import FiniteElement
from ufl.pullback import identity_pullback
from ufl.sobolevspace import H1

domain = Mesh(FiniteElement("Lagrange", tetrahedron, 2, (3,), identity_pullback, H1))

# Define the facet normal
n = FacetNormal(domain)

# Define the boundary integral
c = n[0].dx(0).dx(0) * ds


A = compute_form_data(c)

On a linear grid this is 0.
On a curved grid it is non-zero.
I’ll investigate further.

UFL issue created at: Double derivative of Facet normal broken for higher order meshes · Issue #353 · FEniCS/ufl · GitHub

Dear colleague,
Thanks for the energy spent on this subject.

Has it been merged in the main ? If I update my conda env would I have access to the updated version ?

Best regards,
Pierre.

It has not been merged yet, unfortunately.

It will hopefully be in the 0.10 release of DOLFINx, but will not land on conda prior to that.

All right thanks for the answer !