# Getting Nodes of Boundary in an Ordered Structure

I am implementing a contact model that calculated a penalty term backed on the level of penetration into a boundary face. Right now I have it working in 2D convex objects by pulling the DOFs associated with the boundary and organizing them by the angle to the horizontal:

``````import numpy as np
import ufl
from dolfinx import fem, mesh, plot

domain = mesh.create_rectangle(MPI.COMM_WORLD,[[0.0,0.0],[1,1]],[10,10],mesh.CellType.triangle)
V_vec = fem.VectorFunctionSpace(domain,("CG",2))

# Obtaining the boundary node indexes
edges = mesh.locate_entities_boundary(domain,1,lambda x: np.ones_like(x[0,:]))
edge_dofs = fem.locate_dofs_topological(V_vec,1,edges,1)

X0 = fem.Function(V_vec)         # Local Reference Configuration of the domain
X0.interpolate(lambda x: [x,x])

NN = len(X0.x.array)                 # number of nodes
x_inds = np.arange(0,NN,2)      # indices of the x coordinates
y_inds = np.arange(1,NN+1,2)  # indices of the y coordinates

y_ = X0.x.array[y_inds[edge_dofs]]  # getting x and y coordinates of the edge nodes
x_ = X0.x.array[x_inds[edge_dofs]]
mx = np.mean(x_)
my = np.mean(y_)
angles = np.arctan2((y_ - my),(x_ - mx)) # sorting the nodes by increasing angle

edges_ind = edge_dofs[np.argsort(angles)]
boundary = np.zeros((len(self.edges_ind),2))

# re-organizing the nodes by increasing angle
boundary[:,0] = X0.x.array[x_inds[edges_ind]]
boundary[:,1] = X0.x.array[y_inds[edges_ind]]
``````

With this organized structure, I can loop through all of the nodes and find the penetration beyond a line segment captured by two adjacent nodes. However, this method will not work for concave hulls and does not easily generalize to 3D meshes.

Are there built-in tools that allow for the nodes associated with all boundary facet to be grouped in a way that each facet can be looped over? (ie. a function that returns a structure whose elements are sets of nodes each associated with a single boundary face, and hopefully that are organized in a known way). It seems like there are some tools that do something similar in FEniCS (something like BoundaryMesh), but not yet in FEniCS-x or have changed their implementation.

Thanks!

Consider the following (thanks to @dokken for additional help answering this):

``````import dolfinx
import numpy as np
from mpi4py import MPI

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2)

# Find exterior facets
exterior_facets = dolfinx.mesh.locate_entities_boundary(
mesh, mesh.topology.dim-1, lambda x: np.ones_like(x, dtype=np.int8))
exterior_facets.sort()

# Compute their midpoints
ext_facet_midpoints = dolfinx.mesh.compute_midpoints(
mesh, mesh.topology.dim-1, exterior_facets)

# Print vertices on each facet and its midpoint
f2v = mesh.topology.connectivity(mesh.topology.dim-1, 0)
for f, f_midpoint in zip(exterior_facets, ext_facet_midpoints):
print(f"facet {f} has vertices {vertices} and midpoint {f_midpoint}")

# Find DoFs associated which each facet
V = dolfinx.fem.FunctionSpace(mesh, ("CG", 2))
edl = V.dofmap.dof_layout

c2f = mesh.topology.connectivity(mesh.topology.dim, mesh.topology.dim-1)
f2c = mesh.topology.connectivity(mesh.topology.dim-1, mesh.topology.dim)
for f in exterior_facets:
local_f = np.where(c2f.links(cell) == f)
cell_dofs = V.dofmap.cell_dofs(cell)
single_facet_dofs = cell_dofs[
edl.entity_closure_dofs(mesh.topology.dim-1, local_f)]
print(f"Facet {f} has DoFs: {single_facet_dofs}")
``````

Which outputs

``````facet 0 has vertices [0 1] and midpoint [0.75 0.   0.  ]
facet 2 has vertices [1 2] and midpoint [1.   0.25 0.  ]
facet 5 has vertices [4 0] and midpoint [0.25 0.   0.  ]
facet 7 has vertices [2 5] and midpoint [1.   0.75 0.  ]
facet 10 has vertices [4 6] and midpoint [0.   0.25 0.  ]
facet 12 has vertices [7 5] and midpoint [0.75 1.   0.  ]
facet 14 has vertices [6 8] and midpoint [0.   0.75 0.  ]
facet 15 has vertices [8 7] and midpoint [0.25 1.   0.  ]
Facet 0 has DoFs: [0 1 5]
Facet 2 has DoFs: [1 2 3]
Facet 5 has DoFs: [ 9  0 11]
Facet 7 has DoFs: [ 2 12 13]
Facet 10 has DoFs: [ 9 15 17]
Facet 12 has DoFs: [18 12 19]
Facet 14 has DoFs: [15 22 24]
Facet 15 has DoFs: [22 18 23]
``````

Which should match up with the facet indices

and the DoF map

1 Like

Note that this can be simplified to:

``````exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
``````
3 Likes