Télécharger kmic.eso

Retour à la liste

kmic
C KMIC      SOURCE    CB215821  20/11/25    13:31:36     10792          
      SUBROUTINE KMIC(IKAS,MTABX)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C************************************************************************
C    Operateur KMAC
C
C    OBJET   : Cree un objet de type MATRIK
C
C
C  IKAS=1 KMAB calcul de B     (DIV U)
C  IKAS=2 KMBT calcul de Bt uniquement (GRAD P)
C  IKAS=3 KBBT calcul de B assemblage pour B et Bt
C
C***********************************************************************
 
-INC PPARAM
-INC CCOPTIO
-INC CCGEOME
-INC CCREEL
-INC SMCOORD
-INC SIZFFB
-INC SMCHPOI
      POINTEUR IZCH2.MCHPOI,IZCCH2.MPOVAL
      POINTEUR IZDV.MCHPOI,IZDDV.MPOVAL,IZTU1.MPOVAL,TETAN.MPOVAL
      POINTEUR IZTG1.MCHPOI,IZTGG1.MPOVAL,IZBETA.MPOVAL,VISCO.MPOVAL
-INC SMMATRIK
      POINTEUR IPM.IZAFM
      SEGMENT IMATRS
      ENDSEGMENT
      POINTEUR IPMS.IZAFM,IPS1.IZAFM,IPS2.IZAFM,IPS3.IZAFM
 
      PARAMETER (LRV=64)
      DIMENSION VISC(3)
 
-INC SMLENTI
      POINTEUR IZIPAD.MLENTI,MLENTI1.MLENTI,MLENTI2.MLENTI
-INC SMLMOTS
      POINTEUR LINCO.MLMOTS
-INC SMELEME
      POINTEUR MELEMZ.MELEME,MELEMB.MELEME,MELSTB.MELEME
      POINTEUR MELEM1.MELEME,MELES1.MELEME,MCTREI.MELEME
      POINTEUR IGEOM.MELEME,MELEMA.MELEME
      POINTEUR IZLEMC.MELEME,MELEMS.MELEME,MELEMC.MELEME
      POINTEUR MELEMI.MELEME,MELEMP.MELEME
      POINTEUR MACRO1.MELEME
 
      CHARACTER*8 TYPE,TYPC,NOMZ,NOMP,NOMD,NOM0,NOMK
      CHARACTER*8 NOMPP,NOMDD
      CHARACTER*4 NOM
      INTEGER IPAD,IPAD2,IK
      REAL*8 KAUX,TETA1
      DIMENSION IXV(5)
C
C*************************************************************************
CKMIC
C     write(6,*)' Operateur KMIC MTABX=',MTABX
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(MTABX,'IARG',IARG)
      IF(IARG.GE.3)THEN
C 3ème COEF
      CALL ACMM(MTABX,'ARG3',NOMK)
 
      IF(NOMK.EQ.'SOMM')THEN
*     Cet opérateur calcul un gradient ou div de sommets vers sommets
*     sans contrainte de type multiplicateur de Lagrange.
      CALL DIVS(IKAS,MTABX)
      RETURN
*
      ELSEIF(NOMK.EQ.'NOMU')THEN
*     Cet opérateur calcule la divergence (ou grad) de sommets vers ctp1
*     sans contrainte de cmultiplicateurs de Lagrange.
      CALL DIVC(IKAS,MTABX)
C     write(6,*)' RETOUR KKMIC'
      RETURN
 
      ENDIF
      ENDIF
 
      CALL ACME(MTAB1,'DISCPRES',IDISCP)
 
      CALL ACME(MTAB1,'NAVISTOK',NASTOK)
      IF(NASTOK.EQ.0)THEN
      CALL ZKMIC(IKAS,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 -> CN
C KFORM = 0 -> EFSI  1 -> EF    2 -> VF  3 -> EFMC
C KPRE=3 pression P0   KPRE=4 pression P1  KPRE=2 cas macro 1ère génération
C KPRE=5 pression MSOMMET
 
      NOMUL=1
      IAXI=0
      IK=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)
