--- /dev/null
+## 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
\ No newline at end of file
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff