Unification of one-sided sub-entity integrals + ridge integrals

[jorgensd]

New features

• Introduces ridge integrals (ufl.dr) (tdim - 2) (for rank 0, 1 and 2-tensors)
• Unifies assembly routines in assemble_vector for all exterior entity integrals.
• Adds assemble_matrix and apply_lifting for vertex integrals

Context

Any one-sided sub entity integral boils down to getting the correct (cell, local_entity_index), and assemble over these.
The lower level assembly routines (assemble_exterior_facet, assemble_vertex) does not really differ in this respect. It is only due to the fact that facets might require permutations that codes were not identical.
This became crystal clear when introducing ridge integrals.

Consequences

• The assemble_vertex routines recently introduced are removed in favor of an assemble_exterior_facet function, now renamed as assemble_exterior_entities. This is also used for ridge integrals.

Bugs encountered

• There is a subtle bug in the coefficient allocation of vertex integrals, where the double amount of memory is allocated. Adjusted packing to divide by two for all integral types but cells.

Future work

Mixed dimensional support for codim-2.
For this to work, a few things have to be introduced:

• require_ridge_permutation has to be introduced in FFCx
• Computing ridge-permutation has to be a member of DOLFINx (similar to compute facet permutations)

[schnellerhase]

Nice code duplication removal!

[michalhabera]

The unification reads really well.

One point to consider, for future design, given the increase in type of integral we are seeing:
Vertex integrals are effectively defined by dimension (vertex: dim=0), while all other come with a codimension (cell: codim=0, facets: codim=-1, ridge: codim=-2). The root of this is also in UFL, and surely we want to keep it for user convenience, but we should try to separate the logic early enough, strip all data to a common pattern for both dim/codim and follow unified codepath onwards. I think your PR is the step in right direction as it moves all integral assemble routines to the codimension format. However, we should also avoid IntegralType with duplicate values as much as possible: in 2D ridge integral in == vertex integral.

[jorgensd]

The unification reads really well.

One point to consider, for future design, given the increase in type of integral we are seeing: Vertex integrals are effectively defined by dimension (vertex: dim=0), while all other come with a codimension (cell: codim=0, facets: codim=-1, ridge: codim=-2). The root of this is also in UFL, and surely we want to keep it for user convenience, but we should try to separate the logic early enough, strip all data to a common pattern for both dim/codim and follow unified codepath onwards. I think your PR is the step in right direction as it moves all integral assemble routines to the codimension format. However, we should also avoid IntegralType with duplicate values as much as possible: in 2D ridge integral in == vertex integral.

I think if we ever introduce "peak" integrals, they would effectively replace the vertex integrals, as we should always go down the "co-dimension" route (for the sake of mixed assembly support).

[schnellerhase]

Once we have the ridge integral which aligns with vertex integral for tdim=2, then for tdim=3 peak=vertex so we can just rename vertex->peak after introducing ridge.

[garth-wells]

@garth-wells Renaming exterior_facet to facet cascades deep into FFCx, and should be handled in a separate PR, as there are tons of bad use-cases of ufl strings in there. Especially if we now want to remove the enum altogether, I think we leave it in this PR and change in soon.

Something has gone very badly wrong in the DOLFINX design if this is the case. It should be relatively easy to rename an enum, or change the enum to a pair of integers, in DOLFINx without changing UFCx or FFCx. What specifically causes a 'cascade' into FFCx?

[jorgensd]

@garth-wells Renaming exterior_facet to facet cascades deep into FFCx, and should be handled in a separate PR, as there are tons of bad use-cases of ufl strings in there. Especially if we now want to remove the enum altogether, I think we leave it in this PR and change in soon.

Something has gone very badly wrong in the DOLFINX design if this is the case. It should be relatively easy to rename an enum, or change the enum to a pair of integers, in DOLFINx without changing UFCx or FFCx. What specifically causes a 'cascade' into FFCx?

This is being handled separately in:
FEniCS/ffcx#787

[garth-wells]

@garth-wells Renaming exterior_facet to facet cascades deep into FFCx, and should be handled in a separate PR, as there are tons of bad use-cases of ufl strings in there. Especially if we now want to remove the enum altogether, I think we leave it in this PR and change in soon.

Something has gone very badly wrong in the DOLFINX design if this is the case. It should be relatively easy to rename an enum, or change the enum to a pair of integers, in DOLFINx without changing UFCx or FFCx. What specifically causes a 'cascade' into FFCx?

