Extracting boundary facet connectivity and facet node coordinates from mesh

mesh
,
1

Hi, is there a way of extracting the boundary facet connectivity and the facet node coordinates from a mesh? For example from the domain object below in the linear elasticity example.

This is to pick a facet using another algorithm that needs the surface mesh and then setting that single facet as a zero displacement Diritchlet BC.

Thank you,
Alex

# Scaled variable
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
L = 1
W = 0.2
mu = 1
rho = 1
delta = W / L
gamma = 0.4 * delta**2
beta = 1.25
lambda_ = beta
g = gamma

domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],
                         [20, 6, 6], cell_type=mesh.CellType.hexahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))


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


fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)

u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

ds = ufl.Measure("ds", domain=domain)


def epsilon(u):
    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)

def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, 0, -rho * g)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds

problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
2

See for instance: How to extract the facet co-ordinates at the free end of a cantilever beam - #3 by dokken
as the question seems very similar.

If that post doesn’t answer your question feel free to elaborate here.

3

Hi Jorgen,

Thank you for the prompt reply; I think I have a solution based on the referenced previous post (below). However it still uses the function free_end(). I would like to get all the boundary facets. I can modify the function to return 1.0+x[0]*0.0 but is there a better way of doing it?

Is there a function instead of mesh.locate_entities_boundary() that just returns all boundary facets without a marker function?

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

# Create the mesh
domain = mesh.create_unit_cube(MPI.COMM_WORLD, nx=1, ny=1, nz=1, cell_type=mesh.CellType.tetrahedron)

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

# 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)
print(str(free_end_facets) + '\n')

# get the nodes of all vertices on the free end
free_end_connectivity = mesh.entities_to_geometry(domain, fdim, free_end_facets)
print(str(free_end_connectivity) + '\n')

# Get only unique vertices
free_end_nodes = np.unique(free_end_connectivity.reshape(-1))
print(str(free_end_nodes) + '\n')

# Extract the coordinates of the nodes associated with the facets
free_end_coords = domain.geometry.x[free_end_nodes]
print(str(free_end_coords) + '\n')


# Change connectivity numbering to match order from 0 of extracted nodes
connectivity_copy = np.copy(free_end_connectivity)

for n in np.arange(free_end_nodes.size):
    free_end_connectivity[connectivity_copy == free_end_nodes[n]] = n
    
print(free_end_connectivity)
4

All boundary facets can be found with dolfinx.mesh.exterior_facet_indices(mesh.topology)

5

Great it works!

Replaced part of the code with the following lines and took out free_end() function

# 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.exterior_facet_indices(domain.topology)
print(str(free_end_facets) + '\n')