C MASSE4 SOURCE PV090527 26/04/30 21:15:53 12529 *---------------------------------------------------------------------* * ________________________________ * * | | * * | calcul de la matrice de masse | * * |________________________________| * * * * raccords liquide/massifs,raccords liquide/coque,barre,homogenise * * cerce,joints 2-3d * * * *---------------------------------------------------------------------* * * * entrees : * * ________ * * * * ipmail pointeur sur un segment meleme * * lw dimension du tableau de travail de l'element * * lre nombre de ddl dans la matrice de masse * * ivamat pointeur sur un segment mptval pour le materiau * * nmatt nombre de composante de materiau (imat=1) * * ivacar pointeur sur un segment mptval pour les caracteri- * * stiques * * ncarr nombre de caracteristiques geometriques * * nbpgau nombre de point d'integration pour la masse * * ipmint pointeur sur un segment minte * * ipmin1 pointeur sur un segment minte (aux noeuds) * * nddl nombre de degre de liberte /noeud * * mele numero de l'element fini * * nbpgmi nombre de noeuds /element * * ilump =1 si l'opeateur LUMP est appele * * * * sorties : * * ________ * * * * ipmatr pointeur sur la matrice de masse de la sous-zone * * * *---------------------------------------------------------------------* SUBROUTINE MASSE4(IPMAIL,LW,LRE,IVAMAT,NMATT,IVACAR,NCARR, & NBPGAU,IPMINT,NDDL,MELE,MFR,IPMATR,ILUMP, & ISOUS,IIPDPG,IMOD) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC CCHAMP -INC CCREEL -INC SMRIGID -INC SMCHAML -INC SMELEME -INC SMCOORD -INC SMINTE -INC SMMODEL -INC SMLMOTS -INC TMPTVAL SEGMENT WRK1 REAL*8 REL(LRE,LRE),XE(3,NBBB) ENDSEGMENT SEGMENT WRK2 REAL*8 SHPWRK(6,NBNO),BGENE(NDDL,LRE) ENDSEGMENT SEGMENT WRK3 REAL*8 WORK(LW) ENDSEGMENT SEGMENT WRK4 REAL*8 BPSS(3,3),XEL(3,NBBB) ENDSEGMENT SEGMENT WRK5 REAL*8 XGENE(NCOM,LRN) ENDSEGMENT SEGMENT MVELCH REAL*8 VALMAT(NV1) ENDSEGMENT CHARACTER*8 CMATE CHARACTER*4 lesinc(7),lesdua(7) DATA lesinc/'UX','UY','UZ','RX','RY','RZ','UR'/ DATA lesdua/'FX','FY','FZ','MX','MY','MZ','FR'/ INTEGER KERRE MELEME=IPMAIL NBNN=NUM(/1) NBELEM=NUM(/2) xMATRI=IPMATR NV1=NMATT SEGINI,MVELCH KERRE=0 I195=0 I259=0 WRK1 = 0 WRK2 = 0 WRK3 = 0 WRK4 = 0 WRK5 = 0 * * introduction du point autour duquel se fait le mouvement * de la section en defo plane generalisee * IF (IFOUR.EQ.-3.AND.MFR.NE.35) THEN IREF=(IIPDPG-1)*(IDIM+1) XDPGE=XCOOR(IREF+1) YDPGE=XCOOR(IREF+2) ELSE XDPGE=0.D0 YDPGE=0.D0 ENDIF NHRM=NIFOUR MINTE=IPMINT IMODEL = IMOD CMATE = imodel.CMATEE jmat = 0 iinc = 0 idua = 0 DO imat = 1 , matmod(/2) if (matmod(imat).eq.'IMPEDANCE') then jmat = imat goto 45 endif ENDDO IF (mfr.eq.28) THEN jgn = 8 if (ifour.eq.2) then jgm = 6 segini mlmots iinc = mlmots do igm = 1,jgm mots(igm) = lesinc(igm) enddo segini mlmots idua = mlmots do igm= 1,jgm mots(igm) = lesdua(igm) enddo else if (ifour.lt.0) then jgm = 4 segini mlmots iinc = mlmots mots(1) = lesinc(1) mots(2) = lesinc(2) mots(3) = lesinc(4) mots(4) = lesinc(5) segini mlmots idua = mlmots mots(1) = lesdua(1) mots(2) = lesdua(2) mots(3) = lesdua(4) mots(4) = lesdua(5) else if (ifour.eq.0) then jgm = 3 segini mlmots iinc = mlmots mots(1) = lesinc(7) mots(2) = lesinc(3) mots(3) = lesinc(6) segini mlmots idua = mlmots mots(1) = lesdua(7) mots(2) = lesdua(3) mots(3) = lesdua(6) else if (ifour.eq.1) then * a faire endif ENDIF c numero des etiquettes : c etiquettes de 1 a 98 pour traitement specifique a l element c dans la zone specifique a chaque element commencant par : c 5 continue c element 5 etiquettes 1005 2005 3005 4005 ... c 44 continue c element 44 etiquettes 1044 2044 3044 4044 ... c_______________________________________________________________________ * CABL SEG2 SEG3 TRI3 TRI4 TRI6 TRI7 QUA4 QUA5 QUA8 QUA9 GOTO ( 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 * RAC2 RAC3 CUB8 CU20 PRI6 PR15 LIA3 LIA4 LIA6 LIA8 MULT & , 12, 99, 99, 99, 99, 99, 18, 18, 99, 99, 99 * TET4 TE10 PYR5 PY13 COQ3 DKT POUT LISP FAC3 FAC4 FAC6 & , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 * FAC8 LTR3 LQU4 LCU8 LPR6 LTE4 LPY5 COQ8 TUYA TUFI COQ2 & , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 * POI1 BARR RACO LSU2 COQ4 LISM COF3 RES2 LSU3 LSU4 LICO & , 45, 46, 47, 99, 99, 99, 99, 99, 99, 99, 55 * COQ6 CVS2 CVS3 CVT3 CVT6 CVQ4 CVQ8 THP5 TH13 THP6 TH15 & , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 * THC8 TH20 ICT3 ICQ4 ICT6 ICQ8 ICC8 ICT4 ICP6 IC20 IC10 & , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 * IC15 TRIP QUAP CUBP TETP PRIP TIMO JOI2 JOI3 JOT3 JOI4 & , 99, 99, 99, 99, 99, 99, 99, 85, 99, 87, 88 * JOI6 JOI8 LISC TRIH DST LIC4 CERC TUYO LSE2 LITU HYT3 & , 99, 99, 99, 92, 99, 94, 46, 99, 99, 12, 99 * HYQ4 HYT4 HYP6 HYC8 TRIS QUAS POIS FOR3 JOP3 JOP6 JOP8 & , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 * POL3 POL4 POL5 POL6 POL7 POL8 POL9 PO10 PO11 PO12 PO13 & , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 * PO14 BAR3 BAEX LIA2 QUAH CUBH ROT3 SEF2 TRF3 QUF4 CUF8 & , 99, 46, 124, 99, 126, 127, 99, 99, 99, 99, 99 * PRF6 TEF4 PYF5 MSE3 MTR6 MQU9 MC27 MP18 MT10 MP14 SEF3 & , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 * TRF7 QUF9 CF27 PF21 TF15 PF19 SEG6 TR21 QU36 C216 P126 & , 99, 99, 99, 99, 99, 99, 99, 99, 99, 99, 99 * TE56 PY91 TRH6 & , 99, 99, 157),MELE * GOTO(168,169,170,171,172),MELE-167 * * JOI1 GOTO(265),MELE-264 C--- CAS NON PREVUS ICI------------------------------------------------- 99 CONTINUE MOTERR(1:4)=NOMTP(MELE) MOTERR(5:12)='MASSE4' CALL ERREUR(86) GOTO 510 c_______________________________________________________________________ c secteur de calcul pour les elements de raccord rac2 rac3 litu c liquide massif lineaire cas bidimensionnel c_______________________________________________________________________ 12 CONTINUE IF (ILUMP .EQ. 1) GOTO 99 NBBB=NBNN IF (MELE.NE.98) LW=IDIM SEGINI WRK1,WRK3 DO 3012 IB=1,NBELEM c on cherche les coordonnees de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO(REL,LRE,LRE) c calcul des coefficients de normalisation MPTVAL=IVAMAT IF (MELE.NE.98) THEN DO 5012 IM=1,NMATT MELVAL=IVAL(IM) IBMN=MIN(IB,VELCHE(/2)) VALMAT(IM)=VELCHE(1,IBMN) 5012 CONTINUE RHOREF=VALMAT(1) RLCAR = VALMAT(2) CFPI= RHOREF*RLCAR ELSE c cas de l'element litu DO 7012 IM=1,NMATT MELVAL=IVAL(IM) IBMN=MIN(IB,VELCHE(/2)) WORK(IM+9)=VELCHE(1,IBMN) 7012 CONTINUE ENDIF c lecture des caracteristiques dans work MPTVAL=IVACAR DO 4012 IC=1,NCARR IF (IVAL(IC).NE.0) THEN MELVAL=IVAL(IC) IBMN=MIN(IB,VELCHE(/2)) WORK(IC)=VELCHE(1,IBMN) ELSE WORK(IC)=0.D0 ENDIF 4012 CONTINUE IF (MELE.EQ.98) THEN CALL COUMAS(REL,LRE,WORK,XE,KERRE) ELSE CALL RACMAS(NBPGAU,IFOUR,NIFOUR,IDIM,NBNN,XE,CFPI,WORK, 1 POIGAU,SHPTOT,REL,LRE) ENDIF c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) 3012 CONTINUE GOTO 510 c_______________________________________________________________________ c secteur de calcul pour les elements de raccord lia3 lia4 c liquide massif lineaire cas tridimensionnel c_______________________________________________________________________ 18 CONTINUE IF (ILUMP .EQ. 1) GOTO 99 NBBB=NBNN LW=IDIM SEGINI WRK1,WRK3 DO 3018 IB=1,NBELEM c on cherche les coordonnees de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO(REL,LRE,LRE) c calcul des coefficients de normalisation MPTVAL=IVAMAT DO 5018 IM=1,NMATT MELVAL=IVAL(IM) IBMN=MIN(IB,VELCHE(/2)) VALMAT(IM)=VELCHE(1,IBMN) 5018 CONTINUE RHOREF=VALMAT(1) RLCAR = VALMAT(2) CFPI= RHOREF*RLCAR MPTVAL=IVACAR DO 4018 IC=1,NCARR IF (IVAL(IC).NE.0) THEN MELVAL=IVAL(IC) IBMN=MIN(IB,VELCHE(/2)) WORK(IC)=VELCHE(1,IBMN) ELSE WORK(IC)=0.D0 ENDIF 4018 CONTINUE CALL LIAMAS(NBPGAU,IDIM,NBNN,NDDL,XE,CFPI,WORK,POIGAU, 1 SHPTOT,REL,LRE,IER246) IF(IER246.NE.0) THEN CALL ERREUR(IER246) GOTO 510 ENDIF c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) 3018 CONTINUE GOTO 510 c_______________________________________________________________________ c impedance c_______________________________________________________________________ 45 CONTINUE IF (jmat.gt.0) THEN MPTVAL=IVAMAT MELVAL=IVAL(1) if (ival(/1).gt.1) then melva1 = ival(2) else melva1 = 0 endif jddl = LRE/NBPGAU DO IB = 1,NBELEM JDIAG = 0 XMASS = 0.D0 if (melval.gt.0) IBMN=MIN(IB,VELCHE(/2)) do IG = 1, NBPGAU if (melval.gt.0) then igmn = MIN(IG,VELCHE(/1)) XMASS=VELCHE(IGMN,IBMN) endif XINER = XMASS if (melva1.gt.0) then igmn = MIN(IG,melva1.VELCHE(/1)) XINER = melva1.VELCHE(IGMN,IBMN) endif do idl = 1,jddl JDIAG = JDIAG + 1 RE(JDIAG,JDIAG,ib) = XMASS if (idim.eq.3.and.idl.gt.3) RE(JDIAG,JDIAG,ib) = XINER if (idim.ne.3.and.idl.gt.2) RE(JDIAG,JDIAG,ib) = XINER enddo enddo ENDDO GOTO 510 ENDIF IF (MFR.EQ.26) THEN * MODAL (car goto 510 sous IMPEDANCE) MPTVAL=IVAMAT MELVAL=IVAL(2) DO IB = 1,NBELEM IBMN=MIN(IB,VELCHE(/2)) RE(1,1,ib) = VELCHE(1,IBMN) ENDDO * ELSE IF (MFR.EQ.28) THEN * STATIQUE (car goto 510 sous IMPEDANCE) MPTVAL=IVAMAT DO IB = 1,NBELEM MELVAL=IVAL(1) IBMN=MIN(IB,IELCHE(/2)) idepl=IELCHE(1,IBMN) MELVAL=IVAL(3) IBMN=MIN(IB,IELCHE(/2)) imade=IELCHE(1,IBMN) CALL XTY1(idepl,imade,iinc,idua,X1) re(1,1,ib) = x1 ENDDO ENDIF GOTO 510 c_______________________________________________________________________ c element point (poi1) en defos planes generalisees c_______________________________________________________________________ IF(MELE.EQ.45.AND.IFOUR.NE.-3) GOTO 99 NBBB=NBNN SEGINI WRK1,WRK3 c boucle de calcul pour les differents elements DO 3045 IB=1,NBELEM c on cherche les coordonnees de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) c on recherche rho et la section MPTVAL=IVAMAT MELVAL=IVAL(1) IBMN=MIN(IB,VELCHE(/2)) RR=VELCHE(1,IBMN) MPTVAL=IVACAR MELVAL=IVAL(1) IBMN=MIN(IB,VELCHE(/2)) RR=RR*VELCHE(1,IBMN) c on calcule la matrice de masse CALL PO1MAS(XE,XDPGE,YDPGE,RR,LRE,REL) CALL REMPMT(REL,LRE,RE(1,1,ib)) 3045 CONTINUE GO TO 510 c_______________________________________________________________________ c elements barre et cerce c_______________________________________________________________________ 46 CONTINUE * IF(MELE.EQ.95.AND.IFOUR.NE.0.AND.IFOUR.NE.1) GOTO 99 * NBBB=NBNN SEGINI WRK1,WRK3 IF(MELE.EQ.123) THEN NCOM=NBNN LRN =LRE SEGINI WRK5 ENDIF c boucle de calcul pour les differents elements DO 3046 IB=1,NBELEM c on cherche les coordonnees de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) MPTVAL=IVAMAT MELVAL=IVAL(1) IBMN=MIN(IB,VELCHE(/2)) RR=VELCHE(1,IBMN) MPTVAL=IVACAR MELVAL=IVAL(1) IBMN=MIN(IB,VELCHE(/2)) RR=RR*VELCHE(1,IBMN) c on calcule la matrice de masse IF(MELE.EQ.46) CALL BARMAS(REL,LRE,RR,XE) IF(MELE.EQ.95) CALL CERMAS(REL,LRE,RR,XE) IF(MELE.EQ.123) CALL MASBA3(REL,LRE,RR,XE,XGENE,KERRE) IF (KERRE.NE.0) THEN INTERR(1)=ISOUS INTERR(2)=IB IF(MELE.EQ.123.AND.KERRE.EQ.1) CALL ERREUR(128) GOTO 510 ENDIF IF (ILUMP .EQ. 1) THEN CALL LUMP1(REL,LRE,RE(1,1,ib),IFOUR) ELSE CALL REMPMT(REL,LRE,RE(1,1,ib)) ENDIF 3046 CONTINUE GO TO 510 c_______________________________________________________________________ c JOINT UNIDIMENSIONNEL JOI1 c_______________________________________________________________________ 265 CONTINUE NBBB=NBNN SEGINI WRK1,WRK3,WRK4 IAW1=101 IAW2=IAW1+LRE*LRE IAW3=IAW2+LRE*LRE IAW4=IAW3+LRE*LRE MPTVAL=IVAMAT DO 3265 IB=1,NBELEM DO IC=1,NMATT IF (IVAL(IC).NE.0) THEN MELVAL=IVAL(IC) IBMN=MIN(IB,VELCHE(/2)) WORK(IC)=VELCHE(1,IBMN) ELSE WORK(IC)=0.D0 ENDIF END DO CALL MAPALU(NMATT,WORK,BPSS,IDIM) c on calcule la matrice de masse localement CALL JOIMAS(REL,LRE,WORK,NMATT,IDIM) c on passe en repère global CALL JOIGLO(REL,BPSS,WORK(IAW1),WORK(IAW2), & WORK(IAW3),WORK(IAW4),LRE,IDIM) IF (ILUMP .EQ. 1) THEN CALL LUMP1(REL,LRE,RE(1,1,ib),IFOUR) ELSE CALL REMPMT(REL,LRE,RE(1,1,ib)) ENDIF 3265 CONTINUE GO TO 510 c_______________________________________________________________________ c element barre 3d excentre (baex) c_______________________________________________________________________ 124 CONTINUE NBBB=NBNN NCOM=2 LRN =LRE SEGINI WRK1,WRK3,WRK5 c boucle de calcul pour les differents elements DO 3199 IB=1,NBELEM c on recupere la section de l'element, ses excentrements et son c orientation. les caracteristiques sont rangees dans work c selon l'ordre suivant: sect excz excy vx vy vz MPTVAL=IVACAR DO IC=1,NCARR IF(IVAL(IC).NE.0) THEN MELVAL=IVAL(IC) IBMN=MIN(IB,VELCHE(/2)) WORK(IC)=VELCHE(1,IBMN) ELSE WORK(IC)=0.D0 ENDIF END DO SECT=WORK(1) c xgene stocke la matrice de passage de l'element excentre CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL MAPAEX(XE,NBNN,WORK,AL,XGENE,LRE,KERRE) IF (KERRE.NE.0) THEN INTERR(1)=ISOUS INTERR(2)=IB IF (KERRE.EQ.1) CALL ERREUR(128) ENDIF MPTVAL=IVAMAT MELVAL=IVAL(1) IBMN=MIN(IB,VELCHE(/2)) RR=VELCHE(1,IBMN)*SECT c on calcule la matrice de masse CALL BAMAEX(REL,LRE,RR,AL,XGENE) IF (ILUMP .EQ. 1) THEN CALL LUMP1(REL,LRE,RE(1,1,ib),IFOUR) ELSE CALL REMPMT(REL,LRE,RE(1,1,ib)) ENDIF 3199 CONTINUE GO TO 510 c_______________________________________________________________________ c secteur de calcul pour les elements de raccord c liquide coque cas bidimensionnel c_______________________________________________________________________ 47 CONTINUE IF (ILUMP .EQ. 1) GOTO 99 NBBB=NBNN LW=IDIM SEGINI WRK1,WRK3 DO 3047 IB=1,NBELEM c on cherche les coordonnees de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO(REL,LRE,LRE) c calcul des coefficients de normalisation MPTVAL=IVAMAT DO 5047 IM=1,NMATT MELVAL=IVAL(IM) IBMN=MIN(IB,VELCHE(/2)) VALMAT(IM)=VELCHE(1,IBMN) 5047 CONTINUE RHOREF=VALMAT(1) RLCAR = VALMAT(2) CFPI= RHOREF*RLCAR MPTVAL=IVACAR DO 4047 IC=1,NCARR IF (IVAL(IC).NE.0) THEN MELVAL=IVAL(IC) IBMN=MIN(IB,VELCHE(/2)) WORK(IC)=VELCHE(1,IBMN) ELSE WORK(IC)=0.D0 ENDIF 4047 CONTINUE CALL RACOMA(IFOUR,NIFOUR,IDIM,XE,CFPI,WORK,REL,LRE) c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) 3047 CONTINUE GOTO 510 c_______________________________________________________________________ c secteur de calcul pour les elements de raccord c liquide coque 3 noeuds - cas tridimensionnel c_______________________________________________________________________ 55 CONTINUE IF (ILUMP .EQ. 1) GOTO 99 NBBB=NBNN LW=IDIM SEGINI WRK1,WRK3 DO 3055 IB=1,NBELEM c on cherche les coordonnees de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO(REL,LRE,LRE) c calcul des coefficients de normalisation MPTVAL=IVAMAT DO 5055 IM=1,NMATT MELVAL=IVAL(IM) IBMN=MIN(IB,VELCHE(/2)) VALMAT(IM)=VELCHE(1,IBMN) 5055 CONTINUE RHOREF=VALMAT(1) RLCAR = VALMAT(2) CFPI= RHOREF*RLCAR MPTVAL=IVACAR DO 4055 IC=1,NCARR IF (IVAL(IC).NE.0) THEN MELVAL=IVAL(IC) IBMN=MIN(IB,VELCHE(/2)) WORK(IC)=VELCHE(1,IBMN) ELSE WORK(IC)=0.D0 ENDIF 4055 CONTINUE CALL LICOMA(NBPGAU,IDIM,NBNN,NDDL,XE,CFPI,WORK,POIGAU, 1 QSIGAU,ETAGAU,SHPTOT,REL,LRE,IER246) IF (IER246.NE.0) THEN CALL ERREUR(IER246) GOTO 510 ENDIF c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) 3055 CONTINUE GOTO 510 c_______________________________________________________________________ c secteur de calcul pour les elements joints joi2 c_______________________________________________________________________ 85 CONTINUE IF (ILUMP .EQ. 1) GOTO 99 NBNO=NBNN NBBB=NBNN SEGINI WRK1,WRK2,WRK4 DO 3085 IB=1,NBELEM c on cherche les coordonnees des noeuds de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO (REL,LRE,LRE) c calcul des coordonnees locales de l'element CALL JO2LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) c boucle sur les points de gauss ISDJC=0 DO 4085 IGAU=1,NBPGAU CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 1 1.D0,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) * IF(DJAC.LT.0.) ISDJC=ISDJC+1 IF(DJAC.EQ.0.) I259=IB DJAC=ABS(DJAC)*POIGAU(IGAU) MPTVAL=IVAMAT IF (IVAL(1).NE.0) THEN MELVAL=IVAL(1) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) VALMAT(1)=VELCHE(IGMN,IBMN) ELSE VALMAT(1)=0.D0 ENDIF CCCCCCCCCCC DJAC=DJAC*VALMAT(1)/3.0D0 C IL FAUT DIVISER PAR 4, CE QUI CORRESPOND PLUS EXACTEMENT A DIVISER C LE B PAR 2... DJAC=DJAC*VALMAT(1)/4.0D0 c cas axisymetrique : multiplication par le rayon de courbure IF (IFOUR.EQ.0) THEN RAYON = 0.D0 NUMSUP=NBNO/2 DO 4185 IRAY=1,NUMSUP RAYON=RAYON + ( SHPTOT(1,IRAY,IGAU)*XE(1,IRAY) ) 4185 CONTINUE DJAC=DJAC*RAYON ENDIF CALL NTNST(BGENE,DJAC,LRE,NDDL,REL) 4085 CONTINUE IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) 3085 CONTINUE GOTO 510 c_______________________________________________________________________ c secteur de calcul pour les elements joints jot3 c_______________________________________________________________________ 87 CONTINUE IF (ILUMP .EQ. 1) GOTO 99 NBNO=NBNN NBBB=NBNN SEGINI WRK1,WRK2,WRK4 DO 3087 IB=1,NBELEM c on cherche les coordonnees des noeuds de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO (REL,LRE,LRE) c calcul des coordonnees locales de l'element CALL JT3LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) c boucle sur les points de gauss ISDJC=0 DO 4087 IGAU=1,NBPGAU CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 1 1.D0,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) IF(DJAC.LT.0.) ISDJC=ISDJC+1 IF(DJAC.EQ.0.) I259=IB DJAC=ABS(DJAC)*POIGAU(IGAU) MPTVAL=IVAMAT IF (IVAL(1).NE.0) THEN MELVAL=IVAL(1) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) VALMAT(1)=VELCHE(IGMN,IBMN) ELSE VALMAT(1)=0.D0 ENDIF DJAC=DJAC*VALMAT(1) C IL FAUT DIVISER PAR 4, CE QUI CORRESPOND PLUS EXACTEMENT A DIVISER C LE B PAR 2... DJAC=DJAC/4.0D0 CALL NTNST(BGENE,DJAC,LRE,NDDL,REL) 4087 CONTINUE IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) 3087 CONTINUE GOTO 510 c_______________________________________________________________________ c secteur de calcul pour les elements joints joi4 c_______________________________________________________________________ 88 CONTINUE IF (ILUMP .EQ. 1) GOTO 99 NBNO=NBNN NBBB=NBNN SEGINI WRK1,WRK2,WRK4 DO 3088 IB=1,NBELEM c on cherche les coordonnees des noeuds de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO (REL,LRE,LRE) c calcul des coordonnees locales de l'element CALL JO4LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) c boucle sur les points de gauss ISDJC=0 DO 4088 IGAU=1,NBPGAU CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 1 1.D0,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) IF(DJAC.LT.0.) ISDJC=ISDJC+1 IF(DJAC.EQ.0.) I259=IB DJAC=ABS(DJAC)*POIGAU(IGAU) MPTVAL=IVAMAT IF (IVAL(1).NE.0) THEN MELVAL=IVAL(1) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) VALMAT(1)=VELCHE(IGMN,IBMN) ELSE VALMAT(1)=0.D0 ENDIF DJAC=DJAC*VALMAT(1) C IL FAUT DIVISER PAR 4, CE QUI CORRESPOND PLUS EXACTEMENT A DIVISER C LE B PAR 2... DJAC=DJAC/4.0D0 CALL NTNST(BGENE,DJAC,LRE,NDDL,REL) 4088 CONTINUE IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) 3088 CONTINUE GOTO 510 c_______________________________________________________________________ c secteur de calcul pour les elements joints jgi2 c_______________________________________________________________________ 170 CONTINUE IF (IFOUR.EQ.-3) NDDL=NDDL+1 IF (ILUMP .EQ. 1) GOTO 99 NBNO=NBNN NBBB=NBNN SEGINI WRK1,WRK2,WRK4 IG1=0 DO IB=1,NBELEM c on cherche les coordonnees des noeuds de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO (REL,LRE,LRE) c calcul des coordonnees locales de l'element CALL JO2LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) c boucle sur les points de gauss ISDJC=0 DO IGAU=1,NBPGAU MPTVAL=IVAMAT DO IM=1,NMATT MELVAL=IVAL(IM) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) VALMAT(IM)=VELCHE(IGMN,IBMN) ENDDO EPAIST=VALMAT(2) IF(EPAIST.EQ.0.D0)THEN CALL BJO2GN(IGAU,MFR,IFOUR,NIFOUR,XE,XEL,BPSS,SHPTOT, . SHPWRK,EPAIST,BGENE,DJAC,XDPGE,YDPGE,IERT) IF(IERT.NE.0) IG1=IB ENDIF CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 1 1.D0,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) * IF(DJAC.LT.0.) ISDJC=ISDJC+1 IF(DJAC.EQ.0.) I259=IB DJAC=ABS(DJAC)*POIGAU(IGAU) * c valmat(1)=rho, valmat(2)=epai c /4 correspnnd en fait a diviser les matrices B par 2 CCCCCCCCCCCC DJAC=DJAC*VALMAT(1)*VALMAT(2)/4.0D0 DJAC=DJAC*VALMAT(1)*EPAIST/4.0D0 c cas axisymetrique : multiplication par le rayon de courbure IF (IFOUR.EQ.0) THEN RAYON = 0.D0 NUMSUP=NBNO/2 DO IRAY=1,NUMSUP RAYON=RAYON + ( SHPTOT(1,IRAY,IGAU)*XE(1,IRAY) ) ENDDO DJAC=DJAC*RAYON ENDIF CALL NTNST(BGENE,DJAC,LRE,NDDL,REL) ENDDO IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) ENDDO IF (IG1.NE.0) THEN INTERR(1)=IG1 CALL ERREUR (611) ENDIF GOTO 510 c_______________________________________________________________________ c secteur de calcul pour les elements joints jct3 en 2D cisaillement c_______________________________________________________________________ 168 CONTINUE IF (ILUMP .EQ. 1) GOTO 99 NBNO=NBNN NBBB=NBNN SEGINI WRK1,WRK2,WRK4 DO IB=1,NBELEM c on cherche les coordonnees des noeuds de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO (REL,LRE,LRE) c calcul des coordonnees locales de l'element CALL JT3LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) c boucle sur les points de gauss ISDJC=0 DO IGAU=1,NBPGAU CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 1 1.D0,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) IF(DJAC.LT.0.) ISDJC=ISDJC+1 IF(DJAC.EQ.0.) I259=IB DJAC=ABS(DJAC)*POIGAU(IGAU) MPTVAL=IVAMAT IF (IVAL(1).NE.0) THEN MELVAL=IVAL(1) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) VALMAT(1)=VELCHE(IGMN,IBMN) ELSE VALMAT(1)=0.D0 ENDIF DJAC=DJAC*VALMAT(1) C Il faut diviser par 4, ce qui correspond plus exactement a diviser C le B par 2... DJAC=DJAC/4.0D0 CALL NTNST(BGENE,DJAC,LRE,NDDL,REL) ENDDO IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) ENDDO GOTO 510 c_______________________________________________________________________ c secteur de calcul pour les elements joints jgt3 generalise c_______________________________________________________________________ 171 CONTINUE IF (ILUMP .EQ. 1) GOTO 99 NBNO=NBNN NBBB=NBNN SEGINI WRK1,WRK2,WRK4 IG1=0 DO IB=1,NBELEM c on cherche les coordonnees des noeuds de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO (REL,LRE,LRE) c calcul des coordonnees locales de l'element CALL JT3LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) c boucle sur les points de gauss ISDJC=0 DO IGAU=1,NBPGAU MPTVAL=IVAMAT DO IM=1,NMATT MELVAL=IVAL(IM) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) VALMAT(IM)=VELCHE(IGMN,IBMN) ENDDO EPAIST=VALMAT(2) IF(EPAIST.EQ.0.D0)THEN CALL BJT3G(IGAU,MFR,IFOUR,NIFOUR,XE,XEL,BPSS,SHPTOT, . SHPWRK,EPAIST,BGENE,DJAC,IERT) IF(IERT.NE.0) IG1=IB ENDIF CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 1 1.D0,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) * IF(DJAC.LT.0.) ISDJC=ISDJC+1 IF(DJAC.EQ.0.) I259=IB DJAC=ABS(DJAC)*POIGAU(IGAU) * c valmat(1)=rho, valmat(2)=epai c /4 correspond en fait a diviser les matrices B par 2 DJAC=DJAC*VALMAT(1)*EPAIST/4.0D0 * CALL NTNST(BGENE,DJAC,LRE,NDDL,REL) ENDDO IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) ENDDO IF (IG1.NE.0) THEN INTERR(1)=IG1 CALL ERREUR (611) ENDIF GOTO 510 c_______________________________________________________________________ c secteur de calcul pour les elements joints jci4 en 2D cisaillement c_______________________________________________________________________ 169 CONTINUE IF (ILUMP .EQ. 1) GOTO 99 NBNO=NBNN NBBB=NBNN SEGINI WRK1,WRK2,WRK4 DO IB=1,NBELEM c on cherche les coordonnees des noeuds de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO (REL,LRE,LRE) c calcul des coordonnees locales de l'element CALL JO4LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) c boucle sur les points de gauss ISDJC=0 DO IGAU=1,NBPGAU CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 1 1.D0,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) IF(DJAC.LT.0.) ISDJC=ISDJC+1 IF(DJAC.EQ.0.) I259=IB DJAC=ABS(DJAC)*POIGAU(IGAU) MPTVAL=IVAMAT IF (IVAL(1).NE.0) THEN MELVAL=IVAL(1) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) VALMAT(1)=VELCHE(IGMN,IBMN) ELSE VALMAT(1)=0.D0 ENDIF DJAC=DJAC*VALMAT(1) CCCCCCCCCCC DJAC=DJAC/3.0D0 C Il faut diviser par 4, ce qui correspond plus exactement a diviser C le B par 2... DJAC=DJAC/4.0D0 CALL NTNST(BGENE,DJAC,LRE,NDDL,REL) ENDDO IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) ENDDO GOTO 510 c_______________________________________________________________________ c secteur de calcul pour les elements joints jgi4 generalise c_______________________________________________________________________ 172 CONTINUE IF (ILUMP .EQ. 1) GOTO 99 NBNO=NBNN NBBB=NBNN SEGINI WRK1,WRK2,WRK4 IG1=0 DO IB=1,NBELEM c on cherche les coordonnees des noeuds de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO (REL,LRE,LRE) c calcul des coordonnees locales de l'element CALL JO4LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) c boucle sur les points de gauss ISDJC=0 DO IGAU=1,NBPGAU MPTVAL=IVAMAT DO IM=1,NMATT MELVAL=IVAL(IM) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) VALMAT(IM)=VELCHE(IGMN,IBMN) ENDDO EPAIST=VALMAT(2) IF(EPAIST.EQ.0.D0)THEN CALL BJO4G(IGAU,XE,XEL,BPSS,SHPTOT,SHPWRK,EPAIST, . BGENE,DJAC,IERT) IF(IERT.NE.0) IG1=IB ENDIF CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 1 1.D0,XEL,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) * IF(DJAC.LT.0.) ISDJC=ISDJC+1 IF(DJAC.EQ.0.) I259=IB DJAC=ABS(DJAC)*POIGAU(IGAU) * c valmat(1)=rho, valmat(2)=epai c /4 correspnnd en fait a diviser les matrices B par 2 CCCCCCCCCCCC DJAC=DJAC*VALMAT(1)*VALMAT(2)/4.0D0 DJAC=DJAC*VALMAT(1)*EPAIST/4.0D0 * CALL NTNST(BGENE,DJAC,LRE,NDDL,REL) ENDDO IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) ENDDO IF (IG1.NE.0) THEN INTERR(1)=IG1 CALL ERREUR (611) ENDIF GOTO 510 c_______________________________________________________________________ c secteur de calcul pour les elements homogeneises c (liquide solide) trih c_______________________________________________________________________ 92 CONTINUE IF (ILUMP .EQ. 1) GOTO 99 NBNO=NBNN NBBB=NBNN SEGINI WRK1,WRK2 DO 3092 IB=1,NBELEM c on cherche les coordonnees des noeuds de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO (REL,LRE,LRE) c on cherche les caracteristiques du materiau de l element ib MPTVAL=IVAMAT DO 8092 IM=1,NMATT MELVAL=IVAL(IM) IBMN=MIN(IB,VELCHE(/2)) VALMAT(IM)=VELCHE(1,IBMN) 8092 CONTINUE B11 =VALMAT(1) B22 =VALMAT(2) B12 =VALMAT(3) RHOF =VALMAT(4) RHOS =VALMAT(5) C =VALMAT(6) RHOREF=VALMAT(7) CREF =VALMAT(8) RLCAR =VALMAT(9) c on cherche les caracteristiques geometriques de l element ib MPTVAL=IVACAR IF (IFOUR.EQ.0.OR.IFOUR.EQ.1) THEN MELVAL=IVAL(1) IBMN=MIN(IB,VELCHE(/2)) SECT=VELCHE(1,IBMN) MELVAL=IVAL(2) IBMN=MIN(IB,VELCHE(/2)) SCEL=VELCHE(1,IBMN) MELVAL=IVAL(3) IBMN=MIN(IB,VELCHE(/2)) SFLU=VELCHE(1,IBMN) MELVAL=IVAL(4) IBMN=MIN(IB,VELCHE(/2)) EPS =VELCHE(1,IBMN) ELSE SECT=1.D0 MELVAL=IVAL(1) IBMN=MIN(IB,VELCHE(/2)) SCEL=VELCHE(1,IBMN) MELVAL=IVAL(2) IBMN=MIN(IB,VELCHE(/2)) SFLU=VELCHE(1,IBMN) MELVAL=IVAL(3) IBMN=MIN(IB,VELCHE(/2)) EPS =VELCHE(1,IBMN) MELVAL=IVAL(4) IBMN=MIN(IB,VELCHE(/2)) F11 =VELCHE(1,IBMN) MELVAL=IVAL(5) IBMN=MIN(IB,VELCHE(/2)) F12 =VELCHE(1,IBMN) ENDIF c calcul de la masse m0/eps**2 RHOSS=RHOS*SECT/(EPS*EPS) c calcul des coefficients de normalisation COEFPR=(RHOREF*CREF*CREF)/RLCAR COEFPI=RHOREF*RLCAR VKL12 =-(COEFPR*COEFPI*SFLU)/(RHOF*C*C*SCEL) VKL22 =-(COEFPI*COEFPI)/(RHOF*SCEL) VKL23 = COEFPI/SCEL VKL33 = 1.D0/SCEL IF(IFOUR.EQ.0.OR.IFOUR.EQ.1) THEN VKL23 =COEFPI*0.5D0*(2.D0*SCEL-B11-B22)/SCEL VKL33 =(RHOSS+RHOF*(SFLU-(B11+B22)/2.D0))/SCEL c calcul des termes en pi*pi c integration par nbpgau points de gauss ISDJC=0 DO 4092 IGAU=1,NBPGAU POIGA2=MINTE.POIGAU(IGAU) CALL TRIHM1(IGAU,MELE,MFR,NBNO,XE,SHPTOT,SHPWRK, # IFOUR,NHRM,B11,B22,SFLU,POIGA2,VKL22,LRE,REL,IRRT) IF(IRRT.NE.1) GOTO 7092 c calcul du reste des termes de la matrice masse c integration par nbpgau points de gauss CALL TRIHM2(IGAU,MELE,MFR,NBNO,XE,MINTE.SHPTOT,SHPWRK # ,IFOUR,NHRM,VKL12,VKL23,VKL33,POIGA2,ISDJC,LRE,REL,IRRT) IF(IRRT.NE.1) GOTO 7092 4092 CONTINUE IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB ELSE c boucle sur les points de gauss c cas plan ISDJC=0 DO 6092 IGAU1=1,NBPGAU POIGA1=MINTE.POIGAU(IGAU1) CALL TRIHM3(IGAU1,MELE,NBNO,XE,SHPTOT,SHPWRK # ,RHOSS,RHOF, # B11,B22,B12,F11,F12,SFLU,SCEL,POIGA1,VKL12,VKL22, # VKL23,VKL33,LRE,REL,ISDJC,IRRT) IF(IRRT.NE.1) GOTO 7092 6092 CONTINUE IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB ENDIF c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) 3092 CONTINUE c impression d un eventuel message d erreur 7092 CONTINUE IF(IRRT.EQ.0) THEN MOTERR(1:4)=NOMTP(MELE) CALL ERREUR(420) ELSE IF(IRRT.EQ.2) THEN INTERR(1) = IB CALL ERREUR(405) ENDIF ENDIF IF(IRRT.EQ.3) CALL ERREUR(421) IF(IRRT.EQ.4) CALL ERREUR(422) GOTO 510 c_______________________________________________________________________ c secteur de calcul pour les elements de raccord c liquide coque 4 noeuds - cas tridimensionnel c_______________________________________________________________________ 94 CONTINUE IF (ILUMP .EQ. 1) GOTO 99 NBBB=NBNN LW=IDIM SEGINI WRK1,WRK3 DO 3094 IB=1,NBELEM c on cherche les coordonnees de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO(REL,LRE,LRE) c calcul des coefficients de normalisation MPTVAL=IVAMAT DO 5094 IM=1,NMATT MELVAL=IVAL(IM) IBMN=MIN(IB,VELCHE(/2)) VALMAT(IM)=VELCHE(1,IBMN) 5094 CONTINUE RHOREF=VALMAT(1) RLCAR = VALMAT(2) CFPI= RHOREF*RLCAR MPTVAL=IVACAR DO 4094 IC=1,NCARR IF (IVAL(IC).NE.0) THEN MELVAL=IVAL(IC) IBMN=MIN(IB,VELCHE(/2)) WORK(IC)=VELCHE(1,IBMN) ELSE WORK(IC)=0.D0 ENDIF 4094 CONTINUE CALL LIC4MA(NBPGAU,IDIM,NBNN,NDDL,XE,CFPI,WORK,POIGAU, 1 QSIGAU,ETAGAU,SHPTOT,REL,LRE,IER246) IF(IER246.NE.0) THEN CALL ERREUR(IER246) GOTO 510 ENDIF c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) 3094 CONTINUE GOTO 510 c_______________________________________________________________________ c secteur de calcul pour les elements homogeneises c (liquide solide) quah c_______________________________________________________________________ 126 CONTINUE IF (ILUMP .EQ. 1) GOTO 99 NBNO=NBNN NBBB=NBNN SEGINI WRK1,WRK2 DO 3126 IB=1,NBELEM c on cherche les coordonnees des noeuds de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO (REL,LRE,LRE) c on cherche les caracteristiques du materiau de l element ib MPTVAL=IVAMAT DO 8126 IM=1,NMATT MELVAL=IVAL(IM) IBMN=MIN(IB,VELCHE(/2)) VALMAT(IM)=VELCHE(1,IBMN) 8126 CONTINUE B11 =VALMAT(1) B22 =VALMAT(2) B12 =VALMAT(3) RHOF =VALMAT(4) RHOS =VALMAT(5) C =VALMAT(6) RHOREF=VALMAT(7) CREF =VALMAT(8) RLCAR =VALMAT(9) c on cherche les caracteristiques geometriques de l element ib MPTVAL=IVACAR MELVAL=IVAL(4) IBMN=MIN(IB,VELCHE(/2)) SECT=VELCHE(1,IBMN) MELVAL=IVAL(1) IBMN=MIN(IB,VELCHE(/2)) SCEL=VELCHE(1,IBMN) MELVAL=IVAL(2) IBMN=MIN(IB,VELCHE(/2)) SFLU=VELCHE(1,IBMN) MELVAL=IVAL(3) IBMN=MIN(IB,VELCHE(/2)) EPS =VELCHE(1,IBMN) c calcul de la masse m0/eps**2 RHOSS=RHOS*SECT/(EPS*EPS) c calcul des coefficients de normalisation COEFPR=(RHOREF*CREF*CREF)/RLCAR COEFPI=RHOREF*RLCAR VKL12 =-(COEFPR*COEFPI*SFLU)/(RHOF*C*C*SCEL) VKL22 =-(COEFPI*COEFPI)/(RHOF*SCEL) VKL23 =COEFPI*0.5D0*(2.D0*SCEL-B11-B22)/SCEL VKL33 =(RHOSS+RHOF*(SFLU-(B11+B22)/2.D0))/SCEL c calcul des termes en pi*pi c integration par nbpgau points de gauss ISDJC=0 DO 4126 IGAU=1,NBPGAU POIGA2=MINTE.POIGAU(IGAU) CALL QUAHM1(IGAU,MELE,MFR,NBNO,XE,SHPTOT,SHPWRK,IFOUR # ,NHRM,B11,B22,SFLU,POIGA2,VKL22,LRE,REL,IRRT) IF(IRRT.NE.1) GOTO 7126 c calcul du reste des termes de la matrice masse c integration par nbpgau points de gauss CALL QUAHM2(IGAU,MELE,MFR,NBNO,XE,MINTE.SHPTOT,SHPWRK # ,IFOUR,NHRM,VKL12,VKL23,VKL33,POIGA2,ISDJC,LRE,REL,IRRT) IF(IRRT.NE.1) GOTO 7126 4126 CONTINUE IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) 3126 CONTINUE c impression d un eventuel message d erreur 7126 CONTINUE IF(IRRT.EQ.0) THEN MOTERR(1:4)=NOMTP(MELE) CALL ERREUR(420) ELSE IF(IRRT.EQ.2) THEN INTERR(1) = IB CALL ERREUR(405) ENDIF ENDIF IF(IRRT.EQ.3) CALL ERREUR(421) IF(IRRT.EQ.4) CALL ERREUR(422) GOTO 510 c_______________________________________________________________________ c secteur de calcul pour les elements homogeneises c (liquide solide) cubh c_______________________________________________________________________ 127 CONTINUE IF (ILUMP .EQ. 1) GOTO 99 NBNO=NBNN NBBB=NBNN LW=IDIM SEGINI WRK1,WRK2 DO 3127 IB=1,NBELEM c on cherche les coordonnees des noeuds de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO (REL,LRE,LRE) c on cherche les caracteristiques du materiau de l element ib MPTVAL=IVAMAT DO 8127 IM=1,NMATT MELVAL=IVAL(IM) IBMN=MIN(IB,VELCHE(/2)) VALMAT(IM)=VELCHE(1,IBMN) 8127 CONTINUE B11 =VALMAT(1) B22 =VALMAT(2) B12 =VALMAT(3) RHOF =VALMAT(4) RHOS =VALMAT(5) C =VALMAT(6) RHOREF=VALMAT(7) CREF =VALMAT(8) RLCAR =VALMAT(9) c on cherche les caracteristiques geometriques de l element ib MPTVAL=IVACAR MELVAL=IVAL(1) IBMN=MIN(IB,VELCHE(/2)) SCEL=VELCHE(1,IBMN) MELVAL=IVAL(2) IBMN=MIN(IB,VELCHE(/2)) SFLU=VELCHE(1,IBMN) MELVAL=IVAL(3) IBMN=MIN(IB,VELCHE(/2)) EPS =VELCHE(1,IBMN) MELVAL=IVAL(4) IBMN=MIN(IB,VELCHE(/2)) SECT =VELCHE(1,IBMN) c calcul de la masse m0/eps**2 RHOSS=RHOS*SECT/(EPS*EPS) c calcul des coefficients de normalisation COEFPR=(RHOREF*CREF*CREF)/RLCAR COEFPI=RHOREF*RLCAR VKL12 =-(COEFPR*COEFPI*SFLU)/(RHOF*C*C*SCEL) VKL22 =-(COEFPI*COEFPI)/(RHOF*SCEL) VKL23 = COEFPI/SCEL VKL33 = 1.D0/SCEL ISDJC=0 DO 6127 IGAU1=1,NBPGAU POIGA1=MINTE.POIGAU(IGAU1) CALL CUBHM1(IGAU1,MELE,NBNO,XE,SHPTOT,SHPWRK, # RHOSS,RHOF,B11,B22,B12,SFLU,SCEL,POIGA1,VKL12,VKL22,VKL23, # VKL33,LRE,REL,ISDJC,IRRT) IF(IRRT.NE.1) GOTO 7127 6127 CONTINUE IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) 3127 CONTINUE c impression d un eventuel message d erreur 7127 CONTINUE IF(IRRT.EQ.0) THEN MOTERR(1:4)=NOMTP(MELE) CALL ERREUR(420) ELSE IF(IRRT.EQ.2) THEN INTERR(1) = IB CALL ERREUR(405) ENDIF ENDIF IF(IRRT.EQ.3) CALL ERREUR(421) IF(IRRT.EQ.4) CALL ERREUR(422) GOTO 510 c_______________________________________________________________________ c secteur de calcul pour les elements homogeneises c (liquide solide) trh6 c_______________________________________________________________________ 157 CONTINUE IF (ILUMP .EQ. 1) GOTO 99 NBNO=NBNN NBBB=NBNN SEGINI WRK1,WRK2 DO 3157 IB=1,NBELEM c on cherche les coordonnees des noeuds de l element ib CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO (REL,LRE,LRE) c on cherche les caracteristiques du materiau de l element ib MPTVAL=IVAMAT DO 8157 IM=1,NMATT MELVAL=IVAL(IM) IBMN=MIN(IB,VELCHE(/2)) VALMAT(IM)=VELCHE(1,IBMN) 8157 CONTINUE B11 =VALMAT(1) B22 =VALMAT(2) B12 =VALMAT(3) RHOF =VALMAT(4) RHOS =VALMAT(5) C =VALMAT(6) RHOREF=VALMAT(7) CREF =VALMAT(8) RLCAR =VALMAT(9) E111 =VALMAT(10) E112 =VALMAT(11) E121 =VALMAT(12) E122 =VALMAT(13) E221 =VALMAT(14) E222 =VALMAT(15) c on cherche les caracteristiques geometriques de l element ib MPTVAL=IVACAR IF (IFOUR.EQ.0.OR.IFOUR.EQ.1) THEN MELVAL=IVAL(1) IBMN=MIN(IB,VELCHE(/2)) SECT=VELCHE(1,IBMN) MELVAL=IVAL(2) IBMN=MIN(IB,VELCHE(/2)) SCEL=VELCHE(1,IBMN) MELVAL=IVAL(3) IBMN=MIN(IB,VELCHE(/2)) SFLU=VELCHE(1,IBMN) MELVAL=IVAL(4) IBMN=MIN(IB,VELCHE(/2)) EPS =VELCHE(1,IBMN) ELSE SECT=1.D0 MELVAL=IVAL(1) IBMN=MIN(IB,VELCHE(/2)) SCEL=VELCHE(1,IBMN) MELVAL=IVAL(2) IBMN=MIN(IB,VELCHE(/2)) SFLU=VELCHE(1,IBMN) MELVAL=IVAL(3) IBMN=MIN(IB,VELCHE(/2)) EPS =VELCHE(1,IBMN) MELVAL=IVAL(4) IBMN=MIN(IB,VELCHE(/2)) F11 =VELCHE(1,IBMN) MELVAL=IVAL(5) IBMN=MIN(IB,VELCHE(/2)) F12 =VELCHE(1,IBMN) ENDIF c calcul de la masse m0/eps**2 RHOSS=RHOS*SECT/(EPS*EPS) c calcul des coefficients de normalisation COEFPR=(RHOREF*CREF*CREF)/RLCAR COEFPI=RHOREF*RLCAR VKL12 =-(COEFPR*COEFPI*SFLU)/(RHOF*C*C*SCEL) VKL22 =-(COEFPI*COEFPI)/(RHOF*SCEL) VKL23 = COEFPI/SCEL VKL33 = 1.D0/SCEL VKL41 = EPS*EPS/2.D0/SCEL*(COEFPR*COEFPI) VKL42 = EPS*EPS/2.D0/SCEL*COEFPI*COEFPI VKL43 = EPS*EPS/2.D0/SCEL*COEFPI VKL44 = EPS*EPS/2.D0/SCEL c boucle sur les points de gauss c cas plan ISDJC=0 DO 6157 IGAU1=1,NBPGAU POIGA1=MINTE.POIGAU(IGAU1) CALL TRIHM31(IGAU1,MELE,NBNO,XE,SHPTOT,SHPWRK # ,RHOSS,RHOF, # B11,B22,B12,F11,F12,SFLU,SCEL,POIGA1,VKL12,VKL22, # VKL23,VKL33,VKL42,VKL43,VKL44,E111,E112,E121,E122, # E221,E222,LRE,REL,ISDJC,IRRT) IF (IRRT.NE.1) GOTO 7157 6157 CONTINUE IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB c remplissage de xmatri CALL REMPMT(REL,LRE,RE(1,1,ib)) 3157 CONTINUE c impression d un eventuel message d erreur 7157 CONTINUE IF(IRRT.EQ.0) THEN MOTERR(1:4)=NOMTP(MELE) CALL ERREUR(420) ELSE IF(IRRT.EQ.2) THEN INTERR(1) = IB CALL ERREUR(405) ENDIF ENDIF IF(IRRT.EQ.3) CALL ERREUR(421) IF(IRRT.EQ.4) CALL ERREUR(422) GOTO 510 c_______________________________________________________________________ 510 CONTINUE IF (I195.NE.0) THEN INTERR(1)=I195 CALL ERREUR(195) ENDIF IF (I259.NE.0) THEN INTERR(1)=I259 CALL ERREUR(259) ENDIF SEGSUP,MVELCH SEGSUP,WRK1 IF (WRK2.NE.0) SEGSUP,WRK2 IF (WRK3.NE.0) SEGSUP,WRK3 IF (WRK4.NE.0) SEGSUP,WRK4 IF (WRK5.NE.0) SEGSUP,WRK5 mlmots = iinc if (mlmots.gt.0) segsup mlmots mlmots = idua if (mlmots.gt.0) segsup mlmots RETURN END