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