basis_ordering.icc
Go to the documentation of this file.
#ifndef _RHEOLEF_BASIS_ORDERING_ICC
#define _RHEOLEF_BASIS_ORDERING_ICC
//
// utilities for local ordering, at the element level
// 1) lattice ordering, as (i,j), 0 <= i < n, 0 <= j < n-i
// => used by algorithms that loop locally
// 2) node ordering, by increasing sub-geometry dimension
// first, nodes on vertices,
// second, nodes on interior of edges
// third, nodes on interior of faces
// last, nodes on interior of volume
// => used for FEM basis coefficients
// i.e. for degree-of-freedom ordering
// 3) polynomial degree increasing ordering
// => used for RAW (initial) basis coefficients
//
// author: Pierre.Saramito@imag.fr
//
// date: 2 september 2017
//
#include "rheolef/reference_element.h"
#include "reference_element_aux.icc"
namespace rheolef {
// p(x0..xd) = x0^m0*..xd^md
// compute its degree, as it depends on element shape
// * simplicial elements: sum of m[.]
// * tensorial-product elements: max of m[.]
static
size_t
get_degree (
reference_element hat_K,
const point_basic<size_t>& power_index)
{
size_t d = hat_K.dimension();
size_t deg_p = 0;
for (size_t mu = 0; mu < d; ++mu) {
switch (hat_K.variant()) {
case reference_element::p:
case reference_element::e:
case reference_element::t:
case reference_element::T: deg_p += power_index[mu]; break;
case reference_element::q:
case reference_element::H: deg_p = max(deg_p, power_index[mu]); break;
case reference_element::P: {
if (mu != 2) deg_p += power_index[mu];
else deg_p = max(deg_p, power_index[mu]);
break;
}
default: break;
}
}
return deg_p;
}
// p(x0..xd) = x0^m0*..xd^md
// re-order monoms by polynomial degree
// note: this is useful for all others hierchical basis
// such as Bernstein and Dubiner
static
void
sort_by_degrees (
reference_element hat_K,
size_t degree,
const std::vector<point_basic<size_t> >& power_index,
std::vector<size_t>& ideg2ipow)
{
ideg2ipow.resize (power_index.size());
std::vector<size_t> first_loc_ideg (degree+1);
first_loc_ideg [0] = 0;
for (size_t k = 1; k <= degree; ++k) {
first_loc_ideg [k] = reference_element::n_node(hat_K.variant(), k-1);
}
for (size_t loc_ipow = 0, loc_npow = power_index.size(); loc_ipow < loc_npow; ++loc_ipow) {
size_t deg = get_degree (hat_K, power_index[loc_ipow]);
size_t ideg = first_loc_ideg[deg];
ideg2ipow[ideg] = loc_ipow;
first_loc_ideg[deg]++;
}
}
// p(x0..xd) = x0^m0*..xd^md
// re-order monoms by Lagrange node ordering
// these are numbered by increasing node dimension
// 1rst : nodes on corners
// 2nd : nodes on edges
// 3rd : nodes on faces
// 4rd : nodes on volumes
static
void
sort_by_inodes (
reference_element hat_K,
size_t degree,
const std::vector<point_basic<size_t> >& power_index,
std::vector<size_t>& ipow2inod)
{
ipow2inod.resize (power_index.size());
for (size_t ipow = 0, npow = power_index.size(); ipow < npow; ++ipow) {
ipow2inod[ipow] = ilat2loc_inod (hat_K, degree, power_index[ipow]);
}
}
static
void
invert_permutation (
const std::vector<size_t>& perm,
std::vector<size_t>& iperm)
{
iperm.resize (perm.size());
for (size_t iloc = 0, nloc = perm.size(); iloc < nloc; ++iloc) {
iperm[perm[iloc]] = iloc;
}
}
// ------------------------------------------
// power indexes
// ------------------------------------------
// p(x0..xd) = x0^m0*..xd^md
// => compute all power index monomials m[]
// ordering follows the natural lattice order (sorted by lattice ordering)
static
void
make_power_indexes_permuted (
reference_element hat_K,
size_t degree,
std::vector<size_t> ilat2ipow,
std::vector<point_basic<size_t> >& power_index)
{
typedef reference_element::size_type size_type;
power_index.resize (reference_element::n_node(hat_K.variant(), degree));
switch (hat_K.variant()) {
case reference_element::p: {
power_index[0] = point_basic<size_t>();
break;
}
case reference_element::e: {
for (size_type ilat = 0; ilat <= degree; ilat++) {
power_index[ilat2ipow[ilat]] = point_basic<size_t>(ilat);
}
break;
}
case reference_element::t: {
size_type ilat = 0;
for (size_type j = 0; j <= degree; j++) {
for (size_type i = 0; i+j <= degree; i++) {
power_index[ilat2ipow[ilat++]] = point_basic<size_t>(i,j);
}}
break;
}
case reference_element::q: {
size_type ilat = 0;
for (size_type j = 0; j <= degree; j++) {
for (size_type i = 0; i <= degree; i++) {
power_index[ilat2ipow[ilat++]] = point_basic<size_t>(i,j);
}}
break;
}
case reference_element::T: {
for (size_type k = 0, ilat = 0; k <= degree; k++) {
for (size_type j = 0; j+k <= degree; j++) {
for (size_type i = 0; i+j+k <= degree; i++, ilat++) {
power_index[ilat2ipow[ilat]] = point_basic<size_t>(i,j,k);
}
}
}
break;
}
case reference_element::P: {
for (size_type k = 0, ilat = 0; k <= degree; k++) {
for (size_type j = 0; j <= degree; j++) {
for (size_type i = 0; i+j <= degree; i++, ilat++) {
power_index[ilat2ipow[ilat]] = point_basic<size_t>(i,j,k);
}
}
}
break;
}
case reference_element::H: {
for (size_type k = 0, ilat = 0; k <= degree; k++) {
for (size_type j = 0; j <= degree; j++) {
for (size_type i = 0; i <= degree; i++, ilat++) {
power_index[ilat2ipow[ilat]] = point_basic<size_t>(i,j,k);
}
}
}
break;
}
}
}
// power indexes sorted by lattice ordering:
static
void
make_power_indexes (
reference_element hat_K,
size_t degree,
std::vector<point_basic<size_t> >& power_index)
{
typedef reference_element::size_type size_type;
size_type nloc = reference_element::n_node(hat_K.variant(), degree);
std::vector<size_type> id (nloc);
for (size_type iloc = 0; iloc < nloc; ++iloc) {
id [iloc] = iloc;
}
make_power_indexes_permuted (hat_K, degree, id, power_index);
}
static
void
make_power_indexes_sorted_by_degrees (
reference_element hat_K,
size_t degree,
std::vector<point_basic<size_t> >& power_index)
{
typedef reference_element::size_type size_type;
size_type nloc = reference_element::n_node(hat_K.variant(), degree);
make_power_indexes (hat_K, degree, power_index);
std::vector<size_type> perm (nloc);
sort_by_degrees (hat_K, degree, power_index, perm);
std::vector<size_type> iperm (nloc);
invert_permutation (perm, iperm);
make_power_indexes_permuted (hat_K, degree, iperm, power_index);
}
static
void
make_power_indexes_sorted_by_inodes (
reference_element hat_K,
size_t degree,
std::vector<point_basic<size_t> >& power_index)
{
typedef reference_element::size_type size_type;
size_type nloc = reference_element::n_node(hat_K.variant(), degree);
make_power_indexes (hat_K, degree, power_index);
std::vector<size_type> perm (nloc);
sort_by_inodes (hat_K, degree, power_index, perm);
make_power_indexes_permuted (hat_K, degree, perm, power_index);
}
// ------------------------------------------
// inod2ideg permutation
// ------------------------------------------
// use ilat ordering during the build
static
void
build_inod2ideg (
reference_element hat_K,
size_t degree,
std::vector<size_t>& inod2ideg)
{
std::vector<point_basic<size_t> > power_index;
make_power_indexes (hat_K, degree, power_index);
std::vector<size_t> ideg2ilat;
std::vector<size_t> ilat2inod;
sort_by_inodes (hat_K, degree, power_index, ilat2inod);
sort_by_degrees (hat_K, degree, power_index, ideg2ilat);
std::vector<size_t> ilat2ideg;
std::vector<size_t> inod2ilat;
invert_permutation (ilat2inod, inod2ilat);
invert_permutation (ideg2ilat, ilat2ideg);
inod2ideg.resize (power_index.size());
for (size_t iloc = 0, nloc = inod2ideg.size(); iloc < nloc; ++iloc) {
inod2ideg [iloc] = ilat2ideg [inod2ilat[iloc]];
}
}
}// namespace rheolef
#endif // _RHEOLEF_BASIS_ORDERING_ICC
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::reference_element::H
static const variant_type H
Definition: reference_element.h:81
rheolef::reference_element::q
static const variant_type q
Definition: reference_element.h:78
rheolef::reference_element::e
static const variant_type e
Definition: reference_element.h:76
rheolef::reference_element::p
static const variant_type p
Definition: reference_element.h:75
rheolef::reference_element::size_type
std::vector< int >::size_type size_type
Definition: reference_element.h:71
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::reference_element::P
static const variant_type P
Definition: reference_element.h:80
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
reference_element_aux.icc
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
mkgeo_contraction.mu
mu
Definition: mkgeo_contraction.sh:193
rheolef::space_constant::vector
@ vector
Definition: space_constant.h:137
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37