converters/matlab/load_gmsh.m

1 function mesh = load_gmsh(filename)
2 % Reads a mesh in msh format, version 1 or 2
23 mesh = [];      variable
24 mesh.MIN = zeros(3, 1);
25 mesh.MAX = zeros(3, 1);
48           disp (' This program can only read ASCII mesh file')
62     mesh.nbNod = fscanf(fid, '%d', 1);
63     mesh.POS = zeros(mesh.nbNod, 3);
64     for(I=1:mesh.nbNod)
69 	mesh.MIN = X;
70 	mesh.MAX = X;
72 	if(mesh.MAX(1) < X(1)) mesh.MAX(1) = X(1); end
73 	if(mesh.MAX(2) < X(2)) mesh.MAX(2) = X(2); end
74 	if(mesh.MAX(3) < X(3)) mesh.MAX(3) = X(3); end
75 	if(mesh.MIN(1) > X(1)) mesh.MIN(1) = X(1); end
76 	if(mesh.MIN(2) > X(2)) mesh.MIN(2) = X(2); end
77 	if(mesh.MIN(3) > X(3)) mesh.MIN(3) = X(3); end
79       mesh.POS(I, 1) = X(1);
80       mesh.POS(I, 2) = X(2);
81       mesh.POS(I, 3) = X(3);
87     mesh.nbElm = fscanf(fid, '%d', 1);
96     mesh.ELE_INFOS = zeros(mesh.nbElm, nbinfo);
97     mesh.nbPoints = 0;
98     mesh.nbLines = 0;
99     mesh.nbTriangles = 0;
100     mesh.nbQuads = 0;
101     mesh.nbTets = 0;
102     mesh.nbHexas = 0;
103     mesh.nbPrisms = 0;
104     mesh.nbPyramids = 0;
106     mesh.POINTS = zeros(mesh.nbElm, 2);
107     mesh.LINES = zeros(mesh.nbElm, 3);
108     mesh.TRIANGLES = zeros(mesh.nbElm, 4);
109     mesh.QUADS = zeros(mesh.nbElm, 5);
110     mesh.TETS = zeros(mesh.nbElm, 5);
111     mesh.HEXAS = zeros(mesh.nbElm, 9);
112     mesh.PRISMS = zeros(mesh.nbElm, 7);
113     mesh.PYRAMIDS = zeros(mesh.nbElm, 6);
114     for(I = 1:mesh.nbElm)
115       mesh.ELE_INFOS(I, :) = fscanf(fid, '%d', nbinfo);
117           % take the mesh.ELE_INFOS(I, 5) nodes of the element
118           NODES_ELEM = fscanf(fid, '%d', mesh.ELE_INFOS(I, 5));
121           mesh.ELE_TAGS(I,:) = fscanf(fid, '%d', (mesh.ELE_INFOS(I,3)-1));
122           % take the msh.NODES_PER_TYPE_OF_ELEMENT (mesh.ELE_INFOS(I, 2)) nodes of the element
123           NODES_ELEM = fscanf(fid, '%d', msh.NODES_PER_TYPE_OF_ELEMENT (mesh.ELE_INFOS(I, 2)) );
125       if(mesh.ELE_INFOS(I, 2) == 15) %% point
126 	mesh.nbPoints = mesh.nbPoints + 1;
127     mesh.POINTS(mesh.nbPoints, 1) = IDS(NODES_ELEM(1));
128 	mesh.POINTS(mesh.nbPoints, 2) = mesh.ELE_INFOS(I, tags);
130       if(mesh.ELE_INFOS(I, 2) == 1) %% line
131 	mesh.nbLines = mesh.nbLines + 1;
132 	mesh.LINES(mesh.nbLines, 1) = IDS(NODES_ELEM(1));
133 	mesh.LINES(mesh.nbLines, 2) = IDS(NODES_ELEM(2));
134 	mesh.LINES(mesh.nbLines, 3) = mesh.ELE_INFOS(I, tags);
136       if(mesh.ELE_INFOS(I, 2) == 2) %% triangle
137 	mesh.nbTriangles = mesh.nbTriangles + 1;
138 	mesh.TRIANGLES(mesh.nbTriangles, 1) = IDS(NODES_ELEM(1));
139 	mesh.TRIANGLES(mesh.nbTriangles, 2) = IDS(NODES_ELEM(2));
140 	mesh.TRIANGLES(mesh.nbTriangles, 3) = IDS(NODES_ELEM(3));
141 	mesh.TRIANGLES(mesh.nbTriangles, 4) = mesh.ELE_INFOS(I, tags);
143       if(mesh.ELE_INFOS(I, 2) == 3) %% quadrangle
144 	mesh.nbQuads = mesh.nbQuads + 1;
145 	mesh.QUADS(mesh.nbQuads, 1) = IDS(NODES_ELEM(1));
146 	mesh.QUADS(mesh.nbQuads, 2) = IDS(NODES_ELEM(2));
147 	mesh.QUADS(mesh.nbQuads, 3) = IDS(NODES_ELEM(3));
148 	mesh.QUADS(mesh.nbQuads, 4) = IDS(NODES_ELEM(4));
149 	mesh.QUADS(mesh.nbQuads, 5) = mesh.ELE_INFOS(I, tags);
151       if(mesh.ELE_INFOS(I, 2) == 4) %% tetrahedron
152 	mesh.nbTets = mesh.nbTets + 1;
153 	mesh.TETS(mesh.nbTets, 1) = IDS(NODES_ELEM(1));
154 	mesh.TETS(mesh.nbTets, 2) = IDS(NODES_ELEM(2));
155 	mesh.TETS(mesh.nbTets, 3) = IDS(NODES_ELEM(3));
156 	mesh.TETS(mesh.nbTets, 4) = IDS(NODES_ELEM(4));
157 	mesh.TETS(mesh.nbTets, 5) = mesh.ELE_INFOS(I, tags);
159       if(mesh.ELE_INFOS(I, 2) == 5) %% hexahedron
160 	mesh.nbHexas = mesh.nbHexas + 1;
161 	mesh.HEXAS(mesh.nbHexas, 1) = IDS(NODES_ELEM(1));
162 	mesh.HEXAS(mesh.nbHexas, 2) = IDS(NODES_ELEM(2));
163 	mesh.HEXAS(mesh.nbHexas, 3) = IDS(NODES_ELEM(3));
164 	mesh.HEXAS(mesh.nbHexas, 4) = IDS(NODES_ELEM(4));
165 	mesh.HEXAS(mesh.nbHexas, 5) = IDS(NODES_ELEM(5));
166 	mesh.HEXAS(mesh.nbHexas, 6) = IDS(NODES_ELEM(6));
167 	mesh.HEXAS(mesh.nbHexas, 7) = IDS(NODES_ELEM(7));
168 	mesh.HEXAS(mesh.nbHexas, 8) = IDS(NODES_ELEM(8));
169 	mesh.HEXAS(mesh.nbHexas, 9) = mesh.ELE_INFOS(I, tags);
171       if(mesh.ELE_INFOS(I, 2) == 6) %% prism
172 	mesh.nbPrisms = mesh.nbPrisms + 1;
173 	mesh.PRISMS(mesh.nbPrisms, 1) = IDS(NODES_ELEM(1));
174 	mesh.PRISMS(mesh.nbPrisms, 2) = IDS(NODES_ELEM(2));
175 	mesh.PRISMS(mesh.nbPrisms, 3) = IDS(NODES_ELEM(3));
176 	mesh.PRISMS(mesh.nbPrisms, 4) = IDS(NODES_ELEM(4));
177 	mesh.PRISMS(mesh.nbPrisms, 5) = IDS(NODES_ELEM(5));
178 	mesh.PRISMS(mesh.nbPrisms, 6) = IDS(NODES_ELEM(6));
179 	mesh.PRISMS(mesh.nbPrisms, 7) = mesh.ELE_INFOS(I, tags);
181       if(mesh.ELE_INFOS(I, 2) == 7) %% pyramid
182 	mesh.nbPyramids = mesh.nbPyramids + 1;
183 	mesh.PYRAMIDS(mesh.nbPyramids, 1) = IDS(NODES_ELEM(1));
184 	mesh.PYRAMIDS(mesh.nbPyramids, 2) = IDS(NODES_ELEM(2));
185 	mesh.PYRAMIDS(mesh.nbPyramids, 3) = IDS(NODES_ELEM(3));
186 	mesh.PYRAMIDS(mesh.nbPyramids, 4) = IDS(NODES_ELEM(4));
187 	mesh.PYRAMIDS(mesh.nbPyramids, 5) = IDS(NODES_ELEM(5));
188 	mesh.PYRAMIDS(mesh.nbPyramids, 6) = IDS(NODES_ELEM(6));
189 	mesh.PYRAMIDS(mesh.nbPyramids, 7) = mesh.ELE_INFOS(I, tags);