Rheolef  7.1
an efficient C++ finite element environment
basis_ordering.icc
Go to the documentation of this file.
1 #ifndef _RHEOLEF_BASIS_ORDERING_ICC
2 #define _RHEOLEF_BASIS_ORDERING_ICC
3 //
24 // utilities for local ordering, at the element level
25 // 1) lattice ordering, as (i,j), 0 <= i < n, 0 <= j < n-i
26 // => used by algorithms that loop locally
27 // 2) node ordering, by increasing sub-geometry dimension
28 // first, nodes on vertices,
29 // second, nodes on interior of edges
30 // third, nodes on interior of faces
31 // last, nodes on interior of volume
32 // => used for FEM basis coefficients
33 // i.e. for degree-of-freedom ordering
34 // 3) polynomial degree increasing ordering
35 // => used for RAW (initial) basis coefficients
36 //
37 // author: Pierre.Saramito@imag.fr
38 //
39 // date: 2 september 2017
40 //
41 #include "rheolef/reference_element.h"
43 
44 namespace rheolef {
45 
46 // p(x0..xd) = x0^m0*..xd^md
47 // compute its degree, as it depends on element shape
48 // * simplicial elements: sum of m[.]
49 // * tensorial-product elements: max of m[.]
50 static
51 size_t
52 get_degree (
53  reference_element hat_K,
54  const point_basic<size_t>& power_index)
55 {
56  size_t d = hat_K.dimension();
57  size_t deg_p = 0;
58  for (size_t mu = 0; mu < d; ++mu) {
59  switch (hat_K.variant()) {
63  case reference_element::T: deg_p += power_index[mu]; break;
65  case reference_element::H: deg_p = max(deg_p, power_index[mu]); break;
66  case reference_element::P: {
67  if (mu != 2) deg_p += power_index[mu];
68  else deg_p = max(deg_p, power_index[mu]);
69  break;
70  }
71  default: break;
72  }
73  }
74  return deg_p;
75 }
76 // p(x0..xd) = x0^m0*..xd^md
77 // re-order monoms by polynomial degree
78 // note: this is useful for all others hierchical basis
79 // such as Bernstein and Dubiner
80 static
81 void
82 sort_by_degrees (
83  reference_element hat_K,
84  size_t degree,
85  const std::vector<point_basic<size_t> >& power_index,
86  std::vector<size_t>& ideg2ipow)
87 {
88  ideg2ipow.resize (power_index.size());
89  std::vector<size_t> first_loc_ideg (degree+1);
90  first_loc_ideg [0] = 0;
91  for (size_t k = 1; k <= degree; ++k) {
92  first_loc_ideg [k] = reference_element::n_node(hat_K.variant(), k-1);
93  }
94  for (size_t loc_ipow = 0, loc_npow = power_index.size(); loc_ipow < loc_npow; ++loc_ipow) {
95  size_t deg = get_degree (hat_K, power_index[loc_ipow]);
96  size_t ideg = first_loc_ideg[deg];
97  ideg2ipow[ideg] = loc_ipow;
98  first_loc_ideg[deg]++;
99  }
100 }
101 // p(x0..xd) = x0^m0*..xd^md
102 // re-order monoms by Lagrange node ordering
103 // these are numbered by increasing node dimension
104 // 1rst : nodes on corners
105 // 2nd : nodes on edges
106 // 3rd : nodes on faces
107 // 4rd : nodes on volumes
108 static
109 void
110 sort_by_inodes (
111  reference_element hat_K,
112  size_t degree,
113  const std::vector<point_basic<size_t> >& power_index,
114  std::vector<size_t>& ipow2inod)
115 {
116  ipow2inod.resize (power_index.size());
117  for (size_t ipow = 0, npow = power_index.size(); ipow < npow; ++ipow) {
118  ipow2inod[ipow] = ilat2loc_inod (hat_K, degree, power_index[ipow]);
119  }
120 }
121 static
122 void
123 invert_permutation (
124  const std::vector<size_t>& perm,
125  std::vector<size_t>& iperm)
126 {
127  iperm.resize (perm.size());
128  for (size_t iloc = 0, nloc = perm.size(); iloc < nloc; ++iloc) {
129  iperm[perm[iloc]] = iloc;
130  }
131 }
132 // ------------------------------------------
133 // power indexes
134 // ------------------------------------------
135 // p(x0..xd) = x0^m0*..xd^md
136 // => compute all power index monomials m[]
137 // ordering follows the natural lattice order (sorted by lattice ordering)
138 static
139 void
140 make_power_indexes_permuted (
141  reference_element hat_K,
142  size_t degree,
143  std::vector<size_t> ilat2ipow,
144  std::vector<point_basic<size_t> >& power_index)
145 {
147  power_index.resize (reference_element::n_node(hat_K.variant(), degree));
148  switch (hat_K.variant()) {
149  case reference_element::p: {
150  power_index[0] = point_basic<size_t>();
151  break;
152  }
153  case reference_element::e: {
154  for (size_type ilat = 0; ilat <= degree; ilat++) {
155  power_index[ilat2ipow[ilat]] = point_basic<size_t>(ilat);
156  }
157  break;
158  }
159  case reference_element::t: {
160  size_type ilat = 0;
161  for (size_type j = 0; j <= degree; j++) {
162  for (size_type i = 0; i+j <= degree; i++) {
163  power_index[ilat2ipow[ilat++]] = point_basic<size_t>(i,j);
164  }}
165  break;
166  }
167  case reference_element::q: {
168  size_type ilat = 0;
169  for (size_type j = 0; j <= degree; j++) {
170  for (size_type i = 0; i <= degree; i++) {
171  power_index[ilat2ipow[ilat++]] = point_basic<size_t>(i,j);
172  }}
173  break;
174  }
175  case reference_element::T: {
176  for (size_type k = 0, ilat = 0; k <= degree; k++) {
177  for (size_type j = 0; j+k <= degree; j++) {
178  for (size_type i = 0; i+j+k <= degree; i++, ilat++) {
179  power_index[ilat2ipow[ilat]] = point_basic<size_t>(i,j,k);
180  }
181  }
182  }
183  break;
184  }
185  case reference_element::P: {
186  for (size_type k = 0, ilat = 0; k <= degree; k++) {
187  for (size_type j = 0; j <= degree; j++) {
188  for (size_type i = 0; i+j <= degree; i++, ilat++) {
189  power_index[ilat2ipow[ilat]] = point_basic<size_t>(i,j,k);
190  }
191  }
192  }
193  break;
194  }
195  case reference_element::H: {
196  for (size_type k = 0, ilat = 0; k <= degree; k++) {
197  for (size_type j = 0; j <= degree; j++) {
198  for (size_type i = 0; i <= degree; i++, ilat++) {
199  power_index[ilat2ipow[ilat]] = point_basic<size_t>(i,j,k);
200  }
201  }
202  }
203  break;
204  }
205  }
206 }
207 // power indexes sorted by lattice ordering:
208 static
209 void
210 make_power_indexes (
211  reference_element hat_K,
212  size_t degree,
213  std::vector<point_basic<size_t> >& power_index)
214 {
216  size_type nloc = reference_element::n_node(hat_K.variant(), degree);
217  std::vector<size_type> id (nloc);
218  for (size_type iloc = 0; iloc < nloc; ++iloc) {
219  id [iloc] = iloc;
220  }
221  make_power_indexes_permuted (hat_K, degree, id, power_index);
222 }
223 static
224 void
225 make_power_indexes_sorted_by_degrees (
226  reference_element hat_K,
227  size_t degree,
228  std::vector<point_basic<size_t> >& power_index)
229 {
231  size_type nloc = reference_element::n_node(hat_K.variant(), degree);
232  make_power_indexes (hat_K, degree, power_index);
233  std::vector<size_type> perm (nloc);
234  sort_by_degrees (hat_K, degree, power_index, perm);
235  std::vector<size_type> iperm (nloc);
236  invert_permutation (perm, iperm);
237  make_power_indexes_permuted (hat_K, degree, iperm, power_index);
238 }
239 static
240 void
241 make_power_indexes_sorted_by_inodes (
242  reference_element hat_K,
243  size_t degree,
244  std::vector<point_basic<size_t> >& power_index)
245 {
247  size_type nloc = reference_element::n_node(hat_K.variant(), degree);
248  make_power_indexes (hat_K, degree, power_index);
249  std::vector<size_type> perm (nloc);
250  sort_by_inodes (hat_K, degree, power_index, perm);
251  make_power_indexes_permuted (hat_K, degree, perm, power_index);
252 }
253 // ------------------------------------------
254 // inod2ideg permutation
255 // ------------------------------------------
256 // use ilat ordering during the build
257 static
258 void
259 build_inod2ideg (
260  reference_element hat_K,
261  size_t degree,
262  std::vector<size_t>& inod2ideg)
263 {
264  std::vector<point_basic<size_t> > power_index;
265  make_power_indexes (hat_K, degree, power_index);
266  std::vector<size_t> ideg2ilat;
267  std::vector<size_t> ilat2inod;
268  sort_by_inodes (hat_K, degree, power_index, ilat2inod);
269  sort_by_degrees (hat_K, degree, power_index, ideg2ilat);
270  std::vector<size_t> ilat2ideg;
271  std::vector<size_t> inod2ilat;
272  invert_permutation (ilat2inod, inod2ilat);
273  invert_permutation (ideg2ilat, ilat2ideg);
274  inod2ideg.resize (power_index.size());
275  for (size_t iloc = 0, nloc = inod2ideg.size(); iloc < nloc; ++iloc) {
276  inod2ideg [iloc] = ilat2ideg [inod2ilat[iloc]];
277  }
278 }
279 
280 }// namespace rheolef
281 #endif // _RHEOLEF_BASIS_ORDERING_ICC
rheolef::reference_element::e
static const variant_type e
Definition: reference_element.h:76
rheolef::reference_element::H
static const variant_type H
Definition: reference_element.h:81
rheolef::reference_element::n_node
static size_type n_node(variant_type variant, size_type order)
Definition: reference_element.cc:201
rheolef::reference_element::T
static const variant_type T
Definition: reference_element.h:79
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
reference_element_aux.icc
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
rheolef::reference_element::p
static const variant_type p
Definition: reference_element.h:75
rheolef::reference_element::q
static const variant_type q
Definition: reference_element.h:78
rheolef::reference_element::P
static const variant_type P
Definition: reference_element.h:80
rheolef::space_constant::vector
Definition: space_constant.h:137
mkgeo_contraction.mu
mu
Definition: mkgeo_contraction.sh:193
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::reference_element::size_type
std::vector< int >::size_type size_type
Definition: reference_element.h:71