script from Arthur Levy

utils/converters/matlab/gmsh2pdetoolbox.m

+%GMSH2PDETOOLBOX Reads a mesh in msh format, version 1 or 2 and returns the
+%arrays p, e and t according to the MATLAB pde toolbox syntax. Filename
+%refers to the .msh file. Note that only triangular 2D mesh are processed
+%since these are the only elements allowed in pde toolbox.
+%
+%SEE ALSO load_gmsh load_gmsh4 initmesh
+%
+% Copyright (C) 2014 Arthur Levy. https://github.com/arthurlevy/
+%
+% This fucntion uses the routine load_gmsh4 in load_gmsh2.m: Copyright (C)
+% 2007  JP Moitinho de Almeida (moitinho@civil.ist.utl.pt) and  R
+% Lorphevre(r(point)lorphevre(at)ulg(point)ac(point)be) It is based on
+% load_gmsh.m supplied with gmsh-2.0
+function [p, e, t] = gmsh2pdetoolbox(filename)
+
+%loads the gmsh file using JP Moitinho de Almeida (moitinho@civil.ist.utl.pt)
+% and  R Lorphevre(r(point)lorphevre(at)ulg(point)ac(point)be) code.
+mesh_structure = load_gmsh4(filename,-1);
+
+
+%% NODES
+%nodes positions
+disp('importing nodes');
+p = mesh_structure.POS(:,1:2)';
+
+
+%% EDGES
+disp('importing edges');
+%find the edges (ie dimension 1)
+is_edge = (mesh_structure.ELE_INFOS(:,2)==1);
+%add edge data
+e = mesh_structure.ELE_NODES(is_edge,1:2)';
+e(3,:) = 0;
+e(4,:) = 1;
+%tag is important for applying boundary conditions
+e(5,:) = mesh_structure.ELE_TAGS(is_edge,1)';
+e(6,:) = 1;
+e(7,:) = 0;
+
+%% TRIANGLES
+disp('importing triangles')
+is_triangle = (mesh_structure.ELE_INFOS(:,2)==2);
+t = mesh_structure.ELE_NODES(is_triangle,1:3)';
+t(4,:) = 2;
+
+
+
+
+
+function msh = load_gmsh4(filename, which)
+
+%% Reads a mesh in msh format, version 1 or 2
+%
+% Usage: 
+% To define all variables m.LINES, M.TRIANGLES, etc 
+% (Warning: This creates a very large structure. Do you really need it?)
+%            m = load_gmsh4('a.msh')
+%
+% To define only certain variables (for example TETS and HEXS)
+%            m = load_gmsh4('a.msh', [ 5 6])
+%
+% To define no variables (i.e. if you prefer to use m.ELE_INFOS(i,2))
+%            m = load_gmsh4('a.msh', -1)
+%            m = load_gmsh4('a.msh', [])
+%
+% Copyright (C) 2007  JP Moitinho de Almeida (moitinho@civil.ist.utl.pt)
+% and  R Lorphevre(r(point)lorphevre(at)ulg(point)ac(point)be)
+%
+% based on load_gmsh.m supplied with gmsh-2.0
+%
+% Structure msh always has the following elements:
+%
+% msh.MIN, msh.MAX - Bounding box
+% msh.nbNod - Number of nodes
+% msh.nbElm - Total number of elements
+% msh.nbType(i) - Number of elements of type i (i in 1:19)
+% msh.POS(i,j) - j'th coordinate of node i (j in 1:3, i in 1:msh.nbNod)
+% msh.ELE_INFOS(i,1) - id of element i (i in 1:msh.nbElm)
+% msh.ELE_INFOS(i,2) - type of element i
+% msh.ELE_INFOS(i,3) - number of tags for element i
+% msh.ELE_INFOS(i,4) - Dimension (0D, 1D, 2D or 3D) of element i
+% msh.ELE_TAGS(i,j) - Tags of element i (j in 1:msh.ELE_INFOS(i,3))
+% msh.ELE_NODES(i,j) - Nodes of element i (j in 1:k, with
+%                       k = msh.NODES_PER_TYPE_OF_ELEMENT(msh.ELE_INFOS(i,2)))
+%
+% These elements are created if requested:
+%
+% msh.nbLines - number of 2 node lines
+% msh.LINES(i,1:2) - nodes of line i
+% msh.LINES(i,3) - tag (WHICH ?????) of line i
+%
+% msh.nbTriangles - number of 2 node triangles
+% msh.TRIANGLES(i,1:3) - nodes of triangle i
+% msh.TRIANGLES(i,4) - tag (WHICH ?????) of triangle i
+%
+% Etc
+
+% These definitions need to be updated when new elemens types are added to gmsh
+%
+% msh.Types{i}{1} Number of an element of type i
+% msh.Types{i}{2} Dimension (0D, 1D, 2D or 3D) of element of type i
+% msh.Types{i}{3} Name to add to the structure with elements of type i
+% msh.Types{i}{4} Name to add to the structure with the number of elements of type i
+%
+
+nargchk(1, 2, nargin);
+
+msh.Types = { ...
+    { 2, 1, 'LINES', 'nbLines'}, ... % 1
+    { 3,  2, 'TRIANGLES', 'nbTriangles'}, ...
+    { 4,  2, 'QUADS', 'nbQuads'}, ...  
+    { 4,  3, 'TETS', 'nbTets'}, ...
+    { 8,  3, 'HEXAS', 'nbHexas'}, ... %5
+    { 6,  3, 'PRISMS', 'nbPrisms'}, ...
+    { 5,  3, 'PYRAMIDS', 'nbPyramids'}, ...
+    { 3,  1, 'LINES3', 'nbLines3'}, ...
+    { 6,  2, 'TRIANGLES6', 'nbTriangles6'}, ...
+    { 9,  2, 'QUADS9', 'nbQuads9'}, ... % 10
+    { 10,  3, 'TETS10', 'nbTets10'}, ...
+    { 27,  3, 'HEXAS27', 'nbHexas27'}, ...
+    { 18,  3, 'PRISMS18', 'nbPrisms18'}, ...
+    { 14,  3, 'PYRAMIDS14', 'nbPyramids14'}, ...
+    { 1,  0, 'POINTS', 'nbPoints'}, ... % 15
+    { 8,  3, 'QUADS8', 'nbQuads8'}, ...
+    { 20,  3, 'HEXAS20', 'nbHexas20'}, ...
+    { 15,  3, 'PRISMS15', 'nbPrisms15'}, ...
+    { 13,  3, 'PYRAMIDS13', 'nbPyramids13'}, ...
+};
+                     
+ntypes = length(msh.Types);
+
+if (nargin==1)
+    which = 1:ntypes;
+else
+    if isscalar(which) && which == -1
+        which = [];
+    end
+end
+
+% Could check that "which" is properlly defined....
+
+fid = fopen(filename, 'r');
+fileformat = 0; % Assume wrong file
+
+tline = fgetl(fid);
+if (feof(fid))
+    disp (sprintf('Empty file: %s',  filename));
+    msh = [];
+    return;
+end
+
+%% Read mesh type
+if (strcmp(tline, '$NOD'))
+    fileformat = 1;
+elseif (strcmp(tline, '$MeshFormat'))
+    fileformat = 2;
+    tline = fgetl(fid);
+    if (feof(fid))
+        disp (sprintf('Syntax error (no version) in: %s',  filename));
+        fileformat = 0;
+    end
+    [ form ] = sscanf( tline, '%f %d %d');
+%    if (form(1) ~= 2)
+%        disp (sprintf('Unknown mesh format: %s', tline));
+%        fileformat = 0;
+%    end
+    if (form(2) ~= 0)
+        disp ('Sorry, this program can only read ASCII format');
+        fileformat = 0;
+    end
+    fgetl(fid);    % this should be $EndMeshFormat
+    if (feof(fid))
+        disp (sprintf('Syntax error (no $EndMeshFormat) in: %s',  filename));
+        fileformat = 0;
+    end
+    tline = fgetl(fid);    % this should be $Nodes
+    if (feof(fid))
+        disp (sprintf('Syntax error (no $Nodes) in: %s',  filename));
+        fileformat = 0;
+    end
+end
+
+if (~fileformat)
+    msh = [];
+    return
+end
+
+%% Read nodes
+
+if strcmp(tline, '$NOD') || strcmp(tline, '$Nodes')
+    msh.nbNod = fscanf(fid, '%d', 1);
+    aux = fscanf(fid, '%g', [4 msh.nbNod]);
+    msh.POS = aux(2:4,:)';
+    numids = max(aux(1,:));
+    IDS = zeros(1, numids);
+    IDS(aux(1,:)) = 1:msh.nbNod; % This may not be an identity
+    msh.MAX = max(msh.POS);
+    msh.MIN = min(msh.POS);
+    fgetl(fid); % End previous line
+    fgetl(fid); % Has to be $ENDNOD $EndNodes
+else
+    disp (sprintf('Syntax error (no $Nodes/$NOD) in: %s',  filename));
+    fileformat = 0;
+end
+
+%% Read elements
+tline = fgetl(fid);
+if strcmp(tline,'$ELM') || strcmp(tline, '$Elements')
+    msh.nbElm = fscanf(fid, '%d', 1);
+    % read all info about elements into aux (it is faster!)
+    aux = fscanf(fid, '%g', inf);
+    start = 1;
+    msh.ELE_INFOS = zeros(msh.nbElm, 4); % 1 - id, 2 - type, 3 - ntags, 4 - Dimension
+    msh.ELE_NODES = zeros(msh.nbElm,6); % i - Element, j - ElNodes
+    if (fileformat == 1)
+        ntags = 2;
+    else
+        ntags = 3; % just a prediction
+    end
+    msh.ELE_TAGS = zeros(msh.nbElm, ntags); % format 1: 1 - physical number, 2 - geometrical number
+                                            % format 2: 1 - physical number, 2 - geometrical number, 3 - mesh partition number
+    msh.nbType = zeros(ntypes,1);
+    for I = 1:msh.nbElm
+        if (fileformat == 2)
+            finnish = start + 2; 
+            msh.ELE_INFOS(I, 1:3) = aux(start:finnish); 
+            ntags = aux(finnish); 
+            start = finnish + 1;
+            finnish = start + ntags -1;
+            msh.ELE_TAGS(I, 1:ntags) = aux(start:finnish); 
+            start = finnish + 1;
+        else
+            finnish = start + 1;
+            msh.ELE_INFOS(I, 1:2) = aux(start:finnish); 
+            start = finnish + 1; % the third element is nnodes, which we get from the type
+            msh.ELE_INFOS(I, 3) = 2; 
+            finnish = start + 1;
+            msh.ELE_TAGS(I, 1:2) = aux(start:finnish); 
+            start = finnish + 2;
+        end
+        type = msh.ELE_INFOS(I, 2);
+        msh.nbType(type) = msh.nbType(type) + 1;
+        msh.ELE_INFOS(I, 4) = msh.Types{type}{2};
+        nnodes = msh.Types{type}{1};
+        finnish = start + nnodes - 1;
+        msh.ELE_NODES(I, 1:nnodes) = IDS(aux(start:finnish));
+        start=finnish + 1;
+    end
+    fgetl(fid); % Has to be $ENDELM or $EndElements
+else
+    disp (sprintf('Syntax error (no $Elements/$ELM) in: %s',  filename));
+    fileformat = 0;
+end
+
+if (fileformat == 0)
+    msh = [];
+end
+
+fclose(fid);
+
+%% This is used to create the explicit lists for types of elements
+
+for i = which
+    if (~isempty(msh.Types{i}{3}))
+        cmd = sprintf('msh.%s=msh.nbType(%d);', msh.Types{i}{4}, i);
+        eval(cmd);
+        % Dimension
+        cmd = sprintf('msh.%s=zeros(%d,%d);', msh.Types{i}{3}, msh.nbType(i), msh.Types{i}{1}+1);
+        eval(cmd);
+        % Clear nbType for counting, next loop will recompute it
+        msh.nbType(i) = 0;
+    end
+end
+
+for i = 1:msh.nbElm
+    type = msh.ELE_INFOS(i,2);
+    if (find(which == type))
+        if (~isempty(msh.Types{type}{3}))
+            msh.nbType(type) = msh.nbType(type)+1;
+            aux=[ msh.ELE_NODES(i,1:msh.Types{type}{1}), msh.ELE_TAGS(i,1) ];
+            cmd = sprintf('msh.%s(%d,:)=aux;', msh.Types{type}{3}, msh.nbType(type));
+            eval(cmd);
+        end
+    end
+end
+
+return;