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