Rheolef  7.1
an efficient C++ finite element environment
reference_element_face_transformation.cc
Go to the documentation of this file.
1 
23 
24 namespace rheolef {
25 
26 // Let S and K such that S is the i-th side of K
27 // then map a point hat_x defined in the reference element side hat_S
28 // into tilde_x defined in the reference element tilde_K
29 // note: used to build a quadrature formulae on a side of K
30 
31 template<class T>
34  reference_element tilde_K,
35  const side_information_type& sid,
36  const point_basic<T>& sid_hat_x)
37 {
38  // d=1: tilde_x = tilde_a0
39  // d=2: tilde_x = tilde_a0 + sid_hat_x[0]*(tilde_a1 - tilde_a0)
40  // d=3: tilde_x = tilde_a0 + sid_hat_x[0]*(tilde_a1 - tilde_a0) + sid_hat_x[1]*(tilde_a2 - tilde_a0)
41  // where tilde_a(i) are the vertices of the transformed side tilde_S in tilde_K
43  size_type sid_dim = tilde_K.dimension()-1;
44  switch (sid_dim) {
45  case 0: {
46  check_macro (sid.shift == 0, "unexpected shift="<<sid.shift<<" on a0d side");
47  check_macro (sid.orient > 0, "unexpected orient="<<sid.orient<<" on a 0d side");
48  size_type i0 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, 0);
49  const point& tilde_a0 = tilde_K.vertex (i0);
50  return tilde_a0;
51  }
52  case 1: {
53  check_macro (sid.shift == 0, "unexpected shift="<<sid.shift<<" on an edge");
54  size_type i0 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, 0),
55  i1 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, 1);
56  if (sid.orient < 0) std::swap (i0, i1);
57  const point& tilde_a0 = tilde_K.vertex (i0);
58  const point& tilde_a1 = tilde_K.vertex (i1);
59  return tilde_a0 + sid_hat_x[0]*(tilde_a1 - tilde_a0);
60  }
61  case 2: {
62  reference_element hat_S = tilde_K.side(sid.loc_isid);
63  size_type nv = hat_S.size();
64  size_type i0 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, sid.shift%nv),
65  i1 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, (sid.shift+1)%nv),
66  i2 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, (sid.shift+nv-1)%nv);
67  if (sid.orient < 0) std::swap (i1, i2);
68  const point& tilde_a0 = tilde_K.vertex (i0);
69  const point& tilde_a1 = tilde_K.vertex (i1);
70  const point& tilde_a2 = tilde_K.vertex (i2);
71  if (hat_S.variant() == reference_element::t) {
72  return tilde_a0 + sid_hat_x[0]*(tilde_a1 - tilde_a0)
73  + sid_hat_x[1]*(tilde_a2 - tilde_a0);
74  } else {
75  return tilde_a0 + 0.5*(1+sid_hat_x[0])*(tilde_a1 - tilde_a0)
76  + 0.5*(1+sid_hat_x[1])*(tilde_a2 - tilde_a0);
77  }
78  }
79  default: {
80  fatal_macro ("invalid subgeo_dim="<<sid_dim);
81  return point_basic<T>();
82  }
83  }
84 }
85 template<class T>
88  reference_element tilde_K,
89  const side_information_type& sid,
90  const point_basic<T>& tilde_x)
91 {
92  // d=1: sid_hat_x = 0
93  // d=2: sid_hat_x = dot2d(tilde_a1 - tilde_a0, tilde_x - tilde_a0)/dot2d(tilde_a1 - tilde_a0,tilde_a1 - tilde_a0)
94  // d=3: ?
95  // where tilde_ai are the vertices of the transformed side tilde_S in tilde_K
97  size_type sid_dim = tilde_K.dimension()-1;
98  switch (sid_dim) {
99  case 0: {
100  return point_basic<T>(0);
101  }
102  case 1: {
103  check_macro (sid.shift == 0, "unexpected shift="<<sid.shift<<" on an edge");
104  size_type i0 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, 0),
105  i1 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, 1);
106  if (sid.orient < 0) std::swap (i0, i1);
107  const point& tilde_a0 = tilde_K.vertex (i0);
108  const point& tilde_a1 = tilde_K.vertex (i1);
109  T deno = sqr(tilde_a1[0] - tilde_a0[0])
110  + sqr(tilde_a1[1] - tilde_a0[1]);
111  T num = (tilde_x [0] - tilde_a0[0])
112  *(tilde_a1[0] - tilde_a0[0])
113  + (tilde_x [1] - tilde_a0[1])
114  *(tilde_a1[1] - tilde_a0[1]);
115  return point_basic<T>(num/deno);
116  }
117  case 2: {
118  reference_element hat_S = tilde_K.side(sid.loc_isid);
119  size_type nv = hat_S.size();
120  size_type i0 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, sid.shift%nv),
121  i1 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, (sid.shift+1)%nv),
122  i2 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, (sid.shift+nv-1)%nv);
123  if (sid.orient < 0) std::swap (i1, i2);
124  const point& tilde_a0 = tilde_K.vertex (i0);
125  const point& tilde_a1 = tilde_K.vertex (i1);
126  const point& tilde_a2 = tilde_K.vertex (i2);
127  point_basic<T> hat_x;
128  hat_x[0] = dot(tilde_x - tilde_a0, tilde_a1 - tilde_a0)/norm2(tilde_a1 - tilde_a0);
129  hat_x[1] = dot(tilde_x - tilde_a0, tilde_a2 - tilde_a0)/norm2(tilde_a2 - tilde_a0);
130  return hat_x;
131  }
132  default: {
133  fatal_macro ("invalid subgeo_dim="<<sid_dim);
134  return point_basic<T>();
135  }
136  }
137  return point_basic<T>();
138 }
139 // NOTE: dirty implementation: goes to Float and go-back to size_t
140 // TODO: how to directly do it with size_t ?
141 point_basic<size_t>
143  reference_element tilde_K,
144  const side_information_type& sid,
145  size_t k,
146  const point_basic<size_t>& sid_ilat)
147 {
149  size_type sid_dim = tilde_K.dimension()-1;
150  Float km = (k == 0) ? 1 : k; // manage also k=0 case
151  switch (sid_dim) {
152  case 0: {
153  check_macro (sid.shift == 0, "unexpected shift="<<sid.shift<<" on a0d side");
154  check_macro (sid.orient > 0, "unexpected orient="<<sid.orient<<" on a 0d side");
155  size_type i0 = tilde_K.subgeo_local_vertex (sid_dim, sid.loc_isid, 0);
156  return point_basic<size_type>(i0);
157  }
158  case 1: {
159  point_basic<Float> sid_hat_x (sid_ilat[0]/km);
160  point_basic<Float> hat_x = reference_element_face_transformation (tilde_K, sid, sid_hat_x);
162  if (tilde_K.variant() == reference_element::t) {
163  ilat = point_basic<size_type>(size_type(k*hat_x[0]+0.5), size_type(k*hat_x[1]+0.5));
164  } else { // tilde_K=q
165  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));
166  }
167  return ilat;
168  }
169  case 2: {
170  reference_element hat_S = tilde_K.side(sid.loc_isid);
171  point_basic<Float> sid_hat_x;
172  if (hat_S.variant() == reference_element::t) {
173  sid_hat_x = point_basic<Float>(sid_ilat[0]/km, sid_ilat[1]/km);
174  } else {
175  sid_hat_x = point_basic<Float>(2*sid_ilat[0]/km-1, 2*sid_ilat[1]/km-1);
176  }
177  point_basic<Float> hat_x = reference_element_face_transformation (tilde_K, sid, sid_hat_x);
179  if (tilde_K.variant() == reference_element::T) {
180  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));
181  } else if (tilde_K.variant() == reference_element::H) {
182  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));
183  } else { // tilde_K=P
184  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));
185  }
186  return ilat;
187  }
188  default: {
189  fatal_macro ("invalid subgeo_dim="<<sid_dim);
190  return point_basic<size_type>();
191  }
192  }
193  return point_basic<size_type>();
194 }
195 // ----------------------------------------------------------------------------
196 // instanciation in library
197 // ----------------------------------------------------------------------------
198 
199 #define _RHEOLEF_instanciation(T) \
200 template \
201 point_basic<T> \
202 reference_element_face_transformation ( \
203  reference_element tilde_K, \
204  const side_information_type& sid, \
205  const point_basic<T>& hat_x); \
206 template \
207 point_basic<T> \
208 reference_element_face_inverse_transformation ( \
209  reference_element tilde_K, \
210  const side_information_type& sid, \
211  const point_basic<T>& tilde_x);
212 
214 
215 } // namespace rheolef
field::size_type size_type
Definition: branch.cc:425
see the Float page for the full documentation
see the point page for the full documentation
see the reference_element page for the full documentation
const point_basic< Float > & vertex(size_type iloc) const
static const variant_type H
reference_element side(size_type loc_isid) const
variant_type variant() const
size_type subgeo_local_vertex(size_type subgeo_dim, size_type loc_isid, size_type loc_jsidvert) const
std::vector< int >::size_type size_type
static const variant_type T
static const variant_type t
size_t size_type
Definition: basis_get.cc:76
point_basic< T >
Definition: piola_fem.h:135
Expr1::float_type T
Definition: field_expr.h:261
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
Definition: vec_expr_v2.h:415
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) _RHEOLEF_instanciation_evaluate(Float
T norm2(const vec< T, M > &x)
norm2(x): see the expression page for the full documentation
Definition: vec.h:379
point_basic< T > reference_element_face_transformation(reference_element tilde_K, const side_information_type &sid, const point_basic< T > &sid_hat_x)
geo_element_indirect::orientation_type orient