## 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 , and define the (only) initial triangle as . Then, at each step, take a triangle from the list , add the midpoints of to the list of points and replace the triangle with the four smaller triangles determined by the vertices and midpoints of . 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 or 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?