C ROT3R SOURCE PV090527 26/04/30 21:16:24 12529 ************************************************************************ * * R O T 3 R * --------- * * FONCTION: * --------- * CALCUL DE LA MATRICE DE RESISTANCE POUR L'ELEMENT ROT3 * * PARAMETRES: (E)=ENTREE (S)=SORTIE (+ = CONTENU DANS UN COMMUN) * ----------- * NEF (E) NUMERO DE L'ELEMENT-FINI DANS NOMTP (VOIR CCHAMP) * IPMAIL (E) MAILLAGE ELEMENTAIRE CONSIDERE * IPMODE (E) POINTEUR SUR UN SEGMENT IMODEL * IPCHEM (E) POINTEUR SUR LE CHAMELEM DE CARACTERISTIQUE * IPMATR (E/S) MATRICE DE RIGIDITE ELEMENTAIRE XMATRI ************************************************************************ SUBROUTINE ROT3R(NEF,IPMAIL,IPMODE,IPCHEM,IPMATR) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC CCREEL -INC CCHAMP -INC SMCOORD -INC SMINTE -INC SMMODEL -INC SMRIGID -INC SMELEME -INC SMCHAML -INC TMPTVAL SEGMENT NOTYPE CHARACTER*16 TYPE(NBTYPE) ENDSEGMENT SEGMENT,MMAT1 REAL*8 VALMAT(NMATR) REAL*8 XE(3,NBNN),CEL1(NBNN,NBNN),XE1(3,NBNN) REAL*8 SHP(6,NBNN),GRAD(NDIM,NBNN) ENDSEGMENT REAL*8 COSD1(3),COSD2(3),COSD3(3),YK(2,2) CHARACTER*8 CNM CHARACTER*(NCONCH) CONM PARAMETER (NINF=3) INTEGER INFOS(NINF) IMODEL = IPMODE CONM = imodel.CONMOD MELEME = IPMAIL c* meleme = imodel.IMAMOD NBNN = meleme.NUM(/1) NBELEM = meleme.NUM(/2) * RECUPERATION DES CARACTERISTIQUES D'INTEGRATION DE L'ELEMENT * FINI LIE A NOTRE MAILLAGE if (infmod(/1).lt.4) then write(ioimp,*) 'rot3r infmod(/1)' call erreur(5) endif * INFORMATION SUR L'ELEMENT MINTE = imodel.INFMOD(4) NBPGAU = minte.POIGAU(/1) xMATRI = IPMATR c* SEGACT,xMATRI*MOD NLIGRP = NBNN NLIGRD = NBNN * REMLIR LE TABLEAU INFOS (INFORMATIONS SUR ELEMENT) INFOS(1)=0 INFOS(2)=0 INFOS(3)=NIFOUR * RECHERCHE LES POINTEURS DES SEGMENTS MELVAL QUI CORRESPONDENT * AUX COMPOSANTES DE LA CONDUCTIVITE ET L'EPAISSEUR DES ELEMENT CNM = imodel.CMATEE c* INM = imodel.IMATEE c* INT = imodel.INATUU nbrobl = 0 nbrfac = 0 nomid = 0 MOMATR = nomid NBTYPE = 1 SEGINI,notype TYPE(1) = 'REAL*8' MOTYR8 = notype IF (CNM.EQ.'ISOTROPE') THEN NBROBL=2 SEGINI,NOMID LESOBL(1)='ETA ' LESOBL(2)='EPAI' ELSE IF (CNM.EQ.'ORTHOTRO') THEN NBROBL=5 SEGINI,NOMID LESOBL(1)='ETA1' LESOBL(2)='ETA2' LESOBL(3)='EPAI' LESOBL(4)='V1X ' LESOBL(5)='V1Y ' ELSE CALL ERREUR(251) RETURN ENDIF NMATO = nbrobl NMATF = nbrfac NMATR = NMATO + NMATF MOMATR = nomid IVAMAT = 0 CALL KOMCHA(IPCHEM,IPMAIL,CONM,MOMATR,MOTYR8,1,INFOS,NINF,IVAMAT) IF (IERR.NE.0) GOTO 990 MPTVAL = IVAMAT IF (CNM.EQ.'ISOTROPE')THEN IPMELV = IVAL(2) ELSE IF(CNM.EQ.'ORTHOTRO') THEN IPMELV = IVAL(3) ENDIF CALL QUELCH(IPMELV,ICONS) IF (ICONS.NE.0) THEN CALL ERREUR(566) GOTO 990 ENDIF NDIM = IDIM-1 NFIN = IDIM SEGINI,MMAT1 * BOUCLE (10) SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL DO IEL = 1, NBELEM * ON CHERCHE LES COORDONNEES DES NOEUDS DE L'ELEMENT IEL, * DANS LE REPERE GLOBAL CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE) * CALCUL DES COORDONNEES DES NOEUDS DANS LE REPERE LOCAL DE L' * ELEMENT COQUE COSD1(1) = XE(1,2)-XE(1,1) COSD1(2) = XE(2,2)-XE(2,1) COSD1(3) = XE(3,2)-XE(3,1) COSD2(1) = XE(1,3)-XE(1,1) COSD2(2) = XE(2,3)-XE(2,1) COSD2(3) = XE(3,3)-XE(3,1) COSD3(1) = COSD1(2)*COSD2(3)-COSD1(3)*COSD2(2) COSD3(2) = COSD1(3)*COSD2(1)-COSD1(1)*COSD2(3) COSD3(3) = COSD1(1)*COSD2(2)-COSD1(2)*COSD2(1) COSD1L = SQRT(COSD1(1)*COSD1(1)+COSD1(2)*COSD1(2)+ & COSD1(3)*COSD1(3)) COSD1(1)=COSD1(1)/COSD1L COSD1(2)=COSD1(2)/COSD1L COSD1(3)=COSD1(3)/COSD1L COSD3L = SQRT(COSD3(1)*COSD3(1)+COSD3(2)*COSD3(2)+ & COSD3(3)*COSD3(3)) COSD3(1)=COSD3(1)/COSD3L COSD3(2)=COSD3(2)/COSD3L COSD3(3)=COSD3(3)/COSD3L COSD2(1) = COSD3(2)*COSD1(3)-COSD3(3)*COSD1(2) COSD2(2) = COSD3(3)*COSD1(1)-COSD3(1)*COSD1(3) COSD2(3) = COSD3(1)*COSD1(2)-COSD3(2)*COSD1(1) DO NOE = 1, NBNN XE1(1,NOE) = XE(1,NOE)*COSD1(1) + XE(2,NOE)*COSD1(2) & + XE(3,NOE)*COSD1(3) XE1(2,NOE) = XE(1,NOE)*COSD2(1) + XE(2,NOE)*COSD2(2) & + XE(3,NOE)*COSD2(3) XE1(3,NOE) = 0.D0 ENDDO * MISE A ZERO DU TABLEAU CEL1 CALL ZERO(CEL1,NBNN,NBNN) * BOUCLE (20) SUR LES POINTS DE GAUSS IFOIS = 0 IFOI2 = 0 DO IGAU=1,NBPGAU * CALCUL DE LA MATRCIE GRADIENT DES FONCTIONS DE FORME ET * DU JACOBIEN(DANS LE PLAN), EN UN POINT DE GAUSS DO NOE = 1, NBNN DO I = 1, NFIN SHP(I,NOE) = SHPTOT(I,NOE,IGAU) ENDDO ENDDO * DERIVES DES FONCTIONS DE FORME DANS LA GEOMETRIE REELLE * ET LE JACOBIEN CALL JACOBI(XE1,SHP,NDIM,NBNN,DJAC) IF (DJAC.LT.XZERO) IFOIS=IFOIS+1 IF (ABS(DJAC).LT.XPETIT) IFOI2=IFOI2+1 DJAC=ABS(DJAC)*POIGAU(IGAU) DO NOE = 1, NBNN DO I = 1, NDIM GRAD(I,NOE) = SHP(I+1,NOE) ENDDO ENDDO * ON CHERCHE LES VALEURS DE COMPOSANTES DE LA RESISTIVITE * ET L'EPAISSEUR DE LA COQUE MPTVAL=IVAMAT DO IM=1,NMATR MELVAL=IVAL(IM) IF (MELVAL.NE.0)THEN IBMN=MIN(IEL,VELCHE(/2)) IGMN=MIN(IGAU,VELCHE(/1)) VALMAT(IM)=VELCHE(IGMN,IBMN) ELSE VALMAT(IM)=0.D0 ENDIF ENDDO * MATERIAU ISOTROPE IF (CNM.EQ.'ISOTROPE') THEN EP = VALMAT(2) * L'ELEMENT (IEL) AU POINT DE GAUSS (IGAU) A UNE EPAISSEUR NULLE IF (EP.LE.XPETIT) THEN INTERR(1)=IEL INTERR(2)=IGAU MOTERR(1:4)=NOMTP(NEF) CALL ERREUR(355) GO TO 999 ENDIF XK = VALMAT(1)*DJAC/EP * ON AJOUTE LE PRODUIT K*DJAC*TRANSPOSEE(GRAD)*GRAD * POUR LE POINT DE GAUSS CONSIDERE,A LA MATRICE CEL1 CALL NTNST(GRAD,XK,NBNN,NDIM,CEL1) * CAS ORTHOTROPE ELSE c* IF (CNM.EQ.'ORTHOTRO')THEN EP = VALMAT(3) * L'ELEMENT (IEL) AU POINT DE GAUSS (IGAU) A UNE EPAISSEUR NULLE IF (EP.LE.XPETIT) THEN INTERR(1)=IEL INTERR(2)=IGAU MOTERR(1:4)=NOMTP(NEF) CALL ERREUR(355) GO TO 999 ENDIF XK1 = VALMAT(1) / EP XK2 = VALMAT(2) / EP COSA = VALMAT(4) SINA = -VALMAT(5) * CALCUL DE LA MATRICE DES COEFFICIENTS DE RESISTIVITE DANS LE * PLAN,PAR RAPPORT AU REPERE LOCAL DE L'ELEMENT COS2 = COSA*COSA SIN2 = SINA*SINA SICO = SINA*COSA YK(1,1) = COS2*XK1 + SIN2*XK2 YK(1,2) = SICO*(XK1-XK2) YK(2,1) = YK(1,2) YK(2,2) = SIN2*XK1 + COS2*XK2 * ON AJOUTE LE PRODUIT DJAC*TRANSPOSEE(GRAD)*YK*GRAD * POUR LE POINT DE GAUSS CONSIDERE,A LA MATRICE CEL1 CALL BDBST (GRAD,DJAC,YK,NBNN,NDIM,CEL1) ENDIF ENDDO * FIN BOUCLE (20) SUR LES POINTS DE GAUSS * LE JACOBIEN EST NEGATIF ,MAILLAGE INCORRECT IF (IFOIS.NE.0.AND.IFOIS.NE.NBPGAU) THEN INTERR(1) = IEL CALL ERREUR(195) * CAS OU LE JACOBIEN EST TRES PETIT ELSE IF (IFOI2.EQ.NBPGAU) THEN INTERR(1) = IEL CALL ERREUR (259) ENDIF IF (IERR.NE.0) GO TO 999 * REMPLISSAGE DE XMATRI CALL REMPMT(CEL1,NBNN,RE(1,1,IEL)) ENDDO * FIN BOUCLE (10) SUR LES ELEMENTS * DESACTIVATION DES SEGMENTS 999 CONTINUE SEGSUP,MMAT1 990 CONTINUE mptval = IVAMAT IF (mptval.NE.0) SEGSUP,mptval nomid = MOMATR SEGSUP,nomid notype = MOTYR8 SEGSUP,notype c return END