← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/geo_element/quadrature_t.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/geo_element/quadrature_t.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/quadrature.h"
 
22
#include "rheolef/gauss_jacobi.h"
 
23
using namespace rheolef;
 
24
using namespace std;
 
25
 
 
26
/* ------------------------------------------
 
27
 * Gauss quadrature formulae on the triangle
 
28
 * order r : exact for all polynoms of P_r
 
29
 *
 
30
 * References for optimal formulae :
 
31
 *   1)  G. Dhatt and G. Touzot,
 
32
 *       Une presentation de la methode des elements finis,
 
33
 *       Maloine editeur, Paris
 
34
 *       Deuxieme edition,
 
35
 *       1984,
 
36
 *       page 297
 
37
 *   2)  M. Crouzeix et A. L. Mignot
 
38
 *       Exercices d'analyse des equations differentielles,
 
39
 *       Masson,
 
40
 *       1986,
 
41
 *       p. 61
 
42
 * PROBLEM: direct optimized formulae have hard-coded
 
43
 *   coefficients in double precision.
 
44
 * TODO: compute numeric value at the fly why full 
 
45
 *   precision.
 
46
 * ------------------------------------------ 
 
47
 */
 
48
 
 
49
void
 
50
quadrature_on_geo::init_triangle (quadrature_option_type opt)
 
