Simplify integration domain packing

[jorgensd]

Resolve #3214 and #3213.

This PR changes the packing such that we only pack integral entities for those subdomain id's that exist in the integrals.
For processors with no entities for a given tag, we return the tuple (id, []) indicating that there is an integral of this type on the process, but that it has no entities.

This simplifies the packing of integration entities massively, as we only pack those specified by the form data, not for every value in the meshtags. It also returns a cleaner data-structure that can be used directly by users.

[jorgensd]

@jpdean Thanks for your comments. I've addressed them by simplifying compute_integration_entities.

[jorgensd]

Is there anything specific holding this PR back?
it simplifies code and resolves multiple issues.

[garth-wells]

Is there anything specific holding this PR back? it simplifies code and resolves multiple issues.

Wait for @jpdean to return and add approval.

[jpdean]

I think letting compute_integration_domains just deal with turning a list of mesh entities into data needed for the assembler (instead of worrying about meshtag values) is a good idea. It fits the DOLFINx design philosophy better than the function in main. One downside is that it basically does nothing for cells now, but that's not a big issue. It might make user C++ code a bit more complex though, because they'll have to loop over integration domains and call this function manually. Overall, I prefer the approach in this PR to the one in main.

[jorgensd]

It might make user C++ code a bit more complex though, because they'll have to loop over integration domains and call this function manually. Overall, I prefer the approach in this PR to the one in main.

As I've seen no C++ example where this function is actually used, I don't think we should worry to much.
The current code does way less work for problems where the meshtags are specified for all entities.
For instance, imagine that you have a facet marker where each facet has a unique number (local to process). The function on main will pack integration data for every facet on the process, while the new function will pack data only for those present in subdomain_ids.

[jpdean]

Looks good to me! One thing I'd like to do is simplify this logic if possible:

dolfinx/python/dolfinx/fem/forms.py

Lines 203 to 242 in 159cb1c

	# Make map from integral_type to subdomain id
	subdomain_ids = {type: [] for type in sd.get(domain).keys()}
	for integral in form.integrals():
	if integral.subdomain_data() is not None:
	subdomain_ids[integral.integral_type()].append(integral.subdomain_id())
	# Chain and sort subdomain ids
	for itg_type, marker_ids in subdomain_ids.items():
	flattened_ids = list(chain(marker_ids))
	flattened_ids.sort()
	subdomain_ids[itg_type] = flattened_ids
	def get_integration_domains(integral_type, subdomain, subdomain_ids):
	"""Get integration domains from subdomain data"""
	if subdomain is None:
	return []
	else:
	domains = []
	try:
	if integral_type in (IntegralType.exterior_facet, IntegralType.interior_facet):
	tdim = subdomain.topology.dim
	subdomain._cpp_object.topology.create_connectivity(tdim - 1, tdim)
	subdomain._cpp_object.topology.create_connectivity(tdim, tdim - 1)
	# Compute integration domains only for each subdomain id in the integrals
	# If a process has no integral entities, insert an empty array
	for id in subdomain_ids:
	try:
	integration_entities = _cpp.fem.compute_integration_domains(
	integral_type,
	subdomain._cpp_object.topology,
	subdomain.find(id),
	subdomain.dim,
	)
	domains.append((id, integration_entities))
	except TypeError:
	pass # If subdomain id is "everywhere"
	return [(s[0], np.array(s[1])) for s in domains]
	except AttributeError:
	return [(s[0], np.array(s[1])) for s in subdomain]

It was already complicated in main, and this PR adds more complexity. Happy to look at this in a separate PR though; maybe we should make an issue?

[jorgensd]

Looks good to me! One thing I'd like to do is simplify this logic if possible:

dolfinx/python/dolfinx/fem/forms.py

Lines 203 to 242 in 159cb1c

	# Make map from integral_type to subdomain id
	subdomain_ids = {type: [] for type in sd.get(domain).keys()}
	for integral in form.integrals():
	if integral.subdomain_data() is not None:
	subdomain_ids[integral.integral_type()].append(integral.subdomain_id())
	# Chain and sort subdomain ids
	for itg_type, marker_ids in subdomain_ids.items():
	flattened_ids = list(chain(marker_ids))
	flattened_ids.sort()
	subdomain_ids[itg_type] = flattened_ids
	def get_integration_domains(integral_type, subdomain, subdomain_ids):
	"""Get integration domains from subdomain data"""
	if subdomain is None:
	return []
	else:
	domains = []
	try:
	if integral_type in (IntegralType.exterior_facet, IntegralType.interior_facet):
	tdim = subdomain.topology.dim
	subdomain._cpp_object.topology.create_connectivity(tdim - 1, tdim)
	subdomain._cpp_object.topology.create_connectivity(tdim, tdim - 1)
	# Compute integration domains only for each subdomain id in the integrals
	# If a process has no integral entities, insert an empty array
	for id in subdomain_ids:
	try:
	integration_entities = _cpp.fem.compute_integration_domains(
	integral_type,
	subdomain._cpp_object.topology,
	subdomain.find(id),
	subdomain.dim,
	)
	domains.append((id, integration_entities))
	except TypeError:
	pass # If subdomain id is "everywhere"
	return [(s[0], np.array(s[1])) for s in domains]
	except AttributeError:
	return [(s[0], np.array(s[1])) for s in subdomain]

It was already complicated in main, and this PR adds more complexity. Happy to look at this in a separate PR though; maybe we should make an issue?

We would have to change how ufl exposes integrals to touch upon this (I think).
It should be a separate issue.