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)];

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;
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
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?


Categories: Uncategorized Tags: , ,
  1. No comments yet.
  1. No trackbacks yet.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: