Télécharger la plaquette Cast3M


fr en

Télécharger zclpls.eso

Retour à la liste

Numérotation des lignes : 

zclpls

  1. C ZCLPLS    SOURCE    CHAT      05/01/13    04:21:27     5004
  2.       SUBROUTINE ZCLPLS(HR,RPG,LE,NEL,K0,NPTD,NX,IES,NP,IAXI,IPADL,
  3.      & COEFF,IK1,
  4.      & TN,G,
  5.      & VOLU,COTE,NELZ,DME,DIAM,
  6.      & DT,DTT2,NUEL,DIAEL)
  7. C
  8. C VERSION VECTORISEE (CF XCVTIT POUR PLUS DE DETAILS)
  9. C
  10.       IMPLICIT INTEGER(I-N)
  11.       IMPLICIT REAL*8 (A-H,O-Z)
  12. C NOMBRE MAXI DE POINTS PAR ELEMENT : NPX
  13.       PARAMETER(NPX=9)
  14. C LONGUEUR DES REGISTRES VECTORIELS DE LA MACHINE CIBLE
  15.       PARAMETER(LRV=64)
  16. C***********************************************************************
  17. C
  18. C  CE SP DISCRETISE LE LAPLACIEN D UNE VARIABLE SCALAIRE
  19. C
  20. C    EN 1D SUR L ELEMENT SEG2
  21. C    EN 2D SUR L ELEMENT QUA4       PLAN OU AXI
  22. C    EN 3D SUR L ELEMENT CUB8 PRI6 et TET4
  23. C  LES OPERATEURS SONT SOUS INTEGRES
  24. C
  25. C  COEFF(SCAL DOMA)  LE COEFFICIENT DU LAPLACIEN  ( NU )
  26. C       (SCAL ELEM)
  27. C
  28. C     NUELD(NELZ): NUMERO DE L'ELEMENT NK DE LA ZONE DANS LE DOMAINE
  29. C     (NELZ      : NOMBRE D'ELEMENTS DE LA ZONE)
  30. C     ROC(NELD)  : ROC PAR ELEMENT DONNE SUR TOUT LE DOMAINE
  31. C
  32. C
  33. C***********************************************************************
  34. -INC CCVQUA4
  35. -INC CCREEL
  36. C
  37.       DIMENSION TN(NPTD,NX),GG(8)
  38.       DIMENSION COEFF(*)
  39.       DIMENSION COTE(NELZ,IES),DME(*),VOLU(*),DIAM(*)
  40.  
  41.       DIMENSION IPADL(1),LE(NP,1)
  42.       DIMENSION HR(NEL,NP,IES),RPG(1)
  43.  
  44.       REAL*8    G(NPTD,NX)
  45.       REAL*8    ROC(1)
  46.       DIMENSION NUELD(1)
  47.       DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8)
  48. C
  49. C
  50.  
  51. C* -TABLEAUX ADDITIONNELS POUR LA VECTORISATION   **********************
  52.       DIMENSION AIRE(LRV),ALF (LRV),COEF(LRV),CLSR(LRV),GINC(LRV)
  53.       DIMENSION AL  (LRV),AH  (LRV),AP  (LRV),CFM (LRV),XMA (LRV)
  54.       DIMENSION XMB(LRV),XMD(LRV),XMH(LRV),XMI(LRV),AHL(LRV)
  55.       DIMENSION ALH(LRV),DPR(LRV),TETA(LRV,NPX,3),F(LRV,NPX),BF(LRV,NPX)
  56. C
  57. C --- TABLEAU POUR L'OPTION RAPIDE ---
  58.       DIMENSION ALP(LRV)
  59. C***
  60.       SAVE IPAS,QGGT,Q1,Q2,Q3
  61.       DATA IPAS/0/
  62.  
  63. C
  64. C           INITIALISATIONS DIVERSES
  65. C
  66. C      WRITE(6,*)' ROC=',(ROC(MM),MM=1,NEL)
  67. C      WRITE(6,*)' NUELD=',(NUELD(MM),MM=1,NEL)
  68.  
  69.       NK=K0
  70.       CALL INITD(BF,64,0.D0)
  71.  
  72.       IF(IES.EQ.2)GO TO 20
  73.       IF(IES.EQ.3)GO TO 30
  74. C**********
  75. C                                          *   1D*
  76. C**********
  77. C
  78. C
  79. C K EST LE NUMERO DE L'ELEMENT
  80. C KP PERMET DE SE SITUER A L'INTERIEUR DES TABLEAUX DE LRV ELEMENTS
  81.       NPACK=INT(NEL/FLOAT(LRV))+1
  82.       KPACKD=1
  83.       KPACKF=NPACK
  84. C
  85. C*** BOUCLE SUR LES ELEMENTS ***
  86. C
  87.       DO 70001 KPACK=KPACKD,KPACKF
  88. C     DO 40 K=1,NEL
  89. C
  90. C POUR CHAQUE PAQUET DE LRV ELEMENTS
  91.       KDEB=1+(KPACK-1)*LRV
  92.       KFIN=MIN(NEL,KDEB+LRV-1)
  93.       DO 70002 K=KDEB,KFIN
  94.       KP=K-KDEB+1
  95.       NK=K+K0
  96.       K1=1+(1-IK1)*(NK-1)
  97.       COEF(KP)=COEFF(K1)
  98.       CLSR(KP)=COEFF(K1)
  99.       AIRE(KP)=VOLU(NK)
  100.       XMA(KP)=AIRE(KP)
  101. 70002 CONTINUE
  102.       DO 4 I=1,NP
  103.       DO 4 N=1,NX
  104.       DO 4 K=KDEB,KFIN
  105.       KP=K-KDEB+1
  106.       NF=IPADL(LE(I,K))
  107.       TETA(KP,I,N)=TN(NF,N)
  108.     4 CONTINUE
  109. C
  110. C
  111. C
  112.       DO 70003 K=KDEB,KFIN
  113.       NK=K+K0
  114.       KP=K-KDEB+1
  115.       DT0=DT
  116.       DT2=0.5*XMA(KP)*XMA(KP)/CLSR(KP)
  117.       IF(DT2.LT.DT)DT=DT2
  118.       IF(DT.EQ.DT0)GO TO 152
  119.       DTT2=DT2
  120.       DIAEL=XMA(KP)
  121.       NUEL=NK
  122.  152  CONTINUE
  123. 70003 CONTINUE
  124. C
  125. C
  126.       N=1
  127. C --- CALCUL DE L'INCREMENT GINC ---------------------------
  128.       DO 70104 K=KDEB,KFIN
  129.       KP=K-KDEB+1
  130.       GINC(KP)=COEF(KP)*(TETA(KP,1,N)-TETA(KP,2,N))/AIRE(KP)
  131. 70104 CONTINUE
  132. C --- ACCUMULATION DANS LE TABLEAU G -----------------------
  133.       DO 70004 K=KDEB,KFIN
  134.       KP=K-KDEB+1
  135.       NF=IPADL(LE(1,K))
  136.       G(NF,N)=G(NF,N) +  GINC(KP)
  137.       NF=IPADL(LE(2,K))
  138.       G(NF,N)=G(NF,N) -  GINC(KP)
  139. 70004 CONTINUE
  140. 70001 CONTINUE
  141.   40  CONTINUE
  142.       IPAS=1
  143.       RETURN
  144. C                           ***********
  145. C                           *   2D    *
  146. C                           ***********
  147.   20  CONTINUE
  148. C
  149.  
  150. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  151.  
  152. C                           DIFFERENCES TRIANGLE / QUADRANGLE
  153.       QUA4=0.D0
  154.       IF(NP.EQ.4)QUA4=1.D0
  155.  
  156. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  157.  
  158.  
  159. C
  160. C K EST LE NUMERO DE L'ELEMENT
  161. C KP PERMET DE SE SITUER A L'INTERIEUR DES TABLEAUX DE LRV ELEMENTS
  162.       NPACK=INT(NEL/FLOAT(LRV))+1
  163.       KPACKD=1
  164.       KPACKF=NPACK
  165. C
  166. C*** BOUCLE SUR LES ELEMENTS ***
  167. C
  168.       DO 80001 KPACK=KPACKD,KPACKF
  169. C
  170. C POUR CHAQUE PAQUET DE LRV ELEMENTS
  171.       KDEB=1+(KPACK-1)*LRV
  172.       KFIN=MIN(NEL,KDEB+LRV-1)
  173.       DO 80002 K=KDEB,KFIN
  174.       KP=K-KDEB+1
  175.       NK=K+K0
  176.       K1=1+(1-IK1)*(NK-1)
  177.       XMI(KP)=DIAM(NK)
  178.       XMA(KP)=DME(NK)
  179.       COEF(KP)=COEFF(K1)
  180.       CLSR(KP)=COEFF(K1)
  181.       AIRE(KP)=VOLU(NK)
  182. CC
  183.       AL(KP)=COTE(NK,1)
  184.       AH(KP)=COTE(NK,2)
  185.       XMB(KP)=(AL(KP)+AH(KP))/2.
  186.       XMD(KP)=((AL(KP)*AH(KP))**2)/(AL(KP)**2+AH(KP)**2)
  187.       AHL(KP)=AH(KP)/AL(KP)
  188.       ALH(KP)=AL(KP)/AH(KP)
  189.       DPR(KP)=1.D0
  190.       IF(IAXI.NE.0)DPR(KP)=2.D0*XPI*RPG(K)
  191. 80002 CONTINUE
  192. CC
  193. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  194. C
  195.       DO 5 I=1,NP
  196.       DO 5 N=1,NX
  197.       DO 5 K=KDEB,KFIN
  198.       KP=K-KDEB+1
  199.       NF=IPADL(LE(I,K))
  200.       TETA(KP,I,N)=TN(NF,N)
  201.     5 CONTINUE
  202. C
  203. C
  204. C
  205.       DO 80004 K=KDEB,KFIN
  206.       KP=K-KDEB+1
  207.       NK=K+K0
  208.       DT0=DT
  209.       DT2=0.5*XMD(KP)/CLSR(KP)
  210.       IF(DT2.LT.DT)DT=DT2
  211.       IF(DT.EQ.DT0)GO TO 252
  212.       DTT2=DT2
  213.       DIAEL=XMB(KP)
  214.       NUEL=NK
  215.  252  CONTINUE
  216. 80004 CONTINUE
  217. C
  218. C --- VERSION RAPIDE -----------------------
  219.       CF12=1.D0/12.D0
  220.       DO 80011 K=KDEB,KFIN
  221.       KP=K-KDEB+1
  222.       ALP(KP)=CF12*(AHL(KP)+ALH(KP))*DPR(KP)*QUA4
  223. 80011 CONTINUE
  224. C ------------------------------------------
  225. C
  226.       DO 1 N=1,NX
  227.       DO 80010 I= 1,NP
  228.       DO 80010 K=KDEB,KFIN
  229.       KP=K-KDEB+1
  230.       BF(KP,I)=0.D0
  231. 80010 CONTINUE
  232.  
  233.       DO 1 I=1,NP
  234.       DO 3 J= 1,NP
  235.       DO 3 K=KDEB,KFIN
  236.       KP=K-KDEB+1
  237.       BF(KP,I)=BF(KP,I)+TETA(KP,J,N)*
  238.      &       ((HR(K,I,1)*HR(K,J,1)+HR(K,I,2)*HR(K,J,2))*AIRE(KP)
  239.      &        +ALP(KP)*VGGT(J,I))*COEF(KP)
  240.    3  CONTINUE
  241.  
  242. C --- ACCUMULATION DANS LE TABLEAU G ------------------------------
  243.       DO 80006 K=KDEB,KFIN
  244.       KP=K-KDEB+1
  245.       NF=IPADL(LE(I,K))
  246.       G(NF,N)=G(NF,N)+BF(KP,I)
  247. 80006 CONTINUE
  248.     1 CONTINUE
  249. 80001 CONTINUE
  250.   50  CONTINUE
  251.  
  252. C     WRITE(6,*)' SUB XCLPLS G(1,='
  253. C     WRITE(6,1002) (G(MM,1),MM=1,NPTD)
  254. C     WRITE(6,*)' SUB XCLPLS G(2,='
  255. C     WRITE(6,1002) (G(MM,2),MM=1,NPTD)
  256. C     WRITE(6,*)' FIN **** '
  257.  
  258.       IPAS=1
  259.       RETURN
  260.  30   CONTINUE
  261. C                           ***********
  262. C                           *   3D    *
  263. C                           ***********
  264.  
  265.       IF(IPAS.EQ.0)CALL CALHRH(QGGT,Q1,Q2,Q3,IES)
  266.       CUB8=0.D0
  267.       IF(NP.EQ.8)CUB8=1.D0
  268.  
  269.  1003 FORMAT('   XCVTI ',10I10)
  270. C
  271. C K EST LE NUMERO DE L'ELEMENT
  272. C KP PERMET DE SE SITUER A L'INTERIEUR DES TABLEAUX DE LRV ELEMENTS
  273.       NPACK=INT(NEL/FLOAT(LRV))+1
  274.       KPACKD=1
  275.       KPACKF=NPACK
  276. C
  277. C*** BOUCLE SUR LES ELEMENTS ***
  278. C
  279.       DO 90001 KPACK=KPACKD,KPACKF
  280. C     DO 60 K=1,NEL
  281. C
  282. C POUR CHAQUE PAQUET DE LRV ELEMENTS
  283.       KDEB=1+(KPACK-1)*LRV
  284.       KFIN=MIN(NEL,KDEB+LRV-1)
  285.       DO 90002 K=KDEB,KFIN
  286.       KP=K-KDEB+1
  287.       NK=K+K0
  288.       K1=1+(1-IK1)*(NK-1)
  289.       AL(KP)=COTE(NK,1)
  290.       AH(KP)=COTE(NK,2)
  291.       AP(KP)=COTE(NK,3)
  292.       CFM(KP)=AL(KP)*AH(KP)/AP(KP)+AL(KP)*AP(KP)/AH(KP)+
  293.      &        AP(KP)*AH(KP)/AL(KP)
  294.       XMI(KP)=DIAM(NK)
  295.       XMA(KP)=DME(NK)
  296.       COEF(KP)=COEFF(K1)
  297. C      CLSR(KP)=COEFF(K1)/ROC(NUELD(NK))
  298.       CLSR(KP)=COEFF(K1)
  299.       AIRE(KP)=VOLU(NK)
  300.       XMB(KP)=(XMA(KP)+XMI(KP))*0.5
  301.       XMD(KP)=((XMA(KP)*XMI(KP))**2)/(2.*XMA(KP)**2+XMI(KP)**2)
  302.       XMH(KP)=AIRE(KP)**0.33
  303. 90002 CONTINUE
  304. C
  305.       DO 15 I=1,NP
  306.       DO 15 N=1,NX
  307.       DO 15 K=KDEB,KFIN
  308.       KP=K-KDEB+1
  309.       NF=IPADL(LE(I,K))
  310.       TETA(KP,I,N)=TN(NF,N)
  311.    15 CONTINUE
  312. C     WRITE(6,*)' K ** ',K,' TETA ','    VOLU=',AIRE,'  COEF=',COEF
  313. C    &,'  XMH=',XMH
  314. C     WRITE(6,1002)(TETA(MM,1),MM=1,8)
  315. C
  316. C
  317.       DO 90003 K=KDEB,KFIN
  318.       KP=K-KDEB+1
  319.       NK=K+K0
  320.       DT0=DT
  321.       DT2=0.5*XMD(KP)/CLSR(KP)
  322.       IF(DT2.LT.DT)DT=DT2
  323.       IF(DT.EQ.DT0)GO TO 353
  324.       DTT2=DT2
  325.       DIAEL=XMB(KP)
  326.       NUEL=NK
  327.  353  CONTINUE
  328. 90003 CONTINUE
  329. C
  330.       DO 11 N=1,NX
  331.       DO 90004 I=1,NP
  332.       DO 90004 K=KDEB,KFIN
  333.       KP=K-KDEB+1
  334.       BF(KP,I)=0.
  335. 90004 CONTINUE
  336.       DO 11 I=1,NP
  337.       DO 13 J=1,NP
  338.       DO 13 K=KDEB,KFIN
  339.       KP=K-KDEB+1
  340.       GEO1=CFM(KP)*QGGT(J,I)*CUB8
  341.       BF(KP,I)=BF(KP,I)+TETA(KP,J,N)*
  342.      &(HR(K,I,1)*HR(K,J,1)+HR(K,I,2)*HR(K,J,2)+HR(K,I,3)*HR(K,J,3))
  343.      & *AIRE(KP)+COEF(KP)*XMH(KP)*GEO1
  344.   13  CONTINUE
  345.       DO 90006 K=KDEB,KFIN
  346.       KP=K-KDEB+1
  347.       NF=IPADL(LE(I,K))
  348.       G(NF,N)=G(NF,N)+BF(KP,I)
  349. 90006 CONTINUE
  350.    11 CONTINUE
  351. 90001 CONTINUE
  352. C     WRITE(6,*)' GG(I)='
  353. C     WRITE(6,1002)GG
  354.    60 CONTINUE
  355. C     CALL ARRET(0)
  356.       IPAS=1
  357.       RETURN
  358.  1001 FORMAT('  XCVTI',I10,6E12.5)
  359.  1002 FORMAT(10(1X,1PE11.4))
  360.       END
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  

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