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