C LFROA SOURCE PV090527 26/04/30 21:15:49 12529 SUBROUTINE LFROA(IPOGEO,IPMATR,IPMINT,IVAMAT, 1 IVACAR,MELE,MFR,LRE,NDDL) C C*********************************************************************** C * C Routine appelée par FRVISQ. * C * C Calcule l'amortissement de frontière pour les liquides de face LSE2 * C dans le cas 2D, ou de face LTR3 ou LQU4 dans le cas 3D. * C * C Entrées : * C -------- * C * C IPOGEO : pointeur sur le maillage de l'enveloppe des massifs, * C type MELEME * C IPMATR : pointeur sur le segment IMATRI, chapeau des rigidités * C élémentaires * C IPMINT : pointeur sur le segment d'intégration, type MINTE * C IVAMAT : pointeur sur un segment MPTVAL de données matériau * C IVACAR : pointeur sur un segment MPTVAL de caractéristiques * C (épaisseur dans le cas contraintes planes) * C MELE : numéro de l'élément fini associé à la face du massif * C MFR : numéro de la formulation * C LRE : taille de la matrice d'amortissement à construire * C NDDL : nombre de degrés de liberté * C * C Remplit le segment XMATRI pour chaque élément de la sous-zone. * C*********************************************************************** IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC CCREEL -INC SMRIGID -INC SMELEME -INC SMCOORD -INC SMCHAML -INC SMINTE -INC TMPTVAL SEGMENT,MWORK REAL*8 XE(3,NBNN) REAL*8 REL(LRE,LRE) REAL*8 SHPWRK(6,NBNN),BGENE(NDDL,LRE) REAL*8 VALMAT(NV1) REAL*8 RCLMAT(NDDL,NDDL) ENDSEGMENT MELEME=IPOGEO c* SEGACT MELEME NBNN=NUM(/1) NBELEM=NUM(/2) C MINTE=IPMINT c* SEGACT,MINTE NBPGAU=POIGAU(/1) C DIM3=1.D0 C xMATRI=IPMATR c* SEGACT,xMATRI*MOD c* NLIGRD=LRE c* NLIGRP=LRE NV1=5 SEGINI,MWORK C C boucle sur les éléments C DO 1 IB=1,NBELEM C C on cherche les coordonnées de l'élément IB C CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) CALL ZERO(REL,LRE,LRE) CALL ZERO(RCLMAT,NDDL,NDDL) C C boucle sur les points de Gauss C DO 10 IGAU=1,NBPGAU C C récupération de l'épaisseur C IF (IFOUR.EQ.-2) THEN IF (IVACAR.NE.0) THEN MPTVAL=IVACAR MELVAL=IVAL(1) IF (MELVAL.NE.0) THEN IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) DIM3=VELCHE(IGMN,IBMN) ELSE DIM3=1.D0 ENDIF ENDIF ENDIF C C récupération des données matériau C MPTVAL=IVAMAT DO 11 J=1,5 MELVAL=IVAL(J) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) VALMAT(J)=VELCHE(IGMN,IBMN) 11 CONTINUE C RHO=VALMAT(1) CSON=VALMAT(2) ROREF=VALMAT(3) CREF=VALMAT(4) RLCAR=VALMAT(5) C C coefficient d'amortissement normalisé C RCL=(1.D0/CSON/rho)*ROREF*RLCAR*roref*cref**2/rlcar C C matrice RCLMAT C RCLMAT(1,1)=RCL RCLMAT(2,2)=RCL C C calcul de la matrice N des fonctions de forme C XDPGE=0.D0 YDPGE=0.D0 CALL NMATST(IGAU,MELE,MFR,NBNN,LRE,IFOUR,NIFOUR,NDDL, 1 DIM3,XE,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) C C mise à zéro des composantes sur p C DO 21 J=1,LRE BGENE(2,J)=0.D0 21 CONTINUE C C calcul du jacobien C IF (IDIM.EQ.3) THEN DXDQSI=0.D0 DYDQSI=0.D0 DZDQSI=0.D0 DXDETA=0.D0 DYDETA=0.D0 DZDETA=0.D0 DO 40 I=1,NBNN DXDQSI=DXDQSI+SHPTOT(2,I,IGAU)*XE(1,I) DYDQSI=DYDQSI+SHPTOT(2,I,IGAU)*XE(2,I) DXDETA=DXDETA+SHPTOT(3,I,IGAU)*XE(1,I) DYDETA=DYDETA+SHPTOT(3,I,IGAU)*XE(2,I) DZDQSI=DZDQSI+SHPTOT(2,I,IGAU)*XE(3,I) DZDETA=DZDETA+SHPTOT(3,I,IGAU)*XE(3,I) 40 CONTINUE DJAC=SQRT((DYDQSI*DZDETA-DYDETA*DZDQSI)**2+ 1 (DXDETA*DZDQSI-DXDQSI*DZDETA)**2+ 2 (DXDQSI*DYDETA-DXDETA*DYDQSI)**2) ELSE c* ELSE IF(IDIM.EQ.2) THEN DXDQSI=0.D0 DYDQSI=0.D0 DO I=1,NBNN DXDQSI=DXDQSI+SHPTOT(2,I,IGAU)*XE(1,I) DYDQSI=DYDQSI+SHPTOT(2,I,IGAU)*XE(2,I) ENDDO DJAC=SQRT(DXDQSI**2+DYDQSI**2) ENDIF C C calcul de l'élément de volume C IF (IFOUR.LT.0.OR.IFOUR.EQ.2) THEN R=1.D0 IF (IFOUR.EQ.-2) R=DIM3 ELSE R=0.D0 DO I=1,NBNN R=R+SHPTOT(1,I,IGAU)*XE(1,I) ENDDO IF (IFOUR.EQ.0.OR.(IFOUR.EQ.1 1 .AND.NIFOUR.EQ.0)) THEN R=2*XPI*R ELSE c* ELSEIF (IFOUR.EQ.1.AND.NIFOUR.NE.0) THEN R=XPI*R ENDIF ENDIF C C construction de la matrice d'amortissement C DJAC=DJAC*POIGAU(IGAU)*R CALL BDBST(BGENE,DJAC,RCLMAT,LRE,NDDL,REL) C 10 CONTINUE C C remplissage de XMATRI C CALL REMPMT(REL,LRE,RE(1,1,ib)) * ** la matrice calculée ci dessus serait correcte si l'on avait ** une formulation uniquement en pi. Comme on retient la formulation ** en p et pi et pour ne rien ajouter sur la ligne en pi (si non on ** ne satisfait plus la relation entre p et pi) on aboutit à une matrice ** disymetrique avec des termes uniquement sur le lignes en pi et collones ** en p. D'où la matrice suivante : do 30 i = 1, lre do 31 j = 2, lre ix = mod(j,2) if (ix.eq.0) then re(i,j,ib) = -re( i,(j-1),ib) re( i,(j-1),ib) = 0.d0 endif 31 continue 30 continue 1 CONTINUE SEGSUP,MWORK c* SEGDES MELEME,MINTE,xMATRI RETURN END