Télécharger vnimp.eso

Retour à la liste

vnimp
C VNIMP     SOURCE    FANDEUR   22/01/03    21:15:56     11136          
      SUBROUTINE VNIMP
C************************************************************************
C
C************************************************************************
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
 
-INC PPARAM
-INC CCOPTIO
-INC CCGEOME
-INC SIZFFB
      POINTEUR IZF1.IZFFM,IZH2.IZHR
 
-INC SMCOORD
-INC SMLENTI
-INC SMELEME
      POINTEUR MELEM1.MELEME,MELEMS.MELEME,MELEML.MELEME
      POINTEUR MELENE.MELEME,MELEMQ.MELEME
-INC SMCHPOI
      POINTEUR IZTU1.MPOVAL,IZTU2.MPOVAL,VNORM.MPOVAL
 
-INC SMMATRIK
-INC CCREEL
-INC SMLMOTS
      POINTEUR LINCO.MLMOTS
      CHARACTER*8 TYPE,TYPC
      CHARACTER*(LOCOMP) NOMP,NOMI,NOM,NOM0,NOMZ
      PARAMETER (NTB=1)
      CHARACTER*8 LTAB(NTB)
      DIMENSION KTAB(NTB),IXV(3),RO(1)
      DATA LTAB/'KIZX    '/,RO/1.D0/
C*****************************************************************************
CVNIMP
C       write(6,*)' DEBUT VNIMP '
 
      IAXI=0
      IF(IFOMOD.EQ.0)IAXI=2
 
      DEUPI=1.D0
      IF(IAXI.NE.0)DEUPI=2.D0*XPI
 
      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
 
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
 
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 (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)
      TYPE=' '
      CALL ACMO(MTABX,'DOMZ',TYPE,MMODEL)
      IF(TYPE.NE.'MMODEL')THEN
C On attend un des objets : %m1:8 %m9:16 %m17:24 %m25:32 %m33:40
         MOTERR( 1: 8) = ' MMODEL '
         MOTERR( 8:16) = ' MMODEL '
         MOTERR(17:24) = ' MMODEL '
         MOTERR(25:32) = ' MMODEL '
         MOTERR(33:40) = ' MMODEL '
         CALL ERREUR(471)
         RETURN
      ENDIF
 
C E/  MMODEL : Pointeur de la table contenant l'information cherchée
C  /S IPOINT : Pointeur sur la table DOMAINE
C  /S INEFMD : Type formulation INEFMD=1 LINE,=2 MACRO,=3 QUADRATIQUE
C                               INEFMD=4 LINB
 
      CALL LEKMOD(MMODEL,MTABZ,INEFMD)
 
      IF(INEFMD.EQ.4.AND.KPRE.NE.5)THEN
C% Données incompatibles
      CALL ERREUR(21)
      RETURN
      ENDIF
 
      CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
      CALL LEKTAB(MTABZ,'SOMMET',MELEMS)
      CALL LEKTAB(MTABZ,'MLGVNIMP',MELEML)
      IF(KPRE.EQ.5)CALL LEKTAB(MTABZ,'MSOMMET',MELEML)
      IF(KPRE.EQ.5)CALL LEKTAB(MTABZ,'MMAIL'  ,MELENE)
 
 
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.AND.NBINC.NE.2).AND.KPRE.NE.5)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
      ELSEIF(NBINC.NE.2.AND.KPRE.EQ.5)THEN
C        Indice %m1:8 : contient plus de %i1 %m9:16
         MOTERR( 1:8) = 'LISTINCO'
         INTERR(1) = 2
         MOTERR(9:16) = ' MOTS   '
         CALL ERREUR(799)
         RETURN
      ENDIF
 
      NOMI=LINCO.MOTS(1)
      NOMP='LVNP'
      IF(KPRE.EQ.5)NOMP=LINCO.MOTS(2)
C
C- Récupération de l'inconnue
C
      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)
         NINKO = IZTU1.VPOCHA(/2)
         IF (NINKO.NE.IDIM) THEN
C        Indice %m1:8 : Le %m9:16 n'a pas le bon nombre de composantes
            MOTERR        = 'INC '//NOMI
            MOTERR( 9:16) = 'CHPOINT '
            CALL ERREUR(784)
            RETURN
         ENDIF
