Home > Numerical Analysis, Partial Differential Equations, Programming > Identifying edges and boundary points – 2D Mesh – Matlab

## Identifying edges and boundary points – 2D Mesh – Matlab

A triangulation algorithm often gives as output a list of points, and a list of triangle. Each triangle is represented by the indexes of the points which form it. Sometimes we need extra information on the structure of the triangulation, like the list of edges, or the list of boundary points. I will present below two fast algorithms for doing this.

Finding the list of edges is not hard. The idea is to go through each triangle, and extract all edges. The algorithm proposed below creates the adjacency matrix in the following way: for each triangle ${[i,j,k]}$ we set the elements ${a_{ij},a_{jk},a_{ik}}$ (and their symmetric correspondents) to be equal to ${1}$.

In order to find the points on the boundary (in two dimensions), it is enough to look for the edges which are sides to only one triangle. We can do this using the adjacency matrix. Note that if ${A}$ is the adjacency matrix, then ${A^2=(b_{ik})}$ stores the number of paths of length ${2}$ (two sides) between two points of the triangulation. Note that any edge which is not on the boundary will contain the starting and ending point of two paths of length ${2}$. If ${[i,j,k]}$ is a triangle such that points ${i,j}$ are on the boundary, then ${b_{i,j}=1}$ (there is one path of length ${2}$ going through ${i,k,j}$. We also have ${a_{i,j} = 1}$. Conversely, if ${a_{i,j} = 1}$ and ${b_{i,j} = 1}$ then ${i,j}$ are connected, and there is a unique path of length ${2}$ going from ${i}$ to ${j}$. Thus, ${i,j}$ is an edge on the boundary. Therefore, we just need to identify the indexes ${i}$ such that there exists ${j}$ with ${a_{i,j} b_{i,j} = 1}$.

Below are two short Matlab codes doing these two algorithms. I guess they are close to being optimal, since only sparse and vectorized operations are used.

```
%p is the list of points
%T is the list of triangles, ap is the number of points
%this computes the adjacency matrix
A = min(sparse(T(:,1),T(:,2),1,ap,ap)+sparse(T(:,2),T(:,3),1,ap,ap)+sparse(T(:,3),T(:,1),1,ap,ap),1);
A = min(A+A',1);
% this finds the boundary points, whose indexes are stored in Ibord
B = A^2.*A==1;
Ibord = find(sum(B,2)>0);
```
Advertisements
1. November 2, 2015 at 8:23 am

Quite precise and fast. Good job.

1. No trackbacks yet.