← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to dudley/src/Mesh_write.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_write.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
/*   Dudley: write Mesh */
 
17
 
 
18
/**************************************************************/
 
19
 
 
20
#include "Mesh.h"
 
21
 
 
22
/**************************************************************/
 
23
 
 
24
/*  writes the mesh to the external file fname unsing the Dudley file format: */
 
25
 
 
26
void Dudley_Mesh_write(Dudley_Mesh * in, char *fname)
 
27
{
 
28
    char error_msg[LenErrorMsg_MAX];
 
29
    FILE *f;
 
30
    int NN, i, j, numDim;
 
31
    Dudley_TagMap *tag_map = in->TagMap;
 
32
 
 
33
    if (in->MPIInfo->size > 1)
 
34
    {
 
35
        Dudley_setError(IO_ERROR, "Mesh_write: only single processor runs are supported.");
 
36
        return;
 
37
 
 
38
    }
 
39
    /* open file */
 
40
    f = fopen(fname, "w");
 
41
    if (f == NULL)
 
42
    {
 
43
        sprintf(error_msg, "Mesh_write: Opening file %s for writing failed.", fname);
 
44
        Dudley_setError(IO_ERROR, error_msg);
 
45
        return;
 
46
    }
 
47
 
 
48
    /* write header */
 
49
 
 
50
    fprintf(f, "%s\n", in->Name);
 
51
 
 
52
    /*  write nodes: */
 
53
 
 
54
    if (in->Nodes != NULL)
 
55
    {
 
56
        numDim = Dudley_Mesh_getDim(in);
 
57
        fprintf(f, "%1dD-Nodes %d\n", numDim, in->Nodes->numNodes);
 
58
        for (i = 0; i < in->Nodes->numNodes; i++)
 
59
        {
 
60
            fprintf(f, "%d %d %d", in->Nodes->Id[i], in->Nodes->globalDegreesOfFreedom[i], in->Nodes->Tag[i]);
 
61
            for (j = 0; j < numDim; j++)
 
62
                fprintf(f, " %20.15e", in->Nodes->Coordinates[INDEX2(j, i, numDim)]);
 
63
            fprintf(f, "\n");
 
64
        }
 
65
    }
 
66
    else
 
67
    {
 
68
        fprintf(f, "0D-Nodes 0\n");
 
69
    }
 
70
 
 
71
    /*  write elements: */
 
72
 
 
73
    if (in->Elements != NULL)
 
74
    {
 
75
        fprintf(f, "%s %d\n", in->Elements->ename /*referenceElementSet->referenceElement->Type->Name */ ,
 
76
                in->Elements->numElements);
 
77
        NN = in->Elements->numNodes;
 
78
        for (i = 0; i < in->Elements->numElements; i++)
 
79
        {
 
80
            fprintf(f, "%d %d", in->Elements->Id[i], in->Elements->Tag[i]);
 
81
            for (j = 0; j < NN; j++)
 
82
                fprintf(f, " %d", in->Nodes->Id[in->Elements->Nodes[INDEX2(j, i, NN)]]);
 
83
            fprintf(f, "\n");
 
84
        }
 
85
    }
 
86
    else
 
87
    {
 
88
        fprintf(f, "Tet4 0\n");
 
89
    }
 
90
 
 
91
    /*  write face elements: */
 
92
    if (in->FaceElements != NULL)
 
93
    {
 
94
        fprintf(f, "%s %d\n", in->FaceElements->ename /*referenceElementSet->referenceElement->Type->Name */ ,
 
95
                in->FaceElements->numElements);
 
96
        NN = in->FaceElements->numNodes;
 
97
        for (i = 0; i < in->FaceElements->numElements; i++)
 
98
        {
 
99
            fprintf(f, "%d %d", in->FaceElements->Id[i], in->FaceElements->Tag[i]);
 
100
            for (j = 0; j < NN; j++)
 
101
                fprintf(f, " %d", in->Nodes->Id[in->FaceElements->Nodes[INDEX2(j, i, NN)]]);
 
102
            fprintf(f, "\n");
 
103
        }
 
104
    }
 
105
    else
 
106
    {
 
107
        fprintf(f, "Tri3 0\n");
 
108
    }
 
109
 
 
110
    /*  write points: */
 
111
    if (in->Points != NULL)
 
112
    {
 
113
        fprintf(f, "%s %d\n", in->Points->ename /*referenceElementSet->referenceElement->Type->Name */ ,
 
114
                in->Points->numElements);
 
115
        for (i = 0; i < in->Points->numElements; i++)
 
116
        {
 
117
            fprintf(f, "%d %d %d\n", in->Points->Id[i], in->Points->Tag[i],
 
118
                    in->Nodes->Id[in->Points->Nodes[INDEX2(0, i, 1)]]);
 
119
        }
 
120
    }
 
121
    else
 
122
    {
 
123
        fprintf(f, "Point1 0\n");
 
124
    }
 
125
 
 
126
    /*  write tags: */
 
127
    if (tag_map)
 
128
    {
 
129
        fprintf(f, "Tags\n");
 
130
        while (tag_map)
 
131
        {
 
132
            fprintf(f, "%s %d\n", tag_map->name, tag_map->tag_key);
 
133
            tag_map = tag_map->next;
 
134
        }
 
135
    }
 
136
    fclose(f);
 
137
#ifdef Dudley_TRACE
 
138
    printf("mesh %s has been written to file %s\n", in->Name, fname);
 
139
#endif
 
140
}
 
