Home > Uncategorized > Simple triangle mesh – Matlab code

## Simple triangle mesh – Matlab code

There are cases when you need to define a simple, regular mesh on a nice set, like a triangle, or a polygon. There are software which are capable of creating such meshes, but most of them will require some computation time (Delaunay triangulation, plus optimization, edge swapping, etc.). I propose to you a simple algorithm for defining a mesh on a triangle. Of course, this can be extended to polygons or other cases (I used a variant of it to create a mesh of a spherical surface).

The idea is very simple: start with three points ${P_1,P_2,P_3}$, and define the (only) initial triangle as ${T = [1,\ 2,\ 3]}$. Then, at each step, take a triangle ${t_i}$ from the list ${T}$, add the midpoints of ${t_i}$ to the list of points and replace the triangle ${t_i}$ with the four smaller triangles determined by the vertices and midpoints of ${t_i}$. This is a basic mesh refining strategy. Thus, in short, we start with a triangle and do a number of mesh refinements, until we reach the desired side-length. Of course, at each step we should be careful not to add the midpoint of an edge two times (once for every neighbouring triangle), but in what follows, I will always remove duplicates by bruteforce in the end (not the most economic solution).

I present you below a Matlab code. It is not optimal, so don’t try a great number of refinements. It should work relatively well for ${7}$ or ${8}$ refinements. For improved performance, some vectorization, plus avoiding duplicates strategy should be made. If you want to generalize it to a quadrilateral, or a polygon, simply start from a triangulation of the polygon: give the initial points, and the indexes of the triangulation as starting parameters.

function [p,T] = simple_triangle_mesh(n,p)
% input arguments
% the number of refinements n < 8 (larger may be very slow)
% p the coordinates of the points as a 3x2 matrix
% outputs
% the vertices of the triangles p: a m x 2 matrix
% the triangle matrix: each line represents the indexes of the respective triangle
if nargin < 2
p = [-1 0;
1 0;
0 sqrt(3)];
end

T = [1 2 3];

ap = size(p,1);
for i = 1:n
ct = size(T,1);
T1 = [];
E = [];
for j = 1:ct;
Tr = p(T(j,:),:);
p1 = Tr(1,:);
p2 = Tr(2,:);
p3 = Tr(3,:);

n1 = (p2+p3)/2;
n2 = (p1+p3)/2;
n3 = (p1+p2)/2;

t1 = [ap+1 ap+2 ap+3];
t2 = [T(j,1) ap+2 ap+3];
t3 = [T(j,2) ap+1 ap+3];
t4 = [T(j,3) ap+1 ap+2];
T1 = [T1; t1;t2;t3;t4];
p = [p; n1;n2;n3];

ap = ap+3;
end
T = T1;
dx = norm(n1-n2);
[pt,I,J] = unique(p,'rows');
p = p(I,:);
ap = size(p,1);
T(:,1:3) = J(T(:,1:3));
end ;

%plotting the points
plot(p(:,1),p(:,2),'.');
axis equal

Here is the result for 5 refinements:

Here’s a picture I obtained doing an optimization problem on a mesh obtained with this algorithm. Guess what is represented here?