reference_element.cc

Go to the documentation of this file.

 
#include "reference_element.h"
#include "rheolef/point.h"
 
#include "reference_element_aux.icc"
 
namespace rheolef {
 
namespace edge {
#include "edge.icc"
} // namespace edge
 
namespace triangle {
#include "triangle.icc"
} // namespace triangle
 
namespace quadrangle {
#include "quadrangle.icc"
} // namespace quadrangle
 
namespace tetrahedron {
#include "tetrahedron.icc"
} // namespace tetrahedron
 
namespace prism {
#include "prism.icc"
} // namespace prism
 
namespace hexahedron {
#include "hexahedron.icc"
} // namespace hexahedron
 
#include "reference_element_declare.cc"
 
// ==============================================================================
// generic
// ==============================================================================
void
reference_element::set_name (char name)
{
  _x = variant(name);
  check_macro (_x != max_variant, "undefined reference element `" << name << "'");
}
Float
measure (reference_element hat_K) {
    return hat_K_measure [hat_K.variant()];
}
reference_element::variant_type
reference_element::variant (char name)
{
  return reference_element_variant (_name, name);
} 
reference_element::variant_type
reference_element::variant (size_type n_vertex, size_type dim)
{
  size_type variant = 0;
  for (; variant < max_variant; variant++) {
    if (_dimension[variant] == dim && _n_vertex[variant] == n_vertex) {
      return variant_type(variant);
    }
  }
  error_macro ("undefined "<<dim<<"d reference element with "<<n_vertex << " vertices");
}
reference_element::size_type
reference_element::n_sub_edge (variant_type variant)
{
  switch (variant) {
    case reference_element::p: return 0;
    case reference_element::e: return 0;
    case reference_element::t: return 3;
    case reference_element::q: return 4;
    case reference_element::T: return 6;
    case reference_element::P: return 9;
    case reference_element::H: return 12;
    default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
  }
}
reference_element::size_type
reference_element::n_sub_face (variant_type variant)
{
  switch (variant) {
    case reference_element::p: return 0;
    case reference_element::e: return 0;
    case reference_element::t: return 0;
    case reference_element::q: return 0;
    case reference_element::T: return 4;
    case reference_element::P: return 5;
    case reference_element::H: return 6;
    default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
  }
}
reference_element::size_type
reference_element::n_subgeo (variant_type variant, size_type subgeo_dim)
{
#define _RHEOLEF_reference_element_case(VARIANT)                         \
    case   reference_element::VARIANT:                                \
      return reference_element_##VARIANT::n_subgeo (subgeo_dim);
 
  switch (variant) {
        _RHEOLEF_reference_element_case(p)
        _RHEOLEF_reference_element_case(e)
        _RHEOLEF_reference_element_case(t)
        _RHEOLEF_reference_element_case(q)
        _RHEOLEF_reference_element_case(T)
        _RHEOLEF_reference_element_case(P)
        _RHEOLEF_reference_element_case(H)
    default: error_macro ("unexpected element variant `"<<variant<<"'"); return 0;
  }
#undef  _RHEOLEF_reference_element_case
}
reference_element::size_type
reference_element::subgeo_n_node (
    variant_type variant,
    size_type    order, 
    size_type    subgeo_dim,
    size_type    loc_isid)
{
#define _RHEOLEF_geo_element_auto_case(VARIANT)                         \
    case   reference_element::VARIANT:                                  \
      return reference_element_##VARIANT::subgeo_n_node (order, subgeo_dim, loc_isid);
 
  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
}
reference_element::size_type
reference_element::subgeo_local_node (
    variant_type variant,
    size_type    order, 
    size_type    subgeo_dim,
    size_type    loc_isid,
    size_type    loc_jsidnod)
{
#define _RHEOLEF_geo_element_auto_case(VARIANT)                         \
    case   reference_element::VARIANT:                                  \
      return reference_element_##VARIANT::subgeo_local_node (order, subgeo_dim, loc_isid, loc_jsidnod);
 
  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
}
reference_element::size_type
reference_element::first_inod_by_variant (
    variant_type variant,
    size_type    order, 
    variant_type subgeo_variant)
{
#define _RHEOLEF_geo_element_auto_case(VARIANT)                         \
    case reference_element::VARIANT:                                    \
      return reference_element_##VARIANT##_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
}
reference_element::size_type
reference_element::n_node (variant_type variant, size_type order)
{
  return reference_element_n_node (variant, order);
}
void
reference_element::init_local_nnode_by_variant (
    size_type                                             order, 
    std::array<size_type,reference_element::max_variant>& sz)
{
  reference_element_init_local_nnode_by_variant (order, sz);
}
Float
reference_element::side_measure (size_type loc_isid) const
{
#define _RHEOLEF_geo_element_auto_case(VARIANT)                         \
    case   reference_element::VARIANT:                                  \
      return reference_element_##VARIANT::side_measure (loc_isid);
 
  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
}
void
reference_element::side_normal (size_type loc_isid, point_basic<Float>& hat_n) const
{
#define _RHEOLEF_geo_element_auto_case(VARIANT)                         \
    case   reference_element::VARIANT:                                  \
      reference_element_##VARIANT::side_normal (loc_isid, hat_n); break;
 
  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()<<"'");
  }
