rheolef  6.5
geo_element_v4.cc
Go to the documentation of this file.
1 
2 #include "rheolef/geo_element.h"
3 
4 namespace rheolef {
5 
6 // i/o
7 char
8 skip_blancs_and_tabs (std::istream& is)
9 {
10  if (!is.good()) return 0;
11  char c = is.peek();
12  while (is.good() && (c == ' ' || c == '\t')) {
13  is >> c;
14  } while (is.good() && (c == ' ' || c == '\t'));
15  return c;
16 }
17 std::istream&
18 operator>> (std::istream& is, geo_element& K)
19 {
21  extern char skip_blancs_and_tabs (std::istream&);
22  char c;
23  is >> std::ws >> c;
24  if (!isdigit(c)) {
26  check_macro (variant != reference_element::max_variant, "undefined element variant `"<<c<<"'");
27  size_type order = 1;
28  char c2;
29  is >> std::ws >> c2;
30  if (c2 == 'p') {
31  is >> order;
32  } else {
33  is.unget(); // or is.putback(c); ?
34  }
35  K.reset (variant, order);
36  for (size_type i = 0, n = K.n_node(); i < n; i++) {
37  is >> K[i];
38  }
39  return is;
40  }
41  size_type tmp [3];
42  size_type nvertex = 0;
43  while (is.good() && isdigit(c) && nvertex < 3) {
44  is.unget(); // or is.putback(c); ?
45  is >> tmp [nvertex];
46  nvertex++;
47  c = skip_blancs_and_tabs (is);
48  }
49  size_type variant = reference_element::max_variant;
50  switch (nvertex) {
51  case 1: variant = reference_element::p; break;
52  case 2: variant = reference_element::e; break;
53  case 3: variant = reference_element::t; break;
54  default: error_macro ("unexpected element with "<<nvertex<<" vertices");
55  }
56  K.reset (variant, 1);
57  for (size_type i = 0, n = K.n_node(); i < n; i++) {
58  K[i] = tmp[i];
59  }
60  return is;
61 }
62 std::ostream&
63 operator<< (std::ostream& os, const geo_element& K)
64 {
66  os << K.name() << "\tp" << K.order();
67  for (size_type loc_inod = 0, loc_nnod = K.n_node(); loc_inod < loc_nnod; loc_inod++) {
68  os << " " << K[loc_inod];
69  }
70  return os;
71 }
78 bool
80  const geo_element& S,
81  orientation_type& orient,
82  shift_type& shift) const
83 {
85  "get_orientation_and_shift: elements have different dimensions "<<dimension()<<" and "<<S.dimension());
86  check_macro (S.size() == size(),
87  "get_orientation_and_shift: elements have different sizes "<<size()<<" and "<<S.size());
88  if (S.dimension() == 0) {
89  orient = 1;
90  shift = 0;
91  return true;
92  }
93  if (S.dimension() == 1) {
94  orient = (operator[](0) == S[0]) ? 1 : -1;
95  shift = 0;
96  return true;
97  }
98  if (S.dimension() == 2) {
99  size_type n = size();
100  for (shift = 0; shift < shift_type(n); shift++) {
101  if (operator[](0) == S[shift]) break;
102  }
103  if (shift == shift_type(n)) {
104  orient = 0;
106  return false;
107  }
108  orient = (operator[](1) == S[(shift+1)%n]) ? 1 : -1;
109  return true;
110  }
111  orient = 1;
112  shift = 0;
113  return true;
114 }
117  size_type dis_iv0, size_type dis_iv1) const
118 {
119  return (operator[](0) == dis_iv0) ? 1 : -1;
120 }
121 void
123  size_type dis_iv0, size_type dis_iv1, size_type dis_iv2,
124  orientation_type& orient,
125  shift_type& shift) const
126 {
127  if (operator[](0) == dis_iv0) {
128  orient = (operator[](1) == dis_iv1) ? 1 : -1;
129  shift = 0;
130  } else if (operator[](0) == dis_iv1) {
131  orient = (operator[](1) == dis_iv2) ? 1 : -1;
132  shift = 1;
133  } else {
134  orient = (operator[](1) == dis_iv0) ? 1 : -1;
135  shift = 2;
136  }
137 }
138 void
140  size_type dis_iv0, size_type dis_iv1, size_type dis_iv2, size_type dis_iv3,
141  orientation_type& orient,
142  shift_type& shift) const
143 {
144  if (operator[](0) == dis_iv0) {
145  orient = (operator[](1) == dis_iv1) ? 1 : -1;
146  shift = 0;
147  } else if (operator[](0) == dis_iv1) {
148  orient = (operator[](1) == dis_iv2) ? 1 : -1;
149  shift = 1;
150  } else if (operator[](0) == dis_iv2) {
151  orient = (operator[](1) == dis_iv3) ? 1 : -1;
152  shift = 2;
153  } else {
154  orient = (operator[](1) == dis_iv0) ? 1 : -1;
155  shift = 3;
156  }
157 }
160  const geo_element& S,
161  size_type& loc_isid,
162  size_type& shift) const
163 {
164  loc_isid = shift = 0;
165  check_macro (S.dimension() + 1 == dimension(),
166  "get_side_orientation: side have unexpected dimension "<<S.dimension()<<": "
167  <<dimension()<<" was expected");
168  const geo_element& K = *this;
169  size_type side_dim = S.dimension();
170  for (size_type loc_nsid = K.n_subgeo(side_dim); loc_isid < loc_nsid; loc_isid++) {
171  size_type sid_nloc = K.subgeo_size (side_dim, loc_isid);
172  if (sid_nloc != S.size()) continue;
173  for (shift = 0; shift < sid_nloc; shift++) {
174  size_type loc_jv = K.subgeo_local_vertex (side_dim, loc_isid, shift);
175  if (K[loc_jv] != S[0]) continue;
176  bool matches = true;
177  for (size_type sid_kloc = 1; sid_kloc < sid_nloc; sid_kloc++) {
178  size_type loc_kv = K.subgeo_local_vertex (side_dim, loc_isid, (shift+sid_kloc) % sid_nloc);
179  if (K[loc_kv] != S[sid_kloc]) { matches = false; break; }
180  }
181  if (matches) {
182  if (side_dim == 1 && shift == 1) {
183  shift = 0;
184  return -1;
185  }
186  return 1;
187  }
188  matches = true;
189  for (size_type sid_kloc = 1; sid_kloc < sid_nloc; sid_kloc++) {
190  size_type loc_kv = K.subgeo_local_vertex (side_dim, loc_isid, (shift+sid_nloc-sid_kloc) % sid_nloc);
191  if (K[loc_kv] != S[sid_kloc]) { matches = false; break; }
192  }
193  if (matches) { return -1; }
194  }
195  }
196  fatal_macro ("get_side_orientation: side is not part of the element");
197  return 0; // not reached
198 }
199 void
201  const geo_element& S,
202  side_information_type& sid) const
203 {
204  sid.orient = get_side_informations (S, sid.loc_isid, sid.shift);
205  sid.dim = S.dimension();
206  sid.n_vertex = subgeo_size (sid.dim, sid.loc_isid);
207  sid.hat.set_variant (sid.n_vertex, sid.dim);
208 }
211 {
212  size_type loc_isid, shift;
213  return get_side_informations (S, loc_isid, shift);
214 }
217  orientation_type orient,
218  size_type order,
219  size_type loc_iedg_j)
220 {
221  if (orient > 0) return loc_iedg_j;
222  return order - 2 - loc_iedg_j;
223 }
226  const geo_element& K,
227  size_type loc_iedg,
228  size_type loc_iedg_j,
229  size_type order)
230 {
231  // shift and/or inverse-orient the lattice(i,j)
232  if (K.dimension() < 2 || order <= 2) {
233  return loc_iedg_j;
234  }
235  return fix_edge_indirect (K.edge_indirect(loc_iedg).orientation(), order, loc_iedg_j);
236 }
237 void
239  size_type loc_tri_inod,
240  size_type order,
241  point_basic<size_type>& ij_lattice)
242 {
243  // loc_tri_inod = (n_face_node - (j1-1)*j1/2) + i-1, where j1=order-1-j
244  // i - 1 = 0 = loc_tri_inod - (n_face_node - (j1-1)*j1/2);
245  // <=> j1^2 - j1 - 2*(n_face_node - loc_tri_inod) = 0
246  size_type n_face_node = (order-1)*(order-2)/2;
247  double a = 1;
248  double b = -1;
249  double c = -2*int(n_face_node - loc_tri_inod);
250  double delta = b*b - 4*a*c;
251  check_macro (delta >= 0, "loc_tri_inod2lattice(loc_tri_inod="<<loc_tri_inod<<",order="<<order<<"): unexpected delta="<<delta<<" < 0");
252  double j1_plus = (-b + sqrt(delta))/(2*a);
253  // j1_minus = (-b - sqrt(delta))/(2*a); is negative always
254  size_type j = floor (order - j1_plus);
255  size_type j1 = order - j;
256  size_type i = loc_tri_inod - (n_face_node - (j1-1)*j1/2) + 1 ;
257  ij_lattice = point_basic<size_type>(i,j);
258 }
261  orientation_type orient,
262  shift_type shift,
263  size_type order,
264  size_type loc_itri_j)
265 {
266  point_basic<size_type> ij_lattice;
267  loc_tri_inod2lattice (loc_itri_j, order, ij_lattice);
268  size_type i = ij_lattice[0];
269  size_type j = ij_lattice[1];
270  // 2) then shift and/or inverse-orient (i,j)
271  size_type coord[3];
272  coord [0] = i;
273  coord [1] = j;
274  coord [2] = order - i - j;
275  size_type i_fix = coord [shift%3];
276  size_type j_fix = coord [(shift+1)%3];
277  if (orient < 0) std::swap (i_fix, j_fix);
278  size_type loc_tri_nnod = (order-1)*(order-2)/2;
279  size_type j1_fix = order - j_fix;
280  size_type loc_itri_j_fix = (loc_tri_nnod - (j1_fix-1)*j1_fix/2) + (i_fix-1);
281  return loc_itri_j_fix;
282 }
285  const geo_element& K,
286  size_type loc_ifac,
287  size_type loc_itri_j,
288  size_type order)
289 {
290  // shift and/or inverse-orient the lattice(i,j)
291  if (K.dimension() < 3 ||
292  (K.face_indirect(loc_ifac).orientation() > 0 && K.face_indirect(loc_ifac).shift() == 0)) {
293  return loc_itri_j;
294  }
295  return fix_triangle_indirect (
296  K.face_indirect(loc_ifac).orientation(),
297  K.face_indirect(loc_ifac).shift(),
298  order,
299  loc_itri_j);
300 }
301 void
303  size_type loc_qua_inod,
304  size_type order,
305  point_basic<size_type>& ij_lattice)
306 {
307  ij_lattice[0] = (loc_qua_inod % (order-1)) + 1;
308  ij_lattice[1] = (loc_qua_inod / (order-1)) + 1;
309 }
312  orientation_type orient,
313  shift_type shift,
314  size_type order,
315  size_type loc_iqua_j)
316 {
317  point_basic<size_type> ij_lattice;
318  loc_qua_inod2lattice (loc_iqua_j, order, ij_lattice);
319  size_type i = ij_lattice[0];
320  size_type j = ij_lattice[1];
321  // 2) then shift and/or inverse-orient (i,j)
322  size_type coord [4];
323  coord [0] = i;
324  coord [1] = j;
325  coord [2] = order - i;
326  coord [3] = order - j;
327  size_type i_fix = coord [shift%4];
328  size_type j_fix = coord [(shift+1)%4];
329  if (orient < 0) std::swap (i_fix, j_fix);
330  size_type loc_iqua_j_fix = (order-1)*(j_fix-1) + (i_fix-1);
331  return loc_iqua_j_fix;
332 }
335  const geo_element& K,
336  size_type loc_ifac,
337  size_type loc_iqua_j,
338  size_type order)
339 {
340  // shift and/or inverse-orient the lattice(i,j)
341  if (K.dimension() < 3 || order <= 2 ||
342  (K.face_indirect(loc_ifac).orientation() > 0 && K.face_indirect(loc_ifac).shift() == 0)) {
343  return loc_iqua_j;
344  }
345  return fix_quadrangle_indirect (
346  K.face_indirect(loc_ifac).orientation(),
347  K.face_indirect(loc_ifac).shift(),
348  order,
349  loc_iqua_j);
350 }
351 
352 } // namespace rheolef
353