C KDOM7 SOURCE KK2000 14/04/10 21:15:14 8032 SUBROUTINE KDOM7(IP1,IP2,IP3,IP4,XCEN) C C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : KDOM7 C C DESCRIPTION : Subroutine called by LEKMA0 C Given a QUA4 (IP1,IP2,IP3,IP4), this subroutine C compute its center by averaging the centers of C the two possible configuration it can be C divided in C C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) C C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF C C************************************************************************ C C 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 PPARAM -INC CCOPTIO -INC SMCOORD C DZ0=0.0D0 DZ1=0.0D0 RNOT=0.0D0 C C**** We chose as positive direction the one which C has the maximum modulus IP1IP(1+1) x IP1IP(1+3) C DO I1=1,4 IF(I1 .EQ. 1)THEN C 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 ENDDO C C**** First division in triangles C We consider the triangles IP3, IP1, IP2 and C IP4, IP1, IP3 C 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)THEN C 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 ENDDO C C**** If one of the area is negative or almost zero the division is not C possible (LOGREG=.FALSE.) C IF(LOGREG)THEN VOLT=VOLT1 DO I2=1,IDIM,1 XCEN(I2)=XCEN1(I2) ENDDO ENDIF C C**** Second division in triangles C We consider the triangles IP3, IP2, IP4 and C IP4, IP2, IP1 C C VOLT1=0.0D0 LOGREG=.TRUE. XCEN1(1)=0.0D0 XCEN1(2)=0.0D0 XCEN1(3)=0.0D0 C DO I1=1,2 IF(I1 .EQ. 1)THEN C 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 ENDDO C IF(LOGREG)THEN VOLT=VOLT+VOLT1 DO I2=1,IDIM,1 XCEN(I2)=XCEN(I2)+XCEN1(I2) ENDDO ENDIF C DO I2=1,IDIM,1 XCEN(I2)=XCEN(I2)/VOLT ENDDO C C**** 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)THEN C 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 ENDDO C IF(LOGREG)THEN C WRITE(IOIMP,*) 'Subroutine kdom7.eso' WRITE(IOIMP,*) 'QUA4' C C Erreur -107 0 C Liste des noeuds de l'élément : C CALL ERREUR(-107) WRITE(IOIMP,*) IP1, IP2, IP3, IP4 C C Erreur 845 2 C Maillage donné incorrect ?!!! C CALL ERREUR(845) C ENDIF C RETURN END