#undef _RHEOLEF_geo_element_auto_case
}
const point_basic<Float>&
reference_element::vertex (size_type iloc) const
{
#define _RHEOLEF_geo_element_auto_case(VARIANT)                         \
    case   reference_element::VARIANT:                                  \
      return reference_element_##VARIANT::vertex (iloc);
 
  static point_basic<Float> dummy;
  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 dummy;
  }
#undef _RHEOLEF_geo_element_auto_case
}
// ==============================================================================
// point
// ==============================================================================
reference_element_p::size_type
reference_element_p::n_subgeo (size_type side_dim)
{
  return (side_dim == 0) ? 1 : 0;
}
reference_element_p::size_type
reference_element_p::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
{
  return (side_dim == 0) ? 1 : 0;
}
reference_element_p::size_type
reference_element_p::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
{
  return 0;
}
// edge 0d-lattice for high order elements, i <= 0
// convert to local inod geo_element number for rheolef
reference_element_p::size_type
reference_element_p::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
{
  return reference_element_p_ilat2loc_inod (order, ilat);
}
reference_element_p::size_type
reference_element_p::first_inod_by_variant (
    size_type    order, 
    size_type    subgeo_variant)
{
  return reference_element_p_first_inod_by_variant (order, subgeo_variant);
}
Float
reference_element_p::side_measure (size_type loc_isid)
{
  return 0;
}
void
reference_element_p::side_normal (size_type loc_isid, point_basic<Float>& hat_n)
{
}
const point_basic<Float>&
reference_element_p::vertex (size_type iloc)
{
  static point_basic<Float> dummy;
  return dummy;
}
// ==============================================================================
// edge
// ==============================================================================
reference_element_e::size_type
reference_element_e::n_subgeo (size_type side_dim)
{
  switch (side_dim) {
    case 0:  return 2;
    case 1:  return 1;
    default: return 0;
  }
}
reference_element_e::size_type
reference_element_e::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
{
  switch (side_dim) {
    case 0:  return 1;
    case 1:  return order+1;
    default: return 0;
  }
}
reference_element_e::size_type
reference_element_e::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
{
  switch (side_dim) {
    case 0:  return loc_isid;
    case 1:  return loc_jsidnod;
    default: return 0;
  }
}
// edge 1d-lattice (i) for high order elements, i <= order
// convert to local inod geo_element number for rheolef
reference_element_e::size_type
reference_element_e::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
{
  return reference_element_e_ilat2loc_inod (order, ilat);
}
reference_element_e::size_type
reference_element_e::first_inod_by_variant (
    size_type    order, 
    size_type    subgeo_variant)
{
  return reference_element_e_first_inod_by_variant (order, subgeo_variant);
}
Float
reference_element_e::side_measure (size_type loc_isid)
{
  return 1;
}
void
reference_element_e::side_normal (size_type loc_isid, point_basic<Float>& hat_n)
{
  hat_n[0] = (loc_isid == 0) ? -1 : 1;
}
const point_basic<Float>&
reference_element_e::vertex (size_type iloc)
{
  return edge::vertex[iloc];
}
// ==============================================================================
// triangle
// ==============================================================================
reference_element_t::size_type
reference_element_t::n_subgeo (size_type side_dim)
{
  switch (side_dim) {
    case 0:  return 3;
    case 1:  return 3;
    case 2:  return 1;
    default: return 0;
  }
}
reference_element_t::size_type
reference_element_t::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
{
  switch (side_dim) {
    case 0:  return 1;
    case 1:  return order+1;
    case 2:  return ((order+1)*(order+2))/2;
    default: return 0;
  }
}
reference_element_t::size_type
reference_element_t::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
{
  switch (side_dim) {
    case 0:  return loc_isid;
    case 1:  if (loc_jsidnod < 2) return (loc_isid + loc_jsidnod) % 3;         // edge-node is a vertex
    else                          return (order-1)*loc_isid + loc_jsidnod + 1; // edge-node is edge-internal
    case 2:  return loc_jsidnod;
    default: return 0;
  }
}
// triangle lattice (i,j) for high order elements, i+j <= order
// convert to local inod geo_element number for rheolef
reference_element_t::size_type
reference_element_t::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
{
  return reference_element_t_ilat2loc_inod (order, ilat);
}
reference_element_t::size_type
reference_element_t::first_inod_by_variant (
    size_type    order, 
    size_type    subgeo_variant)
{
  return reference_element_t_first_inod_by_variant (order, subgeo_variant);
}
Float
reference_element_t::side_measure (size_type loc_isid)
{
  return (loc_isid != 1) ? 1 : sqrt(Float(2.));
}
void
reference_element_t::side_normal (size_type loc_isid, point_basic<Float>& hat_n)
{
  switch (loc_isid) {
    case 0:  hat_n = point_basic<Float>( 0,-1); break;
    case 1:  hat_n = point_basic<Float>( 1, 1)/sqrt(Float(2)); break;
    case 2:
    default: hat_n = point_basic<Float>(-1,0); break;
  }
}
const point_basic<Float>&
reference_element_t::vertex (size_type iloc)
{
  return triangle::vertex[iloc];
}
// ==============================================================================
// quadrangle
// ==============================================================================
reference_element_q::size_type
reference_element_q::n_subgeo (size_type side_dim)
{
  switch (side_dim) {
    case 0:  return 4;
    case 1:  return 4;
    case 2:  return 1;
    default: return 0;
  }
}
reference_element_q::size_type
reference_element_q::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
{
  switch (side_dim) {
    case 0:  return 1;
    case 1:  return order+1;
    case 2:  return (order+1)*(order+1);
    default: return 0;
  }
}
reference_element_q::size_type
reference_element_q::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
{
  switch (side_dim) {
    case 0:  return loc_isid;
    case 1:  if (loc_jsidnod < 2) return (loc_isid + loc_jsidnod) % 4;         // edge-node is a vertex
    else                          return (order-1)*loc_isid + loc_jsidnod + 2; // edge-node is edge-internal
    case 2:  return loc_jsidnod;
    default: return 0;
  }
}
// quadrangle lattice (i,j) for high order elements, 0 <= i,j <= order
// convert to local inod geo_element number for rheolef
reference_element_q::size_type
reference_element_q::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
{
  return reference_element_q_ilat2loc_inod (order, ilat);
}
reference_element_q::size_type
reference_element_q::first_inod_by_variant (
    size_type    order, 
    size_type    subgeo_variant)
{
  return reference_element_q_first_inod_by_variant (order, subgeo_variant);
}
Float
reference_element_q::side_measure (size_type loc_isid)
{
  return 2;
}
void
reference_element_q::side_normal (size_type loc_isid, point_basic<Float>& hat_n) 
{
  switch (loc_isid) {
    case 0:  hat_n = point_basic<Float>( 0,-1); break;
    case 1:  hat_n = point_basic<Float>( 1, 0); break;
    case 2:  hat_n = point_basic<Float>( 0, 1); break;
    case 3:
    default: hat_n = point_basic<Float>(-1, 0); break;
  }
}
const point_basic<Float>&
reference_element_q::vertex (size_type iloc)
{
  return quadrangle::vertex[iloc];
}
// ==============================================================================
// tetra
// ==============================================================================
reference_element_T::size_type
reference_element_T::n_subgeo (size_type side_dim)
{
  switch (side_dim) {
    case 0:  return 4;
    case 1:  return 6;
    case 2:  return 4;
    case 3:  return 1;
    default: return 0;
  }
}
reference_element_T::size_type
reference_element_T::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid) 
{
  switch (side_dim) {
    case 0:  return 1;
    case 1:  return order+1;
    case 2:  return (order+1)*(order+2)/2;
    case 3:  return ((order+1)*(order+2)*(order+3))/6;
    default: return 0;
  }
}
reference_element_T::size_type
reference_element_T::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod) 
{
  switch (side_dim) {
    case 0:  return loc_isid;
    case 1:  if (loc_jsidnod < 2) { // edge-node is a vertex
               if (loc_isid < 3) return (loc_isid + loc_jsidnod) % 3;
               else              return loc_jsidnod == 0 ? loc_isid - 3 : 3;
             } else {               // edge-node is internal to the edge
               return loc_isid*(order-1) + loc_jsidnod + 2;
             }
    case 2:  {
               if (loc_jsidnod < 3) { // face-node is a vertex
                 if (loc_isid == 3)    return loc_jsidnod + 1;
                 if (loc_jsidnod == 0) return 0;
                 if (loc_jsidnod == 2) return loc_isid + 1;
                 return ((loc_isid + 1) % 3) + 1;
               }
               // face-node is edge-internal(P2) or face-internal(P3) ; volume-internal(P4) impossible
               size_type last_edge_node_iloc = 3 + 3*(order-1);
               if (loc_jsidnod < last_edge_node_iloc) { // edge-internal
#ifdef TO_CLEAN
                 extern const size_type geo_element_T_fac2edg_idx    [4][3];
                 extern const int       geo_element_T_fac2edg_orient [4][3];
#endif // TO_CLEAN
                 size_type order1 = order - 1; // avoid div by zero compiler error
                 size_type loc_jedg = (loc_jsidnod-3) / order1;
                 size_type loc_kedg = (loc_jsidnod-3) % order1;
                 if (geo_element_T_fac2edg_orient [loc_isid][loc_jedg] < 0) {
                        loc_kedg = order - loc_kedg - 2;
                 }
                 return (order-1)*geo_element_T_fac2edg_idx [loc_isid][loc_jedg] + loc_kedg + 4;
               }
               // face-node is face-internal(P3)
               size_type ij_loc = (loc_jsidnod - last_edge_node_iloc);
               return 4 + 6*(order-1) + loc_isid*(order-1)*(order-2)/2 + ij_loc;
             }
    case 3:  return loc_jsidnod;
    default: return 0;
  }
}
// tetrahedron lattice (i,j,k) for high order elements, i+j+k <= order
// convert to local inod geo_element number for rheolef
reference_element_T::size_type
reference_element_T::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
{
  return reference_element_T_ilat2loc_inod (order, ilat);
}
reference_element_T::size_type
reference_element_T::first_inod_by_variant (
    size_type    order, 
    size_type    subgeo_variant)
{
  return reference_element_T_first_inod_by_variant (order, subgeo_variant);
}
reference_element_T::size_type
reference_element_T::face2edge (size_type loc_iface, size_type loc_iface_jedg)
{
  return geo_element_T_fac2edg_idx [loc_iface][loc_iface_jedg];
}
int
reference_element_T::face2edge_orient (size_type loc_iface, size_type loc_iface_jedg)
{
  return geo_element_T_fac2edg_orient [loc_iface][loc_iface_jedg];
}
Float
reference_element_T::side_measure (size_type loc_isid)
{
  return (loc_isid != 3) ? 0.5 : sqrt(Float(3.))/2.;
}
void
reference_element_T::side_normal (size_type loc_isid, point_basic<Float>& hat_n)
{
  // TODO: generate side_normal automatically from tetrahedron.icc, as it depends upon conventions
  switch (loc_isid) {
    case 0:  hat_n = point_basic<Float>( 0, 0,-1); break;
    case 1:  hat_n = point_basic<Float>(-1, 0, 0); break;
    case 2:  hat_n = point_basic<Float>( 0,-1, 0); break;
    case 3:
    default: hat_n = point_basic<Float>( 1, 1, 1)/sqrt(Float(3)); break;
  }
}
const point_basic<Float>&
reference_element_T::vertex (size_type iloc)
{
  return tetrahedron::vertex[iloc];
}
// ==============================================================================
// prism
// ==============================================================================
reference_element_P::size_type
reference_element_P::n_subgeo (size_type side_dim)
{
  switch (side_dim) {
    case 0:  return 6;
    case 1:  return 9;
    case 2:  return 5;
    case 3:  return 1;
    default: return 0;
  }
}
reference_element_P::size_type
reference_element_P::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
{
  switch (side_dim) {
    case 0:  return 1;
    case 1:  return order+1;
    case 2:  return (loc_isid < 2) ? (order+1)*(order+2)/2 : (order+1)*(order+1);
    case 3:  return ((order+1)*(order+1)*(order+2))/2;
    default: return 0;
  }
}
reference_element_P::size_type
reference_element_P::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
{
  switch (side_dim) {
    case 0:  return loc_isid;
    case 1:  if (loc_jsidnod < 2) { // edge-node is a vertex
               return prism::edge [loc_isid] [loc_jsidnod];
             } else {               // edge-node is internal to the edge
               return 6 + loc_isid*(order-1) + (loc_jsidnod - 2);
             }
    case 2:  {
               size_type n_vert_on_side = (loc_isid < 2) ? 3 : 4;
               if (loc_jsidnod < n_vert_on_side) { // face-node is a vertex
                return prism::face [loc_isid] [loc_jsidnod];
               }
               // face-node is edge-internal(P2) or face-internal(P3) ; volume-internal is impossible
               size_type last_edge_node_iloc = n_vert_on_side + n_vert_on_side*(order-1);
               if (loc_jsidnod < last_edge_node_iloc) { // edge-internal
#ifdef TO_CLEAN
                 extern const size_type geo_element_P_fac2edg_idx    [5][4];
                 extern const int       geo_element_P_fac2edg_orient [5][4];
#endif // TO_CLEAN
                 size_type order1 = order - 1; // avoid div by zero compiler error
                 size_type loc_jedg = (loc_jsidnod - n_vert_on_side) / order1;
                 size_type loc_kedg = (loc_jsidnod - n_vert_on_side) % order1;
                 if (geo_element_P_fac2edg_orient [loc_isid][loc_jedg] < 0) {
                   loc_kedg = order - loc_kedg - 2;
                 }
                 return 6 + (order-1)*geo_element_P_fac2edg_idx [loc_isid][loc_jedg] + loc_kedg;
               }
               // face-node is face-internal(P3)
               size_type ij_loc = (loc_jsidnod - last_edge_node_iloc);
               size_type shift_prev_faces;
               if (loc_isid < 2) {
                 shift_prev_faces = loc_isid*(order-1)*(order-2)/2;
               } else {
                 shift_prev_faces = 2*(order-1)*(order-2)/2 + (loc_isid-2)*(order-1)*(order-1);
               }
               return 6 + 9*(order-1) + shift_prev_faces + ij_loc;
             }
    case 3:  return loc_jsidnod;
    default: return 0;
  }
}
reference_element_P::size_type
reference_element_P::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
{
  return reference_element_P_ilat2loc_inod (order, ilat);
}
reference_element_P::size_type
reference_element_P::first_inod_by_variant (
    size_type    order, 
    size_type    subgeo_variant)
{
  return reference_element_P_first_inod_by_variant (order, subgeo_variant);
}
Float
reference_element_P::side_measure (size_type loc_isid)
{
  if (loc_isid <= 1) return 0.5;              // top and bottom faces
  if (loc_isid == 3) return 2*sqrt(Float(2)); // largest vertical face
  return 2;                                   // two others vertical faces
}
void
reference_element_P::side_normal (size_type loc_isid, point_basic<Float>& hat_n)
{
  fatal_macro ("side_normal: not yet on prism");
}
const point_basic<Float>&
reference_element_P::vertex (size_type iloc)
{
  return prism::vertex[iloc];
}
// ==============================================================================
// hexa
// ==============================================================================
 
