Sources Esope | Cast3M

trihm3
C TRIHM3    SOURCE    CHAT      05/01/13    03:47:17     5004
      SUBROUTINE TRIHM3(IGAU,ITEL,NBNO,XEL,SHPTOT,SHP,RHOS,RHOF,B11,B22,
     #B12,F11,F12,SFLU,SCEL,POIGAU,VKL12,VKL22,VKL23,VKL33,LRE,
     #REL,ISDJC,IRET)
C=======================================================================
C
C    CALCULE LA MATRICE  DE  MASSE     DANS  LE  CAS  PLAN  POUR
C        LA FORMULATION (37) HOMOGENE
C=======================================================================
C     NBDL = NOMBRE DE DDL PAR NOEUD
C  INPUT
C     IGAU=NUMERO DU POINT DE GAUSS
C     ITEL=NUMERO DE L ELEMENT DANS NOMTP
C     NBNO=NOMBRE DE NOEUDS
C     XEL =COORDONNEES  DE L ELEMENT
C     RHOS = MASSE DU  SOLIDE
C     RHOF = MASSE VOLUMIQUE DU  FLUIDE
C     B11,B22,B12 =  PERMEABILITE ACOUSTIQUE DU MILIEU
C     F11,F12  = NORME
C     SFLU = SURFACE  FLUIDE  DANS LA CELLULE ELEMENTAIRE
C     SCEL = SURFACE  DE LA CELLULE ELEMENTAIRE
C     POIGAU=MINTE1.POIGAU(IGAU)
C     VKL12=-((COEFPR*COEFPI)/(RHOF*C**2))*SFLU/SCEL
C     VKL22=-(COEFPI**2)/(RHOF*SCEL)
C     VKL23= COEFPI/SCEL
C     VKL33= 1/SCEL
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  DU SIGNE DU  JACOBIEN
C  ZONE DE TRAVAIL
C     SHP(6,NBNO)=TABLEAU DE TRAVAIL
C OUTPUT
C     REL=MATRICE DE RIGIDITE
C     ISDJC= INDICATEUR  DU SIGNE DU  JACOBIEN
C     IRET=INDICATEUR = 1 : SUCCES
C                     = 0 : ECHEC (ITEL INCOMPATIBLE AVEC LA FORMULATION
C                     = 2 : ECHEC ( JACOBIEN NUL )
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
C
C     SHP(1,I) : FONCTION DE FORME
C     SHP(2,I) : DERIVEE % X DE LA FONCTION DE FORME
C     SHP(3,I) : DERIVEE % Y 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
      IDIM=2
      CALL JACOBI(XEL,SHP,IDIM,NBNO,DJAC)
      IF (DJAC.EQ.0.) GOTO 667
      IRET = 1
      IF (DJAC.LT.0.) ISDJC = ISDJC + 1
      DJAC=ABS(DJAC)*POIGAU
C
C     TERMES EN P*PI
C     NBDL : NOMBRE DE D.D.L PAR NOEUD  ( = 4 )
C
      NBDL = 4
      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) + DJAC*VKL12*SHP(1,IX1)*SHP(1,IY1)
      REL(IX,IY) = REL(IY,IX)
  103 CONTINUE
      IY1=0
  102 CONTINUE
      DO   104 IX=2,LRE-NBDL,NBDL
      IX1=IX +NBDL-1
      DO   105 IY=IX1,LRE ,NBDL
      REL(IX,IY) = REL(IX-1,IY+1)
      REL(IY,IX) = REL(IX,IY)
  105 CONTINUE
  104 CONTINUE
C
C     TERMES EN PI*PI
C
      IX1=0
      IY1=0
      DO   106 IX=2,LRE ,NBDL
      IX1=IX1 + 1
      DO   107 IY=2,IX  ,NBDL
      IY1=IY1 + 1
      REL(IY,IX) = REL(IY,IX) + DJAC*VKL22*((B11*SHP(2,IX1)*SHP(2,IY1)+
     #B22*SHP(3,IX1)*SHP(3,IY1)+B12*(SHP(2,IX1)*SHP(3,IY1)+SHP(3,IX1)*
     #SHP(2,IY1))))
      REL(IX,IY) = REL(IY,IX)
  107 CONTINUE
      IY1=0
  106 CONTINUE
C
C     TERMES EN PI*(UX,UY)
C
      D11=(SCEL-B11)*F12
      D22=(SCEL-B22)*F12
      D12=-B12*F12
      IX1=0
      IY1=0
      DO   108 IX=3,LRE ,NBDL
      IX1=IX1 + 1
      DO   109 IY=2,IX  ,NBDL
      IY1=IY1 + 1
      REL(IY,IX) = REL(IY,IX) + DJAC*VKL23 * ( D11*SHP(2,IY1)  +
     #D12*SHP(3,IY1) ) * SHP(1,IX1)
      REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL23 * ( D12*SHP(2,IY1) +
     #D22*SHP(3,IY1) ) * SHP(1,IX1)
      REL(IX,IY) = REL(IY,IX)
      REL(IX+1,IY) = REL(IY,IX+1)
  109 CONTINUE
      IY1=0
  108 CONTINUE
      IX1=1
      IY1=0
      DO   110 IX=2+NBDL,LRE ,NBDL
      IX1=IX1 + 1
      DO   111 IY=3,IX  ,NBDL
      IY1=IY1 + 1
      REL(IY,IX) = REL(IY,IX) + DJAC*VKL23 * ( D11*SHP(2,IX1) +
     #D12*SHP(3,IX1) ) * SHP(1,IY1)
      REL(IY+1,IX) = REL(IY+1,IX) + DJAC*VKL23 * (D12*SHP(2,IX1)+
     #D22*SHP(3,IX1) ) * SHP(1,IY1)
      REL(IX,IY) = REL(IY,IX)
      REL(IX,IY+1) = REL(IY+1,IX)
  111 CONTINUE
      IY1=0
  110 CONTINUE
C
C     TERMES EN (UX,UY)*(UX,UY)
C
      H11=RHOS + RHOF*(SFLU-B11)*F11
      H22=RHOS + RHOF*(SFLU-B22)*F11
      H12=  (-1.)*RHOF*B12*F11
      IX1=0
      IY1=0
      DO   112 IX=3,LRE ,NBDL
      IX1=IX1 + 1
      DO   113 IY=3,IX  ,NBDL
      IY1=IY1 + 1
      REL(IY,IX) = REL(IY,IX) + DJAC*VKL33*H11*SHP(1,IY1)*SHP(1,IX1)
      REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL33*H12*SHP(1,IY1)*SHP(1,IX1)
      REL(IY+1,IX+1)=REL(IY+1,IX+1)+DJAC*VKL33*H22*SHP(1,IY1)*SHP(1,IX1)
      REL(IX,IY) = REL(IY,IX)
      REL(IY+1,IX) = REL(IY,IX+1)
      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)
  113 CONTINUE
      IY1=0
  112 CONTINUE
      GOTO 666
C
C     MESSAGE D ERREUR : ELEMENT A SURFACE NULLE
C
  667 CONTINUE
      IRET = 2
      GOTO 666
C
  666 CONTINUE
      RETURN
      END