/[escript]/trunk/finley/src/Mesh_saveVTK.c
ViewVC logotype

Contents of /trunk/finley/src/Mesh_saveVTK.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 113 - (show annotations)
Mon Feb 28 07:06:33 2005 UTC (14 years, 7 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_saveVTK.c
File MIME type: text/plain
File size: 17138 byte(s)
*** empty log message ***

1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* writes data and mesh in a vtk file */
6
7 /**************************************************************/
8
9 /* Copyrights by ACcESS, Australia 2004 */
10 /* Author: Paul Cochrane, cochrane@esscc.uq.edu.au */
11
12 /**************************************************************/
13
14 #include "Finley.h"
15 #include "Common.h"
16 #include "Mesh.h"
17 #include "escript/Data/DataC.h"
18 #include "vtkCellType.h" /* copied from vtk source directory !!! */
19
20 void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, escriptDataC* data_p) {
21 /* if there is no mesh we just return */
22 if (mesh_p==NULL) return;
23 Finley_ElementFile* elements=NULL;
24 char elemTypeStr[32];
25 int i, j, k, numVTKNodesPerElement, isCellCentered, nodetype;
26 double* values, rtmp;
27 int nDim = mesh_p->Nodes->numDim;
28
29 /* get a pointer to the relevant mesh component */
30 if (isEmpty(data_p)) {
31 elements=mesh_p->Elements;
32 } else {
33 switch(getFunctionSpaceType(data_p)) {
34 case(FINLEY_DEGREES_OF_FREEDOM):
35 nodetype = FINLEY_DEGREES_OF_FREEDOM;
36 isCellCentered = FALSE;
37 elements = mesh_p->Elements;
38 break;
39 case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
40 Finley_ErrorCode=VALUE_ERROR;
41 sprintf(Finley_ErrorMsg,
42 "Reduced degrees of freedom is not yet "
43 "implemented for saving vtk files\n");
44 return;
45 case(FINLEY_NODES):
46 nodetype=FINLEY_NODES;
47 isCellCentered=FALSE;
48 elements=mesh_p->Elements;
49 break;
50 case(FINLEY_ELEMENTS):
51 isCellCentered=TRUE;
52 elements=mesh_p->Elements;
53 break;
54 case(FINLEY_FACE_ELEMENTS):
55 isCellCentered=TRUE;
56 elements=mesh_p->FaceElements;
57 break;
58 case(FINLEY_POINTS):
59 isCellCentered=TRUE;
60 elements=mesh_p->Points;
61 break;
62 case(FINLEY_CONTACT_ELEMENTS_1):
63 case(FINLEY_CONTACT_ELEMENTS_2):
64 isCellCentered=TRUE;
65 elements=mesh_p->ContactElements;
66 break;
67 default:
68 Finley_ErrorCode=TYPE_ERROR;
69 sprintf(Finley_ErrorMsg,
70 "Finley does not know anything about function space type %d",
71 getFunctionSpaceType(data_p));
72 return;
73 }
74 }
75
76 /* the number of points */
77 int numPoints = mesh_p->Nodes->numNodes;
78
79 /* the number of cells */
80 if (elements == NULL) {
81 Finley_ErrorCode = VALUE_ERROR;
82 sprintf(Finley_ErrorMsg,
83 "elements object is NULL; cannot proceed");
84 return;
85 }
86 int numCells = elements->numElements;
87
88 /* open the file and check handle */
89 FILE * fileHandle_p = fopen(filename_p, "w");
90 if (fileHandle_p==NULL) {
91 Finley_ErrorCode=IO_ERROR;
92 sprintf(Finley_ErrorMsg,
93 "File %s could not be opened for writing.", filename_p);
94 return;
95 }
96 /* xml header */
97 fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");
98 fprintf(fileHandle_p,
99 "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");
100
101 /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */
102 fprintf(fileHandle_p, "<UnstructuredGrid>\n");
103
104 /* is there only one "piece" to the data?? */
105 fprintf(fileHandle_p, "<Piece "
106 "NumberOfPoints=\"%d\" "
107 "NumberOfCells=\"%d\">\n",
108 numPoints, numCells);
109
110 /* now for the points; equivalent to positions section in saveDX() */
111 /* "The points element explicitly defines coordinates for each point
112 * individually. It contains one DataArray element describing an array
113 * with three components per value, each specifying the coordinates of one
114 * point" - from Vtk User's Guide
115 */
116 fprintf(fileHandle_p, "<Points>\n");
117 /*
118 * the reason for this if statement is explained in the long comment below
119 */
120 if (nDim < 3) {
121 fprintf(fileHandle_p, "<DataArray "
122 "NumberOfComponents=\"3\" "
123 "type=\"Float32\" "
124 "format=\"ascii\">\n");
125 } else {
126 fprintf(fileHandle_p, "<DataArray "
127 "NumberOfComponents=\"%d\" "
128 "type=\"Float32\" "
129 "format=\"ascii\">\n",
130 nDim);
131 }
132 for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
133 fprintf(fileHandle_p,
134 "%e", mesh_p->Nodes->Coordinates[INDEX2(0, i, nDim)]);
135 for (j = 1; j < nDim; j++) {
136 fprintf(fileHandle_p,
137 " %f",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
138 /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate
139 * third dimension to handle 2D data (like a sheet of paper). So, if
140 * nDim is 2, we have to append zeros to the array to get this third
141 * dimension, and keep the visualisers happy.
142 * Indeed, if nDim is less than 3, must pad all empty dimensions, so
143 * that the total number of dims is 3.
144 */
145 if (nDim < 3) {
146 for (k=0; k<3-nDim; k++) {
147 fprintf(fileHandle_p, " 0");
148 }
149 }
150 }
151 fprintf(fileHandle_p, "\n");
152 }
153 fprintf(fileHandle_p, "</DataArray>\n");
154 fprintf(fileHandle_p, "</Points>\n");
155
156 /* connections */
157 /* now for the cells */
158 /* "The Cells element defines cells explicitly by specifying point
159 * connectivity and cell types. It contains three DataArray elements. The
160 * first array specifies the point connectivity. All cells' point lists
161 * are concatenated together. The second array specifies th eoffset into
162 * the connectivity array for the end of each cell. The third array
163 * specifies the type of each cell.
164 */
165 /* if no element table is present jump over the connection table */
166 int cellType;
167 if (elements!=NULL) {
168 fprintf(fileHandle_p, "<Cells>\n");
169 ElementTypeId TypeId = elements->ReferenceElement->Type->TypeId;
170 switch(TypeId) {
171 case Point1:
172 cellType = VTK_VERTEX;
173 break;
174 case Line2:
175 cellType = VTK_LINE;
176 break;
177 case Line3:
178 cellType = VTK_QUADRATIC_EDGE;
179 break;
180 case Tri3:
181 cellType = VTK_TRIANGLE;
182 break;
183 case Tri6:
184 cellType = VTK_QUADRATIC_TRIANGLE;
185 break;
186 case Rec4:
187 cellType = VTK_QUAD;
188 break;
189 case Rec8:
190 cellType = VTK_QUADRATIC_QUAD;
191 break;
192 case Tet4:
193 cellType = VTK_TETRA;
194 break;
195 case Tet10:
196 cellType = VTK_QUADRATIC_TETRA;
197 break;
198 case Hex8:
199 cellType = VTK_HEXAHEDRON;
200 break;
201 case Hex20:
202 cellType = VTK_QUADRATIC_HEXAHEDRON;
203 break;
204 case Line2Face:
205 cellType = VTK_VERTEX;
206 break;
207 case Line3Face:
208 cellType = VTK_VERTEX;
209 break;
210 case Tri3Face:
211 cellType = VTK_LINE;
212 break;
213 case Tri6Face:
214 cellType = VTK_QUADRATIC_EDGE;
215 break;
216 case Rec4Face:
217 cellType = VTK_LINE;
218 break;
219 case Rec8Face:
220 cellType = VTK_QUADRATIC_EDGE;
221 break;
222 case Tet4Face:
223 cellType = VTK_TRIANGLE;
224 break;
225 case Tet10Face:
226 cellType = VTK_QUADRATIC_TRIANGLE;
227 break;
228 case Hex8Face:
229 cellType = VTK_QUADRATIC_QUAD;
230 break;
231 case Hex20Face:
232 cellType = VTK_QUADRATIC_QUAD;
233 break;
234 case Point1_Contact:
235 cellType = VTK_VERTEX;
236 break;
237 case Line2_Contact:
238 cellType = VTK_LINE;
239 break;
240 case Line3_Contact:
241 cellType = VTK_QUADRATIC_EDGE;
242 break;
243 case Tri3_Contact:
244 cellType = VTK_TRIANGLE;
245 break;
246 case Tri6_Contact:
247 cellType = VTK_QUADRATIC_TRIANGLE;
248 break;
249 case Rec4_Contact:
250 cellType = VTK_QUAD;
251 break;
252 case Rec8_Contact:
253 cellType = VTK_QUADRATIC_QUAD;
254 break;
255 case Line2Face_Contact:
256 cellType = VTK_VERTEX;
257 break;
258 case Line3Face_Contact:
259 cellType = VTK_VERTEX;
260 break;
261 case Tri3Face_Contact:
262 cellType = VTK_LINE;
263 break;
264 case Tri6Face_Contact:
265 cellType = VTK_QUADRATIC_EDGE;
266 break;
267 case Rec4Face_Contact:
268 cellType = VTK_LINE;
269 break;
270 case Rec8Face_Contact:
271 cellType = VTK_QUADRATIC_EDGE;
272 break;
273 case Tet4Face_Contact:
274 cellType = VTK_TRIANGLE;
275 break;
276 case Tet10Face_Contact:
277 cellType = VTK_QUADRATIC_TRIANGLE;
278 break;
279 case Hex8Face_Contact:
280 cellType = VTK_QUAD;
281 break;
282 case Hex20Face_Contact:
283 cellType = VTK_QUADRATIC_QUAD;
284 break;
285 default:
286 Finley_ErrorCode=VALUE_ERROR;
287 sprintf(Finley_ErrorMsg,
288 "Element type %s is not supported by VTK",
289 elements->ReferenceElement->Type->Name);
290 return;
291 }
292
293 switch(cellType) {
294 case VTK_VERTEX:
295 numVTKNodesPerElement = 1;
296 strcpy(elemTypeStr, "VTK_VERTEX");
297 break;
298 case VTK_LINE:
299 numVTKNodesPerElement = 2;
300 strcpy(elemTypeStr, "VTK_LINE");
301 break;
302 case VTK_TRIANGLE:
303 numVTKNodesPerElement = 3;
304 strcpy(elemTypeStr, "VTK_TRIANGLE");
305 break;
306 case VTK_QUAD:
307 numVTKNodesPerElement = 4;
308 strcpy(elemTypeStr, "VTK_QUAD");
309 break;
310 case VTK_TETRA:
311 numVTKNodesPerElement = 4;
312 strcpy(elemTypeStr, "VTK_TETRA");
313 break;
314 case VTK_HEXAHEDRON:
315 numVTKNodesPerElement = 8;
316 strcpy(elemTypeStr, "VTK_HEXAHEDRON");
317 break;
318 case VTK_QUADRATIC_EDGE:
319 numVTKNodesPerElement = 3;
320 strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
321 break;
322 case VTK_QUADRATIC_TRIANGLE:
323 numVTKNodesPerElement = 6;
324 strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
325 break;
326 case VTK_QUADRATIC_QUAD:
327 numVTKNodesPerElement = 8;
328 strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
329 break;
330 case VTK_QUADRATIC_TETRA:
331 numVTKNodesPerElement = 10;
332 strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
333 break;
334 case VTK_QUADRATIC_HEXAHEDRON:
335 numVTKNodesPerElement = 20;
336 strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
337 break;
338 default:
339 Finley_ErrorCode = VALUE_ERROR;
340 sprintf(Finley_ErrorMsg,
341 "Cell type %d is not supported by VTK", cellType);
342 return;
343 }
344
345 /* write out the DataArray element for the connectivity */
346 fprintf(fileHandle_p, "<DataArray "
347 "Name=\"connectivity\" "
348 "type=\"Int32\" "
349 "format=\"ascii\">\n");
350 int NN = elements->ReferenceElement->Type->numNodes;
351 for (i = 0; i < numCells; i++) {
352 fprintf(fileHandle_p, "%d ",
353 mesh_p->Elements->Nodes[INDEX2(0, i, NN)]);
354 for (j = 1; j < numVTKNodesPerElement; j++) {
355 fprintf(fileHandle_p,"%d ", mesh_p->Elements->Nodes[INDEX2(j, i, NN)]);
356 }
357 }
358 fprintf(fileHandle_p, "\n");
359 fprintf(fileHandle_p, "</DataArray>\n");
360
361 /* write out the DataArray element for the offsets */
362 fprintf(fileHandle_p, "<DataArray "
363 "Name=\"offsets\" "
364 "type=\"Int32\" "
365 "format=\"ascii\">\n");
366 for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement) {
367 fprintf(fileHandle_p, "%d ", i);
368 }
369 fprintf(fileHandle_p, "\n");
370 fprintf(fileHandle_p, "</DataArray>\n");
371
372 /* write out the DataArray element for the types */
373 fprintf(fileHandle_p, "<DataArray "
374 "Name=\"types\" "
375 "type=\"UInt8\" "
376 "format=\"ascii\">\n");
377 for (i=0; i<numCells; i++) {
378 fprintf(fileHandle_p, "%d ", cellType);
379 }
380 fprintf(fileHandle_p, "\n");
381 fprintf(fileHandle_p, "</DataArray>\n");
382
383 /* finish off the <Cells> element */
384 fprintf(fileHandle_p, "</Cells>\n");
385 }
386
387 /* data */
388 if (!isEmpty(data_p)) {
389 int rank = getDataPointRank(data_p);
390 int nComp = getDataPointSize(data_p);
391 /* barf if rank is greater than two */
392 if (rank > 2) {
393 Finley_ErrorCode = VALUE_ERROR;
394 sprintf(Finley_ErrorMsg,
395 "Vtk can't handle objects with rank greater than 2\n"
396 "object rank = %d\n", rank);
397 return;
398 }
399 /* if the rank == 0: --> scalar data
400 * if the rank == 1: --> vector data
401 * if the rank == 2: --> tensor data
402 */
403 char dataNameStr[31], dataTypeStr[63];
404 int nCompReqd; /* the number of components required by vtk */
405 if (rank == 0) {
406 strcpy(dataNameStr, "escript_scalar_data");
407 sprintf(dataTypeStr, "Scalars=\"%s\"", dataNameStr);
408 nCompReqd = 1;
409 }
410 else if (rank == 1) {
411 strcpy(dataNameStr, "escript_vector_data");
412 sprintf(dataTypeStr, "Vectors=\"%s\"", dataNameStr);
413 nCompReqd = 3;
414 }
415 else if (rank == 2) {
416 strcpy(dataNameStr, "escript_tensor_data");
417 sprintf(dataTypeStr, "Tensors=\"%s\"", dataNameStr);
418 nCompReqd = 9;
419 }
420 /* if have cell centred data then use CellData element,
421 * if have node centred data, then use PointData element
422 */
423 if (isCellCentered) {
424 /* now for the cell data */
425 fprintf(fileHandle_p, "<CellData %s>\n", dataTypeStr);
426 fprintf(fileHandle_p,
427 "<DataArray "
428 "Name=\"%s\" "
429 "type=\"Float32\" "
430 "NumberOfComponents=\"%d\" "
431 "format=\"ascii\">\n",
432 dataNameStr, nCompReqd);
433 int numPointsPerSample = elements->ReferenceElement->numQuadNodes;
434 if (numPointsPerSample) {
435 for (i=0; i<numCells; i++) {
436 values = getSampleData(data_p, i);
437 double sampleAvg[nComp];
438 for (k=0; k<nComp; k++) {
439 /* averaging over the number of points in the sample */
440 rtmp = 0.;
441 for (j=0; j<numPointsPerSample; j++) {
442 rtmp += values[INDEX2(k,j,nComp)];
443 }
444 sampleAvg[k] = rtmp/numPointsPerSample;
445 }
446 /* if the number of required components is more than the number
447 * of actual components, pad with zeros
448 */
449 /* probably only need to get shape of first element */
450 int shape = getDataPointShape(data_p, 0);
451 if (shape > 3) {
452 Finley_ErrorCode = VALUE_ERROR;
453 sprintf(Finley_ErrorMsg,
454 "shape should be 1, 2, or 3; I got %d\n", shape);
455 return;
456 }
457 /* write the data different ways for scalar, vector and tensor */
458 if (nCompReqd == 1) {
459 fprintf(fileHandle_p, " %f", sampleAvg[0]);
460 }
461 else if (nCompReqd == 3) {
462 int shape = getDataPointShape(data_p, 0);
463 /* write out the data */
464 for (int m=0; m<shape; m++) {
465 fprintf(fileHandle_p, " %f", sampleAvg[m]);
466 }
467 /* pad with zeros */
468 for (int m=0; m<nCompReqd-shape; m++) {
469 fprintf(fileHandle_p, " 0");
470 }
471 }
472 else if (nCompReqd == 9) {
473 /* tensor data, so have a 3x3 matrix to output as a row
474 * of 9 data points */
475 int count = 0;
476 int elems = 0;
477 for (int m=0; m<shape; m++) {
478 for (int n=0; n<shape; n++) {
479 fprintf(fileHandle_p, " %f", sampleAvg[count]);
480 count++;
481 elems++;
482 }
483 for (int n=0; n<3-shape; n++) {
484 fprintf(fileHandle_p, " 0");
485 elems++;
486 }
487 }
488 for (int m=0; m<nCompReqd-elems; m++) {
489 fprintf(fileHandle_p, " 0");
490 }
491 }
492 fprintf(fileHandle_p, "\n");
493 }
494 }
495 fprintf(fileHandle_p, "</DataArray>\n");
496 fprintf(fileHandle_p, "</CellData>\n");
497 } else {
498 /* now for the point data */
499 fprintf(fileHandle_p, "<PointData %s>\n", dataTypeStr);
500 fprintf(fileHandle_p, "<DataArray "
501 "Name=\"%s\" "
502 "type=\"Float32\" "
503 "NumberOfComponents=\"%d\" "
504 "format=\"ascii\">\n",
505 dataNameStr, nComp);
506 for (i=0; i<mesh_p->Nodes->numNodes; i++) {
507 switch (nodetype) {
508 case(FINLEY_DEGREES_OF_FREEDOM):
509 values = getSampleData(data_p,
510 mesh_p->Nodes->degreeOfFreedom[i]);
511 break;
512 case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
513 if (mesh_p->Nodes->toReduced[i]>=0) {
514 values = getSampleData(data_p,
515 mesh_p->Nodes->reducedDegreeOfFreedom[i]);
516 }
517 break;
518 case(FINLEY_NODES):
519 values = getSampleData(data_p,i);
520 break;
521 }
522 /* write out the data */
523 /* if the number of required components is more than the number
524 * of actual components, pad with zeros
525 */
526 /* probably only need to get shape of first element */
527 int shape = getDataPointShape(data_p, 0);
528 if (shape > 3) {
529 Finley_ErrorCode = VALUE_ERROR;
530 sprintf(Finley_ErrorMsg,
531 "shape should be 1, 2, or 3; I got %d\n", shape);
532 return;
533 }
534 /* write the data different ways for scalar, vector and tensor */
535 if (nCompReqd == 1) {
536 fprintf(fileHandle_p, " %f", values[0]);
537 }
538 else if (nCompReqd == 3) {
539 int shape = getDataPointShape(data_p, 0);
540 /* write out the data */
541 for (int m=0; m<shape; m++) {
542 fprintf(fileHandle_p, " %f", values[m]);
543 }
544 /* pad with zeros */
545 for (int m=0; m<nCompReqd-shape; m++) {
546 fprintf(fileHandle_p, " 0");
547 }
548 }
549 else if (nCompReqd == 9) {
550 /* tensor data, so have a 3x3 matrix to output as a row
551 * of 9 data points */
552 int count = 0;
553 int elems = 0;
554 for (int m=0; m<shape; m++) {
555 for (int n=0; n<shape; n++) {
556 fprintf(fileHandle_p, " %f", values[count]);
557 count++;
558 elems++;
559 }
560 for (int n=0; n<3-shape; n++) {
561 fprintf(fileHandle_p, " 0");
562 elems++;
563 }
564 }
565 for (int m=0; m<nCompReqd-elems; m++) {
566 fprintf(fileHandle_p, " 0");
567 }
568 }
569 fprintf(fileHandle_p, "\n");
570 }
571 fprintf(fileHandle_p, "</DataArray>\n");
572 fprintf(fileHandle_p, "</PointData>\n");
573 }
574 }
575
576 /* finish off the piece */
577 fprintf(fileHandle_p, "</Piece>\n");
578
579 fprintf(fileHandle_p, "</UnstructuredGrid>\n");
580 /* write the xml footer */
581 fprintf(fileHandle_p, "</VTKFile>\n");
582 /* close the file */
583 fclose(fileHandle_p);
584 return;
585 }

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26