octave_packages/fpl-1.2.0/fpl_vtk_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_vtk_write_field (@var{basename}, @
  24 ## @var{mesh}, @var{nodedata},  @var{celldata}, @var{endfile})
  25 ##
  26 ## Output data field in serial XML-VTK UnstructuredGrid format.
  27 ##
  28 ## @var{basename} is a string containing the base-name of the (vtu) file
  29 ## where the data will be saved.
  30 ##
  31 ## @var{mesh} is a PDE-tool like mesh, like the ones generated by the
  32 ## "msh" package.
  33 ##
  34 ## @var{nodedata} and @var{celldata} are (Ndata x 2) cell arrays containing
  35 ## respectively <PointData> and <CellData> representing scalars or
  36 ## vectors:
  37 ## @itemize @minus
  38 ## @item @var{*data}@{:,1@} = variable data;
  39 ## @item @var{*data}@{:,2@} = variable names;
  40 ## @end itemize
  41 ##
  42 ## @var{endfile} should be 0 if you want to add other variables to the
  43 ## same file, 1 otherwise.
  44 ##
  45 ## Example:
  46 ## @example
  47 ## <generate msh1, node centered field nc1, cell centered field cc1>
  48 ## fpl_vtk_write_field("example",msh1,@{nc1, "temperature"@},@
  49 ## @{cc1, "density"@},0);
  50 ## <generate msh2, node centered field nc2>
  51 ## fpl_vtk_write_field("example",msh2,@{nc2, "temperature"@},@
  52 ## @{@},1);
  53 ## @end example
  54 ## will generate a valid XML-VTK UnstructuredGrid file.
  55 ##
  56 ## @seealso{fpl_dx_write_field, fpl_dx_write_series}
  57 ##
  58 ## @end deftypefn
  59
  60 function fpl_vtk_write_field(basename,mesh,nodedata,celldata,endfile)
  61
  62   ## Check input
  63   if nargin!=5
  64     error("fpl_vtk_write_field: wrong number of input");
  65   endif
  66
  67   if !ischar(basename)
  68     error("fpl_vtk_write_field: basename should be a valid string");
  69   elseif !( isstruct(mesh) )
  70     error("fpl_vtk_write_field: mesh should be a valid structure");
  71   elseif !(iscell(nodedata) && iscell(celldata))
  72     error("fpl_vtk_write_field: nodedata and celldata should be a valid cells");
  73   endif
  74
  75   filename = [basename ".vtu"];
  76
  77   if ! exist(filename,"file")
  78     fid = fopen (filename,"w");
  79     ## Header
  80     fprintf (fid,"<?xml version=""1.0""?>\n");
  81     fprintf (fid,"<VTKFile type=""UnstructuredGrid"" version=""0.1"" byte_order=""LittleEndian"">\n");
  82     fprintf (fid,"<UnstructuredGrid>\n");
  83   else
  84     ## FIXME: the following should be performed in a cleaner way! Does a
  85     ## backward fgetl function exist?
  86
  87     ## If file exist, check if it was already closed
  88     fid = fopen (filename,"r");
  89     fseek(fid,-10,SEEK_END);
  90     tst = fgetl(fid);
  91     if strcmp(tst,"</VTKFile>")
  92       error("fpl_vtk_write_field: file %s exist and was already closed",filename);
  93     endif
  94     fclose(fid);
  95     fid = fopen (filename,"a");
  96   endif
  97
  98   p   = mesh.p;
  99   dim = rows(p); # 2D or 3D
 100
 101   if dim == 2
 102     t = mesh.t(1:3,:);
 103   elseif dim == 3
 104     t = mesh.t(1:4,:);
 105   else
 106     error("fpl_vtk_write_field: neither 2D triangle nor 3D tetrahedral mesh");
 107   endif
 108
 109   t -= 1;
 110
 111   nnodes = columns(p);
 112   nelems = columns(t);
 113
 114   ## Header for <Piece>
 115   fprintf (fid,"<Piece NumberOfPoints=""%d"" NumberOfCells=""%d"">\n",nnodes,nelems);
 116   ## Print grid
 117   print_grid(fid,dim,p,nnodes,t,nelems);
 118   ## Print PointData
 119   print_data_points(fid,nodedata,nnodes)
 120   print_cell_data  (fid,celldata,nelems)
 121   ## Footer for <Piece>
 122   fprintf (fid, "</Piece>\n");
 123
 124   if (endfile)
 125     ## Footer
 126     fprintf (fid, "</UnstructuredGrid>\n");
 127     fprintf (fid, "</VTKFile>");
 128   endif
 129
 130   fclose (fid);
 131
 132 endfunction
 133
 134 ## Print Points and Cells Data
 135 function print_grid(fid,dim,p,nnodes,t,nelems)
 136
 137   if dim == 2
 138     p      = [p; zeros(1,nnodes)];
 139     eltype = 5;
 140   else
 141     eltype = 10;
 142   endif
 143
 144   ## VTK-Points (mesh nodes)
 145   fprintf (fid,"<Points>\n");
 146   fprintf (fid,"<DataArray type=""Float64"" Name=""Array"" NumberOfComponents=""3"" format=""ascii"">\n");
 147   fprintf (fid,"%g %g %g\n", p);
 148   fprintf (fid,"</DataArray>\n");
 149   fprintf (fid,"</Points>\n");
 150
 151   ## VTK-Cells (mesh elements)
 152   fprintf (fid, "<Cells>\n");
 153   fprintf (fid, "<DataArray type=""Int32"" Name=""connectivity"" format=""ascii"">\n");
 154   if dim == 2
 155     fprintf (fid, "%d %d %d \n", t);
 156   else
 157     fprintf (fid, "%d %d %d %d \n", t);
 158   endif
 159   fprintf (fid, "</DataArray>\n");
 160   fprintf (fid, "<DataArray type=""Int32"" Name=""offsets"" format=""ascii"">\n");
 161   fprintf (fid, "%d %d %d %d %d %d\n", (dim+1):(dim+1):((dim+1)*nelems));
 162   fprintf (fid, "</DataArray>\n");
 163   fprintf (fid, "<DataArray type=""Int32"" Name=""types"" format=""ascii"">\n");
 164   fprintf (fid, "%d %d %d %d %d %d\n", eltype*ones(nelems,1));
 165   fprintf (fid, "</DataArray>\n");
 166   fprintf (fid, "</Cells>\n");
 167
 168 endfunction
 169
 170 ## Print DataPoints
 171 function print_data_points(fid,nodedata,nnodes)
 172
 173   ## # of data to print in
 174   ## <PointData> field
 175   nvdata = size(nodedata,1);
 176
 177   if (nvdata)
 178     fprintf (fid, "<PointData>\n");
 179     for ii = 1:nvdata
 180       data     = nodedata{ii,1};
 181       dataname = nodedata{ii,2};
 182       nsamples = rows(data);
 183       ncomp    = columns(data);
 184       if nsamples != nnodes
 185         error("fpl_vtk_write_field: wrong number of samples in <PointData> ""%s""",dataname);
 186       endif
 187       fprintf (fid,"<DataArray type=""Float32"" Name=""%s"" ",dataname);
 188       fprintf (fid,"NumberOfComponents=""%d"" format=""ascii"">\n",ncomp);
 189       for jj = 1:nsamples
 190         fprintf (fid,"%g ",data(jj,:));
 191         fprintf (fid,"\n");
 192       endfor
 193       fprintf (fid,"</DataArray>\n");
 194     endfor
 195     fprintf (fid, "</PointData>\n");
 196   endif
 197
 198 endfunction
 199
 200 function print_cell_data(fid,celldata,nelems)
 201
 202   ## # of data to print in
 203   ## <CellData> field
 204   nvdata = size(celldata,1);
 205
 206   if (nvdata)
 207     fprintf (fid, "<CellData>\n");
 208     for ii = 1:nvdata
 209       data     = celldata{ii,1};
 210       dataname = celldata{ii,2};
 211       nsamples = rows(data);
 212       ncomp    = columns(data);
 213       if nsamples != nelems
 214         error("fpl_vtk_write_field: wrong number of samples in <CellData> ""%s""",dataname);
 215       endif
 216       fprintf (fid,"<DataArray type=""Float32"" Name=""%s"" ",dataname);
 217       fprintf (fid,"NumberOfComponents=""%d"" format=""ascii"">\n",ncomp);
 218       for jj = 1:nsamples
 219         fprintf (fid,"%g ",data(jj,:));
 220         fprintf (fid,"\n");
 221       endfor
 222       fprintf (fid,"</DataArray>\n");
 223     endfor
 224     fprintf (fid, "</CellData>\n");
 225   endif
 226
 227 endfunction