1function msh = load_gmsh2(filename, which) 2 3%% Reads a mesh in msh format, version 1 or 2 4 5% Usage: 6% To define all variables m.LINES, M.TRIANGLES, etc 7% (Warning: This creates a very large structure. Do you really need it?) 8% m = load_gmsh2('a.msh') 9% 10% To define only certain variables (for example TETS and HEXS) 11% m = load_gmsh2('a.msh', [ 5 6]) 12% 13% To define no variables (i.e. if you prefer to use m.ELE_INFOS(i,2)) 14% m = load_gmsh2('a.msh', -1) 15% m = load_gmsh2('a.msh', []) 16% 17% Copyright (C) 2007 JP Moitinho de Almeida (moitinho@civil.ist.utl.pt) 18% and R Lorphevre(r(point)lorphevre(at)ulg(point)ac(point)be) 19% 20% based on load_gmsh.m supplied with gmsh-2.0 21 22% Structure msh always has the following elements: 23% 24% msh.MIN, msh.MAX - Bounding box 25% msh.nbNod - Number of nodes 26% msh.nbElm - Total number of elements 27% msh.nbType(i) - Number of elements of type i (i in 1:19) 28% msh.POS(i,j) - j'th coordinate of node i (j in 1:3, i in 1:msh.nbNod) 29% msh.ELE_INFOS(i,1) - id of element i (i in 1:msh.nbElm) 30% msh.ELE_INFOS(i,2) - type of element i 31% msh.ELE_INFOS(i,3) - number of tags for element i 32% msh.ELE_INFOS(i,4) - Dimension (0D, 1D, 2D or 3D) of element i 33% msh.ELE_TAGS(i,j) - Tags of element i (j in 1:msh.ELE_INFOS(i,3)) 34% msh.ELE_NODES(i,j) - Nodes of element i (j in 1:k, with 35% k = msh.NODES_PER_TYPE_OF_ELEMENT(msh.ELE_INFOS(i,2))) 36% 37% These elements are created if requested: 38% 39% msh.nbLines - number of 2 node lines 40% msh.LINES(i,1:2) - nodes of line i 41% msh.LINES(i,3) - tag (WHICH ?????) of line i 42% 43% msh.nbTriangles - number of 2 node triangles 44% msh.TRIANGLES(i,1:3) - nodes of triangle i 45% msh.TRIANGLES(i,4) - tag (WHICH ?????) of triangle i 46% 47% Etc 48 49% These definitions need to be updated when new elemens types are added to gmsh 50% 51% msh.Types{i}{1} Number of an element of type i 52% msh.Types{i}{2} Dimension (0D, 1D, 2D or 3D) of element of type i 53% msh.Types{i}{3} Name to add to the structure with elements of type i 54% msh.Types{i}{4} Name to add to the structure with the number of elements of type i 55% 56 57nargchk(1, 2, nargin); 58 59msh.Types = { ... 60 { 2, 1, 'LINES', 'nbLines'}, ... % 1 61 { 3, 2, 'TRIANGLES', 'nbTriangles'}, ... 62 { 4, 2, 'QUADS', 'nbQuads'}, ... 63 { 4, 3, 'TETS', 'nbTets'}, ... 64 { 8, 3, 'HEXAS', 'nbHexas'}, ... %5 65 { 6, 3, 'PRISMS', 'nbPrisms'}, ... 66 { 5, 3, 'PYRAMIDS', 'nbPyramids'}, ... 67 { 3, 1, 'LINES3', 'nbLines3'}, ... 68 { 6, 2, 'TRIANGLES6', 'nbTriangles6'}, ... 69 { 9, 2, 'QUADS9', 'nbQuads9'}, ... % 10 70 { 10, 3, 'TETS10', 'nbTets10'}, ... 71 { 27, 3, 'HEXAS27', 'nbHexas27'}, ... 72 { 18, 3, 'PRISMS18', 'nbPrisms18'}, ... 73 { 14, 3, 'PYRAMIDS14', 'nbPyramids14'}, ... 74 { 1, 0, 'POINTS', 'nbPoints'}, ... % 15 75 { 8, 3, 'QUADS8', 'nbQuads8'}, ... 76 { 20, 3, 'HEXAS20', 'nbHexas20'}, ... 77 { 15, 3, 'PRISMS15', 'nbPrisms15'}, ... 78 { 13, 3, 'PYRAMIDS13', 'nbPyramids13'}, ... 79}; 80 81ntypes = length(msh.Types); 82 83if (nargin==1) 84 which = 1:ntypes; 85else 86 if isscalar(which) && which == -1 87 which = []; 88 end 89end 90 91% Could check that "which" is properlly defined.... 92 93fid = fopen(filename, 'r'); 94fileformat = 0; % Assume wrong file 95 96tline = fgetl(fid); 97if (feof(fid)) 98 disp (sprintf('Empty file: %s', filename)); 99 msh = []; 100 return; 101end 102 103%% Read mesh type 104if (strcmp(tline, '$NOD')) 105 fileformat = 1; 106elseif (strcmp(tline, '$MeshFormat')) 107 fileformat = 2; 108 tline = fgetl(fid); 109 if (feof(fid)) 110 disp (sprintf('Syntax error (no version) in: %s', filename)); 111 fileformat = 0; 112 end 113 [ form ] = sscanf( tline, '%f %d %d'); 114 if (form(1) < 2 || form(1) > 2.2) 115 disp (sprintf('Unknown mesh format: %s', tline)); 116 fileformat = 0; 117 end 118 if (form(2) ~= 0) 119 disp ('Sorry, this program can only read ASCII format'); 120 fileformat = 0; 121 end 122 fgetl(fid); % this should be $EndMeshFormat 123 if (feof(fid)) 124 disp (sprintf('Syntax error (no $EndMeshFormat) in: %s', filename)); 125 fileformat = 0; 126 end 127 tline = fgetl(fid); % this should be $Nodes or $PhysicalNames 128end 129 130if (~fileformat) 131 msh = []; 132 return 133end 134 135%% Read (and ignore) physical names 136 137if strcmp(tline, '$PhysicalNames') 138 np = fscanf(fid, '%d', 1); 139 tline = fgetl(fid); 140 for I = 1:np+2 141 tline = fgetl(fid); 142 end 143end 144 145%% Read nodes 146 147if strcmp(tline, '$NOD') || strcmp(tline, '$Nodes') 148 msh.nbNod = fscanf(fid, '%d', 1); 149 aux = fscanf(fid, '%g', [4 msh.nbNod]); 150 msh.POS = aux(2:4,:)'; 151 numids = max(aux(1,:)); 152 IDS = zeros(1, numids); 153 IDS(aux(1,:)) = 1:msh.nbNod; % This may not be an identity 154 msh.MAX = max(msh.POS); 155 msh.MIN = min(msh.POS); 156 fgetl(fid); % End previous line 157 fgetl(fid); % Has to be $ENDNOD $EndNodes 158else 159 disp (sprintf('Syntax error (no $Nodes/$NOD) in: %s', filename)); 160 fileformat = 0; 161end 162 163%% Read elements 164tline = fgetl(fid); 165if strcmp(tline,'$ELM') || strcmp(tline, '$Elements') 166 msh.nbElm = fscanf(fid, '%d', 1); 167 % read all info about elements into aux (it is faster!) 168 aux = fscanf(fid, '%g', inf); 169 start = 1; 170 msh.ELE_INFOS = zeros(msh.nbElm, 4); % 1 - id, 2 - type, 3 - ntags, 4 - Dimension 171 msh.ELE_NODES = zeros(msh.nbElm,6); % i - Element, j - ElNodes 172 if (fileformat == 1) 173 ntags = 2; 174 else 175 ntags = 3; % just a prediction 176 end 177 msh.ELE_TAGS = zeros(msh.nbElm, ntags); % format 1: 1 - physical number, 2 - geometrical number 178 % format 2: 1 - physical number, 2 - geometrical number, 3 - mesh partition number 179 msh.nbType = zeros(ntypes,1); 180 for I = 1:msh.nbElm 181 if (fileformat >= 2) 182 finnish = start + 2; 183 msh.ELE_INFOS(I, 1:3) = aux(start:finnish); 184 ntags = aux(finnish); 185 start = finnish + 1; 186 finnish = start + ntags -1; 187 msh.ELE_TAGS(I, 1:ntags) = aux(start:finnish); 188 start = finnish + 1; 189 else 190 finnish = start + 1; 191 msh.ELE_INFOS(I, 1:2) = aux(start:finnish); 192 start = finnish + 1; % the third element is nnodes, which we get from the type 193 msh.ELE_INFOS(I, 3) = 2; 194 finnish = start + 1; 195 msh.ELE_TAGS(I, 1:2) = aux(start:finnish); 196 start = finnish + 2; 197 end 198 type = msh.ELE_INFOS(I, 2); 199 msh.nbType(type) = msh.nbType(type) + 1; 200 msh.ELE_INFOS(I, 4) = msh.Types{type}{2}; 201 nnodes = msh.Types{type}{1}; 202 finnish = start + nnodes - 1; 203 msh.ELE_NODES(I, 1:nnodes) = IDS(aux(start:finnish)); 204 start=finnish + 1; 205 end 206 fgetl(fid); % Has to be $ENDELM or $EndElements 207else 208 disp (sprintf('Syntax error (no $Elements/$ELM) in: %s', filename)); 209 fileformat = 0; 210end 211 212if (fileformat == 0) 213 msh = []; 214end 215 216fclose(fid); 217 218%% This is used to create the explicit lists for types of elements 219 220for i = which 221 if (~isempty(msh.Types{i}{3})) 222 cmd = sprintf('msh.%s=msh.nbType(%d);', msh.Types{i}{4}, i); 223 eval(cmd); 224 % Dimension 225 cmd = sprintf('msh.%s=zeros(%d,%d);', msh.Types{i}{3}, msh.nbType(i), msh.Types{i}{1}+1); 226 eval(cmd); 227 % Clear nbType for counting, next loop will recompute it 228 msh.nbType(i) = 0; 229 end 230end 231 232for i = 1:msh.nbElm 233 type = msh.ELE_INFOS(i,2); 234 if (find(which == type)) 235 if (~isempty(msh.Types{type}{3})) 236 msh.nbType(type) = msh.nbType(type)+1; 237 aux=[ msh.ELE_NODES(i,1:msh.Types{type}{1}), msh.ELE_TAGS(i,1) ]; 238 cmd = sprintf('msh.%s(%d,:)=aux;', msh.Types{type}{3}, msh.nbType(type)); 239 eval(cmd); 240 end 241 end 242end 243 244return; 245