boundaries in built-in meshes

[gdmcbain]

From #776, @gatling-nrl (2021-11-14)

I think what you and I are both really doing in this thread is making the very strong case for labelling points during mesh creation instead of finding them with facets_satsfying.

To begin with, something like MeshLine might have boundaries "left" and "right" corresponding to x == p.min() and x == p.max(), rectangular meshes generated by multiplying two MeshLines together or otherwise might have those and "upper" and "lower", cuboids might also have "front" and "back".

Should this be by universal or optional? (Could there be a reason for wanting an empty .boundaries attribute?)

[kinnala]

I wonder if we first want to reduce the size of boundaries and subdomains by implementing them internally as facet/element flags. Something like:

In [8]: m = MeshTri()                                                                               

In [9]: m.facets                                                                                    
Out[9]: 
array([[0, 0, 1, 1, 2],
       [1, 2, 2, 3, 3]])

In [10]: m.boundary_facets()                                                                        
Out[10]: array([0, 1, 3, 4])

and then m._boundaries = np.array([1, 2, 0, 1, 1], dtype=np.uint8) and m._boundaries_flags = {'all': [1, 2], 'left': [2]}.

Moreover, we need to support preserving boundaries in refined, i.e. #104.

[kinnala]

I guess the size reduction can be done later if we find out that to be an issue.

[gdmcbain]

#104 would be certainly nice (for other reasons too) but isn't necessary here; rectangles and cuboids can be generated using multiplication and .init_tensor without .refine if one wants access to the predefined boundaries (e.g. to avoid agonizing about floating-point comparisons). One could still use .refine before #104 in cases where the predefined boundaries aren't required; e.g.

scikit-fem/docs/examples/ex01.py

Line 10 in 62f72d3

m = MeshTri().refined(4)

[gdmcbain]

So the flags idea would require finding all the intersections of boundaries defined in the current manner (Dict[str, ndarray])? Hmm. Yes it might be a good optimization. But is it orthogonal to the present issue? Is changing the type of the present attribute ._boundaries bad?

[kinnala]

You are correct. I was just worried that if we do supply these by default it will increase the memory footprint of the resulting objects. However, that may not be too large after all.

[kinnala]

I made some tests toward this in #792. It also includes the preservation of tags in MeshTri1._uniform_refine.

[kinnala]

It's actually kind of challenging to add those tags to some of the constructors without relying on facets_satisfying. It requires some knowledge of how the facets are ordered inside Refdom objects.

[gdmcbain]

Ah.
First thought is that it might be reasonable to begin by relying on facets_satisfying. We know that this works correctly, and it's good to begin correctly and worry about optimization later if required. If it is slow or bulky, it might be made optional. Also, hiding the use of the floating-point comparison in the library rather than forcing it into user-level scripts as in docs/examples could be good because floating-point comparison is ugly and freaks people out.
Second thought is that while I'd hoped that an a priori knowledge of which facets would be on which boundary would be available, if not I thought of trying to just use .facets_satisfying on a sequence of very small examples and seeing whether I could spot the pattern. I'm hoping that they might be something like an arithmetic progression? Or would that be fragile?

[kinnala]

https://github.com/kinnala/scikit-fem/blob/master/skfem/refdom.py#L63

This line will tell you that facets for MeshTri1 are created from MeshTri1.p (nowadays called MeshTri1.doflocs internally) by taking (p[t[0]], p[t[1]]), (p[t[1]], p[t[2]]), (p[t[0]], p[t[2]]). However, this process will result in duplicates due to the two equal internal edges coming from two triangles. Hence, we need to run np.unique. Because of calling np.unique, the final order of the facets in the data structure is a bit unpredictable. Well, not unpredictable, but not easy either.

You most likely have to first construct MeshTri1 and afterwards look at the contents of f2t and t2f to figure out the correct facets. However, that's how it's done in MeshTri1._uniform_refine so I guess it's not a problem.

[kinnala]

I think I'll just go with #792. Maybe not start relying on it yet in the examples but only after we have some experiences on using it in practice?

[gdmcbain]

Yeah, O. K.

The other cases that I thought might be worthwhile and not too hard are MeshLine and init_tensor, but I can take a look at those next time i need them.