Flow123d  release_2.2.0-914-gf1a3a4f
intersection.cpp
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file intersection.cpp
15  * @brief
16  */
17 
18 #include "system/exc_common.hh"
23 #include <cmath>
24 //#include "math.h"
25 #include <iostream>
26 #include <armadillo>
28 //#include "mesh/ngh/include/problem.h"
29 
30 using namespace mathfce;
31 
32 namespace ngh {
33 
34 void GetIntersection(const TBisector & B1, const TBisector &B2,
35  TPosition &pos, double &t1, double &t2) {
36  TNSolutions ns;
37  TMatrix A(2);
38  TMVector B(2);
39  TMVector X(2);
40  double U[ 3 ], V[ 3 ];
41 
42  if (AreParallel(B1.GetVector(), B2.GetVector())) {
43  if (B1.Belong(B2.GetPoint())) {
44  pos = same;
45  return;
46  } else {
47  pos = parallel;
48  return;
49  }
50  }
51 
52  B1.GetVector().Get(U);
53  B2.GetVector().Get(V);
54 
55  int r1 = -1;
56  int r2 = -1;
57  for (int i = 0; i < 3; i++) {
58  if (!IsZero(U[ i ])) {
59  r1 = i;
60  }
61  if (!IsZero(V[ i ])) {
62  r2 = i;
63  }
64  }
65 
66  if (r1 == r2) {
67  for (int i = 0; i < 3; i++) {
68  if ((i != r1) && !IsZero(U[ i ])) {
69  r1 = i;
70  break;
71  }
72  if ((i != r2) && !IsZero(V[ i ])) {
73  r2 = i;
74  break;
75  }
76  }
77  }
78 
79  if ((r1 == -1) || (r2 == -1)) {
80  pos = skew;
81  return;
82  }
83 
84  int r3 = -1;
85  for (int i = 0; i < 3; i++) {
86  if ((r1 != i) && (r2 != i)) {
87  r3 = i;
88  break;
89  }
90  }
91 
92  A.Set(1, 1, U[ r1 ]);
93  A.Set(1, 2, -V[ r1 ]);
94  A.Set(2, 1, U[ r2 ]);
95  A.Set(2, 2, -V[ r2 ]);
96  B.Set(1, B2.GetPoint().Get(r1 + 1) - B1.GetPoint().Get(r1 + 1));
97  B.Set(2, B2.GetPoint().Get(r2 + 1) - B1.GetPoint().Get(r2 + 1));
98 
99  ns = Gauss(A, &X, B);
100 
101  if ((ns == no_solution) || (ns == inf_solutions)) {
102  pos = skew;
103  return;
104  }
105 
106  if (IsEqual(U[ r3 ] * X.Get(1) - V[ r3 ] * X.Get(2),
107  B2.GetPoint().Get(r3 + 1) - B1.GetPoint().Get(r3 + 1))) {
108  t1 = X.Get(1);
109  t2 = X.Get(2);
110  pos = intersecting;
111  } else {
112  pos = skew;
113  }
114 }
115 /*
116 void GetIntersection(const TBisector & B1, const TBisector &B2,
117  TPosition &pos, TPoint *P) {
118  double t1, t2;
119  GetIntersection(B1, B2, pos, t1, t2);
120  if (pos != intersecting)
121  return;
122  *P = B1.GetPoint(t1);
123  return;
124 }
125 */
126 /*
127 void GetIntersection(const TAbscissa &A1, const TAbscissa &A2,
128  TPosition &pos, TPoint *P) {
129  double t1, t2;
130  if (!QuickIntersectionTest(A1, A2)) {
131  pos = skew;
132  return;
133  }
134  GetIntersection(A1, A2, pos, t1, t2);
135  if (pos != intersecting)
136  return;
137  *P = A1.GetPoint(t1);
138  return;
139 }
140 */
141 
142 //------------------------UPRAVENO---------------------------
143 // vraci localni souradnice pruniku, prvni souradnice vzhledem k A1 druha souradnice vzhledem k A2
144 
145 void GetIntersection(const TAbscissa &A1, const TAbscissa &A2, IntersectionLocal * &insec) {
146  double t1, t2;
147  TPosition pos;
148 
149  if (!QuickIntersectionTest(A1, A2)) {
150  insec=NULL;
151  return;
152  }
153  GetIntersection(A1, A2, pos, t1, t2);
154  if ( pos == intersecting ) {
155  // test t1 je (0-eps,1+eps) a t2 je z (0-eps,1+eps)
156  if ((t1 > (0 - epsilon)) && (t1 < (1 + epsilon)) && (t2 > (0 - epsilon)) && (t2 < (1 + epsilon))) {
158  vector<double> loc_coord_1(1,t1);
159  vector<double> loc_coord_2(1,t2);
160  insec->add_local_coord(loc_coord_1,loc_coord_2);
161  return;
162  } else {
163  insec = NULL;
164  return;
165  }
166  } else if ( pos == same ) { //A1 a A2 lezi na stejne primce
167  // nalezeni prunikove usecky X1, X2 ; Y1, Y2 ...koncove body usecky XY
168 
169  TAbscissa A2_copy(A2);
170  double dot_product;
171  // TEST - vzajemna orientace A1 a A2 (skalarni soucin vektoru A1, A2)
172  dot_product = Dot(A1.GetVector(), A2.GetVector());
173  if (dot_product < 0 ) { //if opacna orientace A1 a A2
174  TVector AX2((-1)*A2.GetVector());
175  A2_copy.SetVector(AX2);
176  }
177 
178  //vypocet lokalni souradnice koncu A2 vzhledem k A1
179  TVector Diff;
180  double loc_begin_A1, loc_end_A1, loc_begin_A2, loc_end_A2;
181  Diff=A2_copy.GetPoint()-A1.GetPoint();
182  loc_begin_A2=Diff.Length()/(A1.GetVector().Length());
183 
184  Diff=A2_copy.GetPoint(1)-A1.GetPoint();
185  loc_end_A2=Diff.Length()/(A1.GetVector().Length());
186 
187  //vypocet lokalni souradnice koncu A1 vzhledem k A2
188  Diff=A1.GetPoint()-A2_copy.GetPoint();
189  loc_begin_A1=Diff.Length()/(A2_copy.GetVector().Length());
190 
191  Diff=A1.GetPoint(1)-A2_copy.GetPoint();
192  loc_end_A1=Diff.Length()/(A2_copy.GetVector().Length());
193 
194  //X1...loc.souradnice X vzhledem k A1
195  //X2...loc.souradnice X vzhledem k A2
196  //Y1...loc.souradnice Y vzhledem k A1
197  //Y2...loc.souradnice Y vzhledem k A2
198  double X1=0, X2=0, Y1=0, Y2=0;
199  //1.possibility - A2 lezi pred A1
200  if ((loc_begin_A2 < 0) && (loc_end_A2 < 0)) {
201  insec=NULL;
202  return;
203  }
204  //2.possibility - prekryvaji se, A2 lezi pred A1
205  if ((loc_begin_A2 < 0) && (loc_end_A2 >= 0) && (loc_end_A2 < 1)) {
206  X1 = 0;
207  X2 = loc_begin_A1;
208  Y1 = loc_end_A2;
209  Y2 = 1;
210  }
211  //3.possibility - prekryvaji se, A2 lezi mezi koncovymi body A1
212  if ((loc_begin_A2 >= 0) && (loc_end_A2 >= 0) && (loc_end_A2 < 1)) {
213  X1 = loc_begin_A2;
214  X2 = 0;
215  Y1 = loc_end_A2;
216  Y2 = 1;
217  }
218  //4.possibility - prekryvaji se, A2 lezi za A1
219  if ((loc_begin_A2 >= 0) && (loc_begin_A2 < 1) && (loc_end_A2 >= 1)) {
220  X1 = loc_begin_A2;
221  X2 = 0;
222  Y1 = 1;
223  Y2 = loc_end_A1;
224  }
225  //5.possibility - A2 lezi za A1
226  if ((loc_begin_A2 >= 1) && (loc_end_A2 >= 1)) {
227  insec=NULL;
228  return;
229  }
230  //6.possibility - prekryvaji se, A1 lezi mezi koncovymi body A2
231  if ((loc_begin_A2 < 0) && (loc_end_A2 >= 1)) {
232  X1 = 0;
233  X2 = loc_begin_A1;
234  Y1 = 1;
235  Y2 = loc_end_A1;
236  }
237  OLD_ASSERT(X1==0 && X2==0 && Y1==0 && Y2==0,"Uncovered situation.");
238 
239  //set local coords:
241  vector<double> loc_coord_1(1);
242  vector<double> loc_coord_2(1);
243  loc_coord_1[0] = X1;
244  loc_coord_2[0] = X2;
245  insec->add_local_coord(loc_coord_1,loc_coord_2);
246  loc_coord_1[0] = Y1;
247  loc_coord_2[0] = Y2;
248  insec->add_local_coord(loc_coord_1,loc_coord_2);
249  return;
250 
251  } else {
252  insec = NULL;
253  return;
254  }
255  insec = NULL;
256  return;
257 }
258 /*
259 void GetIntersection(const TAbscissa &A, const TBisector &B,
260  TPosition &pos, TPoint *P) {
261  double t1, t2;
262  GetIntersection(A, B, pos, t1, t2);
263  if (pos != intersecting)
264  return;
265  *P = A.GetPoint(t1);
266  return;
267 }
268 */
269 
270 //------------------------UPRAVENO---------------------------
271 // vraci localni souradnice pruniku, prvni souradnice vzhledem k A druha souradnice vzhledem k B
272 // v pripade ze je prunik usecka - poradi bodu podle orientace abscissa A
273 
274 void GetIntersection(const TAbscissa &A, const TBisector &B, IntersectionLocal * &insec) {
275  double t1, t2;
276  TPosition pos;
277  insec = NULL;
278 
279  GetIntersection(A, B, pos, t1, t2);
280  if ( pos == intersecting ) {
281  // test t1 je (0-eps,1+eps)
282  if ((t1 > (0 - epsilon)) && (t1 < (1 + epsilon))) {
284  vector<double> loc_coord_1(1,t1);
285  vector<double> loc_coord_2(1,t2); //t2 na Bisectoru B
286  insec->add_local_coord(loc_coord_1,loc_coord_2);
287  }
288  } else if ( pos == same ) { //A a B lezi na stejne primce
289  // prunik cela usecka =>ulozit koncove body Abscissy A
290 
291  //vypocet lokalni souradnice koncovych bodu A
292  bool x;
293 
294  // get local coords of points of A on B
295  B.GetParameter(A.GetPoint(), t1, x);
296  B.GetParameter(A.GetPoint(1), t2, x);
297 
298  // sort t1, t2
299  if (t2 < t1) {
300  double swap = t1;
301  t1 = t2;
302  t2 = swap;
303  }
304 
305  //set local coords:
306  if ((t2 > (0 - epsilon)) && (t1 < (1 + epsilon))) {
307  double loc_begin_A, loc_end_A;
308  loc_begin_A = (t1 > 0) ? t1 : 0;
309  loc_end_A = (t2 < 1) ? t2 : 1;
310 
312  // first abscissa point
313  vector<double> loc_coord_1(1,0);
314  vector<double> loc_coord_2(1,loc_begin_A);
315  insec->add_local_coord(loc_coord_1,loc_coord_2);
316  // second abscissa point
317  vector<double> loc_coord_1_(1,1);
318  vector<double> loc_coord_2_(1,loc_end_A);
319  insec->add_local_coord(loc_coord_1_,loc_coord_2_);
320  return;
321  }
322  } else {
323  insec=NULL;
324  return;
325  }
326  return;
327 }
328 
329 /*
330 void GetIntersection(const TBisector &B, const TAbscissa &A,
331  TPosition &pos, TPoint *P) {
332  GetIntersection(A, B, pos, P);
333  return;
334 }
335 */
336 void GetIntersection(const TBisector &B, const TAbscissa &A, IntersectionLocal * &insec) {
337  GetIntersection(A, B, insec); //KONTROLA
338  return;
339 }
340 
341 void GetIntersection(const TAbscissa &A1, const TAbscissa &A2,
342  TPosition &pos, double &t1, double &t2) {
343 
344  GetIntersection( (const TBisector &)A1, (const TBisector &)A2, pos, t1, t2);
345  if (pos == intersecting)
346  if (t1 > 1 + epsilon || t1 < 0 - epsilon ||
347  t2 > 1 + epsilon || t2 < 0 - epsilon)
348  pos = skew;
349  return;
350 }
351 
352 void GetIntersection(const TAbscissa &A, const TBisector &B,
353  TPosition &pos, double &t1, double &t2) {
354 
355  GetIntersection( (const TBisector &)A, B, pos, t1, t2);
356  if (pos == intersecting)
357  if (t1 > 1 + epsilon || t1 < 0 - epsilon)
358  pos = skew;
359  return;
360 }
361 
362 void GetIntersection(const TBisector &B, const TAbscissa &A,
363  TPosition &pos, double &t2, double &t1) {
364  GetIntersection(A, B, pos, t1, t2);
365  return;
366 }
367 
368 double Distance(const TBisector & B, const TPoint &P) {
369  double d;
370  TVector U(P, B.GetPoint());
371  d = Cross(U, B.GetVector()).Length() / B.GetVector().Length();
372  return d;
373 }
374 
375 double Distance(const TPlain &P, const TPoint &X) {
376  double dis;
377  dis = fabs(P.GetA() * X.X() + P.GetB() * X.Y() + P.GetC() * X.Z() + P.GetD()) /
378  sqrt(pow(P.GetA(), 2) + pow(P.GetB(), 2) + pow(P.GetC(), 2));
379  return dis;
380 }
381 
382 double Distance(const TPoint &P1, const TPoint &P2) {
383  double dis;
384  dis = sqrt(pow(P1.X() - P2.X(), 2) + pow(P1.Y() - P2.Y(), 2) +
385  pow(P1.Z() - P2.Z(), 2));
386  return dis;
387 }
388 
389 void GetIntersection(const TPlain &P1, const TPlain &P2,
390  TPosition &pos, TBisector *B) {
391  TVector U;
392  TMatrix M(3, 3);
393  TMVector b(3);
394  TMVector x(3);
395  int u1, v1, u2, v2;
396  int i, cit;
397  bool test;
398  if (AreParallel(P1.GetNormal(), P2.GetNormal())) {
399  if (P1.Belong(P2.GetPoint())) {
400  pos = same;
401  return;
402  } else {
403  pos = parallel;
404  return;
405  }
406  }
407  U = Cross(P1.GetNormal(), P2.GetNormal());
408  B->SetVector(U);
409  for (i = 1; i <= 3; i++)
410  b.Set(i, P2.GetPoint().Get(i) - P1.GetPoint().Get(i));
411  test = false;
412  cit = 0;
413  while (!test) {
414  u1 = v1 = u2 = v2 = -1;
415  switch (cit) {
416  case 0: u1 = 1;
417  u2 = 2;
418  v1 = 3;
419  break;
420  case 1: u1 = 1;
421  u2 = 2;
422  v2 = 3;
423  break;
424  case 2: u2 = 1;
425  v1 = 2;
426  v2 = 3;
427  break;
428  case 3: u1 = 1;
429  v1 = 2;
430  v2 = 3;
431  break;
432  default:
433  if (P1.Belong(P2.GetPoint())) {
434  pos = same;
435  return;
436  } else {
437  pos = parallel;
438  return;
439  }
440  // mythrow("The two planes should have intersection bisector",__LINE__, __FUNC__);
441  }
442  for (i = 1; i <= 3; i++) {
443  if (u1 != -1)
444  M.Set(i, u1, P1.GetU().Get(i));
445  if (v1 != -1)
446  M.Set(i, v1, P1.GetV().Get(i));
447  if (u2 != -1)
448  M.Set(i, u2, -P2.GetU().Get(i));
449  if (v2 != -1)
450  M.Set(i, v2, -P2.GetV().Get(i));
451  }
452  if (Gauss(M, &x, b) == one_solution)
453  test = true;
454  cit++;
455  }
456  if (u1 != -1 && v1 != -1)
457  B->SetPoint(P1.GetPoint(x.Get(u1), x.Get(v1)));
458  else
459  B->SetPoint(P2.GetPoint(x.Get(u2), x.Get(v2)));
460  pos = intersecting;
461  return;
462 }
463 
464 void GetIntersection(const TPlain &P, const TBisector &B,
465  TPosition &pos, double &t) {
466  arma::mat aa(3,3);
467  arma::vec xx(3);
468  arma::vec bb(3);
469  //TMatrix M(3, 3);
470  //TMVector b(3);
471  //TMVector x(3);
472 
473  int i;
474  xx.zeros();
475  if (ArePerpendicular(P.GetNormal(), B.GetVector())) {
476  if (P.Belong(B.GetPoint())) {
477  pos = belong;
478  return;
479  } else {
480  pos = parallel;
481  return;
482  }
483  }
484  for (i = 1; i <= 3; i++) {
485  /*M.Set(i, 1, B.GetVector().Get(i));
486  M.Set(i, 2, -P.GetU().Get(i));
487  M.Set(i, 3, -P.GetV().Get(i));
488  b.Set(i, P.GetPoint().Get(i) - B.GetPoint().Get(i));*/
489 
490  aa(i-1,0) = (B.GetVector().Get(i));
491  aa(i-1,1) = (-P.GetU().Get(i));
492  aa(i-1,2) = (-P.GetV().Get(i));
493  bb(i-1) = (P.GetPoint().Get(i) - B.GetPoint().Get(i));
494  }
495  arma::solve(xx, aa, bb);
496  pos = intersecting;
497  t = xx(0);
498  return;
499 }
500 
501 void GetIntersection(const TPlain &P, const TBisector &B,
502  TPosition &pos, TPoint *Pt) {
503  double t;
504  GetIntersection(P, B, pos, t);
505  *Pt = B.GetPoint(t);
506 }
507 
508 void GetIntersection(const TBisector &B, const TPlain &P,
509  TPosition &pos, TPoint *Pt) {
510  GetIntersection(P, B, pos, Pt);
511  return;
512 }
513 
514 void GetIntersection(const TTriangle &T1, const TTriangle &T2,
515  TIntersectionType &it, double &value) { //ZATIM NEPOTREBUJEME =>zakomentovano
516 
517  //******************************************************************************
518  //ONE SHOULD ADD TO THIS FUNCTION POINTS WHICH ARE INSIDE TRIANGLE
519  //******************************************************************************
520 /*
521  if (!QuickIntersectionTest(T1, T2)) {
522  it = none;
523  return;
524  }
525 
526  TPosition pos;
527  TPoint P1,P2;
528  TBisector b(P1,P2);
529  TPolygon pol;
530  double t11, t12, t21, t22, t[ 4 ], t1max, t1min, t2max, t2min;
531  int i, j, cit;
532  // TPoint P[6];
533  GetIntersection(T1.GetPlain(), T2.GetPlain(), pos, &b);
534  if (pos == parallel) {
535  it = none;
536  return;
537  }
538  if (pos == intersecting) {
539  it = line;
540  GetIntersection(b, T1, it, t11, t12);
541  if (it == none)
542  return;
543  GetIntersection(b, T2, it, t21, t22);
544  if (it == none)
545  return;
546  if (t11 > t12) {
547  t1max = t11;
548  t1min = t12;
549  } else {
550  t1max = t12;
551  t1min = t11;
552  }
553  if (t21 > t22) {
554  t2max = t21;
555  t2min = t22;
556  } else {
557  t2max = t22;
558  t2min = t21;
559  }
560  if (t1max < t2min || t2max < t1min) {
561  it = none;
562  return;
563  }
564  t[ 0 ] = t11;
565  t[ 1 ] = t12;
566  t[ 2 ] = t21;
567  t[ 3 ] = t22;
568  SortAsc(t, 4);
569  value = Distance(b.GetPoint(t[ 1 ]), b.GetPoint(t[ 2 ]));
570  return;
571  }
572  if (pos == same) {
573  for (i = 1; i <= 3; i++) {
574  cit = 0;
575  for (j = 1; j <= 3; j++) {
576  GetIntersection(T1.GetAbscissa(i), T2.GetAbscissa(j), pos, t11, t22);
577  if (pos == intersecting) {
578  t[ cit ] = t11;
579  cit++;
580  }
581  }
582  if (cit > 1)
583  SortAsc(t, 2);
584  for (j = 0; j < cit; j++)
585  pol.Add(T1.GetAbscissa(i).GetPoint(t[ j ]));
586  }
587  for (i = 1; i <= 3; i++) {
588  if (T1.IsInner(T2.GetPoint(i)))
589  pol.Add(T2.GetPoint(i));
590  if (T2.IsInner(T1.GetPoint(i)))
591  pol.Add(T1.GetPoint(i));
592  }
593  it = area;
594  value = pol.GetArea();
595  return;
596  }
597 */
598 }
599 
600 //******************************************************************************
601 // The value of it is important for computation
602 //******************************************************************************
603 
604 //------------------------UPRAVENO---------------------------
605 // vraci localni souradnice pruniku, prvni souradnice vzhledem k B druha souradnice vzhledem k T
606 // pro vsechny body pruniku B a T
607 
608 void GetIntersection(const TBisector &B, const TTriangle &T, IntersectionLocal * &insec) {
609  TPosition pos;
610  double t;
611  int cit=0;
612  GetIntersection(T.GetPlain(), B, pos, t);
613  switch (pos) {
614  // POINT INTERSECTION
615  case intersecting: {
616  if (!T.IsInner(B.GetPoint(t))) {
617  insec=NULL;
618  return;
619  }
620 
621  vector<double> loc_coord_1(1);
622  vector<double> loc_coord_2(2);
623 
624  loc_coord_1[0] = t;
625  loc_coord_2[0] = 0; //TODO: values of loc_coord_2 are not computed
626  loc_coord_2[1] = 0;
627 
629  insec->add_local_coord(loc_coord_1, loc_coord_2);
630 
631  return;
632  }
633 
634  // LINE INTERSECTION
635  case belong: {
636  IntersectionLocal* insec_tmp;
637  IntersectionPoint* insec_point_tmp[3];
638 
639  vector<double> loc_tria_tmp(2);
640  // inicializace vektoru pro lokalni souradnice pro pripad useckoveho pruseciku
641  // lokalni souradnice vrholu trojuhelnika
642  vector<double> loc_tria_coord_01(2); //loc. coord 1.vrcholu trojuhelnika
643  vector<double> loc_tria_coord_02(2); //loc. coord 2.vrcholu trojuhelnika
644  vector<double> loc_tria_coord_03(2); //loc. coord 3.vrcholu trojuhelnika
645 
646  loc_tria_coord_01[0]=0;
647  loc_tria_coord_01[1]=0;
648 
649  loc_tria_coord_02[0]=1;
650  loc_tria_coord_02[1]=0;
651 
652  loc_tria_coord_03[0]=0;
653  loc_tria_coord_03[1]=1;
654 
655  //PRVNI STENA TROJUHELNIKU
656  GetIntersection(T.GetAbscissa(1), B, insec_tmp);
657  if (insec_tmp != NULL) {
658  if (insec_tmp->get_type() == IntersectionLocal::point) {
659  // intersection in a point
660  // zde: 1) vybrat IntersectionPoint z insec_tmp
661  // 2) prohodit jeho coord
662  // 3) z druhe coord udelat lokalni souradnici na trojuhelniku
663 
664  // loc. triangle coord (X, 0)
665  loc_tria_tmp[0] = insec_tmp->get_point(0)->el1_coord()[0];
666  loc_tria_tmp[1] = 0;
667  insec_point_tmp[cit] = new IntersectionPoint(insec_tmp->get_point(0)->el2_coord(), loc_tria_tmp);
668  cit++; //citac kolik sten protne
669  //cout<<"\nCit(1.stena)= "<< cit << endl;
670  delete insec_tmp;
671  } else if (insec_tmp->get_type() == IntersectionLocal::line) {
672  // intersection is whole tringle side => ulozit lokalni souradnice strany trojuhelnika do insec
674  insec->add_local_coord(insec_tmp->get_point(0)->el2_coord(), loc_tria_coord_01);
675  insec->add_local_coord(insec_tmp->get_point(1)->el2_coord(), loc_tria_coord_02);
676  delete insec_tmp;
677  return;
678  }
679  }
680 
681  //DRUHA STENA TROJUHELNIKU
682  GetIntersection(T.GetAbscissa(2), B, insec_tmp);
683  if (insec_tmp != NULL) {
684  if (insec_tmp->get_type() == IntersectionLocal::point) {
685  // loc. triangle coord (1-X, X)
686  loc_tria_tmp[0] = 1 - insec_tmp->get_point(0)->el1_coord()[0];
687  loc_tria_tmp[1] = insec_tmp->get_point(0)->el1_coord()[0];
688  insec_point_tmp[cit] = new IntersectionPoint(insec_tmp->get_point(0)->el2_coord(), loc_tria_tmp);
689  cit++; //citac kolik sten protne
690  //cout<<"Cit(2.stena)= "<< cit << endl;
691  delete insec_tmp;
692  } else if (insec_tmp->get_type() == IntersectionLocal::line) {
693  // intersection is whole tringle side => ulozit lokalni souradnice strany trojuhelnika do insec
695  insec->add_local_coord(insec_tmp->get_point(0)->el2_coord(), loc_tria_coord_02);
696  insec->add_local_coord(insec_tmp->get_point(1)->el2_coord(), loc_tria_coord_03);
697  delete insec_tmp;
698  delete insec_point_tmp[0];
699  return;
700  }
701  }
702 
703  //TRETI STENA TROJUHELNIKU
704  GetIntersection(T.GetAbscissa(3), B, insec_tmp);
705  if (insec_tmp != NULL) {
706  if (insec_tmp->get_type() == IntersectionLocal::point) {
707  // loc. triangle coord (0, 1-X)
708  loc_tria_tmp[0] = 0;
709  loc_tria_tmp[1] = 1 - insec_tmp->get_point(0)->el1_coord()[0];
710  insec_point_tmp[cit] = new IntersectionPoint(insec_tmp->get_point(0)->el2_coord(), loc_tria_tmp);
711  cit++; //citac kolik sten protne
712  //cout <<"Cit(3.stena)= "<< cit << endl;
713  delete insec_tmp;
714  } else if (insec_tmp->get_type() == IntersectionLocal::line) {
715  // intersection is whole tringle side => ulozit lokalni souradnice strany trojuhelnika do insec
717  insec->add_local_coord(insec_tmp->get_point(0)->el2_coord(), loc_tria_coord_03);
718  insec->add_local_coord(insec_tmp->get_point(1)->el2_coord(), loc_tria_coord_01);
719  delete insec_tmp;
720  delete insec_point_tmp[0];
721  delete insec_point_tmp[1];
722  return;
723  }
724  }
725 
726  //TEST - pocet protnutych sten trojuhelnika
727  if (cit == 0) {
728  insec = NULL;
729  return;
730  }
731  //lezi pres vrchol a zaroven protne protilehlou stenu => cit==3
732  if (cit == 3) {
733  //cout<<"(cit == 3) => REDUKCE insec_point_tmp!" << endl;
734  if (*(insec_point_tmp[0]) == *(insec_point_tmp[0])) {
735  //cout<<"2 insec_point_tmp[0] jsou stejne" << endl;
736  }
737 
738  for (int i = 0; i < cit; i++) {
739  if (*(insec_point_tmp[i]) == *(insec_point_tmp[(i+1)%3])) {
740 
741  //cout<<"insec_point_tmp[(i+1)%3] :: el1_coord().size(): " << insec_point_tmp[(i+1)%3]->el1_coord().size() << endl;
742  //cout<<"(insec_point_tmp[(i+1)%3]) :: el2_coord().size(): " << (insec_point_tmp[(i+1)%3])->el2_coord().size() << endl;
743 
744  delete insec_point_tmp[(i+1)%3];
745  cit--;
746  //cout<<"PO REDUKCI cit= " << cit << endl;
747 
748  //nutno posunout zbyle prvky pole insec_point_tmp
749  IntersectionPoint* tmp = insec_point_tmp[(i+2)%3];
750  insec_point_tmp[0] = insec_point_tmp[i];
751  insec_point_tmp[1] = tmp;
752  break;
753  }
754  }
755  //return;
756  }
757  if (cit != 2) {
758  cout << "cit = " << cit << endl;
759  THROW( ExcAssertMsg() << EI_Message("Error - pocet bodu pruniku != 2.\n") );
760  return;
761  } else {
762  if (*(insec_point_tmp[0]) == *(insec_point_tmp[1])) { //lezi pres vrchol
764  insec->add_local_point(insec_point_tmp[0]);
765  //cout<<"Insec_point_1 == Insec_point_2" << endl;
766  delete insec_point_tmp[1];
767  } else {
769  insec->add_local_point(insec_point_tmp[0]);
770  insec->add_local_point(insec_point_tmp[1]);
771  //cout<<"(IntersectionType=line) - point_1->el1_coord().size()= " << insec->get_point(0)->el1_coord().size()<< endl;
772  //cout<<"(IntersectionType=line) - point_2->el1_coord().size()= " << insec->get_point(1)->el1_coord().size()<< endl;
773  }
774  return;
775  }
776  } //end case belong
777 
778  // EMPTY INTERSECTION
779  default:
780  insec = NULL;
781  return;
782  }
783 }
784 
785 //******************************************************************************
786 // The value of it is important for computation
787 //******************************************************************************
788 
789 //------------------------UPRAVENO---------------------------
790 // vraci localni souradnice pruniku, prvni souradnice vzhledem k A druha souradnice vzhledem k T
791 // pro vsechny body pruniku A a T
792 
793 void GetIntersection(const TAbscissa &A, const TTriangle &T,
794  IntersectionLocal * & insec) {
795 
796  if (!QuickIntersectionTest(A, T)) {
797  insec = NULL;
798  return;
799  }
800  IntersectionLocal* insec_tmp=NULL;
801  GetIntersection( (const TBisector &)A, T, insec_tmp);
802  if (!insec_tmp) {
803  insec = NULL;
804  return;
805  }
806  if (insec_tmp->get_type() == IntersectionLocal::point) {
807  if (insec_tmp->get_point(0) != NULL) {
808  double t1 = insec_tmp->get_point(0)->el1_coord()[0];
809  if (t1 < (0 - epsilon) || t1 > (1 + epsilon)) {
810  delete insec_tmp;
811  insec = NULL;
812  } else {
813  insec = insec_tmp;
814  }
815  } else{
816  delete insec_tmp;
817  insec = NULL;
818  }
819  } else if(insec_tmp->get_type() == IntersectionLocal::line) {
820  // A1 i A2 ma byt v intervalu (0,1) -> vrati insec
821  // pokud ne tak zkusi zkratit, nebo NULL (delete)
822 
823  const IntersectionPoint* A1;
824  const IntersectionPoint* A2;
825  if (insec_tmp->get_point(0) != NULL) {
826  if (insec_tmp->get_point(0)->el1_coord()[0] > insec_tmp->get_point(1)->el1_coord()[0]) {
827  A2 = insec_tmp->get_point(0);
828  A1 = insec_tmp->get_point(1);
829  } else {
830  A1 = insec_tmp->get_point(0);
831  A2 = insec_tmp->get_point(1);
832  }
833 
834  double A1_t = A1->el1_coord()[0];
835  double A2_t = A2->el1_coord()[0];
836  if (A1_t < 0) A1_t = 0;
837  if (A2_t > 1) A2_t = 1;
838 
839  if (A2_t < A1_t) {
840  delete insec_tmp;
841  insec = NULL;
842  } else {
844  //cout << "A1_t: " << A1_t << endl;
845  //cout << "A2_t: " << A2_t << endl;
846  insec->add_local_point(interpolate(*A1, *A2, A1_t));
847  insec->add_local_point(interpolate(*A1, *A2, A2_t));
848  delete insec_tmp;
849  }
850  } else{
851  delete insec_tmp;
852  insec = NULL;
853  }
854  }
855  return;
856 
857 //-------------------------------------------------------------------
858 /* IntersectionPoint* insec_point_tmp[2];
859  vector<double> loc_tria_tmp(2);
860  double A1_ = insec_tmp->get_point(0)->el1_coord()[0];
861  double A2_ = insec_tmp->get_point(1)->el1_coord()[0]; //test != NULL (lezi pres vrchol!)
862  double *A1 = &A1_;
863  double *A2 = &A2_;
864  bool invert = false;
865  //TEST VZAJEMNE ORIENTACE A1, A2:
866  if (A1_ > A2_) {
867  A1 = &A2_;
868  A2 = &A1_;
869  invert = true;
870  }
871  // pripady:
872  // 1) prusecik uvnitr usecky -> predat nezmenene
873  // 2) prusecik mimo usecku -> delete a vratit NULL
874  // 3) else a) A1 < 0 - zkratit zleva
875  // b) A2 > 1 - zkratit zprava
876 
877  //PRUSECIK UVNITR USECKY
878  if ((*A1 > 0 - epsilon) && (*A1 < 1 + epsilon) && (*A2 > 0 - epsilon) && (*A2 < 1 + epsilon)) {
879  insec = insec_tmp;
880  }
881  //PRUSECIK MIMO USECKU
882  if (((*A1 < 0 - epsilon) && (*A2 < 0 - epsilon)) || ((*A1 > 1 + epsilon) && (*A2 > 1 + epsilon))) {
883  delete insec_tmp;
884  insec = NULL;
885  } else {
886  //1.possibility (A1 < 0) - zkratit zleva
887  if ((*A1 < 0 - epsilon) && (*A2 > 0 - epsilon) && (*A2 < 1 + epsilon)) {
888  if(invert == true) {
889  //loc. coord. r:
890  loc_tria_tmp[0] = ((1 - *A2)/(*A1 - *A2))*(insec_tmp->get_point(1)->el2_coord()[0] - insec_tmp->get_point(0)->el2_coord()[0]) + insec_tmp->get_point(0)->el2_coord()[0];
891  //loc. coord. s:
892  loc_tria_tmp[1] = ((1 - *A2)/(*A1 - *A2))*(insec_tmp->get_point(1)->el2_coord()[1] - insec_tmp->get_point(0)->el2_coord()[1]) + insec_tmp->get_point(0)->el2_coord()[1];
893  } else {
894  //loc. coord. r:
895  loc_tria_tmp[0] = ((1 - *A2)/(*A1 - *A2))*(insec_tmp->get_point(0)->el2_coord()[0] - insec_tmp->get_point(1)->el2_coord()[0]) + insec_tmp->get_point(1)->el2_coord()[0];
896  //loc. coord. s:
897  loc_tria_tmp[1] = ((1 - *A2)/(*A1 - *A2))*(insec_tmp->get_point(0)->el2_coord()[1] - insec_tmp->get_point(1)->el2_coord()[1]) + insec_tmp->get_point(1)->el2_coord()[1];
898  }
899  vector<double> loc_A1(1, 0);
900  vector<double> loc_A2(1, *A2);
901  if(invert == true) {
902  insec_point_tmp[0] = new IntersectionPoint(loc_A2, insec_tmp->get_point(1)->el2_coord());
903  insec_point_tmp[1] = new IntersectionPoint(loc_A1, loc_tria_tmp);
904  } else {
905  insec_point_tmp[0] = new IntersectionPoint(loc_A2, insec_tmp->get_point(0)->el2_coord());
906  insec_point_tmp[1] = new IntersectionPoint(loc_A1, loc_tria_tmp);
907  }
908  insec = new IntersectionLocal(IntersectionLocal::line);
909  insec->add_local_point(insec_point_tmp[0]);
910  insec->add_local_point(insec_point_tmp[1]);
911  delete insec_tmp;
912  }
913  //2.possibility (A2 > 1) - zkratit zprava
914  if ((*A1 > 0 - epsilon) && (*A1 < 1 + epsilon) && (*A2 > 1 + epsilon)) {
915  if(invert == true) {
916  //loc. coord. r:
917  loc_tria_tmp[0] = ((1 - *A1)/(*A2 - *A1))*(insec_tmp->get_point(0)->el2_coord()[0] - insec_tmp->get_point(1)->el2_coord()[0]) + insec_tmp->get_point(1)->el2_coord()[0];
918  //loc. coord. s:
919  loc_tria_tmp[1] = ((1 - *A1)/(*A2 - *A1))*(insec_tmp->get_point(0)->el2_coord()[1] - insec_tmp->get_point(1)->el2_coord()[1]) + insec_tmp->get_point(1)->el2_coord()[1];
920  } else {
921  //loc. coord. r:
922  loc_tria_tmp[0] = ((1 - *A1)/(*A2 - *A1))*(insec_tmp->get_point(1)->el2_coord()[0] - insec_tmp->get_point(0)->el2_coord()[0]) + insec_tmp->get_point(0)->el2_coord()[0];
923  //loc. coord. s:
924  loc_tria_tmp[1] = ((1 - *A1)/(*A2 - *A1))*(insec_tmp->get_point(1)->el2_coord()[1] - insec_tmp->get_point(0)->el2_coord()[1]) + insec_tmp->get_point(0)->el2_coord()[1];
925  }
926  vector<double> loc_A1(1, *A1);
927  vector<double> loc_A2(1, 1);
928  if(invert == true) {
929  insec_point_tmp[0] = new IntersectionPoint(loc_A1, insec_tmp->get_point(1)->el2_coord());
930  insec_point_tmp[1] = new IntersectionPoint(loc_A2, loc_tria_tmp);
931  } else {
932  insec_point_tmp[0] = new IntersectionPoint(loc_A1, insec_tmp->get_point(0)->el2_coord());
933  insec_point_tmp[1] = new IntersectionPoint(loc_A2, loc_tria_tmp);
934  }
935  insec = new IntersectionLocal(IntersectionLocal::line);
936  insec->add_local_point(insec_point_tmp[0]);
937  insec->add_local_point(insec_point_tmp[1]);
938  delete insec_tmp;
939  }
940 
941  //3.possibility - zkratit z obou stran
942  //if ((*A1 < 0 - epsilon) && (*A2 > 1 + epsilon)) {
943  // vector<double> loc_A1(1, 0);
944  // vector<double> loc_A2(1, 1);
945  //}
946  }
947  }
948  return;
949  }
950  return;
951 */
952 }
953 
954 void GetIntersection(const TAbscissa &A, const TTetrahedron &T,
955  TIntersectionType &it, double &coef) {
956 
957  if (!QuickIntersectionTest(A, T)) {
958  it = none;
959  return;
960  }
961 
962  int cit = 0;
963  double tt1, tt2;
964  double tt[2];
965 
966  if (T.IsInner(A.GetPoint())) {
967  tt[cit] = 0;
968  cit++;
969  }
970  if (T.IsInner(A.GetPoint(1))) {
971  tt[cit] = 1;
972  cit++;
973  }
974  if (cit == 2) {
975  coef = fabs(tt[1] - tt[0]);
976  it = line;
977  return;
978  }
979 
980  IntersectionLocal *insec;
981 
982  for (int i = 1; i <= 4; i++) {
983  it = unknown;
984  GetIntersection(A, T.GetTriangle(i), insec);
985  if (insec) {
986  if (insec->get_type() == IntersectionLocal::line) {
987  tt1 = insec->get_point(0)->el1_coord()[0];
988  tt2 = insec->get_point(1)->el1_coord()[0];
989  if (tt1 > tt2) {
990  double swap = tt1;
991  tt1 = tt2;
992  tt2 = swap;
993  }
994 
995  if ((tt2 > (0 - epsilon)) && (tt1 < (1 + epsilon))) {
996  if (tt1 < 0) tt1 = 0;
997  if (tt2 > 1) tt2 = 1;
998  coef = fabs(tt2 - tt1);
999  it = line;
1000  delete insec;
1001  return;
1002  }
1003  }
1004  if (insec->get_type() == IntersectionLocal::point) {
1005  if (cit == 0) {
1006  tt[0] = insec->get_point(0)->el1_coord()[0];
1007  cit++;
1008  } else {
1009  if (IsEqual(tt[0], insec->get_point(0)->el1_coord()[0])) {
1010  delete insec;
1011  insec = nullptr;
1012  continue;
1013  }
1014  if (tt[0] > insec->get_point(0)->el1_coord()[0]) {
1015  tt[1] = tt[0];
1016  tt[0] = insec->get_point(0)->el1_coord()[0];
1017  } else {
1018  tt[1] = insec->get_point(0)->el1_coord()[0];
1019  }
1020  if ((tt[1] > (0 - epsilon)) && (tt[0] < (1 + epsilon))) {
1021  if (tt[0] < 0) tt[0] = 0;
1022  if (tt[1] > 1) tt[1] = 1;
1023  coef = fabs(tt[1] - tt[0]);
1024  it = line;
1025  delete insec;
1026  return;
1027  }
1028  }
1029  }
1030  delete insec;
1031  insec = nullptr;
1032  }
1033  }
1034 
1035  if (insec != nullptr) delete insec;
1036  it = none;
1037 
1038  return;
1039 }
1040 
1041 void GetIntersection(const TTriangle &Tr, const TTetrahedron &Te,
1042  TIntersectionType &it, double &coef) {
1043  //******************************************************************************
1044  //ONE SHOULD ADD TO THIS FUNCTION POINTS WHICH ARE INSIDE TETRAHEDRON
1045  //******************************************************************************
1046  if (!QuickIntersectionTest(Tr, Te)) {
1047  it = none;
1048  return;
1049  }
1050 
1051  TPolygon P;// = new TPolygon();
1052  for (int i = 1; i <= 3; i++) {
1053  if (Te.IsInner(Tr.GetPoint(i))) {
1054  P.Add(Tr.GetPoint(i));
1055  }
1056  }
1057 
1058  if (P.vertexes_count() < 3) {
1059  IntersectionLocal *insec = nullptr;
1060 
1061  for (int i = 1; i <= 3; i++) {
1062  for (int j = 1; j <= 4; j++) {
1063  GetIntersection(Tr.GetAbscissa(i), Te.GetTriangle(j), insec);
1064  if (insec) {
1065  switch (insec->get_type()) {
1066  case IntersectionLocal::point: {
1067  P.Add(Tr.GetAbscissa(i).GetPoint(insec->get_point(0)->el1_coord()[0]));
1068  break;
1069  }
1070  case IntersectionLocal::line: {
1071  P.Add(Tr.GetAbscissa(i).GetPoint(insec->get_point(0)->el1_coord()[0]));
1072  P.Add(Tr.GetAbscissa(i).GetPoint(insec->get_point(1)->el1_coord()[0]));
1073  break;
1074  }
1075  default:
1076  //mythrow((char*) "Runtime error - deny point\n", __LINE__, __FUNC__);
1077  break;
1078  }
1079  delete insec;
1080  }
1081  }
1082  }
1083  for (int i = 1; i <= 6; i++) {
1084  GetIntersection(Te.GetAbscissa(i), Tr, insec);
1085  if (insec) {
1086  switch (insec->get_type()) {
1088  P.Add(Te.GetAbscissa(i).GetPoint(insec->get_point(0)->el1_coord()[0]));
1089  break;
1091  P.Add(Te.GetAbscissa(i).GetPoint(insec->get_point(0)->el1_coord()[0]));
1092  P.Add(Te.GetAbscissa(i).GetPoint(insec->get_point(1)->el1_coord()[0]));
1093  break;
1094  default:
1095  //mythrow((char*) "Runtime error - deny point\n", __LINE__, __FUNC__);
1096  break;
1097  }
1098  delete insec;
1099  }
1100  }
1101  }
1102 
1103  coef = P.GetArea();
1104  if (IsZero(coef)) {
1105  it = none;
1106  } else {
1107  it = area;
1108  }
1109 
1110  return;
1111 }
1112 
1113 //******************************************************************************
1114 // This function returns true if two input geometry can intersect
1115 //******************************************************************************
1116 
1117 template<class A, class B> bool QuickIntersectionTest(const A &a, const B &b) {
1118  for (int i = 1; i <= 3; i++) {
1119  //cout << "i: " << i << "; a: " << a.GetMin(i) << " "<< a.GetMax(i)<< "; b: " << b.GetMin(i) << " "<< b.GetMax(i) << endl;
1120  if (a.GetMin(i) > b.GetMax(i) + epsilon || a.GetMax(i) < b.GetMin(i) - epsilon) {
1121  return false;
1122  }
1123  }
1124  return true;
1125 }
1126 
1127 } // namespace ngh
const TPoint & GetPoint() const
Definition: bisector.cpp:110
enum ngh::Positions TPosition
void Set(int, int, double)
Definition: matrix.cpp:85
const TVector & GetNormal() const
Definition: plain.cpp:115
const TAbscissa & GetAbscissa(int) const
Definition: tetrahedron.cpp:94
bool QuickIntersectionTest(const A &a, const B &b)
bool IsInner(const TPoint &) const
double GetB() const
Definition: plain.cpp:103
#define Length(wts)
Definition: quad.c:31
Definition: abscissa.h:25
double GetC() const
Definition: plain.cpp:107
const IntersectionPoint * get_point(const unsigned int index) const
double Get(int) const
Definition: point.cpp:123
IntersectionType get_type() const
void Add(const TPoint &)
Definition: polygon.cpp:70
void Get(double &, double &, double &) const
Definition: vector.cpp:163
const TPoint & GetPoint(int) const
Definition: triangle.cpp:117
double Y() const
Definition: point.cpp:115
const TVector & GetVector() const
Definition: bisector.cpp:105
double GetA() const
Definition: plain.cpp:99
bool IsZero(double)
Definition: mathfce.cpp:22
const TPoint & GetPoint() const
Definition: plain.cpp:127
TNSolutions
Definition: matrix.h:26
static constexpr bool value
Definition: json.hpp:87
double GetArea()
Definition: polygon.cpp:85
#define OLD_ASSERT(...)
Definition: global_defs.h:131
const std::vector< double > & el2_coord() const
double Length() const
Definition: vector.cpp:159
double Distance(const TBisector &, const TPoint &)
double Get(int)
Definition: matrix.cpp:128
void GetParameter(const TPoint &, double &, bool &) const
Definition: bisector.cpp:124
void add_local_coord(const std::vector< double > &coordin1, const std::vector< double > &coordin2)
bool AreParallel(const TVector &, const TVector &)
Definition: vector.cpp:232
void add_local_point(IntersectionPoint *InPoint)
const TVector & GetV() const
Definition: plain.cpp:123
double X() const
Definition: point.cpp:111
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
void Set(int, double)
Definition: matrix.cpp:137
const std::vector< double > & el1_coord() const
bool Belong(const TPoint &) const
Definition: bisector.cpp:129
void SetPoint(const TPoint &)
Definition: bisector.cpp:92
TNSolutions Gauss(const TMatrix &, TMVector *, const TMVector &)
Definition: matrix.cpp:152
double Z() const
Definition: point.cpp:119
double Dot(const TVector &, const TVector &)
Definition: vector.cpp:199
const double epsilon
Definition: mathfce.h:23
void GetIntersection(const TPlain &P, const TBisector &B, TPosition &pos, double &t)
bool ArePerpendicular(const TVector &, const TVector &)
Definition: vector.cpp:244
enum ngh::Intersections TIntersectionType
const TPlain & GetPlain() const
Definition: triangle.cpp:99
const TAbscissa & GetAbscissa(int) const
Definition: triangle.cpp:103
bool IsInner(const TPoint &) const
Definition: triangle.cpp:205
const TTriangle & GetTriangle(int) const
Definition: tetrahedron.cpp:80
bool Belong(const TPoint &) const
Definition: plain.cpp:137
TVector Cross(const TVector &, const TVector &)
Definition: vector.cpp:182
void SetVector(const TVector &)
Definition: bisector.cpp:96
const TVector & GetU() const
Definition: plain.cpp:119
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
double GetD() const
Definition: plain.cpp:111
int vertexes_count()
Definition: polygon.h:69
IntersectionPoint * interpolate(const IntersectionPoint &A1, const IntersectionPoint &A2, double t)
bool IsEqual(double, double)
Definition: mathfce.cpp:29