Télécharger zzfimp.eso

Retour à la liste

zzfimp
C ZZFIMP    SOURCE    GOUNAND   25/11/12    21:16:04     12399          
      SUBROUTINE ZZFIMP(MTABX,MTAB1)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C
C     SYNTAXE :
C
C        FIMP  coef
C
C       EN 2D
C              elements SEG2  -> Flux
C              elements TRI3  -> Source volumique
C              elements QUA4  -> Source volumique
C       EN 3D
C              elements SEG2  -> Pas de sens !!
C              elements TRI3  -> Flux
C              elements QUA4  -> Flux
C              elements CUB8  -> Source volumique
C              elements PRI6  -> Source volumique
C              elements TET4  -> Source volumique
C
C
C  MTAB1 : Table type EQEX        -> RV
C  MTABZ : Table type DOMAINE     -> Zone definition opérateur
C  MTABD : Table type DOMAINE     -> Zone Totale apres assemblage
C  MTABX : Table type KIZX        -> Description opérateur
C
C  les cartes correspondantes sont commentées C¤.
C  Gare à l'explicite !!!
C
C
C
C
C***********************************************************************
 
 
-INC PPARAM
-INC CCOPTIO
-INC CCGEOME
-INC SMCHAML
-INC SMCOORD
-INC SMLENTI
-INC SMELEME
      POINTEUR MELEM1.MELEME,SPGID.MELEME,SPGZ.MELEME
      POINTEUR MELEMP.MELEME
-INC SMCHPOI
      POINTEUR IZG1.MCHPOI, IZGG1.MPOVAL
      POINTEUR IZTU1.MPOVAL
      POINTEUR MZFLU.MPOVAL
      POINTEUR IZVOL.MPOVAL
-INC SIZFFB
      POINTEUR IZF1.IZFFM
 
-INC SMLMOTS
      POINTEUR LINCO.MLMOTS
      CHARACTER*8 TYPE,TYPC,MTERR
      CHARACTER*(LOCOMP) NOMZ,NOMI,NOM0
      LOGICAL LOGI
      PARAMETER (NTB=1)
      CHARACTER*8 LTAB(NTB)
      DIMENSION KTAB(NTB),IXV(4)
      SAVE IPAS
      DATA LTAB/'KIZX    '/,IPAS/0/
C*****************************************************************************
CFIMP
C?    write(6,*)' Debut FIMP'
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 KFORM = 0 -> SI       1 -> EF       2 -> VF    3 -> EFMC
C IDCEN = 0->rien   1-> CENTREE  2-> SUPGDC  3-> SUPG     4-> TVISQUEU 5-> CNG
C KPOIN = 0->SOMMET 1-> FACE     2-> CENTRE  3-> CENTREP0 4-> CENTREP1
 
      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,'MTRMASS ',MMPG)
      IPG=0
      IF(MMPG.EQ.3)IPG=1
      CALL ACME(KOPTI,'IDCEN',IDCEN)
      KDIM=1
      IF(IDCEN.EQ.2)KDIM=IDIM
      CALL ACME(KOPTI,'IKOMP',IKOMP)
      CALL ACME(KOPTI,'KIMPL',KIMPL)
      CALL ACME(KOPTI,'KFORM',KFORM)
      CALL ACME(KOPTI,'KPOIND',KPOIND)
C     write(6,*)' INCOD=',KPOIND
 
      IF(KFORM.GE.2)THEN
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = 'EF/EFM1 '
            CALL ERREUR(803)
            RETURN
        ENDIF
      CALL ACME(KOPTI,'KMACO',KMACO)
      CALL ACMF(KOPTI,'AIMPL',AIMPL)
      IF (IERR.NE.0) RETURN
 
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
C
C- Récupération des indices CENTRE, FACE, SOMMET et MAILLAGE
C- de la table DOMAINE
C
C?    CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
      CALL LEKTAB(MTABZ,'MAILLAGE',MELEMZ)
      CALL LEKTAB(MTABZ,'SOMMET',MELEMS)
      CALL LEKTAB(MTABZ,'CENTRE',MELEMC)
      CALL LEKTAB(MTABZ,'XXPSOML',MCHELM)
      CALL LEKTAB(MTABZ,'XXVOLUM',MCHVOL)
      CALL LEKTAB(MTABZ,'MACRO',MACRO)
