Getting Nodes of Boundary in an Ordered Structure

mesh
1

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[0],x[1]]) 

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!

2

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[0], 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):
	vertices = f2v.links(f)
	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:
	cell = f2c.links(f)[0]
	local_f = np.where(c2f.links(cell) == f)[0][0]
	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

image
image670×659 7.8 KB

and the DoF map

image
image712×668 9.99 KB

1 Like
3

Note that this can be simplified to:

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