Attachment 'inp2vtu.pl'
Download 1 #!/usr/bin/perl
2 use strict;
3
4 # Copyright (C) 2008 Stefan Tibus
5 #
6 # inp2vtu.pl 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 # inp2vtu.pl 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 magpar; if not, write to the Free Software
18 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19
20 # 2007-11-20 STi: inp2vtu.pl
21 # Supports node data and writes ascii data Vtk-XML files.
22 # Nodes and cells have to be ordered and continuous.
23 #
24 # 2008-02-01 STi: inp2vtu.pl
25 # Added support for cell data.
26
27
28 # Check Number of Arguments
29 # =========================
30 if($#ARGV!=0) {
31 die("Usage: inp2vtu.pl <problem>\n");
32 }
33
34 # Get Arguments
35 # =============
36 my $problem=$ARGV[0];
37
38 # Get List of .inp Files
39 # ======================
40 my @inpfiles=`ls -1 \"$problem.\"????\".inp\"`;
41 # Strip leading/trailing spaces/tabs/newlines
42 foreach my $inpfile (@inpfiles) {
43 $inpfile=~s/^\s*(.*)\s*\n?/$1/;
44 }
45 # Skip last file if it has the number 9999
46 if(@inpfiles[-1]=~/^$problem\.9999\.inp/) {
47 pop(@inpfiles);
48 }
49
50 #DBG#foreach my $inpfile (@inpfiles) {
51 #DBG# print(stderr $inpfile."\n");
52 #DBG#}
53
54 # =================
55 # Create .vtu Files
56 # =================
57 my @vtufiles;
58
59 foreach my $inpfile (@inpfiles) {
60 my $vtufile=$inpfile;
61 # Change and store filename
62 $vtufile=~s/(.*)\.(\d\d\d\d)\.inp/$1_$2.vtu/;
63 push(@vtufiles,$vtufile);
64
65 # skip conversion if VTU file is newer!
66 if ( (stat($vtufile))[9] > (stat($inpfile))[9] ) {
67 print "Skipping $inpfile because $vtufile is newer!\n\n";
68 next;
69 }
70
71 # Read From .inp File
72 # ===================
73 open(INPFILE,"<",$inpfile) || die("Cannot open file \"$inpfile\" for read: $!\n");
74
75 # Read Header
76 # -----------
77 # nnodes, ncells, nnodedata, ncelldata, nmodeldata
78 my $line=<INPFILE>;
79 $line=~s/^\s*(.*)\s*$/$1/;
80 my($nnodes,$ncells,$nnodedata,$ncelldata,$nmodeldata)=split(/\s+/,$line,5);
81 if(($nmodeldata!=0)) {
82 die("Don't know how to handle model data in file \"$inpfile\"!\n");
83 }
84
85 # Read Nodes
86 # ----------
87 print(STDERR $nnodes," nodes\n");
88 my @nodes;
89 for(my $i=0;$i<$nnodes;$i++) {
90 $line=<INPFILE>;
91 $line=~s/^\s*(.*)\s*$/$1/;
92 my($n,$x,$y,$z)=split(/\s+/,$line,4);
93 if($n!=($i+1)) {
94 die("Don't know how to handle unordered nodes in file \"$inpfile\"!\n");
95 }
96 push(@nodes,[$x,$y,$z]);
97 }
98
99 # Read Cells
100 # ----------
101 print(STDERR $ncells," cells\n");
102 my @cells;
103 for(my $i=0;$i<$ncells;$i++) {
104 $line=<INPFILE>;
105 $line=~s/^\s*(.*)\s*$/$1/;
106 my($n,$mat,$type,$remainder)=split(/\s+/,$line,4);
107 if($type!='tet') {
108 die("Don't know how to handle cells other than tetrahedra in file \"$inpfile\"!\n");
109 }
110 if($n!=($i+1)) {
111 die("Don't know how to handle unordered cells in file \"$inpfile\"!\n");
112 }
113 my($n0,$n1,$n2,$n3)=split(/\s+/,$remainder,4);
114 push(@cells,[4,$n0,$n1,$n2,$n3,$mat]);
115 }
116
117 # Read Node Data
118 # --------------
119 my @nodelabels;
120 my @nodedimensions;
121 my @nodefirstelems;
122 my @nodelabeldatatypes;
123 my $nnodelabels;
124 my @nodevectors;
125 if($nnodedata>0) {
126 print(STDERR $nnodedata," nodedata\n");
127 $line=<INPFILE>;
128 $line=~s/^\s*(.*)\s*$/$1/;
129 my($nnodevalues,$remainder)=split(/\s+/,$line,2);
130 my(@snodevalues)=split(/\s+/,$remainder,$nnodevalues);
131 my $totalsize;
132 for(my $i=0;$i<$nnodevalues;$i++) {
133 $totalsize+=$snodevalues[$i];
134 }
135 if($totalsize!=$nnodedata) {
136 die("Node data size mismatch in file \"$inpfile\"!\n");
137 }
138 # Read data labels
139 my $label;
140 my $dimension;
141 my $firstelem=0;
142 my $labeldatatype;
143 for(my $i=0;$i<$nnodevalues;$i++) {
144 $line=<INPFILE>;
145 $line=~s/^\s*(.*)\s*$/$1/;
146 $line=~s/^(.*),.*$/$1/;
147 if($line=~/_?[XxYyZz]$/) {
148 # expect x,y,z
149 $line=~s/^(.*?)_?[XxYyZz]$/$1/;
150 if($line eq $label) {
151 # increase dimensionality of former label
152 $nodedimensions[-1]++;
153 # increase index of first element of next label
154 $firstelem++;
155 next;
156 } else {
157 $label=$line;
158 $dimension=1;
159 }
160 } else {
161 $label=$line;
162 $dimension=1;
163 }
164 push(@nodelabels,$label);
165 push(@nodedimensions,$dimension);
166 push(@nodefirstelems,$firstelem);
167 if($label eq "id") {
168 $labeldatatype="UInt32";
169 } elsif($label eq "proc") {
170 $labeldatatype="UInt32";
171 } elsif($label eq "vert_pid") {
172 $labeldatatype="UInt32";
173 } elsif($label eq "prop") {
174 $labeldatatype="UInt32";
175 } else {
176 $labeldatatype="Float32";
177 }
178 push(@nodelabeldatatypes,$labeldatatype);
179 # increase index of first element of next label
180 $firstelem++;
181 }
182 $nnodelabels=scalar(@nodelabels);
183 # Read data values
184 for(my $i=0;$i<$nnodes;$i++) {
185 $line=<INPFILE>;
186 $line=~s/^\s*(.*)\s*$/$1/;
187 my($n,@vector)=split(/\s+/,$line,($nnodevalues+1));
188 if($n!=($i+1)) {
189 die("Don't know how to handle unordered node data in file \"$inpfile\"!\n");
190 }
191 push(@nodevectors,[@vector]);
192 }
193 }
194
195 # Read Cell Data
196 # --------------
197 my @celllabels;
198 my @celldimensions;
199 my @cellfirstelems;
200 my @celllabeldatatypes;
201 my $ncelllabels;
202 my @cellvectors;
203 if($ncelldata>0) {
204 print(STDERR $ncelldata," celldata\n");
205 $line=<INPFILE>;
206 $line=~s/^\s*(.*)\s*$/$1/;
207 my($ncellvalues,$remainder)=split(/\s+/,$line,2);
208 my(@scellvalues)=split(/\s+/,$remainder,$ncellvalues);
209 my $totalsize;
210 for(my $i=0;$i<$ncellvalues;$i++) {
211 $totalsize+=$scellvalues[$i];
212 }
213 if($totalsize!=$ncelldata) {
214 die("Cell data size mismatch in file \"$inpfile\"!\n");
215 }
216 # Read data labels
217 my $label;
218 my $dimension;
219 my $firstelem=0;
220 my $labeldatatype;
221 for(my $i=0;$i<$ncellvalues;$i++) {
222 $line=<INPFILE>;
223 $line=~s/^\s*(.*)\s*$/$1/;
224 $line=~s/^(.*),.*$/$1/;
225 if($line=~/_?[XxYyZz]$/) {
226 # expect x,y,z
227 $line=~s/^(.*?)_?[XxYyZz]$/$1/;
228 if($line eq $label) {
229 # increase dimensionality of former label
230 $celldimensions[-1]++;
231 # increase index of first element of next label
232 $firstelem++;
233 next;
234 } else {
235 $label=$line;
236 $dimension=1;
237 }
238 } else {
239 $label=$line;
240 $dimension=1;
241 }
242 push(@celllabels,$label);
243 push(@celldimensions,$dimension);
244 push(@cellfirstelems,$firstelem);
245 if($label eq "id") {
246 $labeldatatype="UInt32";
247 } elsif($label eq "proc") {
248 $labeldatatype="UInt32";
249 } elsif($label eq "vert_pid") {
250 $labeldatatype="UInt32";
251 } elsif($label eq "prop") {
252 $labeldatatype="UInt32";
253 } else {
254 $labeldatatype="Float32";
255 }
256 push(@celllabeldatatypes,$labeldatatype);
257 # increase index of first element of next label
258 $firstelem++;
259 }
260 $ncelllabels=scalar(@celllabels);
261 # Read data values
262 for(my $i=0;$i<$ncells;$i++) {
263 $line=<INPFILE>;
264 $line=~s/^\s*(.*)\s*$/$1/;
265 my($n,@vector)=split(/\s+/,$line,($ncellvalues+1));
266 if($n!=($i+1)) {
267 die("Don't know how to handle unordered cell data in file \"$inpfile\"!\n");
268 }
269 push(@cellvectors,[@vector]);
270 }
271 }
272
273 # Close .inp File
274 # ---------------
275 close(INPFILE);
276
277 # Write to .vtu File
278 # ==================
279 open(VTUFILE,">",$vtufile) || die("Cannot open or create file \"$vtufile\" for write: $!\n");
280 binmode(VTUFILE,':raw');
281 # Write Header
282 # ------------
283 print(VTUFILE "<?xml version=\"1.0\"?>\n");
284 print(VTUFILE "<VTKFile type=\"UnstructuredGrid\" byte_order=\"LittleEndian\">\n");
285 # Note: Intel 80x86 is little endian, Sun SPARC and IBM PowerPC are big endian platforms.
286 # See: http://www.netrino.com/Publications/Glossary/Endianness.php
287 # See: http://www.intel.com/design/intarch/papers/endian.pdf
288 print(VTUFILE " <UnstructuredGrid>\n");
289 print(VTUFILE " <Piece NumberOfPoints=\"$nnodes\" NumberOfCells=\"$ncells\">\n");
290
291 # Write Point (Node) Data
292 # -----------------------
293 # User data
294 # Set $nnodelabels=1 if you want first data set only (data for mumag / M for magpar)
295 #$nnodelabels=1;
296 if($nnodelabels>0) {
297 if($nodedimensions[0]>1) {
298 print(VTUFILE " <PointData Vectors=\"$nodelabels[0]\">\n");
299 } else {
300 print(VTUFILE " <PointData>\n");
301 }
302 for(my $i=0;$i<$nnodelabels;$i++) {
303 my $label=$nodelabels[$i];
304 my $dimension=$nodedimensions[$i];
305 my $firstelem=$nodefirstelems[$i];
306 my $labeldatatype=$nodelabeldatatypes[$i];
307 if($dimension==1) {
308 print(VTUFILE " <DataArray Name=\"$label\" NumberOfComponents=\"$dimension\" type=\"$labeldatatype\" format=\"ascii\">\n");
309 for(my $j=0;$j<$nnodes;$j++) {
310 printf(VTUFILE "%.4f\n",$nodevectors[$j][$firstelem]);
311 }
312 } elsif($dimension==2) {
313 print(VTUFILE " <DataArray Name=\"$label\" NumberOfComponents=\"$dimension\" type=\"$labeldatatype\" format=\"ascii\">\n");
314 for(my $j=0;$j<$nnodes;$j++) {
315 printf(VTUFILE "%.4f %.4f\n",$nodevectors[$j][$firstelem],$nodevectors[$j][$firstelem+1]);
316 }
317 } elsif($dimension==3) {
318 print(VTUFILE " <DataArray Name=\"$label\" NumberOfComponents=\"$dimension\" type=\"$labeldatatype\" format=\"ascii\">\n");
319 for(my $j=0;$j<$nnodes;$j++) {
320 printf(VTUFILE "%.4f %.4f %.4f\n",$nodevectors[$j][$firstelem],$nodevectors[$j][$firstelem+1],$nodevectors[$j][$firstelem+2]);
321 }
322 } else {
323 die("Don't know how to handle vectors of dimension > 3 in file \"$inpfile\"!\n");
324 }
325 print(VTUFILE " </DataArray>\n");
326 }
327 print(VTUFILE " </PointData>\n");
328 }
329
330 # Write Cell Data
331 # ---------------
332 # Material
333 print(VTUFILE " <CellData Scalars=\"Material\">\n");
334 print(VTUFILE " <DataArray Name=\"Material\" NumberOfComponents=\"1\" type=\"UInt32\" format=\"ascii\">\n");
335 for (my $i=0;$i<$ncells;$i++) {
336 printf(VTUFILE " %d",$cells[$i][5]);
337 }
338 print(VTUFILE "\n");
339 print(VTUFILE " </DataArray>\n");
340 # User data
341 # Set $ncelllabels=0 if you want to suppress cell data
342 #$ncelllabels=0;
343 if($ncelllabels>0) {
344 for(my $i=0;$i<$ncelllabels;$i++) {
345 my $label=$celllabels[$i];
346 my $dimension=$celldimensions[$i];
347 my $firstelem=$cellfirstelems[$i];
348 my $labeldatatype=$celllabeldatatypes[$i];
349 if($dimension==1) {
350 print(VTUFILE " <DataArray Name=\"$label\" NumberOfComponents=\"$dimension\" type=\"$labeldatatype\" format=\"ascii\">\n");
351 for(my $j=0;$j<$ncells;$j++) {
352 printf(VTUFILE "%.4f\n",$cellvectors[$j][$firstelem]);
353 }
354 } elsif($dimension==2) {
355 print(VTUFILE " <DataArray Name=\"$label\" NumberOfComponents=\"$dimension\" type=\"$labeldatatype\" format=\"ascii\">\n");
356 for(my $j=0;$j<$ncells;$j++) {
357 printf(VTUFILE "%.4f %.4f\n",$cellvectors[$j][$firstelem],$cellvectors[$j][$firstelem+1]);
358 }
359 } elsif($dimension==3) {
360 print(VTUFILE " <DataArray Name=\"$label\" NumberOfComponents=\"$dimension\" type=\"$labeldatatype\" format=\"ascii\">\n");
361 for(my $j=0;$j<$ncells;$j++) {
362 printf(VTUFILE "%.4f %.4f %.4f\n",$cellvectors[$j][$firstelem],$cellvectors[$j][$firstelem+1],$cellvectors[$j][$firstelem+2]);
363 }
364 } else {
365 die("Don't know how to handle vectors of dimension > 3 in file \"$inpfile\"!\n");
366 }
367 print(VTUFILE " </DataArray>\n");
368 }
369 }
370 print(VTUFILE " </CellData>\n");
371
372 # Write Nodes
373 # -----------
374 # Point coordinates
375 print(VTUFILE " <Points>\n");
376 print(VTUFILE " <DataArray NumberOfComponents=\"3\" type=\"Float32\" format=\"ascii\">\n");
377 for(my $i=0;$i<$nnodes;$i++) {
378 printf(VTUFILE "%.4f %.4f %.4f\n",$nodes[$i][0],$nodes[$i][1],$nodes[$i][2]);
379 }
380 print(VTUFILE " </DataArray>\n");
381 print(VTUFILE " </Points>\n");
382
383 # Write Cells
384 # -----------
385 # Cell connectivity (tetrahedrons only)
386 print(VTUFILE " <Cells>\n");
387 print(VTUFILE " <DataArray Name=\"connectivity\" type=\"UInt32\" format=\"ascii\">\n");
388 for(my $i=0;$i<$ncells;$i++) {
389 printf(VTUFILE "%d %d %d %d\n",$cells[$i][1]-1,$cells[$i][2]-1,$cells[$i][3]-1,$cells[$i][4]-1);
390 }
391 print(VTUFILE " </DataArray>\n");
392 # Cell connectivity offsets (for tetrahedrons)
393 print(VTUFILE " <DataArray Name=\"offsets\" type=\"UInt32\" format=\"ascii\">\n");
394 for(my $i=0;$i<$ncells;$i++) {
395 printf(VTUFILE " %d",($i+1)*4);
396 }
397 print(VTUFILE "\n");
398 print(VTUFILE " </DataArray>\n");
399 # Cell types (tetrahedrons only)
400 # Tetrahedron = 10
401 print(VTUFILE " <DataArray Name=\"types\" type=\"UInt8\" format=\"ascii\">\n");
402 for(my $i=0;$i<$ncells;$i++) {
403 printf(VTUFILE " %d",10);
404 }
405 print(VTUFILE "\n");
406 print(VTUFILE " </DataArray>\n");
407 print(VTUFILE " </Cells>\n");
408
409 # Close .vtu File
410 # ---------------
411 print(VTUFILE " </Piece>\n");
412 print(VTUFILE " </UnstructuredGrid>\n");
413 print(VTUFILE "</VTKFile>\n");
414 close(VTUFILE);
415 print($vtufile."\n");
416 }
417
418 # ================
419 # Create .pvd File
420 # ================
421
422 if(scalar(@vtufiles)>1) {
423 my $pvdfile=$problem.".pvd";
424
425 open(PVDFILE,">",$pvdfile) || die("Cannot open or create file \"$pvdfile\" for write: $!\n");
426
427 # Write Header
428 # ============
429 print(PVDFILE "<?xml version=\"1.0\"?>\n");
430 print(PVDFILE "<VTKFile type=\"Collection\">\n");
431 print(PVDFILE " <Collection>\n");
432
433 # Get List of .vtu Files
434 # ======================
435 my @vtufiles=`ls -1 \"${problem}_\"????\".vtu\"`;
436 # Strip leading/trailing spaces/tabs/newlines
437 foreach my $vtufile (@vtufiles) {
438 $vtufile=~s/^\s*(.*)\s*\n?/$1/;
439 }
440
441 # Write List of .vtu Files
442 # ========================
443 my $i=0;
444 foreach my $vtufile (@vtufiles) {
445 printf(PVDFILE " <DataSet timestep=\"%04d\" file=\"%s\"/>\n",$i,$vtufile);
446 $i++;
447 }
448
449 # Close .pvd File
450 # ===============
451 print(PVDFILE " </Collection>\n");
452 print(PVDFILE "</VTKFile>\n");
453 close(PVDFILE);
454 print($pvdfile."\n");
455 }
456
457 # vim:set ft=perl sts=4 sw=4:
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.