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.