Refactor dolfinx::fem::Form to use local indexing of forms, rather than integral ids

[jorgensd]

As shown in: #3735, there are many UFL forms (especially in the case of mixed quadrature rules), where the subdomain_id is not mapped 1-1 with an integration kernel.

This PR refactors the Form class, to rather use the local indexing from the generated code (where integrals have been grouped) as the lookup key for the integral.

Only a minor change occurs in the user-interface, as Form::subdomain_ids no longer maps to the UFL subdomain_ids.
However, this is currently only used for internal looping, as the integral_idx -> subdomain_id mapping happens within the create_form_factory function.

This allows for different quadrature degrees in different parts of a form.

Example:

import basix.ufl
import ufl

cell = basix.CellType.triangle
coord_element = basix.ufl.element("Lagrange", cell, 1, shape=(2,))
domain = ufl.Mesh(coord_element)

degree = 2
el = basix.ufl.element("Lagrange", cell, degree)
V = ufl.FunctionSpace(domain, el)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

dx = ufl.Measure("dx", domain=domain)

a_mass = ufl.inner(u, v)*dx((1, 2), degree=2*degree)
a_stiffness = ufl.inner(ufl.grad(u), ufl.grad(v))*dx((1,), degree=2*(degree-1))

a = a_mass + a_stiffness

This generates two different integration kernels in FFCx, with the integral information as:

uint64_t finite_element_hashes_form_e26395feca75c800a7470ee9a844c8cdc67aff23[2] = {UINT64_C(16933917890882727401), UINT64_C(16933917890882727401)};
int form_integral_offsets_form_e26395feca75c800a7470ee9a844c8cdc67aff23[4] = {0, 3, 3, 3};
static ufcx_integral* form_integrals_form_e26395feca75c800a7470ee9a844c8cdc67aff23[3] = {&integral_0dcf3cd8d98806a1645751b277f55c83190bebb9_triangle, &integral_be233757c48f146683317dc29c6fe8fa75e720a8_triangle, &integral_be233757c48f146683317dc29c6fe8fa75e720a8_triangle};
int form_integral_ids_form_e26395feca75c800a7470ee9a844c8cdc67aff23[3] = {1, 1, 2};

As seen in the latter line here, form_integral_ids_form_* is not a list of unique integers.
In the main branch of DOLFINx, we assume that this id is unique, which is the root cause of #3735.

[garth-wells]

Can you add some gentle background to the PR description?

[jorgensd]

Can you add some gentle background to the PR description?

I've added an example to the description of the PR.

[jorgensd]

@garth-wells any further comments?

[garth-wells]

Yes, will add comments asap.

[jorgensd]

@garth-wells I've addressed your comments. Good to merge?

[michalhabera]

This is more of a critical bug-fix, so should go in asap. Btw. is the local index now needed? Looks like it will always be 0, 1, 2, ..., num_integrals

[jorgensd]

This is more of a critical bug-fix, so should go in asap. Btw. is the local index now needed? Looks like it will always be 0, 1, 2, ..., num_integrals

They will if we integrate over all kernels and associated subdomain ids.
The index is needed to associate a given kernel with a specific integration domain (see:

dolfinx/cpp/dolfinx/fem/assemble_matrix_impl.h

Lines 594 to 598 in d74fcba

	auto fn = a.kernel(IntegralType::cell, i, cell_type_idx);
	assert(fn);
	std::span cells = a.domain(IntegralType::cell, i, cell_type_idx);
	std::span cells0 = a.domain_arg(IntegralType::cell, 0, i, cell_type_idx);
	std::span cells1 = a.domain_arg(IntegralType::cell, 1, i, cell_type_idx);

).
If you do not have that index, you have no sensible way of associating the set of forms with some integral ideas.

I also think this should get merged, as it makes it possible to customize quadrature per integral measure used.

[michalhabera]

@garth-wells Anything more needs addressing? We should merge this bugfix and can improve the data layout further if needed.

All comments resolved. No new changes requested.

[michalhabera]

I've created issue #3876 to discuss a possible redesign of integrals storage. The bug is fixed with this PR, but situation could be much improved.