How to order boundary points/edges

[griff10000]

Is there a way of ordering the set of boundary points of boundary edges, clockwise/anticlockwise?

[griff10000]

I attach a zipped file containing a function sortEdgeList() that orders boundary nodes. It is an amended version of code originally written by Francis Laclé. It manages to order boundaries for spatial domains with one or more concavities, which fail under the simplistic approach of ordering angles from a center point within the domain. However, as it stands, it does not handle boundaries for spatial domains with multiple holes. The zipped file also includes a test file using irregular polynomials. I hope these files a) will be useful for people who need this facility, and b) may help this facility eventually being incorporated into dmsh. boundaryOrderTest.zip

[nschloe]

I don't know if this will ever land in dmsh. As you correctly point out, holes are an issue. I'm not familiar with applications that would need this kind of functionality either. Can you explain what you're doing?

[griff10000]

The application i am working on is the solution of partial differential equations (PDEs) using meshless methods based on radial basis functions (RBFs). In particular, those with irregular spatial domains. The requirement for sorted boundary nodes arises when,

a) Neumann boundary conditions are imposed and outward normals need to be calculated. This involves differentiating along the boundary which requires the nodes to be ordered;

b) the domain is not entirely convex so, for plotting purposes, the triangles in the concave areas need to be identified and eliminated where necessary.

I have identified the python 'graph theory' package networkx which has the ability to order a list of nodes, as follows:

import dmsh
import meshplex
import networkx as nx # For ordering boundary nodes
...
# Generate mesh object
mesh = meshplex.MeshTri(pts, cells)

# Get boundary edge indices
bei = mesh.boundary_edges    # boundary edge indicies
et  = mesh.edges['points']   # total edges
eb  = et[bei,:]              # boundary edges

# Get boundary nodes - ordered
G = nx.Graph()
G.add_edges_from(eb)
b_index = nx.cycle_basis(G,root=3)[0] # nx1 array of ordered nodes

I attache an example that demonstrates the whole process for a circular domain. orderBoundaryNodes.zip

The package networkx apparently works when there are multiple boundaries, but I haven't experimented with this yet.

[griff10000]

I have now tested package networkx with 2 boundaries and it works well. See plot and (updated) zipped file for code orderBoundaryNodes_2.zip

[nschloe]

Nice that you got it to work, although I'm not a big fan of the ordering approach. That's because it doesn't generalize to more complex triangular domains, let alone tetrahedral meshes. A more general approach is to count the adjacent cells for every facet. (This can be done with one np.unique call.) Facets with one adjacent cell are boundary facets. The normal can be computed locally by projecting the non-boundary point to the boundary facet. I would imagine that this approach is a lot faster, and I'm thinking about adding some helper tools to meshplex.

[krober10nd]

Nico’s right this doesn’t generalize to 3d but here’s how you could order them too (based on a walk)

https://github.com/krober10nd/SeismicMesh/blob/06a40992c3073710a2ba34c2613e46d9842cc0ca/SeismicMesh/geometry/utils.py#L328

[nschloe]

For sure, but this doesn't help much with outer normals, right?

[krober10nd]

Hm well if the boundary edges are ordered in a common direction then the boundary normal calculation would be consistent. But this would only work in 2d

I like your approach that landed in meshplex since it’s general to all dimensions

However even then if the normals are all calculated similarly, this doesn’t imply the facets/edges will have some spatial order ? I guess I was just answering the OP

[nschloe]

this doesn’t imply the facets/edges will have some spatial order

No particular order is necessary when determining outer normals since every boundary simplex has exacltly one point which isn't part of the boundary.

[krober10nd]

Sure. As an aside, in 2D a boundary traversal indeed has applications (visualization).

Anyway apparently the OP just needed to calculate boundary normals and did not need to order the edges which is what your approach answers.