C     write(6,*)'KBBT KPRE=',KPRE
      CALL ACMF(KOPTI,'AIMPL',TETA1)
      CALL ACME(KOPTI,'KFORM',KFORM)
C     CALL ACME(KOPTI,'IKOMP',IKOMP)
 
      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
 
      CALL LEKMOD(MMODEL,MTABZ,INEFMD)
 
      CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
      CALL LEKTAB(MTABZ,'SOMMET',MELEMS)
      CALL LEKTAB(MTABZ,'CENTRE',MELEMC)
C??   CALL LEKTAB(MTABZ,'MACRO',MACRO)
      MACRO1=0
      IF(INEFMD.EQ.2)CALL LEKTAB(MTABZ,'MACRO1',MACRO1)
C??   IF(INEFMD.EQ.3)CALL LEKTAB(MTABZ,'QUADRATI',MQUAD)
      IF (IERR.NE.0) RETURN
      IF(INEFMD.EQ.4.AND.KPRE.NE.5)THEN
C% Données incompatibles
      CALL ERREUR(21)
      RETURN
      ENDIF
 
      MELEMI=MELEME
      IF(MACRO1.NE.0.AND.KPRE.NE.2)THEN
      MELEMI=MACRO1
      ENDIF
 
      IF(KPRE.EQ.2.AND.INEFMD.EQ.3)KPRE=3
      IF(INEFMD.EQ.1.AND.KPRE.NE.5)KPRE=2
 
      IF(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)
      ELSEIF(KPRE.EQ.2)THEN
      CALL LEKTAB(MTABZ,'CENTRE',MELEMC)
      MELEMP=MELEMC
      ELSEIF(KPRE.EQ.5)THEN
      CALL LEKTAB(MTABZ,'MSOMMET',MELEMC)
      CALL LEKTAB(MTABZ,'MMAIL  ',MELEMP)
      ENDIF
 
      IF(IDISCP.EQ.0)THEN
      CALL ECME(MTAB1,'DISCPRES',KPRE)
      ELSEIF(IDISCP.NE.KPRE)THEN
C% Données incompatibles
      CALL ERREUR(21)
      RETURN
      ENDIF
 
      IKOMP=1
C     IF(KPRE.EQ.5)IKOMP=0
C     IF(KPRE.EQ.5)IKOMP=1
 
C*************************************************************************
C VERIFICATIONS SUR LES INCONNUES
 
C     write(6,*)' Verification sur les inconnues '
      TYPE='LISTMOTS'
      CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO)
      IF (IERR.NE.0) RETURN
      SEGACT LINCO
      NBINC=LINCO.MOTS(/2)
      IF(NBINC.NE.2)THEN
      WRITE(6,*)'Operateur KMAB '
      WRITE(6,*)'Nombre d''inconnues incorrecte : ',NBINC,' On attend 2'
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
 
C     On recupere PHI n et TETA n pour Cranck-Nicholson
      NOMP=LINCO.MOTS(1)
      TYPE=' '
      CALL ACMO(KINC,NOMP,TYPE,MCHPOI)
      IF(TYPE.NE.'CHPOINT ')THEN
         WRITE(6,*)' Opérateur KMAB :'
         WRITE(6,*)' L objet CHPOINT ',NOMP,
     &        ' n existe pas dans la table'
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(MCHPOI,IZTU1,TYPC,IGEOM0)
      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
C Inconnue Primale
 
C     write(6,*)' Verification inconnue primale '
      CALL KRIPAD(IGEOM0,MLENTI)
      IF(IKAS.EQ.1.OR.IKAS.EQ.3)THEN
      MELEMK=MELEMS
      ELSE
      MELEMK=MELEMC
      ENDIF
 
C     write(6,*)' MELEMK=',melemk
      CALL VERPAD(MLENTI,MELEMK,IRET)
      IF(IRET.NE.0)THEN
      WRITE(6,*)' Opérateur KMAB '
      WRITE(6,*)' La zone ',NOMZ,' n''est pas incluse dans le domaine'
     &        , ' de définition de l''inconnue ',NOMP
      WRITE(6,*)' MELEMK=',melemk,' IGEOM0=',IGEOM0
