Extract coordinates at the mesh boundary and evaluate the function value

General
1

I have a 2d mesh which can be in any shape and I want to extract the coordinate of the mesh vertices (of triangles) at the boundary, then evaluate my function at these points. I can use mesh.geometry.x to print out coordinates of all vertices but I don’t know how to print those only at the boundary.

After having the coordinates, the evaluation can be done following the tutorial. Since I only need the value at each corresponding vertices, I wonder if there is a shorter way of doing this.

Eventually I want to have a list with elements list[i]=[[x,y],[fx(x,y),fy(x,y)]] where [x,y] is the coordinate at the boundary and fx(x,y), fy(x,y) are the vector value evaluated there. My goal here is to use the value (fx, fy) as the boundary coordinate of my new mesh.

Here is a related post unfortunately I couldn’t understand…

Here is an example of preset. The mesh here is bound by a circle but the answer should not rely on the special shape.

import basix.ufl
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import dolfinx.nls.petsc
import dolfinx.log
import dolfinx.plot
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import ufl
import viskex
import matplotlib.pyplot as plt
from petsc4py.PETSc import ScalarType
import pyvista
import sys

abar = 6
beta = 2.
q = 0.
alpha = 1.0 - q/6.
nu_p = 0.8
tau = 0.
R_0 = 10.
Y = 1.
lcar = .27
# circular mesh

gmsh.initialize()
gmsh.model.add("mesh")

p0 = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, lcar)
p1 = gmsh.model.geo.addPoint(0.0, abar, 0.0, lcar)
p2 = gmsh.model.geo.addPoint(0.0, -abar, 0.0, lcar)
c0 = gmsh.model.geo.addCircleArc(p1, p0, p2)
c1 = gmsh.model.geo.addCircleArc(p2, p0, p1)
boundary = gmsh.model.geo.addCurveLoop([c0, c1])

domain = gmsh.model.geo.addPlaneSurface([boundary])

p_center = gmsh.model.geo.addPoint(0, 0, 0, lcar)

gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [c0, c1], 1)
gmsh.model.addPhysicalGroup(0, [p_center], 2)
gmsh.model.addPhysicalGroup(2, [boundary], 0)

gmsh.model.mesh.embed(0, [p_center], 2, domain)

gmsh.model.mesh.generate(2)

mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)
gmsh.finalize()

boundary_facet = boundaries.indices[boundaries.values == 1]
bb_tree = dolfinx.geometry.bb_tree(mesh, mesh.topology.dim)

V = dolfinx.fem.VectorFunctionSpace(mesh, ("Lagrange", 2))

vxvy = ufl.TestFunction(V)
(vx, vy) = vxvy[0], vxvy[1]
duxuy = ufl.TrialFunction(V)
uxuy = dolfinx.fem.Function(V)
(ux, uy) = uxuy[0], uxuy[1]

uxuy.interpolate(lambda x: [x[0]**2 + 1, x[1]])
2

One way to do this is below. The issue is when evaluating the value, if the condition if len(colliding_cells.links(i)) > 0 is not satisfied, the coordinate will be dropped.

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_vertices = dolfinx.mesh.compute_incident_entities(mesh.topology, boundary_facets, mesh.topology.dim-1, 0)

boundary_coordinates = []
for ind in boundary_vertices:
    boundary_coordinates.append(mesh.geometry.x[ind].tolist())

points = np.array(boundary_coordinates).transpose()

cells = []
points_on_proc = []
cell_candidates = dolfinx.geometry.compute_collisions_points(bb_tree, points.T)
colliding_cells = dolfinx.geometry.compute_colliding_cells(mesh, cell_candidates, points.T)
for i, point in enumerate(points.T):
    if len(colliding_cells.links(i)) > 0:
        points_on_proc.append(point)
        cells.append(colliding_cells.links(i)[0])

points_on_proc = np.array(points_on_proc, dtype=np.float64)
actual_boundary_coordinates = uxuy.eval(points_on_proc, cells)

plt.plot(actual_boundary_coordinates[:,0], actual_boundary_coordinates[:,1], '.', color='blue')
plt.plot(points[0], points[1], '.', color='red')
plt.gca().set_aspect('equal')