Sources Esope | Cast3M

matran
C MATRAN    SOURCE    GOUNAND   25/11/12    21:15:40     12399          
      SUBROUTINE MATRAN(MTABD,MCHPOI)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C------------------------------------------------------------------
C PROJET        :  CASTEM 2000
C
C NOM           :  MATRAN
C
C DESCRIPTION   :  Cette subroutine cree un chpoint 'FACE' sur MTABD
C                  contenant la matrice de passage du repaire global
C                  dans le repaire de chaque face (voir CALJQB)
C
C LANGAGE       :  FORTRAN 77 + ESOPE 2000
C
C AUTEUR        :  R. MOREL, DRN/DMT/SEMT/TTMF
C
C-----------------------------------------------------------------
C
C APPELES (E/S) :  LEKTAB
C
C APPELES (Calcul):  CALJBR
C
C-----------------------------------------------------------------
C
C ENTREES       :
C                MTABD  : Table de sous-type domaine
C
C SORTIES       :
C                MCHPOI : Pointeur du CHPOIN contenant les
C                Matrices decrite ci-dessus.
C
C-----------------------------------------------------------------
C
C HISTORIQUE (Anomalies et modifications eventuelles)
C
C HISTORIQUE     : 23.10.98, Creation
C
C-----------------------------------------------------------------
C
 
-INC PPARAM
-INC CCOPTIO
-INC CCGEOME
-INC SMTABLE
      POINTEUR MTABD.MTABLE
-INC SMELEME
      POINTEUR MELEMC.MELEME,MELEMD.MELEME
      POINTEUR MELEMF.MELEME
-INC SMCOORD
-INC SMCHPOI
-INC SMLENTI
-INC SIZFFB
      PARAMETER (NBO=5)
      REAL*8 CFT,NORMX,NORMY,NORMZ,A2J(2,2,1),T,TX,TY,TZ,AIRE
      REAL*8 X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,SCAL,SGN,A3J(3,3,1)
      INTEGER NP,NPG,NES
      CHARACTER*8 TYPE,TYPC,NOM0,LIST1(NBO),LIST2(NBO)
C Les elements complets ne sont pas tous implementes
C On se ramene aux elements facep sans le centre de la face
      DATA LIST1/'SEG3    ','TRI4    ','TRI7    ','QUA5    ',
     &           'QUA9    '/
      DATA LIST2/'SEG2    ','TRI3    ','TRI6    ','QUA4    ',
     &           'QUA8    '/
C
      CALL LEKTAB(MTABD,'FACEL',MELEMD)
      SEGACT MELEMD
      IPT2=MELEMD
      CALL LEKTAB(MTABD,'FACE',MELEMF)
      CALL KRIPAD(MELEMF,MLENTI)
      SEGACT MELEMF
      NC=IDIM*IDIM
      TYPE='FACE'
      CALL CRCHPT(TYPE,MELEMF,NC,1,MCHPOI)
      CALL LICHT(MCHPOI,MPOVAL,TYPC,IGEOM)
      CALL LEKTAB(MTABD,'FACEP',MELEME)
      SEGACT MELEME
C
      SEGACT MCHPOI*MOD
      MSOUPO=MCHPOI.IPCHP(1)
      SEGACT MSOUPO*MOD
      NBSOUS=LISOUS(/1)
      IF(NBSOUS.EQ.0)NBSOUS=1
      IF (IDIM.EQ.2) THEN
C
      NOCOMP(1)='UX  '
      NOCOMP(2)='UY  '
      NOCOMP(3)='RX  '
      NOCOMP(4)='RY  '
      DO 1 L=1,NBSOUS
      IPT1=MELEME
      IF(NBSOUS.NE.1)IPT1=LISOUS(L)
      SEGACT IPT1
      NP=IPT1.NUM(/1)-1
      NEL=IPT1.NUM(/2)
C  Toutes les faces sont des SEG3 donc NP=2
      IF(NP.NE.2)RETURN
