← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/tst/form_2D_D_rz_tst.cc

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2010-06-12 09:08:59 UTC
  • Revision ID: james.westby@ubuntu.com-20100612090859-8gpm2gc7j3ab43et
Tags: upstream-5.89
Import upstream version 5.89

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/tst/form_2D_D_rz_tst.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
//
 
22
// compute:
 
23
//  \int_\Omega 2 D(u):D(v) r dr dz
 
24
//
 
25
// check value on the unit square [0,1]^2
 
26
//
 
27
// component 0 : we check the u_r component : u_z = 0
 
28
//               and u_r takes monomial values : 1, r, ...
 
29
// component 1 : u_r = 0 and u_z is monomial
 
30
//
 
31
#include "rheolef/rheolef.h"
 
32
using namespace std;
 
33
 
 
34
typedef Float (*function_type)(const point&);
 
35
 
 
36
struct list_type {
 
37
        string        name;
 
38
        function_type function;
 
39
};
 
40
Float monom_1  (const point& x) { return 1; }
 
41
Float monom_r  (const point& x) { return x[0]; }
 
42
Float monom_z  (const point& x) { return x[1]; }
 
43
Float monom_r2 (const point& x) { return x[0]*x[0]; }
 
44
Float monom_rz (const point& x) { return x[0]*x[1]; }
 
45
Float monom_z2 (const point& x) { return x[1]*x[1]; }
 
46
 
 
47
// axisymmetric rz
 
48
 
 
49
list_type rz_monom_list [] = {
 
50
        //   monom
 
51
        // id   fct    
 
52
        {"1",   monom_1},
 
53
        {"r",   monom_r},
 
54
        {"z",   monom_z},
 
55
        {"r2",  monom_r2},
 
56
        {"rz",  monom_rz},
 
57
        {"z2",  monom_z2}
 
58
};
 
59
const unsigned int rz_monom_max = sizeof(rz_monom_list)/sizeof(list_type);
 
60
 
 
61
typedef Float T;
 
62
 
 
63
const Float infty = -1;
 
64
 
 
65
Float  expect_table [2][2] [rz_monom_max][rz_monom_max] = {
 
66
  {
 
67
    // a([ur,0],[vr,0])
 
68
    {//         1       r       z      r2      rz       z2
 
69
 
 
70
        {   infty,      2,  infty,      1,      1,  infty},
 
71
        {       2,      2,      1,      2,      1, T(2)/3},
 
72
        {   infty,      1,  infty,    0.5,      1,  infty},
 
73
 
 
74
        {       1,      2,    0.5,    2.5,      1,  T(1)/3},
 
75
        {       1,      1,      1,      1, T(11)/12, T(5)/6},
 
76
        {   infty, T(2)/3,  infty, T(1)/3, T(5)/6,  infty}
 
77
    },
 
78
    // a([ur,0],[0,vz])
 
79
    {//vz=      1       r       z      r2      rz       z2
 
80
        {       0,      0,      0,      0,      0,      0},
 
81
        {       0,      0,      0,      0,      0,      0},
 
82
        {       0,      0.5,    0,      T(2)/3, 0.25,   0},
 
83
        {       0,      0,      0,      0,      0,      0},
 
84
        {       0,      T(1)/3, 0,      0.5,    T(1)/6, 0},
 
85
        {       0,      0.5,    0,      T(2)/3, T(1)/3, 0}
 
86
    }
 
87
  },
 
88
  {
 
89
    // a([0,uz],[vr,0]) : not yet
 
90
    {
 
91
        {-1,    0,      0,      0,      0,      0},
 
92
        {0,     0,      0,      0,      0,      0},
 
93
        {0,     0,      0,      0,      0,      0},
 
94
        {0,     0,      0,      0,      0,      0},
 
95
        {0,     0,      0,      0,      0,      0},
 
96
        {0,     0,      0,      0,      0,      0}
 
97
    },
 
98
    // a([0,uz],[0,vz]) : not yet
 
99
    {
 
100
        {0,     0,      0,      0,      0,      0},
 
101
        {0,     0,      0,      0,      0,      0},
 
102
        {0,     0,      0,      0,      0,      0},
 
103
        {0,     0,      0,      0,      0,      0},
 
104
        {0,     0,      0,      0,      0,      0},
 
105
        {0,     0,      0,      0,      0,      0}
 
106
    }
 
107
  }
 
108
};
 
109
void usage()
 
110
{
 
111
      cerr << "form_2D_D_tst: usage: form_2D_D_tst"
 
112
           << " {-Igeodir}*"
 
113
           << " -|mesh[.geo]"
 
114
           << " [-approx string]"
 
115
           << " [-u-component int]"
 
116
           << " [-v-component int]"
 
117
           << " [-u-monom string]"
 
118
           << " [-v-monom string]"
 
119
           << endl;
 
120
      cerr << "example:\n";
 
121
      cerr << "  form_2D_D_tst square -u-monom rz -v-monom r -u-component 0 -v-component 0\n";
 
122
      exit (1);
 
123
}
 
124
int main(int argc, char**argv)
 
