Flow123d  master-f44eb46
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 
31 
32 
33 
34 template<std::size_t n>
36  std::vector< std::vector<unsigned int> > vec(array_vec.size());
37  for(unsigned int i=0; i<array_vec.size(); i++)
38  for(unsigned int j=0;j<n; j++)
39  vec[i].push_back(array_vec[i][j]);
40  return vec;
41 }
42 
43 
44 
46  = { {0,1} };
47 
49  = { {0,1},
50  {0,2},
51  {1,2}};
52 
54  = { {0,1}, //0
55  {0,2}, //1
56  {0,3}, //2 <-3 (fixed order)
57  {1,2}, //3 <-2
58  {1,3}, //4
59  {2,3}}; //5
60 
61 
63  = { {0},
64  {0}};
65 
67  = { {1,0},
68  {0,2},
69  {2,1}};
70 
71 
72 
73 
74 
75 // Lexicographic order.
77  = { { 0, 1, 2 },
78  { 0, 1, 3 },
79  { 0, 2, 3 },
80  { 1, 2, 3 }};
81 
82 //// Order clockwise, faces opposite to the lines from node_lines.
83 //// !! dependes on S3 inversion
84 //template<> const std::vector<IdxVector<3>> RefElement<3>::node_sides_
85 // = { { { 2, 1, 0 },
86 // { 3, 0, 1 },
87 // { 3, 2, 0 },
88 // { 3, 1, 2 }},
89 // { { 2, 0, 1 },
90 // { 3, 1, 0 },
91 // { 3, 0, 2 },
92 // { 3, 2, 1 }};
93 
94 
95 
96 // lexicographic order
98  = { {0,1,3},
99  {0,2,4},
100  {1,2,5},
101  {3,4,5}};
102 
103 
105  {(1 << 1)} //node 0
106 };
107 
108 
110  {(1 << 1), //node 0
111  (1 << 0)}, //node 1
112  {0, //the element
113  0}
114 };
115 
116 
118  {(1 << 1) | (1 << 2), //node 0
119  (1 << 0) | (1 << 2), //node 1
120  (1 << 0) | (1 << 1)}, //node 2
121  {(1 << 2), //line 0
122  (1 << 1), //line 1
123  (1 << 0)}, //line 2
124  {0, //the element
125  0,
126  0}
127 };
128 
130  {(1 << 1) | (1 << 2) | (1 << 3), //node 0
131  (1 << 0) | (1 << 2) | (1 << 3), //node 1
132  (1 << 0) | (1 << 1) | (1 << 3), //node 2
133  (1 << 0) | (1 << 1) | (1 << 2), //node 3
134  0,
135  0},
136  {(1 << 2) | (1 << 3), //line 0
137  (1 << 1) | (1 << 3), //line 1
138  (1 << 1) | (1 << 2), //line 2
139  (1 << 0) | (1 << 3), //line 3
140  (1 << 0) | (1 << 2), //line 4
141  (1 << 0) | (1 << 1)}, //line 5
142  {1 << 3, //side 0
143  1 << 2, //side 1
144  1 << 1, //side 2
145  1 << 0, //side 3
146  0,
147  0},
148  {0, //the element
149  0,
150  0,
151  0,
152  0,
153  0}
154 };
155 
156 
157 
158 
159 // 0: nodes of nodes
160 // 1: nodes of lines
161 // 2: nodes of sides
162 // 3: nodes of tetrahedron
164  { {0} }
165 };
166 
168  { {0}, {1} },
169  _array_to_vec(line_nodes_)
170 };
171 
173  { {0}, {1}, {2} },
174  _array_to_vec(line_nodes_),
175  { {0,1,2} }
176 };
177 
179  { {0}, {1}, {2}, {3} },
180  _array_to_vec(line_nodes_),
181  _array_to_vec(side_nodes_),
182  { {0,1,2,3} }
183 };
184 
185 
186 template<unsigned int dim>
188 {
189  ASSERT_EQ(lp.n_rows, dim);
190  BaryPoint bp;
191  bp.rows(1, dim ) = lp;
192  bp( 0 ) = 1.0 - arma::sum(lp);
193  return bp;
194 
195  // new armadillo
196  // return arma::join_col( arma::vec::fixed<1>( { 1.0 - arma::sum( local )} ), local);
197 
198 }
199 
200 template<unsigned int dim>
202 {
203  ASSERT_EQ(bp.n_rows, dim+1);
204  LocalPoint lp = bp.rows(1, dim);
205  return lp;
206 }
207 
208 template<unsigned int dim>
209 inline unsigned int RefElement<dim>::oposite_node(unsigned int sid)
210 {
211  return n_sides - sid - 1;
212 }
213 
214 
215 template<unsigned int dim>
216 unsigned int RefElement<dim>::normal_orientation(unsigned int sid)
217 {
218  ASSERT_LT(sid, n_sides).error("Side number is out of range!");
219 
220  return sid % 2;
221 }
222 
223 
224 template<>
225 vec::fixed<1> RefElement<1>::normal_vector(unsigned int sid)
226 {
227  ASSERT_LT(sid, n_sides).error("Side number is out of range!");
228 
229  return node_coords(sid) - node_coords(1-sid);
230 }
231 
232 template<>
233 vec::fixed<2> RefElement<2>::normal_vector(unsigned int sid)
234 {
235  ASSERT_LT(sid, n_sides).error("Side number is out of range!");
236  vec::fixed<2> barycenter, bar_side, n, t;
237 
238  // tangent vector along line
239  t = node_coords(line_nodes_[sid][1]) - node_coords(line_nodes_[sid][0]);
240  // barycenter coordinates
241  barycenter.fill(1./3);
242  // vector from barycenter to the side
243  bar_side = node_coords(line_nodes_[sid][0]) - barycenter;
244  // normal vector to side (modulo sign)
245  n(0) = -t(1);
246  n(1) = t(0);
247  n /= norm(n,2);
248  // check sign of normal vector
249  if (dot(n,bar_side) < 0) n *= -1;
250 
251  return n;
252 }
253 
254 template<>
255 vec::fixed<3> RefElement<3>::normal_vector(unsigned int sid)
256 {
257  ASSERT_LT(sid, n_sides).error("Side number is out of range!");
258  vec::fixed<3> barycenter, bar_side, n, t1, t2;
259 
260  // tangent vectors of side
261  t1 = node_coords(side_nodes_[sid][1]) - node_coords(side_nodes_[sid][0]);
262  t2 = node_coords(side_nodes_[sid][2]) - node_coords(side_nodes_[sid][0]);
263  // baryucenter coordinates
264  barycenter.fill(0.25);
265  // vector from barycenter to the side
266  bar_side = node_coords(side_nodes_[sid][0]) - barycenter;
267  // normal vector (modulo sign)
268  n = cross(t1,t2);
269  n /= norm(n,2);
270  // check sign of normal vector
271  if (dot(n,bar_side) < 0) n = -n;
272 
273  return n;
274 }
275 
276 
277 
278 template<unsigned int dim>
279 auto RefElement<dim>::barycentric_on_face(const BaryPoint &barycentric, unsigned int i_face) -> FaceBaryPoint
280 {
281  ASSERT_EQ(barycentric.n_rows, dim+1);
282  FaceBaryPoint face_barycentric;
283  for(unsigned int i=0; i < dim; i++) {
284 // unsigned int i_sub_node = (i+1)%dim;
285 // unsigned int i_bary = (dim + side_nodes_[i_face][i_sub_node])%(dim+1);
286  unsigned int i_bary = interact_<0,dim-1>(i_face)[i];
287  face_barycentric[i] = barycentric[ i_bary ];
288  }
289  return face_barycentric;
290 }
291 
292 
293 template<unsigned int dim>
294 std::pair<unsigned int, unsigned int> RefElement<dim>::zeros_positions(const BaryPoint &barycentric,
295  double tolerance)
296 {
297  unsigned int zeros = 0;
298  unsigned int n_zeros = 0;
299  for(unsigned int i=0; i < dim+1; i++){
300  if(std::fabs(barycentric[i]) < tolerance)
301  {
302  zeros = zeros | (1 << i);
303  n_zeros++;
304  }
305  }
306 
307  return std::make_pair(n_zeros, zeros);
308 }
309 
310 
311 template<>
312 auto RefElement<0>::clip(const BaryPoint &barycentric) -> BaryPoint
313 {
314  return barycentric;
315 }
316 
317 template<unsigned int dim>
318 auto RefElement<dim>::make_bary_unit_vec()->BarycentricUnitVec
319 {
320  std::vector<arma::vec::fixed<dim+1> > bary_unit_vec(dim+1, arma::zeros(dim+1));
321  for(unsigned int i=0; i<dim; i++) {
322  bary_unit_vec[i][i] = 1.0;
323  bary_unit_vec[i][dim] = -1.0;
324  bary_unit_vec[dim][i] = -1.0 / dim;
325  }
326  bary_unit_vec[dim][dim] = 1.0;
327  return bary_unit_vec;
328 }
329 
330 
331 template<unsigned int dim>
332 auto RefElement<dim>::clip(const BaryPoint &barycentric) -> BaryPoint
333 {
334  static BarycentricUnitVec bary_unit_vec = make_bary_unit_vec();
335  ASSERT_EQ(barycentric.n_rows, dim+1);
336  for(unsigned int i_bary=0; i_bary < dim +1; i_bary ++) {
337  if (barycentric[i_bary] < 0.0) {
338  // index of barycentric coord that is constant on the face i_side
339  // as we use barycentric coords starting with local coordinates:
340  // TODO: rather work only with local coords and/or with canonical barycentric coords
341  unsigned int i_side = (dim - i_bary);
342  // project to face
343  arma::vec projection_to_face(dim+1);
344  //barycentric.print(cout, "input");
345  //cout << "is: " << i_side << endl;
346  //cout << "ibary: " << i_bary << endl;
347  //bary_unit_vec[i_bary].print(cout, "normal");
348  //barycentric.subvec(0, dim-1).print(cout, "bary sub");
349  projection_to_face = barycentric - barycentric[i_bary]*bary_unit_vec[i_bary];
350  //projection_to_face(dim) = 1.0 - arma::sum(projection_to_face.subvec(0, dim-1));
351  //projection_to_face.print(cout, "projection");
352  auto bary_on_face = barycentric_on_face(projection_to_face, i_side);
353  //bary_on_face.print(cout, "b on f");
354  auto sub_clip = RefElement<dim-1>::clip(bary_on_face);
355  //sub_clip.print(cout, "sub clip");
356  return interpolate<dim-1>(sub_clip, i_side);
357  }
358  }
359  return barycentric;
360 
361 }
362 
363 
364 
365 template<unsigned int dim>
366 auto RefElement<dim>::centers_of_subelements(unsigned int sub_dim)->CentersList
367 {
369  if (list.size() == 0) {
370  list.resize(dim+1);
371  for(unsigned int sdim=0; sdim < dim+1; sdim++) {
372  // Temporary solution until we have unified interface to
373  // the internal indexing.
374  // We use the fact that numbering of subelements goes as numbering of
375  // k combinations over nodes.
376 // std::vector<unsigned int> subel_comb(sdim+2);
377  for(auto &sub_el_nodes : nodes_of_subelements[sdim]) {
378  ASSERT_EQ(sub_el_nodes.size(), sdim+1);
379  LocalPoint center = arma::zeros(dim);
380  for( unsigned int i_node : sub_el_nodes)
381  center+=node_coords( i_node );
382  center/=(sdim+1);
383  list[sdim].push_back(center);
384  }
385  }
386  }
387 
388  ASSERT_LE(sub_dim, dim);
389  return list[sub_dim];
390 }
391 
392 
393 template<>
394 double RefElement<1>::side_measure(unsigned int sid)
395 {
396  ASSERT_LT(sid, n_sides).error("Side number is out of range!");
397 
398  return 1;
399 }
400 
401 
402 template<>
403 double RefElement<2>::side_measure(unsigned int sid)
404 {
405  ASSERT_LT(sid, n_sides).error("Side number is out of range!");
406 
407  return norm(node_coords(line_nodes_[sid][1]) - node_coords(line_nodes_[sid][0]),2);
408 }
409 
410 
411 template<>
412 double RefElement<3>::side_measure(unsigned int sid)
413 {
414  ASSERT_LT(sid, n_sides).error("Side number is out of range!");
415 
416  return 0.5*norm(cross(node_coords(side_nodes_[sid][1]) - node_coords(side_nodes_[sid][0]),
417  node_coords(side_nodes_[sid][2]) - node_coords(side_nodes_[sid][0])),2);
418 }
419 
420 
421 
422 template <>
423 unsigned int RefElement<3>::line_between_faces(unsigned int f1, unsigned int f2) {
424  unsigned int i,j;
425  i=j=0;
426  while (side_lines_[f1][i] != side_lines_[f2][j])
427  if (side_lines_[f1][i] < side_lines_[f2][j]) i++;
428  else j++;
429  return side_lines_[f1][i];
430 }
431 
432 
433 /**
434  * Basic line interpolation
435  */
436 template<unsigned int dim>
438  arma::vec::fixed<dim+1> first_coords,
439  arma::vec::fixed<dim+1> second_coords,
440  double first_theta, double second_theta, double theta){
441 
442  arma::vec::fixed<dim+1> bary_interpolated_coords;
443  bary_interpolated_coords = ((theta - first_theta) * second_coords + (second_theta - theta) * first_coords)
444  /(second_theta - first_theta);
445  return bary_interpolated_coords;
446 }
447 
448 template class RefElement<0>;
449 template class RefElement<1>;
450 template class RefElement<2>;
451 template class RefElement<3>;
452 
453 
#define ASSERT_LT(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:301
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
#define ASSERT_LE(a, b)
Definition of comparative assert macro (Less or Equal) only for debug mode.
Definition: asserts.hh:309
static CentersList centers_of_subelements(unsigned int sub_dim)
Definition: ref_element.cc:366
static LocalPoint bary_to_local(const BaryPoint &bp)
Converts from barycentric to local coordinates.
Definition: ref_element.cc:201
static unsigned int normal_orientation(unsigned int sid)
Definition: ref_element.cc:216
Armor::ArmaVec< double, dim > FaceBaryPoint
Definition: ref_element.hh:350
Armor::ArmaVec< double, dim+1 > BaryPoint
Definition: ref_element.hh:349
static FaceBaryPoint barycentric_on_face(const BaryPoint &barycentric, unsigned int i_face)
Definition: ref_element.cc:279
static BaryPoint line_barycentric_interpolation(BaryPoint first_coords, BaryPoint second_coords, double first_theta, double second_theta, double theta)
Definition: ref_element.cc:437
arma::vec::fixed< dim > LocalPoint
Definition: ref_element.hh:342
static BaryPoint local_to_bary(const LocalPoint &lp)
Converts from local to barycentric coordinates.
Definition: ref_element.cc:187
static BarycentricUnitVec make_bary_unit_vec()
Definition: ref_element.cc:318
static std::pair< unsigned int, unsigned int > zeros_positions(const BaryPoint &barycentric, double tolerance=std::numeric_limits< double >::epsilon() *2)
Definition: ref_element.cc:294
static unsigned int line_between_faces(unsigned int f1, unsigned int f2)
static BaryPoint clip(const BaryPoint &barycentric)
Definition: ref_element.cc:332
static double side_measure(unsigned int sid)
static LocalPoint normal_vector(unsigned int sid)
static unsigned int oposite_node(unsigned int sid)
Definition: ref_element.cc:209
static constexpr IdxVector< 3 > S3_node_sides[2][4]
Definition: ref_element.hh:298
static constexpr IdxVector< 3 > S3_node_lines[2][4]
Definition: ref_element.hh:326
static constexpr IdxVector< 2 > S3_line_sides[2][6]
Definition: ref_element.hh:310
Global macros to enhance readability and debugging, general constants.
ArmaVec< double, N > vec
Definition: armor.hh:885
std::vector< std::vector< unsigned int > > _array_to_vec(const std::vector< IdxVector< n >> array_vec)
Definition: ref_element.cc:35
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc.
std::array< unsigned int, Size > IdxVector
Definition: ref_element.hh:276