##  Copyright (C) 2006,2007,2008,2009,2010  Carlo de Falco, Massimiliano Culpo
##
##  This file is part of:
##         FPL - Fem PLotting package for octave
##
##  FPL is free software; you can redistribute it and/or modify
##  it under the terms of the GNU General Public License as published by
##  the Free Software Foundation; either version 2 of the License, or
##  (at your option) any later version.
##
##  FPL is distributed in the hope that it will be useful,
##  but WITHOUT ANY WARRANTY; without even the implied warranty of
##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##  GNU General Public License for more details.
##
##  You should have received a copy of the GNU General Public License
##  along with FPL; If not, see <http://www.gnu.org/licenses/>.
##
##  author: Carlo de Falco     <cdf _AT_ users.sourceforge.net>
##  author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net>

## -*- texinfo -*-
## @deftypefn {Function File} {} fpl_vtk_write_field (@var{basename}, @
## @var{mesh}, @var{nodedata},  @var{celldata}, @var{endfile})
## 
## Output data field in serial XML-VTK UnstructuredGrid format.
##
## @var{basename} is a string containing the base-name of the (vtu) file
## where the data will be saved.
##
## @var{mesh} is a PDE-tool like mesh, like the ones generated by the
## "msh" package.
##
## @var{nodedata} and @var{celldata} are (Ndata x 2) cell arrays containing
## respectively <PointData> and <CellData> representing scalars or
## vectors:
## @itemize @minus
## @item @var{*data}@{:,1@} = variable data;
## @item @var{*data}@{:,2@} = variable names;
## @end itemize
##
## @var{endfile} should be 0 if you want to add other variables to the
## same file, 1 otherwise.
##
## Example:
## @example
## <generate msh1, node centered field nc1, cell centered field cc1>
## fpl_vtk_write_field("example",msh1,@{nc1, "temperature"@},@
## @{cc1, "density"@},0);
## <generate msh2, node centered field nc2>
## fpl_vtk_write_field("example",msh2,@{nc2, "temperature"@},@
## @{@},1);
## @end example
## will generate a valid XML-VTK UnstructuredGrid file.
##
## @seealso{fpl_dx_write_field, fpl_dx_write_series} 
##
## @end deftypefn

function fpl_vtk_write_field(basename,mesh,nodedata,celldata,endfile)

  ## Check input
  if nargin!=5
    error("fpl_vtk_write_field: wrong number of input");
  endif

  if !ischar(basename)
    error("fpl_vtk_write_field: basename should be a valid string");
  elseif !( isstruct(mesh) )
    error("fpl_vtk_write_field: mesh should be a valid structure");
  elseif !(iscell(nodedata) && iscell(celldata))
    error("fpl_vtk_write_field: nodedata and celldata should be a valid cells");
  endif

  filename = [basename ".vtu"];

  if ! exist(filename,"file")
    fid = fopen (filename,"w");
    ## Header
    fprintf (fid,"<?xml version=""1.0""?>\n");
    fprintf (fid,"<VTKFile type=""UnstructuredGrid"" version=""0.1"" byte_order=""LittleEndian"">\n");
    fprintf (fid,"<UnstructuredGrid>\n");
  else
    ## FIXME: the following should be performed in a cleaner way! Does a
    ## backward fgetl function exist?

    ## If file exist, check if it was already closed
    fid = fopen (filename,"r");
    fseek(fid,-10,SEEK_END);
    tst = fgetl(fid);
    if strcmp(tst,"</VTKFile>")
      error("fpl_vtk_write_field: file %s exist and was already closed",filename);
    endif
    fclose(fid);
    fid = fopen (filename,"a");
  endif    

  p   = mesh.p;
  dim = rows(p); # 2D or 3D

  if dim == 2
    t = mesh.t(1:3,:);
  elseif dim == 3
    t = mesh.t(1:4,:);
  else
    error("fpl_vtk_write_field: neither 2D triangle nor 3D tetrahedral mesh");    
  endif
  
  t -= 1;
  
  nnodes = columns(p);
  nelems = columns(t);

  ## Header for <Piece>
  fprintf (fid,"<Piece NumberOfPoints=""%d"" NumberOfCells=""%d"">\n",nnodes,nelems);
  ## Print grid
  print_grid(fid,dim,p,nnodes,t,nelems);
  ## Print PointData
  print_data_points(fid,nodedata,nnodes)
  print_cell_data  (fid,celldata,nelems)
  ## Footer for <Piece>
  fprintf (fid, "</Piece>\n");

  if (endfile)
    ## Footer
    fprintf (fid, "</UnstructuredGrid>\n");
    fprintf (fid, "</VTKFile>");
  endif

  fclose (fid);
	   
