Extracting values along the boundary

I solved the Navier Stokes equation in the following domain:

I created it using the gmsh module and used physical groups to denote the boundary conditions, inflow and outflow locations.

Now suppose I want to extract the velocity and pressure values along the physical groups that I defined earlier.

I do understand how they did this in here:

tol = 0.001 
y = np.linspace(-1 + tol, 1 - tol, 101)
points = np.zeros((3, 101))
points[1] = y
u_values = []
p_values = []

from dolfinx import geometry
bb_tree = geometry.BoundingBoxTree(domain, domain.topology.dim)

cells = []
points_on_proc = []
# Find cells whose bounding-box collide with the the points
cell_candidates = geometry.compute_collisions(bb_tree, points.T)
# Choose one of the cells that contains the point
colliding_cells = geometry.compute_colliding_cells(domain, cell_candidates, points.T)
for i, point in enumerate(points.T):
    if len(colliding_cells.links(i))>0:

points_on_proc = np.array(points_on_proc, dtype=np.float64)
u_values = uh.eval(points_on_proc, cells)
p_values = pressure.eval(points_on_proc, cells)

But in this case, they have manually created the points set. However, I would like to use the names of the physical groups that I used in creating the mesh.

For example, I used:

in_flow_marker  = 100
gmsh.model.addPhysicalGroup(1, [line_7], in_flow_marker)

to denote the left-hand side vertical line. Also, I can observe that the following code will give the corresponding gmsh.points along the inflow boundary:

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)

locate_dofs_topological(Q, fdim, facet_tags.find(in_flow_marker))


So, my question is:
what is the appropriate way to utilize these gmsh.points into the points variable defined at the beginning.
I appreciate your help very much!

For the completeness of my question (and if someone wanted to know how I did that), I’m posting the full code that I used to create the mesh and to solve it.

import numpy as np
import gmsh
from dolfinx.io.gmshio import model_to_mesh
from mpi4py import MPI

import pyvista
from dolfinx import plot

from ufl import (FacetNormal, FiniteElement, Identity,TestFunction, TrialFunction, VectorElement,
                 div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym)
from dolfinx.fem import Constant,Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,locate_dofs_topological
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from dolfinx.io import XDMFFile
from dolfinx.plot import create_vtk_mesh
from petsc4py import PETSc
import matplotlib.pyplot as plt
from dolfinx import geometry

A = np.array([0,0])
B = np.array([5,0])
C = np.array([5,1])
D = np.array([3,1])
E = np.array([2.5,1.6])
F = np.array([2,1])
G = np.array([0,1])

in_flow_marker  = 100
out_flow_marker = 200
wall_marker     = 300

gmsh.option.setNumber("General.Terminal",0) #To hide the mesh output values

mesh_size = 0.1
point_1 = gmsh.model.geo.add_point(A[0],A[1], 0, mesh_size)
point_2 = gmsh.model.geo.add_point(B[0],B[1], 0, mesh_size)
point_3 = gmsh.model.geo.add_point(C[0],C[1], 0, mesh_size)
point_4 = gmsh.model.geo.add_point(D[0],D[1], 0, mesh_size)
point_5 = gmsh.model.geo.add_point(E[0],E[1], 0, mesh_size)
point_6 = gmsh.model.geo.add_point(F[0],F[1], 0, mesh_size)
point_7 = gmsh.model.geo.add_point(G[0],G[1], 0, mesh_size)

line_1 = gmsh.model.geo.add_line(point_1, point_2)
line_2 = gmsh.model.geo.add_line(point_2, point_3)
line_3 = gmsh.model.geo.add_line(point_3, point_4)
line_4 = gmsh.model.geo.add_line(point_4, point_5)
line_5 = gmsh.model.geo.add_line(point_5, point_6)
line_6 = gmsh.model.geo.add_line(point_6, point_7)
line_7 = gmsh.model.geo.add_line(point_7, point_1)

curve_loop    = gmsh.model.geo.add_curve_loop([line_1,line_2,line_3,line_4,line_5,line_6,line_7])
plane_surface = gmsh.model.geo.add_plane_surface([curve_loop])


gmsh.model.addPhysicalGroup(2, [plane_surface], name = "fluid") # You need this for dolfinx

