Loading...
Searching...
No Matches
reference_element_face_transformation.cc
Go to the documentation of this file.
#include "reference_element_face_transformation.h"
namespace rheolef {
// Let S and K such that S is the i-th side of K
// then map a point hat_x defined in the reference element side hat_S
// into tilde_x defined in the reference element tilde_K
// note: used to build a quadrature formulae on a side of K
template<class T>
point_basic<T>
reference_element_face_transformation (
reference_element tilde_K,
const side_information_type& sid,
const point_basic<T>& sid_hat_x)
{
// d=1: tilde_x = tilde_a0
// d=2: tilde_x = tilde_a0 + sid_hat_x[0]*(tilde_a1 - tilde_a0)
// d=3: tilde_x = tilde_a0 + sid_hat_x[0]*(tilde_a1 - tilde_a0) + sid_hat_x[1]*(tilde_a2 - tilde_a0)
// where tilde_a(i) are the vertices of the transformed side tilde_S in tilde_K
typedef reference_element::size_type size_type;
size_type sid_dim = tilde_K.dimension()-1;
switch (sid_dim) {
case 0: {
check_macro (sid.shift == 0, "unexpected shift="<<sid.shift<<" on a0d side");
check_macro (sid.orient > 0, "unexpected orient="<<sid.orient<<" on a 0d side");
size_type i0 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, 0);
const point& tilde_a0 = tilde_K.vertex (i0);
return tilde_a0;
}
case 1: {
check_macro (sid.shift == 0, "unexpected shift="<<sid.shift<<" on an edge");
size_type i0 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, 0),
i1 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, 1);
if (sid.orient < 0) std::swap (i0, i1);
const point& tilde_a0 = tilde_K.vertex (i0);
const point& tilde_a1 = tilde_K.vertex (i1);
return tilde_a0 + sid_hat_x[0]*(tilde_a1 - tilde_a0);
}
case 2: {
reference_element hat_S = tilde_K.side(sid.loc_isid);
size_type nv = hat_S.size();
size_type i0 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, sid.shift%nv),
i1 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, (sid.shift+1)%nv),
i2 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, (sid.shift+nv-1)%nv);
if (sid.orient < 0) std::swap (i1, i2);
const point& tilde_a0 = tilde_K.vertex (i0);
const point& tilde_a1 = tilde_K.vertex (i1);
const point& tilde_a2 = tilde_K.vertex (i2);
if (hat_S.variant() == reference_element::t) {
return tilde_a0 + sid_hat_x[0]*(tilde_a1 - tilde_a0)
+ sid_hat_x[1]*(tilde_a2 - tilde_a0);
} else {
return tilde_a0 + 0.5*(1+sid_hat_x[0])*(tilde_a1 - tilde_a0)
+ 0.5*(1+sid_hat_x[1])*(tilde_a2 - tilde_a0);
}
}
default: {
fatal_macro ("invalid subgeo_dim="<<sid_dim);
return point_basic<T>();
}
}
}
template<class T>
point_basic<T>
reference_element_face_inverse_transformation (
reference_element tilde_K,
const side_information_type& sid,
const point_basic<T>& tilde_x)
{
// d=1: sid_hat_x = 0
// d=2: sid_hat_x = dot2d(tilde_a1 - tilde_a0, tilde_x - tilde_a0)/dot2d(tilde_a1 - tilde_a0,tilde_a1 - tilde_a0)
// d=3: ?
// where tilde_ai are the vertices of the transformed side tilde_S in tilde_K
typedef reference_element::size_type size_type;
size_type sid_dim = tilde_K.dimension()-1;
switch (sid_dim) {
case 0: {
return point_basic<T>(0);
}
case 1: {
check_macro (sid.shift == 0, "unexpected shift="<<sid.shift<<" on an edge");
size_type i0 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, 0),
i1 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, 1);
if (sid.orient < 0) std::swap (i0, i1);
const point& tilde_a0 = tilde_K.vertex (i0);
const point& tilde_a1 = tilde_K.vertex (i1);
T deno = sqr(tilde_a1[0] - tilde_a0[0])
+ sqr(tilde_a1[1] - tilde_a0[1]);
T num = (tilde_x [0] - tilde_a0[0])
*(tilde_a1[0] - tilde_a0[0])
+ (tilde_x [1] - tilde_a0[1])
*(tilde_a1[1] - tilde_a0[1]);
return point_basic<T>(num/deno);
}
case 2: {
reference_element hat_S = tilde_K.side(sid.loc_isid);
size_type nv = hat_S.size();
size_type i0 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, sid.shift%nv),
i1 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, (sid.shift+1)%nv),
i2 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, (sid.shift+nv-1)%nv);
if (sid.orient < 0) std::swap (i1, i2);
const point& tilde_a0 = tilde_K.vertex (i0);
const point& tilde_a1 = tilde_K.vertex (i1);
const point& tilde_a2 = tilde_K.vertex (i2);
point_basic<T> hat_x;
hat_x[0] = dot(tilde_x - tilde_a0, tilde_a1 - tilde_a0)/norm2(tilde_a1 - tilde_a0);
hat_x[1] = dot(tilde_x - tilde_a0, tilde_a2 - tilde_a0)/norm2(tilde_a2 - tilde_a0);
return hat_x;
}
default: {
fatal_macro ("invalid subgeo_dim="<<sid_dim);
return point_basic<T>();
}
}
return point_basic<T>();
}
// NOTE: dirty implementation: goes to Float and go-back to size_t
// TODO: how to directly do it with size_t ?
point_basic<size_t>
reference_element_face_transformation (
reference_element tilde_K,
const side_information_type& sid,
size_t k,
const point_basic<size_t>& sid_ilat)
{
typedef reference_element::size_type size_type;
size_type sid_dim = tilde_K.dimension()-1;
Float km = (k == 0) ? 1 : k; // manage also k=0 case
switch (sid_dim) {
case 0: {
check_macro (sid.shift == 0, "unexpected shift="<<sid.shift<<" on a0d side");
check_macro (sid.orient > 0, "unexpected orient="<<sid.orient<<" on a 0d side");
size_type i0 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, 0);
return point_basic<size_type>(i0);
}
case 1: {
point_basic<Float> sid_hat_x (sid_ilat[0]/km);
point_basic<Float> hat_x = reference_element_face_transformation (tilde_K, sid, sid_hat_x);
point_basic<size_type> ilat;
if (tilde_K.variant() == reference_element::t) {
ilat = point_basic<size_type>(size_type(k*hat_x[0]+0.5), size_type(k*hat_x[1]+0.5));
} else { // tilde_K=q
ilat = point_basic<size_type>(size_type(k*(1+hat_x[0])/2+0.5), size_type(k*(1+hat_x[1])/2+0.5));
}
return ilat;
}
case 2: {
reference_element hat_S = tilde_K.side(sid.loc_isid);
point_basic<Float> sid_hat_x;
if (hat_S.variant() == reference_element::t) {
sid_hat_x = point_basic<Float>(sid_ilat[0]/km, sid_ilat[1]/km);
} else {
sid_hat_x = point_basic<Float>(2*sid_ilat[0]/km-1, 2*sid_ilat[1]/km-1);
}
point_basic<Float> hat_x = reference_element_face_transformation (tilde_K, sid, sid_hat_x);
point_basic<size_type> ilat;
if (tilde_K.variant() == reference_element::T) {
ilat = point_basic<size_type>(size_type(k*hat_x[0]+0.5), size_type(k*hat_x[1]+0.5), size_type(k*hat_x[2]+0.5));
} else if (tilde_K.variant() == reference_element::H) {
ilat = point_basic<size_type>(size_type(k*(1+hat_x[0])/2+0.5), size_type(k*(1+hat_x[1])/2+0.5), size_type(k*(1+hat_x[2])/2+0.5));
} else { // tilde_K=P
ilat = point_basic<size_type>(size_type(k*hat_x[0]+0.5), size_type(k*hat_x[1]+0.5), size_type(k*(1+hat_x[2])/2+0.5));
}
return ilat;
}
default: {
fatal_macro ("invalid subgeo_dim="<<sid_dim);
return point_basic<size_type>();
}
}
return point_basic<size_type>();
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T) \
template \
point_basic<T> \
reference_element_face_transformation ( \
reference_element tilde_K, \
const side_information_type& sid, \
const point_basic<T>& hat_x); \
template \
point_basic<T> \
reference_element_face_inverse_transformation ( \
reference_element tilde_K, \
const side_information_type& sid, \
const point_basic<T>& tilde_x);
_RHEOLEF_instanciation(Float)
} // namespace rheolef
_RHEOLEF_instanciation
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
size_type
field::size_type size_type
Definition branch.cc:430
Float
see the Float page for the full documentation
rheolef::point_basic
Definition point.h:86
rheolef::reference_element
see the reference_element page for the full documentation
Definition reference_element.h:66
rheolef::reference_element::vertex
const point_basic< Float > & vertex(size_type iloc) const
Definition reference_element.cc:251
rheolef::reference_element::H
static const variant_type H
Definition reference_element.h:81
rheolef::reference_element::side
reference_element side(size_type loc_isid) const
Definition reference_element.h:119
rheolef::reference_element::size
size_type size() const
Definition reference_element.h:102
rheolef::reference_element::dimension
size_type dimension() const
Definition reference_element.h:101
rheolef::reference_element::variant
variant_type variant() const
Definition reference_element.h:99
rheolef::reference_element::subgeo_local_vertex
size_type subgeo_local_vertex(size_type subgeo_dim, size_type loc_isid, size_type loc_jsidvert) const
Definition reference_element.h:122
rheolef::reference_element::size_type
std::vector< int >::size_type size_type
Definition reference_element.h:71
rheolef::reference_element::T
static const variant_type T
Definition reference_element.h:79
rheolef::reference_element::t
static const variant_type t
Definition reference_element.h:77
rheolef::point
point_basic< Float > point
Definition point.h:163
rheolef::Float
double Float
see the Float page for the full documentation
Definition Float.h:143
fatal_macro
#define fatal_macro(message)
Definition dis_macros.h:33
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)")
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::reference_element_face_transformation
point_basic< T > reference_element_face_transformation(reference_element tilde_K, const side_information_type &sid, const point_basic< T > &sid_hat_x)
Definition reference_element_face_transformation.cc:33
rheolef::norm2
T norm2(const vec< T, M > &x)
norm2(x): see the expression page for the full documentation
Definition vec.h:379
rheolef::reference_element_face_inverse_transformation
point_basic< T > reference_element_face_inverse_transformation(reference_element tilde_K, const side_information_type &sid, const point_basic< T > &tilde_x)
Definition reference_element_face_transformation.cc:87
rheolef::dot
rheolef::std::enable_if< details::is_vec_expr_v2_arg< Expr1 >::value &&details::is_vec_expr_v2_arg< Expr2 >::value, promote< Expr1::float_type, Expr2::float_type >::type > type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
Definition vec_expr_v2.h:415
reference_element_face_transformation.h
rheolef::side_information_type
Definition reference_element_face_transformation.h:37
rheolef::side_information_type::orient
geo_element_indirect::orientation_type orient
Definition reference_element_face_transformation.h:40
rheolef::side_information_type::loc_isid
size_t loc_isid
Definition reference_element_face_transformation.h:38
rheolef::side_information_type::shift
size_t shift
Definition reference_element_face_transformation.h:39