Sources Esope | Cast3M

trihm2
C TRIHM2    SOURCE    CHAT      05/01/13    03:47:09     5004
      SUBROUTINE TRIHM2(IGAU,ITEL,MFR,NBNO,XEL,SHPTOT,SHP,IFOU,
     #        NHARM,VKL12,VKL23,VKL33,POIGAU,ISDJC,LRE,REL,IRET)
C=======================================================================
C
C    CALCULE LES TERMES  EN  P * PI ,PI * (UR,RT) ,(UR,RT) *(UR,RT)
C        (UT,RR) * (UT,RR) , PI * (UT,RR)     DE   LA     MATRICE
C           MASSE  DANS  LE  CAS  AXISYMETRIQUE OU FOURIER  POUR
C                  LA FORMULATION (37) HOMOGENE
C=======================================================================
C  INPUT
C     IGAU=NUMERO DU POINT DE GAUSS
C     ITEL=NUMERO DE L ELEMENT DANS NOMTP
C     MFR =NUMERO DE LA FORMULATION
C     NBNO=NOMBRE DE NOEUDS
C     XEL =COORDONNEES  DE L ELEMENT
C     IFOU=IFOUR DE CCOPTIO
C     NHARM=NUMERO DU MODE DE FOURIER
C     VKL12=-((COEFPI*COEFPR)/(RHOF*C**2))*SFLU/SCEL
C     VKL23=(BET11+BET22)*COEFPI/(2.*SCEL)
C     VKL33=(RHOS*2.+RHOF*(BET11+BET22))/SCEL
C     POIGAU=MINTE.POIGAU(IGAU)
C     LRE =NOMBRE DE D.D.L DE LA MATRICE  DE  RIGIDITE
C     SHPTOT(6,NBNO,NBGAU)=FONCTIONS DE FORMES ET DERIVEES
C     ISDJC = INDICATEUR SUR LE SIGNE DU JACOBIEN
C  ZONE DE TRAVAIL
C     SHP(5,NBNO)=TABLEAU DE TRAVAIL
C OUTPUT
C     ISDJC = INDICATEUR SUR LE SIGNE DU JACOBIEN
C     REL=MATRICE DE MASSE
C     IRET:INDICATEUR = 1 : SUCCES
C                       0 : ECHEC  (ELEMENT MELE INCOMPATIBLE )
C                       2 : ECHEC  (JACOBIEN  NUL )
C                       3 :ECHEC   (ROUTINE N EST VALABLE QU
C                                          EN AXISYMETRIQUE OU FOURIER )
C=======================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XEL(3,*),SHP(6,*),SHPTOT(6,NBNO,*),REL(LRE,*)
      IF (ITEL.EQ.92)              GOTO 10
C
C     ERREUR : TYPE D' ELEMENT  INCOMPATIBLE AVEC LA FORMULATION
C
      IRET = 0
      GOTO 666
  10  CONTINUE
      IF (IFOU.EQ.0.OR.IFOU.EQ.1)  GOTO 11
C
C     MESSAGE D ERREUR :  ROUTINE  N  EST  VALABLE  QU  EN FOURIER
C                     OU   EN   AXISYMETRIQUE
C
      IRET = 3
      GOTO 666
  11  CONTINUE
C
C     ELEMENTS HOMOGENEISES TRIH   EN AXISYMETRIE OU EN FOURIER
C          NBDL = LRE/NBNO   NOMBRE DE D.D.L PAR NOEUD
C
      NBDL = LRE /NBNO
C
C     SHP(1,I) : FONCTION DE FORME
C     SHP(2,I) : DERIVEE % R DE LA FONCTION DE FORME
C     SHP(3,I) : DERIVEE % Z DE LA FONCTION DE FORME
C
      DO 101 NP=1,NBNO
      SHP(1,NP)=SHPTOT(1,NP,IGAU)
      SHP(2,NP)=SHPTOT(2,NP,IGAU)
      SHP(3,NP)=SHPTOT(3,NP,IGAU)
  101 CONTINUE
      CALL DEVOLU(XEL,SHP,MFR,NBNO,IFOU,NHARM,2,1.D0,RR,DJAC)
      IF (DJAC.EQ.0.) GOTO 667
      IF (DJAC.LT.0.) ISDJC = ISDJC + 1
C
C     SHP(4,I) : FONCTION DE FORME DE UR (DEPLACEMENTS)
C     SHP(5,I) : FONCTION DE FORME DE (DUR/DZ)  (ROTATIONS)
C
      SHP(4,1)=SHP(1,1)*(1.D0-SHP(1,2)*SHP(1,2)-SHP(1,3)*SHP(1,3))  +
     # SHP(1,1)*SHP(1,1)*(SHP(1,2)+SHP(1,3))
      SHP(4,2)=SHP(1,2)*(1.D0-SHP(1,1)*SHP(1,1)-SHP(1,3)*SHP(1,3))  +
     # SHP(1,2)*SHP(1,2)*(SHP(1,1)+SHP(1,3))
      SHP(4,3)=SHP(1,3)*(1.D0-SHP(1,1)*SHP(1,1)-SHP(1,2)*SHP(1,2))  +
     # SHP(1,3)*SHP(1,3)*(SHP(1,1)+SHP(1,2))
C
C     A1=SHP(2,1) , A2=SHP(2,2) , A3 = SHP(2,3)
C
      A1=XEL(2,2)-XEL(2,3)
      A2=XEL(2,3)-XEL(2,1)
      A3=XEL(2,1)-XEL(2,2)
      SHP(5,1)= SHP(1,1)*SHP(1,1)*(A2*SHP(1,3)-A3*SHP(1,2))  +
     # 0.5D0*SHP(1,1)*SHP(1,2)*SHP(1,3)*(A2-A3)
      SHP(5,2)= SHP(1,2)*SHP(1,2)*(A3*SHP(1,1)-A1*SHP(1,3))  +
     # 0.5D0*SHP(1,1)*SHP(1,2)*SHP(1,3)*(A3-A1)
      SHP(5,3)= SHP(1,3)*SHP(1,3)*(A1*SHP(1,2)-A2*SHP(1,1))  +
     # 0.5D0*SHP(1,1)*SHP(1,2)*SHP(1,3)*(A1-A2)
C
C     TERMES EN P * PI
C
      DJAC1 = ABS(DJAC)*POIGAU
      IX1=0
      IY1=0
      DO   102 IX=2,LRE ,NBDL
      IX1=IX1 + 1
      DO   103 IY=1,IX  ,NBDL
      IY1=IY1 + 1
      REL(IY,IX) = REL(IY,IX) + VKL12*DJAC1*SHP(1,IX1)*SHP(1,IY1)
      REL(IX,IY) = REL(IY,IX)
  103 CONTINUE
      IY1=0
  102 CONTINUE
      DO   104 IX=2+NBDL,LRE ,NBDL
      IX2=IX - NBDL
      DO   105 IY=1,IX2 ,NBDL
      REL(IY+1,IX-1) = REL(IY,IX)
      REL(IX-1,IY+1) = REL(IY+1,IX-1)
  105 CONTINUE
  104 CONTINUE
C
C     TERMES EN PI * (UR , RT )
C
      IX1=0
      IY1=0
      DO   106 IX=3,LRE ,NBDL
      IX1=IX1 + 1
      DO   107 IY=2,IX  ,NBDL
      IY1=IY1 + 1
      REL(IY,IX) = REL(IY,IX) + VKL23*DJAC1*SHP(2,IY1)*SHP(4,IX1)
      REL(IY,IX+1) = REL(IY,IX+1) + VKL23*DJAC1*SHP(2,IY1)*SHP(5,IX1)
      REL(IX,IY) = REL(IY,IX)
      REL(IX+1,IY) = REL(IY,IX+1)
  107 CONTINUE
      IY1=0
  106 CONTINUE
      IX1=1
      IY1=0
      DO   108 IX=2+NBDL,LRE ,NBDL
      IX1=IX1 + 1
      DO   109 IY=3,IX  ,NBDL
      IY1=IY1 + 1
      REL(IY,IX) = REL(IY,IX) + VKL23*DJAC1*SHP(2,IX1)*SHP(4,IY1)
      REL(IY+1,IX) = REL(IY+1,IX) + VKL23*DJAC1*SHP(2,IX1)*SHP(5,IY1)
      REL(IX,IY) = REL(IY,IX)
      REL(IX,IY+1) = REL(IY+1,IX)
  109 CONTINUE
      IY1=0
  108 CONTINUE
      IF ( IFOU.EQ.1) THEN
C
C     TERMES EN PI * (UT , RR )
C     NON NULS QU EN FOURIER
C
      DJAC2 = ABS(DJAC)*POIGAU
      VKL25 = -1.D0* VKL23*NHARM
      IX1=0
      IY1=0
      DO   110 IX=5,LRE ,NBDL
      IX1=IX1 + 1
      DO   111 IY=2,IX  ,NBDL
      IY1=IY1 + 1
      REL(IY,IX) = REL(IY,IX) + VKL25*DJAC2*SHP(1,IY1)*SHP(4,IX1)
      REL(IY,IX+1) = REL(IY,IX+1) + VKL25*DJAC2*SHP(1,IY1)*SHP(5,IX1)
      REL(IX,IY) = REL(IY,IX)
      REL(IX+1,IY) = REL(IY,IX+1)
  111 CONTINUE
      IY1=0
  110 CONTINUE
      IX1=1
      IY1=0
      DO   112 IX=2+NBDL,LRE ,NBDL
      IX1=IX1 + 1
      DO   113 IY=5,IX  ,NBDL
      IY1=IY1 + 1
      REL(IY,IX) = REL(IY,IX) + VKL25*DJAC2*SHP(1,IX1)*SHP(4,IY1)
      REL(IY+1,IX) = REL(IY+1,IX) + VKL25*DJAC2*SHP(1,IX1)*SHP(5,IY1)
      REL(IX,IY) = REL(IY,IX)
      REL(IX,IY+1) = REL(IY+1,IX)
  113 CONTINUE
      IY1=0
  112 CONTINUE
      ENDIF
C
C     TERMES EN (UR,RT ) * (UR , RT )
C
      IX1=0
      IY1=0
      DO   114 IX=3,LRE ,NBDL
      IX1=IX1 + 1
      DO   115 IY=3,IX  ,NBDL
      IY1=IY1 + 1
      REL(IY,IX) = REL(IY,IX) + VKL33*DJAC1*SHP(4,IY1)*SHP(4,IX1)
      REL(IY,IX+1) = REL(IY,IX+1) + VKL33*DJAC1*SHP(4,IY1)*SHP(5,IX1)
      REL(IY+1,IX) = REL(IY+1,IX) + VKL33*DJAC1*SHP(5,IY1)*SHP(4,IX1)
      REL(IY+1,IX+1) = REL(IY+1,IX+1)+VKL33*DJAC1*SHP(5,IY1)*SHP(5,IX1)
      REL(IX,IY) = REL(IY,IX)
      REL(IX+1,IY) = REL(IY,IX+1)
      REL(IX,IY+1) = REL(IY+1,IX)
      REL(IX+1,IY+1) = REL(IY+1,IX+1)
  115 CONTINUE
      IY1=0
  114 CONTINUE
      IF ( IFOU.EQ.1) THEN
C
C     TERMES EN (UT,RR ) * (UT , RR )
C     NON NULS QU EN FOURIER
C
      DO   116 IX=3,LRE ,NBDL
      DO   117 IY=3,LRE ,NBDL
      IX2=IX + 2
      IY2=IY + 2
      REL(IX2,IY2)     = REL(IX,IY)
      REL(IX2+1,IY2)   = REL(IX+1,IY)
      REL(IX2,IY2+1)   = REL(IX,IY+1)
      REL(IX2+1,IY2+1) = REL(IX+1,IY+1)
  117 CONTINUE
  116 CONTINUE
      ENDIF
      IRET = 1
      GOTO  666
C
C     MESSAGE D ERREUR : ELEMENT A  SURFACE  NULLE
C
 667  CONTINUE
      IRET = 2
      GOTO  666
C
 666  CONTINUE
      RETURN
      END