← Back to branch summary

~ubuntu-branches/ubuntu/raring/rheolef/raring-proposed

« back to all changes in this revision

Viewing changes to nfem/geo_element/quadrature_Pr.cc

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito, Pierre Saramito, Sylvestre Ledru
  • Date: 2012-05-14 14:02:09 UTC
  • mfrom: (1.1.6)
  • Revision ID: package-import@ubuntu.com-20120514140209-dzbdlidkotyflf9e
Tags: 6.1-1
http://bugs.debian.org/671996
[ Pierre Saramito ]
* New upstream release 6.1 (minor changes):
  - support arbitrarily polynomial order Pk
  - source code supports g++-4.7 (closes: #671996)

[ Sylvestre Ledru ]
* update of the watch file

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/geo_element/quadrature_Pr.cc
20
20
/// =========================================================================
21
21
#include "rheolef/quadrature.h"
22
22
#include "rheolef/gauss_jacobi.h"
23
 
using namespace rheolef;
 
23
namespace rheolef {
24
24
using namespace std;
25
25
 
 
26
template<class T>
26
27
void
27
 
quadrature_on_geo::init_prism  (quadrature_option_type opt)
 
28
quadrature_on_geo<T>::init_prism  (quadrature_option_type opt)
28
29
{
29
30
  quadrature_option_type::family_type f = opt.get_family();
30
31
  // -------------------------------------------------------------------------
33
34
  assert_macro (f == quadrature_option_type::gauss,
34
35
        "unsupported quadrature family \"" << opt.get_family_name() << "\"");
35
36
 
36
 
  typedef Float T;
37
37
  switch (opt.get_order()) {
38
38
   case 0:
39
39
   case 1: {
53
53
    size_type n0 = n_node_gauss(r);
54
54
    size_type n1 = n_node_gauss(r+1);
55
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);
 
56
    vector<T> zeta0(n0), omega0(n0);
 
57
    vector<T> zeta1(n1), omega1(n1);
 
58
    vector<T> zeta2(n2), omega2(n2);
59
59
    gauss_jacobi (n0, 0, 0, zeta0.begin(), omega0.begin());
60
60
    gauss_jacobi (n1, 0, 0, zeta1.begin(), omega1.begin());
61
61
    gauss_jacobi (n2, 0, 0, zeta2.begin(), omega2.begin());
64
64
      for (size_type j = 0; j < n1; j++) {
65
65
        for (size_type k = 0; k < n2; k++) {
66
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;
 
67
          T eta_0 = (1+zeta0[i])*(1-zeta1[j])/4;
 
68
          T eta_1 =              (1+zeta1[j])/2;
 
69
          T eta_2 = zeta2[k];
 
70
          T J     =              (1-zeta1[j])/8;
71
71
          wx (x(eta_0,eta_1,eta_2), J*omega0[i]*omega1[j]*omega2[k]);
72
72
        }
73
73
      }
75
75
   }
76
76
  }
77
77
}
78
 
 
 
78
// ----------------------------------------------------------------------------
 
79
// instanciation in library
 
80
// ----------------------------------------------------------------------------
 
81
#define _RHEOLEF_instanciation(T)                                               \
 
82
template void quadrature_on_geo<T>::init_prism (quadrature_option_type);
 
83
 
 
84
_RHEOLEF_instanciation(Float)
 
85
 
 
86
}// namespace rheolef