← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty

« back to all changes in this revision

Viewing changes to nfem/bin/mkgeo_grid_2d.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_2d.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 square
26
 
// ----------------------------
27
 
# define TRI 3
28
 
# define QUA 4
29
 
 
30
 
Float A = 0;
31
 
Float B = 1;
32
 
Float C = 0;
33
 
Float D = 1;
34
 
int nx = 10;
35
 
 
36
 
# define no(i,j) ((i)*(ny+1)+(j))
37
 
 
38
 
void usage(const char *prog)
39
 
{
40
 
  cerr << "uniform grid nx*ny for the [a,b]x[c,d] rectangle" << endl
41
 
       << "usage: " << endl
42
 
       << prog << " "
43
 
       << "[-t|-q] [nx="<<nx<<"] [ny=nx] "
44
 
       << "[-a float=" << A << "] [-b float=" << B << "] [-c float=" << C << "] [-d float=" << D << "] "
45
 
       << "[-boundary] " << endl
46
 
       << "[-region] " << endl
47
 
       << "[-rz] " << endl
48
 
       << "[-corner] " << endl
49
 
       << "example:" << endl
50
 
       << prog << " -t 10" << endl
51
 
       << prog << " -t 10 20" << endl
52
 
      ;
53
 
  exit (1);
54
 
}
55
 
int main (int argc, char**argv)
56
 
