# Second deriative on a submesh

Dear Community,

I have an issue with the second derivative of a trial function over a submesh. My project is around 10k code lines, and it is not really easy to explain you the context in a simple post.
But basically, I have a differential equation on a global domain where an integral over a surface appears, in this equation p is the main function and q is the secondary function that appears under the integral over the surface, this surface belongs to the global domain, and is also considered as a subdomain. I have another equation relative to p and q that is only defined along the subdomain. So far, I’ve implemented the order 1 and 2 of my operator succesfuly, at the order 3, long story short, I need to have the second derivate of q along the subdomain, and that’s where problems appear.
Here is the system for the 2nd order, I didn’t take the time to write the 3rd one with latex

I’ve build a code as small as possible to reproduce the error, this code doesn’t correspond to the problem, but just reaches the error:

# Needed modules

``````from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio
import dolfinx.mesh as msh
import numpy as np
from dolfinx import plot
from dolfinx.fem import FunctionSpace, form
from ufl import Measure, TrialFunction, TestFunction, grad, inner
``````

# Geometry with gmsh

``````gmsh.initialize()
comm = MPI.COMM_WORLD
model_rank = 0
model = gmsh.model()

side_box = 1
lc = 1e-1
# Definition of the points
p1 = gmsh.model.geo.addPoint(0, 0, 0, lc)
p2 = gmsh.model.geo.addPoint(side_box, 0, 0, lc)
p3 = gmsh.model.geo.addPoint(side_box, side_box, 0, lc)
p4 = gmsh.model.geo.addPoint(0, side_box, 0, lc)

cl1 = [l1, l2, l3, l4]
gmsh.model.geo.synchronize()

gmsh.model.mesh.generate(2)
final_mesh, cell_tags, facet_tags = gmshio.model_to_mesh(model, comm, model_rank)
gmsh.finalize()

tdim = final_mesh.topology.dim
fdim = tdim - 1

submesh, entity_map = msh.create_submesh(final_mesh, fdim, facet_tags.find(1))[0:2]

``````

# Rest of the code

``````deg    = 3
family = "Lagrange"

P = FunctionSpace(final_mesh, (family, deg))
Q = FunctionSpace(submesh, (family, deg))

p, q = TrialFunction(P), TrialFunction(Q)
v, u = TestFunction(P), TestFunction(Q)

dx  = Measure("dx", domain=final_mesh, subdomain_data=cell_tags)
ds  = Measure("ds", domain=final_mesh, subdomain_data=facet_tags)

e0x  = inner(ddqx, u)*ds(1)
a    = form(e0x)
``````

Here is the error (Very long)

# Error

``````---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[4], line 3
2 e0x  = inner(ddqx, u)*ds(1)
----> 3 a    = form(e0x)

File /usr/local/dolfinx-real/lib/python/dist-packages/dolfinx/fem/forms.py:194, in form(form, dtype, form_compiler_options, jit_options, entity_maps)
191         return list(map(lambda sub_form: _create_form(sub_form), form))
192     return form
--> 194 return _create_form(form)

