Error with 1d mesh immersed in 2d with dolfinx

I’m trying dolfinX to generate a line mesh in 2d, where I want to define a form depending on stretching and curvature (double derivative of vertical displacement). I don’t know the meaning of the error message I’m getting, code is as follows

from dolfinx import mesh, fem, io
from mpi4py import MPI
import ufl
import numpy as np
import gmsh

comm = MPI.COMM_WORLD
gdim = 2

omega_0 = 1
boundary0 = 2

gmsh.initialize()
if comm.rank == 0:
    gmsh.model.add("telaio")

    factory = gmsh.model.geo

    h = 0.1

    square_points = [
        factory.addPoint(0.0, 0.0, 0.0, h),
        factory.addPoint(0.0, 1.0, 0.0, h),
    ]

    square_lines = [
        factory.addLine(square_points[0], square_points[1]),
    ]

    square_curve = factory.addCurveLoop(square_lines)

    factory.synchronize()

    gmsh.model.addPhysicalGroup(1, square_lines, omega_0)
    gmsh.model.addPhysicalGroup(0, square_points, boundary0)

    gmsh.model.mesh.generate(1)

partitioner = mesh.create_cell_partitioner(mesh.GhostMode.shared_facet)
msh, ct, ft = io.gmshio.model_to_mesh(
    gmsh.model, comm, 0, gdim=gdim, partitioner=partitioner
)
gmsh.finalize()

dx = ufl.Measure("dx", domain=msh, subdomain_data=ct)

DEGREE = 2
Ue = ufl.FiniteElement("CG", msh.ufl_cell(), DEGREE)
Um = ufl.MixedElement([Ue, Ue])
V_0 = fem.FunctionSpace(msh, Um)

# Test and trial functions
u_0 = ufl.TrialFunction(V_0)
v_0 = ufl.TestFunction(V_0)

w0, v0 = ufl.split(u_0)
w0t, v0t = ufl.split(v_0)

eps0, th0 = ufl.grad(w0)[1], ufl.grad(v0)[1]
eps0t, th0t = ufl.grad(w0t)[1], ufl.grad(v0t)[1]
chi0, chi0t = ufl.grad(th0)[1], ufl.grad(th0t)[1]


a = (
    ufl.inner(chi0, chi0t) * dx
    + ufl.inner(eps0, eps0t) * dx
    )    

a = fem.form(a)

and results in the error

Info    : Meshing 1D...
Info    : Meshing curve 1 (Line)
Info    : Done meshing 1D (Wall 0.000160543s, CPU 0.000146s)
Info    : 11 nodes 12 elements
Traceback (most recent call last):
  File "/root/shared/prova.py", line 69, in <module>
    a = fem.form(a)
  File "/usr/local/dolfinx-real/lib/python/dist-packages/dolfinx/fem/forms.py", line 194, in form
    return _create_form(form)
  File "/usr/local/dolfinx-real/lib/python/dist-packages/dolfinx/fem/forms.py", line 189, in _create_form
    return _form(form)
  File "/usr/local/dolfinx-real/lib/python/dist-packages/dolfinx/fem/forms.py", line 152, in _form
    ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
  File "/usr/local/dolfinx-real/lib/python/dist-packages/dolfinx/jit.py", line 56, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/local/dolfinx-real/lib/python/dist-packages/dolfinx/jit.py", line 204, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 199, in compile_forms
    raise e
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 190, in compile_forms
    impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 260, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/compiler.py", line 102, in compile_ufl_objects
    ir = compute_ir(analysis, object_names, prefix, options, visualise)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/ir/representation.py", line 198, in compute_ir
    irs = [_compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
  File "/usr/local/lib/python3.10/dist-packages/ffcx/ir/representation.py", line 198, in <listcomp>
    irs = [_compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
  File "/usr/local/lib/python3.10/dist-packages/ffcx/ir/representation.py", line 493, in _compute_integral_ir
    integral_ir = compute_integral_ir(itg_data.domain.ufl_cell(), itg_data.integral_type,
  File "/usr/local/lib/python3.10/dist-packages/ffcx/ir/integral.py", line 84, in compute_integral_ir
    mt_table_reference = build_optimized_tables(
  File "/usr/local/lib/python3.10/dist-packages/ffcx/ir/elementtables.py", line 382, in build_optimized_tables
    t = get_ffcx_table_values(quadrature_rule.points, cell,
  File "/usr/local/lib/python3.10/dist-packages/ffcx/ir/elementtables.py", line 144, in get_ffcx_table_values
    tbl = tbl[basix_index(derivative_counts)]
IndexError: index 3 is out of bounds for axis 0 with size 3

1 Like

I’ve made an issue at: Double derivative on manifold · Issue #575 · FEniCS/ffcx · GitHub

Fixed in: Resolve double derivatives for manifolds by jorgensd · Pull Request #576 · FEniCS/ffcx · GitHub, feel free to try the branch

2 Likes

Thanks @dokken, I tried the script above on the nightly docker image and now it runs without errors

1 Like