← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to dudley/src/Mesh_saveDX.c

  • Committer: jfenwick
  • Date: 2010-10-11 01:48:14 UTC
  • Revision ID: svn-v4:77569008-7704-0410-b7a0-a92fef0b09fd:trunk:3259
Merging dudley and scons updates from branches

Show diffs side-by-side

added added

removed removed

Lines of Context:
dudley/src/Mesh_saveDX.c
 
1
 
 
2
/*******************************************************
 
3
*
 
4
* Copyright (c) 2003-2010 by University of Queensland
 
5
* Earth Systems Science Computational Center (ESSCC)
 
6
* http://www.uq.edu.au/esscc
 
7
*
 
8
* Primary Business: Queensland, Australia
 
9
* Licensed under the Open Software License version 3.0
 
10
* http://www.opensource.org/licenses/osl-3.0.php
 
11
*
 
12
*******************************************************/
 
13
 
 
14
/**************************************************************/
 
15
 
 
16
/*   writes data and mesh in an opendx file                      */
 
17
/*   the input data needs to be cell centered or on reducedNodes */
 
18
 
 
19
/**************************************************************/
 
20
 
 
21
#include "Mesh.h"
 
22
#include "Assemble.h"
 
23
 
 
24
/**************************************************************/
 
25
 
 
26
void Dudley_Mesh_saveDX(const char *filename_p, Dudley_Mesh * mesh_p, const dim_t num_data, char **names_p,
 
27
                        escriptDataC * *data_pp)
 
28
{
 
29
    char error_msg[LenErrorMsg_MAX], elemTypeStr[32];
 
30
    /* some tables needed for reordering */
 
31
    int resort[6][9] = {
 
32
        {0, 1},                 /* line */
 
33
        {0, 1, 2},              /* triagle */
 
34
        {0, 1, 2, 3},           /* tetrahedron */
 
35
        {0, 3, 1, 2},           /* quadrilateral */
 
36
        {3, 0, 7, 4, 2, 1, 6, 5},       /* hexahedron */
 
37
    };
 
38
    FILE *fileHandle_p = NULL;
 
39
    int i, j, k, i_data, elementtype, numPoints = 0, nDim, *resortIndex = NULL, p,
 
40
        numDXNodesPerElement = 0, numCells, NN, object_count, rank, nComp, numPointsPerSample;
 
41
    __const double *values;
 
42
    double rtmp;
 
43
    bool_t *isCellCentered = NULL;
 
44
    Dudley_ElementFile *elements = NULL;
 
45
    Dudley_ElementTypeId TypeId;
 
46
    /* open the file  and check handel */
 
47
 
 
48
    /* if there is no mesh we just return */
 
49
    if (mesh_p == NULL)
 
50
        return;
 
51
    isCellCentered = MEMALLOC(num_data, bool_t);
 
52
    if (Dudley_checkPtr(isCellCentered))
 
53
        return;
 
54
 
 
55
    fileHandle_p = fopen(filename_p, "w");
 
56
    if (fileHandle_p == NULL)
 
57
    {
 
58
        sprintf(error_msg, "File %s could not be opened for writing.", filename_p);
 
59
        MEMFREE(isCellCentered);
 
60
        fclose(fileHandle_p);
 
61
        Dudley_setError(IO_ERROR, error_msg);
 
62
        return;
 
63
    }
 
64
    /* find the mesh type to be written */
 
65
    elementtype = DUDLEY_UNKNOWN;
 
66
    for (i_data = 0; i_data < num_data; ++i_data)
 
67
    {
 
68
        if (!isEmpty(data_pp[i_data]))
 
69
        {
 
70
            switch (getFunctionSpaceType(data_pp[i_data]))
 
71
            {
 
72
            case DUDLEY_REDUCED_NODES:
 
73
                if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_ELEMENTS)
 
74
                {
 
75
                    elementtype = DUDLEY_ELEMENTS;
 
76
                }
 
77
                else
 
78
                {
 
79
                    Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
 
80
                    MEMFREE(isCellCentered);
 
81
                    fclose(fileHandle_p);
 
82
                    return;
 
83
                }
 
84
                isCellCentered[i_data] = FALSE;
 
85
                break;
 
86
            case DUDLEY_ELEMENTS:
 
87
            case DUDLEY_REDUCED_ELEMENTS:
 
88
                if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_ELEMENTS)
 
89
                {
 
90
                    elementtype = DUDLEY_ELEMENTS;
 
91
                }
 
92
                else
 
93
                {
 
94
                    Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
 
95
                    MEMFREE(isCellCentered);
 
96
                    fclose(fileHandle_p);
 
97
                    return;
 
98
                }
 
99
                isCellCentered[i_data] = TRUE;
 
100
                break;
 
101
            case DUDLEY_FACE_ELEMENTS:
 
102
            case DUDLEY_REDUCED_FACE_ELEMENTS:
 
103
                if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_FACE_ELEMENTS)
 
104
                {
 
105
                    elementtype = DUDLEY_FACE_ELEMENTS;
 
106
                }
 
107
                else
 
108
                {
 
109
                    Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
 
110
                    MEMFREE(isCellCentered);
 
111
                    fclose(fileHandle_p);
 
112
                    return;
 
113
                }
 
114
                isCellCentered[i_data] = TRUE;
 
115
                break;
 
116
            case DUDLEY_POINTS:
 
117
                if (elementtype == DUDLEY_UNKNOWN || elementtype == DUDLEY_POINTS)
 
118
                {
 
119
                    elementtype = DUDLEY_POINTS;
 
120
                }
 
121
                else
 
122
                {
 
123
                    Dudley_setError(TYPE_ERROR, "saveDX: cannot write given data into single file.");
 
124
                    MEMFREE(isCellCentered);
 
125
                    fclose(fileHandle_p);
 
126
                    return;
 
127
                }
 
128
                isCellCentered[i_data] = TRUE;
 
129
                break;
 
130
            default:
 
131
                sprintf(error_msg, "saveDX: I don't know how to handle function space type %d",
 
132
                        getFunctionSpaceType(data_pp[i_data]));
 
133
                Dudley_setError(TYPE_ERROR, error_msg);
 
134
                MEMFREE(isCellCentered);
 
135
                fclose(fileHandle_p);
 
136
                return;
 
137
            }
 
