C SORE2 SOURCE CHAT 11/03/16 21:31:40 6902 SUBROUTINE SORE2(NEF,IMAIL,IPMODE,IPCHE1,IPCHEM,IPCHE4,IPRIGI) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) ************************************************************************ * * SORE2 * ------ * CREATION DE LA MATRICE DE CONDUCTIVITE N DIV(GRAD T) * ( EFFET SORET) * * FONCTION: * --------- * TRAITEMENT DU CAS DES ELEMENTS-FINIS MASSIFS A INTEGRATION * NUMERIQUE,POUR UN MAILLAGE ELEMENTAIRE * * * PARAMETRES: (E)=ENTREE (S)=SORTIE (+ = CONTENU DANS UN COMMUN) * ----------- * * NEF (E) NUMERO DE L'ELEMENT-FINI DANS NOMTP (VOIR CCHAMP) * IMAIL (E) NUMERO DU MAILLAGE ELEMENTAIRE CONSIDERE,DANS * L'OBJET MODELE * IPMODE (E) POINTEUR SUR UN SEGMENT IMODEL * IPCHEM (E) POINTEUR SUR UN CHAMELEM (GRAD T aux PTS de GAUSS) * IPCHE1 (E) POINTEUR SUR UN CHAMELEM MATERIAU * IPCHE4 (E) POINTEUR SUR UN CHAMELEM FACTEUR DE GRAD T * * IPRIGI (E/S) POINTEUR SUR L'OBJET RESULTAT,DE TYPE RIGIDITE * * VARIABLES: * ---------- * * NBNN NOMBRE DE NOEUDS DANS L'ELEMENT CONSIDERE * NEF NUMERO DE L'ELEMENT FINI DANS NOMTP (VOIR CCHAMP) * NBELEM NOMBRE D'ELEMENTS DANS LE MAILLAGE ELEMENTAIRE * NBPGAU NOMBRE DE POINTS DE GAUSS DANS L'ELEMENT-FINI * NDIM NOMBRE DE LIGNES DE LA MATRICE GRADIENT * CEL(NBNN,NBNN) MATRICE DE CONDUCTIVITE ELEMENTAIRE NON SYMETRIQUE * XE(3,NBNN) COORDONNEES DE L'ELEMENT DANS LE REPERE GLOBAL * SHP(6,NBNN) TABLEAU DE TRAVAIL * GRAD(NDIM,NBNN) MATRICE GRADIENT DES FONCTIONS DE FORME * * * * AUTEUR, DATE DE CREATION: * ------------------------- * * J.M.BAZE AVRIL 97 * * LANGAGE: * -------- * * ESOPE + FORTRAN77 * ************************************************************************ * * MODULES UTILISES: * ----------------- * -INC CCOPTIO -INC CCREEL *- -INC SMCOORD -INC SMINTE C-INC CCHAMP -INC SMMODEL -INC SMRIGID -INC SMELEME -INC SMCHAML * SEGMENT MPTVAL INTEGER IPOS(NS) ,NSOF(NS) INTEGER IVAL(NCOSOU) CHARACTER*16 TYVAL(NCOSOU) ENDSEGMENT * SEGMENT ,MAXE REAL*8 TXR(IDIM,IDIM),XLOC(IDIM,IDIM),XGLOB(IDIM,IDIM) ENDSEGMENT * SEGMENT,MMAT1 REAL*8 CEL(NBNN,NBNN),XE(3,NBNN) REAL*8 SHP(6,NBNN),GRAD(NDIM,NBNN) REAL*8 CMAT(NDIM,NDIM),CMAT1(IDIM,IDIM) ENDSEGMENT * SEGMENT NOTYPE CHARACTER*16 TYPE(NBTYPE) ENDSEGMENT * SEGMENT,MVELCH REAL*8 GDT(IDIM),VALMAT(NV1), GDTL(IDIM) ENDSEGMENT * CHARACTER*8 CNM CHARACTER*(NCONCH) CONM PARAMETER(NINF=3) INTEGER INFOS(NINF) LOGICAL lsupma * * RECUPERATION DES CARACTERISTIQUES GEOMETRIQUES DU MAILLAGE * ELEMENTAIRE * IMODEL=IPMODE SEGACT,IMODEL IPMAIL=IMAMOD MELEME=IMAMOD SEGACT,MELEME NBNN=NUM(/1) NBELEM=NUM(/2) * * RECUPERATION DES CARACTERISTIQUES D'INTEGRATION DE L'ELEMENT * FINI LIE A NOTRE MAILLAGE * CALL TSHAPE(NEF,'GAUSS',IPINTE) IF (IERR.NE.0) THEN SEGDES,IMODEL,MELEME RETURN ENDIF MINTE=IPINTE SEGACT,MINTE NBPGAU=POIGAU(/1) * RECUPERATION DES FONCTIONS DE FORME ET LEURS DERIVEES AU CENTRE * DE L'ELEMENT POUR LE CALCUL DES AXES LOCAUX * NLG=NUMGEO(NEF) CALL RESHPT(1,NBNN,NLG,NEF,0,IPT1,IRT1) MINTE1=IPT1 SEGACT MINTE1 * MRIGID=IPRIGI SEGACT,MRIGID xMATRI=IRIGEL(4,IMAIL) SEGACT,xMATRI*MOD NLIGRP=NBNN NLIGRD=NBNN *-------------------------- * RECHERCHE POINTEUR DU SEGMENTS MELVAL CONTENANT * LA DIFFUSIVITE NFOR=FORMOD(/2) NMAT=MATMOD(/2) CALL NOMATE(FORMOD,NFOR,MATMOD,NMAT,CNM,INM,INT) IF (IERR.NE.0) THEN SEGDES,IMODEL,MELEME SEGDES,MINTE,MINTE1 SEGSUP,MRIGID,xMATRI RETURN ENDIF * * REMLIR LE TABLEAU INFOS (INFORMATIONS SUR ELEMENT) INFOS(1)=0 INFOS(1)=0 INFOS(2)=0 INFOS(3)=NIFOUR * MFR=1 if(lnomid(6).ne.0) then nomid=lnomid(6) segact nomid momatr=nomid nmatr=lesobl(/2) nmatf=lesfac(/2) lsupma=.false. else lsupma=.true. CALL IDMATR(MFR,IMODEL,MOMATR,NMATR,NMATF) endif NMATT=NMATR NV1=NMATT * NBTYPE=1 SEGINI NOTYPE MOTYPE=NOTYPE TYPE(1)='REAL*8' * CALL KOMCHA(IPCHE1,IPMAIL,CONM,MOMATR,MOTYPE,1,INFOS,3,IVAMAT) SEGSUP NOTYPE IF(IERR.NE.0)THEN SEGSUP,xMATRI,MRIGID SEGDES,MELEME SEGDES,IMODEL SEGDES,MINTE,MINTE1 RETURN ENDIF *---------------------------- * recuperation des MELVAL composantes du gradient aux pts de Gauss * et de leurs multiplicateurs MCHEL1=IPCHEM MCHEL2=IPCHE4 SEGACT,MCHEL1,MCHEL2 MCHAM1 = MCHEL1.ICHAML(1) MCHAM2 = MCHEL2.ICHAML(1) SEGACT,MCHAM1,MCHAM2 SEGDES,MCHEL1,MCHEL2 * MELVA1=MCHAM1.IELVAL(1) MELVA2=MCHAM1.IELVAL(2) SEGACT,MELVA1,MELVA2 IF(IDIM.EQ.3) THEN MELVA3=MCHAM1.IELVAL(3) SEGACT,MELVA3 ENDIF MELVA4 =MCHAM2.IELVAL(1) SEGACT,MELVA4 * NDIM=IDIM SEGINI,MMAT1,MVELCH,MAXE * * BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL * DO 10 IEL=1,NBELEM * * MISE A ZERO DU TABLEAU CEL * CALL ZERO(CEL,NBNN,NBNN) * ON CHERCHE LES COORDONNEES DES NOEUDS DE L'ELEMENT IEL, * DANS LE REPERE GLOBAL * CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE) * * CALCUL DES AXES LOCAUX DANS LE CAS ORTHO * IF(INM.EQ.2)THEN NBSH=MINTE1.SHPTOT(/2) CALL RLOCAL(XE,MINTE1.SHPTOT,NBSH,NBNN,TXR) if (nbsh.eq.-1) then call erreur(525) return endif ENDIF * BOUCLE SUR LES POINTS DE GAUSS * IFOIS=0 IFOI2=0 DO 20 IGAU=1,NBPGAU * * CALCUL DE LA MATRICE GRADIENT DES FONCTIONS DE FORME ET * DU JACOBIEN,EN UN POINT DE GAUSS * CALL TCOND5(IGAU,NBNN,NDIM,XE,SHPTOT,SHP,GRAD,DJAC) * SI IFOMOD = 0 axi DJAC CONTITNT DEJA LE MULTIPLICATEUR 2*XPI*R IF (IERR.NE.0) THEN SEGSUP,xMATRI,MRIGID GOTO 99 ENDIF IF (DJAC.LT.XZERO) THEN IFOIS=IFOIS+1 ENDIF IF (ABS(DJAC).LT.XPETIT) THEN IFOI2=IFOI2+1 ENDIF * * ON MULTIPLIE LE JACOBIEN PAR LE POIDS D'INTEGRATION, * POUR LE POINT DE GAUSS CONSIDERE * DJAC=ABS(DJAC)*POIGAU(IGAU) * VALEURS DES COMPOSANTES DES GRADIENTS DO 29 IM=1,IDIM MELVAL=MCHAM1.IELVAL(IM) IBMN=MIN(IEL,VELCHE(/2)) IGMN=MIN(IGAU,VELCHE(/1)) GDT(IM)=VELCHE(IGMN,IBMN) 29 CONTINUE * FACTEUR DE GRAD T IBMN=MIN(IEL,MELVA4.VELCHE(/2)) IGMN=MIN(IGAU,MELVA4.VELCHE(/1)) RMUL = MELVA4.VELCHE(IGMN,IBMN) * * DIFFUSIVITE DANS LE REPERE LOCAL * MPTVAL=IVAMAT DO 30 IM=1,NMATT IF(IVAL(IM).NE.0)THEN MELVAL=IVAL(IM) IBMN=MIN(IEL,VELCHE(/2)) IGMN=MIN(IGAU,VELCHE(/1)) VALMAT(IM)=VELCHE(IGMN,IBMN) ELSE VALMAT(IM)=0.D0 ENDIF 30 CONTINUE * IF (INM.EQ.1) THEN *------------------------ MATERIAU ISOTROPE ---------------------------- * * INTEGRATION DES TERMES N VI B * RK = VALMAT(1)*DJAC*RMUL DO 700 K=1,IDIM XK = GDT(K)*RK DO 300 I=1,NBNN DO 400 J=1,NBNN CEL(J,I) = CEL(J,I) + SHP(1,I)*GRAD(K,J)*XK 400 CONTINUE 300 CONTINUE 700 CONTINUE * ELSE *------------------- MATERIAU ORTHOTROPE ----------------- CALL ZERO(CMAT,NDIM,NDIM) CALL ZERO(CMAT1,IDIM,IDIM) CALL ZERO(XGLOB,IDIM,IDIM) IF(IDIM.EQ.2) THEN *----------BIDIM * MATERIAU ORTHOTROPE XLOC(1,1)=VALMAT(3) XLOC(2,1)=VALMAT(4) XLOC(1,2)=-VALMAT(4) XLOC(2,2)=VALMAT(3) * CALCUL DES COS.DIRECTEURS DES AXES ORTH. /REPERE GLOBAL * XGLOB=TXR*XLOC * DO 40 K=1,IDIM DO 409 J=1,IDIM DO 4099 I=1,IDIM XGLOB(K,J)=TXR(J,I)*XLOC(I,K)+XGLOB(K,J) 4099 CONTINUE 409 CONTINUE 40 CONTINUE * DO 51 I = 1,IDIM CMAT1(I,I) = VALMAT(I) 51 CONTINUE * * RETOUR DANS LE REPERE GLOBAL CALL PRODT(CMAT,CMAT1,XGLOB,IDIM,IDIM) DO 41 I= 1, IDIM GDTL(I) = 0.D0 DO 411 J= 1, IDIM GDTL(I) = GDTL(I) + CMAT(I,J)*GDT(J) 411 CONTINUE 41 CONTINUE * ELSE *----------TRIDIM MATERIAU ORTHOTROPE ------------------- XLOC(1,1)=VALMAT(4) XLOC(2,1)=VALMAT(5) XLOC(3,1)=VALMAT(6) XLOC(1,2)=VALMAT(7) XLOC(2,2)=VALMAT(8) XLOC(3,2)=VALMAT(9) CALL CROSS2 (XLOC(1,1),XLOC(1,2),XLOC(1,3),IRR) DO 45 K=1,IDIM DO 451 J=1,IDIM DO 452 I=1,IDIM XGLOB(K,J)=TXR(J,I)*XLOC(I,K)+XGLOB(K,J) 452 CONTINUE 451 CONTINUE 45 CONTINUE * DO 52 I = 1,IDIM CMAT1(I,I) = VALMAT(I) 52 CONTINUE * * RETOUR DANS LE REPERE GLOBAL CALL PRODT(CMAT,CMAT1,XGLOB,IDIM,IDIM) DO 53 I= 1, IDIM GDTL(I) = 0.D0 DO 531 J= 1, IDIM GDTL(I) = GDTL(I) + CMAT(I,J)*GDT(J) 531 CONTINUE 53 CONTINUE * ENDIF * INTEGRATION DES TERMES N VI B * DO 701 K=1,IDIM XK = GDTL(K)*DJAC*RMUL DO 301 I=1,NBNN DO 401 J=1,NBNN CEL(J,I) = CEL(J,I) + SHP(1,I)*GRAD(K,J)*XK 401 CONTINUE 301 CONTINUE 701 CONTINUE ENDIF * * FIN DE LA BOUCLE SUR LES POINTS DE GAUSS * 20 CONTINUE * END DO IF (IFOIS.NE.0.AND.IFOIS.NE.NBPGAU) THEN * * LE JACOBIEN EST NEGATIF,MAILLAGE INCORRECT * INTERR(1)=IEL CALL ERREUR(195) SEGSUP,xMATRI,MRIGID GO TO 99 ENDIF IF (IFOI2.EQ.NBPGAU) THEN * * CAS OU LE JACOBIEN EST TRES PETIT * INTERR(1)=IEL CALL ERREUR(259) SEGSUP,xMATRI,MRIGID GO TO 99 ENDIF 1100 FORMAT (4E12.5) * SEGINI,XMATRI * IMATTT(IEL)=XMATRI * * REMPLISSAGE DE XMATRI * DO 100 IA=1,NBNN DO 1001 IB=1,NBNN RE(IA,IB,iel)=CEL(IA,IB) 1001 CONTINUE 100 CONTINUE * SEGDES,XMATRI 10 CONTINUE * END DO * * DESACTIVATION DES SEGMENTS * SEGDES,xMATRI,MRIGID 99 CONTINUE SEGSUP,MMAT1,MVELCH 999 CONTINUE SEGDES,MELEME*NOMOD SEGDES,IMODEL*NOMOD SEGDES,MINTE*NOMOD * SEGDES,MELVA4 DO 550 I=1,IDIM MELVAL=MCHAM1.IELVAL(I) SEGDES,MELVAL 550 CONTINUE * MPTVAL = IVAMAT DO 551 I=1,NMATT MELVAL=IVAL(I) SEGDES,MELVAL 551 CONTINUE SEGSUP,MPTVAL nomid=momatr segdes nomid if(lsupma) segsup nomid SEGDES,MCHAM1*NOMOD RETURN END