# Mesh.locate\_entities\_boundary() does not seem to work if importing mesh nodes and cells

**URL:** <https://fenicsproject.discourse.group/t/mesh-locate-entities-boundary-does-not-seem-to-work-if-importing-mesh-nodes-and-cells/17976>\
**Category:** mesh\
**Tags:** boundary-conditions, meshio\
**Created:** [July 13, 2025, 10:16pm UTC](https://fenicsproject.discourse.group/t/mesh-locate-entities-boundary-does-not-seem-to-work-if-importing-mesh-nodes-and-cells/17976 "2025-07-13T22:16:41Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![alex24](https://avatars.discourse-cdn.com/v4/letter/a/4af34b/32.png) [@alex24](https://fenicsproject.discourse.group/u/alex24)\
**Post date:** [July 13, 2025, 10:16pm UTC](https://fenicsproject.discourse.group/t/mesh-locate-entities-boundary-does-not-seem-to-work-if-importing-mesh-nodes-and-cells/17976/1 "2025-07-13T22:16:41Z")

</div>

Hi,

When I generate a mesh using dolfinx, I can do the following and it generates the right boundary condition:

```auto
# Function to mark noslip bc
def noslip_boundary2(x):
    return np.isclose(x[2], z_val)

```

```auto
# No slip boundary condition 2
noslip = fem.Function(V)
facets2 = mesh.locate_entities_boundary(msh, fdim, noslip_boundary2)
dofs2 = fem.locate_dofs_topological((W.sub(0), V), fdim, facets2)
bc1 = fem.dirichletbc(noslip, dofs2, W.sub(0))

```

But when I import a mesh using meshio, using the same above two functions does not work. Am I missing a line to generate boundary facets? How is that done?

How can I generate the cell-node connectivity matrix `vol_mesh.cells[1].data` from msh? (to make a running example)?

```auto
vol_mesh = meshio.read('mesh.vtk')

points = vol_mesh.points
cells = vol_mesh.cells[1].data

c_el = ufl.Mesh(basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(points.shape[1],)))
msh = mesh.create_mesh(MPI.COMM_WORLD, cells, points, c_el)

fdim = msh.topology.dim - 1

msh.topology.create_connectivity(fdim, fdim+1)

```

Thank you,  
Alex

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [July 14, 2025, 3:14am UTC](https://fenicsproject.discourse.group/t/mesh-locate-entities-boundary-does-not-seem-to-work-if-importing-mesh-nodes-and-cells/17976/2 "2025-07-14T03:14:17Z")

</div>

> [@alex24](#):
>
> But when I import a mesh using meshio, using the same above two functions does not work. Am I missing a line to generate boundary facets? How is that done?

Could you be more specific? Does it throw an error or are you getting an unexpected result?

Either way, we would need a reproducible example with the full error trace to help you further

---

<div class="post-metadata">

**Author:** ![alex24](https://avatars.discourse-cdn.com/v4/letter/a/4af34b/32.png) [@alex24](https://fenicsproject.discourse.group/u/alex24)\
**Post date:** [July 14, 2025, 4:43pm UTC](https://fenicsproject.discourse.group/t/mesh-locate-entities-boundary-does-not-seem-to-work-if-importing-mesh-nodes-and-cells/17976/3 "2025-07-14T16:43:51Z")

</div>

Hi, sure, I am importing a vtk mesh:

> **[Dropbox](https://www.dropbox.com/scl/fi/euu0vkue888pa4zmc2et3/mesh.vtk?rlkey=mzqtdrjnz28pz7v0qpidf5lwe&st=8uqm3coq&dl=0)**

I tried the same code using a box generated by dolfinx and the code works. If I use the mesh I get this for the velocity magnitude:

 ![Screenshot from 2025-07-14 12-40-16](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/c/cbae1ce787c9fb9aae1d76583f4fa90a0a53de97.png)

The boundary conditions should be no slip at the cone and top cylinder (not the top cap) boundary facets, and at the bottom cap (which is not working at all)

Thank you for your time,  
Alex

This is the running code:

```auto
vol_mesh = meshio.read('mesh.vtk')

points = vol_mesh.points
cells = vol_mesh.cells[1].data

c_el = ufl.Mesh(basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(points.shape[1],)))
msh = mesh.create_mesh(MPI.COMM_WORLD, cells, points, c_el)

fdim = msh.topology.dim - 1

msh.topology.create_connectivity(fdim, fdim+1)

# Geometry variables
h_array_1 = 0.02 
h_array_2 = 0.035
h_gap = 0.03

rho_air = fem.Constant(msh, 1.225)
mu_air = fem.Constant(msh, 1.8e-3)

# Function to mark noslip
def noslip_boundary1(x):
    return np.logical_and(x[2] < h_array_2, x[2] > -h_array_1)

# Function to mark noslip
def noslip_boundary2(x):
    return np.isclose(x[2], (h_array_2 + h_gap))

# the same as old ufl.VectorElement
# the same as old ufl.FiniteElement
P2 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 2, shape=(msh.geometry.dim, ))
P1 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 1)

V, Q = fem.functionspace(msh, P2), fem.functionspace(msh, P1)

# Create the function space
TH = basix.ufl.mixed_element([P2, P1])
W = fem.functionspace(msh, TH)
W0, _ = W.sub(0).collapse()

# No slip boundary condition 1
noslip = fem.Function(V)
facets1 = mesh.locate_entities_boundary(msh, fdim, noslip_boundary1)
dofs1 = fem.locate_dofs_topological((W.sub(0), V), fdim, facets1)
bc0 = fem.dirichletbc(noslip, dofs1, W.sub(0))

# No slip boundary condition 2
facets2 = mesh.locate_entities_boundary(msh, fdim, noslip_boundary2)
dofs2 = fem.locate_dofs_topological((W.sub(0), V), fdim, facets2)
bc1 = fem.dirichletbc(noslip, dofs2, W.sub(0))

# Since for this problem the pressure is only determined up to a constant, we pin it at a point
zero = fem.Function(Q)
zero.x.array[:] = 0
dofs = fem.locate_dofs_geometrical(
    (W.sub(1), Q), lambda x: np.isclose(x.T, [0.0, 0.0, -h_array_1], atol=0.01).all(axis=1))
bc2 = fem.dirichletbc(zero, dofs, W.sub(1))

def body_force_fun(x):
    F_max = 2.0e-2
    z_offset = 0.02

    r_xy = np.sqrt(x[0]**2 + x[1]**2)

    a = 0.005
    max_val = a * np.exp(-0.5)
    Fr_val = np.exp(-r_xy**2/(2*a**2)) / max_val
        
    b = 0.005
    Fz_val = np.exp(-(x[2]-z_offset)**2/(2*b**2))
    F_mag = F_max * Fr_val * Fz_val
      
    # -ve sign to get clockwise rotation from array top
    return -F_mag * np.vstack((x[1], -x[0], np.zeros_like(x[2])))
        
Force_per_unit_vol = fem.Function(V)
Force_per_unit_vol.interpolate(lambda x: body_force_fun(x))

# Collect Dirichlet boundary conditions
bcs = [bc0, bc1, bc2]

# Define variational problem
w = fem.Function(W)

(u, p) = ufl.split(w)
(v, q) = ufl.TestFunctions(W)

# Tentative velocity step

# https://github.com/bragostin/Fenics/blob/main/HTC_SquareDuct_Steady_Neumann_Minimal_Pub.py

# https://math.stackexchange.com/questions/1952588/
# weak-form-of-steady-navier-stokes-equations-with-special-boundary-condition

F = rho_air*inner(grad(u)*u, v)*dx + \
    mu_air*inner(grad(u), grad(v))*dx + \
    inner(Force_per_unit_vol, v)*dx - \
    inner(p, div(v))*dx - \
    div(u)*q*dx

dw = ufl.TrialFunction(W)
dF = ufl.derivative(F, w, dw)

problem = NonlinearProblem(F, w, bcs=bcs, J=dF)

solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
#solver.max_it = 2
solver.report = True
solver.error_on_nonconvergence = False

# Custom solver options (not needed for code readibility but makes code 9x faster)
opts = PETSc.Options() # type: ignore
ksp = solver.krylov_solver
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

# Compute the solution
log.set_log_level(log.LogLevel.INFO)
start = time.time()
solver.solve(w)
end = time.time()

print('\n')
print(msh.topology.index_map(msh.topology.dim).size_local)
print((end - start) / 60)

# Split the mixed solution and collapse
u = w.sub(0).collapse()
p = w.sub(1).collapse()

# Write the solution to file
with io.XDMFFile(MPI.COMM_WORLD, "pressure.xdmf", "w") as pfile_xdmf:
    p.x.scatter_forward()
    pfile_xdmf.write_mesh(msh)
    pfile_xdmf.write_function(p)

with io.XDMFFile(MPI.COMM_WORLD, "velocity.xdmf", "w") as ufile_xdmf:
    u.x.scatter_forward()
    P1 = basix.ufl.element("Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,))
    u1 = fem.Function(fem.functionspace(msh, P1))
    u1.interpolate(u)
    ufile_xdmf.write_mesh(msh)
    ufile_xdmf.write_function(u1)
    

with io.XDMFFile(MPI.COMM_WORLD, "body_force.xdmf", "w") as ufile_xdmf:
    Force_per_unit_vol.x.scatter_forward()
    P1 = basix.ufl.element("Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,))
    u1 = fem.Function(fem.functionspace(msh, P1))
    u1.interpolate(Force_per_unit_vol)
    ufile_xdmf.write_mesh(msh)
    ufile_xdmf.write_function(u1)

```

---

<div class="post-metadata">

**Author:** ![alex24](https://avatars.discourse-cdn.com/v4/letter/a/4af34b/32.png) [@alex24](https://fenicsproject.discourse.group/u/alex24)\
**Post date:** [July 14, 2025, 6:13pm UTC](https://fenicsproject.discourse.group/t/mesh-locate-entities-boundary-does-not-seem-to-work-if-importing-mesh-nodes-and-cells/17976/4 "2025-07-14T18:13:13Z")

</div>

I sorted it out, I was not considering the meshing program is not placing the facets very accurately on the surfaces, the marking functions should be as follows:

```auto
mark_eps = 2.0e-3

# Function to mark noslip
def noslip_boundary1(x):
    return np.logical_and(x[2] < (h_array_2 + mark_eps), np.sqrt(x[0]**2 + x[1]**2) > (r2-mark_eps))

# Function to mark noslip
def noslip_boundary2(x):
    return np.isclose(x[2], (h_array_2 + h_gap), mark_eps)

```