C     Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
         MOTERR(1: 8) = 'INC '//NOMP
         MOTERR(9:16) = 'CHPOINT '
         CALL ERREUR(788)
      IPAS=0
      RETURN
      ENDIF
 
      SEGSUP MLENTI
 
C*************************************************************************
 
      NOMD=LINCO.MOTS(2)
      TYPE=' '
      CALL ACMO(KINC,NOMD,TYPE,MCHPOI)
      IF(TYPE.NE.'CHPOINT ')THEN
         WRITE(6,*)' Opérateur KMAB :'
         WRITE(6,*)' L objet CHPOINT ',NOMD,
     &        ' n existe pas dans la table'
C        Indice %m1:8 : ne contient pas un objet de type %m9:16
         MOTERR( 1: 8) = 'INC '//NOMD
         MOTERR( 9:16) = 'CHPOINT '
         CALL ERREUR(800)
         RETURN
      ELSE
         CALL LICHT(MCHPOI,TETAN,TYPC,IGEOM0)
      ENDIF
 
      NC=TETAN.VPOCHA(/2)
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
C Inconnue Duale
 
C     write(6,*)' IGEOM0=',igeom0
      CALL KRIPAD(IGEOM0,MLENTI)
      IF(IKAS.EQ.1.OR.IKAS.EQ.3)THEN
      MELEMK=MELEMC
      ELSE
      MELEMK=MELEMS
      ENDIF
 
C     write(6,*)' Verification inconnue duale ',MELEMK
      CALL VERPAD(MLENTI,MELEMK,IRET)
      IF(IRET.NE.0)THEN
      WRITE(6,*)' Opérateur KMAB '
      WRITE(6,*)' La zone ',NOMZ,' n''est pas incluse dans le domaine'
     &         ,' de définition de l''inconnue ',NOMD
      WRITE(6,*)' MELEMK=',melemk,' IGEOM0=',IGEOM0
C     Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
         MOTERR(1: 8) = 'INC '//NOMD
         MOTERR(9:16) = 'CHPOINT '
         CALL ERREUR(788)
      IPAS=0
      RETURN
      ENDIF
 
      SEGSUP MLENTI
 
C*************************************************************************
C Lecture du ou des coefficients
C Type du coefficient :
C IK1=0 CHPOINT IK1=1 scalaire  IK1=2 vecteur
 
C     write(6,*)' Verification sur les coefficients '
      CALL ACME(MTABX,'IARG',IARG)
C     write(6,*)' IARG=',iarg
 
      IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.2)THEN
      TYPE=' '
      CALL ACMO(MTABZ,'MELSTB',TYPE,MELSTB)
      SEGACT MELSTB
      IF(IDIM.EQ.2)NBELEM=MELSTB.NUM(/2)/4
      IF(IDIM.EQ.3)NBELEM=MELSTB.NUM(/2)/8
      NBNN=MELSTB.NUM(/1)
      NBSOUS=0
      NBREF=0
      SEGINI MELEMA
      MELEMA.ITYPEL=MELSTB.ITYPEL
 
      NKPE=4
      IF(IDIM.EQ.3)NKPE=8
      do 4878 k=1,nbelem
      mi=(k-1)*NKPE+1
      do 4879 i=1,nbnn
      MELEMA.num(i,k)=melstb.num(i,mi)
 4879 continue
C     write(6,*)k,(MELEMA.num(i,k),i=1,nbnn)
 4878 continue
 
      TYPE=' '
      CALL ACMO(MTABZ,'MCHPOC',TYPE,MCHPOC)
      TYPE=' '
      CALL ACMO(MTABZ,'CENTRE',TYPE,MCTREI)
      ENDIF
 
 
C 1er COEF
 
      IXV(1)=MELEMC
      IXV(2)=1
      IXV(3)=0
      IXV(4)=MELEMS
      IXV(5)=-MELEMS
      IRET=5
      CALL LEKCOF('Opérateur KMAB :',
     & MTABX,KINC,1,IXV,IZTG1,IZTGG1,NPT1,NC1,IK1,IRET)
      IF(IRET.EQ.0)RETURN
      IF(IK1.GE.4)CALL KRIPAD(MELEMS,MLENTI)
 
 
      BETA0=1.D-1
      IF(IARG.EQ.2)THEN
