Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
geo_element.cc
Go to the documentation of this file.
1
21
22#include "rheolef/geo_element.h"
23
24#include "geo_element_aux.icc"
25
26namespace rheolef {
27
28// --------------------------------------------------------------------------
29// i/o
30// --------------------------------------------------------------------------
31char
32skip_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}
41void
42geo_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}
92void
93geo_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
113bool
115 const geo_element& S,
116 orientation_type& orient,
117 shift_type& shift) const
118{
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
162void
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
172void
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}
229void
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,
257 size_type order)
258{
259 return geo_element_fix_edge_indirect (K, loc_iedg, loc_iedg_j, order);
260}
261// --------------
262// 2) triangles
263// --------------
264void
266 size_type loc_tri_inod,
267 size_type order,
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,
276 size_type order,
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,
286 size_type order)
287{
288 return geo_element_fix_triangle_indirect (K, loc_ifac, loc_itri_j, order);
289}
290// --------------
291// 3) quadrangles
292// --------------
293void
295 size_type loc_qua_inod,
296 size_type order,
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,
305 size_type order,
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,
315 size_type order)
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,
325 size_type order)
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
geo_element_indirect::orientation_type orientation_type
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
size_type & operator[](size_type loc_inod)
size_type n_node() const
size_type subgeo_local_vertex(size_type subgeo_dim, size_type i_subgeo, size_type i_subgeo_vertex) const
static size_type fix_edge_indirect(const geo_element &K, size_type loc_iedg, size_type loc_iedg_j, size_type order)
size_type size() const
static size_type fix_triangle_indirect(const geo_element &K, size_type loc_itri, size_type loc_itri_j, size_type order)
reference_element::size_type size_type
void put(std::ostream &is) const
size_type dimension() const
void get(std::istream &os)
size_type subgeo_size(size_type subgeo_dim, size_type loc_isid) const
orientation_type get_side_orientation(const geo_element &S) const
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)
variant_type variant() const
static void loc_qua_inod2lattice(size_type loc_qua_inod, size_type order, point_basic< size_type > &ij_lattice)
orientation_type get_edge_orientation(size_type dis_iv0, size_type dis_iv1) const
static void loc_tri_inod2lattice(size_type loc_tri_inod, size_type order, point_basic< size_type > &ij_lattice)
geo_element_indirect::shift_type shift_type
orientation_type get_side_informations(const geo_element &S, size_type &loc_isid, size_type &shift) const
size_type n_subgeo(size_type subgeo_dim) const
static size_type fix_quadrangle_indirect(const geo_element &K, size_type loc_iqua, size_type loc_iqua_j, size_type order)
virtual void reset(variant_type variant, size_type order)=0
reference_element::variant_type variant_type
size_type order() const
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)
geo_element_indirect::orientation_type orient