← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/basis/Pk_symbolic.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/Pk_symbolic.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
// P2 approximation
 
23
//
 
24
#include "rheolef/rheostream.h"
 
25
#include "basis_symbolic.h"
 
26
#include "basis_symbolic.h"
 
27
using namespace rheolef;
 
28
using namespace std;
 
29
using namespace GiNaC;
 
30
 
 
31
class Pk_symbolic : public basis_symbolic_nodal
 
32
{
 
33
    typedef basis_symbolic_nodal::size_type size_type;
 
34
public:
 
35
    Pk_symbolic (size_type k);
 
36
};
 
37
Pk_symbolic::Pk_symbolic (size_type degree)
 
38
: basis_symbolic_nodal("P"+itos(degree),degree)
 
39
{
 
40
  typedef point_basic<size_type> ilat;
 
41
  std::vector<point_basic<ex> >  hat_xnod;
 
42
  // ------------------------------------
 
43
  on('p') << node(0) << poly (1) << end;
 
44
  // ------------------------------------
 
45
  hat_xnod.resize (reference_element::n_node(reference_element::e, degree));
 
46
  for (size_type i = 0; i <= degree; i++) { 
 
47
    on('e') << poly (pow(x,i));
 
48
    size_type loc_idof = reference_element_e::ilat2loc_inod (degree, ilat(i));
 
49
    hat_xnod [loc_idof] = point_basic<ex> ((ex(i))/ex(degree));
 
50
  }
 
51
  for (size_type loc_idof = 0, loc_ndof = hat_xnod.size(); loc_idof < loc_ndof; loc_idof++) { 
 
52
    on('e') << hat_xnod [loc_idof];
 
53
  }
 
54
  on('e') << end;
 
55
  // ------------------------------------
 
56
  hat_xnod.resize (reference_element::n_node(reference_element::t, degree));
 
57
  for (size_type j = 0; j <= degree; j++) { 
 
58
    for (size_type i = 0; i+j <= degree; i++) { 
 
59
      on('t') << poly (pow(x,i)*pow(y,j));
 
60
      size_type loc_idof = reference_element_t::ilat2loc_inod (degree, ilat(i,j));
 
61
      hat_xnod [loc_idof] = point_basic<ex> ((ex(i))/ex(degree), ex(j)/ex(degree));
 
62
    }
 
63
  }
 
64
  for (size_type loc_idof = 0, loc_ndof = hat_xnod.size(); loc_idof < loc_ndof; loc_idof++) { 
 
65
    on('t') << hat_xnod [loc_idof];
 
66
  }
 
67
  on('t') << end;
 
68
  // ------------------------------------
 
69
  hat_xnod.resize (reference_element::n_node(reference_element::q, degree));
 
70
  for (size_type j = 0; j <= degree; j++) { 
 
71
    for (size_type i = 0; i <= degree; i++) { 
 
72
      on('q') << poly (pow(x,i)*pow(y,j));
 
73
      size_type loc_idof = reference_element_q::ilat2loc_inod (degree, ilat(i,j));
 
74
      hat_xnod [loc_idof] = point_basic<ex> (-1+2*(ex(i))/ex(degree), -1+2*ex(j)/ex(degree));
 
75
    }
 
76
  }
 
77
  for (size_type loc_idof = 0, loc_ndof = hat_xnod.size(); loc_idof < loc_ndof; loc_idof++) { 
 
78
    on('q') << hat_xnod [loc_idof];
 
79
  }
 
80
  on('q') << end;
 
81
  // -----------------------------------------
 
82
  if (degree >= 5) return; // too big: not yet
 
83
  // ------------------------------------
 
84
  hat_xnod.resize (reference_element::n_node(reference_element::T, degree));
 
85
  for (size_type k = 0; k <= degree; k++) {
 
86
    for (size_type j = 0; j+k <= degree; j++) {
 
87
      for (size_type i = 0; i+j+k <= degree; i++) { 
 
88
        on('T') << poly (pow(x,i)*pow(y,j)*pow(z,k));
 
89
        size_type loc_idof = reference_element_T::ilat2loc_inod (degree, ilat(i,j,k));
 
90
        hat_xnod [loc_idof] = point_basic<ex> ((ex(i))/ex(degree), ex(j)/ex(degree), ex(k)/ex(degree));
 
91
      }
 
92
    }
 
93
  }
 
94
  for (size_type loc_idof = 0, loc_ndof = hat_xnod.size(); loc_idof < loc_ndof; loc_idof++) { 
 
95
    on('T') << hat_xnod [loc_idof];
 
96
  }
 
97
  on('T') << end;
 
98
  // -----------------------------------------
 
99
  hat_xnod.resize (reference_element::n_node(reference_element::P, degree));
 
100
  for (size_type k = 0; k <= degree; k++) { 
 
101
    for (size_type j = 0; j <= degree; j++) { 
 
102
      for (size_type i = 0; i+j <= degree; i++) { 
 
103
        on('P') << poly (pow(x,i)*pow(y,j)*pow(z,k));
 
104
        size_type loc_idof = reference_element_P::ilat2loc_inod (degree, ilat(i,j,k));
 
105
        hat_xnod [loc_idof] = point_basic<ex> ((ex(i))/ex(degree), ex(j)/ex(degree), -1+2*ex(k)/ex(degree));
 
106
      }
 
107
    }
 
108
  }
 
109
  for (size_type loc_idof = 0, loc_ndof = hat_xnod.size(); loc_idof < loc_ndof; loc_idof++) { 
 
110
    on('P') << hat_xnod [loc_idof];
 
111
  }
 
112
  on('P') << end;
 
113
  // ------------------------------------
 
114
  hat_xnod.resize (reference_element::n_node(reference_element::H, degree));
 
115
  for (size_type k = 0; k <= degree; k++) { 
 
116
    for (size_type j = 0; j <= degree; j++) { 
 
117
      for (size_type i = 0; i <= degree; i++) { 
 
118
        on('H') << poly (pow(x,i)*pow(y,j)*pow(z,k));
 
119
        size_type loc_idof = reference_element_H::ilat2loc_inod (degree, ilat(i,j,k));
 
120
        hat_xnod [loc_idof] = point_basic<ex> (-1+2*(ex(i))/ex(degree), -1+2*ex(j)/ex(degree), -1+2*ex(k)/ex(degree));
 
121
      }
 
122
    }
 
123
  }
 
124
  for (size_type loc_idof = 0, loc_ndof = hat_xnod.size(); loc_idof < loc_ndof; loc_idof++) { 
 
125
    on('H') << hat_xnod [loc_idof];
 
126
  }
 
127
  on('H') << end;
 
128
}
 
129
int main (int argc, char **argv) {
 
130
        if (argc < 3) {
 
131
           cerr << "Pk_symbolic: usage: Pk_symbolic [-h|-c] degree" << endl;
 
132
           exit (1);
 
133
        }
 
134
        size_t degree = atoi (argv[2]);
 
135
        Pk_symbolic Pk (degree);
 
136
        Pk.put_cxx_main (argc,argv);
 
137
}