1 /****************************************************************************/
2 /* This file is part of FreeFEM. */
3 /* */
4 /* FreeFEM is free software: you can redistribute it and/or modify */
5 /* it under the terms of the GNU Lesser General Public License as */
6 /* published by the Free Software Foundation, either version 3 of */
7 /* the License, or (at your option) any later version. */
8 /* */
9 /* FreeFEM is distributed in the hope that it will be useful, */
10 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
11 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
12 /* GNU Lesser General Public License for more details. */
13 /* */
14 /* You should have received a copy of the GNU Lesser General Public License */
15 /* along with FreeFEM. If not, see <http://www.gnu.org/licenses/>. */
16 /****************************************************************************/
17 // SUMMARY : ...
18 // LICENSE : LGPLv3
19 // ORG : LJLL Universite Pierre et Marie Curie, Paris, FRANCE
20 // AUTHORS : ...
21 // E-MAIL : ...
22
23 /* clang-format off */
24 //ff-c++-LIBRARY-dep: [flann]
25 //ff-c++-cpp-dep:
26 /* clang-format on */
27
28 #include <ff++.hpp>
29 #include <AFunction_ext.hpp>
30 using namespace Fem2D;
31
32 static bool debug = false;
33
34 class R2close {
35 public:
36 double *data; // minimal points+ delta
37 typedef double *Point;
38 int n, nx, offset;
39 Point *P;
40 const double EPSILON;
41 double x0, y0, x1, y1, coef; // boundin box
R2close()42 R2close( ) : data(0), n(0), nx(1000000), P(new Point[nx]), EPSILON(1e-6), offset(0) {
43 InitialiserListe( );
44 }
45
R2close(double * dd,int mx,double eps=1e-6,int offsett=1)46 R2close(double *dd, int mx, double eps = 1e-6, int offsett = 1)
47 : data(dd), n(0), nx(mx), P(new Point[nx]), EPSILON(eps), offset(offsett) {
48 InitialiserListe( );
49 }
50
operator [](int i) const51 Point operator[](int i) const { return P[i]; }
52
53 int N, m; // le nombre de code = ncase()%n
54 int *head; // pointeur tableau de dimension m contenant les tete de liste pour chaque code
55 int *next; // pointeur tableau de dimension nx donnant le point suivant de mÍme code
56 static const int NotaPoint = -1;
InitialiserListe()57 void InitialiserListe( ) {
58 if (data) {
59 x0 = data[0];
60 y0 = data[1];
61 x1 = data[2];
62 y1 = data[3];
63 } else {
64 x0 = 0;
65 y0 = 1;
66 x1 = 0;
67 y1 = 1;
68 }
69
70 coef = 1. / max(x1 - x0, y1 - y0);
71 if (verbosity > 10) {
72 cout << " bounding box ClosePoints Pmin=[" << x0 << ", " << y0 << "], Pmax=[ " << x1
73 << " " << y1 << "] "
74 << "eps= " << EPSILON << endl;
75 }
76
77 N = max(sqrt(nx), 10.);
78 m = max(nx / 10, 100);
79 next = new int[nx];
80 head = new int[m];
81
82 for (int i = 0; i < m; ++i) {
83 head[i] = NotaPoint;
84 }
85 }
86
AddSimple(double * p)87 int AddSimple(double *p) {
88 double x = p[0], y = p[offset];
89
90 assert(n < nx);
91 P[n] = p;
92
93 int k = ncase(x, y) % m;
94 assert(k >= 0);
95 next[n] = head[k];
96 head[k] = n;
97 if (debug) {
98 cout << " AddSimple " << n << " <- " << k << " / " << x << " " << y << " / " << offset
99 << endl;
100 }
101
102 return n++;
103 }
104
105 private:
Exist(double * p) const106 Point *Exist(double *p) const {
107 double x = p[0], y = p[offset];
108
109 for (int i = 0; i < n; ++i) {
110 if (Equivalent(x, y, P[i])) {
111 return P + i;
112 }
113 }
114
115 return 0;
116 }
117
118 public:
Equivalent(double x0,double y0,const Point Q) const119 bool Equivalent(double x0, double y0, const Point Q) const {
120 return ((x0 - Q[0]) * (x0 - Q[0]) + (y0 - Q[offset]) * (y0 - Q[offset])) < EPSILON * EPSILON;
121 }
122
Exist(double x,double y,int k) const123 Point *Exist(double x, double y, int k) const {
124 for (int i = head[k % m]; i != NotaPoint; i = next[i]) {
125 if (Equivalent(x, y, P[i])) {
126 return P + i;
127 }
128 }
129
130 return 0;
131 }
132
Add(double * p)133 int Add(double *p) {
134 Point *q = Exist(p);
135
136 if (q) {
137 return q - P; // numero du point
138 } else {
139 return AddSimple(p);
140 }
141 }
142
FindAll(double x,double y,int * p)143 int FindAll(double x, double y, int *p) {
144 int nbp = 0;
145
146 if (debug) {
147 cout << " Find " << x << " " << y << " " << EPSILON << " " << ncase(x, y) << ": ";
148 }
149
150 double eps = EPSILON / 2;
151 Point *q = 0;
152 int kk[9] = {}, k = 0, nc;
153
154 for (int i = -1; i < 2; ++i) {
155 for (int j = -1; j < 2; ++j) {
156 nc = ncase(x + eps * i, y + eps * j);
157 if (nc >= 0) {
158 for (int i = 0; i < k; ++i) { // remove double cas e..
159 if (kk[i] == nc) {
160 nc = -1;
161 break;
162 }
163 }
164 }
165
166 if (nc >= 0) {
167 kk[k++] = nc;
168 }
169 }
170 }
171
172 if (k > 4) {
173 cout << " ClosePoints: FindAll Bug ??? : " << k << " : ";
174
175 for (int i = 0; i < k; ++i) {
176 cout << " " << kk[i];
177 }
178
179 cout << endl;
180 }
181
182 assert(k <= 4);
183
184 for (int i = 0; i < k; ++i) {
185 int k = kk[i];
186
187 for (int i = head[k % m]; i != NotaPoint; i = next[i]) {
188 if (Equivalent(x, y, P[i])) {
189 p[nbp++] = i;
190 }
191 }
192 }
193
194 return nbp;
195 }
196
Find(double x,double y)197 Point *Find(double x, double y) {
198 if (debug) {
199 cout << " Find " << x << " " << y << " " << EPSILON << " " << ncase(x, y) << ": ";
200 }
201
202 // double x = p[0], y=p[offset];
203 double eps = EPSILON / 2;
204 Point *q = 0;
205 int kk[9] = {}, k = 0, nc;
206
207 for (int i = -1; i < 2; ++i) {
208 for (int j = -1; j < 2; ++j) {
209 nc = ncase(x + eps * i, y + eps * j);
210 // cout <<x+eps*i << " " << y+eps*j << " " << nc << " . ";
211 if (nc >= 0) {
212 for (int i = 0; i < k; ++i) { // remove double cas e..
213 if (kk[i] == nc) {
214 nc = -1;
215 break;
216 }
217 }
218 }
219
220 if (nc >= 0) {
221 kk[k++] = nc;
222 }
223 }
224 }
225
226 if (k > 4) {
227 cout << " ClosePoints: Bug ??? : " << k << " : ";
228
229 for (int i = 0; i < k; ++i) {
230 cout << " " << kk[i];
231 }
232
233 cout << endl;
234 }
235
236 assert(k <= 4);
237
238 for (int i = 0; i < k; ++i) {
239 q = Exist(x, y, kk[i]);
240 if (debug) {
241 cout << " " << kk[i];
242 }
243
244 if (q) {
245 break;
246 }
247 }
248
249 if (debug) {
250 cout << " q= " << q << endl;
251 }
252
253 if (q) {
254 return q; // numero du point
255 } else {
256 return 0;
257 }
258 }
259
Find(double * p)260 Point *Find(double *p) { return Find(*p, *(p + offset)); }
261
AddOpt(double * p)262 int AddOpt(double *p) {
263 Point *q = Find(p);
264
265 if (q) {
266 return q - P; // numero du point
267 } else {
268 return AddSimple(p);
269 }
270 }
271
ncase(double x,double y)272 long ncase(double x, double y) {
273 if (x < x0 || x >= x1 || y < y0 || y >= y1) {
274 return -1; // dehors
275 } else {
276 return long((x - x0) / EPSILON / 2) +
277 long((y - y0) / EPSILON / 2) * N; // indice de la case contenant (x,y).
278 }
279 }
280
~R2close()281 ~R2close( ) {
282 delete[] P;
283 delete[] head;
284 delete[] next;
285 }
286
287 // pas de copie de cette structure
288
289 private:
290 R2close(const R2close &);
291 R2close &operator=(const R2close &);
292 };
293
dist2(int n,double * p,double * q)294 double dist2(int n, double *p, double *q) {
295 double s = 0;
296
297 for (int i = 0; i < n; ++n) {
298 s += (q[i] - p[i]) * (q[i] - p[i]);
299 }
300
301 return s;
302 }
303
304 long Hcode(int n, double eps, double *p, double *p0);
305
CloseTo2(Stack stack,double const & eps,KNM_<double> const & P,KNM_<double> const & Q)306 KN< long > *CloseTo2(Stack stack, double const &eps, KNM_< double > const &P,
307 KNM_< double > const &Q) {
308 long Pm0 = P.M( );
309 long Qm0 = Q.M( );
310 int po10 = P.step * P.shapei.step;
311 double x0 = P(0, ':').min( );
312 double y0 = P(1, ':').min( );
313 double x1 = P(0, ':').max( );
314 double y1 = P(1, ':').max( );
315
316 // add cc
317 double dd = max(x1 - x0, y1 - y0) * 0.01;
318
319 if (dd == 0) {
320 dd = max(abs(x0), abs(y0)) * 1e-8;
321 }
322
323 if (dd == 0) {
324 dd = 1e-8;
325 }
326
327 double data[] = {x0 - dd, y0 - dd, x1 + dd, y1 + dd};
328 R2close S(data, Pm0, eps, po10);
329
330 for (int i = 0; i < Pm0; ++i) {
331 if (verbosity > 19) {
332 cout << i << " :: " << P(0, i) << " " << P(1, i) << endl;
333 }
334
335 S.AddSimple(&P(0, i));
336 }
337
338 KN< long > *pr = new KN< long >(Qm0);
339
340 for (int i = 0; i < Qm0; ++i) {
341 R2close::Point *p = S.Find(Q(0, i), Q(1, i));
342 if (p) {
343 (*pr)[i] = p - S.P;
344 } else {
345 (*pr)[i] = -1;
346 }
347 }
348
349 return Add2StackOfPtr2FreeRC(stack, pr);
350 }
351
CloseTo2t(Stack stack,double const & eps,KNM_<double> const & P,KNM_<double> const & Q)352 KN< long > *CloseTo2t(Stack stack, double const &eps, KNM_< double > const &P,
353 KNM_< double > const &Q) {
354 return CloseTo2(stack, eps, P.t( ), Q.t( ));
355 }
356
CloseTo(Stack stack,double const & eps,KNM_<double> const & P,KNM<double> * const & q,bool tq,bool inv=0)357 KN< long > *CloseTo(Stack stack, double const &eps, KNM_< double > const &P,
358 KNM< double > *const &q, bool tq, bool inv = 0) {
359 long n0 = P.N( );
360 long m0 = P.M( );
361
362 if (verbosity > 2) {
363 cout << " -ClosePoints Size array; n0 " << n0 << " m0 " << m0 << endl;
364 }
365
366 ffassert(n0 == 2);
367 KNM< double > &Qo = *q;
368 double *p0 = &(P(0, 0));
369 int offset10 = P.step * P.shapei.step;
370 int offset01 = P.step * P.shapej.step;
371 if (verbosity > 10) {
372 cout << " offset of 0 1 : " << offset01 << endl;
373 cout << " offset of 1 0 : " << offset10 << endl;
374 }
375
376 // MeshPoint &mp = *MeshPointStack(stack); // the struct to get x, y, normal, value
377 double x0 = P(0, ':').min( );
378 double y0 = P(1, ':').min( );
379 double x1 = P(0, ':').max( );
380 double y1 = P(1, ':').max( );
381
382 // add cc
383 double dd = max(x1 - x0, y1 - y0) * 0.01;
384 double data[] = {x0 - dd, y0 - dd, x1 + dd, y1 + dd};
385 KN< long > *pr = 0;
386 if (inv) {
387 pr = new KN< long >(m0);
388 }
389
390 R2close S(data, m0, eps, offset10);
391
392 for (int i = 0; i < m0; ++i) {
393 if (verbosity > 19) {
394 cout << i << " :: " << P(0, i) << " " << P(1, i) << endl;
395 }
396
397 int j = S.AddOpt(&P(0, i));
398 if (pr) {
399 (*pr)[i] = j;
400 }
401 }
402
403 if (pr == 0) {
404 pr = new KN< long >(S.n);
405 }
406
407 if (q) {
408 if (tq) {
409 Qo.resize(S.n, n0);
410 } else {
411 Qo.resize(n0, S.n);
412 }
413
414 KNM_< double > Q(tq ? Qo.t( ) : Qo);
415
416 for (int i = 0; i < S.n; ++i) {
417 double *pi = S[i];
418 if (!inv) {
419 (*pr)[i] = (pi - p0) / offset01;
420 }
421
422 for (int j = 0; j < n0; ++j) {
423 Q(j, i) = pi[j * offset10];
424 }
425 }
426 } else if (!inv) {
427 for (int i = 0; i < S.n; ++i) {
428 double *pi = S[i];
429 (*pr)[i] = (pi - p0) / offset01;
430 }
431 }
432
433 if (verbosity > 2) {
434 cout << " - ClosePoint: nb of common points " << m0 - S.n;
435 }
436
437 return Add2StackOfPtr2FreeRC(stack, pr);
438 }
439
440 template< bool inv >
CloseTo(Stack stack,double const & eps,KNM<double> * const & p,KNM<double> * const & q)441 KN< long > *CloseTo(Stack stack, double const &eps, KNM< double > *const &p,
442 KNM< double > *const &q) {
443 KNM_< double > P(*p);
444 return CloseTo(stack, eps, P, q, false, inv);
445 }
446
447 template< bool inv >
CloseTo(Stack stack,double const & eps,Transpose<KNM<double> * > const & p,KNM<double> * const & q)448 KN< long > *CloseTo(Stack stack, double const &eps, Transpose< KNM< double > * > const &p,
449 KNM< double > *const &q) {
450 KNM_< double > P(p.t->t( ));
451 return CloseTo(stack, eps, P, q, false, inv);
452 }
453
454 template< bool inv >
CloseTo(Stack stack,double const & eps,Transpose<KNM<double> * > const & p,Transpose<KNM<double> * > const & q)455 KN< long > *CloseTo(Stack stack, double const &eps, Transpose< KNM< double > * > const &p,
456 Transpose< KNM< double > * > const &q) {
457 KNM_< double > P(p.t->t( ));
458 return CloseTo(stack, eps, P, q, true, inv);
459 }
460
461 template< bool inv >
CloseTo(Stack stack,double const & eps,KNM<double> * const & p,Transpose<KNM<double> * > const & q)462 KN< long > *CloseTo(Stack stack, double const &eps, KNM< double > *const &p,
463 Transpose< KNM< double > * > const &q) {
464 KNM_< double > P(*p);
465 return CloseTo(stack, eps, P, q.t, true, inv);
466 }
467
468 template< bool inv >
CloseTo(Stack stack,double const & eps,KNM<double> * const & p)469 KN< long > *CloseTo(Stack stack, double const &eps, KNM< double > *const &p) {
470 KNM_< double > P(*p);
471 return CloseTo(stack, eps, P, 0, false, inv);
472 }
473
474 template< bool inv >
CloseTo(Stack stack,double const & eps,KNM_<double> const & p)475 KN< long > *CloseTo(Stack stack, double const &eps, KNM_< double > const &p) {
476 KNM_< double > P(p);
477 if (verbosity > 5) {
478 cout << " CloseTo KNM_ " << P.N( ) << " " << P.M( ) << endl;
479 }
480
481 return CloseTo(stack, eps, P, 0, false, inv);
482 }
483
484 template< bool inv >
CloseTo(Stack stack,double const & eps,Transpose<KNM<double> * > const & p)485 KN< long > *CloseTo(Stack stack, double const &eps, Transpose< KNM< double > * > const &p) {
486 KNM_< double > P(p.t->t( ));
487 return CloseTo(stack, eps, P, 0, false, inv);
488 }
489
490 template< bool inv >
CloseTo(Stack stack,double const & eps,pmesh const & pTh,KNM<double> * const & pq)491 KN< long > *CloseTo(Stack stack, double const &eps, pmesh const &pTh, KNM< double > *const &pq) {
492 ffassert(pTh && pq);
493 const Mesh &Th = *pTh;
494 const KNM< double > Q = *pq;
495 int np = Q.N( );
496 assert(Q.M( ) >= 2);
497 KN< long > *pr = new KN< long >(inv ? Th.nv : np);
498 KN< int > b(Th.nv);
499 b = 0;
500
501 for (int i = 0; i < Th.nv; i++) {
502 if (Th(i).lab != 0) {
503 b[i] = 1;
504 }
505 }
506
507 for (int i = 0; i < Th.neb; i++) {
508 b[Th(Th.bedges[i][0])] = 1;
509 b[Th(Th.bedges[i][1])] = 1;
510 }
511
512 *pr = -1L;
513 R2 Pmin, Pmax;
514 Th.BoundingBox(Pmin, Pmax);
515 int nv = b.sum( );
516 if (verbosity > 9) {
517 cout << " Th.nv " << Th.nv << " " << nv << "/ " << Pmin << " " << Pmax << endl;
518 }
519
520 FQuadTree *quadtree =
521 new FQuadTree(pTh, Pmin, Pmax, nv); // put all the old vertices in the quadtree
522
523 // copy of the old vertices
524 for (int i = 0; i < Th.nv; i++) {
525 if (b[i]) {
526 cout << i << " " << Th.vertices[i] << endl;
527 quadtree->Add(Th.vertices[i]);
528 }
529 }
530
531 cout << " After quadterr" << endl;
532
533 for (int j = 0; j < np; ++j) {
534 R2 P(Q(j, 0), Q(j, 1));
535 Vertex *pV = quadtree->ToClose(P, eps);
536 if (pV) {
537 pV = quadtree->NearestVertex(P);
538 long k = Th(pV);
539 if (inv) {
540 (*pr)[k] = j;
541 } else {
542 (*pr)[j] = k;
543 }
544 }
545 }
546
547 delete quadtree;
548
549 return Add2StackOfPtr2FreeRC(stack, pr);
550 }
551
InterAB_Disq(R2 A,R2 B,R2 C,double r)552 bool InterAB_Disq(R2 A, R2 B, R2 C, double r) {
553 double ACB = det(A, C, B);
554
555 if (ACB < 0) {
556 return false;
557 }
558
559 R2 AB(A, B);
560 R2 AC(A, C), CB(C, B);
561 double LAB = AB.norme( );
562 double H = ACB / 4 / LAB;
563 if (H > r) {
564 return false;
565 }
566
567 double AC2 = AC.norme2( ), r2 = r * r;
568 if (AC2 < r2) {
569 return true;
570 }
571
572 double CB2 = CB.norme2( );
573 if (CB2 < r2) {
574 return true;
575 }
576
577 double ABAC = (AB, AC), ABCB = (AB, CB);
578 return ABAC > 0 && ABCB > 0;
579 }
580
Add(KN<long> & I,int i)581 void Add(KN< long > &I, int i) {
582 int n = I.N( );
583 int m = -I[n - 1];
584
585 if (m <= 0) {
586 m = n;
587 n = m * 2;
588 I.resize(n);
589 I[n - 1] = -m - 1;
590 } else {
591 --m;
592 }
593
594 if (debug) {
595 cout << " add " << m << " " << i << " " << n << endl;
596 }
597
598 I[m] = i;
599 if (m < n - 1) {
600 --I[n - 1];
601 }
602 }
603
Clean(KN<long> & I)604 void Clean(KN< long > &I) {
605 int n = I.N( );
606 int m = -I[n - 1];
607
608 if (m >= 0) {
609 I.resize(m - 1);
610 }
611 }
612
Voisinage(KNM_<double> const & P,KNM_<double> const & Q,double const & eps,KN<KN<long>> * const & IJ)613 long Voisinage(KNM_< double > const &P, KNM_< double > const &Q, double const &eps,
614 KN< KN< long > > *const &IJ) {
615 debug = (verbosity > 999);
616 int np = P.N( );
617 int nq = Q.N( );
618 int mp = P.M( );
619 int mq = Q.M( );
620 double *p = &P(0, 0);
621 int offset01 = P.step * P.shapej.step;
622 int offset10 = P.step * P.shapei.step;
623 ffassert(mp == 2);
624 ffassert(mq == 2);
625 KN< int > lp(np);
626
627 IJ->resize(nq);
628
629 for (int i = 0; i < nq; i++) {
630 (*IJ)[i].resize(2);
631 (*IJ)[i][0] = -1;
632 (*IJ)[i][1] = -1;
633 }
634
635 if (verbosity > 99) {
636 cout << " offset01 " << offset01 << " " << offset10 << " p" << p << " " << np << " " << P.M( )
637 << endl;
638 }
639
640 // store - the size in last value of IJ[i][j]
641 double data[4];
642 data[0] = data[2] = p[0];
643 data[1] = data[3] = p[offset01];
644
645 for (int i = 0, k = 0; i < np; ++i, k += offset10) {
646 data[0] = min(data[0], p[k]);
647 data[2] = max(data[2], p[k]);
648 data[1] = min(data[1], p[k + offset01]);
649 data[3] = max(data[3], p[k + offset01]);
650 }
651
652 double eps2 = eps + eps;
653 data[0] -= eps2;
654 data[2] += eps2;
655 data[1] -= eps2;
656 data[3] += eps2;
657
658 R2close SP(data, np, eps, offset01);
659
660 for (int i = 0; i < np; ++i) {
661 SP.AddSimple(&P(i, 0));
662 }
663
664 for (int j = 0; j < nq; ++j) {
665 int nlp = SP.FindAll(Q(j, 0), Q(j, 1), lp);
666
667 for (int k = 0; k < nlp; ++k) {
668 int i = lp[k];
669 if (verbosity > 99) {
670 cout << " Add to i=" << i << " -> j " << j << endl;
671 }
672
673 Add((*IJ)[i], j);
674 }
675 }
676
677 for (int j = 0; j < nq; ++j) {
678 Clean((*IJ)[j]);
679 }
680
681 debug = 0;
682 return 0;
683 }
684
685 #ifdef WITH_flann
686 #include <flann/flann.hpp>
ff_flann_search(KNM_<double> const & P,KNM_<double> const & Q,double const & eps,KN<KN<long>> * const & pIJ)687 long ff_flann_search(KNM_< double > const &P, KNM_< double > const &Q, double const &eps,
688 KN< KN< long > > *const &pIJ) {
689 KN< KN< long > > &IJ = *pIJ;
690 int mp = P.M( );
691 int mq = Q.M( );
692 int np = P.N( );
693 int nq = Q.N( );
694 double *p = &P(0, 0);
695 double *q = &Q(0, 0);
696 int offset01 = P.step * P.shapej.step;
697 int offset10 = P.step * P.shapei.step;
698 int qoffset01 = Q.step * Q.shapej.step;
699 int qoffset10 = Q.step * Q.shapei.step;
700
701 cout << np << " " << nq << " po 01,10 =:" << offset01 << ", " << offset10 << ", " << &P(0, 1) - p
702 << " qo 01,10 =:" << qoffset01 << ", " << qoffset10 << ", " << &Q(0, 1) - q << endl;
703 cout << np << " " << mp << endl;
704 cout << nq << " " << mq << endl;
705
706 ffassert(mp == mq && offset10 == 1 && qoffset10 == 1);
707
708 IJ.resize(nq);
709 flann::Matrix< double > dataset(p, np, mp);
710 flann::Matrix< double > query(q, nq, mq);
711 std::vector< std::vector< int > > indices;
712 std::vector< std::vector< double > > dists;
713 flann::SearchParams params(128);
714 params.eps = eps;
715 flann::Index< flann::L2< double > > index(dataset, flann::KDTreeIndexParams(4));
716 index.buildIndex( );
717 index.radiusSearch(query, indices, dists, eps, params);
718
719 for (int j = 0; j < indices.size( ); ++j) {
720 IJ[j].resize(indices[j].size( ));
721
722 for (int k = 0; k < indices[j].size( ); ++k) {
723 IJ[j][k] = indices[j][k];
724 }
725 }
726
727 for (int j = 0; j < indices.size( ); ++j) {
728 int k = j * mq;
729 cout << j << " [ " << q[k] << " " << q[k + 1] << "] : ";
730
731 for (int k = 0; k < indices[j].size( ); ++k) {
732 cout << indices[j][k] << " " << dists[j][k] << ", ";
733 }
734
735 cout << endl;
736 }
737
738 return 0;
739 }
740
741 #endif
742 using bamg::BinaryRand;
WalkInTriangle(const Mesh & Th,int it,double * lambda,R2 PF)743 int WalkInTriangle(const Mesh &Th, int it, double *lambda, R2 PF) {
744 const Triangle &T(Th[it]);
745 const R2 Q[3] = {(const R2)T[0], (const R2)T[1], (const R2)T[2]};
746
747 R2 P = lambda[0] * Q[0] + lambda[1] * Q[1] + lambda[2] * Q[2];
748
749 // couleur(15);MoveTo( P); LineTo( PF);
750 R l[3];
751 l[0] = Area2(PF, Q[1], Q[2]);
752 l[1] = Area2(Q[0], PF, Q[2]);
753 l[2] = Area2(Q[0], Q[1], PF);
754 R Det = l[0] + l[1] + l[2];
755 l[0] /= Det;
756 l[1] /= Det;
757 l[2] /= Det;
758 const R eps = 1e-5;
759 int neg[3], k = 0;
760 int kk = -1;
761 if (l[0] > -eps && l[1] > -eps && l[2] > -eps) {
762
763 lambda[0] = l[0];
764 lambda[1] = l[1];
765 lambda[2] = l[2];
766 } else {
767
768 if (l[0] < eps && lambda[0] != l[0]) neg[k++] = 0;
769 if (l[1] < eps && lambda[1] != l[1]) neg[k++] = 1;
770 if (l[2] < eps && lambda[2] != l[2]) neg[k++] = 2;
771 R eps1 = T.area * 1.e-5;
772
773 if (k == 2) // 2
774 {
775 // let j be the vertex beetween the 2 edges
776 int j = 3 - neg[0] - neg[1];
777 R S = Area2(P, PF, Q[j]);
778
779 if (S > eps1)
780 kk = (j + 1) % 3;
781 else if (S < -eps1)
782 kk = (j + 2) % 3;
783 else if (BinaryRand( ))
784 kk = (j + 1) % 3;
785 else
786 kk = (j + 2) % 3;
787
788 } else if (k == 1)
789 kk = neg[0];
790 if (kk >= 0) {
791 R d = lambda[kk] - l[kk];
792
793 throwassert(d);
794 R coef = lambda[kk] / d;
795 R coef1 = 1 - coef;
796 lambda[0] = lambda[0] * coef1 + coef * l[0];
797 lambda[1] = lambda[1] * coef1 + coef * l[1];
798 lambda[2] = lambda[2] * coef1 + coef * l[2];
799 lambda[kk] = 0;
800 }
801 }
802 int jj = 0;
803 R lmx = lambda[0];
804 if (lmx < lambda[1]) jj = 1, lmx = lambda[1];
805 if (lmx < lambda[2]) jj = 2, lmx = lambda[2];
806 if (lambda[0] < 0) lambda[jj] += lambda[0], lambda[0] = 0;
807 if (lambda[1] < 0) lambda[jj] += lambda[1], lambda[1] = 0;
808 if (lambda[2] < 0) lambda[jj] += lambda[2], lambda[2] = 0;
809 return kk;
810 }
Cut(const Mesh & Th,R2 P,R2 & PF)811 BoundaryEdge *Cut(const Mesh &Th, R2 P, R2 &PF) {
812 R2 PHat;
813 bool outside;
814 R l[3];
815 const Triangle *K = Th.Find(P, PHat, outside);
816 ffassert(!outside);
817 int it = Th(K);
818 PHat.toBary(l);
819 int k = 0;
820 int j;
821 while ((j = WalkInTriangle(Th, it, l, PF)) >= 0) {
822 int jj = j;
823 throwassert(l[j] == 0);
824 R a = l[(j + 1) % 3], b = l[(j + 2) % 3];
825 int itt = Th.ElementAdj(it, j);
826 if (itt == it || itt < 0) {
827 PF = Th[it](R2(l + 1)); // point de sortie
828 int i1 = Th(it, (jj + 1) % 3), i2 = Th(it, (jj + 2) % 3);
829 return Th.TheBoundaryEdge(i1, i2);
830 }
831 it = itt;
832 l[j] = 0;
833 l[(j + 1) % 3] = b;
834 l[(j + 2) % 3] = a;
835 ffassert(k++ < 1000);
836 }
837 return 0;
838 }
BorderIntersect(pmesh const & pTh,KN_<double> const & IX,KN_<double> const & IY,KN_<double> const & OX,KN_<double> const & OY,KN_<long> const & cL)839 long BorderIntersect(pmesh const &pTh, KN_< double > const &IX, KN_< double > const &IY,
840 KN_< double > const &OX, KN_< double > const &OY, KN_< long > const &cL) {
841 const Mesh &Th = *pTh;
842 KN_< double > ox = OX, oy = OY;
843 KN_< long > L = cL;
844 int np = IX.N( );
845 ffassert(OX.N( ) == np && OY.N( ) == np && IY.N( ) == np && cL.N( ) == np); // verif taille
846 long ni = 0;
847 for (int i = 0; i < np; ++i) {
848 R2 P(IX[i], IY[i]);
849 R2 Q(OX[i], OY[i]);
850 BoundaryEdge *e = Cut(Th, P, Q);
851
852 if (e) {
853 L[i] = e->lab;
854 OX[i] = Q.x, OY[i] = Q.y;
855 } else
856 L[i] = notaregion;
857 }
858 return ni;
859 }
860
init()861 static void init( ) {
862 #ifdef WITH_flann
863 Global.Add("radiusSearch", "(",
864 new OneOperator4_< long, KNM_< double >, KNM_< double >, double, KN< KN< long > > * >(
865 ff_flann_search));
866 #endif
867 Global.Add("Voisinage", "(",
868 new OneOperator4_< long, KNM_< double >, KNM_< double >, double, KN< KN< long > > * >(
869 Voisinage));
870 Global.Add("neighborhood", "(",
871 new OneOperator4_< long, KNM_< double >, KNM_< double >, double, KN< KN< long > > * >(
872 Voisinage));
873 Global.Add("ClosePoints2", "(",
874 new OneOperator3s_< KN< long > *, double, KNM_< double >, KNM_< double > >(CloseTo2));
875
876 // numbering ..
877 Global.Add("ClosePoints", "(",
878 new OneOperator2s_< KN< long > *, double, KNM_< double > >(CloseTo< false >));
879 Global.Add("ClosePoints", "(",
880 new OneOperator3s_< KN< long > *, double, pmesh, KNM< double > * >(CloseTo< false >));
881 Global.Add("ClosePoints1", "(",
882 new OneOperator2s_< KN< long > *, double, KNM_< double > >(CloseTo< true >));
883 Global.Add("ClosePoints1", "(",
884 new OneOperator3s_< KN< long > *, double, pmesh, KNM< double > * >(CloseTo< true >));
885 Global.Add("BorderIntersect", "(",
886 new OneOperator6_< long, pmesh, KN_< double >, KN_< double >, KN_< double >,
887 KN_< double >, KN_< long > >(BorderIntersect));
888 }
889
890 LOADFUNC(init);
891