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

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