← Back to branch summary

~ubuntu-branches/ubuntu/trusty/silo-llnl/trusty

« back to all changes in this revision

Viewing changes to tests/efcentering.c

  • Committer: Bazaar Package Importer
  • Author(s): Alastair McKinstry
  • Date: 2011-01-02 00:03:01 UTC
  • Revision ID: james.westby@ubuntu.com-20110102000301-9s2hfsjrkguz6h4r
Tags: upstream-4.8
Import upstream version 4.8

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/efcentering.c
 
1
/*
 
2
Copyright (c) 1994 - 2010, Lawrence Livermore National Security, LLC.
 
3
LLNL-CODE-425250.
 
4
All rights reserved.
 
5
 
 
6
This file is part of Silo. For details, see silo.llnl.gov.
 
7
 
 
8
Redistribution and use in source and binary forms, with or without
 
9
modification, are permitted provided that the following conditions
 
10
are met:
 
11
 
 
12
   * Redistributions of source code must retain the above copyright
 
13
     notice, this list of conditions and the disclaimer below.
 
14
   * Redistributions in binary form must reproduce the above copyright
 
15
     notice, this list of conditions and the disclaimer (as noted
 
16
     below) in the documentation and/or other materials provided with
 
17
     the distribution.
 
18
   * Neither the name of the LLNS/LLNL nor the names of its
 
19
     contributors may be used to endorse or promote products derived
 
20
     from this software without specific prior written permission.
 
21
 
 
22
THIS SOFTWARE  IS PROVIDED BY  THE COPYRIGHT HOLDERS  AND CONTRIBUTORS
 
23
"AS  IS" AND  ANY EXPRESS  OR IMPLIED  WARRANTIES, INCLUDING,  BUT NOT
 
24
LIMITED TO, THE IMPLIED  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 
25
A  PARTICULAR  PURPOSE ARE  DISCLAIMED.  IN  NO  EVENT SHALL  LAWRENCE
 
26
LIVERMORE  NATIONAL SECURITY, LLC,  THE U.S.  DEPARTMENT OF  ENERGY OR
 
27
CONTRIBUTORS BE LIABLE FOR  ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 
28
EXEMPLARY, OR  CONSEQUENTIAL DAMAGES  (INCLUDING, BUT NOT  LIMITED TO,
 
29
PROCUREMENT OF  SUBSTITUTE GOODS  OR SERVICES; LOSS  OF USE,  DATA, OR
 
30
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 
31
LIABILITY, WHETHER  IN CONTRACT, STRICT LIABILITY,  OR TORT (INCLUDING
 
32
NEGLIGENCE OR  OTHERWISE) ARISING IN  ANY WAY OUT  OF THE USE  OF THIS
 
33
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
34
 
 
35
This work was produced at Lawrence Livermore National Laboratory under
 
36
Contract  No.   DE-AC52-07NA27344 with  the  DOE.  Neither the  United
 
37
States Government  nor Lawrence  Livermore National Security,  LLC nor
 
38
any of  their employees,  makes any warranty,  express or  implied, or
 
39
assumes   any   liability   or   responsibility  for   the   accuracy,
 
40
completeness, or usefulness of any information, apparatus, product, or
 
41
process  disclosed, or  represents  that its  use  would not  infringe
 
42
privately-owned   rights.  Any  reference   herein  to   any  specific
 
43
commercial products,  process, or  services by trade  name, trademark,
 
44
manufacturer or otherwise does not necessarily constitute or imply its
 
45
endorsement,  recommendation,   or  favoring  by   the  United  States
 
46
Government or Lawrence Livermore National Security, LLC. The views and
 
47
opinions  of authors  expressed  herein do  not  necessarily state  or
 
48
reflect those  of the United  States Government or  Lawrence Livermore
 
49
National  Security, LLC,  and shall  not  be used  for advertising  or
 
50
product endorsement purposes.
 
51
*/
 
52
#include <silo.h>
 
53
#include <stdio.h>
 
54
#include <stdlib.h>
 