C?    IF(CALL LEKTAB(MTABZ,'MACRO1',MACRO1)
      CALL LEKTAB(MTABZ,'QUADRATI',MQUAD)
      IF (IERR.NE.0) RETURN
 
      CALL LICHT(MCHVOL,IZVOL,TYPC,IGEOM)
 
      SEGACT MCHELM
 
C
C- Vérification des compatiblités Formulation/SPG
C- Identification du spg de l'inconnue
C- SPGZ=spg inco duale de la zone; MELEME=connectivité associée ;
C
 
C     write(6,*)' KPOIND,KFORM=',KPOIND,KFORM
      IF(KFORM.EQ.0)THEN
        IF(KPOIND.EQ.99)THEN
        KPOIND=0
        SPGZ  =MELEMS
C       MELEME=MELEMS
        MELEME=MELEMZ
        ELSEIF(KPOIND.EQ.2)THEN
        SPGZ  =MELEMC
        MELEME=MELEMC
        MELEMP=MELEMC
        ELSEIF(KPOIND.EQ.0)THEN
        SPGZ  =MELEMS
C       MELEME=MELEMS
        MELEME=MELEMZ
        ELSEIF(KPOIND.NE.2.AND.KPOIND.NE.0)THEN
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = '  EF    '
            CALL ERREUR(803)
            RETURN
        ENDIF
 
      ELSEIF(KFORM.EQ.1)THEN
C       write(6,*)' KFORM=',kform,' KPOIND=',kpoind
        IF(KPOIND.EQ.99)THEN
        KPOIND=0
        SPGZ  =MELEMS
        MELEME=MELEMZ
        ELSEIF(KPOIND.EQ.0)THEN
        SPGZ  =MELEMS
        MELEME=MELEMZ
        ELSEIF(KPOIND.EQ.2)THEN
        SPGZ  =MELEMC
        MELEME=MELEMZ
        MELEMP=MELEMC
        ELSEIF(KPOIND.EQ.3)THEN
          MTERR='EF CTRP0'
          IF(MACRO.EQ.0.AND.MQUAD.EQ.0)GO TO 90
          CALL LEKTAB(MTABZ,'CENTREP0',MELEMC)
          SPGZ  =MELEMC
          MELEME=MELEMZ
          IF(MACRO.NE.0)CALL LEKTAB(MTABZ,'MACRO1',MELEME)
          MELEMP=MELEMC
        ELSEIF(KPOIND.EQ.4)THEN
          MTERR='EF CTRP1'
          IF(MACRO.EQ.0.AND.MQUAD.EQ.0)GO TO 90
          CALL LEKTAB(MTABZ,'CENTREP1',MELEMC)
          SPGZ  =MELEMC
          MELEME=MELEMZ
          IF(MACRO.NE.0)CALL LEKTAB(MTABZ,'MACRO1',MELEME)
          CALL LEKTAB(MTABZ,'ELTP1NC ',MELEMP)
        ELSEIF(KPOIND.NE.2.AND.KPOIND.NE.0.AND.KPOIND.NE.3
     &    .AND.KPOIND.NE.4)THEN
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = '  EF    '
            CALL ERREUR(803)
            RETURN
        ENDIF
 
      ELSEIF(KFORM.EQ.2)THEN
        IF(KPOIND.EQ.99)THEN
        KPOIND=2
        SPGZ  =MELEMC
        MELEME=MELEMC
        MELEMP=MELEMC
        ELSEIF(KPOIND.EQ.2)THEN
        SPGZ  =MELEMC
        MELEME=MELEMC
        MELEMP=MELEMC
        ELSEIF(KPOIND.NE.2)THEN
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = '  VF    '
            CALL ERREUR(803)
            RETURN
        ENDIF
 
      ELSEIF(KFORM.EQ.3)THEN
        MTERR='EFMC'
        GO TO 90
        IF(KPOIND.EQ.99.OR.KPOIND.EQ.2)THEN
        KPOIND=2
        SPGZ  =MELEMC
        MELEME=MELEMC
        ELSEIF(KPOIND.EQ.3.AND.(MACRO.NE.0.OR.MQUAD.NE.0))THEN
        CALL LEKTAB(MTABZ,'CENTREP0',MELEMC)
        SPGZ  =MELEMC
        MELEME=MELEMC
        ELSEIF(KPOIND.EQ.4.AND.(MACRO.NE.0.OR.MQUAD.NE.0))THEN
        CALL LEKTAB(MTABZ,'CENTREP1',MELEMC)
        CALL LEKTAB(MTABZ,'ELTP1NC',MELEMQ)
        SPGZ  =MELEMC
        MELEME=MELEMQ
C? SG 2025/11/10 MELEMO jamais utilise ?
C           IF(MACRO.NE.0)MELEMO=MACRO1
C           IF(MQUAD.NE.0)MELEMO=MELEMZ
C       ELSEIF(KPOIND.NE.2.AND.KPOIND.NE.3.AND.KPOIND.NE.4)THEN
        ELSE
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = '  EFMC  '
            CALL ERREUR(803)
            RETURN
        ENDIF
 
      ENDIF
C*************************************************************************
C VERIFICATIONS SUR LES INCONNUES
C
C- Récupération du nombre d'inconnues et du nom de l'inconnue NOMI
C
      TYPE='LISTMOTS'
      CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO)
      IF (IERR.NE.0) RETURN
      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)
 
      IF(KIMPL.EQ.0)THEN
      CALL LEKTAB(MTAB1,'KIZG',KIZG)
      IF(KIZG.EQ.0)THEN
      CALL CRTABL(KIZG)
      CALL ECMM(KIZG,'SOUSTYPE','KIZG')
      CALL ECMO(MTAB1,'KIZG','TABLE   ',KIZG)
      ENDIF
      ENDIF
 
