Rheolef  7.2
an efficient C++ finite element environment
geo_element.cc
Go to the documentation of this file.
1 
22 #include "rheolef/geo_element.h"
23 
24 #include "geo_element_aux.icc"
25 
26 namespace rheolef {
27 
28 // --------------------------------------------------------------------------
29 // i/o
30 // --------------------------------------------------------------------------
31 char
32 skip_blancs_and_tabs (std::istream& is)
33 {
34  if (!is.good()) return 0;
35  char c = is.peek();
36  while (is.good() && (c == ' ' || c == '\t')) {
37  is >> c;
38  } while (is.good() && (c == ' ' || c == '\t'));
39  return c;
40 }
41 void
42 geo_element::get (std::istream& is)
43 {
44  geo_element& K = *this; // for more readability
46  extern char skip_blancs_and_tabs (std::istream&);
47  char c;
48  is >> std::ws >> c;
49  if (!isdigit(c)) {
50  // have an element type: 'e', 't' etc
52  check_macro (variant != reference_element::max_variant, "undefined element variant `"<<c<<"'");
53  size_type order = 1;
54  char c2;
55  is >> std::ws >> c2;
56  if (c2 == 'p') {
57  // have a high order spec as "p2"
58  is >> order;
59  } else {
60  is.unget(); // or is.putback(c); ?
61  }
62  K.reset (variant, order);
63  for (size_type i = 0, n = K.n_node(); i < n; i++) {
64  is >> K[i];
65  }
66  return;
67  }
68  // here, starts with a digit: triangle (e.g. "21 22 23"), edge (e.g. "12 13") or point (e.g. "11")
69  // end of element input: with end-of-line or non-digit input char
70  // Note: with 4 vertices, there is an ambiguity between quadrangle and thetraedra
71  size_type tmp [3];
72  size_type nvertex = 0;
73  while (is.good() && isdigit(c) && nvertex < 3) {
74  is.unget(); // or is.putback(c); ?
75  is >> tmp [nvertex];
76  nvertex++;
77  c = skip_blancs_and_tabs (is);
78  }
80  switch (nvertex) {
81  case 1: variant = reference_element::p; break;
82  case 2: variant = reference_element::e; break;
83  case 3: variant = reference_element::t; break;
84  default: error_macro ("unexpected element with "<<nvertex<<" vertices");
85  }
86  K.reset (variant, 1);
87  for (size_type i = 0, n = K.n_node(); i < n; i++) {
88  K[i] = tmp[i];
89  }
90  return;
91 }
92 void
93 geo_element::put (std::ostream& os) const
94 {
96  os << name() << "\t";
97  if (order() > 1) {
98  os << "p" << order() << " ";
99  }
100  for (size_type loc_inod = 0, loc_nnod = n_node(); loc_inod < loc_nnod; loc_inod++) {
101  if (loc_inod != 0) os << " ";
102  os << operator[](loc_inod);
103  }
104 }
105 // --------------------------------------------------------------------------
107 
113 bool
115  const geo_element& S,
116  orientation_type& orient,
117  shift_type& shift) const
118 {
119  check_macro (S.dimension() == dimension(),
120  "get_orientation_and_shift: elements have different dimensions "<<dimension()<<" and "<<S.dimension());
121  check_macro (S.size() == size(),
122  "get_orientation_and_shift: elements have different sizes "<<size()<<" and "<<S.size());
123  if (S.dimension() == 0) {
124  orient = 1;
125  shift = 0;
126  return true;
127  }
128  if (S.dimension() == 1) {
129  orient = (operator[](0) == S[0]) ? 1 : -1;
130  shift = 0;
131  return true;
132  }
133  if (S.dimension() == 2) {
134  size_type n = size();
135  // 3d face: triangle or quadrangle
136  for (shift = 0; shift < shift_type(n); shift++) {
137  if (operator[](0) == S[shift]) break;
138  }
139  if (shift == shift_type(n)) {
140  orient = 0;
141  shift = std::numeric_limits<shift_type>::max();
142  return false;
143  }
144  orient = (operator[](1) == S[(shift+1)%n]) ? 1 : -1;
145  return true;
146  }
147  // S.dimension() == 3: TODO volumic domain: tetra can be rotated ?
148  orient = 1;
149  shift = 0;
150  return true;
151 }
152 // if edge (dis_iv0,dis_iv1) has the same orientation as current edge
153 // assume that current geo_element is an edge that contains the same vertices
156  size_type dis_iv0, size_type dis_iv1) const
157 {
158  return (operator[](0) == dis_iv0) ? 1 : -1;
159 }
160 // for a triangular side (dis_iv0,dis_iv1,dis_iv2)
161 // assume that current geo_element is a triangle that contains the same vertices
162 void
164  size_type dis_iv0, size_type dis_iv1, size_type dis_iv2,
165  orientation_type& orient,
166  shift_type& shift) const
167 {
168  geo_element_get_orientation_and_shift (*this, dis_iv0, dis_iv1, dis_iv2, orient, shift);
169 }
170 // for a quadrangular side (dis_iv0,dis_iv1,dis_iv2,dis_iv3)
171 // assume that current geo_element is a quadrangular that contains the same vertices
172 void
174  size_type dis_iv0, size_type dis_iv1, size_type dis_iv2, size_type dis_iv3,
175  orientation_type& orient,
176  shift_type& shift) const
177 {
178  geo_element_get_orientation_and_shift (*this, dis_iv0, dis_iv1, dis_iv2, dis_iv3, orient, shift);
179 }
180 // let K=(*this) one element and S one side of K
181 // S could have the good or the opposite orientation on K:
182 // gives its sign, local side index and node shift (in 3D)
183 // Note: assume that S is a valid side of K
186  const geo_element& S,
187  size_type& loc_isid,
188  size_type& shift) const
189 {
190  loc_isid = shift = 0;
191  check_macro (S.dimension() + 1 == dimension(),
192  "get_side_orientation: side have unexpected dimension "<<S.dimension()<<": "
193  <<dimension()<<"-1 was expected");
194  const geo_element& K = *this;
195  size_type side_dim = S.dimension();
196  for (size_type loc_nsid = K.n_subgeo(side_dim); loc_isid < loc_nsid; loc_isid++) {
197  size_type sid_nloc = K.subgeo_size (side_dim, loc_isid);
198  if (sid_nloc != S.size()) continue;
199  for (shift = 0; shift < sid_nloc; shift++) {
200  size_type loc_jv = K.subgeo_local_vertex (side_dim, loc_isid, shift);
201  if (K[loc_jv] != S[0]) continue;
202  // one node matches with S[0] on K: loc_isid and shift
203  // check others nodes in a first rotation direction:
204  bool matches = true;
205  for (size_type sid_kloc = 1; sid_kloc < sid_nloc; sid_kloc++) {
206  size_type loc_kv = K.subgeo_local_vertex (side_dim, loc_isid, (shift+sid_kloc) % sid_nloc);
207  if (K[loc_kv] != S[sid_kloc]) { matches = false; break; }
208  }
209  if (matches) {
210  if (side_dim == 1 && shift == 1) {
211  // shift=1 for an edge means a change of orientation of this edge (for DG)
212  shift = 0;
213  return -1;
214  }
215  return 1;
216  }
217  // check others nodes in the opposite rotation direction:
218  matches = true;
219  for (size_type sid_kloc = 1; sid_kloc < sid_nloc; sid_kloc++) {
220  size_type loc_kv = K.subgeo_local_vertex (side_dim, loc_isid, (shift+sid_nloc-sid_kloc) % sid_nloc);
221  if (K[loc_kv] != S[sid_kloc]) { matches = false; break; }
222  }
223  if (matches) { return -1; }
224  }
225  }
226  fatal_macro ("get_side_orientation: side is not part of the element");
227  return 0; // not reached
228 }
229 void
231  const geo_element& S,
232  side_information_type& sid) const
233 {
234  sid.orient = get_side_informations (S, sid.loc_isid, sid.shift);
235  sid.dim = S.dimension();
236  sid.n_vertex = subgeo_size (sid.dim, sid.loc_isid);
237  sid.hat.set_variant (sid.n_vertex, sid.dim);
238 }
241 {
242  size_type loc_isid, shift;
243  return get_side_informations (S, loc_isid, shift);
244 }
245 // =========================================================================
246 // fix rotation and orientation on a 2d edge or 3d face
247 // =========================================================================
248 
249 // --------------
250 // 1) edges
251 // --------------
254  const geo_element& K,
255  size_type loc_iedg,
256  size_type loc_iedg_j,
258 {
259  return geo_element_fix_edge_indirect (K, loc_iedg, loc_iedg_j, order);
260 }
261 // --------------
262 // 2) triangles
263 // --------------
264 void
266  size_type loc_tri_inod,
268  point_basic<size_type>& ij_lattice)
269 {
270  geo_element_loc_tri_inod2lattice (loc_tri_inod, order, ij_lattice);
271 }
274  orientation_type orient,
275  shift_type shift,
277  size_type loc_itri_j)
278 {
279  return geo_element_fix_triangle_indirect (orient, shift, order, loc_itri_j);
280 }
283  const geo_element& K,
284  size_type loc_ifac,
285  size_type loc_itri_j,
287 {
288  return geo_element_fix_triangle_indirect (K, loc_ifac, loc_itri_j, order);
289 }
290 // --------------
291 // 3) quadrangles
292 // --------------
293 void
295  size_type loc_qua_inod,
297  point_basic<size_type>& ij_lattice)
298 {
299  geo_element_loc_qua_inod2lattice (loc_qua_inod, order, ij_lattice);
300 }
303  orientation_type orient,
304  shift_type shift,
306  size_type loc_iqua_j)
307 {
308  return geo_element_fix_quadrangle_indirect (orient, shift, order, loc_iqua_j);
309 }
312  const geo_element& K,
313  size_type loc_ifac,
314  size_type loc_iqua_j,
316 {
317  return geo_element_fix_quadrangle_indirect (K, loc_ifac, loc_iqua_j, order);
318 }
321  const geo_element& K,
322  size_type subgeo_variant,
323  size_type loc_ige,
324  size_type loc_comp_idof_on_subgeo,
326 {
327  switch (subgeo_variant) {
328  case reference_element::e: return geo_element::fix_edge_indirect (K, loc_ige, loc_comp_idof_on_subgeo, order);
329  case reference_element::t: return geo_element::fix_triangle_indirect (K, loc_ige, loc_comp_idof_on_subgeo, order);
330  case reference_element::q: return geo_element::fix_quadrangle_indirect (K, loc_ige, loc_comp_idof_on_subgeo, order);
331  default: return loc_comp_idof_on_subgeo;
332  }
333 }
334 
335 } // namespace rheolef
see the geo_element page for the full documentation
Definition: geo_element.h:102
size_type & operator[](size_type loc_inod)
Definition: geo_element.h:183
geo_element_indirect::orientation_type orientation_type
Definition: geo_element.h:134
bool get_orientation_and_shift(const geo_element &S, orientation_type &orient, shift_type &shift) const
return orientation and shift between *this element and S
Definition: geo_element.cc:114
size_type n_node() const
Definition: geo_element.h:170
size_type subgeo_local_vertex(size_type subgeo_dim, size_type i_subgeo, size_type i_subgeo_vertex) const
Definition: geo_element.h:223
static size_type fix_edge_indirect(const geo_element &K, size_type loc_iedg, size_type loc_iedg_j, size_type order)
Definition: geo_element.cc:253
size_type size() const
Definition: geo_element.h:168
static size_type fix_triangle_indirect(const geo_element &K, size_type loc_itri, size_type loc_itri_j, size_type order)
Definition: geo_element.cc:282
reference_element::size_type size_type
Definition: geo_element.h:125
void put(std::ostream &is) const
Definition: geo_element.cc:93
size_type dimension() const
Definition: geo_element.h:167
void get(std::istream &os)
Definition: geo_element.cc:42
size_type subgeo_size(size_type subgeo_dim, size_type loc_isid) const
Definition: geo_element.h:221
orientation_type get_side_orientation(const geo_element &S) const
Definition: geo_element.cc:240
static size_type fix_indirect(const geo_element &K, size_type subgeo_variant, size_type loc_ige, size_type loc_comp_idof_on_subgeo, size_type order)
Definition: geo_element.cc:320
variant_type variant() const
Definition: geo_element.h:161
static void loc_qua_inod2lattice(size_type loc_qua_inod, size_type order, point_basic< size_type > &ij_lattice)
Definition: geo_element.cc:294
orientation_type get_edge_orientation(size_type dis_iv0, size_type dis_iv1) const
Definition: geo_element.cc:155
static void loc_tri_inod2lattice(size_type loc_tri_inod, size_type order, point_basic< size_type > &ij_lattice)
Definition: geo_element.cc:265
geo_element_indirect::shift_type shift_type
Definition: geo_element.h:135
orientation_type get_side_informations(const geo_element &S, size_type &loc_isid, size_type &shift) const
Definition: geo_element.cc:185
size_type n_subgeo(size_type subgeo_dim) const
Definition: geo_element.h:212
static size_type fix_quadrangle_indirect(const geo_element &K, size_type loc_iqua, size_type loc_iqua_j, size_type order)
Definition: geo_element.cc:311
virtual void reset(variant_type variant, size_type order)=0
char name() const
Definition: geo_element.h:169
reference_element::variant_type variant_type
Definition: geo_element.h:126
size_type order() const
Definition: geo_element.h:162
static const variant_type q
static const variant_type e
void set_variant(variant_type x)
static const variant_type max_variant
static const variant_type p
variant_type variant() const
static const variant_type t
#define error_macro(message)
Definition: dis_macros.h:49
#define fatal_macro(message)
Definition: dis_macros.h:33
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.
char skip_blancs_and_tabs(std::istream &is)
Definition: geo_element.cc:32
geo_element_indirect::orientation_type orient