Numérotation des lignes :

C KDOM7     SOURCE    KK2000    14/04/10    21:15:14     8032      SUBROUTINE KDOM7(IP1,IP2,IP3,IP4,XCEN)CC************************************************************************CC PROJET            :  CASTEM 2000CC NOM               :  KDOM7CC DESCRIPTION       :  Subroutine called by LEKMA0C                      Given a QUA4 (IP1,IP2,IP3,IP4), this subroutineC                      compute its center by averaging the centers ofC                      the two possible configuration it can beC                      divided inCC LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)CC AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/LTMFCC************************************************************************CC       IMPLICIT INTEGER(I-N)      IMPLICIT REAL*8(A-H,O-Z)       INTEGER IP1,IP2,IP3, IP4, I1,I2, IPCEN, IPL, IPR      REAL*8 VOL, XCEN(3),VOLT,DX0,DY0,DZ0,DX1,DY1,DZ1,CEL     &     ,RNOX0,RNOY0,RNOZ0,RNOX1,RNOY1,RNOZ1,XCEN1(3),VOLT1,PSCA     &     ,RNOT      LOGICAL LOGREG-INC CCOPTIO-INC SMCOORDC      DZ0=0.0D0      DZ1=0.0D0      RNOT=0.0D0CC**** We chose as positive direction the one whichC     has the maximum modulus IP1IP(1+1) x IP1IP(1+3)C      DO I1=1,4         IF(I1 .EQ. 1)THENC           We consider the triangles IP1, IP2, IP4            IPCEN=IP1            IPR=IP2            IPL=IP4         ELSEIF(I1 .EQ. 2)THEN            IPCEN=IP2            IPR=IP3            IPL=IP1         ELSEIF(I1 .EQ. 3)THEN            IPCEN=IP3            IPR=IP4            IPL=IP2         ELSEIF(I1 .EQ. 4)THEN            IPCEN=IP4            IPR=IP1            IPL=IP3         ENDIF         DX0=XCOOR((IPR-1)*(IDIM+1)+1)-     &        XCOOR((IPCEN-1)*(IDIM+1)+1)         DY0=XCOOR((IPR-1)*(IDIM+1)+2)-     &        XCOOR((IPCEN-1)*(IDIM+1)+2)         IF(IDIM. EQ. 3) DZ0=XCOOR((IPR-1)*(IDIM+1)+3)-     &        XCOOR((IPCEN-1)*(IDIM+1)+3)         DX1=XCOOR((IPL-1)*(IDIM+1)+1)-     &        XCOOR((IPCEN-1)*(IDIM+1)+1)         DY1=XCOOR((IPL-1)*(IDIM+1)+2)-     &        XCOOR((IPCEN-1)*(IDIM+1)+2)         IF(IDIM. EQ. 3) DZ1=XCOOR((IPL-1)*(IDIM+1)+3)-     &        XCOOR((IPCEN-1)*(IDIM+1)+3)         RNOX1=(DY0*DZ1-DY1*DZ0)         RNOY1=(DZ0*DX1-DZ1*DX0)         RNOZ1=(DX0*DY1-DX1*DY0)         VOL=((RNOX1**2)+(RNOY1**2)+(RNOZ1**2))**0.5D0         VOL=VOL/0.5D0         IF(VOL .GT. RNOT)THEN            RNOT=VOL            RNOX0=RNOX1            RNOY0=RNOY1            RNOZ0=RNOZ1         ENDIF      ENDDOCC**** First division in trianglesC     We consider the triangles IP3, IP1, IP2 andC     IP4, IP1, IP3C      XCEN1(1)=0.0D0      XCEN1(2)=0.0D0      XCEN1(3)=0.0D0      XCEN(1)=0.0D0      XCEN(2)=0.0D0      XCEN(3)=0.0D0      VOLT=0.0D0      VOLT1=0.0D0      LOGREG=.TRUE.C      DO I1=1,2         IF(I1 .EQ. 1)THENC           We consider the triangles IP3, IP1, IP2            IPCEN=IP1            IPR=IP2            IPL=IP3         ELSEIF(I1 .EQ. 2)THEN            IPCEN=IP1            IPR=IP3            IPL=IP4         ENDIF         DX0=XCOOR((IPR-1)*(IDIM+1)+1)-     &        XCOOR((IPCEN-1)*(IDIM+1)+1)         DY0=XCOOR((IPR-1)*(IDIM+1)+2)-     &        XCOOR((IPCEN-1)*(IDIM+1)+2)         IF(IDIM. EQ. 3) DZ0=XCOOR((IPR-1)*(IDIM+1)+3)-     &        XCOOR((IPCEN-1)*(IDIM+1)+3)         DX1=XCOOR((IPL-1)*(IDIM+1)+1)-     &        XCOOR((IPCEN-1)*(IDIM+1)+1)         DY1=XCOOR((IPL-1)*(IDIM+1)+2)-     &        XCOOR((IPCEN-1)*(IDIM+1)+2)         IF(IDIM. EQ. 3) DZ1=XCOOR((IPL-1)*(IDIM+1)+3)-     &        XCOOR((IPCEN-1)*(IDIM+1)+3)         RNOX1=(DY0*DZ1-DY1*DZ0)         RNOY1=(DZ0*DX1-DZ1*DX0)         RNOZ1=(DX0*DY1-DX1*DY0)         VOL=((RNOX1**2)+(RNOY1**2)+(RNOZ1**2))**0.5D0         VOL=VOL*0.5D0         PSCA=(RNOX0*RNOX1)+(RNOY0*RNOY1)+(RNOZ0*RNOZ1)         IF((VOL .LE. (RNOT*1.0D-6)) .OR. (PSCA.LE.0))THEN            LOGREG=.FALSE.         ENDIF         VOLT1=VOLT1+VOL         DO I2=1,IDIM,1            CEL=(XCOOR((IPCEN-1)*(IDIM+1)+I2)+     &           XCOOR((IPL-1)*(IDIM+1)+I2)+XCOOR((IPR-1)*(IDIM+1)+I2))     &           /3.0D0            XCEN1(I2)=XCEN1(I2)+(VOL*CEL)         ENDDO      ENDDOCC**** If one of the area is negative or almost zero the division is notC     possible (LOGREG=.FALSE.)C      IF(LOGREG)THEN         VOLT=VOLT1         DO I2=1,IDIM,1            XCEN(I2)=XCEN1(I2)         ENDDO      ENDIFCC**** Second division in trianglesC     We consider the triangles IP3, IP2, IP4 andC     IP4, IP2, IP1CC      VOLT1=0.0D0      LOGREG=.TRUE.      XCEN1(1)=0.0D0      XCEN1(2)=0.0D0      XCEN1(3)=0.0D0C      DO I1=1,2         IF(I1 .EQ. 1)THENC           We consider the triangles            IPCEN=IP2            IPR=IP3            IPL=IP4         ELSEIF(I1 .EQ. 2)THEN            IPCEN=IP2            IPR=IP4            IPL=IP1         ENDIF         DX0=XCOOR((IPR-1)*(IDIM+1)+1)-     &        XCOOR((IPCEN-1)*(IDIM+1)+1)         DY0=XCOOR((IPR-1)*(IDIM+1)+2)-     &        XCOOR((IPCEN-1)*(IDIM+1)+2)         IF(IDIM. EQ. 3) DZ0=XCOOR((IPR-1)*(IDIM+1)+3)-     &        XCOOR((IPCEN-1)*(IDIM+1)+3)         DX1=XCOOR((IPL-1)*(IDIM+1)+1)-     &        XCOOR((IPCEN-1)*(IDIM+1)+1)         DY1=XCOOR((IPL-1)*(IDIM+1)+2)-     &        XCOOR((IPCEN-1)*(IDIM+1)+2)         IF(IDIM. EQ. 3) DZ1=XCOOR((IPL-1)*(IDIM+1)+3)-     &        XCOOR((IPCEN-1)*(IDIM+1)+3)         RNOX1=(DY0*DZ1-DY1*DZ0)         RNOY1=(DZ0*DX1-DZ1*DX0)         RNOZ1=(DX0*DY1-DX1*DY0)         VOL=((RNOX1**2)+(RNOY1**2)+(RNOZ1**2))**0.5D0         VOL=VOL*0.5D0         PSCA=(RNOX0*RNOX1)+(RNOY0*RNOY1)+(RNOZ0*RNOZ1)         IF((VOL .LE. (RNOT*1.0D-6)) .OR. (PSCA.LE.0))THEN            LOGREG=.FALSE.         ENDIF         VOLT1=VOLT1+VOL         DO I2=1,IDIM,1            CEL=(XCOOR((IPCEN-1)*(IDIM+1)+I2)+     &           XCOOR((IPL-1)*(IDIM+1)+I2)+XCOOR((IPR-1)*(IDIM+1)+I2))     &           /3.0D0            XCEN1(I2)=XCEN1(I2)+(VOL*CEL)         ENDDO      ENDDOC      IF(LOGREG)THEN         VOLT=VOLT+VOLT1         DO I2=1,IDIM,1            XCEN(I2)=XCEN(I2)+XCEN1(I2)         ENDDO      ENDIFC      DO I2=1,IDIM,1         XCEN(I2)=XCEN(I2)/VOLT      ENDDOCC**** Finally, we check whether the center of gravity is "inside"C     the domain, i.e. whether the normals CENIPi x CENIP(i+1)C     have the same direction of (RNOX0,RNOY0,RNOZ0)C      LOGREG=.FALSE.      DO I1=1,4         IF(I1 .EQ. 1)THENC           We consider the triangles CEN, IP1, IP2            IPR=IP1            IPL=IP2         ELSEIF(I1 .EQ. 2)THEN            IPR=IP2            IPL=IP3         ELSEIF(I1 .EQ. 3)THEN            IPR=IP3            IPL=IP4         ELSEIF(I1 .EQ. 4)THEN            IPR=IP4            IPL=IP1         ENDIF         DX0=XCOOR((IPR-1)*(IDIM+1)+1)-     &        XCEN(1)         DY0=XCOOR((IPR-1)*(IDIM+1)+2)-     &        XCEN(2)         IF(IDIM. EQ. 3) DZ0=XCOOR((IPR-1)*(IDIM+1)+3)-     &        XCEN(3)         DX1=XCOOR((IPL-1)*(IDIM+1)+1)-     &        XCEN(1)         DY1=XCOOR((IPL-1)*(IDIM+1)+2)-     &        XCEN(2)         IF(IDIM. EQ. 3) DZ1=XCOOR((IPL-1)*(IDIM+1)+3)-     &        XCEN(3)         RNOX1=(DY0*DZ1-DY1*DZ0)         RNOY1=(DZ0*DX1-DZ1*DX0)         RNOZ1=(DX0*DY1-DX1*DY0)         PSCA=(RNOX0*RNOX1)+(RNOY0*RNOY1)+(RNOZ0*RNOZ1)         IF(PSCA .LE. 0.0D0)THEN            LOGREG=.TRUE.         ENDIF      ENDDOC      IF(LOGREG)THENC         WRITE(IOIMP,*) 'Subroutine kdom7.eso'         WRITE(IOIMP,*) 'QUA4'CC        Erreur -107 0C        Liste des noeuds de l'élément :C         CALL ERREUR(-107)         WRITE(IOIMP,*) IP1, IP2, IP3, IP4CC        Erreur 845 2C        Maillage donné incorrect ?!!!C         CALL ERREUR(845)C      ENDIFC      RETURN      END

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