src/tool/Utilitaire/utb3b0.F

   1       subroutine utb3b0 ( hetnoe, coonoe,
   2      >                    numcoi, coinpt, coinnn,
   3      >                    somare,
   4      >                    hettri, aretri, voltri,
   5      >                    np2are,
   6      >                    nbpbco, mess08, mess54,
   7      >                    ulbila, ulsort, langue, codret )
   8 c ______________________________________________________________________
   9 c
  10 c                             H O M A R D
  11 c
  12 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
  13 c
  14 c Version originale enregistree le 18 juin 1996 sous le numero 96036
  15 c aupres des huissiers de justice Simart et Lavoir a Clamart
  16 c Version 11.2 enregistree le 13 fevrier 2015 sous le numero 2015/014
  17 c aupres des huissiers de justice
  18 c Lavoir, Silinski & Cherqui-Abrahmi a Clamart
  19 c
  20 c    HOMARD est une marque deposee d'Electricite de France
  21 c
  22 c Copyright EDF 1996
  23 c Copyright EDF 1998
  24 c Copyright EDF 2002
  25 c Copyright EDF 2020
  26 c ______________________________________________________________________
  27 c
  28 c    UTilitaire - Bilan - option 3 - phase B0
  29 c    --           -              -         --
  30 c ______________________________________________________________________
  31 c
  32 c but : controle l'interpenetration des triangles.
  33 c ______________________________________________________________________
  34 c .        .     .        .                                            .
  35 c .  nom   . e/s . taille .           description                      .
  36 c .____________________________________________________________________.
  37 c . hetnoe . e   . nbnoto . historique de l'etat des noeuds            .
  38 c . coonoe . e   . nbnoto . coordonnees des noeuds                     .
  39 c .        .     . * sdim .                                            .
  40 c . numcoi . e   . nbnoto . numero de la coincidence du noeud          .
  41 c . coinpt . e   .   *    . pointeur de la i-eme coincidence dans coinn.
  42 c . coinnn . e   .   *    . liste des noeuds coincidents               .
  43 c . somare . e   .2*nbarto. numeros des extremites d'arete             .
  44 c . hettri . e   . nbtrto . historique de l'etat des triangles         .
  45 c . aretri . e   .nbtrto*3. numeros des 3 aretes des triangles         .
  46 c . voltri . e   .2*nbtrto. numeros des 2 volumes par triangle         .
  47 c .        .     .        . voltri(i,k) definit le i-eme voisin de k   .
  48 c .        .     .        .   0 : pas de voisin                        .
  49 c .        .     .        . j>0 : tetraedre j                          .
  50 c .        .     .        . j<0 : pyramide/pentaedre dans pypetr(1/2,j).
  51 c . np2are . e   . nbarto . noeud milieux des aretes                   .
  52 c . nbpbco . es  .  -1:7  . nombre de problemes de coincidences        .
  53 c . mess54 . e   .nblang,*. messages                                   .
  54 c . ulbila . e   .   1    . unite logique d'ecriture du bilan          .
  55 c . ulsort . e   .   1    . unite logique de la sortie generale        .
  56 c . langue . e   .    1   . langue des messages                        .
  57 c .        .     .        . 1 : francais, 2 : anglais                  .
  58 c . codret .  s  .    1   . code de retour des modules                 .
  59 c .        .     .        . 0 : pas de probleme                        .
  60 c .        .     .        . 1 : probleme                               .
  61 c .____________________________________________________________________.
  62 c
  63 c====
  64 c 0. declarations et dimensionnement
  65 c====
  66 c
  67 c 0.1. ==> generalites
  68 c
  69       implicit none
  70       save
  71 c
  72       character*6 nompro
  73       parameter ( nompro = 'UTB3B0' )
  74 c
  75       integer typenh
  76       parameter ( typenh = 2 )
  77 c
  78 #include "nblang.h"
  79 c
  80 c 0.2. ==> communs
  81 c
  82 #include "nombno.h"
  83 #include "nombar.h"
  84 #include "nombtr.h"
  85 #include "nombte.h"
  86 #include "nombpy.h"
  87 #include "nombpe.h"
  88 #include "envca1.h"
  89 #include "precis.h"
  90 #include "impr02.h"
  91 c
  92 c 0.3. ==> arguments
  93 c
  94       double precision coonoe(nbnoto,sdim)
  95 c
  96       integer hetnoe(nbnoto)
  97       integer numcoi(nbnoto), coinpt(*), coinnn(*)
  98       integer somare(2,nbarto)
  99       integer aretri(nbtrto,3), hettri(nbtrto), voltri(2,nbtrto)
 100       integer np2are(nbarto)
 101       integer nbpbco(-1:7)
 102 c
 103       character*08 mess08(nblang,*)
 104       character*54 mess54(nblang,*)
 105 c
 106       integer ulbila
 107       integer ulsort, langue, codret
 108 c
 109 c 0.4. ==> variables locales
 110 c
 111       integer iaux, jaux
 112       integer letria, lenoeu
 113       integer nucoin, ptcoin, ptcode, ptcofi
 114       integer sommet(6), nbsomm
 115       integer a1, a2, a3
 116       integer sa1a2, sa2a3, sa3a1
 117 #ifdef _DEBUG_HOMARD_
 118       integer glop
 119 #endif
 120 c
 121       double precision v1(3), v2(3), v3(3), v4(3), v6(4,3), vn(3)
 122       double precision xmax, xmin, ymax, ymin, zmax, zmin
 123       double precision prosca
 124       double precision daux1, daux2
 125 c
 126       logical logaux(7)
 127 c
 128       integer nbmess
 129       parameter (nbmess = 10 )
 130       character*80 texte(nblang,nbmess)
 131 c
 132 c 0.5. ==> initialisations
 133 c ______________________________________________________________________
 134 c
 135 c====
 136 c 1. initialisations
 137 c====
 138 c
 139 c 1.1. ==> messages
 140 c
 141 #include "impr01.h"
 142 c
 143 #ifdef _DEBUG_HOMARD_
 144       write (ulsort,texte(langue,1)) 'Entree', nompro
 145       call dmflsh (iaux)
 146 #endif
 147 c
 148 #include "utb300.h"
 149 c
 150 #include "utb301.h"
 151 c
 152 c 1.2. ==> constantes
 153 c
 154       codret = 0
 155 c
 156       if ( degre.eq.1 ) then
 157         nbsomm = 3
 158       else
 159         nbsomm = 6
 160       endif
 161 c
 162       if ( sdim.eq.2 ) then
 163         v1(3) = 0.d0
 164         v2(3) = 0.d0
 165         v3(3) = 0.d0
 166         vn(3) = 0.d0
 167       endif
 168 c
 169       if ( nbteto.eq.0 .and. nbpyto.eq.0 .and. nbpeto.eq.0 ) then
 170         logaux(2) = .true.
 171       else
 172         logaux(2) = .false.
 173       endif
 174 c
 175 c====
 176 c 2. controle de la non-interpenetration des triangles
 177 c    remarque :
 178 c    La verification est sujette a caution car le test sur la
 179 c    coplanarite est un test sur une egalite de reels ...
 180 c====
 181 c
 182       do 20 , letria = 1 , nbtrto
 183 c
 184 #ifdef _DEBUG_HOMARD_
 185         if ( letria.lt.0 ) then
 186           glop = 1
 187         else
 188           glop = 0
 189         endif
 190 #endif
 191 c
 192 c 2.1. ==> Filtrage
 193 c    1. on ne s'interesse qu'aux actifs car les autres sont
 194 c    censes avoir ete controles aux iterations anterieures
 195 c    2. on ne s'interesse qu'aux triangles de region 2D, car ceux qui
 196 c    bordent des volumes seront vus par la suite.
 197 c
 198         if ( mod(hettri(letria),10).eq.0 ) then
 199           if ( logaux(2) ) then
 200             logaux(1) = .true.
 201           else
 202             if ( voltri(1,letria).eq.0 ) then
 203               logaux(1) = .true.
 204             else
 205               logaux(1) = .false.
 206             endif
 207           endif
 208         else
 209           logaux(1) = .false.
 210         endif
 211 c
 212         if ( logaux(1) ) then
 213 c
 214           if ( nbpbco(typenh).eq.-1 ) then
 215 #ifdef _DEBUG_HOMARD_
 216       write (ulsort,texte(langue,4)) mess14(langue,3,typenh)
 217 #endif
 218             nbpbco(typenh) = 0
 219           endif
 220 c
 221 c 2.3. ==> les aretes et les sommets de ce triangle
 222 c
 223 #ifdef _DEBUG_HOMARD_
 224         if ( glop.ne.0 ) then
 225         write (ulsort,*) '.. ', mess14(langue,2,typenh), letria
 226         endif
 227 #endif
 228 c
 229           a1  = aretri(letria,1)
 230           a2  = aretri(letria,2)
 231           a3  = aretri(letria,3)
 232 c
 233           call utsotr ( somare, a1, a2, a3,
 234      >                  sa1a2, sa2a3, sa3a1 )
 235 c
 236           sommet(1) = sa1a2
 237           sommet(2) = sa2a3
 238           sommet(3) = sa3a1
 239 c
 240           if ( degre.eq.2 ) then
 241             sommet(4) = np2are(a1)
 242             sommet(5) = np2are(a2)
 243             sommet(6) = np2are(a3)
 244           endif
 245 c
 246           v1(1) = coonoe(sommet(1),1)
 247           v1(2) = coonoe(sommet(1),2)
 248           if ( sdim.eq.3 ) then
 249             v1(3) = coonoe(sommet(1),3)
 250           endif
 251 c
 252           v2(1) = coonoe(sommet(2),1)
 253           v2(2) = coonoe(sommet(2),2)
 254           if ( sdim.eq.3 ) then
 255             v2(3) = coonoe(sommet(2),3)
 256           endif
 257 c
 258           v3(1) = coonoe(sommet(3),1)
 259           v3(2) = coonoe(sommet(3),2)
 260           if ( sdim.eq.3 ) then
 261             v3(3) = coonoe(sommet(3),3)
 262           endif
 263 c
 264           xmin = min(v1(1),v2(1),v3(1))
 265           xmax = max(v1(1),v2(1),v3(1))
 266           ymin = min(v1(2),v2(2),v3(2))
 267           ymax = max(v1(2),v2(2),v3(2))
 268           zmin = min(v1(3),v2(3),v3(3))
 269           zmax = max(v1(3),v2(3),v3(3))
 270 c
 271 c         v6(1,.) est le produit vectoriel n1n2 x n1n3.
 272 c
 273           v6(1,1) = (v2(2)-v1(2)) * (v3(3)-v1(3))
 274      >            - (v2(3)-v1(3)) * (v3(2)-v1(2))
 275           v6(1,2) = (v2(3)-v1(3)) * (v3(1)-v1(1))
 276      >            - (v2(1)-v1(1)) * (v3(3)-v1(3))
 277           v6(1,3) = (v2(1)-v1(1)) * (v3(2)-v1(2))
 278      >            - (v2(2)-v1(2)) * (v3(1)-v1(1))
 279 c
 280 c         v6(2,.) est le produit vectoriel n2n3 x n2n1.
 281 c
 282           v6(2,1) = (v3(2)-v2(2)) * (v1(3)-v2(3))
 283      >            - (v3(3)-v2(3)) * (v1(2)-v2(2))
 284           v6(2,2) = (v3(3)-v2(3)) * (v1(1)-v2(1))
 285      >            - (v3(1)-v2(1)) * (v1(3)-v2(3))
 286           v6(2,3) = (v3(1)-v2(1)) * (v1(2)-v2(2))
 287      >            - (v3(2)-v2(2)) * (v1(1)-v2(1))
 288 c
 289 c         v6(3,.) est le produit vectoriel n3n1 x n3n2.
 290 c
 291           v6(3,1) = (v1(2)-v3(2)) * (v2(3)-v3(3))
 292      >            - (v1(3)-v3(3)) * (v2(2)-v3(2))
 293           v6(3,2) = (v1(3)-v3(3)) * (v2(1)-v3(1))
 294      >            - (v1(1)-v3(1)) * (v2(3)-v3(3))
 295           v6(3,3) = (v1(1)-v3(1)) * (v2(2)-v3(2))
 296      >            - (v1(2)-v3(2)) * (v2(1)-v3(1))
 297 c
 298 c 2.3. ==> on passe en revue tous les autres sommets qui ne sont pas des
 299 c          sommets isoles.
 300 c       . on ne s'interesse qu'a ceux qui sont strictement dans le
 301 c         parallelepide enveloppe du triangle
 302 c       . on commence par eliminer les trois noeuds du triangle
 303 c       . ensuite, on elimine les noeuds coincidents
 304 c       . on recherche si le noeud est a l'interieur du triangle
 305 c
 306 c         Remarque hyper importante : il ne faut faire les affectations
 307 c         de vn(2) et vn(3) que si c'est utile car elles coutent
 308 c         tres cheres (30% du temps total !)
 309 c         En revanche, inutile de deplier davantage les tests
 310 c         Remarque hyper importante : il vaut mieux mettre en dernier
 311 c         le test sur l'identite de lenoeu avec les noeuds du triangle
 312 c         car on gagne aussi 40% !
 313 c
 314           do 23 , lenoeu = numip1, numap1
 315 c
 316             logaux(7) = .false.
 317 c
 318             vn(1) = coonoe(lenoeu,1)
 319             if ( vn(1).ge.xmin .and. vn(1).le.xmax ) then
 320               vn(2) = coonoe(lenoeu,2)
 321               if ( vn(2).ge.ymin .and. vn(2).le.ymax ) then
 322                 if ( sdim.eq.3 ) then
 323                   vn(3) = coonoe(lenoeu,3)
 324                   if ( vn(3).ge.zmin .and. vn(3).le.zmax ) then
 325                     logaux(7) = .true.
 326                   endif
 327                 else
 328                   logaux(7) = .true.
 329                 endif
 330               endif
 331             endif
 332 c
 333             if ( logaux(7) ) then
 334               do 231 , iaux = 1 , nbsomm
 335                 if ( lenoeu.eq.sommet(iaux) ) then
 336                   goto 23
 337                 endif
 338   231         continue
 339             endif
 340 c
 341 c 2.3.2. ==> le noeud est-il coincident avec un des sommets ?
 342 c
 343             if ( logaux(7) ) then
 344 c
 345               if ( nbpbco(-1).gt.0 ) then
 346 c
 347                 nucoin = numcoi(lenoeu)
 348 c
 349                 if ( nucoin.ne.0 ) then
 350 c
 351                   ptcode = coinpt(nucoin)+1
 352                   ptcofi = coinpt(nucoin+1)
 353                   do 232 , ptcoin = ptcode, ptcofi
 354                     jaux = coinnn(ptcoin)
 355                     do 2321 , iaux = 1 , nbsomm
 356                       if ( jaux.eq.sommet(iaux) ) then
 357                         goto 23
 358                       endif
 359  2321               continue
 360   232             continue
 361 c
 362                 endif
 363 c
 364               endif
 365 c
 366             endif
 367 c
 368 c 2.3.3. ==> exclusivement les noeuds p1
 369 c
 370             if ( logaux(7) ) then
 371 c
 372               if ( hetnoe(lenoeu).ne.1  .and.
 373      >             hetnoe(lenoeu).ne.21 .and.
 374      >             hetnoe(lenoeu).ne.51 ) then
 375 c
 376                 logaux(7) = .false.
 377 c
 378               endif
 379 c
 380             endif
 381 c
 382 c 2.3.2. ==> le noeud est-il dans le plan du triangle ?
 383 c            on calcule les deux vecteurs normaux aux triangles
 384 c            (s1,s2,s3) et (s1,s2,n) : vecteurs v6 et v4
 385 c            on cherche s'ils sont paralleles, c'est-a-dire de produit
 386 c            vectoriel nul
 387 c
 388             if ( logaux(7) ) then
 389 c
 390 c            v4 est le produit vectoriel n1n2 x n1n.
 391 c
 392               v4(1) = (v2(2)-v1(2)) * (vn(3)-v1(3))
 393      >              - (v2(3)-v1(3)) * (vn(2)-v1(2))
 394               v4(2) = (v2(3)-v1(3)) * (vn(1)-v1(1))
 395      >              - (v2(1)-v1(1)) * (vn(3)-v1(3))
 396               v4(3) = (v2(1)-v1(1)) * (vn(2)-v1(2))
 397      >              - (v2(2)-v1(2)) * (vn(1)-v1(1))
 398 c
 399 c            daux2 represente la norme du produit vectoriel v4 x v6(1,.)
 400 c
 401               daux2 = abs ( v6(1,2) * v4(3) - v6(1,3) * v4(2) )
 402      >              + abs ( v6(1,3) * v4(1) - v6(1,1) * v4(3) )
 403      >              + abs ( v6(1,1) * v4(2) - v6(1,2) * v4(1) )
 404 c
 405 c 2.3.3. ==> dans ce plan, n est-il dedans ?
 406 c            cela est vrai si le noeud et un sommet sont du meme cote
 407 c            de l'arete formee par les deux autres sommets. pour cela,
 408 c            on regarde si les produits vectoriels (ab,ac) et (ab,an)
 409 c            sont de meme orientation pour les trois permutations
 410 c            circulaires sur (a,b,c), c'est-a-dire si le produit
 411 c            scalaire des deux produits vectoriels est positif.
 412 c            on teste le caractere strictement positif du produit
 413 c            scalaire, pour pouvoir pieger les cas ou le noeud est sur
 414 c            une arete.
 415 c
 416               if ( daux2.le.epsima ) then
 417 c
 418 c 2.3.3.1. ==> critere absolu ou relatif
 419 c
 420                 daux1 = 0.d0
 421 c
 422 c 2.3.3.2. ==> arete (s1,s2)
 423 c
 424                 prosca = v4(1)*v6(1,1) + v4(2)*v6(1,2) + v4(3)*v6(1,3)
 425 c
 426                 if ( prosca.lt.daux1 ) then
 427                   goto 23
 428                 endif
 429 c
 430 c 2.3.3.3. ==> arete (s2,s3)
 431 c
 432 c            v4 est le produit vectoriel n2n3 x s2n.
 433 c
 434                 v4(1) = (v3(2)-v2(2)) * (vn(3)-v2(3))
 435      >                - (v3(3)-v2(3)) * (vn(2)-v2(2))
 436                 v4(2) = (v3(3)-v2(3)) * (vn(1)-v2(1))
 437      >                - (v3(1)-v2(1)) * (vn(3)-v2(3))
 438                 v4(3) = (v3(1)-v2(1)) * (vn(2)-v2(2))
 439      >                - (v3(2)-v2(2)) * (vn(1)-v2(1))
 440 c
 441                 prosca = v4(1)*v6(2,1) + v4(2)*v6(2,2) + v4(3)*v6(2,3)
 442 c
 443                 if ( prosca.lt.daux1 ) then
 444                   goto 23
 445                 endif
 446 c
 447 c 2.3.3.4. ==> arete (s3,s1)
 448 c
 449 c            v4 est le produit vectoriel n3n1 x s3n.
 450 c
 451                 v4(1) = (v1(2)-v3(2)) * (vn(3)-v3(3))
 452      >                - (v1(3)-v3(3)) * (vn(2)-v3(2))
 453                 v4(2) = (v1(3)-v3(3)) * (vn(1)-v3(1))
 454      >                - (v1(1)-v3(1)) * (vn(3)-v3(3))
 455                 v4(3) = (v1(1)-v3(1)) * (vn(2)-v3(2))
 456      >                - (v1(2)-v3(2)) * (vn(1)-v3(1))
 457 c
 458                 prosca = v4(1)*v6(3,1) + v4(2)*v6(3,2) + v4(3)*v6(3,3)
 459 c
 460                 if ( prosca.lt.daux1 ) then
 461                   goto 23
 462                 endif
 463 c
 464 c 2.3.5. ==> si logaux(7) est encore vrai, c'est que le noeud est
 465 c            a l'interieur du triangle ... malaise ...
 466 c
 467                 iaux = letria
 468 c
 469 #include "utb302.h"
 470 c
 471                 if ( sdim.eq.2 ) then
 472                   write (ulbila,14202) sommet(1), v1(1), v1(2)
 473                   write (ulbila,14202) sommet(2), v2(1), v2(2)
 474                   write (ulbila,14202) sommet(3), v3(1), v3(2)
 475                 else
 476                   write (ulbila,14203) sommet(1), v1(1), v1(2), v1(3)
 477                   write (ulbila,14203) sommet(2), v2(1), v2(2), v2(3)
 478                   write (ulbila,14203) sommet(3), v3(1), v3(2), v3(3)
 479                 endif
 480 c
 481                 write (ulbila,10200)
 482 c
 483               endif
 484 c
 485             endif
 486 c
 487    23     continue
 488 c
 489         endif
 490 c
 491    20 continue
 492 c
 493       end