51
{
 
52
  // -------------------------------------------------------------------------
 
53
  // special case : superconvergent patch recovery nodes & weights
 
54
  // -------------------------------------------------------------------------
 
55
  quadrature_option_type::family_type f = opt.get_family();
 
56
  if (f == quadrature_option_type::superconvergent) {
 
57
    switch (opt.get_order()) {
 
58
      case 0:
 
59
      case 1:
 
60
        f =  quadrature_option_type::gauss;     
 
61
        break;
 
62
      case 2:
 
63
        f =  quadrature_option_type::middle_edge;       
 
64
        break;
 
65
      default:
 
66
        error_macro ("unsupported superconvergent("<<opt.get_order()<<")");
 
67
    }
 
68
  }
 
69
  // -------------------------------------------------------------------------
 
70
  // special case : Middle Edge formulae
 
71
  // e.g. when using special option in riesz_representer
 
72
  // -------------------------------------------------------------------------
 
73
  if (f == quadrature_option_type::middle_edge) {
 
74
    switch (opt.get_order()) {
 
75
     case 0 :
 
76
     case 1 :
 
77
     case 2 :
 
78
            // middle edge formulae:
 
79
            wx(x(Float(0.5), Float(0.0)), 1./6);
 
80
            wx(x(Float(0.0), Float(0.5)), 1./6);
 
81
            wx(x(Float(0.5), Float(0.5)), 1./6);
 
82
            break;
 
83
     default:
 
84
            error_macro ("unsupported Middle-Edge("<<opt.get_order()<<")");
 
85
    }
 
86
    return;
 
87
  }
 
88
  // -------------------------------------------------------------------------
 
89
  // special case : Gauss-Lobatto points
 
90
  // e.g. when using special option in riesz() function
 
91
  // -------------------------------------------------------------------------
 
92
  if (f == quadrature_option_type::gauss_lobatto) {
 
93
    switch (opt.get_order()) {
 
94
     case 0 :
 
95
     case 1 :
 
96
            // trapezes:
 
97
            wx(x(Float(0), Float(0)), 1./6);
 
98
            wx(x(Float(1), Float(0)), 1./6);
 
99
            wx(x(Float(0), Float(1)), 1./6);
 
100
            break;
 
101
     case 2 :
 
102
            // simpson:
 
103
            wx(x(Float(0), Float(0)), 3/Float(120));
 
104
            wx(x(Float(1), Float(0)), 3/Float(120));
 
105
            wx(x(Float(0), Float(1)), 3/Float(120));
 
106
            wx(x(Float(0.5), Float(0.0)), 8/Float(120));
 
107
            wx(x(Float(0.0), Float(0.5)), 8/Float(120));
 
108
            wx(x(Float(0.5), Float(0.5)), 8/Float(120));
 
109
            wx(x(1/Float(3), 1/Float(3)), 27/Float(120));
 
110
            break;
 
111
     default:
 
112
            error_macro ("unsupported Gauss-Lobatto("<<opt.get_order()<<")");
 
113
    }
 
114
    return;
 
115
  }
 
116
  // -------------------------------------------------------------------------
 
117
  // default & general case : Gauss points
 
118
  // -------------------------------------------------------------------------
 
119
  assert_macro (f == quadrature_option_type::gauss,
 
120
        "unsupported quadrature family \"" << opt.get_family_name() << "\"");
 
121
 
 
122
  switch (opt.get_order()) {
 
123
   case 0 :
 
124
   case 1 :
 
125
            // better than the general case: 
 
126
            //  r=0 => 2 nodes
 
127
            //  r=1 => 4 nodes
 
128
            //  here : 1 node
 
129
            wx(x(1/Float(3), 1/Float(3)), 1./2);
 
130
            break;
 
131
   case 2 :
 
132
            // better than the general case: 
 
133
            //  r=2 => 6 nodes
 
134
            //  here : 3 node
 
135
            wx(x(1/Float(6), 1/Float(6)), 1/Float(6));
 
136
            wx(x(4/Float(6), 1/Float(6)), 1/Float(6));
 
137
            wx(x(1/Float(6), 4/Float(6)), 1/Float(6));
 
138
            break;
 
139
#if !(defined(_RHEOLEF_HAVE_CLN) || defined(_RHEOLEF_HAVE_DOUBLEDOUBLE) || defined(_RHEOLEF_HAVE_LONG_DOUBLE))
 
140
            // numerical values are inlined in double precision : 15 digits !
 
141
   case 3 :
 
142
   case 4 : {
 
143
            // better than the general case: 
 
144
            //  r=3 => 9 nodes
 
145
            //  r=4 => 12 nodes
 
146
            //  here : 6 node
 
147
            Float a  = 0.445948490915965;
 
148
            Float b  = 0.091576213509771; 
 
149
            Float wa = 0.111690794839005;
 
150
            Float wb = 0.054975871827661;
 
151
            wx(x(a, a),     wa);
 
152
            wx(x(1-2*a, a), wa);
 
153
            wx(x(a, 1-2*a), wa);
 
154
            wx(x(b, b),     wb);
 
155
            wx(x(1-2*b, b), wb);
 
156
            wx(x(b, 1-2*b), wb);
 
157
            break;
 
158
          }
 
159
   case 5 :
 
160
   case 6 : {
 
161
            // better than the general case: 
 
162
            //  r=5 => 16 nodes
 
163
            //  r=6 => 20 nodes
 
164
            //  here : 12 node
 
165
            Float a  = 0.063089014491502;
 
166
            Float b  = 0.249286745170910;
 
167
            Float c  = 0.310352451033785;
 
168
            Float d  = 0.053145049844816;
 
169
            Float wa = 0.025422453185103;
 
170
            Float wb = 0.058393137863189;
 
171
            Float wc = 0.041425537809187;
 
172
            wx(x(a, a),         wa);
 
173
            wx(x(1-2*a, a),     wa);
 
174
            wx(x(a, 1-2*a),     wa);
 
175
            wx(x(b, b),         wb);
 
176
            wx(x(1-2*b, b),     wb);
 
177
            wx(x(b, 1-2*b),     wb);
 
178
            wx(x(c, d),         wc);
 
179
            wx(x(d, c),         wc);
 
180
            wx(x(1-(c+d), c),   wc);
 
181
            wx(x(1-(c+d), d),   wc);
 
182
            wx(x(c, 1-(c+d)),   wc);
 
183
            wx(x(d, 1-(c+d)),   wc);
 
184
            break;
 
185
          }
 
186
#endif // BIGFLOAT
 
187
   default: {
 
188
    // Gauss-Legendre quadrature formulae 
 
189
    //  where Legendre = Jacobi(alpha=0,beta=0) polynoms
 
190
    size_t r = opt.get_order();
 
191
    size_type n0 = n_node_gauss(r);
 
192
    size_type n1 = n_node_gauss(r+1);
 
193
    vector<Float> zeta0(n0), omega0(n0);
 
194
    vector<Float> zeta1(n1), omega1(n1);
 
195
    gauss_jacobi (n0, 0, 0, zeta0.begin(), omega0.begin());
 
196
    gauss_jacobi (n1, 0, 0, zeta1.begin(), omega1.begin());
 
197
 
 
198
    for (size_type i = 0; i < n0; i++) {
 
199
      for (size_type j = 0; j < n1; j++) {
 
200
        // we transform the square into the triangle 
 
201
        Float eta_0 = (1+zeta0[i])*(1-zeta1[j])/4;
 
202
        Float eta_1 =              (1+zeta1[j])/2;
 
203
        Float J     =              (1-zeta1[j])/8;
 
204
        wx (x(eta_0,eta_1), J*omega0[i]*omega1[j]);
 
205
      }
 
206
    }
 
207
   }
 
208
  }
 
209
}