Télécharger la plaquette Cast3M


fr en

Télécharger ffaxca.eso

Retour à la liste

Numérotation des lignes : 

ffaxca

  1. C FFAXCA    SOURCE    CB215821  25/04/22    21:15:06     12245          
  2.       SUBROUTINE FFAXCA (ithr,IDONNEES)
  3. C
  4.       IMPLICIT INTEGER(I-N)
  5.       IMPLICIT REAL*8 (A-H,O-Z)
  6. C
  7.       LOGICAL ICOQ,IFACINF
  8. C----------------------------------------------------------------------
  9. C     Calcul des facteurs de forme en axisymetrique (appele par FACAXI)
  10. C----------------------------------------------------------------------
  11. -INC CCREEL
  12. -INC SMELEME
  13. -INC SMMODEL
  14.  
  15. -INC PPARAM
  16. -INC CCOPTIO
  17. -INC SMCOORD
  18.       POINTEUR MYMOD.MMODEL
  19. C
  20. C     SEGMENT POUR LA PARALLELISATION
  21.       SEGMENT SPARAL
  22.         INTEGER NBTHRD
  23.         INTEGER IAFAIR(NBEL2)
  24.         INTEGER IMYMOD,ISEGEL,KNPAX,KNGAX,KKACHE,IIFACFO
  25.         INTEGER KITYP,KNELT1,IWRKTH
  26.         REAL*8  XEXTINC,XRAD
  27.         LOGICAL BICOQ
  28.       ENDSEGMENT
  29. C
  30. C----------------------------------------------------------------------
  31.       SEGMENT , INFOEL
  32.         LOGICAL KCOQ(N1),KQUAD(N1)
  33.       ENDSEGMENT
  34. C----------------------------------------------------------------------
  35. C     FACTEURS DE FORME
  36. C           NNBEL1 = NOMBRE DE LIGNES + 1
  37. C           NBEL2  = NOMBRE DE COLONNES
  38. C     LFACT(NNBEL1) POINTE SUR LE TABLEAU DES SURFACES
  39. C
  40.       SEGMENT IFACFO
  41.         INTEGER   LFACT(NNBEL1)
  42.       ENDSEGMENT
  43.       SEGMENT LFAC
  44.         REAL*8    FACT(NBEL2)
  45.       ENDSEGMENT
  46.       POINTEUR PSUR.LFAC, PLIG.LFAC, PCOL.LFAC
  47. C----------------------------------------------------------------------
  48. C     coordonnees des obstructeurs
  49.       SEGMENT SFOBS
  50.         REAL*8 OBS(2,NTOBS)
  51.       ENDSEGMENT
  52. C----------------------------------------------------------------------
  53.       SEGMENT STRAV
  54.         REAL*8 A1(NA,NA),A2(NA,NA),A3(NA,NA),AA2(NA,NA)
  55.         REAL*8 U1(NA1),U2(NA1),UU2(NA1)
  56.       ENDSEGMENT
  57. C----------------------------------------------------------------------
  58.       SEGMENT SEGTH
  59.         INTEGER  SSFOBS(NTHRD)
  60.         INTEGER  SSTRAV(NTHRD)
  61.       ENDSEGMENT
  62. C
  63. C  reprise de la valeur de facaxi
  64.       ECOQ=1.D-3
  65.       NS = 2 
  66.       EPS = 1D-5
  67.       KIMP = IIMPI
  68.       NES = IDIM
  69. C
  70.       SPARAL = IDONNEES
  71.       NBTHR  = SPARAL.NBTHRD
  72.       MYMOD  = SPARAL.IMYMOD
  73.       INFOEL = SPARAL.ISEGEL
  74.       NPAX   = SPARAL.KNPAX
  75.       NGAX   = SPARAL.KNGAX
  76.       KACHE  = SPARAL.KKACHE
  77.       EXTINC = SPARAL.XEXTINC
  78.       IFACFO = SPARAL.IIFACFO
  79.       ITYP   = SPARAL.KITYP
  80.       RAD    = SPARAL.XRAD
  81.       ICOQ   = SPARAL.BICOQ
  82.       NELT1  = SPARAL.KNELT1
  83.       SEGTH  = SPARAL.IWRKTH
  84.       SFOBS  = SEGTH.SSFOBS(ithr)
  85.       STRAV  = SEGTH.SSTRAV(ithr)
  86. CC    WRITE(*,*) 'ith SFOBS STRAV',ithr,sfobs,strav
  87. C
  88.       NNBEL1 = LFACT(/1)
  89.       PSUR = LFACT(NNBEL1)
  90.       NTYP = MYMOD.KMODEL(/1)
  91. C
  92.       IMODEL = MYMOD.KMODEL(ITYP)
  93.       IPT1 = IMAMOD
  94.       NSGEO1 = IPT1.NUM(/1)
  95.       NSPA1=1
  96.       IF(ICOQ.AND.KQUAD(ITYP)) NSPA1=2
  97.       NEL1 = IPT1.NUM(/2)
  98. C
  99. C     Decoupage pour le travail d'ecriture en parallele
  100.       IRES=MOD(NEL1,NBTHR)
  101.       IF (IRES.EQ.0) THEN
  102.         ILON = NEL1 / NBTHR
  103.         IDEB = (ithr -1)* ILON + 1
  104.       ELSE
  105.         IF (ithr.LE.IRES) THEN
  106.           ILON = (NEL1 / NBTHR) + 1
  107.           IDEB = (ithr -1)* ILON + 1
  108.         ELSE
  109.           ILON = NEL1 / NBTHR
  110.           IDEB = (IRES * (ILON+1)) + (ithr-IRES-1)* ILON + 1
  111.         ENDIF
  112.       ENDIF
  113.       IFIN = IDEB + ILON - 1
  114. C
  115.       DO 110 I1 = IDEB,IFIN
  116.  
  117.         IF (ICOQ.AND.KCOQ(ITYP)) THEN
  118.           K1 = (2*I1-1) + NELT1
  119.         ELSE
  120.           K1 = I1+ NELT1
  121.         ENDIF
  122.  
  123.         PLIG = LFACT(K1)
  124.  
  125. C*** traitement de la face K1 ***************************************
  126.  
  127.         DO 111 IS = 1,NSGEO1,NSPA1
  128.           LS=IS
  129.           IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2
  130.           IREF = (IDIM+1)*(IPT1.NUM(IS,I1)-1)
  131.           IF(KIMP.GE.4) THEN
  132.             WRITE (6,*) 'Noeud numero',IPT1.NUM(IS,I1),'ref :',IREF
  133.           ENDIF
  134.           DO 112 K = 1,NES
  135.             A1(K,LS) = XCOOR(IREF+K)/RAD
  136.  112      CONTINUE
  137.  111    CONTINUE
  138.  
  139.         CALL KNORM2(NES,A1,NS,U1,S1)
  140.         S1=XPI*S1*ABS(A1(1,1)+A1(1,2))
  141.         PSUR.FACT(K1)=S1
  142.         ZMIN1= MIN(A1(2,1),A1(2,2))
  143.         ZMAX1= MAX(A1(2,1),A1(2,2))
  144.         RMAX1= MAX(A1(1,1),A1(1,2))
  145. C
  146. C>>     BOUCLE FACE 2
  147. C       -----------------------------------------------------------
  148. C
  149.         IFACINF = .TRUE.
  150.   10    CONTINUE
  151. C
  152.         NELT2= 0
  153.         DO 200 JTYP = 1,NTYP
  154. C
  155.           IMODEL = MYMOD.KMODEL(JTYP)
  156.           IPT2 = IMAMOD
  157.           NSGEO2 = IPT2.NUM(/1)
  158.           NSPA2=1
  159.           IF(ICOQ.AND.KQUAD(JTYP)) NSPA2=2
  160.           NEL2 = IPT2.NUM(/2)
  161. C
  162.           DO 210 I2 = 1,NEL2
  163.             IF (ICOQ.AND.KCOQ(JTYP))THEN
  164.               K2 = 2*I2-1  + NELT2
  165.             ELSE
  166.               K2 = I2 + NELT2
  167.             ENDIF
  168.  
  169. C... face K2 .....................................................
  170.  
  171.             KVU=1
  172.             FF=0.D0
  173.  
  174. C.. UTILISATION DE LA RECIPROCITE
  175. C
  176.             IF(K1.GT.K2) THEN
  177.               IF (NBTHR.NE.1) THEN
  178.                 IAFAIR(K1) = K2
  179.               ELSE
  180.                 S2=PSUR.FACT(K2)
  181.                 PCOL=LFACT(K2)
  182.                 PLIG.FACT(K2)=(S2/S1)*PCOL.FACT(K1)
  183.               ENDIF
  184.             ELSE
  185. C.. TEST K1=K2 ET VISIBILITE A PRIORI
  186.  
  187. C------------------------------------------------------------------
  188.  
  189. C>> K1=K2 : ELIMINATION DES CAS TRIVIAUX OU F(K1,K2)=0
  190. C
  191.               IF(K1.EQ.K2) THEN
  192.                 IF(KIMP.GE.4) THEN
  193.                   WRITE(6,*) ' A1 ',A1(1,1),A1(2,1)
  194.                   WRITE(6,*) ' A1 ',A1(1,2),A1(2,2)
  195.                   WRITE(6,*) ' U1 ',U1(1),U1(2)
  196.                 ENDIF
  197.                 IF (ABS(U1(1)).LT.EPS) KVU=0
  198.                 IF (U1(1).GE.EPS) KVU=0
  199.  
  200.                 IF(KVU.NE.0) THEN
  201.                   DO 214 K = 1,NES
  202.                     A2(K,1) = A1(K,1)
  203.                     A2(K,2) = A1(K,2)
  204.  214              CONTINUE
  205.                 ENDIF
  206.  
  207. C>> K1::K2 : ELIMINATION DES CAS DE NON-VISIBLITE
  208. C            UTILISATION DE LA VISIBILITE A PRIORI DANS LE PLAN R-Z
  209. C            ON TESTE LA FACE K2 ET SON SYMETRIQUE (-R,Z)
  210.  
  211.               ELSE
  212.  
  213.                 DO 211 IS = 1,NSGEO2,NSPA2
  214.                   LS=IS
  215.                   IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2
  216.                   IREF = (IDIM+1)*(IPT2.NUM(IS,I2)-1)
  217.                   DO 212 K = 1,NES
  218.                     A2(K,LS) = XCOOR(IREF+K)/RAD
  219.  212              CONTINUE
  220.  211            CONTINUE
  221.                 CALL KNORM2(NES,A2,NS,U2,S2)
  222.                 S2=XPI*S2*ABS(A2(1,1)+A2(1,2))
  223.  
  224. C.. calcul du symetrique
  225.                 AA2(1,1)= -A2(1,2)
  226.                 AA2(2,1)=  A2(2,2)
  227.                 AA2(1,2)= -A2(1,1)
  228.                 AA2(2,2)=  A2(2,1)
  229.                 CALL KNORM2(NES,AA2,NS,UU2,SS2)
  230.  
  231. C.. visibilite a priori
  232.                 CALL KPRIOR(NES,NS,NS,A1,A2,U1,U2,IVU)
  233.                 CALL KPRIOR(NES,NS,NS,A1,AA2,U1,UU2,IVUS)
  234.  
  235.                 IF(KIMP.GE.4)WRITE(6,*)' K1 K2 IVU IVUS ',K1,K2,IVU,IVUS
  236.                 IF(IVU.EQ.0.AND.IVUS.EQ.0) KVU=0
  237.  
  238.               ENDIF
  239. C   ---------------------------------------------------------------
  240. C<< FIN TEST K1=K2 ET VISIBIITE A PRIORI
  241.               IF(KIMP.GE.4.AND.KVU.NE.0) THEN
  242.                 WRITE(6,*) ' FACES VISIBLES ',K1,K2
  243.               ENDIF
  244.  
  245. C>> CALCUL
  246. C... option convexe
  247.               IF(KACHE.EQ.0) THEN
  248.                 IF(KVU.NE.0) THEN
  249.                   CALL KAXC(A1,A2,NPAX,NGAX,FF,KIMP,EXTINC,RAD)
  250.                   IF(KIMP.GE.4) WRITE(6,*) ' FF CONVEXE ',K1,K2,FF/S1
  251.  
  252.                 ENDIF
  253.               ELSE
  254. C.. option cache
  255.                 IF(KVU.NE.0) THEN
  256.  
  257. C>> RECHERCHE DES OBSTRUCTEURS POTENTIELS.--------------------------
  258.  
  259.                   NOBS=0
  260.  
  261.                   ZMIN2= MIN(A2(2,1),A2(2,2))
  262.                   ZMAX2= MAX(A2(2,1),A2(2,2))
  263.                   RMAX2= MAX(A2(1,1),A2(1,2))
  264.                   ZTMIN= MIN(ZMIN1,ZMIN2)
  265.                   ZTMAX= MAX(ZMAX1,ZMAX2)
  266.                   RMAX = MAX(RMAX1,RMAX2)
  267.  
  268. C>> BOUCLE FACE 3---------------------------------------------------
  269.                   NELT3= 0
  270.                   DO 300 KTYP = 1,NTYP
  271.                     IMODEL = MYMOD.KMODEL(KTYP)
  272.                     IPT3 = IMAMOD
  273.                     NSGEO3 = IPT3.NUM(/1)
  274.                     NSPA3=1
  275.                     IF(ICOQ.AND.KQUAD(KTYP)) NSPA3=2
  276.                     NEL3 = IPT3.NUM(/2)
  277.                     DO 310 I3 = 1,NEL3
  278.                       K3 = I3+ NELT3
  279.  
  280.                       IF(K3.NE.K1.AND.K3.NE.K2) THEN
  281.  
  282.                         DO 311 IS = 1,NSGEO3,NSPA3
  283.                           LS=IS
  284.                           IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2
  285.                           IREF = (IDIM+1)*(IPT3.NUM(IS,I3)-1)
  286.                           DO 312 K = 1,NES
  287.                             A3(K,LS) = XCOOR(IREF+K)/RAD
  288.  312                      CONTINUE
  289.  311                    CONTINUE
  290.  
  291.                         IF( MAX(A3(2,1),A3(2,2)).LE.ZTMIN) THEN
  292.                           KOBS=0
  293.                         ELSEIF(MIN(A3(2,1),A3(2,2)).GE.ZTMAX) THEN
  294.                           KOBS=0
  295.                         ELSEIF(MIN(A3(1,1),A3(1,2)).GE.RMAX) THEN
  296.                           KOBS=0
  297.                         ELSE
  298.                           KOBS=1
  299.                           OBS(1,2*NOBS+1)=A3(1,1)
  300.                           OBS(2,2*NOBS+1)=A3(2,1)
  301.                           OBS(1,2*NOBS+2)=A3(1,2)
  302.                           OBS(2,2*NOBS+2)=A3(2,2)
  303.                           NOBS=NOBS+1
  304.                         ENDIF
  305.                       ENDIF
  306.  310                CONTINUE
  307. C
  308.                     NELT3=NELT3+NEL3
  309. C
  310.  300              CONTINUE
  311.  
  312. C<< FIN RECHERCHE OBSTRUCTEURS--------------------------------------
  313.  
  314.                   IF(KIMP.GE.4) THEN
  315.                     WRITE(6,*) ' FACES K1 K2 ',K1,K2,' NOBS ',NOBS
  316.                   ENDIF
  317.  
  318.                   CALL KAXK(A1,A2,OBS,NOBS,NTOBS,NPAX,NGAX,FF
  319.      $                            ,KIMP,EXTINC,RAD)
  320.                   IF(KIMP.GE.4) WRITE(6,*) ' FF  ',K1,K2,FF/S1
  321.  
  322.                 ENDIF
  323.               ENDIF
  324. C<< FIN CALCUL
  325. C   ---------------------------------------------------------------
  326.  
  327. C             WRITE(6,*) ' K1=',K1,' K2 PLIG ',K2,PLIG
  328.               PLIG.FACT(K2) =  FF/S1
  329.  
  330.             ENDIF
  331. C... fin face K2 ..................................................
  332.  
  333. C... face inverse de K2 (cas des coques) ..........................
  334.  
  335.             IF (ICOQ.AND.KCOQ(JTYP)) THEN
  336.               K2=K2 + 1
  337.               KVU=1
  338.               FF=0.D0
  339.  
  340. C UTILISATION DE LA RECIPROCITE
  341. C
  342.               IF(K1.GT.K2) THEN
  343.                 S2=PSUR.FACT(K2)
  344.                 PCOL=LFACT(K2)
  345.                 PLIG.FACT(K2)=(S2/S1)*PCOL.FACT(K1)
  346.               ELSE
  347.  
  348. C>> TEST K1=K2 ET VISIBILITE A PRIORI
  349.  
  350. C------------------------------------------------------------------
  351.  
  352. C>> K1=K2 : ELIMINATION DES CAS TRIVIAUX OU F(K1,K2)=0
  353. C
  354.                 IF(K1.EQ.K2) THEN
  355.                   IF(KIMP.GE.4) THEN
  356.                     WRITE(6,*) ' A1 ',A1(1,1),A1(2,1)
  357.                     WRITE(6,*) ' A1 ',A1(1,2),A1(2,2)
  358.                     WRITE(6,*) ' U1 ',U1(1),U1(2)
  359.                   ENDIF
  360.                   IF (ABS(U1(1)).LT.EPS) KVU=0
  361.                   IF (U1(1).GE.EPS) KVU=0
  362.  
  363.                   IF(KVU.NE.0) THEN
  364.                     DO 514 K = 1,NES
  365.                       A2(K,1) = A1(K,1)
  366.                       A2(K,2) = A1(K,2)
  367.  514                CONTINUE
  368.                   ENDIF
  369.  
  370. C>> K1::K2 : ELIMINATION DES CAS DE NON-VISIBLITE
  371. C            UTILISATION DE LA VISIBILITE A PRIORI DANS LE PLAN R-Z
  372. C            ON TESTE LA FACE K2 ET SON SYMETRIQUE (-R,Z)
  373.  
  374.  
  375.                 ELSE
  376.  
  377.                   DO 512 K = 1,NES
  378.                     U2(K)=-U2(K)
  379.                     XX1= A2(K,1)
  380.                     A2(K,1) = A2(K,2) + ECOQ*U2(K)
  381.                     A2(K,2) = XX1     + ECOQ*U2(K)
  382.  512              CONTINUE
  383.  
  384. C... calcul du symetrique
  385.                   AA2(1,1)= -A2(1,2)
  386.                   AA2(2,1)=  A2(2,2)
  387.                   AA2(1,2)= -A2(1,1)
  388.                   AA2(2,2)=  A2(2,1)
  389.                   CALL KNORM2(NES,AA2,NS,UU2,SS2)
  390.  
  391. C... visibilite a priori
  392.                   CALL KPRIOR(NES,NS,NS,A1,A2,U1,U2,IVU)
  393.                   CALL KPRIOR(NES,NS,NS,A1,AA2,U1,UU2,IVUS)
  394.  
  395.                   IF(KIMP.GE.4)WRITE(6,*) ' K1 K2 IVU IVUS ',K1,K2,IVU,
  396.      &                                      IVUS
  397.                   IF(IVU.EQ.0.AND.IVUS.EQ.0) KVU=0
  398.  
  399.                 ENDIF
  400. C   ---------------------------------------------------------------
  401. C<< FIN TEST K1=K2 ET VISIBIITE A PRIORI
  402.                 IF(KIMP.GE.4.AND.KVU.NE.0) THEN
  403.                   WRITE(6,*) ' FACES VISIBLES ',K1,K2
  404.                 ENDIF
  405.  
  406. C>> CALCUL --------------------------------------------------------
  407. C.. option convexe
  408.                 IF(KACHE.EQ.0) THEN
  409.                   IF(KVU.NE.0) THEN
  410.                     CALL KAXC(A1,A2,NPAX,NGAX,FF,KIMP,EXTINC,RAD)
  411.                     IF(KIMP.GE.4) WRITE(6,*) ' FF CONVEXE ',K1,K2,FF/S1
  412.                   ENDIF
  413. C.. option cache
  414.                 ELSE
  415.                   IF(KVU.NE.0) THEN
  416.  
  417. C>> les obstructeurs potentiels sont les memes que pour la face K2
  418.  
  419.                     IF(KIMP.GE.4) THEN
  420.                       WRITE(6,*) ' FACES K1 K2 ',K1,K2,' NOBS ',NOBS
  421.                     ENDIF
  422.  
  423.                     CALL KAXK(A1,A2,OBS,NOBS,NTOBS,NPAX,NGAX
  424.      $                               ,FF,KIMP,EXTINC,RAD)
  425.                     IF(KIMP.GE.4) WRITE(6,*) ' FF  ',K1,K2,FF/S1
  426.                   ENDIF
  427.                 ENDIF
  428. C<< FIN CALCUL ------------------------------------------------------
  429.  
  430. C               WRITE(6,*) ' K1 K2 PLIG ',K1,K2,PLIG
  431.                 PLIG.FACT(K2) =  FF/S1
  432.  
  433.               ENDIF
  434.             ENDIF
  435.  
  436. C... fin face inverse de K2 (cas des coques) .........................
  437.  
  438.  210      CONTINUE
  439.  
  440.           IF (ICOQ.AND.KCOQ(JTYP)) THEN
  441.             NELT2 = NELT2 + 2 * NEL2
  442.           ELSE
  443.             NELT2 = NELT2 + NEL2
  444.           ENDIF
  445.  
  446.  200    CONTINUE
  447.  
  448. C*** traitement de la face inverse de K1*****************************
  449.         IF (ICOQ.AND.KCOQ(ITYP).AND.IFACINF) THEN
  450.           K1=K1+1
  451.           PLIG = LFACT(K1)
  452.  
  453.           DO 122 K = 1,NES
  454.             U1(K) =-U1(K)
  455.             XX1=A1(K,1)
  456.             A1(K,1) = A1(K,2) + ECOQ * U1(K)
  457.             A1(K,2) = XX1     + ECOQ * U1(K)
  458.  122      CONTINUE
  459.  
  460.           PSUR.FACT(K1)=S1
  461.           IFACINF = .FALSE.
  462.           GOTO 10
  463.  
  464.         ENDIF
  465. C*** fin traitement de la face inverse de K1*************************
  466.  
  467.  110  CONTINUE
  468. C
  469.       RETURN
  470.       END
  471. C///////////
  472. C mettre apres KPRIOR
  473. *       WRITE (6,*) 'K1 =',K1,' et K2 =',K2
  474. *       WRITE (6,*) 'NES =',NES,'NS =',NS
  475. *       WRITE (6,*) 'A1(*,1) :',(A1(K,1),K=1,IDIM)
  476. *       WRITE (6,*) 'A1(*,2) :',(A1(K,2),K=1,IDIM)
  477. *       WRITE (6,*) 'U1      :',(U1(K),K=1,IDIM+1)
  478. *       WRITE (6,*) 'A2(*,1) :',(A2(K,1),K=1,IDIM)
  479. *       WRITE (6,*) 'A2(*,2) :',(A2(K,2),K=1,IDIM)
  480. *       WRITE (6,*) 'U2      :',(U2(K),K=1,IDIM+1)
  481. *       WRITE (6,*) 'IVU =',IVU
  482. *       WRITE (6,*) 'AA2(*,1) :',(AA2(K,1),K=1,IDIM)
  483. *       WRITE (6,*) 'AA2(*,2) :',(AA2(K,2),K=1,IDIM)
  484. *       WRITE (6,*) 'UU2      :',(UU2(K),K=1,IDIM+1)
  485. *       WRITE (6,*) 'IVUS =',IVUS
  486. C///////////
  487.  
  488.  
  489.  
  490.  

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