141
 
 
142
void Dudley_PrintMesh_Info(Dudley_Mesh * in, bool_t full)
 
143
{
 
144
    int NN, i, j, numDim;
 
145
    Dudley_TagMap *tag_map = in->TagMap;
 
146
 
 
147
    fprintf(stdout, "Dudley_PrintMesh_Info running on CPU %d of %d\n", in->MPIInfo->rank, in->MPIInfo->size);
 
148
    fprintf(stdout, "\tMesh name '%s'\n", in->Name);
 
149
    fprintf(stdout, "\tApproximation order %d\n", in->approximationOrder);
 
150
    fprintf(stdout, "\tReduced Approximation order %d\n", in->reducedApproximationOrder);
 
151
    fprintf(stdout, "\tIntegration order %d\n", in->integrationOrder);
 
152
    fprintf(stdout, "\tReduced Integration order %d\n", in->reducedIntegrationOrder);
 
153
 
 
154
    /* write nodes: */
 
155
    if (in->Nodes != NULL)
 
156
    {
 
157
        numDim = Dudley_Mesh_getDim(in);
 
158
        fprintf(stdout, "\tNodes: %1dD-Nodes %d\n", numDim, in->Nodes->numNodes);
 
159
        if (full)
 
160
        {
 
161
            fprintf(stdout, "\t     Id   Tag  gDOF   gNI grDfI  grNI:  Coordinates\n");
 
162
            for (i = 0; i < in->Nodes->numNodes; i++)
 
163
            {
 
164
                fprintf(stdout, "\t  %5d %5d %5d %5d %5d %5d: ", in->Nodes->Id[i], in->Nodes->Tag[i],
 
165
                        in->Nodes->globalDegreesOfFreedom[i], in->Nodes->globalNodesIndex[i],
 
166
                        in->Nodes->globalReducedDOFIndex[i], in->Nodes->globalReducedNodesIndex[i]);
 
167
                for (j = 0; j < numDim; j++)
 
168
                    fprintf(stdout, " %20.15e", in->Nodes->Coordinates[INDEX2(j, i, numDim)]);
 
169
                fprintf(stdout, "\n");
 
170
            }
 
171
        }
 
172
    }
 
173
    else
 
174
    {
 
175
        fprintf(stdout, "\tNodes: 0D-Nodes 0\n");
 
176
    }
 
177
 
 
178
    /* write elements: */
 
179
    if (in->Elements != NULL)
 
180
    {
 
181
        int mine = 0, overlap = 0;
 
182
        for (i = 0; i < in->Elements->numElements; i++)
 
183
        {
 
184
            if (in->Elements->Owner[i] == in->MPIInfo->rank)
 
185
                mine++;
 
186
            else
 
187
                overlap++;
 
188
        }
 
189
        fprintf(stdout, "\tElements: %s %d (TypeId=%d) owner=%d overlap=%d\n",
 
190
                in->Elements->ename /*referenceElementSet->referenceElement->Type->Name */ , in->Elements->numElements,
 
191
                in->Elements->etype /*referenceElementSet->referenceElement->Type->TypeId */ , mine, overlap);
 
192
        NN = in->Elements->numNodes;
 
193
        if (full)
 
194
        {
 
195
            fprintf(stdout, "\t     Id   Tag Owner Color:  Nodes\n");
 
196
            for (i = 0; i < in->Elements->numElements; i++)
 
197
            {
 
198
                fprintf(stdout, "\t  %5d %5d %5d %5d: ", in->Elements->Id[i], in->Elements->Tag[i],
 
199
                        in->Elements->Owner[i], in->Elements->Color[i]);
 
200
                for (j = 0; j < NN; j++)
 
201
                    fprintf(stdout, " %5d", in->Nodes->Id[in->Elements->Nodes[INDEX2(j, i, NN)]]);
 
202
                fprintf(stdout, "\n");
 
203
            }
 
204
        }
 
205
    }
 
206
    else
 
207
    {
 
208
        fprintf(stdout, "\tElements: Tet4 0\n");
 
209
    }
 
210
 
 
211
    /* write face elements: */
 
212
    if (in->FaceElements != NULL)
 
213
    {
 
214
        int mine = 0, overlap = 0;
 
215
        for (i = 0; i < in->FaceElements->numElements; i++)
 
216
        {
 
217
            if (in->FaceElements->Owner[i] == in->MPIInfo->rank)
 
218
                mine++;
 
219
            else
 
220
                overlap++;
 
221
        }
 
222
        fprintf(stdout, "\tFace elements: %s %d (TypeId=%d) owner=%d overlap=%d\n",
 
223
                in->FaceElements->ename /*referenceElementSet->referenceElement->Type->Name */ ,
 
224
                in->FaceElements->numElements,
 
225
                in->FaceElements->etype /*->referenceElementSet->referenceElement->Type->TypeId*/ , mine, overlap);
 
226
        NN = in->FaceElements->numNodes;
 
227
        if (full)
 
228
        {
 
229
            fprintf(stdout, "\t     Id   Tag Owner Color:  Nodes\n");
 
230
            for (i = 0; i < in->FaceElements->numElements; i++)
 
231
            {
 
232
                fprintf(stdout, "\t  %5d %5d %5d %5d: ", in->FaceElements->Id[i], in->FaceElements->Tag[i],
 
233
                        in->FaceElements->Owner[i], in->FaceElements->Color[i]);
 
234
                for (j = 0; j < NN; j++)
 
235
                    fprintf(stdout, " %5d", in->Nodes->Id[in->FaceElements->Nodes[INDEX2(j, i, NN)]]);
 
236
                fprintf(stdout, "\n");
 
237
            }
 
238
        }
 
239
    }
 
240
    else
 
241
    {
 
242
        fprintf(stdout, "\tFace elements: Tri3 0\n");
 
243
    }
 
244
 
 
245
    /* write points: */
 
246
    if (in->Points != NULL)
 
247
    {
 
248
        int mine = 0, overlap = 0;
 
249
        for (i = 0; i < in->Points->numElements; i++)
 
250
        {
 
251
            if (in->Points->Owner[i] == in->MPIInfo->rank)
 
252
                mine++;
 
253
            else
 
254
                overlap++;
 
255
        }
 
256
        fprintf(stdout, "\tPoints: %s %d (TypeId=%d) owner=%d overlap=%d\n",
 
257
                in->Points->ename /*->referenceElementSet->referenceElement->Type->Name*/ , in->Points->numElements,
 
258
                in->Points->etype /*referenceElementSet->referenceElement->Type->TypeId */ , mine, overlap);
 
259
        if (full)
 
260
        {
 
261
            fprintf(stdout, "\t     Id   Tag Owner Color:  Nodes\n");
 
262
            for (i = 0; i < in->Points->numElements; i++)
 
263
            {
 
264
                fprintf(stdout, "\t  %5d %5d %5d %5d %5d\n", in->Points->Id[i], in->Points->Tag[i],
 
265
                        in->Points->Owner[i], in->Points->Color[i], in->Nodes->Id[in->Points->Nodes[INDEX2(0, i, 1)]]);
 
266
            }
 
267
        }
 
268
    }
 
269
    else
 
270
    {
 
271
        fprintf(stdout, "\tPoints: Point1 0\n");
 
272
    }
 
273
 
 
274
    /* write tags: */
 
275
    if (tag_map)
 
276
    {
 
277
        fprintf(stdout, "\tTags:\n");
 
278
        while (tag_map)
 
279
        {
 
280
            fprintf(stdout, "\t  %5d %s\n", tag_map->tag_key, tag_map->name);
 
281
            tag_map = tag_map->next;
 
282
        }
 
283
    }
 
284
}