reference_element_H::size_type
reference_element_H::n_subgeo (size_type side_dim)
{
  switch (side_dim) {
    case 0:  return 8;
    case 1:  return 12;
    case 2:  return 6;
    case 3:  return 1;
    default: return 0;
  }
}
reference_element_H::size_type
reference_element_H::subgeo_n_node (size_type order, size_type side_dim, size_type loc_isid)
{
  switch (side_dim) {
    case 0:  return 1;
    case 1:  return order+1;
    case 2:  return (order+1)*(order+1);
    case 3:  return (order+1)*(order+1)*(order+1);
    default: return 0;
  }
}
reference_element_H::size_type
reference_element_H::subgeo_local_node (size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
{
  switch (side_dim) {
    case 0:  return loc_isid;
    case 1:  if (loc_jsidnod < 2) { // edge-node is a vertex
                    if (loc_isid < 4) return (loc_isid + loc_jsidnod) % 4;
               else if (loc_isid < 8) return loc_isid - 4 + loc_jsidnod*4;
               else                   return loc_jsidnod == 0 ? loc_isid - 4 : (loc_isid - 7)%4 + 4;
             } else {               // edge-node is internal to the edge
               return loc_isid*(order-1) + loc_jsidnod + 6; // TODO : BUG? +8 a la place de +6 ??
             }
    case 2:  {
               if (loc_jsidnod < 4) { // face-node is a vertex
                      return hexahedron::face [loc_isid][loc_jsidnod];
               } 
               // face-node is edge-internal(P2) or face-internal(P2) ; volume-internal(P2) impossible
               size_type last_edge_node_iloc = 4 + 4*(order-1);
               if (loc_jsidnod < last_edge_node_iloc) { // edge-internal
#ifdef TO_CLEAN
                 extern const size_t geo_element_H_fac2edg_idx    [6][4];
                 extern const int    geo_element_H_fac2edg_orient [6][4];
#endif // TO_CLEAN
                 size_type order1 = order - 1; // avoid div by zero compiler error
                 size_type loc_jedg = (loc_jsidnod-4) / order1;
                 size_type loc_kedg = (loc_jsidnod-4) % order1;
                 if (geo_element_H_fac2edg_orient [loc_isid][loc_jedg] < 0) {
                   loc_kedg = order - loc_kedg - 2;
                 }
                 return (order-1)*geo_element_H_fac2edg_idx [loc_isid][loc_jedg] + loc_kedg + 8;
               }
               // face-node is face-internal(P2)
               return 8 + 12*(order-1) + (order-1)*(order-1)*loc_isid + (loc_jsidnod - last_edge_node_iloc);
             }
    case 3:  return loc_jsidnod;
    default: return 0;
  }
}
// edge 0d-lattice for high order elements, i <= 0
// convert to local inod geo_element number for rheolef
reference_element_H::size_type
reference_element_H::ilat2loc_inod (size_type order, const point_basic<size_type>& ilat)
{
  return reference_element_H_ilat2loc_inod (order, ilat);
}
reference_element_H::size_type
reference_element_H::first_inod_by_variant (
    size_type    order, 
    size_type    subgeo_variant)
{
  return reference_element_H_first_inod_by_variant (order, subgeo_variant);
}
Float
reference_element_H::side_measure (size_type loc_isid)
{
  return 2;
}
void
reference_element_H::side_normal (size_type loc_isid, point_basic<Float>& hat_n)
{
  switch (loc_isid) {
    case 0:  hat_n = point_basic<Float>( 0, 0,-1); break;
    case 1:  hat_n = point_basic<Float>(-1, 0, 0); break;
    case 2:  hat_n = point_basic<Float>( 0,-1, 0); break;
    case 3:  hat_n = point_basic<Float>( 0, 0, 1); break;
    case 4:  hat_n = point_basic<Float>( 1, 0, 0); break;
    case 5: 
    default: hat_n = point_basic<Float>( 0, 1, 0); break;
  }
}
const point_basic<Float>&
reference_element_H::vertex (size_type iloc)
{
  return hexahedron::vertex[iloc];
}
// ==============================================================================
 
} // namespace rheolef