octave_packages/fpl-1.2.0/FPL2vtkoutputdata.m

   1 ## Copyright (C) 2008 Carlo de Falco
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 2 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; if not, write to the Free Software
  15 ## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  16
  17
  18 ## -*- texinfo -*-
  19 ## @deftypefn {Function File} {} FPL2vtkoutputdata ( @var{filename}, @var{p}, @var{t}, @var{nodedata}, @var{celldata}, @var{header}, @var{vtkver})
  20 ##
  21 ## Save data in VTK ASCII format.
  22 ##
  23 ## @itemize @minus
  24 ##  @item  @var{filename} = name of file to save into
  25 ##  @item  @var{p}, @var{t} = mesh node coordinates and connectivity
  26 ##  @item  @var{name} = name of a mesh variable
  27 ##  @item  @var{nodedata}/@var{celldata} = node/cell centered data
  28 ##  fields (2xNfields cell array), @var{*data}@{:,1@} = variable names;
  29 ##  @var{*data}@{:,2@} = variable data;
  30 ##  @item @var{header} comment to add in the file header
  31 ##  @item @var{vtkver} format version (default is 3.0)
  32 ## @end itemize
  33 ##
  34 ## @seealso{FPL2dxoutputdata}
  35 ## @end deftypefn
  36
  37 function FPL2vtkoutputdata (filename, p, t, nodedata, celldata, header, vtkver)
  38
  39   fid = fopen (filename, "w");
  40   if ( fid )
  41
  42     ## version ID
  43     if (!exist("vtkver"))
  44       vtkver = [3 0];
  45     endif
  46
  47     fprintf (fid, "# vtk DataFile Version %d.%d\n", vtkver(1), vtkver(2));
  48
  49     ## header
  50     if (!exist("header"))
  51       header = "";
  52     elseif (length(header) > 255)
  53       header (255:end) = [];
  54     endif
  55
  56     fprintf (fid, "%s\n", header);
  57
  58     ## File format: only ASCII supported for the moment
  59     fprintf (fid, "ASCII\n");
  60
  61     ## Mesh: only triangles suported
  62     fprintf (fid, "DATASET UNSTRUCTURED_GRID\n");
  63
  64     Nnodes = columns(p);
  65     fprintf (fid,"POINTS %d double\n", Nnodes);
  66     fprintf (fid,"%g %g 0\n", p);
  67
  68     Nelem = columns(t);
  69     T     = zeros(4,Nelem);
  70     T(1,:)= 3;
  71     T(2:4,:) = t(1:3,:) -1;
  72     fprintf (fid,"CELLS %d %d\n", Nelem, Nelem*4);
  73     fprintf (fid,"%d %d %d %d\n", T);
  74     fprintf (fid,"CELL_TYPES %d \n", Nelem);
  75     fprintf (fid,"%d\n", 5*ones(Nelem,1));
  76
  77     ## node data
  78     if (exist("nodedata"))
  79       nfields = rows(nodedata);
  80       if nfields
  81         fprintf (fid,"POINT_DATA %d\n", Nnodes);
  82         for ifield = 1:nfields
  83           V = nodedata {ifield, 2};
  84           N = nodedata {ifield, 1};
  85           if (isvector (V))
  86             fprintf (fid,"SCALARS %s double\nLOOKUP_TABLE default\n", N);
  87             fprintf (fid,"%g\n", V);
  88           else
  89             V(:,3) = 0;
  90             fprintf (fid,"VECTORS %s double\nLOOKUP_TABLE default\n", N);
  91             fprintf (fid,"%g %g %g\n", V);
  92           endif
  93         endfor
  94       endif
  95     endif
  96
  97     ## cell data
  98     if (exist("celldata"))
  99       nfields = rows(celldata);
 100       if nfields
 101         fprintf (fid,"CELL_DATA %d\n", Nelem);
 102         for ifield = 1:nfields
 103           V = celldata {ifield, 2};
 104           N = celldata {ifield, 1};
 105           if (isvector (V))
 106             fprintf (fid,"SCALARS %s double\nLOOKUP_TABLE default\n", N);
 107             fprintf (fid,"%g\n", V);
 108           else
 109             V(:,3) = 0;
 110             fprintf (fid,"VECTORS %s double\nLOOKUP_TABLE default\n", N);
 111             fprintf (fid,"%g %g %g\n", V);
 112           endif
 113         endfor
 114       endif
 115     endif
 116
 117     ## cleanup
 118     fclose (fid);
 119   else
 120     error(["FPL2vtkoutputdata: could not open file " filename]);
 121   endif
 122 endfunction
 123
 124 %!test
 125 %! msh.p =[ 0   0   0   0   1   1   1   1   2   2   2   2   3   3   3   3
 126 %!          0   1   2   3   0   1   2   3   0   1   2   3   0   1   2   3];
 127 %! msh.e =[1    5    9   13   14   15    4    8   12    1    2    3
 128 %!     5    9   13   14   15   16    8   12   16    2    3    4
 129 %!     0    0    0    0    0    0    0    0    0    0    0    0
 130 %!     0    0    0    0    0    0    0    0    0    0    0    0
 131 %!     1    1    1    2    2    2    3    3    3    4    4    4
 132 %!     0    0    0    0    0    0    0    0    0    0    0    0
 133 %!     1    1    1    1    1    1    1    1    1    1    1    1];
 134 %! msh.t =[ 1    2    3    5    6    7    9   10   11    1    2    3    5    6    7    9   10   11
 135 %!     5    6    7    9   10   11   13   14   15    6    7    8   10   11   12   14   15   16
 136 %!     6    7    8   10   11   12   14   15   16    2    3    4    6    7    8   10   11   12
 137 %!     1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1];
 138 %! testfile = "# vtk DataFile Version 3.0\n\nASCII\nDATASET UNSTRUCTURED_GRID\nPOINTS 16 double\n0 0 0\n0 1 0\n0 2 0\n0 3 0\n1 0 0\n1 1 0\n1 2 0\n1 3 0\n2 0 0\n2 1 0\n2 2 0\n2 3 0\n3 0 0\n3 1 0\n3 2 0\n3 3 0\nCELLS 18 72\n3 0 4 5\n3 1 5 6\n3 2 6 7\n3 4 8 9\n3 5 9 10\n3 6 10 11\n3 8 12 13\n3 9 13 14\n3 10 14 15\n3 0 5 1\n3 1 6 2\n3 2 7 3\n3 4 9 5\n3 5 10 6\n3 6 11 7\n3 8 13 9\n3 9 14 10\n3 10 15 11\nCELL_TYPES 18 \n5\n5\n5\n5\n5\n5\n5\n5\n5\n5\n5\n5\n5\n5\n5\n5\n5\n5\nPOINT_DATA 16\nSCALARS u double\nLOOKUP_TABLE default\n0\n0\n0\n0\n1\n1\n1\n1\n2\n2\n2\n2\n3\n3\n3\n3\n";
 139 %! FPL2vtkoutputdata ("__FPL2vtkoutputdata__test__.vtk", msh.p, msh.t, {"u", msh.p(1,:).'});
 140 %! fid = fopen("__FPL2vtkoutputdata__test__.vtk","r");
 141 %! s = fread(fid);
 142 %! fclose(fid);
 143 %! assert (char(s'), testfile);
 144 %! delete __FPL2vtkoutputdata__test__.vtk