src/tool/Utilitaire/utb3c0.F

   1       subroutine utb3c0 ( hetnoe, coonoe,
   2      >                    numcoi, coinpt, coinnn,
   3      >                    somare,
   4      >                    hetqua, arequa, volqua,
   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 C0
  29 c    --           -              -         --
  30 c ______________________________________________________________________
  31 c
  32 c but : controle l'interpenetration des quadrangles.
  33 c      remarque : cela suppose que les quadrangles sont plans
  34 c ______________________________________________________________________
  35 c .        .     .        .                                            .
  36 c .  nom   . e/s . taille .           description                      .
  37 c .____________________________________________________________________.
  38 c . hetnoe . e   . nbnoto . historique de l'etat des noeuds            .
  39 c . coonoe . e   . nbnoto . coordonnees des noeuds                     .
  40 c .        .     . * sdim .                                            .
  41 c . numcoi . e   . nbnoto . numero de la coincidence du noeud          .
  42 c . coinpt . e   .   *    . pointeur de la i-eme coincidence dans coinn.
  43 c . coinnn . e   .   *    . liste des noeuds coincidents               .
  44 c . somare . e   .2*nbarto. numeros des extremites d'arete             .
  45 c . hetqua . e   . nbtrto . historique de l'etat des quadrangles       .
  46 c . arequa . e   .nbtrto*4. numeros des 3 aretes des quadrangles       .
  47 c . volqua . e   .2*nbquto. numeros des 2 volumes par quadrangle       .
  48 c .        .     .        . volqua(i,k) definit le i-eme voisin de k   .
  49 c .        .     .        .   0 : pas de voisin                        .
  50 c .        .     .        . j>0 : hexaedre j                           .
  51 c .        .     .        . j<0 : pyramide/pentaedre dans pypequ(1/2,j).
  52 c . np2are . e   . nbarto . noeud milieux des aretes                   .
  53 c . nbpbco . es  .  -1:7  . nombre de problemes de coincidences        .
  54 c . mess54 . e   .nblang,*. messages                                   .
  55 c . ulbila . e   .   1    . unite logique d'ecriture du bilan          .
  56 c . ulsort . e   .   1    . unite logique de la sortie generale        .
  57 c . langue . e   .    1   . langue des messages                        .
  58 c .        .     .        . 1 : francais, 2 : anglais                  .
  59 c . codret .  s  .    1   . code de retour des modules                 .
  60 c .        .     .        . 0 : pas de probleme                        .
  61 c .        .     .        . 1 : probleme                               .
  62 c .____________________________________________________________________.
  63 c
  64 c====
  65 c 0. declarations et dimensionnement
  66 c====
  67 c
  68 c 0.1. ==> generalites
  69 c
  70       implicit none
  71       save
  72 c
  73       character*6 nompro
  74       parameter ( nompro = 'UTB3C0' )
  75 c
  76       integer typenh
  77       parameter ( typenh = 4 )
  78 c
  79 #include "nblang.h"
  80 c
  81 c 0.2. ==> communs
  82 c
  83 #include "nombno.h"
  84 #include "nombar.h"
  85 #include "nombqu.h"
  86 #include "nombhe.h"
  87 #include "nombpy.h"
  88 #include "nombpe.h"
  89 #include "precis.h"
  90 #include "envca1.h"
  91 #include "impr02.h"
  92 c
  93 c 0.3. ==> arguments
  94 c
  95       double precision coonoe(nbnoto,sdim)
  96 c
  97       integer hetnoe(nbnoto)
  98       integer numcoi(nbnoto), coinpt(*), coinnn(*)
  99       integer somare(2,nbarto)
 100       integer hetqua(nbquto), arequa(nbquto,4), volqua(2,nbquto)
 101       integer np2are(nbarto)
 102       integer nbpbco(-1:7)
 103 c
 104       character*08 mess08(nblang,*)
 105       character*54 mess54(nblang,*)
 106 c
 107       integer ulbila
 108       integer ulsort, langue, codret
 109 c
 110 c 0.4. ==> variables locales
 111 c
 112       integer iaux, jaux
 113       integer lequad, lenoeu
 114       integer nucoin, ptcoin, ptcode, ptcofi
 115       integer sommet(8), nbsomm
 116       integer a1, a2, a3, a4
 117       integer sa1a2, sa2a3, sa3a4, sa4a1
 118 #ifdef _DEBUG_HOMARD_
 119       integer glop
 120 #endif
 121 c
 122       double precision v1(3), v2(3), v3(3), v4(3), vn(3)
 123       double precision v5(2,3), v6(4,3)
 124       double precision v12(3)
 125       double precision xmax, xmin, ymax, ymin, zmax, zmin
 126       double precision prosca
 127       double precision daux1, daux2
 128 c
 129       logical logaux(7)
 130 c
 131       integer nbmess
 132       parameter (nbmess = 10 )
 133       character*80 texte(nblang,nbmess)
 134 c
 135 c 0.5. ==> initialisations
 136 c ______________________________________________________________________
 137 c
 138 c====
 139 c 1. initialisations
 140 c====
 141 c
 142 c 1.1. ==> messages
 143 c
 144 #include "impr01.h"
 145 c
 146 #ifdef _DEBUG_HOMARD_
 147       write (ulsort,texte(langue,1)) 'Entree', nompro
 148       call dmflsh (iaux)
 149 #endif
 150 c
 151 #include "utb300.h"
 152 c
 153 #include "utb301.h"
 154 c
 155 c 1.2. ==> constantes
 156 c
 157       codret = 0
 158 c
 159       if ( degre.eq.1 ) then
 160         nbsomm = 4
 161       else
 162         nbsomm = 8
 163       endif
 164 c
 165       if ( sdim.eq.2 ) then
 166         v1(3) = 0.d0
 167         v2(3) = 0.d0
 168         v3(3) = 0.d0
 169         v4(3) = 0.d0
 170         vn(3) = 0.d0
 171       endif
 172 c
 173       if ( nbheto.eq.0 .and. nbpyto.eq.0 .and. nbpeto.eq.0 ) then
 174         logaux(2) = .true.
 175       else
 176         logaux(2) = .false.
 177       endif
 178 c
 179 c====
 180 c 2. controle de la non-interpenetration des quadrangles
 181 c    remarques :
 182 c    3. La verification est sujette a caution car le test sur la
 183 c    coplanarite est un test sur une egalite de reels ...
 184 c====
 185 cgn      call gtdems (112)
 186 c
 187       do 20 , lequad = 1 , nbquto
 188 c
 189 #ifdef _DEBUG_HOMARD_
 190         if ( lequad.lt.0 ) then
 191           glop = 1
 192         else
 193           glop = 0
 194         endif
 195 #endif
 196 c
 197 c 2.1. ==> Filtrage
 198 c    1. on ne s'interesse qu'aux actifs car les autres sont
 199 c    censes avoir ete controles aux iterations anterieures
 200 c    2. on ne s'interesse qu'aux quadrangles de region 2D, car ceux qui
 201 c    bordent des volumes seront vus par la suite.
 202 c
 203         if ( mod(hetqua(lequad),100).eq.0 ) then
 204           if ( logaux(2) ) then
 205             logaux(1) = .true.
 206           else
 207             if ( volqua(1,lequad).eq.0 ) then
 208               logaux(1) = .true.
 209             else
 210               logaux(1) = .false.
 211             endif
 212           endif
 213         else
 214           logaux(1) = .false.
 215         endif
 216 c
 217         if ( logaux(1) ) then
 218 c
 219           if ( nbpbco(typenh).eq.-1 ) then
 220 #ifdef _DEBUG_HOMARD_
 221       write (ulsort,texte(langue,4)) mess14(langue,3,typenh)
 222 #endif
 223             nbpbco(typenh) = 0
 224           endif
 225 c
 226 c 2.3. ==> les aretes et les sommets de ce quadrangle
 227 c
 228 #ifdef _DEBUG_HOMARD_
 229         if ( glop.ne.0 ) then
 230         write (ulsort,*) '.. ', mess14(langue,2,typenh), lequad
 231         endif
 232 #endif
 233 c
 234           a1  = arequa(lequad,1)
 235           a2  = arequa(lequad,2)
 236           a3  = arequa(lequad,3)
 237           a4  = arequa(lequad,4)
 238 c
 239           call utsoqu ( somare, a1, a2, a3, a4,
 240      >                  sa1a2, sa2a3, sa3a4, sa4a1 )
 241 c
 242           sommet(1) = sa1a2
 243           sommet(2) = sa2a3
 244           sommet(3) = sa3a4
 245           sommet(4) = sa4a1
 246 c
 247           if ( degre.eq.2 ) then
 248             sommet(5) = np2are(a1)
 249             sommet(6) = np2are(a2)
 250             sommet(7) = np2are(a3)
 251             sommet(8) = np2are(a4)
 252           endif
 253 c
 254 c 2.3. ==> le parallelepipede enveloppe
 255 c
 256           v1(1) = coonoe(sommet(1),1)
 257           v1(2) = coonoe(sommet(1),2)
 258           if ( sdim.eq.3 ) then
 259             v1(3) = coonoe(sommet(1),3)
 260           endif
 261 c
 262           v2(1) = coonoe(sommet(2),1)
 263           v2(2) = coonoe(sommet(2),2)
 264           if ( sdim.eq.3 ) then
 265             v2(3) = coonoe(sommet(2),3)
 266           endif
 267 c
 268           v3(1) = coonoe(sommet(3),1)
 269           v3(2) = coonoe(sommet(3),2)
 270           if ( sdim.eq.3 ) then
 271             v3(3) = coonoe(sommet(3),3)
 272           endif
 273 c
 274           v4(1) = coonoe(sommet(4),1)
 275           v4(2) = coonoe(sommet(4),2)
 276           if ( sdim.eq.3 ) then
 277             v4(3) = coonoe(sommet(4),3)
 278           endif
 279 c
 280           xmin = min(v1(1),v2(1),v3(1),v4(1))
 281           xmax = max(v1(1),v2(1),v3(1),v4(1))
 282           ymin = min(v1(2),v2(2),v3(2),v4(2))
 283           ymax = max(v1(2),v2(2),v3(2),v4(2))
 284           zmin = min(v1(3),v2(3),v3(3),v4(3))
 285           zmax = max(v1(3),v2(3),v3(3),v4(3))
 286 cgn          if ( lequad.eq.1270 ) then
 287 cgn        print *,xmin,xmax,ymin,ymax,zmin,zmax
 288 cgn          endif
 289 c
 290           logaux(3) = .true.
 291           logaux(4) = .true.
 292 c
 293 c 2.3. ==> on verifie que le quadrangle est plan
 294 c          on calcule les deux vecteurs normaux aux deux triangles
 295 c          inclus dans le quadrangle, (s1,s2,s4) et (s3,s2,s4).
 296 c          on cherche s'ils sont paralleles, c'est-a-dire de produit
 297 c          vectoriel nul
 298 c
 299 c         v6(1,.) est le produit vectoriel n1n2 x n1n4.
 300 c
 301           v6(1,1) = (v2(2)-v1(2)) * (v4(3)-v1(3))
 302      >            - (v2(3)-v1(3)) * (v4(2)-v1(2))
 303           v6(1,2) = (v2(3)-v1(3)) * (v4(1)-v1(1))
 304      >            - (v2(1)-v1(1)) * (v4(3)-v1(3))
 305           v6(1,3) = (v2(1)-v1(1)) * (v4(2)-v1(2))
 306      >            - (v2(2)-v1(2)) * (v4(1)-v1(1))
 307 c
 308 c         v6(3,.) est le produit vectoriel n3n4 x n3n2.
 309 c
 310           v6(3,1) = (v4(2)-v3(2)) * (v2(3)-v3(3))
 311      >            - (v4(3)-v3(3)) * (v2(2)-v3(2))
 312           v6(3,2) = (v4(3)-v3(3)) * (v2(1)-v3(1))
 313      >            - (v4(1)-v3(1)) * (v2(3)-v3(3))
 314           v6(3,3) = (v4(1)-v3(1)) * (v2(2)-v3(2))
 315      >            - (v4(2)-v3(2)) * (v2(1)-v3(1))
 316 c
 317 c         daux2 est le carre de la norme du produit vectoriel
 318 c         v6(1,.) x v6(3,.)
 319 c
 320           daux2 = abs ( v6(3,2)*v6(1,3) - v6(3,3)*v6(1,2) )
 321      >          + abs ( v6(3,3)*v6(1,1) - v6(3,1)*v6(1,3) )
 322      >          + abs ( v6(3,1)*v6(1,2) - v6(3,2)*v6(1,1) )
 323 c
 324           if ( daux2.gt.epsima ) then
 325 cgn            print *,'Quadrangle ',lequad,' non plan.'
 326             goto 20
 327           endif
 328 c
 329 c 2.3. ==> on passe en revue tous les autres sommets qui ne sont pas des
 330 c          sommets isoles.
 331 c       . on ne s'interesse qu'a ceux qui sont strictement dans le
 332 c         parallelepide enveloppe du quadrangle
 333 c       . on commence par eliminer les quatre noeuds du quadrangle
 334 c       . ensuite, on elimine les noeuds coincidents
 335 c       . on elimine les noeuds qui ne sont pas dans le plan du
 336 c         quadrangle
 337 c       . on recherche si le noeud est a l'interieur du quadrangle
 338 c
 339 c         Remarque hyper importante : il ne faut faire les affectations
 340 c         de vn(2) et vn(3) que si c'est utile car elles coutent
 341 c         tres cheres (30% du temps total !)
 342 c         Remarque hyper importante : il vaut mieux mettre en dernier
 343 c         le test sur l'identite de lenoeu avec les noeuds du quadrangle
 344 c         car on gagne aussi 40% !
 345 c         En revanche, inutile de deplier davantage les tests
 346 c
 347           do 23 , lenoeu = numip1, numap1
 348 c
 349             logaux(7) = .false.
 350 c
 351             vn(1) = coonoe(lenoeu,1)
 352             if ( vn(1).ge.xmin .and. vn(1).le.xmax ) then
 353               vn(2) = coonoe(lenoeu,2)
 354               if ( vn(2).ge.ymin .and. vn(2).le.ymax ) then
 355                 if ( sdim.eq.3 ) then
 356                   vn(3) = coonoe(lenoeu,3)
 357                   if ( vn(3).ge.zmin .and. vn(3).le.zmax ) then
 358                     logaux(7) = .true.
 359                   endif
 360                 else
 361                   logaux(7) = .true.
 362                 endif
 363               endif
 364             endif
 365 c
 366             if ( logaux(7) ) then
 367               do 231 , iaux = 1 , nbsomm
 368                 if ( lenoeu.eq.sommet(iaux) ) then
 369                   goto 23
 370                 endif
 371   231         continue
 372             endif
 373 c
 374 c 2.3.2. ==> le noeud est-il coincident avec un des sommets ?
 375 c
 376             if ( logaux(7) ) then
 377 c
 378               if ( nbpbco(-1).gt.0 ) then
 379 c
 380                 nucoin = numcoi(lenoeu)
 381 c
 382                 if ( nucoin.ne.0 ) then
 383 c
 384                   ptcode = coinpt(nucoin)+1
 385                   ptcofi = coinpt(nucoin+1)
 386                   do 232 , ptcoin = ptcode, ptcofi
 387                     jaux = coinnn(ptcoin)
 388                     do 2321 , iaux = 1 , nbsomm
 389                       if ( jaux.eq.sommet(iaux) ) then
 390                         goto 23
 391                       endif
 392  2321               continue
 393   232             continue
 394 c
 395                 endif
 396 c
 397               endif
 398 c
 399             endif
 400 c
 401 c 2.3.3. ==> exclusivement les noeuds p1
 402 c
 403             if ( logaux(7) ) then
 404 c
 405               if ( hetnoe(lenoeu).ne.1  .and.
 406      >             hetnoe(lenoeu).ne.21 .and.
 407      >             hetnoe(lenoeu).ne.51 ) then
 408 c
 409                 logaux(7) = .false.
 410 c
 411               endif
 412 c
 413             endif
 414 c
 415 c 2.3.3. ==> En 3D, le noeud est-il dans le plan du quadrangle ?
 416 c            on calcule les deux vecteurs normaux aux quadrangles
 417 c            (s1,s2,s3,s4) et (s1,s2,s3,n). En fait, il suffit de
 418 c            s'interesser aux vecteurs normaux aux sous-triangles
 419 c            (s1,s2,s4) et (s1,s2,n), vecteurs v6 et v5
 420 c            on cherche s'ils sont paralleles, c'est-a-dire de produit
 421 c            vectoriel nul
 422 c
 423             if ( logaux(7) ) then
 424 c
 425 c            v5(1,.) est le produit vectoriel n1n2 x n1n.
 426 c
 427               if ( sdim.eq.3 ) then
 428 c
 429                 if ( logaux(3) ) then
 430 c
 431                   v12(1) = v2(1)-v1(1)
 432                   v12(2) = v2(2)-v1(2)
 433                   v12(3) = v2(3)-v1(3)
 434 c
 435                   logaux(3) = .false.
 436 c
 437                 endif
 438 c
 439                 v5(1,1) = v12(2) * (vn(3)-v1(3))
 440      >                  - v12(3) * (vn(2)-v1(2))
 441                 v5(1,2) = v12(3) * (vn(1)-v1(1))
 442      >                  - v12(1) * (vn(3)-v1(3))
 443                 v5(1,3) = v12(1) * (vn(2)-v1(2))
 444      >                  - v12(2) * (vn(1)-v1(1))
 445 c
 446 c            daux2 represente la norme du produit vectoriel v5 x v6(1,.)
 447 c
 448                 daux2 = abs ( v6(1,2)*v5(1,3) - v6(1,3)*v5(1,2) )
 449      >                + abs ( v6(1,3)*v5(1,1) - v6(1,1)*v5(1,3) )
 450      >                + abs ( v6(1,1)*v5(1,2) - v6(1,2)*v5(1,1) )
 451 c
 452               else
 453 c
 454                 daux2 = -1.d0
 455 c
 456               endif
 457 c
 458 c 2.3.4. ==> dans ce plan, n est-il dedans ?
 459 c            cela est vrai si les 4 triangles batis sur ce noeud sont
 460 c            orientes de la meme facon. Pour cela, on regarde si les
 461 c            quatre produits vectoriels (na,nb) sont de meme orientation
 462 c            pour les quatre permutations circulaires sur (a,b,c,d),
 463 c            c'est-a-dire si le produit scalaire des produits vectoriels
 464 c            est de meme signe.
 465 c
 466 c           n1 x------------------x n2
 467 c              x  .               x .
 468 c              x     .            x  .
 469 c              x        .         x   .
 470 c              x           .      x    .
 471 c              x              .   x     .
 472 c              x                 .x      .
 473 c              x                  x .     .
 474 c           n4 x------------------x    .   .
 475 c                             .   n3  .  . .
 476 c                                            . n
 477 c
 478 c
 479               if ( daux2.le.epsima ) then
 480 c
 481 c 2.3.4.1. ==> critere absolu ou relatif
 482 c
 483                 daux1 = 0.d0
 484 c
 485 c 2.3.4.2. ==> triangle (n,n1,n2)
 486 c
 487 c            v5(1,.) est le produit vectoriel nn1 x nn2
 488 c
 489                 v5(1,1) = (v1(2)-vn(2)) * (v2(3)-vn(3))
 490      >                  - (v1(3)-vn(3)) * (v2(2)-vn(2))
 491                 v5(1,2) = (v1(3)-vn(3)) * (v2(1)-vn(1))
 492      >                  - (v1(1)-vn(1)) * (v2(3)-vn(3))
 493                 v5(1,3) = (v1(1)-vn(1)) * (v2(2)-vn(2))
 494      >                  - (v1(2)-vn(2)) * (v2(1)-vn(1))
 495 c
 496 c 2.3.4.3. ==> triangle (n,n2,n3)
 497 c
 498 c            v5(2,.) est le produit vectoriel nn2 x nn3
 499 c
 500                 v5(2,1) = (v2(2)-vn(2)) * (v3(3)-vn(3))
 501      >                  - (v2(3)-vn(3)) * (v3(2)-vn(2))
 502                 v5(2,2) = (v2(3)-vn(3)) * (v3(1)-vn(1))
 503      >                  - (v2(1)-vn(1)) * (v3(3)-vn(3))
 504                 v5(2,3) = (v2(1)-vn(1)) * (v3(2)-vn(2))
 505      >                  - (v2(2)-vn(2)) * (v3(1)-vn(1))
 506 c
 507 c            si v5(1,.) est le vecteur nul, cela signifie que les noeuds
 508 c            n, n1 et n2 sont alignes. Il ne faut donc pas se poser la
 509 c            question du triangle (n,n1,n2). On remplace alors (n,n1,n2)
 510 c            par (n,n2,n3) comme triangle de reference, c'est-a-dire
 511 c            v5(1,.) par v5(2,.).
 512 c            On aura remarque que dans ces conditions le vecteur v5(2,.)
 513 c            ne peut pas etre nul, sinon cela signifierait que les 3
 514 c            noeuds n1, n2 et n3 sont alignes.
 515 c
 516                 if ( ( v5(1,1)*v5(1,1) + v5(1,2)*v5(1,2)
 517      >               + v5(1,3)*v5(1,3) ) .le. daux1 ) then
 518                   v5(1,1) = v5(2,1)
 519                   v5(1,2) = v5(2,2)
 520                   v5(1,3) = v5(2,3)
 521                 else
 522                   prosca =
 523      >             v5(1,1)*v5(2,1) + v5(1,2)*v5(2,2) + v5(1,3)*v5(2,3)
 524 cgn            if ( lequad.eq.1270 .and. lenoeu.eq.7119 ) then
 525 cgn            print *,v5(1,1),v5(1,2),v5(1,3)
 526 cgn            print *,v5(2,1),v5(2,2),v5(2,3)
 527 cgn            print *,'==> prosca = ',prosca
 528 cgn            endif
 529 c
 530                   if ( prosca.lt.daux1 ) then
 531                     goto 23
 532                   endif
 533                 endif
 534 c
 535 c 2.3.4.4. ==> triangle (n,n3,n4)
 536 c
 537 c            v5(2,.) est le produit vectoriel nn3 x nn4
 538 c
 539                 v5(2,1) = (v3(2)-vn(2)) * (v4(3)-vn(3))
 540      >                  - (v3(3)-vn(3)) * (v4(2)-vn(2))
 541                 v5(2,2) = (v3(3)-vn(3)) * (v4(1)-vn(1))
 542      >                  - (v3(1)-vn(1)) * (v4(3)-vn(3))
 543                 v5(2,3) = (v3(1)-vn(1)) * (v4(2)-vn(2))
 544      >                  - (v3(2)-vn(2)) * (v4(1)-vn(1))
 545 c
 546                 prosca =
 547      >             v5(1,1)*v5(2,1) + v5(1,2)*v5(2,2) + v5(1,3)*v5(2,3)
 548 c
 549 cgn            if ( lequad.eq.1270 .and. lenoeu.eq.7119 ) then
 550 cgn            print *,v5(2,1),v5(2,2),v5(2,3)
 551 cgn            print *,'==> prosca = ',prosca
 552 cgn            endif
 553                 if ( prosca.lt.daux1 ) then
 554                   goto 23
 555                 endif
 556 c
 557 c 2.3.4.5. ==> triangle (n,n4,n1)
 558 c
 559 c            v5(2,.) est le produit vectoriel nn4 x nn1
 560 c
 561                 v5(2,1) = (v4(2)-vn(2)) * (v1(3)-vn(3))
 562      >                  - (v4(3)-vn(3)) * (v1(2)-vn(2))
 563                 v5(2,2) = (v4(3)-vn(3)) * (v1(1)-vn(1))
 564      >                  - (v4(1)-vn(1)) * (v1(3)-vn(3))
 565                 v5(2,3) = (v4(1)-vn(1)) * (v1(2)-vn(2))
 566      >                  - (v4(2)-vn(2)) * (v1(1)-vn(1))
 567 c
 568                 prosca =
 569      >             v5(1,1)*v5(2,1) + v5(1,2)*v5(2,2) + v5(1,3)*v5(2,3)
 570 cgn            if ( lequad.eq.1270 .and. lenoeu.eq.7119 ) then
 571 cgn            print *,v5(2,1),v5(2,2),v5(2,3)
 572 cgn            print *,'==> prosca = ',prosca
 573 cgn            endif
 574 c
 575                 if ( prosca.lt.daux1 ) then
 576                   goto 23
 577                 endif
 578 c
 579 c 2.3.5. ==> si logaux(7) est encore vrai, c'est que le noeud est
 580 c            a l'interieur du quadrangle ... malaise ...
 581 c
 582                 iaux = lequad
 583 c
 584 #include "utb302.h"
 585 c
 586                 if ( sdim.eq.2 ) then
 587                   write (ulbila,14202) sommet(1), v1(1), v1(2)
 588                   write (ulbila,14202) sommet(2), v2(1), v2(2)
 589                   write (ulbila,14202) sommet(3), v3(1), v3(2)
 590                   write (ulbila,14202) sommet(4), v4(1), v4(2)
 591                 else
 592                   write (ulbila,14203) sommet(1), v1(1), v1(2), v1(3)
 593                   write (ulbila,14203) sommet(2), v2(1), v2(2), v2(3)
 594                   write (ulbila,14203) sommet(3), v3(1), v3(2), v3(3)
 595                   write (ulbila,14203) sommet(4), v4(1), v4(2), v4(3)
 596                 endif
 597 c
 598                 write (ulbila,10200)
 599 c
 600               endif
 601 c
 602             endif
 603 c
 604    23     continue
 605 c
 606         endif
 607 c
 608    20 continue
 609 cgn      call gtfims (112)
 610 c
 611       end