Replace IntegralType enum with something that is more expressive

[jorgensd]

Replaces dolfinx::fem::IntegralType-enum with a more expressive integral type.
The only thing IntegralType requires to know of the DOLFINx side is:

• What codimension we want to pack integration entities over
• How many integration entities/coefficients we need to pack per kernel

Default behavior is num_cells=1. Thus:

• fem::IntegralType::cell -> fem::IntegralType(0)
• fem::IntegralType::exterior_facet -> fem::IntegralType(1)
• fem::IntegralType::interior_facet -> fem::IntegralType(1, 2)
• fem::IntegralType::ridge -> fem::IntegralType(2)
• fem::IntegralType::vertex -> fem::IntegralType(-1)

Future integral type that should replace vertex is peak, and it should be fem::IntegralType(3).

Depends on FEniCS/ffcx#787 (due to renaming of exterior_facet to facet to be consistent with
vertex and ridge.

As I have kept IntegralType as input to all existing functions, they only user-facing API change is intializing the IntegralType.

Also fixes assigning environment variables through config file on windows CI.

This idea will also extend to Integrals where one have more than one or two cells that should be packed for a kernel.

[jorgensd]

@garth-wells CI (via Spack, experimental) is blocking here.

[jorgensd]

@garth-wells I've now replaced the dolfinx::fem::IntegralType with a struct that currently holds codim and num_cells, which expresses all integral types we currently have (dx, ds, dS, dP), future ones as dr (ridge) and even integral types we haven't looked at yet, such as interior peak or ridge integrals, as well as interior_facet on branched manifolds.

[jorgensd]

Thinking even more about this, I think the integral type should also contain a default_domain enum, which should be all or exterior to distinguish between dP and ds when tdim is 1.

currently these two measure would give different behavior for the two measures, which is unintuitive, at least to me. This would affect compute_integration_domains only. @schnellerhase ; a penny for your thoughts:)

thinking even more about it, ds,dP and dr (facet, vertex, ridge) integrals should only take entities on the exterior boundary as default. Anything else will lead to non-deterministic results when having discontinuous data in such cells, and calling dP or dr, as the first connected cell is not the same when executed with a different number of processes.

if a user wants to integrate over all vertices, they should compute the integration data themselves to ensure correctness, and attach it as subdomain data.

[jorgensd]

@michalhabera @schnellerhase
I think we need to redesign the vertex measure right away, as it has introduced some very weird behavior on 1D grids:

• dP on a 1D grid takes all owned facets, restricted to one side of a cell (arbitrary)
• ds on a 1D grid takes all owned exterior facets (consistent)
• dS on a 1D grid takes all owned facets that are not exterior (consistent).

I think dP shouldn't work on a 1D grid.

[schnellerhase]

@michalhabera @schnellerhase I think we need to redesign the vertex measure right away, as it has introduced some very weird behavior on 1D grids:

• dP on a 1D grid takes all owned facets, restricted to one side of a cell (arbitrary)
• ds on a 1D grid takes all owned exterior facets (consistent)
• dS on a 1D grid takes all owned facets that are not exterior (consistent).

I think dP shouldn't work on a 1D grid.

I agree these are a lot of special cases, but they are not combinable at the moment. Consider a branching mesh (truss layout), here dS will not work as expected on the joints, and the classification of boundary facets for ds is no longer correct, so the assembly with measures ds + dS will not be equivalent to dP. It covers a different use case.

thinking even more about it, ds,dP and dr (facet, vertex, ridge) integrals should only take entities on the exterior boundary as default. Anything else will lead to non-deterministic results when having discontinuous data in such cells, and calling dP or dr, as the first connected cell is not the same when executed with a different number of processes.

Integration of discontinuous data should not be performed with these integral types. An attempt of facilitating this is to be found at https://github.com/FEniCS/ffcx/blob/main/ffcx/analysis.py#L205-L209 for dP, which however does not fully protect against it - is there functionality in place to check for integrand continuity at the moment?

[schnellerhase]

I don't think IntegralType carrying the num_cells information is the right abstraction. This is connectivity information from the mesh, and changes across the domain. As we have not figured out how to do patch assemblies for now, best to not add interface support for it.

[garth-wells]

I don't think IntegralType carrying the num_cells information is the right abstraction. This is connectivity information from the mesh, and changes across the domain. As we have not figured out how to do patch assemblies for now, best to not add interface support for it.

Agree - just need the dimension and special case handling for interior facets for now.

[jorgensd]

I don't think IntegralType carrying the num_cells information is the right abstraction. This is connectivity information from the mesh, and changes across the domain. As we have not figured out how to do patch assemblies for now, best to not add interface support for it.

Using num_cells or a Boolean flag doesn’t change the implementation much.

This depends on how you want to generate your kernels (in the future). If you think along the lines of how we do interior facet, the kernel needs to know how many cells (for packing coordinate data, coefficients etc) it is taking in. For instance, if I could specify restriction (0,1,2,….,N) in UFL, rather than +,-, this is the information that would be required.

One of my key points is that dP (vertex) does all the same assumptions as my ridge implementation (things should be one sided), which is the same when manually packing ds integral domains, which one can use to integrate over any facet of your mesh.