138
        }
 
139
    }
 
140
    /* select number of points and the mesh component */
 
141
    numPoints = mesh_p->Nodes->reducedNodesMapping->numTargets;
 
142
    nDim = mesh_p->Nodes->numDim;
 
143
    if (elementtype == DUDLEY_UNKNOWN)
 
144
        elementtype = DUDLEY_ELEMENTS;
 
145
    switch (elementtype)
 
146
    {
 
147
    case DUDLEY_ELEMENTS:
 
148
        elements = mesh_p->Elements;
 
149
        break;
 
150
    case DUDLEY_FACE_ELEMENTS:
 
151
        elements = mesh_p->FaceElements;
 
152
        break;
 
153
    case DUDLEY_POINTS:
 
154
        elements = mesh_p->Points;
 
155
        break;
 
156
    }
 
157
    if (elements == NULL)
 
158
    {
 
159
        Dudley_setError(SYSTEM_ERROR, "saveDX: undefined element file");
 
160
        MEMFREE(isCellCentered);
 
161
        fclose(fileHandle_p);
 
162
        return;
 
163
    }
 
164
 
 
165
    /* map dudley element type to DX element type */
 
166
    TypeId = elements->etype;
 
167
    numDXNodesPerElement = 0;
 
168
    numCells = elements->numElements;
 
169
    if (TypeId == Dudley_Line2)
 
170
    {
 
171
        numDXNodesPerElement = 2;
 
172
        resortIndex = resort[0];
 
173
        strcpy(elemTypeStr, "lines");
 
174
    }
 
175
    else if (TypeId == Dudley_Tri3)
 
176
    {
 
177
        numDXNodesPerElement = 3;
 
178
        resortIndex = resort[1];
 
179
        strcpy(elemTypeStr, "triangles");
 
180
/*   } else if (TypeId==Rec4 || TypeId==Rec8 || TypeId==Rec9 || TypeId==Rec12 || TypeId==Rec16 || TypeId==Rec9Macro ) {
 
181
     numDXNodesPerElement = 4;
 
182
     resortIndex=resort[3];
 
183
     strcpy(elemTypeStr, "quads");*/
 
184
    }
 
185
    else if (TypeId == Dudley_Tet4)
 
186
    {
 
187
        numDXNodesPerElement = 4;
 
188
        resortIndex = resort[2];
 
189
        strcpy(elemTypeStr, "tetrahedra");
 
190
/*   } else if (TypeId==Hex8 || TypeId==Hex20 || TypeId==Hex32 ) {
 
191
     numDXNodesPerElement = 8;
 
192
     resortIndex=resort[4];
 
193
     strcpy(elemTypeStr, "cubes");*/
 
194
    }
 
195
    else
 
196
    {
 
197
        sprintf(error_msg, "saveDX: Element type %s is not supported by DX", elements->ename);
 
198
        Dudley_setError(VALUE_ERROR, error_msg);
 
199
        MEMFREE(isCellCentered);
 
200
        fclose(fileHandle_p);
 
201
        return;
 
202
    }
 
203
 
 
204
    /* positions */
 
205
    fprintf(fileHandle_p, "object 1 class array type float rank 1 shape %d items %d data follows\n", nDim, numPoints);
 
206
    for (i = 0; i < numPoints; i++)
 
207
    {
 
208
        p = mesh_p->Nodes->reducedNodesMapping->map[i];
 
209
        for (j = 0; j < nDim; j++)
 
210
            fprintf(fileHandle_p, " %g", mesh_p->Nodes->Coordinates[INDEX2(j, p, nDim)]);
 
211
        fprintf(fileHandle_p, "\n");
 
212
    }
 
213
    /* connection table */
 
214
    NN = elements->numNodes;
 
