Télécharger la plaquette Cast3M


fr en

Télécharger facaxi.eso

Retour à la liste

Numérotation des lignes : 

  1. C FACAXI    SOURCE    CB215821  17/11/30    21:16:07     9639           
  2.       SUBROUTINE FACAXI (MYMOD,INFOEL,NPAX,NGAX,KACHE,INOR,ICHFAC
  3.      &     ,EXTINC)
  4.       IMPLICIT INTEGER(I-N)
  5.       IMPLICIT REAL*8 (A-H,O-Z)
  6. C----------------------------------------------------------------------
  7. C SP appele par FFORME
  8. C
  9. C     Calcul des facteurs de forme en axisymetrique
  10. C   entée:
  11. C     MYMOD  : pointeur de l'objet modèle
  12. C     INFOEL : utile pour les coques ou quadratiques
  13. C     NPAX   : nombre de points d integration (disc.reguliere)
  14. C     NGAX   : nombre de points de Gauss
  15. C     KACHE  : 0 si option convexe, sinon option cache
  16. C     INOR   : 1 si normalisation, 0 sinon
  17. C     EXTINC : coefficient d'extinction (si cavite absorbante)
  18. C   sortie:
  19. C     ICHFAC : pointeur sur l'objet MCHAML resultat
  20. C
  21. C----------------------------------------------------------------------
  22. C     traitement des coques par dedoublement des elements a partir
  23. C     de la normale
  24. C                          ->
  25. C       A (inverse) = A - e*n  (e=1e-3)
  26. C     cas des boucles 1 sur k1 et 2 sur k2
  27. C     mais pas de la boucle 3 obstructeurs
  28. C
  29. C     bcl face k1
  30. C      ** face k1  **
  31. C            bcl face k2
  32. C            .. face k2 ..
  33. C                  bcl 3 obstructeurs
  34. C            .. si coq: inverse face k2 ..
  35. C                  les obstructeurs sont les memes que pour k2
  36. C            fin bcl face k2
  37. C
  38. C      ** si coq : inverse face k1 **
  39. C            bcl face k2
  40. C            .. face k2 ..
  41. C                  bcl 3 obstructeurs
  42. C            .. si coq: inverse face k2 ..
  43. C                  les obstructeurs sont les memes que pour k2
  44. C            fin bcl face k2
  45. C     fin bcl face k1
  46. C----------------------------------------------------------------------
  47. C
  48.       LOGICAL ICOQ
  49. C
  50. -INC CCREEL
  51. -INC SMELEME
  52. -INC SMMODEL
  53. -INC CCOPTIO
  54. -INC SMCOORD
  55.       POINTEUR MYMOD.MMODEL
  56.  
  57. C----------------------------------------------------------------------
  58.       SEGMENT , INFOEL
  59.       LOGICAL KCOQ(N1),KQUAD(N1)
  60.       ENDSEGMENT
  61. C----------------------------------------------------------------------
  62. C     FACTEURS DE FORME
  63. C           NNBEL1 = NOMBRE DE LIGNES + 1
  64. C           NBEL2  = NOMBRE DE COLONNES
  65. C     LFACT(NNBEL1) POINTE SUR LE TABLEAU DES SURFACES
  66. C
  67.       SEGMENT IFACFO
  68.       INTEGER   LFACT(NNBEL1)
  69.       ENDSEGMENT
  70.       SEGMENT LFAC
  71.       REAL*8    FACT(NBEL2)
  72.       ENDSEGMENT
  73.       POINTEUR PSUR.LFAC, PLIG.LFAC
  74.       POINTEUR PCOL.LFAC
  75. C----------------------------------------------------------------------
  76. C     coordonnees des obstructeurs
  77.       SEGMENT SFOBS
  78.       REAL*8 OBS(2,NTOBS)
  79.       ENDSEGMENT
  80. C----------------------------------------------------------------------
  81.       SEGMENT STRAV
  82.       REAL*8 A1(NA,NA),A2(NA,NA),A3(NA,NA),AA2(NA,NA)
  83.       REAL*8 U1(NA1),U2(NA1),UU2(NA1)
  84.       ENDSEGMENT
  85. C----------------------------------------------------------------------
  86. C
  87.       EPS = 1D-5
  88.  
  89.       KIMP = IIMPI
  90.       NES = IDIM
  91.  
  92. C... critere de dedoublement des coques
  93.       ECOQ=1.D-3
  94.       IF (INFOEL.EQ.0) THEN
  95.          ICOQ = .FALSE.
  96.       ELSE
  97.          ICOQ = .TRUE.
  98.          SEGACT INFOEL
  99.       ENDIF
  100. C...  quadratique
  101.       NSPA1=1
  102.       NSPA2=1
  103.       NSPA3=1
  104.       NS=2
  105.  
  106.       RAD = 0
  107. C
  108. C     Calcul du nombre de faces NFACE et du rayon RAD
  109.       SEGACT MYMOD
  110.       NTYP = MYMOD.KMODEL(/1)
  111.       NFACE = 0
  112.       DO 10 ITYP=1,NTYP
  113.          IMODEL = MYMOD.KMODEL(ITYP)
  114.          SEGACT IMODEL
  115.          IPT1 = IMAMOD
  116.          SEGDES IMODEL
  117.          SEGACT IPT1
  118.          NEL = IPT1.NUM(/2)
  119.          NSGEO  = IPT1.NUM(/1)
  120. C         Recherche du max sur Ox
  121.          DO 5 IEL = 1,NEL
  122.             DO 6 IPT = 1,NSGEO
  123.                IREF = (IDIM+1)*(IPT1.NUM(IPT,IEL)-1)
  124.                VALX = XCOOR(IREF+1)
  125.                IF (VALX.GT.RAD)  RAD = VALX
  126.  6          CONTINUE
  127.  5       CONTINUE
  128.          IF (ICOQ.AND.KCOQ(ITYP)) THEN
  129.             NFACE = NFACE + 2 * NEL
  130.          ELSE
  131.             NFACE = NFACE + NEL
  132.          ENDIF
  133.  
  134.  10   CONTINUE
  135.  
  136. C
  137.  
  138.       IF (KIMP.GE.3) THEN
  139.          WRITE( 6,*) ' DIMENSIONNEMENT : RAD =',RAD
  140.          WRITE( 6,*) ' NOMBRE TOTAL DE FACES ',NFACE
  141.       ENDIF
  142. C
  143. C>>>  INITIALISATION SFOBS et STRAV
  144. C     -----------------------------
  145.       IF(KACHE.NE.0) THEN
  146.          NTOBS = 2*NFACE
  147.          SEGINI SFOBS
  148.       ENDIF
  149.       NA=2
  150.       NA1=3
  151.       SEGINI STRAV
  152.  
  153. C
  154. C>>>  INITIALISATION OBJET FACFOR
  155. C     ---------------------------
  156. C
  157.       NNBEL1=NFACE+1
  158.       NBEL2=NFACE
  159.       SEGINI IFACFO
  160.       DO 900 LS=1,NNBEL1
  161.          SEGINI PLIG
  162.          LFACT(LS)=PLIG
  163.          SEGACT PLIG*MOD
  164.  900  CONTINUE
  165.       PSUR = LFACT(NNBEL1)
  166.       SEGACT PSUR*MOD
  167.  
  168. C     -------------------------------------------------------------
  169. C>>>  CALCUL
  170.  
  171. C
  172. C>>   BOUCLE FACE 1
  173. C
  174.  
  175.       NELT1= 0
  176.       DO 100 ITYP = 1,NTYP
  177.          IMODEL = MYMOD.KMODEL(ITYP)
  178.          SEGACT IMODEL
  179.          IPT1 = IMAMOD
  180.          SEGDES IMODEL
  181.          IF(KIMP.GE.4) THEN
  182.             WRITE (6,*) 'Maillage :',IPT1
  183.          ENDIF
  184.  
  185.          NSGEO1 = IPT1.NUM(/1)
  186.          NSPA1=1
  187.          IF(ICOQ.AND.KQUAD(ITYP)) NSPA1=2
  188.          NEL1 = IPT1.NUM(/2)
  189.          DO 110 I1 = 1,NEL1
  190.  
  191.             IF (ICOQ.AND.KCOQ(ITYP)) THEN
  192.                K1 = (2*I1-1) + NELT1
  193.             ELSE
  194.                K1 = I1+ NELT1
  195.             ENDIF
  196.  
  197. C*** traitement de la face K1 ***************************************
  198.  
  199.             PLIG=LFACT(K1)
  200. C
  201.             DO 111 IS = 1,NSGEO1,NSPA1
  202.                LS=IS
  203.                IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2
  204.                IREF = (IDIM+1)*(IPT1.NUM(IS,I1)-1)
  205.                IF(KIMP.GE.4) THEN
  206.                   WRITE (6,*) 'Noeud numéro',IPT1.NUM(IS,I1),'ref :'
  207.      $                 ,IREF
  208.                ENDIF
  209.                DO 112 K = 1,NES
  210. C**          A1(K,IS) = XCOOR(IREF+K)/RAD
  211.                   A1(K,LS) = XCOOR(IREF+K)/RAD
  212.  112           CONTINUE
  213.  111        CONTINUE
  214.             CALL KNORM2(NES,A1,NS,U1,S1)
  215.             S1=XPI*S1*ABS(A1(1,1)+A1(1,2))
  216.             PSUR.FACT(K1)=S1
  217.             ZMIN1=MIN(A1(2,1),A1(2,2))
  218.             ZMAX1= MAX(A1(2,1),A1(2,2))
  219.             RMAX1= MAX(A1(1,1),A1(1,2))
  220.  
  221.  
  222. C>>    BOUCLE FACE 2
  223. C     -------------------------------------------------------------
  224.  
  225.             NELT2= 0
  226.             DO 200 JTYP = 1,NTYP
  227.                IMODEL = MYMOD.KMODEL(JTYP)
  228.                SEGACT IMODEL
  229.                IPT2 = IMAMOD
  230.                SEGDES IMODEL
  231.  
  232.                NSGEO2 = IPT2.NUM(/1)
  233.                NSPA2=1
  234.                IF(ICOQ.AND.KQUAD(JTYP)) NSPA2=2
  235.                NEL2 = IPT2.NUM(/2)
  236.                DO 210 I2 = 1,NEL2
  237.  
  238.                   IF (ICOQ.AND.KCOQ(JTYP))THEN
  239.                      K2 = 2*I2-1  + NELT2
  240.                   ELSE
  241.                      K2 = I2 + NELT2
  242.                   ENDIF
  243.  
  244. C... face K2 .....................................................
  245.  
  246.                   KVU=1
  247.                   FF=0.D0
  248.  
  249. C.. UTILISATION DE LA RECIPROCITE
  250. C
  251.                   IF(K1.GT.K2) THEN
  252.                      S2=PSUR.FACT(K2)
  253.                      PCOL=LFACT(K2)
  254.  
  255.                      PLIG.FACT(K2)=(S2/S1)*PCOL.FACT(K1)
  256.  
  257.                   ELSE
  258. C.. TEST K1=K2 ET VISIBILITE A PRIORI
  259.  
  260. C------------------------------------------------------------------
  261.  
  262. C>> K1=K2 : ELIMINATION DES CAS TRIVIAUX OU F(K1,K2)=0
  263. C
  264.                      IF(K1.EQ.K2) THEN
  265.                         IF(KIMP.GE.4) THEN
  266.                            WRITE(6,*) ' A1 ',A1(1,1),A1(2,1)
  267.                            WRITE(6,*) ' A1 ',A1(1,2),A1(2,2)
  268.                            WRITE(6,*) ' U1 ',U1(1),U1(2)
  269.                         ENDIF
  270.                         IF (ABS(U1(1)).LT.EPS) KVU=0
  271.                         IF (U1(1).GE.EPS) KVU=0
  272.  
  273.                         IF(KVU.NE.0) THEN
  274.                            DO 214 K = 1,NES
  275.                               A2(K,1) = A1(K,1)
  276.                               A2(K,2) = A1(K,2)
  277.  214                       CONTINUE
  278.                         ENDIF
  279.  
  280. C>> K1::K2 : ELIMINATION DES CAS DE NON-VISIBLITE
  281. C            UTILISATION DE LA VISIBILITE A PRIORI DANS LE PLAN R-Z
  282. C            ON TESTE LA FACE K2 ET SON SYMETRIQUE (-R,Z)
  283.  
  284.                      ELSE
  285.                         DO 211 IS = 1,NSGEO2,NSPA2
  286.                            LS=IS
  287.                            IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2
  288.                            IREF = (IDIM+1)*(IPT2.NUM(IS,I2)-1)
  289.                            DO 212 K = 1,NES
  290. C**        A2(K,IS) = XCOOR(IREF+K)/RAD
  291.                               A2(K,LS) = XCOOR(IREF+K)/RAD
  292.  212                       CONTINUE
  293.  211                    CONTINUE
  294.                         CALL KNORM2(NES,A2,NS,U2,S2)
  295.                         S2=XPI*S2*ABS(A2(1,1)+A2(1,2))
  296.  
  297. C.. calcul du symetrique
  298.                         AA2(1,1)= -A2(1,2)
  299.                         AA2(2,1)=  A2(2,2)
  300.                         AA2(1,2)= -A2(1,1)
  301.                         AA2(2,2)=  A2(2,1)
  302.                         CALL KNORM2(NES,AA2,NS,UU2,SS2)
  303.  
  304. C.. visibilite a priori
  305.                         CALL KPRIOR(NES,NS,NS,A1,A2,U1,U2,IVU)
  306.                         CALL KPRIOR(NES,NS,NS,A1,AA2,U1,UU2,IVUS)
  307.  
  308.                         IF(KIMP.GE.4)WRITE(6,*) ' K1 K2 IVU IVUS ',K1,K2
  309.      $                       ,IVU,IVUS
  310.                         IF(IVU.EQ.0.AND.IVUS.EQ.0) KVU=0
  311.  
  312.                      ENDIF
  313. C   ---------------------------------------------------------------
  314. C<< FIN TEST K1=K2 ET VISIBIITE A PRIORI
  315.                      IF(KIMP.GE.4.AND.KVU.NE.0) THEN
  316.                         WRITE(6,*) ' FACES VISIBLES ',K1,K2
  317.                      ENDIF
  318.  
  319. C>> CALCUL
  320. C... option convexe
  321.                      IF(KACHE.EQ.0) THEN
  322.                         IF(KVU.NE.0) THEN
  323.  
  324.                            CALL KAXC(A1,A2,NPAX,NGAX,FF,KIMP,EXTINC,RAD)
  325.                            IF(KIMP.GE.4) WRITE(6,*) ' FF CONVEXE ',K1,K2
  326.      $                          ,FF/S1
  327.  
  328.                         ENDIF
  329. C.. option cache
  330.                      ELSE
  331.                         IF(KVU.NE.0) THEN
  332.  
  333. C>> RECHERCHE DES OBSTRUCTEURS POTENTIELS.--------------------------
  334.  
  335.                            NOBS=0
  336.  
  337. C>> BOUCLE FACE 3---------------------------------------------------
  338.                            NELT3= 0
  339.                            DO 300 KTYP = 1,NTYP
  340.                               IMODEL = MYMOD.KMODEL(KTYP)
  341.                               SEGACT IMODEL
  342.                               IPT3 = IMAMOD
  343.                               SEGDES IMODEL
  344.  
  345.                               NSGEO3 = IPT3.NUM(/1)
  346.                               NSPA3=1
  347.                               IF(ICOQ.AND.KQUAD(KTYP)) NSPA3=2
  348.                               NEL3 = IPT3.NUM(/2)
  349.                               DO 310 I3 = 1,NEL3
  350.                                  K3 = I3+ NELT3
  351.  
  352.                                  IF(K3.NE.K1.AND.K3.NE.K2) THEN
  353.  
  354.                                     DO 311 IS = 1,NSGEO3,NSPA3
  355.                                        LS=IS
  356.                                        IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1
  357.      $                                      )/2
  358.                                        IREF = (IDIM+1)*(IPT3.NUM(IS,I3)
  359.      $                                      -1)
  360.                                        DO 312 K = 1,NES
  361. C**        A3(K,IS) = XCOOR(IREF+K)/RAD
  362.                                           A3(K,LS) = XCOOR(IREF+K)/RAD
  363.  312                                   CONTINUE
  364.  311                                CONTINUE
  365.  
  366.                                     ZMIN2=MIN(A2(2,1),A2(2,2))
  367.                                     ZMAX2= MAX(A2(2,1),A2(2,2))
  368.                                     RMAX2= MAX(A2(1,1),A2(1,2))
  369.                                     ZTMIN=MIN(ZMIN1,ZMIN2)
  370.                                     ZTMAX= MAX(ZMAX1,ZMAX2)
  371.                                     RMAX = MAX(RMAX1,RMAX2)
  372.  
  373.                                     IF( MAX(A3(2,1),A3(2,2)).LE.ZTMIN)
  374.      $                                   THEN
  375.                                        KOBS=0
  376.                                     ELSEIF(MIN(A3(2,1),A3(2,2)).GE
  377.      $                                      .ZTMAX) THEN
  378.                                        KOBS=0
  379.                                     ELSEIF(MIN(A3(1,1),A3(1,2)).GE
  380.      $                                      .RMAX) THEN
  381.                                        KOBS=0
  382.                                     ELSE
  383.                                        KOBS=1
  384.                                        OBS(1,2*NOBS+1)=A3(1,1)
  385.                                        OBS(2,2*NOBS+1)=A3(2,1)
  386.                                        OBS(1,2*NOBS+2)=A3(1,2)
  387.                                        OBS(2,2*NOBS+2)=A3(2,2)
  388.                                        NOBS=NOBS+1
  389.  
  390.                                     ENDIF
  391.                                  ENDIF
  392.  
  393.  310                          CONTINUE
  394. C
  395.                               NELT3=NELT3+NEL3
  396. C
  397.  300                       CONTINUE
  398.  
  399. C<< FIN RECHERCHE OBSTRUCTEURS--------------------------------------
  400.  
  401.                            IF(KIMP.GE.4) THEN
  402.                               WRITE(6,*) ' FACES K1 K2 ',K1,K2,' NOBS '
  403.      $                             ,NOBS
  404.                            ENDIF
  405.  
  406.                            CALL KAXK(A1,A2,OBS,NOBS,NTOBS,NPAX,NGAX,FF
  407.      $                          ,KIMP,EXTINC,RAD)
  408.                            IF(KIMP.GE.4) WRITE(6,*) ' FF  ',K1,K2,FF/S1
  409.  
  410.                         ENDIF
  411.  
  412.                      ENDIF
  413. C<< FIN CALCUL
  414. C   ---------------------------------------------------------------
  415.  
  416. C       WRITE(6,*) ' K1 K2 PLIG ',K1,K2,PLIG
  417.                      PLIG.FACT(K2) =  FF/S1
  418.  
  419.                   ENDIF
  420. C... fin face K2 ..................................................
  421.  
  422.  
  423. C... face inverse de K2 (cas des coques) ..........................
  424.  
  425.                   IF (ICOQ.AND.KCOQ(JTYP)) THEN
  426.  
  427.                      K2=K2 + 1
  428.  
  429.                      KVU=1
  430.                      FF=0.D0
  431.  
  432. C UTILISATION DE LA RECIPROCITE
  433. C
  434.                      IF(K1.GT.K2) THEN
  435.                         S2=PSUR.FACT(K2)
  436.                         PCOL=LFACT(K2)
  437.  
  438.                         PLIG.FACT(K2)=(S2/S1)*PCOL.FACT(K1)
  439.  
  440.                      ELSE
  441.  
  442. C>> TEST K1=K2 ET VISIBILITE A PRIORI
  443.  
  444. C------------------------------------------------------------------
  445.  
  446. C>> K1=K2 : ELIMINATION DES CAS TRIVIAUX OU F(K1,K2)=0
  447. C
  448.                         IF(K1.EQ.K2) THEN
  449.                            IF(KIMP.GE.4) THEN
  450.                               WRITE(6,*) ' A1 ',A1(1,1),A1(2,1)
  451.                               WRITE(6,*) ' A1 ',A1(1,2),A1(2,2)
  452.                               WRITE(6,*) ' U1 ',U1(1),U1(2)
  453.                            ENDIF
  454.                            IF (ABS(U1(1)).LT.EPS) KVU=0
  455.                            IF (U1(1).GE.EPS) KVU=0
  456.  
  457.                            IF(KVU.NE.0) THEN
  458.                               DO 514 K = 1,NES
  459.                                  A2(K,1) = A1(K,1)
  460.                                  A2(K,2) = A1(K,2)
  461.  514                          CONTINUE
  462.                            ENDIF
  463.  
  464. C>> K1::K2 : ELIMINATION DES CAS DE NON-VISIBLITE
  465. C            UTILISATION DE LA VISIBILITE A PRIORI DANS LE PLAN R-Z
  466. C            ON TESTE LA FACE K2 ET SON SYMETRIQUE (-R,Z)
  467.  
  468.  
  469.                         ELSE
  470.  
  471.                            DO 512 K = 1,NES
  472.                               U2(K)=-U2(K)
  473.                               XX1= A2(K,1)
  474.                               A2(K,1) = A2(K,2) + ECOQ*U2(K)
  475.                               A2(K,2) = XX1     + ECOQ*U2(K)
  476.  512                       CONTINUE
  477.  
  478. C... calcul du symetrique
  479.                            AA2(1,1)= -A2(1,2)
  480.                            AA2(2,1)=  A2(2,2)
  481.                            AA2(1,2)= -A2(1,1)
  482.                            AA2(2,2)=  A2(2,1)
  483.                            CALL KNORM2(NES,AA2,NS,UU2,SS2)
  484.  
  485. C... visibilité a priori
  486.                            CALL KPRIOR(NES,NS,NS,A1,A2,U1,U2,IVU)
  487.                            CALL KPRIOR(NES,NS,NS,A1,AA2,U1,UU2,IVUS)
  488.  
  489.                            IF(KIMP.GE.4)WRITE(6,*) ' K1 K2 IVU IVUS ',K1
  490.      $                          ,K2,IVU,IVUS
  491.                            IF(IVU.EQ.0.AND.IVUS.EQ.0) KVU=0
  492.  
  493.                         ENDIF
  494. C   ---------------------------------------------------------------
  495. C<< FIN TEST K1=K2 ET VISIBIITE A PRIORI
  496.                         IF(KIMP.GE.4.AND.KVU.NE.0) THEN
  497.                            WRITE(6,*) ' FACES VISIBLES ',K1,K2
  498.                         ENDIF
  499.  
  500. C>> CALCUL --------------------------------------------------------
  501. C.. option convexe
  502.                         IF(KACHE.EQ.0) THEN
  503.                            IF(KVU.NE.0) THEN
  504.  
  505.                               CALL KAXC(A1,A2,NPAX,NGAX,FF,KIMP,EXTINC
  506.      $                             ,RAD)
  507.                               IF(KIMP.GE.4) WRITE(6,*) ' FF CONVEXE ',K1
  508.      $                             ,K2,FF/S1
  509.  
  510.                            ENDIF
  511. C.. option cache
  512.                         ELSE
  513.                            IF(KVU.NE.0) THEN
  514.  
  515. C>> les obstructeurs potentiels sont les memes que pour la face K2
  516.  
  517.                               IF(KIMP.GE.4) THEN
  518.                                  WRITE(6,*) ' FACES K1 K2 ',K1,K2
  519.      $                                ,' NOBS ',NOBS
  520.                               ENDIF
  521.  
  522.                               CALL KAXK(A1,A2,OBS,NOBS,NTOBS,NPAX,NGAX
  523.      $                             ,FF,KIMP,EXTINC,RAD)
  524.  
  525.  
  526.                               IF(KIMP.GE.4) WRITE(6,*) ' FF  ',K1,K2,FF
  527.      $                             /S1
  528.  
  529.                            ENDIF
  530.  
  531.                         ENDIF
  532. C<< FIN CALCUL ------------------------------------------------------
  533.  
  534. C       WRITE(6,*) ' K1 K2 PLIG ',K1,K2,PLIG
  535.                         PLIG.FACT(K2) =  FF/S1
  536.  
  537.                      ENDIF
  538.  
  539.                   ENDIF
  540.  
  541. C... fin face inverse de K2 (cas des coques) .........................
  542.  
  543.  210           CONTINUE
  544.  
  545.                IF (ICOQ.AND.KCOQ(JTYP)) THEN
  546.                   NELT2 = NELT2 + 2 * NEL2
  547.                ELSE
  548.                   NELT2 = NELT2 + NEL2
  549.                ENDIF
  550.  
  551.  200        CONTINUE
  552.  
  553. C     ---FIN BOUCLE FACE 2-----------------------------------------
  554.  
  555.  
  556.  
  557.  
  558. C*** fin traitement de la face K1 ***********************************
  559.  
  560. C*** traitement de la face inverse de K1*****************************
  561.             IF (ICOQ.AND.KCOQ(ITYP)) THEN
  562.  
  563.                K1=K1+1
  564.                PLIG=LFACT(K1)
  565.  
  566.                DO 122 K = 1,NES
  567.                   U1(K) =-U1(K)
  568.                   XX1=A1(K,1)
  569.                   A1(K,1) = A1(K,2) + ECOQ * U1(K)
  570.                   A1(K,2) = XX1     + ECOQ * U1(K)
  571.  122           CONTINUE
  572.  
  573.                PSUR.FACT(K1)=S1
  574.  
  575. C>>    BOUCLE FACE 2
  576. C     -------------------------------------------------------------
  577.  
  578.                NELT2= 0
  579.                DO 2000 JTYP = 1,NTYP
  580.                   IMODEL = MYMOD.KMODEL(JTYP)
  581.                   SEGACT IMODEL
  582.                   IPT2 = IMAMOD
  583.                   SEGDES IMODEL
  584.  
  585.                   NSGEO2 = IPT2.NUM(/1)
  586.                   NSPA2=1
  587.                   IF(ICOQ.AND.KQUAD(JTYP)) NSPA2=2
  588.                   NEL2 = IPT2.NUM(/2)
  589.                   DO 2100 I2 = 1,NEL2
  590.  
  591.                      IF (ICOQ.AND.KCOQ(JTYP))THEN
  592.                         K2 = 2*I2-1  + NELT2
  593.                      ELSE
  594.                         K2 = I2 + NELT2
  595.                      ENDIF
  596.  
  597. C... face K2 .....................................................
  598.  
  599.                      KVU=1
  600.                      FF=0.D0
  601.  
  602. C.. UTILISATION DE LA RECIPROCITE
  603. C
  604.                      IF(K1.GT.K2) THEN
  605.                         S2=PSUR.FACT(K2)
  606.                         PCOL=LFACT(K2)
  607.  
  608.                         PLIG.FACT(K2)=(S2/S1)*PCOL.FACT(K1)
  609.  
  610.                      ELSE
  611. C.. TEST K1=K2 ET VISIBILITE A PRIORI
  612.  
  613. C------------------------------------------------------------------
  614.  
  615. C>> K1=K2 : ELIMINATION DES CAS TRIVIAUX OU F(K1,K2)=0
  616. C
  617.                         IF(K1.EQ.K2) THEN
  618.                            IF(KIMP.GE.4) THEN
  619.                               WRITE(6,*) ' A1 ',A1(1,1),A1(2,1)
  620.                               WRITE(6,*) ' A1 ',A1(1,2),A1(2,2)
  621.                               WRITE(6,*) ' U1 ',U1(1),U1(2)
  622.                            ENDIF
  623.                            IF (ABS(U1(1)).LT.EPS) KVU=0
  624.                            IF (U1(1).GE.EPS) KVU=0
  625.  
  626.                            IF(KVU.NE.0) THEN
  627.                               DO 2140 K = 1,NES
  628.                                  A2(K,1) = A1(K,1)
  629.                                  A2(K,2) = A1(K,2)
  630.  2140                         CONTINUE
  631.                            ENDIF
  632.  
  633. C>> K1::K2 : ELIMINATION DES CAS DE NON-VISIBLITE
  634. C            UTILISATION DE LA VISIBILITE A PRIORI DANS LE PLAN R-Z
  635. C            ON TESTE LA FACE K2 ET SON SYMETRIQUE (-R,Z)
  636.  
  637.                         ELSE
  638.                            DO 2110 IS = 1,NSGEO2,NSPA2
  639.                               LS=IS
  640.                               IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2
  641.                               IREF = (IDIM+1)*(IPT2.NUM(IS,I2)-1)
  642.                               DO 2120 K = 1,NES
  643. C**        A2(K,IS) = XCOOR(IREF+K)/RAD
  644.                                  A2(K,LS) = XCOOR(IREF+K)/RAD
  645.  2120                         CONTINUE
  646.  2110                      CONTINUE
  647.                            CALL KNORM2(NES,A2,NS,U2,S2)
  648.                            S2=XPI*S2*ABS(A2(1,1)+A2(1,2))
  649.  
  650. C.. calcul du symetrique
  651.                            AA2(1,1)= -A2(1,2)
  652.                            AA2(2,1)=  A2(2,2)
  653.                            AA2(1,2)= -A2(1,1)
  654.                            AA2(2,2)=  A2(2,1)
  655.                            CALL KNORM2(NES,AA2,NS,UU2,SS2)
  656.  
  657. C.. visibilite a priori
  658.                            CALL KPRIOR(NES,NS,NS,A1,A2,U1,U2,IVU)
  659.                            CALL KPRIOR(NES,NS,NS,A1,AA2,U1,UU2,IVUS)
  660.  
  661.                            IF(KIMP.GE.4)WRITE(6,*) ' K1 K2 IVU IVUS ',K1
  662.      $                          ,K2,IVU,IVUS
  663.                            IF(IVU.EQ.0.AND.IVUS.EQ.0) KVU=0
  664.  
  665.                         ENDIF
  666. C   ---------------------------------------------------------------
  667. C<< FIN TEST K1=K2 ET VISIBIITE A PRIORI
  668.                         IF(KIMP.GE.4.AND.KVU.NE.0) THEN
  669.                            WRITE(6,*) ' FACES VISIBLES ',K1,K2
  670.                         ENDIF
  671.  
  672. C>> CALCUL
  673. C... option convexe
  674.                         IF(KACHE.EQ.0) THEN
  675.                            IF(KVU.NE.0) THEN
  676.  
  677.                               CALL KAXC(A1,A2,NPAX,NGAX,FF,KIMP,EXTINC
  678.      $                             ,RAD)
  679.                               IF(KIMP.GE.4) WRITE(6,*) ' FF CONVEXE ',K1
  680.      $                             ,K2,FF/S1
  681.  
  682.                            ENDIF
  683. C.. option cache
  684.                         ELSE
  685.                            IF(KVU.NE.0) THEN
  686.  
  687. C>> RECHERCHE DES OBSTRUCTEURS POTENTIELS.--------------------------
  688.  
  689.                               NOBS=0
  690.  
  691. C>> BOUCLE FACE 3---------------------------------------------------
  692.                               NELT3= 0
  693.                               DO 3000 KTYP = 1,NTYP
  694.                                  IMODEL = MYMOD.KMODEL(KTYP)
  695.                                  SEGACT IMODEL
  696.                                  IPT3 = IMAMOD
  697.                                  SEGDES IMODEL
  698.  
  699.                                  NSGEO3 = IPT3.NUM(/1)
  700.                                  NSPA3=1
  701.                                  IF(ICOQ.AND.KQUAD(KTYP)) NSPA3=2
  702.                                  NEL3 = IPT3.NUM(/2)
  703.                                  DO 3100 I3 = 1,NEL3
  704.                                     K3 = I3+ NELT3
  705.  
  706.                                     IF(K3.NE.K1.AND.K3.NE.K2) THEN
  707.  
  708.                                        DO 3110 IS = 1,NSGEO3,NSPA3
  709.                                           LS=IS
  710.                                           IF(ICOQ.AND.KQUAD(ITYP)) LS
  711.      $                                         =(IS+1)/2
  712.                                           IREF = (IDIM+1)*(IPT3.NUM(IS
  713.      $                                         ,I3)-1)
  714.                                           DO 3120 K = 1,NES
  715. C**        A3(K,IS) = XCOOR(IREF+K)/RAD
  716.                                              A3(K,LS) = XCOOR(IREF+K)
  717.      $                                            /RAD
  718.  3120                                     CONTINUE
  719.  3110                                  CONTINUE
  720.  
  721.                                        ZMIN2=MIN(A2(2,1),A2(2,2))
  722.                                        ZMAX2= MAX(A2(2,1),A2(2,2))
  723.                                        RMAX2= MAX(A2(1,1),A2(1,2))
  724.                                        ZTMIN=MIN(ZMIN1,ZMIN2)
  725.                                        ZTMAX= MAX(ZMAX1,ZMAX2)
  726.                                        RMAX = MAX(RMAX1,RMAX2)
  727.  
  728.                                        IF( MAX(A3(2,1),A3(2,2)).LE
  729.      $                                      .ZTMIN) THEN
  730.                                           KOBS=0
  731.                                        ELSEIF(MIN(A3(2,1),A3(2,2)).GE
  732.      $                                         .ZTMAX) THEN
  733.                                           KOBS=0
  734.                                        ELSEIF(MIN(A3(1,1),A3(1,2)).GE
  735.      $                                         .RMAX) THEN
  736.                                           KOBS=0
  737.                                        ELSE
  738.                                           KOBS=1
  739.                                           OBS(1,2*NOBS+1)=A3(1,1)
  740.                                           OBS(2,2*NOBS+1)=A3(2,1)
  741.                                           OBS(1,2*NOBS+2)=A3(1,2)
  742.                                           OBS(2,2*NOBS+2)=A3(2,2)
  743.                                           NOBS=NOBS+1
  744.  
  745.                                        ENDIF
  746.                                     ENDIF
  747.  
  748.  3100                            CONTINUE
  749. C
  750.                                  NELT3=NELT3+NEL3
  751. C
  752.  3000                         CONTINUE
  753.  
  754. C<< FIN RECHERCHE OBSTRUCTEURS--------------------------------------
  755.  
  756.                               IF(KIMP.GE.4) THEN
  757.                                  WRITE(6,*) ' FACES K1 K2 ',K1,K2
  758.      $                                ,' NOBS ',NOBS
  759.                               ENDIF
  760.  
  761.                               CALL KAXK(A1,A2,OBS,NOBS,NTOBS,NPAX,NGAX
  762.      $                             ,FF,KIMP,EXTINC,RAD)
  763.                               IF(KIMP.GE.4) WRITE(6,*) ' FF  ',K1,K2,FF
  764.      $                             /S1
  765.  
  766.                            ENDIF
  767.  
  768.                         ENDIF
  769. C<< FIN CALCUL
  770. C   ---------------------------------------------------------------
  771.  
  772. C       WRITE(6,*) ' K1 K2 PLIG ',K1,K2,PLIG
  773.                         PLIG.FACT(K2) =  FF/S1
  774.  
  775.                      ENDIF
  776. C... fin face K2 ..................................................
  777.  
  778.  
  779. C... face inverse de K2 (cas des coques) ..........................
  780.  
  781.                      IF (ICOQ.AND.KCOQ(JTYP)) THEN
  782.  
  783.                         K2=K2 + 1
  784.  
  785.                         KVU=1
  786.                         FF=0.D0
  787.  
  788. C UTILISATION DE LA RECIPROCITE
  789. C
  790.                         IF(K1.GT.K2) THEN
  791.                            S2=PSUR.FACT(K2)
  792.                            PCOL=LFACT(K2)
  793.  
  794.                            PLIG.FACT(K2)=(S2/S1)*PCOL.FACT(K1)
  795.  
  796.                         ELSE
  797. C**
  798.  
  799. C>> TEST K1=K2 ET VISIBILITE A PRIORI
  800.  
  801. C------------------------------------------------------------------
  802.  
  803. C>> K1=K2 : ELIMINATION DES CAS TRIVIAUX OU F(K1,K2)=0
  804. C
  805.                            IF(K1.EQ.K2) THEN
  806.                               IF(KIMP.GE.4) THEN
  807.                                  WRITE(6,*) ' A1 ',A1(1,1),A1(2,1)
  808.                                  WRITE(6,*) ' A1 ',A1(1,2),A1(2,2)
  809.                                  WRITE(6,*) ' U1 ',U1(1),U1(2)
  810.                               ENDIF
  811.                               IF (ABS(U1(1)).LT.EPS) KVU=0
  812.                               IF (U1(1).GE.EPS) KVU=0
  813.  
  814.                               IF(KVU.NE.0) THEN
  815.                                  DO 5140 K = 1,NES
  816.                                     A2(K,1) = A1(K,1)
  817.                                     A2(K,2) = A1(K,2)
  818.  5140                            CONTINUE
  819.                               ENDIF
  820.  
  821. C>> K1::K2 : ELIMINATION DES CAS DE NON-VISIBLITE
  822. C            UTILISATION DE LA VISIBILITE A PRIORI DANS LE PLAN R-Z
  823. C            ON TESTE LA FACE K2 ET SON SYMETRIQUE (-R,Z)
  824.  
  825.  
  826.                            ELSE
  827.  
  828.                               DO 5120 K = 1,NES
  829.                                  U2(K)=-U2(K)
  830.                                  XX1=A2(K,1)
  831.                                  A2(K,1) = A2(K,2) + ECOQ*U2(K)
  832.                                  A2(K,2) = XX1     + ECOQ*U2(K)
  833.  5120                         CONTINUE
  834.  
  835. C... calcul du symetrique
  836.                               AA2(1,1)= -A2(1,2)
  837.                               AA2(2,1)=  A2(2,2)
  838.                               AA2(1,2)= -A2(1,1)
  839.                               AA2(2,2)=  A2(2,1)
  840.                               CALL KNORM2(NES,AA2,NS,UU2,SS2)
  841.  
  842. C... visibilité a priori
  843.                               CALL KPRIOR(NES,NS,NS,A1,A2,U1,U2,IVU)
  844.                               CALL KPRIOR(NES,NS,NS,A1,AA2,U1,UU2,IVUS)
  845.  
  846.                               IF(KIMP.GE.4)WRITE(6,*) ' K1 K2 IVU IVUS '
  847.      $                             ,K1,K2,IVU,IVUS
  848.                               IF(IVU.EQ.0.AND.IVUS.EQ.0) KVU=0
  849.  
  850.                            ENDIF
  851. C   ---------------------------------------------------------------
  852. C<< FIN TEST K1=K2 ET VISIBIITE A PRIORI
  853.                            IF(KIMP.GE.4.AND.KVU.NE.0) THEN
  854.                               WRITE(6,*) ' FACES VISIBLES ',K1,K2
  855.                            ENDIF
  856.  
  857. C>> CALCUL --------------------------------------------------------
  858. C.. option convexe
  859.                            IF(KACHE.EQ.0) THEN
  860.                               IF(KVU.NE.0) THEN
  861.  
  862.                                  CALL KAXC(A1,A2,NPAX,NGAX,FF,KIMP
  863.      $                                ,EXTINC,RAD)
  864.                                  IF(KIMP.GE.4) WRITE(6,*) ' FF CONVEXE '
  865.      $                                ,K1,K2,FF/S1
  866.  
  867.                               ENDIF
  868. C.. option cache
  869.                            ELSE
  870.                               IF(KVU.NE.0) THEN
  871.  
  872. C>> les obstructeurs potentiels sont les memes que pour la face K2
  873.  
  874.                                  IF(KIMP.GE.4) THEN
  875.                                     WRITE(6,*) ' FACES K1 K2 ',K1,K2
  876.      $                                   ,' NOBS ',NOBS
  877.                                  ENDIF
  878.  
  879.                                  CALL KAXK(A1,A2,OBS,NOBS,NTOBS,NPAX
  880.      $                                ,NGAX,FF,KIMP,EXTINC,RAD)
  881.  
  882.                                  IF(KIMP.GE.4) WRITE(6,*) ' FF  ',K1,K2
  883.      $                                ,FF/S1
  884.  
  885.                               ENDIF
  886.  
  887.                            ENDIF
  888. C<< FIN CALCUL ------------------------------------------------------
  889.  
  890. C       WRITE(6,*) ' K1 K2 PLIG ',K1,K2,PLIG
  891.                            PLIG.FACT(K2) =  FF/S1
  892.  
  893.                         ENDIF
  894.  
  895.                      ENDIF
  896.  
  897. C... fin face inverse de K2 (cas des coques) .........................
  898.  
  899.  2100             CONTINUE
  900.  
  901.                   IF (ICOQ.AND.KCOQ(JTYP)) THEN
  902.                      NELT2 = NELT2 + 2 * NEL2
  903.                   ELSE
  904.                      NELT2 = NELT2 + NEL2
  905.                   ENDIF
  906.  
  907.  2000          CONTINUE
  908.  
  909. C     ---FIN BOUCLE FACE 2-----------------------------------------
  910.  
  911.             ENDIF
  912. C*** fin traitement de la face inverse de K1*************************
  913.  
  914.  110     CONTINUE
  915.  
  916.  
  917.          IF (ICOQ.AND.KCOQ(ITYP)) THEN
  918.             NELT1 = NELT1 + 2 * NEL1
  919.          ELSE
  920.             NELT1 = NELT1 + NEL1
  921.          ENDIF
  922.  
  923.  100  CONTINUE
  924.  
  925. C     ---FIN BOUCLE FACE 1-----------------------------------------
  926. C
  927.  
  928.  
  929.  
  930. C     Désactivation des maillages élémentaires
  931.       DO 920 ITYP = 1,NTYP
  932.          IMODEL = MYMOD.KMODEL(ITYP)
  933.          SEGACT IMODEL
  934.          IPT1 = IMAMOD
  935.          SEGDES IPT1
  936.          SEGDES IMODEL
  937.  920  CONTINUE
  938.       SEGDES MYMOD
  939.  
  940. C>> SURFACES DIMENSIONNEES
  941.  
  942.       DO 910 LS=1,NFACE
  943.          PSUR.FACT(LS)=PSUR.FACT(LS)*RAD*RAD
  944.  910  CONTINUE
  945. C
  946.       IF(KACHE.NE.0) SEGSUP SFOBS
  947. C
  948.  
  949. C>>>  NORMALISATION,CACUL DES BILANS ET IMPRESSION
  950. C     --------------------------------------------
  951.  
  952.       IF(EXTINC.GT.1D-3) THEN
  953.          INOR=0
  954.       ENDIF
  955.  
  956.       CALL KFN(IFACFO,INOR,KIMP)
  957.  
  958. C     Traduction puis suppression de IFACFO
  959.  
  960.       IF (KIMP.GE.3) THEN
  961.          CALL PRFACF(IFACFO)
  962.       ENDIF
  963.  
  964.       LTITR=1
  965.       CALL FFMCHA(MYMOD,INFOEL,IFACFO,ICHFAC,LTITR)
  966.  
  967.       SEGACT IFACFO
  968.       DO 930 NNEL = 1,LFACT(/1)
  969.          LFAC = LFACT(NNEL)
  970.          SEGSUP LFAC
  971.  930  CONTINUE
  972.       SEGSUP IFACFO
  973.  
  974.       SEGSUP STRAV
  975.  
  976.       RETURN
  977.       END
  978. C///////////
  979. C mettre apres KPRIOR
  980. *       WRITE (6,*) 'K1 =',K1,' et K2 =',K2
  981. *       WRITE (6,*) 'NES =',NES,'NS =',NS
  982. *       WRITE (6,*) 'A1(*,1) :',(A1(K,1),K=1,IDIM)
  983. *       WRITE (6,*) 'A1(*,2) :',(A1(K,2),K=1,IDIM)
  984. *       WRITE (6,*) 'U1      :',(U1(K),K=1,IDIM+1)
  985. *       WRITE (6,*) 'A2(*,1) :',(A2(K,1),K=1,IDIM)
  986. *       WRITE (6,*) 'A2(*,2) :',(A2(K,2),K=1,IDIM)
  987. *       WRITE (6,*) 'U2      :',(U2(K),K=1,IDIM+1)
  988. *       WRITE (6,*) 'IVU =',IVU
  989. *       WRITE (6,*) 'AA2(*,1) :',(AA2(K,1),K=1,IDIM)
  990. *       WRITE (6,*) 'AA2(*,2) :',(AA2(K,2),K=1,IDIM)
  991. *       WRITE (6,*) 'UU2      :',(UU2(K),K=1,IDIM+1)
  992. *       WRITE (6,*) 'IVUS =',IVUS
  993. C///////////
  994.  
  995.  
  996.  
  997.  
  998.  
  999.  
  1000.  
  1001.  
  1002.  
  1003.  
  1004.  
  1005.  
  1006.  
  1007.  
  1008.  
  1009.  
  1010.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales