C PORMAO SOURCE CHAT 05/01/13 02:17:57 5004 SUBROUTINE PORMAO(VELA,MATE,IFOU,IDIM,TXR,XLOC, . XGLOB,D1HOOK,ROTHOO,DDHOOK,LHOOK,COBMA,XMOB,KCAS,IRET) C----------------------------------------------------------------------- C C MATRICE DE HOOK DES ELEMENTS MASSIFS ORTHOTROPES OU ANISOTROPES C SPECIAL POUR MILIEU POREUX C C ENTREES C VELA() = materiau dans un tableau de travail C MATE = Nom du materiau C IFOU = num{ro d'harmonique de fourier: IFOUR de CCOPTIO C IDIM = DEFINIT SI ON EST EN 2D OU 3D C TXR = COS-DIRECTEURS DES AXES LOCAUX /REPERE GLOBAL C LHOOK = TAILLE DE LA MATRICE DE HOOKE C KCAS = 1 SI ON VEUT LA MATRICE POUR ELLE MEME C = 2 SI ON VEUT LA MATRICE POUR L'INVERSER ENSUITE C XLOC | C XGLOB |= TABLEAU DE TRAVAIL C D1HOOK| C C SORTIES C DDHOOK(LHOOK,LHOOK) = matrice de hooke C IRET = 1 si option existante 0 SINON C C INSPIRE DE DOHMAO C----------------------------------------------------------------------- C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(UNDEMI=.5D0,XZER=0.0D0,XK=1.2D0) PARAMETER(UN=1.D0,DEUX=2.D0) CHARACTER*8 MATE C DIMENSION VELA(*),DDHOOK(LHOOK,*),TXR(IDIM,*) DIMENSION XLOC(IDIM,*),D1HOOK(LHOOK,*) DIMENSION XGLOB(IDIM,*),ROTHOO(LHOOK,*) DIMENSION COBMA(*),COBAUX(10) C IRET=1 CALL ZERO(D1HOOK,LHOOK,LHOOK) CALL ZERO(XGLOB,IDIM,IDIM) CALL ZERO(ROTHOO,LHOOK,LHOOK) CALL ZERO(DDHOOK,LHOOK,LHOOK) CALL ZERO(COBMA,LHOOK,1) CALL ZERO(COBAUX,LHOOK,1) C IF(MATE.EQ.'ORTHOTRO')THEN C MATRICE DE HOOK / AXES ORTHO IF(IFOU.EQ.-2.AND.KCAS.EQ.1)THEN C CONTRAINTE PLANE ET KCAS = 1 YG1=VELA(1) YG2=VELA(2) XNU12=VELA(3) G12=VELA(4) COBAUX(1)=VELA(7) COBAUX(2)=VELA(8) XMOB =VELA(9) XNU21=(YG2/YG1)*XNU12 AUX=UN-XNU12*XNU21 D1HOOK(1,1)=YG1/AUX D1HOOK(2,1)=XNU21*(YG1/AUX) D1HOOK(1,2)=D1HOOK(2,1) D1HOOK(2,2)=YG2/AUX D1HOOK(4,4)=G12 ELSE IF(IFOU.EQ.-2.AND.KCAS.EQ.2)THEN C CONTRAINTE PLANE ET KCAS=2 YG1=VELA(1) YG2=VELA(2) XNU12=VELA(3) YG3=VELA(10) XNU23=VELA(11) XNU13=VELA(12) G12=VELA(4) COBAUX(1)=VELA(7) COBAUX(2)=VELA(8) XMOB =VELA(9) XNU21=(YG2/YG1)*XNU12 XNU32=(YG3/YG2)*XNU23 XNU31=(YG3/YG1)*XNU13 AUX=(UN-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31 . -DEUX*XNU21*XNU32*XNU13) AUX1=AUX/YG1 AUX2=AUX/YG2 AUX3=AUX/YG3 D1HOOK(1,1)=(UN-XNU23*XNU32)/AUX1 D1HOOK(1,2)=(XNU21+XNU31*XNU23)/AUX1 D1HOOK(2,1)=D1HOOK(1,2) D1HOOK(1,3)=(XNU31+XNU21*XNU32)/AUX1 D1HOOK(3,1)=D1HOOK(1,3) D1HOOK(2,2)=(UN-XNU13*XNU31)/AUX2 D1HOOK(2,3)=(XNU32+XNU12*XNU31)/AUX2 D1HOOK(3,2)=D1HOOK(2,3) D1HOOK(3,3)=(UN-XNU12*XNU21)/AUX3 D1HOOK(4,4)=G12 ELSEIF(IFOU.EQ.-1.OR.IFOU.EQ.0) THEN C DEFORMATION PLANE,AXISYMETRIE YG1=VELA(1) YG2=VELA(2) YG3=VELA(3) XNU12=VELA(4) XNU23=VELA(5) XNU13=VELA(6) G12=VELA(7) COBAUX(1)=VELA(10) COBAUX(2)=VELA(11) COBAUX(3)=VELA(12) XMOB =VELA(13) XNU21=(YG2/YG1)*XNU12 XNU32=(YG3/YG2)*XNU23 XNU31=(YG3/YG1)*XNU13 AUX=(UN-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31 . -DEUX*XNU21*XNU32*XNU13) AUX1=AUX/YG1 AUX2=AUX/YG2 AUX3=AUX/YG3 D1HOOK(1,1)=(UN-XNU23*XNU32)/AUX1 D1HOOK(1,2)=(XNU21+XNU31*XNU23)/AUX1 D1HOOK(2,1)=D1HOOK(1,2) D1HOOK(1,3)=(XNU31+XNU21*XNU32)/AUX1 D1HOOK(3,1)=D1HOOK(1,3) D1HOOK(2,2)=(UN-XNU13*XNU31)/AUX2 D1HOOK(2,3)=(XNU32+XNU12*XNU31)/AUX2 D1HOOK(3,2)=D1HOOK(2,3) D1HOOK(3,3)=(UN-XNU12*XNU21)/AUX3 D1HOOK(4,4)=G12 ELSEIF(IFOU.EQ.1.OR.IFOU.EQ.2)THEN C AXISYMETRIE DE FOURIER ET 3D YG1=VELA(1) YG2=VELA(2) YG3=VELA(3) XNU12=VELA(4) XNU23=VELA(5) XNU13=VELA(6) G12=VELA(7) G23=VELA(8) G13=VELA(9) * IF(IFOU.EQ.1) THEN COBAUX(1)=VELA(12) COBAUX(2)=VELA(13) COBAUX(3)=VELA(14) XMOB =VELA(15) ELSE IF (IFOU.EQ.2) THEN COBAUX(1)=VELA(16) COBAUX(2)=VELA(17) COBAUX(3)=VELA(18) XMOB =VELA(19) ENDIF * XNU21=(YG2/YG1)*XNU12 XNU32=(YG3/YG2)*XNU23 XNU31=(YG3/YG1)*XNU13 AUX=(UN-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31 . -DEUX*XNU21*XNU32*XNU13) AUX1=AUX/YG1 AUX2=AUX/YG2 AUX3=AUX/YG3 D1HOOK(1,1)=(UN-XNU23*XNU32)/AUX1 D1HOOK(1,2)=(XNU21+XNU31*XNU23)/AUX1 D1HOOK(2,1)=D1HOOK(1,2) D1HOOK(1,3)=(XNU31+XNU21*XNU32)/AUX1 D1HOOK(3,1)=D1HOOK(1,3) D1HOOK(2,2)=(UN-XNU13*XNU31)/AUX2 D1HOOK(2,3)=(XNU32+XNU12*XNU31)/AUX2 D1HOOK(3,2)=D1HOOK(2,3) D1HOOK(3,3)=(UN-XNU12*XNU21)/AUX3 D1HOOK(4,4)=G12 D1HOOK(5,5)=G13 D1HOOK(6,6)=G23 ENDIF C DEFINITION DES AXES ORTHO./AXES LOCAUX IF(IDIM.EQ.2)THEN IF(IFOU.EQ.-2)THEN XLOC(1,1)=VELA(5) XLOC(2,1)=VELA(6) ELSEIF(IFOU.EQ.-1.OR.IFOU.EQ.0)THEN XLOC(1,1)=VELA(8) XLOC(2,1)=VELA(9) ELSEIF(IFOU.EQ.1)THEN XLOC(1,1)=VELA(10) XLOC(2,1)=VELA(11) ENDIF XLOC(1,2)=-XLOC(2,1) XLOC(2,2)=XLOC(1,1) ELSEIF(IDIM.EQ.3)THEN XLOC(1,1)=VELA(10) XLOC(2,1)=VELA(11) XLOC(3,1)=VELA(12) XLOC(1,2)=VELA(13) XLOC(2,2)=VELA(14) XLOC(3,2)=VELA(15) C CALL CROSS2(XLOC(1,1),XLOC(1,2),XLOC(1,3),IRR) ENDIF ELSEIF(MATE.EQ.'ANISOTRO')THEN C C MATRICE DE HOOK / AXES ORTHO IF(IFOU.EQ.-2.AND.KCAS.EQ.1)THEN C CONTRAINTE PLANE ET KCAS = 1 D11=VELA(1) D21=VELA(2) D22=VELA(3) D31=VELA(13) D32=VELA(14) D33=VELA(15) IF (D33.EQ.0.D0) D33=D11*1.E0-6 D41=VELA(4) D42=VELA(5) D43=VELA(16) D44=VELA(6) D1HOOK(1,1)=D11 - ((D31**2)/D33) D1HOOK(2,1)=D21 - ((D31*D32)/D33) D1HOOK(1,2)=D1HOOK(2,1) D1HOOK(2,2)=D22 - ((D32**2)/D33) D1HOOK(4,1)=D41 - ((D31*D43)/D33) D1HOOK(1,4)=D1HOOK(4,1) D1HOOK(4,2)=D42 - ((D32*D43)/D33) D1HOOK(2,4)=D1HOOK(4,2) D1HOOK(4,4)=D44 - ((D43**2)/D33) COBAUX(1)=VELA( 9) COBAUX(2)=VELA(10) COBAUX(3)=0.D0 COBAUX(4)=VELA(11) XMOB =VELA(12) ELSE IF(IFOU.EQ.-2.AND.KCAS.EQ.2)THEN C CONTRAINTE PLANE ET KCAS = 2 D1HOOK(1,1)=VELA(1) D1HOOK(2,1)=VELA(2) D1HOOK(1,2)=D1HOOK(2,1) D1HOOK(2,2)=VELA(3) D1HOOK(3,1)=VELA(13) D1HOOK(1,3)=D1HOOK(3,1) D1HOOK(3,2)=VELA(14) D1HOOK(2,3)=D1HOOK(3,2) D1HOOK(3,3)=VELA(15) D1HOOK(4,1)=VELA(4) D1HOOK(1,4)=D1HOOK(4,1) D1HOOK(4,2)=VELA(5) D1HOOK(2,4)=D1HOOK(4,2) D1HOOK(4,3)=VELA(16) D1HOOK(3,4)=D1HOOK(4,3) D1HOOK(4,4)=VELA(6) COBAUX(1)=VELA( 9) COBAUX(2)=VELA(10) COBAUX(3)=0.D0 COBAUX(4)=VELA(11) XMOB =VELA(12) ELSEIF(IFOU.EQ.-1.OR.IFOU.EQ.0) THEN C DEFORMATION PLANE,AXISYMETRIE D1HOOK(1,1)=VELA(1) D1HOOK(2,1)=VELA(2) D1HOOK(1,2)=D1HOOK(2,1) D1HOOK(2,2)=VELA(3) D1HOOK(3,1)=VELA(4) D1HOOK(1,3)=D1HOOK(3,1) D1HOOK(3,2)=VELA(5) D1HOOK(2,3)=D1HOOK(3,2) D1HOOK(3,3)=VELA(6) D1HOOK(4,1)=VELA(7) D1HOOK(1,4)=D1HOOK(4,1) D1HOOK(4,2)=VELA(8) D1HOOK(2,4)=D1HOOK(4,2) D1HOOK(4,3)=VELA(9) D1HOOK(3,4)=D1HOOK(4,3) D1HOOK(4,4)=VELA(10) COBAUX(1)=VELA(13) COBAUX(2)=VELA(14) COBAUX(3)=VELA(16) COBAUX(4)=VELA(15) XMOB =VELA(17) ELSEIF(IFOU.EQ.1.OR.IFOU.EQ.2)THEN C AXISYMETRIE DE FOURIER , 3D D1HOOK(1,1)=VELA(1) D1HOOK(2,1)=VELA(2) D1HOOK(1,2)=D1HOOK(2,1) D1HOOK(2,2)=VELA(3) D1HOOK(3,1)=VELA(4) D1HOOK(1,3)=D1HOOK(3,1) D1HOOK(3,2)=VELA(5) D1HOOK(2,3)=D1HOOK(3,2) D1HOOK(3,3)=VELA(6) D1HOOK(4,1)=VELA(7) D1HOOK(1,4)=D1HOOK(4,1) D1HOOK(4,2)=VELA(8) D1HOOK(2,4)=D1HOOK(4,2) D1HOOK(4,3)=VELA(9) D1HOOK(3,4)=D1HOOK(4,3) D1HOOK(4,4)=VELA(10) IF(IFOU.EQ.1)THEN D1HOOK(5,5)=VELA(11) D1HOOK(6,5)=VELA(12) D1HOOK(5,6)=D1HOOK(6,5) D1HOOK(6,6)=VELA(13) COBAUX(1)=VELA(16) COBAUX(2)=VELA(17) COBAUX(3)=VELA(19) COBAUX(4)=VELA(18) XMOB =VELA(20) ELSE D1HOOK(5,1)=VELA(11) D1HOOK(1,5)=D1HOOK(5,1) D1HOOK(5,2)=VELA(12) D1HOOK(2,5)=D1HOOK(5,2) D1HOOK(5,3)=VELA(13) D1HOOK(3,5)=D1HOOK(5,3) D1HOOK(5,4)=VELA(14) D1HOOK(4,5)=D1HOOK(5,4) D1HOOK(5,5)=VELA(15) D1HOOK(6,1)=VELA(16) D1HOOK(1,6)=D1HOOK(6,1) D1HOOK(6,2)=VELA(17) D1HOOK(2,6)=D1HOOK(6,2) D1HOOK(6,3)=VELA(18) D1HOOK(3,6)=D1HOOK(6,3) D1HOOK(6,4)=VELA(19) D1HOOK(4,6)=D1HOOK(6,4) D1HOOK(6,5)=VELA(20) D1HOOK(5,6)=D1HOOK(6,5) D1HOOK(6,6)=VELA(21) COBAUX(1)=VELA(28) COBAUX(2)=VELA(29) COBAUX(3)=VELA(30) COBAUX(4)=VELA(31) COBAUX(5)=VELA(32) COBAUX(6)=VELA(33) XMOB =VELA(34) ENDIF ENDIF C DEFINITION DES AXES ORTHO./AXES LOCAUX IF(IDIM.EQ.2)THEN IF(IFOU.EQ.-2)THEN XLOC(1,1)=VELA(7) XLOC(2,1)=VELA(8) ELSEIF(IFOU.EQ.-1.OR.IFOU.EQ.0)THEN XLOC(1,1)=VELA(11) XLOC(2,1)=VELA(12) ELSEIF(IFOU.EQ.1)THEN XLOC(1,1)=VELA(14) XLOC(2,1)=VELA(15) ENDIF XLOC(1,2)=-XLOC(2,1) XLOC(2,2)=XLOC(1,1) ELSEIF(IDIM.EQ.3)THEN XLOC(1,1)=VELA(22) XLOC(2,1)=VELA(23) XLOC(3,1)=VELA(24) XLOC(1,2)=VELA(25) XLOC(2,2)=VELA(26) XLOC(3,2)=VELA(27) C CALL CROSS2(XLOC(1,1),XLOC(1,2),XLOC(1,3),IRR) ENDIF ELSEIF(MATE.EQ.'UNIDIREC')THEN C MATRICE DE HOOKE /AXES D'ARMAURE D1HOOK(1,1)=VELA(1) COBAUX(1)=VELA(4) XMOB =VELA(5) C DEFINITION DES AXES ORTHO./AXES LOCAUX IF(IDIM.EQ.2)THEN XLOC(1,1)=VELA(2) XLOC(2,1)=VELA(3) XLOC(1,2)=-XLOC(2,1) XLOC(2,2)=XLOC(1,1) ELSEIF(IDIM.EQ.3)THEN XLOC(1,1)=VELA(2) XLOC(2,1)=VELA(3) XLOC(3,1)=VELA(4) XLOC(1,2)=VELA(5) XLOC(2,2)=VELA(6) XLOC(3,2)=VELA(7) C CALL CROSS2(XLOC(1,1),XLOC(1,2),XLOC(1,3),IRR) ENDIF ELSE IRET=0 ENDIF C C DEFINITION DES AXES ORTHO./AXES GLOBAUX C IF(IRET.EQ.1)THEN C DO 90 K=1,IDIM DO 90 J=1,IDIM DO 90 I=1,IDIM XGLOB(K,J)=TXR(J,I)*XLOC(I,K)+XGLOB(K,J) 90 CONTINUE C MATRICE DE TRANSFORMATION IF(IDIM.EQ.2)THEN ROTHOO(1,1)=XGLOB(1,1)*XGLOB(1,1) ROTHOO(1,2)=XGLOB(1,2)*XGLOB(1,2) ROTHOO(1,4)=XGLOB(1,1)*XGLOB(1,2) ROTHOO(2,1)=XGLOB(2,1)*XGLOB(2,1) ROTHOO(2,2)=XGLOB(2,2)*XGLOB(2,2) ROTHOO(2,4)=XGLOB(2,1)*XGLOB(2,2) ROTHOO(3,3)=UN ROTHOO(4,1)=DEUX*XGLOB(1,1)*XGLOB(2,1) ROTHOO(4,2)=DEUX*XGLOB(1,2)*XGLOB(2,2) ROTHOO(4,4)=XGLOB(1,2)*XGLOB(2,1)+XGLOB(1,1)*XGLOB(2,2) IF(IFOU.EQ.1)THEN ROTHOO(5,5)=XGLOB(1,1) ROTHOO(5,6)=XGLOB(1,2) ROTHOO(6,5)=XGLOB(2,1) ROTHOO(6,6)=XGLOB(2,2) ENDIF ELSEIF(IDIM.EQ.3)THEN DO 100 IC=1,3 DO 100 IL=1,3 ROTHOO(IL,IC)=XGLOB(IL,IC)*XGLOB(IL,IC) 100 CONTINUE C DO 110 IL=1,3 ROTHOO(IL,4)=XGLOB(IL,1)*XGLOB(IL,2) ROTHOO(IL,5)=XGLOB(IL,2)*XGLOB(IL,3) ROTHOO(IL,6)=XGLOB(IL,1)*XGLOB(IL,3) 110 CONTINUE C DO 120 IC=1,3 ROTHOO(4,IC)=DEUX*XGLOB(1,IC)*XGLOB(2,IC) ROTHOO(5,IC)=DEUX*XGLOB(2,IC)*XGLOB(3,IC) ROTHOO(6,IC)=DEUX*XGLOB(1,IC)*XGLOB(3,IC) 120 CONTINUE C DO 130 IL=4,6 IL1=IL-3 IL2=IL1+1 IF(IL2.GT.3)IL2=IL2-3 DO 130 IC=4,6 IC1=IC-3 IC2=IC1+1 IF(IC2.GT.3)IC2=IC2-3 ROTHOO(IL,IC)=XGLOB(IL1,IC1)*XGLOB(IL2,IC2)+ . XGLOB(IL1,IC2)*XGLOB(IL2,IC1) 130 CONTINUE DO 135 IC=1,6 AA=ROTHOO(6,IC) ROTHOO(6,IC)=ROTHOO(5,IC) ROTHOO(5,IC)=AA 135 CONTINUE DO 136 IL=1,6 AA=ROTHOO(IL,6) ROTHOO(IL,6)=ROTHOO(IL,5) ROTHOO(IL,5)=AA 136 CONTINUE ENDIF C C TRANSFORMATION DE LA MATRICE DE HOOKE ET DE COBAUX C CALL PRODT(DDHOOK,D1HOOK,ROTHOO,LHOOK,LHOOK) * DO 140 IL=1,LHOOK DO 140 IC=1,LHOOK COBMA(IL)=COBMA(IL)+ROTHOO(IC,IL)*COBAUX(IC) 140 CONTINUE ENDIF RETURN END