Loading...
Searching...
No Matches
reference_element_aux.icc
Go to the documentation of this file.
#ifndef _RHEOLEF_REFERENCE_ELEMENT_AUX_ICC
#define _RHEOLEF_REFERENCE_ELEMENT_AUX_ICC
// minimal auto-contained rheolef helpers for msh2geo
// used in:
// - reference_element.cc
// - msh2geo.cc
// in order to avoid code duplication
//
namespace rheolef {
// used by msh2geo.cc that do not include reference_element.h
#ifndef _RHEOLEF_REFERENCE_ELEMENT_H
constexpr size_t reference_element__p = 0;
constexpr size_t reference_element__e = 1;
constexpr size_t reference_element__t = 2;
constexpr size_t reference_element__q = 3;
constexpr size_t reference_element__T = 4;
constexpr size_t reference_element__P = 5;
constexpr size_t reference_element__H = 6;
constexpr size_t reference_element__max_variant = 7;
#else // _RHEOLEF_REFERENCE_ELEMENT_H
constexpr size_t reference_element__p = reference_element::p;
constexpr size_t reference_element__e = reference_element::e;
constexpr size_t reference_element__t = reference_element::t;
constexpr size_t reference_element__q = reference_element::q;
constexpr size_t reference_element__T = reference_element::T;
constexpr size_t reference_element__P = reference_element::P;
constexpr size_t reference_element__H = reference_element::H;
constexpr size_t reference_element__max_variant = reference_element::max_variant;
#endif // _RHEOLEF_REFERENCE_ELEMENT_H
static size_t _reference_element_dimension_by_variant[reference_element__max_variant] = {0, 1, 2, 2, 3, 3, 3};
static size_t _reference_element_n_edge_by_variant [reference_element__max_variant] = {0, 1, 3, 4, 6, 9, 12};
static size_t _reference_element_n_face_by_variant [reference_element__max_variant] = {0, 0, 1, 1, 4, 5, 6};
static
size_t
reference_element_variant (const char* name_table, char name)
{
for (size_t i = 0; i < reference_element__max_variant; i++) {
if (name_table [i] == name) return i;
}
return reference_element__max_variant;
}
static
size_t
reference_element_dimension_by_name(char element_name) {
switch (element_name) {
case 'p': return 0;
case 'e': return 1;
case 't':
case 'q': return 2;
case 'T':
case 'P':
case 'H': return 3;
default: error_macro("unexpected element name `"<<element_name<<"' (ascii="<<int(element_name)<<")");
return 0;
}
}
static
size_t
reference_element_dimension_by_variant (size_t variant)
{
return _reference_element_dimension_by_variant[variant];
}
static
size_t
reference_element_variant (char element_name)
{
switch (element_name) {
case 'p': return reference_element__p;
case 'e': return reference_element__e;
case 't': return reference_element__t;
case 'q': return reference_element__q;
case 'T': return reference_element__T;
case 'P': return reference_element__P;
case 'H': return reference_element__H;
default: error_macro("unexpected element name `"<<element_name<<"'\"=char("<<int(element_name)<<")");
return 0;
}
}
static
size_t
reference_element_n_vertex (size_t variant)
{
switch (variant) {
case reference_element__p: return 1;
case reference_element__e: return 2;
case reference_element__t: return 3;
case reference_element__q: return 4;
case reference_element__T: return 4;
case reference_element__P: return 6;
case reference_element__H: return 8;
default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
}
}
static
size_t
reference_element_n_node (size_t variant, size_t order)
{
switch (variant) {
case reference_element__p: return 1;
case reference_element__e: return order+1;
case reference_element__t: return (order+1)*(order+2)/2;
case reference_element__q: return (order+1)*(order+1);
case reference_element__T: return (order+1)*(order+2)*(order+3)/6;
case reference_element__P: return (order+1)*(order+1)*(order+2)/2;
case reference_element__H: return (order+1)*(order+1)*(order+1);
default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
}
}
static
void
reference_element_init_local_nnode_by_variant (
size_t order,
std::array<size_t,reference_element__max_variant>& sz)
{
check_macro (order > 0, "invalid order="<<order<<": expect order > 0");
sz [reference_element__p] = 1;
sz [reference_element__e] = order-1;
sz [reference_element__t] = (order-1)*(order-2)/2;
sz [reference_element__q] = (order-1)*(order-1);
sz [reference_element__T] = (order-1)*(order-2)*(order-3)/6;
sz [reference_element__P] = (order-1)*(order-1)*(order-2)/2;
sz [reference_element__H] = (order-1)*(order-1)*(order-1);
}
// -----------------------------------------------------------------------------
// reference_element::first_inod_by_dimension
// -----------------------------------------------------------------------------
static const size_t
_first_variant_by_dimension [5] = {
reference_element__p,
reference_element__e,
reference_element__t,
reference_element__T,
reference_element__max_variant
};
static
size_t
reference_element_first_variant_by_dimension (size_t dim)
{
return _first_variant_by_dimension[dim];
}
static
size_t
reference_element_last_variant_by_dimension (size_t dim)
{
return _first_variant_by_dimension[dim+1];
}
// -----------------------------------------------------------------------------
// reference_element::first_inod_by_variant
// -----------------------------------------------------------------------------
static
size_t
reference_element_p_first_inod_by_variant (
size_t order,
size_t subgeo_variant)
{
if (subgeo_variant == reference_element__p) return 0;
return 1;
}
static
size_t
reference_element_e_first_inod_by_variant (
size_t order,
size_t subgeo_variant)
{
if (subgeo_variant == reference_element__p) return 0;
if (subgeo_variant == reference_element__e) return 2;
return order+1;
}
static
size_t
reference_element_t_first_inod_by_variant (
size_t order,
size_t subgeo_variant)
{
if (subgeo_variant == reference_element__p) return 0;
if (subgeo_variant == reference_element__e) return 3;
if (subgeo_variant == reference_element__t) return 3 + 3*(order-1);
return (order+1)*(order+2)/2;
}
static
size_t
reference_element_q_first_inod_by_variant (
size_t order,
size_t subgeo_variant)
{
if (subgeo_variant == reference_element__p) return 0;
if (subgeo_variant == reference_element__e) return 4;
if (subgeo_variant == reference_element__t) return 4 + 4*(order-1);
if (subgeo_variant == reference_element__q) return 4 + 4*(order-1);
return (order+1)*(order+1);
}
static
size_t
reference_element_T_first_inod_by_variant (
size_t order,
size_t subgeo_variant)
{
if (subgeo_variant == reference_element__p) return 0;
if (subgeo_variant == reference_element__e) return 4;
if (subgeo_variant == reference_element__t) return 4 + 6*(order-1);
if (subgeo_variant == reference_element__q) return 4 + 6*(order-1) + 4*((order-1)*(order-2)/2);
if (subgeo_variant == reference_element__T) return 4 + 6*(order-1) + 4*((order-1)*(order-2)/2);
return (order+1)*(order+2)*(order+3)/6;
}
static
size_t
reference_element_P_first_inod_by_variant (
size_t order,
size_t subgeo_variant)
{
if (subgeo_variant == reference_element__p) return 0;
if (subgeo_variant == reference_element__e) return 6;
if (subgeo_variant == reference_element__t) return 6 + 9*(order-1);
if (subgeo_variant == reference_element__q) return 6 + 9*(order-1) + 2*((order-1)*(order-2)/2);
if (subgeo_variant == reference_element__T) return 6 + 9*(order-1) + 2*((order-1)*(order-2)/2) + 3*(order-1)*(order-1);
if (subgeo_variant == reference_element__P) return 6 + 9*(order-1) + 2*((order-1)*(order-2)/2) + 3*(order-1)*(order-1);
return ((order+1)*(order+1)*(order+2))/2;
}
static
size_t
reference_element_H_first_inod_by_variant (
size_t order,
size_t subgeo_variant)
{
if (subgeo_variant == reference_element__p) return 0;
if (subgeo_variant == reference_element__e) return 8;
if (subgeo_variant == reference_element__t) return 8 + 12*(order-1);
if (subgeo_variant == reference_element__q) return 8 + 12*(order-1);
if (subgeo_variant == reference_element__T) return 8 + 12*(order-1) + 6*(order-1)*(order-1);
if (subgeo_variant == reference_element__P) return 8 + 12*(order-1) + 6*(order-1)*(order-1);
if (subgeo_variant == reference_element__H) return 8 + 12*(order-1) + 6*(order-1)*(order-1);
return (order+1)*(order+1)*(order+1);
}
static
size_t
reference_element_first_inod_by_variant (
size_t variant,
size_t order,
size_t subgeo_variant)
{
#define _RHEOLEF_geo_element_auto_case(NAME) \
case reference_element__##NAME: \
return reference_element_##NAME##_first_inod_by_variant (order, subgeo_variant);
switch (variant) {
_RHEOLEF_geo_element_auto_case(p)
_RHEOLEF_geo_element_auto_case(e)
_RHEOLEF_geo_element_auto_case(t)
_RHEOLEF_geo_element_auto_case(q)
_RHEOLEF_geo_element_auto_case(T)
_RHEOLEF_geo_element_auto_case(P)
_RHEOLEF_geo_element_auto_case(H)
default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
}
#undef _RHEOLEF_geo_element_auto_case
}
static
size_t
reference_element_last_inod_by_variant (
size_t variant,
size_t order,
size_t subgeo_variant)
{
return reference_element_first_inod_by_variant (variant, order, subgeo_variant+1);
}
// -----------------------------------------------------------------------------
// reference_element::ilat2loc_inod
// -----------------------------------------------------------------------------
static
size_t
reference_element_p_ilat2loc_inod (size_t order, const point_basic<size_t>& ilat)
{
return 0;
}
static
size_t
reference_element_e_ilat2loc_inod (size_t order, const point_basic<size_t>& ilat)
{
size_t i = ilat[0];
if (i == 0) return 0;
if (i == order) return 1;
return 2 + (i - 1);
}
static
size_t
reference_element_t_ilat2loc_inod (size_t order, const point_basic<size_t>& ilat)
{
size_t i = ilat[0];
size_t j = ilat[1];
if (j == 0) { // horizontal edge [0,1]
if (i == 0) return 0;
if (i == order) return 1;
return 3 + 0*(order-1) + (i - 1);
}
if (i + j == order) { // oblique edge ]1,2]
if (j == order) return 2;
return 3 + 1*(order-1) + (j - 1);
}
size_t j1 = order - j;
if (i == 0) { // vertical edge ]2,0[
return 3 + 2*(order-1) + (j1 - 1);
}
// internal face ]0,1[x]0,2[: i, then j
size_t n_face_node = (order-1)*(order-2)/2;
size_t ij = (n_face_node - (j1-1)*j1/2) + (i - 1);
return 3 + 3*(order-1) + ij;
}
static
size_t
reference_element_q_ilat2loc_inod (size_t order, const point_basic<size_t>& ilat)
{
size_t i = ilat[0];
size_t j = ilat[1];
if (j == 0) { // horizontal edge [0,1]
if (i == 0) return 0;
if (i == order) return 1;
return 4 + 0*(order-1) + (i - 1);
}
if (i == order) { // vertical edge ]1,2]
if (j == order) return 2;
return 4 + 1*(order-1) + (j - 1);
}
size_t i1 = order - i;
if (j == order) { // horizontal edge ]2,3]
if (i == 0) return 3;
return 4 + 2*(order-1) + (i1 - 1);
}
size_t j1 = order - j;
if (i == 0) { // vertical edge ]3,0[
return 4 + 3*(order-1) + (j1 - 1);
}
// internal face ]0,1[x]0,3[: i, then j
size_t ij = (order-1)*(j-1) + (i-1);
size_t n_bdry_node = 4 + 4*(order-1);
return n_bdry_node + ij;
}
static
size_t
reference_element_T_ilat2loc_inod (size_t order, const point_basic<size_t>& ilat)
{
size_t i = ilat[0];
size_t j = ilat[1];
size_t k = ilat[2];
size_t n_tot_edge_node = 6*(order-1);
size_t n_face_node = (order-1)*(order-2)/2;
if (k == 0) {
// horizontal face 021 = [0,2]x[0,1]
if (j == 0) { // left edge [0,1]
if (i == 0) return 0;
if (i == order) return 1;
return 4 + 0*(order-1) + (i - 1);
}
if (i + j == order) { // oblique edge ]1,2[
if (j == order) return 2;
return 4 + 1*(order-1) + (j - 1);
}
size_t j1 = order - j;
if (i == 0) { // right edge ]0,2[
return 4 + 2*(order-1) + (j1 - 1);
}
// internal face 021 = ]0,2[x]0,1[: j, then i
size_t i1 = order - i;
size_t ji = (n_face_node - i1*(i1-1)/2) + (j - 1);
return 4 + n_tot_edge_node + 0*n_face_node + ji;
}
if (i == 0) {
// back face 032 =]0,3]x]0,2[
if (j == 0) { // vertical edge ]0,3]
if (k == order) return 3;
return 4 + 3*(order-1) + (k - 1);
}
if (j + k == order) { // oblique edge ]2,3[
return 4 + 5*(order-1) + (k - 1);
}
// internal face 032 =]0,3]x]0,2[, k, then j
size_t j1 = order - j;
size_t kj = (n_face_node - j1*(j1-1)/2) + (k - 1);
return 4 + n_tot_edge_node + 1*n_face_node + kj;
}
if (j == 0) {
// left face 013 = ]0,1[x]0,3[
if (i + k == order) { // oblique edge ]1,3[
return 4 + 4*(order-1) + (k - 1);
}
// internal face 013 = ]0,1[x]0,3[: i, then k
size_t k1 = order - k;
size_t ik = (n_face_node - k1*(k1-1)/2) + (i - 1);
return 4 + n_tot_edge_node + 2*n_face_node + ik;
}
if (i + j + k == order) {
// internal to oblique face 123 = ]1,2[x]1,3], j, then k
size_t k1 = order - k;
size_t jk = (n_face_node - k1*(k1-1)/2) + (j - 1);
return 4 + n_tot_edge_node + 3*n_face_node + jk;
}
// internal to volume: i, then j, then k
size_t n_tot_face_node = 4*n_face_node;
size_t n_tot_node = (order+1)*(order+2)*(order+3)/6;
size_t n_volume_node = (order-1)*(order-2)*(order-3)/6;
size_t k1 = order - k;
size_t j1 = order - j - k;
size_t n_level_k_node = (k1-1)*(k1-2)/2;
size_t ijk = (n_volume_node - k1*(k1-1)*(k1-2)/6)
+ (n_level_k_node - j1*(j1-1)/2)
+ (i - 1);
return 4 + n_tot_edge_node + n_tot_face_node + ijk;
}
static
size_t
reference_element_P_ilat2loc_inod (size_t order, const point_basic<size_t>& ilat)
{
size_t i = ilat[0];
size_t j = ilat[1];
size_t k = ilat[2];
size_t i1 = order - i;
size_t j1 = order - j;
size_t k1 = order - k;
size_t n_node_edge = order-1;
size_t n_tot_edge_node = 9*(order-1);
size_t n_node_face_tri = (order-1)*(order-2)/2;
size_t n_node_face_qua = (order-1)*(order-1);
if (k == 0) {
// horizontal face 021 = [0,1]x[0,2]
if (j == 0) { // left edge [0,1]
if (i == 0) return 0;
if (i == order) return 1;
return 6 + 0*(order-1) + (i - 1);
}
if (i + j == order) { // oblique edge ]1,2[
if (j == order) return 2;
return 6 + 1*(order-1) + (j - 1);
}
if (i == 0) { // right edge ]0,2[
return 6 + 2*(order-1) + (j1 - 1);
}
// internal face 021 = ]0,2[x]0,1[: j, then i
size_t ji = (n_node_face_tri - i1*(i1-1)/2) + (j - 1);
return 6 + n_tot_edge_node + 0*n_node_face_tri + ji;
}
if (k == order) {
// horizontal face 354 = [3,4]x[3,5]
if (j == 0) { // left edge [3,4]
if (i == 0) return 3;
if (i == order) return 4;
return 6 + 6*(order-1) + (i - 1);
}
if (i + j == order) { // oblique edge ]4,5[
if (j == order) return 5;
return 6 + 7*(order-1) + (j - 1);
}
if (i == 0) { // right edge ]3,5[
return 6 + 8*(order-1) + (j1 - 1);
}
// internal face 354 = ]3,4[x]3,5[: i, then j
size_t ij = (n_node_face_tri - j1*(j1-1)/2) + (i - 1);
return 6 + n_tot_edge_node + 1*n_node_face_tri + ij;
}
if (j == 0) {
// left face 0143 = [0,1]x]0,3[
if (i == 0) {
// vertical edge ]0,3[
return 6 + 3*(order-1) + (k - 1);
}
if (i == order) {
// vertical edge ]1,4[
return 6 + 4*(order-1) + (k - 1);
}
// internal face 01243 = ]0,1[x]0,3[: i then k
size_t ik = n_node_edge*(k-1) + (i-1);
return 6 + n_tot_edge_node + 2*n_node_face_tri + 0*n_node_face_qua + ik;
}
if (i == 0) {
// back face 0352 = ]0,2]x]0,3[
if (j == order) {
// vertical edge ]2,5[
return 6 + 5*(order-1) + (k - 1);
}
// internal face 0352 = ]0,2[x]0,3[: k then j
size_t kj = n_node_edge*(j-1) + (k-1);
return 6 + n_tot_edge_node + 2*n_node_face_tri + 2*n_node_face_qua + kj;
}
if (i + j == order) {
// internal face 1254 = ]1,2[x]1,4[: j then k
size_t jk = n_node_edge*(k-1) + (j-1);
return 6 + n_tot_edge_node + 2*n_node_face_tri + 1*n_node_face_qua + jk;
}
// internal volume ]0,1[x]0,2[x][0,3[: i, then j, then k
size_t n_tot_face_node = 2*n_node_face_tri + 3*n_node_face_qua;
size_t ij = (n_node_face_tri - (j1-1)*j1/2) + (i - 1);
size_t ijk = n_node_face_tri*(k-1) + ij;
return 6 + n_tot_edge_node + n_tot_face_node + ijk;
}
static
size_t
reference_element_H_ilat2loc_inod (size_t order, const point_basic<size_t>& ilat)
{
size_t i = ilat[0];
size_t j = ilat[1];
size_t k = ilat[2];
size_t i1 = order - i;
size_t j1 = order - j;
size_t k1 = order - k;
size_t n_node_edge = order-1;
size_t n_node_face = (order-1)*(order-1);
if (k == 0) { // bottom face [0,3]x[0,1]
if (j == 0) { // bottom-left edge [0,1]
if (i == 0) return 0;
if (i == order) return 1;
return 8 + 0*n_node_edge + (i - 1);
}
if (j == order) { // bottom-right edge [3,2]
if (i == 0) return 3;
if (i == order) return 2;
return 8 + 2*n_node_edge + (i1 - 1);
}
if (i == 0) { // bottom-back edge ]3,0[
return 8 + 3*n_node_edge + (j1 - 1);
}
if (i == order) { // bottom-front edge ]1,2[
return 8 + 1*n_node_edge + (j - 1);
}
// internal face ]0,3[x]0,1[: j then i (TODO: not checked since Pk(H) only with k <= 2 yet
size_t ji = n_node_edge*(i-1) + (j-1);
return 8 + 12*n_node_edge + 0*n_node_face + ji;
}
if (k == order) { // top face [4,5]x[4,7]
if (j == 0) { // top-left edge [4,5]
if (i == 0) return 4;
if (i == order) return 5;
return 8 + 8*n_node_edge + (i - 1);
}
if (j == order) { // top-right edge [7,6]
if (i == 0) return 7;
if (i == order) return 6;
return 8 + 10*n_node_edge + (i1 - 1);
}
if (i == 0) { // top-back edge ]7,4[
return 8 + 11*n_node_edge + (j1 - 1);
}
if (i == order) { // top-front edge ]5,6[
return 8 + 9*n_node_edge + (j - 1);
}
// internal face ]4,5[x]4,7[: i then j
size_t ij = n_node_edge*(j-1) + (i-1);
return 8 + 12*n_node_edge + 3*n_node_face + ij;
}
if (j == 0) { // left face ]0,1[x[0,4]
if (i == 0) { // left-back edge [0,4]
return 8 + 4*n_node_edge + (k - 1);
}
if (i == order) { // left-front edge [1,5]
return 8 + 5*n_node_edge + (k - 1);
}
// internal face ]0,1[x]0,4[: i then k
size_t ik = n_node_edge*(k-1) + (i-1);
return 8 + 12*n_node_edge + 2*n_node_face + ik;
}
if (j == order) { // right face ]2,3[x[2,6]
if (i == 0) { // right-back edge [3,7]
return 8 + 7*n_node_edge + (k - 1);
}
if (i == order) { // right-front edge [2,6]
return 8 + 6*n_node_edge + (k - 1);
}
// internal face ]2,3[x]2,6[: i1 then k
size_t i1k = n_node_edge*(k-1) + (i1-1);
return 8 + 12*n_node_edge + 5*n_node_face + i1k;
}
if (i == 0) { // internal back face ]0,4[x[0,3]: k then j
size_t kj = n_node_edge*(j-1) + (k-1);
return 8 + 12*n_node_edge + 1*n_node_face + kj;
}
if (i == order) { // internal front face ]1,2[x[1,5]: j then k
size_t jk = n_node_edge*(k-1) + (j-1);
return 8 + 12*n_node_edge + 4*n_node_face + jk;
}
// internal volume: i then j then k
size_t n_node_bdry = 8 + 12*n_node_edge + 6*n_node_face;
size_t ijk = (order-1)*((order-1)*(k-1) + (j-1)) + (i-1);
return n_node_bdry + ijk;
}
} // namespace rheolef
#endif // _RHEOLEF_REFERENCE_ELEMENT_AUX_ICC
rheolef::point_basic
Definition point.h:86
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::max_variant
static const variant_type max_variant
Definition reference_element.h:82
rheolef::reference_element::p
static const variant_type p
Definition reference_element.h:75
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::reference_element__H
constexpr size_t reference_element__H
Definition reference_element_aux.icc:39
rheolef::reference_element__e
constexpr size_t reference_element__e
Definition reference_element_aux.icc:34
rheolef::reference_element__P
constexpr size_t reference_element__P
Definition reference_element_aux.icc:38
rheolef::_first_variant_by_dimension
static const size_t _first_variant_by_dimension[5]
Definition reference_element_aux.icc:151
rheolef::_reference_element_n_face_by_variant
static size_t _reference_element_n_face_by_variant[reference_element__max_variant]
Definition reference_element_aux.icc:54
rheolef::reference_element__T
constexpr size_t reference_element__T
Definition reference_element_aux.icc:37
rheolef::reference_element__p
constexpr size_t reference_element__p
Definition reference_element_aux.icc:33
rheolef::reference_element__t
constexpr size_t reference_element__t
Definition reference_element_aux.icc:35
rheolef::reference_element__q
constexpr size_t reference_element__q
Definition reference_element_aux.icc:36
rheolef::reference_element__max_variant
constexpr size_t reference_element__max_variant
Definition reference_element_aux.icc:40
rheolef::_reference_element_n_edge_by_variant
static size_t _reference_element_n_edge_by_variant[reference_element__max_variant]
Definition reference_element_aux.icc:53
rheolef::_reference_element_dimension_by_variant
static size_t _reference_element_dimension_by_variant[reference_element__max_variant]
Definition reference_element_aux.icc:52
error_macro
#define error_macro(message)
Definition dis_macros.h:49
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.dim
int dim
Definition mkgeo_ball.sh:307
mkgeo_ball.variant
int variant
Definition mkgeo_ball.sh:149
mkgeo_ball.order
int order
Definition mkgeo_ball.sh:343
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
_RHEOLEF_geo_element_auto_case
#define _RHEOLEF_geo_element_auto_case(VARIANT)