Sources Esope | Cast3M

pormao
C PORMAO    SOURCE    CB215821  17/07/21    21:15:26     9513           
      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(XZER=0.0D0)
      PARAMETER(UN=1.D0,DEUX=2.D0)
      CHARACTER*8 MATE
C
      DIMENSION VELA(*),DDHOOK(LHOOK,*),TXR(IDIM,*)
      DIMENSION XLOC(3,3),D1HOOK(LHOOK,*)
      DIMENSION XGLOB(3,3),ROTHOO(LHOOK,*)
      DIMENSION COBMA(*),COBAUX(10)
C
      IRET=1
      CALL ZERO(D1HOOK,LHOOK,LHOOK)
      CALL ZERO(XGLOB,3,3)
      CALL ZERO(XLOC,3,3)
      CALL ZERO(ROTHOO,LHOOK,LHOOK)
      CALL ZERO(DDHOOK,LHOOK,LHOOK)
      CALL ZERO(COBMA,LHOOK,1)
      CALL ZERO(COBAUX,LHOOK,1)
C
C  1 - MATERIAU ORTHOTROPE
C =========================
      IF(MATE.EQ.'ORTHOTRO')THEN
C =====
C  1.1 - Definition de la matrice de HOOKE dans les axes d'ORTHOTROPIE
C =====
C= -> Contraintes planes et KCAS=1
        IF(IFOU.EQ.-2.AND.KCAS.EQ.1)THEN
          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
C= -> Contraintes planes et KCAS=2
        ELSEIF(IFOU.EQ.-2.AND.KCAS.EQ.2)THEN
          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
C= -> Deformations planes, Axisymetrie
        ELSEIF(IFOU.EQ.-1.OR.IFOU.EQ.0) THEN
          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
C= -> Serie de Fourier et 3D
        ELSEIF(IFOU.EQ.1.OR.IFOU.EQ.2)THEN
          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)
            ELSEIF (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 =====
C  1.2 - Definition de la matrice de passage des axes d'ORTHOTROPIE
C        aux axes LOCAUX au point considere (dimensions 2 et 3)
C =====
cbp     IF (IDIM.EQ.2) THEN
        IF (IDIM.EQ.2.AND.IFOU.NE.1) 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)
c          ELSEIF(IFOU.EQ.1)THEN
c           XLOC(1,1)=VELA(10)
c           XLOC(2,1)=VELA(11)
          ENDIF
          XLOC(1,2)=-XLOC(2,1)
          XLOC(2,2)=XLOC(1,1)
cbp     ELSEIF(IDIM.EQ.3)THEN
        ELSE
          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)
          CALL CROSS2(XLOC(1,1),XLOC(1,2),XLOC(1,3),IRR)
        ENDIF
 
C  2 - MATERIAU ANISOTROPE
C =========================
      ELSEIF(MATE.EQ.'ANISOTRO')THEN
C =====
C  2.1 - Definition de la matrice de HOOKE dans les axes d'ANISOTROPIE
C =====
C= -> Contraintes planes et KCAS=1
        IF(IFOU.EQ.-2.AND.KCAS.EQ.1)THEN
          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)
C= -> Contraintes planes et KCAS=2
        ELSEIF(IFOU.EQ.-2.AND.KCAS.EQ.2)THEN
          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)
C= -> Deformations planes et Axisymetrie
        ELSEIF(IFOU.EQ.-1.OR.IFOU.EQ.0) THEN
          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)
C= -> Serie de Fourier et 3D
        ELSEIF(IFOU.EQ.1.OR.IFOU.EQ.2)THEN
          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 =====
C  2.2 - Definition de la matrice de passage des axes d'ANSIOTROPIE
C        aux axes LOCAUX au point considere (dimensions 2 et 3)
C =====
        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)
          CALL CROSS2(XLOC(1,1),XLOC(1,2),XLOC(1,3),IRR)
        ENDIF
 
C  3 - MATERIAU UNIDIRECTIONNEL
C ==============================        
      ELSEIF(MATE.EQ.'UNIDIREC')THEN
C =====
C  3.1 - Definition de la matrice de HOOKE dans les axes d'ANISOTROPIE
C        (par rapport aux axes d'armatures)
C =====
        D1HOOK(1,1)=VELA(1)
        COBAUX(1)=VELA(4)
        XMOB    =VELA(5)
C =====
C  3.2 - Definition de la matrice de passage des axes d'ANISOTROPIE
C        aux axes LOCAUX au point considere (dimensions 2 et 3)
C =====
        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)
          CALL CROSS2(XLOC(1,1),XLOC(1,2),XLOC(1,3),IRR)
        ENDIF
 
C  4 - MATERIAUX NON PREVUS
C ==========================
      ELSE
        IRET=0
        RETURN
      ENDIF
C
C  5 - DEFINITION DE LA MATRICE DE PASSAGE DES AXES
C      D'ORTHO/ANISOTROPIE AUX AXES DU REPERE GLOBAL
C ===================================================
      IF (IDIM.EQ.1) RETURN
C
      IDIM2=IDIM
      IF(IFOU.EQ.1) IDIM2=3
      DO j=1,IDIM
        DO k=1,IDIM2
          cc=XZER
          DO i=1,IDIM
            cc=cc+TXR(j,i)*XLOC(i,k)
          ENDDO
          XGLOB(k,j)=cc
        ENDDO
      ENDDO
cbp   en 2D Fourier, vrai TXR = [TXR(2x2) [0] ; [0] 1]   
      IF (IFOU.EQ.1) THEN
        XGLOB(1,3)=XLOC(3,1)
        XGLOB(2,3)=XLOC(3,2)
        XGLOB(3,3)=XLOC(3,3)        
      ENDIF
 
C     MATRICE DE TRANSFORMATION
cbp      IF (IDIM.EQ.2) THEN
      IF (IDIM.EQ.2.AND.IFOU.NE.1) 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)
c          IF(IFOU.EQ.1)THEN
c           ROTHOO(5,5)=XGLOB(1,1)
c           ROTHOO(5,6)=XGLOB(1,2)
c           ROTHOO(6,5)=XGLOB(2,1)
c           ROTHOO(6,6)=XGLOB(2,2)
c          ENDIF
cbp        ELSEIF(IDIM.EQ.3)THEN
      ELSE
         DO 100 IC=1,3
         DO 100 IL=1,3
           ROTHOO(IL,IC)=XGLOB(IL,IC)*XGLOB(IL,IC)
 100     CONTINUE
         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
         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
         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  6 - CALCUL DE LA MATRICE DE HOOKE DANS LE REPERE GLOBAL
C      ET DE COBAUX
C =========================================================
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
 
 
      RETURN
      END