Flow123d  release_3.0.0-1008-gca43bb7
compute_intersection.cc
Go to the documentation of this file.
1 /*
2  * Author: viktor
3  */
4 
6 #include "mesh/ref_element.hh"
7 #include "mesh/side_impl.hh"
8 #include "mesh/mesh.h"
9 #include "mesh/accessors.hh"
10 #include "system/system.hh"
11 
12 #include "plucker.hh"
14 #include "intersection_aux.hh"
15 
16 using namespace std;
17 
18 /*************************************************************************************************************
19  * COMPUTE INTERSECTION FOR: 1D AND 2D
20  ************************************************************************************************************/
22 : computed_(false)
23 {
24  plucker_coordinates_abscissa_ = nullptr;
25  plucker_coordinates_triangle_.resize(3, nullptr);
26  plucker_products_.resize(3, nullptr);
27 };
28 
29 
31  ElementAccessor<3> triangle, Mesh *mesh)
32 : computed_(false)
33 {
34  ASSERT_DBG(abscissa->dim() == 1);
35  ASSERT_DBG(triangle->dim() == 2);
36  // in this constructor, we suppose this is the final object -> we create all data members
37  plucker_coordinates_abscissa_ = new Plucker(abscissa.node(0),
38  abscissa.node(1), true);
39  scale_line_=plucker_coordinates_abscissa_->scale();
40 
41  plucker_coordinates_triangle_.resize(3);
42  plucker_products_.resize(3);
43  scale_triangle_=std::numeric_limits<double>::max();
44  for(unsigned int side = 0; side < 3; side++){
45  plucker_coordinates_triangle_[side] = new Plucker(triangle.node(RefElement<2>::interact(Interaction<0,1>(side))[0]),
46  triangle.node(RefElement<2>::interact(Interaction<0,1>(side))[1]),
47  true);
48  scale_triangle_ = std::min( scale_triangle_, plucker_coordinates_triangle_[side]->scale());
49 
50  // allocate and compute new Plucker products
51  plucker_products_[side] = new double((*plucker_coordinates_abscissa_)*(*plucker_coordinates_triangle_[side]));
52  }
53 };
54 
56 {
57  if(plucker_coordinates_abscissa_ != nullptr)
58  delete plucker_coordinates_abscissa_;
59 
60  for(unsigned int side = 0; side < RefElement<2>::n_sides; side++){
61  if(plucker_products_[side] != nullptr)
62  delete plucker_products_[side];
63  if(plucker_coordinates_triangle_[side] != nullptr)
64  delete plucker_coordinates_triangle_[side];
65  }
66 }
67 
68 
70  // unset all pointers
71  for(unsigned int side = 0; side < RefElement<2>::n_sides; side++)
72  {
73  plucker_products_[side] = nullptr;
74  plucker_coordinates_triangle_[side] = nullptr;
75  }
76  plucker_coordinates_abscissa_ = nullptr;
77 };
78 
80 
81  // if not already computed, compute plucker coordinates of abscissa
82  plucker_coordinates_abscissa_->compute();
83  scale_line_=plucker_coordinates_abscissa_->scale();
84 
85 // DBGMSG("Abscissa:\n");
86 // (*abscissa_)[0].point_coordinates().print();
87 // (*abscissa_)[1].point_coordinates().print();
88  scale_triangle_=std::numeric_limits<double>::max();
89  // if not already computed, compute plucker coordinates of triangle sides
90  for(unsigned int side = 0; side < RefElement<2>::n_sides; side++){
91  plucker_coordinates_triangle_[side]->compute();
92  scale_triangle_ = std::min( scale_triangle_, plucker_coordinates_triangle_[side]->scale());
93 
94  ASSERT_DBG(plucker_products_[side]).error("Undefined plucker product.");
95  if(*plucker_products_[side] == plucker_empty){
96  *plucker_products_[side] = (*plucker_coordinates_abscissa_)*(*plucker_coordinates_triangle_[side]);
97  }
98 // DBGMSG("Plucker product = %e\n", *(plucker_products_[side]));
99  }
100 };
101 
102 
104 {
105  double tol = geometry_epsilon * scale_line_;
106 
107  double t = RefElement<1>::bary_to_local(IP.local_bcoords_A())(0);
108 
109  int sign = (t < tol ? -2 : (t > 1+tol ? +2 : 0) ); // outside or inside
110  if (std::abs(t) <= tol){ // left end
111  sign = -1;
112  IP.set_topology_A(0,0);
113  }
114  else if (std::abs(1-t) <= tol){ // right end
115  sign = +1;
116  IP.set_topology_A(1,0);
117  }
118 
119  return sign;
120 }
121 
122 
124 {
125  // compute local barycentric coordinates of IP: see formula (3) on pg. 12 in BP VF
126  // local alfa = w2/sum; local beta = w1/sum; => local barycentric coordinates in the triangle
127  // where sum = w0+w1+w2
128 
129  // plucker products with correct signs according to ref_element, ordered by sides
130 // double w[3] = {signed_plucker_product(0),
131 // signed_plucker_product(1),
132 // signed_plucker_product(2)};
133 
134  // local barycentric coordinates of IP, depends on barycentric coordinates order !!!
135 // arma::vec3 local_triangle({w[2],w[1],w[0]});
136 // DBGMSG("Plucker product sum = %f\n",w[0]+w[1]+w[2]);
137 // local_triangle = local_triangle / (w[0]+w[1]+w[2]);
138 // double w_sum = local[0] + local[1] + local[2];
139 // DBGMSG("Plucker product sum = %e %e %e\n",w_sum, 1-rounding_epsilon, 1+rounding_epsilon);
140 
141  //assert inaccurate barycentric coordinates
142  ASSERT_DBG(fabs(1.0 - local[0] - local[1] - local[2]) < geometry_epsilon)(local[0]+local[1]+local[2])
143  (local[0])(local[1])(local[2]);
144 
145 
146 
147  arma::vec3 local_triangle({local[2],local[1],local[0]});
148 
149  // local coordinate T on the line
150  // for i-th coordinate it holds: (from formula (4) on pg. 12 in BP VF)
151  // T = localAbscissa= (- A(i) + ( 1 - alfa - beta ) * V0(i) + alfa * V1(i) + beta * V2 (i)) / U(i)
152  // let's choose [max,i] = max {U(i)}
153  arma::vec3 u = plucker_coordinates_abscissa_->get_u_vector();
154 
155  //find max in u in abs value:
156  unsigned int i = 0; // index of maximum in u
157  double max = u[0]; // maximum in u
158  if(fabs((double)u[1]) > fabs(max)){ max = u[1]; i = 1;}
159  if(fabs((double)u[2]) > fabs(max)){ max = u[2]; i = 2;}
160 
161  // global coordinates in triangle
162  double isect_coord_i =
163  local_triangle[0] * plucker_coordinates_triangle_[0]->point(RefElement<2>::normal_orientation(0))[i] +
164  local_triangle[1] * plucker_coordinates_triangle_[0]->point(1-RefElement<2>::normal_orientation(0))[i] +
165  local_triangle[2] * plucker_coordinates_triangle_[1]->point(RefElement<2>::normal_orientation(1))[i];
166 
167  //theta on abscissa
168  double t = (-plucker_coordinates_abscissa_->point(0)[i] + isect_coord_i)/max;
169  //DebugOut() << print_var(t) << print_var(isect_coord_i) << print_var(max);
170  arma::vec2 local_abscissa = {1-t, t};
171 
172  /*
173  DBGMSG("Coordinates: line; local and global triangle\n");
174  theta.print();
175  local_triangle.print();
176  global_triangle.print();
177  */
178 
179  IP.set_topology_A(0, 1);
180  IP.set_topology_B(0, 2);
181  IP.set_coordinates(local_abscissa,local_triangle);
182 
183  return true;
184 }
185 
186 
188 {
189  compute_plucker_products();
190  computed_ = true;
191 
192  //DebugOut() << "LINE: \n" << (*abscissa_)[0].point_coordinates()
193  // << (*abscissa_)[1].point_coordinates();
194 
195 
196  // convert plucker products to local coords
197  arma::vec3 w = {signed_plucker_product(0),
198  signed_plucker_product(1),
199  signed_plucker_product(2)};
200  double w_sum = w[0] + w[1] + w[2];
201 
202  unsigned int n_positive = 0;
203  unsigned int n_negative = 0;
204  unsigned int zero_idx_sum =0;
205  //DebugOut() << print_var(std::fabs(w_sum));
206  //DebugOut() << print_var(scale_line_);
207  //DebugOut() << print_var(scale_triangle_);
208 
209  double scaled_epsilon = geometry_epsilon*scale_line_*scale_triangle_*scale_triangle_;
210  if(std::fabs(w_sum) > scaled_epsilon) {
211  w = w / w_sum;
212  for (unsigned int i=0; i < 3; i++) {
213  //DebugOut() << print_var(i) << print_var(w[i]);
214  if (w[i] > geometry_epsilon) n_positive++;
215  else if ( w[i] > -geometry_epsilon) zero_idx_sum+=i;
216  else n_negative++;
217  }
218 
219  } else {
220  /* case 'w_sum == 0':
221  * 1] all products are zero => n_negative=0 and n_positive=0 => it is degenerate case (coplanar case)
222  * 2] at least two products are nonzero AND some of them must be negative => no intersection
223  * (it happens when line is parallel with the triangle but not coplanar; unit test line_triangle09.msh)
224  * See the IF conditions below.
225  */
226 
227  for (unsigned int i=0; i < 3; i++) {
228  //DebugOut().fmt("i: {} w[i]: {:g}", i, w[i]);
229  if (w[i] > scaled_epsilon || w[i] < -scaled_epsilon) n_negative++;
230  }
231  // n_positive == 0
232  //DebugOut() << print_var(n_negative);
233  }
234 
235 // w.print("w");
236 
237  // any negative barycentric coordinate means, no intersection
238  if (n_negative>0) return IntersectionResult::none;
239 
240  // test whether any plucker products is non-zero
241  if (n_positive > 0) {
242 
243  compute_plucker(IP, w);
244  // edge of triangle
245  unsigned int non_zero_idx=0;
246  if (n_positive == 2) {
247  // one zero product, intersection on the zero edge
248  // the zero edge index is equal to zero_idx_sum
249  IP.set_topology_B(zero_idx_sum, 1);
250  non_zero_idx = (zero_idx_sum + 1) % 3;
251  }
252  else if (n_positive == 1) {
253  // two zero products, intersection in vertex oposite to single non-zero edge
254  // index of the non-zero edge is 3-zero_idx_sum
255  IP.set_topology_B(RefElement<2>::oposite_node(3-zero_idx_sum), 0);
256  non_zero_idx = 3-zero_idx_sum;
257  }
258  //DebugOut() << print_var(non_zero_idx) << print_var(signed_plucker_product(non_zero_idx));
259 
260  IntersectionResult result = signed_plucker_product(non_zero_idx) > 0 ?
262  IP.set_result(result);
263 
264  return result;
265 
266  } else {
267  ASSERT_DBG(IP.topology_equal(IPAux12())); // check empty IP (none)
268  IP.set_result(IntersectionResult::degenerate); // set deg. result
270  }
271 };
272 
274 {
275  IPAux12 IP;
276  IntersectionResult result = compute(IP);
277 // DBGVAR((int)result);
278  // skip empty cases
279  if(result == IntersectionResult::none) return 0;
280 
281  // standard case with a single intersection corner
282  if(result < IntersectionResult::degenerate){
283 // DBGCOUT(<< "12d plucker case\n");
284  int sign = check_abscissa_topology(IP);
285 // DBGVAR(sign);
286  if(std::abs(sign) > 1) return 0;
287 
288  IP12s.push_back(IP);
289  return IP12s.size();
290  }
291  else{
293 // DBGCOUT(<< "12d degenerate case, all products are zero\n");
294  return compute_final_in_plane(IP12s);
295  }
296 }
297 
298 
300  IPAux12& IP)
301 {
302 // DBGMSG("PluckerProduct[%d]: %f\n",side, *plucker_products_[side]);
303 
304  // We solve following equation for parameters s,t:
305  /* A + sU = C + tV = intersection point
306  * sU - tV = C - A
307  *
308  * which is by components:
309  * (u1 -v1) (s) = (c1-a1)
310  * (u2 -v2) (t) = (c2-a2)
311  * (u3 -v3) = (c3-a3)
312  *
313  * these are 3 equations for variables s,t
314  * see (4.3) on pg. 19 in DP VF
315  *
316  * We will solve this using Crammer's rule for the maximal subdeterminant det_ij of the matrix.
317  * s = detX_ij / det_ij
318  * t = detY_ij / det_ij
319  */
320 
321  // starting point of abscissa
322  arma::vec3 A = plucker_coordinates_abscissa_->point(0);
323  // direction vector of abscissa
324  arma::vec3 U = plucker_coordinates_abscissa_->get_u_vector();
325  // vertex of triangle side
326  arma::vec3 C = plucker_coordinates_triangle_[side]->point(RefElement<2>::normal_orientation(side));
327  // direction vector of triangle side
328  arma::vec3 V = plucker_coordinates_triangle_[side]->get_u_vector();
329  // right hand side
330  arma::vec3 K = C - A;
331  // subdeterminants det_ij of the system equal minus normal vector to common plane of U and V
332  // det_12 = (-UxV)[1]
333  // det_13 = -(-UxV)[2]
334  // det_23 = (-UxV)[3]
335  arma::vec3 Det = -arma::cross(U,V);
336 
337 
338 // A.print("A");
339 // U.print("U");
340 // C.print("C");
341 // V.print("V");
342 // K.print("K");
343 // Det.print();
344 // cout <<endl;
345 
346  unsigned int max_index = 0;
347  double maximum = fabs(Det[0]);
348  if(fabs((double)Det[1]) > maximum){
349  maximum = fabs(Det[1]);
350  max_index = 1;
351  }
352  if(fabs((double)Det[2]) > maximum){
353  maximum = fabs(Det[2]);
354  max_index = 2;
355  }
356 // DBGVAR(maximum);
357  //abscissa is parallel to triangle side
358  //TODO What tolerance should we use here and why?
359 // if(std::abs(maximum) <= std::sqrt(geometry_epsilon)) return false;
360  if(std::abs(maximum) <= geometry_epsilon) return false;
361 
362  // map maximum index in {-UxV} to i,j of subdeterminants
363  // i j
364  // max_index 0: 1 2
365  // 1: 2 0 (switch due to sign change)
366  // 2: 0 1
367  unsigned int i = (max_index+1)%3,
368  j = (max_index+2)%3;
369 
370  double DetX = -K[i]*V[j] + K[j]*V[i];
371  double DetY = -K[i]*U[j] + K[j]*U[i];
372 
373  double s = DetX/Det[max_index]; //parameter on abscissa
374  double t = DetY/Det[max_index]; //parameter on triangle side
375 
376  // change sign according to side orientation
377  if(RefElement<2>::normal_orientation(side)) t=-t;
378 
379 // DBGVAR(s);
380 // DBGVAR(t);
381 
382  //TODO correct tolerance with scale; compute scale without plucker coords
383  double tol = geometry_epsilon;
384  // if IP is inside of triangle side
385  if(t >= -tol && t <= 1+tol){
386 
387  IP.set_result(IntersectionResult::degenerate); // set orientation as a pathologic case ( > 1)
388 
389  // possibly set triangle vertex {0,1,2}
391  else if(fabs(1-t) <= tol) { IP.set_topology_B(RefElement<2>::interact(Interaction<0,1>(side))[1-RefElement<2>::normal_orientation(side)],0);}
392  else { IP.set_topology_B(side,1);} // no vertex, side side, dim = 1
393 
394  arma::vec2 local_abscissa({1-s, s});
395  arma::vec3 local_triangle({0,0,0});
396 
397  // set local triangle barycentric coordinates according to nodes of the triangle side:
398  local_triangle[RefElement<2>::interact(Interaction<0,1>(side))[RefElement<2>::normal_orientation(side)]] = 1 - t;
400 // local_abscissa.print();
401 // local_triangle.print();
402 
403  IP.set_coordinates(local_abscissa,local_triangle);
404  return true; // IP found
405  }
406 
407  return false; // IP NOT found
408 }
409 
410 
412 {
413  std::vector<int> sign;
414 
415  for(unsigned int i = 0; i < 3;i++){
416 // DBGCOUT( << "side [" << i << "]\n");
417  if(IP12s.size() == 2) break;
418 
419  IPAux12 IP;
420  if (compute_degenerate(i,IP))
421  {
422  uint s = check_abscissa_topology(IP);
423 // cout << IP;
424  // check whether we have found the same IP (e.g. vertex case)
425  if(IP12s.size() > 0 && IP.topology_equal(IP12s.back())) continue;
426 
427  sign.push_back(s);
428  IP12s.push_back(IP);
429  }
430  }
431  if(IP12s.size() == 0) return IP12s.size();
432 
433  // in the case, that line goes through vertex, but outside tetrahedron (touching vertex)
434  if(IP12s.size() == 1){
435  if(std::abs(sign[0]) > 1) IP12s.pop_back(); // outside abscissa
436  return IP12s.size();
437  }
438 
439  // 2 IPs CASE:
440  ASSERT_EQ_DBG(2, IP12s.size());
441 
442 // DBGVAR(sign[0]);
443 // DebugOut() << IP12s[0];
444 // DBGVAR(sign[1]);
445 // DebugOut() << IP12s[1];
446 
447  // intersection outside of abscissa => NO intersection
448  if (sign[0] == sign[1] && std::abs(sign[0]) > 1) {
449  IP12s.clear();
450  return 0;
451  }
452 
453  // order IPs according to the abscissa parameter
454  if(IP12s[0].local_bcoords_A()[1] > IP12s[1].local_bcoords_A()[1]){
455  std::swap(IP12s[0], IP12s[1]);
456  std::swap(sign[0], sign[1]);
457  }
458 
459  // possibly cut IPs to abscissa ends and interpolate tetrahedron bcoords
460  const int ip_sign[] = {-2, +2}; // states to cut
461  for( unsigned int ip=0; ip<2; ip++) {
462  // cut every IP to its end of the abscissa
463  if (sign[ip] == ip_sign[ip]) {
464  sign[ip] /=2; // -2 to -1; +2 to +1
465  correct_triangle_ip_topology(double(ip), ip, IP12s);
466  IP12s[ip].set_topology_A(ip, 0);
467  }
468  }
469 
470  // if IPs are the same, then throw the second one away
471  if(IP12s[0].topology_equal(IP12s[1])){
472 // DebugOut() << "IP same, pop up\n";
473  IP12s.pop_back();
474  }
475 
476  return IP12s.size();
477 }
478 
480  double t, unsigned int ip, std::vector<IPAux12> &ips)
481 {
483  ips[0].local_bcoords_B(), ips[1].local_bcoords_B(),
484  ips[0].local_bcoords_A()[1], ips[1].local_bcoords_A()[1], t);
485  arma::vec2 local_abscissa({1 - t, t}); // abscissa local barycentric coords
486  ips[ip].set_coordinates(local_abscissa, local_tria);
487 
488  // create mask for zeros in barycentric coordinates
489  // coords (*,*,*,*) -> byte bitwise xxxx
490  // only least significant one byte used from the integer
491  std::pair<unsigned int, unsigned int> zeros =
492  RefElement<2>::zeros_positions(ips[ip].local_bcoords_B(), geometry_epsilon);
493 
494  /**
495  * TODO:
496  * 1. Try to avoid setting topology from coords. Try to use just topology information.
497  * 2. If can not be done. Use interact method to setup a map mapping 16 possible zeros positions to appropriate topology,
498  * remove topology_idx from RefElement.
499  */
500 
501 // DebugOut() << "zeros: <" << zeros.first << ", " << zeros.second << ">\n";
502  switch(zeros.first)
503  {
504  default: ips[ip].set_topology_B(0,2); //inside tetrahedon
505  break;
506  case 1: ips[ip].set_topology_B(RefElement<2>::topology_idx<1>(zeros.second),1);
507  break;
508  case 2: ips[ip].set_topology_B(RefElement<2>::topology_idx<0>(zeros.second),0);
509  break;
510  }
511 };
512 
513 
515  os << "\tPluckerCoordinates Abscissa[0]";
516  if(plucker_coordinates_abscissa_ == nullptr){
517  os << "NULL" << endl;
518  }else{
519  os << plucker_coordinates_abscissa_;
520  }
521  for(unsigned int i = 0; i < 3;i++){
522  os << "\tPluckerCoordinates Triangle[" << i << "]";
523  if(plucker_coordinates_triangle_[i] == nullptr){
524  os << "NULL" << endl;
525  }else{
526  os << plucker_coordinates_triangle_[i];
527  }
528  }
529 };
530 
531 
532 
533 
534 /*************************************************************************************************************
535  * COMPUTE INTERSECTION FOR: 2D AND 2D
536  ************************************************************************************************************/
538 {
539  plucker_coordinates_.resize(2*RefElement<2>::n_sides, nullptr);
540  plucker_products_.resize(3*RefElement<2>::n_sides, nullptr);
541 };
542 
544  ElementAccessor<3> triaB,
545  Mesh *mesh)
546 {
547  ASSERT_DBG(triaA->dim() == 2);
548  ASSERT_DBG(triaB->dim() == 2);
549  plucker_coordinates_.resize(2*RefElement<2>::n_sides);
550  plucker_products_.resize(3*RefElement<2>::n_sides);
551 
552  for(unsigned int side = 0; side < RefElement<2>::n_sides; side++){
553  plucker_coordinates_[side] = new Plucker(triaA.node(RefElement<2>::interact(Interaction<0,1>(side))[0]),
555  plucker_coordinates_[RefElement<2>::n_sides+side]
556  = new Plucker(triaB.node(RefElement<2>::interact(Interaction<0,1>(side))[0]),
558  }
559 
560  // compute Plucker products for each pair triangle A side and triangle B side
561  for(unsigned int p = 0; p < 3*RefElement<2>::n_sides; p++){
562  plucker_products_[p] = new double(plucker_empty);
563  }
564 };
565 
567 {
568  // unset pointers:
569  for(unsigned int side = 0; side < 2*RefElement<2>::n_sides; side++)
570  CI12[side].clear_all();
571 
572  // then delete objects:
573  for(unsigned int side = 0; side < 2*RefElement<2>::n_sides; side++){
574  if(plucker_coordinates_[side] != nullptr)
575  delete plucker_coordinates_[side];
576  }
577 
578  for(unsigned int p = 0; p < 3*RefElement<2>::n_sides; p++){
579  if(plucker_products_[p] != nullptr)
580  delete plucker_products_[p];
581  }
582 }
583 
585 {
586  // unset all pointers
587  for(unsigned int side = 0; side < 2*RefElement<2>::n_sides; side++)
588  {
589  plucker_coordinates_[side] = nullptr;
590  }
591  for(unsigned int p = 0; p < 3*RefElement<2>::n_sides; p++){
592  plucker_products_[p] = nullptr;
593  }
594 }
595 
597 
598 // DBGMSG("init\n");
599  for(unsigned int i = 0; i < RefElement<2>::n_sides; i++){
600  // set side A vs triangle B
601  for(unsigned int j = 0; j < RefElement<2>::n_sides; j++)
602  CI12[i].set_pc_triangle(plucker_coordinates_[3+j], j); // set triangle B
603  // set side of triangle A
604  CI12[i].set_pc_abscissa(plucker_coordinates_[i]);
605 
606  // set side B vs triangle A
607  for(unsigned int j = 0; j < RefElement<2>::n_sides; j++)
608  CI12[RefElement<2>::n_sides + i].set_pc_triangle(plucker_coordinates_[j], j); // set triangle A
609  // set side of triangle B
610  CI12[RefElement<2>::n_sides + i].set_pc_abscissa(plucker_coordinates_[RefElement<2>::n_sides + i]);
611 
612  // set plucker products
613  for(unsigned int j = 0; j < RefElement<2>::n_sides; j++)
614  {
615  //for A[i]_B set pp. A[i] x B[j]
616  CI12[i].set_plucker_product(get_plucker_product(i,j),j);
617  //for B[i]_A set pp. A[j] x B[i]
618  CI12[RefElement<2>::n_sides + i].set_plucker_product(get_plucker_product(j,i),j);
619  }
620 
621  }
622 };
623 
624 
626 {
627  // final intersection points
628  std::vector<IPAux22> &IP22s = intersection.points();
629  // temporary vector for lower dimensional IPs
630  std::vector<IPAux12> IP12s;
631  IP12s.reserve(2);
632  unsigned int ip_coutner = 0;
633 
634  // loop over CIs (side vs triangle): [A0_B, A1_B, A2_B, B0_A, B1_A, B2_A].
635  for(unsigned int i = 0; i < 2*RefElement<2>::n_sides && ip_coutner < 2; i++){
636  if(!CI12[i].is_computed()) // if not computed yet
637  {
638 // DBGVAR(i);
639  if(CI12[i].compute_final(IP12s) == 0) continue;
640 
641  unsigned int triangle_side = i%3; //i goes from 0 to 5 -> i%3 = 0,1,2,0,1,2
642 
643  for(IPAux12 &IP : IP12s)
644  {
645  IPAux22 IP22(IP.switch_objects(), triangle_side); // swicth dim 12 -> 21; interpolate dim 21 -> 22
646 
647  if(i < 3){
648  //switch back to keep order of triangles [A,B]
649  IP22 = IP22.switch_objects();
650 
651  if( IP22.dim_A() == 0 ) // if IP is vertex of triangle
652  {
653 // DBGCOUT("set_node A\n");
654  // we are on line of the triangle A, and IP.idx_A contains local node of the line
655  // we know vertex index
656 
657  // set flag on all sides of tetrahedron connected by the node
658  for(unsigned int s=0; s < RefElement<2>::n_sides_per_node; s++)
659  CI12[RefElement<2>::interact(Interaction<1,0>(IP22.idx_A()))[s]].set_computed();
660  }
661  if( IP22.dim_B() == 0 ) // if IP is vertex of triangle
662  {
663 // DBGCOUT("set_node B\n");
664  // set flag on both sides of triangle connected by the node
665  for(unsigned int s=0; s < RefElement<2>::n_sides_per_node; s++)
666  CI12[3 + RefElement<2>::interact(Interaction<1,0>(IP22.idx_B()))[s]].set_computed();
667  }
668  else if( IP22.dim_B() == 1 ) // if IP is vertex of triangle
669  {
670 // DBGCOUT("set line B\n");
671  // set flag on both sides of triangle connected by the node
672  CI12[3 + IP22.idx_B()].set_computed();
673  }
674  }
675  else if( IP22.dim_B() == 0 ) // if IP is vertex of triangle B (triangles switched!!! A <-> B)
676  {
677  //we do not need to look back to triangle A, because if IP was at vertex, we would know already
678 // DBGCOUT("set_node B\n");
679  // we are on line of the triangle, and IP.idx_A contains local node of the line
680  // we know vertex index
681 
682  // set flag on both sides of triangle connected by the node
683  for(unsigned int s=0; s < RefElement<2>::n_sides_per_node; s++)
684  CI12[3 + RefElement<2>::interact(Interaction<1,0>(IP22.idx_B()))[s]].set_computed();
685  }
686 // DBGCOUT( << IP22);
687  ip_coutner++;
688  IP22s.push_back(IP22);
689  }
690  IP12s.clear();
691  }
692  }
693 
694  return ip_coutner;
695 };
696 
697 // void ComputeIntersection< 2, 2>::correct_triangle_ip_topology(IntersectionPointAux<2,2>& ip)
698 // {
699 // // create mask for zeros in barycentric coordinates
700 // // coords (*,*,*,*) -> byte bitwise xxxx
701 // // only least significant one byte used from the integer
702 // unsigned int zeros = 0;
703 // unsigned int n_zeros = 0;
704 // for(char i=0; i < 3; i++){
705 // if(std::fabs(ip.local_bcoords_B()[i]) < geometry_epsilon)
706 // {
707 // zeros = zeros | (1 << (2-i));
708 // n_zeros++;
709 // }
710 // }
711 //
712 // switch(n_zeros)
713 // {
714 // default: ip.set_topology_B(0,2); //inside triangle
715 // break;
716 // case 1: ip.set_topology_B(RefElement<2>::topology_idx<1>(zeros),2);
717 // break;
718 // case 2: ip.set_topology_B(RefElement<2>::topology_idx<0>(zeros),1);
719 // break;
720 // }
721 // };
722 
723 
725 
726  for(unsigned int i = 0; i < RefElement<2>::n_lines; i++){
727  os << "\tPluckerCoordinates Triangle A[" << i << "]";
728  if(plucker_coordinates_[i] == nullptr){
729  os << "NULL" << endl;
730  }else{
731  os << plucker_coordinates_[i];
732  }
733  }
734  for(unsigned int i = 0; i < RefElement<2>::n_lines; i++){
735  os << "\tPluckerCoordinates Triangle B[" << i << "]";
736  if(plucker_coordinates_[RefElement<2>::n_lines + i] == nullptr){
737  os << "NULL" << endl;
738  }else{
739  os << plucker_coordinates_[RefElement<2>::n_lines + i];
740  }
741  }
742 };
743 
745  os << "ComputeIntersection<2,2> Plucker Coordinates Tree:" << endl;
746  print_plucker_coordinates(os);
747  for(unsigned int i = 0; i < 6;i++){
748  os << "ComputeIntersection<1,2>["<< i <<"] Plucker Coordinates:" << endl;
749  CI12[i].print_plucker_coordinates(os);
750  }
751 };
752 
753 
754 /*************************************************************************************************************
755  * COMPUTE INTERSECTION FOR: 1D AND 3D
756  ************************************************************************************************************/
758 {
759  plucker_coordinates_abscissa_ = nullptr;
760  plucker_coordinates_tetrahedron.resize(6, nullptr);
761  plucker_products_.resize(6, nullptr);
762 };
763 
765  ElementAccessor<3> tetrahedron,
766  Mesh *mesh)
767 {
768  ASSERT_DBG(abscissa->dim() == 1);
769  ASSERT_DBG(tetrahedron->dim() == 3);
770 
771  plucker_coordinates_abscissa_ = new Plucker(abscissa.node(0), abscissa.node(1));
772  plucker_coordinates_tetrahedron.resize(6);
773  plucker_products_.resize(6);
774 
775  for(unsigned int line = 0; line < RefElement<3>::n_lines; line++){
776  plucker_coordinates_tetrahedron[line] = new Plucker(tetrahedron.node(RefElement<3>::interact(Interaction<0,1>(line))[0]),
777  tetrahedron.node(RefElement<3>::interact(Interaction<0,1>(line))[1]));
778  // compute Plucker products (abscissa X tetrahedron line)
779  plucker_products_[line] = new double(plucker_empty);
780  }
781 };
782 
784 {
785  // unset pointers:
786  for(unsigned int side = 0; side < RefElement<3>::n_sides; side++)
787  CI12[side].clear_all();
788 
789  // then delete objects:
790  if(plucker_coordinates_abscissa_ != nullptr)
791  delete plucker_coordinates_abscissa_;
792 
793  for(unsigned int line = 0; line < RefElement<3>::n_lines; line++){
794  if(plucker_products_[line] != nullptr)
795  delete plucker_products_[line];
796  if(plucker_coordinates_tetrahedron[line] != nullptr)
797  delete plucker_coordinates_tetrahedron[line];
798  }
799 }
800 
802 {
803  // unset all pointers
804  for(unsigned int side = 0; side < RefElement<3>::n_lines; side++)
805  {
806  plucker_products_[side] = nullptr;
807  plucker_coordinates_tetrahedron[side] = nullptr;
808  }
809  plucker_coordinates_abscissa_ = nullptr;
810 }
811 
813 
814  for(unsigned int side = 0; side < RefElement<3>::n_sides; side++){
815  for(unsigned int line = 0; line < RefElement<3>::n_lines_per_side; line++){
816  CI12[side].set_pc_triangle(plucker_coordinates_tetrahedron[
817  RefElement<3>::interact(Interaction<1,2>(side))[line]], line);
818 
819  CI12[side].set_plucker_product(
820  plucker_products_[RefElement<3>::interact(Interaction<1,2>(side))[line]],
821  line);
822  }
823  CI12[side].set_pc_abscissa(plucker_coordinates_abscissa_);
824  }
825 };
826 
827 
829 {
830  return compute(intersection.i_points_);
831 }
832 
834  ASSERT_EQ_DBG(0, IP13s.size());
835  std::vector<int> sign;
836 
837  // loop over faces of tetrahedron
838  for(unsigned int face = 0; face < RefElement<3>::n_sides && IP13s.size() < 2; face++){
839 
840  if (CI12[face].is_computed()) continue;
841 
843  IntersectionResult result = CI12[face].compute(IP12);
844  //DebugOut() << print_var(face) << print_var(int(result)) << "1d-3d";
845 
846  if (int(result) < int(IntersectionResult::degenerate) ) {
847  // fix topology of A (abscissa) in IP12 and save sign
848  sign.push_back(CI12[face].check_abscissa_topology(IP12));
849 
850  IPAux IP13(IP12, face);
851  //DebugOut() << IP13;
852  IP13s.push_back(IP13);
853 
854  // set the 'computed' flag on the connected sides by IP
855  if(IP13.dim_B() == 0) // IP is vertex of triangle
856  {
857  // set flag on all sides of tetrahedron connected by the node
858  for(unsigned int node_face : RefElement<3>::interact(Interaction<2,0>(IP13.idx_B())))
859  CI12[node_face].set_computed();
860  }
861  else if(IP13.dim_B() == 1) // IP is on edge of triangle
862  {
863  for(unsigned int edge_face : RefElement<3>::interact(Interaction<2,1>(IP13.idx_B())))
864  CI12[edge_face].set_computed();
865  }
866  }
867  }
868  if (IP13s.size() == 0) return IP13s.size();
869 
870  // in the case, that line goes through vertex, but outside tetrahedron (touching vertex)
871  if(IP13s.size() == 1){
872  if(std::abs(sign[0]) > 1) IP13s.pop_back(); // outside abscissa
873  return IP13s.size();
874  }
875 
876  // 2 IPs CASE:
877  ASSERT_EQ_DBG(2, IP13s.size());
878 
879  // intersection outside of abscissa => NO intersection
880  if (sign[0] == sign[1] && std::abs(sign[0]) > 1) {
881  IP13s.clear();
882  return 0;
883  }
884 
885  // order IPs according to the abscissa parameter
886  if(IP13s[0].local_bcoords_A()[1] > IP13s[1].local_bcoords_A()[1]){
887  std::swap(IP13s[0], IP13s[1]);
888  std::swap(sign[0], sign[1]);
889  }
890 
891  // possibly cut IPs to abscissa ends and interpolate tetrahedron bcoords
892  const int ip_sign[] = {-2, +2}; // states to cut
893  for( unsigned int ip=0; ip<2; ip++) {
894  // cut every IP to its end of the abscissa
895  if (sign[ip] == ip_sign[ip]) {
896  sign[ip] /=2; // -2 to -1; +2 to +1
897  correct_tetrahedron_ip_topology(double(ip), ip, IP13s);
898  IP13s[ip].set_topology_A(ip, 0);
899  }
900  }
901 
902  // if IPs are the same, then throw the second one away
903  if(IP13s[0].topology_equal(IP13s[1])){
904  IP13s.pop_back();
905  }
906 
907  return IP13s.size();
908 };
909 
911  double t, unsigned int ip, std::vector<IPAux> &ips)
912 {
913  arma::vec4 local_tetra = RefElement<3>::line_barycentric_interpolation(
914  ips[0].local_bcoords_B(), ips[1].local_bcoords_B(),
915  ips[0].local_bcoords_A()[1], ips[1].local_bcoords_A()[1], t);
916  arma::vec2 local_abscissa({1 - t, t}); // abscissa local barycentric coords
917  ips[ip].set_coordinates(local_abscissa, local_tetra);
918 
919  // create mask for zeros in barycentric coordinates
920  // coords (*,*,*,*) -> byte bitwise xxxx
921  // only least significant one byte used from the integer
922  unsigned int zeros = 0;
923  unsigned int n_zeros = 0;
924  for(char i=0; i < 4; i++){
925  if(std::fabs(ips[ip].local_bcoords_B()[i]) < geometry_epsilon)
926  {
927  zeros = zeros | (1 << i);
928  n_zeros++;
929  }
930  }
931 
932  /**
933  * TODO:
934  * 1. Try to avoid setting topology from coords. Try to use just topology information.
935  * 2. If can not be done. Use interact method to setup a map mapping 16 possible zeros positions to appropriate topology,
936  * remove topology_idx from RefElement.
937  */
938 
939  switch(n_zeros)
940  {
941  default: ips[ip].set_topology_B(0,3); //inside tetrahedon
942  break;
943  case 1: ips[ip].set_topology_B(RefElement<3>::topology_idx<2>(zeros),2);
944  break;
945  case 2: ips[ip].set_topology_B(RefElement<3>::topology_idx<1>(zeros),1);
946  break;
947  case 3: ips[ip].set_topology_B(RefElement<3>::topology_idx<0>(zeros),0);
948  break;
949  }
950 };
951 
952 
954  os << "\tPluckerCoordinates Abscissa[0]";
955  if(plucker_coordinates_abscissa_ == nullptr){
956  os << "NULL" << endl;
957  }else{
958  os << plucker_coordinates_abscissa_;
959  }
960 
961  for(unsigned int i = 0; i < 6;i++){
962  os << "\tPluckerCoordinates Tetrahedron[" << i << "]";
963  if(plucker_coordinates_tetrahedron[i] == nullptr){
964  os << "NULL" << endl;
965  }else{
966  os << plucker_coordinates_tetrahedron[i];
967  }
968  }
969 };
970 
972  os << "ComputeIntersection<1,3> Plucker Coordinates Tree:" << endl;
973  print_plucker_coordinates(os);
974  for(unsigned int i = 0; i < 4;i++){
975  os << "ComputeIntersection<1,2>["<< i <<"] Plucker Coordinates:" << endl;
976  CI12[i].print_plucker_coordinates(os);
977  }
978 };
979 
980 
981 
982 /*************************************************************************************************************
983  * COMPUTE INTERSECTION FOR: 2D AND 3D
984  ************************************************************************************************************/
986 : no_idx(100),
987 s3_dim_starts({0, 4, 10, 14}), // vertices, edges, faces, volume
988 s2_dim_starts({15, 18, 21}), // vertices, sides, surface
989 object_next(22, no_idx), // 4 vertices, 6 edges, 4 faces, 1 volume, 3 corners, 3 sides, 1 surface; total 22
990 on_faces(_on_faces())
991  {
992 
993  plucker_coordinates_triangle_.resize(3, nullptr);
994  plucker_coordinates_tetrahedron.resize(6, nullptr);
995  plucker_products_.resize(3*6, nullptr);
996 };
997 
998 
1000  ElementAccessor<3> tetrahedron,
1001  Mesh *mesh)
1003 {
1004  mesh_ = mesh;
1005  plucker_coordinates_triangle_.resize(3);
1006  plucker_coordinates_tetrahedron.resize(6);
1007 
1008  // set CI object for 1D-2D intersection 'tetrahedron edge - triangle'
1009  for(unsigned int i = 0; i < RefElement<3>::n_lines; i++){
1010  plucker_coordinates_tetrahedron[i] = new Plucker(tetrahedron.node(RefElement<3>::interact(Interaction<0,1>(i))[0]),
1011  tetrahedron.node(RefElement<3>::interact(Interaction<0,1>(i))[1]));
1012  }
1013  // set CI object for 1D-3D intersection 'triangle side - tetrahedron'
1014  for(unsigned int i = 0; i < RefElement<2>::n_lines;i++){
1015  plucker_coordinates_triangle_[i] = new Plucker(triangle.node(RefElement<2>::interact(Interaction<0,1>(i))[0]),
1016  triangle.node(RefElement<2>::interact(Interaction<0,1>(i))[1]));
1017  }
1018 
1019  // compute Plucker products (triangle side X tetrahedron line)
1020  // order: triangle sides X tetrahedron lines:
1021  // TS[0] X TL[0..6]; TS[1] X TL[0..6]; TS[1] X TL[0..6]
1022  unsigned int np = RefElement<2>::n_sides * RefElement<3>::n_lines;
1023  plucker_products_.resize(np, nullptr);
1024  for(unsigned int line = 0; line < np; line++){
1025  plucker_products_[line] = new double(plucker_empty);
1026 
1027  }
1028 };
1029 
1031 {
1032  // unset pointers:
1033  for(unsigned int triangle_side = 0; triangle_side < RefElement<2>::n_sides; triangle_side++)
1034  CI13[triangle_side].clear_all();
1035 
1036  for(unsigned int line = 0; line < RefElement<3>::n_lines; line++)
1037  CI12[line].clear_all();
1038 
1039  // then delete objects:
1040  unsigned int np = RefElement<2>::n_sides * RefElement<3>::n_lines;
1041  for(unsigned int line = 0; line < np; line++){
1042  if(plucker_products_[line] != nullptr)
1043  delete plucker_products_[line];
1044  }
1045 
1046  for(unsigned int i = 0; i < RefElement<3>::n_lines;i++){
1047  if(plucker_coordinates_tetrahedron[i] != nullptr)
1048  delete plucker_coordinates_tetrahedron[i];
1049  }
1050  for(unsigned int i = 0; i < RefElement<2>::n_sides;i++){
1051  if(plucker_coordinates_triangle_[i] != nullptr)
1052  delete plucker_coordinates_triangle_[i];
1053  }
1054 };
1055 
1057 
1058  // set pointers to Plucker coordinates for 1D-2D
1059  // set pointers to Plucker coordinates for 1D-3D
1060  // distribute Plucker products - CI13
1061  for(unsigned int triangle_side = 0; triangle_side < RefElement<2>::n_sides; triangle_side++){
1062  for(unsigned int line = 0; line < RefElement<3>::n_lines; line++){
1063  CI13[triangle_side].set_plucker_product(
1064  plucker_products_[triangle_side * RefElement<3>::n_lines + line],
1065  line);
1066  CI12[line].set_plucker_product(
1067  plucker_products_[triangle_side * RefElement<3>::n_lines + line],
1068  triangle_side);
1069 
1070  CI13[triangle_side].set_pc_tetrahedron(plucker_coordinates_tetrahedron[line],line);
1071  CI12[line].set_pc_triangle(plucker_coordinates_triangle_[triangle_side],triangle_side);
1072  }
1073  CI13[triangle_side].set_pc_abscissa(plucker_coordinates_triangle_[triangle_side]);
1074  CI13[triangle_side].init();
1075  }
1076 
1077  // set pointers to Plucker coordinates for 1D-2D
1078  for(unsigned int line = 0; line < RefElement<3>::n_lines; line++)
1079  CI12[line].set_pc_abscissa(plucker_coordinates_tetrahedron[line]);
1080 
1081 };
1082 
1083 
1084 
1086  ASSERT_LT_DBG(i_obj, object_next.size());
1087  unsigned int ip = object_next[i_obj];
1088  if (ip == no_idx) return false;
1089  ASSERT_LT_DBG(ip, IP_next.size());
1090  return IP_next[ip] == i_obj;
1091 }
1092 
1093 /**
1094  * Set links: obj_before -> IP -> obj_after
1095  * if obj_after have null successor, set obj_after -> IP (backlink)
1096  */
1097 void ComputeIntersection<2,3>::set_links(uint obj_before_ip, uint ip_idx, uint obj_after_ip) {
1098  if (have_backlink(obj_after_ip)) {
1099  // target object is already target of other IP, so it must be source object
1100  std::swap(obj_before_ip, obj_after_ip);
1101  }
1102  //DebugOut().fmt("before: {} ip: {} after: {}\n", obj_before_ip, ip_idx, obj_after_ip );
1103  ASSERT_DBG( ! have_backlink(obj_after_ip) )
1104  (mesh_->element_accessor(intersection_->component_ele_idx()).idx())
1105  (mesh_->element_accessor(intersection_->bulk_ele_idx()).idx())
1106  (obj_before_ip)(ip_idx)(obj_after_ip); // at least one could be target object
1107  object_next[obj_before_ip] = ip_idx;
1108  IP_next.push_back( obj_after_ip);
1109  if (object_next[obj_after_ip] == no_idx) {
1110  object_next[obj_after_ip] = ip_idx;
1111  }
1112 }
1113 
1114 
1115 
1117 {
1118  intersection_= &intersection;
1119  //DebugOut().fmt("2d ele: {} 3d ele: {}\n",
1120  // intersection.component_ele_idx(),
1121  // intersection.bulk_ele_idx());
1122 
1123  IP23_list.clear();
1124  IP_next.clear();
1125  std::fill(object_next.begin(), object_next.end(), no_idx);
1126  std::vector<IPAux13> IP13s;
1127 
1128  std::array<bool, 6> edge_touch={false,false, false, false, false, false};
1129  unsigned int object_before_ip, object_after_ip;
1130 
1131  //unsigned int last_triangle_vertex=30; // no vertex at last IP
1132 
1133  // pass through the ccwise oriented sides in ccwise oriented order
1134  // How to make this in independent way?
1135  // Move this into RefElement?
1136  std::vector<unsigned int> side_cycle_orientation = { 0, 0, 1};
1137  std::vector<unsigned int> cycle_sides = {0, 2, 1};
1138 
1139  // TODO:
1140  // better mechanism for detecting vertex duplicities, do no depend on cyclic order of sides
1141  // still need cyclic orientation
1142  for(unsigned int _i_side = 0; _i_side < RefElement<2>::n_lines; _i_side++) { // go through triangle lines
1143  unsigned int i_side = cycle_sides[_i_side];
1144  IP13s.clear();
1145  CI13[ i_side ].compute(IP13s);
1146  ASSERT_DBG(IP13s.size() < 3);
1147  if (IP13s.size() == 0) continue;
1148  for(unsigned int _ip=0; _ip < IP13s.size(); _ip++) {
1149  //int ip_topo_position = _ip*2-1; // -1 (incoming ip), +1 (outcoming ip), 0 both
1150 
1151  // fix order of IPs
1152  unsigned int ip = (side_cycle_orientation[_i_side] + _ip) % IP13s.size();
1153 
1154 // DebugOut().fmt("rside: {} cside: {} rip: {} cip: {}", _i_side, i_side, _ip, ip);
1155 
1156  // convert from 13 to 23 IP
1157  IntersectionPointAux<3,1> IP31 = IP13s[ip].switch_objects(); // switch idx_A and idx_B and coords
1158  IntersectionPointAux<3,2> IP32(IP31, i_side); // interpolation uses local_bcoords_B and given idx_B
1159  IPAux23 IP23 = IP32.switch_objects(); // switch idx_A and idx_B and coords back
1160  //DebugOut() << IP;
1161 
1162  // Tracking info
1163  unsigned int tetra_object = s3_dim_starts[IP23.dim_B()] + IP23.idx_B();
1164  unsigned int side_object = s2_dim_starts[1] + i_side;
1165 
1166  object_before_ip = tetra_object;
1167  object_after_ip = side_object;
1168 
1169 
1170  // IP is vertex of triangle,
1171  if( IP23.dim_A() == 0 && IP23.dim_B() == 3)
1172  {
1173  // we are on line of the triangle, and IP.idx_A contains local node of the line
1174  // E-E, we know vertex index
1175  object_before_ip = s2_dim_starts[0]+IP23.idx_A();
1176  }// else current_triangle_vertex=3+IP23_list.size(); // no vertex, and unique
1177 
1178  // side of triangle touching S3, in vertex or in edge
1179  if (IP13s.size() == 1 ) {
1180  if (IP23.dim_B() == 0) {
1181  continue; // skip, S3 vertices are better detected in phase 2
1182  }
1183  if (IP23.dim_A() == 0) { // vertex of triangle
1184  object_before_ip = tetra_object;
1185  object_after_ip = s2_dim_starts[0]+IP23.idx_A();
1186 
1187  // source vertex of the side vector (oriented CCwise)
1188  //if ( (IP.idx_A()+side_cycle_orientation[_i_side])%2 == 0)
1189  // std::swap(object_before_ip, object_after_ip);
1190  } else {
1191  // touch in edge
1192 
1193  //continue;
1194  ASSERT_EQ_DBG(IP23.dim_B(), 1);
1195  edge_touch[IP23.idx_B()]=true;
1196  std::swap(object_before_ip, object_after_ip);
1197  }
1198  }
1199 
1200  IP23_list.push_back(IP23);
1201 
1202  unsigned int ip_idx = IP23_list.size()-1;
1203  //DebugOut().fmt("before: {} after: {} ip: {}\n", object_before_ip, object_after_ip, ip_idx);
1204  ASSERT_EQ_DBG(IP23_list.size(), IP_next.size()+1);
1205  set_links(object_before_ip, ip_idx, object_after_ip);
1206  }
1207 
1208  }
1209 
1210  // now we have at most single true degenerate IP in IP23
1211  // TODO:
1212  // - deal with degenerate IPs in the edge-trinagle phase
1213  // - remove degenerate_list (just count degen points)
1214  // - remove check for duplicities in final list copy
1215  // - add more comment to individual cases in order to be sure that any case in particular branch is
1216  // treated right
1217 
1218 
1219  //TODO: cannot the two cycles be merged in one now?
1220  // and instead of processed_edge use CI12[tetra_edge].set_computed();
1221  // and I think we don't have to generate dummy IP12s
1222  // so we don't have to have IP12s vector here
1223  IP12s_.clear();
1224  // S3 Edge - S2 intersections; collect all signs, make dummy intersections
1225  for(unsigned int tetra_edge = 0; tetra_edge < 6; tetra_edge++) {
1226  IPAux12 IP12;
1227  IntersectionResult result = CI12[tetra_edge].compute(IP12);
1228 // DebugOut() << print_var(tetra_edge) << print_var(int(result));
1229  // in degenerate case: IP12 is empty with degenerate result
1230  IP12s_.push_back(IP12);
1231  }
1232  vector<uint> processed_edge(6, 0);
1233  FacePair face_pair;
1234  for(unsigned int tetra_edge = 0; tetra_edge < 6;tetra_edge++) {
1235  if (! processed_edge[tetra_edge]) {
1236 // DBGVAR(tetra_edge);
1237  IPAux12 &IP12 = IP12s_[tetra_edge];
1238  if(IP12.result() >= IntersectionResult::degenerate) continue;
1239 
1240  int sign = CI12[tetra_edge].check_abscissa_topology(IP12);
1241 // DBGVAR(sign);
1242  if(std::abs(sign) > 1) continue;
1243 
1244  IPAux23 IP23(IP12s_[tetra_edge].switch_objects(), tetra_edge);
1245 
1246  const uint edge_dim = IP23.dim_B();
1247  const uint i_edge = IP23.idx_B();
1248  ASSERT_LT_DBG(edge_dim, 2);
1249 
1250  if ( edge_dim == 0) {
1251  face_pair = vertex_faces(i_edge);
1252  // mark edges coincident with the vertex
1253  for( uint ie : RefElement<3>::interact(Interaction<1,0>(i_edge)) )
1254  processed_edge[ie] = 1;
1255  }
1256  else
1257  face_pair = edge_faces(i_edge);
1258 
1259  //DebugOut() << print_var(face_pair[0])<< print_var(face_pair[1]);
1260 
1261  IP23_list.push_back(IP23);
1262  unsigned int ip_idx = IP23_list.size()-1;
1263 
1264  unsigned int s3_object = s3_dim_starts[edge_dim] + i_edge;
1265 
1266  //DebugOut() << print_var(edge_touch[i_edge]) << print_var(have_backlink(s3_object));
1267  if (IP23.dim_A() < 2
1268  && (! edge_touch[i_edge])
1269  && object_next[s3_object] != no_idx) { // boundary of S2, these ICs are duplicate
1270 
1271  if ( have_backlink(s3_object) ) {
1272  set_links(s3_object, ip_idx, face_pair[1]);
1273  } else {
1274  set_links(face_pair[0], ip_idx, s3_object);
1275  }
1276 
1277  } else { // interior of S2, just use the face pair
1278  //DebugOut() << print_var(face_pair[0])<< print_var(face_pair[1]);
1279  set_links(face_pair[0], ip_idx, face_pair[1]);
1280 
1281  if ( have_backlink(s3_object) ) {
1282  object_next[s3_object]=ip_idx;
1283  }
1284  }
1285  }
1286  }
1287 
1288 
1289  // Return IPs in correct order and remove duplicates
1290  ASSERT_EQ(0, intersection.size());
1291 
1292  if (IP23_list.size() == 0) return; // empty intersection
1293 
1294  // detect first IP, this needs to be done only in the case of
1295  // point or line intersections, where IPs links do not form closed cycle
1296  // Possibly we do this only if we detect such case through previous phases.
1297  vector<char> have_predecessor(IP23_list.size(), 0);
1298  for(auto obj : IP_next) {
1299  ASSERT_LT_DBG(obj, object_next.size());
1300  unsigned int ip = object_next[obj];
1301  if (ip < IP_next.size()) have_predecessor[ip]=1;
1302  }
1303  unsigned int ip_init=0;
1304  for(unsigned int i=0; i< IP23_list.size(); i++) if (! have_predecessor[i]) ip_init=i;
1305 
1306  arma::uvec::fixed<RefElement<3>::n_sides> ips_face_counter;
1307  ips_face_counter.zeros();
1308 
1309  // regular case, merge duplicit IPs
1310  unsigned int ip=ip_init;
1311  ASSERT_EQ_DBG(IP_next.size(), IP23_list.size());
1312  intersection.points().push_back(IP23_list[ip]);
1313  //DebugOut() << print_var(ip) << IP23_list[ip];
1314  while (1) {
1315  //DebugOut() << print_var(ip) << IP23_list[ip];
1316 
1317  unsigned int object = IP_next[ip];
1318  //IP_next[ip]=no_idx;
1319  ASSERT_LT_DBG(object, object_next.size());
1320  ip = object_next[object];
1321  object_next[object]=no_idx;
1322  if ((ip == no_idx)) break;
1323  ASSERT_LT_DBG(ip, IP_next.size());
1324 
1325  if ( ! IP23_list[ip].topology_equal(intersection.points().back()) ) {
1326  IPAux23 &IP = IP23_list[ip];
1327  //DebugOut() << print_var(ip) << IP23_list[ip];
1328  intersection.points().push_back(IP);
1329  if(IP.dim_B() < 3)
1330  ips_face_counter += on_faces[IP.dim_B()][IP.idx_B()];
1331  }
1332  }
1333 
1334  if (intersection.points().size() == 1) return;
1335 
1336  if (IP23_list[ip_init].topology_equal(intersection.points().back()) )
1337  intersection.points().pop_back();
1338 
1339  if (intersection.points().size() > 2){
1340  for(uint i=0; i< ips_face_counter.n_elem; i++){
1341  if(ips_face_counter(i) >= intersection.points().size()){ // all IPs lie in a single face
1342 // DBGCOUT(<< "all IPs in face: " << i <<"\n");
1343 // ips_face_counter.print("face_counter");
1344  intersection.set_ips_in_face(i);
1345  }
1346  }
1347  }
1348 }
1349 
1351 {
1352  auto &line_faces=RefElement<3>::interact(Interaction<2,1>(i_edge));
1353  unsigned int ip_ori = (unsigned int)(IP12s_[i_edge].result());
1354  ASSERT_DBG(ip_ori < 2); // no degenerate case
1355 
1356  // RefElement returns edge faces in clockwise order (edge pointing to us)
1357  // negative ip sign (ori 0) = faces counter-clockwise
1358  // positive ip sign (ori 1) = faces clockwise
1359  return { s3_dim_starts[2] + line_faces[1-ip_ori], s3_dim_starts[2] + line_faces[ip_ori] };
1360 }
1361 
1363 {
1364  // vertex edges clockwise
1365  const IdxVector<3> &vtx_edges = RefElement<3>::interact(Interaction<1,0>(i_vertex));
1366  std::array<unsigned int, 3> n_ori, sum_idx;
1367  n_ori.fill(0);
1368  sum_idx.fill(0);
1369  for(unsigned int ie=0; ie <3; ie++) {
1370  unsigned int edge_ip_ori = (unsigned int)(IP12s_[ vtx_edges[ie]].result());
1371  if (RefElement<3>::interact(Interaction<0,1>(vtx_edges[ie]))[0] != i_vertex
1372  && edge_ip_ori!= int(IntersectionResult::degenerate) )
1373  edge_ip_ori = (edge_ip_ori +1)%2;
1374  //ASSERT_LT_DBG(edge_ip_ori, 3)(ie);
1375  if (edge_ip_ori == 3) edge_ip_ori=2; // none as degenerate
1376  n_ori[edge_ip_ori]++;
1377  sum_idx[edge_ip_ori]+=ie;
1378  }
1379  unsigned int n_degen = n_ori[ int(IntersectionResult::degenerate) ];
1380  unsigned int sum_degen = sum_idx[ int(IntersectionResult::degenerate) ];
1381  unsigned int n_positive = n_ori[ int(IntersectionResult::positive) ];
1382  unsigned int n_negative= n_ori[ int(IntersectionResult::negative) ];
1383  //DebugOut().fmt("nd: {} sd: {} np: {} nn: {}", n_degen, sum_degen, n_positive, n_negative);
1384  if ( n_degen == 2 ) {
1385  // S2 plane match a face of S3, we treat degenerated edges as the faces
1386  // incident with the single regualr edge.
1387 
1388  unsigned int i_edge = 3 - sum_degen; // regular edge index
1389  FacePair pair = edge_faces(vtx_edges[i_edge]);
1390  auto &vtx_faces = RefElement<3>::interact(Interaction<2,0>(i_vertex));
1391  // replace faces by edges
1392  if (pair[0] == s3_dim_starts[2] + vtx_faces[(i_edge+1)%3])
1393  return { s3_dim_starts[1] + (i_edge+2)%3, s3_dim_starts[1] + (i_edge+1)%3 };
1394  else
1395  return { s3_dim_starts[1] + (i_edge+1)%3, s3_dim_starts[1] + (i_edge+2)%3 };
1396 
1397  } else if (n_degen == 1) {
1398  // One edge in S2 plane.
1399  unsigned int i_edge = sum_degen;
1400  ASSERT( n_positive + n_negative == 2);
1401  if ( n_positive == 1) {
1402  // opposite signs, S2 plane cuts S3
1403 
1404  FacePair pair = edge_faces(vtx_edges[(i_edge+1)%3]);
1405  unsigned int face = RefElement<3>::interact(Interaction<2,0>(i_vertex))[i_edge];
1406  // assign edges to faces
1407  //DebugOut().fmt("vtx: {} edg: {} face: {}", i_vertex, i_edge, face);
1408  if (pair[0] == s3_dim_starts[2] + face)
1409  return { s3_dim_starts[2] + face, s3_dim_starts[1] + vtx_edges[i_edge]};
1410  else
1411  return { s3_dim_starts[1] + vtx_edges[i_edge], s3_dim_starts[2] + face };
1412  } else {
1413  // same signs; S2 plane touch S3 vertex and a single edge
1414  //DebugOut() << "Touch in edge.";
1415  // same signs S2 plane touchs S3
1416  ASSERT(n_positive == 0 || n_positive== 2);
1417  return { s3_dim_starts[0]+i_vertex, s3_dim_starts[1] + vtx_edges[i_edge]};
1418  }
1419 
1420 
1421  } else {
1422  ASSERT(n_degen == 0);
1423  ASSERT( n_positive + n_negative == 3);
1424 
1425  if (n_positive == 1) {
1426  unsigned int i_edge = sum_idx[ int(IntersectionResult::positive) ];
1427  return edge_faces(vtx_edges[i_edge]);
1428  } else if (n_negative == 1) {
1429  unsigned int i_edge = sum_idx[ int(IntersectionResult::negative) ];
1430  return edge_faces(vtx_edges[i_edge]);
1431  } else {
1432  // S2 touch vertex of S3 in
1433  ASSERT( n_positive == 0 || n_positive == 3);
1434  return { s3_dim_starts[0]+i_vertex, s3_dim_starts[0]+i_vertex};
1435  }
1436  }
1437 }
1438 
1439 
1441 {
1443 
1444  on_faces.resize(3);
1445  arma::uvec::fixed<RefElement<3>::n_sides> v; v.zeros();
1446 
1447  on_faces[0].resize(RefElement<3>::n_nodes, v);
1448  for(uint i=0; i<RefElement<3>::n_nodes; i++)
1449  for(uint j=0; j<RefElement<3>::n_sides_per_node; j++)
1450  on_faces[0][i](RefElement<3>::interact(Interaction<2,0>(i))[j]) = 1;
1451 
1452  on_faces[1].resize(RefElement<3>::n_lines, v);
1453  for(uint i=0; i<RefElement<3>::n_lines; i++)
1454  for(uint j=0; j<RefElement<3>::n_sides_per_line; j++)
1455  on_faces[1][i](RefElement<3>::interact(Interaction<2,1>(i))[j]) = 1;
1456 
1457  on_faces[2].resize(RefElement<3>::n_sides, v);
1458  for(uint i=0; i<RefElement<3>::n_sides; i++)
1459  on_faces[2][i](i) = 1;
1460 
1461 // DBGCOUT("Print on_faces:\n");
1462 // for(uint d=0; d<on_faces.size(); d++)
1463 // for(uint i=0; i<on_faces[d].size(); i++){
1464 // for(uint j=0; j<on_faces[d][i].n_elem; j++)
1465 // cout << on_faces[d][i](j) << " ";
1466 // cout << endl;
1467 // }
1468  return on_faces;
1469 }
1470 
1471 
1473  for(unsigned int i = 0; i < 3;i++){
1474  os << "\tPluckerCoordinates Triangle[" << i << "]";
1475  if(plucker_coordinates_triangle_[i] == nullptr){
1476  os << "NULL" << endl;
1477  }else{
1478  os << plucker_coordinates_triangle_[i];
1479  }
1480  }
1481  for(unsigned int i = 0; i < 6;i++){
1482  os << "\tPluckerCoordinates Tetrahedron[" << i << "]";
1483  if(plucker_coordinates_tetrahedron[i] == nullptr){
1484  os << "NULL" << endl;
1485  }else{
1486  os << plucker_coordinates_tetrahedron[i];
1487  }
1488  }
1489 };
1490 
1492  os << "ComputeIntersection<2,3> Plucker Coordinates Tree:" << endl;
1493  print_plucker_coordinates(os);
1494  for(unsigned int i = 0; i < 6;i++){
1495  os << "ComputeIntersection<1,2>["<< i <<"] Plucker Coordinates:" << endl;
1496  CI12[i].print_plucker_coordinates(os);
1497  }
1498  for(unsigned int i = 0; i < 3;i++){
1499  CI13[i].print_plucker_coordinates_tree(os);
1500  }
1501 };
1502 
unsigned int dim_A() const
Returns dimension of object A.
#define ASSERT_EQ_DBG(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:331
void set_topology_B(unsigned int idx, unsigned int dim_B)
Sets the topology of object B in Simplex<M>.
std::vector< IntersectionPointAux< dimA, dimB > > i_points_
Vector of internal intersection points.
unsigned int uint
Plucker coordinates class.
IntersectionResult result() const
Result: 0 - negative sign, 1 - positive sign, 2 - degenerate (zero for all sides), 3 - none.
static BaryPoint line_barycentric_interpolation(BaryPoint first_coords, BaryPoint second_coords, double first_theta, double second_theta, double theta)
Definition: ref_element.cc:552
Fundamental simplicial intersections.
void set_topology_A(unsigned int idx, unsigned int dim_A)
Sets the topology of object A in Simplex<N>.
Definition: mesh.h:76
unsigned int size() const
Returns number of intersection points.
static const double geometry_epsilon
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
#define ASSERT(expr)
Allow use shorter versions of macro names if these names is not used with external library...
Definition: asserts.hh:346
const arma::vec::fixed< N+1 > & local_bcoords_A() const
Returns barycentric coordinates in the Simplex<N>.
static unsigned int normal_orientation(unsigned int sid)
Definition: ref_element.cc:315
unsigned int idx_A() const
Returns the index of Simplex<N>.
unsigned int dim_B() const
Returns dimension of object B.
static LocalPoint bary_to_local(const BaryPoint &bp)
Converts from barycentric to local coordinates.
Definition: ref_element.cc:300
unsigned int dim() const
Definition: elements.h:124
void set_coordinates(const arma::vec::fixed< N+1 > &lcA, const arma::vec::fixed< M+1 > &lcB)
Setter for coordinates.
static const IdxVector< (InDim >OutDim?InDim+1:dim-InDim) > interact(TInteraction< OutDim, InDim > interaction)
std::vector< IntersectionPointAux< dimA, dimB > > & points()
Returns intersection points by a reference.
void swap(nlohmann::json &j1, nlohmann::json &j2) noexcept(is_nothrow_move_constructible< nlohmann::json >::value andis_nothrow_move_assignable< nlohmann::json >::value)
exchanges the values of two JSON objects
Definition: json.hpp:8688
std::array< unsigned int, Size > IdxVector
Definition: ref_element.hh:146
bool topology_equal(const IntersectionPointAux< N, M > &other) const
Returns true, if other intersection point has the same topology.
IntersectionPointAux< M, N > switch_objects() const
Switches the object A and B.
#define ASSERT_DBG(expr)
Definition: asserts.hh:349
static const double plucker_empty
Auxiliary value for Plucker product. If equal this value, it is supposed not to be set yet...
void set_ips_in_face(unsigned int face_idx)
Internal auxiliary class represents an intersection point of simplex<N> and simplex<M>.
Class RefElement defines numbering of vertices, sides, calculation of normal vectors etc...
unsigned int idx_B() const
void set_result(IntersectionResult result)
Setter orientation flag.
Internal class representing intersection point.
Plucker coordinates representing line given by points A,B.
Definition: plucker.hh:44
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual)
Definition: asserts.hh:327
double scale() const
Definition: plucker.hh:72
Internal class representing intersection object.
#define ASSERT_LT_DBG(a, b)
Definition of comparative assert macro (Less Than) only for debug mode.
Definition: asserts.hh:299
const Node * node(unsigned int ni) const
Definition: accessors.hh:145