# Compute_colliding_cells are far from my points

I try to add a jet boundary inside a flow around cylinder example.The cylinder’s centre is （0.3，0.2，0）and its radius is 0.05,the domain’s shape is [0.41 * 2.2].I try to add two jet boundaries to the top and bottom，which means the coordinate of the jet is（0.3，0.25, 0）and（0.3，0.15, 0）.But jets didn’t be created as expected：

Here are the codes I creat points：

``````def positions(jet_position):
jet_coordinates = np.zeros((3,len(jet_positions)))
jet_x = [i[0]for i in jet_positions]
jet_x = [i[0]for i in jet_positions]
jet_coordinates[0] = jet_x
jet_coordinates[1] = jet_y
return jet_coordinates
jet_positions = [(0.3,0.15,0),(0.3,0.25,0)]
jet_coordinates = positions(jet_positions)
``````

Then compute the colliding cells

``````from dolfinx.fem.geometry import BoundingBoxTree, compute_collisions, compute_colliding_cells

bb_tree = geometry.BoundingBoxTree(mesh,mesh.geometry.dim)
jet_cells = []
jet_points = []
jet_cell_candidates = geometry.compute_collisions(bb_tree, jet_coordinates.T)
jet_colliding_cells = geometry.compute_colliding_cells(mesh, jet_cell_candidates, jet_coordinates.T)
for i, jet_point in enumerate(jet_coordinate.T):
jet_points.append(jet_point)
``````

``````class JetVelocity():
def __init__(t, jet_x, jet_y):
self.t_jet = t
self.jet_x = jet_x
self.jet_y = jet_y
def __call__(self, x):
values = np.zeros((2,x.shape[1]),dtype=PETSC.ScalarType)
values[0] = self.jet_x
values[1] = self.jet_y
return values
u_jet = Function(V)
jet_velocity = JetVelocity(t,jet_x,jet_y)
u_jet.interpolate(jet_velocity)
bcu_jets = diricheletbc(u_jet, locate_dofs_topological(V,fdim,jet_cells))
``````

Are there any problems?Any advices would be appreciated

As you have fragmented your code into snippets, it is to me very unclear what you want to do here.
Are you tryinbg to apply the `JetVelocity` to a subset of cells, which should be those intersecting with the jet-coordinates)?

Note that you use `fdim` (which I assume is `mesh.topology.dim-1`, and not `mesh.topology.dim` in

which could explain your issues. However, without a minimal reproducible example, this is just an educated guess.

1 Like

Hi @dokken ，thanks for your reply.Your guess is what I want to do.And here are the reproducible example：

``````import gmsh
import numpy as np
import sys
from mpi4py import MPI

from petsc4py import PETSc

from dolfinx import geometry
from dolfinx.fem import (Constant, Function, FunctionSpace,
assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
create_vector, create_matrix, set_bc)
from dolfinx.geometry import BoundingBoxTree, compute_collisions, compute_colliding_cells
from dolfinx.io import (XDMFFile, distribute_entity_data, gmshio)

from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,

class simulation:
def __init__(self, c_x, c_y, o_x, o_y, r, r2,
obstacle_type='cylinder',
save_mesh='false',
save='false',
visualize='false'):

self.L = 2.2
self.H = 0.41
self.c_x = c_x
self.c_y = c_y
self.r = r
self.gdim = 2
self.model_rank = 0
self.mesh_comm = MPI.COMM_WORLD
self.save = save
self.o_x = o_x
self.o_y = o_y
self.r2 = r2
self.obstacle_type = obstacle_type
self.save_mesh = save_mesh
self.visualize = visualize

def generate_mesh(self):
gmsh.initialize()

rectangle = gmsh.model.occ.addRectangle(0, 0, 0, self.L, self.H, tag=1)
obstacle1 = gmsh.model.occ.addDisk(self.c_x, self.c_y, 0, self.r, self.r)
pre_fluid = gmsh.model.occ.cut([(self.gdim, rectangle)], [(self.gdim, obstacle1)])
if self.r2 == 0:
fluid = pre_fluid
else:
if self.obstacle_type == 'cylinder':
obstacle2 = gmsh.model.occ.addDisk(self.o_x, self.o_y, 0, self.r2, self.r2)
elif self.obstacle_type == 'rectangular':
obstacle2 = gmsh.model.occ.addRectangle(self.o_x, self.o_y, 0, self.r2, self.r2)
fluid = gmsh.model.occ.cut([(self.gdim, rectangle)], [(self.gdim, obstacle2)])

gmsh.model.occ.synchronize()
fluid_marker = 1

volumes = gmsh.model.getEntities(dim=self.gdim)
assert (len(volumes) == 1)
gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")

self.inlet_marker, self.outlet_marker, self.wall_marker, self.obstacle_marker = 2, 3, 4, 5
inflow, outflow, walls, obstacle = [], [], [], []

boundaries = gmsh.model.getBoundary(volumes, oriented=False)
for boundary in boundaries:
center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
if np.allclose(center_of_mass, [0, self.H / 2, 0]):
inflow.append(boundary[1])
elif np.allclose(center_of_mass, [self.L, self.H / 2, 0]):
outflow.append(boundary[1])
elif np.allclose(center_of_mass, [self.L / 2, self.H, 0]) or np.allclose(center_of_mass, [self.L / 2, 0, 0]):
walls.append(boundary[1])
else:
obstacle.append(boundary[1])
gmsh.model.setPhysicalName(1, self.wall_marker, "Walls")
gmsh.model.setPhysicalName(1, self.inlet_marker, "Inlet")
gmsh.model.setPhysicalName(1, self.outlet_marker, "Outlet")
gmsh.model.setPhysicalName(1, self.obstacle_marker, "Obstacle")

res_min = self.r / 3
gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", obstacle)
gmsh.model.mesh.field.setNumber(threshold_field, "IField", distance_field)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMin", res_min)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMax", 0.25 * self.H)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", self.r)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", 2 * self.H)
gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field])
gmsh.model.mesh.field.setAsBackgroundMesh(min_field)

