Flow123d  JS_before_hm-1712-g4d968e368
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  * - functions from DEAL.ii:
53  * bool is_inside_unit_cell( point )
54  * line_to_cell_vertices(line, vertex) vertex index on line to index on whole element
55  * face_to_cell_vertices(face, vertex, Orientation), Orientation should be some class describing permutation of shared face to element's face/side
56  * face_to_cel_lines
57  * standard_to_real_face_vertex(vertex, Orientation), maps vertex to permuted face
58  * real_to_standard_face_vertex - inverse
59  * ... same for line; we should need also something like standard_to_real_line_vertex; Deal dosn;t support Orientation changes for 2D element faces
60  * Point unit_cell_vertex(vertex) - coordinates
61  * project_to_unit_cell
62  * distance_to_unit_cell
63  * d_linear_shape_function
64  *
65  * - can not change numbering of element sides due to DarcyFlow, which use hardwired side numbering in construction of basis functions
66  * - any change of side numbering requires also change in flow/old_bcd.cc
67  *
68  *
69  */
70 
71 #ifndef REF_ELEMENT_HH_
72 #define REF_ELEMENT_HH_
73 
74 #include <vector> // for vector
75 #include <array>
76 #include <armadillo>
77 #include "system/armor.hh"
78 #include "system/asserts.hh"
79 
80 
81 /*
82  * Ordering of nodes and sides in reference elements
83  * =================================================
84  *
85  * TODO we want the following (22.10.):
86  *
87  * 1D element (line segment) 2D element (triangle) 3D element (tetrahedron)
88  *
89  * z
90  * .
91  * ,/
92  * /
93  * 3
94  * y ,/|`\
95  * ^ ,/ | `\
96  * | ,/ '. `\
97  * 2 ,/ | `\
98  * |`\ ,/ | `\
99  * | `\ 0-----------'.--------1 --> x
100  * | `\ `\. | ,/
101  * | `\ `\. | ,/
102  * | `\ `\. '. ,/
103  * 0----------1 --> x 0----------1 --> x `\. |/
104  * `2
105  * `\.
106  * `y
107  *
108  * side id node ids side id node ids side id node ids normal
109  * 0 0 0 0,1 0 0,1,2 OUT
110  * 1 1 1 0,2 1 0,1,3 IN
111  * 2 1,2 2 0,2,3 OUT
112  * 3 1,2,3 IN
113  *
114  *
115  * nodes coordinates:
116  * 0 [0] 0 [0,0] 0 [0,0,0]
117  * 1 [1] 1 [1,0] 1 [1,0,0]
118  * 2 [0,1] 2 [0,1,0]
119  * 3 [0,0,1]
120  *
121  * barycentric coordinates of nodes:
122  * 0 [1,0] 0 [1,0,0] 0 [1,0,0,0]
123  * 1 [0,1] 1 [0,1,0] 1 [0,1,0,0]
124  * 2 [0,0,1] 2 [0,0,1,0]
125  * 3 [0,0,0,1]
126  */
127 
128 /** Auxilliary class representing vector of indices (unsigned int).
129  * @tparam Size is the fixed size of the vector.
130  */
131 /*
132 template<unsigned int Size>
133 class IdxVector{
134  unsigned int data_[Size]; ///< Array with indices.
135 
136  public:
137  /// Constructor taking in array of indices.
138  IdxVector(std::array<unsigned int,Size> data_in);
139  /// Constructor enabling creating object with initializer list {...}.
140  IdxVector(std::initializer_list<unsigned int> data_in);
141  /// Getter for index @p idx.
142  unsigned int operator[](unsigned int idx) const;
143 };
144 */
145 
146 template<std::size_t Size>
147 using IdxVector = std::array<unsigned int, Size>;
148 
149 
150 /** Auxilliary structure that is used to pass template arguments into interact function of RefElement:
151  * RefElement<dim>::interact( Interaction<OutDim,InDim>(i) )
152  *
153  * This enables automatic deduction of dimensional template arguments.
154  * @see @p RefElement<dim>::interact
155  */
156 template <unsigned int OutDim, unsigned int InDim>
157 struct Interaction {
158  Interaction(unsigned int i) : i_(i) {}
159  unsigned int i_;
160 };
161 
162 template<unsigned int dim>
163 class RefElement
164 {
165 public:
166  typedef arma::vec::fixed<dim> LocalPoint;
167  /**
168  * Barycentric coordinates.
169  *
170  * e.g. coordinates (a,b,c) on triangle with vertices X, Y, Z
171  * represents a point: a*X+b*Y+c*Z
172  */
173  typedef Armor::ArmaVec<double, dim+1> BaryPoint;
174  typedef Armor::ArmaVec<double, dim> FaceBaryPoint;
175 
176  DECLARE_EXCEPTION( ExcInvalidPermutation, << "Side permutation not found.\n" );
177 
178  /**
179  * Return coordinates of given node.
180  * @see the class documentation @p RefElement
181  * @param nid Node number.
182  * NOTE: Implementation is dependent on current node and side numbering.
183  */
184  static LocalPoint node_coords(unsigned int nid);
185 
186  /**
187  * Compute normal vector to a given side.
188  * @param sid Side number.
189  */
190  static LocalPoint normal_vector(unsigned int sid);
191 
192 
193  /**
194  * If the given barycentric coordinate is in the ref. element, return unchanged.
195  * If the given barycentric coordinate is out of the ref. element,
196  * project it on the surface of the ref. element.
197  */
198  static BaryPoint clip(const BaryPoint &barycentric);
199 
200  /** Returns orientation of the normal of side @p sid. 0 -> OUT, 1 -> IN.
201  * NOTE: Implementation is dependent on current node and side numbering.
202  */
203  static unsigned int normal_orientation(unsigned int sid);
204 
205  static double side_measure(unsigned int sid);
206 
207  /**
208  * Returns index of the node that is oposite to side of given index @p sid.
209  * Note: It is dependent on current node and side numbering.
210  * @param sid Side number.
211  * NOTE: Implementation is dependent on current node and side numbering.
212  */
213  static unsigned int oposite_node(unsigned int sid);
214 
215  /**
216  * Return index of 1D line, shared by two faces @p f1 and @p f2 of the reference tetrahedron.
217  * Implemented only for @p dim == 3.
218  */
219  static unsigned int line_between_faces(unsigned int f1, unsigned int f2);
220 
221 
222  static const unsigned int n_sides = dim + 1; ///< Number of sides.
223  static const unsigned int n_nodes = dim + 1; ///< Number of nodes.
224  static const unsigned int n_nodes_per_side = dim; ///< Number of nodes on one side.
225  static const unsigned int n_lines_per_node = dim; ///< Number of lines with one common node.
226  static const unsigned int n_nodes_per_line = 2; ///< Number of nodes in one line.
227  static const unsigned int n_sides_per_line = 2; ///< Number of sides with one common line. @p dim == 3.
228  static const unsigned int n_sides_per_node = dim; ///< Number of sides with one common line.
229 
230  /// Number of lines on boundary of one side.
231  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
232 
233  /// Number of lines, i.e. @p object of dimension @p dim-2 on the boundary of the reference element.
234  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
235 
236 
237 // /**
238 // * Node numbers for each side.
239 // */
240 // static const unsigned int side_nodes[n_sides][n_nodes_per_side];
241 //
242 // /**
243 // * Indices of 1D lines of the 2D sides of an tetrahedron. Nonempty only for @p dim==3.
244 // */
245 // static const unsigned int side_lines[n_sides][n_lines_per_side];
246 //
247 // /**
248 // * Nodes of 1D lines of the tetrahedron.
249 // */
250 // static const unsigned int line_nodes[n_lines][2];
251 //
252 // /**
253 // * Indices of sides for each line. Nonempty only for @p dim==3 and @p dim==2.
254 // */
255 // static const unsigned int line_sides[n_lines][2];
256 
257 
258  static const std::vector< std::vector< std::vector<unsigned int> > > nodes_of_subelements;
259 
260  /**
261  * Number of permutations of nodes on sides.
262  * dim value
263  * -----------
264  * 1 1
265  * 2 2
266  * 3 6
267  */
268  static constexpr unsigned int n_side_permutations = (dim+1)*(2*dim*dim-5*dim+6)/6;
269 
270  /**
271  * Permutations of nodes on sides.
272  * [n_side_permutations][n_nodes_per_side]
273  */
274  static const std::vector< std::vector<unsigned int> > side_permutations;
275 
276  /**
277  * For a given permutation @p p of nodes finds its index within @p side_permutations.
278  * @param p Permutation of nodes.
279  */
280  static unsigned int permutation_index(unsigned int p[n_nodes_per_side]);
281 
282  /** @brief Converts from local to barycentric coordinates.
283  * @param lp point in local coordinates (x,y)
284  * @return point in barycentric coordinates (1-x-y, x, y)
285  */
286  static BaryPoint local_to_bary(const LocalPoint& lp);
287 
288  /** @brief Converts from barycentric to local coordinates.
289  * @param bp point in barycentric coordinates
290  * @return point in local coordinates
291  */
292  static LocalPoint bary_to_local(const BaryPoint& bp);
293 
294  typedef std::vector<BaryPoint> BarycentricUnitVec;
295 
296  /**
297  * Used in the clip method.
298  */
299  static BarycentricUnitVec make_bary_unit_vec();
300 
301  /**
302  * For given barycentric coordinates on the ref element returns barycentric coordinates
303  * on the ref. element of given face. Assumes that the input point is on the face.
304  * Barycentric order: (complanatory, local_coords )
305  */
306  static FaceBaryPoint barycentric_on_face(const BaryPoint &barycentric, unsigned int i_face);
307 
308 
309  typedef const std::vector<LocalPoint> & CentersList;
310  static CentersList centers_of_subelements(unsigned int sub_dim);
311 
312  /**
313  * Return (1) number of zeros and (2) positions of zeros in barycentric coordinates.
314  * @p tolerance serves for testing zero values of @p barycentric coordinates.
315  */
316  static std::pair<unsigned int, unsigned int> zeros_positions(const BaryPoint &barycentric,
317  double tolerance = std::numeric_limits<double>::epsilon()*2);
318 
319  /**
320  * According to positions of zeros in barycentric coordinates, it gives the index of subdim-simplex
321  * in the reference element. Number of zeros must be equal to (3-subdim).
322  * e.g.:
323  * if 1 zeros, return index of side (subdim 2)
324  * if 2 zeros, return index of edge (subdim 1)
325  * if 3 zeros, return index of vertex (subdim 0)
326  */
327  template<unsigned int subdim> static unsigned int topology_idx(unsigned int zeros_positions);
328 
329  /** Function returns number of subdim-simplices inside dim-simplex.
330  * The aim is covering all the n_**** members with a single function.
331  * TODO: think of generalization for n_****_per_**** members, like function @p interact:
332  * template<unsigned int subdimA, unsigned int subdimB> static unsigned int count();
333  */
334  template<unsigned int subdim> static unsigned int count();
335 
336  /**
337  * @param sid - index of a sub-simplex in a simplex
338  * return an array of barycentric coordinates on <dim> simplex from <subdim> simplex
339  * for example: simplex<3> - ABCD and its subsubsimplex<1> AD (line index: 3)
340  * 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)
341  * this method creates array ((1,0,0,0),(0,0,0,1))
342  */
343  template<unsigned int subdim> static arma::mat::fixed<dim+1,subdim+1> bary_coords(unsigned int sid);
344 
345  /** Interpolate barycentric coords to a higher dimension of a simplex.
346  * @param coord - barycentric coords of a point on a sub-simplex
347  * @param sub_simplex_idx - id of sub-simplex on a simplex
348  */
349  template<unsigned int subdim> static BaryPoint interpolate(arma::vec::fixed<subdim+1> coord, int sub_simplex_idx);
350 
351 
352  /**
353  * Basic line interpolation.
354  */
355  static BaryPoint line_barycentric_interpolation(BaryPoint first_coords,
356  BaryPoint second_coords,
357  double first_theta, double second_theta, double theta);
358 
359  /**
360  * Usage:
361  * RefElement<3>::interact(Interaction<2,0>(1))
362  * (means: In tetrahedron <3>, give indices of sides <2>, connected by node <0> with index 1)
363  * RefElement<3>::interact(Interaction<2,0>(1))[1]
364  * (as above, but give only the side with index 1)
365  *
366  * Template usage: RefElement<dim>::interact(Interaction<OutDim, InDim>(i))[j]
367  * (means: on dim-dimensional reference element, go on InDim-dimensional subelement with index i,
368  * which connects OutDim-dimnesional subelements and select the one with index j)
369  *
370  * This method serves as an interface to topology information of the reference element.
371  * It returns indices of OutDim-dimensional object
372  * of InDim-dimnesional object of given index
373  * in dim-dimnesional reference element.
374  * @tparam interaction - auxilliary object carying the index and the template arguments OutDim and InDim
375  * @tparam OutDim - output dimension (give me node-0, line-1, side-2), <= dim
376  * @tparam InDim - input dimension (for node-0, line-1, side-2), <= dim
377  * @return vector of indices of OutDim-dimensional subelements represented by @p IdxVector object.
378  *
379  * possible calls:
380  * dim OutDim InDim return
381  * 1,2,3 0 1 InDim+1 - give me indices of nodes of line of given index
382  * 3 0 2 InDim+1 - give me indices of nodes of a side (triangle) of given index
383  * 3 1 2 InDim+1 - give me indices of lines of side (triangle) of given index
384  *
385  * 1,2,3 1 0 dim-InDim - give me indices of lines with common node of given index
386  * 3 2 0 dim-InDim - give me indices of sides (triangles) with common node of given index
387  * 3 2 1 dim-InDim - give me indices of sides (triangles) with common line of given index
388  *
389  */
390  template < template <unsigned int OutDim, unsigned int InDim> class TInteraction, unsigned int OutDim, unsigned int InDim>
391  static const IdxVector< (InDim>OutDim ? InDim+1 : dim-InDim) > interact( TInteraction<OutDim,InDim> interaction );
392 
393 
394 private:
395  /// Internal part of the interact function.
396  template<unsigned int OutDim, unsigned int InDim>
397  static const IdxVector< (InDim>OutDim ? InDim+1 : dim-InDim) > interact_(unsigned int index);
398 
399  static const std::vector<IdxVector<n_nodes_per_line>> line_nodes_; ///< [n_lines] For given line, returns its nodes indices.
400  static const std::vector<IdxVector<n_lines_per_node>> node_lines_; ///< [n_nodes] For given node, returns lines indices.
401  static const std::vector<IdxVector<n_nodes_per_side>> side_nodes_; ///< [n_sides] For given side, returns nodes indices. For @p dim == 3.
402  static const std::vector<IdxVector<n_sides_per_node>> node_sides_; ///< [n_nodes] For given node, returns sides indices. For @p dim == 3.
403  static const std::vector<IdxVector<n_sides_per_line>> line_sides_; ///< [n_lines] For given line, returns sides indices. For @p dim == 3.
404  static const std::vector<IdxVector<n_lines_per_side>> side_lines_; ///< [n_sides] For given side, returns lines indices. For @p dim == 3.
405 
406  //TODO: implement for 1d and 2d
407  /**
408  * Consider an n-face (node, edge, face, bulk) with dimension `subdim` and
409  * index within subdimension `idx`. Barycentric coordinates of all points
410  * on the n-face have unique pattern of zero coordinates.
411  *
412  * topology_zeros_[subdim][idx] is a bitfield with '1' where the pattern have zeros.
413  */
414  static const IdxVector<(n_lines > n_nodes) ? n_lines : n_nodes> topology_zeros_[dim+1];
415 };
416 
417 
418 // Declarations of explicit specialization of static memebers.
419 
420 template<> const IdxVector<1> RefElement<0>::topology_zeros_[];
421 template<> const IdxVector<2> RefElement<1>::topology_zeros_[];
422 template<> const IdxVector<3> RefElement<2>::topology_zeros_[];
423 template<> const IdxVector<6> RefElement<3>::topology_zeros_[];
424 
425 
426 template<> const std::vector<IdxVector<2>> RefElement<1>::line_nodes_;
427 template<> const std::vector<IdxVector<2>> RefElement<2>::line_nodes_;
428 template<> const std::vector<IdxVector<2>> RefElement<3>::line_nodes_;
429 template<> const std::vector<IdxVector<1>> RefElement<1>::node_lines_;
430 template<> const std::vector<IdxVector<2>> RefElement<2>::node_lines_;
431 
432 template<> const std::vector<IdxVector<3>> RefElement<3>::node_lines_;
433 template<> const std::vector<IdxVector<3>> RefElement<3>::side_nodes_;
434 template<> const std::vector<IdxVector<3>> RefElement<3>::node_sides_;
435 template<> const std::vector<IdxVector<2>> RefElement<3>::line_sides_;
436 template<> const std::vector<IdxVector<3>> RefElement<3>::side_lines_;
437 
438 
439 template<> const std::vector< std::vector<unsigned int> > RefElement<0>::side_permutations;
440 template<> const std::vector< std::vector<unsigned int> > RefElement<1>::side_permutations;
441 template<> const std::vector< std::vector<unsigned int> > RefElement<2>::side_permutations;
442 template<> const std::vector< std::vector<unsigned int> > RefElement<3>::side_permutations;
443 
444 template<> const IdxVector<1> RefElement<0>::topology_zeros_[];
445 template<> const IdxVector<2> RefElement<1>::topology_zeros_[];
446 template<> const IdxVector<3> RefElement<2>::topology_zeros_[];
447 template<> const IdxVector<6> RefElement<3>::topology_zeros_[];
448 
449 
450 
451 
452 
453 
454 // 0: nodes of nodes
455 // 1: nodes of lines
456 // 2: nodes of sides
457 // 3: nodes of tetrahedron
458 template<> const std::vector< std::vector< std::vector<unsigned int> > > RefElement<0>::nodes_of_subelements;
459 template<> const std::vector< std::vector< std::vector<unsigned int> > > RefElement<1>::nodes_of_subelements;
460 template<> const std::vector< std::vector< std::vector<unsigned int> > > RefElement<2>::nodes_of_subelements;
461 template<> const std::vector< std::vector< std::vector<unsigned int> > > RefElement<3>::nodes_of_subelements;
462 
463 
464 
465 
466 
467 
468 
469 
470 
471 /************************* template implementation ****************************/
472 
473 template<unsigned int dim>
474 template<unsigned int subdim> inline
475 arma::mat::fixed<dim+1,subdim+1> RefElement<dim>::bary_coords(unsigned int sid){
476  ASSERT_LT_DBG(subdim, dim).error("Dimension mismatch!");
477  arma::mat::fixed<dim+1,subdim+1> bary_c;
478 
479  for(unsigned int i = 0; i < subdim+1; i++){
480  unsigned int nid = interact_<0,subdim>(sid)[i];
481  bary_c.col(i).zeros();
482  bary_c.col(i)(nid) = 1;
483  }
484 
485  return bary_c;
486 }
487 
488 
489 template<unsigned int dim> inline
490 arma::vec::fixed<dim> RefElement<dim>::node_coords(unsigned int nid)
491 {
492  ASSERT_LT_DBG(nid, n_nodes).error("Node number is out of range!");
493 
494  arma::vec::fixed<dim> p;
495  p.zeros();
496 
497  if (nid > 0)
498  p(nid-1) = 1;
499 
500  return p;
501 }
502 
503 
504 template<unsigned int dim>
505 template<unsigned int subdim>
506 auto RefElement<dim>::interpolate(arma::vec::fixed<subdim+1> coord, int sub_simplex_idx) -> BaryPoint
507 {
508  return RefElement<dim>::bary_coords<subdim>(sub_simplex_idx)*coord;
509 }
510 /*
511 template <unsigned int Size>
512 IdxVector<Size>::IdxVector(std::array<unsigned int,Size> data_in)
513 : data_(data_in){}
514 
515 template <unsigned int Size>
516 IdxVector<Size>::IdxVector(std::initializer_list<unsigned int> data_in)
517 {
518  ASSERT_EQ_DBG(data_in.size(), Size).error("Incorrect data size.");
519  std::copy(data_in.begin(), data_in.end(), data_);
520 }
521 
522 template <unsigned int Size>
523 inline unsigned int IdxVector<Size>::operator[](unsigned int idx) const
524 { ASSERT_LT_DBG(idx, Size).error("Index out of bounds.");
525  return data_[idx]; }
526 
527 */
528 
529 template<> template<> inline unsigned int RefElement<3>::count<0>()
530 { return n_nodes; }
531 template<> template<> inline unsigned int RefElement<3>::count<1>()
532 { return n_lines; }
533 template<> template<> inline unsigned int RefElement<3>::count<2>()
534 { return n_sides; }
535 template<> template<> inline unsigned int RefElement<3>::count<3>()
536 { return 1; }
537 template<> template<> inline unsigned int RefElement<2>::count<0>()
538 { return n_nodes; }
539 template<> template<> inline unsigned int RefElement<2>::count<1>()
540 { return n_lines; }
541 template<> template<> inline unsigned int RefElement<2>::count<2>()
542 { return 1; }
543 template<> template<> inline unsigned int RefElement<2>::count<3>()
544 { return 0; }
545 template<> template<> inline unsigned int RefElement<1>::count<0>()
546 { return n_nodes; }
547 template<> template<> inline unsigned int RefElement<1>::count<1>()
548 { return 1; }
549 template<> template<> inline unsigned int RefElement<1>::count<2>()
550 { return 0; }
551 template<> template<> inline unsigned int RefElement<1>::count<3>()
552 { return 0; }
553 template<> template<> inline unsigned int RefElement<0>::count<0>()
554 { return 1; }
555 template<> template<> inline unsigned int RefElement<0>::count<1>()
556 { return 0; }
557 template<> template<> inline unsigned int RefElement<0>::count<2>()
558 { return 0; }
559 template<> template<> inline unsigned int RefElement<0>::count<3>()
560 { return 0; }
561 
562 template<unsigned int dim>
563 template<unsigned int subdim>
564 unsigned int RefElement<dim>::topology_idx(unsigned int zeros_positions)
565 {
566  for(unsigned int i=0; i < RefElement<dim>::count<subdim>(); i++){
567  if(zeros_positions == topology_zeros_[subdim][i]) return i;
568  }
569  ASSERT(0).error("Undefined zero pattern.");
570  return -1;
571 }
572 
573 
574 /// This function is for "side_nodes" - for given side, give me nodes (0->0, 1->1).
575 template<> template<> inline const IdxVector<1> RefElement<1>::interact_<0,0>(unsigned int i)
576 { ASSERT_LT_DBG(i, RefElement<1>::n_nodes).error("Index out of bounds.");
577  return IdxVector<1>({i});}
578 
579 /// For line i {0}, give me indices of its nodes.
580 template<> template<> inline const IdxVector<2> RefElement<1>::interact_<0,1>(unsigned int i)
581 { ASSERT_LT_DBG(i, RefElement<1>::n_lines).error("Index out of bounds.");
582  return line_nodes_[i];}
583 
584 /// For line i {0,1,2}, give me indices of its nodes.
585 template<> template<> inline const IdxVector<2> RefElement<2>::interact_<0,1>(unsigned int i)
586 { ASSERT_LT_DBG(i, RefElement<2>::n_lines).error("Index out of bounds.");
587  return line_nodes_[i];}
588 
589 /// For line i {0,1,2,3,4,5}, give me indices of its nodes.
590 template<> template<> inline const IdxVector<2> RefElement<3>::interact_<0,1>(unsigned int i)
591 { ASSERT_LT_DBG(i, RefElement<3>::n_lines).error("Index out of bounds.");
592  return line_nodes_[i];}
593 
594 /// For node i {0,1}, give me indices of lines.
595 template<> template<> inline const IdxVector<1> RefElement<1>::interact_<1,0>(unsigned int i)
596 { ASSERT_LT_DBG(i, RefElement<1>::n_nodes).error("Index out of bounds.");
597  return node_lines_[i];}
598 
599 /// For node i {0,1,2}, give me indices of lines.
600 template<> template<> inline const IdxVector<2> RefElement<2>::interact_<1,0>(unsigned int i)
601 { ASSERT_LT_DBG(i, RefElement<2>::n_nodes).error("Index out of bounds.");
602  return node_lines_[i];}
603 
604 /// For node i {0,1,2,3}, give me indices of lines.
605 template<> template<> inline const IdxVector<3> RefElement<3>::interact_<1,0>(unsigned int i)
606 { ASSERT_LT_DBG(i, RefElement<3>::n_nodes).error("Index out of bounds.");
607  return node_lines_[i];}
608 
609 /// For side i {0,1,2}, give me indices of its nodes.
610 template<> template<> inline const IdxVector<3> RefElement<3>::interact_<0,2>(unsigned int i)
611 { ASSERT_LT_DBG(i, RefElement<3>::n_sides).error("Index out of bounds.");
612  return side_nodes_[i];}
613 
614 /// For node i {0,1,2,3}, give me indices of sides.
615 template<> template<> inline const IdxVector<3> RefElement<3>::interact_<2,0>(unsigned int i)
616 { ASSERT_LT_DBG(i, RefElement<3>::n_sides).error("Index out of bounds.");
617  return node_sides_[i];}
618 
619 /// For line i {0,1,2,3}, give me indices of sides.
620 template<> template<> inline const IdxVector<2> RefElement<3>::interact_<2,1>(unsigned int i)
621 { ASSERT_LT_DBG(i, RefElement<3>::n_lines).error("Index out of bounds.");
622  return line_sides_[i];}
623 
624 /// For side i {0,1,2}, give me indices of its lines.
625 template<> template<> inline const IdxVector<3> RefElement<3>::interact_<1,2>(unsigned int i)
626 { ASSERT_LT_DBG(i, RefElement<3>::n_sides).error("Index out of bounds.");
627  return side_lines_[i];}
628 
629 template<unsigned int dim> template<unsigned int OutDim, unsigned int InDim>
630 inline const IdxVector< (InDim>OutDim ? InDim+1 : dim-InDim) > RefElement<dim>::interact_(unsigned int i)
631 {
632  ASSERT(false)(dim)(OutDim)(InDim)(i).error("Not implemented.");
633  //ASSERT_LT_DBG(OutDim, dim);
634  //ASSERT_LT_DBG(InDim, dim);
635  return IdxVector< (InDim>OutDim ? InDim+1 : dim-InDim) >(); // just to avoid warning for missing return
636 }
637 
638 
639 template<unsigned int dim>
640 template < template <unsigned int OutDim, unsigned int InDim> class TInteraction, unsigned int OutDim, unsigned int InDim>
641 inline const IdxVector< (InDim>OutDim ? InDim+1 : dim-InDim) > RefElement<dim>::interact( TInteraction<OutDim,InDim> interaction )
642 {
643  return interact_<OutDim,InDim>(interaction.i_);
644 }
645 
646 #endif /* REF_ELEMENT_HH_ */
RefElement::normal_vector
static LocalPoint normal_vector(unsigned int sid)
RefElement::permutation_index
static unsigned int permutation_index(unsigned int p[n_nodes_per_side])
Definition: ref_element.cc:552
RefElement::FaceBaryPoint
Armor::ArmaVec< double, dim > FaceBaryPoint
Definition: ref_element.hh:174
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:571
RefElement
Definition: ref_element.hh:163
RefElement::normal_orientation
static unsigned int normal_orientation(unsigned int sid)
Definition: ref_element.cc:334
RefElement::BaryPoint
Armor::ArmaVec< double, dim+1 > BaryPoint
Definition: ref_element.hh:173
RefElement::LocalPoint
arma::vec::fixed< dim > LocalPoint
Definition: ref_element.hh:166
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:227
RefElement::bary_coords
static arma::mat::fixed< dim+1, subdim+1 > bary_coords(unsigned int sid)
Definition: ref_element.hh:475
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:401
Interaction::i_
unsigned int i_
Definition: ref_element.hh:159
RefElement::n_nodes
static const unsigned int n_nodes
Number of nodes.
Definition: ref_element.hh:223
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:412
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:404
RefElement::n_sides
static const unsigned int n_sides
Number of sides.
Definition: ref_element.hh:222
RefElement::make_bary_unit_vec
static BarycentricUnitVec make_bary_unit_vec()
Definition: ref_element.cc:436
RefElement::centers_of_subelements
static CentersList centers_of_subelements(unsigned int sub_dim)
Definition: ref_element.cc:484
RefElement::nodes_of_subelements
static const std::vector< std::vector< std::vector< unsigned int > > > nodes_of_subelements
Definition: ref_element.hh:258
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:402
RefElement::bary_to_local
static LocalPoint bary_to_local(const BaryPoint &bp)
Converts from barycentric to local coordinates.
Definition: ref_element.cc:319
RefElement::n_nodes_per_line
static const unsigned int n_nodes_per_line
Number of nodes in one line.
Definition: ref_element.hh:226
RefElement::n_sides_per_node
static const unsigned int n_sides_per_node
Number of sides with one common line.
Definition: ref_element.hh:228
RefElement::clip
static BaryPoint clip(const BaryPoint &barycentric)
Definition: ref_element.cc:450
IdxVector
std::array< unsigned int, Size > IdxVector
Definition: ref_element.hh:147
RefElement::topology_zeros_
static const IdxVector<(n_lines > n_nodes) ? n_lines :n_nodes > topology_zeros_[dim+1]
Definition: ref_element.hh:414
RefElement::oposite_node
static unsigned int oposite_node(unsigned int sid)
Definition: ref_element.cc:327
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:309
tolerance
Definition: isotherm.hh:67
RefElement::DECLARE_EXCEPTION
DECLARE_EXCEPTION(ExcInvalidPermutation,<< "Side permutation not found.\n")
Interaction
Definition: ref_element.hh:157
RefElement::barycentric_on_face
static FaceBaryPoint barycentric_on_face(const BaryPoint &barycentric, unsigned int i_face)
Definition: ref_element.cc:397
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:399
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:234
RefElement::node_coords
static LocalPoint node_coords(unsigned int nid)
Definition: ref_element.hh:490
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:231
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:158
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:305
RefElement::BarycentricUnitVec
std::vector< BaryPoint > BarycentricUnitVec
Definition: ref_element.hh:294
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:403
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:225
RefElement::n_side_permutations
static constexpr unsigned int n_side_permutations
Definition: ref_element.hh:268
RefElement::n_nodes_per_side
static const unsigned int n_nodes_per_side
Number of nodes on one side.
Definition: ref_element.hh:224
Armor::ArmaVec
typename arma::Col< Type >::template fixed< nr > ArmaVec
Definition: armor.hh:505
RefElement::side_permutations
static const std::vector< std::vector< unsigned int > > side_permutations
Definition: ref_element.hh:274
RefElement::topology_idx
static unsigned int topology_idx(unsigned int zeros_positions)
Definition: ref_element.hh:564
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:400
Generated by   doxygen 1.8.17