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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (13 years, 10 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_saveVTK.c
File MIME type: text/plain
File size: 23534 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

1 jgs 150 /*
2     ******************************************************************************
3     * *
4     * COPYRIGHT ACcESS 2003,2004,2005 - All Rights Reserved *
5     * *
6     * This software is the property of ACcESS. No part of this code *
7     * may be copied in any form or by any means without the expressed written *
8     * consent of ACcESS. Copying, use or modification of this software *
9     * by any unauthorised person is illegal unless that person has a software *
10     * license agreement with ACcESS. *
11     * *
12     ******************************************************************************
13     */
14 jgs 110
15 jgs 150
16 jgs 110 /**************************************************************/
17    
18     /* writes data and mesh in a vtk file */
19    
20     /**************************************************************/
21    
22     /* Author: Paul Cochrane, cochrane@esscc.uq.edu.au */
23 jgs 150 /* Version: $Id$ */
24 jgs 110
25     /**************************************************************/
26    
27     #include "Mesh.h"
28 jgs 113 #include "vtkCellType.h" /* copied from vtk source directory !!! */
29 jgs 110
30 jgs 150 /**************************************************************/
31    
32 jgs 110 void Finley_Mesh_saveVTK(const char * filename_p, Finley_Mesh *mesh_p, escriptDataC* data_p) {
33 jgs 150 char error_msg[LenErrorMsg_MAX];
34 jgs 110 /* if there is no mesh we just return */
35     if (mesh_p==NULL) return;
36     Finley_ElementFile* elements=NULL;
37     char elemTypeStr[32];
38 jgs 123 int i, j, k, numVTKNodesPerElement, isCellCentered=FALSE, nodetype=FINLEY_DEGREES_OF_FREEDOM;
39 jgs 150 index_t j2;
40 jgs 113 double* values, rtmp;
41 jgs 110 int nDim = mesh_p->Nodes->numDim;
42 jgs 150 int numPoints=0;
43    
44     if (nDim==1) {
45     Finley_setError(VALUE_ERROR,"saveVTK: 1-dimensional domains are not supported.");
46     return;
47     }
48 jgs 113 /* get a pointer to the relevant mesh component */
49     if (isEmpty(data_p)) {
50 jgs 150 numPoints = mesh_p->Nodes->numNodes;
51 jgs 113 elements=mesh_p->Elements;
52     } else {
53     switch(getFunctionSpaceType(data_p)) {
54     case(FINLEY_DEGREES_OF_FREEDOM):
55 jgs 150 numPoints = mesh_p->Nodes->numNodes;
56 jgs 113 nodetype = FINLEY_DEGREES_OF_FREEDOM;
57     isCellCentered = FALSE;
58     elements = mesh_p->Elements;
59     break;
60     case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
61 jgs 150 numPoints = mesh_p->Nodes->reducedNumNodes;
62     nodetype =FINLEY_REDUCED_DEGREES_OF_FREEDOM;
63     isCellCentered = FALSE;
64     elements = mesh_p->Elements;
65     break;
66 jgs 113 case(FINLEY_NODES):
67 jgs 150 numPoints = mesh_p->Nodes->numNodes;
68 jgs 113 nodetype=FINLEY_NODES;
69     isCellCentered=FALSE;
70     elements=mesh_p->Elements;
71     break;
72     case(FINLEY_ELEMENTS):
73 jgs 150 numPoints = mesh_p->Nodes->numNodes;
74     nodetype=FINLEY_NODES;
75 jgs 113 isCellCentered=TRUE;
76     elements=mesh_p->Elements;
77     break;
78     case(FINLEY_FACE_ELEMENTS):
79 jgs 150 numPoints = mesh_p->Nodes->numNodes;
80     nodetype=FINLEY_NODES;
81 jgs 113 isCellCentered=TRUE;
82     elements=mesh_p->FaceElements;
83     break;
84     case(FINLEY_POINTS):
85 jgs 150 numPoints = mesh_p->Nodes->numNodes;
86     nodetype=FINLEY_NODES;
87 jgs 113 isCellCentered=TRUE;
88     elements=mesh_p->Points;
89     break;
90     case(FINLEY_CONTACT_ELEMENTS_1):
91     case(FINLEY_CONTACT_ELEMENTS_2):
92 jgs 150 numPoints = mesh_p->Nodes->numNodes;
93     nodetype=FINLEY_NODES;
94 jgs 113 isCellCentered=TRUE;
95     elements=mesh_p->ContactElements;
96     break;
97     default:
98 jgs 150 sprintf(error_msg,"saveVTK: Finley does not know anything about function space type %d",getFunctionSpaceType(data_p));
99     Finley_setError(TYPE_ERROR,error_msg);
100 jgs 113 return;
101     }
102     }
103    
104     /* the number of points */
105    
106     /* the number of cells */
107     if (elements == NULL) {
108 jgs 150 Finley_setError(VALUE_ERROR,"saveVTK: elements object is NULL; cannot proceed");
109 jgs 113 return;
110     }
111     int numCells = elements->numElements;
112     /* open the file and check handle */
113 jgs 110 FILE * fileHandle_p = fopen(filename_p, "w");
114     if (fileHandle_p==NULL) {
115 jgs 150 sprintf(error_msg, "saveVTK: File %s could not be opened for writing.", filename_p);
116     Finley_setError(IO_ERROR,error_msg);
117 jgs 113 return;
118 jgs 110 }
119     /* xml header */
120     fprintf(fileHandle_p, "<?xml version=\"1.0\"?>\n");
121 jgs 150 fprintf(fileHandle_p, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\">\n");
122 jgs 110
123     /* finley uses an unstructured mesh, so UnstructuredGrid *should* work */
124     fprintf(fileHandle_p, "<UnstructuredGrid>\n");
125    
126     /* is there only one "piece" to the data?? */
127 jgs 113 fprintf(fileHandle_p, "<Piece "
128 jgs 150 "NumberOfPoints=\"%d\" "
129     "NumberOfCells=\"%d\">\n",
130     numPoints, numCells);
131 jgs 113 /* now for the points; equivalent to positions section in saveDX() */
132     /* "The points element explicitly defines coordinates for each point
133     * individually. It contains one DataArray element describing an array
134     * with three components per value, each specifying the coordinates of one
135     * point" - from Vtk User's Guide
136     */
137 jgs 110 fprintf(fileHandle_p, "<Points>\n");
138 jgs 113 /*
139     * the reason for this if statement is explained in the long comment below
140     */
141     if (nDim < 3) {
142     fprintf(fileHandle_p, "<DataArray "
143     "NumberOfComponents=\"3\" "
144     "type=\"Float32\" "
145     "format=\"ascii\">\n");
146     } else {
147     fprintf(fileHandle_p, "<DataArray "
148     "NumberOfComponents=\"%d\" "
149     "type=\"Float32\" "
150     "format=\"ascii\">\n",
151     nDim);
152     }
153 jgs 150 /* vtk/mayavi doesn't like 2D data, it likes 3D data with a degenerate
154     * third dimension to handle 2D data (like a sheet of paper). So, if
155     * nDim is 2, we have to append zeros to the array to get this third
156     * dimension, and keep the visualisers happy.
157     * Indeed, if nDim is less than 3, must pad all empty dimensions, so
158     * that the total number of dims is 3.
159     */
160     if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
161     for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
162     if (mesh_p->Nodes->toReduced[i]>=0) {
163     for (j = 0; j < nDim; j++)
164     fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
165     for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
166     fprintf(fileHandle_p, "\n");
167     }
168     }
169     } else {
170     for (i = 0; i < mesh_p->Nodes->numNodes; i++) {
171     for (j = 0; j < nDim; j++)
172     fprintf(fileHandle_p, " %e",mesh_p->Nodes->Coordinates[INDEX2(j, i, nDim)]);
173     for (k=0; k<3-nDim; k++) fprintf(fileHandle_p, " %e",0.);
174     fprintf(fileHandle_p, "\n");
175     }
176 jgs 113 }
177 jgs 110 fprintf(fileHandle_p, "</DataArray>\n");
178     fprintf(fileHandle_p, "</Points>\n");
179    
180 jgs 113 /* connections */
181 jgs 110 /* now for the cells */
182 jgs 113 /* "The Cells element defines cells explicitly by specifying point
183     * connectivity and cell types. It contains three DataArray elements. The
184     * first array specifies the point connectivity. All cells' point lists
185     * are concatenated together. The second array specifies th eoffset into
186     * the connectivity array for the end of each cell. The third array
187     * specifies the type of each cell.
188     */
189     /* if no element table is present jump over the connection table */
190     if (elements!=NULL) {
191 jgs 150 int cellType;
192     ElementTypeId TypeId;
193    
194     if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
195     TypeId = elements->LinearReferenceElement->Type->TypeId;
196     } else {
197     TypeId = elements->ReferenceElement->Type->TypeId;
198     }
199    
200 jgs 113 switch(TypeId) {
201     case Point1:
202     cellType = VTK_VERTEX;
203     break;
204     case Line2:
205     cellType = VTK_LINE;
206     break;
207     case Line3:
208     cellType = VTK_QUADRATIC_EDGE;
209     break;
210     case Tri3:
211     cellType = VTK_TRIANGLE;
212     break;
213     case Tri6:
214     cellType = VTK_QUADRATIC_TRIANGLE;
215     break;
216     case Rec4:
217     cellType = VTK_QUAD;
218     break;
219     case Rec8:
220     cellType = VTK_QUADRATIC_QUAD;
221     break;
222     case Tet4:
223     cellType = VTK_TETRA;
224     break;
225     case Tet10:
226     cellType = VTK_QUADRATIC_TETRA;
227     break;
228     case Hex8:
229     cellType = VTK_HEXAHEDRON;
230     break;
231     case Hex20:
232     cellType = VTK_QUADRATIC_HEXAHEDRON;
233     break;
234     case Line2Face:
235     cellType = VTK_VERTEX;
236     break;
237     case Line3Face:
238     cellType = VTK_VERTEX;
239     break;
240     case Tri3Face:
241     cellType = VTK_LINE;
242     break;
243     case Tri6Face:
244     cellType = VTK_QUADRATIC_EDGE;
245     break;
246     case Rec4Face:
247     cellType = VTK_LINE;
248     break;
249     case Rec8Face:
250     cellType = VTK_QUADRATIC_EDGE;
251     break;
252     case Tet4Face:
253     cellType = VTK_TRIANGLE;
254     break;
255     case Tet10Face:
256     cellType = VTK_QUADRATIC_TRIANGLE;
257     break;
258     case Hex8Face:
259 jgs 150 cellType = VTK_QUAD;
260 jgs 113 break;
261     case Hex20Face:
262     cellType = VTK_QUADRATIC_QUAD;
263     break;
264     case Point1_Contact:
265     cellType = VTK_VERTEX;
266     break;
267     case Line2_Contact:
268     cellType = VTK_LINE;
269     break;
270     case Line3_Contact:
271     cellType = VTK_QUADRATIC_EDGE;
272     break;
273     case Tri3_Contact:
274 jgs 150 cellType = VTK_LINE;
275 jgs 113 break;
276     case Tri6_Contact:
277     cellType = VTK_QUADRATIC_TRIANGLE;
278     break;
279     case Rec4_Contact:
280     cellType = VTK_QUAD;
281     break;
282     case Rec8_Contact:
283     cellType = VTK_QUADRATIC_QUAD;
284     break;
285     case Line2Face_Contact:
286     cellType = VTK_VERTEX;
287     break;
288     case Line3Face_Contact:
289     cellType = VTK_VERTEX;
290     break;
291     case Tri3Face_Contact:
292     cellType = VTK_LINE;
293     break;
294     case Tri6Face_Contact:
295     cellType = VTK_QUADRATIC_EDGE;
296     break;
297     case Rec4Face_Contact:
298     cellType = VTK_LINE;
299     break;
300     case Rec8Face_Contact:
301     cellType = VTK_QUADRATIC_EDGE;
302     break;
303     case Tet4Face_Contact:
304     cellType = VTK_TRIANGLE;
305     break;
306     case Tet10Face_Contact:
307     cellType = VTK_QUADRATIC_TRIANGLE;
308     break;
309     case Hex8Face_Contact:
310     cellType = VTK_QUAD;
311     break;
312     case Hex20Face_Contact:
313     cellType = VTK_QUADRATIC_QUAD;
314     break;
315     default:
316 jgs 150 sprintf(error_msg, "saveVTK: Element type %s is not supported by VTK",elements->ReferenceElement->Type->Name);
317     Finley_setError(VALUE_ERROR,error_msg);
318 jgs 113 return;
319     }
320 jgs 110
321 jgs 113 switch(cellType) {
322     case VTK_VERTEX:
323     numVTKNodesPerElement = 1;
324     strcpy(elemTypeStr, "VTK_VERTEX");
325     break;
326     case VTK_LINE:
327     numVTKNodesPerElement = 2;
328     strcpy(elemTypeStr, "VTK_LINE");
329     break;
330     case VTK_TRIANGLE:
331     numVTKNodesPerElement = 3;
332     strcpy(elemTypeStr, "VTK_TRIANGLE");
333     break;
334     case VTK_QUAD:
335     numVTKNodesPerElement = 4;
336     strcpy(elemTypeStr, "VTK_QUAD");
337     break;
338     case VTK_TETRA:
339     numVTKNodesPerElement = 4;
340     strcpy(elemTypeStr, "VTK_TETRA");
341     break;
342     case VTK_HEXAHEDRON:
343     numVTKNodesPerElement = 8;
344     strcpy(elemTypeStr, "VTK_HEXAHEDRON");
345     break;
346     case VTK_QUADRATIC_EDGE:
347     numVTKNodesPerElement = 3;
348     strcpy(elemTypeStr, "VTK_QUADRATIC_EDGE");
349     break;
350     case VTK_QUADRATIC_TRIANGLE:
351     numVTKNodesPerElement = 6;
352     strcpy(elemTypeStr, "VTK_QUADRATIC_TRIANGLE");
353     break;
354     case VTK_QUADRATIC_QUAD:
355     numVTKNodesPerElement = 8;
356     strcpy(elemTypeStr, "VTK_QUADRATIC_QUAD");
357     break;
358     case VTK_QUADRATIC_TETRA:
359     numVTKNodesPerElement = 10;
360     strcpy(elemTypeStr, "VTK_QUADRATIC_TETRA");
361     break;
362     case VTK_QUADRATIC_HEXAHEDRON:
363     numVTKNodesPerElement = 20;
364     strcpy(elemTypeStr, "VTK_QUADRATIC_HEXAHEDRON");
365     break;
366     default:
367 jgs 150 sprintf(error_msg,"saveVTK: Cell type %d is not supported by VTK", cellType);
368     Finley_setError(VALUE_ERROR,error_msg);
369 jgs 113 return;
370     }
371 jgs 110
372 jgs 113 /* write out the DataArray element for the connectivity */
373 jgs 150 int NN = elements->ReferenceElement->Type->numNodes;
374     fprintf(fileHandle_p, "<Cells>\n");
375 jgs 113 fprintf(fileHandle_p, "<DataArray "
376     "Name=\"connectivity\" "
377     "type=\"Int32\" "
378     "format=\"ascii\">\n");
379 jgs 150 if (nodetype==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
380 jgs 147 for (i = 0; i < numCells; i++) {
381 jgs 150 for (j = 0; j < numVTKNodesPerElement; j++) {
382     j2=elements->ReferenceElement->Type->linearNodes[j];
383     fprintf(fileHandle_p,"%d ",mesh_p->Nodes->toReduced[elements->Nodes[INDEX2(j2, i, NN)]]);
384     }
385     fprintf(fileHandle_p, "\n");
386 jgs 147 }
387 jgs 150 } else if (VTK_QUADRATIC_HEXAHEDRON==cellType) {
388     for (i = 0; i < numCells; i++) {
389     fprintf(fileHandle_p,"%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
390     elements->Nodes[INDEX2(0, i, NN)],
391     elements->Nodes[INDEX2(1, i, NN)],
392     elements->Nodes[INDEX2(2, i, NN)],
393     elements->Nodes[INDEX2(3, i, NN)],
394     elements->Nodes[INDEX2(4, i, NN)],
395     elements->Nodes[INDEX2(5, i, NN)],
396     elements->Nodes[INDEX2(6, i, NN)],
397     elements->Nodes[INDEX2(7, i, NN)],
398     elements->Nodes[INDEX2(8, i, NN)],
399     elements->Nodes[INDEX2(9, i, NN)],
400     elements->Nodes[INDEX2(10, i, NN)],
401     elements->Nodes[INDEX2(11, i, NN)],
402     elements->Nodes[INDEX2(16, i, NN)],
403     elements->Nodes[INDEX2(17, i, NN)],
404     elements->Nodes[INDEX2(18, i, NN)],
405     elements->Nodes[INDEX2(19, i, NN)],
406     elements->Nodes[INDEX2(12, i, NN)],
407     elements->Nodes[INDEX2(13, i, NN)],
408     elements->Nodes[INDEX2(14, i, NN)],
409     elements->Nodes[INDEX2(15, i, NN)]);
410     }
411     } else if (numVTKNodesPerElement!=NN) {
412     for (i = 0; i < numCells; i++) {
413     for (j = 0; j < numVTKNodesPerElement; j++) {
414     j2=elements->ReferenceElement->Type->geoNodes[j];
415     fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j2, i, NN)]);
416     }
417     fprintf(fileHandle_p, "\n");
418     }
419     } else {
420     for (i = 0; i < numCells; i++) {
421     for (j = 0; j < numVTKNodesPerElement; j++) fprintf(fileHandle_p,"%d ", elements->Nodes[INDEX2(j, i, NN)]);
422     fprintf(fileHandle_p, "\n");
423     }
424     }
425 jgs 113 fprintf(fileHandle_p, "</DataArray>\n");
426    
427     /* write out the DataArray element for the offsets */
428     fprintf(fileHandle_p, "<DataArray "
429     "Name=\"offsets\" "
430     "type=\"Int32\" "
431     "format=\"ascii\">\n");
432 jgs 150 for (i=numVTKNodesPerElement; i<=numCells*numVTKNodesPerElement; i+=numVTKNodesPerElement)
433     fprintf(fileHandle_p, "%d\n", i);
434 jgs 113 fprintf(fileHandle_p, "</DataArray>\n");
435    
436     /* write out the DataArray element for the types */
437     fprintf(fileHandle_p, "<DataArray "
438     "Name=\"types\" "
439     "type=\"UInt8\" "
440     "format=\"ascii\">\n");
441 jgs 150 for (i=0; i<numCells; i++)
442     fprintf(fileHandle_p, "%d\n", cellType);
443 jgs 113 fprintf(fileHandle_p, "</DataArray>\n");
444    
445     /* finish off the <Cells> element */
446     fprintf(fileHandle_p, "</Cells>\n");
447 jgs 110 }
448    
449     /* data */
450     if (!isEmpty(data_p)) {
451 jgs 113 int rank = getDataPointRank(data_p);
452     int nComp = getDataPointSize(data_p);
453     /* barf if rank is greater than two */
454     if (rank > 2) {
455 jgs 150 sprintf(error_msg, "saveVTK: Vtk can't handle objects with rank greater than 2. object rank = %d\n", rank);
456     Finley_setError(VALUE_ERROR,error_msg);
457 jgs 113 return;
458     }
459     /* if the rank == 0: --> scalar data
460     * if the rank == 1: --> vector data
461     * if the rank == 2: --> tensor data
462     */
463     char dataNameStr[31], dataTypeStr[63];
464 jgs 123 int nCompReqd=1; /* the number of components required by vtk */
465 jgs 113 if (rank == 0) {
466 jgs 150 strcpy(dataNameStr, "scalar");
467 jgs 113 sprintf(dataTypeStr, "Scalars=\"%s\"", dataNameStr);
468     nCompReqd = 1;
469     }
470     else if (rank == 1) {
471 jgs 150 strcpy(dataNameStr, "vector");
472 jgs 113 sprintf(dataTypeStr, "Vectors=\"%s\"", dataNameStr);
473     nCompReqd = 3;
474     }
475     else if (rank == 2) {
476 jgs 150 strcpy(dataNameStr, "tensor");
477 jgs 113 sprintf(dataTypeStr, "Tensors=\"%s\"", dataNameStr);
478     nCompReqd = 9;
479     }
480     /* if have cell centred data then use CellData element,
481     * if have node centred data, then use PointData element
482     */
483     if (isCellCentered) {
484     /* now for the cell data */
485     fprintf(fileHandle_p, "<CellData %s>\n", dataTypeStr);
486     fprintf(fileHandle_p,
487     "<DataArray "
488     "Name=\"%s\" "
489     "type=\"Float32\" "
490     "NumberOfComponents=\"%d\" "
491     "format=\"ascii\">\n",
492     dataNameStr, nCompReqd);
493     int numPointsPerSample = elements->ReferenceElement->numQuadNodes;
494     if (numPointsPerSample) {
495 jgs 150 int shape = getDataPointShape(data_p, 0);
496     if (shape > 3) {
497     sprintf(error_msg,"saveVTK: shape should be 1, 2, or 3; I got %d\n", shape);
498     Finley_setError(VALUE_ERROR,error_msg);
499     return;
500     }
501 jgs 113 for (i=0; i<numCells; i++) {
502     values = getSampleData(data_p, i);
503     double sampleAvg[nComp];
504     for (k=0; k<nComp; k++) {
505     /* averaging over the number of points in the sample */
506     rtmp = 0.;
507 jgs 150 for (j=0; j<numPointsPerSample; j++)
508 jgs 113 rtmp += values[INDEX2(k,j,nComp)];
509     sampleAvg[k] = rtmp/numPointsPerSample;
510     }
511     /* if the number of required components is more than the number
512     * of actual components, pad with zeros
513     */
514     /* probably only need to get shape of first element */
515     /* write the data different ways for scalar, vector and tensor */
516     if (nCompReqd == 1) {
517 jgs 150 fprintf(fileHandle_p, " %e", sampleAvg[0]);
518     } else if (nCompReqd == 3) {
519 jgs 113 /* write out the data */
520     for (int m=0; m<shape; m++) {
521 jgs 150 fprintf(fileHandle_p, " %e", sampleAvg[m]);
522 jgs 113 }
523     /* pad with zeros */
524 jgs 150 for (int m=0; m<nCompReqd-shape; m++)
525     fprintf(fileHandle_p, " %e", 0.);
526     } else if (nCompReqd == 9) {
527 jgs 113 /* tensor data, so have a 3x3 matrix to output as a row
528     * of 9 data points */
529     int count = 0;
530     for (int m=0; m<shape; m++) {
531     for (int n=0; n<shape; n++) {
532 jgs 150 fprintf(fileHandle_p, " %e", sampleAvg[count]);
533 jgs 113 count++;
534     }
535 jgs 150 for (int n=0; n<3-shape; n++)
536     fprintf(fileHandle_p, " %e", 0.);
537 jgs 113 }
538 jgs 150 for (int m=0; m<3-shape; m++) {
539     for (int n=0; n<3; n++)
540     fprintf(fileHandle_p, " %e", 0.);
541     }
542 jgs 113 }
543     fprintf(fileHandle_p, "\n");
544     }
545 jgs 110 }
546 jgs 113 fprintf(fileHandle_p, "</DataArray>\n");
547     fprintf(fileHandle_p, "</CellData>\n");
548     } else {
549     /* now for the point data */
550     fprintf(fileHandle_p, "<PointData %s>\n", dataTypeStr);
551     fprintf(fileHandle_p, "<DataArray "
552     "Name=\"%s\" "
553     "type=\"Float32\" "
554     "NumberOfComponents=\"%d\" "
555     "format=\"ascii\">\n",
556 jgs 121 dataNameStr, nCompReqd);
557 jgs 150 /* write out the data */
558     /* if the number of required components is more than the number
559     * of actual components, pad with zeros
560     */
561     bool_t do_write=TRUE;
562     int shape = getDataPointShape(data_p, 0);
563     if (shape > 3) {
564     sprintf(error_msg,"shape should be 1, 2, or 3; I got %d\n", shape);
565     Finley_setError(VALUE_ERROR,error_msg);
566     return;
567     }
568 jgs 113 for (i=0; i<mesh_p->Nodes->numNodes; i++) {
569     switch (nodetype) {
570 jgs 150 case(FINLEY_DEGREES_OF_FREEDOM):
571     values = getSampleData(data_p,mesh_p->Nodes->degreeOfFreedom[i]);
572 jgs 113 break;
573     case(FINLEY_REDUCED_DEGREES_OF_FREEDOM):
574     if (mesh_p->Nodes->toReduced[i]>=0) {
575 jgs 150 do_write=TRUE;
576     values = getSampleData(data_p,mesh_p->Nodes->reducedDegreeOfFreedom[i]);
577     } else {
578     do_write=FALSE;
579     }
580 jgs 113 break;
581     case(FINLEY_NODES):
582     values = getSampleData(data_p,i);
583     break;
584     }
585     /* write the data different ways for scalar, vector and tensor */
586 jgs 150 if (do_write) {
587     if (nCompReqd == 1) {
588     fprintf(fileHandle_p, " %e", values[0]);
589 jgs 113 }
590 jgs 150 else if (nCompReqd == 3) {
591     /* write out the data */
592     for (int m=0; m<shape; m++)
593     fprintf(fileHandle_p, " %e", values[m]);
594     /* pad with zeros */
595     for (int m=0; m<nCompReqd-shape; m++)
596     fprintf(fileHandle_p, " %e", 0.);
597 jgs 113 }
598 jgs 150 else if (nCompReqd == 9) {
599     /* tensor data, so have a 3x3 matrix to output as a row
600     * of 9 data points */
601     int count = 0;
602     for (int m=0; m<shape; m++) {
603     for (int n=0; n<shape; n++) {
604     fprintf(fileHandle_p, " %e", values[count]);
605     count++;
606     }
607     for (int n=0; n<3-shape; n++)
608     fprintf(fileHandle_p, " %e", 0.);
609     }
610     for (int m=0; m<3-shape; m++) {
611     for (int n=0; n<3; n++)
612     fprintf(fileHandle_p, " %e", 0.);
613     }
614     }
615     fprintf(fileHandle_p, "\n");
616     }
617 jgs 110 }
618 jgs 113 fprintf(fileHandle_p, "</DataArray>\n");
619     fprintf(fileHandle_p, "</PointData>\n");
620     }
621 jgs 110 }
622    
623 jgs 113 /* finish off the piece */
624     fprintf(fileHandle_p, "</Piece>\n");
625 jgs 110
626     fprintf(fileHandle_p, "</UnstructuredGrid>\n");
627     /* write the xml footer */
628     fprintf(fileHandle_p, "</VTKFile>\n");
629     /* close the file */
630     fclose(fileHandle_p);
631     return;
632     }
633 jgs 121
634     /*
635 jgs 147 * Revision 1.6 2005/08/12 01:45:43 jgs
636     * erge of development branch dev-02 back to main trunk on 2005-08-12
637     *
638 jgs 150 * Revision 1.5.2.4 2005/09/09 08:15:17 gross
639     * bugs in vtk and dx writer fixed
640     *
641     * Revision 1.5.2.3 2005/09/08 08:28:39 gross
642     * some cleanup in savevtk
643     *
644     * Revision 1.5.2.2 2005/09/07 06:26:20 gross
645     * the solver from finley are put into the standalone package paso now
646     *
647 jgs 147 * Revision 1.5.2.1 2005/08/10 06:14:37 gross
648     * QUADRATIC HEXAHEDRON elements fixed
649     *
650 jgs 123 * Revision 1.5 2005/07/08 04:07:55 jgs
651     * Merge of development branch back to main trunk on 2005-07-08
652     *
653 jgs 121 * Revision 1.4 2005/05/06 04:26:15 jgs
654     * Merge of development branch back to main trunk on 2005-05-06
655     *
656 jgs 123 * Revision 1.1.2.7 2005/06/29 02:34:54 gross
657     * some changes towards 64 integers in finley
658     *
659 jgs 121 * Revision 1.1.2.6 2005/05/06 01:17:19 cochrane
660     * Fixed incorrect reporting of number of components in PointData arrays for
661     * vector data.
662     *
663     * Revision 1.1.2.5 2005/05/05 05:38:44 cochrane
664     * Improved formatting of VTK file output.
665     *
666     * Revision 1.1.2.4 2005/02/22 10:03:54 cochrane
667     * Implementation of writing of vtk xml files from finley. This function will
668     * require more testing, but on the cases that I have tried (and with the help
669     * of Lutz and mayavi), it looks like it's producing the correct output. Testing
670     * with more realistic data would be good. I'm at least confident enough to
671     * commit my changes.
672     *
673     * Revision 1.1.2.3 2005/02/17 05:53:26 gross
674     * some bug in saveDX fixed: in fact the bug was in
675     * DataC/getDataPointShape
676     *
677     * Revision 1.1.2.2 2005/02/10 01:34:22 cochrane
678     * Quick fix to make sure that saveVTK compiles so that finley is still buildable. Apologies to those this has affected.
679     *
680     * Revision 1.1.2.1 2005/02/09 06:53:15 cochrane
681     * Initial import to repository. This is the file to implement saving finley/escript meshes out to vtk formatted files. It is basically just a hack of the opendx equivalent, with a lot of the opendx stuff still in the file, so it doesn't actually work just yet, but it probably needs to be added to the cvs.
682     *
683     * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
684     * initial import of project esys2
685     *
686     * Revision 1.1 2004/07/27 08:27:11 gross
687     * Finley: saveDX added: now it is possible to write data on boundary and contact elements
688     *
689     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26