File /usr/local/dolfinx-real/lib/python/dist-packages/dolfinx/fem/forms.py:189, in form.<locals>._create_form(form)
186 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
187 return form argument"""
188 if isinstance(form, ufl.Form):
--> 189     return _form(form)
190 elif isinstance(form, collections.abc.Iterable):
191     return list(map(lambda sub_form: _create_form(sub_form), form))

File /usr/local/dolfinx-real/lib/python/dist-packages/dolfinx/fem/forms.py:152, in form.<locals>._form(form)
150 if mesh is None:
151     raise RuntimeError("Expecting to find a Mesh in the form.")
--> 152 ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
153                                        form_compiler_options=form_compiler_options,
154                                        jit_options=jit_options)
156 # For each argument in form extract its function space
157 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]

File /usr/local/dolfinx-real/lib/python/dist-packages/dolfinx/jit.py:56, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
51 @functools.wraps(local_jit)
52 def mpi_jit(comm, *args, **kwargs):
53
54     # Just call JIT compiler when running in serial
55     if comm.size == 1:
---> 56         return local_jit(*args, **kwargs)
58     # Default status (0 == ok, 1 == fail)
59     status = 0

File /usr/local/dolfinx-real/lib/python/dist-packages/dolfinx/jit.py:204, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
202 # Switch on type and compile, returning cffi object
203 if isinstance(ufl_object, ufl.Form):
--> 204     r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
205 elif isinstance(ufl_object, ufl.FiniteElementBase):
206     r = ffcx.codegeneration.jit.compile_elements([ufl_object], options=p_ffcx, **p_jit)

File /usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py:199, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
197     except Exception:
198         pass
--> 199     raise e
201 obj, module = _load_objects(cache_dir, module_name, form_names)
202 return obj, module, (decl, impl)

File /usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py:190, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
187     for name in form_names:
188         decl += form_template.format(name=name)
--> 190     impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
191                             cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
192 except Exception as e:
193     try:
194         # remove c file so that it will not timeout next time

File /usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py:260, in _compile_objects(decl, ufl_objects, object_names, module_name, options, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
256 import ffcx.compiler
258 # JIT uses module_name as prefix, which is needed to make names of all struct/function
259 # unique across modules
--> 260 _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
262 ffibuilder = cffi.FFI()
263 ffibuilder.set_source(module_name, code_body, include_dirs=[ffcx.codegeneration.get_include_path()],
264                       extra_compile_args=cffi_extra_compile_args, libraries=cffi_libraries)

File /usr/local/lib/python3.10/dist-packages/ffcx/compiler.py:102, in compile_ufl_objects(ufl_objects, object_names, prefix, options, visualise)
100 # Stage 2: intermediate representation
101 cpu_time = time()
--> 102 ir = compute_ir(analysis, object_names, prefix, options, visualise)
103 _print_timing(2, time() - cpu_time)
105 # Stage 3: code generation

File /usr/local/lib/python3.10/dist-packages/ffcx/ir/representation.py:198, in compute_ir(analysis, object_names, prefix, options, visualise)
192 ir_elements = [_compute_element_ir(e, analysis.element_numbers, finite_element_names)
193                for e in analysis.unique_elements]
195 ir_dofmaps = [_compute_dofmap_ir(e, analysis.element_numbers, dofmap_names)
196               for e in analysis.unique_elements]
--> 198 irs = [_compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
199                             options, visualise)
200        for (i, fd) in enumerate(analysis.form_data)]
201 ir_integrals = list(itertools.chain(*irs))
203 ir_forms = [_compute_form_ir(fd, i, prefix, form_names, integral_names, analysis.element_numbers,
204                              finite_element_names, dofmap_names, object_names)
205             for (i, fd) in enumerate(analysis.form_data)]

File /usr/local/lib/python3.10/dist-packages/ffcx/ir/representation.py:198, in <listcomp>(.0)
192 ir_elements = [_compute_element_ir(e, analysis.element_numbers, finite_element_names)
193                for e in analysis.unique_elements]
195 ir_dofmaps = [_compute_dofmap_ir(e, analysis.element_numbers, dofmap_names)
196               for e in analysis.unique_elements]
--> 198 irs = [_compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
199                             options, visualise)
200        for (i, fd) in enumerate(analysis.form_data)]
201 ir_integrals = list(itertools.chain(*irs))
203 ir_forms = [_compute_form_ir(fd, i, prefix, form_names, integral_names, analysis.element_numbers,
204                              finite_element_names, dofmap_names, object_names)
205             for (i, fd) in enumerate(analysis.form_data)]

File /usr/local/lib/python3.10/dist-packages/ffcx/ir/representation.py:493, in _compute_integral_ir(form_data, form_index, element_numbers, integral_names, finite_element_names, options, visualise)
490 integrands = {rule: integral.integrand() for rule, integral in sorted_integrals.items()}
492 # Build more specific intermediate representation
--> 493 integral_ir = compute_integral_ir(itg_data.domain.ufl_cell(), itg_data.integral_type,
494                                   ir["entitytype"], integrands, ir["tensor_shape"],
495                                   ir["mixed_dim"], options, visualise)
497 ir.update(integral_ir)
499 # Fetch name

File /usr/local/lib/python3.10/dist-packages/ffcx/ir/integral.py:84, in compute_integral_ir(cell, integral_type, entitytype, integrands, argument_shape, is_mixed_dim, p, visualise)
74 # Build terminal_data from V here before factorization. Then we
75 # can use it to derive table properties for all modified
76 # terminals, and then use that to rebuild the scalar graph more
77 # efficiently before argument factorization. We can build
78 # terminal_data again after factorization if that's necessary.
80 initial_terminals = {i: analyse_modified_terminal(v['expression'])
81                      for i, v in S.nodes.items()
82                      if is_modified_terminal(v['expression'])}
---> 84 mt_table_reference = build_optimized_tables(
86     cell,
87     integral_type,
88     entitytype,
89     initial_terminals.values(),
90     ir["unique_tables"],
91     is_mixed_dim,
92     rtol=p["table_rtol"],
93     atol=p["table_atol"])
95 # Fetch unique tables for this quadrature rule
96 table_types = {v.name: v.ttype for v in mt_table_reference.values()}

File /usr/local/lib/python3.10/dist-packages/ffcx/ir/elementtables.py:351, in build_optimized_tables(quadrature_rule, cell, integral_type, entitytype, modified_terminals, existing_tables, is_mixed_dim, rtol, atol)
349 new_table = []
350 for ref in range(2):
--> 351     new_table.append(get_ffcx_table_values(
353         integral_type, element, avg, entitytype, local_derivatives, flat_component, is_mixed_dim))
355 t = new_table[0]
356 t['array'] = np.vstack([td['array'] for td in new_table])

File /usr/local/lib/python3.10/dist-packages/ffcx/ir/elementtables.py:144, in get_ffcx_table_values(points, cell, integral_type, element, avg, entitytype, derivative_counts, flat_component, is_mixed_dim)
141         raise RuntimeError(f"Domain cell and ufl element cell not compatible for {integral_type} integral")
143     tbl = component_element.tabulate(deriv_order, entity_points)
--> 144     tbl = tbl[basix_index(derivative_counts)]
145     component_tables.append(tbl)
147 if avg in ("cell", "facet"):
148     # Compute numeric integral of the each component table

IndexError: index 4 is out of bounds for axis 0 with size 3
``````

Within my project, in 3D with a 2D subdomain, the last line of the error is

``````IndexError: index 7 is out of bounds for axis 0 with size 6
``````

I’ve been very deep in FEniCSx files in order to understand what’s going on and how to solve it. But I have no clue…
Once again, the second order derivative of p (And all derivatives of p) doesn’t provide any error, the first derivative of q doesn’t neither.

Thank’s a billion for your help,
Regards,
Pierre.

Hi Folks,
Did anyone had chance to get into this?

Regards,

Dear Pierre,

As far as I can tell, you are looking for the support sketched out in: Add support for mixed-domain problems of codimension 1 by jpdean · Pull Request #3224 · FEniCS/dolfinx · GitHub
by @jpdean .
I.e. a mixed assembly of spaces defined on submeshes of different co-dimensions.