Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
reference_element_face_transformation.cc
Go to the documentation of this file.
1
21
23
24namespace 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
31template<class T>
32point_basic<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}
85template<class T>
86point_basic<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 ?
141point_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) \
200template \
201point_basic<T> \
202reference_element_face_transformation ( \
203 reference_element tilde_K, \
204 const side_information_type& sid, \
205 const point_basic<T>& hat_x); \
206template \
207point_basic<T> \
208reference_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
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
field::size_type size_type
Definition branch.cc:430
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
#define fatal_macro(message)
Definition dis_macros.h:33
Expr1::float_type T
Definition field_expr.h:230
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.
point_basic< T > reference_element_face_transformation(reference_element tilde_K, const side_information_type &sid, const point_basic< T > &sid_hat_x)
T norm2(const vec< T, M > &x)
norm2(x): see the expression page for the full documentation
Definition vec.h:379
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
point_basic< T > reference_element_face_inverse_transformation(reference_element tilde_K, const side_information_type &sid, const point_basic< T > &tilde_x)
geo_element_indirect::orientation_type orient