function out=fe_gmsh(varargin);

% Attempt to write an interface to GMSH
%
% fe_gmsh('addline',model,sel);
%  Nodes are stored in model.Node
%  Other geometric info GM=stack_get(model,'geom','GMSH','getdata');
%
% mo2=fe_gmsh('write FileName.geo [-lc (length)] [-run options]',model)
% mode=fe_gmsh('read FileName.msh')   
%
% Information about GMSH can be found at
%       http://www.geuz.org/gmsh/

%       Etienne Balmes
%       Copyright (c) 2001-2005 by INRIA and SDTools, All Rights Reserved.
%       Use under OpenFEM trademark.html license and LGPL.txt library license
%       For revision information use fe_gmsh('cvs')


[CAM,Cam]=comstr(varargin{1},1);carg=2;

% ------------------------------------------------------------------------
if comstr(Cam,'add') [CAM,Cam] = comstr(CAM,4);

model=varargin{carg};carg=carg+1;

GM=stack_get(model,'geom','GMSH','getdata');
if isempty(GM)
 GM=struct('Line',[],'Circle',[],'LineLoop',[],'PlaneSurface',[]);
end

% provides center, one edge node, normal - - - - - - - - - - - - - - - - - -
if comstr(Cam,'fullcircle') 

 r2=varargin{carg};carg=carg+1;
 if ~isequal(size(r2),[3 3]); error('You must provide two nodes and a normal');end

 r3=r2(1,:)-r2(2,:);r4=cross(r3,r2(3,:));
 r2=r2([1 1 1 1 1],:)+[0 0 0;r3;r4;-r3;-r4];

 [model.Node,i1]=feutil('addnode',model.Node,r2);
 GM.circle(end+[1:4],1:3)=i1([2 1 3;3 1 4;4 1 5;5 1 2]);

% Define a surface between two circles - - - - - - - - - - - - - - - - - -
% provide : center, node_radius_1, node_radius_2, normal
elseif comstr(Cam,'disk')


 circle=stack_get(model,'geom','circle','getdata');
 line=stack_get(model,'geom','line','getdata');
 LineLoop=stack_get(model,'geom','LineLoop','getdata');

 r2=varargin{carg};carg=carg+1;
 if ~isequal(size(r2),[4 3]); error('You must provide 3 nodes and a normal');end

 r3=r2(2:3,:)-r2([1 1],:); % two radius
 for j1=1:3 
  r3(end+1,:)=cross(r2(4,:),r3(end-1,:));
  r3(end+1,:)=cross(r2(4,:),r3(end-1,:));
 end
 r3=r2(ones(size(r3,1),1),:)+r3;
 [model.Node,i1]=feutil('addnode',model.Node,[r2(1,:);r3]);
 circle(end+[1:8],1:3)=i1([2 1 4;4 1 6;6 1 8;8 1 2;3 1 5;5 1 7;7 1 9; ...
    9 1 3]);
 %line(end+[1:4],1:2)=i1([2 3;4 5;6 7;8 9]);
 LineLoop(end+[1:2],:)=size(circle,1)-8+[1 2 3 4;5 6 7 8];

 model=stack_set(model,'geom','circle',circle);
 model=stack_set(model,'geom','line',line);
 model=stack_set(model,'geom','LineLoop',LineLoop);
 out=model;


% Add a straightline between two nodes - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'line')

 r2=varargin{carg};carg=carg+1;

 if ischar(r2) % Line is specified as a selection
  elt=feutil(horzcat('selelt',r2),model);

  r3=[];r1=fe_gmsh('lineloops',elt);
  for j1=1:length(r1) % list of continuous lines
   r2=r1{j1}(:);    
   ind=size(GM.Line,1)+[1:length(r2)-1];
   GM.Line(ind,1:2)=[r2(1:end-1) r2([2:end])];
   if r2(1)==r2(end) % LoopLine
    GM.LineLoop(end+1,1:length(ind))=ind(:)';
    r3(end+1)=size(GM.LineLoop,1);
   end
  end % list of continuous lines
  if ~isempty(r3); GM.PlaneSurface(end+1,1:length(r3))=r3(:)';end
 elseif isequal(size(r2),[2 3]); % Line specified using two nodes
   [model.Node,i1]=feutil('addnode',model.Node,r2);
   GM.Line(end+1,1:2)=i1(:)';
 else; error('Not a valid line input format');
 end