215
    fprintf(fileHandle_p, "object 2 class array type int rank 1 shape %d items %d data follows\n", numDXNodesPerElement,
 
216
            numCells);
 
217
    for (i = 0; i < numCells; i++)
 
218
    {
 
219
        for (j = 0; j < numDXNodesPerElement; j++)
 
220
            fprintf(fileHandle_p, " %d",
 
221
                    mesh_p->Nodes->reducedNodesMapping->target[elements->Nodes[INDEX2(resortIndex[j], i, NN)]]);
 
222
        fprintf(fileHandle_p, "\n");
 
223
    }
 
224
    fprintf(fileHandle_p, "attribute \"element type\" string \"%s\"\n", elemTypeStr);
 
225
    fprintf(fileHandle_p, "attribute \"ref\" string \"positions\"\n");
 
226
 
 
227
    /* data */
 
228
    object_count = 2;
 
229
    for (i_data = 0; i_data < num_data; ++i_data)
 
230
    {
 
231
        if (!isEmpty(data_pp[i_data]))
 
232
        {
 
233
            object_count++;
 
234
            rank = getDataPointRank(data_pp[i_data]);
 
235
            nComp = getDataPointSize(data_pp[i_data]);
 
236
            fprintf(fileHandle_p, "object %d class array type float rank %d ", object_count, rank);
 
237
            if (0 < rank)
 
238
            {
 
239
                fprintf(fileHandle_p, "shape ");
 
240
                for (i = 0; i < rank; i++)
 
241
                    fprintf(fileHandle_p, "%d ", getDataPointShape(data_pp[i_data], i));
 
242
            }
 
243
            if (isCellCentered[i_data])
 
244
            {
 
245
                if (Dudley_Assemble_reducedIntegrationOrder(data_pp[i_data]))
 
246
                {
 
247
                    numPointsPerSample = (elements->numLocalDim == 0) ? 0 : 1;
 
248
                }
 
249
                else
 
250
                {
 
251
                    numPointsPerSample = (elements->numLocalDim == 0) ? 0 : (elements->numLocalDim + 1);
 
252
                }
 
253
                if (numPointsPerSample > 0)
 
254
                {
 
255
                    fprintf(fileHandle_p, "items %d data follows\n", numCells);
 
256
                    for (i = 0; i < elements->numElements; i++)
 
257
                    {
 
258
                        values = getSampleDataRO(data_pp[i_data], i);
 
259
                        for (k = 0; k < nComp; k++)
 
260
                        {
 
261
                            if (isExpanded(data_pp[i_data]))
 
262
                            {
 
263
                                rtmp = 0.;
 
264
                                for (j = 0; j < numPointsPerSample; j++)
 
265
                                    rtmp += values[INDEX2(k, j, nComp)];
 
266
                                fprintf(fileHandle_p, " %g", rtmp / numPointsPerSample);
 
267
                            }
 
268
                            else
 
269
                            {
 
270
                                fprintf(fileHandle_p, " %g", values[k]);
 
271
                            }
 
272
                        }
 
273
                        fprintf(fileHandle_p, "\n");
 
274
                    }
 
275
                    fprintf(fileHandle_p, "attribute \"dep\" string \"connections\"\n");
 
276
                }
 
277
            }
 
278
            else
 
279
            {
 
280
                fprintf(fileHandle_p, "items %d data follows\n", numPoints);
 
281
                for (i = 0; i < numPoints; i++)
 
282
                {
 
283
                    values = getSampleDataRO(data_pp[i_data], i);
 
284
                    for (k = 0; k < nComp; k++)
 
285
                        fprintf(fileHandle_p, " %g", values[k]);
 
286
                    fprintf(fileHandle_p, "\n");
 
287
                }
 
288
                fprintf(fileHandle_p, "attribute \"dep\" string \"positions\"\n");
 
289
            }
 
290
        }
 
291
    }
 
292
 
 
293
    /* and finish it up */
 
294
    if (num_data == 0)
 
295
    {
 
296
        fprintf(fileHandle_p, "object %d class field\n", object_count + 1);
 
297
        fprintf(fileHandle_p, "component \"positions\" value 1\n");
 
298
        fprintf(fileHandle_p, "component \"connections\" value 2\n");
 
299
    }
 
300
    else
 
301
    {
 
302
        object_count = 2;
 
303
        for (i_data = 0; i_data < num_data; ++i_data)
 
304
        {
 
305
            if (!isEmpty(data_pp[i_data]))
 
306
            {
 
307
                object_count++;
 
308
                fprintf(fileHandle_p, "object \"%s\" class field\n", names_p[i_data]);
 
309
                fprintf(fileHandle_p, "component \"positions\" value 1\n");
 
310
                fprintf(fileHandle_p, "component \"connections\" value 2\n");
 
311
                fprintf(fileHandle_p, "component \"data\" value %d\n", object_count);
 
312
            }
 
313
        }
 
314
    }
 
315
    fprintf(fileHandle_p, "end\n");
 
316
    /* close the file */
 
317
    fclose(fileHandle_p);
 
318
    MEMFREE(isCellCentered);
 
319
    return;
 
320
}