I have a list of vertical cut planes, and I have a polygonal mesh ( it's a 2D+0.5D mesh, something like a 2D mesh with an extra dimension, $z$ attached to each point). One can assume that the mesh contains vertices $V$ , edges $E$ and face $F$.The mesh is a real life 3D terrain.Now, given a list of vertical cut planes, I want to find all the intersection points ( to form a polyline) between the cut planes and the polygonal mesh, in order to build a vertical view of the terrain for each of the cut plane.Currently, my approach is that I use the Aabb tree to test each of the mesh triangles against the cut section, and then compute the intersection points subsequently for the cut section and the triangle that are within each other's Aabb tree range. But this approach doesn't make use of the mesh structure of the terrain, and I think we can significantly speed up the computation if we can make use of the mesh structure.Is there a more efficient algorithm I can use for this purpose? 
How to efficiently determine the intersection of a vertical cutting plane with a mesh
computational geometry;mesh traversal
You have an 2D mesh of height values. I'm going to assume you create the faces of the surface for the $i$,$j$th square using two triangles of $\left\langle\left(i,j\right),\left(i+1,j\right),\left(i+1,j+1\right)\right\rangle$ and $\left\langle\left(i,j\right),\left(i,j+1\right),\left(i+1,j+1\right)\right\rangle$.Your vertical cut plain can be defined as a 2D line called $l$. $l$ can be represented by using 2 known points on the edge of the mesh or by using an 2D line equation $i = mj+c$, $l=\left\langle{m,c}\right\rangle$.The intersection points are the intersection with any of the edges of triangles along that 2D line $l$.All the points heights are worked out by using linear interpolation between the heights given by the two array coordinated connected by the triangle edge being intersected. Distance along the cut plane is given by Pythagorean theorem.#!/usr/bin/env pythonimport mathfrom random import choice, uniformmesh = [[10,5,8],[12,7,5],[4,7,3],[4,6,4],[8,3,1]]print meshiMax = len(mesh)-1jMax = len(mesh[0])-1# pick random vertical cut planeif choice([True,False]):  p1 = (uniform(0,iMax), choice([0,jMax]))else:  p1 = (choice([0,iMax]), uniform(0,jMax))while True:  if choice([True,False]):    p2 = (uniform(0,iMax), choice([0,jMax]))  else:    p2 = (choice([0,iMax]), uniform(0,jMax))  if p1[0] != p2[0] and  p1[1] != p2[1]:    breakprint P1,p1print P2,p2if p1[1] == p2[1]:  m = float('Inf') # Not handling this, needs fix  c = Noneelse:  m = float(p2[0]-p1[0])/(p2[1]-p1[1])  c = p1[0]-m*p1[1]# i=mj+c formularprint i={0}j+{1}.format(m,c)def floor(a):  return int(math.floor(a))def ceil(a):  return int(math.ceil(a))iL = floor(min(p1[0], p2[0]))iU = ceil(max(p1[0], p2[0]))jL = floor(min(p1[1], p2[1]))jU = ceil(max(p1[1], p2[1]))print points to plot unordered, d is distance from p1 (pythagoras), h is height as given by linear interpolationdef linearInterpolation(h0,h1,r):  return (h1-h0)*r+h0def pythag(a,b):  return math.sqrt(a**2 + b**2)for i in xrange(iL, iU+1):  j = (i-c)/m  if 0 <= j <= jMax:    j1, j2 = floor(j), ceil(j)    h = linearInterpolation(mesh[i][j1], mesh[i][j2], j - j1)    d = pythag(p1[0]-i, p1[1]-j)    print (d,h)for j in xrange(jL, jU+1):  i = m*j+c  if 0 <= i <= iMax:    i1, i2 = floor(i), ceil(i)    h = linearInterpolation(mesh[i1][j], mesh[i2][j], i - i1)    d = pythag(p1[0]-i, p1[1]-j)    print (d,h)for v in xrange(-iU-1, jU):  i = (c - v)/(1.0-m)  j = m*i+c  if 0 <= i <= iMax and 0 <= j <= jMax:    i1, i2 = floor(i), ceil(i)    j1, j2 = floor(j), ceil(j)    r = pythag(i - i1, j - j1) / math.sqrt(2)    h = linearInterpolation(mesh[i1][j1], mesh[i2][j2], r)    d = pythag(p1[0]-i, p1[1]-j)    print (d,h)