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