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