C --> 1 ere Inconnue
 
      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,SPGID)
      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(SPGID,MLENTI)
      IF(IPAS.EQ.0)THEN
      CALL VERPAD(MLENTI,SPGZ,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)
      IPAS=0
      RETURN
      ENDIF
      ENDIF
 
C*************************************************************************
C Lecture du coefficient
C Type du coefficient :
C IK1=0 CHPOINT IK1=1 scalaire  IK1=2 vecteur
C     write(6,*)' Lecture des coefficients '
 
      CALL ACME(MTABX,'IARG',IARG)
      IF(IARG.NE.1)THEN
C           Indice %m1:8 : nombre d'arguments incorrect
            MOTERR(1:8) = 'IARG    '
            CALL ERREUR(804)
      RETURN
      ENDIF
 
      IRET=3
      IXV(1)=MELEMC
      IXV(2)=1
      IXV(3)=0
C?    IXV(4)=MELEMS
C     write(6,*)' MELEMQ=',melemq
      CALL LEKCOF('Opérateur FIMP :',
     & MTABX,KINC,1,IXV,MFLU,MZFLU,NPT1,NC1,IK1,IRET)
      IF(IRET.EQ.0)RETURN
 
C Fin lecture Arguments **************************************************
 
 
      IF(KIMPL.EQ.0)THEN
      TYPE=' '
      CALL ACMO(KIZG,NOMI,TYPE,IZG1)
      IF(TYPE.NE.'CHPOINT ')THEN
      NC=IZTU1.VPOCHA(/2)
      TYPE='SOMMET'
      CALL CRCHPT(TYPE,SPGID,NC,2,IZG1)
C     SEGACT IZG1
C     MSOUPO=IZG1.IPCHP(1)
C     SEGACT MSOUPO
C     NOCOMP(1)=NOMI
C     SEGDES MSOUPO
      CALL ECMO(KIZG,NOMI,'CHPOINT ',IZG1)
      ENDIF
 
      ELSE
      NC=IZTU1.VPOCHA(/2)
      TYPE='SOMMET'
