serialize oriented boundaries

[gdmcbain]

#865 added support for orientation in FacetBasis, but how do we serialize the ‘boundaries’ (as subsets of facets of a Mesh are known)?

• @kinnala 2022-01-27

For some reason I cannot figure out how to put those orientations into _encode_cell_data. :/

Edit: Maybe I will just try encoding the corresponding element indices instead.

…

I tried doing negative indices, but then I realized that the encoding is ultimately based on having matching indices in t2f and there are no negative indices there.

I was thinking that if instead of encoding 0 or 1 per facet ix, we somehow binary encode the triangle Mesh.f2t[0, ix] or Mesh.f2t[1, ix], and then in decode try reversing that process. Because the same triangle might belong to multiple oriented facets it is not completely obvious how it should work.

…

Try to do it and you'll notice where the issue is I think. :-)

[gdmcbain]

My first idea is not exactly to encode the orientation by signing the index of the facet since this is confusing for facet zero but instead to index the facets twice, once for the default orientation and then once flipped, like so:

scikit-fem/skfem/mesh/mesh.py

Lines 227 to 233 in d02e967

	def encode_boundary(boundary: Union[ndarray, OrientedBoundary]) -> ndarray:
	if hasattr(boundary, "ori"):
	ori = boundary.ori.astype(bool)
	cells = np.vstack([np.isin(self.t2f, boundary[o]) for o in [ori, ~ori]])
	else:
	cells = np.isin(self.t2f, boundary)
	return (1 << np.arange(cells.shape[0])) @ cells

[gdmcbain]

Decoding is straightforward: after unpacking the binary, the upper half holds the already oriented facets and the lower those that need flipping.

scikit-fem/skfem/mesh/mesh.py

Lines 260 to 284 in abb4578

	oriented = data[0].max() > 2 ** self.t2f.shape[0] - 1
	indices = (
	1 << np.arange(self.t2f.shape[0] * (1 + oriented))[:, None]
	& data[0].astype(np.int64)
	> 0
	)
	if oriented:
	oriented_facets = (
	self.t2f[indices[: self.t2f.shape[0]]],
	self.t2f[indices[self.t2f.shape[0] :]],
	)
	facets, index = np.unique(
	np.concatenate(oriented_facets), return_index=True
	)
	boundaries[subnames[2]] = OrientedBoundary(
	facets,
	np.concatenate(
	[
	np.ones_like(oriented_facets[0]),
	np.zeros_like(oriented_facets[1]),
	]
	)[index],
	)
	else:
	boundaries[subnames[2]] = np.sort(self.t2f[indices])

[gdmcbain]

No, this is silly: the facets as indexed in t2f are already oriented since the column pertains to the cell on the inside of the facet. There is no need to double the rows of t2f! See below.

[gdmcbain]

There was a draft of this in #874 but it was abandoned in favour of #875 and then ultimately #876.

[gdmcbain]

So a simpler alternative to #874 was mooted over there.

assuming we only want facets that have traces from inside cells, which I think is correct. Then every oriented facet is just a specific facet of a specific cell and we just need to put True in the corresponding slot of a boolean array shaped like .t2f and bit-pack the column of facets for each cell, which is really what we have been doing already since #681. Let's try that before merging this.

The present PR works for an arbitrary .ori but what I'm wondering now is whether that's only meaningful for internal facets which separate two cells rather than boundary facets which only pertain to one cell and so only have one meaningful orientation (hence our not bothering for years to orient .boundaries since they were mostly strictly boundaries).

Consider the example from the test.

scikit-fem/tests/test_basis.py

Lines 606 to 609 in 33a2513

	# m = MeshTri().refined(4)
	# m = m.with_boundaries({
	# 'mid': m.facets_around([5]),
	# })

We have for boundary = self.boundaries["mid"] that boundary.ori, boundary is:

array([0, 0, 1]), OrientedBoundary([  8,   9, 431])

which work in this example because self.f2t[[boundary.ori, boundary]] is

array([5, 5, 5])

These are the columns of self.t2f that we're to look in for the facets [ 8, 9, 431], since self.t2f[:, self.f2t[[boundary.ori, boundary]]] is

array([[  8,   8,   8],
       [431, 431, 431],
       [  9,   9,   9]])

I think this is

np.nonzero(boundary == self.t2f[:, self.f2t[[boundary.ori, boundary]]])[1]

i.e. here [0, 2, 1], but really we just want to put the boolean array

boundary == self.t2f[:, self.f2t[[boundary.ori, boundary]]]

into the columns of a False array shaped like .t2f for the appropriate cells, so the intermediate representation is:

ir = np.zeros_like(self.t2f).astype(np.bool)
ir[[np.nonzero(boundary == self.t2f[:, self.f2t[[boundary.ori, boundary]]])[1], self.f2t[[boundary.ori, boundary]]]] = True

since then ir.nonzero() is [[0, 1, 2], [5, 5, 5]], which is the indices of the facets of the cells and of the cells. This is bit-packed as

(1 << np.arange(self.refdom.nfacets)) @ ir

[gdmcbain]

        def encode_boundary(
            boundary: Union[ndarray, OrientedBoundary]
        ) -> ndarray:
            b = (
                boundary
                if isinstance(boundary, OrientedBoundary)
                else OrientedBoundary(boundary, np.zeros_like(boundary))
            )
            ir = np.zeros_like(self.t2f)
            columns = self.f2t[(b.ori, b)]
            ir[(np.nonzero(b == self.t2f[:, columns])[0], columns)] = 1
            return (1 << np.arange(self.refdom.nfacets)) @ ir

(or should ir be COO?).

For decoding,

facets = self.t2f[((1 << np.arange(self.refdom.nfacets))[:, None] & data[0]).astype(np.bool)]

[gdmcbain]

I think this is going to work and be better than #874 but I've run out of time for this week.

[gdmcbain]

I started on this in #875 but struggled with the unexpectedly pathological boundaries "left" and "right" in test_meshio_cycle_boundaries when internal_facets=False.

scikit-fem/tests/test_mesh.py

Lines 488 to 490 in 64da34c

	def test_meshio_cycle_boundaries(internal_facets, m):
	m = m.with_boundaries(_test_lambda, internal_facets)

It looks like this flag has been flipped and renamed boundaries_only.

scikit-fem/skfem/mesh/mesh.py

Lines 149 to 153 in 64da34c

	def with_boundaries(self,
	boundaries: Dict[str, Union[Callable[[ndarray],
	ndarray],
	ndarray]],
	boundaries_only: bool = True):

So the failing cases are those with internal facets which are thereby collectively unorientable to the extent that there are facets whose outsides are insides of other facets. I think this can be accounted for in the presently proposed t2f-based scheme but it did break an implicit assumption in the implementation.

Are there any applications for this kind of unorientable "boundary"? Anyway it does make a good rigorous test.

[gdmcbain]

Implemented satisfactorily and merged in #876.