/
triangulate.m
84 lines (74 loc) · 2 KB
/
triangulate.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
% triangulate - interface to 'triangle' library.
%
% [vertex, face] = compute_constrained_delaunay(vertex, edges, holes);
%
% Triangulates a set of points with boundary, using a constrainted Delaunay
% triangulation. It uses the Triangle library of Jonathan Richard Shewchuk:
% http://www-2.cs.cmu.edu/~quake/triangle.html
%
% Copyright (c) 2004 Gabriel Peyré
function [nodes, triangles] = triangulate(nodes, boundary, holes)
if size(nodes,1)<size(nodes,2)
nodes = nodes';
end
if size(boundary,1)<size(boundary,2)
boundary = boundary';
end
if size(holes,1)<size(holes,2)
holes = holes';
end
nodes = [(1:size(nodes,1))', nodes];
boundary = [(1:size(boundary,1))', boundary];
holes = [(1:size(holes,1))', holes];
n = size(nodes, 1);
m = size(boundary, 1);
h = size(holes, 1);
% Write out domain in PSLG format.
fid = fopen('mesh.poly', 'w');
if fid<0
error('Unable to create file mesh.poly.');
end
fprintf(fid, '%d 2 0 0\n', n);
for i = 1:n
fprintf(fid, '%d %f %f\n', nodes(i,1), nodes(i,2), nodes(i,3));
end
fprintf(fid, '%d 0\n', m);
for j = 1:m
fprintf(fid, '%d %d %d\n', boundary(j,1), ...
boundary(j,2), ...
boundary(j,3));
end
fprintf(fid, '%d\n', h);
for k = 1:h
fprintf(fid, '%f %f %f\n', holes(k,1), holes(k,2), holes(k,3));
end
fclose(fid);
% Triangulate domain.
!triangle -pq mesh.poly > /dev/tmp
% Read in new nodes.
fid = fopen('mesh.1.node', 'r');
if fid<0
error('No file mesh.1.node produced.');
end
n = fscanf(fid, '%d', 1);
dummy = fscanf(fid, '%d', 3);
nodes = [];
for i = 1:n
point = [ fscanf(fid, '%d', 1)' fscanf(fid, '%f', 2)' ];
nodes = [ nodes; point ];
fscanf(fid, '%d', 1);
end
fclose(fid);
% Read in triangle mesh.
fid = fopen('mesh.1.ele', 'r');
m = fscanf(fid, '%d', 1);
dummy = fscanf(fid, '%d', 2);
triangles = [];
for j = 1:m
triangles = [ triangles; fscanf(fid, '%d', 4)' ];
end
fclose(fid);
triangles = triangles(:,2:end)';
nodes = nodes(:,2:end)';
% Delete temporary files.
!rm -f mesh.poly mesh.1.poly mesh.1.node mesh.1.ele