55
#include <std.c>
 
56
 
 
57
void
 
58
add_edge(int nid1, int nid2, int *nedges, int *maxedges, int **edges)
 
59
{
 
60
    if (*nedges == *maxedges)
 
61
    {
 
62
        *maxedges = (*maxedges) * 2 + 1;
 
63
        *edges = (int *) realloc(*edges, *maxedges * 2 * sizeof(int));
 
64
    }
 
65
    (*edges)[2*(*nedges)+0] = nid1;
 
66
    (*edges)[2*(*nedges)+1] = nid2;
 
67
    *nedges = (*nedges) + 1;
 
68
}
 
69
 
 
70
int
 
71
have_edge(int nid1, int nid2, int nedges, const int *edges)
 
72
{
 
73
    int i;
 
74
    for (i = 0; i < nedges; i++)
 
75
    {
 
76
        if ((edges[2*i+0] == nid1 &&
 
77
             edges[2*i+1] == nid2) ||
 
78
            (edges[2*i+0] == nid2 &&
 
79
             edges[2*i+1] == nid1))
 
80
            return 1;
 
81
    }
 
82
    return 0;
 
83
}
 
84
 
 
85
#define HANDLE_EDGE(A,B)                                                \
 
86
    if (!have_edge(nodelist[nlidx+A],nodelist[nlidx+B],*nedges,*edges)) \
 
87
        add_edge(nodelist[nlidx+A],nodelist[nlidx+B],nedges,&maxedges,edges)
 
88
 
 
89
void
 
90
build_edgelist(int nzones, int ndims, const int *nodelist,
 
91
    int lnodelist, int origin, int lo_off, int hi_off,
 
92
    const int *st, const int *sz, const int *sc, int nshapes,
 
93
    int *nedges, int **edges)
 
94
{
 
95
    int nlidx = 0;
 
96
    int maxedges = 0;
 
97
    int shape;
 
98
 
 
99
    *nedges = 0;
 
100
    *edges = 0;
 
101
    for (shape = 0; shape < nshapes; shape++)
 
102
    {
 
103
        int zonetype = st[shape];
 
104
        int zonesize = sz[shape];
 
105
        int zone;
 
106
        for (zone = 0; zone < sc[shape]; zone++)
 
107
        {
 
108
            switch (zonetype)
 
109
            {
 
110
                case DB_ZONETYPE_HEX:
 
111
                {
 
112
                    HANDLE_EDGE(0,1);
 
113
                    HANDLE_EDGE(0,3);
 
114
                    HANDLE_EDGE(0,4);
 
115
                    HANDLE_EDGE(1,2);
 
116
                    HANDLE_EDGE(1,5);
 
117
                    HANDLE_EDGE(2,3);
 
118
                    HANDLE_EDGE(2,6);
 
119
                    HANDLE_EDGE(3,7);
 
120
                    HANDLE_EDGE(4,5);
 
121
                    HANDLE_EDGE(4,7);
 
122
                    HANDLE_EDGE(5,6);
 
123
                    HANDLE_EDGE(6,7);
 
124
                    break;
 
125
                }
 
126
            }
 
127
            nlidx += zonesize;
 
128
        }
 
129
    }
 
130
}
 
131
 
 
132
int compare_nodes(const void *nap, const void *nbp)
 
133
{
 
134
    int na = *((int*)nap);
 
135
    int nb = *((int*)nbp);
 
136
    if (na < nb) return -1;
 
137
    else if (na > nb) return 1;
 
138
    return 0;
 
139
}
 
140
 
 
141
void
 
142
add_face(int nid1, int nid2, int nid3, int nid4, int *nfaces, int *maxfaces, int **faces)
 
