gmsh_stl2tet

GMSH_STL2TET creates a tetrahedral mesh from an stl file
 mdl = gmsh_stl2tet(stlfile, maxh, extra, nopt) where:
        mdl - EIDORS fwd_model struct
    stlfile - either:
                 - a path to an stl file
               OR
                 - an EIDORS fwd_model with .nodes and .boundary or .elems
               OR
                 - a struct with .vertices and .faces
  mesh_size - controls mesh size. Either:
                 - scalar: maximum edge length 
                   (default: [] matches surface)
               OR
                 - [SizeMax, DistMax]: Linear grading from a fine surface,
                   achieving SizeMax at DistMax distance from surface.
      extra - extra command line options to gmsh (default: [])

       nopt - number of extra optimization runs (default: 0)

 If stlfile is a struct, stl_write will be called first and an STL file
 written in a temporary location. 

 CASHING: Calls are cashed iff stlfile is a struct.

 NOTE: Only one surface per file is allowed.

 See also CALL_GMSH, STL_WRITE

0001 function mdl = gmsh_stl2tet(stlfile, mesh_size, extra, nopt)
0002 %GMSH_STL2TET creates a tetrahedral mesh from an stl file
0003 % mdl = gmsh_stl2tet(stlfile, maxh, extra, nopt) where:
0004 %        mdl - EIDORS fwd_model struct
0005 %    stlfile - either:
0006 %                 - a path to an stl file
0007 %               OR
0008 %                 - an EIDORS fwd_model with .nodes and .boundary or .elems
0009 %               OR
0010 %                 - a struct with .vertices and .faces
0011 %  mesh_size - controls mesh size. Either:
0012 %                 - scalar: maximum edge length
0013 %                   (default: [] matches surface)
0014 %               OR
0015 %                 - [SizeMax, DistMax]: Linear grading from a fine surface,
0016 %                   achieving SizeMax at DistMax distance from surface.
0017 %      extra - extra command line options to gmsh (default: [])
0018 %
0019 %       nopt - number of extra optimization runs (default: 0)
0020 %
0021 % If stlfile is a struct, stl_write will be called first and an STL file
0022 % written in a temporary location.
0023 %
0024 % CASHING: Calls are cashed iff stlfile is a struct.
0025 %
0026 % NOTE: Only one surface per file is allowed.
0027 %
0028 % See also CALL_GMSH, STL_WRITE
0029 
0030 % (C) Bartlomiej Grychtol, 2012-2021.
0031 % $Id: gmsh_stl2tet.m 6983 2024-11-14 20:21:00Z bgrychtol $
0032 
0033 if nargin < 4, nopt = 0; end
0034 if nargin < 3, extra = []; end
0035 if nargin < 2, mesh_size = []; end
0036 
0037 if isstruct(stlfile)
0038     try
0039        name = stlfile.name;
0040     catch
0041        name = 'gmsh_stl2tet';
0042     end
0043     try
0044        opt.cache_obj{1} = stlfile.nodes;
0045        if isfield(stlfile, 'boundary')
0046           opt.cache_obj{2} = stlfile.boundary;
0047        else 
0048           opt.cache_obj{2} = stlfile.elems;
0049        end
0050     catch
0051        opt.cache_obj = {stlfile.vertices, stlfile.faces};
0052     end
0053     opt.cache_obj{3} = mesh_size;
0054     opt.cache_obj{4} = extra;
0055     opt.cache_obj{5} = nopt;
0056     mdl = eidors_cache(@do_gmsh_stl2tet,{stlfile, mesh_size, extra, nopt}, opt);
0057 else
0058     [~, name, ext] = fileparts(stlfile); name = [name ext];
0059     mdl = do_gmsh_stl2tet(stlfile, mesh_size, extra, nopt);
0060 end
0061 
0062 mdl = eidors_obj('fwd_model', mdl, 'name', name);
0063 
0064 
0065 function mdl = do_gmsh_stl2tet(stlfile, mesh_size, extra, nopt)
0066 if isstruct(stlfile)
0067     stem = tempname;
0068     stl_write(stlfile, [stem, '.stl'])
0069     stlname = [stem '.stl'];
0070 else
0071     stem = strrep(stlfile,'.stl','');
0072     stlname = stlfile;
0073 end
0074 %TODO: Some of this could be exposed as options (Algorithm, Optimize, ...)
0075 fid = fopen([stem '.geo'],'w');
0076 fprintf(fid,'General.Terminal=1;\n');
0077 fprintf(fid,'Merge "%s";\n',stlname);
0078 fprintf(fid,'Surface Loop(1) = {1};\n');
0079 fprintf(fid,'Volume(2) = {1};\n');
0080 fprintf(fid,'Physical Volume(3) = {2};\n');
0081 
0082 if numel(mesh_size) == 1
0083    fprintf(fid,'Mesh.MeshSizeMax = %f;\n', mesh_size);
0084 elseif numel(mesh_size) == 2
0085    % Only works for HPX - experimental!
0086    fprintf(fid, 'Field[1] = Extend;\n');
0087    fprintf(fid, 'Field[1].SurfacesList = {1};\n');
0088    fprintf(fid, 'Field[1].SizeMax = %g;\n', mesh_size(1));
0089    fprintf(fid, 'Field[1].DistMax = %g;\n', mesh_size(2));
0090    fprintf(fid, 'Field[1].Power = 1;\n');
0091    fprintf(fid, 'Background Field = 1;\n');
0092    fprintf(fid, 'Mesh.MeshSizeExtendFromBoundary = 0;\n');
0093    fprintf(fid, 'Mesh.Algorithm3D=10;\n'); % HPX (multi-threaded delaunay)
0094 end
0095 if numel(mesh_size) < 2
0096    fprintf(fid,'Mesh.Algorithm3D=4;\n'); % 4=frontal (netgen)
0097 end
0098 fprintf(fid,'Mesh.OptimizeNetgen=1;\n');
0099 fprintf(fid,'Mesh 3;\n');
0100 % It seems that if Gmsh chooses to optimize, it does Gmsh, followed by
0101 % Netgen (because OptimizeNetgen=1). This can be repeated.
0102 for i = 1:nopt
0103    fprintf(fid, 'OptimizeMesh "Gmsh";\n');
0104    fprintf(fid, 'OptimizeMesh "Netgen";\n');
0105 end
0106 % Gmsh optimization is quick and experience suggests it's worth running
0107 % after Netgen (again).
0108 fprintf(fid, 'OptimizeMesh "Gmsh";\n');
0109 fclose(fid);
0110 
0111 call_gmsh([stem '.geo'], extra);
0112 
0113 mdl = gmsh_mk_fwd_model([stem '.msh'],[],[],[]);
0114 
0115 delete([stem '.geo']);
0116 delete([stem '.msh']);
0117 if isstruct(stlfile)
0118     delete(stlname);
0119 end