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()
gmsh.model.add("test")
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)
l1 = gmsh.model.geo.addLine(p1, p2)
l2 = gmsh.model.geo.addLine(p2, p3)
l3 = gmsh.model.geo.addLine(p3, p4)
l4 = gmsh.model.geo.addLine(p4, p1)
cl1 = [l1, l2, l3, l4]
s1 = gmsh.model.geo.addPlaneSurface([gmsh.model.geo.addCurveLoop(cl1)])
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [l1], tag=1)
gmsh.model.addPhysicalGroup(2, [s1], tag=1)
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)
ddqx = grad(grad(q)[0])[0]
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
1 ddqx = grad(grad(q)[0])[0]
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(
85 quadrature_rule,
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(
352 permute_quadrature_interval(quadrature_rule.points, ref), cell,
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.