C 2ème COEF
      IXV(1)=0
      IXV(2)=1
      IXV(3)=0
      CALL LEKCOF('Opérateur KMAB :',
     & MTABX,KINC,2,IXV,IZTG2,IZBETA,NPT2,NC2,IK2,IRET)
      IF(IRET.EQ.0)RETURN
      BETA0=IZBETA.VPOCHA(1,1)
C?    IF(IZBETA.VPOCHA(1,1).LT.0.D0)THEN
C% Données incompatibles
C?    CALL ERREUR(21)
C?    RETURN
C?    ENDIF
      ENDIF
 
      IF(IARG.EQ.3)THEN
C 3ème COEF
      CALL ACMM(MTABX,'ARG3',NOMK)
C     write(6,*)' NOMK=',nomk
      IF(NOMK.EQ.'CAVIFERM')IKOMP=0
      IF(NOMK.EQ.'NOMUL   ')NOMUL=0
      ENDIF
 
      NOMP=LINCO.MOTS(1)
      NOMD=LINCO.MOTS(2)(1:4)
 
      NRIGE=7
      NKID =9
      NKMT =7
      NMATRI=1
      IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.2)NMATRI=2
      IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.4.AND.IDIM.EQ.2)NMATRI=2
      IF(KPRE.EQ.5)NMATRI=NMATRI+1
C     write(6,*)' NMATRI=',nmatri,' KPRE=',kpre
      SEGINI MATRIK
 
C CAS Stabilisation via MACRO CENTRE
      IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.2)THEN
      IF(IKAS.EQ.3)BETA0=-BETA0
C     write(6,*)' CAS Stabilisation via MACRO CENTRE ',BETA0,IKAS
      NK=0
      I2=2
      NBME=1
      NBSOUS=1
      SEGINI IMATRI
      IRIGEL(4,i2)=IMATRI
      KSPGP=MCTREI
      KSPGD=MCTREI
      IRIGEL(1,i2)=MELEMA
      IRIGEL(2,i2)=MELEMA
      IRIGEL(7,i2)=0
      CALL LICHT(MCHPOC,MPOVAL,TYPC,IGEOM)
 
C     write(6,*)' MELSTB=',melstb,' MCTREI=',MCTREI
 
      SEGACT MELSTB
 
      NBSOUS=MELSTB.LISOUS(/1)
      IF(NBSOUS.NE.0)THEN
      CALL ERREUR(5)
      ENDIF
 
      NBEL=MELEMA.NUM(/2)
      NBCI=MELSTB.NUM(/2)
      NP  =MELSTB.NUM(/1)
      MP  =NP
 
C     write(6,*)' nbel,nbci,np,mp=',nbel,nbci,np,mp
      SEGINI IZAFM
      LIZAFM(1,1)=IZAFM
      LISPRI(1)=NOMD
      LISDUA(1)=NOMD
 
      CALL KRIPAD(MCTREI,MLENT1)
 
C     write(6,*)' IK1=',ik1
      IF(IK1.LT.4)THEN
      DO 33 K=1,NBEL
 
      NK=NK+1
      KC=1+(1-IK1)*(NK-1)
C     write(6,*)' K=',K,' NK=',nk,np,kc
 
      DO 32 J=1,NP
      K1=MLENT1.LECT(MELEMA.NUM(J,K))
      II=J
      DO 321 I=1,NP
      U=VPOCHA(K1,I)*BETA0*IZTGG1.VPOCHA(KC,1)
      IF(I.EQ.1)U=ABS(VPOCHA(K1,I))*BETA0*IZTGG1.VPOCHA(KC,1)
      IF(II.LE.NP)THEN
      AM(K,II,J)=U
      ELSE
      AM(K,II-NP,J)=U
      ENDIF
      II=II+1
 321  CONTINUE
 32   CONTINUE
 33   CONTINUE
