function gmesh = load_gmsh ( filename ) [fid,ierr] = mopen ( filename, 'r'); if ierr <> 0 then; mclose(fid); error('no such file: '+filename) end gmesh = []; gmesh.MIN = zeros(3,1); gmesh.MAX = zeros(3,1); while 1 endoffile = 0; while (1) // jump over irrelevant text tline = mgetl(fid,1); if meof(fid), endoffile=1,break, end if (part(tline,1) == '$'), break, end end if (endoffile == 1), break, end if (part(tline,1:4) == '$NOD' | part(tline,1:4) == '$Nod') disp('reading nodes') gmesh.nbNod = mfscanf (1, fid, '%d'); gmesh.POS = zeros (gmesh.nbNod, 3); for(I=1:gmesh.nbNod) iNod = mfscanf(1,fid,'%d'); X = mfscanf(3,fid,'%g'); IDS(iNod) = I; if (I == 1) gmesh.MIN = X; gmesh.MAX = X; else if (gmesh.MAX(1) < X(1)), gmesh.MAX(1)=X(1);end if (gmesh.MAX(2) < X(2)), gmesh.MAX(2)=X(2);end if (gmesh.MAX(3) < X(3)), gmesh.MAX(3)=X(3);end if (gmesh.MIN(1) > X(1)), gmesh.MIN(1)=X(1);end if (gmesh.MIN(2) > X(2)), gmesh.MIN(2)=X(2);end if (gmesh.MIN(3) > X(3)), gmesh.MIN(3)=X(3);end end gmesh.POS(I,1) = X(1); gmesh.POS(I,2) = X(2); gmesh.POS(I,3) = X(3); end // Here we have to read two lines, instead of one, in order to continue. // The reason may be due to an unread end-of-line character from the previous // scan. tline = mgetl(fid,2); // read 2 dummy lines disp('nodes have been read') elseif(part(tline,1:4)=='$ELM' | part(tline,1:4)=='$Ele') disp('reading elements') gmesh.nbElm = mfscanf (1,fid, '%d'); gmesh.ELE_INFOS = zeros (gmesh.nbElm,5); gmesh.nbLines = 0; gmesh.nbPoints = 0; gmesh.nbTriangles = 0; gmesh.nbQuads = 0; gmesh.nbPrism=0; gmesh.nbPyr=0; gmesh.nbTet = 0; gmesh.nbHex = 0; gmesh.nbQTriangles = 0; gmesh.nbQQuads = 0; gmesh.POINTS=[]; gmesh.LINES=[]; gmesh.TRIANGLES=[]; gmesh.QUADS=[]; gmesh.TETS=[]; gmesh.HEXA=[]; gmesh.PRISM=[]; gmesh.PYRAMID=[]; gmesh.QTRIANGLES=[]; gmesh.QQUADS=[]; for (I=1:gmesh.nbElm) if (part(tline,1:4)=='$Ele') // file is in GMSH version 2.0 format gmesh.ELE_INFOS(I,1:3) = mfscanf(3,fid,'%d')'; n_of_tags = gmesh.ELE_INFOS(I,3); gmesh.ELE_INFOS(I,4+(1:n_of_tags)) = mfscanf(n_of_tags,fid,'%d')'; // the number of points select gmesh.ELE_INFOS(I,2) case 15 then, np=1; case 1 then, np=2; case 2 then, np=3; case 3 then, np=4; case 4 then, np=4; case 5 then, np=8; case 6 then, np=6; case 7 then, np=5; case 9 then, np=6; case 10 then, np=9; end NODES_ELEM = mfscanf(np,fid,'%d')'; else // version 1.0 format gmesh.ELE_INFOS(I,:) = mfscanf(5,fid,'%d')'; NODES_ELEM = mfscanf(gmesh.ELE_INFOS(I,5),fid,'%d')'; end select gmesh.ELE_INFOS(I,2) case 15 then // point gmesh.nbPoints = gmesh.nbPoints + 1; gmesh.POINTS =[gmesh.POINTS; IDS(NODES_ELEM (1))' I]; case 1 then // beam gmesh.nbLines = gmesh.nbLines + 1; gmesh.LINES =[gmesh.LINES; IDS(NODES_ELEM (1:2))' I]; case 2 then // triangle gmesh.nbTriangles = gmesh.nbTriangles + 1; gmesh.TRIANGLES = [gmesh.TRIANGLES; IDS(NODES_ELEM (1:3))' I]; case 3 then // quadrangle gmesh.nbQuads = gmesh.nbQuads + 1; gmesh.QUADS =[gmesh.QUADS; IDS(NODES_ELEM (1:4))' I]; case 4 then // tetrahedron (4 node) gmesh.nbTet = gmesh.nbTet + 1; gmesh.TETS = [gmesh.TETS; IDS(NODES_ELEM (1:4))' I]; case 5 then // hexahedron (8 nodes) gmesh.nbHex = gmesh.nbHex + 1; gmesh.HEXA =[gmesh.HEXA; IDS(NODES_ELEM (1:8))' I]; case 6 then // prism (6 nodes) gmesh.nbPrism = gmesh.nbPrism + 1; gmesh.PRISM =[gmesh.PRISM; IDS(NODES_ELEM (1:6))' I]; case 7 then // pyramid (5 nodes) gmesh.nbPyr = gmesh.nbPyr + 1; gmesh.PYRAMID =[gmesh.PYRAMID; IDS(NODES_ELEM (1:5))' I]; case 9 then // second order triangle (6 nodes) gmesh.nbQTriangles = gmesh.nbQTriangles + 1; gmesh.QTRIANGLES = [gmesh.QTRIANGLES; IDS(NODES_ELEM (1:6))' I]; case 10 then // second order quadrangle (9 nodes) gmesh.nbQQuads = gmesh.nbQQuads + 1; gmesh.QQUADS = [gmesh.QQUADS; IDS(NODES_ELEM (1:9))' I]; else disp(' ') warning('Unknown element type ' + string(gmesh.ELE_INFOS(I,2)) + ' !') end end disp('elements have been read') tline = mgetl(fid,1); end // end // while mclose (fid); endfunction
Generated by dwww version 1.15 on Thu Jun 27 09:42:43 CEST 2024.