C SQTP SOURCE CHAT 11/03/16 21:31:46 6902 SUBROUTINE SQTP C----------------------------------------------------------------------- C Calcul de la contribution au systeme en trace de charge d'une force C volumique dans le cas de la résolution des équations de Darcy par EFMH C avec le modèle DARCY. C----------------------------------------------------------------------- C C--------------------------- C Phrase d'appel (GIBIANE) : C--------------------------- C C CHP2 = 'SQTP' MMODEL RIG1 RIG2 CHP1 ; C C------------------------ C Operandes et resultat : C------------------------ C C MMODEL : Objet modele décrivant la formulation. C TABLE1 : TABLE DOMAINE contenant les maillages et les connectivités. C RIG1 : Matrices masses hybrides elementaires de sous-type MASSE. C RIG2 : Matrices elementaires de type HYBTP C CHP1 : CHPO centre de composante FX,FY(,FZ), C densite de force moyenne par élément. C C----------------------------------------------------------------------- C C Langage : ESOPE + FORTRAN77 C C Auteurs : 02/96 L.V.BENET - Cas permanent C C----------------------------------------------------------------------- IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) * -INC CCOPTIO -INC SMCHAML -INC SMCHPOI -INC SMELEME -INC SMMODEL -INC SMRIGID -INC SMTABLE * SEGMENT IPMAHY INTEGER MAHYBR(NSOUS) * CHARACTER*4 NOMTOT(IDIM) ENDSEGMENT * LOGICAL LOGRE,LOGIN CHARACTER*4 NOMTO3(3) CHARACTER*4 NOMTO2(2) CHARACTER*6 CHAR6 CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE,LETYPE,LETYP1,LETYP2 * * Initialisations * ISOU1 = 0 IHN1 = 0 ITPN1 = 0 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 * * Lecture de la TABLE domaine * IPTABL = 0 CALL LEKMOD(IPMODE,IPTABL,IRET) IF (IERR.NE.0) RETURN CHARIN = 'MAILLAGE' TYPOBJ = 'MAILLAGE' CALL LEKTAB(IPTABL,CHARIN,IOBRE) IF (IERR.NE.0) RETURN IPGEOM = IOBRE 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 IPFACE = IOBRE * * Lecture de RIGIDITE * CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET) IF (IERR.NE.0) RETURN RI1 = IPRIGI * * Lecture de RIGIDITE * CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET) IF (IERR.NE.0) RETURN RI2 = IPRIGI * * Lecture du champoint FORCE * CALL LIROBJ('CHPOINT',IFORC,1,IRET) IF (IERR.NE.0) RETURN * * récup. MCHAML orientation normale * CALL LEKTAB(IPTABL,'XXNORMAE',IORIE) IF (IERR.NE.0) RETURN * * récup. CHPO orientation normale * CALL LEKTAB(IPTABL,'XXNORMAF',INORM) IF (IERR.NE.0) RETURN * * récup. CHPO surface des faces * CALL LEKTAB(IPTABL,'XXSURFAC',ISURF) IF (IERR.NE.0) RETURN * * Test du CHAMPOINT FORCE de spg CENTRE * NBCOMP = IDIM INDIC = 1 IF(IDIM.EQ.2)THEN NOMTO2(1) = 'FX' NOMTO2(2) = 'FY' CALL QUEPOI(IFORC,ICENTR,INDIC,NBCOMP,NOMTO2) IF (IERR.NE.0) RETURN ELSE IF(IDIM.EQ.3)THEN NOMTO3(1) = 'FX' NOMTO3(2) = 'FY' NOMTO3(3) = 'FZ' CALL QUEPOI(IFORC,ICENTR,INDIC,NBCOMP,NOMTO3) IF (IERR.NE.0) RETURN ENDIF * * Test de la formulation * SEGACT MMODEL NSOUS = KMODEL(/1) SEGINI IPMAHY IDARCY = 0 IPT1=IPGEOM IPT2=IPGEOM DO 10 ISOUS=1,NSOUS IF(NSOUS.GT.1)THEN SEGACT IPT2 IPT1=IPT2.LISOUS(ISOUS) SEGDES IPT2 ENDIF IMODEL = KMODEL(ISOUS) SEGACT IMODEL LETYPE = FORMOD(1) IF (LETYPE.EQ.'DARCY') THEN IDARCY = IDARCY + 1 MAHYBR(ISOUS) = IPT1 ENDIF SEGDES IMODEL 10 CONTINUE SEGDES MMODEL IF (IDARCY.EQ.0) THEN MOTERR = LETYPE CALL ERREUR(193) GOTO 100 ENDIF * * Recuperation des pointeurs ELTFA pour les zones ou DARCY est defini * MELEME = IELTFA SEGACT MELEME LZONES = LISOUS(/1) IF (LZONES.EQ.0) LZONES = 1 IPT1 = IPGEOM SEGACT IPT1 DO 30 ISOUS=1,NSOUS IMAMEL = MAHYBR(ISOUS) IF (IMAMEL.NE.0) THEN DO 20 ISZ=1,LZONES IF (LZONES.EQ.1) THEN IPT2 = IPT1 IPT3 = MELEME ELSE IPT2 = IPT1.LISOUS(ISZ) IPT3 = LISOUS(ISZ) ENDIF IF (IPT2.EQ.IMAMEL) THEN MAHYBR(ISOUS) = IPT3 GOTO 30 ENDIF 20 CONTINUE IF (IMAMEL.EQ.MAHYBR(ISOUS)) THEN INTERR(1) = ISOUS CALL ERREUR(698) SEGDES IPT1 SEGDES MELEME GOTO 100 ENDIF ENDIF 30 CONTINUE SEGDES IPT1 SEGDES MELEME * * Test du sous-type des matrices de rigiditée MASSE et HYBTP récupérées * SEGACT RI1 LETYP1 = RI1.MTYMAT SEGACT RI2 LETYP2 = RI2.MTYMAT IF (LETYP1.NE.'MASSE') THEN IF (LETYP2.NE.'MASSE') THEN MOTERR(1:8) = 'RIGIDITE' MOTERR(9:16) = 'MASSE' CALL ERREUR(79) SEGDES RI1 SEGDES RI2 GOTO 100 ELSE IPFORC=RI2 IF (LETYP1.NE.'HYBTP') THEN MOTERR(1:8) = 'RIGIDITE' MOTERR(9:16) = 'HYBTP' CALL ERREUR(79) SEGDES RI1 SEGDES RI2 GOTO 100 ELSE IPMATP=RI1 ENDIF ENDIF ELSE IPFORC=RI1 IF (LETYP2.NE.'HYBTP') THEN MOTERR(1:8) = 'RIGIDITE' MOTERR(9:16) = 'HYBTP' CALL ERREUR(79) SEGDES RI1 SEGDES RI2 GOTO 100 ELSE IPMATP=RI2 ENDIF ENDIF * * Controle des pointeurs de MELEME des deux rigidités * DO 40 ISOUS=1,NSOUS IMAMEL = MAHYBR(ISOUS) IF (IMAMEL.NE.0) THEN IPTEST = RI1.IRIGEL(1,ISOUS) IF (MAHYBR(ISOUS).NE.IPTEST) THEN MOTERR(1:8) = LETYP1 MOTERR(9:16) = 'ELTFA ' INTERR(1) = ISOUS CALL ERREUR(698) SEGDES RI1 GOTO 100 ENDIF IPTEST = RI2.IRIGEL(1,ISOUS) IF (MAHYBR(ISOUS).NE.IPTEST) THEN MOTERR(1:8) = LETYP2 MOTERR(9:16) = 'ELTFA ' INTERR(1) = ISOUS CALL ERREUR(698) SEGDES RI2 GOTO 100 ENDIF ENDIF 40 CONTINUE SEGDES RI1 SEGDES RI2 * * Construction de MCHPOI * SEGDES IPMAHY CALL SQTP1(IPMAHY,IPFACE,IPFORC,IPMATP,IFORC, S IORIE,INORM,ISURF,MCHPOI) * CALL ECROBJ('CHPOINT',MCHPOI) * * Ménage * 100 CONTINUE SEGSUP IPMAHY RETURN END