function [] = opfem2VTK(varargin)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  convert openFEM format to VTK
%
%  opfem2VTK(rootvtkfile,model)
%
%  with following mandatory fields:
%
%     model.Node
%     model.Elt
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rootvtkfile=varargin{1};
model = varargin{2};
Nodes = model.Node(:,5:7);
elt = model.Elt;
ig = find(elt(:,1)==Inf);
ng=length(ig);
ig=[ig;size(elt,1)+1];

  % NEW FILE
  obj = [rootvtkfile,'.vtk'];
  fidend = fopen(obj,'w');
  
  % HEADER  
  fprintf(fidend,['# vtk DataFile Version 2.0','\n']);
  fprintf(fidend,'%s\n',rootvtkfile);
  fprintf(fidend,['ASCII','\n']);
  fprintf(fidend,['DATASET UNSTRUCTURED_GRID','\n']);
  
  % POINTS
  %disp('Write Points')
  fprintf(fidend,'POINTS ');
  fprintf(fidend,['%i',' '],size(Nodes,1));
  fprintf(fidend,['float','\n']);
  fprintf(fidend,'%13.6f %13.6f %13.6f\n',Nodes');
  
  
  cell_type=zeros(ng,1);
  cell_type=zeros(ng,1);
  ne=zeros(ng,1);
  elemF=cell(ng,1);
  T=cell(ng,1);
  n_cell=0;
  size_cell = 0;
  % ELEMENTS IDENTIFICATION
  for k=1:ng
    elemF{k} = char(elt(ig(k),2:end));
    ne(k)=size(elt(ig(k)+1:ig(k+1)-1,:),1);
    n_cell = n_cell+ne(k);
    if (comstr(elemF{k},'tria3') | comstr(elemF{k},'t3p') | comstr(elemF{k},'dktp'))
        size_cell = size_cell + ne(k)*4;
        T{k} = [3*ones(ne(k),1) elt(ig(k)+1:ig(k+1)-1,1:3)-1];
        cell_type(k)=5;
    elseif (comstr(elemF{k},'quad4') | comstr(elemF{k},'q4p') | comstr(elemF{k},'mitc4'))
        size_cell = size_cell + ne(k)*5;
        T{k} = [4*ones(ne(k),1) elt(ig(k)+1:ig(k+1)-1,1:4)-1];
        cell_type(k)=9;
    elseif comstr(elemF{k},'tetra4')
        size_cell = size_cell + ne(k)*5;
        T{k} = [4*ones(ne(k),1) elt(ig(k)+1:ig(k+1)-1,1:4)-1];
        cell_type(k)=10;
    elseif comstr(elemF{k},'hexa8')
        size_cell = size_cell + ne(k)*9;
        T{k} = [8*ones(ne(k),1) elt(ig(k)+1:ig(k+1)-1,1:8)-1];
        cell_type(k)=12;
    elseif comstr(elemF{k},'penta6')
        size_cell = size_cell + ne(k)*7;
        T{k} = [6*ones(ne(k),1) elt(ig(k)+1:ig(k+1)-1,1:6)-1];
        cell_type(k)=13;
    elseif comstr(elemF{k},'mitc6')
        size_cell = size_cell + ne(k)*7;
        T{k} = [6*ones(ne(k),1) elt(ig(k)+1:ig(k+1)-1,1:6)-1];
        cell_type(k)=22;
    elseif (comstr(elemF{k},'q8p') | comstr(elemF{k},'quadb'))
        size_cell = size_cell + ne(k)*9;
        T{k} = [8*ones(ne(k),1) elt(ig(k)+1:ig(k+1)-1,1:8)-1];
        cell_type(k)=23;
    elseif comstr(elemF{k},'tetra10')
        size_cell = size_cell + ne(k)*11;
        T{k} = [10*ones(ne(k),1) elt(ig(k)+1:ig(k+1)-1,1:10)-1];
        cell_type(k)=24;
    elseif comstr(elemF{k},'hexa20')
        size_cell = size_cell + ne(k)*21;
        T{k} = [20*ones(ne(k),1) elt(ig(k)+1:ig(k+1)-1,1:12)-1 elt(ig(k)+1:ig(k+1)-1,17:20)-1 elt(ig(k)+1:ig(k+1)-1,13:16)-1];
        cell_type(k)=25;
    else
        error(sprintf('Sorry, %s not yet implemented :-) !',elemF{k}))
    end
  end
  
  ref = [];
  fprintf(fidend,'CELLS ');
  fprintf(fidend,['%i',' '],n_cell);
  fprintf(fidend,['%i','\n'],size_cell);
  
  % ELEMENTS WRITING
  for k=1:ng
	if cell_type(k)==5
		fprintf(fidend,'%i %6i %6i %6i\n',T{k}');
	elseif (cell_type(k)==9 |cell_type(k)==10)
		fprintf(fidend,'%i %6i %6i %6i %6i\n',T{k}');
  	elseif (cell_type(k)==12 | cell_type(k)==23)
		fprintf(fidend,'%i %6i %6i %6i %6i %6i %6i %6i %6i\n',T{k}');
        elseif (cell_type(k)==13 | cell_type(k)==22)
                fprintf(fidend,'%i %6i %6i %6i %6i %6i %6i\n',T{k}');
  	elseif cell_type(k)==24
		fprintf(fidend,'%i %6i %6i %6i %6i %6i %6i %6i %6i %6i %6i\n',T{k}');
  	elseif cell_type(k)==25
		fprintf(fidend,'%i %6i %6i %6i %6i %6i %6i %6i %6i %6i %6i %6i %6i %6i %6i %6i %6i %6i %6i %6i %6i\n',T{k}');
	end
  end	

  fprintf(fidend,'CELL_TYPES ');
  fprintf(fidend,'%i\n',n_cell);
  
  for k=1:ng
    L = cell_type(k)*ones(1,ne(k));
    %disp('Write Cells types')
    fprintf(fidend,'%i\n', L');
  end

if nargin > 2
    nodevar =[];
    eltvar=[];
    for ivar=3:nargin
        if size(varargin{ivar}.Data,1) == size(Nodes,1)
            nodevar = [nodevar ivar];
        elseif size(varargin{ivar}.Data,1) == n_cell
            eltvar = [eltvar ivar];
        end
    end
    if length(nodevar)
        nodearg =[];
        for inod=1:length(nodevar)
            nodearg = [nodearg ',varargin{nodevar(' num2str(inod) ')}'];
        end
        eval(['writeVTKnodedata(fidend' nodearg ')']);
    end
    if length(eltvar)
        eltarg = [];
        for ielt=1:length(eltvar)
            eltarg = [eltarg ',varargin{eltvar(' num2str(ielt) ')}'];
        end 
        eval(['writeVTKeltdata(fidend' eltarg ')']);
    end
  
end

fclose('all');

% ----------------------------------------------------------------

function [] = writeVTKnodedata(varargin);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    append data at node to an existant vtk file
%
%    io_writeVTKnodedata(vtkrootfile,var1,var2,...varn)
%
%    with var_i as follow:
%
%            var_i.Name = 'namegiventodata'
%            var_i.Data = datavalues   
%                    nnode rows & one column for scalars
%                    nnode rows & three columns for vectors 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%rootMeshFile = varargin{1};

  % APPEND FILE
  %obj = [rootMeshFile,'.vtk'];
  %fidend = fopen(obj,'a');
  fidend = varargin{1};


 fprintf(fidend,'\nPOINT_DATA ');
 fprintf(fidend,'%i\n',size(varargin{2}.Data,1));


for k=2:nargin
  if size(varargin{k}.Data,2)>1
     fprintf(fidend,'\nVECTORS %s float\n',varargin{k}.Name);
     fprintf(fidend,'%e %e %e\n',varargin{k}.Data');
  else
     fprintf(fidend,'SCALARS %s float\n',varargin{k}.Name);
     fprintf(fidend,'LOOKUP_TABLE default\n');
     fprintf(fidend,'%e\n',varargin{k}.Data');
  end
end

%fclose('all');


% -------------------------------------------------------------------

function [] = writeVTKeltdata(varargin);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    append data at node to an existant vtk file
%
%    io_writeVTKeltdata(vtkrootfile,var1,var2,...varn)
%
%    with var_i as follow:
%
%            var_i.Name = 'namegiventodata'
%            var_i.Data = datavalues   
%                    nelt rows & one column for scalars
%                    nelt rows & three columns for vectors 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%rootMeshFile = varargin{1};

  % APPEND FILE
  %obj = [rootMeshFile,'.vtk'];
  %fidend = fopen(obj,'a');
  fidend = varargin{1};


 fprintf(fidend,'\nCELL_DATA ');
 fprintf(fidend,'%i\n',size(varargin{2}.Data,1));


for k=2:nargin
  if size(varargin{k}.Data,2)>1
     fprintf(fidend,'\nVECTORS %s float\n',varargin{k}.Name);
     fprintf(fidend,'%e %e %e\n',varargin{k}.Data');
  else
     fprintf(fidend,'SCALARS %s float\n',varargin{k}.Name);
     fprintf(fidend,'LOOKUP_TABLE default\n');
     fprintf(fidend,'%e\n',varargin{k}.Data');
  end
end

%fclose('all');