% Add a straightline between two nodes - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'line')

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
else sdtw('''Add%s'' unknow',CAM); % subcommand selection - - - - - - - 
end
out=stack_set(model,'geom','GMSH',GM);

% ------------------------------------------------------------------------
elseif comstr(Cam,'lineloops')

elt=varargin{carg};carg=carg+1;

[ElemF,st2,ElemP]= feutil('getelemf',elt(1,:),1);
elt(find(~finite(elt(:,1))),:)=[];
i1=strfind(Cam,'-elt'); if ~isempty(i1); 
 RunOpt.out='elt';
else RunOpt.out='loops';
end
if strcmp(ElemP,'beam3') % allow for mid-node export
      RunOpt.edge=sort(elt(:,1:2),2); RunOpt.edge3=elt(:,1:3);
end

  j2=[];r1=cell(1,0);
  conn=sparse(elt(:,1:2),[1:size(elt,1)]'*[1 1],ones(size(elt,1),1)*[1 2]);
  while any(elt(:));
   if isempty(j2);  % new LineLoop
     j2=min(find(any(conn)));
     r1{end+1}=elt(j2,1:2);elt(j2,:)=0;conn(:,j2)=0;
   else  % append at the end
    r3=conn(r1{end}(1),:);j2=[];
    if any(r3); 
     [j1,j2,r3]=find(r3); r1{end}=[elt(j2,3-r3) r1{end}];
     conn(:,j2)=0;elt(j2,:)=0;
    end
    r3=conn(r1{end}(end),:);
    if any(r3); 
     [j1,j2,r3]=find(r3); r1{end}=[r1{end} elt(j2,3-r3)];
     conn(:,j2)=0;elt(j2,:)=0;
    end
   end
  end % while any elt

if strcmp(RunOpt.out,'elt') % return oriented elements
 out=[]; for j1=1:length(r1)
  out(end+1,1:length(ElemP)+1)=[Inf abs(ElemP)];
  r2=[r1{j1}(1:end-1)' r1{j1}(2:end)'];
  if strcmp(ElemP,'beam3')
    [r3,i3,i4]=intersect(sort(r2,2),RunOpt.edge,'rows');
    i2=[];i2(i3)=1:length(i3);isequal(sort(r3(i2,:),2),sort(r2,2))
    r3=RunOpt.edge3(i4,:);r3=r3(i2,:);r3(:,1:2)=r2; r2=r3;
  end
  out(end+[1:size(r2,1)],1:size(r2,2))=r2;
 end
else
 out=r1;
end

% ------------------------------------------------------------------------
elseif comstr(Cam,'write'); [CAM,Cam] = comstr(CAM,6);


model=varargin{carg};carg=carg+1;
GM=stack_get(model,'geom','GMSH','getdata')
RunOpt=[];

lc=[];i1=findstr(Cam,'-lc'); if ~isempty(i1)
 [lc,st,Cam] = comstr(CAM(i1:end),'-lc','%g');
 [CAM,Cam]=comstr([CAM(1:i1-1) st],1);
end
i1=findstr(Cam,'-run'); if ~isempty(i1);
 RunOpt=struct('Run',CAM(i1+4:end));[CAM,Cam]=comstr([CAM(1:i1-1)],1);
end

if isempty(lc); lc=getpref('OpenFEM','DefaultElementLength',.1);end
if isempty(Cam);fid=1;else;fname=CAM; fid=fopen(fname,'w');end

ind=find(model.Node(:,4)==0);if ~isempty(ind); model.Node(ind,4)=lc;end
% Write Nodes

fprintf(fid,'Point(%i)={%.10g,%.10g,%.10g,%.10g};\n', ...
 model.Node(:,[1 5 6 7 4])');

% Write Circles
i0=0;
r1=GM.Circle;
if ~isempty(r1)
 r1=[i0+[1:size(r1,1)]' r1]; i0=r1(end,1);
 fprintf(fid,'Circle(%i)={%i,%i,%i};\n',r1');
end
% Write Lines
r1=GM.Line;
if ~isempty(r1)
 r1=[i0+[1:size(r1,1)]' r1];i0=r1(end,1);
 fprintf(fid,'Line(%i)={%i,%i};\n',r1');
end
% Write LineLoops
r1=GM.LineLoop;GM.iLineLoop=i0;
if ~isempty(r1)
 r1=[i0+[1:size(r1,1)]' r1];i0=r1(end,1);
 for j1=1:size(r1,1)
  i2=r1(j1,find(r1(j1,:)));
  fprintf(fid,'Line Loop(%i)={%i%s};\n',i2(1:2),sprintf(',%i',i2(3:end)));
 end
end
% Write PlaneSurfaces
r1=GM.PlaneSurface;
if ~isempty(r1)
 r1=[i0+[1:size(r1,1)]' r1+GM.iLineLoop];i0=r1(end,1);
 for j1=1:size(r1,1)
  i2=r1(j1,:);i2=i2(find(i2));
  if length(i2)==2;st='';else;st=sprintf(',%i',i2(3:end));end
  fprintf(fid,'Plane Surface(%i)={%i%s};\n',i2(1:2),st);
 end
end

if fid~=1; fclose(fid);end
if isfield(RunOpt,'Run'); 
 RunOpt.Run=sprintf('!"%s" "%s" %s', ...
   getpref('OpenFEM','gmsh','gmsh.exe'),fname,RunOpt.Run);
 eval(RunOpt.Run);
 [wd,fname,ext]=fileparts(fname);if isempty(wd);wd=pwd;end
 out=fe_gmsh(sprintf('read %s',fullfile(wd,horzcat(fname,'.msh'))));
end

% ------------------------------------------------------------------------
elseif comstr(Cam,'read') [CAM,Cam] = comstr(CAM,5);

[wd,fname,ext]=fileparts(CAM);
fid=fopen(CAM,'r');

if comstr(ext,'.msh') % Reading .MSH files

st=fgetl(fid);
model=struct('Node',[],'Elt',[]);

while isstr(st);

st1=sscanf(st,'%s',1);
switch st1
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
case '$NOD' % Reading nodes

 i1=fscanf(fid,'%i',1); % number of nodes
 r1=fscanf(fid,'%i %g %g %g',[4 i1])';
 r1(:,5)=0;
 model.Node=r1(:,[1 5 5 5 2 3 4]); 

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
case '$ELM' % Reading elements
 %elm-number elm-type reg-phys reg-elem number-of-nodes node-number-list
 ne=fscanf(fid,'%i',1); % number of elements
 i1=fscanf(fid,'%i');
 while ~isempty(i1)
  switch i1(2)
  case 1;  ElemF='beam1';   i3=7;   ind=[1:2];
  case 2;  ElemF='tria3';   i3=8;   ind=[1:3];
  case 3;  ElemF='quad4';   i3=9;   ind=[1:4];
  case 4;  ElemF='tetra4';  i3=9;   ind=[1:4];
  case 5;  ElemF='hexa8';   i3=13;  ind=[1:8];
  case 6;  ElemF='penta6';  i3=11;  ind=[1:6];
  case 7;  ElemF='pyra5';   i3=10;  ind=[1:5];
  case 8;  ElemF='beam1';   i3=8;   ind=[1:2];
  case 9;  ElemF='tria6';   i3=11;  ind=[1:6];
  case 10; ElemF='quad9';   i3=14;  ind=[1:9];
  case 11; ElemF='tetra10'; i3=15;  ind=[1:10];
  case 12; ElemF='hexa20';  i3=25;  ind=[1:20];
  case 13; ElemF='penta15'; i3=20;  ind=[1:15];
  case 14; ElemF='pyra13';  i3=18;  ind=[1:13];
  case 15; ElemF='mass2';   i3=6;   ind=[1];
  otherwise error(1)
  end

  i2=floor(length(i1)/i3);i2=reshape(i1(1:i2*i3),i3,i2)';
  i3=find(diff(i2(:,2))~=0); if ~isempty(i3);i2=i2(1:min(i3),:); end
  i1(1:prod(size(i2)))=[];
  model.Elt(end+1,1:length(ElemF)+1)=[Inf abs(ElemF)];
  model.Elt(end+[1:size(i2,1)],1:size(i2,2)-3)=i2(:,[6:end 3 4]);
 end



% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
otherwise; 
  if ~comstr(st1,'$END')&~isempty(st1); fprintf('''%s'' ignored\n',st1);end
end
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
st=fgetl(fid);

end % loop on file

else sdtw('''%s'' unknow format',ext); % subcommand selection - - - - - - - 
end

if fid~=1; fclose(fid);end
if nargout==0; feplot(model);else; out=model;end

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'cvs')

 out='$Revision: 1.8 $  $Date: 2006/01/27 14:28:50 $';

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
else ; sdtw('''%s'' unknow',CAM); % subcommand selection - - - - - - - 
end % function
