Flow123d  release_3.0.0-1211-g84d0b74
ref_element.cc
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.cc
15  * @brief Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc.
16  * @author Jan Stebel
17  */
18 
19 #include "system/global_defs.h"
20 #include "system/system.hh"
21 #include "mesh/ref_element.hh"
22 
23 
24 
25 using namespace arma;
26 using namespace std;
27 
28 template<std::size_t n>
31  for(unsigned int i=0; i<m; i++)
32  for(unsigned int j=0;j<n; j++)
33  vec[i].push_back(array[i][j]);
34  return vec;
35 }
36 
37 
38 // template<unsigned int n>
39 // std::vector< std::vector<unsigned int> > _array_to_vec( const unsigned int array[][n], unsigned int m) {
40 // std::vector< std::vector<unsigned int> > vec(m);
41 // for(unsigned int i=0; i<m; i++)
42 // for(unsigned int j=0;j<n; j++)
43 // vec[i].push_back(array[i][j]);
44 // return vec;
45 // }
46 
47 
48 
49 
51  {0,1}
52 };
53 
55  {0,1},
56  {0,2},
57  {1,2}
58 };
59 
61  {0,1}, //0
62  {0,2}, //1
63  {0,3}, //2 <-3 (fixed order)
64  {1,2}, //3 <-2
65  {1,3}, //4
66  {2,3} //5
67 };
68 
69 
71  {0},
72  {0}
73 };
74 
76  {1,0},
77  {0,2},
78  {2,1}
79 };
80 
81 // Order clockwise looking over the vertex to center; smallest index first
83  {0,1,2},
84  {0,4,3},
85  {1,3,5},
86  {2,5,4},
87 };
88 
89 
90 
91 // Lexicographic order.
93  { 0, 1, 2 },
94  { 0, 1, 3 },
95  { 0, 2, 3 },
96  { 1, 2, 3 }
97 };
98 
99 // Order clockwise, faces opposite to the lines from node_lines.
101  { 2, 1, 0 },
102  { 3, 0, 1 },
103  { 3, 2, 0 },
104  { 3, 1, 2 }
105 };
106 
107 // faces ordered clock wise with respect to edge shifted to center of tetrahedron
109  {0,1},
110  {2,0},
111  {1,2},
112  {0,3},
113  {3,1},
114  {2,3}
115 };
116 
117 
119  {0,1,3},
120  {0,2,4},
121  {1,2,5},
122  {3,4,5}
123 };
124 
125 
126 template<> const unsigned int RefElement<0>::side_permutations[][n_nodes_per_side] = { };
127 
128 template<> const unsigned int RefElement<1>::side_permutations[][n_nodes_per_side] = { { 0 } };
129 
130 template<> const unsigned int RefElement<2>::side_permutations[][n_nodes_per_side] = { { 0, 1 }, { 1, 0 } };
131 
132 template<> const unsigned int RefElement<3>::side_permutations[][n_nodes_per_side] = {
133  { 0, 1, 2 },
134  { 0, 2, 1 },
135  { 1, 0, 2 },
136  { 1, 2, 0 },
137  { 2, 0, 1 },
138  { 2, 1, 0 }
139 };
140 
141 
143  {(1 << 1)} //node 0
144 };
145 
146 
148  {(1 << 1), //node 0
149  (1 << 0)}, //node 1
150  {0, //the element
151  0}
152 };
153 
154 
156  {(1 << 1) | (1 << 2), //node 0
157  (1 << 0) | (1 << 2), //node 1
158  (1 << 0) | (1 << 1)}, //node 2
159  {(1 << 2), //line 0
160  (1 << 1), //line 1
161  (1 << 0)}, //line 2
162  {0, //the element
163  0,
164  0}
165 };
166 
168  {(1 << 1) | (1 << 2) | (1 << 3), //node 0
169  (1 << 0) | (1 << 2) | (1 << 3), //node 1
170  (1 << 0) | (1 << 1) | (1 << 3), //node 2
171  (1 << 0) | (1 << 1) | (1 << 2), //node 3
172  0,
173  0},
174  {(1 << 2) | (1 << 3), //line 0
175  (1 << 1) | (1 << 3), //line 1
176  (1 << 1) | (1 << 2), //line 2
177  (1 << 0) | (1 << 3), //line 3
178  (1 << 0) | (1 << 2), //line 4
179  (1 << 0) | (1 << 1)}, //line 5
180  {1 << 3, //side 0
181  1 << 2, //side 1
182  1 << 1, //side 2
183  1 << 0, //side 3
184  0,
185  0},
186  {0, //the element
187  0,
188  0,
189  0,
190  0,
191  0}
192 };
193 
194 
195 
196 
197 // template<> const unsigned int RefElement<1>::side_nodes[][1] = {
198 // { 0 },
199 // { 1 }
200 // };
201 //
202 // template<> const unsigned int RefElement<2>::side_nodes[][2] = {
203 // { 0, 1 },
204 // { 0, 2 },
205 // { 1, 2 }
206 // };
207 //
208 // template<> const unsigned int RefElement<3>::side_nodes[][3] = {
209 // { 0, 1, 2 },
210 // { 0, 1, 3 },
211 // { 0, 2, 3 },
212 // { 1, 2, 3 }
213 // };
214 //
215 //
216 //
217 // template<> const unsigned int RefElement<3>::side_lines[][3] = {
218 // {0,1,2},
219 // {0,3,4},
220 // {1,3,5},
221 // {2,4,5}
222 // };
223 //
224 //
225 // template<> const unsigned int RefElement<1>::line_nodes[][2] = {
226 // {0,1}
227 // };
228 //
229 // template<> const unsigned int RefElement<2>::line_nodes[][2] = {
230 // {0,1},
231 // {0,2},
232 // {1,2}
233 // };
234 //
235 // template<> const unsigned int RefElement<3>::line_nodes[][2] = {
236 // {0,1},
237 // {0,2},
238 // {1,2},
239 // {0,3},
240 // {1,3},
241 // {2,3}
242 // };
243 //
244 //
245 // /**
246 // * Indexes of sides for each line - with right orientation
247 // */
248 //
249 // template<> const unsigned int RefElement<3>::line_sides[][2] = {
250 // {0,1},
251 // {2,0},
252 // {0,3},
253 // {1,2},
254 // {3,1},
255 // {2,3}
256 // };
257 //
258 // /**
259 // * Indexes of sides for each line - for Simplex<2>, with right orientation
260 // */
261 // template<> const unsigned int RefElement<2>::line_sides[][2] = {
262 // {1,0},
263 // {0,2},
264 // {2,1}
265 // };
266 
267 
268 // 0: nodes of nodes
269 // 1: nodes of lines
270 // 2: nodes of sides
271 // 3: nodes of tetrahedron
273  { {0} }
274 };
275 
277  { {0},{1} },
278  _array_to_vec(line_nodes_, n_lines)
279 };
280 
282  { {0}, {1}, {2} },
283  _array_to_vec(line_nodes_, n_lines),
284  { {0,1,2} }
285 };
286 
288  { {0}, {1}, {2}, {3} },
289  _array_to_vec(line_nodes_, n_lines),
290  _array_to_vec(side_nodes_, n_sides),
291  { {0,1,2,3} }
292 };
293 
294 
295 template<unsigned int dim>
297 {
298  ASSERT_EQ_DBG(lp.n_rows, dim);
299  BaryPoint bp;
300  bp.rows(1, dim ) = lp;
301  bp( 0 ) = 1.0 - arma::sum(lp);
302  return bp;
303 
304  // new armadillo
305  // return arma::join_col( arma::vec::fixed<1>( { 1.0 - arma::sum( local )} ), local);
306 
307 }
308 
309 template<unsigned int dim>
311 {
312  ASSERT_EQ_DBG(bp.n_rows, dim+1);
313  LocalPoint lp = bp.rows(1, dim);
314  return lp;
315 }
316 
317 template<unsigned int dim>
318 inline unsigned int RefElement<dim>::oposite_node(unsigned int sid)
319 {
320  return n_sides - sid - 1;
321 }
322 
323 
324 template<unsigned int dim>
325 unsigned int RefElement<dim>::normal_orientation(unsigned int sid)
326 {
327  ASSERT_LT_DBG(sid, n_sides).error("Side number is out of range!");
328 
329  return sid % 2;
330 }
331 
332 
333 template<>
334 vec::fixed<1> RefElement<1>::normal_vector(unsigned int sid)
335 {
336  ASSERT_LT_DBG(sid, n_sides).error("Side number is out of range!");
337 
338  return node_coords(sid) - node_coords(1-sid);
339 }
340 
341 template<>
342 vec::fixed<2> RefElement<2>::normal_vector(unsigned int sid)
343 {
344  ASSERT_LT_DBG(sid, n_sides).error("Side number is out of range!");
345  vec::fixed<2> barycenter, bar_side, n, t;
346 
347  // tangent vector along line
348  t = node_coords(line_nodes_[sid][1]) - node_coords(line_nodes_[sid][0]);
349  // barycenter coordinates
350  barycenter.fill(1./3);
351  // vector from barycenter to the side
352  bar_side = node_coords(line_nodes_[sid][0]) - barycenter;
353  // normal vector to side (modulo sign)
354  n(0) = -t(1);
355  n(1) = t(0);
356  n /= norm(n,2);
357  // check sign of normal vector
358  if (dot(n,bar_side) < 0) n *= -1;
359 
360  return n;
361 }
362 
363 template<>
364 vec::fixed<3> RefElement<3>::normal_vector(unsigned int sid)
365 {
366  ASSERT_LT_DBG(sid, n_sides).error("Side number is out of range!");
367  vec::fixed<3> barycenter, bar_side, n, t1, t2;
368 
369  // tangent vectors of side
370  t1 = node_coords(side_nodes_[sid][1]) - node_coords(side_nodes_[sid][0]);
371  t2 = node_coords(side_nodes_[sid][2]) - node_coords(side_nodes_[sid][0]);
372  // baryucenter coordinates
373  barycenter.fill(0.25);
374  // vector from barycenter to the side
375  bar_side = node_coords(side_nodes_[sid][0]) - barycenter;
376  // normal vector (modulo sign)
377  n = cross(t1,t2);
378  n /= norm(n,2);
379  // check sign of normal vector
380  if (dot(n,bar_side) < 0) n = -n;
381 
382  return n;
383 }
384 
385 
386 
387 template<unsigned int dim>
388 auto RefElement<dim>::barycentric_on_face(const BaryPoint &barycentric, unsigned int i_face) -> FaceBaryPoint
389 {
390  ASSERT_EQ_DBG(barycentric.n_rows, dim+1);
391  FaceBaryPoint face_barycentric;
392  for(unsigned int i=0; i < dim; i++) {
393 // unsigned int i_sub_node = (i+1)%dim;
394 // unsigned int i_bary = (dim + side_nodes_[i_face][i_sub_node])%(dim+1);
395  unsigned int i_bary = interact_<0,dim-1>(i_face)[i];
396  face_barycentric[i] = barycentric[ i_bary ];
397  }
398  return face_barycentric;
399 }
400 
401 
402 template<unsigned int dim>
403 std::pair<unsigned int, unsigned int> RefElement<dim>::zeros_positions(const BaryPoint &barycentric,
404  double tolerance)
405 {
406  unsigned int zeros = 0;
407  unsigned int n_zeros = 0;
408  for(char i=0; i < dim+1; i++){
409  if(std::fabs(barycentric[i]) < tolerance)
410  {
411  zeros = zeros | (1 << i);
412  n_zeros++;
413  }
414  }
415 
416  return std::make_pair(n_zeros, zeros);
417 }
418 
419 
420 template<>
421 auto RefElement<0>::clip(const BaryPoint &barycentric) -> BaryPoint
422 {
423  return barycentric;
424 }
425 
426 template<unsigned int dim>
427 auto RefElement<dim>::make_bary_unit_vec()->BarycentricUnitVec
428 {
429  std::vector<arma::vec::fixed<dim+1> > bary_unit_vec(dim+1, arma::zeros(dim+1));
430  for(unsigned int i=0; i<dim; i++) {
431  bary_unit_vec[i][i] = 1.0;
432  bary_unit_vec[i][dim] = -1.0;
433  bary_unit_vec[dim][i] = -1.0 / dim;
434  }
435  bary_unit_vec[dim][dim] = 1.0;
436  return bary_unit_vec;
437 }
438 
439 
440 template<unsigned int dim>
441 auto RefElement<dim>::clip(const BaryPoint &barycentric) -> BaryPoint
442 {
443  static BarycentricUnitVec bary_unit_vec = make_bary_unit_vec();
444  ASSERT_EQ_DBG(barycentric.n_rows, dim+1);
445  for(unsigned int i_bary=0; i_bary < dim +1; i_bary ++) {
446  if (barycentric[i_bary] < 0.0) {
447  // index of barycentric coord that is constant on the face i_side
448  // as we use barycentric coords starting with local coordinates:
449  // TODO: rather work only with local coords and/or with canonical barycentric coords
450  unsigned int i_side = (dim - i_bary);
451  // project to face
452  arma::vec projection_to_face(dim+1);
453  //barycentric.print(cout, "input");
454  //cout << "is: " << i_side << endl;
455  //cout << "ibary: " << i_bary << endl;
456  //bary_unit_vec[i_bary].print(cout, "normal");
457  //barycentric.subvec(0, dim-1).print(cout, "bary sub");
458  projection_to_face = barycentric - barycentric[i_bary]*bary_unit_vec[i_bary];
459  //projection_to_face(dim) = 1.0 - arma::sum(projection_to_face.subvec(0, dim-1));
460  //projection_to_face.print(cout, "projection");
461  auto bary_on_face = barycentric_on_face(projection_to_face, i_side);
462  //bary_on_face.print(cout, "b on f");
463  auto sub_clip = RefElement<dim-1>::clip(bary_on_face);
464  //sub_clip.print(cout, "sub clip");
465  return interpolate<dim-1>(sub_clip, i_side);
466  }
467  }
468  return barycentric;
469 
470 }
471 
472 
473 
474 template<unsigned int dim>
475 auto RefElement<dim>::centers_of_subelements(unsigned int sub_dim)->CentersList
476 {
478  if (list.size() == 0) {
479  list.resize(dim+1);
480  for(unsigned int sdim=0; sdim < dim+1; sdim++) {
481  // Temporary solution until we have unified interface to
482  // the internal indexing.
483  // We use the fact that numbering of subelements goes as numbering of
484  // k combinations over nodes.
485 // std::vector<unsigned int> subel_comb(sdim+2);
486  for(auto &sub_el_nodes : nodes_of_subelements[sdim]) {
487  ASSERT_EQ_DBG(sub_el_nodes.size(), sdim+1);
488  LocalPoint center = arma::zeros(dim);
489  for( unsigned int i_node : sub_el_nodes)
490  center+=node_coords( i_node );
491  center/=(sdim+1);
492  list[sdim].push_back(center);
493  }
494  }
495  }
496 
497  ASSERT_LE_DBG(sub_dim, dim);
498  return list[sub_dim];
499 }
500 
501 
502 template<>
503 double RefElement<1>::side_measure(unsigned int sid)
504 {
505  ASSERT_LT_DBG(sid, n_sides).error("Side number is out of range!");
506 
507  return 1;
508 }
509 
510 
511 template<>
512 double RefElement<2>::side_measure(unsigned int sid)
513 {
514  ASSERT_LT_DBG(sid, n_sides).error("Side number is out of range!");
515 
516  return norm(node_coords(line_nodes_[sid][1]) - node_coords(line_nodes_[sid][0]),2);
517 }
518 
519 
520 template<>
521 double RefElement<3>::side_measure(unsigned int sid)
522 {
523  ASSERT_LT_DBG(sid, n_sides).error("Side number is out of range!");
524 
525  return 0.5*norm(cross(node_coords(side_nodes_[sid][1]) - node_coords(side_nodes_[sid][0]),
526  node_coords(side_nodes_[sid][2]) - node_coords(side_nodes_[sid][0])),2);
527 }
528 
529 
530 
531 template <>
532 unsigned int RefElement<3>::line_between_faces(unsigned int f1, unsigned int f2) {
533  unsigned int i,j;
534  i=j=0;
535  while (side_lines_[f1][i] != side_lines_[f2][j])
536  if (side_lines_[f1][i] < side_lines_[f2][j]) i++;
537  else j++;
538  return side_lines_[f1][i];
539 }
540 
541 
542 template<unsigned int dim>
543 unsigned int RefElement<dim>::permutation_index(unsigned int p[n_nodes_per_side])
544 {
545  unsigned int index;
546  for (index = 0; index < n_side_permutations; index++)
547  if (equal(p, p + n_nodes_per_side, side_permutations[index]))
548  return index;
549 
550  xprintf(PrgErr, "Side permutation not found.\n");
551 
552  // The following line is present in order to suppress compilers warning
553  // about missing return value.
554  return 0;
555 }
556 
557 
558 /**
559  * Basic line interpolation
560  */
561 template<unsigned int dim>
563  arma::vec::fixed<dim+1> first_coords,
564  arma::vec::fixed<dim+1> second_coords,
565  double first_theta, double second_theta, double theta){
566 
567  arma::vec::fixed<dim+1> bary_interpolated_coords;
568  bary_interpolated_coords = ((theta - first_theta) * second_coords + (second_theta - theta) * first_coords)
569  /(second_theta - first_theta);
570  return bary_interpolated_coords;
571 };
572 
573 template class RefElement<0>;
574 template class RefElement<1>;
575 template class RefElement<2>;
576 template class RefElement<3>;
577 
578 
static LocalPoint normal_vector(unsigned int sid)
#define ASSERT_EQ_DBG(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:331
std::vector< std::vector< unsigned int > > _array_to_vec(const IdxVector< n > array[], unsigned int m)
Definition: ref_element.cc:29
static BaryPoint line_barycentric_interpolation(BaryPoint first_coords, BaryPoint second_coords, double first_theta, double second_theta, double theta)
Definition: ref_element.cc:562
static unsigned int permutation_index(unsigned int p[n_nodes_per_side])
Definition: ref_element.cc:543
arma::vec::fixed< dim > LocalPoint
Definition: ref_element.hh:165
static unsigned int line_between_faces(unsigned int f1, unsigned int f2)
Mat< double, N, 1 > vec
Definition: armor.hh:211
static std::pair< unsigned int, unsigned int > zeros_positions(const BaryPoint &barycentric, double tolerance=std::numeric_limits< double >::epsilon()*2)
Definition: ref_element.cc:403
static unsigned int normal_orientation(unsigned int sid)
Definition: ref_element.cc:325
static unsigned int oposite_node(unsigned int sid)
Definition: ref_element.cc:318
static CentersList centers_of_subelements(unsigned int sub_dim)
Definition: ref_element.cc:475
Array< double > array
Definition: armor.hh:216
static LocalPoint bary_to_local(const BaryPoint &bp)
Converts from barycentric to local coordinates.
Definition: ref_element.cc:310
Global macros to enhance readability and debugging, general constants.
static FaceBaryPoint barycentric_on_face(const BaryPoint &barycentric, unsigned int i_face)
Definition: ref_element.cc:388
#define xprintf(...)
Definition: system.hh:92
Type dot(const Mat< Type, nRows, nCols > &a, const Mat< Type, nRows, nCols > &b)
Definition: armor.hh:176
#define ASSERT_LE_DBG(a, b)
Definition of comparative assert macro (Less or Equal) only for debug mode.
Definition: asserts.hh:307
std::array< unsigned int, Size > IdxVector
Definition: ref_element.hh:146
static double side_measure(unsigned int sid)
Definition: system.hh:64
arma::vec::fixed< dim > FaceBaryPoint
Definition: ref_element.hh:173
static BaryPoint local_to_bary(const LocalPoint &lp)
Converts from local to barycentric coordinates.
Definition: ref_element.cc:296
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc...
static BarycentricUnitVec make_bary_unit_vec()
Definition: ref_element.cc:427
static BaryPoint clip(const BaryPoint &barycentric)
Definition: ref_element.cc:441
arma::vec::fixed< dim+1 > BaryPoint
Definition: ref_element.hh:172
#define ASSERT_LT_DBG(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:299