Sources Esope | Cast3M

kdxdy
C KDXDY     SOURCE    GOUNAND   25/11/12    21:15:18     12399          
      SUBROUTINE KDXDY
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C*************************************************************************
C
C     OBJET   : Cree un CHAMPOINT CENTRE contenant les informations
C               décrites ci dessous, sur les éléments du domaine.
C
C     SYNTAXE : CHPC = KDXDY OBJDOM ;
C
C                      OBJDOM : TABLE de SOUSTYPE DOMAINE
C
C     AVIS IMPORTANT ;
C
C     Ces informations sont destinées aux opérateurs de discrétisation
C     et sont rangées dans un ordre particulier.
C
C Les dimensions du tableau sont : (DIME NBEL)  DIME=11 en 2D
C                                               DIME=13 en 3D
C et il contient :
C
C  pour un SEG2 la longueur de l'element (XML) et la matrice P(2,2) en 2D
C                                                            P(3,3) en 3D
C  pour un TRI3 XML,XMH,AJ1/XML,AJ2/XML, et la matrice P(2,2) en 2D
C                                                      P(3,3) en 3D
C  pour un QUA4 IDEM
C  pour un CUB8 XML XMH XME et la matrice P
C  pour un PRI6 IDEM
C
C***********************************************************************
C   MATRICE P
C       LA MATRICE DE ROTATION DU REPERE GLOBALE VERS LE REPERE LOCAL
C   DEFINI PAR DEUX OU TROIS POINTS PRIS DANS XYZ SUIVANT QU'ON EST
C   EN 2D OU EN 3D
C   ON PREND P1 P2 ET PNP
C
C     U                TEL QUE T SOIT DIRIGE SUIVANT P1P2 V TOURNE VERS
C     .     .V          P1PNP ET U = T VECTORIEL V
C     .   .
C     . .                             __       __
C (P1). . . . .T (P2)                | tx ty tz  |
C                                    |           |
C  ON A ALORS WL= P WG          P  = | vx vy vz  |
C                                    |           |
C                                    | ux uy uz  |
C                                    |__       __|
C*************************************************************************
 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
*-
-INC CCGEOME
-INC SMTABLE
-INC SMELEME
      POINTEUR MELEMC.MELEME
-INC SMCOORD
-INC SMCHPOI
-INC SIZFFB
      PARAMETER (NLTS=10)
      CHARACTER*8 LISTS(NLTS)
C***
      DIMENSION AAJ(3,3,9),XC(3),T(13),DELTAX(3)
      PARAMETER (NTB=1)
      CHARACTER*8 LTAB(NTB),TYPE,NOM0,TYPC
      DIMENSION KTAB(NTB)
      DATA LTAB/'DOMAINE '/
      DATA LISTS/'QUA4    ','CUB8    ','TRI3    ','PRI6    ',
     &'SEG2    ','TET4    ','PYR5    ','TRI7    ','QUA9    ',
     &'SEG3    '/
C***
      NTO=NTB
      CALL LITABS(LTAB,KTAB,NTB,NTO,IRET)
      IF(IRET.EQ.0)RETURN
      MTABD=KTAB(1)
 
      TYPE=' '
      CALL ACMO(MTABD,'CENTRE',TYPE,MELEMC)
 
      NC=IDIM
 
      TYPE='CENTRE'
      CALL CRCHPT(TYPE,MELEMC,NC,1,MCHPOI)
      CALL LICHTM(MCHPOI,MPOVAL,TYPC,IGEOM)
 
      TYPE=' '
      CALL ACMO(MTABD,'MAILLAGE',TYPE,MELEME)
      IF(TYPE.NE.'MAILLAGE')GO TO 90
      SEGACT MELEME
 
      IAXI=0
      IF(IFOMOD.EQ.0)IAXI=2
      NPGI=1
 
      NBSOUS=LISOUS(/1)
      IF(NBSOUS.EQ.0)NBSOUS=1
 
      K0=0
      DO 1 L=1,NBSOUS
      IPT1=MELEME
      IF(NBSOUS.NE.1)IPT1=LISOUS(L)
      SEGACT IPT1
      NP=IPT1.NUM(/1)
      NEL=IPT1.NUM(/2)
      NOM0=NOMS(IPT1.ITYPEL)//'    '
C     write(6,*)' KDXDY : NOM0=',nom0,' NEL=',nel
      CALL OPTLI(IP,LISTS,NOM0,NLTS)
      IF(IP.EQ.0)THEN
      WRITE(6,*)' Sub KDXDY :'
      WRITE(6,*)' Il y a des types d''éléments pour lesquels on ne sait'
     &,'pas faire'
      GO TO 1
      ENDIF
      CALL KALPBG(NOM0,'FONFORM ',IZFFM)
      SEGACT IZFFM*MOD
      IZHR=KZHR(1)
      SEGACT IZHR*MOD
      NES=GR(/1)
      NPG=GR(/3)
 
      DO 10 K=1,NEL
      KK=K0+K
      NPGR=0
      IF(IAXI.NE.0)NPGR=NPG
C
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
 
      CALL CALJAJ(GR,XYZ,NES,IDIM,NP,NPGI,AAJ)
C
      CALL CALJBC(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,AIRE)
 
      IF(AIRE.LE.0.D0)AIRE=ABS(AIRE)
 
 
C             CALCUL DES LONGUEURS DES COTES DE L'ELEMENT
C                    UNIQUEMENT POUR LES SEG2 QUA4 ET CUB8
 
      XML=0.
      DO 231 N=1,IDIM
      XML=XML+(XYZ(N,1)-XYZ(N,2))*(XYZ(N,1)-XYZ(N,2))
 231  CONTINUE
      XML=SQRT(XML)
 
      DO 234 N=1,IDIM
      DELTAX(N)=ABS(XYZ(N,NP)-XYZ(N,1))
      DO 233 I=2,NP
      DELTAX(N)=DELTAX(N)+ABS(XYZ(N,I)-XYZ(N,I-1))
 233  CONTINUE
      VPOCHA(KK,N)=(DELTAX(N)/FLOAT(NP))+XPETIT
      IF(NOM0.EQ.'QUA4    ')THEN
      VPOCHA(KK,N)=VPOCHA(KK,N)*FLOAT(IDIM)
      ENDIF
 234  CONTINUE
 
C
C
      IF(XML.EQ.0.)
     &WRITE(6,*)' ELEMENT DEGENERE ( AYANT UN COTE DE LONGUEUR NULLE) '
C
C
 10   CONTINUE
      K0=K0+NEL
      SEGDES IPT1
 1    CONTINUE
      SEGDES MELEME,MPOVAL,MCHPOI
C
      CALL ECROBJ('CHPOINT ',MCHPOI)
 
      RETURN
 
 90   CONTINUE
      WRITE(6,*)' Interruption anormale de KDXDY'
      RETURN
 1001 FORMAT(20(1X,I5))
 1002 FORMAT(10(1X,1PE11.4))
      END