X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Ffpl-1.2.0%2Ffpl_dx_write_field.m;fp=octave_packages%2Ffpl-1.2.0%2Ffpl_dx_write_field.m;h=ca91cf405795ae6f471271f17c200606d78e715e;hp=0000000000000000000000000000000000000000;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hpb=1705066eceaaea976f010f669ce8e972f3734b05
diff --git a/octave_packages/fpl-1.2.0/fpl_dx_write_field.m b/octave_packages/fpl-1.2.0/fpl_dx_write_field.m
new file mode 100644
index 0000000..ca91cf4
--- /dev/null
+++ b/octave_packages/fpl-1.2.0/fpl_dx_write_field.m
@@ -0,0 +1,222 @@
+## 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 .
+##
+## author: Carlo de Falco
+## author: Massimiliano Culpo
+
+## -*- 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
+##
+## 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
+##
+## 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