125
{
 
126
    //
 
127
    // load geometry and options
 
128
    //
 
129
    geo omega;  
 
130
    string coord_sys = "rz";
 
131
    string approx = "P2";
 
132
    bool mesh_done = false;
 
133
    size_t i_comp_u = 0;
 
134
    size_t i_comp_v = 0;
 
135
    string        u_monom_id = "";
 
136
    function_type u_monom_function = 0;
 
137
    unsigned int  u_monom_idx = 0;
 
138
    string        v_monom_id = "";
 
139
    function_type v_monom_function = 0;
 
140
    unsigned int  v_monom_idx = 0;
 
141
 
 
142
    if (argc <= 1) usage() ;
 
143
 
 
144
    for (int i = 1; i < argc; i++ ) {
 
145
 
 
146
      if (argv [i][0] == '-' && argv [i][1] == 'I')
 
147
          append_dir_to_rheo_path (argv[i]+2) ;
 
148
      else if (strcmp(argv[i], "-approx") == 0)
 
149
          approx = argv[++i];
 
150
      else if (strcmp(argv[i], "-u-component") == 0)
 
151
          i_comp_u = atoi(argv[++i]);
 
152
      else if (strcmp(argv[i], "-v-component") == 0)
 
153
          i_comp_v = atoi(argv[++i]);
 
154
      else if (strcmp(argv[i], "-u-monom") == 0) {
 
155
          u_monom_id = argv[++i];
 
156
          for (unsigned int i = 0; i < rz_monom_max; i++) {
 
157
            if (u_monom_id == rz_monom_list [i].name) {
 
158
              u_monom_function = rz_monom_list [i].function;
 
159
              u_monom_idx = i;
 
160
            }
 
161
          }
 
162
      } else if (strcmp(argv[i], "-v-monom") == 0) {
 
163
          v_monom_id = argv[++i];
 
164
          for (unsigned int i = 0; i < rz_monom_max; i++) {
 
165
            if (v_monom_id == rz_monom_list [i].name) {
 
166
              v_monom_function = rz_monom_list [i].function;
 
167
              v_monom_idx = i;
 
168
            }
 
169
          }
 
170
      } else if (strcmp(argv[i], "-") == 0) {
 
171
          // input geo on standard input
 
172
          if (mesh_done) usage() ;
 
173
          cerr << " ! load geo on stdin" << endl ;
 
174
          cin >> omega ;
 
175
          mesh_done = true ;
 
176
      } else {
 
177
          // input geo on file
 
178
          if (mesh_done) usage() ;
 
179
          cerr << " ! load " << argv[i] << endl ;
 
180
          omega = geo(argv[i]);
 
181
          mesh_done = true ;
 
182
      }
 
183
    }
 
184
    if (!mesh_done) {
 
185
        cerr << "form_2D_D_tst: no mesh specified" << endl;
 
186
        usage() ;
 
187
    }
 
188
    omega.set_coordinate_system (coord_sys);
 
189
    if (!u_monom_function) {
 
190
        error_macro("invalid u-monom identifier: " << u_monom_id);
 
191
    }
 
192
    if (!v_monom_function) {
 
193
        error_macro("invalid v-monom identifier: " << v_monom_id);
 
194
    }
 
195
    cerr << "syscoord = " << omega.coordinate_system() << endl;
 
196
    cerr << "u_monom = " << u_monom_id << endl;
 
197
    cerr << "v_monom = " << v_monom_id << endl;
 
198
    space V0h(omega, approx);
 
199
    space Vh = V0h*V0h;
 
200
    field u_monom (Vh, Float(0));
 
201
    u_monom[i_comp_u] = interpolate (V0h, u_monom_function);
 
202
    field v_monom (Vh, Float(0));
 
203
    v_monom[i_comp_v] = interpolate (V0h, v_monom_function);
 
204
    form  a (Vh,Vh,"2D_D");
 
205
    Float result = a(u_monom, v_monom);
 
206
    string u_vec_monom_id;
 
207
    string v_vec_monom_id;
 
208
    if (i_comp_u == 0) u_vec_monom_id = u_monom_id + ",0";
 
209
    else               u_vec_monom_id = "0," + u_monom_id;
 
210
    if (i_comp_v == 0) v_vec_monom_id = v_monom_id + ",0";
 
211
    else               v_vec_monom_id = "0," + v_monom_id;
 
212
    cerr << setprecision(numeric_limits<Float>::digits10)
 
213
         << "a(omega," << approx << "[" << u_vec_monom_id << "]," << approx << "[" << v_vec_monom_id << "]) = " << double(result) << endl;
 
214
    Float expect = expect_table [i_comp_u][i_comp_v][u_monom_idx][v_monom_idx];
 
215
    Float tol = 1e-10;
 
216
 
 
217
    if (expect == infty) {
 
218
        cerr << "ok (expect infty)" << endl;
 
219
        return 0;
 
220
    } else if (fabs(result - expect) <= tol) {
 
221
        cerr << "ok" << endl;
 
222
        return 0;
 
223
    } else {
 
224
        cerr << "but it was expected " << expect << endl;
 
225
        cerr << "and error = " << fabs(result - expect) << endl;
 
226
        cerr << "and tol = " << tol << endl;
 
227
        return 1;
 
228
    }
 
229
}