143
{
 
144
    int i;
 
145
    int fnid[4];
 
146
    fnid[0] = nid1; 
 
147
    fnid[1] = nid2; 
 
148
    fnid[2] = nid3; 
 
149
    fnid[3] = nid4; 
 
150
    qsort(fnid, sizeof(int), 4, compare_nodes);
 
151
    if (*nfaces == *maxfaces)
 
152
    {
 
153
        *maxfaces = (*maxfaces) * 2 + 1;
 
154
        *faces = (int *) realloc(*faces, *maxfaces * 4 * sizeof(int));
 
155
    }
 
156
    for (i = 0; i < 4; i++)
 
157
        (*faces)[4*(*nfaces)+i] = fnid[i];
 
158
    *nfaces = (*nfaces) + 1;
 
159
}
 
160
 
 
161
int
 
162
have_face(int nid1, int nid2, int nid3, int nid4, int nfaces, const int *faces)
 
163
{
 
164
    int i,j;
 
165
    int fnid[4];
 
166
    fnid[0] = nid1; 
 
167
    fnid[1] = nid2; 
 
168
    fnid[2] = nid3; 
 
169
    fnid[3] = nid4; 
 
170
    qsort(fnid, sizeof(int), 4, compare_nodes);
 
171
    for (i = 0; i < nfaces; i++)
 
172
    {
 
173
        int allmatch = 1;
 
174
        for (j = 0; j < 4; j++)
 
175
        {
 
176
            if (faces[4*i+j] != fnid[j])
 
177
            {
 
178
                allmatch = 0;
 
179
                break;
 
180
            }
 
181
        }
 
182
        if (allmatch) return 1;
 
183
    }
 
184
    return 0;
 
185
}
 
186
 
 
187
#define HANDLE_FACE(A,B,C,D)                                            \
 
188
    if (!have_face(nodelist[nlidx+A],nodelist[nlidx+B],                 \
 
189
                   nodelist[nlidx+C],nodelist[nlidx+D],*nfaces,*faces)) \
 
190
        add_face(nodelist[nlidx+A],nodelist[nlidx+B],                   \
 
191
                 nodelist[nlidx+C],nodelist[nlidx+D],nfaces,&maxfaces,faces)
 
192
 
 
193
void
 
194
build_facelist(int nzones, int ndims, const int *nodelist,
 
195
    int lnodelist, int origin, int lo_off, int hi_off,
 
196
    const int *st, const int *sz, const int *sc, int nshapes,
 
197
    int *nfaces, int **faces)
 
198
{
 
199
    int nlidx = 0;
 
200
    int maxfaces = 0;
 
201
    int shape;
 
202
 
 
203
    *nfaces = 0;
 
204
    *faces = 0;
 
205
    for (shape = 0; shape < nshapes; shape++)
 
206
    {
 
207
        int zonetype = st[shape];
 
208
        int zonesize = sz[shape];
 
209
        int zone;
 
210
        for (zone = 0; zone < sc[shape]; zone++)
 
211
        {
 
212
            switch (zonetype)
 
213
            {
 
214
                case DB_ZONETYPE_HEX:
 
215
                {
 
216
                    HANDLE_FACE(0,1,5,4);
 
217
                    HANDLE_FACE(0,3,2,1);
 
218
                    HANDLE_FACE(0,4,7,3);
 
219
                    HANDLE_FACE(1,2,6,5);
 
220
                    HANDLE_FACE(2,3,7,6);
 
221
                    HANDLE_FACE(4,5,6,7);
 
222
                    break;
 
223
                }
 
224
            }
 
225
            nlidx += zonesize;
 
226
        }
 
227
    }
 
228
}
 
229
 
 
230
int
 
231
main(int argc, char *argv[])
 
