function md = importgmsh(filename)

%some checks
if ~exist(filename),
	error(['expread error message: file ' filename ' not found!']);
end

%open file
fid=fopen(filename,'r');

%Get Mesh format
A=fscanf(fid,'%s',1);
if ~strcmp(A,'$MeshFormat'), 
	error(['Expecting $MeshFormat (' A ')']);
end
A=fscanf(fid,'%f %i %i',[1 3]);
A=fscanf(fid,'%s',1);
if ~strcmp(A,'$EndMeshFormat'), 
	error(['Expecting $EndMeshFormat (' A ')']);
end

%Nodes
A=fscanf(fid,'%s',1);
if ~strcmp(A,'$Nodes'), 
	error(['Expecting $Nodes (' A ')']);
end
nbv=fscanf(fid,'%i',1);
disp(['Number of nodes: ' num2str(nbv) ]);
A=fscanf(fid,'%i %f %f %f',[4 nbv]);
x = A(2,:)';
y = A(3,:)';
A=fscanf(fid,'%s',1);
if ~strcmp(A,'$EndNodes'), 
	error(['Expecting $EndNodes (' A ')']);
end

%Elements
A=fscanf(fid,'%s',1);
if ~strcmp(A,'$Elements'), 
	error(['Expecting $Elements (' A ')']);
end
nbt=fscanf(fid,'%i',1);
disp(['Number of elements: ' num2str(nbt) ]);
counter = 0;
index   = zeros(0,3);
segments       = zeros(0,2);
segmentmarkers = zeros(0,1);
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;
	end
	counter = counter + 1;
end

%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;
end

%close file
fclose(fid);

%Create model
md=meshconvert(model,index,x,y);
md.mesh=mesh2dvertical(md.mesh);
md.mesh.segmentmarkers=segmentmarkers;
md.mesh.segments=segments;
md.mesh.vertexonbed=zeros(md.mesh.numberofvertices,1);
md.mesh.vertexonbed(find(vertexflags(md.mesh,1)))=1;
md.mesh.vertexonsurface=zeros(md.mesh.numberofvertices,1);
md.mesh.vertexonsurface(find(vertexflags(md.mesh,3)))=1;
