fpma3d
C FPMA3D SOURCE OF166741 25/02/06 21:15:05 12146 C____________________________________________________________________ C CALCULE LES FORCES DE PRESSIONS SUR LES FACES D ELEMENTS C MASSIFS TRIDIMENSIONNELS C C ENTREES : C --------- C C IPTVPR POINTEUR SUR UN MELVAL CONTENANT LES PRESSIONS APPLIQUEES C 0 SI ON A DONNE UNE PRESSION CONSTANTE C IPMAIL POINTEUR SUR UN OBJET GEOMETRIQUE C IPTINT POINTEUR SUR UN MINTE CONTENANT LES POINTS D INTEGRATION C ACTIF EN ENTREE ET EN SORTIE SANS MODIFICATION C IVAFOR POINTEUR SUR UN MPTVAL ET LES MELVAL CONTENANT LES FORCES C NODALES RESUL C C JACQUELINE BROCHARD AVRIL 85 C C PASSAGE AUX NOUVEAUX CHAMELEM PAR JM CAMPENON LE 17 09 90 C C______________________________________________________________________ & ,netn1,ietn1) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC SMCHAML -INC SMELEME -INC SMINTE -INC SMCOORD segment netn(notn) segment ietn(letn) SEGMENT WORK REAL*8 XE(3,NBNN) ENDSEGMENT SEGMENT MPTVAL INTEGER IPOS(NS),NSOF(NS) INTEGER IVAL(NCOSOU) CHARACTER*16 TYVAL(NCOSOU) ENDSEGMENT idimp1 = IDIM+1 * prob optimiseur il faut initialiser melva1 MELVA1=IVAFOR IF (IPTVPR.NE.0) THEN MELVA1=IPTVPR c* SEGACT,MELVA1 <- ACTIF EN E/S IG11 = MELVA1.VELCHE(/1) IB12 = MELVA1.VELCHE(/2) ENDIF MINTE=IPTINT C* SEGACT,MINTE <- ACTIF EN E/S NBPGAU=POIGAU(/1) MELEME = IPMAIL c* SEGACT,MELEME <- ACTIF EN E/S NBNN = meleme.NUM(/1) NBELEM = meleme.NUM(/2) SEGINI,WORK netn = netn1 ietn = ietn1 ipt1 = ipmaim IF (IPT1.GT.0) THEN if (netn.eq.0 .or. ietn.eq.0) then write(ioimp,*) 'FPMA3D : incompatibilite netn, ietn & IPMAIM' endif c* SEGACT,IPT1 <- ACTIF en E/S nbnn1 = ipt1.num(/1) nbel1 = ipt1.num(/2) ELSE if (netn.gt.0 .or. ietn.gt.0) then write(ioimp,*) 'FPMA2D : incompatibilite netn, ietn & IPMAIM' endif ENDIF C C BOUCLE SUR LES ELEMENTS C DO IB = 1, NBELEM XFLOT = +1.D0 IF (netn.GT.0) THEN DO inf = 1, nbnn ip = meleme.num(inf,ib) ideb = netn(ip)+1 ifin = netn(ip+1) do itn = ideb, ifin IEM = ietn(itn) jne = 0 do i = 1, nbnn ino = num(i,ib) do i1 = 1, nbnn1 if (ino.eq.ipt1.num(i1,IEM)) jne=jne+1 enddo enddo if (jne.eq.nbnn) goto 170 enddo ENDDO GOTO 9900 170 continue XG = 0.D0 YG = 0.D0 ZG = 0.D0 DO I = 1, NBNN1 ino = (IPT1.NUM(I,IEM)-1)*idimp1 XG=XG+XCOOR(ino+1) YG=YG+XCOOR(ino+2) ZG=ZG+XCOOR(ino+3) ENDDO XG=XG / NBNN1 YG=YG / NBNN1 ZG=ZG / NBNN1 XK=0.D0 YK=0.D0 ZK=0.D0 DO i = 1,NBNN XK=XK+XE(1,I) YK=YK+XE(2,I) ZK=ZK+XE(3,I) ENDDO XK=XK/NBNN YK=YK/NBNN ZK=ZK/NBNN V_1 = XG - XK V_2 = YG - YK V_3 = ZG - ZK r_z = 1.D0 / SQRT(V_1*V_1+V_2*V_2+V_3*V_3) V_1 = V_1 * r_z V_2 = V_2 * r_z V_3 = V_3 * r_z ENDIF C C BOUCLE SUR LES POINTS DE GAUSS C DO IGAU = 1, NBPGAU VNQSI1 = 0.D0 VNQSI2 = 0.D0 VNQSI3 = 0.D0 VNETA1 = 0.D0 VNETA2 = 0.D0 VNETA3 = 0.D0 DO I = 1, NBNN VNQSI1 = VNQSI1+SHPTOT(2,I,IGAU)*XE(1,I) VNQSI2 = VNQSI2+SHPTOT(2,I,IGAU)*XE(2,I) VNQSI3 = VNQSI3+SHPTOT(2,I,IGAU)*XE(3,I) VNETA1 = VNETA1+SHPTOT(3,I,IGAU)*XE(1,I) VNETA2 = VNETA2+SHPTOT(3,I,IGAU)*XE(2,I) VNETA3 = VNETA3+SHPTOT(3,I,IGAU)*XE(3,I) ENDDO VNN1 = VNQSI2*VNETA3-VNQSI3*VNETA2 VNN2 = VNQSI3*VNETA1-VNQSI1*VNETA3 VNN3 = VNQSI1*VNETA2-VNQSI2*VNETA1 if (igau.eq.1.and.netn.gt.0) then vnnn = 1.D0 / SQRT(vnn1*vnn1+vnn2*vnn2+vnn3*vnn3) endif r_z = POIGAU(IGAU) * XFLOT IF (IPTVPR.NE.0) THEN IGMN=MIN(IGAU,IG11) IBMN=MIN(IB ,IB12) r_z = r_z * MELVA1.VELCHE(IGMN,IBMN) ELSE r_z = r_z * XP ENDIF T3 = r_z * VNN3 MPTVAL=IVAFOR MELVAL=IVAL(1) DO i=1,NBNN ENDDO MELVAL=IVAL(2) DO i=1,NBNN ENDDO MELVAL=IVAL(3) DO i=1,NBNN VELCHE(i,IB) = VELCHE(i,IB) + SHPTOT(1,i,IGAU)*T3 ENDDO ENDDO ENDDO 9900 CONTINUE SEGSUP,WORK c* SEGDES,MINTE <- ACTIF en E/S c* SEGDES,MELEME <- ACTIF en E/S c* IF (IPTVPR.NE.0) SEGDES,MELVA1 <- ACTIF en E/S RETURN END
© Cast3M 2003 - Tous droits réservés.
Mentions légales