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
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
|
function [newnode, newelem, newelem0] = surfboolean(node, elem, varargin)
%
% [newnode,newelem,newelem0]=surfboolean(node1,elem1,op2,node2,elem2,op3,node3,elem3,...)
%
% merge two or more triangular meshes and resolve intersecting elements
%
% author: Qianqian Fang <q.fang at neu.edu>
%
% input:
% node: node coordinates, dimension (nn,3)
% elem: triangle surfaces (ne,3)
% op: a string of a boolean operator, possible op values include
% 'union' or 'or': the outter surface of the union of the enclosed space
% 'inter' or 'and': the surface of the domain contained by both meshes
% 'diff' or '-': the surface of the domain in mesh 1 excluding that of
% mesh 2
% 'all', 'resolve', or 'xor' or '+': the output contains 4 subsurfaces, identified by the 4th
% column of newelem:
% 1: mesh 1 outside of mesh 2
% 2: mesh 2 outside of mesh 1
% 3: mesh 1 inside of mesh 2
% 4: mesh 2 inside of mesh 1
% you can use newelem(find(mod(newelem(:,4),2)==1),:) to
% get mesh 1 cut by mesh 2, or newelem(find(mod(newelem(:,4),2)==0),:)
% to get mesh 2 cut by mesh 1;
% 'first': combine 1 and 3 from the output of 'all'
% 'second': combine 2 and 4 from the output of 'all'
% 'self': test for self-intersections; only the first mesh is
% tested; other inputs are ignored.
% 'remesh': remesh the first input mesh (2nd input mesh is ignored)
% 'decouple': separate two shells and make sure there is no intersection;
% the input surfaces must be closed and ordered from outer to inner
%
% if a dash '-' is preceding a command, including
% -diff, -union, -isct, -resolve, -first, -second, -xor
% a 'remesh' operation is applied to the output mesh to remove poorly shaped
% triangles
%
% output:
% newnode: the node coordinates after boolean operations, dimension (nn,3)
% newelem: tetrahedral element or surfaces after boolean operations (nn,4) or (nhn,5)
% newelem0: when the operator is 'self', return the intersecting
% element list in terms of the input node list (experimental)
%
% example:
%
% [node1,face1,elem1]=meshabox([0 0 0],[10 10 10],1,1);
% [node2,face2,elem2]=meshabox([0 0 0]+5,[10 10 10]+5,1,1);
% [newnode,newface]=surfboolean(node1,face1,'union',node2,face2);
% plotmesh(newnode,newface);
% [newnode,newface]=surfboolean(node1,face1,'-union',node2,face2);
% figure; plotmesh(newnode,newface);
% [newnode,newface]=surfboolean(node1,face1,'diff',node2,face2);
% figure;plotmesh(newnode,newface,'x>5');
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%
len = length(varargin);
newnode = node;
newelem = elem;
if (len > 0 && mod(len, 3) ~= 0)
error('you must give operator, node and element in a triplet form');
end
try
exename = evalin('caller', 'ISO2MESH_SURFBOOLEAN');
catch
try
exename = evalin('base', 'ISO2MESH_SURFBOOLEAN');
catch
exename = 'cork';
end
end
exesuff = fallbackexeext(getexeext, exename);
randseed = hex2dec('623F9A9E'); % "U+623F U+9A9E"
if (~isempty(getvarfrom({'caller', 'base'}, 'ISO2MESH_RANDSEED')))
randseed = getvarfrom({'caller', 'base'}, 'ISO2MESH_RANDSEED');
end
for i = 1:3:len
op = varargin{i};
no = varargin{i + 1};
el = varargin{i + 2};
opstr = op;
if (strcmp(op, 'or'))
opstr = 'union';
end
if (strcmp(op, 'xor'))
opstr = 'all';
end
if (strcmp(op, 'and'))
opstr = 'isct';
end
if (strcmp(op, '-'))
opstr = 'diff';
end
if (strcmp(op, 'self'))
opstr = 'solid';
end
if (strcmp(opstr, 'all'))
opstr = 'resolve';
end
tempsuff = 'off';
deletemeshfile(mwpath(['pre_surfbool*.' tempsuff]));
deletemeshfile(mwpath('post_surfbool.off'));
if (strcmp(opstr, 'all'))
deletemeshfile(mwpath('s1out2.off'));
deletemeshfile(mwpath('s1in2.off'));
deletemeshfile(mwpath('s2out1.off'));
deletemeshfile(mwpath('s2in1.off'));
end
if (strcmp(op, 'decouple'))
if (exist('node1', 'var') == 0)
node1 = node;
elem1 = elem;
newnode(:, 4) = 1;
newelem(:, 4) = 1;
end
opstr = ' --decouple-inin 1 --shells 2'; % -q
saveoff(node1(:, 1:3), elem1(:, 1:3), mwpath('pre_decouple1.off'));
if (size(no, 2) ~= 3)
opstr = ['-q --shells ' num2str(no)];
cmd = sprintf('cd "%s" && "%s%s" "%s" %s', mwpath, mcpath('meshfix'), exesuff, ...
mwpath('pre_decouple1.off'), opstr);
else
saveoff(no(:, 1:3), el(:, 1:3), mwpath('pre_decouple2.off'));
cmd = sprintf('cd "%s" && "%s%s" "%s" "%s" %s', mwpath, mcpath('meshfix'), exesuff, ...
mwpath('pre_decouple1.off'), mwpath('pre_decouple2.off'), opstr);
end
elseif (strcmp(op, 'decoupleout'))
if (exist('node1', 'var') == 0)
node1 = node;
elem1 = elem;
newnode(:, 4) = 1;
newelem(:, 4) = 1;
end
opstr = [' --decouple-outout 1 --shells 2']; % -q
saveoff(node1(:, 1:3), elem1(:, 1:3), mwpath('pre_decouple1.off'));
if (size(no, 2) ~= 3)
opstr = ['-q --shells ' num2str(no)];
cmd = sprintf('cd "%s" && "%s%s" "%s" %s', mwpath, mcpath('meshfix'), exesuff, ...
mwpath('pre_decouple1.off'), opstr);
else
saveoff(no(:, 1:3), el(:, 1:3), mwpath('pre_decouple2.off'));
cmd = sprintf('cd "%s" && "%s%s" "%s" "%s" %s', mwpath, mcpath('meshfix'), exesuff, ...
mwpath('pre_decouple1.off'), mwpath('pre_decouple2.off'), opstr);
end
elseif (strcmp(op, 'separate'))
if (exist('node1', 'var') == 0)
node1 = node;
elem1 = elem;
newnode(:, 4) = 1;
newelem(:, 4) = 1;
end
opstr = ' --shells 2'; % -q
saveoff(node1(:, 1:3), elem1(:, 1:3), mwpath('pre_decouple1.off'));
if (size(no, 2) ~= 3)
opstr = ['-q --shells ' num2str(no)];
cmd = sprintf('cd "%s" && "%s%s" "%s" %s', mwpath, mcpath('meshfix'), exesuff, ...
mwpath('pre_decouple1.off'), opstr);
else
saveoff(no(:, 1:3), el(:, 1:3), mwpath('pre_decouple2.off'));
cmd = sprintf('cd "%s" && "%s%s" "%s" "%s" %s', mwpath, mcpath('meshfix'), exesuff, ...
mwpath('pre_decouple1.off'), mwpath('pre_decouple2.off'), opstr);
end
elseif (strcmp(op, 'remesh'))
saveoff(newnode(:, 1:3), newelem(:, 1:3), mwpath(['pre_surfbool1.' tempsuff]));
cmd = sprintf('cd "%s" && "%s%s" %s%s "%s" "%s" -%d', mwpath, mcpath(exename), exesuff, '-', ...
opstr, mwpath(['pre_surfbool1.' tempsuff]), mwpath('post_surfbool.off'), randseed);
else
saveoff(newnode(:, 1:3), newelem(:, 1:3), mwpath(['pre_surfbool1.' tempsuff]));
saveoff(no(:, 1:3), el(:, 1:3), mwpath(['pre_surfbool2.' tempsuff]));
cmd = sprintf('cd "%s" && "%s%s" %s%s "%s" "%s" "%s" -%d', mwpath, mcpath(exename), exesuff, '-', ...
opstr, mwpath(['pre_surfbool1.' tempsuff]), mwpath(['pre_surfbool2.' tempsuff]), mwpath('post_surfbool.off'), randseed);
end
[status, outstr] = system(cmd);
if (status ~= 0 && strcmp(op, 'self') == 0)
error('surface boolean command failed:\n%s\nERROR: %s\n', cmd, outstr);
end
if (strcmp(op, 'self'))
if (isempty(strfind(outstr, 'NOT SOLID')))
fprintf(1, 'no self-intersection was found!\n');
newnode = [];
newelem = [];
newelem0 = [];
else
fprintf(1, 'input mesh is self-intersecting\n');
newnode = 1;
newelem = [];
newelem0 = 1;
end
return
end
if (strcmp(opstr, 'all'))
% tag the 4 piceses of meshes, this tag do not propagate to the next boolean operation
[nnode, nelem] = readoff(mwpath('s1out2.off'));
newelem = [nelem ones(size(nelem, 1), 1)];
newnode = [nnode ones(size(nnode, 1), 1)];
[nnode, nelem] = readoff(mwpath('s1in2.off'));
newelem = [newelem; nelem + size(newnode, 1) 3 * ones(size(nelem, 1), 1)];
newnode = [newnode; nnode 3 * ones(size(nnode, 1), 1)];
[nnode, nelem] = readoff(mwpath('s2out1.off'));
newelem = [newelem; nelem + size(newnode, 1) 2 * ones(size(nelem, 1), 1)];
newnode = [newnode; nnode 2 * ones(size(nnode, 1), 1)];
[nnode, nelem] = readoff(mwpath('s2in1.off'));
newelem = [newelem; nelem + size(newnode, 1) 4 * ones(size(nelem, 1), 1)];
newnode = [newnode; nnode 4 * ones(size(nnode, 1), 1)];
elseif (strcmp(op, 'decouple'))
[newnode, newelem] = readoff(mwpath('pre_decouple1_fixed.off')); % [node1,elem1]
% newelem=[newelem;elem1+size(newnode,1) (i+1)*ones(size(elem1,1),1)];
% newnode=[newnode;node1 (i+1)*ones(size(node1,1),1)];
elseif (strcmp(op, 'separate'))
[newnode, newelem] = readoff(mwpath('pre_decouple1_fixed.off'));
elseif (strcmp(op, 'decoupleout'))
[newnode, newelem] = readoff(mwpath('pre_decouple1_fixed.off')); % [node1,elem1]
else
[newnode, newelem] = readoff(mwpath('post_surfbool.off'));
if (strcmp(op, 'self'))
fprintf(1, 'a total of %d self-intersecting elements were found\n', size(newelem, 1));
if (nargout >= 3)
[found, newelem0] = ismember(newnode, node, 'rows');
if (~all(found))
error('self intersecting elements contain new nodes');
end
newelem0 = newelem0(newelem);
end
return
end
end
end
|