232
{
 
233
    DBfile         *dbfile = NULL;
 
234
    char           *coordnames[3];
 
235
    float          *coords[3];
 
236
    float           x[64], y[64], z[64];
 
237
    int             shapesize[1];
 
238
    int             shapecnt[1];
 
239
    DBfacelist     *facelist = NULL;
 
240
    int             matnos[1], matlist[1], dims[3];
 
241
    int             i, j, k, len;
 
242
    float           evar2d[2*16], evar3d[3*64], fvar3d[3*64];
 
243
    int             driver = DB_PDB;
 
244
    char            *filename = "efcentering.silo";
 
245
    int             layer, zone;
 
246
    int             nodelist2[9*4] = {0,1,5,4,    1,2,6,5,    2,3,7,6,
 
247
                                      4,5,9,8,    5,6,10,9,   6,7,11,10,
 
248
                                      8,9,13,12,  9,10,14,13, 10,11,15,14};
 
249
    int st2 = DB_ZONETYPE_QUAD;
 
250
    int ss2 = 4;
 
251
    int sc2 = 9;
 
252
    int nodelist3[27*8];
 
253
    int st3 = DB_ZONETYPE_HEX;
 
254
    int ss3 = 8;
 
255
    int sc3 = 27;
 
256
 
 
257
    int nedges;
 
258
    int *edges;
 
259
    int nfaces;
 
260
    int *faces;
 
261
    int ndims;
 
262
    int show_all_errors = FALSE;
 
263
 
 
264
    /* Parse command-line */
 
265
    for (i=1; i<argc; i++) {
 
266
        if (!strncmp(argv[i], "DB_PDB",6)) {
 
267
            driver = StringToDriver(argv[i]);
 
268
        } else if (!strncmp(argv[i], "DB_HDF5", 7)) {
 
269
            driver = StringToDriver(argv[i]);
 
270
        } else if (!strcmp(argv[i], "show-all-errors")) {
 
271
            show_all_errors = 1;
 
272
        } else if (argv[i][0] != '\0') {
 
273
            fprintf(stderr, "%s: ignored argument `%s'\n", argv[0], argv[i]);
 
274
        }
 
275
    }
 
276
    
 
277
    DBShowErrors(show_all_errors?DB_ALL_AND_DRVR:DB_ABORT, NULL);
 
278
    dbfile = DBCreate(filename, DB_CLOBBER, DB_LOCAL, "edge and face centered data", driver);
 
279
 
 
280
    coordnames[0] = "xcoords";
 
281
    coordnames[1] = "ycoords";
 
282
    coordnames[2] = "zcoords";
 
283
 
 
284
    dims[0] = 4;
 
285
    dims[1] = 4;
 
286
    dims[2] = 4;
 
287
    for (k = 0; k < 4; k++)
 
288
    {
 
289
        for (j = 0; j < 4; j++)
 
290
        {
 
291
            for (i = 0; i < 4; i++)
 
292
            {
 
293
                x[k*4*4+j*4+i] = (float) i;
 
294
                y[k*4*4+j*4+i] = (float) j;
 
295
                z[k*4*4+j*4+i] = (float) k;
 
296
                evar2d[0*16+j*4+i] = (float) i;
 
297
                evar2d[1*16+j*4+i] = (float) j;
 
298
                evar3d[0*64+k*4*4+j*4+i] = (float) i;
 
299
                evar3d[1*64+k*4*4+j*4+i] = (float) j;
 
300
                evar3d[2*64+k*4*4+j*4+i] = (float) k;
 
301
                fvar3d[0*64+k*4*4+j*4+i] = (float) 10*i;
 
302
                fvar3d[1*64+k*4*4+j*4+i] = (float) 100*j;
 
303
                fvar3d[2*64+k*4*4+j*4+i] = (float) 1000*k;
 
304
            }
 
305
        }
 
306
    }
 
307
 
 
308
    coords[0] = x;
 
309
    coords[1] = y;
 
310
    coords[2] = z;
 
311
 
 
312
    /* build 3d zonelist by layering 2d zonelist */
 
313
    for (layer = 0; layer < 3; layer++)
 
314
    {
 
315
        for (zone = 0; zone < 9; zone++)
 
316
        {
 
317
            nodelist3[layer*9*8+zone*8+0] = nodelist2[zone*4+0]+(layer+1)*16;
 
318
            nodelist3[layer*9*8+zone*8+1] = nodelist2[zone*4+0]+layer*16;
 
319
            nodelist3[layer*9*8+zone*8+2] = nodelist2[zone*4+1]+layer*16;
 
320
            nodelist3[layer*9*8+zone*8+3] = nodelist2[zone*4+1]+(layer+1)*16;
 
321
            nodelist3[layer*9*8+zone*8+4] = nodelist2[zone*4+3]+(layer+1)*16;
 
322
            nodelist3[layer*9*8+zone*8+5] = nodelist2[zone*4+3]+layer*16;
 
323
            nodelist3[layer*9*8+zone*8+6] = nodelist2[zone*4+2]+layer*16;
 
324
            nodelist3[layer*9*8+zone*8+7] = nodelist2[zone*4+2]+(layer+1)*16;
 
325
        }
 
326
    }
 
327
 
 
328
    DBPutQuadmesh(dbfile, "qmesh2", coordnames, coords, dims, 2, DB_FLOAT, DB_NONCOLLINEAR, 0);
 
329
    DBPutQuadmesh(dbfile, "qmesh3", coordnames, coords, dims, 3, DB_FLOAT, DB_NONCOLLINEAR, 0);
 
330
    DBPutQuadvar1(dbfile, "qevar2", "qmesh2", evar2d, dims, 2, 0, 0, DB_FLOAT, DB_EDGECENT, 0);
 
331
    DBPutQuadvar1(dbfile, "qevar3", "qmesh3", evar3d, dims, 3, 0, 0, DB_FLOAT, DB_EDGECENT, 0);
 
332
    DBPutQuadvar1(dbfile, "qfvar3", "qmesh3", fvar3d, dims, 3, 0, 0, DB_FLOAT, DB_FACECENT, 0);
 
333
 
 
334
    DBPutUcdmesh(dbfile, "umesh2", 2, coordnames, coords, 16, 9,  "um2zl", 0, DB_FLOAT, 0);
 
335
    DBPutUcdmesh(dbfile, "umesh3", 3, coordnames, coords, 64, 27, "um3zl", 0, DB_FLOAT, 0);
 
336
    DBPutZonelist2(dbfile, "um2zl", 9, 2, nodelist2, ss2*sc2, 0, 0, 0, &st2, &ss2, &sc2, 1, 0);
 
337
    DBPutZonelist2(dbfile, "um3zl", 27, 3, nodelist3, ss3*sc3, 0, 0, 0, &st3, &ss3, &sc3, 1, 0);
 
338
 
 
339
    /* Only reason we build an edgelist is so we know the number of unique edges in the mesh */
 
340
    build_edgelist(27, 3, nodelist3, ss3*sc3, 0, 0, 0, &st3, &ss3, &sc3, 1, &nedges, &edges);
 
341
    for (i = 0; i < nedges; i++)
 
342
        evar3d[i] = i;
 
343
    DBPutUcdvar1(dbfile, "uevar3", "umesh3", evar3d, nedges, 0, 0, DB_FLOAT, DB_EDGECENT, 0);
 
344
    ndims = 2;
 
345
    dims[0] = nedges;
 
346
    dims[1] = 2;
 
347
    DBWrite(dbfile, "edges", edges, dims, ndims, DB_INT);
 
348
    free(edges);
 
349
 
 
350
    /* Only reason we build a facelist is so we know the number of unique faces in the mesh */
 
351
    build_facelist(27, 3, nodelist3, ss3*sc3, 0, 0, 0, &st3, &ss3, &sc3, 1, &nfaces, &faces);
 
352
    for (i = 0; i < nfaces; i++)
 
353
        fvar3d[i] = i;
 
354
    DBPutUcdvar1(dbfile, "ufvar3", "umesh3", fvar3d, nfaces, 0, 0, DB_FLOAT, DB_FACECENT, 0);
 
355
    dims[0] = nfaces;
 
356
    dims[1] = 4;
 
357
    DBWrite(dbfile, "faces", faces, dims, ndims, DB_INT);
 
358
    free(faces);
 
359
 
 
360
    DBClose(dbfile);
 
361
 
 
362
    CleanupDriverStuff();
 
363
    return (0);
 
364
}