C
C
C  On considere la face sans son centre comme un elt fini SEG2
C
      TYPE = NOMS(IPT1.ITYPEL)//'    '
      CALL OPTLI(IP,LIST1,TYPE,NBO)
      IF (IP .EQ. 0) CALL ERREUR(5)
      TYPE = LIST2(IP)
      CALL KALPBG(TYPE,'FONFORM0',IZFFM)
      SEGACT IZFFM*MOD
      IZHR=KZHR(1)
      SEGACT IZHR*MOD
      NES=1
      NPG=1
C Boucle sur toutes les faces pour le paquet L
      DO 10 K=1,NEL
      NF=LECT(IPT1.NUM(NP+1,K))
C     REMPLISSAGE DE XYZ
C     ------------------
C
      DO 12 I=1,NP
      J=IPT1.NUM(I,K)
      DO 121 N=1,IDIM
      XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N)
 121  CONTINUE
 12   CONTINUE
C
C   CALJBR calcul le jacobien du passage de l'elt de ref.
C   a l'elt. reel
C   A2J=Jacobien AIRE=detA2J
      CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,
     &NPG,IAXI,AIRE,A2J,SGEN)
C
      NORMX=A2J(1,2,1)
      NORMY=A2J(2,2,1)
C On verifie que n est dans le meme sens que vec(13) de FACEL
C Calcul de vec(13)
C FACEL est repere par IPT2
         J1=IPT2.NUM(1,NF)
         J2=IPT2.NUM(2,NF)
         JJ2=LECT(J2)
         J3=IPT2.NUM(3,NF)
C
         IF (J1.eq.J3) THEN
            X1=XCOOR((J1-1)*(IDIM+1)+1)
            Y1=XCOOR((J1-1)*(IDIM+1)+2)
            X2=XCOOR((J2-1)*(IDIM+1)+1)
            Y2=XCOOR((J2-1)*(IDIM+1)+2)
            SCAL=(X2-X1)*NORMX+(Y2-Y1)*NORMY
            SGN=SIGN(1.D0,SCAL)
C
         ELSE
            X1=XCOOR((J1-1)*(IDIM+1)+1)
            Y1=XCOOR((J1-1)*(IDIM+1)+2)
            X3=XCOOR((J3-1)*(IDIM+1)+1)
            Y3=XCOOR((J3-1)*(IDIM+1)+2)
            SCAL=(X3-X1)*NORMX+(Y3-Y1)*NORMY
            SGN=SIGN(1.D0,SCAL)
         ENDIF
         MPOVAL.VPOCHA(JJ2,1)=A2J(1,1,1)
         MPOVAL.VPOCHA(JJ2,2)=A2J(2,1,1)
         MPOVAL.VPOCHA(JJ2,3)=SGN*NORMX
         MPOVAL.VPOCHA(JJ2,4)=SGN*NORMY
C
C
 10   CONTINUE
C
      SEGDES IPT1
 1    CONTINUE
C
C
      ELSE
C
      NOCOMP(1)='UX  '
      NOCOMP(2)='UY  '
      NOCOMP(3)='UZ  '
      NOCOMP(4)='RX  '
      NOCOMP(5)='RY  '
      NOCOMP(6)='RZ  '
      NOCOMP(7)='MX  '
      NOCOMP(8)='MY  '
      NOCOMP(9)='MZ  '
      DO 21 L=1,NBSOUS
      IPT1=MELEME
      IF(NBSOUS.NE.1)IPT1=LISOUS(L)
      SEGACT IPT1
      NP=IPT1.NUM(/1)-1
      NEL=IPT1.NUM(/2)
C  Les elts complets ne sont pas implementes, on utilise
C  l'elt. face sans le centre
C
      TYPE = NOMS(IPT1.ITYPEL)//'    '
      CALL OPTLI(IP,LIST1,TYPE,NBO)
      IF (IP .EQ. 0) CALL ERREUR(5)
      TYPE = LIST2(IP)
C  Calcul des fonctions de forme sur l'elt. de reference
C
      CALL KALPBG(TYPE,'FONFORM0',IZFFM)
      SEGACT IZFFM*MOD
      IZHR=KZHR(1)
      SEGACT IZHR*MOD
      NES=2
      NPG=1
C  Boucle sur les faces du paquet L
      DO 210 K=1,NEL
      NF=LECT(IPT1.NUM(NP+1,K))
C
C     REMPLISSAGE DE XYZ
C     ------------------
C
C     write(6,*)'num',(ipt1.num(ii,k),ii=1,np+1),' NF=',nf
       DO 212 I=1,NP
       J=IPT1.NUM(I,K)
       DO 2121 N=1,IDIM
       XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N)
 2121 CONTINUE
 212   CONTINUE
C
C
C  calcul du jacobien idem dim 2
C
       CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,
     & NPG,IAXI,AIRE,A3J,SGEN)
       NORMX=A3J(1,3,1)
       NORMY=A3J(2,3,1)
       NORMZ=A3J(3,3,1)
C
C On verifie que n est dans le meme sens que vec(13) de FACEL
C Calcul de vec(13)
C FACEL est repere par IPT2
         J1=IPT2.NUM(1,NF)
         J2=IPT2.NUM(2,NF)
         JJ2=LECT(J2)
         J3=IPT2.NUM(3,NF)
C
         IF (J1.eq.J3) THEN
            X1=XCOOR((J1-1)*(IDIM+1)+1)
            Y1=XCOOR((J1-1)*(IDIM+1)+2)
            Z1=XCOOR((J1-1)*(IDIM+1)+3)
            X2=XCOOR((J2-1)*(IDIM+1)+1)
            Y2=XCOOR((J2-1)*(IDIM+1)+2)
            Z2=XCOOR((J2-1)*(IDIM+1)+3)
C
            SCAL=(X2-X1)*NORMX+(Y2-Y1)*NORMY+(Z2-Z1)*NORMZ
            SGN=SIGN(1.D0,SCAL)
C
         ELSE
            X1=XCOOR((J1-1)*(IDIM+1)+1)
            Y1=XCOOR((J1-1)*(IDIM+1)+2)
            Z1=XCOOR((J1-1)*(IDIM+1)+3)
            X3=XCOOR((J3-1)*(IDIM+1)+1)
            Y3=XCOOR((J3-1)*(IDIM+1)+2)
            Z3=XCOOR((J3-1)*(IDIM+1)+3)
            SCAL=(X3-X1)*NORMX+(Y3-Y1)*NORMY+(Z3-Z1)*NORMZ
            SGN=SIGN(1.D0,SCAL)
         ENDIF
C
         TX=A3J(1,1,1)
         TY=A3J(2,1,1)
         TZ=A3J(3,1,1)
         T=SQRT(TX*TX+TY*TY+TZ*TZ)
         TX=TX / T
         TY=TY / T
         TZ=TZ / T
         MPOVAL.VPOCHA(JJ2,1)=TX
         MPOVAL.VPOCHA(JJ2,2)=TY
         MPOVAL.VPOCHA(JJ2,3)=TZ
         MPOVAL.VPOCHA(JJ2,4)=NORMY*TZ-NORMZ*TY
         MPOVAL.VPOCHA(JJ2,5)=NORMZ*TX-NORMX*TZ
         MPOVAL.VPOCHA(JJ2,6)=NORMX*TY-NORMY*TX
         MPOVAL.VPOCHA(JJ2,7)=SGN*NORMX
         MPOVAL.VPOCHA(JJ2,8)=SGN*NORMY
         MPOVAL.VPOCHA(JJ2,9)=SGN*NORMZ
C
 210  CONTINUE
C
      SEGDES IPT1
 21   CONTINUE
      ENDIF
      CALL ECMO(MTABD,'MATROT','CHPOINT ',MCHPOI)
      SEGSUP MLENTI
      SEGDES MELEMD
      SEGDES MELEME
      SEGDES MELEMF
      SEGDES MCHPOI
      SEGDES MSOUPO
      SEGDES MPOVAL
      RETURN
      END