How to extract the facet co-ordinates at the free end of a cantilever beam

[dokken]

@prusso630

# Find facets on the boundary where x = L (free end)
fdim = domain.topology.dim - 1
free_end_facets = mesh.locate_entities_boundary(domain, fdim, free_end)

# Get the vertices (nodes) connected to the facets
facet_to_vertex_connectivity = domain.topology.connectivity(fdim, 0)

# Extract the vertex indices for the free end facets
free_end_vertex_indices = np.unique(np.concatenate([facet_to_vertex_connectivity.links(facet) for facet in free_end_facets]))

# Extract the coordinates of the nodes associated with the facets
free_end_coords = domain.geometry.x[free_end_vertex_indices]

This is not an always correct solution (especially in parallel).
You should use dolfinx.mesh.entities_to_geometry and dolfinx.plot.vtk_mesh:

import pyvista as pv
from mpi4py import MPI
from dolfinx import mesh, plot
import numpy as np

# Geometry and material properties
L = 120.0
W = 4.0
H = 8.0

# Create the mesh
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, H])],
                         [20, 6, 6], cell_type=mesh.CellType.hexahedron)

def free_end(x):
    return np.isclose(x[0], L)

# Find facets on the boundary where x = L (free end)
fdim = domain.topology.dim - 1
domain.topology.create_connectivity(fdim, fdim+1)
free_end_facets = mesh.locate_entities_boundary(domain, fdim, free_end)

# get the nodes of all vertices on the free end
free_end_nodes = mesh.entities_to_geometry(domain, fdim, free_end_facets)

# Get only unique vertices
free_end_nodes = np.unique(free_end_nodes.reshape(-1))



# --- PyVista Visualization ---
# Create a PyVista mesh object for the wireframe
points = domain.geometry.x  # All points in the mesh

# Get the cell (hexahedron) connectivity for PyVista
# Get the PyVista cell connectivity
pv_mesh = plot.vtk_mesh(domain, dim=domain.topology.dim)


# Create the PyVista mesh for the entire domain
domain.topology.create_connectivity(fdim+1, 0)
grid = pv.UnstructuredGrid(*pv_mesh)

# Start a PyVista plot
plotter = pv.Plotter()

# Add the wireframe of the mesh (domain)
plotter.add_mesh(grid, color='lightblue', show_edges=True, style='wireframe', opacity=0.5)

# Highlight the free end nodes (in red)
plotter.add_mesh(pv.PolyData(points[free_end_nodes]), color='red', point_size=10, render_points_as_spheres=True)

# Show the plot
plotter.show()

print('done...')