C     write(6,*)' FIN BCL 33 '
      ELSE
 
      DO 34 K=1,NBEL
 
      NK=NK+1
      UA=1.D0
 
      DO 35 J=1,NP
      K1=MLENT1.LECT(MELEMA.NUM(J,K))
      II=J
      DO 351 I=1,NP
      U=VPOCHA(K1,I)*BETA0*UA
      IF(I.EQ.1)U=ABS(VPOCHA(K1,I))*BETA0*UA
      IF(II.LE.NP)THEN
      AM(K,II,J)=U
      ELSE
      AM(K,II-NP,J)=U
      ENDIF
      II=II+1
 351  CONTINUE
 35   CONTINUE
 34   CONTINUE
 
 
      ENDIF
      SEGSUP MLENT1
      ENDIF
C     write(6,*)' Apres 34 '
 
C CAS Stabilisation via MACRO CENTREP1
      DEUPI=1.D0
      IF(IAXI.NE.0)DEUPI=2.D0*XPI
      IF(IKAS.EQ.3)DEUPI=-DEUPI
 
      IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.4.AND.IDIM.EQ.2)THEN
      I2=2
      NBME=1
      NBSOUS=1
      SEGINI IMATRI
      IRIGEL(4,i2)=IMATRI
      KSPGP=MELEMC
      KSPGD=MELEMC
      IRIGEL(1,i2)=MELEMP
      IRIGEL(2,i2)=MELEMP
      IRIGEL(7,i2)=0
 
      SEGACT MELEMP
      NBSOUS=MELEMP.LISOUS(/1)
C     write(6,*)' MELEMP2 NBSOUS=',NBSOUS,' MELEMP=',MELEMP
      IF(NBSOUS.NE.0)THEN
      CALL ERREUR(5)
      ENDIF
 
      SEGACT MELEMI
 
      NBEL=MELEMP.NUM(/2)
      NP  =MELEMP.NUM(/1)
      MP  =NP
      NOM0=NOMS(MELEMP.ITYPEL)//'    '
 
      SEGINI IZAFM
      LIZAFM(1,1)=IZAFM
      LISPRI(1)=NOMD
      LISDUA(1)=NOMD
 
      CALL KALPBG(NOM0,'FONFORM ',IZFFM)
 
      SEGACT IZFFM*MOD
      IZHR=KZHR(1)
      SEGACT IZHR*MOD
      NES=GR(/1)
      NPG=GR(/3)
 
      NK=0
      DO 430 LM=1,MAX(1,MELEMI.LISOUS(/1))
      IPT1=MELEMI
      IF(MELEMI.LISOUS(/1).NE.0)IPT1=MELEMI.LISOUS(LM)
      SEGACT IPT1
      ITYP=IPT1.ITYPEL
      BETA=0.D0
      IF(NOMS(ITYP).EQ.'TRI6')BETA=BETA0
C     write(6,*)' NOMS(ITYP)=',NOMS(ITYP),BETA
      NBEL=IPT1.NUM(/2)
 
