Sources Esope | Cast3M

kdom11
C KDOM11    SOURCE    PV        21/10/31    21:15:01     11158          
      SUBROUTINE KDOM11(MELF,MELFL,MELFP,MCHPSU,MCHPNO,MCHPMR)
C
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  KDOM11
C
C DESCRIPTION       :  Subroutine called by KDOM10 in the case of EULER
C                      model.
C                      We create the CHPOINTS
C                      MCHPSU =  surfaces (faces dimension)
C                      MCHPNO =  normals (oriented from the first
C                                to the second point of FACEL)
C                      MCHPMR =  rotation matrix
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/LTMF
C
C************************************************************************
C
C INPUTS :
C
C MELF   : meleme 'FACE'
C MELFL  : meleme 'FACEL'
C MELFP  : meleme 'FACEP'
C
C OUTPUTS :
C
C MCHPSU : mchpoi 'XXSURFAC'
C MCHPNO : mchpoi 'XXNORMAF'
C MCHPMR : mchpoi 'MATROT'
C
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C
      INTEGER IGEOM, MCHPSU, MCHPNO,  MCHPMR
     &     ,JGN, JGM, NBSOUS, ISOUS, NP, NEL, IEL, IPO, NLCF
     &     , NF, IFAC
      REAL*8 XF,YF,X1,Y1, DSURFX, DSURFY, XSURF, YSURF, SURF, A
     &     ,ORIENT
     &     ,ZF, DXP, DYP, DZP, DXPM1, DYPM1, DZPM1, DSURFZ, ZSURF
     &     ,RMX, RMY, RMZ, RRX, RRY, RRZ
      CHARACTER*8 TYPE
 
 
-INC PPARAM
-INC CCOPTIO
-INC SMELEME
      POINTEUR MELFL.MELEME,MELFP.MELEME,MELF.MELEME
-INC SMCHPOI
      POINTEUR MPOVSU.MPOVAL, MPOVNO.MPOVAL, MPOVMR.MPOVAL
-INC SMLENTI
-INC SMLMOTS
-INC SMCOORD
-INC CCREEL
C
C**** Corresp. FACE
C
      CALL KRIPAD(MELF,MLENTI)
C     SEGINI MLENTI
C
C**** Champoint surfaces
C
      JGN=4
      JGM=1
      SEGINI MLMOTS
      MLMOTS.MOTS(1)='SCAL'
      TYPE='FACE    '
      CALL KRCHP1(TYPE,MELF,MCHPSU,MLMOTS)
      IF(IERR.NE.0)GOTO 9999
      CALL LICHT(MCHPSU,MPOVSU,TYPE,IGEOM)
      IF(IERR.NE.0)GOTO 9999
C     SEGACT MPOVSU
      SEGSUP MLMOTS
C
C**** Champoint normales
C
      JGN=4
      JGM=IDIM
      SEGINI MLMOTS
      MLMOTS.MOTS(1)='UX'
      MLMOTS.MOTS(2)='UY'
      IF(IDIM .EQ. 3) MLMOTS.MOTS(3)='UZ'
      TYPE='FACE    '
      CALL KRCHP1(TYPE,MELF,MCHPNO,MLMOTS)
      IF(IERR.NE.0)GOTO 9999
      CALL LICHT(MCHPNO,MPOVNO,TYPE,IGEOM)
      IF(IERR.NE.0)GOTO 9999
C     SEGACT MPOVNO
      SEGSUP MLMOTS
C
C**** Champoint matrice de rotation
C
      JGN=4
      JGM=IDIM*IDIM
      SEGINI MLMOTS
      IF(IDIM.EQ.2)THEN
         MLMOTS.MOTS(1)='RX'
         MLMOTS.MOTS(2)='RY'
         MLMOTS.MOTS(3)='MX'
         MLMOTS.MOTS(4)='MY'
*        Normale en M
*        vect(M,U) = z
      ELSE
         MLMOTS.MOTS(1)='UX'
         MLMOTS.MOTS(2)='UY'
         MLMOTS.MOTS(3)='UZ'
         MLMOTS.MOTS(4)='RX'
         MLMOTS.MOTS(5)='RY'
         MLMOTS.MOTS(6)='RZ'
         MLMOTS.MOTS(7)='MX'
         MLMOTS.MOTS(8)='MY'
         MLMOTS.MOTS(9)='MZ'
*        Normale en M
*        vect(M,R) = U
      ENDIF
      CALL KRCHP1(TYPE,MELF,MCHPMR,MLMOTS)
      IF(IERR.NE.0)GOTO 9999
      CALL LICHT(MCHPMR,MPOVMR,TYPE,IGEOM)
      IF(IERR.NE.0)GOTO 9999
