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
|
function [node, elem, face] = cgals2m(v, f, opt, maxvol, varargin)
%
% [node,elem,face]=cgals2m(v,f,opt,maxvol)
%
% wrapper for CGAL 3D mesher (CGAL 3.5 and newer)
% convert a triangular surface to tetrahedral mesh
%
% http://www.cgal.org/Manual/3.5/doc_html/cgal_manual/Mesh_3/Chapter_main.html
%
% author: Qianqian Fang (q.fang at neu.edu)
%
% input:
% v: the node coordinate list of a surface mesh (nn x 3)
% f: the face element list of a surface mesh (be x 3)
% opt: parameters for CGAL mesher, if opt is a structure, then
% opt.radbound: defines the maximum surface element size
% opt.angbound: defines the miminum angle of a surface triangle
% opt.distbound: defines the maximum distance between the
% center of the surface bounding circle and center of the
% element bounding sphere
% opt.reratio: maximum radius-edge ratio
% if opt is a scalar, it only specifies radbound.
% maxvol: target maximum tetrahedral elem volume
%
% output:
% node: output, node coordinates of the tetrahedral mesh
% elem: output, element list of the tetrahedral mesh, the last
% column is the region id
% face: output, mesh surface element list of the tetrahedral mesh
% the last column denotes the boundary ID
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%
fprintf(1, 'creating surface and tetrahedral mesh from a polyhedral surface ...\n');
exesuff = fallbackexeext(getexeext, 'cgalpoly');
ang = 30;
ssize = 6;
approx = 0.5;
reratio = 3;
flags = varargin2struct(varargin{:});
if (~isstruct(opt))
ssize = opt;
end
if (isstruct(opt) && length(opt) == 1) % does not support settings for multiple labels
ssize = jsonopt('radbound', ssize, opt);
ang = jsonopt('angbound', ang, opt);
approx = jsonopt('distbound', approx, opt);
reratio = jsonopt('reratio', reratio, opt);
end
if (jsonopt('DoRepair', 0, flags) == 1)
[v, f] = meshcheckrepair(v, f);
end
saveoff(v, f, mwpath('pre_cgalpoly.off'));
deletemeshfile(mwpath('post_cgalpoly.mesh'));
randseed = hex2dec('623F9A9E'); % "U+623F U+9A9E"
if (~isempty(getvarfrom({'caller', 'base'}, 'ISO2MESH_RANDSEED')))
randseed = getvarfrom({'caller', 'base'}, 'ISO2MESH_RANDSEED');
end
cmd = sprintf('"%s%s" "%s" "%s" %.16f %.16f %.16f %.16f %.16f %d', mcpath('cgalpoly'), exesuff, ...
mwpath('pre_cgalpoly.off'), mwpath('post_cgalpoly.mesh'), ang, ssize, ...
approx, reratio, maxvol, randseed);
status = system(cmd);
if (status)
error('cgalpoly command failed');
end
if (~exist(mwpath('post_cgalpoly.mesh'), 'file'))
error('output file was not found, failure was encountered when running command: \n%s\n', cmd);
end
[node, elem, face] = readmedit(mwpath('post_cgalpoly.mesh'));
fprintf(1, 'node number:\t%d\ntriangles:\t%d\ntetrahedra:\t%d\nregions:\t%d\n', ...
size(node, 1), size(face, 1), size(elem, 1), length(unique(elem(:, end))));
fprintf(1, 'surface and volume meshes complete\n');
|
