Flow123d  master-f44eb46
ref_element.hh
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file ref_element.hh
15  * @brief Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc.
16  * @author Jan Stebel
17  * @todo
18  *
19  * TODO: reconsider following whether it is actual...
20  * - design interface in such a way, that we can change numbering
21  * - design numbering and orientations on ref element that is consistent (orientation and numbering od 2d el. match sides of 3d),
22  * and possibly allows permute vertices of elements so that sides sharing an edge match numbering and orientation (we dos'nt need permute faces)
23  *
24  * Proposal(prefers combinatoric order) :
25  * 1D - orientation V0 -> V1
26  *
27  * 2D - edges: E0: V0 -> V1,
28  * E1: V0 -> V2
29  * E2: V1 -> V2
30  * This maximize number of edge orientations achievable by edge permutations
31  * edge numbering edge orientation( in original numbering)
32  * 0 1 2 + + +
33  * 0 2 1 - + +
34  * 1 0 2 + + -
35  * 1 2 0 - - +
36  * 2 0 1 + - -
37  * 2 1 0 - - -
38  *
39  * vertices edges normal (out = +)
40  * 3D - sides: S0: 0 1 2 E0 E1 E3 -
41  * S1: 0 1 3 E0 E2 E4 +
42  * S2: 0 2 3 E1 E2 E5 -
43  * S3: 1 2 3 E3 E4 E5 -
44  *
45  * edges: E0: V0 -> V1 x direction
46  * E1: V0 -> V2 y direction
47  * E2: V0 -> V3 z direction
48  * E3: V1 -> V2
49  * E4: V1 -> V3
50  * E5: V2 -> V3
51 
52  *
53  * - functions from DEAL.ii:
54  * bool is_inside_unit_cell( point )
55  * line_to_cell_vertices(line, vertex) vertex index on line to index on whole element
56  * face_to_cell_vertices(face, vertex, Orientation), Orientation should be some class describing permutation of shared face to element's face/side
57  * face_to_cel_lines
58  * standard_to_real_face_vertex(vertex, Orientation), maps vertex to permuted face
59  * real_to_standard_face_vertex - inverse
60  * ... same for line; we should need also something like standard_to_real_line_vertex; Deal dosn;t support Orientation changes for 2D element faces
61  * Point unit_cell_vertex(vertex) - coordinates
62  * project_to_unit_cell
63  * distance_to_unit_cell
64  * d_linear_shape_function
65  *
66  * - can not change numbering of element sides due to DarcyFlow, which use hardwired side numbering in construction of basis functions
67  * - any change of side numbering requires also change in flow/old_bcd.cc
68  *
69  *
70  */
71 
72 #ifndef REF_ELEMENT_HH_
73 #define REF_ELEMENT_HH_
74 
75 #include <vector> // for vector
76 #include <array>
77 #include <armadillo>
78 #include "system/armor.hh"
79 #include "system/asserts.hh"
80 
81 
82 /*
83  * Ordering of nodes and sides in reference elements
84  * =================================================
85  *
86  * TODO we want the following (22.10.):
87  *
88  * 1D element (line segment) 2D element (triangle) 3D element (tetrahedron)
89  *
90  * z
91  * .
92  * ,/
93  * /
94  * 3
95  * y ,/|`\
96  * ^ ,/ | `\
97  * | ,/ '. `\
98  * 2 ,/ | `\
99  * |`\ ,/ | `\
100  * | `\ 0-----------'.--------1 --> x
101  * | `\ `\. | ,/
102  * | `\ `\. | ,/
103  * | `\ `\. '. ,/
104  * 0----------1 --> x 0----------1 --> x `\. |/
105  * `2
106  * `\.
107  * `y
108  *
109  * side id node ids side id node ids side id node ids normal
110  * 0 0 0 0,1 0 0,1,2 OUT
111  * 1 1 1 0,2 1 0,1,3 IN
112  * 2 1,2 2 0,2,3 OUT
113  * 3 1,2,3 IN
114  *
115  *
116  * nodes coordinates:
117  * 0 [0] 0 [0,0] 0 [0,0,0]
118  * 1 [1] 1 [1,0] 1 [1,0,0]
119  * 2 [0,1] 2 [0,1,0]
120  * 3 [0,0,1]
121  *
122  * barycentric coordinates of nodes:
123  * 0 [1,0] 0 [1,0,0] 0 [1,0,0,0]
124  * 1 [0,1] 1 [0,1,0] 1 [0,1,0,0]
125  * 2 [0,0,1] 2 [0,0,1,0]
126  * 3 [0,0,0,1]
127  *
128  *
129  *
130  * Element node permutation for matching sides
131  *
132  * 1. 2D case: Can add mesh elements in "layers" so that at most toe edges are prescribed and one is free.
133  * possible cases, two edges A and B considered counter clockwise around a node:
134  * A out, B out : vertex 0 to node, regular orientation (normal up)
135  * A out, B in : vertex 1 to node, regular orientation
136  * A in , B out : vertex 1 to node, inverted orientation (normal down)
137  * A in , B in : vertex 2 to node regular orientation
138  *
139  * Result: can be matched using two element orientations, element inversion in 1/3 of cases
140  *
141  * 2. 3d elements around an oriented edge
142  * a. possible face configurations
143  * U: side 0 attached to the edge, positive; face vertices: 0 (edge 0), 1 (edge 1), 2
144  * V: side 1 attached to the edge, negative; face vertices: 0 (edge 0), 2 (edge 1), 1
145  * W: side 2 attached to the edge, positive; face vertices: 1 (edge 0), 2 (edge 1), 0
146  *
147  * b. faces of the element edges, configuration of faces attached to the element edge; edge pointing up, looking from outside position, face A on right, face B on left
148  * E0: A (side 0) in U, B (side 1) in U
149  * E1: A (side 2) in U, B (side 0) in V
150  * E2: A (side 1) in V, B (side 2) in V
151  * E3. A (side 0) in W, B (side 3) in U
152  * E4: A (side 3) in V, B (side 1) in W
153  * E5: A (side 2) in W, B (side 3) in W
154  * not presented combination:
155  * A in U, B in W = E3 inverted
156  * A in V, B in U = E1 inverted
157  * A in W, b in V = E4 inverted
158  *
159  * Result can permute element veritices to arrange elements around an edge, possibly element inversion in 1/3 of cases
160  *
161  * 3. 3d elements around a node, up to 3 given faces, one free face
162  * denote V0, V1, V2 vertices of the free face, V3 common vertex of given faces
163  * Denote edges of this element ABCDEF with +- orientation
164  * a. orientation of edges of the vertices
165  * V0: E0+ >V1, E1+ >V2, E2+ >V3 +++
166  * V1: E0- >V0, E3+ >V2, E4+ >V3 -++
167  * V2: E1- >V0, E3- >V1, E5+ >V3 --+
168  * V3: E2- >V0, E4- >V1, E5- >V2 ---
169  *
170 * vertices edges normal (out = +)
171  * 3D - sides: S0: 0 1 2 E0 E1 E3 -
172  * S1: 0 1 3 E0 E2 E4 +
173  * S2: 0 2 3 E1 E2 E5 -
174  * S3: 1 2 3 E3 E4 E5 -
175  *
176  * edges: A E0: V0 -> V1 x direction
177  * B E1: V0 -> V2 y direction
178  * C E2: V0 -> V3 z direction
179  * D E3: V1 -> V2
180  * E E4: V1 -> V3
181  * F E5: V2 -> V3
182 
183  CE -> A
184  EF -> D
185  FC -> B
186  *
187  * with edges DEF, other three edges ABC
188  * configuration given by orientation of edges ABC and possibly by orientation of edges DEF
189  * ABC orientation + down, - up; DEF orientation positive if match side 0, i.e. D (AB), E (AC), F (BC)
190  * faces denoted fD,fE,fF, their configurations considered with respect to the edges DEF respectively
191  * front face denoted: ff
192  * vertices V0, V1, V2 top vertices with edges D (V0V1) E (V0V2) F (V1V2); V3 bottom
193  * invalid = corner case, can not continue must either modify one of the faces, or introduce unmatching face
194  * impossible = combination of edges that can not happen
195  * edge cases arranged for the vertex 3
196  *
197  * CEF ABD face config: fD fE fF ff
198  * +++ +++ U U U
199  * +++ ++- U U U-
200  * +++ +-+ invalid
201  * +++ +-- U U U
202  * +++ -++ U U U
203  * +++ -+- invalid U U U-
204  * +++ --+
205  * +++ --- U U U
206  *
207  * ++- +++ impossible U U U
208  * ++- ++- impossible U U U-
209  * ++- +-+ impossible
210  * ++- +-- U U U
211  * ++- -++ impossible U U U
212  * ++- -+- impossible U U U-
213  * ++- --+ impossible
214  * ++- --- U U U
215 
216  * +-+ +++ impossible U U U
217  * +-+ ++- impossible U U U-
218  * +-+ +-+ impossible
219  * +-+ +-- U U U
220  * +-+ -++ impossible U U U
221  * +-+ -+- impossible U U U-
222  * +-+ --+ impossible
223  * +-+ --- U U U
224 
225  * +-- +++ impossible U U U
226  * +-- ++- impossible U U U-
227  * +-- +-+ impossible
228  * +-- +-- U U U
229  * +-- -++ impossible U U U
230  * +-- -+- impossible U U U-
231  * +-- --+ impossible
232  * +-- --- U U U
233 
234  * -++ +++ impossible U U U
235  * -++ ++- impossible U U U-
236  * -++ +-+ impossible
237  * -++ +-- U U U
238  * -++ -++ impossible U U U
239  * -++ -+- impossible U U U-
240  * -++ --+ impossible
241  * -++ --- U U U
242 
243  * -+- +++ impossible U U U
244  * -+- ++- impossible U U U-
245  * -+- +-+ impossible
246  * -+- +-- U U U
247  * -+- -++ impossible U U U
248  * -+- -+- impossible U U U-
249  * -+- --+ impossible
250  * -+- --- U U U
251 
252  * --+ +++ impossible U U U
253  * --+ ++- impossible U U U-
254  * --+ +-+ impossible
255  * --+ +-- U U U
256  * --+ -++ impossible U U U
257  * --+ -+- impossible U U U-
258  * --+ --+ impossible
259  * --+ --- U U U
260 
261  * --- +++ U U U
262  * --- ++- U U U-
263  * --- +-+ invalid
264  * --- +-- U U U
265  * --- -++ U U U
266  * --- -+- invalid U U U-
267  * --- --+
268  * --- --- U U U
269  *
270  *
271  *
272  */
273 
274 
275 template<std::size_t Size>
276 using IdxVector = std::array<unsigned int, Size>;
277 
278 
279 /** Auxilliary structure that is used to pass template arguments into interact function of RefElement:
280  * RefElement<dim>::interact( Interaction<OutDim,InDim>(i) )
281  *
282  * This enables automatic deduction of dimensional template arguments.
283  * @see @p RefElement<dim>::interact
284  */
285 template <unsigned int OutDim, unsigned int InDim>
286 struct Interaction {
287  Interaction(unsigned int i) : i_(i) {}
288  unsigned int i_;
289 };
290 
291 
292 
294 public:
295  // Order clockwise, faces opposite to the lines from node_lines.
296  // !! dependes on S3 inversion
297  static constexpr IdxVector<3> S3_node_sides [2][4]
298  = { { { 2, 1, 0 },
299  { 3, 0, 1 },
300  { 3, 2, 0 },
301  { 3, 1, 2 }},
302  { { 2, 0, 1 },
303  { 3, 1, 0 },
304  { 3, 0, 2 },
305  { 3, 2, 1 }}};
306 
307  // faces adjecent to given edge, first is the right face when looking form outside with the edge pointing up.
308  // !! dependes on S3 inversion
309  static constexpr IdxVector<2> S3_line_sides [2][6]
310  = { { {0,1},
311  {2,0},
312  {1,2},
313  {0,3},
314  {3,1},
315  {2,3}},
316  { {1,0},
317  {0,2},
318  {2,1},
319  {3,0},
320  {1,3},
321  {3,2}}};
322 
323  // Order clockwise looking over the vertex to center; smallest index first
324  // !! dependes on S3 inversion
325  static constexpr IdxVector<3> S3_node_lines [2][4]
326  = { { {0,1,2},
327  {0,4,3},
328  {1,3,5},
329  {2,5,4}},
330  { {0,2,1},
331  {0,3,4},
332  {1,5,3},
333  {2,4,5}}};
334 
335 };
336 
337 
338 template<unsigned int dim>
340 {
341 public:
342  typedef arma::vec::fixed<dim> LocalPoint;
343  /**
344  * Barycentric coordinates.
345  *
346  * e.g. coordinates (a,b,c) on triangle with vertices X, Y, Z
347  * represents a point: a*X+b*Y+c*Z
348  */
351 
352  DECLARE_EXCEPTION( ExcInvalidPermutation, << "Side permutation not found.\n" );
353 
354  /**
355  * Return coordinates of given node.
356  * @see the class documentation @p RefElement
357  * @param nid Node number.
358  * NOTE: Implementation is dependent on current node and side numbering.
359  */
360  static LocalPoint node_coords(unsigned int nid);
361 
362  /**
363  * Compute normal vector to a given side.
364  * @param sid Side number.
365  */
366  static LocalPoint normal_vector(unsigned int sid);
367 
368 
369  /**
370  * If the given barycentric coordinate is in the ref. element, return unchanged.
371  * If the given barycentric coordinate is out of the ref. element,
372  * project it on the surface of the ref. element.
373  */
374  static BaryPoint clip(const BaryPoint &barycentric);
375 
376  /** Returns orientation of the normal of side @p sid. 0 -> OUT, 1 -> IN.
377  * NOTE: Implementation is dependent on current node and side numbering.
378  */
379  static unsigned int normal_orientation(unsigned int sid);
380 
381  static double side_measure(unsigned int sid);
382 
383  /**
384  * Returns index of the node that is oposite to side of given index @p sid.
385  * Note: It is dependent on current node and side numbering.
386  * @param sid Side number.
387  * NOTE: Implementation is dependent on current node and side numbering.
388  */
389  static unsigned int oposite_node(unsigned int sid);
390 
391  /**
392  * Return index of 1D line, shared by two faces @p f1 and @p f2 of the reference tetrahedron.
393  * Implemented only for @p dim == 3.
394  */
395  static unsigned int line_between_faces(unsigned int f1, unsigned int f2);
396 
397 
398  static const unsigned int n_sides = dim + 1; ///< Number of sides.
399  static const unsigned int n_nodes = dim + 1; ///< Number of nodes.
400  static const unsigned int n_nodes_per_side = dim; ///< Number of nodes on one side.
401  static const unsigned int n_lines_per_node = dim; ///< Number of lines with one common node.
402  static const unsigned int n_nodes_per_line = 2; ///< Number of nodes in one line.
403  static const unsigned int n_sides_per_line = 2; ///< Number of sides with one common line. @p dim == 3.
404  static const unsigned int n_sides_per_node = dim; ///< Number of sides with one common line.
405 
406  /// Number of lines on boundary of one side.
407  static const unsigned int n_lines_per_side = (unsigned int)((dim * (dim - 1)) / 2);//( dim == 3 ? 3 : 0);// Kombinační číslo dim nad dvěma
408 
409  /// Number of lines, i.e. @p object of dimension @p dim-2 on the boundary of the reference element.
410  static const unsigned int n_lines = (unsigned int)((dim * (dim + 1)) / 2); //( dim == 3 ? 6 : dim == 2 ? 3 : dim == 1 ? 1 : 0); součet posloupnosti
411 
412 
414 
415 
416  /** @brief Converts from local to barycentric coordinates.
417  * @param lp point in local coordinates (x,y)
418  * @return point in barycentric coordinates (1-x-y, x, y)
419  */
420  static BaryPoint local_to_bary(const LocalPoint& lp);
421 
422  /** @brief Converts from barycentric to local coordinates.
423  * @param bp point in barycentric coordinates
424  * @return point in local coordinates
425  */
426  static LocalPoint bary_to_local(const BaryPoint& bp);
427 
429 
430  /**
431  * Used in the clip method.
432  */
434 
435  /**
436  * For given barycentric coordinates on the ref element returns barycentric coordinates
437  * on the ref. element of given face. Assumes that the input point is on the face.
438  * Barycentric order: (complanatory, local_coords )
439  */
440  static FaceBaryPoint barycentric_on_face(const BaryPoint &barycentric, unsigned int i_face);
441 
442 
444  static CentersList centers_of_subelements(unsigned int sub_dim);
445 
446  /**
447  * Return (1) number of zeros and (2) positions of zeros in barycentric coordinates.
448  * @p tolerance serves for testing zero values of @p barycentric coordinates.
449  */
450  static std::pair<unsigned int, unsigned int> zeros_positions(const BaryPoint &barycentric,
451  double tolerance = std::numeric_limits<double>::epsilon()*2);
452 
453  /**
454  * According to positions of zeros in barycentric coordinates, it gives the index of subdim-simplex
455  * in the reference element. Number of zeros must be equal to (3-subdim).
456  * e.g.:
457  * if 1 zeros, return index of side (subdim 2)
458  * if 2 zeros, return index of edge (subdim 1)
459  * if 3 zeros, return index of vertex (subdim 0)
460  */
461  template<unsigned int subdim> static unsigned int topology_idx(unsigned int zeros_positions);
462 
463  /** Function returns number of subdim-simplices inside dim-simplex.
464  * The aim is covering all the n_**** members with a single function.
465  * TODO: think of generalization for n_****_per_**** members, like function @p interact:
466  * template<unsigned int subdimA, unsigned int subdimB> static unsigned int count();
467  */
468  template<unsigned int subdim> static unsigned int count();
469 
470  /**
471  * @param sid - index of a sub-simplex in a simplex
472  * return an array of barycentric coordinates on <dim> simplex from <subdim> simplex
473  * for example: simplex<3> - ABCD and its subsubsimplex<1> AD (line index: 3)
474  * AD has barycoords for A (1,0), for D (0,1), but A in ABCD is (1,0,0,0) and D is (0,0,0,1)
475  * this method creates array ((1,0,0,0),(0,0,0,1))
476  */
477  template<unsigned int subdim> static arma::mat::fixed<dim+1,subdim+1> bary_coords(unsigned int sid);
478 
479  /** Interpolate barycentric coords to a higher dimension of a simplex.
480  * @param coord - barycentric coords of a point on a sub-simplex
481  * @param sub_simplex_idx - id of sub-simplex on a simplex
482  */
483  template<unsigned int subdim> static BaryPoint interpolate(arma::vec::fixed<subdim+1> coord, int sub_simplex_idx);
484 
485 
486  /**
487  * Basic line interpolation.
488  */
490  BaryPoint second_coords,
491  double first_theta, double second_theta, double theta);
492 
493  /**
494  * Usage:
495  * RefElement<3>::interact(Interaction<2,0>(1))
496  * (means: In tetrahedron <3>, give indices of sides <2>, connected by node <0> with index 1)
497  * RefElement<3>::interact(Interaction<2,0>(1))[1]
498  * (as above, but give only the side with index 1)
499  *
500  * Template usage: RefElement<dim>::interact(Interaction<OutDim, InDim>(i))[j]
501  * (means: on dim-dimensional reference element, go on InDim-dimensional subelement with index i,
502  * which connects OutDim-dimnesional subelements and select the one with index j)
503  *
504  * This method serves as an interface to topology information of the reference element.
505  * It returns indices of OutDim-dimensional object
506  * of InDim-dimnesional object of given index
507  * in dim-dimnesional reference element.
508  * @tparam interaction - auxilliary object carying the index and the template arguments OutDim and InDim
509  * @tparam OutDim - output dimension (give me node-0, line-1, side-2), <= dim
510  * @tparam InDim - input dimension (for node-0, line-1, side-2), <= dim
511  * @return vector of indices of OutDim-dimensional subelements represented by @p IdxVector object.
512  *
513  * possible calls:
514  * dim OutDim InDim return
515  * 1,2,3 0 1 InDim+1 - give me indices of nodes of line of given index
516  * 3 0 2 InDim+1 - give me indices of nodes of a side (triangle) of given index
517  * 3 1 2 InDim+1 - give me indices of lines of side (triangle) of given index
518  *
519  * 1,2,3 1 0 dim-InDim - give me indices of lines with common node of given index
520  * 3 2 0 dim-InDim - give me indices of sides (triangles) with common node of given index
521  * 3 2 1 dim-InDim - give me indices of sides (triangles) with common line of given index
522  *
523  */
524  template < template <unsigned int OutDim, unsigned int InDim> class TInteraction, unsigned int OutDim, unsigned int InDim>
525  static const IdxVector< (InDim>OutDim ? InDim+1 : dim-InDim) > interact( TInteraction<OutDim,InDim> interaction, bool inv = false );
526 
527 
528 private:
529  /// Internal part of the interact function.
530  template<unsigned int OutDim, unsigned int InDim>
531  static const IdxVector< (InDim>OutDim ? InDim+1 : dim-InDim) > interact_(unsigned int index, bool inv = false);
532 
533  static const std::vector<IdxVector<n_nodes_per_line>> line_nodes_; ///< [n_lines] For given line, returns its nodes indices.
534  static const std::vector<IdxVector<n_lines_per_node>> node_lines_; ///< [n_nodes] For given node, returns lines indices.
535  static const std::vector<IdxVector<n_nodes_per_side>> side_nodes_; ///< [n_sides] For given side, returns nodes indices. For @p dim == 3.
536  static const std::vector<IdxVector<n_sides_per_node>> node_sides_; ///< [n_nodes] For given node, returns sides indices. For @p dim == 3.
537  static const std::vector<IdxVector<n_sides_per_line>> line_sides_; ///< [n_lines] For given line, returns sides indices. For @p dim == 3.
538  static const std::vector<IdxVector<n_lines_per_side>> side_lines_; ///< [n_sides] For given side, returns lines indices. For @p dim == 3.
539 
540  //TODO: implement for 1d and 2d
541  /**
542  * Consider an n-face (node, edge, face, bulk) with dimension `subdim` and
543  * index within subdimension `idx`. Barycentric coordinates of all points
544  * on the n-face have unique pattern of zero coordinates.
545  *
546  * topology_zeros_[subdim][idx] is a bitfield with '1' where the pattern have zeros.
547  */
548  static const IdxVector<(n_lines > n_nodes) ? n_lines : n_nodes> topology_zeros_[dim+1];
549 };
550 
551 
552 // Declarations of explicit specialization of static memebers.
553 
558 
559 
565 
571 
572 
577 
578 
579 
580 
581 
582 
583 // 0: nodes of nodes
584 // 1: nodes of lines
585 // 2: nodes of sides
586 // 3: nodes of tetrahedron
591 
592 
593 
594 
595 
596 
597 
598 
599 
600 /************************* template implementation ****************************/
601 
602 template<unsigned int dim>
603 template<unsigned int subdim> inline
604 arma::mat::fixed<dim+1,subdim+1> RefElement<dim>::bary_coords(unsigned int sid){
605  ASSERT_LT(subdim, dim).error("Dimension mismatch!");
606  arma::mat::fixed<dim+1,subdim+1> bary_c;
607 
608  for(unsigned int i = 0; i < subdim+1; i++){
609  unsigned int nid = interact_<0,subdim>(sid)[i];
610  bary_c.col(i).zeros();
611  bary_c.col(i)(nid) = 1;
612  }
613 
614  return bary_c;
615 }
616 
617 
618 template<unsigned int dim> inline
619 arma::vec::fixed<dim> RefElement<dim>::node_coords(unsigned int nid)
620 {
621  ASSERT_LT(nid, n_nodes).error("Node number is out of range!");
622 
623  arma::vec::fixed<dim> p;
624  p.zeros();
625 
626  if (nid > 0)
627  p(nid-1) = 1;
628 
629  return p;
630 }
631 
632 
633 template<unsigned int dim>
634 template<unsigned int subdim>
635 auto RefElement<dim>::interpolate(arma::vec::fixed<subdim+1> coord, int sub_simplex_idx) -> BaryPoint
636 {
637  return RefElement<dim>::bary_coords<subdim>(sub_simplex_idx)*coord;
638 }
639 
640 
641 template<> template<> inline unsigned int RefElement<3>::count<0>()
642 { return n_nodes; }
643 template<> template<> inline unsigned int RefElement<3>::count<1>()
644 { return n_lines; }
645 template<> template<> inline unsigned int RefElement<3>::count<2>()
646 { return n_sides; }
647 template<> template<> inline unsigned int RefElement<3>::count<3>()
648 { return 1; }
649 template<> template<> inline unsigned int RefElement<2>::count<0>()
650 { return n_nodes; }
651 template<> template<> inline unsigned int RefElement<2>::count<1>()
652 { return n_lines; }
653 template<> template<> inline unsigned int RefElement<2>::count<2>()
654 { return 1; }
655 template<> template<> inline unsigned int RefElement<2>::count<3>()
656 { return 0; }
657 template<> template<> inline unsigned int RefElement<1>::count<0>()
658 { return n_nodes; }
659 template<> template<> inline unsigned int RefElement<1>::count<1>()
660 { return 1; }
661 template<> template<> inline unsigned int RefElement<1>::count<2>()
662 { return 0; }
663 template<> template<> inline unsigned int RefElement<1>::count<3>()
664 { return 0; }
665 template<> template<> inline unsigned int RefElement<0>::count<0>()
666 { return 1; }
667 template<> template<> inline unsigned int RefElement<0>::count<1>()
668 { return 0; }
669 template<> template<> inline unsigned int RefElement<0>::count<2>()
670 { return 0; }
671 template<> template<> inline unsigned int RefElement<0>::count<3>()
672 { return 0; }
673 
674 template<unsigned int dim>
675 template<unsigned int subdim>
676 unsigned int RefElement<dim>::topology_idx(unsigned int zeros_positions)
677 {
678  for(unsigned int i=0; i < RefElement<dim>::count<subdim>(); i++){
679  if(zeros_positions == topology_zeros_[subdim][i]) return i;
680  }
681  ASSERT_PERMANENT(0).error("Undefined zero pattern.");
682  return -1;
683 }
684 
685 
686 /// This function is for "side_nodes" - for given side, give me nodes (0->0, 1->1).
687 template<> template<> inline const IdxVector<1> RefElement<1>::interact_<0,0>(unsigned int i, bool inv)
688 { ASSERT_LT(i, RefElement<1>::n_nodes)(inv).error("Index out of bounds.");
689  return IdxVector<1>({i});}
690 
691 /// For line i {0}, give me indices of its nodes.
692 template<> template<> inline const IdxVector<2> RefElement<1>::interact_<0,1>(unsigned int i, bool inv)
693 { ASSERT_LT(i, RefElement<1>::n_lines)(inv).error("Index out of bounds.");
694  return line_nodes_[i];}
695 
696 /// For line i {0,1,2}, give me indices of its nodes.
697 template<> template<> inline const IdxVector<2> RefElement<2>::interact_<0,1>(unsigned int i, bool inv)
698 { ASSERT_LT(i, RefElement<2>::n_lines)(inv).error("Index out of bounds.");
699  return line_nodes_[i];}
700 
701 /// For line i {0,1,2,3,4,5}, give me indices of its nodes.
702 template<> template<> inline const IdxVector<2> RefElement<3>::interact_<0,1>(unsigned int i, bool inv)
703 { ASSERT_LT(i, RefElement<3>::n_lines)(inv).error("Index out of bounds.");
704  return line_nodes_[i];}
705 
706 /// For node i {0,1}, give me indices of lines.
707 template<> template<> inline const IdxVector<1> RefElement<1>::interact_<1,0>(unsigned int i, bool inv)
708 { ASSERT_LT(i, RefElement<1>::n_nodes)(inv).error("Index out of bounds.");
709  return node_lines_[i];}
710 
711 /// For node i {0,1,2}, give me indices of lines.
712 template<> template<> inline const IdxVector<2> RefElement<2>::interact_<1,0>(unsigned int i, bool inv)
713 { ASSERT_LT(i, RefElement<2>::n_nodes)(inv).error("Index out of bounds.");
714  return node_lines_[i];}
715 
716 /// For node i {0,1,2,3}, give me indices of lines.
717 template<> template<> inline const IdxVector<3> RefElement<3>::interact_<1,0>(unsigned int i, bool inv)
718 { ASSERT_LT(i, RefElement<3>::n_nodes).error("Index out of bounds.");
719  return _AuxInteract::S3_node_lines[inv][i];}
720 
721 /// For side i {0,1,2}, give me indices of its nodes.
722 template<> template<> inline const IdxVector<3> RefElement<3>::interact_<0,2>(unsigned int i, bool inv)
723 { ASSERT_LT(i, RefElement<3>::n_sides)(inv).error("Index out of bounds.");
724  return side_nodes_[i];}
725 
726 /// For node i {0,1,2,3}, give me indices of sides.
727 template<> template<> inline const IdxVector<3> RefElement<3>::interact_<2,0>(unsigned int i, bool inv)
728 { ASSERT_LT(i, RefElement<3>::n_sides).error("Index out of bounds.");
729  return _AuxInteract::S3_node_sides[inv][i];}
730 
731 /// For line i {0,1,2,3}, give me indices of sides.
732 template<> template<> inline const IdxVector<2> RefElement<3>::interact_<2,1>(unsigned int i, bool inv)
733 { ASSERT_LT(i, RefElement<3>::n_lines).error("Index out of bounds.");
734  return _AuxInteract::S3_line_sides[inv][i];}
735 
736 /// For side i {0,1,2}, give me indices of its lines.
737 template<> template<> inline const IdxVector<3> RefElement<3>::interact_<1,2>(unsigned int i, bool inv)
738 { ASSERT_LT(i, RefElement<3>::n_sides)(inv).error("Index out of bounds.");
739  return side_lines_[i];}
740 
741 template<unsigned int dim> template<unsigned int OutDim, unsigned int InDim>
742 inline const IdxVector< (InDim>OutDim ? InDim+1 : dim-InDim) > RefElement<dim>::interact_(unsigned int i, bool inv)
743 {
744  ASSERT_PERMANENT(false)(dim)(OutDim)(InDim)(i)(inv).error("Not implemented.");
745  //ASSERT_LT(OutDim, dim);
746  //ASSERT_LT(InDim, dim);
747  return IdxVector< (InDim>OutDim ? InDim+1 : dim-InDim) >(); // just to avoid warning for missing return
748 }
749 
750 
751 template<unsigned int dim>
752 template < template <unsigned int OutDim, unsigned int InDim> class TInteraction, unsigned int OutDim, unsigned int InDim>
753 inline const IdxVector< (InDim>OutDim ? InDim+1 : dim-InDim) > RefElement<dim>::interact( TInteraction<OutDim,InDim> interaction , bool inv)
754 {
755  return interact_<OutDim,InDim>(interaction.i_, inv);
756 }
757 
758 #endif /* REF_ELEMENT_HH_ */
Definitions of ASSERTS.
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:348
#define ASSERT_LT(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:301
static const std::vector< IdxVector< n_lines_per_side > > side_lines_
[n_sides] For given side, returns lines indices. For dim == 3.
Definition: ref_element.hh:538
static CentersList centers_of_subelements(unsigned int sub_dim)
Definition: ref_element.cc:366
static const std::vector< IdxVector< n_nodes_per_line > > line_nodes_
[n_lines] For given line, returns its nodes indices.
Definition: ref_element.hh:533
static LocalPoint bary_to_local(const BaryPoint &bp)
Converts from barycentric to local coordinates.
Definition: ref_element.cc:201
static unsigned int normal_orientation(unsigned int sid)
Definition: ref_element.cc:216
static const unsigned int n_lines_per_node
Number of lines with one common node.
Definition: ref_element.hh:401
Armor::ArmaVec< double, dim > FaceBaryPoint
Definition: ref_element.hh:350
Armor::ArmaVec< double, dim+1 > BaryPoint
Definition: ref_element.hh:349
static FaceBaryPoint barycentric_on_face(const BaryPoint &barycentric, unsigned int i_face)
Definition: ref_element.cc:279
const std::vector< LocalPoint > & CentersList
Definition: ref_element.hh:443
std::vector< BaryPoint > BarycentricUnitVec
Definition: ref_element.hh:428
static BaryPoint line_barycentric_interpolation(BaryPoint first_coords, BaryPoint second_coords, double first_theta, double second_theta, double theta)
Definition: ref_element.cc:437
static const unsigned int n_sides
Number of sides.
Definition: ref_element.hh:398
static const std::vector< IdxVector< n_sides_per_line > > line_sides_
[n_lines] For given line, returns sides indices. For dim == 3.
Definition: ref_element.hh:537
static const IdxVector<(n_lines > n_nodes) ? n_lines :n_nodes > topology_zeros_[dim+1]
Definition: ref_element.hh:548
static const std::vector< IdxVector< n_lines_per_node > > node_lines_
[n_nodes] For given node, returns lines indices.
Definition: ref_element.hh:534
arma::vec::fixed< dim > LocalPoint
Definition: ref_element.hh:342
static BaryPoint local_to_bary(const LocalPoint &lp)
Converts from local to barycentric coordinates.
Definition: ref_element.cc:187
static const unsigned int n_lines_per_side
Number of lines on boundary of one side.
Definition: ref_element.hh:407
static BarycentricUnitVec make_bary_unit_vec()
Definition: ref_element.cc:318
static std::pair< unsigned int, unsigned int > zeros_positions(const BaryPoint &barycentric, double tolerance=std::numeric_limits< double >::epsilon() *2)
Definition: ref_element.cc:294
static const unsigned int n_nodes
Number of nodes.
Definition: ref_element.hh:399
DECLARE_EXCEPTION(ExcInvalidPermutation,<< "Side permutation not found.\n")
static const unsigned int n_nodes_per_side
Number of nodes on one side.
Definition: ref_element.hh:400
static const IdxVector<(InDim >OutDim ? InDim+1 :dim-InDim) > interact_(unsigned int index, bool inv=false)
Internal part of the interact function.
static arma::mat::fixed< dim+1, subdim+1 > bary_coords(unsigned int sid)
Definition: ref_element.hh:604
static const unsigned int n_sides_per_node
Number of sides with one common line.
Definition: ref_element.hh:404
static const unsigned int n_nodes_per_line
Number of nodes in one line.
Definition: ref_element.hh:402
static unsigned int line_between_faces(unsigned int f1, unsigned int f2)
static unsigned int topology_idx(unsigned int zeros_positions)
Definition: ref_element.hh:676
static const std::vector< IdxVector< n_sides_per_node > > node_sides_
[n_nodes] For given node, returns sides indices. For dim == 3.
Definition: ref_element.hh:536
static const unsigned int n_sides_per_line
Number of sides with one common line. dim == 3.
Definition: ref_element.hh:403
static BaryPoint interpolate(arma::vec::fixed< subdim+1 > coord, int sub_simplex_idx)
static LocalPoint node_coords(unsigned int nid)
Definition: ref_element.hh:619
static BaryPoint clip(const BaryPoint &barycentric)
Definition: ref_element.cc:332
static double side_measure(unsigned int sid)
static const unsigned int n_lines
Number of lines, i.e. object of dimension dim-2 on the boundary of the reference element.
Definition: ref_element.hh:410
static unsigned int count()
static LocalPoint normal_vector(unsigned int sid)
static const IdxVector<(InDim >OutDim ? InDim+1 :dim-InDim) > interact(TInteraction< OutDim, InDim > interaction, bool inv=false)
static unsigned int oposite_node(unsigned int sid)
Definition: ref_element.cc:209
static const std::vector< std::vector< std::vector< unsigned int > > > nodes_of_subelements
Definition: ref_element.hh:413
static const std::vector< IdxVector< n_nodes_per_side > > side_nodes_
[n_sides] For given side, returns nodes indices. For dim == 3.
Definition: ref_element.hh:535
static constexpr IdxVector< 3 > S3_node_sides[2][4]
Definition: ref_element.hh:298
static constexpr IdxVector< 3 > S3_node_lines[2][4]
Definition: ref_element.hh:326
static constexpr IdxVector< 2 > S3_line_sides[2][6]
Definition: ref_element.hh:310
typename arma::Col< Type >::template fixed< nr > ArmaVec
Definition: armor.hh:505
std::array< unsigned int, Size > IdxVector
Definition: ref_element.hh:276
Interaction(unsigned int i)
Definition: ref_element.hh:287
unsigned int i_
Definition: ref_element.hh:288