C     SEGACT MPOVMR
C
      SEGACT MELFP
      NBSOUS=MELFP.LISOUS(/1)
      IF(NBSOUS.EQ.0)NBSOUS=1
      SEGACT MELFL
C**********************
C CAS DE LA DIMENSION 2
C**********************
      surf=0.d0
      xsurf=0.d0
      ysurf=0.d0
      zsurf=0.d0
      IF (IDIM.EQ.2) THEN
C
         DO ISOUS=1,NBSOUS,1
            IPT1 = MELFP
            IF(NBSOUS .NE. 1)THEN
               IPT1 = MELFP.LISOUS(ISOUS)
               SEGACT IPT1
            ENDIF
C
            NP=IPT1.NUM(/1)-1
            NEL=IPT1.NUM(/2)
            IF(NP .NE. 2)THEN
               WRITE(IOIMP,*) 'Subroutine knrf.eso'
               CALL ERREUR(5)
            ENDIF
C
            DO IEL=1,NEL,1
               NF=IPT1.NUM(NP+1,IEL)
               XF = XCOOR((NF-1)*(IDIM+1)+1)
               YF = XCOOR((NF-1)*(IDIM+1)+2)
*
               IPO=IPT1.NUM(1,IEL)
               X1=XCOOR((IPO-1)*(IDIM+1)+1)
               Y1=XCOOR((IPO-1)*(IDIM+1)+2)
               DSURFX=Y1-YF
               DSURFY=XF-X1
               XSURF=DSURFX
               YSURF=DSURFY
               SURF=((DSURFX*DSURFX)+(DSURFY*DSURFY))**0.5D0
*
               IPO=IPT1.NUM(2,IEL)
               X1=XCOOR((IPO-1)*(IDIM+1)+1)
               Y1=XCOOR((IPO-1)*(IDIM+1)+2)
               DSURFX=YF-Y1
               DSURFY=X1-XF
               XSURF=XSURF+DSURFX
               YSURF=YSURF+DSURFY
               SURF=SURF+(((DSURFX*DSURFX)+(DSURFY*DSURFY))**0.5D0)
C
               NLCF=MLENTI.LECT(NF)
               MPOVSU.VPOCHA(NLCF,1)=SURF
C
C************* Orientation selon FACEL
C
               IPO=MELFL.NUM(1,NLCF)
               X1=XCOOR((IPO-1)*(IDIM+1)+1)
               Y1=XCOOR((IPO-1)*(IDIM+1)+2)
               DSURFX=XF-X1
               DSURFY=YF-Y1
               ORIENT=SIGN(1.0D0,((DSURFX*XSURF)+(DSURFY*YSURF)))
C
               MPOVNO.VPOCHA(NLCF,1)=XSURF/SURF*ORIENT
               MPOVNO.VPOCHA(NLCF,2)=YSURF/SURF*ORIENT
C
               MPOVMR.VPOCHA(NLCF,3)=XSURF/SURF*ORIENT
               MPOVMR.VPOCHA(NLCF,4)=YSURF/SURF*ORIENT
               MPOVMR.VPOCHA(NLCF,1)=-1*YSURF/SURF*ORIENT
               MPOVMR.VPOCHA(NLCF,2)=XSURF/SURF*ORIENT
C
            ENDDO
            IF(NBSOUS .NE. 1) SEGDES IPT1
         ENDDO
C**********************
C CAS DE LA DIMENSION 3
C**********************
      ELSE
C
         DO ISOUS=1,NBSOUS,1
            IPT1 = MELFP
            IF(NBSOUS .NE. 1)THEN
               IPT1 = MELFP.LISOUS(ISOUS)
               SEGACT IPT1
            ENDIF
C
            NP=IPT1.NUM(/1)-1
            NEL=IPT1.NUM(/2)
C
            DO IEL=1,NEL,1
               NF=IPT1.NUM(NP+1,IEL)
               XF = XCOOR((NF-1)*(IDIM+1)+1)
               YF = XCOOR((NF-1)*(IDIM+1)+2)
               ZF = XCOOR((NF-1)*(IDIM+1)+3)
