src/tool/Suivi_Frontiere/sfseno.F

   1       subroutine sfseno ( cono, numlig, unst2x, epsid2,
   2      >                    seglig, somseg, geocoo, abscur,
   3      >                    seg, acnoeu )
   4 c ______________________________________________________________________
   5 c
   6 c                             H O M A R D
   7 c
   8 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
   9 c
  10 c Version originale enregistree le 18 juin 1996 sous le numero 96036
  11 c aupres des huissiers de justice Simart et Lavoir a Clamart
  12 c Version 11.2 enregistree le 13 fevrier 2015 sous le numero 2015/014
  13 c aupres des huissiers de justice
  14 c Lavoir, Silinski & Cherqui-Abrahmi a Clamart
  15 c
  16 c    HOMARD est une marque deposee d'Electricite de France
  17 c
  18 c Copyright EDF 1996
  19 c Copyright EDF 1998
  20 c Copyright EDF 2002
  21 c Copyright EDF 2020
  22 c ______________________________________________________________________
  23 c
  24 c   Suivi de Frontiere - SEgment - NOeud
  25 c   -        -           --        --
  26 c ______________________________________________________________________
  27 c
  28 c but : recherche du segment de la ligne auquel appartient le noeud
  29 c ______________________________________________________________________
  30 c .        .     .        .                                            .
  31 c .  nom   . e/s . taille .           description                      .
  32 c .____________________________________________________________________.
  33 c . cono   . es  .  sdim  . coordonnees du noeud                       .
  34 c . numlig . e   .   1    . numero de ligne sur laquelle on cherche    .
  35 c . unst2x . e   .   1    . inverse de la taille maximale au carre     .
  36 c . epsid2 . e   .   1    . precision relative pour carre de distance  .
  37 c . seglig . e   .0:sfnbli. pointeur dans le tableau somseg : les      .
  38 c .        .     .        . segments de la ligne i sont aux places de  .
  39 c .        .     .        . seglig(i-1)+1 a seglig(i)-1 inclus         .
  40 c . somseg . e   . sfnbse . liste des sommets des lignes separees par  .
  41 c                           des 0                                      .
  42 c . geocoo . e   .sfnbso**. coordonnees des sommets de la frontiere    .
  43 c . seg    .  s  .   1    . numero de segment trouve                   .
  44 c                           - : noeud sur extremite de ligne           .
  45 c                           + : noeud hors extremites                  .
  46 c . acnoeu .  s  .   1    . abscisse curviligne du noeud sur la ligne  .
  47 c ______________________________________________________________________
  48 c
  49 c====
  50 c 0. declarations et dimensionnement
  51 c====
  52 c
  53 c 0.1. ==> generalites
  54 c
  55       implicit none
  56       save
  57 c
  58 c 0.2. ==> communs
  59 c
  60 #include "front1.h"
  61 #include "infini.h"
  62 #include "envca1.h"
  63 c
  64 c 0.3. ==> arguments
  65 c
  66       double precision unst2x, epsid2
  67       double precision geocoo(sfnbso,*), cono(sdim)
  68       integer seglig(0:sfnbli), somseg(sfnbse)
  69       integer numlig, seg
  70       double precision acnoeu
  71       double precision abscur(sfnbse)
  72 c
  73 c 0.4. ==> variables locales
  74 c
  75       integer iaux, jaux, kaux
  76       integer seg1, seg2, segmin
  77 c
  78       double precision cooa(3), coob(3)
  79       double precision daux, daux1, daux2
  80 c
  81 #ifdef _DEBUG_HOMARD_
  82       integer ulsort
  83       parameter (ulsort=1)
  84       integer glop
  85       common / tutu / glop
  86 #endif
  87 c
  88 c 0.5. ==> initialisations
  89 c ______________________________________________________________________
  90 c
  91 #include "impr03.h"
  92 c
  93 c====
  94 c 1. initialisation
  95 c====
  96 c
  97 #ifdef _DEBUG_HOMARD_
  98       if ( glop.ne.0 ) then
  99       write(ulsort,*) 'Entree dans SFSENO'
 100       write(ulsort,90004) 'unst2x',unst2x
 101       write(ulsort,90004) 'epsid2',epsid2
 102       write(ulsort,90004) 'cono', (cono(iaux),iaux = 1 , sdim)
 103       write(ulsort,90002) 'numlig', numlig
 104       endif
 105 #endif
 106 c
 107 c 1.1. ==> segment pas encore trouve
 108 c
 109       seg = 0
 110 c
 111 c 1.2. ==> bornes de la ligne dans la numerotation des segments
 112 c
 113       seg1 = seglig(numlig-1) + 1
 114       seg2 = seglig(numlig  ) - 2
 115 #ifdef _DEBUG_HOMARD_
 116       if ( glop.ne.0 ) then
 117       write(ulsort,90006) 'Pointeurs seglig pour la ligne', numlig,
 118      >' de',seg1,' a ',seg2+2
 119       write(ulsort,90002) '==> Segments extremites',seg1,seg2
 120       write(ulsort,90002) '==> Sommets extremites ',
 121      >                    somseg(seg1),somseg(seg2+1)
 122         endif
 123 #endif
 124 c
 125 c====
 126 c 2. Le noeud est-il une extremite de la ligne ?
 127 c    Remarque : on commence par la fin de la ligne pour pouvoir
 128 c               initialiser correctement l'etape suivante en cas de
 129 c               non coincidence
 130 c    Remarque : comme ce sont les coordonnees du debut du segment qui
 131 c               sont stockees, il faut examiner le segment fictif
 132 c               (dernier+1) pour trouver les coordonnees de sa fin
 133 c====
 134 c
 135       do 20 , jaux = 1 , 2
 136 c
 137         if ( jaux.eq.1 ) then
 138           kaux = seg2+1
 139         else
 140           kaux = seg1
 141         endif
 142 c
 143         daux = 0.d0
 144         do 21 , iaux = 1 , sdim
 145           coob(iaux) = geocoo(somseg(kaux),iaux)
 146           daux = daux + (coob(iaux)-cono(iaux))**2
 147    21   continue
 148 c
 149 #ifdef _DEBUG_HOMARD_
 150       if ( glop.ne.0 ) then
 151         if ( jaux.eq.1 ) then
 152           write(ulsort,90002) 'segment (jaux=1)', seg2
 153         else
 154           write(ulsort,90002) 'segment (jaux=2)', seg1
 155         endif
 156         write(ulsort,90024) 'Noeud', somseg(kaux),
 157      >  (geocoo(somseg(kaux),iaux),iaux=1,sdim)
 158         write(ulsort,90024) '==> carre distance a l''extremite',
 159      >   1+mod(jaux,2),daux
 160         endif
 161 #endif
 162         if ( daux*unst2x.le.epsid2 ) then
 163           seg = -kaux
 164           acnoeu = abscur(kaux)
 165 #ifdef _DEBUG_HOMARD_
 166         if ( glop.ne.0 ) then
 167           write(ulsort,90015) '==> Le noeud est l''extremite',
 168      >      1+mod(jaux,2),' de la ligne', numlig
 169         endif
 170 #endif
 171           goto 50
 172 c
 173         endif
 174 c
 175    20 continue
 176 c
 177 c====
 178 c 3. Le noeud n'est pas une extremite de la ligne
 179 c    On va chercher le segment dont la premiere extremite est
 180 c    la plus proche du noeud courant.
 181 c    Si une de ces extremites est le noeud courant, on le note
 182 c    et on sort.
 183 c====
 184 c
 185 #ifdef _DEBUG_HOMARD_
 186       if ( glop.ne.0 ) then
 187       write(ulsort,90002) 'Le noeud n''est pas extremite de la ligne',
 188      >                     numlig
 189       endif
 190 #endif
 191 c
 192       daux1 = vinfpo
 193 c
 194       do 30 , kaux = seg1 , seg2
 195 c
 196 #ifdef _DEBUG_HOMARD_
 197       if ( glop.ne.0 ) then
 198         if ( kaux.eq.seg1 .or. kaux.le.-1 ) then
 199         write(ulsort,90006) '    Segment [',
 200      >                       somseg(kaux),' ',somseg(kaux+1),'] '
 201         endif
 202         endif
 203 #endif
 204 c
 205 c 3.1. ==> Calcul de la distance au point A, debut du segment
 206 c
 207         daux = 0.d0
 208         do 31 , iaux = 1 , sdim
 209           daux = daux + (geocoo(somseg(kaux),iaux)-cono(iaux))**2
 210    31   continue
 211 cgn        if ( kaux.eq.seg1 .or. kaux.le.-1 ) then
 212 cgn          write(ulsort,90004) '      Carre de la distance', daux
 213 cgn        endif
 214 c
 215         if ( daux*unst2x.le.epsid2 ) then
 216 #ifdef _DEBUG_HOMARD_
 217       if ( glop.ne.0 ) then
 218           write(ulsort,90004) 'Carre de la distance', daux
 219           write(ulsort,90002)
 220      >    '==> Le noeud est au debut du segment', kaux
 221         endif
 222 #endif
 223           seg = kaux
 224           acnoeu = abscur(kaux)
 225           goto 50
 226 c
 227         endif
 228 c
 229 c 3.2. ==> Memorisation du minimum
 230 c
 231 cgn        write(ulsort,90004) '     Distance ', daux
 232         if ( daux.le.daux1 ) then
 233 c
 234           daux1 = daux
 235           segmin = kaux
 236 cgn        write(ulsort,*) '     Minimum pour le segment ',segmin,
 237 cgn     >               ' : [',somseg(segmin),',',somseg(segmin+1),'] '
 238 c
 239         endif
 240 c
 241    30 continue
 242 c
 243 #ifdef _DEBUG_HOMARD_
 244       if ( glop.ne.0 ) then
 245         write(ulsort,90004) 'Distance minimale', daux1
 246         write(ulsort,90006) '    atteinte sur le segment ',segmin,
 247      >                  '[',somseg(segmin),' ',somseg(segmin+1),'] '
 248       endif
 249 #endif
 250 c
 251 c====
 252 c 4. Le noeud n'est pas une extremite d'un segment
 253 c    On precise les points A et B
 254 c====
 255 c
 256 #ifdef _DEBUG_HOMARD_
 257       if ( glop.ne.0 ) then
 258       write(ulsort,90002) 'Le premier segment est seg1', seg1
 259       write(ulsort,90002) 'Le dernier segment est seg2', seg2
 260       write(ulsort,90002)
 261      >  'Le segment le plus proche du noeud est segmin', segmin
 262       endif
 263 #endif
 264 c
 265 c 4.1. ==> On a trouve le segment dont le debut, A, est le sommet le
 266 c          plus proche de N.
 267 c          Si c'est le premier segment, il n'y a pas d'equivoque :
 268 c          N sera place entre A et B
 269 c          Sinon, il faut preciser entre les deux segments concernes,
 270 c          segmin-1 et segmin. On postule que le segment le plus proche
 271 c          est celui dont la deuxieme extremite est la plus proche de N.
 272 c
 273 c              N
 274 c                       A
 275 c                       .
 276 c        segmin-1  .            . segmin
 277 c             .                         .
 278 c      C .                                      . B
 279 c
 280       if ( segmin.ne.seg1 ) then
 281 c
 282         daux  = 0.d0
 283         daux2 = 0.d0
 284         do 41 , iaux = 1 , sdim
 285           daux  = daux  +
 286      >           (geocoo(somseg(segmin-1),iaux)-cono(iaux))**2
 287           daux2 = daux2 +
 288      >           (geocoo(somseg(segmin+1),iaux)-cono(iaux))**2
 289    41   continue
 290 c
 291 #ifdef _DEBUG_HOMARD_
 292       if ( glop.ne.0 ) then
 293        write(ulsort,90024) '(C-A) Carre de la distance au debut de',
 294      >          segmin-1, daux
 295        write(ulsort,90024) '(A-B) Carre de la distance au debut de',
 296      >          segmin  , daux1
 297        write(ulsort,90024) '(B- ) Carre de la distance au debut de',
 298      >          segmin+1, daux2
 299        write(ulsort,90004) '    N',(cono(iaux),iaux=1,sdim)
 300        write(ulsort,90004) '    C',
 301      >  (geocoo(somseg(segmin-1),iaux),iaux=1,sdim), abscur(segmin-1)
 302        write(ulsort,90004) '    A',
 303      >  (geocoo(somseg(segmin),iaux),iaux=1,sdim), abscur(segmin)
 304        write(ulsort,90004) '    B',
 305      >  (geocoo(somseg(segmin+1),iaux),iaux=1,sdim), abscur(segmin+1)
 306       endif
 307 #endif
 308 c
 309         if ( daux.lt.daux2 ) then
 310           seg = segmin - 1
 311         else
 312           seg = segmin
 313         endif
 314 c
 315       else
 316 c
 317         seg = segmin
 318 c
 319       endif
 320 c
 321 c 4.2. ==> Definition des points A et B
 322 c
 323 #ifdef _DEBUG_HOMARD_
 324       if ( glop.ne.0 ) then
 325       write(ulsort,90002) '... ==> Le noeud est lie au segment', seg
 326       endif
 327 #endif
 328 c
 329       do 42 , iaux = 1 , sdim
 330         cooa(iaux) = geocoo(somseg(seg  ),iaux)
 331         coob(iaux) = geocoo(somseg(seg+1),iaux)
 332    42 continue
 333 c
 334 #ifdef _DEBUG_HOMARD_
 335       if ( glop.ne.0 ) then
 336       write(ulsort,90004) '   de A', (cooa(iaux),iaux=1,sdim),
 337      >                   abscur(seg)
 338       write(ulsort,90004) '    a B', (coob(iaux),iaux=1,sdim),
 339      >                   abscur(seg+1)
 340       endif
 341 #endif
 342 c
 343 c 4.3. ==> Il faut recoller N sur la frontiere le cas echeant,
 344 c          puis redefinir son abscisse curviligne.
 345 c
 346 c 4.3.1. ==> Calcul du produit scalaire AB.AN
 347 c
 348       daux = 0.d0
 349       do 431 , iaux = 1 , sdim
 350         daux = daux +
 351      >       ( coob(iaux) - cooa(iaux) ) *
 352      >       ( cono(iaux) - cooa(iaux) )
 353   431   continue
 354 #ifdef _DEBUG_HOMARD_
 355       if ( glop.ne.0 ) then
 356       write(ulsort,90004) 'Produit scalaire AB.AN', daux
 357       endif
 358 #endif
 359 c
 360 c 4.3.2. ==> Positionnement
 361 c 4.3.2.1. ==> Si le produit scalaire est negatif, c'est que P est
 362 c              "en arriere" de A. On ramene P sur A.
 363 c
 364       if ( daux.le.0 ) then
 365 c
 366         do 4321 , iaux = 1 , sdim
 367           cono(iaux) = cooa(iaux)
 368  4321   continue
 369         acnoeu = abscur(seg)
 370 #ifdef _DEBUG_HOMARD_
 371       if ( glop.ne.0 ) then
 372         write(ulsort,90004) '.. ==> Le noeud a ete replace sur A'
 373       endif
 374 #endif
 375 c
 376       else
 377 c
 378 c         daux1 : distance AB
 379         daux1 = abscur(seg+1)-abscur(seg)
 380         daux = daux / (daux1**2)
 381 #ifdef _DEBUG_HOMARD_
 382       if ( glop.ne.0 ) then
 383         write(ulsort,90004) 'Distance AB', daux1
 384         write(ulsort,90004) 'Carre de la distance AB  ', daux1**2
 385       endif
 386 #endif
 387 c
 388 c 4.3.2.2. ==> Si le produit scalaire est superieur au carre de la
 389 c              distance AB, c'est que P est "en avant" de B.
 390 c              On ramenera P sur B.
 391 c
 392         if ( daux.ge.1.d0 ) then
 393 c
 394           do 4322 , iaux = 1 , sdim
 395             cono(iaux) = coob(iaux)
 396  4322     continue
 397           acnoeu = abscur(seg+1)
 398 #ifdef _DEBUG_HOMARD_
 399       if ( glop.ne.0 ) then
 400         write(ulsort,90004) '.. ==> Le noeud a ete replace sur B'
 401       endif
 402 #endif
 403 c
 404 c 4.3.2.3. ==> N est "entre" A et B.
 405 c              On decompose le vecteur AN en une partie le long du
 406 c              segment, alpha.AB, et une partie orthogonale.
 407 c              Cela revient a projeter N orthogonalement au segment.
 408 c
 409         else
 410 c
 411           do 4323 , iaux = 1 , sdim
 412             cono(iaux) = cooa(iaux) + daux*(coob(iaux)-cooa(iaux))
 413  4323     continue
 414           acnoeu = abscur(seg) + daux*daux1
 415 #ifdef _DEBUG_HOMARD_
 416       if ( glop.ne.0 ) then
 417         write(ulsort,90004) '.. ==> Le noeud a ete place entre A et B'
 418       endif
 419 #endif
 420         endif
 421 c
 422       endif
 423 c
 424 #ifdef _DEBUG_HOMARD_
 425       if ( glop.ne.0 ) then
 426         write(ulsort,90004) '. Nouveau noeud', (cono(iaux),iaux=1,sdim)
 427       endif
 428 #endif
 429 c
 430 c====
 431 c 5. Sortie
 432 c====
 433 c
 434 #ifdef _DEBUG_HOMARD_
 435       write(ulsort,*) '5. sortie'
 436 #endif
 437 c
 438    50 continue
 439 c
 440 #ifdef _DEBUG_HOMARD_
 441       if ( glop.ne.0 ) then
 442       write(ulsort,*) 'Bilan :'
 443       write(ulsort,90002) '==> Segment ', seg
 444       write(ulsort,90004) '==> Abcisse curviligne', acnoeu
 445       write(ulsort,*) 'Sortie de SFSENO'
 446       endif
 447 #endif
 448 c
 449       end