Changeset 17719

Timestamp:

04/14/14 13:34:20 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: importgmsh.m now compatible with 3d meshes

File:

• issm/trunk-jpl/src/m/contrib/dassflow/importgmsh.m (modified) (5 diffs)

issm/trunk-jpl/src/m/contrib/dassflow/importgmsh.m ¶

@@ -1 @@
-function md = importgmsh(filename)
+function md = importgmsh(filename,dim)
 
 %some checks
@@ -30 +30 @@
 x = A(2,:)';
 y = A(3,:)';
+z = A(4,:)';
+
 A=fscanf(fid,'%s',1);
 if ~strcmp(A,'$EndNodes'), 
@@ -43 +45 @@
 disp(['Number of elements: ' num2str(nbt) ]);
 counter = 0;
-index   = zeros(0,3);
-segments       = zeros(0,2);
-segmentmarkers = zeros(0,1);
+if (dim==2),
+        index   = zeros(0,3);
+        segments       = zeros(0,2);
+        segmentmarkers = zeros(0,1);
+elseif (dim==3),
+        index   = zeros(0,4);
+        segments       = zeros(0,3);
+        segmentmarkers = zeros(0,1);
+else
+        error('not supported');
+end
+
 while(counter<nbt);
         id = fscanf(fid,'%i',1);
         ty = fscanf(fid,'%i',1);
-        if(ty>2) error('Type not supported'); end
         nbf = fscanf(fid,'%i',1);
         flags = fscanf(fid,'%i',nbf);
 
-        if ty==1,
-                A=fscanf(fid,'%i %i',2);
-                segments(end+1,:)=A;
-                if(flags(1)==5 & flags(2)==3), segmentmarkers(end+1)=3; 
-                elseif(flags(1)==1 & flags(2)==4), segmentmarkers(end+1)=4;
-                elseif(flags(1)==2 & flags(2)==1), segmentmarkers(end+1)=1;
-                elseif(flags(1)==4 & flags(2)==2), segmentmarkers(end+1)=2;
-                else error(['flags ' num2str(flags') ' not supported']);
-                end
-        else
-                A=fscanf(fid,'%i %i %i',3);
-                index(end+1,:)=A;
+        switch(ty)
+                case 1, %segments
+                        A=fscanf(fid,'%i %i',2);
+                        if (dim==2),  %Actual element
+                                segments(end+1,:)=A;
+                                if    (flags(1)==5 & flags(2)==3), segmentmarkers(end+1)=3; 
+                                elseif(flags(1)==1 & flags(2)==4), segmentmarkers(end+1)=4;
+                                elseif(flags(1)==2 & flags(2)==1), segmentmarkers(end+1)=1;
+                                elseif(flags(1)==4 & flags(2)==2), segmentmarkers(end+1)=2;
+                                else error(['flags ' num2str(flags') ' not supported']);
+                                end
+                        else
+                                error('not supported');
+                        end
+                case 2, %tria
+                        A=fscanf(fid,'%i %i %i',3);
+                        if (dim==2), %Actual element
+                                index(end+1,:)=A;
+                        else         %Boundary element
+                                segments(end+1,:)=A;
+                                if    (flags(1)==3), segmentmarkers(end+1)=3; 
+                                elseif(flags(1)==4), segmentmarkers(end+1)=4;
+                                elseif(flags(1)==1), segmentmarkers(end+1)=1;
+                                elseif(flags(1)==2), segmentmarkers(end+1)=2;
+                                elseif(flags(1)==5), segmentmarkers(end+1)=5;
+                                elseif(flags(1)==6), segmentmarkers(end+1)=6;
+                                else error(['flags ' num2str(flags') ' not supported']);
+                                end
+                        end
+                case 4, %tetra
+                        A=fscanf(fid,'%i %i %i %i',4);
+                        if (dim==3), %Actual element
+                                index(end+1,:)=A;
+                        else
+                                error('not supported');
+                        end
+                otherwise,
+                        error(['Type ' num2str(ty) ' not supported']);
         end
         counter = counter + 1;
@@ -70 +106 @@
 
 %recreate segments
-nbs = size(segments,1);
-segments = [segments zeros(nbs,1)];
-for i=1:nbs,
-        E = find(sum(ismember(index,segments(i,:)),2)>1);
-        segments(i,3)=E;
+if dim==2,
+        nbs = size(segments,1);
+        segments = [segments zeros(nbs,1)];
+        for i=1:nbs,
+                E = find(sum(ismember(index,segments(i,:)),2)>1);
+                segments(i,3)=E;
+        end
+else
+        nbs = size(segments,1);
+        segments = [segments zeros(nbs,1)];
+        for i=1:nbs,
+                E = find(sum(ismember(index,segments(i,:)),2)>2);
+                segments(i,4)=E;
+        end
 end
 
@@ -81 +126 @@
 
 %Create model
-md=meshconvert(model,index,x,y);
-md.mesh=mesh2dvertical(md.mesh);
-md.mesh.segmentmarkers=segmentmarkers;
-md.mesh.segments=segments;
-md.mesh.vertexonbase=zeros(md.mesh.numberofvertices,1);
-md.mesh.vertexonbase(find(vertexflags(md.mesh,1)))=1;
-md.mesh.vertexonsurface=zeros(md.mesh.numberofvertices,1);
-md.mesh.vertexonsurface(find(vertexflags(md.mesh,3)))=1;
+if 0, %2d triangles
+        md=meshconvert(model,index,x,y);
+        md.mesh=mesh2dvertical(md.mesh);
+        md.mesh.segmentmarkers=segmentmarkers;
+        md.mesh.segments=segments;
+        md.mesh.vertexonbase=zeros(md.mesh.numberofvertices,1);
+        md.mesh.vertexonbase(find(vertexflags(md.mesh,1)))=1;
+        md.mesh.vertexonsurface=zeros(md.mesh.numberofvertices,1);
+        md.mesh.vertexonsurface(find(vertexflags(md.mesh,3)))=1;
+else
+        md=model();
+        md.mesh=mesh3dprisms();
+        md.mesh.x = x;
+        md.mesh.y = y;
+        md.mesh.z = z;
+        md.mesh.elements = index;
+        md.mesh.numberofelements=size(md.mesh.elements,1);
+        md.mesh.numberofvertices=length(md.mesh.x);
+        md.mesh.vertexonbase=zeros(md.mesh.numberofvertices,1);
+        md.mesh.vertexonbase(find(segments(find(segmentmarkers==1),1:3)))=1;
+        md.mesh.vertexonsurface=zeros(md.mesh.numberofvertices,1);
+        md.mesh.vertexonsurface(find(segments(find(segmentmarkers==3),1:3)))=1;
+end