← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/ptst/form_circle_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_circle_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 weight(x) dx = mes(\Omega) with weight
 
24
//
 
25
// check value on the unit ball C(0,1) in IR^d
 
26
//
 
27
#include "rheolef/rheolef.h"
 
28
using namespace rheolef;
 
29
using namespace std;
 
30
 
 
31
/*
 
32
  I_{m,n}
 
33
  = int_Omega x^m y^n dx
 
34
  = (int_{r=0}^1 r^{m+n+1} dr)*(int_0^{2*pi} cos(theta)^m sin(theta)^n d theta)
 
35
  et :
 
36
   int_{r=0}^1 r^{m+n+1} dr = 1/(m+n+2)
 
37
  => I_{m,n}   = J_{m,n}/(m+n+2)
 
38
  J_{m,n} = int_0^{2*pi} cos(theta)^m sin(theta)^n d theta
 
39
        cf Bro-1985, p. 557 & 561 : recurrence
 
40
  J_{m,n} = 0 si m ou n impair
 
41
  J_{2p,2q} = (2q-1)/(2p+2q)*J_{2p,2q-2}
 
42
  J_{2p,0} = (2p-1)/(2p)*J_{2p-2,0}
 
43
  J_{0,0} = 2*pi
 
44
                  (2p)!
 
45
  => J_{2p,0} = ------------ * (2*pi)        pour q = 0
 
46
                (2^p * p!)^2
 
47
                 (2p)! (2q-1)! pi
 
48
  => J_{2p,2q} = --------------------------  pour q >=1
 
49
                 4^(p+q-1) p! (q-1)! (p+q)!
 
50
*/
 
51
Float factorial (size_t n)
 
52
{
 
53
  Float res = 1;
 
54
  for (size_t i = 1; i <= n; i++) res *= i;
 
55
  return res;
 
56
}
 
57
Float
 
58
integrate_2d_mass_xy (size_t nx, size_t ny, bool is_boundary)
 
59
{
 
60
  static Float pi = acos(-1.0);
 
61
  if ((nx % 2 == 1) || (ny % 2 == 1)) return 0;
 
62
  size_t p = nx/2;
 
63
  size_t q = ny/2;
 
64
  Float value;
 
65
  if (q == 0) {
 
66
    value = 2*pi*factorial(2*p)/sqr(pow(2.0,p)*factorial(p));
 
67
  } else {
 
68
    value = pi*factorial(2*p)*factorial(2*q-1)/(pow(4.0,p+q-1)*factorial(p)*factorial(q-1)*factorial(p+q));
 
69
  }
 
70
  if (!is_boundary) value = value/(nx+ny+2); // integrate also in r : interior of the circle
 
71
  return value;
 
72
}
 
73
// handle also grad_grad form:
 
74
Float
 
75
integrate_xy (size_t d, string name, size_t nx, size_t ny, size_t mx, size_t my, bool is_boundary)
 
76
{
 
77
  if (d == 2) { // circle
 
78
    if (name == "mass") return integrate_2d_mass_xy (nx+mx, ny+my, is_boundary);
 
79
    check_macro (name == "grad_grad", "unexpected form `"<<name<<"'");
 
80
    if (is_boundary) return 0;
 
81
    // here: form="grad_grad" on a filled circle
 
82
    Float value = 0;
 
83
    if (nx+mx >= 2) value += nx*mx*integrate_2d_mass_xy (nx+mx-2, ny+my,   is_boundary);
 
84
    if (ny+my >= 2) value += ny*my*integrate_2d_mass_xy (nx+mx,   ny+my-2, is_boundary);
 
85
    return value;
 
86
  } else { // sphere: only mass and nx=ny=0 yet
 
87
    check_macro (name == "mass", "form `"<<name<<"' not yet");
 
88
    check_macro (nx + ny + mx + my == 0, "non null nx etc: not yet");
 
89
    static Float pi = acos(-1.0);
 
90
    Float value = 4*pi;
 
91
    return value;
 
92
  }
 
93
}
 
94
void usage()
 
95
{
 
96
      derr << "usage: form_mass_circle_tst"
 
97
           << " -|mesh[.geo]"
 
98
           << " {-Igeodir}*"
 
99
           << " [-form name=mass]"
 
100
           << " -app approx"
 
101
           << " [-nx int]"
 
102
           << " [-ny int]"
 
103
           << " [-nz int]"
 
104
           << " [-mx int]"
 
105
           << " [-my int]"
 
106
           << " [-mz int]"
 
107
           << " [-tol float]"
 
108
           << endl
 
109
           << "example:" << endl
 
110
           << "  mkgeo_ball -order 3 -t 10 > ball.geo" << endl
 
111
           << "  form_mass_circle_tst ball -app P3 -nx 3" << endl;
 
112
      exit (1);
 
113
}
 
