savemsh

 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)

0001 function savemsh(node,elem,fname,rname)
0002 %
0003 % savemsh(node,elem,fname,rname)
0004 %
0005 % save a tetrahedral mesh to GMSH mesh format
0006 %
0007 % author: Riccardo Scorretti (riccardo.scorretti<at> univ-lyon1.fr)
0008 % date: 2013/07/22
0009 %
0010 % input:
0011 %      node: input, node list, dimension (nn,3)
0012 %      elem: input, tetrahedral mesh element list, dimension (ne,4) or (ne,5) for multi-region meshes
0013 %      fname: output file name
0014 %      rname: name of the regions, cell-array of strings (optional)
0015 %
0016 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0017 %
0018 
0019 if nargin < 4 , rname = {} ; end
0020 if size(elem,2) < 5 , elem(:,5) = 1 ; end
0021 fid = fopen(fname,'wt');
0022 if(fid==-1)
0023     error('You do not have permission to save mesh files.');
0024 end
0025 
0026 % Check that all the elements are correctly oriented
0027 elem(:,1:4)=meshreorient(node,elem(:,1:4));
0028 
0029 nbNodes = size (node,1);
0030 reg = unique (elem(:,5));
0031 nbRegion = length (reg);
0032 nbElements = size (elem, 1);
0033 
0034 % Create the skeleton of the mesh structure
0035 M.Info.version = [];
0036 
0037 M.Nodes.nb = 0;
0038 M.Nodes.x = [];
0039 M.Nodes.y = [];
0040 M.Nodes.z = [];
0041 
0042 M.Elements.nb = 0;
0043 M.Elements.type = zeros (0, 0, 'uint8');
0044 M.Elements.tableOfNodes = zeros (0, 0, 'uint32');
0045 M.Elements.region = zeros (0, 0, 'uint16');
0046 
0047 M.Regions.nb = 0;
0048 M.Regions.name = {};
0049 M.Regions.dimension = [];
0050 
0051 % Build the table of nodes
0052 
0053 M.Nodes.nb = nbNodes;
0054 M.Nodes.x = node(:,1);
0055 M.Nodes.y = node(:,2);
0056 M.Nodes.z = node(:,3);
0057 clear node
0058 
0059 % Build the table of elements
0060 
0061 M.Elements.nb = nbElements;
0062 M.Elements.type = uint8(4*ones(nbElements, 1));
0063 M.Elements.tableOfNodes = uint32(elem(:,1:4));
0064 M.Elements.region = uint16(elem(:,5));
0065 clear elem
0066 
0067 % Build the table of regions
0068 
0069 M.Regions.nb = max(reg);
0070 for k = 1 : nbRegion
0071    if length(rname) < k , rname{k} = sprintf('region_%d', k) ; end
0072    M.Regions.name{reg(k)} = sprintf ('%s', rname{k});
0073    M.Regions.dimension(reg(k)) = 3;
0074 end
0075 
0076 % Writhe the header
0077 fprintf (fid, '$MeshFormat\n2.2 0 8\n$EndMeshFormat\n');
0078 
0079 % Write the physical names
0080 if M.Regions.nb > 0 
0081    fprintf (fid, '$PhysicalNames\n');
0082    fprintf (fid, '%d\n', M.Regions.nb);
0083    for r = 1 : M.Regions.nb
0084       name = M.Regions.name{r};
0085       if isempty (name)
0086          name = sprintf ('Region_%d', r);
0087       end
0088       fprintf (fid, '%d %d "%s"\n', M.Regions.dimension(r), r, name);
0089    end
0090    fprintf (fid, '$EndPhysicalNames\n');
0091 end
0092 
0093 % Write the nodes
0094 fprintf (fid, '$Nodes\n');
0095 fprintf (fid, '%d\n', size(M.Nodes.x,1));
0096 buffer = [ 1:M.Nodes.nb ; M.Nodes.x' ; M.Nodes.y' ; M.Nodes.z' ];
0097 fprintf (fid, '%d %10.10f %10.10f %10.10f\n', buffer);
0098 fprintf (fid, '$EndNodes\n');
0099 
0100 % Write the elements
0101 %
0102 % In order to accelerate the printing, the elements are printed in groups of (blockSize) elements, and
0103 % are grouped by (homogeneous) type. This variable sets the size of each group.
0104 %
0105 blockSize = 100000;
0106 
0107 fprintf (fid, '$Elements\n');
0108 fprintf (fid, '%d\n', M.Elements.nb);
0109 for h = 1 : blockSize : ceil(length(M.Elements.type)/blockSize)*blockSize
0110    e = h : min(length(M.Elements.type) , h+blockSize-1);  % = elements being considered
0111    type = unique (M.Elements.type(e));                    % = types of elements found in this group
0112    
0113    %
0114    % Process each type of element separately
0115    %
0116    for k = 1 : length(type)
0117       if type(k) == 0, continue; end
0118       et = e(find(M.Elements.type(e) == type(k)));  % = elements of the group of the same type
0119       
0120       %
0121       % Determine the format for printing the elements
0122       %
0123       elementFormat = '%d %d %d %d %d %d\n';
0124       for n = 1 : 4
0125          elementFormat = [ elementFormat '%d ' ];
0126       end
0127       elementFormat = [ elementFormat '\n' ];
0128       
0129       %
0130       % Collect in a buffer all the data of the elements of index (et)
0131       %
0132       buffer = zeros (10 , length(et));
0133       buffer(1,:) = et;
0134       buffer(2,:) = type(k);
0135       buffer(3,:) = 3;
0136       buffer(4,:) = M.Elements.region(et);
0137       buffer(5,:) = M.Elements.region(et);
0138       buffer(6,:) = 0;
0139       for n = 1 : 4
0140          buffer(6+n,:) = M.Elements.tableOfNodes(et,n);
0141       end
0142       
0143       %
0144       % Print all the homogeneous elements in the group with a single instruction
0145       %
0146       fprintf (fid, elementFormat, buffer);
0147    end
0148 end
0149 fprintf (fid, '$EndElements\n');
0150 
0151 
0152 fclose(fid);