For a 2D unstructured mesh and a parametric curve (here, consider a uniform quadratic B-spline), I wish to obtain the x,y coordinates of all intersection points between the mesh cell edges and the curve.
The form of the curve is \mathbf{r}(\hat{u}_i) = [\mathbf{p}_i, \mathbf{p}_{i+1}, \mathbf{p}_{i+2}]\mathbf{C}_i\array{[\hat{u}_i^2 &\hat{u}_i&1]^T}, where each \mathbf{p} is a control point and \hat{u}_i\in[0,1] is the parameter for each Bezier curve segment.
Two assumptions: 1) the curve passes through the first and last control points, and 2) these two control points live on the mesh boundary. I would guess a good approach is to start with the first control point, find the host cell, and then do a breadth-first search of neighbors until I reach the last control point. However, I haven’t been able to figure out how to implement this beyond the 1st cell… MWE follows - apologies it is rather lengthy but it is mostly plotting code
Performance is important as I need to do this with multiple curves on meshes as part of an iterative optimization procedure. Probably most cells will not be intersected by any one curve, so is there a way to only consider cells that intersect the bbox of the control polygon? Still using fenics 2019 because I need dolfin adjoint.
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
# Plot a mesh
fig, ax = plt.subplots(figsize=[10,10])
mesh = UnitSquareMesh(4, 4)
coords = mesh.coordinates()
ax.triplot(coords[:,0], coords[:,1], triangles=mesh.cells(), linewidth=0.5, color='k')
ax.axis('equal')
# label cells
for cell_no in range(mesh.num_cells()):
centroid = Cell(mesh, cell_no).midpoint()
ax.text(*centroid.array()[:2], str(cell_no))
# consider quadratic uniform bspline with natural knot vector {0,1,2...}
# some bspline control points, three bezier curve segments
p0 = [0,0.45]
p1 = [0.225, 0.9]
p2 = [0.425, 0.15]
p3 = [0.625, 0.125]
p4 = [1.0, 0.35]
p = np.array([p0, p1, p2, p3, p4])
# plot control points
ax.plot(list(p[:,0]), list(p[:,1]), '--', lw=0.5, marker='o', color='r')
# BSpline coefficients
# 1st segment
C0 = np.matrix([[1,-2, 1],
[-1.5, 2, 0],
[0.5, 0, 0]])
# middle segement
Ci = np.matrix([[0.5,-1, 0.5],
[-1, 1, 0.5],
[0.5, 0, 0]])
# last segment
Cm = np.matrix([[0.5,-1, 0.5],
[-1.5, 1, 0.5],
[1, 0, 0]])
def u(u):
return np.matrix([u**2, u, 1])
# plot the bspline for visual reference
uu = np.linspace(0, 1, 50)
U = np.array([uu**2, uu, np.ones(50)])
r = [np.array(p[0:3]).T*C0*U,
np.array(p[1:4]).T*Ci*U,
np.array(p[2:]).T*Cm*U]
for ri in r:
ax.plot((np.asarray(ri[0])).flatten(), (np.asarray(ri[1])).flatten(), lw=1)
x_t1 = []
y_t1 = []
# only consider first curve for simplicity
# 1st bezier curve coefficients
C1 = np.asarray(np.array(p[0:3]).T*C0)
C1x = np.asarray(C1[0]) # 1st curve coefficients, x = C1x(u)
C1y = np.asarray(C1[1]) # 1st curve coefficients, y = C1y(u)
# Start with first control point
tree = mesh.bounding_box_tree()
cell_index = tree.compute_first_entity_collision(Point(*p0))
# iterate over cell edges and test for intersection
for edge in edges(Cell(mesh, cell_index)):
pt1, pt2 = mesh.coordinates()[list(edge.entities(0))] # get edge nodal coords
xx = [pt1[0], pt2[0]] # the two nodal x coords
yy = [pt1[1], pt2[1]] # the two nodal y coords
xx.sort()
yy.sort()
if abs(pt2[0] - pt1[0]) <= 1e-6: # case of vertical edge, x = c
x_intercept = pt2[0]
roots = np.roots(C1x - np.array([0, 0, x_intercept]).T)
for root in roots:
xy = np.asarray(np.array(p[0:3]).T*C0*u(root).T) # cartesian coords of found root
if (root>=0 and root <=1) and (xy[1] >= yy[0] and xy[1]<= yy[1]): # check root real and if point on edge
y_t1.append(root)
else: # cell edge is line of form y = mx+b
m = (pt2[1] - pt1[1])/(pt2[0]-pt1[0])
b = -pt1[0]*m + pt1[1]
# plug in parametric x, y coords in to y=mx + b and solve
roots = np.roots(C1y - m*C1x - np.array([0, 0, b]))
for root in roots:
if root >= 0 and root <=1: # check if root is real
# check if point on edge
xy = np.asarray(np.array(p[0:3]).T*C0*u(root).T) # cartesian coords of found root
if (xy[0] >= xx[0] and xy[0]<= xx[1]) and (xy[1] >= yy[1] and xy[1]<= yy[1]):
x_t1.append(root)
for root in x_t1:
pt = np.asarray(np.array(p[0:3]).T*C0*u(root).T)
ax.plot(pt[0], pt[1], marker='s', color='g', fillstyle='none')
for root in y_t1:
pt = np.asarray(np.array(p[0:3]).T*C0*u(root).T)
ax.plot(pt[0], pt[1], marker='s', color='k', fillstyle='none')