gmsh.option.setNumber("Mesh.Algorithm", 8)
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.generate(self.gdim)
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")
mesh, _, ft = gmshio.model_to_mesh(gmsh.model, self.mesh_comm, self.model_rank, gdim=self.gdim)
ft.name = "facet markers"
return mesh, ft

def compute(self, mesh,  ft, jet):
jet_position = jet.get("position")
jet_x = jet.get("jet_x")
jet_y = jet.get("jet_y")
global u_probes_t, p_probes_t
bb_tree = geometry.BoundingBoxTree(mesh, mesh.topology.dim)
jet_cells = []
jet_points = []
jet_cell_candidates = geometry.compute_collisions(bb_tree, jet_position.T)
jet_colliding_cells = geometry.compute_colliding_cells(mesh, jet_cell_candidates, jet_position.T)
for i, jet_point in enumerate(jet_position.T):
jet_points.append(jet_point)

t = 0
T = 7  # Final time
dt = 1 / 1600  # Time step size
num_steps = 100 # It is just for testing
k = Constant(mesh, PETSc.ScalarType(dt))
mu = Constant(mesh, PETSc.ScalarType(0.001))  # Dynamic viscosity
rho = Constant(mesh, PETSc.ScalarType(1))  # Density

v_cg2 = VectorElement("CG", mesh.ufl_cell(), 2)
s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)

fdim = mesh.topology.dim - 1

# Define boundary conditions
class InletVelocity():
def __init__(self, t):
self.t = t

def __call__(self, x):
values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = 4 * 1.5 * np.sin(self.t * np.pi / 8) * x[1] * (0.41 - x[1]) / (0.41 ** 2)
return values

class JetVelocity():
def __init__(self, t, jet_x, jet_y):
self.t_jet = t
self.jet_x = jet_x
self.jet_y = jet_y

def __call__(self, x):
values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = self.jet_x
values[1] = self.jet_y
return values
# Inlet
u_inlet = Function(V)
u_jet = Function(V)
inlet_velocity = InletVelocity(t)
jet_velocity = JetVelocity(t, jet_x, jet_y)
u_inlet.interpolate(inlet_velocity)
u_jet.interpolate(jet_velocity)
bcu_inflow = dirichletbc(u_inlet, locate_dofs_topological(V, fdim, ft.find(self.inlet_marker)))
bcu_jets = dirichletbc(u_jet, locate_dofs_topological(V, fdim, jet_cells))
# Walls
u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bcu_walls = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(self.wall_marker)), V)
# Obstacle
bcu_obstacle = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(self.obstacle_marker)), V)
bcu = [bcu_inflow, bcu_obstacle, bcu_walls, bcu_jets]
# Outlet
bcp_outlet = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(Q, fdim, ft.find(self.outlet_marker)), Q)
bcp = [bcp_outlet]

