← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/basis/quadrature_Pr.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/basis/quadrature_Pr.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
 
void
27
 
quadrature_on_geo::init_prism  (quadrature_option_type opt)
28
 
{
29
 
  quadrature_option_type::family_type f = opt.get_family();
30
 
  // -------------------------------------------------------------------------
31
 
  // general case: Gauss
32
 
  // -------------------------------------------------------------------------
33
 
  assert_macro (f == quadrature_option_type::gauss,
34
 
        "unsupported quadrature family \"" << opt.get_family_name() << "\"");
35
 
 
36
 
  typedef Float T;
37
 
  switch (opt.get_order()) {
38
 
   case 0:
39
 
   case 1: {
40
 
        // central point:
41
 
        // better than the general case
42
 
        // r=0 : 2 nodes
43
 
        // r=1 : 8 nodes
44
 
        // here: 1 node
45
 
        T a = T(1)/3;
46
 
        wx(x(a,a,0), 1);
47
 
        break;
48
 
   }
49
 
   default: {
50
 
    // Gauss-Legendre quadrature formulae 
51
 
    //  where Legendre = Jacobi(alpha=0,beta=0) polynoms
52
 
    size_t r = opt.get_order();
53
 
    size_type n0 = n_node_gauss(r);
54
 
    size_type n1 = n_node_gauss(r+1);
55
 
    size_type n2 = n_node_gauss(r);
56
 
    vector<Float> zeta0(n0), omega0(n0);
57
 
    vector<Float> zeta1(n1), omega1(n1);
58
 
    vector<Float> zeta2(n2), omega2(n2);
59
 
    gauss_jacobi (n0, 0, 0, zeta0.begin(), omega0.begin());
60
 
    gauss_jacobi (n1, 0, 0, zeta1.begin(), omega1.begin());
61
 
    gauss_jacobi (n2, 0, 0, zeta2.begin(), omega2.begin());
62
 
 
63
 
    for (size_type i = 0; i < n0; i++) {
64
 
      for (size_type j = 0; j < n1; j++) {
65
 
        for (size_type k = 0; k < n2; k++) {
66
 
          // we transform the square into the triangle 
67
 
          Float eta_0 = (1+zeta0[i])*(1-zeta1[j])/4;
68
 
          Float eta_1 =              (1+zeta1[j])/2;
69
 
          Float eta_2 = zeta2[k];
70
 
          Float J     =              (1-zeta1[j])/8;
71
 
          wx (x(eta_0,eta_1,eta_2), J*omega0[i]*omega1[j]*omega2[k]);
72
 
        }
73
 
      }
74
 
    }
75
 
   }
76
 
  }
77
 
}
78