This is being handled separately in: FEniCS/ffcx#787

But what's the issue? I can't see it. If we have such tight coupling between DOLFINx and FFCx/UFCx there is a design mistake that needs to be fixed.

[jorgensd]

@garth-wells Renaming exterior_facet to facet cascades deep into FFCx, and should be handled in a separate PR, as there are tons of bad use-cases of ufl strings in there. Especially if we now want to remove the enum altogether, I think we leave it in this PR and change in soon.

Something has gone very badly wrong in the DOLFINX design if this is the case. It should be relatively easy to rename an enum, or change the enum to a pair of integers, in DOLFINx without changing UFCx or FFCx. What specifically causes a 'cascade' into FFCx?

This is being handled separately in: FEniCS/ffcx#787

But what's the issue? I can't see it. If we have such tight coupling between DOLFINx and FFCx/UFCx there is a design mistake that needs to be fixed.

Currently in DOLFINx we look things up by strings, and assume that they match with things in FFCx and UFL.
This is bad design and something I have taken care of in:

• Replace IntegralType enum with something that is more expressive #3919
• Introduce integral type and try to separate ufl concepts from FFCx concepts ffcx#787

[michalhabera]

What are you checking here? If the IntegralType happens to be codim=1, i.e. if the entities we integrate over happen to be facets?

[jorgensd]

This was initially an if statement, then a switch then moved back after Garth’s comment.

in general any integral of codim>0 needs permutations (except codim=tdim). At the time of writing, we only have facet permutations implemented.
It matters for mesh-submesh coupling, ie this PR does not introduce 3D-1D mixed problems, as I need ridge permutations, which is in another PR (#3904)

[michalhabera]

Should the condition be what you said then, that codim > 0 (or better codim=1, and warning for codim > 1)? I don't mind either of the if, switch or ternary, I am asking about the logic.

[jorgensd]

With what we currently have implemented, it should be equality.

if I get ridge permutations into DOLFINx, it will be an if codim>0. Garth asked to not implement for features that are not there yet, which is why I removed the switch and went for a tenary logic for what we support.

[michalhabera]

I am happy with this PR, I only see minor potential for naming improvement, as after the PR we will have:

1. assemble_entities
   assembles over provided entities (regardless of their dim/codim) restricted to 0th adjacent cell, could as well be called assemble_entities_cell0?.

And two more specialized cases:

2. assemble_cells
   assembles over cells restricted to 0th (but the only) adjacent cell (itself),
3. assemble_interior_facets
   assembles over facets (codim=1) with restriction being based on exactly first two adjacent cells.

If we'd like to unify the special cases, with some work we can use assemble_entities_cell0 to handle assemble_cells -- the only complication now I think is different data format of entities array, and also generalize assemble_interior_facets to assemble_entities_cell0cell1. This proposed naming is very explicit in what the function does.

This PR does not add any complexity, it actually exposes the shortcomings of the current design where the assumptions on: 1. what do we integrate over and 2. what adjacent entities is it restricted to, are simply buried in the code, not clear from the API.

[schnellerhase]

Suggested change

	/// @brief Execute kernel over entities of codimension > 1 and accumulate result
	/// @brief Execute kernel over entities of codimension ≥ 1 and accumulate result

[schnellerhase]

Suggested change

	/// @brief Execute kernel over entities of codimension > 1 and accumulate result
	/// @brief Execute kernel over entities of codimension ≥ 1 and accumulate result

[schnellerhase]

Suggested change

	/// @brief Execute kernel over entities of codimension > 1 and accumulate result
	/// @brief Execute kernel over entities of codimension ≥ 1 and accumulate result

[schnellerhase]

Suggested change

	/// @param[in] perms Facet permutation integer. Empty if facet
	/// @param[in] perms Entity permutation integer. Empty if entity

[cdaversin]

Glad to see the ridge integrals are coming their way! The unification makes real sense to me. In addition to the significant duplicated code removal, the "re-naming" of the assembly functions with entities instead of facets makes the code more generic and less confusing when working with various dimension / co-dimensions.

I agree with @michalhabera 's point on moving assembly routines to "the codimension format". I think @jorgensd 's PR#3919 with a dedicated class IntegralType is a smart way to define the integral types from their codimension and number of adjacent cells. This class is simple/light but flexible to include future type of integrals (some which can be difficult to name / define as a string)

This is a big step towards codimension > 1 :)