1 // -*- Mode : c++ -*-
2 //
3 // -*- Mode : c++ -*-
4 //
5 // SUMMARY : mshptg mesh generator ( form mshptg.f / inria / emc2 software)
6 // USAGE : LGPL
7 // ORG : LJLL Universite Pierre et Marie Curi, Paris, FRANCE
8 // AUTHOR : Frederic Hecht
9 // E-MAIL : frederic.hecht@ann.jussieu.fr
10 //
11 //
12
13 /*
14
15 This file is part of Freefem++
16
17 Freefem++ is free software; you can redistribute it and/or modify
18 it under the terms of the GNU Lesser General Public License as published by
19 the Free Software Foundation; either version 2.1 of the License, or
20 (at your option) any later version.
21
22 Freefem++ is distributed in the hope that it will be useful,
23 but WITHOUT ANY WARRANTY; without even the implied warranty of
24 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 GNU Lesser General Public License for more details.
26
27 You should have received a copy of the GNU Lesser General Public License
28 along with Freefem++; if not, write to the Free Software
29 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
30 */
31
32
33 // bof bof FH je ne sais pas pourquoi nnnBAMG_LONG_LONG
34 // nov 2005 j'ai corrige les problemes
35 // il y avait un gros bug de le trie
36 // overflow long car on calculais x*x + y*y
37 // pour que cela tienne dans un long long
38 // x*x + y*y < 2^63-1 => x et y < 2^31-1
39 // ---
40 #include <iostream>
41 using namespace std;
42
43 #ifdef BAMG_LONG_LONG
44 #define LONG8 long long
45 #define DLONG8LONG8 1.e17
46 #define MAXCOOR ((double) 1073741823)// 2^30-1=1073741823. )
47 #else
48 #define LONG8 long
49 #define DLONG8LONG8 1e9
50 #define MAXCOOR ((double) 32767. )
51 #endif
52 #define DET8(a11,a12,a21,a22) ( (LONG8) (a11) * (LONG8) (a22) - (LONG8) (a21)* (LONG8) (a12))
53 #include <cmath>
54 #include <cstdio>
55 #include <cstdlib>
56 namespace Fem2D {
57 #define i_nint(x) ((long) x >= 0 ? x + 0.5 : x - 0.5)
58 #define amin(a, b) ((a) < (b) ? (a) : (b))
59 #define amax(a, b) ((a) < (b) ? (b) : (a))
60 #define aabs(a) ((a) < 0 ? -(a) : (a))
61
62 static long nothing = -1073741824L;
63
64
65
66 /*****************************************************************************
67 I. Fabrique une triangulation (de nbsmax sommets maxi.) a partir de:
68 --------------------------------------------------------------------
69 -cr: tableau de taille 2*nbs (abscisse & ordonnees des points en entree)
70 -arete: tableau de taille 2*nba (depart & arrivee des aretes)
71 (decrire le contour dans le sens direct: domaine a
72 gauche des aretes exterieures)
73 -h: tableau de taille nbs (precision relative aux points en entree)
74
75
76 II. Alloue et affecte une structure t de type triangulation:
77 ------------------------------------------------------------
78 -t.np: nombre de sommets en sortie
79 -t.nt: nombre de triangles en sortie
80 -t.rp[i].x:} abscisse
81 t.rp[i].y:} & ordonnee du ieme point du maillage
82 -t.tr[i][j]: numero du jeme point du ieme triangle du maillage
83
84
85 III. Renvoie le type d'erreur:
86 ------------------------------
87 - 0 : pas d'erreur
88 -(-1): pas assez de memoire
89 - >0 : erreur mshptg8_
90 mshptg erreurs
91 1: mshptg := 'le nb de point est < 3 or > que le nb max de point';
92 2: mshptg := 'il y des points confondus ';
93 3: mshptg := 'tout les points sont align»s ';
94 7: mshptg := 'bug dans l''algorithme pour retrouver la frontiŸre (Internal Bug mshfrt)';
95 8: mshptg := 'Une arte a forcer traverse trop de triangles (Internal Bug : mshfr1)';
96 9: mshptg := 'la frontiŸre est crois»e ';
97 10: mshptg := 'il y a un point qui est sur une arte frontiŸre';
98 11: mshptg := 'Un sous domaine est reference par aucun triangle ';
99 20: mshptg := 'mshopt: 3 points confondus (Internal Bug) ';
100 21: mshptg := 'la pile de mshopt est trop petit (Internal Bug)';
101
102
103
104 Remarque: penser a liberer les pointeurs t.rp et t.tr!!!
105 --------- un tableau de taille N est indexe de 0 a N-1
106 les numeros affectes aux points et aux triangles commencent a 0
107 ******************************************************************************/
108 /*
109 long
110 MakeTriangulation (triangulation * t, long nbs, long nbsmax, long nba,
111 double *crbdy, double *hbdy, long *arete, int *ngbdy, long *sd, long nbsd, int* flag, int fflag)
112 {
113 int i, j;
114 long err = 0, nbt, nbsold = nbs;
115 long *c = NULL;
116 long *tri = NULL;
117 long *nu = NULL;
118 long *reft = NULL;
119 int *ngg = NULL;
120 double *cr = NULL;
121 double *h = NULL;
122
123 nbt = 2 * nbsmax;
124 nu = new long[6*nbt];
125 c = new long[2*nbsmax];
126 ngg = new int[nbsmax];
127 tri = new long[(4 * nbsmax + 2 * nbsd)];
128 reft = new long[nbt];
129 cr = new double[(2 * nbsmax + 2)];
130 h = new double[nbsmax];
131
132 for (i = 0; i < 2 * nba; i++)
133 arete[i]++;
134 for (i = 0; i < nbs; i++)
135 {
136 ngg[i] = ngbdy[i];
137 cr[2 * i] = crbdy[2 * i];
138 cr[2 * i + 1] = crbdy[2 * i + 1];
139 h[i] = hbdy[i];
140 }
141 for (i = nbs; i < nbsmax; i++)
142 ngg[i] = 0;
143
144 mshptg8_ (cr, h, c, nu, &nbs, nbsmax, tri, arete, nba, (long *) sd, nbsd, reft, &nbt, .25, .75, &err);
145 for (i = 0; i < 2 * nba; i++)
146 arete[i]--;
147 if (err)
148 goto L1;
149 if (*flag)
150 {
151 delete [] t->rp;t->rp = NULL;
152 delete [] t->tr;t->tr = NULL;
153 delete [] t->ng;t->ng = NULL;
154 delete [] t->ngt;t->ngt = NULL;
155 }
156 t->rp =new rpoint[nbs];
157 t->tr = new triangle[nbt];
158 t->ng = new int[nbs];
159 t->ngt = new int[nbt];
160
161 *flag = 1;
162 t->np = nbs;
163 t->nt = nbt;
164 for (i = 0; i < nbt; i++)
165 {
166 for (j = 0; j < 3; j++)
167 t->tr[i][j] = nu[3 * i + j] - 1;
168 t->ngt[i] = reft[i] - 1;
169 }
170 for (i = 0; i < nbs; i++)
171 {
172 t->rp[i].x = cr[2 * i];
173 t->rp[i].y = cr[2 * i + 1];
174 if (i < nbsold)
175 t->ng[i] = ngg[i];
176 else
177 t->ng[i] = 0;
178 }
179 renum (t);
180 if(!fflag) removeBdyT (t);
181 L1:
182 delete [] nu;nu = NULL;
183 delete [] cr;cr = NULL;
184 delete [] c;c = NULL;
185 delete [] tri;tri = NULL;
186 delete [] reft;reft = NULL;
187 delete [] ngg;ngg = NULL;
188 delete [] h;h = NULL;
189 return err;
190 }
191
192
193
194 int
195 verifietr (double *cr, int ns)
196 {
197 int i;
198 double xh, diameter = 1.F;
199
200 if (ns == 0)
201 return -1;
202 if (ns > 1)
203 {
204 diameter = 0.F;
205 for (i = 0; i < ns; i++)
206 diameter = amax (diameter, aabs (cr[2 * i] - cr[0]) + aabs (cr[2 * i + 1] - cr[1]));
207 }
208 else
209 diameter = 0.001F;
210 for (i = 0; i < ns; i++)
211 {
212 xh = aabs (cr[2 * i] - cr[2 * ns]) + aabs (cr[2 * i + 1] - cr[2 * ns + 1]);
213 if (xh < 1e-5 * diameter)
214 return i;
215 }
216 return -1;
217 }
218
219 */
220
221 int mshrgl_ (double *c, long *nrfs, long *nbs, long *nu, long *w1,
222 long *w, double omega, long itermx, double eps);
223 int mshopt_ (long *c, long *nu, long *t, long a, long *err);
224 void mshvoi_ (long *nu, long *w1, long *w, long *nbt, long *nbs);
225 int msha1p_ (long *t, long *s, long *c, long *nu, long *reft, long *tete, long *nbt,
226 long *err);
227 int mshtri_ (double *cr, long *c, long *nbs, long *tri, LONG8 *nu, double *trfri, long *err);
228 int mshcxi_ (long *c, long *nu, long *tri, long *nbs, long *tete, long *err);
229 int mshfrt_ (long *c, long *nu, long *nbs, long *arete, long nba, long *sd,
230 long nbsd, long *reft, long *w, long *err);
231 int mshgpt_ (long *c, double *cr, long *nu, double *h, long *reft, long *nbs,
232 long nbsmx, long *nbt, double coef, double puis, double *trfri, long *err);
233 long mshlcl_ (long *c, long *nu, long *tete, long *s);
234 int mshtr1_ (LONG8 *criter, long *record, long *n);
235 int mshcvx_ (long direct, long *c, long *nu, long *pfold, long *err);
236 int mshfr1_ (long *c, long *nu, long *it1, long *ita, long *is1, long *s2, long *err);
237 int mshfr2_ (long *c, long *nu, long *lst, long *nbac, long *t, long *ta,
238 long *ss1, long *ss2, long *err);
239
240 int
mshptg8_(double * cr,double * h,long * c,long * nu,long * nbs,long nbsmx,long * tri,long * arete,long nba,long * sd,long nbsd,long * reft,long * nbt,double coef,double puis,long * err)241 mshptg8_ (double *cr, double *h, long *c, long *nu, long *nbs, long nbsmx, long *tri,
242 long *arete, long nba, long *sd,
243 long nbsd, long *reft, long *nbt, double coef, double puis, long *err)
244 {
245 /* System generated locals */
246 long i_1;
247
248 /* Local variables */
249 static long tete, i, j, k, t;
250 static double trfri[4];
251 static long nbsgrn, nbtgrn;
252
253 /* ----------------------------------------------------------------------- */
254 /* but: construire une triangulation a partir d'un ensemble de */
255 /* points et d'un maillage frontalier */
256 /* ----------------------------------------------------------------------- */
257 /* entre : */
258 /* ------- */
259 /* cr(2,nbsmx) tableau des coordonnees des nbs points donnes */
260
261
262 /* h (nbsmx) tableau du h local voulu autour de chaque point */
263 /* donnes */
264 /* nbs nombre de points donnes */
265 /* nbsmx nombre de points maximal a cree */
266 /* si nbs = nbsmx ont ne cree pas de points */
267 /* si nbs < nbsmx => erreur */
268 /* arete(2,nba) tableau des aretes du maillage a forcer */
269 /* exemple :la frontiere */
270 /* nba le nombre d'aretes du maillage */
271 /* sd(2,nbsd) tableau definisant les nbsd sous domaine */
272 /* (reference des triangles gerener) */
273 /* abs(sd(1,i)) = numero d'une l'arete */
274 /* si sd(1,i) est positive alors le sous domaine */
275 /* est a gauche de l'arete sinon il est a droite */
276 /* sd(2,i) donne le numero du sous doimaine */
277
278 /* puis coefficient de generation des points */
279 /* .1 => on propage plus loin les rafinement */
280 /* donnes par h */
281 /* .25 => valeur conseillee */
282 /* coef coefficient sur le test arret */
283 /* le valeur conseillee est .75 */
284 /* remarque le nombre de point genere est en */
285 /* O(coef**2) */
286
287 /* tableaux de travail: */
288 /* -------------------- */
289 /* c(2,nbsmx) tableau d'entiers (copie de coordonnees) */
290 /* tri(ltri) tableau d'entiers */
291 /* out : */
292 /* ----- */
293 /* nbs nombre de points donnes + generes */
294 /* nbt nombre de triangles generes */
295 /* cr(1:2,nbs) coordonnees des sommets donnes + generes */
296 /* nu(1:3,nbt) sommets des triangles (tableau des connections) */
297 /* telle que les sommets tourne dans le sens direct */
298 /* reft(1:nbt) numero de sous domaine de chaque triangle */
299 /* err si err = 0 alors pas de probleme */
300 /* sinon nbt = 0 et pas de triangulation */
301 /* dimension des tableaux */
302 /* ---------------------- */
303 /* definition des parameters */
304 /* nbtmx = 2*(nbs-1) , ltri = max(4*nbs+2*nbsd,nba) */
305 /* long : nu(6*nbtmx) , reft(nbtmx) , c(2*nbsmx) , tri(ltri) */
306 /* long : arete(2,nba), sd(2,nbsd) */
307 /* double : cr(2*nbsmx) , h(nbsmx) */
308
309 /* ---------------------------------------------------------------------- */
310 /* programmeur F.Hecht, France */
311 /* version 1.0 mars 1986 */
312 /* ----------------------------------------------------------------------- */
313
314 /* Parameter adjustments */
315 --reft;
316 sd -= 3;
317 arete -= 3;
318 --tri;
319 --nu;
320 c -= 3;
321 --h;
322 cr -= 3;
323
324 /* Function Body */
325 *err = 0;
326 *nbt = 0;
327 if (*nbs < 3 || nbsmx < *nbs)
328 {
329 *err = 1;
330 printf("mshptg bug in number of points %ld > %ld == max nb points \n",*nbs,nbsmx);
331 return 0;
332 }
333 /* ------------------------- */
334 /* preparation des donnees */
335 /* ------------------------- */
336 mshtri_ (&cr[3], &c[3], nbs, &tri[1], (LONG8*) (void *) &tri[*nbs + 1], trfri, err);
337 if (*err != 0)
338 {
339 return 0;
340 }
341 /* -------------------------------- */
342 /* maillage de l enveloppe convexe */
343 /* -------------------------------- */
344 mshcxi_ (&c[3], &nu[1], &tri[1], nbs, &tete, err);
345 /* ----------------------------------------------------------------------- */
346 /* definition de tableau nu(1:6,2*nbs-2) */
347 /* ----------------------------------------------------------------------- */
348 /* nu(*,ie) definit soit un element ,soit un sommet frontiere */
349 /* si nu(5:6,ie) = (0,0) alors ie est un sommet frontiere */
350 /* avec nu(1,ie) = numero du sommet */
351 /* nu(2,ie) = 8*t + a */
352 /* ou t est le numero du triangle ayant l'arete */
353 /* frontiere (a) dont le premier sommet est nu(1,ie) */
354 /* nu(3,ie) = pointeur dans nu sur sommet frontiere precedent */
355 /* nu(4,ie) = pointeur dans nu sur sommet frontiere suivant */
356 /* sinon ie est un element : */
357 /* nu(1:3,ie) numero des 3 sommets du triangle ie tournant dans */
358 /* le sens direct */
359 /* nu(4:6,ie) = (d4,d5,d6) donnee des 3 aretes ai */
360 /* ai est forme des sommets nu(i-3,ie),nu(mod(i,3)+1,ie) */
361 /* si di < 0 alors arete i est frontiere et -di est pointeur */
362 /* sur 1er sommet frontiere de i */
363 /* sinon arete est interne et di = 8*ta + ata */
364 /* ou ta est le numero du triangle adjacent a l'arete */
365 /* et ata est le numero de l'arete dans ta */
366 /* ------------------------------------------------------------------------ */
367 if (*err != 0)
368 {
369 return 0;
370 }
371
372 i_1 = *nbs;
373 for (i = 1; i <= i_1; ++i)
374 {
375 tri[i] = 0;
376 }
377 i = tete;
378 L20:
379 j = nu[(i - 1) * 6 + 4];
380 tri[nu[(i - 1) * 6 + 1]] = nu[(j - 1) * 6 + 1];
381 i = j;
382 if (i != tete)
383 {
384 goto L20;
385 }
386 /* ----------------------------- */
387 /* traitement frontiere */
388 /* ----------------------------- */
389 k = 0;
390 mshfrt_ (&c[3], &nu[1], nbs, &arete[3], nba, &sd[3], nbsd, &reft[1], &tri[
391 1], err);
392 if (*err != 0)
393 {
394 return 0;
395 }
396 /* ------------------------------------------------------------------- */
397 /* on a modifie nu les sommets frontiere n'ont plus de sens */
398 /* ainsi que les pointeurs sur ces elements */
399 /* ------------------------------------------------------------------- */
400 nbsgrn = *nbs;
401 mshgpt_ (&c[3], &cr[3], &nu[1], &h[1], &reft[1], &nbsgrn, nbsmx, &nbtgrn,
402 coef, puis, trfri, err);
403 if (*err != 0)
404 {
405 return 0;
406 }
407 /* construction du tableau nu(1:3,1:nbt) */
408 /* ------------------------------------------ */
409 *nbt = 0;
410 k = 0;
411 j = 1;
412 i_1 = nbtgrn;
413 for (t = 1; t <= i_1; ++t)
414 {
415 if (nu[j + 5] != 0)
416 {
417 ++(*nbt);
418 reft[*nbt] = reft[t];
419 for (i = 0; i <= 2; ++i)
420 {
421 ++k;
422 nu[k] = nu[j + i];
423 }
424 }
425 j += 6;
426 }
427 /* dans nu il y a (s1(t),s2(t),s3(t),t=1,nbt) */
428 /* ou s1 s2 s3 sont les 3 sommets de t */
429 /* ------------------------------------------------ */
430 i_1 = *nbs;
431 for (i = 1; i <= i_1; ++i)
432 {
433 tri[i] = 1;
434 }
435 i_1 = nbsgrn;
436 for (i = *nbs + 1; i <= i_1; ++i)
437 {
438 tri[i] = 0;
439 }
440 mshvoi_ (&nu[1], &tri[nbsgrn + 1], &nu[*nbt * 3 + 1], nbt, &nbsgrn);
441 mshrgl_ (&cr[3], &tri[1], &nbsgrn, &nu[1], &tri[nbsgrn + 1], &nu[*nbt * 3
442 + 1], 1.4F, 20L, .005F);
443 *nbs = nbsgrn;
444 return 1;
445 } /* mshptg_ */
446
447 /* ********************************************************************** */
448 void
mshvoi_(long * nu,long * w1,long * w,long * nbt,long * nbs)449 mshvoi_ (long *nu, long *w1, long *w, long *nbt, long *nbs)
450 {
451 /* System generated locals */
452 long i_1;
453
454 /* Local variables */
455 static long i, is;
456
457
458 /* RECHERCHE DU VOISINAGE */
459 /* ----------------------- */
460 /* Parameter adjustments */
461 --w;
462 --nu;
463
464 /* Function Body */
465 i_1 = *nbs;
466 for (i = 1; i <= i_1; ++i)
467 {
468 w1[i] = 0;
469 }
470 i_1 = *nbt * 3;
471 for (i = 1; i <= i_1; ++i)
472 {
473 ++w1[nu[i]];
474 }
475 w1[0] = 0;
476 i_1 = *nbs;
477 for (i = 1; i <= i_1; ++i)
478 {
479 w1[i] = w1[i - 1] + w1[i];
480 }
481 i_1 = *nbt * 3;
482 for (i = 1; i <= i_1; ++i)
483 {
484 is = nu[i] - 1;
485 ++w1[is];
486 w[w1[is]] = i;
487 }
488 for (i = *nbs; i >= 1; --i)
489 {
490 w1[i] = w1[i - 1];
491 }
492 w1[0] = 0;
493 } /* mshvoi_ */
494
495 /* ********************************************************************** */
496 int
mshrgl_(double * c,long * nrfs,long * nbs,long * nu,long * w1,long * w,double omega,long itermx,double eps)497 mshrgl_ (double *c, long *nrfs, long *nbs, long *nu, long *w1,
498 long *w, double omega, long itermx, double eps)
499 {
500 /* System generated locals */
501 long i_1, i_2, i_3;
502 double r_1, r_2;
503
504 /* Local variables */
505 static double depx, depy;
506 static long iter;
507 static double xmin, ymin, xmax, ymax;
508 static long i, k, i1, i2, ic;
509 static double bx, by, dx;
510 static long is;
511 static double err;
512
513
514 /* REGULARISATION PAR MOYENNE BARYCENTRIQUE */
515 /* ----------------------------------------- */
516 /* Parameter adjustments */
517 --w;
518 --nu;
519 --nrfs;
520 c -= 3;
521
522 /* Function Body */
523 xmin = c[3];
524 ymin = c[4];
525 xmax = c[3];
526 ymax = c[4];
527 i_1 = *nbs;
528 for (ic = 1; ic <= i_1; ++ic)
529 {
530 /* Computing MIN */
531 r_1 = c[(ic << 1) + 1];
532 xmin = amin (r_1, xmin);
533 /* Computing MIN */
534 r_1 = c[(ic << 1) + 2];
535 ymin = amin (r_1, ymin);
536 /* Computing MAX */
537 r_1 = c[(ic << 1) + 1];
538 xmax = amax (r_1, xmax);
539 /* Computing MAX */
540 r_1 = c[(ic << 1) + 2];
541 ymax = amax (r_1, ymax);
542 }
543 /* Computing MAX */
544 r_1 = xmax - xmin, r_2 = ymax - ymin;
545 dx = amax (r_1, r_2);
546 i_1 = itermx;
547 for (iter = 1; iter <= i_1; ++iter)
548 {
549 err = (double) 0.;
550 i2 = w1[0];
551 i_2 = *nbs;
552 for (is = 1; is <= i_2; ++is)
553 {
554 i1 = i2 + 1;
555 i2 = w1[is];
556 if (i2 >= i1 && nrfs[is] == 0)
557 {
558 bx = (double) 0.;
559 by = (double) 0.;
560 i_3 = i2;
561 for (i = i1; i <= i_3; ++i)
562 {
563 if (w[i] % 3 == 0)
564 {
565 k = w[i] - 2;
566 }
567 else
568 {
569 k = w[i] + 1;
570 }
571 bx += c[(nu[k] << 1) + 1];
572 by += c[(nu[k] << 1) + 2];
573 }
574 bx /= i2 - i1 + 1;
575 by /= i2 - i1 + 1;
576 depx = omega * (c[(is << 1) + 1] - bx);
577 depy = omega * (c[(is << 1) + 2] - by);
578 c[(is << 1) + 1] -= depx;
579 c[(is << 1) + 2] -= depy;
580 /* Computing MAX */
581 r_1 = err, r_2 = aabs (depx), r_1 = amax (r_1, r_2), r_2 = aabs (depy);
582 err = amax (r_1, r_2);
583 }
584 }
585 /* -------------------------------- */
586 if (err <= eps * dx)
587 {
588 return 0;
589 }
590 }
591 return 1;
592 } /* mshrgl_ */
593
594 /* ********************************************************************** */
595 int
mshgpt_(long * c,double * cr,long * nu,double * h,long * reft,long * nbs,long nbsmx,long * nbt,double coef,double puis,double * trfri,long * err)596 mshgpt_ (long *c, double *cr, long *nu, double *h, long *reft, long *nbs,
597 long nbsmx, long *nbt, double coef, double puis, double *trfri, long *err)
598 {
599 /* System generated locals */
600 long i_1;
601 double r_1, r_2, r_3;
602 double d_1, d_2, d_3, d_4, d_5, d_6, d_7, d_8;
603
604 /* Local variables */
605 static double aire;
606 static long tete;
607 static long t;
608 static double x, y;
609 static long itera;
610 static double h1, h2, h3;
611 static long s1, s2, s3;
612 static double hs;
613 static long ix, iy, nbsold;
614 static double det, pui;
615
616 /* Parameter adjustments */
617 --trfri;
618 --reft;
619 --h;
620 nu -= 7;
621 cr -= 3;
622 c -= 3;
623
624 /* Function Body */
625 pui = puis;
626 *nbt = (*nbs << 1) - 2;
627 if (*nbs >= nbsmx)
628 { /* ADD FH avril 2007 */
629 i_1 = *nbt;
630 for (t = 1; t <= i_1; ++t)
631 if (nu[t * 6 + 6] != 0)
632 {
633 mshopt_ (&c[3], &nu[7], &t, 4L, err);
634 mshopt_ (&c[3], &nu[7], &t, 5L, err);
635 mshopt_ (&c[3], &nu[7], &t, 6L, err);
636 }
637 /* FIN ADD FH vril 2007 */
638 return 0;
639 }
640 tete = 0;
641 /* initialisation de la liste des triangles libre */
642 i_1 = *nbt;
643 for (t = 1; t <= i_1; ++t)
644 {
645 if (nu[t * 6 + 6] == 0)
646 {
647 nu[t * 6 + 1] = tete;
648 tete = t;
649 }
650 }
651 itera = 0;
652 L20:
653 ++itera;
654 nbsold = *nbs;
655 i_1 = *nbt;
656 for (t = 1; t <= i_1; ++t)
657 {
658 if (nu[t * 6 + 6] != 0)
659 {
660 s1 = nu[t * 6 + 1];
661 s2 = nu[t * 6 + 2];
662 s3 = nu[t * 6 + 3];
663 /* calcul de 2 fois l'aire du triangle */
664 det = (cr[(s2 << 1) + 1] - cr[(s1 << 1) + 1])
665 * (cr[(s3 << 1) + 2] - cr[(s1 << 1) + 2])
666 - (cr[(s2 << 1) + 2] - cr[(s1 << 1) + 2])
667 * (cr[(s3 << 1) + 1] - cr[(s1 << 1) + 1]);
668 aire = det * coef;
669 if (puis > (double) 0.)
670 {
671 d_2 = (double) h[s1];
672 d_3 = (double) pui;
673 d_4 = (double) h[s2];
674 d_5 = (double) pui;
675 d_6 = (double) h[s3];
676 d_7 = (double) pui;
677 d_1 = (pow (d_2, d_3) + pow (d_4, d_5) + pow (d_6, d_7)) / (double) 3.;
678 d_8 = (double) (1.F / pui);
679 hs = (double)pow (d_1, d_8);
680 }
681 else if (puis > (double) -1.)
682 {
683 d_1 = (double) (h[s1] * h[s2] * h[s3]);
684 hs = (double)pow (d_1, 1. / 3);
685 }
686 else if (puis > (double) -2.)
687 {
688 hs = h[s1] * (double) 3. *h[s2] * h[s3] / (h[s1] * h[s2] + h[
689 s1] * h[s3] + h[s2] * h[s3]);
690 }
691 else
692 {
693 /* Computing 2nd power */
694 r_1 = (double)sqrt (h[s1] * h[s2]);
695 /* Computing 2nd power */
696 r_2 = h[s1] * h[s3];
697 /* Computing 2nd power */
698 r_3 = h[s2] * h[s3];
699 hs = (double)sqrt (3.0) * (h[s1] * h[s2] * h[s3]
700 / (r_1 * r_1) + r_2 * r_2 + r_3 * r_3);
701 }
702 if (aire > hs * hs)
703 {
704 h1 = (double) 1.;
705 h2 = (double) 1.;
706 h3 = (double) 1.;
707 x = (cr[(s1 << 1) + 1] * h1 + cr[(s2 << 1) + 1] * h2 + cr[(s3
708 << 1) + 1] * h3) / (h1 + h2 + h3);
709 y = (cr[(s1 << 1) + 2] * h1 + cr[(s2 << 1) + 2] * h2 + cr[(s3
710 << 1) + 2] * h3) / (h1 + h2 + h3);
711 r_1 = trfri[1] * (x - trfri[2]);
712 ix = (long) i_nint (r_1);
713 r_1 = trfri[1] * (y - trfri[3]) - trfri[4];
714 iy = (long) i_nint (r_1);
715 if ((c[(s2 << 1) + 1] - ix) * (c[(s3 << 1) + 2] - iy) - (c[(
716 s2 << 1) + 2] - iy) * (c[(s3 << 1) + 1] - ix) <= 0 ||
717 (ix - c[(s1 << 1) + 1]) * (c[(s3 << 1) + 2] - c[(s1 <<
718 1) + 2]) - (iy - c[(s1 << 1) + 2]) * (c[(s3 << 1) +
719 1] - c[(s1 << 1) + 1]) <= 0 || (c[(s2 << 1) + 1] - c[(
720 s1 << 1) + 1]) * (iy - c[(s1 << 1) + 2]) - (c[(s2 <<
721 1) + 2] - c[(s1 << 1) + 2]) * (ix - c[(s1 << 1) + 1])
722 <= 0)
723 {
724 }
725 else
726 {
727 if (*nbs >= nbsmx)
728 {
729 return 0;
730 }
731 ++(*nbs);
732 c[(*nbs << 1) + 1] = ix;
733 c[(*nbs << 1) + 2] = iy;
734 cr[(*nbs << 1) + 1] = ix / trfri[1] + trfri[2];
735 cr[(*nbs << 1) + 2] = (iy + trfri[4]) / trfri[1] + trfri[
736 3];
737 h[*nbs] = hs;
738 msha1p_ (&t, nbs, &c[3], &nu[7], &reft[1], &tete, nbt, err);
739 if (*err != 0)
740 {
741 return 0;
742 }
743 }
744 }
745 }
746 }
747 if (nbsold != *nbs)
748 {
749 goto L20;
750 }
751 return 1;
752 } /* mshgpt_ */
753
754 /* ********************************************************************** */
755 int
msha1p_(long * t,long * s,long * c,long * nu,long * reft,long * tete,long * nbt,long * err)756 msha1p_ (long *t, long *s, long *c, long *nu, long *reft, long *tete, long *nbt,
757 long *err)
758 {
759 static long t1, t2, t3, ia2, ia3;
760 static long ta2, ta3;
761
762 // extern int mshopt_();
763 static long tta;
764
765 /* Parameter adjustments */
766 --reft;
767 nu -= 7;
768 c -= 3;
769
770 /* Function Body */
771 t1 = *t;
772 if (*tete == 0)
773 {
774 ++(*nbt);
775 t2 = *nbt;
776 }
777 else
778 {
779 t2 = *tete;
780 *tete = nu[*tete * 6 + 1];
781 }
782 if (*tete == 0)
783 {
784 ++(*nbt);
785 t3 = *nbt;
786 }
787 else
788 {
789 t3 = *tete;
790 *tete = nu[*tete * 6 + 1];
791 }
792 nu[t2 * 6 + 1] = *s;
793 nu[t2 * 6 + 2] = nu[*t * 6 + 2];
794 nu[t2 * 6 + 3] = nu[*t * 6 + 3];
795 nu[t2 * 6 + 4] = (t1 << 3) + 5;
796 nu[t2 * 6 + 5] = nu[*t * 6 + 5];
797 nu[t2 * 6 + 6] = (t3 << 3) + 5;
798 nu[t3 * 6 + 1] = nu[*t * 6 + 1];
799 nu[t3 * 6 + 2] = *s;
800 nu[t3 * 6 + 3] = nu[*t * 6 + 3];
801 nu[t3 * 6 + 4] = (t1 << 3) + 6;
802 nu[t3 * 6 + 5] = (t2 << 3) + 6;
803 nu[t3 * 6 + 6] = nu[*t * 6 + 6];
804 tta = nu[*t * 6 + 5];
805 if (tta > 0)
806 {
807 ta2 = tta / 8;
808 ia2 = tta - (ta2 << 3);
809 nu[ia2 + ta2 * 6] = (t2 << 3) + 5;
810 }
811 tta = nu[*t * 6 + 6];
812 if (tta > 0)
813 {
814 ta3 = tta / 8;
815 ia3 = tta - (ta3 << 3);
816 nu[ia3 + ta3 * 6] = (t3 << 3) + 6;
817 }
818 nu[t1 * 6 + 3] = *s;
819 nu[t1 * 6 + 5] = (t2 << 3) + 4;
820 nu[t1 * 6 + 6] = (t3 << 3) + 4;
821 reft[t2] = reft[*t];
822 reft[t3] = reft[*t];
823 mshopt_ (&c[3], &nu[7], &t1, 4L, err);
824 if (*err != 0)
825 {
826 return 0;
827 }
828 mshopt_ (&c[3], &nu[7], &t2, 5L, err);
829 if (*err != 0)
830 {
831 return 0;
832 }
833 mshopt_ (&c[3], &nu[7], &t3, 6L, err);
834 if (*err != 0)
835 {
836 return 0;
837 }
838 return 1;
839 } /* msha1p_ */
840
841 /* ********************************************************************** */
842 long
mshlcl_(long * c,long * nu,long * tete,long * s)843 mshlcl_ (long *c, long *nu, long *tete, long *s)
844 {
845 /* System generated locals */
846 long ret_val;
847
848 /* Local variables */
849 static long init;
850 static long x, y, pt, ppt;
851 LONG8 det;
852
853 /* Parameter adjustments */
854 nu -= 7;
855 c -= 3;
856
857 /* Function Body */
858 x = c[(*s << 1) + 1];
859 y = c[(*s << 1) + 2];
860 init = 1;
861 pt = *tete;
862 L10:
863 ppt = pt;
864 pt = nu[pt * 6 + 4];
865 if (pt != *tete)
866 {
867 det = (LONG8) x * (LONG8) c[(nu[pt * 6 + 1] << 1) + 2]
868 - (LONG8) y * (LONG8) c[(nu[pt * 6 + 1] << 1) + 1];
869 if (det < 0)
870 {
871 init = 0;
872 goto L10;
873 }
874 else if (init && det == 0)
875 {
876 goto L10;
877 }
878 }
879 ret_val = ppt;
880 return ret_val;
881 } /* mshlcl_ */
882
883 /* ********************************************************************** */
884 int
mshtri_(double * cr,long * c,long * nbs,long * tri,LONG8 * nu,double * trfri,long * err)885 mshtri_ (double *cr, long *c, long *nbs, long *tri, LONG8 *nu, double *trfri, long *err)
886 {
887 /* System generated locals */
888 long i_1, i_2, i_3;
889 double r_1;
890
891 /* Local variables */
892 static long ierr, trik;
893 static double xmin, ymin, xmax, ymax;
894 static long i, j, k, ic, jc;
895
896 // extern int mshtr1_();
897 static long ip, xx;
898 static double aa1, aa2;
899 static long iii, tri3;
900 LONG8 det;
901
902 /* Parameter adjustments */
903 --trfri;
904 --nu;
905 --tri;
906 c -= 3;
907 cr -= 3;
908
909 /* Function Body */
910 ierr = 0;
911 iii = 1;
912 xmin = cr[3];
913 ymin = cr[4];
914 xmax = cr[3];
915 ymax = cr[4];
916 i_1 = *nbs;
917 for (ic = 1; ic <= i_1; ++ic)
918 {
919 /* Computing MIN */
920 r_1 = cr[(ic << 1) + 1];
921 xmin = amin (r_1, xmin);
922 /* Computing MIN */
923 r_1 = cr[(ic << 1) + 2];
924 ymin = amin (r_1, ymin);
925 /* Computing MAX */
926 r_1 = cr[(ic << 1) + 1];
927 xmax = amax (r_1, xmax);
928 /* Computing MAX */
929 r_1 = cr[(ic << 1) + 2];
930 ymax = amax (r_1, ymax);
931 tri[ic] = ic;
932 if (cr[(ic << 1) + 1] < cr[(iii << 1) + 1])
933 {
934 iii = ic;
935 }
936 }
937 aa1 = MAXCOOR / (xmax - xmin);
938 aa2 = MAXCOOR / (ymax - ymin);
939 aa1 = amin (aa1, aa2);
940 aa2 = aa1 * (cr[(iii << 1) + 2] - ymin);
941 trfri[1] = aa1;
942 trfri[2] = cr[(iii << 1) + 1];
943 trfri[3] = ymin;
944 trfri[4] = aa2;
945 i_1 = *nbs;
946 for (ic = 1; ic <= i_1; ++ic)
947 {
948 r_1 = aa1 * (cr[(ic << 1) + 1] - cr[(iii << 1) + 1]);
949 c[(ic << 1) + 1] = (long) i_nint (r_1);
950 r_1 = aa1 * (cr[(ic << 1) + 2] - ymin) - aa2;
951 c[(ic << 1) + 2] = (long) i_nint (r_1);
952 /* Computing 2nd power */
953 i_2 = c[(ic << 1) + 1];
954 /* Computing 2nd power */
955 i_3 = c[(ic << 1) + 2];
956 nu[ic] = (LONG8) i_2 * (LONG8) i_2 + (LONG8) i_3 * (LONG8) i_3;
957 }
958 /* ---------------------------------------------------------- */
959 mshtr1_ (&nu[1], &tri[1], nbs);
960 ip = 1;
961 xx = nu[ip];
962 i_1 = *nbs;
963 for (jc = 1; jc <= i_1; ++jc)
964 {
965 if (nu[jc] > xx)
966 {
967 i_2 = jc - ip;
968 mshtr1_ (&nu[ip], &tri[ip], &i_2);
969 i_2 = jc - 2;
970 for (i = ip; i <= i_2; ++i)
971 {
972 if (nu[i] == nu[i + 1])
973 {
974 ++ierr;
975 if (ierr <10)
976 printf(" The points %ld and %ld are too close \n",tri[i],tri[i+1]);
977 }
978 }
979 xx = nu[jc];
980 ip = jc;
981 }
982 ic = tri[jc];
983 nu[jc] = c[(ic << 1) + 2];
984 }
985 i_1 = *nbs - ip;
986 mshtr1_ (&nu[ip], &tri[ip], &i_1);
987 i_1 = jc - 2;
988 for (i = ip; i <= i_1; ++i)
989 {
990 if (nu[i] == nu[i + 1])
991 {
992 if (ierr <10)
993 printf(" The points %ld and %ld are to close \n",tri[i],tri[i+1]);
994 ++ierr;
995 }
996 }
997 if (ierr != 0)
998 {
999 *err = 2;
1000 printf("mshptg bug 2 \n");
1001
1002 return 0;
1003 }
1004 k = 2;
1005 L50:
1006 if (k <= *nbs)
1007 {
1008 ++k;
1009 det = (LONG8) c[(tri[2] << 1) + 1] * (LONG8) c[(tri[k] << 1) + 2]
1010 - (LONG8) c[(tri[2] << 1) + 2] * (LONG8) c[(tri[k] << 1) + 1];
1011 if (det == 0)
1012 {
1013 goto L50;
1014 }
1015 }
1016 else
1017 {
1018 *err = 3;
1019 return 0;
1020 }
1021 /* k est le premier point non aligne */
1022 trik = tri[k];
1023 for (j = k - 1; j >= 3; --j)
1024 {
1025 tri[j + 1] = tri[j];
1026 }
1027 tri[3] = trik;
1028 if (det < 0)
1029 {
1030 /* on inverse les points 2 3 tries */
1031 tri3 = tri[3];
1032 tri[3] = tri[2];
1033 tri[2] = tri3;
1034 }
1035 return 1;
1036 } /* mshtri_ */
1037
1038 /* ********************************************************************** */
1039 int
mshtr1_(LONG8 * criter,long * record,long * n)1040 mshtr1_ (LONG8 *criter, long *record, long *n)
1041 {
1042 /* System generated locals */
1043 long i_1;
1044
1045 /* Local variables */
1046 long i, j, l, r, rec;
1047 LONG8 crit;
1048 /* trie selon les valeurs de criter croissantes */
1049 /* record suit le reordonnancement */
1050
1051
1052 /* Parameter adjustments */
1053 --record;
1054 --criter;
1055
1056 /* Function Body */
1057 if (*n <= 1)
1058 {
1059 return 0;
1060 }
1061 l = *n / 2 + 1;
1062 r = *n;
1063 L2:
1064 if (l <= 1)
1065 {
1066 goto L20;
1067 }
1068 --l;
1069 rec = record[l];
1070 crit = criter[l];
1071 goto L3;
1072 L20:
1073 rec = record[r];
1074 crit = criter[r];
1075 record[r] = record[1];
1076 criter[r] = criter[1];
1077 --r;
1078 if (r == 1)
1079 {
1080 goto L999;
1081 }
1082 L3:
1083 j = l;
1084 L4:
1085 i = j;
1086 j <<= 1;
1087 if ((i_1 = j - r) < 0)
1088 {
1089 goto L5;
1090 }
1091 else if (i_1 == 0)
1092 {
1093 goto L6;
1094 }
1095 else
1096 {
1097 goto L8;
1098 }
1099 L5:
1100 if (criter[j] < criter[j + 1])
1101 {
1102 ++j;
1103 }
1104 L6:
1105 if (crit >= criter[j])
1106 {
1107 goto L8;
1108 }
1109 record[i] = record[j];
1110 criter[i] = criter[j];
1111 goto L4;
1112 L8:
1113 record[i] = rec;
1114 criter[i] = crit;
1115 goto L2;
1116 L999:
1117 record[1] = rec;
1118 criter[1] = crit;
1119 return 0;
1120 } /* mshtr1_ */
1121
1122 /* ********************************************************************** */
1123 int
mshcvx_(long direct,long * c,long * nu,long * pfold,long * err)1124 mshcvx_ (long direct, long *c, long *nu, long *pfold, long *err)
1125 {
1126 static long t, a4, a5, i1, i2, i3, i4, i5, i6, s1, s2, s3, t4, t5,
1127 pf, pp, ps;
1128
1129 // extern int mshopt_();
1130 static long tt4, tt5, ppf, psf;
1131 LONG8 det;
1132
1133 /* Parameter adjustments */
1134 nu -= 7;
1135 c -= 3;
1136
1137 /* Function Body */
1138 if (direct)
1139 {
1140 pp = 3;
1141 ps = 4;
1142 i1 = 1;
1143 i2 = 3;
1144 i3 = 2;
1145 i4 = 6;
1146 i5 = 5;
1147 i6 = 4;
1148 }
1149 else
1150 {
1151 pp = 4;
1152 ps = 3;
1153 i1 = 1;
1154 i2 = 2;
1155 i3 = 3;
1156 i4 = 4;
1157 i5 = 5;
1158 i6 = 6;
1159 }
1160 L10:
1161 ppf = *pfold;
1162 pf = nu[ps + *pfold * 6];
1163 psf = nu[ps + pf * 6];
1164 s1 = nu[ppf * 6 + 1];
1165 s2 = nu[pf * 6 + 1];
1166 s3 = nu[psf * 6 + 1];
1167 det = (LONG8) (c[(s2 << 1) + 1] - c[(s1 << 1) + 1])
1168 * (LONG8) (c[(s3 << 1) + 2] - c[(s1 << 1) + 2])
1169 - (LONG8) (c[(s2 << 1) + 2] - c[(s1 << 1) + 2])
1170 * (LONG8) (c[(s3 << 1) + 1] - c[(s1 << 1) + 1]);
1171
1172 if ((!(direct) && det > 0) || (direct && det < 0))
1173 {
1174 /* on ajoute un triangle t et on detruit une arete */
1175 /* ----------------------------------------------- */
1176 if (direct)
1177 {
1178 tt4 = nu[ppf * 6 + 2];
1179 tt5 = nu[pf * 6 + 2];
1180 }
1181 else
1182 {
1183 tt4 = nu[pf * 6 + 2];
1184 tt5 = nu[psf * 6 + 2];
1185 }
1186 t4 = tt4 / 8;
1187 t5 = tt5 / 8;
1188 a4 = tt4 - (t4 << 3);
1189 a5 = tt5 - (t5 << 3);
1190 /* destruction de l'arete frontiere en pf */
1191 /* -------------------------------------- */
1192 nu[ps + ppf * 6] = psf;
1193 nu[pp + psf * 6] = ppf;
1194 /* on remplace l'arete frontiere par l'element genere */
1195 /* --------------------------------------------------- */
1196 t = pf;
1197 /* update de l'arete non detruite */
1198 /* ------------------------------ */
1199 if (direct)
1200 {
1201 nu[ppf * 6 + 2] = (t << 3) + i6;
1202 }
1203 else
1204 {
1205 nu[psf * 6 + 2] = (t << 3) + i6;
1206 }
1207 /* on cree l'element */
1208 /* ----------------- */
1209 nu[i1 + t * 6] = s1;
1210 nu[i2 + t * 6] = s2;
1211 nu[i3 + t * 6] = s3;
1212 nu[i4 + t * 6] = (t4 << 3) + a4;
1213 nu[i5 + t * 6] = (t5 << 3) + a5;
1214 if (direct)
1215 {
1216 nu[i6 + t * 6] = -ppf;
1217 }
1218 else
1219 {
1220 nu[i6 + t * 6] = -psf;
1221 }
1222 nu[a4 + t4 * 6] = (t << 3) + i4;
1223 nu[a5 + t5 * 6] = (t << 3) + i5;
1224 mshopt_ (&c[3], &nu[7], &t5, a5, err);
1225 if (*err != 0)
1226 {
1227 return 0;
1228 }
1229 goto L10;
1230 }
1231 return 1;
1232 } /* mshcvx_ */
1233
1234 /* ********************************************************************** */
1235 int
mshcxi_(long * c,long * nu,long * tri,long * nbs,long * tete,long * err)1236 mshcxi_ (long *c, long *nu, long *tri, long *nbs, long *tete, long *err)
1237 {
1238 /* System generated locals */
1239 long i_1;
1240
1241 /* Local variables */
1242 static long sfree, ttaf, i, j, s, t, pf;
1243
1244 // extern long mshlcl_();
1245 // extern int mshcvx_(), mshopt_();
1246 static long iaf, taf, npf, ppf, psf;
1247
1248 /* initialisation de la sfree liste dans nu */
1249 /* Parameter adjustments */
1250 --tri;
1251 nu -= 7;
1252 c -= 3;
1253
1254 /* Function Body */
1255 i_1 = *nbs + *nbs - 2;
1256 for (i = 1; i <= i_1; ++i)
1257 {
1258 nu[i * 6 + 1] = i + 1;
1259 for (j = 2; j <= 6; ++j)
1260 {
1261 nu[j + i * 6] = 0;
1262 }
1263 }
1264 nu[(*nbs + *nbs - 2) * 6 + 1] = 0;
1265 sfree = 1;
1266 /* initialisation du premier triangle */
1267 t = sfree;
1268 sfree = nu[sfree * 6 + 1];
1269 /* initialisation de la liste frontiere */
1270 *tete = sfree;
1271 pf = sfree;
1272 for (i = 1; i <= 3; ++i)
1273 {
1274 nu[i + t * 6] = tri[i];
1275 nu[i + 3 + t * 6] = -pf;
1276 ppf = pf;
1277 sfree = nu[pf * 6 + 1];
1278 pf = sfree;
1279 if (i == 3)
1280 {
1281 pf = *tete;
1282 }
1283 nu[ppf * 6 + 1] = tri[i];
1284 nu[ppf * 6 + 2] = i + 3 + (t << 3);
1285 nu[ppf * 6 + 4] = pf;
1286 nu[pf * 6 + 3] = ppf;
1287 }
1288 i_1 = *nbs;
1289 for (i = 4; i <= i_1; ++i)
1290 {
1291 s = tri[i];
1292 pf = mshlcl_ (&c[3], &nu[7], tete, &s);
1293 /* creation d'un nouveau triangle et modification de la
1294 frontiere */
1295 /*
1296 --------------------------------------------------------------
1297 */
1298 t = sfree;
1299 sfree = nu[sfree * 6 + 1];
1300 npf = sfree;
1301 sfree = nu[sfree * 6 + 1];
1302 ppf = nu[pf * 6 + 3];
1303 psf = nu[pf * 6 + 4];
1304 ttaf = nu[pf * 6 + 2];
1305 taf = ttaf / 8;
1306 iaf = ttaf - (taf << 3);
1307
1308 /* npf */
1309 /* 1 x s --- */
1310 /* / \ --- */
1311 /* 4 / \ 6 --- vide --- */
1312 /* / t \ --- */
1313 /* 2 / 5 \ 3 --- */
1314 /* ------ --<---x---------x---------x- frontiere--<--- */
1315 /* psf \ iaf / pf --- */
1316 /* \ taf / --- omega --- */
1317 /* \ / --- */
1318 /* \ / --- */
1319 /* x --- */
1320 /* --- */
1321 /* generation de l'element t */
1322 nu[t * 6 + 1] = s;
1323 nu[t * 6 + 2] = nu[psf * 6 + 1];
1324 nu[t * 6 + 3] = nu[pf * 6 + 1];
1325 nu[t * 6 + 4] = -npf;
1326 nu[t * 6 + 5] = (taf << 3) + iaf;
1327 nu[t * 6 + 6] = -pf;
1328 nu[iaf + taf * 6] = (t << 3) + 5;
1329 /* update de la liste frontiere */
1330 nu[npf * 6 + 4] = psf;
1331 nu[pf * 6 + 4] = npf;
1332 nu[npf * 6 + 3] = pf;
1333 nu[psf * 6 + 3] = npf;
1334 nu[npf * 6 + 1] = s;
1335 nu[npf * 6 + 2] = (t << 3) + 4;
1336 nu[pf * 6 + 2] = (t << 3) + 6;
1337 mshopt_ (&c[3], &nu[7], &t, 5L, err);
1338 if (*err != 0)
1339 {
1340 return 0;
1341 }
1342 mshcvx_ (1, &c[3], &nu[7], &npf, err);
1343 if (*err != 0)
1344 {
1345 return 0;
1346 }
1347 mshcvx_ (0, &c[3], &nu[7], &npf, err);
1348 if (*err != 0)
1349 {
1350 return 0;
1351 }
1352 }
1353 return 1;
1354 } /* mshcxi_ */
1355
1356 /* ********************************************************************** */
1357 int
mshopt_(long * c,long * nu,long * t,long a,long * err)1358 mshopt_ (long *c, long *nu, long *t, long a, long *err)
1359 {
1360 /* Initialized data */
1361
1362 static long mod3[3] =
1363 {2, 3, 1};
1364
1365 /* System generated locals */
1366 long i_1;
1367 double d_1;
1368
1369 /* Local variables */
1370 static long pile[4096] /* was [2][256] */ ;
1371 static double reel1, reel2;
1372 static double reel8;
1373 static long i, a1, a2, s1, t1, t2, s2, s3, s4, aa, i11, i12, i13,
1374 i21, i22, i23, tt;
1375 static long tt1;
1376 LONG8 cos1, cos2, sin1, sin2, sgn;
1377 int kerr21=0;
1378 /* Parameter adjustments */
1379 nu -= 7;
1380 c -= 3;
1381
1382 /* Function Body */
1383 i = 1;
1384 pile[(i << 1) - 2] = *t;
1385 pile[(i << 1) - 1] = a;
1386 L10:
1387 if (i > 0)
1388 {
1389 t1 = pile[(i << 1) - 2];
1390 a1 = pile[(i << 1) - 1];
1391 --i;
1392 if (t1 <= 0)
1393 {
1394 goto L10;
1395 }
1396 tt1 = nu[a1 + t1 * 6];
1397 if (tt1 <= 0)
1398 {
1399 goto L10;
1400 }
1401 t2 = tt1 / 8;
1402 a2 = tt1 - (t2 << 3);
1403 i11 = a1 - 3;
1404 i12 = mod3[i11 - 1];
1405 i13 = mod3[i12 - 1];
1406 i21 = a2 - 3;
1407 i22 = mod3[i21 - 1];
1408 i23 = mod3[i22 - 1];
1409 s1 = nu[i13 + t1 * 6];
1410 s2 = nu[i11 + t1 * 6];
1411 s3 = nu[i12 + t1 * 6];
1412 s4 = nu[i23 + t2 * 6];
1413 sin1 = (LONG8) (c[(s3 << 1) + 2] - c[(s1 << 1) + 2])
1414 * (LONG8) (c[(s2 << 1) + 1] - c[(s1 << 1) + 1])
1415 - (LONG8) (c[(s3 << 1) + 1] - c[(s1 << 1) + 1])
1416 * (LONG8) (c[(s2 << 1) + 2] - c[(s1 << 1) + 2]);
1417
1418 cos1 = (LONG8) (c[(s3 << 1) + 1] - c[(s1 << 1) + 1])
1419 * (LONG8) (c[(s3 << 1) + 1] - c[(s2 << 1) + 1])
1420 + (LONG8) (c[(s3 << 1) + 2] - c[(s1 << 1) + 2])
1421 * (LONG8) (c[(s3 << 1) + 2] - c[(s2 << 1) + 2]);
1422 if (sin1 == 0 && cos1 == 0)
1423 {
1424 *err = 20;
1425 return 0;
1426 }
1427 /* b est la cotangente de angle (s1,s3,s2) */
1428
1429 sin2 = (LONG8) (c[(s4 << 1) + 1] - c[(s1 << 1) + 1])
1430 * (LONG8) (c[(s2 << 1) + 2] - c[(s1 << 1) + 2])
1431 - (LONG8) (c[(s4 << 1) + 2] - c[(s1 << 1) + 2])
1432 * (LONG8) (c[(s2 << 1) + 1] - c[(s1 << 1) + 1]);
1433 // FH correct 6/11/2005 forgotted cast
1434 cos2 = (LONG8) (c[(s4 << 1) + 1] - c[(s2 << 1) + 1])
1435 * (LONG8) (c[(s4 << 1) + 1] - c[(s1 << 1) + 1])
1436 + (LONG8) (c[(s4 << 1) + 2] - c[(s2 << 1) + 2])
1437 * (LONG8) (c[(s4 << 1) + 2] - c[(s1 << 1) + 2]);
1438
1439 reel1 = (double) cos2 *(double) sin1;
1440 reel2 = (double) cos1 *(double) sin2;
1441
1442 if (aabs (reel1) + aabs (reel2) >= (double ) DLONG8LONG8)
1443 {
1444 reel8 = (double) cos2 *(double) sin1 + (double) cos1 *(double) sin2;
1445
1446 /* Computing MIN */
1447 d_1 = amax (reel8, -1.);
1448 reel8 = amin (d_1, 1.);
1449 sgn = (LONG8) reel8;
1450 }
1451 else
1452 {
1453 sgn = cos2 * sin1 + cos1 * sin2;
1454 }
1455 /* Computing MIN */
1456 i_1 = amin(amax (sgn, -1),1);
1457 // cout << sgn << " " << i_1 << endl;
1458 if ( i_1 * sin1 >= 0)
1459 {
1460 goto L10;
1461 }
1462 /* on inverse le quadrilatere */
1463 /* update des sommets */
1464 /* ------------------------- */
1465 nu[i12 + t1 * 6] = s4;
1466 nu[i22 + t2 * 6] = s1;
1467 /* update des aretes a1,a2 */
1468 /* ------------------------------- */
1469 tt1 = nu[i22 + 3 + t2 * 6];
1470 nu[a1 + t1 * 6] = tt1;
1471 if (tt1 > 0)
1472 {
1473 tt = tt1 / 8;
1474 aa = tt1 - (tt << 3);
1475 nu[aa + tt * 6] = a1 + (t1 << 3);
1476 }
1477 else if (tt1 != nothing)
1478 {
1479 nu[-tt1 * 6 + 2] = a1 + (t1 << 3);
1480 }
1481 tt1 = nu[i12 + 3 + t1 * 6];
1482 nu[a2 + t2 * 6] = tt1;
1483 if (tt1 > 0)
1484 {
1485 tt = tt1 / 8;
1486 aa = tt1 - (tt << 3);
1487 nu[aa + tt * 6] = a2 + (t2 << 3);
1488 }
1489 else if (tt1 != nothing)
1490 {
1491 nu[-tt1 * 6 + 2] = a2 + (t2 << 3);
1492 }
1493 nu[i12 + 3 + t1 * 6] = i22 + 3 + (t2 << 3);
1494 nu[i22 + 3 + t2 * 6] = i12 + 3 + (t1 << 3);
1495 if (i + 4 > 2048) /* Modif F. H june 2015 */
1496 {
1497 if(kerr21)
1498 cout << " Bizarre mshptg err 21 pile too small (continue) "<< endl;
1499 if(kerr21++<10000)
1500 goto L10;
1501 *err = 21;
1502 return 0;
1503 }
1504 ++i;
1505 pile[(i << 1) - 2] = t1;
1506 pile[(i << 1) - 1] = a1;
1507 ++i;
1508 pile[(i << 1) - 2] = t2;
1509 pile[(i << 1) - 1] = a2;
1510 ++i;
1511 pile[(i << 1) - 2] = t1;
1512 pile[(i << 1) - 1] = i13 + 3;
1513 ++i;
1514 pile[(i << 1) - 2] = t2;
1515 pile[(i << 1) - 1] = i23 + 3;
1516 goto L10;
1517 }
1518 return 1;
1519 } /* mshopt_ */
1520
1521 /* ********************************************************************** */
1522 int
mshfrt_(long * c,long * nu,long * nbs,long * arete,long nba,long * sd,long nbsd,long * reft,long * w,long * err)1523 mshfrt_ (long *c, long *nu, long *nbs, long *arete, long nba, long *sd,
1524 long nbsd, long *reft, long *w, long *err)
1525 {
1526 /* Initialized data */
1527
1528 static long p3[3] =
1529 {2, 3, 1};
1530
1531 /* System generated locals */
1532 long i_1, i_2, i_3;
1533
1534 /* Local variables */
1535 static long nbac, ifrt, a, i, t, itera, s1, s2;
1536
1537 // extern int mshfr1_();
1538 static long ie, ap, ta, is, nbacpp;
1539 static long is1, ss1, s2t, s3t, isd, jsd, nbt, err1;
1540 LONG8 det2, det3;
1541
1542 /* Parameter adjustments */
1543 --w;
1544 --reft;
1545 sd -= 3;
1546 arete -= 3;
1547 nu -= 7;
1548 c -= 3;
1549
1550 /* Function Body */
1551 if (nba == 0)
1552 {
1553 return 0;
1554 }
1555 ifrt = 0;
1556 nbt = *nbs + *nbs - 2;
1557 i_1 = *nbs;
1558 for (i = 1; i <= i_1; ++i)
1559 {
1560 reft[i] = 0;
1561 }
1562 i_1 = nba;
1563 for (i = 1; i <= i_1; ++i)
1564 {
1565 reft[arete[(i << 1) + 1]] = nothing;
1566 reft[arete[(i << 1) + 2]] = nothing;
1567 }
1568 nbac = 0;
1569 i_1 = nba;
1570 for (a = 1; a <= i_1; ++a)
1571 {
1572 /* Computing MIN */
1573 i_2 = arete[(a << 1) + 1], i_3 = arete[(a << 1) + 2];
1574 s1 = amin (i_2, i_3);
1575 /* Computing MAX */
1576 i_2 = arete[(a << 1) + 1], i_3 = arete[(a << 1) + 2];
1577 s2 = amax (i_2, i_3);
1578 if (s1 == s2)
1579 {
1580 ++nbac;
1581 }
1582 else
1583 {
1584 i = reft[s1];
1585 L25:
1586 if (i != nothing)
1587 {
1588 /* Computing MAX */
1589 i_2 = arete[(i << 1) + 1], i_3 = arete[(i << 1) + 2];
1590 if (s2 == amax (i_2, i_3))
1591 {
1592 ++nbac;
1593 }
1594 else
1595 {
1596 i = w[i];
1597 goto L25;
1598 }
1599 }
1600 else
1601 {
1602 w[a] = reft[s1];
1603 reft[s1] = a;
1604 }
1605 }
1606 }
1607 nbacpp = 1;
1608 itera = 0;
1609 err1 = 0;
1610 L50:
1611 ++itera;
1612 if (err1 != 0)
1613 {
1614 *err = err1;
1615 return 0;
1616 }
1617 if (nbac < nba)
1618 {
1619 if (nbacpp == 0)
1620 {
1621 i_1 = *nbs;
1622 for (i = 1; i <= i_1; ++i)
1623 {
1624 a = reft[i];
1625 L60:
1626 if (a > 0)
1627 {
1628 s1 = arete[(i << 1) + 1];
1629 s2 = arete[(i << 1) + 2];
1630 a = w[a];
1631 goto L60;
1632 }
1633 }
1634 *err = 7;
1635 return 0;
1636 }
1637 /* --------------------------------------------------------------------- */
1638 /* on s'occupe des aretes a forcer */
1639 /* --------------------------------------------------------------------- */
1640 nbacpp = 0;
1641 i_1 = nbt;
1642 for (ie = 1; ie <= i_1; ++ie)
1643 {
1644 if (nu[ie * 6 + 5] != 0)
1645 {
1646 for (is = 1; is <= 3; ++is)
1647 {
1648 s1 = nu[is + ie * 6];
1649 s2t = nu[p3[is - 1] + ie * 6];
1650 ss1 = amin (s1, s2t);
1651 ap = 0;
1652 a = reft[ss1];
1653 L80:
1654 if (a > 0)
1655 {
1656 /* Computing MAX */
1657 i_2 = arete[(a << 1) + 1], i_3 = arete[(a << 1) + 2];
1658 s2 = amax (i_2, i_3);
1659 t = ie;
1660 ta = 0;
1661 if (s2 == amax (s1, s2t))
1662 {
1663 if (nu[is + 3 + ie * 6] > 0)
1664 {
1665 ta = nu[is + 3 + ie * 6] / 8;
1666 i = nu[is + 3 + ie * 6] - (ta << 3);
1667 nu[i + ta * 6] = nothing;
1668 }
1669 nu[is + 3 + ie * 6] = nothing;
1670 goto L100;
1671 }
1672 ap = a;
1673 a = w[a];
1674 goto L80;
1675 }
1676 if (itera == 1)
1677 {
1678 goto L110;
1679 }
1680 ss1 = s1;
1681 ap = 0;
1682 a = reft[ss1];
1683 L90:
1684 if (a > 0)
1685 {
1686 /* Computing MAX */
1687 i_2 = arete[(a << 1) + 1], i_3 = arete[(a << 1) + 2];
1688 s2 = amax (i_2, i_3);
1689 t = ie;
1690 ta = 0;
1691 /* recherche si l' element coupe l''arete a */
1692 is1 = is;
1693 s3t = nu[p3[p3[is - 1] - 1] + t * 6];
1694 det2 = (LONG8) (c[(s2t << 1) + 1] - c[(s1 << 1) + 1]) * (LONG8) (c[( s2 << 1) + 2] - c[(s1 << 1) + 2])
1695 - (LONG8) (c[(s2t << 1) + 2] - c[(s1 << 1) + 2]) * (LONG8) (c[( s2 << 1) + 1] - c[(s1 << 1) + 1]);
1696 det3 = (LONG8) (c[(s3t << 1) + 1] - c[(s1 << 1) + 1]) * (LONG8) (c[( s2 << 1) + 2] - c[(s1 << 1) + 2])
1697 - (LONG8) (c[(s3t << 1) + 2] - c[(s1 << 1) + 2]) * (LONG8) (c[( s2 << 1) + 1] - c[(s1 << 1) + 1]);
1698 if (det2 > 0 && det3 < 0)
1699 {
1700 mshfr1_ (&c[3], &nu[7], &t, &ta, &is1, &s2, err);
1701 if (*err != 0)
1702 {
1703 printf(" Error %ld inforce edge %ld = %ld;%ld with vertices %ld %ld \n",*err,a,i_2, i_3, is1,s2 );
1704 return 0;
1705 }
1706 goto L100;
1707 }
1708 else if (det2 == 0 && det3 < 0 && reft[s2t] == 0)
1709 {
1710 err1 = 10;
1711 printf(" det = %ld %ld %ld %ld %ld == %ld \n ",
1712 (c[(s2t << 1) + 1] - c[(s1 << 1) + 1]), (c[( s2 << 1) + 2] - c[(s1 << 1) + 2]) ,
1713 (c[(s2t << 1) + 2] - c[(s1 << 1) + 2]), (c[( s2 << 1) + 1] - c[(s1 << 1) + 1]) ,
1714 (c[(s2t << 1) + 1] - c[(s1 << 1) + 1]) * (c[( s2 << 1) + 2] - c[(s1 << 1) + 2]) ,
1715 (c[(s2t << 1) + 2] - c[(s1 << 1) + 2]) * (c[( s2 << 1) + 1] - c[(s1 << 1) + 1])
1716 );
1717
1718
1719 printf("bug 2, mshptg: point %ld is on boundary edge %ld %ld \n",s2t,i_2,i_3);
1720 }
1721 else if (det2 > 0 && det3 == 0 && reft[s3t] == 0)
1722 {
1723 err1 = 10;
1724 printf(" det = %ld %ld %ld %ld %ld %ld \n ",
1725 (c[(s3t << 1) + 1] - c[(s1 << 1) + 1]), (c[( s2 << 1) + 2] - c[(s1 << 1) + 2]) ,
1726 (c[(s3t << 1) + 2] - c[(s1 << 1) + 2]), (c[( s2 << 1) + 1] - c[(s1 << 1) + 1]) ,
1727 (c[(s3t << 1) + 1] - c[(s1 << 1) + 1]) * (c[( s2 << 1) + 2] - c[(s1 << 1) + 2]) ,
1728 (c[(s3t << 1) + 2] - c[(s1 << 1) + 2]) * (c[( s2 << 1) + 1] - c[(s1 << 1) + 1])
1729 );
1730 printf("bug 2, mshptg: point %ld is on boundary %ld %ld\n",s3t,i_2,i_3);
1731 }
1732 ap = a;
1733 a = w[a];
1734 goto L90;
1735 }
1736 goto L110;
1737 L100:
1738 ++nbacpp;
1739 if (ap == 0)
1740 {
1741 reft[ss1] = w[a];
1742 }
1743 else
1744 {
1745 w[ap] = w[a];
1746 }
1747 if (nbac + nbacpp == nba)
1748 {
1749 goto L130;
1750 }
1751 L110:
1752 ;
1753 }
1754 }
1755 }
1756 nbac += nbacpp;
1757 goto L50;
1758 }
1759 L130:
1760 /* ----------------------------------------------------------------------- */
1761 /* prise en compte des sous domaines */
1762 /* ----------------------------------------------------------------------- */
1763 /* add FH si pas de nbsd --- jan 2004 */
1764
1765 if (nbsd<0) {
1766 i_1 = nbt;
1767 for (ie = 1; ie <= i_1; ++ie)
1768 {
1769 if (nu[ie * 6 + 1] > 0)
1770 {
1771 nu[ie * 6 + 1] = -nu[ie * 6 + 1];
1772 }
1773 reft[ie]=1;
1774 }
1775
1776 goto L205; // pas triangle retire
1777 }
1778 if (nbsd == 0) {
1779 long i__1 = nbt,nbssd,exter,i__,headt,nst,j;
1780 for (t = 1; t <= i__1; ++t) {
1781 reft[t] = -1073741824;
1782 }
1783 nbssd = 0;
1784 /* if(impre.gt.1) print *,'nbsd .eq. 0 => recherche de ssd' */
1785 i__1 = nbt;
1786 for (t = 1; t <= i__1; ++t) {
1787 if (nu[t * 6 + 1] > 0 && nu[t * 6 + 5] != 0) {
1788 ++nbssd;
1789 exter = 0;
1790 i__ = 2;
1791 w[i__ - 1] = t;
1792 w[i__] = 3;
1793 reft[t] = 0;
1794 headt = t;
1795 nu[t * 6 + 1] = -nu[t * 6 + 1];
1796 nst = 1;
1797 /* print *,' ssd ',nbssd */
1798 L131:
1799 if (i__ > 0) {
1800 ++w[i__];
1801 if (w[i__] <= 6) {
1802 ta = nu[w[i__] + w[i__ - 1] * 6];
1803 if (ta <= 0) {
1804 if (ta != -1073741824) {
1805 exter = 1;
1806 }
1807 } else if (ta > 0) {
1808 ta /= 8;
1809 if (nu[ta * 6 + 1] > 0) {
1810 /* print *,ta */
1811 nu[ta * 6 + 1] = -nu[ta * 6 + 1];
1812 ++nst;
1813 reft[ta] = headt;
1814 headt = ta;
1815 w[i__ + 1] = ta;
1816 w[i__ + 2] = 3;
1817 i__ += 2;
1818 }
1819 }
1820 } else {
1821 i__ += -2;
1822 }
1823 goto L131;
1824 }
1825
1826 if (exter) {
1827 --nbssd;
1828 i__ = headt;
1829 L133:
1830 if (i__ > 0) {
1831 j = reft[i__];
1832 /* print *,i */
1833 reft[i__] = 0;
1834 nu[i__ * 6 + 1] = 0;
1835 i__ = j;
1836 goto L133;
1837 }
1838 } else {
1839 i__ = headt;
1840 L136:
1841 if (i__ > 0) {
1842 j = reft[i__];
1843 reft[i__] = nbssd;
1844 i__ = j;
1845 goto L136;
1846 }
1847 }
1848 }
1849 }
1850 goto L205;
1851 }
1852 /* fin ajoute FH jan 2004 */
1853
1854
1855 i_1 = *nbs + nbsd + nbsd;
1856 for (i = 1; i <= i_1; ++i)
1857 {
1858 w[i] = 0;
1859 }
1860 i_1 = nbsd;
1861 for (i = 1; i <= i_1; ++i)
1862 {
1863 a = (i_2 = sd[(i << 1) + 1], aabs (i_2));
1864 /* Computing MIN */
1865 i_2 = arete[(a << 1) + 1], i_3 = arete[(a << 1) + 2];
1866 s1 = amin (i_2, i_3);
1867 w[i + i] = w[s1 + nbsd + nbsd];
1868 w[s1 + nbsd + nbsd] = i;
1869 }
1870 i_1 = nbt;
1871 for (t = 1; t <= i_1; ++t)
1872 {
1873 reft[t] = nothing;
1874 if (nu[t * 6 + 6] != 0)
1875 {
1876 for (i = 1; i <= 3; ++i)
1877 {
1878 /* Computing MIN */
1879 i_2 = nu[i + t * 6], i_3 = nu[p3[i - 1] + t * 6];
1880 ss1 = amin (i_2, i_3);
1881 jsd = nbsd + nbsd + ss1;
1882 L160:
1883 isd = w[jsd];
1884 if (isd > 0)
1885 {
1886 a = sd[(isd << 1) + 1];
1887 if (a > 0)
1888 {
1889 if (nu[i + t * 6] == arete[(a << 1) + 1] && nu[p3[i -
1890 1] + t * 6] == arete[(a << 1) + 2])
1891 {
1892 reft[t] = sd[(isd << 1) + 2];
1893 w[isd + isd - 1] = t;
1894 w[jsd] = w[isd + isd];
1895 goto L170;
1896 }
1897 }
1898 else if (a < 0)
1899 {
1900 if (nu[i + t * 6] == arete[(-a << 1) + 2] && nu[p3[i
1901 - 1] + t * 6] == arete[(-a << 1) + 1])
1902 {
1903 reft[t] = sd[(isd << 1) + 2];
1904 w[isd + isd - 1] = t;
1905 w[jsd] = w[isd + isd];
1906 goto L170;
1907 }
1908 }
1909 else
1910 {
1911 cout << " Err sous domaine " << isd << "ref par d'arete +/- = " << a << " n'est dans aucun triangle " << endl;
1912 *err = 11;
1913 }
1914 jsd = isd + isd;
1915 goto L160;
1916 }
1917 L170:
1918 ;
1919 }
1920 }
1921 }
1922 i_1 = nbsd;
1923 for (isd = 1; isd <= i_1; ++isd)
1924 {
1925 if (w[isd + isd - 1] == 0)
1926 {
1927 cout << " Err sous domaine " << isd << "ref par d'arete +/- = " << sd[(isd << 1) + 1] << " n'est dans aucun triangle " << endl;
1928 *err = 11;
1929 }
1930 else
1931 {
1932 w[isd + isd] = 3;
1933 }
1934 }
1935 if (*err != 0)
1936 {
1937 return 0;
1938 }
1939 i = nbsd + nbsd;
1940 L200:
1941 if (i > 0)
1942 {
1943 ++w[i];
1944 if (w[i] <= 6)
1945 {
1946 ta = nu[w[i] + w[i - 1] * 6];
1947 if (ta > 0)
1948 {
1949 ta /= 8;
1950 if (nu[ta * 6 + 1] > 0)
1951 {
1952 nu[ta * 6 + 1] = -nu[ta * 6 + 1];
1953 if (reft[ta] != reft[w[i - 1]])
1954 {
1955 if (reft[ta] != nothing)
1956 {
1957 }
1958 else
1959 {
1960 reft[ta] = reft[w[i - 1]];
1961 }
1962 w[i + 1] = ta;
1963 w[i + 2] = 3;
1964 i += 2;
1965 }
1966 }
1967 }
1968 }
1969 else
1970 {
1971 ta = w[i-1] ;
1972 if( nu[ta * 6 + 1]>=0) nu[ta * 6 + 1]= -nu[ta * 6 + 1]; // pour etre sur que le traingle est bien marque
1973
1974 i += -2;
1975 }
1976 goto L200;
1977 }
1978 L205:
1979 i_1 = nbt;
1980 for (ie = 1; ie <= i_1; ++ie)
1981 {
1982 if (nu[ie * 6 + 1] < 0)
1983 {
1984 nu[ie * 6 + 1] = -nu[ie * 6 + 1];
1985 }
1986 else
1987 {
1988 for (i = 1; i <= 6; ++i)
1989 {
1990 nu[i + ie * 6] = 0;
1991 }
1992 }
1993 }
1994 return 1;
1995 } /* mshfrt_ */
1996
1997 /* ********************************************************************** */
1998 int
mshfr1_(long * c,long * nu,long * it1,long * ita,long * is1,long * s2,long * err)1999 mshfr1_ (long *c, long *nu, long *it1, long *ita, long *is1, long *s2, long *err)
2000 {
2001 /* Initialized data */
2002
2003 static long p3[5] =
2004 {2, 3, 1, 2, 3};
2005
2006 static long nbac, t, x, y, l1, l2, l3, s1, s3;
2007
2008 // extern int mshfr2_();
2009 static long la, ta;
2010 static long s2t, s3t, lst[768] /* was [3][256] */ ;
2011 LONG8 det;
2012
2013 /* Parameter adjustments */
2014 nu -= 7;
2015 c -= 3;
2016
2017 /* Function Body */
2018 t = *it1;
2019 s1 = nu[*is1 + t * 6];
2020 x = c[(*s2 << 1) + 1] - c[(s1 << 1) + 1];
2021 y = c[(*s2 << 1) + 2] - c[(s1 << 1) + 2];
2022 nbac = 0;
2023 l1 = *is1;
2024 l2 = p3[l1 - 1];
2025 l3 = p3[l2 - 1];
2026 s2t = nu[l2 + t * 6];
2027 s3t = nu[l3 + t * 6];
2028 la = l2 + 3;
2029 L20:
2030 ++nbac;
2031 if (nbac > 256)
2032 {
2033 *err = 8;
2034 return 0;
2035 }
2036 lst[nbac * 3 - 2] = t;
2037 lst[nbac * 3 - 1] = la;
2038 ta = nu[la + t * 6];
2039 if (ta <= 0)
2040 {
2041 *err = 9;
2042 long sa= nu[ (la+1)%3 + t * 6], sb=nu[(la+2)%3 + t * 6];
2043 printf("mshptg: bug boundary crossing %ld,%ld , %ld %ld \n", s1, *s2 , sa,sb);
2044 return 0;
2045 }
2046 t = ta / 8;
2047 la = ta - (t << 3);
2048 s3 = nu[p3[la - 3] + t * 6];
2049 if (s3 != *s2)
2050 {
2051 det = (LONG8) x * (LONG8) (c[(s3 << 1) + 2] - c[(s1 << 1) + 2])
2052 - (LONG8) y * (LONG8) (c[(s3 << 1) + 1] - c[(s1 << 1) + 1]);
2053 if (det > 0)
2054 {
2055 la = p3[la - 4] + 3;
2056 }
2057 else if (det < 0)
2058 {
2059 la = p3[la - 3] + 3;
2060 }
2061 else
2062 {
2063 printf("mshptg: bug the point %ld is on boundary \n", s3 );
2064 *err = 10+s3*10;
2065 return 0;
2066 }
2067 goto L20;
2068 }
2069 mshfr2_ (&c[3], &nu[7], lst, &nbac, it1, ita, &s1, s2, err);
2070 return 0;
2071 } /* mshfr1_ */
2072
2073 /* ********************************************************************** */
2074 int
mshfr2_(long * c,long * nu,long * lst,long * nbac,long * t,long * ta,long * ss1,long * ss2,long * err)2075 mshfr2_ (long *c, long *nu, long *lst, long *nbac, long *t, long *ta,
2076 long *ss1, long *ss2, long *err)
2077 {
2078 /* Initialized data */
2079
2080 static long mod3[3] =
2081 {2, 3, 1};
2082
2083 /* System generated locals */
2084 long i_1;
2085
2086 /* Local variables */
2087 static long i, x, y, a1, a2, pplst, s1, pslst, ptlst, s2, s3, s4,
2088 ttlst, t1, t2, aa, i11, i12, i13, i21, i22, i23, x41,
2089 y41, tt;
2090
2091 // extern int mshopt_();
2092 static long tt1, aas;
2093 LONG8 det1, det2, det3, det4;
2094
2095 /* Parameter adjustments */
2096 lst -= 4;
2097 nu -= 7;
2098 c -= 3;
2099
2100 /* Function Body */
2101 x = c[(*ss1 << 1) + 1] - c[(*ss2 << 1) + 1];
2102 y = c[(*ss1 << 1) + 2] - c[(*ss2 << 1) + 2];
2103 i_1 = *nbac - 1;
2104 for (i = 1; i <= i_1; ++i)
2105 {
2106 lst[i * 3 + 1] = i + 1;
2107 }
2108 lst[*nbac * 3 + 1] = 0;
2109 ttlst = 1;
2110 L20:
2111 ptlst = ttlst;
2112 pplst = 0;
2113 L30:
2114 if (ptlst > 0)
2115 {
2116 t1 = lst[ptlst * 3 + 2];
2117 a1 = lst[ptlst * 3 + 3];
2118 tt1 = nu[a1 + t1 * 6];
2119 t2 = tt1 / 8;
2120 a2 = tt1 - (t2 << 3);
2121 i11 = a1 - 3;
2122 i12 = mod3[i11 - 1];
2123 i13 = mod3[i12 - 1];
2124 i21 = a2 - 3;
2125 i22 = mod3[i21 - 1];
2126 i23 = mod3[i22 - 1];
2127 s1 = nu[i13 + t1 * 6];
2128 s2 = nu[i11 + t1 * 6];
2129 s3 = nu[i12 + t1 * 6];
2130 s4 = nu[i23 + t2 * 6];
2131 x41 = c[(s4 << 1) + 1] - c[(s1 << 1) + 1];
2132 y41 = c[(s4 << 1) + 2] - c[(s1 << 1) + 2];
2133 det2 = (LONG8) (c[(s2 << 1) + 1] - c[(s1 << 1) + 1]) * (LONG8) y41
2134 - (LONG8) (c[(s2 << 1) + 2] - c[(s1 << 1) + 2]) * (LONG8) x41;
2135
2136 det3 = (LONG8) (c[(s3 << 1) + 1] - c[(s1 << 1) + 1]) * (LONG8) y41
2137 - (LONG8) (c[(s3 << 1) + 2] - c[(s1 << 1) + 2]) * (LONG8) x41;
2138 if (det2 > 0 && det3 < 0)
2139 {
2140 /* le quadrilataire est convexe on le retourne */
2141 /* update des sommets */
2142 /* ------------------------- */
2143 nu[i12 + t1 * 6] = s4;
2144 nu[i22 + t2 * 6] = s1;
2145 /* update du pointeur suivant */
2146 /* ----------------------------------- */
2147 pslst = lst[ptlst * 3 + 1];
2148 if (pslst > 0)
2149 {
2150 aas = lst[pslst * 3 + 3];
2151 if (aas == i22 + 3)
2152 {
2153 lst[pslst * 3 + 2] = t1;
2154 lst[pslst * 3 + 3] = i11 + 3;
2155 }
2156 }
2157 /* update des aretes a1,a2 */
2158 /* ------------------------------- */
2159 tt1 = nu[i22 + 3 + t2 * 6];
2160 nu[a1 + t1 * 6] = tt1;
2161 if (tt1 > 0)
2162 {
2163 tt = tt1 / 8;
2164 aa = tt1 - (tt << 3);
2165 nu[aa + tt * 6] = a1 + (t1 << 3);
2166 }
2167 else if (tt1 != nothing)
2168 {
2169 nu[-tt1 * 6 + 2] = a1 + (t1 << 3);
2170 }
2171 tt1 = nu[i12 + 3 + t1 * 6];
2172 nu[a2 + t2 * 6] = tt1;
2173 if (tt1 > 0)
2174 {
2175 tt = tt1 / 8;
2176 aa = tt1 - (tt << 3);
2177 nu[aa + tt * 6] = a2 + (t2 << 3);
2178 }
2179 else if (tt1 != nothing)
2180 {
2181 nu[-tt1 * 6 + 2] = a2 + (t2 << 3);
2182 }
2183 nu[i12 + 3 + t1 * 6] = i22 + 3 + (t2 << 3);
2184 nu[i22 + 3 + t2 * 6] = i12 + 3 + (t1 << 3);
2185
2186 det1 = (LONG8) (c[(s1 << 1) + 1] - c[(*ss1 << 1) + 1]) * (LONG8) y
2187 - (LONG8) (c[(s1 << 1) + 2] - c[(*ss1 << 1) + 2]) * (LONG8) x;
2188
2189 det4 = (LONG8) (c[(s4 << 1) + 1] - c[(*ss1 << 1) + 1]) * (LONG8) y
2190 - (LONG8) (c[(s4 << 1) + 2] - c[(*ss1 << 1) + 2]) * (LONG8) x;
2191 if (det1 < 0 && det4 > 0)
2192 {
2193 /* le sommets s4 est dans omega */
2194 lst[ptlst * 3 + 2] = t2;
2195 lst[ptlst * 3 + 3] = i22 + 3;
2196 }
2197 else if (det1 > 0 && det4 < 0)
2198 {
2199 /* le sommets s1 est dans omega */
2200 lst[ptlst * 3 + 2] = t1;
2201 lst[ptlst * 3 + 3] = i12 + 3;
2202 }
2203 else
2204 {
2205 if (pplst == 0)
2206 {
2207 ttlst = lst[ptlst * 3 + 1];
2208 ptlst = ttlst;
2209 }
2210 else
2211 {
2212 ptlst = lst[ptlst * 3 + 1];
2213 lst[pplst * 3 + 1] = ptlst;
2214 }
2215 goto L30;
2216 }
2217 }
2218 pplst = ptlst;
2219 ptlst = lst[ptlst * 3 + 1];
2220 goto L30;
2221 }
2222 if (ttlst != 0)
2223 {
2224 goto L20;
2225 }
2226 nu[i12 + 3 + t1 * 6] = nothing;
2227 nu[i22 + 3 + t2 * 6] = nothing;
2228 *t = t2;
2229 *ta = t1;
2230 i_1 = *nbac;
2231 for (i = 1; i <= i_1; ++i)
2232 {
2233 mshopt_ (&c[3], &nu[7], &lst[i * 3 + 2], 4L, err);
2234 mshopt_ (&c[3], &nu[7], &lst[i * 3 + 2], 5L, err);
2235 mshopt_ (&c[3], &nu[7], &lst[i * 3 + 2], 6L, err);
2236 }
2237 return 1;
2238 } /* mshfr2_ */
2239
2240
2241 }
2242