← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/geo_element/quadrature_e.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_e.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
/*
 
28
 * quadrature formulae
 
29
 * refs: G. Dhatt and G. Touzot,
 
30
 *       Une presentation de la methode des elements finis,
 
31
 *       Maloine editeur, Paris
 
32
 *       Deuxieme edition,
 
33
 *       1984,
 
34
 *       page 288
 
35
 */
 
36
 
 
37
void
 
38
quadrature_on_geo::init_edge (quadrature_option_type opt)
 
39
{
 
40
  // -------------------------------------------------------------------------
 
41
  // special case : superconvergent patch recovery nodes & weights
 
42
  // -------------------------------------------------------------------------
 
43
  quadrature_option_type::family_type f = opt.get_family();
 
44
  if (f == quadrature_option_type::superconvergent) {
 
45
      f =  quadrature_option_type::gauss;
 
46
  }
 
47
  // -------------------------------------------------------------------------
 
48
  // special case : Gauss-Lobatto points
 
49
  // e.g. when using special option in riesz_representer
 
50
  // -------------------------------------------------------------------------
 
51
  if (f == quadrature_option_type::gauss_lobatto) {
 
52
    switch (opt.get_order()) {
 
53
     case 0 :
 
54
     case 1 :
 
55
            // trapezes:
 
56
            wx(x(Float(0)), 0.5);
 
57
            wx(x(Float(1)), 0.5);
 
58
            break;
 
59
     case 2 :
 
60
     case 3 :
 
61
            // http://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules
 
62
            wx(x(Float(0)),   Float(1)/Float(6));
 
63
            wx(x(Float(0.5)), Float(4)/Float(6));
 
64
            wx(x(Float(1)),   Float(1)/Float(6));
 
65
            break;
 
66
     case 4 :
 
67
     case 5 :
 
68
            // http://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules
 
69
            wx(x(Float(0)),               Float(1)/Float(12));
 
70
            wx(x(0.5-0.5/sqrt(Float(5))), Float(5)/Float(12));
 
71
            wx(x(0.5+0.5/sqrt(Float(5))), Float(5)/Float(12));
 
72
            wx(x(Float(1)),               Float(1)/Float(12));
 
73
            break;
 
74
     case 6 :
 
75
     case 7 :
 
76
            // http://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules
 
77
            wx(x(Float(0)),                        Float( 9)/Float(180));
 
78
            wx(x(0.5-0.5/sqrt(Float(3)/Float(7))), Float(49)/Float(180));
 
79
            wx(x(Float(0.5)),                      Float(64)/Float(180));
 
80
            wx(x(0.5+0.5/sqrt(Float(3)/Float(7))), Float(49)/Float(180));
 
81
            wx(x(Float(1)),                        Float( 9)/Float(180));
 
82
            break;
 
83
     default:
 
84
            error_macro ("unsupported Gauss-Lobatto("<<opt.get_order()<<")");
 
85
    }
 
86
    return;
 
87
  }
 
88
  // -------------------------------------------------------------------------
 
89
  // default & general case : Gauss points
 
90
  // -------------------------------------------------------------------------
 
91
  assert_macro (f == quadrature_option_type::gauss, 
 
92
        "unsupported quadrature family \"" << opt.get_family_name() << "\"");
 
93
 
 
94
  // Gauss-Legendre quadrature formulae 
 
95
  //  where Legendre = Jacobi(alpha=0,beta=0) polynoms
 
96
 
 
97
  size_type n = n_node_gauss(opt.get_order());
 
98
  vector<Float> zeta(n), omega(n);
 
99
  gauss_jacobi (n, 0, 0, zeta.begin(), omega.begin());
 
100
  for (size_type i = 0; i < n; i++)
 
101
    wx (x((1+zeta[i])/2), omega[i]/2);
 
102
}