Sources Esope | Cast3M

kcot
C KCOT      SOURCE    GOUNAND   25/11/12    21:15:16     12399          
      SUBROUTINE KCOT
      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 = KCOT 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 CCGEOME
-INC SMTABLE
-INC SMELEME
      POINTEUR MELEMC.MELEME
-INC SMCOORD
-INC SMCHPOI
-INC SIZFFB
      PARAMETER (NLTS=8)
      CHARACTER*8 LISTS(NLTS)
C***
      REAL*8 AAJ(3,3,9),U,XC(3),T(13),XQ(81)
      PARAMETER (NTB=1)
      CHARACTER*8 LTAB(NTB),TYPE,NOM0,TYPC
      DIMENSION KTAB(NTB)
      DATA LTAB/'DOMAINE '/
      DATA LISTS/'QUA4    ','CUB8    ','TRI3    ','PRI6    ','SEG2    ',
     &'SEG3    ','QUA9    ','TRI7    '/
C***
C     write(6,*)' DBUT KCOT '
      NTO=NTB
      CALL LITABS(LTAB,KTAB,NTB,NTO,IRET)
      IF (IERR.NE.0) RETURN
      MTABD=KTAB(1)
 
      CALL INITD(AAJ,81,0.D0)
      TYPE=' '
      CALL ACMO(MTABD,'CENTRE',TYPE,MELEMC)
      IF(TYPE.NE.'MAILLAGE')GO TO 90
 
      IF(IDIM.EQ.2)NC=11
      IF(IDIM.EQ.3)NC=13
 
      TYPE='CENTRE'
      CALL CRCHPT(TYPE,MELEMC,NC,1,MCHPOI)
      CALL LICHT(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)//'    '
      CALL OPTLI(IP,LISTS,NOM0,NLTS)
      IF(IP.EQ.0)THEN
      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
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
 
      IF(IP.EQ.6)THEN
      IN=0
      DO 13 I=1,NP,2
         J=IPT1.NUM(I,K)
         DO 131 N=1,IDIM
            IN=IN+1
            XQ(IN)=XCOOR((J-1)*(IDIM+1)+N)
 131     CONTINUE
 13   CONTINUE
      ELSEIF(IP.GT.6)THEN
      IN=0
      DO 14 I=1,NP-1,2
         J=IPT1.NUM(I,K)
         DO 141 N=1,IDIM
            IN=IN+1
            XQ(IN)=XCOOR((J-1)*(IDIM+1)+N)
 141     CONTINUE
 14   CONTINUE
      ENDIF
 
      IF(IP.LE.5)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)AIRE=ABS(AIRE)
 
 
C             CALCUL DES LONGUEURS DES COTES DE L'ELEMENT
C                    UNIQUEMENT POUR LES SEG2 QUA4 ET CUB8
 
      XML=0.D0
      DO 231 N=1,IDIM
      XML=XML+(XYZ(N,1)-XYZ(N,2))*(XYZ(N,1)-XYZ(N,2))
 231  CONTINUE
      XML=SQRT(XML)
C
C
      IF(XML.EQ.0.D0)
     &WRITE(6,*)' ELEMENT DEGENERE ( AYANT UN COTE DE LONGUEUR NULLE) '
C
C
      IF(NOM0.EQ.'SEG2    ')THEN
      T(1)=XML
      CALL CALJQB(XYZ,T(5),IDIM,IDIM)
      GOTO 244
      ELSEIF(NOM0.EQ.'SEG3    ')THEN
      T(1)=XML
      CALL CALJQB(XQ,T(5),IDIM,IDIM)
      GOTO 244
      ENDIF
C
C
      XMH=0.D0
      DO 232 N=1,IDIM
      XMH=XMH+(XYZ(N,2)-XYZ(N,3))*(XYZ(N,2)-XYZ(N,3))
 232  CONTINUE
      XMH=SQRT(XMH)
C     WRITE(6,*)'........XML XMH NP NES ',XML,XMH,NP,NES
 
C
C
      IF(XMH.EQ.0.D0)
     &WRITE(6,*)' ELEMENT DEGENERE ( AYANT UN COTE DE LONGUEUR NULLE) '
C
C
 
      IF(NOM0.EQ.'QUA4    '.OR.NOM0.EQ.'TRI3    ')THEN
      T(1)=XML
      T(2)=XMH
      T(3)=AAJ(1,1,1)/XML
      T(4)=AAJ(2,1,1)/XML
      CALL CALJQB(XYZ,T(5),IDIM,IDIM)
      GOTO 244
      ELSEIF(NOM0.EQ.'QUA9    '.OR.NOM0.EQ.'TRI7    ')THEN
      T(1)=XML
      T(2)=XMH
      T(3)=AAJ(1,1,1)/XML
      T(4)=AAJ(2,1,1)/XML
      CALL CALJQB(XQ,T(5),IDIM,IDIM)
      GOTO 244
      ENDIF
 
      IF(NOM0.EQ.'CUB8    '.OR.NOM0.EQ.'PRI6    ')THEN
 
      XME=0.D0
      IF(NOM0.EQ.'CUB8    ')THEN
      DO 242 N=1,IDIM
      XME=XME+(XYZ(N,1)-XYZ(N,5))*(XYZ(N,1)-XYZ(N,5))
 242  CONTINUE
      ENDIF
 
      IF(NOM0.EQ.'PRI6    ')THEN
      DO 243 N=1,IDIM
      XME=XME+(XYZ(N,1)-XYZ(N,4))*(XYZ(N,1)-XYZ(N,4))
 243  CONTINUE
      ENDIF
 
      XME=SQRT(XME)
C     WRITE(6,*)'........XML XMH XME NP NES ',XML,XMH,XME,NP,NES
C
      IF(XME.EQ.0.D0)
     &WRITE(6,*)' ELEMENT DEGENERE ( AYANT UN COTE DE LONGUEUR NULLE) '
C
C
      CALL CALJQB(XYZ,T(5),IDIM,4)
      T(1)=XML
      T(2)=XMH
      T(3)=XME
      GO TO 244
      ENDIF
 
 244  CONTINUE
 
      CALL RSETIV(VPOCHA,NC,KK,T,NC)
 
C     write(6,*)' KK=',KK
C     write(6,1002)T
 10   CONTINUE
      K0=K0+NEL
      SEGDES IPT1
 1    CONTINUE
      SEGDES MELEME,MPOVAL,MCHPOI
C
      CALL ECROBJ('CHPOINT ',MCHPOI)
 
C     write(6,*)' FIN  KCOT '
      RETURN
 
 90   CONTINUE
      WRITE(6,*)' Interruption anormale de KCOT'
C        Option %m1:8 incompatible avec les données
      CALL ERREUR(803)
      RETURN
 1001 FORMAT(20(1X,I5))
 1002 FORMAT(10(1X,1PE11.4))
      END