endfunction

## Print Points and Cells Data
function print_grid(fid,dim,p,nnodes,t,nelems)
  
  if dim == 2
    p      = [p; zeros(1,nnodes)];
    eltype = 5;
  else
    eltype = 10;
  endif
  
  ## VTK-Points (mesh nodes)
  fprintf (fid,"<Points>\n");
  fprintf (fid,"<DataArray type=""Float64"" Name=""Array"" NumberOfComponents=""3"" format=""ascii"">\n");  
  fprintf (fid,"%g %g %g\n", p);
  fprintf (fid,"</DataArray>\n");
  fprintf (fid,"</Points>\n");

  ## VTK-Cells (mesh elements)
  fprintf (fid, "<Cells>\n");
  fprintf (fid, "<DataArray type=""Int32"" Name=""connectivity"" format=""ascii"">\n");
  if dim == 2
    fprintf (fid, "%d %d %d \n", t);
  else
    fprintf (fid, "%d %d %d %d \n", t);
  endif
  fprintf (fid, "</DataArray>\n");
  fprintf (fid, "<DataArray type=""Int32"" Name=""offsets"" format=""ascii"">\n");
  fprintf (fid, "%d %d %d %d %d %d\n", (dim+1):(dim+1):((dim+1)*nelems));
  fprintf (fid, "</DataArray>\n");
  fprintf (fid, "<DataArray type=""Int32"" Name=""types"" format=""ascii"">\n");
  fprintf (fid, "%d %d %d %d %d %d\n", eltype*ones(nelems,1));
  fprintf (fid, "</DataArray>\n");
  fprintf (fid, "</Cells>\n");

endfunction

## Print DataPoints
function print_data_points(fid,nodedata,nnodes)
  
  ## # of data to print in 
  ## <PointData> field
  nvdata = size(nodedata,1);  
  
  if (nvdata)
    fprintf (fid, "<PointData>\n");   
    for ii = 1:nvdata    
      data     = nodedata{ii,1};
      dataname = nodedata{ii,2};
      nsamples = rows(data);
      ncomp    = columns(data);
      if nsamples != nnodes
	error("fpl_vtk_write_field: wrong number of samples in <PointData> ""%s""",dataname);
      endif
      fprintf (fid,"<DataArray type=""Float32"" Name=""%s"" ",dataname);
      fprintf (fid,"NumberOfComponents=""%d"" format=""ascii"">\n",ncomp);
      for jj = 1:nsamples
	fprintf (fid,"%g ",data(jj,:));
	fprintf (fid,"\n");
      endfor
      fprintf (fid,"</DataArray>\n"); 
    endfor
    fprintf (fid, "</PointData>\n");
  endif

endfunction

function print_cell_data(fid,celldata,nelems)
  
  ## # of data to print in 
  ## <CellData> field
  nvdata = size(celldata,1); 

  if (nvdata)
    fprintf (fid, "<CellData>\n");
    for ii = 1:nvdata
      data     = celldata{ii,1};
      dataname = celldata{ii,2};
      nsamples = rows(data);
      ncomp    = columns(data);
      if nsamples != nelems
	error("fpl_vtk_write_field: wrong number of samples in <CellData> ""%s""",dataname);
      endif
      fprintf (fid,"<DataArray type=""Float32"" Name=""%s"" ",dataname);
      fprintf (fid,"NumberOfComponents=""%d"" format=""ascii"">\n",ncomp);
      for jj = 1:nsamples
	fprintf (fid,"%g ",data(jj,:));
	fprintf (fid,"\n");
      endfor
      fprintf (fid,"</DataArray>\n");
    endfor
    fprintf (fid, "</CellData>\n"); 
  endif

endfunction