C CUBHM1 SOURCE CHAT 05/01/12 22:32:31 5004 C CUBHM1 SOURCE LOFT 90/01/26 21:08:44 SUBROUTINE CUBHM1(IGAU,ITEL,NBNO,XEL,SHPTOT,SHP,RHOS,RHOF,B11,B22, #B12,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 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.127) 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 SHP(4,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) SHP(4,NP)=SHPTOT(4,NP,IGAU) 101 CONTINUE IDIM=3 CALL JACOBI(XEL,SHP,IDIM,NBNO,DJAC) C C A1 = XEL(3,2) - XEL(3,1) A1 = ABS(A1) A2 = XEL(3,3) - XEL(3,1) A2 = ABS(A2) A3 = XEL(3,4) - XEL(3,1) A3 = ABS(A3) A4 = XEL(3,5) - XEL(3,1) A4 = ABS(A4) A5 = XEL(3,6) - XEL(3,1) A5 = ABS(A5) A6 = XEL(3,7) - XEL(3,1) A6 = ABS(A6) A7 = XEL(3,8) - XEL(3,1) A7 = ABS(A7) C A3 = (0.25D0)*(A1+A2+A3+A4+A5+A6+A7) C IF (DJAC.EQ.0.) GOTO 667 IRET = 1 IF (DJAC.LT.0.) ISDJC = ISDJC + 1 DJAC=ABS(DJAC)*POIGAU C WRITE(6,*)'DJAC=',DJAC C C C FONCTIONS DE FORME EN Z C ZH=SHP(1,5)+SHP(1,6)+SHP(1,7)+SHP(1,8) C Z1=1.D0-3.D0*ZH*ZH+2.D0*ZH*ZH*ZH Z2=3.D0*(ZH**2)-2.D0*ZH*ZH*ZH Z3=(A3*ZH)*(1.D0-2.D0*ZH+ZH*ZH) Z4=(A3*ZH)*(ZH*ZH-ZH) C C S1=SHP(1,1)+SHP(1,4)+SHP(1,5)+SHP(1,8) S2=SHP(1,2)+SHP(1,3)+SHP(1,6)+SHP(1,7) C T1=SHP(1,1)+SHP(1,2)+SHP(1,5)+SHP(1,6) T2=SHP(1,3)+SHP(1,4)+SHP(1,7)+SHP(1,8) C C C FONCTIONS DE FORME POUR LES DEPLACEMENTS C SHP(5,1)=S1*T1*Z1 SHP(5,2)=S2*T1*Z1 SHP(5,3)=S2*T2*Z1 SHP(5,4)=S1*T2*Z1 SHP(5,5)=S1*T1*Z2 SHP(5,6)=S2*T1*Z2 SHP(5,7)=S2*T2*Z2 SHP(5,8)=S1*T2*Z2 SHP(5,9)=S1*T1*Z3 SHP(5,10)=S2*T1*Z3 SHP(5,11)=S2*T2*Z3 SHP(5,12)=S1*T2*Z3 SHP(5,13)=S1*T1*Z4 SHP(5,14)=S2*T1*Z4 SHP(5,15)=S2*T2*Z4 SHP(5,16)=S1*T2*Z4 C C TERMES EN P*PI C NBDL : NOMBRE DE D.D.L PAR NOEUD ( = 6 ) C NBDL = 6 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 B33=SFLU IX1=0 IY1=0 DO 106 IX=2,LRE ,NBDL C WRITE(6,*)'IX= ',IX IX1=IX1 + 1 C WRITE(6,*)'IX1= ',IX1 DO 107 IY=2,IX ,NBDL C WRITE(6,*)'IY= ',IY IY1=IY1 + 1 C WRITE(6,*)'IY1= ',IY1 REL(IY,IX) = REL(IY,IX) + DJAC*VKL22*(B11*SHP(2,IX1)*SHP(2,IY1) #+B22*SHP(3,IX1)*SHP(3,IY1)+B33*SHP(4,IY1)*SHP(4,IX1)) REL(IX,IY) = REL(IY,IX) 107 CONTINUE IY1=0 106 CONTINUE C C TERMES EN PI*(UX,RY) C D11 = (SCEL - B11) D22 = (SCEL - B22) D33=SFLU C 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)*SHP(5,IX1)) REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL23 * (D11*SHP(2,IY1) #*SHP(5,IX1+8)) 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)*SHP(5,IY1) REL(IY+1,IX) = REL(IY+1,IX) + DJAC*VKL23* #D11*SHP(2,IX1)*SHP(5,IY1+8) REL(IX,IY) = REL(IY,IX) REL(IX,IY+1) = REL(IY+1,IX) 111 CONTINUE IY1=0 110 CONTINUE C C TERMES EN PI*(UY,RX) C IX1=0 IY1=0 DO 114 IX=5,LRE ,NBDL IX1=IX1 + 1 DO 115 IY=2,IX ,NBDL IY1=IY1 + 1 REL(IY,IX) = REL(IY,IX) + DJAC*VKL23 * D22*SHP(3,IY1) #* SHP(5,IX1) REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL23 * #D22*SHP(3,IY1) * SHP(5,IX1+8)*(-1.D0) REL(IX,IY) = REL(IY,IX) REL(IX+1,IY) = REL(IY,IX+1) 115 CONTINUE IY1=0 114 CONTINUE IX1=1 IY1=0 DO 116 IX=2+NBDL,LRE ,NBDL IX1=IX1 + 1 DO 117 IY=5,IX ,NBDL IY1=IY1 + 1 REL(IY,IX) = REL(IY,IX) + DJAC*VKL23 * D22*SHP(3,IX1) #* SHP(5,IY1) REL(IY+1,IX) = REL(IY+1,IX) + DJAC*VKL23 * #D22*SHP(3,IX1) * SHP(5,IY1+8)*(-1.D0) REL(IX,IY) = REL(IY,IX) REL(IX,IY+1) = REL(IY+1,IX) 117 CONTINUE IY1=0 116 CONTINUE C C TERMES EN (UX,RY)*(UX,RY) C H11=RHOS + RHOF*(SFLU-B11) H22=RHOS + RHOF*(SFLU-B22) 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(5,IY1)*SHP(5,IX1) REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL33*H11*SHP(5,IY1)*SHP(5,IX1+8) REL(IY+1,IX) = REL(IY+1,IX)+DJAC*VKL33*H11*SHP(5,IY1+8)*SHP(5,IX1) REL(IY+1,IX+1)=REL(IY+1,IX+1)+DJAC*VKL33*H22*SHP(5,IY1+8) #*SHP(5,IX1+8) 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) 113 CONTINUE IY1=0 112 CONTINUE C C TERMES EN (UY,RX)*(UY,RX) C H11=RHOS + RHOF*(SFLU-B11) H22=RHOS + RHOF*(SFLU-B22) IX1=0 IY1=0 DO 118 IX=5,LRE ,NBDL IX1=IX1 + 1 DO 119 IY=5,IX ,NBDL IY1=IY1 + 1 REL(IY,IX) = REL(IY,IX) + DJAC*VKL33*H11*SHP(5,IY1)*SHP(5,IX1) REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL33*H11*SHP(5,IY1)*SHP(5,IX1+8) #*(-1.D0) REL(IY+1,IX) = REL(IY+1,IX)-DJAC*VKL33*H11*SHP(5,IY1+8)*SHP(5,IX1) REL(IY+1,IX+1)=REL(IY+1,IX+1)+DJAC*VKL33*H22*SHP(5,IY1+8) #*SHP(5,IX1+8) 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) 119 CONTINUE IY1=0 118 CONTINUE GOTO 666 C C MESSAGE D ERREUR : ELEMENT A SURFACE NULLE C 667 CONTINUE IRET = 2 GOTO 666 C 666 CONTINUE RETURN END