114
struct weight : std::unary_function<point,Float> {
 
115
  Float operator() (const point& x) const {
 
116
    if (nx == 0 && ny == 0) return 1;
 
117
    Float factx, facty;
 
118
    if (nx % 2 == 0) factx = pow(fabs(x[0]),nx);
 
119
    else             factx = x[0]*pow(fabs(x[0]),nx-1);
 
120
    if (ny % 2 == 0) facty = pow(fabs(x[1]),ny);
 
121
    else             facty = x[1]*pow(fabs(x[1]),ny-1);
 
122
    return factx*facty;
 
123
  }
 
124
  weight(size_t nx1=0, size_t ny1=0) : nx(nx1), ny(ny1) {}
 
125
  size_t nx, ny;
 
126
};
 
127
int main(int argc, char**argv)
 
128
{
 
129
    //
 
130
    // load geometry and options
 
131
    //
 
132
    environment rheolef (argc,argv);
 
133
    geo omega;  
 
134
    string form_name = "mass";
 
135
    string approx1 = "P1";
 
136
    string approx2 = "";
 
137
    size_t nx = 0;
 
138
    size_t ny = 0;
 
139
    size_t nz = 0;
 
140
    size_t mx = 0;
 
141
    size_t my = 0;
 
142
    size_t mz = 0;
 
143
    bool mesh_done = false;
 
144
    Float tol = 1e-10;
 
145
 
 
146
    if (argc <= 1) usage() ;
 
147
 
 
148
    for (int i = 1; i < argc; i++ ) {
 
149
 
 
150
      if (argv [i][0] == '-' && argv [i][1] == 'I') { append_dir_to_rheo_path (argv[i]+2) ; }
 
151
      else if (strcmp(argv[i], "-form") == 0) { form_name = argv[++i]; }
 
152
      else if (strcmp(argv[i], "-app") == 0) { approx1 = argv[++i]; }
 
153
      else if (strcmp(argv[i], "-tol") == 0) { tol = atof(argv[++i]); }
 
154
      else if (strcmp(argv[i], "-nx") == 0) { nx = atoi(argv[++i]); }
 
155
      else if (strcmp(argv[i], "-ny") == 0) { ny = atoi(argv[++i]); }
 
156
      else if (strcmp(argv[i], "-nz") == 0) { nz = atoi(argv[++i]); }
 
157
      else if (strcmp(argv[i], "-mx") == 0) { mx = atoi(argv[++i]); }
 
158
      else if (strcmp(argv[i], "-my") == 0) { my = atoi(argv[++i]); }
 
159
      else if (strcmp(argv[i], "-mz") == 0) { mz = atoi(argv[++i]); }
 
160
      else if (strcmp(argv[i], "-") == 0) {
 
161
          // input geo on standard input
 
162
          if (mesh_done) usage() ;
 
163
          derr << "! load geo on stdin" << endl ;
 
164
          din >> omega ;
 
165
          mesh_done = true ;
 
166
      } else {
 
167
          // input geo on file
 
168
          if (mesh_done) usage() ;
 
169
          //derr << "! load " << argv[i] << endl ;
 
170
          omega = geo(argv[i]);
 
171
          mesh_done = true ;
 
172
      }
 
173
    }
 
174
    warning_macro ("form = " << form_name);
 
175
    warning_macro ("nx = " << nx);
 
176
    warning_macro ("ny = " << ny);
 
177
    warning_macro ("mx = " << mx);
 
178
    warning_macro ("my = " << my);
 
179
    if (!mesh_done) usage() ;
 
180
    if (approx2 == "") {
 
181
        approx2 = approx1;
 
182
    }
 
183
    space V1(omega, approx1);
 
184
    space V2(omega, approx2);
 
185
    field w1h = interpolate (V1, weight(nx,ny));
 
186
    field w2h = interpolate (V2, weight(mx,my));
 
187
    form  a(V1,V2,form_name);
 
188
#ifdef TO_CLEAN
 
189
    odiststream out; out.open("a","mm"); out << a.uu;
 
190
    odiststream zzz; zzz.open("w1h","field"); zzz << w1h;
 
191
    form  mass(V1,V2,"mass");
 
192
    solver sm (mass.uu);
 
193
    field aw1h (V2);
 
194
    aw1h.u = sm.solve ((a*w1h).u);
 
195
    aw1h["boundary"] = -2;
 
196
    odiststream yyy; yyy.open("aw1h","field"); yyy << aw1h;
 
197
#endif // TO_CLEAN
 
198
    Float mes_omega = a(w1h, w2h);
 
199
    dout << setprecision(numeric_limits<Float>::digits10)
 
200
         << "mes(omega," << approx1 << "," << approx2 << ") = " << double(mes_omega) << endl;
 
201
    int status = 0;
 
202
    // TODO: ny+my
 
203
    bool is_boundary = (omega.dimension() > omega.map_dimension());
 
204
    Float expect = integrate_xy (omega.dimension(),form_name,nx, ny, mx, my, is_boundary);
 
205
    dout << "exact  = " << expect << endl
 
206
         << "error  = " << fabs(mes_omega - expect) << endl;
 
207
    return 0;
 
208
}