C? pour plutard    CALL CRCHPT(TYPE,SPGZ,NC,2,IZG1)
      CALL CRCHPT(TYPE,SPGID,NC,2,IZG1)
      SEGACT IZG1
      MSOUPO=IZG1.IPCHP(1)
      SEGACT MSOUPO*MOD
      NOCOMP(1)=NOMI
      ENDIF
 
      CALL LICHT(IZG1,IZGG1,TYPC,IGEOM)
 
      IF(IGEOM.NE.SPGID)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
 
      SEGACT MELEME
      NBSOUS=LISOUS(/1)
      IF(NBSOUS.EQ.0)NBSOUS=1
      NUTOEL=0
 
      NPTD=IZTU1.VPOCHA(/1)
      IES=IDIM
 
      DO 1 L=1,NBSOUS
      IPT1=MELEME
      IF(NBSOUS.NE.1)IPT1=LISOUS(L)
      SEGACT IPT1
 
      MCHAML=ICHAML(L)
      SEGACT MCHAML
      MELVAL=IELVAL(1)
      SEGACT MELVAL
 
      NP  =IPT1.NUM(/1)
      NBEL=IPT1.NUM(/2)
 
      IF(KPOIND.EQ.0)THEN
 
      CALL ZXFIMP(NBEL,NUTOEL,NP,LECT,IPT1.NUM,
     & VELCHE,IZGG1.VPOCHA,MZFLU.VPOCHA,IK1)
 
      ELSE
      IPT2=MELEMP
      IF(NBSOUS.NE.1)IPT2=LISOUS(L)
      SEGACT IPT2
 
      IF(MQUAD.NE.0)THEN
        IF(KPOIND.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//'PRP0'
        IF(KPOIND.EQ.3)NOM0=NOMS(IPT1.ITYPEL)//'PRP0'
        IF(KPOIND.EQ.4)NOM0=NOMS(IPT1.ITYPEL)//'PRP1'
      ELSEIF(MACRO.NE.0)THEN
        IF(KPOIND.EQ.2)NOM0=NOMS(IPT1.ITYPEL)
        IF(KPOIND.EQ.3)NOM0=NOMS(IPT1.ITYPEL)//'MCP0'
        IF(KPOIND.EQ.4)NOM0=NOMS(IPT1.ITYPEL)//'MCP1'
      ELSE
        IF(KPOIND.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)
 
      NK=NUTOEL
      DO 21 K=1,NBEL
 
      NK=NK+1
      JC=(1-IK1)*(NK-1)+1
      DO 109 I=1,NP
      J=IPT1.NUM(I,K)
      DO 1091 N=1,IDIM
      XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N)
 1091 CONTINUE
 109  CONTINUE
 
      CALL CALJBC(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,AIRE)
 
      DO 39 M=1,MP1
      M1=LECT(IPT2.NUM(M,K))
      U=0.D0
      DO 33 LL=1,NPG
      U=U+IZF1.FN(M,LL)*MZFLU.VPOCHA(JC,1)*PGSQ(LL)
 33   CONTINUE
      IZGG1.VPOCHA(M1,1)=IZGG1.VPOCHA(M1,1)-U
 39   CONTINUE
 
 21   CONTINUE
 
      ENDIF
 
      SEGDES IPT1,MCHAML,MELVAL
      NUTOEL=NUTOEL+NBEL
 
 1    CONTINUE
C     SEGDES MZFLU,IZVOL
      SEGDES MZFLU
 
      IF(KIMPL.NE.0)THEN
      TYPE=' '
      CALL ACMO(MTAB1,'SMBR',TYPE,MCHPO2)
      IF(TYPE.NE.'CHPOINT')THEN
      CALL ECMO(MTAB1,'SMBR','CHPOINT',IZG1)
      ELSE
      CALL ECROBJ('CHPOINT',MCHPO2)
      CALL ECROBJ('CHPOINT',IZG1)
      CALL PRFUSE
      CALL LIROBJ('CHPOINT',MCHPOI,1,IRET)
      CALL ECMO(MTAB1,'SMBR','CHPOINT',MCHPOI)
      ENDIF
 
      ENDIF
 
      SEGDES MELEME
 
      SEGDES IZTU1
      SEGDES IZG1,IZGG1
      SEGDES LINCO
      SEGSUP MLENTI
 
      IPAS=1
      RETURN
 90   CONTINUE
C        Option %m1:8 incompatible avec les données
            MOTERR( 1: 8) = MTERR
            CALL ERREUR(803)
            RETURN
 
 1002 FORMAT(10(1X,1PE11.4))
      END
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales