C LISPS3 SOURCE CHAT 05/01/13 01:23:47 5004 SUBROUTINE LISPS3(XE,EPA1,FISS1,V1,EPA2,FISS2,V2,D,XDDL,XEL,BPSS, 1 XDDLOC,NPOINT,XSTRS,I69,I70,I195,I157) C======================================================================= C C EBERSOLT AVRIL 87 ELEMENT LISM C UHLMANN OCTOBRE 94 C C ENTREES C ENTREES C XE(3,4) = COORDONNEES DE LA POUTRE LINESPRING C EPA1 = EPAISSEUR NOEUDS 1 4 C EPA2 = EPAISSEUR NOEUDS 2 3 C FISS1 = PROFONDEUR DE LA FISSURE NOEUDS 1 4 C FISS2 = PROFONDEUR DE LA FISSURE NOEUDS 2 3 C V1(3) = VECTEUR ORIENTANT LES NOEUDS 1 4 C V2(3) = VECTEUR ORIENTANT LES NOEUDS 2 3 C D(2,2) = MATRICE DE HOOKE C XDDL(24) = D.D.L. DU LINESPRING C NPOINT = NOMBRE DE POINTS DE CONTRAINTES C TABLEAU DE TRAVAIL C XEL(3,3) = COORDONNEES LOCALES C BPSS(3,3) = MATRICE DE PASSAGE C XDDLOC(24) = D.D.L. LOCAUX C SORTIES C XSTRS(NPOINT*6)= CONTRAINTES LINESPRING C I69 = FISSURE TOTALEMENT TRAVERSANTE CONTRAINTES MISES A 0 C I70 = INDICERNABILITE DES 2 LEVRES C I343 = PROFONDEUR DE FISSURE NULLE C I157 = LES 2 LEVRES SONT TROP ELOIGNEES C C======================================================================= C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) C Include contenant quelques constantes dont XPI : -INC CCREEL PARAMETER(UNDEMI=.5D0,SIX=6.D0) PARAMETER(DOUZE=12.D0,TRSIX=36.D0) PARAMETER(EPS=2.D-3,PENA=1.D6,EPSINV=1.D-3,PENB=1.D2) PARAMETER(X774=.774596669241483D0) C DIMENSION XE(3,*),D(2,*),XSTRS(*),V1(*),V2(*),BPSS(3,*),XEL(3,*) DIMENSION S(3),XDDL(*),XDDLOC(*) C IF(NPOINT.EQ.1) THEN S(1)= XZERO ELSE IF(NPOINT.EQ.3) THEN S(1)=-X774 S(2)= XZERO S(3)= X774 ENDIF C C ON RECUPERE LES FISS AUX POINTS DE GAUSS IL FAUT LES CALCULER C AUX EXTREMITES C FISS1 = (FISS1*(UNDEMI +UNDEMI/X774))+(FISS2*(UNDEMI-UNDEMI/X774)) FISS2 = (FISS1*(UNDEMI -UNDEMI/X774))+(FISS2*(UNDEMI+UNDEMI/X774)) C NPOIN6=6*NPOINT C C MISE A ZERO DES CONTRAINTES DES INDICATEURS D ERREURS ET DES D.D. C DO 90 IA=1,NPOIN6 XSTRS(IA)=XZERO 90 CONTINUE I69 =0 I70 =0 C I343=0 I157=0 I195=0 C IF(FISS1.LT.XZERO) THEN C I343=1 FISS1=XZERO ENDIF C IF(FISS2.LT.XZERO) THEN C I343=1 FISS2=XZERO ENDIF C C EXTRACTION DE LA MATRICE DE PASSAGE C DO 100 IA=1,3 XEL(IA,1)=XE(IA,1) XEL(IA,2)=XE(IA,2) XEL(IA,3)=XE(IA,1)+(V1(IA)+V2(IA))*UNDEMI 100 CONTINUE C C----------------------------------------------- CALL VPAST(XEL,BPSS) CALL MATVEC(XDDL,XDDLOC,BPSS,8) C----------------------------------------------- C système local: x,y et z C Noeud 1 (_NO1): C XNO1=(XE(1,1)*BPSS(1,1))+(XE(2,1)*BPSS(1,2))+(XE(3,1)*BPSS(1,3)) C YNO1=(XE(1,1)*BPSS(2,1))+(XE(2,1)*BPSS(2,2))+(XE(3,1)*BPSS(2,3)) C ZNO1=(XE(1,1)*BPSS(3,1))+(XE(2,1)*BPSS(3,2))+(XE(3,1)*BPSS(3,3)) C Noeud2 (_NO2): C XNO2=(XE(1,2)*BPSS(1,1))+(XE(2,2)*BPSS(1,2))+(XE(3,2)*BPSS(1,3)) C YNO2=(XE(1,2)*BPSS(2,1))+(XE(2,2)*BPSS(2,2))+(XE(3,2)*BPSS(2,3)) ZNO2=(XE(1,2)*BPSS(3,1))+(XE(2,2)*BPSS(3,2))+(XE(3,2)*BPSS(3,3)) C Noeud 3 (_NO3): C XNO3=(XE(1,3)*BPSS(1,1))+(XE(2,3)*BPSS(1,2))+(XE(3,3)*BPSS(1,3)) C YNO3=(XE(1,3)*BPSS(2,1))+(XE(2,3)*BPSS(2,2))+(XE(3,3)*BPSS(2,3)) ZNO3=(XE(1,3)*BPSS(3,1))+(XE(2,3)*BPSS(3,2))+(XE(3,3)*BPSS(3,3)) C C Direction z: Différence entre noeud 2 et noeud 3: DZ23 = ZNO2 - ZNO3 C ________________________________________________________________ C DJA1=XZERO DJA2=XZERO DO 105 IA=1,3 DJA1=DJA1+(XE(IA,1)-XE(IA,4))*BPSS(3,IA) DJA2=DJA2+(XE(IA,2)-XE(IA,3))*BPSS(3,IA) 105 CONTINUE DJAC=DJA1*DJA2 IF(DJAC.LT.-1.D-20) I195=1 C HAUT = LARGEUR ENTRE LES NOEUDS 1,4 ET 2,3 C HAUT C = LARGEUR ENTRE LES NOEUDS 1,4 ET 2,3 HAUT=XZERO XLARG1=XZERO XLARG2=XZERO DO 110 IA=1,3 HAUT =(XE(IA,2)-XE(IA,1))*(XE(IA,2)-XE(IA,1))+HAUT XLARG1=(XE(IA,4)-XE(IA,1))*(XE(IA,4)-XE(IA,1))+XLARG1 XLARG2=(XE(IA,3)-XE(IA,2))*(XE(IA,3)-XE(IA,2))+XLARG2 110 CONTINUE HAUT =SQRT(HAUT) XLARG1=SQRT(XLARG1) XLARG2=SQRT(XLARG2) EPS1=XLARG1/HAUT EPS2=XLARG2/HAUT IF(EPS1.GT.EPS.OR.EPS2.GT.EPS) I157=1 DJA1=DJA1/HAUT DJA2=DJA2/HAUT IF(DJA1.LT.1.D-6.AND.DJA2.LT.1.D-6) I70=1 C C ASURW = A / W NOTATION CHEISSOUX C W=(EPA1+EPA2)*UNDEMI ASURW=(FISS1+FISS2)/W ASUR1=FISS1/W ASUR2=FISS2/W IF(ASUR1.GT..98D0.AND.ASUR2.GT..98D0) I69=1 IF(I69.EQ.1) GOTO 666 C C TRANSFORMATION DE D D(1,1)=D(1,1)*W D(2,2)=D(2,2)*W*W*W/DOUZE C PEWM=D(1,1)*PENA PEWF=D(2,2)*PENA PEWMB=D(1,1)*PENB PEWFB=D(2,2)*PENA IF (ASURW.GT.EPSINV) THEN C C BOUCLE SUR LES POINTS DE CONTRAINTES C DO 200 IA=1,NPOINT H1=UNDEMI-UNDEMI*S(IA) H2=UNDEMI+UNDEMI*S(IA) ASURW=(H1*FISS1+H2*FISS2)/W A = H1*FISS1+H2*FISS2 C C ON RECUPERE LES COEFFICIENTS ALPHAS ET F I C-------------------------------------------------- CALL LISPAL(ASURW,ALMM,ALMF,ALFF,DELTA) CALL LISPFI(ASURW,FM,FF) C CALL INTER2D C-------------------------------------------------- DELTA=D(1,1)/DELTA C C CALCUL DES COEFFICIENTS R1 R2 R3 R4 C R1= DELTA*ALFF/W R2=-DELTA*ALMF/SIX R3=-DELTA*ALMF/SIX R4= DELTA*ALMM*W/TRSIX C C REMPLISSAGE DES CONTRAINTES C IF (DZ23.GE.0) THEN EPXZ=(XDDLOC(1 )-XDDLOC(19))*H1+(XDDLOC(7 )-XDDLOC(13)) $ *H2 EPYZ=(XDDLOC(2 )-XDDLOC(20))*H1+(XDDLOC(8 )-XDDLOC(14)) $ *H2 DD =(XDDLOC(3 )-XDDLOC(21))*H1+(XDDLOC(9 )-XDDLOC(15)) $ *H2 RT =(XDDLOC(4 )-XDDLOC(22))*H1+(XDDLOC(10)-XDDLOC(16)) $ *H2 RTZZ=(XDDLOC(6 )-XDDLOC(24))*H1+(XDDLOC(12)-XDDLOC(18)) $ *H2 ELSE EPXZ=(XDDLOC(19)-XDDLOC(1 ))*H1+(XDDLOC(13)-XDDLOC(7 )) $ *H2 EPYZ=(XDDLOC(20)-XDDLOC(2 ))*H1+(XDDLOC(14)-XDDLOC(8 )) $ *H2 DD =(XDDLOC(21)-XDDLOC(3 ))*H1+(XDDLOC(15)-XDDLOC(9 )) $ *H2 RT =(XDDLOC(22)-XDDLOC(4 ))*H1+(XDDLOC(16)-XDDLOC(10)) $ *H2 RTZZ=(XDDLOC(24)-XDDLOC(6 ))*H1+(XDDLOC(18)-XDDLOC(12)) $ *H2 ENDIF C ___________________________________________________________ C IX=6*(IA-1) C XSTRS(IX+1)=DD*R1+RT*R2 XSTRS(IX+2)=PEWMB*EPXZ XSTRS(IX+3)=PEWMB*EPYZ XSTRS(IX+4)=DD*R3+RT*R4 XSTRS(IX+5)=PEWFB*RTZZ C X1=XSTRS(IX+1)/W X4=XSTRS(IX+4)*SIX/(W*W) C ________________________________________ C C CALCUL DE K I ET DE J C XXX=XPI*A XXX=SQRT(XXX) XKIE =(X1*FM+X4*FF)*XXX XSTRS(IX+6)= XKIE C 200 CONTINUE C C SI FISSURE INEXISTANTE C C ON FAIT UN C.S.T. POUR L ELEMENT DU BOUT C ELSE IF(ASURW.LE.EPSINV) THEN C DO 150 IA=1,NPOINT H1=UNDEMI-UNDEMI*S(IA) H2=UNDEMI+UNDEMI*S(IA) EPXZ =(XDDLOC(1 )-XDDLOC(19))*H1+(XDDLOC(7 )-XDDLOC(13))*H2 EPYZ =(XDDLOC(2 )-XDDLOC(20))*H1+(XDDLOC(8 )-XDDLOC(14))*H2 C DD1 = XDDLOC(3 )-XDDLOC(21) DD2 = XDDLOC(9 )-XDDLOC(15) DD = DD1*H1 + DD2*H2 C RT1 = XDDLOC(4 )-XDDLOC(22) RT2 = XDDLOC(10)-XDDLOC(16) RT = RT1*H1 + RT2*H2 C RTZZ =(XDDLOC(6 )-XDDLOC(24))*H1+(XDDLOC(12)-XDDLOC(18))*H2 IX=6*(IA-1) XSTRS(IX+1)=PEWM*DD XSTRS(IX+2)=PEWMB*EPXZ XSTRS(IX+3)=PEWMB*EPYZ XSTRS(IX+4)=PEWF*RT XSTRS(IX+5)=PEWFB*RTZZ XSTRS(IX+6)=XZERO 150 CONTINUE ENDIF 666 CONTINUE RETURN END