*
               IPO=IPT1.NUM(NP,IEL)
               DXP = XCOOR((IPO-1)*(IDIM+1)+1) - XF
               DYP = XCOOR((IPO-1)*(IDIM+1)+2) - YF
               DZP = XCOOR((IPO-1)*(IDIM+1)+3) - ZF
               XSURF=0.0D0
               YSURF=0.0D0
               ZSURF=0.0D0
               DO IFAC=1,NP,1
                  DXPM1 = DXP
                  DYPM1 = DYP
                  DZPM1 = DZP
                  IPO=IPT1.NUM(IFAC,IEL)
                  DXP = XCOOR((IPO-1)*(IDIM+1)+1) - XF
                  DYP = XCOOR((IPO-1)*(IDIM+1)+2) - YF
                  DZP = XCOOR((IPO-1)*(IDIM+1)+3) - ZF
                  DSURFX = 0.5D0 * ((DYPM1 * DZP) - (DZPM1 * DYP))
                  DSURFY = 0.5D0 * ((DZPM1 * DXP) - (DXPM1 * DZP))
                  DSURFZ = 0.5D0 * ((DXPM1 * DYP) - (DYPM1 * DXP))
                  XSURF=XSURF+DSURFX
                  YSURF=YSURF+DSURFY
                  ZSURF=ZSURF+DSURFZ
               ENDDO
*
               SURF=(XSURF*XSURF)+(YSURF*YSURF)+(ZSURF*ZSURF)
               surf=max(surf,xpetit)
               SURF=sqrt(SURF)
               NLCF=MLENTI.LECT(NF)
               MPOVSU.VPOCHA(NLCF,1)=SURF
C
C************* Orientation selon FACEL
C
               IPO=MELFL.NUM(1,NLCF)
               DSURFX=XF-XCOOR((IPO-1)*(IDIM+1)+1)
               DSURFY=YF-XCOOR((IPO-1)*(IDIM+1)+2)
               DSURFZ=ZF-XCOOR((IPO-1)*(IDIM+1)+3)
               ORIENT=SIGN(1.0D0,((DSURFX*XSURF)+(DSURFY*YSURF)+(DSURFZ
     &              *ZSURF)))
C
               MPOVNO.VPOCHA(NLCF,1)=XSURF/SURF*ORIENT
               MPOVNO.VPOCHA(NLCF,2)=YSURF/SURF*ORIENT
               MPOVNO.VPOCHA(NLCF,3)=ZSURF/SURF*ORIENT
C
C************* MATROT
C
C              Normal
C
               MPOVMR.VPOCHA(NLCF,7)=XSURF/SURF*ORIENT
               MPOVMR.VPOCHA(NLCF,8)=YSURF/SURF*ORIENT
               MPOVMR.VPOCHA(NLCF,9)=ZSURF/SURF*ORIENT
 
C              First direction (RX,RY,RZ) is normal to the
C              normal and FP (P = first point of FACEP)
C
               IPO=IPT1.NUM(1,IEL)
               DXP = XCOOR((IPO-1)*(IDIM+1)+1) - XF
               DYP = XCOOR((IPO-1)*(IDIM+1)+2) - YF
               DZP = XCOOR((IPO-1)*(IDIM+1)+3) - ZF
               DSURFX = (ZSURF * DYP) - (YSURF * DZP)
               DSURFY = (XSURF * DZP) - (ZSURF * DXP)
               DSURFZ = (YSURF * DXP) - (XSURF * DYP)
C
C              DZP=modulus of the RX,RY,RZ
C
               DZP=(((DSURFX*DSURFX)+(DSURFY*DSURFY)+
     &                 (DSURFZ*DSURFZ))**0.5D0)
C
               MPOVMR.VPOCHA(NLCF,4)=DSURFX/DZP
               MPOVMR.VPOCHA(NLCF,5)=DSURFY/DZP
               MPOVMR.VPOCHA(NLCF,6)=DSURFZ/DZP
C
C              (UX,UY,UZ,RX,RY,RZ,MX,MY,MZ)=(1,2,3,4,5,6,7,8,9)
C              M,R,U is a right-hand normal frame
C
               RMX=MPOVMR.VPOCHA(NLCF,7)
               RMY=MPOVMR.VPOCHA(NLCF,8)
               RMZ=MPOVMR.VPOCHA(NLCF,9)
               RRX=MPOVMR.VPOCHA(NLCF,4)
               RRY=MPOVMR.VPOCHA(NLCF,5)
               RRZ=MPOVMR.VPOCHA(NLCF,6)
C
               MPOVMR.VPOCHA(NLCF,1)=(RMY*RRZ) - (RMZ*RRY)
               MPOVMR.VPOCHA(NLCF,2)=(RMZ*RRX) - (RMX*RRZ)
               MPOVMR.VPOCHA(NLCF,3)=(RMX*RRY) - (RMY*RRX)
            ENDDO
            IF(NBSOUS .NE. 1) SEGDES IPT1
         ENDDO
      ENDIF
C
      SEGDES MPOVSU
      SEGDES MPOVNO
      SEGDES MPOVMR
      SEGDES MELFP
C
      SEGDES MELFL
      SEGSUP MLENTI
C
 9999 RETURN
      END