C YDUDW SOURCE GOUNAND 23/08/22 21:15:04 11725 SUBROUTINE YDUDW IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C*********************************************************************** C C CET OPERATEUR DISCRETISE L'OPERATEUR DE PENALISATION DIV(U)=EPS*P C EN 2D SUR LES ELEMENTS QUA4 TRI3 TRI7 et QUA9 PLAN OU AXI C EN 3D SUR LES ELEMENTS CUB8 ET PRI6 C C SYNTAXE : C --------- C C DUDW(EPS) INCO UN : C C COEFFICIENTS : C -------------- C C C EPS (SCAL DOMA) PARAMETRE DE PENALISATION C (SCAL ELEM) C C INCONNUES : C ----------- C C UN CHAMP DE VITESSE C C C OPTIONS : POROSITE DIV(PU)=EPS*P C SOURCE DIV(U)-Q=EPS*P C C************************************************************************ -INC PPARAM -INC CCOPTIO -INC CCGEOME -INC SMCOORD -INC SMLENTI -INC SMELEME POINTEUR MELEM1.MELEME,MELEMS.MELEME,MELEMI.MELEME,MELEMP.MELEME POINTEUR MELEMC.MELEME -INC SMCHPOI POINTEUR IZGG1.MPOVAL,SRCE.MPOVAL,IZTGG1.MPOVAL -INC SIZFFB POINTEUR IZF1.IZFFM -INC SMMATRIK -INC SMLMOTS POINTEUR LINCO.MLMOTS SEGMENT TRAV REAL*8 AF(NP,NP,NINC,NINC),AS(NP,NINC),CT(MP1,NP,NINC),PQ(MP1) ENDSEGMENT CHARACTER*8 NOMZ,TYPE,TYPC CHARACTER*(LOCOMP) NOM0,NOM1,NOM2,NOM3,NOMI DIMENSION IXV(3) PARAMETER (NTB=1) CHARACTER*8 LTAB(NTB) DIMENSION KTAB(NTB) DATA LTAB/'KIZX '/ C***************************************************************************** CDUDW C write(6,*)' Operateur DUDW ' CALL LITABS(LTAB,KTAB,NTB,1,IRET) IF(IERR.NE.0)RETURN MTABX=KTAB(1) C C- Récupération de la table EQEX (pointeur MTAB1) C CALL LEKTAB(MTABX,'EQEX',MTAB1) IF(MTAB1.EQ.0)THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' EQEX ' MOTERR( 9:16) = ' EQEX ' MOTERR(17:24) = ' KIZX ' CALL ERREUR(786) RETURN ENDIF CALL ACME(MTAB1,'NAVISTOK',NASTOK) IF(NASTOK.EQ.0)THEN CALL ZDUDW(MTABX,MTAB1) RETURN ENDIF C C- Récupération de la table INCO (pointeur KINC) C CALL LEKTAB(MTAB1,'INCO',KINC) IF(KINC.EQ.0)THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' INCO ' MOTERR( 9:16) = ' INCO ' MOTERR(17:24) = ' EQEX ' CALL ERREUR(786) RETURN ENDIF C***************************************************************************** C OPTIONS C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> SEMI C KFORM = 0 -> SI 1 -> EF 2 -> VF 3 -> EFMC C IDCEN = 0-> rien 1-> CENTREE 2-> SUPGDC 3-> SUPG 4-> TVISQUEU 5-> CNG IAXI=0 IF(IFOMOD.EQ.0)IAXI=2 C C- Récupération de la table des options KOPT (pointeur KOPTI) C CALL LEKTAB(MTABX,'KOPT',KOPTI) IF (KOPTI.EQ.0) THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' KOPT ' MOTERR( 9:16) = ' KOPT ' MOTERR(17:24) = ' KIZX ' CALL ERREUR(786) RETURN ENDIF CALL ACME(KOPTI,'KIMPL',KIMPL) CALL ACME(KOPTI,'KPOIN',KPRE) CALL ACME(KOPTI,'KFORM',KFORM) IF(KFORM.NE.0.AND.KFORM.NE.1)THEN C Option %m1:8 incompatible avec les données MOTERR( 1: 8) = 'EF/EFM1 ' CALL ERREUR(803) RETURN ENDIF C write(6,*)' Apres les options ' C***************************************************************************** C C- Récupération de la table DOMAINE associée au domaine local C CALL ACMM(MTABX,'NOMZONE',NOMZ) CALL LEKTAB(MTABX,'DOMZ',MTABZ) IF(MTABZ.EQ.0)THEN C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24 MOTERR( 1: 8) = ' DOMZ ' MOTERR( 9:16) = ' DOMZ ' MOTERR(17:24) = ' KIZX ' CALL ERREUR(786) RETURN ENDIF CALL LEKTAB(MTABZ,'SOMMET',MELEMS) CALL LEKTAB(MTABZ,'MAILLAGE',MELEME) CALL LEKTAB(MTABZ,'MACRO1',MACRO1) MELEMI=MACRO1 CALL LEKTAB(MTABZ,'QUADRATI',MQUAD) IF (IERR.NE.0) RETURN IF(MQUAD.EQ.0.AND.MACRO1.EQ.0.AND.KPRE.NE.2)THEN WRITE(6,*)'Operateur DUDW ' WRITE(6,*)'Type d''éléments non prévu' RETURN ENDIF IF(KPRE.EQ.2)THEN CALL LEKTAB(MTABZ,'CENTRE',MELEMC) MELEMP=MELEMC MELEMI=MELEME ELSEIF(KPRE.EQ.3)THEN CALL LEKTAB(MTABZ,'CENTREP0',MELEMC) MELEMP=MELEMC ELSEIF(KPRE.EQ.4)THEN CALL LEKTAB(MTABZ,'CENTREP1',MELEMC) CALL LEKTAB(MTABZ,'ELTP1NC ',MELEMP) ENDIF SEGACT MELEME SEGACT MELEMP C************************************************************************* C VERIFICATIONS SUR LES INCONNUES C write(6,*)' Verification des inconnues ' TYPE='LISTMOTS' CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO) SEGACT LINCO NBINC=LINCO.MOTS(/2) IF(NBINC.NE.1)THEN C Indice %m1:8 : contient plus de %i1 %m9:16 MOTERR( 1:8) = 'LISTINCO' INTERR(1) = 1 MOTERR(9:16) = ' MOTS ' CALL ERREUR(799) RETURN ENDIF NOMI=LINCO.MOTS(1) TYPE=' ' CALL ACMO(KINC,NOMI,TYPE,MCHPOI) IF(TYPE.NE.'CHPOINT ')THEN C Indice %m1:8 : ne contient pas un objet de type %m9:16 MOTERR( 1: 8) = 'INC '//NOMI MOTERR( 9:16) = 'CHPOINT ' CALL ERREUR(800) RETURN ELSE CALL LICHT(MCHPOI,IZTU1,TYPC,MELEM1) ENDIF C************************************************************************* C Le domaine de definition est donne par le SPG de la premiere inconnue C Les inconnues suivantes devront posseder ce meme pointeur C On verifie que les points de la zone sont tous inclus dans ce SPG CALL KRIPAD(MELEM1,MLENTI) CALL VERPAD(MLENTI,MELEME,IRET) IF(IRET.NE.0)THEN C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique MOTERR(1: 8) = 'INC '//NOMI MOTERR(9:16) = 'CHPOINT ' CALL ERREUR(788) RETURN ENDIF NINC=IDIM MLENT1=MLENTI C***************************************************************************** C Lecture du ou des coefficients C write(6,*)' Lecture des coefficients' CALL ACME(MTABX,'IARG',IARG) C? IXV(1)=MELEMC IXV(1)=0 IXV(2)=1 IXV(3)=0 CALL LEKCOF('Opérateur DUDW :', & MTABX,KINC,1,IXV,IZTG1,IZTGG1,NPT1,NC1,IK1,IRET) IF(IRET.EQ.0)RETURN KITT=2 INDGS=0 SRCE=IZTGG1 MPOVA1 =IZTGG1 MPOVA3 =IZTGG1 IF(IARG.EQ.2)THEN IXV(1)=MELEMC IXV(2)=1 IXV(3)=0 CALL LEKCOF('Opérateur DUDW :', & MTABX,KINC,2,IXV,IZS,SRCE,NPTS,NCS,IKS,IRET) IF(IRET.EQ.0)RETURN INDGS=1 CALL KRIPAD(MELEMC,MLENT1) ENDIF C************************************************************************* IF(KFORM.EQ.0)THEN C CAS FORMULATION EF SI (GRESHO) IF(KIMPL.NE.0)THEN GO TO 90 ENDIF WRITE(6,*)' Operateur DUDW ' WRITE(6,*)' Cas Formulation EF SI ' WRITE(6,*)' Cas invalide ' GO TO 90 C*********************************************************************** ELSE IF(KFORM.EQ.1)THEN C CAS FORMULATION EF NUTOEL=0 SEGACT MELEMI NBSOUS=MELEMI.LISOUS(/1) IF(NBSOUS.EQ.0)NBSOUS=1 NRIGE=7 NKID =9 NKMT =7 NMATRI=1 SEGINI MATRIK IRIGEL(1,1)=MELEMI IRIGEL(2,1)=MELEMI NBOP=0 NBME=NINC*NINC C NBME=2 NBELC=0 SEGINI IMATRI IRIGEL(4,1)=IMATRI KSPGP=MELEMS KSPGD=MELEMS IF(NBME.EQ.4)THEN NOM1='1'//NOMI(1:7) NOM2='2'//NOMI(1:7) LISPRI(1)=NOM1 LISDUA(1)=NOM1 LISPRI(2)=NOM2 LISDUA(2)=NOM1 LISPRI(3)=NOM1 LISDUA(3)=NOM2 LISPRI(4)=NOM2 LISDUA(4)=NOM2 ELSEIF(NBME.EQ.9)THEN NOM1='1'//NOMI(1:7) NOM2='2'//NOMI(1:7) NOM3='3'//NOMI(1:7) LISPRI(1)=NOM1 LISDUA(1)=NOM1 LISPRI(2)=NOM2 LISDUA(2)=NOM1 LISPRI(3)=NOM3 LISDUA(3)=NOM1 LISPRI(4)=NOM1 LISDUA(4)=NOM2 LISPRI(5)=NOM2 LISDUA(5)=NOM2 LISPRI(6)=NOM3 LISDUA(6)=NOM2 LISPRI(7)=NOM1 LISDUA(7)=NOM3 LISPRI(8)=NOM2 LISDUA(8)=NOM3 LISPRI(9)=NOM3 LISDUA(9)=NOM3 ELSEIF(NBME.EQ.2)THEN LISPRI(1)='1'//NOMI(1:7) LISDUA(1)='1'//NOMI(1:7) LISPRI(2)='2'//NOMI(1:7) LISDUA(2)='2'//NOMI(1:7) ELSE WRITE(6,*)' Operateur DUDW ' WRITE(6,*)' Cas invalide' GO TO 90 ENDIF IF(INDGS.NE.0)THEN NAT=2 NSOUPO=1 SEGACT MELEMS N=MELEMS.NUM(/2) NC=NINC SEGINI MCHPO1,MSOUP1,MPOVA1 MCHPO1.IFOPOI=IFOUR MCHPO1.MOCHDE=TITREE MCHPO1.MTYPOI='SMBR' MCHPO1.JATTRI(1)=2 MCHPO1.IPCHP(1)=MSOUP1 MSOUP1.NOCOMP(1)=NOM1 MSOUP1.NOCOMP(2)=NOM2 IF(NINC.EQ.3)MSOUP1.NOCOMP(3)=NOM3 MSOUP1.IGEOC=MELEMS MSOUP1.IPOVAL=MPOVA1 SEGACT MELEMC N=MELEMC.NUM(/2) NC=1 SEGINI MCHPO3,MSOUP3,MPOVA3 MCHPO3.IFOPOI=IFOUR MCHPO3.MOCHDE=TITREE MCHPO3.MTYPOI='SMBR' MCHPO3.JATTRI(1)=2 MCHPO3.IPCHP(1)=MSOUP3 MSOUP3.NOCOMP(1)='SCAL' MSOUP3.IGEOC=MELEMC MSOUP3.IPOVAL=MPOVA3 ENDIF DO 101 L=1,NBSOUS IPT1=MELEMI IPT2=MELEMI IF(NBSOUS.NE.1)IPT1=MELEMI.LISOUS(L) SEGACT IPT1 IF(INDGS.NE.0)THEN IPT2=MELEMP IF(NBSOUS.NE.1)IPT2=MELEMP.LISOUS(L) SEGACT IPT2 ENDIF IF(MQUAD.NE.0)THEN IF(KPRE.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//'PRP0' IF(KPRE.EQ.3)NOM0=NOMS(IPT1.ITYPEL)//'PRP0' IF(KPRE.EQ.4)NOM0=NOMS(IPT1.ITYPEL)//'PRP1' ELSEIF(MACRO1.NE.0)THEN IF(KPRE.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//' ' IF(KPRE.EQ.3)NOM0=NOMS(IPT1.ITYPEL)//'MCP0' IF(KPRE.EQ.4)NOM0=NOMS(IPT1.ITYPEL)//'MCP1' ELSE IF(KPRE.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//' ' ENDIF CALL KALPBG(NOM0,'FONFORM ',IZFFM) SEGACT IZFFM*MOD IZHR=KZHR(1) SEGACT IZHR*MOD NES=GR(/1) NPG=GR(/3) IZF1=KTP(1) SEGACT IZF1*MOD MP1=IZF1.FN(/1) NP = IPT1.NUM(/1) MP = NP NBEL=IPT1.NUM(/2) SEGINI IPM1,IPM2,IPM3,IPM4 LIZAFM(L,1)=IPM1 C LIZAFM(L,2)=IPM4 LIZAFM(L,2)=IPM2 LIZAFM(L,3)=IPM3 LIZAFM(L,4)=IPM4 IPM5=IPM1 IPM6=IPM1 IPM7=IPM1 IPM8=IPM1 IPM9=IPM1 IF(NBME.EQ.9)THEN SEGINI IPM5,IPM6,IPM7,IPM8,IPM9 LIZAFM(L,5)=IPM5 LIZAFM(L,6)=IPM6 LIZAFM(L,7)=IPM7 LIZAFM(L,8)=IPM8 LIZAFM(L,9)=IPM9 ENDIF C Pour l'instant les sources et puits de masse ne sont pas actifs INDG2=0 KI2=0 KJ2=0 SEGINI TRAV NPT=MPOVA1.VPOCHA(/1) SEGACT,MCOORD CALL XDUDW(FN,IZF1.FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ, & NES,IDIM,NP,MP1,NPG,IAXI,NINC, & IZTGG1.VPOCHA,IK1,SRCE.VPOCHA,INDGS,IKS, & IPT1.NUM,NBEL,NUTOEL,XCOOR,AF,AS,CT,PQ, & IPM1.AM,IPM2.AM,IPM3.AM,IPM4.AM,IPM5.AM,IPM6.AM,IPM7.AM, & IPM8.AM,IPM9.AM,MPOVA1.VPOCHA,NPT,LECT,IPT2.NUM,MLENT1.LECT, & MPOVA3.VPOCHA) SEGDES,MCOORD SEGSUP TRAV NUTOEL=NUTOEL+NBEL 101 CONTINUE SEGDES IMATRI,MATRIK CALL ECROBJ('MATRIK',MATRIK) C? CALL ECMO(MTABX,'MATELM','MATRIK',MATRIK) NAT=2 NSOUPO=0 SEGINI MCHPOI JATTRI(1)=2 CALL ECROBJ('CHPOINT',MCHPOI) IF(INDGS.NE.0)THEN TYPE=' ' CALL ACMO(MTAB1,'SMBR',TYPE,MCHPO2) IF(TYPE.NE.'CHPOINT')THEN CALL ECMO(MTAB1,'SMBR','CHPOINT',MCHPO1) ELSE CALL ECROBJ('CHPOINT',MCHPO2) CALL ECROBJ('CHPOINT',MCHPO1) CALL PRFUSE CALL LIROBJ('CHPOINT',MCHPOI,1,IRET) CALL ECMO(MTAB1,'SMBR','CHPOINT',MCHPOI) ENDIF TYPE=' ' C CALL ACMO(MTAB1,'DUDWSRCE',TYPE,MCHPO4) C IF(TYPE.NE.'CHPOINT')THEN CALL ECMO(MTAB1,'DUDWSRCE','CHPOINT',MCHPO3) C ELSE C CALL ECROBJ('CHPOINT',MCHPO3) C CALL ECROBJ('CHPOINT',MCHPO4) C CALL PRFUSE C CALL LIROBJ('CHPOINT',MCHPOI,1,IRET) C CALL ECMO(MTAB1,'DUDWSRCE','CHPOINT',MCHPOI) C ENDIF ENDIF SEGSUP MLENTI IF(INDGS.NE.0)SEGSUP MLENT1 RETURN C************************************************************************* ELSE IF(KFORM.EQ.2)THEN C CAS FORMULATION VF RETURN ENDIF C************************************************************************* 1002 FORMAT(10(1X,1PE11.4)) 90 CONTINUE * WRITE(6,*)' Interuption anormale de DUDW' * 34 2 * Operateur incompatible avec les options utilisees CALL ERREUR(34) RETURN END