C On fait pointer ces deux tableaux sur le champ U inconu (tjs présent) pour
C eviter de les enlever lors de l'appel FORTRAN si les options sont absentes
      ENDIF
 
      IF(KPRE.EQ.5)THEN
 
      TYPE=' '
      CALL ACMO(KINC,NOMP,TYPE,MCHPO2)
      IF(TYPE.NE.'CHPOINT ')THEN
C        Indice %m1:8 : ne contient pas un objet de type %m9:16
         MOTERR( 1: 8) = 'INC '//NOMP
         MOTERR( 9:16) = 'CHPOINT '
         CALL ERREUR(800)
         RETURN
      ELSE
      CALL LICHT(MCHPO2,IZTU2,TYPC,MELEM2)
         NINK2 = IZTU2.VPOCHA(/2)
         IF (NINK2.NE.1) THEN
C        Indice %m1:8 : Le %m9:16 n'a pas le bon nombre de composantes
            MOTERR( 1: 8) = 'INC '//NOMP
            MOTERR( 9:16) = 'CHPOINT '
            CALL ERREUR(784)
            RETURN
         ENDIF
C On fait pointer ces deux tableaux sur le champ U inconu (tjs présent) pour
C eviter de les enlever lors de l'appel FORTRAN si les options sont absentes
      ENDIF
 
      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,IPADI)
 
      CALL VERPAD(IPADI,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
 
C*****************************************************************************
 
 
 
C*************************************************************************
C Lecture des coefficients
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.LT.2)THEN
C           Indice %m1:8 : nombre d'arguments incorrect
            MOTERR(1:8) = 'IARG    '
            CALL ERREUR(804)
      RETURN
      ENDIF
 
      TYPE=' '
      CALL ACMO(MTABX,'ARG1',TYPE,MMODER)
      IF(TYPE.NE.'MMODEL')THEN
C On attend un des objets : %m1:8 %m9:16 %m17:24 %m25:32 %m33:40
         MOTERR( 1: 8) = ' MMODEL '
         MOTERR( 8:16) = ' MMODEL '
         MOTERR(17:24) = ' MMODEL '
         MOTERR(25:32) = ' MMODEL '
         MOTERR(33:40) = ' MMODEL '
         CALL ERREUR(471)
         RETURN
      ENDIF
 
C E/  MMODEL : Pointeur de la table contenant l'information cherchée
C  /S IPOINT : Pointeur sur la table DOMAINE
C  /S INEFMD : Type formulation INEFMD=1 LINE,=2 MACRO,=3 QUADRATIQUE
C                               INEFMD=4 LINB
 
      CALL LEKMOD(MMODER,MTABR,INEFDR)
      IF(INEFDR.NE.INEFMD)THEN
C% Données incompatibles
      CALL ERREUR(21)
      RETURN
      ENDIF
      CALL LEKTAB(MTABR,'SOMMET',MELEMR)
      CALL KRIPAD(MELEMR,MLENTI)
      CALL VERPAD(MLENTI,MELEMS,IRET)
      IF(IRET.NE.0)THEN
C     Le support du %m1:8 est incompatible avec celui-ci
         MOTERR(1: 8) = 'MODELE'
         CALL ERREUR(708)
      RETURN
      ENDIF
      SEGSUP MLENTI
 
 
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      IXV(1)=MELEMS
C     IXV(1)=0
      IXV(2)=1
      IXV(3)=0
      CALL LEKCOF('Opérateur VNIMP :',
     & MTABX,KINC,2,IXV,IZTG1,VNORM,NPT1,NC1,IK1,IRET)
      IF(IRET.EQ.0)RETURN
 
C     write(6,*)' Operateur VNIMP : Fin lecture Arguments '
C Fin lecture Arguments ************************************************
 
      IF(KPRE.NE.5)THEN
