← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/pbasis/Pk_numbering.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/pbasis/Pk_numbering.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 "Pk_numbering.h"
 
22
#include "rheolef/geo.h"
 
23
#include "rheolef/rheostream.h"
 
24
namespace rheolef {
 
25
using namespace std;
 
26
 
 
27
template <class T, class M>
 
28
numbering_Pk<T,M>::numbering_Pk (std::string name)
 
29
 : _name(name), _degree()
 
30
{
 
31
  if (_name.length() > 0 && (_name[0] == 'P' || _name[0] == 'R')) {
 
32
    // TODO: check also that name fits "Pk" where is an k integer
 
33
    _degree = atoi(name.c_str()+1);
 
34
  } else if (_name.length() > 0) { // missing 'P' !
 
35
    error_macro ("invalid polynomial name `"<<_name<<"' for the Pk family");
 
36
  } else {
 
37
    // empty _name : default cstor
 
38
    _degree = 0;
 
39
  }
 
40
  base::_basis = basis (_name);
 
41
  reference_element::init_local_nnode_by_variant (_degree, _loc_ndof_by_variant);
 
42
}
 
43
template <class T, class M>
 
44
void
 
45
numbering_Pk<T,M>::set_degree (size_type degree) const
 
46
{
 
47
  _name   = "P" + itos(degree);
 
48
  _degree = degree;
 
49
  base::_basis = basis (_name);
 
50
  reference_element::init_local_nnode_by_variant (_degree, _loc_ndof_by_variant);
 
51
}
 
52
template <class T, class M>
 
53
typename numbering_Pk<T,M>::size_type
 
54
numbering_Pk<T,M>::degree () const
 
55
{
 
56
  return _degree;
 
57
}
 
58
template <class T, class M>
 
59
std::string
 
60
numbering_Pk<T,M>::name() const
 
61
{
 
62
  return _name;
 
63
}
 
64
template <class T, class M>
 
65
bool
 
66
numbering_Pk<T,M>::is_continuous () const
 
67
{
 
68
  return true;
 
69
}
 
70
template <class T, class M>
 
71
typename numbering_Pk<T,M>::size_type
 
72
numbering_Pk<T,M>::ndof (
 
73
        const geo_size& gs, 
 
74
        size_type       map_dim) const
 
75
{
 
76
  size_type ndof = 0;
 
77
  for (size_type variant = 0; variant < reference_element::max_variant; variant++) { 
 
78
    ndof += gs.ownership_by_variant [variant].size() * _loc_ndof_by_variant [variant];
 
79
  }
 
80
  return ndof;
 
81
}
 
82
template <class T, class M>
 
83
typename numbering_Pk<T,M>::size_type
 
84
numbering_Pk<T,M>::dis_ndof (
 
85
        const geo_size& gs, 
 
86
        size_type       map_dim) const
 
87
{
 
88
  size_type dis_ndof = 0;
 
89
  for (size_type variant = 0; variant < reference_element::max_variant; variant++) { 
 
90
    dis_ndof += gs.ownership_by_variant [variant].dis_size() * _loc_ndof_by_variant [variant];
 
91
  }
 
92
  return dis_ndof;
 
93
}
 
94
template <class T, class M>
 
95
void
 
96
numbering_Pk<T,M>::dis_idof (
 
97
        const geo_size&       gs, 
 
98
        const geo_element&    K, 
 
99
        std::vector<size_type>& dis_idof1) const
 
100
{
 
101
  dis_idof1.resize (reference_element::n_node (K.variant(), _degree));
 
102
  // PARANO mode:
 
103
  std::fill (dis_idof1.begin(), dis_idof1.end(), std::numeric_limits<size_type>::max());
 
104
 
 
105
  for (size_type subgeo_variant = 0, variant = K.variant(); subgeo_variant <= variant; subgeo_variant++) {
 
106
    size_type subgeo_dim   = reference_element::dimension (subgeo_variant);
 
107
    size_type loc_sub_ndof = _loc_ndof_by_variant [subgeo_variant];
 
108
    for (size_type first_loc_idof = reference_element::first_inod_by_variant (variant, _degree, subgeo_variant),
 
109
                    last_loc_idof = reference_element:: last_inod_by_variant (variant, _degree, subgeo_variant), 
 
110
                         loc_idof = first_loc_idof;
 
111
                         loc_idof <  last_loc_idof; loc_idof++) {
 
112
      // 1) local loc_igev on subgeo
 
113
      size_type loc_igev   = (loc_idof - first_loc_idof) / loc_sub_ndof;
 
114
      size_type loc_igev_j = (loc_idof - first_loc_idof) % loc_sub_ndof;
 
115
      // 2) then compute dis_ige, global distributed by dimension; computation depends upon the subgeo dimension:
 
116
      size_type dis_ige;
 
117
      switch (subgeo_dim) {
 
118
        case 0: {
 
119
          // convert node numbering to vertex numbering, for geo order > 1
 
120
          size_type loc_inod = loc_idof; // only one dof per vertex
 
121
          size_type dis_inod = K [loc_inod];
 
122
          size_type iproc           = gs.node_ownership.find_owner (dis_inod);
 
123
          size_type first_dis_inod  = gs.node_ownership.first_index(iproc);
 
124
          assert_macro (dis_inod >= first_dis_inod, "invalid index");
 
125
          size_type inod = dis_inod - first_dis_inod;
 
126
          size_type first_dis_iv    = gs.ownership_by_variant [reference_element::p].first_index(iproc);
 
127
          size_type iv = inod;
 
128
          size_type dis_iv = first_dis_iv + iv;
 
129
          dis_ige = dis_iv;
 
130
          break;
 
131
        }
 
132
        case 1: {
 
133
          loc_igev_j = geo_element::fix_edge_indirect (K, loc_igev, loc_igev_j, _degree);
 
134
          size_type loc_ige = loc_igev;
 
135
          dis_ige = K.edge(loc_ige);
 
136
          break;
 
137
        }
 
138
        case 2: {
 
139
          size_type loc_ige = loc_igev;
 
140
          if (subgeo_variant == reference_element::t) {
 
141
             loc_igev_j = geo_element::fix_triangle_indirect (K, loc_igev, loc_igev_j, _degree);
 
142
          } else {
 
143
             size_type loc_ntri = (K.variant() == reference_element::P) ? 2 : 0;
 
144
             loc_ige += loc_ntri;
 
145
             loc_igev_j = geo_element::fix_quadrangle_indirect (K, loc_ige, loc_igev_j, _degree);
 
146
          }
 
147
          dis_ige = K.face(loc_ige);
 
148
          break;
 
149
        }
 
150
        case 3:
 
151
        default: {
 
152
          dis_ige = K.dis_ie();
 
153
        }
 
154
      }
 
155
      // 3) then compute dis_igev, distributed by variant
 
156
      size_type iproc         = gs.ownership_by_dimension [subgeo_dim].find_owner (dis_ige);
 
157
      size_type first_dis_ige = gs.ownership_by_dimension [subgeo_dim].first_index(iproc);
 
158
      assert_macro (dis_ige >= first_dis_ige, "invalid index");
 
159
      size_type ige = dis_ige - first_dis_ige;
 
160
      size_type igev = ige;
 
161
      for (size_type prev_subgeo_variant = reference_element::first_variant_by_dimension(subgeo_dim);
 
162
                     prev_subgeo_variant < subgeo_variant; 
 
163
                     prev_subgeo_variant++) {
 
164
        size_type nprev = gs.ownership_by_variant [prev_subgeo_variant].size(iproc);
 
165
        assert_macro (ige >= nprev, "invalid index");
 
166
        igev -= nprev;
 
167
      } 
 
168
      size_type first_dis_igev = gs.ownership_by_variant [subgeo_variant].first_index(iproc);
 
169
      size_type dis_igev = first_dis_igev + igev;
 
170
      // 3) finally compute dis_idof
 
171
      size_type dis_idof = 0;
 
172
      for (size_type prev_subgeo_variant = 0;
 
173
                     prev_subgeo_variant < subgeo_variant; 
 
174
                     prev_subgeo_variant++) {
 
175
        dis_idof += gs.ownership_by_variant [prev_subgeo_variant]. last_index(iproc) * _loc_ndof_by_variant [prev_subgeo_variant];
 
176
      }
 
177
      dis_idof +=  dis_igev*_loc_ndof_by_variant [subgeo_variant] + loc_igev_j;
 
178
      for (size_type next_subgeo_variant = subgeo_variant+1;
 
179
                     next_subgeo_variant < reference_element::max_variant; 
 
180
                     next_subgeo_variant++) {
 
181
        dis_idof += gs.ownership_by_variant [next_subgeo_variant].first_index(iproc) * _loc_ndof_by_variant [next_subgeo_variant];
 
182
      }
 
183
      assert_macro (loc_idof < dis_idof1.size(), "invalid index");
 
184
      dis_idof1 [loc_idof] = dis_idof;
 
185
    }
 
186
  }
 
187
}
 
188
// -------------------------------------------------------------
 
189
// permut to/from ios dof numbering, for distributed environment
 
190
// -------------------------------------------------------------
 
191
template <class T, class M>
 
192
static 
 
193
void
 
194
Pk_set_ios_permutations (
 
195
    const geo_basic<T,M>& omega,
 
196
    size_t                             degree,
 
197
    array<distributor::size_type,M>&   idof2ios_dis_idof,
 
198
    array<distributor::size_type,M>&   ios_idof2dis_idof)
 
199
{
 
200
}
 
201
#ifdef _RHEOLEF_HAVE_MPI
 
202
template <class T>
 
203
static 
 
204
void
 
205
Pk_set_ios_permutations (
 
206
    const geo_basic<T,distributed>&              omega,
 
207
    size_t                                       degree,
 
208
    array<distributor::size_type,distributed>&   idof2ios_dis_idof,
 
209
    array<distributor::size_type,distributed>&   ios_idof2dis_idof)
 
210
{
 
211
  typedef size_t size_type;
 
212
  boost::array<size_type,reference_element::max_variant> loc_ndof_by_variant ;
 
213
  reference_element::init_local_nnode_by_variant (degree, loc_ndof_by_variant);
 
214
  omega.set_ios_permutation (loc_ndof_by_variant, idof2ios_dis_idof);
 
215
  size_type dis_ndof = idof2ios_dis_idof.ownership().dis_size();
 
216
  distributor ios_dof_ownership (dis_ndof, idof2ios_dis_idof.comm(), distributor::decide);
 
217
  ios_idof2dis_idof.resize (ios_dof_ownership, std::numeric_limits<size_type>::max());
 
218
  idof2ios_dis_idof.reverse_permutation (ios_idof2dis_idof);
 
219
}
 
220
#endif // _RHEOLEF_HAVE_MPI
 
221
template <class T, class M>
 
222
void
 
223
numbering_Pk<T,M>::set_ios_permutations (
 
224
    const geo_basic<T,M>& omega,
 
225
    array<size_type,M>&   idof2ios_dis_idof,
 
226
    array<size_type,M>&   ios_idof2dis_idof) const
 
227
{
 
228
    Pk_set_ios_permutations (omega, degree(), idof2ios_dis_idof, ios_idof2dis_idof);
 
229
}
 
230
// ----------------------------------------------------------------------------
 
231
// instanciation in library
 
232
// ----------------------------------------------------------------------------
 
233
 
 
234
template class numbering_Pk<Float,sequential>;
 
235
 
 
236
#ifdef _RHEOLEF_HAVE_MPI
 
237
template class numbering_Pk<Float,distributed>;
 
238
#endif // _RHEOLEF_HAVE_MPI
 
239
 
 
240
} // namespace rheolef