← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/ptst/form_div_rz_tst.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/ptst/form_div_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 div(u) r dr dz
 
24
//
 
25
// check value on the unit square [0,1]^2
 
26
//
 
27
//           d u_r   u_r   d u_z
 
28
//  div(u) = ----- + --- + -----
 
29
//            d r     r     d z
 
30
//
 
31
// component 0 : we check the u_r component : u_z = 0
 
32
//               and u_r takes monomial values : 1, r, ...
 
33
// component 1 : u_r = 0 and u_z is monomial
 
34
//
 
35
#include "rheolef.h"
 
36
using namespace rheolef;
 
37
using namespace std;
 
38
 
 
39
typedef Float (*function_type)(const point&);
 
40
 
 
41
struct list_type {
 
42
        string        name;
 
43
        function_type function;
 
44
        Float  expect [2];
 
45
};
 
46
Float monom_1  (const point& x) { return 1; }
 
47
Float monom_r  (const point& x) { return x[0]; }
 
48
Float monom_z  (const point& x) { return x[1]; }
 
49
Float monom_r2 (const point& x) { return x[0]*x[0]; }
 
50
Float monom_rz (const point& x) { return x[0]*x[1]; }
 
51
Float monom_z2 (const point& x) { return x[1]*x[1]; }
 
52
 
 
53
// axisymmetric rz
 
54
 
 
55
list_type rz_monom_list [] = {
 
56
        //   monom                    component
 
57
        // id   fct             u_r test     u_z test
 
58
        {"1",   monom_1,        {1,             0}              },
 
59
        {"r",   monom_r,        {1,             0}              },
 
60
        {"z",   monom_z,        {0.5,           0.5}            },
 
61
        {"r2",  monom_r2,       {1.,            0}              },
 
62
        {"rz",  monom_rz,       {0.5,           Float(1)/3}     },
 
63
        {"z2",  monom_z2,       {Float(1)/3,    0.5}            }
 
64
};
 
65
const unsigned int rz_monom_max = sizeof(rz_monom_list)/sizeof(list_type);
 
66
 
 
67
void usage()
 
68
{
 
69
      derr << "form_div_tst: usage: form_div_tst"
 
70
           << " {-Igeodir}*"
 
71
           << " -|mesh[.geo]"
 
72
           << " [-component int]"
 
73
           << " [-monom string]"
 
74
           << endl
 
75
           << "example:\n"
 
76
           << "  form_div_tst square -monom rz -component 0\n";
 
77
      exit (1);
 
78
}
 
79
int main(int argc, char**argv)
 
80
{
 
81
    environment rheolef (argc, argv);
 
82
    //
 
83
    // load geometry and options
 
84
    //
 
85
    geo omega;  
 
86
    string approx1 = "P2";
 
87
    string approx2 = "P1";
 
88
    string monom_id = "";
 
89
    size_t i_comp = 0;
 
90
    bool mesh_done = false;
 
91
    function_type monom_function = 0;
 
92
    unsigned int monom_idx = 0;
 
93
    string sys_coord_name = "rz";
 
94
 
 
95
    if (argc <= 1) usage() ;
 
96
 
 
97
    for (int i = 1; i < argc; i++ ) {
 
98
 
 
99
      if (argv [i][0] == '-' && argv [i][1] == 'I')
 
100
          append_dir_to_rheo_path (argv[i]+2) ;
 
101
      else if (strcmp(argv[i], "-component") == 0)
 
102
          i_comp = atoi(argv[++i]);
 
103
      else if (strcmp(argv[i], "-monom") == 0) {
 
104
          monom_id = argv[++i];
 
105
          for (unsigned int i = 0; i < rz_monom_max; i++) {
 
106
            if (monom_id == rz_monom_list [i].name) {
 
107
              monom_function = rz_monom_list [i].function;
 
108
              monom_idx = i;
 
109
            }
 
110
          }
 
111
      } else if (strcmp(argv[i], "-") == 0) {
 
112
          // input geo on standard input
 
113
          if (mesh_done) usage() ;
 
114
          derr << " ! load geo on stdin" << endl ;
 
115
          din >> omega ;
 
116
          mesh_done = true ;
 
117
      } else {
 
118
          // input geo on file
 
119
          if (mesh_done) usage() ;
 
120
          derr << " ! load " << argv[i] << endl ;
 
121
          omega = geo(argv[i]);
 
122
          mesh_done = true ;
 
123
      }
 
124
    }
 
125
    if (!mesh_done) usage() ;
 
126
    omega.set_coordinate_system (sys_coord_name);
 
127
    if (!monom_function) {
 
128
        error_macro("invalid monom identifier: " << monom_id);
 
129
    }
 
130
    derr << "syscoord = " << omega.coordinate_system() << endl
 
131
         << "monom = " << monom_id << endl;
 
132
    space V0h(omega, approx1);
 
133
    space Vh (omega, approx1, "vector");
 
134
    space Qh(omega, approx2);
 
135
    field monom (Vh, Float(0));
 
136
    monom[i_comp] = interpolate (V0h, monom_function);
 
137
    field one(Qh, Float(1));
 
138
    form  b (Vh,Qh,"div");
 
139
    Float result = b(monom, one);
 
140
    string vec_monom_id;
 
141
    if (i_comp == 0) vec_monom_id = monom_id + ",0";
 
142
    else             vec_monom_id = "0," + monom_id;
 
143
    derr << setprecision(numeric_limits<Float>::digits10)
 
144
         << "b(omega," << approx1 << "[" << vec_monom_id << "]," << approx2 << "[1]) = " << double(result) << endl;
 
145
    Float expect = rz_monom_list[monom_idx].expect[i_comp];
 
146
    Float tol = 1e-10;
 
147
 
 
148
    if (fabs(result - expect) <= tol) {
 
149
        derr << "ok" << endl;
 
150
        return 0;
 
151
    } else {
 
152
        derr << "but it was expected " << expect << endl
 
153
             << "and error = " << fabs(result - expect) << endl
 
154
             << "and tol = " << tol << endl;
 
155
        return 1;
 
156
    }
 
157
}