Sources Esope | Cast3M

zdudw
C ZDUDW     SOURCE    CB215821  23/01/25    21:15:41     11573          
      SUBROUTINE ZDUDW(MTABX,MTAB1)
      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 TYPE,TYPC
      CHARACTER*(LOCOMP) NOM0,NOM1,NOM2,NOM3,NOMZ,NOMI
      DIMENSION IXV(3)
      PARAMETER (NTB=1)
      DIMENSION KTAB(NTB)
C*****************************************************************************
CDUDW
C     write(6,*)' Operateur DUDW '
 
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)
      MACRO1=0
      CALL LEKTAB(MTABZ,'MACRO1',MELEMI)
      IF(MELEMI.NE.0)MACRO1=1
      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= '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= '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:LOCOMP-1)
      NOM2='2'//NOMI(1:LOCOMP-1)
      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:LOCOMP-1)
      NOM2='2'//NOMI(1:LOCOMP-1)
      NOM3='3'//NOMI(1:LOCOMP-1)
      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:LOCOMP-1)
      LISDUA(1)='1'//NOMI(1:LOCOMP-1)
      LISPRI(2)='2'//NOMI(1:LOCOMP-1)
      LISDUA(2)='2'//NOMI(1:LOCOMP-1)
 
      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)
 
      SEGSUP TRAV
      SEGDES IPT1
      NUTOEL=NUTOEL+NBEL
 101  CONTINUE
      SEGDES IMATRI
      SEGDES MELEME,MATRIK
      IF(IK1.EQ.1)THEN
      SEGSUP IZTGG1
      ELSE
      SEGDES IZTGG1
      ENDIF
      CALL ECMO(MTABX,'MATELM','MATRIK',MATRIK)
 
 
      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'
      RETURN
      END