src/tool/AV_Conversion/vcms22.F

   1       subroutine vcms22 ( maextr,
   2      >                    nbnoto, nbelem,
   3      >                    nbse2d, nbtr2d, nbqu2d, nbele2,
   4      >                    nu2dno, coonoe,
   5      >                    fameel, typele, noeele,
   6      >                    fame2d, type2d, noee2d,
   7      >                    faminf, famsup,
   8      >                    ulsort, langue, codret )
   9 c ______________________________________________________________________
  10 c
  11 c                             H O M A R D
  12 c
  13 c Outil de Maillage Adaptatif par Raffinement et Deraffinement d'EDF R&D
  14 c
  15 c Version originale enregistree le 18 juin 1996 sous le numero 96036
  16 c aupres des huissiers de justice Simart et Lavoir a Clamart
  17 c Version 11.2 enregistree le 13 fevrier 2015 sous le numero 2015/014
  18 c aupres des huissiers de justice
  19 c Lavoir, Silinski & Cherqui-Abrahmi a Clamart
  20 c
  21 c    HOMARD est une marque deposee d'Electricite de France
  22 c
  23 c Copyright EDF 1996
  24 c Copyright EDF 1998
  25 c Copyright EDF 2002
  26 c Copyright EDF 2020
  27 c ______________________________________________________________________
  28 c
  29 c    aVant adaptation - Conversion de Maillage -
  30 c     -                 -             -
  31 c                       Saturne 2D - phase 2 - Neptune 2D
  32 c                       -       -          -
  33 c ______________________________________________________________________
  34 c .        .     .        .                                            .
  35 c .  nom   . e/s . taille .           description                      .
  36 c .____________________________________________________________________.
  37 c . maextr . e   .   1    . maillage extrude                           .
  38 c .        .     .        . 0 : non (defaut)                           .
  39 c .        .     .        . 1 : selon X                                .
  40 c .        .     .        . 2 : selon Y                                .
  41 c .        .     .        . 3 : selon Z (cas de Saturne ou Neptune)    .
  42 c . nbnoto . e   .   1    . nombre de noeuds du maillage externe       .
  43 c . nbse2d . e   .   1    . nombre de segments du maillage 2d          .
  44 c . nbtr2d . e   .   1    . nombre de triangles du maillage 2d         .
  45 c . nbqu2d . e   .   1    . nombre de quadrangles du maillage 2d       .
  46 c . nbelem . e   .   1    . nombre d'elements du maillage externe      .
  47 c . nu2dno . e   . nbnoto . numero du calcul des noeuds saturne/neptune.
  48 c . coonoe . e   . nbnoto . coordonnees des noeuds                     .
  49 c .        .     . * sdim .                                            .
  50 c . fameel . e   . nbelem . famille med des elements                   .
  51 c . typele . e   . nbelem . type des elements pour le code de calcul   .
  52 c . noeele . e   . nbelem . noeuds des elements                        .
  53 c .        .     .*nbmane .                                            .
  54 c . fame2d .  s  . nbele2 . famille med des elements du maillage 2d    .
  55 c . type2d .  s  . nbele2 . type des elements du maillage 2d           .
  56 c . noee2d .  s  . nbele2 . noeuds des elements du maillage 2d         .
  57 c .        .     .*nbman2 .                                            .
  58 c . faminf .  s  .   1    . famille med des quad de la face inferieure .
  59 c . famsup .  s  .   1    . famille med des quad de la face superieure .
  60 c . ulsort . e   .   1    . numero d'unite logique de la liste standard.
  61 c . langue . e   .    1   . langue des messages                        .
  62 c .        .     .        . 1 : francais, 2 : anglais                  .
  63 c . codret . es  .    1   . code de retour des modules                 .
  64 c .        .     .        . 0 : pas de probleme                        .
  65 c ______________________________________________________________________
  66 c
  67 c====
  68 c 0. declarations et dimensionnement
  69 c====
  70 c
  71 c 0.1. ==> generalites
  72 c
  73       implicit none
  74       save
  75 c
  76       character*6 nompro
  77       parameter ( nompro = 'VCMS22' )
  78 c
  79 #include "nblang.h"
  80 #include "consta.h"
  81 #include "consts.h"
  82 #include "fractc.h"
  83 c
  84 c 0.2. ==> communs
  85 c
  86 #include "envex1.h"
  87 c
  88 #include "infini.h"
  89 #include "meddc0.h"
  90 #include "impr02.h"
  91 #include "op0123.h"
  92 c
  93 c 0.3. ==> arguments
  94 c
  95       integer maextr
  96       integer nbnoto
  97       integer nbse2d, nbtr2d, nbqu2d, nbele2
  98       integer nbelem
  99       integer faminf, famsup
 100       integer nu2dno(nbnoto)
 101       integer fameel(nbelem), typele(nbelem), noeele(nbelem,*)
 102       integer fame2d(nbele2), type2d(nbele2), noee2d(nbele2,*)
 103 c
 104       double precision coonoe(nbnoto,*)
 105 c
 106       integer ulsort, langue, codret
 107 c
 108 c 0.4. ==> variables locales
 109 c
 110       integer iaux, jaux, kaux, laux
 111       integer iaux1, iaux2
 112       integer tbiaux(4)
 113       integer el, nuel2d
 114       integer numloc(4)
 115 c
 116       character*1 saux01
 117 c
 118       double precision xymil(2)
 119       double precision v1(2), vn(2)
 120       double precision daux1, daux2
 121       double precision epsilo
 122       double precision daux(0:4)
 123 c
 124       integer nbmess
 125       parameter ( nbmess = 20 )
 126       character*80 texte(nblang,nbmess)
 127 c
 128 c 0.5. ==> initialisations
 129 c ______________________________________________________________________
 130 c
 131 c====
 132 c 1. messages
 133 c====
 134 c
 135 #include "impr01.h"
 136 c
 137 #ifdef _DEBUG_HOMARD_
 138       write (ulsort,texte(langue,1)) 'Entree', nompro
 139       call dmflsh (iaux)
 140 #endif
 141 c
 142       texte(1,4) = '(''Maille numero :'',i10,'', de noeuds '',8i10)'
 143       texte(1,5) = '(i1,'' noeud(s) sont dans le plan cooinf.'')'
 144       texte(1,6) = '(''Pour un '',a,'', il en faudrait '',a)'
 145       texte(1,7) = '(''Famille de la face '',a,'' : '',i6)'
 146       texte(1,8) = '(''Famille du '',a,i10,'' : '',i6)'
 147       texte(1,9) =
 148      >'(''Nombre de '',a,'' attendus pour le maillage 2D :'',i10)'
 149       texte(1,10) =
 150      >'(''Nombre de '',a,'' trouves pour le maillage 2D  :'',i10)'
 151       texte(1,11) = '(''Element '',i10,'' ('',a,''), numloc = '',4i10)'
 152 c
 153       texte(2,4) = '(''Mesh # :'',i10,'', with nodes '',8i10)'
 154       texte(2,5) = '(i1,'' node(s) are in cooinf plane.'')'
 155       texte(2,6) = '(''For '',a,'', '',a,'' were expected.'')'
 156       texte(2,7) = '(''Family for '',a,'' face : '',i6)'
 157       texte(2,8) = '(''Family for '',a,'' #'',i10,'' : '',i6)'
 158       texte(2,9) =
 159      > '(''Expected number of '',a,'' for the 2D mesh :'',i10)'
 160       texte(2,10) =
 161      > '(''Found number of '',a,'' for the 2D mesh    :'',i10)'
 162       texte(2,11) = '(''Element '',i10,'' ('',a,''), numloc = '',4i10)'
 163 c
 164 #include "impr03.h"
 165 c
 166       codret = 0
 167 c
 168       if ( maextr.eq.1 ) then
 169         saux01 = 'X'
 170         iaux1 = 2
 171         iaux2 = 3
 172       elseif ( maextr.eq.2 ) then
 173         saux01 = 'Y'
 174         iaux1 = 1
 175         iaux2 = 3
 176       elseif ( maextr.eq.3 ) then
 177         saux01 = 'Z'
 178         iaux1 = 1
 179         iaux2 = 2
 180       else
 181         codret = 1
 182       endif
 183 c
 184       nuel2d = 0
 185 #ifdef _DEBUG_HOMARD_
 186       write (ulsort,90002) 'maextr', maextr
 187       write (ulsort,90002) 'nbelem', nbelem
 188       write (ulsort,90002) 'nbele2', nbele2
 189       write (ulsort,90002) 'nbse2d', nbse2d
 190       write (ulsort,90002) 'nbqu2d', nbqu2d
 191       write (ulsort,90002) 'nbtr2d', nbtr2d
 192 #endif
 193 c
 194       epsilo = 1.d-10 * pi
 195 c
 196 c====
 197 c 2. transformations des quadrangles en segments
 198 c    Attention a les mettre avant les quadrangles pour respecter
 199 c    les conventions ...
 200 c====
 201 #ifdef _DEBUG_HOMARD_
 202       write (ulsort,90002) '2. quad -> segm ; codret', codret
 203 #endif
 204 c
 205       faminf = 1
 206       famsup = 1
 207 c
 208       do 21 , el = 1 , nbelem
 209 c
 210         if ( codret.eq.0 ) then
 211 c
 212         if ( typele(el).eq.edqua4 ) then
 213 c
 214 c 2.1. ==> recherche des noeuds dans le plan cooinf
 215 c
 216           jaux = 0
 217           do 211 , iaux = 1 , 4
 218             if ( nu2dno(noeele(el,iaux)).ne.0 ) then
 219               jaux = jaux + 1
 220               numloc(jaux) = noeele(el,iaux)
 221               if ( jaux.eq.1 ) then
 222                 kaux = iaux
 223               else
 224                 laux = iaux
 225               endif
 226             endif
 227   211     continue
 228 c
 229 #ifdef _DEBUG_HOMARD_
 230       if ( jaux.gt.0 .and. el.lt.1 ) then
 231       write (ulsort,texte(langue,11)) el, 'quad',
 232      >                               (numloc(iaux),iaux=1,jaux)
 233       endif
 234 #endif
 235 c
 236 c 2.2. ==> si exactement 2 noeuds sont dans le plan, on cree le segment
 237 c          Le segment est cree avec les noeuds dans l'ordre
 238 c          d'apparition dans la connectivite du quadrangle. Cela
 239 c          permettra a la reconstitution du quadrangle apres adaptation
 240 c          de retrouver la meme orientation de la face.
 241 c          attention au retournemnt eventuel ...
 242 c
 243           if ( jaux.eq.2 ) then
 244 c
 245             nuel2d = nuel2d + 1
 246             if ( kaux.eq.1 .and. laux.eq.4 ) then
 247               iaux = numloc(2)
 248               numloc(2) = numloc(1)
 249               numloc(1) = iaux
 250             endif
 251             noee2d(nuel2d,1) = nu2dno(numloc(1))
 252             noee2d(nuel2d,2) = nu2dno(numloc(2))
 253             fame2d(nuel2d) = fameel(el)
 254             type2d(nuel2d) = edseg2
 255 c
 256 c 2.3. ==> si exactement 4 noeuds sont dans ce plan, c'est la face
 257 c          inferieure. On memorise le numero de famille
 258 c
 259           elseif ( jaux.eq.4 ) then
 260 c
 261             if ( faminf.eq.1 ) then
 262               faminf = fameel(el)
 263             else
 264               if ( fameel(el).ne.faminf ) then
 265                 write (ulsort,texte(langue,7)) 'inf', faminf
 266                 write (ulsort,texte(langue,8))
 267      >          mess14(langue,1,4), el, fameel(el)
 268                 codret = 23
 269               endif
 270             endif
 271 c
 272 c 2.4. ==> si exactement 4 noeuds sont dans l'autre plan, c'est la face
 273 c          superieure. On memorise le numero de famille
 274 c
 275           elseif ( jaux.eq.0 ) then
 276 c
 277             if ( famsup.eq.1 ) then
 278               famsup = fameel(el)
 279             else
 280               if ( fameel(el).ne.famsup ) then
 281                 write (ulsort,texte(langue,7)) 'sup', famsup
 282                 write (ulsort,texte(langue,8))
 283      >          mess14(langue,1,4), el, fameel(el)
 284                 codret = 24
 285               endif
 286             endif
 287 c
 288 c 2.5. ==> sinon, c'est louche ...
 289 c
 290           else
 291 c
 292             write (ulsort,texte(langue,4)) el,(noeele(el,iaux),iaux=1,4)
 293             write (ulsort,texte(langue,5)) jaux
 294             write (ulsort,texte(langue,6)) mess14(langue,1,4),
 295      >                                     '0, 2 ou 4'
 296 c
 297           endif
 298 c
 299         endif
 300 c
 301         endif
 302 c
 303    21 continue
 304 c
 305       if ( codret.eq.0 ) then
 306 c
 307 cgn      write (ulsort,90002) 'nuel2d', nuel2d
 308 cgn      write (ulsort,90002) 'nbse2d', nbse2d
 309       if ( nuel2d.ne.nbse2d ) then
 310         write (ulsort,texte(langue,9)) mess14(langue,3,1), nbse2d
 311         write (ulsort,texte(langue,10))
 312      >    mess14(langue,3,1), nuel2d-nbse2d
 313         codret = 222
 314       endif
 315 c
 316       endif
 317 c
 318 #ifdef _DEBUG_HOMARD_
 319       write (ulsort,texte(langue,7)) 'inf', faminf
 320       write (ulsort,texte(langue,7)) 'sup', famsup
 321 #endif
 322 c
 323 c====
 324 c 3. transformations des hexaedres en quadrangles
 325 c====
 326 #ifdef _DEBUG_HOMARD_
 327       write (ulsort,90002) '3. hexa -> quad ; codret', codret
 328 #endif
 329 c
 330       if ( nbqu2d.ne.0 ) then
 331 c
 332       do 31 , el = 1 , nbelem
 333 c
 334         if ( codret.eq.0 ) then
 335 c
 336         if ( typele(el).eq.edhex8 ) then
 337 c
 338 c 3.1. ==> recherche des noeuds dans le plan cooinf
 339 c
 340 cgn       write (ulsort,90012) 'noeele',el,(noeele(el,iaux),iaux=1,8)
 341           jaux = 0
 342           do 311 , iaux = 1 , 8
 343             if ( nu2dno(noeele(el,iaux)).ne.0 ) then
 344               jaux = jaux + 1
 345               numloc(jaux) = noeele(el,iaux)
 346             endif
 347   311     continue
 348 c
 349 #ifdef _DEBUG_HOMARD_
 350       write (ulsort,90002) 'jaux', jaux
 351       if ( nuel2d.eq.nbse2d.or. nuel2d.eq.988.or.nuel2d.eq.987 ) then
 352       write (ulsort,texte(langue,11)) el, 'hexa',
 353      >                               (numloc(iaux),iaux=1,jaux)
 354       write (ulsort,90004) 'X',(coonoe(numloc(iaux),1),iaux = 1,jaux)
 355       write (ulsort,90004) 'Y',(coonoe(numloc(iaux),2),iaux = 1,jaux)
 356       write (ulsort,90004) 'Z',(coonoe(numloc(iaux),3),iaux = 1,jaux)
 357       endif
 358 #endif
 359 c
 360 c 3.2. ==> si exactement 4 noeuds sont dans ce plan, creation du
 361 c          quadrangle asssocie
 362 c          attention a bien ranger les noeuds pour ne pas croiser !
 363 c          pour cela, on les positionne en fonction de leur milieu :
 364 c
 365 c                                 Y
 366 c                                 .
 367 c                                 .
 368 c                         4...............3
 369 c                         .       .       .
 370 c                         .       .       .
 371 c                         .       .       .
 372 c                         ........M...........> X
 373 c                         .       .       .
 374 c                         .       .       .
 375 c                         .       .       .
 376 c                         1...............2
 377 c
 378 c
 379           if ( jaux.eq.4 ) then
 380 c
 381 c           Le milieu du quadrangle
 382 c
 383             xymil(1) = 0.d0
 384             xymil(2) = 0.d0
 385             do 313 , iaux = 1 , 4
 386               xymil(1) = xymil(1) + coonoe(numloc(iaux),iaux1)
 387               xymil(2) = xymil(2) + coonoe(numloc(iaux),iaux2)
 388   313       continue
 389             xymil(1) = unsqu * xymil(1)
 390             xymil(2) = unsqu * xymil(2)
 391 cgn            write (ulsort,90004) 'xymil', xymil
 392 c
 393 c           Le vecteur entre le milieu et le premier noeud
 394 c
 395             v1(1) = coonoe(numloc(1),iaux1) - xymil(1)
 396             v1(2) = coonoe(numloc(1),iaux2) - xymil(2)
 397             daux1 = sqrt(v1(1)**2+v1(2)**2)
 398             v1(1) = v1(1)/daux1
 399             v1(2) = v1(2)/daux1
 400 cgn            write (ulsort,90004) 'v1', v1
 401 c
 402 c           Les angles entre le segment (milieu,noeud 1) et
 403 c           (milieu,noeud suivant)
 404 c
 405             do 314 , iaux = 2 , 4
 406               vn(1) = coonoe(numloc(iaux),iaux1) - xymil(1)
 407               vn(2) = coonoe(numloc(iaux),iaux2) - xymil(2)
 408               daux1 = sqrt(vn(1)**2+vn(2)**2)
 409               vn(1) = vn(1)/daux1
 410               vn(2) = vn(2)/daux1
 411 cgn            write (ulsort,*) ' '
 412 cgn            write (ulsort,90114) 'vn', iaux,vn
 413               daux1 = v1(1)*vn(1) + v1(2)*vn(2)
 414               daux2 = v1(1)*vn(2) - v1(2)*vn(1)
 415 cgn        write (ulsort,90114) 'p scal', iaux, daux1
 416 cgn        write (ulsort,90114) 'p vect', iaux, daux2
 417 c
 418               if ( (daux1+1.d0).le.zeroma ) then
 419                 daux(iaux) = pi
 420               else
 421                 daux(iaux) = acos(daux1)
 422               endif
 423               if ( daux2.le.0.d0 ) then
 424                 daux(iaux) = deuxpi - daux(iaux)
 425               endif
 426   314       continue
 427 #ifdef _DEBUG_HOMARD_
 428       if ( nuel2d.lt.1 ) then
 429         write (ulsort,90004) 'angles', daux(2), daux(3), daux(4)
 430       endif
 431 #endif
 432 c
 433 c           Classement des angles
 434 c
 435             daux1 = min(daux(2), daux(3), daux(4))
 436             daux2 = max(daux(2), daux(3), daux(4))
 437 cgn        write (ulsort,90004) 'angles min', daux1
 438 cgn        write (ulsort,90004) 'angles max', daux2
 439             tbiaux(1) = 1
 440             do 315 , iaux = 2 , 4
 441               if ( abs(daux(iaux)-daux1).le.epsilo ) then
 442                 tbiaux(2) = iaux
 443               endif
 444               if ( abs(daux(iaux)-daux2).le.epsilo ) then
 445                 tbiaux(4) = iaux
 446               endif
 447   315       continue
 448 c
 449             tbiaux(3) = fp0123(tbiaux(2)-1, tbiaux(4)-1) + 1
 450 cgn      write (ulsort,90002) 'tbiaux final  ', tbiaux
 451 c
 452 c           Transfert de la connectivite
 453 c
 454             nuel2d = nuel2d + 1
 455             do 316 , iaux = 1 , 4
 456               noee2d(nuel2d,iaux) = nu2dno(numloc(tbiaux(iaux)))
 457   316       continue
 458             fame2d(nuel2d) = fameel(el)
 459             type2d(nuel2d) = edqua4
 460 #ifdef _DEBUG_HOMARD_
 461       if ( nuel2d.lt.1 ) then
 462         write (ulsort,90015) 'typele(',nuel2d,') = ',type2d(nuel2d)
 463         do 32221 , iaux = 1 , 4
 464           write (ulsort,90007)
 465      >    '    noeele',nuel2d,iaux,noee2d(nuel2d,iaux)
 466 32221   continue
 467       endif
 468 #endif
 469 c
 470           else
 471             write (ulsort,texte(langue,4)) el,(noeele(el,iaux),iaux=1,8)
 472             write (ulsort,texte(langue,5)) jaux
 473             write (ulsort,texte(langue,6)) mess14(langue,1,9), '4'
 474 c             a changer un jour
 475             codret = 31
 476           endif
 477 c
 478         endif
 479 c
 480         endif
 481 c
 482    31 continue
 483 c
 484       if ( codret.eq.0 ) then
 485 c
 486 cgn      write (ulsort,90002) 'nuel2d', nuel2d
 487 cgn      write (ulsort,90002) 'nbqu2d', nbqu2d
 488 cgn      write (ulsort,90002) 'nuel2d-nbse2d', nuel2d-nbse2d
 489       if ( (nuel2d-nbse2d).ne.nbqu2d ) then
 490         write (ulsort,texte(langue,9))  mess14(langue,3,4), nbqu2d
 491         write (ulsort,texte(langue,10)) mess14(langue,3,4),
 492      >                                  nuel2d-nbse2d
 493         codret = 333
 494       endif
 495 c
 496       endif
 497 c
 498       endif
 499 c
 500 c====
 501 c 4. transformations des pentaedres en triangles
 502 c====
 503 #ifdef _DEBUG_HOMARD_
 504       write (ulsort,90002) '4. pent -> tria ; codret', codret
 505 #endif
 506 c
 507       if ( nbtr2d.ne.0 ) then
 508 c
 509       if ( codret.eq.0 ) then
 510 c
 511       do 41 , el = 1 , nbelem
 512 c
 513         if ( typele(el).eq.edpen6 ) then
 514 c
 515 c 4.1. ==> recherche des noeuds dans le plan cooinf
 516 c
 517           jaux = 0
 518           do 411 , iaux = 1 , 6
 519             if ( nu2dno(noeele(el,iaux)).ne.0 ) then
 520               jaux = jaux + 1
 521               numloc(jaux) = noeele(el,iaux)
 522             endif
 523   411     continue
 524 c
 525 #ifdef _DEBUG_HOMARD_
 526       if ( jaux.gt.0 ) then
 527       write (ulsort,texte(langue,11)) el, 'pent',
 528      >                               (numloc(iaux),iaux=1,jaux)
 529       endif
 530 #endif
 531 c
 532 c 4.2. ==> si exactement 3 noeuds sont dans ce plan, creation du
 533 c          triangle asssocie
 534 c
 535           if ( jaux.eq.3 ) then
 536 c
 537             nuel2d = nuel2d + 1
 538             noee2d(nuel2d,1) = nu2dno(numloc(3))
 539             noee2d(nuel2d,2) = nu2dno(numloc(2))
 540             noee2d(nuel2d,3) = nu2dno(numloc(1))
 541             fame2d(nuel2d) = fameel(el)
 542             type2d(nuel2d) = edtri3
 543 c
 544           else
 545             write (ulsort,texte(langue,4)) el,(noeele(el,iaux),iaux=1,8)
 546             write (ulsort,texte(langue,5)) jaux
 547             write (ulsort,texte(langue,6)) mess14(langue,1,9), '4'
 548 c             a changer un jour
 549             codret = 42
 550           endif
 551 c
 552         endif
 553 c
 554    41 continue
 555 c
 556       if ( codret.eq.0 ) then
 557 c
 558 cgn      write (ulsort,90002) 'nuel2d', nuel2d
 559 cgn      write (ulsort,90002) 'nbtr2d', nbtr2d
 560 cgn      write (ulsort,90002) 'nuel2d-nbse2d-nbqu2d', nuel2d-nbse2d-nbqu2d
 561       if ( (nuel2d-nbse2d-nbqu2d).ne.nbtr2d ) then
 562         write (ulsort,texte(langue,9))  mess14(langue,3,2), nbtr2d
 563         write (ulsort,texte(langue,10)) mess14(langue,3,2),
 564      >                                  nuel2d-nbse2d-nbqu2d
 565         codret = 444
 566       endif
 567 c
 568       endif
 569 c
 570       endif
 571 c
 572       endif
 573 c
 574 c====
 575 c 5. la fin
 576 c====
 577 c
 578 #ifdef _DEBUG_HOMARD_
 579       write (ulsort,90002) '5. fin ; codret', codret
 580 #endif
 581 c
 582       if ( codret.ne.0 ) then
 583 c
 584 #include "envex2.h"
 585 c
 586       write (ulsort,texte(langue,1)) 'Sortie', nompro
 587       write (ulsort,texte(langue,2)) codret
 588       endif
 589 c
 590 #ifdef _DEBUG_HOMARD_
 591       write (ulsort,texte(langue,1)) 'Sortie', nompro
 592       call dmflsh (iaux)
 593 #endif
 594 c
 595       end