1#ifndef _RHEOLEF_GEO_ELEMENT_CURVED_BALL_H
2#define _RHEOLEF_GEO_ELEMENT_CURVED_BALL_H
64 T x = 0.5*(hat_x[0]+1);
65 T y = 0.5*(hat_x[1]+1);
141 T la = 1 - hat_x[0] - hat_x[1];
144 return 0.5*( la/(1-lb)*
edge(0,lb) + lb/(1-la)*
edge(0,1-la)
145 + lb/(1-lc)*
edge(1,lc) + lc/(1-lb)*
edge(1,1-lb)
146 + lc/(1-la)*
edge(2,la) + la/(1-lc)*
edge(2,1-lc)
165 std::array<point_basic<T>,3>
node;
194 for (
size_t loc_ifac = 0; loc_ifac < 4; loc_ifac++) {
200 for (
size_t loc_ifac_jedg = 0; loc_ifac_jedg < 3; loc_ifac_jedg++) {
209 for (
size_t loc_ifac = 0; loc_ifac < 4; loc_ifac++) {
210 for (
size_t loc_ifac_jedg = 0; loc_ifac_jedg < 3; loc_ifac_jedg++) {
212 if (loc_jedg != loc_iedg_curved)
continue;
235 l[0] = 1 - hat_x[0] - hat_x[1] - hat_x[2];
241 std::array<size_t,3> S;
242 std::array<bool,4> is_on_S;
244 for (
size_t loc_ifac = 0; loc_ifac < 4; loc_ifac++) {
245 is_on_S.fill (
false);
246 for (
size_t loc_ifac_jv = 0; loc_ifac_jv < 3; loc_ifac_jv++) {
248 is_on_S [S[loc_ifac_jv]] =
true;
251 size_t jv_3 = size_t(-1);
252 for (
size_t loc_jv = 0; loc_jv < 4; loc_jv++) {
253 if (!is_on_S [loc_jv]) jv_3 = loc_jv;
255 T sum = l[S[0]] + l[S[1]] + l[S[2]];
257 x_all_faces += (1 - l[jv_3])*
x_face (loc_ifac, hat_xi);
261 for (
size_t loc_iedg = 0; loc_iedg < 6; loc_iedg++) {
265 size_t ic = size_t(-1);
266 size_t id = size_t(-1);
267 for (
size_t loc_iv = 0; loc_iv < 4; loc_iv++) {
268 if (loc_iv == ia || loc_iv == ib)
continue;
269 if (ic ==
size_t(-1)) { ic = loc_iv;
continue; }
272 assert_macro (ic !=
size_t(-1) &&
id !=
size_t(-1),
"unexpected");
273 T hat_xi = l[ib]/(l[ia] + l[ib]);
274 x_all_edges += (1 - l[ic] - l[id])*
x_edge (loc_iedg, hat_xi);
278 for (
size_t loc_iv = 0; loc_iv < 4; loc_iv++) {
279 x_all_vertices += l[loc_iv]*
node[loc_iv];
281 return x_all_faces - x_all_edges + x_all_vertices;
294 size_t loc_ifac_jedg_bdry = size_t(-1);
295 for (
size_t loc_ifac_jedg = 0; loc_ifac_jedg < 3; loc_ifac_jedg++) {
298 loc_ifac_jedg_bdry = loc_ifac_jedg;
302 assert_macro (loc_ifac_jedg_bdry !=
size_t(-1),
"unexpected");
306 point_basic<T> x = (1 - hat_x[0] - hat_x[1])*a + hat_x[0]*b + hat_x[1]*c;
340 std::array<point_basic<T>,4>
node;
see the edge page for the full documentation
point_basic< T > operator()(const point_basic< T > &z) const
curved_ball_H(const point_basic< T > &a0, const point_basic< T > &b0, const point_basic< T > &c0, const point_basic< T > &d0, const point_basic< T > &e0, const point_basic< T > &f0, const point_basic< T > &g0, const point_basic< T > &h0, const point_basic< T > ¢er0=point_basic< T >(0, 0), const T &radius0=1)
curved_ball_T(const point_basic< T > &a0, const point_basic< T > &b0, const point_basic< T > &c0, const point_basic< T > &d0, const point_basic< T > ¢er0=point_basic< T >(0, 0), const T &radius0=1)
point_basic< T > project_on_boundary(const point_basic< T > &x) const
point_basic< T > f_bcd(const point_basic< T > &x) const
std::array< bool, 4 > is_curved_fac
point_basic< T > f_cd(const point_basic< T > &x) const
std::array< bool, 4 > is_bdry_fac
point_basic< T > f_acd(const point_basic< T > &x) const
point_basic< T > f_ab(const point_basic< T > &x) const
std::array< bool, 6 > is_bdry_edg
point_basic< T > f_ac(const point_basic< T > &x) const
point_basic< T > operator()(const point_basic< T > &hat_x) const
point_basic< T > f_bc(const point_basic< T > &x) const
point_basic< T > f_bd(const point_basic< T > &x) const
void set_boundary_face(size_t loc_ifac_curved)
point_basic< T > f_abd(const point_basic< T > &x) const
point_basic< T > x_edge(size_t loc_iedg, const T &hat_x) const
point_basic< T > f_abc(const point_basic< T > &x) const
point_basic< T > x_face(size_t loc_ifac, const point_basic< T > &hat_x) const
void set_boundary_edge(size_t loc_iedg_curved)
point_basic< T > f_ad(const point_basic< T > &x) const
std::array< point_basic< T >, 4 > node
point_basic< T > f_small_carre(const T &x, const T &y) const
point_basic< T > operator()(const point_basic< T > &hat_x) const
curved_ball_q(const point_basic< T > &a0, const point_basic< T > &b0, const point_basic< T > &c0, const point_basic< T > &d0, const point_basic< T > ¢er0=point_basic< T >(0, 0), const T &radius0=1)
point_basic< T > f_ab(const T &s) const
point_basic< T > f_bc(const T &s) const
point_basic< T > f_ad(const T &s) const
point_basic< T > f_dc(const T &s) const
point_basic< T > project_on_boundary(const point_basic< T > &x) const
curved_ball_t(const point_basic< T > &a0, const point_basic< T > &b0, const point_basic< T > &c0, size_t loc_curved_iedg, const point_basic< T > ¢er0=point_basic< T >(0, 0), const T &radius0=1)
std::array< bool, 3 > is_bdry_edg
point_basic< T > operator()(const point_basic< T > &hat_x) const
point_basic< T > edge(size_t loc_iedg, const T &hat_x) const
std::array< point_basic< T >, 3 > node
static size_type subgeo_local_node(size_type order, size_type side_dim, size_type loc_isid, size_type loc_jsidnod)
static size_type face2edge(size_type loc_iface, size_type loc_iface_jedg)
#define assert_macro(ok_condition, message)
#define error_macro(message)
This file is part of Rheolef.
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
point - d-dimensional physical point or vector