{
57
 
    char* prog = argv[0];
58
 
    bool boundary = false;
59
 
    bool region   = false;
60
 
    bool corner   = false;
61
 
    bool axisymmetric   = false;
62
 
    if (argc == 1) usage(prog);
63
 
 
64
 
    // machine-dependent precision
65
 
    int digits10 = numeric_limits<Float>::digits10;
66
 
    cout << setprecision(digits10);
67
 
 
68
 
    // default triangle
69
 
    int type = TRI;
70
 
    
71
 
    // default nb edge on x axis
72
 
    int has_reset_nx = 0;
73
 
 
74
 
    // default nb edge on y axis
75
 
    int ny = 10;
76
 
 
77
 
    int i,j ;
78
 
 
79
 
    // parse command line 
80
 
    for (i = 1; i < argc; i++) {
81
 
 
82
 
        switch (argv[i][0]) {
83
 
 
84
 
        case '-' : {
85
 
            if (strcmp (argv[i], "-boundary") == 0) { boundary = true; }
86
 
            else if (strcmp (argv[i], "-region") == 0)   { region = true; }
87
 
            else if (strcmp (argv[i], "-rz") == 0)   { axisymmetric = true; }
88
 
            else if (strcmp (argv[i], "-corner") == 0)   { corner = true; }
89
 
            else switch (argv[i][1]) {
90
 
 
91
 
                case 't' : type = TRI; break;
92
 
                case 'q' : type = QUA; break;
93
 
 
94
 
                case 'a' : A = atof (argv[++i]); break;
95
 
                case 'b' : B = atof (argv[++i]); break;
96
 
                case 'c' : C = atof (argv[++i]); break;
97
 
                case 'd' : D = atof (argv[++i]); break;
98
 
 
99
 
                default :  usage(prog);
100
 
            }
101
 
            break;
102
 
        }
103
 
        default :
104
 
            if (has_reset_nx) {
105
 
                ny = atoi(argv[i]);
106
 
            } else {
107
 
                nx = ny = atoi(argv[i]);
108
 
                has_reset_nx = 1;
109
 
            }
110
 
            break;
111
 
        }
112
 
    }
113
 
    if (A >= B || C >= D) {
114
 
        cerr << prog << ": [a,b]x[c,d] may be non-empty.\n";
115
 
        exit (1);
116
 
    }
117
 
    if (region && nx % 2 != 0) {
118
 
        cerr << prog << ": region: nx may be an even number.\n";
119
 
        exit (1);
120
 
    }
121
 
    if (type == TRI) cerr << "triangular" ; else cerr << "quadrangular" ;
122
 
    cerr << " mesh " << nx << " x " << ny << "\n" ; 
123
 
    
124
 
    int n_vert = (nx+1)*(ny+1);
125
 
    int n_elt;
126
 
    int n_edg;
127
 
    if (type == TRI) {
128
 
        n_elt = 2*nx*ny;
129
 
        n_edg = 3*nx*ny+nx+ny; // = (nx+1)*ny + (ny+1)*nx + nx*ny
130
 
    } else {
131
 
        n_elt = nx*ny;
132
 
        n_edg = 2*nx*ny+nx+ny; // = (nx+1)*ny + (ny+1)*nx
133
 
    }
134
 
    // header
135
 
    size_t version = (axisymmetric ? 3 : 2);
136
 
    cout << "#!geo\n\nmesh\n"
137
 
         << version << " 2 " << n_vert << " " << n_elt << " " << n_edg << "\n";
138
 
    if (axisymmetric)
139
 
        cout << "coordinate_system rz" << endl;
140
 
    cout << endl;
141
 
    
142
 
    // geometry
143
 
    for( i= 0 ; i <= nx ; i++ )
144
 
      for( j= 0 ; j <= ny ; j++ )
145
 
        cout << (B-A)*i/Float(nx)+A << " " << (D-C)*j/Float(ny)+C << "\n" ;
146
 
    cout << "\n" ;
147
 
 
148
 
    // connectivity
149
 
    //    elements
150
 
    for (i = 0; i < nx; i++) {
151
 
      for (j = 0; j < ny; j++) {
152
 
        if (type == TRI)
153
 
          cout << "t\t"   << no(i,j) << " " << no(i+1,j)   << " " << no(i+1,j+1) 
154
 
               << "\nt\t" << no(i,j) << " " << no(i+1,j+1) << " " << no(i,j+1) << "\n" ;
155
 
        else
156
 
          cout << "q\t" << no(i,j)     << " " << no(i+1,j) << " " 
157
 
               <<          no(i+1,j+1) << " " << no(i,j+1) << "\n" ;
158
 
        }
159
 
    }
160
 
    cout << "\n" ;
161
 
    //   edges
162
 
    //      vertical
163
 
    for( int i_ve=0 ; i_ve < nx+1 ; i_ve++ )
164
 
      for( int j_ve=0 ; j_ve < ny ; j_ve++ )
165
 
        cout << no(i_ve, j_ve) << " " << no(i_ve, j_ve+1) << "\n" ;
166
 
    //      horizontal
167
 
    for( int i_he=0 ; i_he < nx ; i_he++ )
168
 
      for( int j_he=0 ; j_he < ny+1 ; j_he++ )
169
 
        cout << no(i_he, j_he) << " " << no(i_he+1, j_he) << "\n" ;
170
 
    //      diagonal (TRI only)
171
 
    if (type==TRI)
172
 
      for ( int i_elt = 0 ; i_elt < nx ; i_elt++ )
173
 
        for ( int j_elt = 0 ; j_elt < ny ; j_elt++ )
174
 
          cout << no(i_elt, j_elt) <<  " " << no(i_elt+1, j_elt+1) << "\n" ;
175
 
    cout << "\n" ;
176
 
 
177
 
    // domains = 4 sides
178
 
    
179
 
    if (boundary) {
180
 
        cout << "domain\nboundary\n1 1 " << 2*(nx+ny) << "\n" ;
181
 
    } else {
182
 
        cout << "domain\nbottom\n1 1 " << nx << "\n" ;
183
 
    }
184
 
    for (i = 0; i < nx; i++)
185
 
      cout << no(i,0) << " " << no(i+1,0) << "\n" ;
186
 
 
187
 
    if (!boundary) cout << "\ndomain\nright\n1 1 " << ny << "\n" ;
188
 
    for (j = 0; j < ny; j++)
189
 
      cout << no(nx,j) << " " << no(nx,j+1) << "\n" ;
190
 
 
191
 
    if (!boundary) cout << "\ndomain\ntop\n1 1 " << nx << "\n" ;
192
 
    for (i = 0; i < nx; i++)
193
 
      cout << no(i+1,ny) << " " << no(i,ny) << "\n" ;
194
 
 
195
 
    if (!boundary) cout << "\ndomain\nleft\n1 1 " << ny << "\n" ;
196
 
    for (j = 0; j < ny; j++)
197
 
      cout << no(0,j+1) << " " << no(0,j) << "\n" ;
198
 
    cout << "\n" ;
199
 
 
200
 
    if (region) {
201
 
            int n_face = (type == TRI) ? nx*ny : nx*ny/2;
202
 
            cout << "\ndomain\nwest\n1 2 " << n_face << "\n" ;
203
 
            for (i = 0; i < nx/2; i++) {
204
 
              for (j = 0; j < ny; j++) {
205
 
                if (type == TRI)
206
 
                  cout << "t\t"   << no(i,j) << " " << no(i+1,j)   << " " << no(i+1,j+1) 
207
 
                       << "\nt\t" << no(i,j) << " " << no(i+1,j+1) << " " << no(i,j+1) << "\n" ;
208
 
                else
209
 
                  cout << "q\t" << no(i,j)     << " " << no(i+1,j) << " " 
210
 
                       <<          no(i+1,j+1) << " " << no(i,j+1) << "\n" ;
211
 
                }
212
 
            }
213
 
            cout << "\n" ;
214
 
            cout << "\ndomain\neast\n1 2 " << n_face << "\n" ;
215
 
            for (i = nx/2; i < nx; i++) {
216
 
              for (j = 0; j < ny; j++) {
217
 
                if (type == TRI)
218
 
                  cout << "t\t"   << no(i,j) << " " << no(i+1,j)   << " " << no(i+1,j+1) 
219
 
                       << "\nt\t" << no(i,j) << " " << no(i+1,j+1) << " " << no(i,j+1) << "\n" ;
220
 
                else
221
 
                  cout << "q\t" << no(i,j)     << " " << no(i+1,j) << " " 
222
 
                       <<          no(i+1,j+1) << " " << no(i,j+1) << "\n" ;
223
 
                }
224
 
            }
225
 
    }
226
 
    // corners = 4 0d-domains
227
 
    if (corner) {
228
 
            cout << "\ndomain\nleft_bottom\n1 0 1\n"  << no( 0, 0) << endl
229
 
                 << "\ndomain\nright_bottom\n1 0 1\n" << no(nx, 0) << endl
230
 
                 << "\ndomain\nleft_top\n1 0 1\n"     << no(0, ny) << endl
231
 
                 << "\ndomain\nright_top\n1 0 1\n"    << no(nx,ny) << endl
232
 
                 << endl;
233
 
    }
234
 
    return 0;
235
 
}