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 : Signed distance function
18 // LICENSE : LGPLv3
19 // ORG : LJLL Universite Pierre et Marie Curie, Paris, FRANCE
20 // AUTHORS : Frederic Hecht
21 // E-MAIL : frederic.hecht@sorbonne-universite.fr
22
23 /* clang-format off */
24 //ff-c++-LIBRARY-dep:
25 //ff-c++-cpp-dep:
26 /* clang-format on */
27
28 #ifndef WITH_NO_INIT
29 #include "ff++.hpp"
30 #include "AFunction_ext.hpp"
31 #endif
32
33 using namespace std;
34
35 #include <set>
36 #include <vector>
37 #include <map>
38 #include <algorithm>
39 #include <queue>
40 #include <limits>
41 static int debug = 0, debug1 = 0;
42 using namespace Fem2D;
43
44 template< class Rd >
zero(double a,const Rd & A,double b,const Rd & B)45 Rd zero(double a, const Rd &A, double b, const Rd &B) { // zero affine
46 double la = b / (b - a), lb = a / (a - b);
47
48 return la * A + lb * B;
49 }
50
51 template< class Rd >
distmin(const Rd & A,const Rd & B,const Rd & Q)52 double distmin(const Rd &A, const Rd &B, const Rd &Q) {
53 double d = 0;
54 Rd AB(A, B), AQ(A, Q);
55 double ab2 = (AB, AB);
56 double lc = (AQ, AB) / ab2;
57 Rd CQ = AQ - lc * AB; // ( CQ , AB) = 0
58 Rd C = A + lc * AB; // or Q - CQ
59
60 if (lc < 0) {
61 d = Norme2(AQ);
62 } else if (lc > 1.) {
63 d = Norme2(Rd(B, Q));
64 } else {
65 d = Norme2(CQ);
66 }
67
68 if (verbosity > 9999) {
69 cout << " distmin: d =" << d << " /" << lc << " :: " << A << " " << B << " " << Q << " C" << C
70 << endl;
71 }
72
73 return d;
74 }
75
76 template< class Rd >
distmin(const Rd & A,const Rd & B,const Rd & Q,double aq,double bq)77 double distmin(const Rd &A, const Rd &B, const Rd &Q, double aq, double bq) {
78 double d = 0;
79 Rd AB(A, B), AQ(A, Q);
80 double ab2 = (AB, AB);
81 double lp = (AQ, AB) / ab2;
82 Rd PQ = AQ - lp * AB; // ( PQ , AB) = 0
83
84 if (lp < 0) {
85 d = aq;
86 } else if (lp > 1.) {
87 d = bq;
88 } else {
89 d = Norme2(PQ);
90 }
91
92 if (verbosity > 9999) {
93 cout << " distmin:AB/Q: d =" << d << " /" << lp << " :: A " << A << " B " << B << " Q " << Q
94 << " P " << (A + lp * AB) << endl;
95 }
96
97 return d;
98 }
99
100 template< class Rd >
distmin(const Rd & A,double a,const Rd & B,double b,const Rd & Q,double aq,double bq)101 double distmin(const Rd &A, double a, const Rd &B, double b, const Rd &Q, double aq, double bq) {
102 // compule the min_P in [A,B] of p + || PQ ||
103 // where p = affine interpolation of a,b on [A,B]
104 // P = (1-l) A + l B ; p = (1-l) a + l b for l in [0,1]
105 double dmin = min(a + aq, b + bq);
106 double ab = b - a;
107 Rd AB(A, B), AQ(A, Q);
108 double ab2 = (AB, AB);
109 Rd H = ab * AB / ab2; // H = d p/ dP
110 double h2 = (H, H);
111 int cas = 0;
112
113 if (h2 < 1) {
114 cas = 1;
115 double lc = (AQ, AB) / ab2;
116 Rd CQ = AQ - lc * AB; // ( CQ , AB) = 0
117 Rd C = A + lc * AB; // or Q - CQ
118 assert(abs((CQ, AB)) < 1e-6);
119 double r2 = (CQ, CQ) / ab2;
120 double lpm = lc + copysign(sqrt(r2 * h2 / (1 - h2)), -ab);
121 // lpm in [0,1]
122 if (verbosity > 999) {
123 Rd M = A + lpm * AB;
124 cout << " lgm " << lpm << " r= " << sqrt(r2) << " M= " << M << " Q =" << Q
125 << " ::" << a + lpm * ab << " " << ab << endl;
126 }
127
128 lpm = max(0., min(1., lpm));
129
130 if (lpm > 0. && lpm < 1) {
131 cas = 2;
132 Rd M = A + lpm * AB, MQ(M, Q);
133 dmin = a + lpm * ab + sqrt((MQ, MQ));
134 /* // check ???? min ....
135 * lpm += 0.001;
136 * MQ=Rd(A +lpm*AB,Q);
137 * double dmin1 = a + lpm*ab + sqrt((MQ,MQ));
138 * lpm -= 0.002;
139 * MQ=Rd(A +lpm*AB,Q);
140 * double dmin2 = a + lpm*ab + sqrt((MQ,MQ));
141 * if(verbosity>99) cout << " ### "<< dmin1 << " " << dmin << " " << dmin2 << " " << endl;
142 * ffassert(dmin1 >= dmin1&& dmin2 >= dmin );
143 */
144 }
145 }
146
147 if (verbosity > 99) {
148 cout << " distmin/ AaBaQ " << A << " " << a << " / " << B << " " << b << " / " << Q
149 << " / dmin= " << dmin << " cas =" << cas << endl;
150 }
151
152 return dmin;
153 }
154
distmin(const R3 & A,double a,const R3 & B,double b,const R3 & C,double c,const R3 & Q,double aq,double bq,double cq)155 double distmin(const R3 &A, double a, const R3 &B, double b, const R3 &C, double c, const R3 &Q,
156 double aq, double bq, double cq) {
157 // let fA the affine function on ABC
158 const int in = 1;
159 int cas = 0, flat = 0;
160 double dmin = min(min(a + aq, b + bq), c + cq);
161 R3 AB(A, B), AC(A, C), AQ(A, Q);
162 double ab2 = (AB, AB), acab = (AC, AB), ac2 = (AC, AC);
163 double aqab = (AQ, AB), aqac = (AQ, AC);
164 double pdet = ab2 * ac2 - acab * acab;
165 double pb = (aqab * ac2 - aqac * acab) / pdet;
166 double pc = (ab2 * aqac - aqab * acab) / pdet;
167 double pa = 1 - pb - pc;
168 R3 P = pa * A + pb * B + pc * C; // proj of Q on plan ABC.
169 R3 PQ(P, Q);
170 // PG = grad(fA)
171 double fab = b - a, fac = c - a;
172
173 if (abs(fab) + abs(fac) < 1e-16) { // const function fA (flat)
174 flat = 1;
175 if (a >= 0. && b >= 0. && c >= 0.) {
176 dmin = a + Norme2(PQ);
177 cas = in;
178 }
179 } else {
180 // Calcul de d FA / dP =
181 R3 AZ = fab * AC - fac * AB;
182 // fA in constant on line AZ
183 R3 AG = AZ ^ PQ; // otho of AZ in ABC :
184 // AG = gb AB + gc AC ;
185 double agab = (AG, AB), agac = (AG, AC);
186 double gab = (agab * ac2 - agac * acab) / pdet;
187 double gac = (ab2 * agac - agab * acab) / pdet;
188 R3 AGG = gab * AB + gac * AC;
189 ffassert(Norme2(AGG - AG) < 1e-6); // verif ..
190 double fag = fab * gab + fac * gac;
191 R3 H = AG / fag;
192 double r2 = (PQ, PQ);
193 double h2 = (H, H);
194 double lpm = -sqrt(r2 * h2 / (1 - h2));
195 double gb = gab / fag;
196 double gc = gac / fag;
197 double ga = -gb - gc;
198 double am = pa + lpm * ga, bm = pb + lpm * gb, cm = pc + lpm * gc;
199 if (am >= 0 && bm >= 0 && cm > 0.) {
200 R3 M = am * A + bm * B + cm * C;
201 R3 QM(Q, M);
202 dmin = am * a + bm * b + cm * c + Norme2(QM);
203 cas = in;
204 if (debug1) {
205 // verif .
206 {
207 double df1 = ga * a + gb * b + gc * c;
208 cout << "df1 " << df1 << " " << fag << endl;
209
210 ffassert(abs(df1 - 1) < 1e-6);
211 ffassert(abs((H, PQ)) < 1e-6);
212 ffassert(abs((AZ, PQ)) < 1e-6);
213 ffassert(abs((AG, PQ)) < 1e-6);
214 }
215 {
216 am += 0.001, bm -= 0.001;
217 R3 M = am * A + bm * B + cm * C, QM(Q, M);
218 double dmin1 = am * a + bm * b + cm * c + Norme2(QM);
219 ffassert(dmin <= dmin1);
220 }
221 {
222 bm += 0.001, cm -= 0.001;
223 R3 M = am * A + bm * B + cm * C, QM(Q, M);
224 double dmin2 = am * a + bm * b + cm * c + Norme2(QM);
225 ffassert(dmin <= dmin2);
226 }
227 }
228 }
229 }
230
231 if (cas != in) {
232 if (flat) { // externe => test 3 bord ...
233 double dminab = a + distmin(A, B, Q, aq, bq);
234 double dminac = a + distmin(A, C, Q, aq, cq);
235 double dminbc = a + distmin(B, C, Q, bq, cq);
236 dmin = min(min(dminab, dminac), min(dminbc, dmin));
237 } else { // externe => test 3 bord ...
238 double dminab = distmin(A, a, B, b, Q, aq, bq);
239 double dminac = distmin(A, a, C, c, Q, aq, cq);
240 double dminbc = distmin(B, b, C, c, Q, bq, cq);
241 dmin = min(min(dminab, dminac), min(dminbc, dmin));
242 }
243 }
244
245 if (debug) {
246 cout << " AaBbCc/q " << dmin << " " << cas << flat << endl;
247 }
248
249 return dmin;
250 }
251
distmin(const R3 & A,double a,const R3 & B,double b,const R3 & C,double c,const R3 & Q)252 double distmin(const R3 &A, double a, const R3 &B, double b, const R3 &C, double c, const R3 &Q) {
253 R3 AQ(A, Q), BQ(B, Q), CQ(C, Q);
254 double aq = Norme2(AQ), bq = Norme2(BQ), cq = Norme2(CQ);
255
256 return distmin(A, a, B, b, C, c, Q, aq, bq, cq);
257 }
258
259 template< class Rd >
distmin(const Rd & A,double a,const Rd & B,double b,const Rd & Q)260 double distmin(const Rd &A, double a, const Rd &B, double b, const Rd &Q) {
261 double aq = Norme2(Rd(A, Q));
262 double bq = Norme2(Rd(B, Q));
263
264 return distmin(A, a, B, b, Q, aq, bq);
265 }
266
distmin(const R3 & A,const R3 & B,const R3 & C,const R3 & Q)267 double distmin(const R3 &A, const R3 &B, const R3 &C, const R3 &Q) {
268 double d = 0;
269 R3 AB(A, B), AC(A, C), AQ(A, Q);
270 double ab2 = (AB, AB), acab = (AC, AB), ac2 = (AC, AC);
271 double aqab = (AQ, AB), aqac = (AQ, AC);
272 double det = ab2 * ac2 - acab * acab;
273 double b = (aqab * ac2 - aqac * acab) / det;
274 double c = (ab2 * aqac - aqab * acab) / det;
275 double a = 1 - b - c;
276 R3 P = a * A + b * B + c * C; // proj of Q on plan ABC.
277
278 if (debug) {
279 cout << " distmin ABC/q " << a << " " << b << " " << c << endl;
280 }
281
282 R3 PQ(P, Q);
283 assert(abs((PQ, AB)) < 1e-7);
284 assert(abs((PQ, AC)) < 1e-7);
285 if (a >= 0. && b >= 0. && c >= 0.) {
286 d = Norme2(PQ);
287 } else {
288 double d1 = distmin(A, B, Q);
289 double d2 = distmin(B, C, Q);
290 double d3 = distmin(C, A, Q);
291 d = min(min(d1, d2), d3);
292 }
293
294 return d;
295 }
296
297 typedef Mesh3::Element Tet;
298
DistanceIso0(const Tet & K,double * f,double * fK)299 int DistanceIso0(const Tet &K, double *f, double *fK) {
300 int ret = 1;
301 // fk value f
302 const double eps = 1e-16;
303 R3 P[10];
304 int np = 0;
305
306 if (abs(f[0]) < eps) {
307 f[0] = 0.;
308 P[np++] = K[0];
309 }
310
311 if (abs(f[1]) < eps) {
312 f[1] = 0.;
313 P[np++] = K[1];
314 }
315
316 if (abs(f[2]) < eps) {
317 f[2] = 0.;
318 P[np++] = K[2];
319 }
320
321 if (abs(f[3]) < eps) {
322 f[3] = 0.;
323 P[np++] = K[3];
324 }
325
326 for (int e = 0; e < 6; ++e) {
327 int i1 = Tet::nvedge[e][0];
328 int i2 = Tet::nvedge[e][1];
329
330 if ((f[i1] < 0. && f[i2] > 0.) or (f[i1] > 0. && f[i2] < 0.)) {
331 P[np++] = zero< R3 >(f[i1], K[i1], f[i2], K[i2]);
332 }
333 }
334
335 if (np && debug) {
336 cout << " np " << np << " " << P[0] << " " << P[1] << " :: " << f[0] << " " << f[1] << " "
337 << f[2] << " " << f[3] << endl;
338 }
339
340 if (np == 0) {
341 ret = 0;
342 } else if (np == 1) {
343 for (int j = 0; j < 4; ++j) {
344 fK[j] = Norme2(R3(P[0], K[j]));
345 }
346 } else if (np == 2) {
347 for (int j = 0; j < 4; ++j) {
348 fK[j] = distmin(P[0], P[1], (R3)K[j]);
349 }
350 } else if (np == 3 || np == 4) {
351 for (int j = 0; j < 4; ++j) {
352 fK[j] = distmin(P[0], P[1], P[2], (R3)K[j]);
353 }
354 } else {
355 fK[0] = fK[1] = fK[2] = fK[3] = 0.; //
356 }
357
358 if (debug) {
359 cout << ret << " 3d DistanceIso0 " << np << " " << fK[0] << " " << fK[1] << fK[2] << " "
360 << fK[3] << endl;
361 }
362
363 return ret;
364 }
365
DistanceIso0(const Triangle & K,double * f,double * fK)366 int DistanceIso0(const Triangle &K, double *f, double *fK) {
367 // fk value f
368
369 const double eps = 1e-16;
370 R2 P[6];
371 int np = 0;
372
373 if (abs(f[0]) < eps) {
374 f[0] = 0.;
375 }
376
377 if (abs(f[1]) < eps) {
378 f[1] = 0.;
379 }
380
381 if (abs(f[2]) < eps) {
382 f[2] = 0.;
383 }
384
385 int ke[6];
386
387 for (int e = 0; e < 3; ++e) {
388 int i1 = (e + 1) % 3;
389 int i2 = (e + 2) % 3;
390 if (f[i1] == 0) {
391 ke[np] = i1, P[np++] = K[i1];
392 } else if ((f[i1] < 0. && f[i2] > 0.) or (f[i1] > 0. && f[i2] < 0.)) {
393 ke[np] = e;
394 P[np++] = zero< R2 >(f[i1], K[i1], f[i2], K[i2]);
395 }
396 }
397
398 if (np && debug) {
399 cout << " np " << np << " " << P[0] << " " << P[1] << " :: " << f[0] << " " << f[1] << " "
400 << f[2] << endl;
401 }
402
403 if (np == 0) {
404 return 0;
405 } else if (np == 1) {
406 for (int j = 0; j < 3; ++j) {
407 fK[j] = Norme2(R2(P[0], K[j]));
408 }
409 } else if (np == 2) {
410 for (int j = 0; j < 3; ++j) {
411 fK[j] = distmin(P[0], P[1], (R2)K[j]);
412 }
413 } else {
414 fK[0] = fK[1] = fK[2] = 0.; //
415 }
416
417 if (debug) {
418 cout << np << " DistanceIso0 np="
419 << " " << fK[0] << " " << fK[1] << " " << fK[2] << endl;
420 }
421
422 return np;
423 }
424
CheckDist(double a,double b)425 double CheckDist(double a, double b) {
426 for (int i = 0; i < 30; ++i) {
427 double a = 1, b = 1.1, c = 1.5;
428 R3 A(-0.5, 0.001, 0.002), B(0.5, -0.001, 0.0001), C(0.0001, 1., -0.0003), Q(i * 0.1, 0.001, 1.);
429 double d = distmin(A, a, B, b, C, c, Q);
430 cout << " d = " << i << " == " << d << endl;
431 }
432
433 return 0;
434 }
435
436 // fonction determinant les points d'intersection
437
438 class Distance2d_Op : public E_F0mps {
439 public:
440 Expression eTh, eff, exx;
441 static const int n_name_param = 1; //
442 static basicAC_F0::name_and_type name_param[];
443 Expression nargs[n_name_param];
444
arg(int i,Stack stack,double a) const445 double arg(int i, Stack stack, double a) const {
446 return nargs[i] ? GetAny< double >((*nargs[i])(stack)) : a;
447 }
448
arg(int i,Stack stack,long a) const449 long arg(int i, Stack stack, long a) const {
450 return nargs[i] ? GetAny< long >((*nargs[i])(stack)) : a;
451 }
452
arg(int i,Stack stack,KN<long> * a) const453 KN< long > *arg(int i, Stack stack, KN< long > *a) const {
454 return nargs[i] ? GetAny< KN< long > * >((*nargs[i])(stack)) : a;
455 }
456
arg(int i,Stack stack,string * a) const457 string *arg(int i, Stack stack, string *a) const {
458 return nargs[i] ? GetAny< string * >((*nargs[i])(stack)) : a;
459 }
460
461 public:
Distance2d_Op(const basicAC_F0 & args,Expression tth,Expression fff,Expression xxx)462 Distance2d_Op(const basicAC_F0 &args, Expression tth, Expression fff, Expression xxx)
463 : eTh(tth), eff(fff), exx(xxx) {
464 args.SetNameParam(n_name_param, name_param, nargs);
465 }
466
467 AnyType operator( )(Stack stack) const;
468 };
469
470 basicAC_F0::name_and_type Distance2d_Op::name_param[] = {{"distmax", &typeid(double)}};
471
Add(const Mesh & Th,int kk,int ee,double * fv)472 pair< double, long > Add(const Mesh &Th, int kk, int ee, double *fv) {
473 const Triangle &K = Th[kk];
474 int a = (ee + 1) % 3;
475 int b = (ee + 2) % 3;
476 int q = ee;
477 int ia = Th(kk, a), ib = Th(kk, b), iq = Th(kk, q);
478 double fq = distmin< R2 >(K[a], fv[ia], K[b], fv[ib], K[q]);
479
480 if (debug) {
481 cout << iq << " ** add " << kk << " " << ee << " ; " << fq << " :: " << fv[ia] << " " << fv[ib]
482 << " || " << fv[iq] << endl;
483 }
484
485 return pair< double, long >(fq, kk * 3 + ee);
486 }
487
Add(const Mesh3 & Th,int kk,int ee,double * fv)488 pair< double, long > Add(const Mesh3 &Th, int kk, int ee, double *fv) {
489 typedef Mesh3::Element Tet;
490 const Tet &K = Th[kk];
491 int a = Tet::nvface[ee][0];
492 int b = Tet::nvface[ee][1];
493 int c = Tet::nvface[ee][2];
494 int q = ee;
495 int ia = Th(kk, a), ib = Th(kk, b), ic = Th(kk, c), iq = Th(kk, q);
496 double fq = distmin(K[a], fv[ia], K[b], fv[ib], K[c], fv[ic], K[q]);
497 if (debug) {
498 cout << " ** add " << kk << " " << ee << " ; " << fq << " :: " << fv[ia] << " " << fv[ib] << " "
499 << fv[ic] << " || " << fv[iq] << endl;
500 }
501
502 return pair< double, long >(fq, kk * 4 + ee);
503 }
504
DistanceIso0(const Mesh & Th,int k,double * f,double * fv)505 int DistanceIso0(const Mesh &Th, int k, double *f, double *fv) {
506 typedef Mesh::Element Elem;
507 const int nbve = 3;
508 const Elem &K = Th[k];
509 int iK[nbve] = {Th(k, 0), Th(k, 1), Th(k, 2)};
510 R fk[nbve] = {f[iK[0]], f[iK[1]], f[iK[2]]};
511 // cut here ..
512 double FK[nbve] = {fv[iK[0]], fv[iK[1]], fv[iK[2]]};
513 int cas = DistanceIso0(K, fk, FK);
514 if (cas > 1) { // OK iso cut triangle
515 //
516
517 fv[iK[0]] = min(fv[iK[0]], FK[0]);
518 fv[iK[1]] = min(fv[iK[1]], FK[1]);
519 fv[iK[2]] = min(fv[iK[2]], FK[2]);
520 if (debug) {
521 cout << " DistanceIso0 set K" << cas << " " << iK[0] << " " << iK[1] << " " << iK[2] << " "
522 << fv[iK[0]] << " " << fv[iK[1]] << " " << fv[iK[2]] << endl;
523 }
524 }
525
526 return cas > 1;
527 }
528
DistanceIso0(const Mesh3 & Th,int k,double * f,double * fv)529 int DistanceIso0(const Mesh3 &Th, int k, double *f, double *fv) {
530 typedef Mesh3::Element Elem;
531 const int nbve = 4;
532 const Elem &K = Th[k];
533 int iK[nbve] = {Th(k, 0), Th(k, 1), Th(k, 2), Th(k, 3)};
534 R fk[nbve] = {f[iK[0]], f[iK[1]], f[iK[2]], f[iK[3]]};
535 // cut here ..
536 double FK[nbve] = {fv[iK[0]], fv[iK[1]], fv[iK[2]], fv[iK[3]]};
537 int cas = DistanceIso0(K, fk, FK);
538 if (cas > 0) { // OK iso cut triangle
539 fv[iK[0]] = min(fv[iK[0]], FK[0]);
540 fv[iK[1]] = min(fv[iK[1]], FK[1]);
541 fv[iK[2]] = min(fv[iK[2]], FK[2]);
542 fv[iK[3]] = min(fv[iK[3]], FK[3]);
543 }
544
545 return cas;
546 }
547
548 template< class Mesh >
Distance(Stack stack,const Mesh * pTh,Expression eff,KN<double> * pxx,double dmax)549 AnyType Distance(Stack stack, const Mesh *pTh, Expression eff, KN< double > *pxx, double dmax) {
550 typedef typename Mesh::Element Elem;
551 const int nbve = Mesh::Rd::d + 1;
552 debug = 0;
553 if (verbosity > 99) {
554 debug = 1;
555 }
556
557 double unset = -std::numeric_limits< double >::max( );
558 double distinf = std::numeric_limits< double >::max( );
559 MeshPoint *mp(MeshPointStack(stack)), mps = *mp;
560 double isovalue = 0.;
561 ffassert(pTh);
562 const Mesh &Th(*pTh);
563 long nbv = Th.nv; // nombre de sommet
564 long nbt = Th.nt; // nombre de triangles
565 // long nbe=Th.neb; // nombre d'aretes fontiere
566 typedef KN< double > Rn;
567 typedef KN< long > Zn;
568 typedef pair< double, long > KEY; // Distance max , vertex i/triangle k= 3*k+i
569 if (verbosity > 2) {
570 cout << " distance max = " << dmax << " nt =" << nbt << endl;
571 }
572
573 Rn tK(nbt), tV(nbv), f(nbv);
574 pxx->resize(nbv);
575 Rn &fv = *pxx;
576 Zn mK(nbt);
577 mK = 0L;
578 f = unset;
579 fv = distinf;
580 typedef std::priority_queue< KEY, std::vector< KEY >, std::greater< KEY > > myPQ;
581 myPQ pqs;
582 vector< long > markT(Th.nt, 1);
583
584 for (int it = 0; it < Th.nt; ++it) {
585 for (int iv = 0; iv < nbve; ++iv) {
586 int i = Th(it, iv);
587 if (f[i] == unset) {
588 mp->setP(pTh, it, iv);
589 f[i] = GetAny< double >((*eff)(stack)) - isovalue;
590 }
591 }
592 }
593
594 long err = 0, nt0 = 0;
595
596 for (long k = 0; k < Th.nt; ++k) {
597 int cas = DistanceIso0(Th, k, f, fv);
598 if (cas) {
599 ++nt0;
600 markT[k] = 0;
601 } // init
602 }
603
604 for (long k = 0; k < Th.nt; ++k) {
605 if (markT[k] == 0) {
606 for (int e = 0; e < nbve; ++e) {
607 int ee = e, kk = Th.ElementAdj(k, ee);
608 if (kk >= 0 && (markT[kk] != 0)) {
609 pqs.push(Add(Th, kk, ee, fv));
610 }
611 }
612 }
613 }
614
615 //
616 if (verbosity > 3) {
617 cout << " Distance: nb elemets in queue after init" << pqs.size( ) << " / nb set " << nt0
618 << endl;
619 }
620
621 double fm;
622
623 while (!pqs.empty( )) {
624 KEY t = pqs.top( );
625 pqs.pop( );
626 fm = t.first;
627 if (fm > dmax) {
628 break;
629 }
630
631 int e = t.second % nbve;
632 int k = t.second / nbve;
633 int iq = Th(k, e);
634
635 if (debug) {
636 cout << iq << " " << fm << " -- k=" << k << " e=" << e << " fv:" << fv[iq] << " / "
637 << markT[k] << endl;
638 }
639
640 if (markT[k] != 0) { // already done, skeep
641 markT[k] = 0;
642 fv[iq] = min(fm, fv[iq]);
643
644 for (int e = 0; e < nbve; ++e) {
645 int ee = e, kk = Th.ElementAdj(k, ee);
646 if (kk >= 0 && (markT[kk] != 0)) {
647 pqs.push(Add(Th, kk, ee, fv));
648 }
649 }
650 }
651 }
652
653 // clean
654 err = 0;
655
656 for (int i = 0; i < nbv; ++i) {
657 if (fv[i] == distinf) {
658 fv[i] = dmax;
659 err++;
660 }
661
662 fv[i] = copysign(fv[i], f[i]);
663 }
664
665 if (err && fm < dmax) {
666 if (verbosity) {
667 cout << " Distance 2d : we miss some point " << err << endl;
668 }
669 }
670
671 return err;
672 }
673
operator ( )(Stack stack) const674 AnyType Distance2d_Op::operator( )(Stack stack) const {
675 double distinf = std::numeric_limits< double >::max( );
676 double dmax = arg(0, stack, distinf);
677
678 KN< double > *pxx = 0;
679 pxx = GetAny< KN< double > * >((*exx)(stack));
680 const Mesh *pTh = GetAny< const Mesh * >((*eTh)(stack));
681
682 return Distance< Mesh >(stack, pTh, eff, pxx, dmax);
683 }
684
685 class Distance3d_Op : public E_F0mps {
686 public:
687 Expression eTh, eff, exx;
688 static const int n_name_param = 1;
689 static basicAC_F0::name_and_type name_param[];
690 Expression nargs[n_name_param];
691
arg(int i,Stack stack,double a) const692 double arg(int i, Stack stack, double a) const {
693 return nargs[i] ? GetAny< double >((*nargs[i])(stack)) : a;
694 }
695
arg(int i,Stack stack,long a) const696 long arg(int i, Stack stack, long a) const {
697 return nargs[i] ? GetAny< long >((*nargs[i])(stack)) : a;
698 }
699
arg(int i,Stack stack,KN<long> * a) const700 KN< long > *arg(int i, Stack stack, KN< long > *a) const {
701 return nargs[i] ? GetAny< KN< long > * >((*nargs[i])(stack)) : a;
702 }
703
arg(int i,Stack stack,string * a) const704 string *arg(int i, Stack stack, string *a) const {
705 return nargs[i] ? GetAny< string * >((*nargs[i])(stack)) : a;
706 }
707
708 public:
Distance3d_Op(const basicAC_F0 & args,Expression tth,Expression fff,Expression xxx)709 Distance3d_Op(const basicAC_F0 &args, Expression tth, Expression fff, Expression xxx)
710 : eTh(tth), eff(fff), exx(xxx) {
711 args.SetNameParam(n_name_param, name_param, nargs);
712 }
713
714 AnyType operator( )(Stack stack) const;
715 };
716
717 basicAC_F0::name_and_type Distance3d_Op::name_param[] = {{"distmax", &typeid(double)}};
max(double a,double b,double c)718 inline double max(double a, double b, double c) { return max(max(a, b), c); }
719
min(double a,double b,double c)720 inline double min(double a, double b, double c) { return min(min(a, b), c); }
721
operator ( )(Stack stack) const722 AnyType Distance3d_Op::operator( )(Stack stack) const {
723 double distinf = std::numeric_limits< double >::max( );
724 double dmax = arg(0, stack, distinf);
725
726 KN< double > *pxx = 0;
727 pxx = GetAny< KN< double > * >((*exx)(stack));
728 const Mesh3 *pTh = GetAny< const Mesh3 * >((*eTh)(stack));
729
730 return Distance< Mesh3 >(stack, pTh, eff, pxx, dmax);
731 }
732
733 class Distance2d_P1 : public OneOperator {
734 public:
735 typedef const Mesh *pmesh;
Distance2d_P1()736 Distance2d_P1( )
737 : OneOperator(atype< long >( ), atype< pmesh >( ), atype< double >( ),
738 atype< KN< double > * >( )) {}
739
code(const basicAC_F0 & args) const740 E_F0 *code(const basicAC_F0 &args) const {
741 return new Distance2d_Op(args, t[0]->CastTo(args[0]), t[1]->CastTo(args[1]),
742 t[2]->CastTo(args[2]));
743 }
744 };
745
746 class Distance3d_P1 : public OneOperator {
747 public:
748 typedef const Mesh3 *pmesh3;
Distance3d_P1()749 Distance3d_P1( )
750 : OneOperator(atype< long >( ), atype< pmesh3 >( ), atype< double >( ),
751 atype< KN< double > * >( )) {}
752
code(const basicAC_F0 & args) const753 E_F0 *code(const basicAC_F0 &args) const {
754 return new Distance3d_Op(args, t[0]->CastTo(args[0]), t[1]->CastTo(args[1]),
755 t[2]->CastTo(args[2]));
756 }
757 };
758
finit()759 static void finit( ) {
760 typedef const Mesh *pmesh;
761 typedef const Mesh3 *pmesh3;
762
763 Global.Add("distance", "(", new Distance2d_P1);
764 Global.Add("distance", "(", new Distance3d_P1);
765 Global.Add("checkdist", "(", new OneOperator2< double, double, double >(CheckDist));
766 }
767
768 LOADFUNC(finit); // une variable globale qui serat construite au chargement dynamique
769