Loading...
Searching...
No Matches
geo_element_aux.icc
Go to the documentation of this file.
#ifndef _RHEOLEF_GEO_ELEMENT_V4_AUX_ICC
#define _RHEOLEF_GEO_ELEMENT_V4_AUX_ICC
// shared code, used in:
// - geo_element_v4.cc
// - msh2geo.cc
// in order to avoid code duplication
//
namespace rheolef {
// =============================================================================
// compute orientation & shift
// =============================================================================
// for a triangular side (dis_iv0,dis_iv1,dis_iv2)
// assume that current geo_element is a triangle that contains the same vertices
template <typename GeoElement, typename OrientationType, typename ShiftType>
static
void
geo_element_get_orientation_and_shift (
const GeoElement& S,
size_t dis_iv0, size_t dis_iv1, size_t dis_iv2,
OrientationType& orient,
ShiftType& shift)
{
if (S[0] == dis_iv0) {
orient = (S[1] == dis_iv1) ? 1 : -1;
shift = 0;
} else if (S[0] == dis_iv1) {
orient = (S[1] == dis_iv2) ? 1 : -1;
shift = 1;
} else {
orient = (S[1] == dis_iv0) ? 1 : -1;
shift = 2;
}
}
// for a quadrangular side (dis_iv0,dis_iv1,dis_iv2,dis_iv3)
// assume that current geo_element is a quadrangular that contains the same vertices
template <typename GeoElement, typename OrientationType, typename ShiftType>
static
void
geo_element_get_orientation_and_shift (
const GeoElement& S,
size_t dis_iv0, size_t dis_iv1, size_t dis_iv2, size_t dis_iv3,
OrientationType& orient,
ShiftType& shift)
{
if (S[0] == dis_iv0) {
orient = (S[1] == dis_iv1) ? 1 : -1;
shift = 0;
} else if (S[0] == dis_iv1) {
orient = (S[1] == dis_iv2) ? 1 : -1;
shift = 1;
} else if (S[0] == dis_iv2) {
orient = (S[1] == dis_iv3) ? 1 : -1;
shift = 2;
} else {
orient = (S[1] == dis_iv0) ? 1 : -1;
shift = 3;
}
}
// =============================================================================
// fix rotation and orientation on a 2d edge or 3d face
// =============================================================================
// --------------
// 1) edge
// --------------
static
size_t
geo_element_fix_edge_indirect (
const geo_element& K,
size_t loc_iedg,
size_t loc_iedg_j,
size_t order)
{
if (K.dimension() < 2 || order <= 2 || K.edge_indirect(loc_iedg).orientation() > 0) return loc_iedg_j;
return order - 2 - loc_iedg_j;
}
// --------------
// 2) triangles
// --------------
static
void
geo_element_loc_tri_inod2lattice (
size_t loc_tri_inod,
size_t order,
point_basic<size_t>& ij_lattice)
{
// search (i,j) integers such that
// loc_tri_inod = (n_face_node - (j1-1)*j1/2) + i-1, where j1=order-1-j
// first, search j1 integer such that for i=1:
// i - 1 = 0 = loc_tri_inod - (n_face_node - (j1-1)*j1/2);
// <=> j1^2 - j1 - 2*(n_face_node - loc_tri_inod) = 0
size_t n_face_node = (order-1)*(order-2)/2;
double a = 1;
double b = -1;
double c = -2*int(n_face_node - loc_tri_inod);
double delta = b*b - 4*a*c;
check_macro (delta >= 0, "loc_tri_inod2lattice(loc_tri_inod="<<loc_tri_inod<<",order="<<order<<"): unexpected delta="<<delta<<" < 0");
double j1_plus = (-b + sqrt(delta))/(2*a);
// j1_minus = (-b - sqrt(delta))/(2*a); is negative always
size_t j = floor (order - j1_plus);
size_t j1 = order - j;
size_t i = loc_tri_inod - (n_face_node - (j1-1)*j1/2) + 1 ;
ij_lattice = point_basic<size_t>(i,j);
}
static
size_t
geo_element_fix_triangle_indirect (
geo_element_indirect::orientation_type orient,
geo_element_indirect::shift_type shift,
size_t order,
size_t loc_itri_j)
{
// 1) compute lattice coord (i,j) on this face, from current loc_itri_j
point_basic<size_t> ij_lattice;
geo_element_loc_tri_inod2lattice (loc_itri_j, order, ij_lattice);
size_t i = ij_lattice[0];
size_t j = ij_lattice[1];
// 2) then shift and/or inverse-orient (i,j)
size_t coord[3];
coord [0] = i;
coord [1] = j;
coord [2] = order - i - j;
size_t i_fix = coord [shift%3];
size_t j_fix = coord [(shift+1)%3];
if (orient < 0) std::swap (i_fix, j_fix);
// 3) re-compute the fixed loc_itri_j face index
size_t loc_tri_nnod = (order-1)*(order-2)/2;
size_t j1_fix = order - j_fix;
size_t loc_itri_j_fix = (loc_tri_nnod - (j1_fix-1)*j1_fix/2) + (i_fix-1);
return loc_itri_j_fix;
}
static
size_t
geo_element_fix_triangle_indirect (
const geo_element& K,
size_t loc_ifac,
size_t loc_itri_j,
size_t order)
{
// shift and/or inverse-orient the lattice(i,j)
if (K.dimension() < 3 ||
(K.face_indirect(loc_ifac).orientation() > 0 && K.face_indirect(loc_ifac).shift() == 0)) {
return loc_itri_j;
}
return geo_element_fix_triangle_indirect (
K.face_indirect(loc_ifac).orientation(),
K.face_indirect(loc_ifac).shift(),
order,
loc_itri_j);
}
// --------------
// 3) quadrangles
// --------------
static
void
geo_element_loc_qua_inod2lattice (
size_t loc_qua_inod,
size_t order,
point_basic<size_t>& ij_lattice)
{
// search (i,j) integers in [0:order-1[ such that
// loc_qua_inod = (order-1)*(j-1) + (i-1);
ij_lattice[0] = (loc_qua_inod % (order-1)) + 1;
ij_lattice[1] = (loc_qua_inod / (order-1)) + 1;
}
static
size_t
geo_element_fix_quadrangle_indirect (
geo_element_indirect::orientation_type orient,
geo_element_indirect::shift_type shift,
size_t order,
size_t loc_iqua_j)
{
// 1) compute lattice coord (i,j) on this face, from current loc_itri_j
point_basic<size_t> ij_lattice;
geo_element_loc_qua_inod2lattice (loc_iqua_j, order, ij_lattice);
size_t i = ij_lattice[0];
size_t j = ij_lattice[1];
// 2) then shift and/or inverse-orient (i,j)
size_t coord [4];
coord [0] = i;
coord [1] = j;
coord [2] = order - i;
coord [3] = order - j;
size_t i_fix = coord [shift%4];
size_t j_fix = coord [(shift+1)%4];
if (orient < 0) std::swap (i_fix, j_fix);
// 3) re-compute the fixed loc_iqua_j face index
size_t loc_iqua_j_fix = (order-1)*(j_fix-1) + (i_fix-1);
return loc_iqua_j_fix;
}
static
size_t
geo_element_fix_quadrangle_indirect (
const geo_element& K,
size_t loc_ifac,
size_t loc_iqua_j,
size_t order)
{
// shift and/or inverse-orient the lattice(i,j)
if (K.dimension() < 3 || order <= 2 ||
(K.face_indirect(loc_ifac).orientation() > 0 && K.face_indirect(loc_ifac).shift() == 0)) {
return loc_iqua_j;
}
return geo_element_fix_quadrangle_indirect (
K.face_indirect(loc_ifac).orientation(),
K.face_indirect(loc_ifac).shift(),
order,
loc_iqua_j);
}
} // namespace rheolef
#endif // _RHEOLEF_GEO_ELEMENT_V4_AUX_ICC
rheolef::geo_element_indirect::shift_type
short int shift_type
Definition geo_element_indirect.h:39
rheolef::geo_element_indirect::orientation_type
short int orientation_type
Definition geo_element_indirect.h:38
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
delta
Float delta(Float f, Float g)
Definition mosolov_error_yield_surface.cc:29
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39