gmsh.model.addPhysicalGroup(1, [line_1,line_3,line_4,line_5,line_6], wall_marker)
gmsh.model.addPhysicalGroup(1, [line_7], in_flow_marker)
gmsh.model.addPhysicalGroup(1, [line_2], out_flow_marker)


mesh, cell_tags, facet_tags = model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0,gdim=2)


topology, cell_types, geometry_for_plotting = plot.create_vtk_mesh(mesh, 2)
#If you have the word geometry in place of geometry_for_plotting, it might conflict with the boundingBoxTree statement in the get coordinate function
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry_for_plotting)


plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)

if not pyvista.OFF_SCREEN:
    figure = plotter.screenshot("fundamentals_mesh.png")

t = 0
T = 10
num_steps = 500
dt = T/num_steps

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)

u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
fdim = mesh.topology.dim - 1

class inlet_pressure_class():
    def __init__(self, t):
        self.t = t
    def __call__(self, x):
        return 1 +0*x[0]
inlet_pressure = inlet_pressure_class(t)
u_inlet = Function(Q)
bc_inflow = dirichletbc(u_inlet, locate_dofs_topological(Q, fdim, facet_tags.find(in_flow_marker)))

u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bc_noslip = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, facet_tags.find(wall_marker)), V)

class outflow_pressure_class():
    def __init__(self, t):
        self.t = t
    def __call__(self, x):
        return 0 +0*x[0]
outlet_pressure = outflow_pressure_class(t)
u_outlet = Function(Q)
bc_outflow = dirichletbc(u_outlet, locate_dofs_topological(Q, fdim, facet_tags.find(out_flow_marker)))

bcu = [bc_noslip]
bcp = [bc_inflow, bc_outflow]

u_n = Function(V)
u_n.name = "u_n"
U = 0.5 * (u_n + u)
n = FacetNormal(mesh)
f = Constant(mesh, PETSc.ScalarType((0,0)))
k = Constant(mesh, PETSc.ScalarType(dt))
mu = Constant(mesh, PETSc.ScalarType(1))
rho = Constant(mesh, PETSc.ScalarType(1))

# Define strain-rate tensor
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(u.geometric_dimension())

# Define the variational problem for the first step
p_n = Function(Q)
p_n.name = "p_n"
F1 = rho*dot((u - u_n) / k, v)*dx
F1 += rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx
F1 += inner(sigma(U, p_n), epsilon(v))*dx
F1 += dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds
F1 -= dot(f, v)*dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))

A1 = assemble_matrix(a1, bcs=bcu)
b1 = create_vector(L1)

# Define variational problem for step 2
u_ = Function(V)
a2 = form(dot(nabla_grad(p), nabla_grad(q))*dx)
L2 = form(dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx)
A2 = assemble_matrix(a2, bcs=bcp)
b2 = create_vector(L2)

# Define variational problem for step 3
p_ = Function(Q)
a3 = form(dot(u, v)*dx)
L3 = form(dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx)
A3 = assemble_matrix(a3)
b3 = create_vector(L3)

# Solver for step 1
solver1 = PETSc.KSP().create(mesh.comm)
pc1 = solver1.getPC()

# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
pc2 = solver2.getPC()

# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
pc3 = solver3.getPC()

xdmf = XDMFFile(mesh.comm, "poiseuille.xdmf", "w")
xdmf.write_function(u_n, t)
xdmf.write_function(p_n, t)

for i in range(num_steps):
    # Update current time step
    t += dt

    # Step 1: Tentative veolcity step
    with b1.localForm() as loc_1:
    assemble_vector(b1, L1)
    apply_lifting(b1, [a1], [bcu])
    b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b1, bcu)
    solver1.solve(b1, u_.vector)
    # Step 2: Pressure corrrection step
    with b2.localForm() as loc_2:
    assemble_vector(b2, L2)
    apply_lifting(b2, [a2], [bcp])
    b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b2, bcp)
    solver2.solve(b2, p_.vector)

    # Step 3: Velocity correction step
    with b3.localForm() as loc_3:
    assemble_vector(b3, L3)
    b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    solver3.solve(b3, u_.vector)
    # Update variable with solution form this time step
    u_n.x.array[:] = u_.x.array[:]
    p_n.x.array[:] = p_.x.array[:]

    # Write solutions to file
    xdmf.write_function(u_n, t)
    xdmf.write_function(p_n, t)