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

General
1 

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

# 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)

#fdin == 2
fdim = domain.topology.dim - 1

# Get the facets @ the free end...
free_end_facets = mesh.locate_entities_boundary(domain, fdim, free_end)
'''
array([1960, 1961, 1962, 1963, 2072, 2073, 2074, 2075, 2172, 2173, 2174,
       2175, 2176, 2258, 2259, 2260, 2261, 2262, 2263, 2327, 2328, 2329,
       2330, 2331, 2332, 2333, 2376, 2377, 2378, 2379, 2406, 2407, 2408,
       2424, 2425, 2433], dtype=int32)
'''

# Extract the coordinates of the free end facets
free_end_coords = domain.geometry.x[free_end_facets] #<<--
'''
Traceback (most recent call last):
  File "/home/prusso/dolfinx-beam_moment/main.py", line 31, in <module>
    free_end_coords = domain.geometry.x[free_end_facets] #<<--
                      ~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^
IndexError: index 1960 is out of bounds for axis 0 with size 1029
'''

There are free_end_facets to be had however when extracting the nodes from domain.geometry.x it seems that those indices are being reported as out of range. Why is that?

2

I was able to get so far as a sample script that can do this. I am posting just as a propesed solution. If anyone has anything left to add, please post after me:


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

# 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
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]

# --- 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
# Each hexahedron has 8 vertices, so prepend each cell with the number 8
cells = domain.topology.connectivity(domain.topology.dim, 0).array
n_cells = cells.shape[0] // 8

# Properly format the connectivity array for PyVista (prepend 8 for each hexahedron)
cells_for_pv = np.hstack([np.full((n_cells, 1), 8), cells.reshape((n_cells, 8))]).flatten()

# Create the PyVista mesh for the entire domain
grid = pv.UnstructuredGrid(cells_for_pv, np.full(n_cells, 12), points)  # 12 is the VTK_HEXAHEDRON ID

# 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(free_end_coords), color='red', point_size=10, render_points_as_spheres=True)

# Show the plot
plotter.show()

print('done...')

Screenshot from 2024-09-01 17-34-35
Screenshot from 2024-09-01 17-34-351028×813 38.4 KB

3

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...')
1 Like