C ROT3M SOURCE PV090527 26/04/30 21:16:24 12529 ************************************************************************ * * R O T 3 M * --------- * * FONCTION: * --------- * CALCUL DE LA MATRICE DE MUTUELLES POUR L'ELEMENT ROT3 * * 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 LE CHAMELEM DE CARACTERISTIQUE * IPSUPJ (E) POINTEUR SUR LE MAILLAGE SUPPORT DES COURANTS DE FOUCAULT * 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 * 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 BIDIM. * VALMAT(NMATR) TABLEAU DE TRAVAIL * * AUTEUR, DATE DE CREATION: * ------------------------- * * YANN STEPHAN , FEVRIER 1997 (COPIE DE ROT3R) * ************************************************************************ SUBROUTINE ROT3M(NEF,IMAIL,IPMODE,IPCHEM,IPSUPJ,XMATRI, $ ICPR,ICPR2) 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 MAXT REAL*8 RA(NBNN,NBNN) ENDSEGMENT SEGMENT,MMAT1 REAL*8 VALMAT(NMATR) REAL*8 XE(3,NBNN),XE1(3,NBNN) REAL*8 SHP(6,NBNN),GRAD(NDIM,NBNN,NBPGAU) REAL*8 COSD1(3),COSD2(3) ENDSEGMENT POINTEUR MMAT2.MMAT1,MMATX.MMAT1 SEGMENT NOTYPE CHARACTER*16 TYPE(NBTYPE) ENDSEGMENT SEGMENT SGAUSS REAL*8 XGAUSS(3,NBPGAU) REAL*8 DX(NBPGAU) ENDSEGMENT POINTEUR SGX.SGAUSS,SGY.SGAUSS SEGMENT ICPR(NA) SEGMENT ICPR2(NA) CHARACTER*8 CNM CHARACTER*(NCONCH) CONM PARAMETER (NINF=3) INTEGER INFOS(NINF) LOGICAL SELF,NEAR * PERMEABILITE DU VIDE SUR 4PI DATA PM0S4P/1.D-7/ IMODEL = IPMODE CONM = imodel.CONMOD IPMAIL = imodel.IMAMOD CNM = imodel.CMATEE * RECUPERATION DES CARACTERISTIQUES GEOMETRIQUES DU MAILLAGE * ELEMENTAIRE MELEME = IPMAIL NBNN = meleme.NUM(/1) NBELEM = meleme.NUM(/2) * * REMLIR LE TABLEAU INFOS (INFORMATIONS SUR ELEMENT) INFOS(1)=0 INFOS(2)=0 INFOS(3)=NIFOUR * RECUPERATION DES CARACTERISTIQUES D'INTEGRATION * POUR LA MATRICE MUTUELLE (RIGIDITE) DE L'ELEMENT * FINI LIE A NOTRE MAILLAGE if (infmod(/1).lt.5) then write(ioimp,*) 'ROT3 - INFMOD(/1) < 5' call erreur(5) endif MINTE = imodel.INFMOD(5) NBPGAU = minte.POIGAU(/1) * * RECHERCHE LES POINTEURS DES SEGMENTS MELVAL QUI CORRESPONDENT * A LA PERMEABILITE ET L'EPAISSEUR DES ELEMENTS * nomid = 0 nbrobl = 0 nbrfac = 0 IF (CNM.EQ.'ISOTROPE'.OR.CNM.EQ.'ORTHOTRO') THEN NBROBL=1 SEGINI,NOMID LESOBL(1)='PERM' ELSE CALL ERREUR(251) RETURN ENDIF NMATO = nbrobl NMATF = nbrfac NMATR = NMATO+NMATF MOMATR = nomid NBTYPE=1 SEGINI,NOTYPE TYPE(1) ='REAL*8' MOTYPE = NOTYPE CALL KOMCHA(IPCHEM,IPMAIL,CONM,MOMATR,MOTYPE,1,INFOS,3,IVAMAT) IF (IERR.NE.0) RETURN MPTVAL = IVAMAT SEGINI SGX,SGY NDIM = IDIM NDIM2 = IDIM-1 SEGINI,MAXT,MMAT1,MMAT2 MMATX = MMAT1 * * CALCUL DE LA DISTANCE DE REFERENCE * DREF = 0.D0 DO IEL = 1, NBELEM CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE) CALL MAXLT3(NBNN,XE,DARET) DREF = MAX(DREF,DARET) ENDDO IPT1 = IPSUPJ NBSOUJ = IPT1.LISOUS(/1) NBSOU1 = MAX(1,NBSOUJ) * BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL DO 10 IEL=1,NBELEM MMAT1=MMATX * * 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 * CALL COQLOC(NBNN,XE,COSD1,COSD2,XE1) * * ON CALCULE LES FONCTIONS DE FORME AUX POINTS DE GAUSS * CALL ELGAUS(MINTE,MMAT1,SGX,IFOIS,IFOI2) * * LE JACOBIEN EST NEGATIF ,MAILLAGE INCORRECT IF(IFOIS.NE.0.AND.IFOIS.NE.NBPGAU)THEN INTERR(1)=IEL CALL ERREUR(195) GO TO 999 * * CAS OU LE JACOBIEN EST TRES PETIT ELSEIF(IFOI2.EQ.NBPGAU)THEN INTERR(1)=IEL CALL ERREUR (259) GO TO 999 ENDIF * * ON BOUCLE SUR LE MAILLAGE SUPPORT DE COURANTS IPT1 = IPSUPJ IPT2 = IPT1 DO 110 ISOUJ = 1, NBSOU1 IF (NBSOUJ.NE.0) IPT2 = IPT1.LISOUS(ISOUJ) NBELJ=IPT2.NUM(/2) NBNNJ=IPT2.NUM(/1) NBNNT=NBNN+NBNNJ NLIGRP=NBNN NLIGRD=NBNN DO 111 IELJ=1,NBELJ DO IX=1,NBNN DO JX=1,NBNN RA(JX,IX) = 0.D0 ENDDO ENDDO * NEAR=.FALSE. * MMAT1=MMAT2 SGAUSS = SGY * * ON CHERCHE LES COORDONNEES DES NOEUDS DE L'ELEMENT IEL, * DANS LE REPERE GLOBAL * CALL DOXE(XCOOR,IDIM,NBNN,IPT2.NUM,IELJ,XE) * * CALCUL DES COORDONNEES DES NOEUDS DANS LE REPERE LOCAL DE L' * ELEMENT COQUE * CALL COQLOC(NBNN,XE,COSD1,COSD2,XE1) * CALL ELGAUS(MINTE,MMAT1,SGAUSS,JFOIS,JFOI2) * IF (JFOIS.NE.0.AND.JFOIS.NE.NBPGAU)THEN * * LE JACOBIEN EST NEGATIF ,MAILLAGE INCORRECT INTERR(1)=IEL CALL ERREUR(195) GO TO 999 ELSEIF(JFOI2.EQ.NBPGAU)THEN * * CAS OU LE JACOBIEN EST TRES PETIT * INTERR(1)=IEL CALL ERREUR (259) GO TO 999 ENDIF * * CALCUL DE LA DISTANCE ENTRE LES DEUX ELEMENTS * CALL S2ELT3(NBNN,NUM,IEL,IPT2.NUM,IELJ,NEAR,SELF) CALL D2ELT3(NBNN,XE,MMATX.XE,DT3) NEAR=NEAR.OR.(DT3.LE.DREF) * * BOUCLE SUR LES POINTS DE GAUSS MAILLAGE 1 * DO 22 IGAU=1,NBPGAU * MMAT1=MMATX SGAUSS=SGX * * ON CHERCHE LES VALEURS DE LA PERMEABILITE * MPTVAL=IVAMAT DO 30 IM=1,NMATR 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 ENDIF 30 CONTINUE PERM=VALMAT(1) XK=PERM*PM0S4P*DX(IGAU) * * BOUCLE SUR LES POINTS DE GAUSS MAILLAGE 2 * DO 23 JGAU=1,NBPGAU * MMAT1 = MMAT2 SGAUSS = SGY * IF (SELF) THEN IF (JGAU.GT.1) GOTO 23 CALL SELFT3(SGX.XGAUSS(1,IGAU),NBNN,XE,QQ) YK=XK*QQ ELSE IF (NEAR) THEN IF(JGAU.GT.1) GO TO 23 CALL NEART3(SGX.XGAUSS(1,IGAU),NBNN,MMATX.XE,XE,QQ) YK=XK*QQ ELSE DIST=0.D0 DO I=1,IDIM DIST=DIST+(SGX.XGAUSS(I,IGAU)-XGAUSS(I,JGAU))**2 ENDDO DIST=SQRT(DIST) YK=XK*DX(JGAU)/DIST ENDIF * * ON AJOUTE LE PRODUIT K*DJAC*TRANSPOSEE(GRADX)*GRADY * POUR LE POINT DE GAUSS CONSIDERE,A LA MATRICE RE * MMAT1=MMATX CALL WXWYSR(GRAD(1,1,IGAU),MMAT2.GRAD(1,1,JGAU), & YK,NBNN,NBNNJ,IDIM,RA) 23 CONTINUE 22 CONTINUE * realisation de l'assemblage DO IX = 1, IPT2.NUM(/1) IA = IPT2.NUM(IX,IELJ) IB = ICPR2(IA) DO JX = 1, NUM(/1) IC = NUM(JX,IEL) ID = ICPR(IC) RE(IB,ID,1) = RA(IX,JX) + RE(IB,ID,1) ENDDO ENDDO * 111 CONTINUE 110 CONTINUE * 10 CONTINUE * END DO * on symetrise la matrice DO IX = 1, RE(/1) DO JX = 1, IX XP = 0.5D0 * ( RE(IX,JX,1) + RE(JX,IX,1) ) RE(IX,JX,1)=XP RE(JX,IX,1)=XP ENDDO ENDDO * * DESACTIVATION DES SEGMENTS * 999 CONTINUE MMAT1=MMATX SEGSUP,MMAT1,MMAT2,MAXT MPTVAL=IVAMAT SEGSUP,MPTVAL NOMID=MOMATR SEGSUP,NOMID,NOTYPE SEGSUP,SGX,SGY C RETURN END