C RTENS4 SOURCE CB215821 26/08/24 21:18:20 12622 SUBROUTINE RTENS4(IPCHE1,IFOMEM,IMOT,IPTV2,IELEME,IVAVEC,IVACOM, & IVARES,IDEFO,IINTE,MELE,NPINT,NVEC,V1,V2,W2,W3, & CENTR1,CENTR2,AXEI1,ICAS,IER1) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) *-----------------------------------------------------------------------* * Operateur RTENS : cas de la formulation coque avec * * cisaillement transverse (DST, COQ4) * * * * IPCHE1 (e) pointeur sur un MCHAML de caracteristiques * * = 0 si isotropie * * IFOMEM (e) = IFOUR de CCOPTIO * * IMOT (e) indique le type de repere desire (cf RTENS) * * IPTV2 (e) pointeur sur le 2nd point repere * * IELEME (e) pointeur sur le segment MELEME (actif) * * IVAVEC (e/s) pointeur sur un segment MPTVAL (actif) * * IVACOM (e/s) pointeur sur un segment MPTVAL (actif) * * IVARES (e/s) pointeur sur un segment MPTVAL (actif) * * IDEFO (e) =1 : tenseur de deformations (contraintes sinon) * * IINTE (e) pointeur sur le segment MINTE (actif) * * MELE (e) numero de l'element-fini dans NOMTP * * NPINT (e) nombre de points d'integration (coques) * * NVEC (e) nombre de composantes du futur MCHAML * * V1 (e) coordonnees et norme du 1er vecteur * * V2 (e) coordonnees et norme du 2nd vecteur * * W2 (e) coordonnees d'un 1er vecteur de travail * * W3 (e) coordonnees d'un 2nd vecteur de travail * * CENTR1 (e) coordonnees du 1er point repere * * CENTR2 (e) coordonnees du 2nd point repere * * AXEI1 (e) coordonnees du vecteur de l'axe de symetrie * * ICAS (e) distingue les differentes syntaxes d'appel * * IER1 (s) code d'erreur pour desactivation dans RTENS * * D.R.-M. le 18/3/94 * *-----------------------------------------------------------------------* -INC PPARAM -INC CCOPTIO -INC CCHAMP -INC SMCHAML -INC SMINTE -INC SMCOORD -INC SMELEME -INC TMPTVAL * * MWRK1,3 initialises dans RTENS4 * SEGMENT MWRK1 REAL*8 XEL(3,NBNN),XEL2(3,NBNN) ENDSEGMENT * SEGMENT MWRK3 REAL*8 A(NDIM,NDIM),R(NDIM,NDIM),RT(NDIM,NDIM),TRAV(NDIM,NDIM) ENDSEGMENT * DIMENSION BPSS(3,3),T(9),VECWRK(3),V1(4),V2(4),W2(3),W3(3) DIMENSION CENTR1(3),CENTR2(3),AXEI1(3),VECX(3),VECY(3) DIMENSION UR(3),UTHETA(3),UPHI(3),UN(3),UT(3),XIGAU(3) * NDIM=3 IER1 = 0 MELEME = IELEME NBNN = NUM(/1) NBELEM = NUM(/2) MINTE = IINTE NBPGAU = POIGAU(/1) * SEGINI MWRK3 IF (ICAS.NE.1) SEGINI MWRK1 CALL ZERO(R,3,3) * * Boucle sur les elements * DO 1101 IB=1,NBELEM * IF (ICAS.NE.1) THEN CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XEL) IF (MELE.EQ.49) THEN CALL CQ4LOC (XEL,XEL2,BPSS,IRRT,0) ELSE IF (MELE.EQ.93) THEN CALL VPAST(XEL,BPSS) ENDIF ENDIF * IF (ICAS.EQ.2) THEN * * On veut le tenseur dans un repere defini par * un ou deux vecteurs * VL11=BPSS(1,1)*V1(1)+BPSS(1,2)*V1(2)+BPSS(1,3)*V1(3) VL12=BPSS(2,1)*V1(1)+BPSS(2,2)*V1(2)+BPSS(2,3)*V1(3) VL1N=SQRT(VL11**2+VL12**2) IF (VL1N.EQ.0.) THEN CALL ERREUR(344) IER1 = 1 GOTO 1100 ENDIF IF (IPTV2.NE.0) THEN VL21=BPSS(1,1)*V2(1)+BPSS(1,2)*V2(2)+BPSS(1,3)*V2(3) VL22=BPSS(2,1)*V2(1)+BPSS(2,2)*V2(2)+BPSS(2,3)*V2(3) VL2N=SQRT(VL21**2+VL22**2) IF (VL2N.EQ.0.) THEN CALL ERREUR(344) IER1 = 1 GOTO 1100 ENDIF WL33=(VL11*VL22-VL12*VL21)/(VL1N*VL2N) WL21=(-WL33*VL12)/VL1N WL22=(WL33*VL11)/VL1N ENDIF * * Construction de la matrice de rotation qui fait * passer du repere local a la coque au nouveau repere * local defini a partir de la projection des vecteurs * V1 et V2 sur la coque * IF (IPTV2.EQ.0) THEN R(1,1)= VL11/VL1N R(2,1)= VL12/VL1N R(1,2)= -VL12/VL1N R(2,2)= VL11/VL1N R(3,3)= 1.D0 ELSE R(1,1)= VL11/VL1N R(2,1)= VL12/VL1N R(1,2)= WL21 R(2,2)= WL22 R(3,3)= WL33 ENDIF SIGFLX=R(3,3) CALL TRSPOD(R,NDIM,NDIM,RT) ENDIF * * Boucle sur les points de Gauss * DO 1090 IGAU=1,NBPGAU * IF (ICAS.EQ.1) THEN * * On veut le tenseur dans le repere d'orthotropie * MPTVAL=IVAVEC * MELVAL=IVAL(1) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) R(1,1)=VELCHE(IGMN,IBMN) * MELVAL=IVAL(2) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) R(2,1)=VELCHE(IGMN,IBMN) * RN=R(1,1)*R(1,1)+R(2,1)*R(2,1) IF (RN.EQ.0.D0) THEN CALL ERREUR(344) IER1 = 1 GOTO 1100 ENDIF * R(1,2)=-R(2,1) R(2,2)= R(1,1) R(3,3)= 1.D0 SIGFLX= 1.D0 CALL TRSPOD(R,NDIM,NDIM,RT) ENDIF * IF (ICAS.EQ.3) THEN * * Matrice de passage entre le repere local de la coque * et la projection sur celle-ci du repere global choisi * CALL RTENS5(IMOT,4,IGAU,NDIM,V1,CENTR1,CENTR2,BPSS, & SHPTOT,XEL,NBNN,NBPGAU,R,SIGFLX,IER1) IF (IER1.NE.0) THEN IF (IER1.EQ.1) CALL ERREUR(344) IF (IER1.EQ.2) CALL ERREUR(642) GOTO 1100 ENDIF CALL TRSPOD(R,NDIM,NDIM,RT) ENDIF * * Rotation du tenseur pour les termes de membrane * CALL ZERO(A,3,3) * MPTVAL=IVACOM MELVAL=IVAL(1) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) A(1,1) = VELCHE(IGMN,IBMN) * MELVAL=IVAL(2) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) A(2,2) = VELCHE(IGMN,IBMN) * MELVAL=IVAL(3) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) A(1,2) = VELCHE(IGMN,IBMN) * MELVAL=IVAL(7) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) A(1,3) = VELCHE(IGMN,IBMN) * MELVAL=IVAL(8) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) A(2,3) = VELCHE(IGMN,IBMN) * MELVAL=IVAL(4) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) AUX4=VELCHE(IGMN,IBMN) * MELVAL=IVAL(5) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) AUX5 = VELCHE(IGMN,IBMN) * MELVAL=IVAL(6) IGMN=MIN(IGAU,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) AUX6 = VELCHE(IGMN,IBMN) * IF (IDEFO.EQ.1) A(1,2)=A(1,2)/2.D0 IF (IDEFO.EQ.1) A(1,3)=A(1,3)/2.D0 IF (IDEFO.EQ.1) A(2,3)=A(2,3)/2.D0 A(2,1)=A(1,2) A(3,1)=A(1,3) A(3,2)=A(2,3) * t * >>> Changement de repere : A = R A R <<< * CALL MULMAT(TRAV,A,R,3,3,3) CALL MULMAT(A,RT,TRAV,3,3,3) * MPTVAL=IVARES * MELVAL=IVAL(1) VELCHE(IGAU,IB)=A(1,1) * MELVAL=IVAL(2) VELCHE(IGAU,IB)=A(2,2) * IF (IDEFO.EQ.1) A(1,2)=A(1,2)*2.D0 IF (IDEFO.EQ.1) A(1,3)=A(1,3)*2.D0 IF (IDEFO.EQ.1) A(2,3)=A(2,3)*2.D0 * MELVAL=IVAL(3) VELCHE(IGAU,IB) = A(1,2) * MELVAL=IVAL(7) VELCHE(IGAU,IB)=A(1,3) * MELVAL=IVAL(8) VELCHE(IGAU,IB)=A(2,3) * * Rotation du tenseur pour les termes de flexion * CALL ZERO(A,3,3) * A(1,1)=SIGFLX*AUX4 A(2,2)=SIGFLX*AUX5 A(1,2)=SIGFLX*AUX6 IF (IDEFO.EQ.1) A(1,2)=A(1,2)/2.D0 A(2,1)=A(1,2) * t * >>> Changement de repere : A = R A R <<< * CALL MULMAT(TRAV,A,R,3,3,3) CALL MULMAT(A,RT,TRAV,3,3,3) * MELVAL=IVAL(4) VELCHE(IGAU,IB)=A(1,1) * MELVAL=IVAL(5) VELCHE(IGAU,IB)=A(2,2) * IF (IDEFO.EQ.1) A(1,2)=A(1,2)*2.D0 * MELVAL=IVAL(6) VELCHE(IGAU,IB)=A(1,2) * * Fin des deux boucles * 1090 CONTINUE 1101 CONTINUE 1100 CONTINUE SEGSUP MWRK3 IF (ICAS.NE.1) SEGSUP MWRK1 RETURN END