u = TrialFunction(V)
v = TestFunction(V)
u_ = Function(V)
u_.name = "u"
u_s = Function(V)
u_n = Function(V)
u_n1 = Function(V)
p = TrialFunction(Q)
q = TestFunction(Q)
p_ = Function(Q)
p_.name = "p"
phi = Function(Q)

f = Constant(mesh, PETSc.ScalarType((0, 0)))
F1 = rho / k * dot(u - u_n, v) * dx
F1 += inner(dot(1.5 * u_n - 0.5 * u_n1, 0.5 * nabla_grad(u + u_n)), v) * dx
F1 += 0.5 * mu * inner(grad(u + u_n), grad(v)) * dx - dot(p_, div(v)) * dx
F1 += dot(f, v) * dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))
A1 = create_matrix(a1)
b1 = create_vector(L1)

L2 = form(-rho / k * dot(div(u_s), q) * dx)
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)

a3 = form(rho * dot(u, v) * dx)
L3 = form(rho * dot(u_s, v) * dx - k * dot(nabla_grad(phi), v) * dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)
# +
# Solver for step 1
solver1 = PETSc.KSP().create(mesh.comm)
solver1.setOperators(A1)
solver1.setType(PETSc.KSP.Type.BCGS)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.JACOBI)

# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.MINRES)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.HYPRE)
pc2.setHYPREType("boomeramg")

# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)

xdmf = XDMFFile(mesh.comm, "flow_around_cylinder.xdmf", "w")
xdmf.write_mesh(mesh)
for i in range(num_steps):
# Update current time step
t += dt
# Update inlet velocity
inlet_velocity.t = t
u_inlet.interpolate(inlet_velocity)

# Step 1: Tentative velocity step
A1.zeroEntries()
assemble_matrix(A1, a1, bcs=bcu)
A1.assemble()
with b1.localForm() as loc:
loc.set(0)
assemble_vector(b1, L1)
apply_lifting(b1, [a1], [bcu])
set_bc(b1, bcu)
solver1.solve(b1, u_s.vector)
u_s.x.scatter_forward()

# Step 2: Pressure corrrection step
with b2.localForm() as loc:
loc.set(0)
assemble_vector(b2, L2)
apply_lifting(b2, [a2], [bcp])
set_bc(b2, bcp)
solver2.solve(b2, phi.vector)
phi.x.scatter_forward()

p_.vector.axpy(1, phi.vector)
p_.x.scatter_forward()

# Step 3: Velocity correction step
with b3.localForm() as loc:
loc.set(0)
assemble_vector(b3, L3)
solver3.solve(b3, u_.vector)
u_.x.scatter_forward()
xdmf.write_function(u_, t)
xdmf.write_function(p_, t)
xdmf.close()
gmsh.clear()

if __name__ == '__main__':
jet_positions = [(0.3, 0.15, 0), (0.3, 0.25, 0)]
def positions(jet_positions):
jet_coordinates = np.zeros((3, len(jet_positions)))
jet_x = [i[0] for i in jet_positions]
jet_y = [i[1] for i in jet_positions]
jet_coordinates[0] = jet_x
jet_coordinates[1] = jet_y
return jet_coordinates

jet_coordinates = positions(jet_positions)

jet = {
"jet_x": 1,
"jet_y": 1,
"position": jet_coordinates
}
simu = simulation(c_x=0.3, c_y=0.2, o_x=0.3, o_y=0.2, r=0.05, r2=0)
mesh, ft = simu.generate_mesh()
simu.compute(mesh=mesh, ft=ft, jet=jet)

``````

As I said, this should not be fdim, but mesh.topology.dim

Hi @dokken ，could you explain that why `dim` in `bcu_inflow` is 1, and `dim` in `bcu_jet` is 2 ?

Because you are using `facet_tags` (`ft`) in `bcu_inflow`, which consists of entities of dimension `mesh.topology.dim-1`. As you are sending in cells in `bcu_jets`, you need to use `mesh.topology.dim`.

1 Like