Mercurial > forge

comparison extra/fpl/deprecated/fpl_dx_write_field.m @ 12671:20e8aca47b2c octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
prepare for release
author cdf
date Mon, 17 Aug 2015 10:19:39 +0000
parents extra/fpl/inst/fpl_dx_write_field.m@2748190086ad
children
comparison
equal deleted inserted replaced
12670:d68da2f2417b 12671:20e8aca47b2c
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 3 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