diff --git a/Mesh/BoundaryLayer.cpp b/Mesh/BoundaryLayer.cpp
index 6bac73c7670a10e712ae298fc4b5faa9d5c216ee..a29694c8854596f3ee59e4e14c2cacd24e647181 100644
--- a/Mesh/BoundaryLayer.cpp
+++ b/Mesh/BoundaryLayer.cpp
@@ -1,4 +1,4 @@
-// $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());
diff --git a/utils/converters/matlab/load_gmsh.m b/utils/converters/matlab/load_gmsh.m
index f1a9aa7dc2131b0dd329b62b0785f9ee8939876e..bba09258b3396d5999c26842419858fb17dc9848 100644
--- a/utils/converters/matlab/load_gmsh.m
+++ b/utils/converters/matlab/load_gmsh.m
@@ -1,18 +1,63 @@
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 Lorph�vre (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);
diff --git a/utils/converters/matlab/load_gmsh2.m b/utils/converters/matlab/load_gmsh2.m
new file mode 100644
index 0000000000000000000000000000000000000000..09db4bb37e60bf685750bc2a9a75314fe3b22b22
--- /dev/null
+++ b/utils/converters/matlab/load_gmsh2.m
@@ -0,0 +1,148 @@
+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;