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
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
|
function savemsh(node, elem, fname, rname)
%
% savemsh(node,elem,fname,rname)
%
% save a tetrahedral mesh to GMSH mesh format
%
% author: Riccardo Scorretti (riccardo.scorretti<at> univ-lyon1.fr)
% date: 2013/07/22
%
% input:
% node: input, node list, dimension (nn,3)
% elem: input, tetrahedral mesh element list, dimension (ne,4) or (ne,5) for multi-region meshes
% fname: output file name
% rname: name of the regions, cell-array of strings (optional)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%
if nargin < 4
rname = {};
end
if size(elem, 2) < 5
elem(:, 5) = 1;
end
fid = fopen(fname, 'wt');
if (fid == -1)
error('You do not have permission to save mesh files.');
end
% Check that all the elements are correctly oriented
elem(:, 1:4) = meshreorient(node, elem(:, 1:4));
nbNodes = size (node, 1);
reg = unique (elem(:, 5));
reg(reg <= 0) = max(reg) + 1 - reg(reg <= 0); % convert label 0 to max(reg)+1, -1 to max(reg)+2, and so on
nbRegion = length (reg);
nbElements = size (elem, 1);
% Create the skeleton of the mesh structure
M.Info.version = [];
M.Nodes.nb = 0;
M.Nodes.x = [];
M.Nodes.y = [];
M.Nodes.z = [];
M.Elements.nb = 0;
M.Elements.type = zeros (0, 0, 'uint8');
M.Elements.tableOfNodes = zeros (0, 0, 'uint32');
M.Elements.region = zeros (0, 0, 'uint16');
M.Regions.nb = 0;
M.Regions.name = {};
M.Regions.dimension = [];
% Build the table of nodes
M.Nodes.nb = nbNodes;
M.Nodes.x = node(:, 1);
M.Nodes.y = node(:, 2);
M.Nodes.z = node(:, 3);
clear node;
% Build the table of elements
M.Elements.nb = nbElements;
M.Elements.type = uint8(4 * ones(nbElements, 1));
M.Elements.tableOfNodes = uint32(elem(:, 1:4));
M.Elements.region = uint16(elem(:, 5));
clear elem;
% Build the table of regions
M.Regions.nb = max(reg);
for k = 1:nbRegion
if length(rname) < k
rname{k} = sprintf('region_%d', k);
end
M.Regions.name{reg(k)} = sprintf ('%s', rname{k});
M.Regions.dimension(reg(k)) = 3;
end
% Writhe the header
fprintf (fid, '$MeshFormat\n2.2 0 8\n$EndMeshFormat\n');
% Write the physical names
if M.Regions.nb > 0
fprintf (fid, '$PhysicalNames\n');
fprintf (fid, '%d\n', M.Regions.nb);
for r = 1:M.Regions.nb
name = M.Regions.name{r};
if isempty (name)
name = sprintf ('Region_%d', r);
end
fprintf (fid, '%d %d "%s"\n', M.Regions.dimension(r), r, name);
end
fprintf (fid, '$EndPhysicalNames\n');
end
% Write the nodes
fprintf (fid, '$Nodes\n');
fprintf (fid, '%d\n', size(M.Nodes.x, 1));
buffer = [1:M.Nodes.nb; M.Nodes.x'; M.Nodes.y'; M.Nodes.z'];
fprintf (fid, '%d %10.10f %10.10f %10.10f\n', buffer);
fprintf (fid, '$EndNodes\n');
% Write the elements
%
% In order to accelerate the printing, the elements are printed in groups of (blockSize) elements, and
% are grouped by (homogeneous) type. This variable sets the size of each group.
%
blockSize = 100000;
fprintf (fid, '$Elements\n');
fprintf (fid, '%d\n', M.Elements.nb);
for h = 1:blockSize:ceil(length(M.Elements.type) / blockSize) * blockSize
e = h:min(length(M.Elements.type), h + blockSize - 1); % = elements being considered
type = unique (M.Elements.type(e)); % = types of elements found in this group
%
% Process each type of element separately
%
for k = 1:length(type)
if type(k) == 0
continue
end
et = e(find(M.Elements.type(e) == type(k))); % = elements of the group of the same type
%
% Determine the format for printing the elements
%
elementFormat = '%d %d %d %d %d %d\n';
for n = 1:4
elementFormat = [elementFormat '%d '];
end
elementFormat = [elementFormat '\n'];
%
% Collect in a buffer all the data of the elements of index (et)
%
buffer = zeros (10, length(et));
buffer(1, :) = et;
buffer(2, :) = type(k);
buffer(3, :) = 3;
buffer(4, :) = M.Elements.region(et);
buffer(5, :) = M.Elements.region(et);
buffer(6, :) = 0;
for n = 1:4
buffer(6 + n, :) = M.Elements.tableOfNodes(et, n);
end
%
% Print all the homogeneous elements in the group with a single instruction
%
fprintf (fid, elementFormat, buffer);
end
end
fprintf (fid, '$EndElements\n');
fclose(fid);
|