C Cas coefficient FLOTTANT ou CENTRE
      IF(IK1.LT.4)THEN
      DO 43 KK=1,NBEL
      NK=NK+1
      DO    I=1,NP
      J=MELEMP.NUM(I,NK)
      DO    N=1,IDIM
      XYZ(N,I)=XCOOR((J-1)*(IDIM+1)    +N)
      ENDDO
      ENDDO
 
      CALL CALJBR
     &(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
      KC=1+(1-IK1)*(NK-1)
 
      DO    I=1,NP
      DO    J=1,NP
      U=0.D0
      DO 433 L=1,NPG
      UH=0.D0
      DO 4333 KH=1,IDIM
      UH=UH+HR(KH,J,L)*HR(KH,I,L)
 4333 CONTINUE
      U=U-UH*PGSQ(L)*DEUPI*BETA*IZTGG1.VPOCHA(KC,1)
 433  CONTINUE
      AM(NK,I,J)=U
      ENDDO
      ENDDO
 
 43   CONTINUE
 
      ELSEIF(IK1.EQ.4.OR.IK1.EQ.5)THEN
C Cas coefficient SOMMET
 
      DO 45 KK=1,NBEL
      NK=NK+1
      DO    I=1,NP
      J=MELEMP.NUM(I,NK)
      DO    N=1,IDIM
      XYZ(N,I)=XCOOR((J-1)*(IDIM+1)    +N)
      ENDDO
      ENDDO
 
      CALL CALJBR
     &(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
 
      DO    I=1,NP
      DO    J=1,NP
      U=0.D0
      DO 453 L=1,NPG
      UH=0.D0
      DO 4533 KH=1,IDIM
      UH=UH+HR(KH,J,L)*HR(KH,I,L)
 4533 CONTINUE
 
      UA=0.D0
      DO     JJ=1,NP
      J1=LECT(IPT1.NUM(JJ,NK))
      DO     KH=1,IDIM
      KG=(IK1-4)*(KH-1)+1
      UA=UA+FN(JJ,L)*IZTGG1.VPOCHA(J1,KG)
      ENDDO
      ENDDO
      UA=UA/KG
 
      U=U-UH*PGSQ(L)*DEUPI*BETA*UA
 453  CONTINUE
      AM(NK,I,J)=U
      ENDDO
      ENDDO
 
 45   CONTINUE
      ENDIF
 
      SEGDES IPT1
 430  CONTINUE
 
C      write(6,*)'C Fin Stabilisation de toutes sortes'
      ENDIF
C Fin Stabilisation de toutes sortes
 
      NBME=IDIM
      CALL KRIPAD(MELEMI,MLENTI1)
      SEGACT MELEMI
      NBSOUS=MELEMI.LISOUS(/1)
      IF(NBSOUS.EQ.0)NBSOUS=1
      SEGINI IMATRI
 
      IF(IKAS.EQ.2)THEN
      KSPGD=MELEMS
      KSPGP=MELEMC
      IRIGEL(2,1)=MELEMI
      IRIGEL(1,1)=MELEMP
      ELSE
      KSPGP=MELEMS
      KSPGD=MELEMC
      IRIGEL(1,1)=MELEMI
      IRIGEL(2,1)=MELEMP
      ENDIF
      SEGACT MELEMP
 
C     write(6,*)' ds kmac melemp=',IRIGEL(1,1)
C     write(6,*)' ds kmac melemd=',IRIGEL(2,1)
 
      IRIGEL(4,1)=IMATRI
      IF(NOMUL.EQ.0)THEN
      IRIGEL(7,1)=2
      ELSE
      IF(IKAS.EQ.1)IRIGEL(7,1)=-3
      IF(IKAS.EQ.2)IRIGEL(7,1)=3
      IF(IKAS.EQ.3)IRIGEL(7,1)=4
      ENDIF
 
      SEGACT MELEMP
 
      K0=0
      DO 11 L=1,MAX(1,MELEMI.LISOUS(/1))
      IPT1=MELEMI
      IPT2=MELEMP
      IF(MELEMI.LISOUS(/1).NE.0)IPT1=MELEMI.LISOUS(L)
      IF(MELEMP.LISOUS(/1).NE.0)IPT2=MELEMP.LISOUS(L)
      SEGACT IPT1,IPT2
      NBEL=IPT1.NUM(/2)
 
      IF(IKAS.EQ.2)THEN
      MP=IPT1.NUM(/1)
C avt NP=MELEMP.NUM(/1)
      NP=IPT2.NUM(/1)
      ELSE
      NP=IPT1.NUM(/1)
C avt MP=MELEMP.NUM(/1)
      MP=IPT2.NUM(/1)
      ENDIF
 
      DO 12 I=1,NBME
      SEGINI IZAFM
      LIZAFM(L,I)=IZAFM
      IF(IKAS.EQ.2)THEN
      WRITE(NOM,FMT='(I1,A3)')I,NOMD(1:3)
      LISDUA(I)=NOM//'    '
      LISPRI(I)=NOMP
      ELSE
      WRITE(NOM,FMT='(I1,A3)')I,NOMP(1:3)
      LISPRI(I)=NOM//'    '
      LISDUA(I)=NOMD
      ENDIF
 12   CONTINUE
      IPM1=LIZAFM(L,1)
      IPM2=LIZAFM(L,2)
      IPM3=LIZAFM(L,2)
      IF(IDIM.EQ.3)IPM3=LIZAFM(L,3)
 
      IF(IK1.LT.4)THEN
C     write(6,*)' Appel KPRUSS'
      CALL KPRUSS(IPT1,IPT2,
     &IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1,IK1,K0,IKOMP)
C     write(6,*)'Retour KPRUSS'
      ELSE
C     write(6,*)' Appel KPRJSS'
      CALL KPRJSS(IPT1,IPT2,IPM1,IPM2,IPM3,
     & IAXI,IKAS,INEFMD,KPRE,IZTGG1,MLENTI,IK1,IKOMP)
      ENDIF
 
      K0=K0+NBEL
      SEGDES IPM1,IPM2,IPM3
      SEGDES IPT1
 11   CONTINUE
      SEGDES MELEMI
      IF(IK1.EQ.4)SEGSUP MLENTI
 
 
C%%%%%%%%%%%%%%%%% CAS DES PRESSIONS CONTINUES %%%%%%%%%%%%%%%%%%%%%%%
 
      IF(KPRE.EQ.5)THEN
      CALL LEKTAB(MTABZ,'MMAIL',MELEME)
      CALL LEKTAB(MTABZ,'MSOMMET',MELEMS)
 
      SEGACT MELEME
      NBSOU1=LISOUS(/1)
      IF(NBSOU1.EQ.0)NBSOU1=1
 
      IRIGEL(1,2)=MELEME
      IRIGEL(2,2)=MELEME
      IRIGEL(7,2)=0
 
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
      IHV=0
      AIMPL=1.D0
C --- CALCUL DU NOMBRE DE PAQUETS DE LRV ELEMENTS ---
      CALL KPLRVM(MELEME,LRV,NBSOUS)
      NBME=1
 
      SEGINI IMATRI
      IRIGEL(4,2)=IMATRI
      KSPGP=MELEMS
      KSPGD=MELEMS
      LISPRI(1)=NOMD//'    '
      LISDUA(1)=NOMD//'    '
 
      LL=0
      NUTOEL=0
      SEGACT MELEME
 
      DO 101 L=1,MAX(1,LISOUS(/1))
      IPT1=MELEME
      IF(LISOUS(/1).NE.0)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)
      MP = NP
      NBELEM=IPT1.NUM(/2)
 
      NNN=MOD(NBELEM,LRV)
      IF(NNN.EQ.0) NPACK=NBELEM/LRV
      IF(NNN.NE.0) NPACK=1+(NBELEM-NNN)/LRV
 
      DO 701 KPACK=1,NPACK
 
C --- POUR CHAQUE PAQUET DE LRV ELEMENTS ou moins
         KDEB=1+(KPACK-1)*LRV
         KFIN=MIN(NBELEM,KDEB+LRV-1)
 
      NBEL=KFIN-KDEB+1
      LL=LL+1
 
      SEGINI IPM1,IPS1
C     write(6,*)' IPM1,LL=',IPM1,LL
      LIZAFM(LL,1)=IPM1
 
      IPM2=IPM1
      IPM3=IPM1
      IPS2=IPS1
      IPS3=IPS1
 
      KITT=2
      KJTT=IK1
      VISCO=IZTGG1
      CALL INITD(VISC,3, 1.D0)
 
      CALL XLAPL(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,
     & VISC,VISC,VISC,KITT,KJTT,IK1,
     & IPT1.NUM(1,KDEB),NBEL,NUTOEL,XCOOR,AIMPL,IKOMP,
     & IPM1.AM,IPM2.AM,IPM3.AM,
     & IPS1.AM,IPS2.AM,IPS3.AM,
     & NINKO,IHV,IARG,VISC)
 
      NUTOEL=NUTOEL+NBEL
      SEGDES IPM1
 
 701  CONTINUE
      SEGDES IPT1
 101  CONTINUE
 
 
 
 
 
 
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
      ENDIF
 
      CALL ECROBJ('MATRIK',MATRIK)
 
            NAT=2
            NSOUPO=0
            SEGINI MCHPOI
            JATTRI(1)=2
            CALL ECROBJ('CHPOINT',MCHPOI)
 
C     write(6,*)' Fin operateur KMIC'
      SEGDES IMATRI,MATRIK
      RETURN
 1001 FORMAT(20(1X,I5))
 1002 FORMAT(10(1X,1PE11.4))
      END
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

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