Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
basis_ordering.icc
Go to the documentation of this file.
1#ifndef _RHEOLEF_BASIS_ORDERING_ICC
2#define _RHEOLEF_BASIS_ORDERING_ICC
23//
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
44namespace 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[.]
50static
51size_t
52get_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;
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
80static
81void
82sort_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
108static
109void
110sort_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}
121static
122void
123invert_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)
138static
139void
140make_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()) {
150 power_index[0] = point_basic<size_t>();
151 break;
152 }
154 for (size_type ilat = 0; ilat <= degree; ilat++) {
155 power_index[ilat2ipow[ilat]] = point_basic<size_t>(ilat);
156 }
157 break;
158 }
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 }
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 }
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 }
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 }
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:
208static
209void
210make_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}
223static
224void
225make_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}
239static
240void
241make_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
257static
258void
259build_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
field::size_type size_type
Definition branch.cc:430
static const variant_type H
static const variant_type q
static const variant_type e
static const variant_type p
std::vector< int >::size_type size_type
static size_type n_node(variant_type variant, size_type order)
static const variant_type T
static const variant_type P
static const variant_type t
This file is part of Rheolef.