- fixed remeshing of boundary layers when geo entities are constructed
using bnd layer points (no need to reload the file from scratch anymore
to avoid accumulating the transformations)

- patch load_gmsh.m for .msh version 2 (patch from Raoul)

- new load_gmsh2.m from the mailing list

Mesh/BoundaryLayer.cpp

-// $Id: BoundaryLayer.cpp,v 1.5 2007-09-10 04:47:04 geuzaine Exp $
+// $Id: BoundaryLayer.cpp,v 1.6 2007-10-04 13:04:37 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -47,14 +47,12 @@ static void addExtrudeNormals(std::vector<T*> &elements)
 
 int Mesh2DWithBoundaryLayers(GModel *m)
 {
-  ExtrudeParams *ep = 0;
-
   std::set<GFace*> sourceFaces, otherFaces;
   std::set<GEdge*> sourceEdges, otherEdges;
   for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
     GFace *gf = *it;
     if(gf->geomType() == GEntity::BoundaryLayerSurface){
-      ep = gf->meshAttributes.extrude;
+      ExtrudeParams *ep = gf->meshAttributes.extrude;
       if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == COPIED_ENTITY){
 	GFace *from = m->faceByTag(std::abs(ep->geo.Source));
 	if(!from){
@@ -96,17 +94,29 @@ int Mesh2DWithBoundaryLayers(GModel *m)
   ExtrudeParams::normals->normalize();
 
   // set the position of boundary layer points using the smooth normal
-  // field
-  for(GModel::viter it = m->firstVertex(); it != m->lastVertex(); it++){
-    GVertex *gv = *it;
-    if(gv->geomType() == GEntity::BoundaryLayerPoint){
-      GPoint p = gv->point();
-      ep->Extrude(ep->mesh.NbLayer - 1, ep->mesh.NbElmLayer[ep->mesh.NbLayer - 1],
-		  p.x(), p.y(), p.z());
-      gv->setPosition(p);
+  // field 
+  for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++){
+    GEdge *ge = *it;
+    if(ge->geomType() == GEntity::BoundaryLayerCurve){
+      ExtrudeParams *ep = ge->meshAttributes.extrude;
+      if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == EXTRUDED_ENTITY){
+	GVertex *vsrc, *vdest;
+	if(ge->getBeginVertex()->geomType() == GEntity::BoundaryLayerPoint){
+	  vsrc = ge->getEndVertex();
+	  vdest = ge->getBeginVertex();
+	}
+	else{
+	  vsrc = ge->getBeginVertex();
+	  vdest = ge->getEndVertex();
+	}
+	GPoint p = vsrc->point();
+	ep->Extrude(ep->mesh.NbLayer - 1, ep->mesh.NbElmLayer[ep->mesh.NbLayer - 1],
+		    p.x(), p.y(), p.z());
+	vdest->setPosition(p);
+      }
     }
   }
-
+  
   // remesh non-source edges (since they might have been modified by
   // the change in boundary layer points)
   std::for_each(otherFaces.begin(), otherFaces.end(), deMeshGFace());

utils/converters/matlab/load_gmsh.m

 function mesh = load_gmsh(filename)
-% Reads a mesh in msh format, version 1
+% Reads a mesh in msh format, version 1 or 2
+
+% Copyright (C) 10/2007 R Lorphvre (r.lorphevre@ulg.ac.be) 
+
+% Based on load_gmsh supplied with gmsh-2.0 and load_gmsh2 from JP
+% Moitinho de Almeida (moitinho@civil.ist.utl.pt)
+
+% number of nodes in function of the element type
+msh.NODES_PER_TYPE_OF_ELEMENT = [ 
+    2  3  4  4  8  6  5  3  6  9 10 27 18 14  1  8 20 15 13 ];
+
+% The format 2 don't sort the elements by reg phys but by
+% reg-elm. If this classification is important for your program,
+% use this (after calling this function):
+%
+% [OldRowNumber, NewRowNumber] = sort(OldMatrix(:,SortColumn)); 
+% NewMatrix = OldMatrix(NewRowNumber,:);
+%
+% Change the name of OldMatrix and NewMatrix with the name of yours
+% SortColumn by the number of the last column
+
 mesh = [];    
 mesh.MIN = zeros(3, 1);
 mesh.MAX = zeros(3, 1);
 fid = fopen(filename, 'r');
+disp (' ')
 while 1
   endoffile = 0;
   while 1
     tline = fgetl(fid);
     if feof(fid), endoffile=1, break, end
-    if tline(1) == '$', break, end  
+    if (tline(1) == '$' )
+      if tline(2) == 'N' && tline(3) == 'O'
+        fileformat = 1 ;
+        break  
+      end
+      if tline(2) == 'M' && tline(3) == 'e'
+        fileformat = 2;
+        tline = fgetl(fid);
+        disp('Mesh Type') 
+        disp (tline)
+        tline = fgetl(fid);
+        if (tline(1) == '$' && tline(2) == 'E'&& tline(3) == 'n')
+          tline = fgetl(fid);
+          break
+        else
+          disp (' This program can only read ASCII mesh file')
+          disp (' of format 1 or 2 from GMSH, try again?')
+          endoffile=1
+          break
+        end    
+      end
+      if tline(2) == 'E' && (tline(3) == 'L' || tline(3) == 'l' )
+        break  
+      end
+    end
   end
   if endoffile == 1, break, end
-  if tline(2) == 'N' && tline(3) == 'O'
+  if tline(2) == 'N' && ( tline(3) == 'O' || tline(3) == 'o' )
     disp('reading nodes')
     mesh.nbNod = fscanf(fid, '%d', 1);
     mesh.POS = zeros(mesh.nbNod, 3);
@@ -37,10 +82,18 @@ while 1
     end
     tline = fgetl(fid);
     disp('nodes have been read')
-  elseif tline(2) == 'E' && tline(3) == 'L'
+  elseif tline(2) == 'E' && ( tline(3) == 'L' || tline(3) == 'l')
     disp('reading elements')
     mesh.nbElm = fscanf(fid, '%d', 1);
-    mesh.ELE_INFOS = zeros(mesh.nbElm, 5);
+    if (fileformat == 1)
+        nbinfo = 5;
+        tags = 3;
+    end
+    if (fileformat == 2)
+        nbinfo = 4;
+        tags = 4;
+    end   
+    mesh.ELE_INFOS = zeros(mesh.nbElm, nbinfo); 
     mesh.nbPoints = 0;
     mesh.nbLines = 0;
     mesh.nbTriangles = 0;
@@ -59,25 +112,33 @@ while 1
     mesh.PRISMS = zeros(mesh.nbElm, 7);
     mesh.PYRAMIDS = zeros(mesh.nbElm, 6);
     for(I = 1:mesh.nbElm)
-      mesh.ELE_INFOS(I, :) = fscanf(fid, '%d', 5);
-      NODES_ELEM = fscanf(fid, '%d', mesh.ELE_INFOS(I, 5));
+      mesh.ELE_INFOS(I, :) = fscanf(fid, '%d', nbinfo);
+      if (fileformat == 1)
+          % take the mesh.ELE_INFOS(I, 5) nodes of the element
+          NODES_ELEM = fscanf(fid, '%d', mesh.ELE_INFOS(I, 5));  
+      end
+      if (fileformat == 2)
+          mesh.ELE_TAGS(I,:) = fscanf(fid, '%d', (mesh.ELE_INFOS(I,3)-1));
+          % take the msh.NODES_PER_TYPE_OF_ELEMENT (mesh.ELE_INFOS(I, 2)) nodes of the element
+          NODES_ELEM = fscanf(fid, '%d', msh.NODES_PER_TYPE_OF_ELEMENT (mesh.ELE_INFOS(I, 2)) ); 
+      end
       if(mesh.ELE_INFOS(I, 2) == 15) %% point
 	mesh.nbPoints = mesh.nbPoints + 1;
-	mesh.POINTS(mesh.nbPoints, 1) = IDS(NODES_ELEM(1));
-	mesh.POINTS(mesh.nbPoints, 2) = mesh.ELE_INFOS(I, 3);
+    mesh.POINTS(mesh.nbPoints, 1) = IDS(NODES_ELEM(1));
+	mesh.POINTS(mesh.nbPoints, 2) = mesh.ELE_INFOS(I, tags);
       end
       if(mesh.ELE_INFOS(I, 2) == 1) %% line
 	mesh.nbLines = mesh.nbLines + 1;
 	mesh.LINES(mesh.nbLines, 1) = IDS(NODES_ELEM(1));
 	mesh.LINES(mesh.nbLines, 2) = IDS(NODES_ELEM(2));   
-	mesh.LINES(mesh.nbLines, 3) = mesh.ELE_INFOS(I, 3);
+	mesh.LINES(mesh.nbLines, 3) = mesh.ELE_INFOS(I, tags);
       end
       if(mesh.ELE_INFOS(I, 2) == 2) %% triangle
 	mesh.nbTriangles = mesh.nbTriangles + 1;
 	mesh.TRIANGLES(mesh.nbTriangles, 1) = IDS(NODES_ELEM(1));
 	mesh.TRIANGLES(mesh.nbTriangles, 2) = IDS(NODES_ELEM(2));
 	mesh.TRIANGLES(mesh.nbTriangles, 3) = IDS(NODES_ELEM(3));
-	mesh.TRIANGLES(mesh.nbTriangles, 4) = mesh.ELE_INFOS(I, 3);
+	mesh.TRIANGLES(mesh.nbTriangles, 4) = mesh.ELE_INFOS(I, tags);
       end
       if(mesh.ELE_INFOS(I, 2) == 3) %% quadrangle
 	mesh.nbQuads = mesh.nbQuads + 1;
@@ -85,7 +146,7 @@ while 1
 	mesh.QUADS(mesh.nbQuads, 2) = IDS(NODES_ELEM(2));
 	mesh.QUADS(mesh.nbQuads, 3) = IDS(NODES_ELEM(3));
 	mesh.QUADS(mesh.nbQuads, 4) = IDS(NODES_ELEM(4));
-	mesh.QUADS(mesh.nbQuads, 5) = mesh.ELE_INFOS(I, 3);
+	mesh.QUADS(mesh.nbQuads, 5) = mesh.ELE_INFOS(I, tags);
       end
       if(mesh.ELE_INFOS(I, 2) == 4) %% tetrahedron
 	mesh.nbTets = mesh.nbTets + 1;
@@ -93,7 +154,7 @@ while 1
 	mesh.TETS(mesh.nbTets, 2) = IDS(NODES_ELEM(2));
 	mesh.TETS(mesh.nbTets, 3) = IDS(NODES_ELEM(3));
 	mesh.TETS(mesh.nbTets, 4) = IDS(NODES_ELEM(4));
-	mesh.TETS(mesh.nbTets, 5) = mesh.ELE_INFOS(I, 3);
+	mesh.TETS(mesh.nbTets, 5) = mesh.ELE_INFOS(I, tags);
       end
       if(mesh.ELE_INFOS(I, 2) == 5) %% hexahedron
 	mesh.nbHexas = mesh.nbHexas + 1;
@@ -105,7 +166,7 @@ while 1
 	mesh.HEXAS(mesh.nbHexas, 6) = IDS(NODES_ELEM(6));
 	mesh.HEXAS(mesh.nbHexas, 7) = IDS(NODES_ELEM(7));
 	mesh.HEXAS(mesh.nbHexas, 8) = IDS(NODES_ELEM(8));
-	mesh.HEXAS(mesh.nbHexas, 9) = mesh.ELE_INFOS(I, 3);
+	mesh.HEXAS(mesh.nbHexas, 9) = mesh.ELE_INFOS(I, tags);
       end
       if(mesh.ELE_INFOS(I, 2) == 6) %% prism
 	mesh.nbPrisms = mesh.nbPrisms + 1;
@@ -115,7 +176,7 @@ while 1
 	mesh.PRISMS(mesh.nbPrisms, 4) = IDS(NODES_ELEM(4));
 	mesh.PRISMS(mesh.nbPrisms, 5) = IDS(NODES_ELEM(5));
 	mesh.PRISMS(mesh.nbPrisms, 6) = IDS(NODES_ELEM(6));
-	mesh.PRISMS(mesh.nbPrisms, 7) = mesh.ELE_INFOS(I, 3);
+	mesh.PRISMS(mesh.nbPrisms, 7) = mesh.ELE_INFOS(I, tags);
       end
       if(mesh.ELE_INFOS(I, 2) == 7) %% pyramid
 	mesh.nbPyramids = mesh.nbPyramids + 1;
@@ -125,11 +186,13 @@ while 1
 	mesh.PYRAMIDS(mesh.nbPyramids, 4) = IDS(NODES_ELEM(4));
 	mesh.PYRAMIDS(mesh.nbPyramids, 5) = IDS(NODES_ELEM(5));
 	mesh.PYRAMIDS(mesh.nbPyramids, 6) = IDS(NODES_ELEM(6));
-	mesh.PYRAMIDS(mesh.nbPyramids, 7) = mesh.ELE_INFOS(I, 3);
+	mesh.PYRAMIDS(mesh.nbPyramids, 7) = mesh.ELE_INFOS(I, tags);
       end
     end
     tline = fgetl(fid);
     disp('elements have been read')
+    
   end
 end
+
 fclose(fid);

utils/converters/matlab/load_gmsh2.m

+function msh = load_gmsh2(filename)
+
+%% Reads a mesh in msh format, version 1 or 2
+
+% Copyright (C) 2007  JP Moitinho de Almeida (moitinho@civil.ist.utl.pt)
+% based on load_gmsh.m supplied with gmsh-2.0
+
+% Structure msh has the following elements:
+
+% msh.NODES_PER_TYPE_OF_ELEMENT - Number of nodes for each type of element
+% 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)))
+
+msh.MIN = [ inf; inf; inf ];
+msh.MAX = [ -inf; -inf; -inf ];
+
+%
+% These definitions need to be updated when new elemens types are added to gmsh
+%
+%   0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  1 
+%   1  2  3  4  5  6  7  8  9  0  1  2  3  4  5  6  7  8  9 
+msh.NODES_PER_TYPE_OF_ELEMENT = [ 
+    2  3  4  4  8  6  5  3  6  9 10 27 18 14  1  8 20 15 13 ];
+%
+%       0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2
+%       1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4
+Dim = [ 1 2 2 3 3 3 3 1 2 2 3 3 3 3 0 3 3 3 3 ];
+
+ntypes = length(msh.NODES_PER_TYPE_OF_ELEMENT);
+
+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);
+    msh.POS = zeros(msh.nbNod, 3);
+    for I=1:msh.nbNod
+        % iNod =
+        fscanf(fid, '%d', 1);
+        X = fscanf(fid, '%g', 3);
+        % IDS(iNod) = I;                      % Is it safe to assume that the nodes are sequentially presented?
+        msh.MAX = max(msh.MAX, X);
+        msh.MIN = min(msh.MIN, X);
+        msh.POS(I,:) = X;
+    end
+    fgetl(fid); % End previous line
+    fgetl(fid); % Has to be $ENDNOD $EndNodes
+end
+
+%% Read elements
+tline = fgetl(fid);
+if strcmp(tline,'$ELM') || strcmp(tline, '$Elements')
+    msh.nbElm = fscanf(fid, '%d', 1);
+    msh.ELE_INFOS = zeros(msh.nbElm, 4); % 1 - id, 2 - type, 3 -tags, 4 - Dimension
+    msh.ELE_NODES = zeros(msh.nbElm,1); % i - Element, j - ElNodes
+    if eq(fileformat, 1)
+        ntags = 2;
+    else
+        ntags = 3; % just a prediction
+    end
+    msh.ELE_TAGS = zeros(msh.nbElm, ntags);
+    msh.nbType = zeros(ntypes,1);
+    for I = 1:msh.nbElm
+        if (fileformat == 2)
+            msh.ELE_INFOS(I, 1:3) = fscanf(fid, '%d', 3);
+            ntags = msh.ELE_INFOS(I,3);
+            msh.ELE_TAGS(I, 1:ntags) = fscanf(fid, '%d', ntags);
+        else
+            aux = fscanf(fid, '%d', 5);
+            msh.ELE_INFOS(I, 1:2) = aux(1:2);
+            msh.ELE_INFOS(I, 3) = 2;
+            msh.ELE_TAGS(I, :) = aux(3:4);
+        end
+        type = msh.ELE_INFOS(I, 2);
+        msh.nbType(type) = msh.nbType(type) + 1;
+        msh.ELE_INFOS(I, 4) = Dim(type);
+        nnodes = msh.NODES_PER_TYPE_OF_ELEMENT(type);
+        NODES_ELEM = fscanf(fid, '%d', nnodes);
+        msh.ELE_NODES(I, 1:nnodes) = NODES_ELEM;
+    end
+    fgetl(fid); % End previous line
+    fgetl(fid); % Has to be $ENDELM or $EndElements
+else
+    disp('Error reading elements');
+    fileformat = 0;
+end
+
+if (fileformat == 0)
+    msh = [];
+end
+
+fclose(fid);
+return;