C DMMU SOURCE CHAT 11/03/16 21:19:16 6902 SUBROUTINE DMMU C----------------------------------------------------------------------- C -1 t C Calcul du produit D (M -U) B CHPface C -1 t C ou du produit D M B CHPface C le résultat est un CHPcentre C----------------------------------------------------------------------- C C--------------------------- C Phrase d'appel (GIBIANE) : C--------------------------- C C CHP3 = 'DMMU' MMODEL RIG1 (CHP1) CHP2 ; C C------------------------ C Operandes et resultat : C------------------------ C C MMODEL : modèle formulation DARCY. C RIG1 : Matrices hybrides elementaires de DARCY crees par MHYB. C CHP1 : flux aux faces C CHP2 : CHPO face à plusieur composantes. C CHP3 : CHPO centre à plusieur composantes dont le support géometrique C est le maillage centre du domaine les noms des composantes C sont ceux de CHP2 C C----------------------------------------------------------------------- C C Langage : ESOPE + FORTRAN77 C C C----------------------------------------------------------------------- IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) * -INC CCOPTIO -INC SMCOORD -INC SMCHPOI -INC SMCHAML -INC SMELEME -INC SMRIGID -INC SMTABLE -INC SMMODEL * SEGMENT ICCPR INTEGER ICPR(NNGOT) ENDSEGMENT SEGMENT IORGA INTEGER IAA(ITES) ENDSEGMENT * LOGICAL LOGRE,LOGIN CHARACTER*8 TAPIND,TYPOBJ,CHARRE,LETYPE CHARACTER*4 NOMTOT(1) * * Initialisations * IVALIN = 0 XVALIN = 0.D0 LOGIN = .TRUE. IOBIN = 0 TAPIND = 'MOT ' * * * Lecture du MMODEL * CALL LIROBJ('MMODEL',IPMODE,1,IRET) IF(IERR.NE.0)RETURN MMODEL=IPMODE SEGACT MMODEL N1=KMODEL(/1) DO 7 I=1,N1 IMODEL=KMODEL(I) SEGACT IMODEL IF(FORMOD(1).NE.'DARCY')THEN MOTERR(1:16) = 'DARCY ' CALL ERREUR(719) RETURN ENDIF 7 CONTINUE C on récupère la table DOMAINE IPTABL = 0 CALL LEKMOD(MMODEL,IPTABL,IRET) IF (IERR.NE.0) RETURN TYPOBJ = 'MAILLAGE' CALL LEKTAB(IPTABL,'ELTFA',IOBRE) IF (IERR.NE.0) RETURN IELTFA = IOBRE CALL LEKTAB(IPTABL,'CENTRE',IOBRE) IF (IERR.NE.0) RETURN ICENTR = IOBRE CALL LEKTAB(IPTABL,'FACE',IOBRE) IF (IERR.NE.0) RETURN IFACE = IOBRE * * Lecture de RIGIDITE * CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET) IF (IERR.NE.0) RETURN MRIGID = IPRIGI * * * * Test du sous-type de la matrice de rigiditée récupérée * SEGACT MRIGID LETYPE = MTYMAT IF (LETYPE.NE.'DARCY') THEN MOTERR(1:8) = 'RIGIDITE' MOTERR(9:16) = 'DARCY ' CALL ERREUR(79) SEGDES MRIGID GOTO 100 ENDIF * * Controle des pointeurs de MELEME de la rigidité * NRIGEL=IRIGEL(/2) C write(6,*)' NRIGEL----- ' ,nrigel MELEME=IELTFA SEGACT MELEME NBSOUS=LISOUS(/1) IF(NBSOUS.EQ.0)THEN IF((NRIGEL.NE.1).OR.(IRIGEL(1,1).NE.MELEME))THEN MOTERR(1:8) = 'DARCY ' MOTERR(9:16) = 'ELTFA ' INTERR(1) = 1 CALL ERREUR(698) SEGDES MRIGID GOTO 100 ENDIF ELSE IF(NRIGEL.NE.NBSOUS)THEN MOTERR(1:8) = 'DARCY ' MOTERR(9:16) = 'ELTFA ' INTERR(1) = 1 CALL ERREUR(698) SEGDES MRIGID GOTO 100 ENDIF DO 10 ISOUS=1,NBSOUS IF (LISOUS(ISOUS).NE.IRIGEL(1,ISOUS)) THEN MOTERR(1:8) = 'DARCY ' MOTERR(9:16) = 'ELTFA ' INTERR(1) = ISOUS CALL ERREUR(698) SEGDES MRIGID GOTO 100 ENDIF 10 CONTINUE ENDIF * * LECTURE DES FLUX OU * DU CHPOIN FACE A MULTIPLIER * CALL LIROBJ('CHPOINT',IPCHF,1,IRET) IF(IRET.NE.1) RETURN MCHPO1=IPCHF NOMTOT(1)=' ' INDIC=0 NBCOMP=-1 CALL QUEPOI(MCHPO1,IFACE,INDIC,NBCOMP,NOMTOT) IF(IERR.NE.0)RETURN SEGACT MCHPO1 MSOUP1=MCHPO1.IPCHP(1) SEGACT MSOUP1 IF(MSOUP1.NOCOMP(1).EQ.'FLUX')THEN MPOVA1=MSOUP1.IPOVAL SEGACT MPOVA1 * * CALL LIROBJ('CHPOINT',IPCHF,1,IRET) IF(IRET.NE.1) RETURN MCHPO2=IPCHF NOMTOT(1)=' ' INDIC=0 NBCOMP=-1 CALL QUEPOI(MCHPO2,IFACE,INDIC,NBCOMP,NOMTOT) IF(IERR.NE.0)RETURN SEGACT MCHPO2 MSOUP2=MCHPO2.IPCHP(1) SEGACT MSOUP2 MPOVA2=MSOUP2.IPOVAL SEGACT MPOVA2 ELSE MCHPO2=MCHPO1 MCHPO1=0 MSOUP2=MSOUP1 MPOVA2=MSOUP2.IPOVAL SEGACT MPOVA2 CALL LIROBJ('CHPOINT',IPCHF,0,IRET) IF(IRET.EQ.1) THEN MCHPO1=IPCHF NOMTOT(1)=' ' INDIC=0 NBCOMP=-1 CALL QUEPOI(MCHPO1,IFACE,INDIC,NBCOMP,NOMTOT) IF(IERR.NE.0)RETURN SEGACT MCHPO1 MSOUP1=MCHPO1.IPCHP(1) SEGACT MSOUP1 IF(MSOUP1.NOCOMP(1).EQ.'FLUX')THEN MPOVA1=MSOUP1.IPOVAL SEGACT MPOVA1 ELSE CALL ERREUR(21) RETURN ENDIF ENDIF ENDIF * * Construction de MCHPOI * * IPT1=ICENTR SEGACT IPT1 NPN=IPT1.NUM(/2) NSOUPO=1 NAT=1 SEGINI MCHPOI MTYPOI=' ' MOCHDE=' CHPOIN CREE PAR DMMU ' IFOPOI=IFOUR JATTRI(1)=2 NC=NBCOMP SEGINI MSOUPO IPCHP(1)=MSOUPO DO 5 L=1,NBCOMP NOCOMP(L)=MSOUP2.NOCOMP(L) NOHARM(L)=MSOUP2.NOHARM(L) 5 CONTINUE IGEOC=ICENTR N=NPN SEGINI MPOVAL IPOVAL=MPOVAL NB=N*NC CALL INITD(VPOCHA,NB,0.D0) * * Creation du tableau ICPR * IK = 0 NNGOT = XCOOR(/1)/(IDIM+1) SEGINI ICCPR MELEME = IFACE SEGACT MELEME N2 = NUM(/2) DO 15 I2=1,N2 K = NUM(1,I2) IF (ICPR(K).EQ.0) THEN IK = IK + 1 ICPR(K) = IK ENDIF 15 CONTINUE ITES=IK SEGDES MELEME SEGINI IORGA MELEME = MSOUP2.IGEOC SEGACT MELEME N2 = NUM(/2) DO 25 I2=1,N2 K = NUM(1,I2) IF (ICPR(K).EQ.0) THEN INTERR(1) = K MOTERR(1:8) = 'FACE ' CALL ERREUR(64) SEGSUP ICCPR, IORGA RETURN ELSE IAA(ICPR(K)) = I2 ENDIF 25 CONTINUE SEGDES MELEME CALL LEKTAB(IPTABL,'XXNORMAE',IORIE) MCHELM=IORIE SEGACT MCHELM C C Calcul du produit C MELEME=IELTFA SEGACT MELEME IPT1=MELEME ITELEM=0 IF(MCHPO1.NE.0)THEN C C Cas avec vitesses C DO 50 ISOUS=1,NRIGEL MCHAML= ICHAML(ISOUS) SEGACT MCHAML MELVAL=IELVAL(1) SEGACT MELVAL IF(NRIGEL.NE.1)IPT1=LISOUS(ISOUS) SEGACT IPT1 xMATRI=IRIGEL(4,ISOUS) SEGACT xMATRI NELRIG=re(/3) C write(6,*)' ISOUS NELRIG ',ISOUS,NELRIG DO 60 IEL=1,NELRIG ITELEM=ITELEM+1 * XMATRI=IMATTT(IEL) * SEGACT XMATRI NLIGRD=RE(/1) NLIGRP=RE(/2) C write(6,*) ' NLIGRD NLIGRP ',NLIGRD,NLIGRP DO 40 I=1,NLIGRD RLIGN=0.D0 DO 30 J=1,NLIGRP RLIGN=RLIGN+RE(I,J,iel) 30 CONTINUE IPOPTS=ICPR(IPT1.NUM(I,IEL)) III=IAA(IPOPTS) C write(6,*)' dmmu ',itelem,i,rlign,MPOVA1.VPOCHA(IPOPTS,1) C * , VELCHE(I,IEL) RLIGN=RLIGN-MPOVA1.VPOCHA(IPOPTS,1)*VELCHE(I,IEL) DO 20 K=1,NBCOMP VPOCHA(ITELEM,K)=VPOCHA(ITELEM,K)+RLIGN* * MPOVA2.VPOCHA(III,K) 20 CONTINUE 40 CONTINUE * SEGDES XMATRI 60 CONTINUE SEGDES xMATRI,MCHAML,MELVAL 50 CONTINUE ELSE C C Cas ou il n'y a pas de vitesse C On a dupliqué toute la boucle 50 pour ne pas faire C NRIGEL*NELRIG*NLIGRD fois le meme test C DO 51 ISOUS=1,NRIGEL MCHAML= ICHAML(ISOUS) SEGACT MCHAML MELVAL=IELVAL(1) SEGACT MELVAL IF(NRIGEL.NE.1)IPT1=LISOUS(ISOUS) SEGACT IPT1 xMATRI=IRIGEL(4,ISOUS) SEGACT xMATRI NELRIG=re(/3) DO 61 IEL=1,NELRIG ITELEM=ITELEM+1 * XMATRI=IMATTT(IEL) * SEGACT XMATRI NLIGRD=RE(/1) NLIGRP=RE(/2) DO 41 I=1,NLIGRD RLIGN=0.D0 DO 31 J=1,NLIGRP RLIGN=RLIGN+RE(I,J,iel) 31 CONTINUE IPOPTS=ICPR(IPT1.NUM(I,IEL)) III=IAA(IPOPTS) DO 21 K=1,NBCOMP VPOCHA(ITELEM,K)=VPOCHA(ITELEM,K)+RLIGN* * MPOVA2.VPOCHA(III,K) 21 CONTINUE 41 CONTINUE * SEGDES XMATRI 61 CONTINUE SEGDES xMATRI,MCHAML,MELVAL 51 CONTINUE ENDIF SEGDES MRIGID CALL ECROBJ('CHPOINT',MCHPOI) * * Ménage * SEGDES MELEME SEGSUP ICCPR ,IORGA SEGDES MPOVAL,MSOUPO,MCHPOI SEGDES MPOVA2,MSOUP2,MCHPO2 IF(MCHPO1.NE.0) THEN SEGDES MPOVA1,MSOUP1,MCHPO1 ENDIF 100 CONTINUE RETURN END