C
C Pressions discontinues
C
      CALL LEKTAB(MTABZ,'NORMALEV',MNORM)
      CALL LICHT(MNORM,MPOVA2,TYPC,MELEML)
 
      NAT=2
      NSOUPO=1
      SEGACT MELEML
      N=MELEML.NUM(/2)
      NC=1
      SEGINI MCHPO1,MSOUP1,MPOVA1
      MCHPO1.IFOPOI=IFOUR
      MCHPO1.MOCHDE='                                                  '
      MCHPO1.MTYPOI='SMBR'
      MCHPO1.JATTRI(1)=2
      MCHPO1.IPCHP(1)=MSOUP1
      MSOUP1.NOCOMP(1)='LVNP'
      MSOUP1.IGEOC=MELEML
      MSOUP1.IPOVAL=MPOVA1
 
      NRIGE=7
      NKID =9
      NKMT =7
      NMATRI=1
      SEGINI MATRIK
 
      NBSOUS=1
 
      IRIGEL(1,1)=MELEMS
      IRIGEL(2,1)=MELEML
      IRIGEL(7,1)=4
      IF(KFORM.EQ.2)IRIGEL(7,1)=9
      NBME=NINKO
      SEGINI IMATRI
      IRIGEL(4,1)=IMATRI
      KSPGP=MELEMS
      KSPGD=MELEML
      DO 102 I=1,NBME
      WRITE(NOM,FMT='(I1)')I
      NOM=NOM(1:1)//NOMI(1:LOCOMP-1)
      LISPRI(I)= NOM
      LISDUA(I)='LVNP'
 102  CONTINUE
 
 
      NUTOEL=0
      NBEL =MELEML.NUM(/2)
      NP =1
      MP =1
 
      SEGINI IPM1,IPM2
      LIZAFM(1,1)=IPM1
      LIZAFM(1,2)=IPM2
      IF(NBME.EQ.3)THEN
      SEGINI IPM3
      LIZAFM(1,3)=IPM3
      ENDIF
 
      DO 301 K=1,NBEL
      IPM1.AM(K,1,1)= MPOVA2.VPOCHA(K,1)
      IPM2.AM(K,1,1)= MPOVA2.VPOCHA(K,2)
      IF(IDIM.EQ.3)IPM3.AM(K,1,1)= MPOVA2.VPOCHA(K,3)
 301  CONTINUE
 
      ELSEIF(KPRE.EQ.5)THEN
