Hello everyone,
I am currently experiencing difficulties with form creation in an HDG scheme for compressible Euler equations when the mesh is one-dimensional with trace elements defined at the extremities. More specifically, there appears to be an issue during the block form creation. This error also manifests when I run the Poisson equation tutorial and replace the 3D mesh with a 1D mesh. Would anyone have advice to tell me what I’m missing?
Here is the code:
from mpi4py import MPI
import numpy as np
import ufl
from petsc4py import PETSc
import dolfinx
from dolfinx import fem, mesh
from dolfinx.cpp.mesh import cell_num_entities
def compute_cell_boundary_facets(msh):
"""Compute the integration entities for integrals around the
boundaries of all cells in msh.
Parameters:
msh: The mesh.
Returns:
Facets to integrate over, identified by ``(cell, local facet
index)`` pairs.
"""
tdim = msh.topology.dim
fdim = tdim - 1
n_f = cell_num_entities(msh.topology.cell_type, fdim)
n_c = msh.topology.index_map(tdim).size_local
return np.vstack((np.repeat(np.arange(n_c), n_f), np.tile(np.arange(n_f), n_c))).T.flatten()
# Number of elements
n = 8
# Create the mesh
comm = MPI.COMM_WORLD
rank = comm.rank
dtype = PETSc.ScalarType
#Not working
msh = mesh.create_unit_interval(comm, n)
#Working
# msh = mesh.create_unit_square(comm, n, n)
# msh = mesh.create_unit_cube(comm, n, n, n)
# We need to create a broken Lagrange space defined over the facets of the
# mesh. To do so, we require a sub-mesh of the all facets. We begin by
# creating a list of all of the facets in the mesh
tdim = msh.topology.dim
fdim = tdim - 1
msh.topology.create_entities(fdim)
facet_imap = msh.topology.index_map(fdim)
num_facets = facet_imap.size_local + facet_imap.num_ghosts
facets = np.arange(num_facets, dtype=np.int32)
# Create the sub-mesh
# NOTE Despite all facets being present in the submesh, the entity map
# isn't necessarily the identity in parallel
facet_mesh, facet_mesh_to_mesh = mesh.create_submesh(msh, fdim, facets)[:2]
# Define function spaces
k = 0 # Polynomial order
V = fem.functionspace(msh, ("Discontinuous Lagrange", k))
Vbar = fem.functionspace(facet_mesh, ("Discontinuous Lagrange", k))
# Trial and test functions in mixed space
W = ufl.MixedFunctionSpace(V, Vbar)
u, ubar = ufl.TrialFunctions(W)
v, vbar = ufl.TestFunctions(W)
# Define integration measures
# Cell
dx_c = ufl.Measure("dx", domain=msh)
# Cell boundaries
# We need to define an integration measure to integrate around the
# boundary of each cell. The integration entities can be computed
# using the following convenience function.
cell_boundary_facets = compute_cell_boundary_facets(msh)
cell_boundaries = 1 # A tag
# Create the measure
ds_c = ufl.Measure("ds", subdomain_data=[(cell_boundaries, cell_boundary_facets)], domain=msh)
# Create a cell integral measure over the facet mesh
dx_f = ufl.Measure("dx", domain=facet_mesh)
# We write the mixed domain forms as integrals over msh. Hence, we must
# provide a map from facets in msh to cells in facet_mesh. This is the
# 'inverse' of facet_mesh_to_mesh, which we compute as follows:
mesh_to_facet_mesh = np.full(num_facets, -1)
mesh_to_facet_mesh[facet_mesh_to_mesh] = np.arange(len(facet_mesh_to_mesh))
entity_maps = {facet_mesh: mesh_to_facet_mesh}
# Define forms
h = ufl.CellDiameter(msh)
n = ufl.FacetNormal(msh)
gamma = 16.0 * k**2 / h # Scaled penalty parameter
x = ufl.SpatialCoordinate(msh)
c = 1.0 + 0.1 * ufl.sin(ufl.pi * x[0])
a = (
ufl.inner(c * ufl.grad(u), ufl.grad(v)) * dx_c
- ufl.inner(c * (u - ubar), ufl.dot(ufl.grad(v), n)) * ds_c(cell_boundaries)
- ufl.inner(ufl.dot(ufl.grad(u), n), c * (v - vbar)) * ds_c(cell_boundaries)
+ gamma * ufl.inner(c * (u - ubar), v - vbar) * ds_c(cell_boundaries)
)
# Define block structure
a_blocked = ufl.extract_blocks(a)
a_blocked_form = dolfinx.fem.form(a_blocked, entity_maps=entity_maps)
Which raises the error:
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
File /nfs/home/bouteillerp/Téléchargements/demo_hdg(1).py:107
105 # Define block structure
106 a_blocked = ufl.extract_blocks(a)
--> 107 a_blocked_form = dolfinx.fem.form(a_blocked, entity_maps=entity_maps)
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/fem/forms.py:337, in form(form, dtype, form_compiler_options, jit_options, entity_maps)
334 else:
335 return form
--> 337 return _create_form(form)
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/fem/forms.py:333, in form.<locals>._create_form(form)
331 return _form(form)
332 elif isinstance(form, collections.abc.Iterable):
--> 333 return list(map(lambda sub_form: _create_form(sub_form), form))
334 else:
335 return form
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/fem/forms.py:333, in form.<locals>._create_form.<locals>.<lambda>(sub_form)
331 return _form(form)
332 elif isinstance(form, collections.abc.Iterable):
--> 333 return list(map(lambda sub_form: _create_form(sub_form), form))
334 else:
335 return form
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/fem/forms.py:333, in form.<locals>._create_form(form)
331 return _form(form)
332 elif isinstance(form, collections.abc.Iterable):
--> 333 return list(map(lambda sub_form: _create_form(sub_form), form))
334 else:
335 return form
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/fem/forms.py:333, in form.<locals>._create_form.<locals>.<lambda>(sub_form)
331 return _form(form)
332 elif isinstance(form, collections.abc.Iterable):
--> 333 return list(map(lambda sub_form: _create_form(sub_form), form))
334 else:
335 return form
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/fem/forms.py:331, in form.<locals>._create_form(form)
329 return None
330 else:
--> 331 return _form(form)
332 elif isinstance(form, collections.abc.Iterable):
333 return list(map(lambda sub_form: _create_form(sub_form), form))
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/fem/forms.py:254, in form.<locals>._form(form)
252 if mesh is None:
253 raise RuntimeError("Expecting to find a Mesh in the form.")
--> 254 ufcx_form, module, code = jit.ffcx_jit(
255 mesh.comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options
256 )
258 # For each argument in form extract its function space
259 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/jit.py:62, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
58 @functools.wraps(local_jit)
59 def mpi_jit(comm, *args, **kwargs):
60 # Just call JIT compiler when running in serial
61 if comm.size == 1:
---> 62 return local_jit(*args, **kwargs)
64 # Default status (0 == ok, 1 == fail)
65 status = 0
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/jit.py:212, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
210 # Switch on type and compile, returning cffi object
211 if isinstance(ufl_object, ufl.Form):
--> 212 r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
213 elif isinstance(ufl_object, ufl.Mesh):
214 r = ffcx.codegeneration.jit.compile_coordinate_maps([ufl_object], options=p_ffcx, **p_jit)
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/codegeneration/jit.py:225, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
223 except Exception:
224 pass
--> 225 raise e
227 obj, module = _load_objects(cache_dir, module_name, form_names)
228 return obj, module, (decl, impl)
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/codegeneration/jit.py:205, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
202 for name in form_names:
203 decl += form_template.format(name=name)
--> 205 impl = _compile_objects(
206 decl,
207 forms,
208 form_names,
209 module_name,
210 p,
211 cache_dir,
212 cffi_extra_compile_args,
213 cffi_verbose,
214 cffi_debug,
215 cffi_libraries,
216 visualise=visualise,
217 )
218 except Exception as e:
219 try:
220 # remove c file so that it will not timeout next time
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/codegeneration/jit.py:330, in _compile_objects(decl, ufl_objects, object_names, module_name, options, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
326 libraries = _libraries + cffi_libraries if cffi_libraries is not None else _libraries
328 # JIT uses module_name as prefix, which is needed to make names of all struct/function
329 # unique across modules
--> 330 _, code_body = ffcx.compiler.compile_ufl_objects(
331 ufl_objects, prefix=module_name, options=options, visualise=visualise
332 )
334 # Raise error immediately prior to compilation if no support for C99
335 # _Complex. Doing this here allows FFCx to be used for complex codegen on
336 # Windows.
337 if sys.platform.startswith("win32"):
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/compiler.py:113, in compile_ufl_objects(ufl_objects, options, object_names, prefix, visualise)
111 # Stage 2: intermediate representation
112 cpu_time = time()
--> 113 ir = compute_ir(analysis, _object_names, _prefix, options, visualise)
114 _print_timing(2, time() - cpu_time)
116 # Stage 3: code generation
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/ir/representation.py:135, in compute_ir(analysis, object_names, prefix, options, visualise)
129 for itg_index, itg_data in enumerate(fd.integral_data):
130 integral_names[(fd_index, itg_index)] = naming.integral_name(
131 fd.original_form, itg_data.integral_type, fd_index, itg_data.subdomain_id, prefix
132 )
134 irs = [
--> 135 _compute_integral_ir(
136 fd,
137 i,
138 analysis.element_numbers,
139 integral_names,
140 finite_element_hashes,
141 options,
142 visualise,
143 )
144 for (i, fd) in enumerate(analysis.form_data)
145 ]
146 ir_integrals = list(itertools.chain(*irs))
148 ir_forms = [
149 _compute_form_ir(
150 fd,
(...)
159 for (i, fd) in enumerate(analysis.form_data)
160 ]
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/ir/representation.py:396, in _compute_integral_ir(form_data, form_index, element_numbers, integral_names, finite_element_hashes, options, visualise)
391 integrand_map: dict[QuadratureRule, ufl.core.expr.Expr] = {
392 rule: integral.integrand() for rule, integral in sorted_integrals.items()
393 }
395 # Build more specific intermediate representation
--> 396 integral_ir = compute_integral_ir(
397 itg_data.domain.ufl_cell(),
398 itg_data.integral_type,
399 expression_ir["entity_type"],
400 integrand_map,
401 expression_ir["tensor_shape"],
402 options,
403 visualise,
404 )
406 expression_ir.update(integral_ir)
408 # Fetch name
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/ir/integral.py:91, in compute_integral_ir(cell, integral_type, entity_type, integrands, argument_shape, p, visualise)
88 if domain.topological_dimension() != cell.topological_dimension():
89 is_mixed_dim = True
---> 91 mt_table_reference = build_optimized_tables(
92 quadrature_rule,
93 cell,
94 integral_type,
95 entity_type,
96 initial_terminals.values(),
97 ir["unique_tables"],
98 use_sum_factorization=p["sum_factorization"],
99 is_mixed_dim=is_mixed_dim,
100 rtol=p["table_rtol"],
101 atol=p["table_atol"],
102 )
104 # Fetch unique tables for this quadrature rule
105 table_types = {v.name: v.ttype for v in mt_table_reference.values()}
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/ir/elementtables.py:444, in build_optimized_tables(quadrature_rule, cell, integral_type, entity_type, modified_terminals, existing_tables, use_sum_factorization, is_mixed_dim, rtol, atol)
442 t["array"] = np.vstack([td["array"] for td in new_table])
443 else:
--> 444 t = get_ffcx_table_values(
445 quadrature_rule.points,
446 cell,
447 integral_type,
448 element,
449 avg,
450 entity_type,
451 local_derivatives,
452 flat_component,
453 codim,
454 )
455 # Clean up table
456 tbl = clamp_table_small_numbers(t["array"], rtol=rtol, atol=atol)
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/ir/elementtables.py:153, in get_ffcx_table_values(points, cell, integral_type, element, avg, entity_type, derivative_counts, flat_component, codim)
151 raise RuntimeError("Codimension > 1 isn't supported.")
152 tbl = component_element.tabulate(deriv_order, entity_points)
--> 153 tbl = tbl[basix_index(derivative_counts)]
154 component_tables.append(tbl)
156 if avg in ("cell", "facet"):
157 # Compute numeric integral of the each component table
File /nfs/home/bouteillerp/miniforge3/envs/fenicsx-env/lib/python3.13/site-packages/ffcx/element_interface.py:17, in basix_index(indices)
15 def basix_index(indices: tuple[int]) -> int:
16 """Get the Basix index of a derivative."""
---> 17 return basix.index(*indices)
TypeError: index() missing 1 required positional argument: 'p'
Thanks for your advice, sincerely,
Paul