octave_packages/fpl-1.2.0/fpl_dx_write_field.m

   1 ##  Copyright (C) 2006,2007,2008,2009,2010  Carlo de Falco, Massimiliano Culpo
   2 ##
   3 ##  This file is part of:
   4 ##         FPL - Fem PLotting package for octave
   5 ##
   6 ##  FPL is free software; you can redistribute it and/or modify
   7 ##  it under the terms of the GNU General Public License as published by
   8 ##  the Free Software Foundation; either version 2 of the License, or
   9 ##  (at your option) any later version.
  10 ##
  11 ##  FPL is distributed in the hope that it will be useful,
  12 ##  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 ##  GNU General Public License for more details.
  15 ##
  16 ##  You should have received a copy of the GNU General Public License
  17 ##  along with FPL; If not, see <http://www.gnu.org/licenses/>.
  18 ##
  19 ##  author: Carlo de Falco     <cdf _AT_ users.sourceforge.net>
  20 ##  author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net>
  21
  22 ## -*- texinfo -*-
  23 ## @deftypefn {Function File} {} fpl_dx_writefield (@var{basename}, @
  24 ## @var{mesh}, @var{u},  @var{attr_name}, @var{attr_rank},  @
  25 ## @var{attr_shape}, @var{endfile})
  26 ##
  27 ## Output data field in ASCII Open-DX format.
  28 ##
  29 ## @var{basename} is a string containing the base-name of the dx file where the
  30 ## data will be saved.
  31 ##
  32 ## @var{mesh} is a PDE-tool like mesh, like the ones generated by the
  33 ## "msh" package.
  34 ##
  35 ## @var{u} is the field to be saved. It should represent scalar, vector
  36 ## or tensor of doubles.
  37 ##
  38 ## @var{attr_name} is a descriptive name for the field @var{u}, while
  39 ## @var{attr_rank} is the rank of the field (0 for scalar, 1 for vector,
  40 ## etc.) and @var{attr_shape} is the number of components of the field
  41 ## (assumed 1 for scalar).
  42 ##
  43 ## @var{endfile} should be 0 if you want to add other variables to the
  44 ## same file, 1 otherwise.
  45 ##
  46 ## Notice that when appending fields to an already existing file:
  47 ##
  48 ## @itemize
  49 ## @item @var{mesh} will not be printed to @var{filename}, but it will
  50 ## be only used to determine if the field is piece-wise constant or
  51 ## piece-wise linear
  52 ## @item @var{u} is not checked for consistency against the @var{mesh}
  53 ## already printed in @var{filename}
  54 ## @end itemize
  55 ##
  56 ## Example 1 (wrong usage):
  57 ## @example
  58 ## <generate msh1/u1 msh2/u2 in some way>
  59 ## fpl_dx_write_field("field.dx",msh1,u1,"density",1,0,0);
  60 ## fpl_dx_write_field("field.dx",msh2,u2,"temperature",1,0,1);
  61 ## @end example
  62 ## generate a file that fails at OpenDX run-time.
  63 ##
  64 ## Example 2:
  65 ## @example
  66 ## <generate msh1 and two fields u1-u2 in some way>
  67 ## fpl_dx_write_field("field",msh1,u1,"density",1,0,0);
  68 ## fpl_dx_write_field("field",msh1,u2,"temperature",1,0,1);
  69 ## @end example
  70 ## will generate a valid OpenDX ASCII data file.
  71 ##
  72 ## @seealso{fpl_dx_write_series}
  73 ##
  74 ## @end deftypefn
  75
  76 function fpl_dx_write_field(basename,mesh,u,attr_name,attr_rank,attr_shape,endfile)
  77
  78   ## Check input
  79   if nargin!=7
  80     error("fpl_dx_write_field: wrong number of input");
  81   endif
  82
  83   if !ischar(basename)
  84     error("fpl_dx_write_field: basename should be a valid string");
  85   elseif !( isstruct(mesh) )
  86     error("fpl_dx_write_field: mesh should be a valid structure");
  87   elseif !ismatrix(u)
  88     error("fpl_dx_write_field: u should be a valid matrix");
  89   elseif !ischar(attr_name)
  90     error("fpl_dx_write_field: attr_name should be a valid string");
  91   elseif !isscalar(attr_rank)
  92     error("fpl_dx_write_field: attr_rank should be a valid scalar");
  93   elseif !isscalar(attr_shape)
  94     error("fpl_dx_write_field: attr_shape should be a valid scalar");
  95   elseif !isscalar(endfile)
  96     error("fpl_dx_write_field: endfile should be a valid scalar");
  97   endif
  98
  99   filename = [basename ".dx"];
 100
 101   if ! exist(filename,"file")
 102     ## If file does not exist, create it
 103     fid    = fopen (filename,"w");
 104     create = 1;
 105   else
 106     ## FIXME: the following should be performed in a cleaner way! Does a
 107     ## backward fgetl function exist?
 108
 109     ## If file exist, check if it was already closed
 110     fid = fopen (filename,"r");
 111     fseek(fid,-4,SEEK_END);
 112     tst = fgetl(fid);
 113     if strcmp(tst,"end")
 114       error("fpl_dx_write_field: file %s exist and was already closed",filename);
 115     endif
 116     fclose(fid);
 117     fid    = fopen(filename,"a");
 118     create = 0;
 119   endif
 120
 121   p   = mesh.p';
 122   dim = columns(p); # 2D or 3D
 123
 124   if dim == 2
 125     t = mesh.t(1:3,:)';
 126   elseif dim == 3
 127     t = mesh.t(1:4,:)';
 128   else
 129     error("fpl_dx_write_field: neither 2D triangle nor 3D tetrahedral mesh");
 130   endif
 131
 132   nnodes = rows(p);
 133   nelems = rows(t);
 134   ndatas = rows(u);
 135
 136   if ndatas == nnodes
 137     dep = "positions";
 138   elseif ndatas == nelems
 139     dep = "connections";
 140   else
 141     error("fpl_dx_write_field: neither position nor connection data type")
 142   endif
 143
 144   if create
 145     ## If the file has just been created, print mesh information
 146     print_grid(fid,dim,p,nnodes,t,nelems);
 147   endif
 148   ## Otherwise assume the mesh is consistent with the one in the file
 149   ## and print only field information
 150   print_data(fid,u,ndatas,dep,attr_name,attr_rank,attr_shape);
 151
 152   if(endfile)
 153     fprintf(fid,"\nend\n");
 154   endif
 155   fclose (fid);
 156
 157 endfunction
 158
 159 ## fprint a 2Dtrg or 3Dtet mesh
 160 function print_grid(fid,dim,p,nnodes,t,nelems)
 161
 162   fprintf(fid,"object ""pos""\n");
 163   fprintf(fid,"class array type float rank 1 shape %d items %d data follows",dim,nnodes);
 164
 165   for ii = 1:nnodes
 166     fprintf(fid,"\n");
 167     fprintf(fid,"    %1.7e",p(ii,:));
 168   endfor
 169
 170   ## In DX format nodes are
 171   ## numbered starting from zero,
 172   ## instead we want to number
 173   ## them starting from 1!
 174   ## Here we restore the DX
 175   ## format
 176   if (min(min(t))==1)
 177     t -= 1;
 178   elseif(min(min(t))~=0)
 179     error("fpl_dx_write_field: check triangle structure")
 180   endif
 181
 182   fprintf(fid,"\n\nobject ""con""\n");
 183   fprintf(fid,"class array type int rank 1 shape %d items %d data follows",dim+1,nelems);
 184   for ii = 1:nelems
 185     fprintf(fid,"\n");
 186     fprintf(fid,"      %d",t(ii,:));
 187   endfor
 188
 189   fprintf(fid,"\n");
 190   if dim == 2
 191     fprintf(fid,"attribute ""element type"" string ""triangles""\n");
 192   elseif dim == 3
 193     fprintf(fid,"\nattribute ""element type"" string ""tetrahedra""\n");
 194   endif
 195   fprintf(fid,"attribute ""ref"" string ""positions""\n\n");
 196
 197 endfunction
 198
 199 ## fprint data on a trg grid
 200 function print_data(fid,u,ndatas,dep,attr_name,attr_rank,attr_shape)
 201
 202   if ((attr_rank == 0) && (min(size(u))==1))
 203     fprintf(fid,"object ""%s.data""\n",attr_name);
 204     fprintf(fid,"class array type double rank 0 items %d data follows",ndatas);
 205     fprintf(fid,"\n %1.7e",u);
 206   else
 207     fprintf(fid,"object ""%s.data""\n",attr_name);
 208     fprintf(fid,"class array type double rank %d shape %d items %d data follows",attr_rank,attr_shape,ndatas);
 209     for ii=1:ndatas
 210       fprintf(fid,"\n");
 211       fprintf(fid,"    %1.7e",u(ii,:));
 212     endfor
 213   endif
 214
 215   fprintf(fid,"\n");
 216   fprintf(fid,"attribute ""dep"" string ""%s"" \n\n",dep);
 217   fprintf(fid,"object ""%s"" class field\n",attr_name);
 218   fprintf(fid,"component ""positions"" value ""pos""\n");
 219   fprintf(fid,"component ""connections"" value ""con""\n");
 220   fprintf(fid,"component ""data"" value ""%s.data""\n",attr_name);
 221   fprintf(fid,"\n");
 222 endfunction