Flow123d  release_3.0.0-922-g9fcbba1
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 
127 template<> const unsigned int RefElement<1>::side_permutations[][n_nodes_per_side] = { { 0 } };
128 
129 template<> const unsigned int RefElement<2>::side_permutations[][n_nodes_per_side] = { { 0, 1 }, { 1, 0 } };
130 
131 template<> const unsigned int RefElement<3>::side_permutations[][n_nodes_per_side] = {
132  { 0, 1, 2 },
133  { 0, 2, 1 },
134  { 1, 0, 2 },
135  { 1, 2, 0 },
136  { 2, 0, 1 },
137  { 2, 1, 0 }
138 };
139 
140 
142  {(1 << 1), //node 0
143  (1 << 0)}, //node 1
144  {0, //the element
145  0}
146 };
147 
148 
150  {(1 << 1) | (1 << 2), //node 0
151  (1 << 0) | (1 << 2), //node 1
152  (1 << 0) | (1 << 1)}, //node 2
153  {(1 << 2), //line 0
154  (1 << 1), //line 1
155  (1 << 0)}, //line 2
156  {0, //the element
157  0,
158  0}
159 };
160 
162  {(1 << 1) | (1 << 2) | (1 << 3), //node 0
163  (1 << 0) | (1 << 2) | (1 << 3), //node 1
164  (1 << 0) | (1 << 1) | (1 << 3), //node 2
165  (1 << 0) | (1 << 1) | (1 << 2), //node 3
166  0,
167  0},
168  {(1 << 2) | (1 << 3), //line 0
169  (1 << 1) | (1 << 3), //line 1
170  (1 << 1) | (1 << 2), //line 2
171  (1 << 0) | (1 << 3), //line 3
172  (1 << 0) | (1 << 2), //line 4
173  (1 << 0) | (1 << 1)}, //line 5
174  {1 << 3, //side 0
175  1 << 2, //side 1
176  1 << 1, //side 2
177  1 << 0, //side 3
178  0,
179  0},
180  {0, //the element
181  0,
182  0,
183  0,
184  0,
185  0}
186 };
187 
188 
189 
190 
191 // template<> const unsigned int RefElement<1>::side_nodes[][1] = {
192 // { 0 },
193 // { 1 }
194 // };
195 //
196 // template<> const unsigned int RefElement<2>::side_nodes[][2] = {
197 // { 0, 1 },
198 // { 0, 2 },
199 // { 1, 2 }
200 // };
201 //
202 // template<> const unsigned int RefElement<3>::side_nodes[][3] = {
203 // { 0, 1, 2 },
204 // { 0, 1, 3 },
205 // { 0, 2, 3 },
206 // { 1, 2, 3 }
207 // };
208 //
209 //
210 //
211 // template<> const unsigned int RefElement<3>::side_lines[][3] = {
212 // {0,1,2},
213 // {0,3,4},
214 // {1,3,5},
215 // {2,4,5}
216 // };
217 //
218 //
219 // template<> const unsigned int RefElement<1>::line_nodes[][2] = {
220 // {0,1}
221 // };
222 //
223 // template<> const unsigned int RefElement<2>::line_nodes[][2] = {
224 // {0,1},
225 // {0,2},
226 // {1,2}
227 // };
228 //
229 // template<> const unsigned int RefElement<3>::line_nodes[][2] = {
230 // {0,1},
231 // {0,2},
232 // {1,2},
233 // {0,3},
234 // {1,3},
235 // {2,3}
236 // };
237 //
238 //
239 // /**
240 // * Indexes of sides for each line - with right orientation
241 // */
242 //
243 // template<> const unsigned int RefElement<3>::line_sides[][2] = {
244 // {0,1},
245 // {2,0},
246 // {0,3},
247 // {1,2},
248 // {3,1},
249 // {2,3}
250 // };
251 //
252 // /**
253 // * Indexes of sides for each line - for Simplex<2>, with right orientation
254 // */
255 // template<> const unsigned int RefElement<2>::line_sides[][2] = {
256 // {1,0},
257 // {0,2},
258 // {2,1}
259 // };
260 
261 
262 // 0: nodes of nodes
263 // 1: nodes of lines
264 // 2: nodes of sides
265 // 3: nodes of tetrahedron
267  { {0},{1} },
268  _array_to_vec(line_nodes_, n_lines)
269 };
270 
272  { {0}, {1}, {2} },
273  _array_to_vec(line_nodes_, n_lines),
274  { {0,1,2} }
275 };
276 
278  { {0}, {1}, {2}, {3} },
279  _array_to_vec(line_nodes_, n_lines),
280  _array_to_vec(side_nodes_, n_sides),
281  { {0,1,2,3} }
282 };
283 
284 
285 template<unsigned int dim>
287 {
288  ASSERT_EQ_DBG(lp.n_rows, dim);
289  BaryPoint bp;
290  bp.rows(1, dim ) = lp;
291  bp( 0 ) = 1.0 - arma::sum(lp);
292  return bp;
293 
294  // new armadillo
295  // return arma::join_col( arma::vec::fixed<1>( { 1.0 - arma::sum( local )} ), local);
296 
297 }
298 
299 template<unsigned int dim>
301 {
302  ASSERT_EQ_DBG(bp.n_rows, dim+1);
303  LocalPoint lp = bp.rows(1, dim);
304  return lp;
305 }
306 
307 template<unsigned int dim>
308 inline unsigned int RefElement<dim>::oposite_node(unsigned int sid)
309 {
310  return n_sides - sid - 1;
311 }
312 
313 
314 template<unsigned int dim>
315 unsigned int RefElement<dim>::normal_orientation(unsigned int sid)
316 {
317  ASSERT_LT_DBG(sid, n_sides).error("Side number is out of range!");
318 
319  return sid % 2;
320 }
321 
322 
323 template<>
324 vec::fixed<1> RefElement<1>::normal_vector(unsigned int sid)
325 {
326  ASSERT_LT_DBG(sid, n_sides).error("Side number is out of range!");
327 
328  return node_coords(sid) - node_coords(1-sid);
329 }
330 
331 template<>
332 vec::fixed<2> RefElement<2>::normal_vector(unsigned int sid)
333 {
334  ASSERT_LT_DBG(sid, n_sides).error("Side number is out of range!");
335  vec::fixed<2> barycenter, bar_side, n, t;
336 
337  // tangent vector along line
338  t = node_coords(line_nodes_[sid][1]) - node_coords(line_nodes_[sid][0]);
339  // barycenter coordinates
340  barycenter.fill(1./3);
341  // vector from barycenter to the side
342  bar_side = node_coords(line_nodes_[sid][0]) - barycenter;
343  // normal vector to side (modulo sign)
344  n(0) = -t(1);
345  n(1) = t(0);
346  n /= norm(n,2);
347  // check sign of normal vector
348  if (dot(n,bar_side) < 0) n *= -1;
349 
350  return n;
351 }
352 
353 template<>
354 vec::fixed<3> RefElement<3>::normal_vector(unsigned int sid)
355 {
356  ASSERT_LT_DBG(sid, n_sides).error("Side number is out of range!");
357  vec::fixed<3> barycenter, bar_side, n, t1, t2;
358 
359  // tangent vectors of side
360  t1 = node_coords(side_nodes_[sid][1]) - node_coords(side_nodes_[sid][0]);
361  t2 = node_coords(side_nodes_[sid][2]) - node_coords(side_nodes_[sid][0]);
362  // baryucenter coordinates
363  barycenter.fill(0.25);
364  // vector from barycenter to the side
365  bar_side = node_coords(side_nodes_[sid][0]) - barycenter;
366  // normal vector (modulo sign)
367  n = cross(t1,t2);
368  n /= norm(n,2);
369  // check sign of normal vector
370  if (dot(n,bar_side) < 0) n = -n;
371 
372  return n;
373 }
374 
375 
376 
377 template<unsigned int dim>
378 auto RefElement<dim>::barycentric_on_face(const BaryPoint &barycentric, unsigned int i_face) -> FaceBaryPoint
379 {
380  ASSERT_EQ_DBG(barycentric.n_rows, dim+1);
381  FaceBaryPoint face_barycentric;
382  for(unsigned int i=0; i < dim; i++) {
383 // unsigned int i_sub_node = (i+1)%dim;
384 // unsigned int i_bary = (dim + side_nodes_[i_face][i_sub_node])%(dim+1);
385  unsigned int i_bary = interact_<0,dim-1>(i_face)[i];
386  face_barycentric[i] = barycentric[ i_bary ];
387  }
388  return face_barycentric;
389 }
390 
391 
392 template<unsigned int dim>
393 std::pair<unsigned int, unsigned int> RefElement<dim>::zeros_positions(const BaryPoint &barycentric,
394  double tolerance)
395 {
396  unsigned int zeros = 0;
397  unsigned int n_zeros = 0;
398  for(char i=0; i < dim+1; i++){
399  if(std::fabs(barycentric[i]) < tolerance)
400  {
401  zeros = zeros | (1 << i);
402  n_zeros++;
403  }
404  }
405 
406  return std::make_pair(n_zeros, zeros);
407 }
408 
409 
410 template<>
411 auto RefElement<0>::clip(const BaryPoint &barycentric) -> BaryPoint
412 {
413  return barycentric;
414 }
415 
416 template<unsigned int dim>
417 auto RefElement<dim>::make_bary_unit_vec()->BarycentricUnitVec
418 {
419  std::vector<arma::vec::fixed<dim+1> > bary_unit_vec(dim+1, arma::zeros(dim+1));
420  for(unsigned int i=0; i<dim; i++) {
421  bary_unit_vec[i][i] = 1.0;
422  bary_unit_vec[i][dim] = -1.0;
423  bary_unit_vec[dim][i] = -1.0 / dim;
424  }
425  bary_unit_vec[dim][dim] = 1.0;
426  return bary_unit_vec;
427 }
428 
429 
430 template<unsigned int dim>
431 auto RefElement<dim>::clip(const BaryPoint &barycentric) -> BaryPoint
432 {
433  static BarycentricUnitVec bary_unit_vec = make_bary_unit_vec();
434  ASSERT_EQ_DBG(barycentric.n_rows, dim+1);
435  for(unsigned int i_bary=0; i_bary < dim +1; i_bary ++) {
436  if (barycentric[i_bary] < 0.0) {
437  // index of barycentric coord that is constant on the face i_side
438  // as we use barycentric coords starting with local coordinates:
439  // TODO: rather work only with local coords and/or with canonical barycentric coords
440  unsigned int i_side = (dim - i_bary);
441  // project to face
442  arma::vec projection_to_face(dim+1);
443  //barycentric.print(cout, "input");
444  //cout << "is: " << i_side << endl;
445  //cout << "ibary: " << i_bary << endl;
446  //bary_unit_vec[i_bary].print(cout, "normal");
447  //barycentric.subvec(0, dim-1).print(cout, "bary sub");
448  projection_to_face = barycentric - barycentric[i_bary]*bary_unit_vec[i_bary];
449  //projection_to_face(dim) = 1.0 - arma::sum(projection_to_face.subvec(0, dim-1));
450  //projection_to_face.print(cout, "projection");
451  auto bary_on_face = barycentric_on_face(projection_to_face, i_side);
452  //bary_on_face.print(cout, "b on f");
453  auto sub_clip = RefElement<dim-1>::clip(bary_on_face);
454  //sub_clip.print(cout, "sub clip");
455  return interpolate<dim-1>(sub_clip, i_side);
456  }
457  }
458  return barycentric;
459 
460 }
461 
462 
463 
464 template<unsigned int dim>
465 auto RefElement<dim>::centers_of_subelements(unsigned int sub_dim)->CentersList
466 {
468  if (list.size() == 0) {
469  list.resize(dim+1);
470  for(unsigned int sdim=0; sdim < dim+1; sdim++) {
471  // Temporary solution until we have unified interface to
472  // the internal indexing.
473  // We use the fact that numbering of subelements goes as numbering of
474  // k combinations over nodes.
475 // std::vector<unsigned int> subel_comb(sdim+2);
476  for(auto &sub_el_nodes : nodes_of_subelements[sdim]) {
477  ASSERT_EQ_DBG(sub_el_nodes.size(), sdim+1);
478  LocalPoint center = arma::zeros(dim);
479  for( unsigned int i_node : sub_el_nodes)
480  center+=node_coords( i_node );
481  center/=(sdim+1);
482  list[sdim].push_back(center);
483  }
484  }
485  }
486 
487  ASSERT_LE_DBG(sub_dim, dim);
488  return list[sub_dim];
489 }
490 
491 
492 template<>
493 double RefElement<1>::side_measure(unsigned int sid)
494 {
495  ASSERT_LT_DBG(sid, n_sides).error("Side number is out of range!");
496 
497  return 1;
498 }
499 
500 
501 template<>
502 double RefElement<2>::side_measure(unsigned int sid)
503 {
504  ASSERT_LT_DBG(sid, n_sides).error("Side number is out of range!");
505 
506  return norm(node_coords(line_nodes_[sid][1]) - node_coords(line_nodes_[sid][0]),2);
507 }
508 
509 
510 template<>
511 double RefElement<3>::side_measure(unsigned int sid)
512 {
513  ASSERT_LT_DBG(sid, n_sides).error("Side number is out of range!");
514 
515  return 0.5*norm(cross(node_coords(side_nodes_[sid][1]) - node_coords(side_nodes_[sid][0]),
516  node_coords(side_nodes_[sid][2]) - node_coords(side_nodes_[sid][0])),2);
517 }
518 
519 
520 
521 template <>
522 unsigned int RefElement<3>::line_between_faces(unsigned int f1, unsigned int f2) {
523  unsigned int i,j;
524  i=j=0;
525  while (side_lines_[f1][i] != side_lines_[f2][j])
526  if (side_lines_[f1][i] < side_lines_[f2][j]) i++;
527  else j++;
528  return side_lines_[f1][i];
529 }
530 
531 
532 template<unsigned int dim>
533 unsigned int RefElement<dim>::permutation_index(unsigned int p[n_nodes_per_side])
534 {
535  unsigned int index;
536  for (index = 0; index < n_side_permutations; index++)
537  if (equal(p, p + n_nodes_per_side, side_permutations[index]))
538  return index;
539 
540  xprintf(PrgErr, "Side permutation not found.\n");
541 
542  // The following line is present in order to suppress compilers warning
543  // about missing return value.
544  return 0;
545 }
546 
547 
548 /**
549  * Basic line interpolation
550  */
551 template<unsigned int dim>
553  arma::vec::fixed<dim+1> first_coords,
554  arma::vec::fixed<dim+1> second_coords,
555  double first_theta, double second_theta, double theta){
556 
557  arma::vec::fixed<dim+1> bary_interpolated_coords;
558  bary_interpolated_coords = ((theta - first_theta) * second_coords + (second_theta - theta) * first_coords)
559  /(second_theta - first_theta);
560  return bary_interpolated_coords;
561 };
562 
563 template class RefElement<1>;
564 template class RefElement<2>;
565 template class RefElement<3>;
566 
567 
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:552
static unsigned int permutation_index(unsigned int p[n_nodes_per_side])
Definition: ref_element.cc:533
arma::vec::fixed< dim > LocalPoint
Definition: ref_element.hh:165
static unsigned int line_between_faces(unsigned int f1, unsigned int f2)
static std::pair< unsigned int, unsigned int > zeros_positions(const BaryPoint &barycentric, double tolerance=std::numeric_limits< double >::epsilon()*2)
Definition: ref_element.cc:393
static unsigned int normal_orientation(unsigned int sid)
Definition: ref_element.cc:315
static unsigned int oposite_node(unsigned int sid)
Definition: ref_element.cc:308
static CentersList centers_of_subelements(unsigned int sub_dim)
Definition: ref_element.cc:465
static LocalPoint bary_to_local(const BaryPoint &bp)
Converts from barycentric to local coordinates.
Definition: ref_element.cc:300
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:378
#define xprintf(...)
Definition: system.hh:92
#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:286
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc...
static BarycentricUnitVec make_bary_unit_vec()
Definition: ref_element.cc:417
static BaryPoint clip(const BaryPoint &barycentric)
Definition: ref_element.cc:431
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