C
C Cas Pressions continue
C
C     write(6,*)'Pressions continue'
 
      IF(INEFMD.EQ.2)CALL LEKTAB(MTABZ,'MACRO',MELEME)
      CALL LEKTAB(MTABR,'OENVELOP',MCHPOR)
      CALL LEKTAB(MTABZ,'QUAF    ',MELEMQ)
 
      CALL LICHT(MCHPOR,MPOVA3,TYPC,IGEOM)
      CALL KRIPAD(IGEOM,MLENT3)
 
      NAT=2
      NSOUPO=1
      SEGACT MELEML
      N=MELEML.NUM(/2)
      NC=1
      SEGINI MCHPO1,MSOUP1,MPOVA1
      MCHPO1.IFOPOI=IFOUR
      MCHPO1.MOCHDE='                                                  '
      MCHPO1.MTYPOI='SMBR'
      MCHPO1.JATTRI(1)=2
      MCHPO1.IPCHP(1)=MSOUP1
      MSOUP1.NOCOMP(1)=NOMP
      MSOUP1.IGEOC=MELEML
      MSOUP1.IPOVAL=MPOVA1
 
      NRIGE=7
      NKID =9
      NKMT =7
      NMATRI=1
      SEGINI MATRIK
 
      SEGACT MELEME,MELENE
      NBSOUS=LISOUS(/1)
 
      IRIGEL(1,1)=MELEME
      IRIGEL(2,1)=MELENE
      IRIGEL(7,1)=-3
      NBME=NINKO
      NBSOUS=MAX(1,LISOUS(/1))
      SEGINI IMATRI
      IRIGEL(4,1)=IMATRI
      KSPGP=MELEMS
      KSPGD=MELEML
      DO 202 I=1,NBME
      WRITE(NOM,FMT='(I1)')I
      NOM=NOM(1:1)//NOMI(1:LOCOMP-1)
      LISPRI(I)=NOM
      LISDUA(I)=NOMP
 202  CONTINUE
 
      SEGACT MELEMQ
      NUTOEL=0
      DO 203 L=1,MAX(1,LISOUS(/1))
      IPT1=MELEME
      IPT2=MELENE
      IPT3=MELEMQ
      IF(LISOUS(/1).NE.0)IPT1=LISOUS(L)
      IF(MELENE.LISOUS(/1).NE.0)IPT2=MELENE.LISOUS(L)
      IF(MELEMQ.LISOUS(/1).NE.0)IPT3=MELEMQ.LISOUS(L)
      SEGACT IPT1,IPT2,IPT3
 
      NBEL =IPT1.NUM(/2)
      NP =IPT1.NUM(/1)
      MP =IPT2.NUM(/1)
      NQ =IPT3.NUM(/1)
      IF(NQ.EQ.3)THEN
      IFA=2
      ELSE
      IFA=NQ
      ENDIF
 
 
      SEGINI IPM1,IPM2
      LIZAFM(L,1)=IPM1
      LIZAFM(L,2)=IPM2
      IF(NBME.EQ.3)THEN
      SEGINI IPM3
      LIZAFM(L,3)=IPM3
      ENDIF
 
      NOM0=NOMS(IPT1.ITYPEL)
      IF(INEFMD.EQ.1)NOM0=NOMS(IPT1.ITYPEL)//'P1P1'
      IF(INEFMD.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//'MCF1'
      IF(INEFMD.EQ.3)NOM0=NOMS(IPT1.ITYPEL)//'PFP1'
      IF(INEFMD.EQ.4)NOM0=NOMS(IPT1.ITYPEL)//'P1P1'
 
      CALL KALPBG(NOM0,'FONFORM ',IZFFM)
      SEGACT IZFFM*MOD
      IZHR=KZHR(1)
      IZH2=KZHR(2)
      IZF1=KTP(1)
      SEGACT IZHR*MOD,IZF1*MOD
      NES=GR(/1)
      NPG=GR(/3)
 
      DO 201 K=1,NBEL
 
      DO 20 I=1,NP
            J1 = IPT1.NUM(I,K)
            DO 10 N=1,IDIM
               XYZ(N,I) = XCOOR((J1-1)*(IDIM+1)+N)
  10        CONTINUE
  20     CONTINUE
        CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,
     &        NES,IDIM,NP,NPG,IAXI,AIRE,AJ,ASGN)
 
      BSGN=0.D0
      DO 21 N=1,IDIM
      BSGN=BSGN+
     &(AJ(N,IDIM,1)*MPOVA3.VPOCHA(MLENT3.LECT(IPT3.NUM(IFA,K)),N))
 21   CONTINUE
      CSGN=-1.D0
      IF(BSGN.LT.0.D0)CSGN=1.D0
 
      DO  I=1,NP
      DO  J=1,MP
      UX=0.D0
      UY=0.D0
      UZ=0.D0
      UT=0.D0
      DO 304 LG=1,NPG
      UX=UX+FN(I,LG)*IZF1.FN(J,LG)*AJ(1,IDIM,LG)*PGSQ(Lg)*DEUPI*RPG(Lg)
      UY=UY+FN(I,LG)*IZF1.FN(J,LG)*AJ(2,IDIM,LG)*PGSQ(Lg)*DEUPI*RPG(Lg)
      IF(IDIM.EQ.3)
     &UZ=UZ+FN(I,LG)*IZF1.FN(J,LG)*AJ(3,IDIM,LG)*PGSQ(Lg)*DEUPI*RPG(Lg)
 304  CONTINUE
      IPM1.AM(K,I,J)= UX*CSGN
      IPM2.AM(K,I,J)= UY*CSGN
      IF(IDIM.EQ.3)IPM3.AM(K,I,J)= UZ*CSGN
      ENDDO
      ENDDO
 201  CONTINUE
      SEGSUP IZFFM,IZHR,IZH2,IZF1
      SEGDES IPT1,IPT2
 203  CONTINUE
 
      SEGSUP MLENT3
 
      ENDIF
 
      CALL ECROBJ('MATRIK',MATRIK)
      CALL ECROBJ('CHPOINT',MCHPO1)
 
      RETURN
 1001 FORMAT(20(1X,I5))
 1002 FORMAT(10(1X,1PE11.4))
      END
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

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