Flow123d  JS_before_hm-1754-g1847fd3ed
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 /** Auxilliary class representing vector of indices (unsigned int).
275  * @tparam Size is the fixed size of the vector.
276  */
277 /*
278 template<unsigned int Size>
279 class IdxVector{
280  unsigned int data_[Size]; ///< Array with indices.
281 
282  public:
283  /// Constructor taking in array of indices.
284  IdxVector(std::array<unsigned int,Size> data_in);
285  /// Constructor enabling creating object with initializer list {...}.
286  IdxVector(std::initializer_list<unsigned int> data_in);
287  /// Getter for index @p idx.
288  unsigned int operator[](unsigned int idx) const;
289 };
290 */
291 
292 template<std::size_t Size>
293 using IdxVector = std::array<unsigned int, Size>;
294 
295 
296 /** Auxilliary structure that is used to pass template arguments into interact function of RefElement:
297  * RefElement<dim>::interact( Interaction<OutDim,InDim>(i) )
298  *
299  * This enables automatic deduction of dimensional template arguments.
300  * @see @p RefElement<dim>::interact
301  */
302 template <unsigned int OutDim, unsigned int InDim>
303 struct Interaction {
304  Interaction(unsigned int i) : i_(i) {}
305  unsigned int i_;
306 };
307 
308 template<unsigned int dim>
309 class RefElement
310 {
311 public:
312  typedef arma::vec::fixed<dim> LocalPoint;
313  /**
314  * Barycentric coordinates.
315  *
316  * e.g. coordinates (a,b,c) on triangle with vertices X, Y, Z
317  * represents a point: a*X+b*Y+c*Z
318  */
319  typedef Armor::ArmaVec<double, dim+1> BaryPoint;
320  typedef Armor::ArmaVec<double, dim> FaceBaryPoint;
321 
322  DECLARE_EXCEPTION( ExcInvalidPermutation, << "Side permutation not found.\n" );
323 
324  /**
325  * Return coordinates of given node.
326  * @see the class documentation @p RefElement
327  * @param nid Node number.
328  * NOTE: Implementation is dependent on current node and side numbering.
329  */
330  static LocalPoint node_coords(unsigned int nid);
331 
332  /**
333  * Compute normal vector to a given side.
334  * @param sid Side number.
335  */
336  static LocalPoint normal_vector(unsigned int sid);
337 
338 
339  /**
340  * If the given barycentric coordinate is in the ref. element, return unchanged.
341  * If the given barycentric coordinate is out of the ref. element,
342  * project it on the surface of the ref. element.
343  */
344  static BaryPoint clip(const BaryPoint &barycentric);
345 
346  /** Returns orientation of the normal of side @p sid. 0 -> OUT, 1 -> IN.
347  * NOTE: Implementation is dependent on current node and side numbering.
348  */
349  static unsigned int normal_orientation(unsigned int sid);
350 
351  static double side_measure(unsigned int sid);
352 
353  /**
354  * Returns index of the node that is oposite to side of given index @p sid.
355  * Note: It is dependent on current node and side numbering.
356  * @param sid Side number.
357  * NOTE: Implementation is dependent on current node and side numbering.
358  */
359  static unsigned int oposite_node(unsigned int sid);
360 
361  /**
362  * Return index of 1D line, shared by two faces @p f1 and @p f2 of the reference tetrahedron.
363  * Implemented only for @p dim == 3.
364  */
365  static unsigned int line_between_faces(unsigned int f1, unsigned int f2);
366 
367 
368  static const unsigned int n_sides = dim + 1; ///< Number of sides.
369  static const unsigned int n_nodes = dim + 1; ///< Number of nodes.
370  static const unsigned int n_nodes_per_side = dim; ///< Number of nodes on one side.
371  static const unsigned int n_lines_per_node = dim; ///< Number of lines with one common node.
372  static const unsigned int n_nodes_per_line = 2; ///< Number of nodes in one line.
373  static const unsigned int n_sides_per_line = 2; ///< Number of sides with one common line. @p dim == 3.
374  static const unsigned int n_sides_per_node = dim; ///< Number of sides with one common line.
375 
376  /// Number of lines on boundary of one side.
377  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
378 
379  /// Number of lines, i.e. @p object of dimension @p dim-2 on the boundary of the reference element.
380  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
381 
382 
383 // /**
384 // * Node numbers for each side.
385 // */
386 // static const unsigned int side_nodes[n_sides][n_nodes_per_side];
387 //
388 // /**
389 // * Indices of 1D lines of the 2D sides of an tetrahedron. Nonempty only for @p dim==3.
390 // */
391 // static const unsigned int side_lines[n_sides][n_lines_per_side];
392 //
393 // /**
394 // * Nodes of 1D lines of the tetrahedron.
395 // */
396 // static const unsigned int line_nodes[n_lines][2];
397 //
398 // /**
399 // * Indices of sides for each line. Nonempty only for @p dim==3 and @p dim==2.
400 // */
401 // static const unsigned int line_sides[n_lines][2];
402 
403 
404  static const std::vector< std::vector< std::vector<unsigned int> > > nodes_of_subelements;
405 
406 
407  /** @brief Converts from local to barycentric coordinates.
408  * @param lp point in local coordinates (x,y)
409  * @return point in barycentric coordinates (1-x-y, x, y)
410  */
411  static BaryPoint local_to_bary(const LocalPoint& lp);
412 
413  /** @brief Converts from barycentric to local coordinates.
414  * @param bp point in barycentric coordinates
415  * @return point in local coordinates
416  */
417  static LocalPoint bary_to_local(const BaryPoint& bp);
418 
419  typedef std::vector<BaryPoint> BarycentricUnitVec;
420 
421  /**
422  * Used in the clip method.
423  */
424  static BarycentricUnitVec make_bary_unit_vec();
425 
426  /**
427  * For given barycentric coordinates on the ref element returns barycentric coordinates
428  * on the ref. element of given face. Assumes that the input point is on the face.
429  * Barycentric order: (complanatory, local_coords )
430  */
431  static FaceBaryPoint barycentric_on_face(const BaryPoint &barycentric, unsigned int i_face);
432 
433 
434  typedef const std::vector<LocalPoint> & CentersList;
435  static CentersList centers_of_subelements(unsigned int sub_dim);
436 
437  /**
438  * Return (1) number of zeros and (2) positions of zeros in barycentric coordinates.
439  * @p tolerance serves for testing zero values of @p barycentric coordinates.
440  */
441  static std::pair<unsigned int, unsigned int> zeros_positions(const BaryPoint &barycentric,
442  double tolerance = std::numeric_limits<double>::epsilon()*2);
443 
444  /**
445  * According to positions of zeros in barycentric coordinates, it gives the index of subdim-simplex
446  * in the reference element. Number of zeros must be equal to (3-subdim).
447  * e.g.:
448  * if 1 zeros, return index of side (subdim 2)
449  * if 2 zeros, return index of edge (subdim 1)
450  * if 3 zeros, return index of vertex (subdim 0)
451  */
452  template<unsigned int subdim> static unsigned int topology_idx(unsigned int zeros_positions);
453 
454  /** Function returns number of subdim-simplices inside dim-simplex.
455  * The aim is covering all the n_**** members with a single function.
456  * TODO: think of generalization for n_****_per_**** members, like function @p interact:
457  * template<unsigned int subdimA, unsigned int subdimB> static unsigned int count();
458  */
459  template<unsigned int subdim> static unsigned int count();
460 
461  /**
462  * @param sid - index of a sub-simplex in a simplex
463  * return an array of barycentric coordinates on <dim> simplex from <subdim> simplex
464  * for example: simplex<3> - ABCD and its subsubsimplex<1> AD (line index: 3)
465  * 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)
466  * this method creates array ((1,0,0,0),(0,0,0,1))
467  */
468  template<unsigned int subdim> static arma::mat::fixed<dim+1,subdim+1> bary_coords(unsigned int sid);
469 
470  /** Interpolate barycentric coords to a higher dimension of a simplex.
471  * @param coord - barycentric coords of a point on a sub-simplex
472  * @param sub_simplex_idx - id of sub-simplex on a simplex
473  */
474  template<unsigned int subdim> static BaryPoint interpolate(arma::vec::fixed<subdim+1> coord, int sub_simplex_idx);
475 
476 
477  /**
478  * Basic line interpolation.
479  */
480  static BaryPoint line_barycentric_interpolation(BaryPoint first_coords,
481  BaryPoint second_coords,
482  double first_theta, double second_theta, double theta);
483 
484  /**
485  * Usage:
486  * RefElement<3>::interact(Interaction<2,0>(1))
487  * (means: In tetrahedron <3>, give indices of sides <2>, connected by node <0> with index 1)
488  * RefElement<3>::interact(Interaction<2,0>(1))[1]
489  * (as above, but give only the side with index 1)
490  *
491  * Template usage: RefElement<dim>::interact(Interaction<OutDim, InDim>(i))[j]
492  * (means: on dim-dimensional reference element, go on InDim-dimensional subelement with index i,
493  * which connects OutDim-dimnesional subelements and select the one with index j)
494  *
495  * This method serves as an interface to topology information of the reference element.
496  * It returns indices of OutDim-dimensional object
497  * of InDim-dimnesional object of given index
498  * in dim-dimnesional reference element.
499  * @tparam interaction - auxilliary object carying the index and the template arguments OutDim and InDim
500  * @tparam OutDim - output dimension (give me node-0, line-1, side-2), <= dim
501  * @tparam InDim - input dimension (for node-0, line-1, side-2), <= dim
502  * @return vector of indices of OutDim-dimensional subelements represented by @p IdxVector object.
503  *
504  * possible calls:
505  * dim OutDim InDim return
506  * 1,2,3 0 1 InDim+1 - give me indices of nodes of line of given index
507  * 3 0 2 InDim+1 - give me indices of nodes of a side (triangle) of given index
508  * 3 1 2 InDim+1 - give me indices of lines of side (triangle) of given index
509  *
510  * 1,2,3 1 0 dim-InDim - give me indices of lines with common node of given index
511  * 3 2 0 dim-InDim - give me indices of sides (triangles) with common node of given index
512  * 3 2 1 dim-InDim - give me indices of sides (triangles) with common line of given index
513  *
514  */
515  template < template <unsigned int OutDim, unsigned int InDim> class TInteraction, unsigned int OutDim, unsigned int InDim>
516  static const IdxVector< (InDim>OutDim ? InDim+1 : dim-InDim) > interact( TInteraction<OutDim,InDim> interaction );
517 
518 
519 private:
520  /// Internal part of the interact function.
521  template<unsigned int OutDim, unsigned int InDim>
522  static const IdxVector< (InDim>OutDim ? InDim+1 : dim-InDim) > interact_(unsigned int index);
523 
524  static const std::vector<IdxVector<n_nodes_per_line>> line_nodes_; ///< [n_lines] For given line, returns its nodes indices.
525  static const std::vector<IdxVector<n_lines_per_node>> node_lines_; ///< [n_nodes] For given node, returns lines indices.
526  static const std::vector<IdxVector<n_nodes_per_side>> side_nodes_; ///< [n_sides] For given side, returns nodes indices. For @p dim == 3.
527  static const std::vector<IdxVector<n_sides_per_node>> node_sides_; ///< [n_nodes] For given node, returns sides indices. For @p dim == 3.
528  static const std::vector<IdxVector<n_sides_per_line>> line_sides_; ///< [n_lines] For given line, returns sides indices. For @p dim == 3.
529  static const std::vector<IdxVector<n_lines_per_side>> side_lines_; ///< [n_sides] For given side, returns lines indices. For @p dim == 3.
530 
531  //TODO: implement for 1d and 2d
532  /**
533  * Consider an n-face (node, edge, face, bulk) with dimension `subdim` and
534  * index within subdimension `idx`. Barycentric coordinates of all points
535  * on the n-face have unique pattern of zero coordinates.
536  *
537  * topology_zeros_[subdim][idx] is a bitfield with '1' where the pattern have zeros.
538  */
539  static const IdxVector<(n_lines > n_nodes) ? n_lines : n_nodes> topology_zeros_[dim+1];
540 };
541 
542 
543 // Declarations of explicit specialization of static memebers.
544 
545 template<> const IdxVector<1> RefElement<0>::topology_zeros_[];
546 template<> const IdxVector<2> RefElement<1>::topology_zeros_[];
547 template<> const IdxVector<3> RefElement<2>::topology_zeros_[];
548 template<> const IdxVector<6> RefElement<3>::topology_zeros_[];
549 
550 
551 template<> const std::vector<IdxVector<2>> RefElement<1>::line_nodes_;
552 template<> const std::vector<IdxVector<2>> RefElement<2>::line_nodes_;
553 template<> const std::vector<IdxVector<2>> RefElement<3>::line_nodes_;
554 template<> const std::vector<IdxVector<1>> RefElement<1>::node_lines_;
555 template<> const std::vector<IdxVector<2>> RefElement<2>::node_lines_;
556 
557 template<> const std::vector<IdxVector<3>> RefElement<3>::node_lines_;
558 template<> const std::vector<IdxVector<3>> RefElement<3>::side_nodes_;
559 template<> const std::vector<IdxVector<3>> RefElement<3>::node_sides_;
560 template<> const std::vector<IdxVector<2>> RefElement<3>::line_sides_;
561 template<> const std::vector<IdxVector<3>> RefElement<3>::side_lines_;
562 
563 
564 template<> const IdxVector<1> RefElement<0>::topology_zeros_[];
565 template<> const IdxVector<2> RefElement<1>::topology_zeros_[];
566 template<> const IdxVector<3> RefElement<2>::topology_zeros_[];
567 template<> const IdxVector<6> RefElement<3>::topology_zeros_[];
568 
569 
570 
571 
572 
573 
574 // 0: nodes of nodes
575 // 1: nodes of lines
576 // 2: nodes of sides
577 // 3: nodes of tetrahedron
578 template<> const std::vector< std::vector< std::vector<unsigned int> > > RefElement<0>::nodes_of_subelements;
579 template<> const std::vector< std::vector< std::vector<unsigned int> > > RefElement<1>::nodes_of_subelements;
580 template<> const std::vector< std::vector< std::vector<unsigned int> > > RefElement<2>::nodes_of_subelements;
581 template<> const std::vector< std::vector< std::vector<unsigned int> > > RefElement<3>::nodes_of_subelements;
582 
583 
584 
585 
586 
587 
588 
589 
590 
591 /************************* template implementation ****************************/
592 
593 template<unsigned int dim>
594 template<unsigned int subdim> inline
595 arma::mat::fixed<dim+1,subdim+1> RefElement<dim>::bary_coords(unsigned int sid){
596  ASSERT_LT_DBG(subdim, dim).error("Dimension mismatch!");
597  arma::mat::fixed<dim+1,subdim+1> bary_c;
598 
599  for(unsigned int i = 0; i < subdim+1; i++){
600  unsigned int nid = interact_<0,subdim>(sid)[i];
601  bary_c.col(i).zeros();
602  bary_c.col(i)(nid) = 1;
603  }
604 
605  return bary_c;
606 }
607 
608 
609 template<unsigned int dim> inline
610 arma::vec::fixed<dim> RefElement<dim>::node_coords(unsigned int nid)
611 {
612  ASSERT_LT_DBG(nid, n_nodes).error("Node number is out of range!");
613 
614  arma::vec::fixed<dim> p;
615  p.zeros();
616 
617  if (nid > 0)
618  p(nid-1) = 1;
619 
620  return p;
621 }
622 
623 
624 template<unsigned int dim>
625 template<unsigned int subdim>
626 auto RefElement<dim>::interpolate(arma::vec::fixed<subdim+1> coord, int sub_simplex_idx) -> BaryPoint
627 {
628  return RefElement<dim>::bary_coords<subdim>(sub_simplex_idx)*coord;
629 }
630 /*
631 template <unsigned int Size>
632 IdxVector<Size>::IdxVector(std::array<unsigned int,Size> data_in)
633 : data_(data_in){}
634 
635 template <unsigned int Size>
636 IdxVector<Size>::IdxVector(std::initializer_list<unsigned int> data_in)
637 {
638  ASSERT_EQ_DBG(data_in.size(), Size).error("Incorrect data size.");
639  std::copy(data_in.begin(), data_in.end(), data_);
640 }
641 
642 template <unsigned int Size>
643 inline unsigned int IdxVector<Size>::operator[](unsigned int idx) const
644 { ASSERT_LT_DBG(idx, Size).error("Index out of bounds.");
645  return data_[idx]; }
646 
647 */
648 
649 template<> template<> inline unsigned int RefElement<3>::count<0>()
650 { return n_nodes; }
651 template<> template<> inline unsigned int RefElement<3>::count<1>()
652 { return n_lines; }
653 template<> template<> inline unsigned int RefElement<3>::count<2>()
654 { return n_sides; }
655 template<> template<> inline unsigned int RefElement<3>::count<3>()
656 { return 1; }
657 template<> template<> inline unsigned int RefElement<2>::count<0>()
658 { return n_nodes; }
659 template<> template<> inline unsigned int RefElement<2>::count<1>()
660 { return n_lines; }
661 template<> template<> inline unsigned int RefElement<2>::count<2>()
662 { return 1; }
663 template<> template<> inline unsigned int RefElement<2>::count<3>()
664 { return 0; }
665 template<> template<> inline unsigned int RefElement<1>::count<0>()
666 { return n_nodes; }
667 template<> template<> inline unsigned int RefElement<1>::count<1>()
668 { return 1; }
669 template<> template<> inline unsigned int RefElement<1>::count<2>()
670 { return 0; }
671 template<> template<> inline unsigned int RefElement<1>::count<3>()
672 { return 0; }
673 template<> template<> inline unsigned int RefElement<0>::count<0>()
674 { return 1; }
675 template<> template<> inline unsigned int RefElement<0>::count<1>()
676 { return 0; }
677 template<> template<> inline unsigned int RefElement<0>::count<2>()
678 { return 0; }
679 template<> template<> inline unsigned int RefElement<0>::count<3>()
680 { return 0; }
681 
682 template<unsigned int dim>
683 template<unsigned int subdim>
684 unsigned int RefElement<dim>::topology_idx(unsigned int zeros_positions)
685 {
686  for(unsigned int i=0; i < RefElement<dim>::count<subdim>(); i++){
687  if(zeros_positions == topology_zeros_[subdim][i]) return i;
688  }
689  ASSERT(0).error("Undefined zero pattern.");
690  return -1;
691 }
692 
693 
694 /// This function is for "side_nodes" - for given side, give me nodes (0->0, 1->1).
695 template<> template<> inline const IdxVector<1> RefElement<1>::interact_<0,0>(unsigned int i)
696 { ASSERT_LT_DBG(i, RefElement<1>::n_nodes).error("Index out of bounds.");
697  return IdxVector<1>({i});}
698 
699 /// For line i {0}, give me indices of its nodes.
700 template<> template<> inline const IdxVector<2> RefElement<1>::interact_<0,1>(unsigned int i)
701 { ASSERT_LT_DBG(i, RefElement<1>::n_lines).error("Index out of bounds.");
702  return line_nodes_[i];}
703 
704 /// For line i {0,1,2}, give me indices of its nodes.
705 template<> template<> inline const IdxVector<2> RefElement<2>::interact_<0,1>(unsigned int i)
706 { ASSERT_LT_DBG(i, RefElement<2>::n_lines).error("Index out of bounds.");
707  return line_nodes_[i];}
708 
709 /// For line i {0,1,2,3,4,5}, give me indices of its nodes.
710 template<> template<> inline const IdxVector<2> RefElement<3>::interact_<0,1>(unsigned int i)
711 { ASSERT_LT_DBG(i, RefElement<3>::n_lines).error("Index out of bounds.");
712  return line_nodes_[i];}
713 
714 /// For node i {0,1}, give me indices of lines.
715 template<> template<> inline const IdxVector<1> RefElement<1>::interact_<1,0>(unsigned int i)
716 { ASSERT_LT_DBG(i, RefElement<1>::n_nodes).error("Index out of bounds.");
717  return node_lines_[i];}
718 
719 /// For node i {0,1,2}, give me indices of lines.
720 template<> template<> inline const IdxVector<2> RefElement<2>::interact_<1,0>(unsigned int i)
721 { ASSERT_LT_DBG(i, RefElement<2>::n_nodes).error("Index out of bounds.");
722  return node_lines_[i];}
723 
724 /// For node i {0,1,2,3}, give me indices of lines.
725 template<> template<> inline const IdxVector<3> RefElement<3>::interact_<1,0>(unsigned int i)
726 { ASSERT_LT_DBG(i, RefElement<3>::n_nodes).error("Index out of bounds.");
727  return node_lines_[i];}
728 
729 /// For side i {0,1,2}, give me indices of its nodes.
730 template<> template<> inline const IdxVector<3> RefElement<3>::interact_<0,2>(unsigned int i)
731 { ASSERT_LT_DBG(i, RefElement<3>::n_sides).error("Index out of bounds.");
732  return side_nodes_[i];}
733 
734 /// For node i {0,1,2,3}, give me indices of sides.
735 template<> template<> inline const IdxVector<3> RefElement<3>::interact_<2,0>(unsigned int i)
736 { ASSERT_LT_DBG(i, RefElement<3>::n_sides).error("Index out of bounds.");
737  return node_sides_[i];}
738 
739 /// For line i {0,1,2,3}, give me indices of sides.
740 template<> template<> inline const IdxVector<2> RefElement<3>::interact_<2,1>(unsigned int i)
741 { ASSERT_LT_DBG(i, RefElement<3>::n_lines).error("Index out of bounds.");
742  return line_sides_[i];}
743 
744 /// For side i {0,1,2}, give me indices of its lines.
745 template<> template<> inline const IdxVector<3> RefElement<3>::interact_<1,2>(unsigned int i)
746 { ASSERT_LT_DBG(i, RefElement<3>::n_sides).error("Index out of bounds.");
747  return side_lines_[i];}
748 
749 template<unsigned int dim> template<unsigned int OutDim, unsigned int InDim>
750 inline const IdxVector< (InDim>OutDim ? InDim+1 : dim-InDim) > RefElement<dim>::interact_(unsigned int i)
751 {
752  ASSERT(false)(dim)(OutDim)(InDim)(i).error("Not implemented.");
753  //ASSERT_LT_DBG(OutDim, dim);
754  //ASSERT_LT_DBG(InDim, dim);
755  return IdxVector< (InDim>OutDim ? InDim+1 : dim-InDim) >(); // just to avoid warning for missing return
756 }
757 
758 
759 template<unsigned int dim>
760 template < template <unsigned int OutDim, unsigned int InDim> class TInteraction, unsigned int OutDim, unsigned int InDim>
761 inline const IdxVector< (InDim>OutDim ? InDim+1 : dim-InDim) > RefElement<dim>::interact( TInteraction<OutDim,InDim> interaction )
762 {
763  return interact_<OutDim,InDim>(interaction.i_);
764 }
765 
766 #endif /* REF_ELEMENT_HH_ */
RefElement::normal_vector
static LocalPoint normal_vector(unsigned int sid)
RefElement::FaceBaryPoint
Armor::ArmaVec< double, dim > FaceBaryPoint
Definition: ref_element.hh:320
armor.hh
RefElement::line_barycentric_interpolation
static BaryPoint line_barycentric_interpolation(BaryPoint first_coords, BaryPoint second_coords, double first_theta, double second_theta, double theta)
Definition: ref_element.cc:524
RefElement
Definition: ref_element.hh:309
RefElement::normal_orientation
static unsigned int normal_orientation(unsigned int sid)
Definition: ref_element.cc:301
RefElement::BaryPoint
Armor::ArmaVec< double, dim+1 > BaryPoint
Definition: ref_element.hh:319
RefElement::LocalPoint
arma::vec::fixed< dim > LocalPoint
Definition: ref_element.hh:312
ASSERT
#define ASSERT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:347
RefElement::n_sides_per_line
static const unsigned int n_sides_per_line
Number of sides with one common line. dim == 3.
Definition: ref_element.hh:373
RefElement::bary_coords
static arma::mat::fixed< dim+1, subdim+1 > bary_coords(unsigned int sid)
Definition: ref_element.hh:595
RefElement::side_nodes_
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:526
Interaction::i_
unsigned int i_
Definition: ref_element.hh:305
RefElement::n_nodes
static const unsigned int n_nodes
Number of nodes.
Definition: ref_element.hh:369
asserts.hh
Definitions of ASSERTS.
RefElement::interpolate
static BaryPoint interpolate(arma::vec::fixed< subdim+1 > coord, int sub_simplex_idx)
std::vector
Definition: doxy_dummy_defs.hh:7
RefElement::zeros_positions
static std::pair< unsigned int, unsigned int > zeros_positions(const BaryPoint &barycentric, double tolerance=std::numeric_limits< double >::epsilon() *2)
Definition: ref_element.cc:379
RefElement::count
static unsigned int count()
RefElement::side_lines_
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:529
RefElement::n_sides
static const unsigned int n_sides
Number of sides.
Definition: ref_element.hh:368
RefElement::make_bary_unit_vec
static BarycentricUnitVec make_bary_unit_vec()
Definition: ref_element.cc:403
RefElement::centers_of_subelements
static CentersList centers_of_subelements(unsigned int sub_dim)
Definition: ref_element.cc:451
RefElement::nodes_of_subelements
static const std::vector< std::vector< std::vector< unsigned int > > > nodes_of_subelements
Definition: ref_element.hh:404
RefElement::interact_
static const IdxVector<(InDim >OutDim ? InDim+1 :dim-InDim) > interact_(unsigned int index)
Internal part of the interact function.
RefElement::node_sides_
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:527
RefElement::bary_to_local
static LocalPoint bary_to_local(const BaryPoint &bp)
Converts from barycentric to local coordinates.
Definition: ref_element.cc:286
RefElement::n_nodes_per_line
static const unsigned int n_nodes_per_line
Number of nodes in one line.
Definition: ref_element.hh:372
RefElement::n_sides_per_node
static const unsigned int n_sides_per_node
Number of sides with one common line.
Definition: ref_element.hh:374
RefElement::clip
static BaryPoint clip(const BaryPoint &barycentric)
Definition: ref_element.cc:417
IdxVector
std::array< unsigned int, Size > IdxVector
Definition: ref_element.hh:293
RefElement::topology_zeros_
static const IdxVector<(n_lines > n_nodes) ? n_lines :n_nodes > topology_zeros_[dim+1]
Definition: ref_element.hh:539
RefElement::oposite_node
static unsigned int oposite_node(unsigned int sid)
Definition: ref_element.cc:294
RefElement::interact
static const IdxVector<(InDim >OutDim ? InDim+1 :dim-InDim) > interact(TInteraction< OutDim, InDim > interaction)
RefElement::CentersList
const typedef std::vector< LocalPoint > & CentersList
Definition: ref_element.hh:434
tolerance
Definition: isotherm.hh:67
RefElement::DECLARE_EXCEPTION
DECLARE_EXCEPTION(ExcInvalidPermutation,<< "Side permutation not found.\n")
Interaction
Definition: ref_element.hh:303
RefElement::barycentric_on_face
static FaceBaryPoint barycentric_on_face(const BaryPoint &barycentric, unsigned int i_face)
Definition: ref_element.cc:364
RefElement::line_nodes_
static const std::vector< IdxVector< n_nodes_per_line > > line_nodes_
[n_lines] For given line, returns its nodes indices.
Definition: ref_element.hh:524
RefElement::n_lines
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:380
RefElement::node_coords
static LocalPoint node_coords(unsigned int nid)
Definition: ref_element.hh:610
RefElement::n_lines_per_side
static const unsigned int n_lines_per_side
Number of lines on boundary of one side.
Definition: ref_element.hh:377
RefElement::line_between_faces
static unsigned int line_between_faces(unsigned int f1, unsigned int f2)
Interaction::Interaction
Interaction(unsigned int i)
Definition: ref_element.hh:304
ASSERT_LT_DBG
#define ASSERT_LT_DBG(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:300
RefElement::local_to_bary
static BaryPoint local_to_bary(const LocalPoint &lp)
Converts from local to barycentric coordinates.
Definition: ref_element.cc:272
RefElement::BarycentricUnitVec
std::vector< BaryPoint > BarycentricUnitVec
Definition: ref_element.hh:419
RefElement::line_sides_
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:528
RefElement::side_measure
static double side_measure(unsigned int sid)
RefElement::n_lines_per_node
static const unsigned int n_lines_per_node
Number of lines with one common node.
Definition: ref_element.hh:371
RefElement::n_nodes_per_side
static const unsigned int n_nodes_per_side
Number of nodes on one side.
Definition: ref_element.hh:370
Armor::ArmaVec
typename arma::Col< Type >::template fixed< nr > ArmaVec
Definition: armor.hh:505
RefElement::topology_idx
static unsigned int topology_idx(unsigned int zeros_positions)
Definition: ref_element.hh:684
RefElement::node_lines_
static const std::vector< IdxVector< n_lines_per_node > > node_lines_
[n_nodes] For given node, returns lines indices.
Definition: ref_element.hh:525
Generated by   doxygen 1.8.17