Télécharger la plaquette Cast3M


fr en

Télécharger yclpls.eso

Retour à la liste

Numérotation des lignes : 

yclpls

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

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