← Back to branch summary

~ubuntu-branches/ubuntu/raring/rheolef/raring-proposed

« back to all changes in this revision

Viewing changes to nfem/bin/mkgeo_grid_3d.cc

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito
  • Date: 2012-04-06 09:12:21 UTC
  • mfrom: (1.1.5)
  • Revision ID: package-import@ubuntu.com-20120406091221-m58me99p1nxqui49
Tags: 6.0-1
http://bugs.debian.org/667356
http://bugs.debian.org/661689
* New upstream release 6.0 (major changes):
  - massively distributed and parallel support
  - full FEM characteristic method (Lagrange-Gakerkin method) support
  - enhanced users documentation 
  - source code supports g++-4.7 (closes: #667356)
* debian/control: dependencies for MPI distributed solvers added
* debian/rules: build commands simplified
* debian/librheolef-dev.install: man1/* to man9/* added
* debian/changelog: package description rewritted (closes: #661689)

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/bin/mkgeo_grid_3d.cc
1
 
///
2
 
/// This file is part of Rheolef.
3
 
///
4
 
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
5
 
///
6
 
/// Rheolef is free software; you can redistribute it and/or modify
7
 
/// it under the terms of the GNU General Public License as published by
8
 
/// the Free Software Foundation; either version 2 of the License, or
9
 
/// (at your option) any later version.
10
 
///
11
 
/// Rheolef is distributed in the hope that it will be useful,
12
 
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
13
 
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14
 
/// GNU General Public License for more details.
15
 
///
16
 
/// You should have received a copy of the GNU General Public License
17
 
/// along with Rheolef; if not, write to the Free Software
18
 
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
19
 
///
20
 
/// =========================================================================
21
 
# include <rheolef/compiler.h>
22
 
using namespace rheolef;
23
 
using namespace std;
24
 
// 
25
 
// generate a mesh for a cube
26
 
// ----------------------------
27
 
# define TETRA 'T'
28
 
# define PENTA 'P'
29
 
# define HEXA  'H'
30
 
 
31
 
Float A = 0;
32
 
Float B = 1;
33
 
Float C = 0;
34
 
Float D = 1;
35
 
Float F = 0;
36
 
Float G = 1;
37
 
 
38
 
# define no(i,j,k) ((i)+(nx+1)*((j)+(ny+1)*(k)))
39
 
 
40
 
void usage (const char* prog)
41
 
{
42
 
        cerr << "uniform grid nx*ny*nz for the [a,b]x[c,d]x[e,f] parallelotope" << endl
43
 
             << "usage: " << endl
44
 
             << prog << " "
45
 
             << "[-T|-P|-H]"
46
 
             << "[nx=10][ny=nx][nz=ny] "
47
 
             << "[-a float="<< A <<"] "
48
 
             << "[-b float="<< B <<"] "
49
 
             << "[-c float="<< C <<"] "
50
 
             << "[-d float="<< D <<"] "
51
 
             << "[-f float="<< F <<"] "
52
 
             << "[-g float="<< G <<"] "
53
 
             << endl
54
 
             << "example: " << endl
55
 
             << prog << " 10 " << endl
56
 
             << prog << " 10 20 " << endl
57
 
             << prog << " 10 20 15" << endl
58
 
             ;
59
 
        exit (1);
60
 
}
61
 
void print_tetra (int p0, int p1, int p2, int p3)
62
 
{
63
 
      cout << "T\t"
64
 
           << p0 << " " << p1 << " " 
65
 
           << p2 << " " << p3 << endl;
66
 
}
67
 
void print_penta (int type, int p0, int p1, int p2, int p3, int p4, int p5)
68
 
{
69
 
    if (type == PENTA)
70
 
      cout << "P\t"
71
 
           << p0 << " " << p1 << " " 
72
 
           << p2 << " " << p3 << " " 
73
 
           << p4 << " " << p5 << endl;
74
 
    else {
75
 
        print_tetra (p0, p1, p2, p3);
76
 
        print_tetra (p4, p1, p3, p5);
77
 
        print_tetra (p5, p1, p3, p2);
78
 
    }
79
 
}
80
 
void print_triangle (int p0, int p1, int p2) {
81
 
      cout << "t\t" << p0 << " " << p1 << " " << p2 << endl;
82
 
}
83
 
void print_rectangle (int p0, int p1, int p2, int p3) {
84
 
      cout << "q\t" << p0 << " " << p1 << " " << p2 << " " << p3 << endl;
85
 
}
86
 
void print_hexa (int type, int p0, int p1, int p2, int p3, int p4, int p5, int p6, int p7)
87
 
{
88
 
    if (type == HEXA)
89
 
      cout << "H\t"
90
 
           << p0 << " " << p1 << " " 
91
 
           << p2 << " " << p3 << " " 
92
 
           << p4 << " " << p5 << " " 
93
 
           << p6 << " " << p7 << endl;
94
 
    else if (type == PENTA) {
95
 
            print_penta (type, p0, p1, p2, p4, p5, p6);
96
 
            print_penta (type, p7, p6, p4, p3, p2, p0);
97
 
    } else {
98
 
            print_tetra (p1, p2, p0, p5);
99
 
            print_tetra (p6, p5, p4, p2);
100
 
            print_tetra (p4, p2, p5, p0);
101
 
            
102
 
            print_tetra (p7, p6, p4, p3);
103
 
            print_tetra (p0, p4, p2, p3);
104
 
            print_tetra (p6, p2, p4, p3);
105
 
    }
106
 
}
107
 
void print_domain (const char* name, int nface) {
108
 
    cout << "domain"        << endl
109
 
         << name            << endl
110
 
         << "1 2 " << nface << endl; 
111
 
}
112
 
int
113
 
main (int argc, char**argv)
114
 
{
115
 
    char* prog = argv[0];
116
 
    bool boundary = false;
117
 
    bool region = false;
118
 
    bool corner = false;
119
 
    if (argc == 1) usage(prog);
120
 
 
121
 
    // machine-dependent precision
122
 
    int digits10 = numeric_limits<Float>::digits10;
123
 
    cout << setprecision(digits10);
124
 
 
125
 
    // default triangle
126
 
    int type = TETRA;
127
 
    
128
 
    // default nb edge on x axis
129
 
    int nx = 10;
130
 
    int has_reset_nx = 0;
131
 
 
132
 
    // default nb edge on y axis
133
 
    int ny = 10;
134
 
    int has_reset_ny = 0;
135
 
 
136
 
    // default nb edge on z axis
137
 
    int nz = 10;
138
 
 
139
 
    int i,j,k;
140
 
 
141
 
    // parse command line 
142
 
    for (i = 1; i < argc; i++) {
143
 
 
144
 
        switch (argv[i][0]) {
145
 
 
146
 
          case '-' : {
147
 
            if (strcmp (argv[i], "-boundary") == 0)    { boundary = true; } 
148
 
            else if (strcmp (argv[i], "-region") == 0) { region = true; } 
149
 
            else if (strcmp (argv[i], "-corner") == 0) { corner = true; } 
150
 
            else switch (argv[i][1]) {
151
 
                case TETRA : 
152
 
                case PENTA :
153
 
                case HEXA  : {
154
 
                    type = argv[i][1];
155
 
                    break;
156
 
                }
157
 
                default : {
158
 
                    int j = i+1;
159
 
                    if (j >= argc) usage(prog);
160
 
                    switch (argv[i][1]) {
161
 
                    case 'a'   : A = atof(argv[j]);
162
 
                             break;
163
 
                    case 'b'   : B = atof(argv[j]);
164
 
                             break;
165
 
                    case 'c'   : C = atof(argv[j]);
166
 
                             break;
167
 
                    case 'd'   : D = atof(argv[j]);
168
 
                             break;
169
 
                    case 'f'   : F = atof(argv[j]);
170
 
                             break;
171
 
                    case 'g'   : G = atof(argv[j]);
172
 
                             break;
173
 
                    default    : usage(prog);
174
 
                             break;
175
 
                    }
176
 
                    ++i;
177
 
                    break;
178
 
                }
179
 
            }
180
 
            break;
181
 
          }
182
 
          default : {
183
 
            if (has_reset_ny) {
184
 
                nz = atoi(argv[i]);
185
 
            } else if (has_reset_nx) {
186
 
                nz = ny = atoi(argv[i]);
187
 
                has_reset_ny = 1;
188
 
            } else {
189
 
                nz = ny = nx = atoi(argv[i]);
190
 
                has_reset_nx = 1;
191
 
            }
192
 
            break;
193
 
          }
194
 
        }
195
 
    }
196
 
    if (A >= B || C >= D || F >= G) {
197
 
        cerr << prog << ": [a,b]x[c,d]x[f,g] may be non-empty.\n";
198
 
        exit (1);
199
 
    }
200
 
    if (region && nx % 2 != 0) {
201
 
        cerr << prog << ": region: nx may be an even number.\n";
202
 
        exit (1);
203
 
    }
204
 
    int np = (nx+1)*(ny+1)*(nz+1);
205
 
    int ne;
206
 
    if (type == TETRA)
207
 
        ne = 6*nx*ny*nz;
208
 
    else if (type == PENTA)
209
 
        ne = 2*nx*ny*nz;
210
 
    else
211
 
        ne = nx*ny*nz;
212
 
    
213
 
    // header
214
 
    cout << "mesh" << endl
215
 
         << "1 3 " << np << " " << ne << endl
216
 
         << endl;
217
 
 
218
 
 
219
 
    // geometry
220
 
    for (k = 0; k <= nz; k++)
221
 
      for (j = 0; j <= ny; j++)
222
 
        for (i = 0; i <= nx; i++)
223
 
           cout << A+(B-A)*i/Float(nx) << " "
224
 
                << C+(D-C)*j/Float(ny) << " "
225
 
                << F+(G-F)*k/Float(nz) << endl;
226
 
    // connectivity
227
 
    for (i = 0; i < nx; i++) 
228
 
      for (j = 0; j < ny; j++) 
229
 
        for (k = 0; k < nz; k++) 
230
 
              print_hexa (type, 
231
 
                no(i,j,k),   no(i+1,j,k),   no(i+1,j+1,k),   no(i,j+1,k),
232
 
                no(i,j,k+1), no(i+1,j,k+1), no(i+1,j+1,k+1), no(i,j+1,k+1));
233
 
    cout << endl;
234
 
 
235
 
    // domains = 6 sides
236
 
    if (boundary) {
237
 
        if (type == HEXA)  print_domain("boundary",  6*nx*ny);
238
 
        if (type == PENTA) print_domain("boundary",  8*nx*ny);
239
 
        if (type == TETRA) print_domain("boundary", 12*nx*ny);
240
 
    } else {
241
 
        print_domain("bottom", (type != HEXA) ? 2*nx*ny : nx*ny);
242
 
    }
243
 
    for (i = 0; i < nx; i++) {
244
 
        for (j = 0; j < ny; j++) {
245
 
            if (type != HEXA) {
246
 
                print_triangle (no(i,j,0), no(i+1,j+1,0), no(i+1,j,0));
247
 
                print_triangle (no(i,j,0), no(i,j+1,0), no(i+1,j+1,0));
248
 
            } else {
249
 
                print_rectangle (no(i,j,0), no(i,j+1,0), no(i+1,j+1,0), no(i+1,j,0));
250
 
            }
251
 
        }
252
 
    }
253
 
    if (!boundary) {
254
 
        cout << endl;
255
 
        print_domain("top", (type != HEXA) ? 2*nx*ny : nx*ny);
256
 
    }
257
 
    for (i = 0; i < nx; i++) {
258
 
        for (j = 0; j < ny; j++) {
259
 
            if (type != HEXA) {
260
 
                print_triangle(no(i,j,nz), no(i+1,j,nz), no(i+1,j+1,nz));
261
 
                print_triangle(no(i,j,nz), no(i+1,j+1,nz), no(i,j+1,nz));
262
 
            } else {
263
 
                print_rectangle (no(i,j,nz), no(i+1,j,nz), no(i+1,j+1,nz), no(i,j+1,nz));
264
 
            }
265
 
        }
266
 
    }
267
 
    if (!boundary) {
268
 
        cout << endl;
269
 
        print_domain("left", (type == TETRA) ? 2*nx*ny : nx*ny);
270
 
    }
271
 
    for (i = 0; i < nx; i++) {
272
 
        for (k = 0; k < nz; k++) {
273
 
            if (type == TETRA) {
274
 
                print_triangle(no(i,0,k), no(i+1,0,k), no(i+1,0,k+1));
275
 
                print_triangle(no(i,0,k), no(i+1,0,k+1), no(i,0,k+1));
276
 
            } else {
277
 
                print_rectangle (no(i,0,k), no(i+1,0,k), no(i+1,0,k+1), no(i,0,k+1));
278
 
            }
279
 
        }
280
 
    }
281
 
    if (!boundary) {
282
 
        cout << endl;
283
 
        print_domain("front", (type == TETRA) ? 2*nx*ny : nx*ny);
284
 
    }
285
 
    for (j = 0; j < ny; j++) {
286
 
        for (k = 0; k < nz; k++)  {
287
 
            if (type == TETRA) {
288
 
                print_triangle(no(nx,j,k), no(nx,j+1,k), no(nx,j,k+1)); 
289
 
                print_triangle(no(nx,j+1,k), no(nx,j+1,k+1), no(nx,j,k+1));
290
 
            } else {
291
 
                print_rectangle (no(nx,j,k), no(nx,j+1,k), no(nx,j+1,k+1), no(nx,j,k+1));
292
 
            }
293
 
        }
294
 
    }
295
 
    if (!boundary) {
296
 
        cout << endl;
297
 
        print_domain("right", (type == TETRA) ? 2*nx*ny : nx*ny);
298
 
    }
299
 
    for (i = 0; i < nx; i++) {
300
 
        for (k = 0; k < nz; k++) {
301
 
            if (type == TETRA) {
302
 
                print_triangle(no(i,ny,k), no(i+1,ny,k+1), no(i+1,ny,k)); 
303
 
                print_triangle(no(i,ny,k), no(i,ny,k+1), no(i+1,ny,k+1));
304
 
            } else {
305
 
                print_rectangle (no(i,ny,k), no(i,ny,k+1), no(i+1,ny,k+1), no(i+1,ny,k));
306
 
            }
307
 
        }
308
 
    }
309
 
    if (!boundary) {
310
 
        cout << endl;
311
 
        print_domain("back", (type == TETRA) ? 2*nx*ny : nx*ny);
312
 
    }
313
 
    for (j = 0; j < ny; j++) {
314
 
        for (k = 0; k < nz; k++)  {
315
 
            if (type == TETRA) {
316
 
                print_triangle(no(0,j,k), no(0,j,k+1), no(0,j+1,k));
317
 
                print_triangle(no(0,j+1,k), no(0,j,k+1), no(0,j+1,k+1));
318
 
            } else {
319
 
                print_rectangle (no(0,j,k), no(0,j,k+1), no(0,j+1,k+1), no(0,j+1,k));
320
 
            }
321
 
        }
322
 
    }
323
 
    cout << endl;
324
 
 
325
 
    if (region) {
326
 
            int nvol;
327
 
            if (type == TETRA)
328
 
                nvol = 6*nx*ny*nz/2;
329
 
            else if (type == PENTA)
330
 
                nvol = 2*nx*ny*nz/2;
331
 
            else
332
 
                nvol = nx*ny*nz/2;
333
 
        
334
 
            cout << "\ndomain\nwest\n1 3 " << nvol << "\n" ;
335
 
            for (i = 0; i < nx/2; i++) 
336
 
              for (j = 0; j < ny; j++) 
337
 
                for (k = 0; k < nz; k++) 
338
 
                      print_hexa (type, 
339
 
                        no(i,j,k),   no(i+1,j,k),   no(i+1,j+1,k),   no(i,j+1,k),
340
 
                        no(i,j,k+1), no(i+1,j,k+1), no(i+1,j+1,k+1), no(i,j+1,k+1));
341
 
            cout << endl;
342
 
        
343
 
            cout << "\ndomain\neast\n1 3 " << nvol << "\n" ;
344
 
            for (i = nx/2; i < nx; i++) 
345
 
              for (j = 0; j < ny; j++) 
346
 
                for (k = 0; k < nz; k++) 
347
 
                      print_hexa (type, 
348
 
                        no(i,j,k),   no(i+1,j,k),   no(i+1,j+1,k),   no(i,j+1,k),
349
 
                        no(i,j,k+1), no(i+1,j,k+1), no(i+1,j+1,k+1), no(i,j+1,k+1));
350
 
            cout << endl;
351
 
    }
352
 
    // corners = 4 0d-domains
353
 
    if (corner) {
354
 
            cout << "\ndomain\nleft_bottom_back\n1 0 1\n"   << no( 0, 0, 0) << endl
355
 
                 << "\ndomain\nright_bottom_back\n1 0 1\n"  << no(nx, 0, 0) << endl
356
 
                 << "\ndomain\nleft_top_back\n1 0 1\n"      << no(0, ny, 0) << endl
357
 
                 << "\ndomain\nright_top_back\n1 0 1\n"     << no(nx,ny, 0) << endl
358
 
                 << "\ndomain\nleft_bottom_front\n1 0 1\n"  << no( 0, 0,nz) << endl
359
 
                 << "\ndomain\nright_bottom_front\n1 0 1\n" << no(nx, 0,nz) << endl
360
 
                 << "\ndomain\nleft_top_front\n1 0 1\n"     << no(0, ny,nz) << endl
361
 
                 << "\ndomain\nright_top_front\n1 0 1\n"    << no(nx,ny,nz) << endl
362
 
                 << endl;
363
 
    }
364
 
    return 0;
365
 
}
366