1function gmesh = load_gmsh ( filename ) 2 [fid,ierr] = mopen ( filename, 'r'); 3 if ierr <> 0 then; 4 mclose(fid); 5 error('no such file: '+filename) 6 end 7 gmesh = []; 8 gmesh.MIN = zeros(3,1); 9 gmesh.MAX = zeros(3,1); 10 11 while 1 12 endoffile = 0; 13 14 while (1) 15 // jump over irrelevant text 16 tline = mgetl(fid,1); 17 if meof(fid), endoffile=1,break, end 18 if (part(tline,1) == '$'), break, end 19 end 20 if (endoffile == 1), break, end 21 if (part(tline,1:4) == '$NOD' | part(tline,1:4) == '$Nod') 22 disp('reading nodes') 23 gmesh.nbNod = mfscanf (1, fid, '%d'); 24 gmesh.POS = zeros (gmesh.nbNod, 3); 25 for(I=1:gmesh.nbNod) 26 iNod = mfscanf(1,fid,'%d'); 27 X = mfscanf(3,fid,'%g'); 28 IDS(iNod) = I; 29 if (I == 1) 30 gmesh.MIN = X; 31 gmesh.MAX = X; 32 else 33 if (gmesh.MAX(1) < X(1)), gmesh.MAX(1)=X(1);end 34 if (gmesh.MAX(2) < X(2)), gmesh.MAX(2)=X(2);end 35 if (gmesh.MAX(3) < X(3)), gmesh.MAX(3)=X(3);end 36 if (gmesh.MIN(1) > X(1)), gmesh.MIN(1)=X(1);end 37 if (gmesh.MIN(2) > X(2)), gmesh.MIN(2)=X(2);end 38 if (gmesh.MIN(3) > X(3)), gmesh.MIN(3)=X(3);end 39 end 40 gmesh.POS(I,1) = X(1); 41 gmesh.POS(I,2) = X(2); 42 gmesh.POS(I,3) = X(3); 43 end 44 // Here we have to read two lines, instead of one, in order to continue. 45 // The reason may be due to an unread end-of-line character from the previous 46 // scan. 47 tline = mgetl(fid,2); // read 2 dummy lines 48 disp('nodes have been read') 49 50 elseif(part(tline,1:4)=='$ELM' | part(tline,1:4)=='$Ele') 51 disp('reading elements') 52 gmesh.nbElm = mfscanf (1,fid, '%d'); 53 gmesh.ELE_INFOS = zeros (gmesh.nbElm,5); 54 gmesh.nbLines = 0; 55 gmesh.nbPoints = 0; 56 gmesh.nbTriangles = 0; 57 gmesh.nbQuads = 0; 58 gmesh.nbPrism=0; 59 gmesh.nbPyr=0; 60 gmesh.nbTet = 0; 61 gmesh.nbHex = 0; 62 gmesh.nbQTriangles = 0; 63 gmesh.nbQQuads = 0; 64 gmesh.POINTS=[]; 65 gmesh.LINES=[]; 66 gmesh.TRIANGLES=[]; 67 gmesh.QUADS=[]; 68 gmesh.TETS=[]; 69 gmesh.HEXA=[]; 70 gmesh.PRISM=[]; 71 gmesh.PYRAMID=[]; 72 gmesh.QTRIANGLES=[]; 73 gmesh.QQUADS=[]; 74 for (I=1:gmesh.nbElm) 75 if (part(tline,1:4)=='$Ele') 76 // file is in GMSH version 2.0 format 77 gmesh.ELE_INFOS(I,1:3) = mfscanf(3,fid,'%d')'; 78 n_of_tags = gmesh.ELE_INFOS(I,3); 79 gmesh.ELE_INFOS(I,4+(1:n_of_tags)) = mfscanf(n_of_tags,fid,'%d')'; 80 // the number of points 81 select gmesh.ELE_INFOS(I,2) 82 case 15 then, np=1; 83 case 1 then, np=2; 84 case 2 then, np=3; 85 case 3 then, np=4; 86 case 4 then, np=4; 87 case 5 then, np=8; 88 case 6 then, np=6; 89 case 7 then, np=5; 90 case 9 then, np=6; 91 case 10 then, np=9; 92 end 93 NODES_ELEM = mfscanf(np,fid,'%d')'; 94 else 95 // version 1.0 format 96 gmesh.ELE_INFOS(I,:) = mfscanf(5,fid,'%d')'; 97 NODES_ELEM = mfscanf(gmesh.ELE_INFOS(I,5),fid,'%d')'; 98 end 99 select gmesh.ELE_INFOS(I,2) 100 case 15 then // point 101 gmesh.nbPoints = gmesh.nbPoints + 1; 102 gmesh.POINTS =[gmesh.POINTS; IDS(NODES_ELEM (1))' I]; 103 case 1 then // beam 104 gmesh.nbLines = gmesh.nbLines + 1; 105 gmesh.LINES =[gmesh.LINES; IDS(NODES_ELEM (1:2))' I]; 106 case 2 then // triangle 107 gmesh.nbTriangles = gmesh.nbTriangles + 1; 108 gmesh.TRIANGLES = [gmesh.TRIANGLES; IDS(NODES_ELEM (1:3))' I]; 109 case 3 then // quadrangle 110 gmesh.nbQuads = gmesh.nbQuads + 1; 111 gmesh.QUADS =[gmesh.QUADS; IDS(NODES_ELEM (1:4))' I]; 112 case 4 then // tetrahedron (4 node) 113 gmesh.nbTet = gmesh.nbTet + 1; 114 gmesh.TETS = [gmesh.TETS; IDS(NODES_ELEM (1:4))' I]; 115 case 5 then // hexahedron (8 nodes) 116 gmesh.nbHex = gmesh.nbHex + 1; 117 gmesh.HEXA =[gmesh.HEXA; IDS(NODES_ELEM (1:8))' I]; 118 case 6 then // prism (6 nodes) 119 gmesh.nbPrism = gmesh.nbPrism + 1; 120 gmesh.PRISM =[gmesh.PRISM; IDS(NODES_ELEM (1:6))' I]; 121 case 7 then // pyramid (5 nodes) 122 gmesh.nbPyr = gmesh.nbPyr + 1; 123 gmesh.PYRAMID =[gmesh.PYRAMID; IDS(NODES_ELEM (1:5))' I]; 124 case 9 then // second order triangle (6 nodes) 125 gmesh.nbQTriangles = gmesh.nbQTriangles + 1; 126 gmesh.QTRIANGLES = [gmesh.QTRIANGLES; IDS(NODES_ELEM (1:6))' I]; 127 case 10 then // second order quadrangle (9 nodes) 128 gmesh.nbQQuads = gmesh.nbQQuads + 1; 129 gmesh.QQUADS = [gmesh.QQUADS; IDS(NODES_ELEM (1:9))' I]; 130 else 131 disp(' ') 132 warning('Unknown element type ' + string(gmesh.ELE_INFOS(I,2)) + ' !') 133 end 134 end 135 disp('elements have been read') 136 tline = mgetl(fid,1); 137 end // 138 end // while 139mclose (fid); 140endfunction 141