##  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 3 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_dx_writefield (@var{basename}, @
## @var{mesh}, @var{u},  @var{attr_name}, @var{attr_rank},  @
## @var{attr_shape}, @var{endfile})
##
## Output data field in ASCII Open-DX format.
##
## @var{basename} is a string containing the base-name of the dx file where the
## data will be saved.
##
## @var{mesh} is a PDE-tool like mesh, like the ones generated by the
## "msh" package.
##
## @var{u} is the field to be saved. It should represent scalar, vector
## or tensor of doubles.
##
## @var{attr_name} is a descriptive name for the field @var{u}, while
## @var{attr_rank} is the rank of the field (0 for scalar, 1 for vector,
## etc.) and @var{attr_shape} is the number of components of the field
## (assumed 1 for scalar).
##
## @var{endfile} should be 0 if you want to add other variables to the
## same file, 1 otherwise.
##
## Notice that when appending fields to an already existing file:
##
## @itemize
## @item @var{mesh} will not be printed to @var{filename}, but it will
## be only used to determine if the field is piece-wise constant or
## piece-wise linear
## @item @var{u} is not checked for consistency against the @var{mesh}
## already printed in @var{filename}
## @end itemize
##
## Example 1 (wrong usage):
## @example
## <generate msh1/u1 msh2/u2 in some way>
## fpl_dx_write_field("field.dx",msh1,u1,"density",1,0,0);
## fpl_dx_write_field("field.dx",msh2,u2,"temperature",1,0,1);
## @end example
## generate a file that fails at OpenDX run-time.
##
## Example 2:
## @example
## <generate msh1 and two fields u1-u2 in some way>
## fpl_dx_write_field("field",msh1,u1,"density",1,0,0);
## fpl_dx_write_field("field",msh1,u2,"temperature",1,0,1);
## @end example
## will generate a valid OpenDX ASCII data file.
##
## @seealso{fpl_dx_write_series} 
##
## @end deftypefn

function fpl_dx_write_field(basename,mesh,u,attr_name,attr_rank,attr_shape,endfile)

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

  if !ischar(basename)
    error("fpl_dx_write_field: basename should be a valid string");
  elseif !( isstruct(mesh) )
    error("fpl_dx_write_field: mesh should be a valid structure");
  elseif !ismatrix(u)
    error("fpl_dx_write_field: u should be a valid matrix");
  elseif !ischar(attr_name)
    error("fpl_dx_write_field: attr_name should be a valid string");
  elseif !isscalar(attr_rank)
    error("fpl_dx_write_field: attr_rank should be a valid scalar");
  elseif !isscalar(attr_shape)
    error("fpl_dx_write_field: attr_shape should be a valid scalar");
  elseif !isscalar(endfile)
    error("fpl_dx_write_field: endfile should be a valid scalar");
  endif

  filename = [basename ".dx"];
    
  if ! exist(filename,"file")
    ## If file does not exist, create it
    fid    = fopen (filename,"w");
    create = 1;
  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,-4,SEEK_END);
    tst = fgetl(fid);
    if strcmp(tst,"end")
      error("fpl_dx_write_field: file %s exist and was already closed",filename);
    endif
    fclose(fid);
    fid    = fopen(filename,"a");
    create = 0;
  endif
  
  p   = mesh.p';
  dim = columns(p); # 2D or 3D

  if dim == 2
    t = mesh.t(1:3,:)';
  elseif dim == 3
    t = mesh.t(1:4,:)';
  else
    error("fpl_dx_write_field: neither 2D triangle nor 3D tetrahedral mesh");    
  endif
  
  nnodes = rows(p);
  nelems = rows(t);
  ndatas = rows(u);

  if ndatas == nnodes
    dep = "positions";
  elseif ndatas == nelems
    dep = "connections";
  else
    error("fpl_dx_write_field: neither position nor connection data type")
  endif
  
  if create 
    ## If the file has just been created, print mesh information
    print_grid(fid,dim,p,nnodes,t,nelems);
  endif
  ## Otherwise assume the mesh is consistent with the one in the file
  ## and print only field information
  print_data(fid,u,ndatas,dep,attr_name,attr_rank,attr_shape);

  if(endfile)
    fprintf(fid,"\nend\n");
  endif
  fclose (fid);

endfunction

## fprint a 2Dtrg or 3Dtet mesh
function print_grid(fid,dim,p,nnodes,t,nelems)

  fprintf(fid,"object ""pos""\n");
  fprintf(fid,"class array type float rank 1 shape %d items %d data follows",dim,nnodes);

  for ii = 1:nnodes
    fprintf(fid,"\n");
    fprintf(fid,"    %1.7e",p(ii,:));
  endfor

  ## In DX format nodes are 
  ## numbered starting from zero,
  ## instead we want to number
  ## them starting from 1!
  ## Here we restore the DX
  ## format		  
  if (min(min(t))==1)
    t -= 1;
  elseif(min(min(t))~=0)
    error("fpl_dx_write_field: check triangle structure")
  endif                    

  fprintf(fid,"\n\nobject ""con""\n");
  fprintf(fid,"class array type int rank 1 shape %d items %d data follows",dim+1,nelems);
  for ii = 1:nelems
    fprintf(fid,"\n");
    fprintf(fid,"      %d",t(ii,:));
  endfor

  fprintf(fid,"\n");
  if dim == 2
    fprintf(fid,"attribute ""element type"" string ""triangles""\n");
  elseif dim == 3
    fprintf(fid,"\nattribute ""element type"" string ""tetrahedra""\n");
  endif
  fprintf(fid,"attribute ""ref"" string ""positions""\n\n");

endfunction

## fprint data on a trg grid
function print_data(fid,u,ndatas,dep,attr_name,attr_rank,attr_shape)
  
  if ((attr_rank == 0) && (min(size(u))==1))
    fprintf(fid,"object ""%s.data""\n",attr_name);
    fprintf(fid,"class array type double rank 0 items %d data follows",ndatas);
    fprintf(fid,"\n %1.7e",u);
  else
    fprintf(fid,"object ""%s.data""\n",attr_name);
    fprintf(fid,"class array type double rank %d shape %d items %d data follows",attr_rank,attr_shape,ndatas);
    for ii=1:ndatas
      fprintf(fid,"\n");
      fprintf(fid,"    %1.7e",u(ii,:));
    endfor    
  endif
  
  fprintf(fid,"\n");
  fprintf(fid,"attribute ""dep"" string ""%s"" \n\n",dep);
  fprintf(fid,"object ""%s"" class field\n",attr_name);
  fprintf(fid,"component ""positions"" value ""pos""\n");
  fprintf(fid,"component ""connections"" value ""con""\n");
  fprintf(fid,"component ""data"" value ""%s.data""\n",attr_name);
  fprintf(fid,"\n");
endfunction