Sources Esope | Cast3M

vtimp
C VTIMP     SOURCE    FANDEUR   22/01/03    21:15:56     11136          
      SUBROUTINE VTIMP
C************************************************************************
C   CALCUL DE LA MATRICE M ( P*DIV(U) )  --> AM(NP,IES,NEL)
C
C   IKAS=1  V NORMALE IMPOSEE
C   IKAS=-1 V TANGENTE IMPOSEE
C
C************************************************************************
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
 
-INC PPARAM
-INC CCOPTIO
-INC CCGEOME
-INC SIZFFB
      POINTEUR IZF1.IZFFM
 
-INC SMCOORD
-INC SMLENTI
      POINTEUR IPADI.MLENTI
-INC SMELEME
      POINTEUR MELEM1.MELEME,MELEMC.MELEME,MELEMS.MELEME,MELEML.MELEME
-INC SMCHPOI
      POINTEUR IZTU1.MPOVAL,VTANG.MPOVAL
 
-INC SMMATRIK
 
-INC SMLMOTS
      POINTEUR LINCO.MLMOTS
      CHARACTER*8 NOMZ,NOMI,TYPE,TYPC,NOM0,NOMA,NOM
      PARAMETER (NTB=1)
      CHARACTER*8 LTAB(NTB)
      DIMENSION KTAB(NTB),IXV(3),RO(1)
      DATA LTAB/'KIZX    '/,RO/1.D0/
C*****************************************************************************
CVTIMP
C     write(6,*)' DEBUT VTIMP '
 
      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
 
      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
      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
 
      CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
      CALL LEKTAB(MTABZ,'SOMMET',MELEMS)
      CALL LEKTAB(MTABZ,'CENTRE',MELEMC)
      CALL LEKTAB(MTABZ,'MLGVTIMP',MELEML)
 
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)
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( 1: 8) = '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( 1: 8) = '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
 
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(1: 8) = '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.NE.1)THEN
C           Indice %m1:8 : nombre d'arguments incorrect
            MOTERR(1:8) = 'IARG    '
            CALL ERREUR(804)
      RETURN
      ENDIF
 
      IXV(1)=MELEMS
      IXV(2)=1
      IXV(3)=0
      CALL LEKCOF('Opérateur VTIMP :',
     & MTABX,KINC,1,IXV,IZTG1,VTANG,NPT1,NC1,IK1,IRET)
      IF(IRET.EQ.0)RETURN
 
 
C     write(6,*)' Operateur VTIMP : Fin lecture Arguments '
C Fin lecture Arguments ************************************************
 
      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)='LVTP'
      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
      NBME=NINKO
      SEGINI IMATRI
      IRIGEL(4,1)=IMATRI
      KSPGP=MELEMS
      KSPGD=MELEML
      DO 102 I=1,NBME
      WRITE(NOM,FMT='(I1,A7)')I,NOMI(1:7)
      LISPRI(I)=NOM(1:4)//'    '
      LISDUA(I)='LVTP'//'    '
 102  CONTINUE
 
      SEGACT MELEMS
      NP =1
      MP =1
      NBEL=MELEMS.NUM(/2)
 
      SEGINI IPM1,IPM2
      LIZAFM(1,1)=IPM1
      LIZAFM(1,2)=IPM2
      IF(NBME.EQ.3)THEN
      SEGINI IPM3
      LIZAFM(1,3)=IPM3
      ENDIF
 
      CALL KRIPAD(MELEMS,MLENTI)
 
      SEGACT MELEME
 
      NBSOU1=LISOUS(/1)
      IF(NBSOU1.EQ.0)NBSOU1=1
 
      NUTOEL=0
      DO 101 L=1,NBSOU1
      IPT1=MELEME
      IF(NBSOU1.NE.1)IPT1=LISOUS(L)
      SEGACT IPT1
      NOM0=NOMS(IPT1.ITYPEL)//'    '
      CALL KALPBG(NOM0,'FONFORM ',IZFFM)
      SEGACT IZFFM*MOD
      IZHR=KZHR(1)
      SEGACT IZHR*MOD
      NES=GR(/1)
      NPG=GR(/3)
      NP=IPT1.NUM(/1)
      NBELEM=IPT1.NUM(/2)
 
      DO 301 K=1,NBELEM
      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,SGN)
 
 
      U0=1.D0
      DO 302 I=1,NP
      K1=LECT(IPT1.NUM(I,K))
      UX=0.D0
      UY=0.D0
      UZ=0.D0
      UT=0.D0
      DO 304 LG=1,NPG
      UX=UX+FN(I,LG)*AJ(1,1,LG)
      UY=UY+FN(I,LG)*AJ(2,1,LG)
      IF(IDIM.EQ.3)UZ=UZ+FN(I,LG)*AJ(3,1,LG)
      UNL=0.D0
      DO 321 J=1,NP
      K2=LECT(IPT1.NUM(J,K))
      NK=(K2-1)*(1-IK1)+1
      UNL=UNL+VTANG.VPOCHA(NK,1)*FN(J,L)
 321  CONTINUE
      UT=UT+FN(I,LG)*UNL
 304  CONTINUE
      IPM1.AM(K1,1,1)=IPM1.AM(K1,1,1)-UX
      IPM2.AM(K1,1,1)=IPM2.AM(K1,1,1)-UY
      IF(IDIM.EQ.3)IPM3.AM(K1,1,1)=IPM3.AM(K1,1,1)-UZ
      MPOVA1.VPOCHA(K1,1)=MPOVA1.VPOCHA(K1,1)-UT
 302  CONTINUE
 301  CONTINUE
 
 101  CONTINUE
 
 
 
      CALL ECROBJ('MATRIK',MATRIK)
      CALL ECROBJ('CHPOINT',MCHPO1)
 
      RETURN
 1001 FORMAT(20(1X,I5))
 1002 FORMAT(10(1X,1PE11.4))
      END