C ZCLPLS SOURCE CHAT 05/01/13 04:21:27 5004 SUBROUTINE ZCLPLS(HR,RPG,LE,NEL,K0,NPTD,NX,IES,NP,IAXI,IPADL, & COEFF,IK1, & TN,G, & VOLU,COTE,NELZ,DME,DIAM, & DT,DTT2,NUEL,DIAEL) C C VERSION VECTORISEE (CF XCVTIT POUR PLUS DE DETAILS) C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C NOMBRE MAXI DE POINTS PAR ELEMENT : NPX PARAMETER(NPX=9) C LONGUEUR DES REGISTRES VECTORIELS DE LA MACHINE CIBLE PARAMETER(LRV=64) C*********************************************************************** C C CE SP DISCRETISE LE LAPLACIEN D UNE VARIABLE SCALAIRE C C EN 1D SUR L ELEMENT SEG2 C EN 2D SUR L ELEMENT QUA4 PLAN OU AXI C EN 3D SUR L ELEMENT CUB8 PRI6 et TET4 C LES OPERATEURS SONT SOUS INTEGRES C C COEFF(SCAL DOMA) LE COEFFICIENT DU LAPLACIEN ( NU ) C (SCAL ELEM) C C NUELD(NELZ): NUMERO DE L'ELEMENT NK DE LA ZONE DANS LE DOMAINE C (NELZ : NOMBRE D'ELEMENTS DE LA ZONE) C ROC(NELD) : ROC PAR ELEMENT DONNE SUR TOUT LE DOMAINE C C C*********************************************************************** -INC CCVQUA4 -INC CCREEL C DIMENSION TN(NPTD,NX),GG(8) DIMENSION COEFF(*) DIMENSION COTE(NELZ,IES),DME(*),VOLU(*),DIAM(*) DIMENSION IPADL(1),LE(NP,1) DIMENSION HR(NEL,NP,IES),RPG(1) REAL*8 G(NPTD,NX) REAL*8 ROC(1) DIMENSION NUELD(1) DIMENSION QGGT(8,8),Q1(8,8),Q2(8,8),Q3(8,8) C C C* -TABLEAUX ADDITIONNELS POUR LA VECTORISATION ********************** DIMENSION AIRE(LRV),ALF (LRV),COEF(LRV),CLSR(LRV),GINC(LRV) DIMENSION AL (LRV),AH (LRV),AP (LRV),CFM (LRV),XMA (LRV) DIMENSION XMB(LRV),XMD(LRV),XMH(LRV),XMI(LRV),AHL(LRV) DIMENSION ALH(LRV),DPR(LRV),TETA(LRV,NPX,3),F(LRV,NPX),BF(LRV,NPX) C C --- TABLEAU POUR L'OPTION RAPIDE --- DIMENSION ALP(LRV) C*** SAVE IPAS,QGGT,Q1,Q2,Q3 DATA IPAS/0/ C C INITIALISATIONS DIVERSES C C WRITE(6,*)' ROC=',(ROC(MM),MM=1,NEL) C WRITE(6,*)' NUELD=',(NUELD(MM),MM=1,NEL) NK=K0 CALL INITD(BF,64,0.D0) IF(IES.EQ.2)GO TO 20 IF(IES.EQ.3)GO TO 30 C********** C * 1D* C********** C C C K EST LE NUMERO DE L'ELEMENT C KP PERMET DE SE SITUER A L'INTERIEUR DES TABLEAUX DE LRV ELEMENTS NPACK=INT(NEL/FLOAT(LRV))+1 KPACKD=1 KPACKF=NPACK C C*** BOUCLE SUR LES ELEMENTS *** C DO 70001 KPACK=KPACKD,KPACKF C DO 40 K=1,NEL C C POUR CHAQUE PAQUET DE LRV ELEMENTS KDEB=1+(KPACK-1)*LRV KFIN=MIN(NEL,KDEB+LRV-1) DO 70002 K=KDEB,KFIN KP=K-KDEB+1 NK=K+K0 K1=1+(1-IK1)*(NK-1) COEF(KP)=COEFF(K1) CLSR(KP)=COEFF(K1) AIRE(KP)=VOLU(NK) XMA(KP)=AIRE(KP) 70002 CONTINUE DO 4 I=1,NP DO 4 N=1,NX DO 4 K=KDEB,KFIN KP=K-KDEB+1 NF=IPADL(LE(I,K)) TETA(KP,I,N)=TN(NF,N) 4 CONTINUE C C C DO 70003 K=KDEB,KFIN NK=K+K0 KP=K-KDEB+1 DT0=DT DT2=0.5*XMA(KP)*XMA(KP)/CLSR(KP) IF(DT2.LT.DT)DT=DT2 IF(DT.EQ.DT0)GO TO 152 DTT2=DT2 DIAEL=XMA(KP) NUEL=NK 152 CONTINUE 70003 CONTINUE C C N=1 C --- CALCUL DE L'INCREMENT GINC --------------------------- DO 70104 K=KDEB,KFIN KP=K-KDEB+1 GINC(KP)=COEF(KP)*(TETA(KP,1,N)-TETA(KP,2,N))/AIRE(KP) 70104 CONTINUE C --- ACCUMULATION DANS LE TABLEAU G ----------------------- DO 70004 K=KDEB,KFIN KP=K-KDEB+1 NF=IPADL(LE(1,K)) G(NF,N)=G(NF,N) + GINC(KP) NF=IPADL(LE(2,K)) G(NF,N)=G(NF,N) - GINC(KP) 70004 CONTINUE 70001 CONTINUE 40 CONTINUE IPAS=1 RETURN C *********** C * 2D * C *********** 20 CONTINUE C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DIFFERENCES TRIANGLE / QUADRANGLE QUA4=0.D0 IF(NP.EQ.4)QUA4=1.D0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C K EST LE NUMERO DE L'ELEMENT C KP PERMET DE SE SITUER A L'INTERIEUR DES TABLEAUX DE LRV ELEMENTS NPACK=INT(NEL/FLOAT(LRV))+1 KPACKD=1 KPACKF=NPACK C C*** BOUCLE SUR LES ELEMENTS *** C DO 80001 KPACK=KPACKD,KPACKF C C POUR CHAQUE PAQUET DE LRV ELEMENTS KDEB=1+(KPACK-1)*LRV KFIN=MIN(NEL,KDEB+LRV-1) DO 80002 K=KDEB,KFIN KP=K-KDEB+1 NK=K+K0 K1=1+(1-IK1)*(NK-1) XMI(KP)=DIAM(NK) XMA(KP)=DME(NK) COEF(KP)=COEFF(K1) CLSR(KP)=COEFF(K1) AIRE(KP)=VOLU(NK) CC AL(KP)=COTE(NK,1) AH(KP)=COTE(NK,2) XMB(KP)=(AL(KP)+AH(KP))/2. XMD(KP)=((AL(KP)*AH(KP))**2)/(AL(KP)**2+AH(KP)**2) AHL(KP)=AH(KP)/AL(KP) ALH(KP)=AL(KP)/AH(KP) DPR(KP)=1.D0 IF(IAXI.NE.0)DPR(KP)=2.D0*XPI*RPG(K) 80002 CONTINUE CC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 5 I=1,NP DO 5 N=1,NX DO 5 K=KDEB,KFIN KP=K-KDEB+1 NF=IPADL(LE(I,K)) TETA(KP,I,N)=TN(NF,N) 5 CONTINUE C C C DO 80004 K=KDEB,KFIN KP=K-KDEB+1 NK=K+K0 DT0=DT DT2=0.5*XMD(KP)/CLSR(KP) IF(DT2.LT.DT)DT=DT2 IF(DT.EQ.DT0)GO TO 252 DTT2=DT2 DIAEL=XMB(KP) NUEL=NK 252 CONTINUE 80004 CONTINUE C C --- VERSION RAPIDE ----------------------- CF12=1.D0/12.D0 DO 80011 K=KDEB,KFIN KP=K-KDEB+1 ALP(KP)=CF12*(AHL(KP)+ALH(KP))*DPR(KP)*QUA4 80011 CONTINUE C ------------------------------------------ C DO 1 N=1,NX DO 80010 I= 1,NP DO 80010 K=KDEB,KFIN KP=K-KDEB+1 BF(KP,I)=0.D0 80010 CONTINUE DO 1 I=1,NP DO 3 J= 1,NP DO 3 K=KDEB,KFIN KP=K-KDEB+1 BF(KP,I)=BF(KP,I)+TETA(KP,J,N)* & ((HR(K,I,1)*HR(K,J,1)+HR(K,I,2)*HR(K,J,2))*AIRE(KP) & +ALP(KP)*VGGT(J,I))*COEF(KP) 3 CONTINUE C --- ACCUMULATION DANS LE TABLEAU G ------------------------------ DO 80006 K=KDEB,KFIN KP=K-KDEB+1 NF=IPADL(LE(I,K)) G(NF,N)=G(NF,N)+BF(KP,I) 80006 CONTINUE 1 CONTINUE 80001 CONTINUE 50 CONTINUE C WRITE(6,*)' SUB XCLPLS G(1,=' C WRITE(6,1002) (G(MM,1),MM=1,NPTD) C WRITE(6,*)' SUB XCLPLS G(2,=' C WRITE(6,1002) (G(MM,2),MM=1,NPTD) C WRITE(6,*)' FIN **** ' IPAS=1 RETURN 30 CONTINUE C *********** C * 3D * C *********** IF(IPAS.EQ.0)CALL CALHRH(QGGT,Q1,Q2,Q3,IES) CUB8=0.D0 IF(NP.EQ.8)CUB8=1.D0 1003 FORMAT(' XCVTI ',10I10) C C K EST LE NUMERO DE L'ELEMENT C KP PERMET DE SE SITUER A L'INTERIEUR DES TABLEAUX DE LRV ELEMENTS NPACK=INT(NEL/FLOAT(LRV))+1 KPACKD=1 KPACKF=NPACK C C*** BOUCLE SUR LES ELEMENTS *** C DO 90001 KPACK=KPACKD,KPACKF C DO 60 K=1,NEL C C POUR CHAQUE PAQUET DE LRV ELEMENTS KDEB=1+(KPACK-1)*LRV KFIN=MIN(NEL,KDEB+LRV-1) DO 90002 K=KDEB,KFIN KP=K-KDEB+1 NK=K+K0 K1=1+(1-IK1)*(NK-1) AL(KP)=COTE(NK,1) AH(KP)=COTE(NK,2) AP(KP)=COTE(NK,3) CFM(KP)=AL(KP)*AH(KP)/AP(KP)+AL(KP)*AP(KP)/AH(KP)+ & AP(KP)*AH(KP)/AL(KP) XMI(KP)=DIAM(NK) XMA(KP)=DME(NK) COEF(KP)=COEFF(K1) C CLSR(KP)=COEFF(K1)/ROC(NUELD(NK)) CLSR(KP)=COEFF(K1) AIRE(KP)=VOLU(NK) XMB(KP)=(XMA(KP)+XMI(KP))*0.5 XMD(KP)=((XMA(KP)*XMI(KP))**2)/(2.*XMA(KP)**2+XMI(KP)**2) XMH(KP)=AIRE(KP)**0.33 90002 CONTINUE C DO 15 I=1,NP DO 15 N=1,NX DO 15 K=KDEB,KFIN KP=K-KDEB+1 NF=IPADL(LE(I,K)) TETA(KP,I,N)=TN(NF,N) 15 CONTINUE C WRITE(6,*)' K ** ',K,' TETA ',' VOLU=',AIRE,' COEF=',COEF C &,' XMH=',XMH C WRITE(6,1002)(TETA(MM,1),MM=1,8) C C DO 90003 K=KDEB,KFIN KP=K-KDEB+1 NK=K+K0 DT0=DT DT2=0.5*XMD(KP)/CLSR(KP) IF(DT2.LT.DT)DT=DT2 IF(DT.EQ.DT0)GO TO 353 DTT2=DT2 DIAEL=XMB(KP) NUEL=NK 353 CONTINUE 90003 CONTINUE C DO 11 N=1,NX DO 90004 I=1,NP DO 90004 K=KDEB,KFIN KP=K-KDEB+1 BF(KP,I)=0. 90004 CONTINUE DO 11 I=1,NP DO 13 J=1,NP DO 13 K=KDEB,KFIN KP=K-KDEB+1 GEO1=CFM(KP)*QGGT(J,I)*CUB8 BF(KP,I)=BF(KP,I)+TETA(KP,J,N)* &(HR(K,I,1)*HR(K,J,1)+HR(K,I,2)*HR(K,J,2)+HR(K,I,3)*HR(K,J,3)) & *AIRE(KP)+COEF(KP)*XMH(KP)*GEO1 13 CONTINUE DO 90006 K=KDEB,KFIN KP=K-KDEB+1 NF=IPADL(LE(I,K)) G(NF,N)=G(NF,N)+BF(KP,I) 90006 CONTINUE 11 CONTINUE 90001 CONTINUE C WRITE(6,*)' GG(I)=' C WRITE(6,1002)GG 60 CONTINUE C CALL ARRET(0) IPAS=1 RETURN 1001 FORMAT(' XCVTI',I10,6E12.5) 1002 FORMAT(10(1X,1PE11.4)) END