[Gmsh] load_gmsh.sci file

Christophe Geuzaine cag32 at case.edu
Sat May 13 23:28:49 CEST 2006


Gürsoy Turan wrote:
> Earlier I had send the scilab version of the load_gmsh.m file. Back then it 
> only worked for GMSH version 1.0 formatted files. Now, I included version 2.0
> Bye.

Thanks!



> 
> 
> ------------------------------------------------------------------------
> 
> 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
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> gmsh mailing list
> gmsh at geuz.org
> http://www.geuz.org/mailman/listinfo/gmsh


-- 
Christophe Geuzaine
Assistant Professor, Case Western Reserve University, Mathematics
http://www.case.edu/artsci/math/geuzaine