C KDOM2 SOURCE KK2000 14/04/10 21:15:11 8032 SUBROUTINE KDOM2(MELEMQ) C C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : KDOM2 C C DESCRIPTION : Subroutine called by KDOM1 C In the case of EULER model, we change the C position of the QUAF points 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 MELEMQ, NBSOUS, ISOUS, NBEL, IEL, NN1, NN2, NN3, I1 & ,NN4,NN5,NN6,NN7,NN8,NN9,NN10,NN11,NN12,NN13,NN14,NN15 & ,NN16,NN17,NN18,NN19,NN20,NN21,NN22,NN23,NN24,NN25,NN26 & ,NN27 REAL*8 XCEN(3),VOL,X1(3),X2(3),X3(3),X4(3),X5(3),X6(3) & ,VOL1,VOL2,VOL3,VOL4,VOL5,VOL6 LOGICAL LOSEG3, LOTRI7, LOQUA9, LOTE15, LOPY19, LOPR21, LOCU27 C -INC PPARAM -INC CCOPTIO -INC SMELEME -INC SMCOORD C C Elements allowed are SEG3, TRI7, QUA9, TE15, 'PY19', 'PR21', 'CU27' C ITYPEL 3 7 11 35 36 34 33 C They can be referred just once C LOSEG3=.FALSE. LOTRI7=.FALSE. LOQUA9=.FALSE. LOTE15=.FALSE. LOPY19=.FALSE. LOPR21=.FALSE. LOCU27=.FALSE. C MELEME=MELEMQ SEGACT MELEME NBSOUS=MELEME.LISOUS(/1) IF(NBSOUS .EQ. 0) NBSOUS=1 DO ISOUS=1,NBSOUS,1 IF(NBSOUS .NE. 1)THEN IPT1=MELEME.LISOUS(ISOUS) SEGACT IPT1 ELSE IPT1=MELEME ENDIF NBEL=IPT1.NUM(/2) C IF(IPT1.ITYPEL .EQ. 3)THEN IF(LOSEG3)THEN C Elt already referred WRITE(IOIMP,*) 'Subroutine kdom2' WRITE(IOIMP,*) 'Mesh type not recognized' CALL ERREUR(5) ENDIF LOSEG3=.TRUE. C C********** SEG3 C DO IEL=1,NBEL,1 NN1=IPT1.NUM(1,IEL) NN2=IPT1.NUM(2,IEL) NN3=IPT1.NUM(3,IEL) DO I1=1,IDIM XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN1-1)*(IDIM+1)+I1)) ENDDO ENDDO ELSEIF(IPT1.ITYPEL .EQ. 7)THEN C C********** TRI7 C IF(LOTRI7)THEN C Elt already referred WRITE(IOIMP,*) 'Subroutine kdom2' WRITE(IOIMP,*) 'Mesh type not recognized' CALL ERREUR(5) ENDIF LOTRI7=.TRUE. C DO IEL=1,NBEL,1 NN1=IPT1.NUM(1,IEL) NN2=IPT1.NUM(2,IEL) NN3=IPT1.NUM(3,IEL) NN4=IPT1.NUM(4,IEL) NN5=IPT1.NUM(5,IEL) NN6=IPT1.NUM(6,IEL) NN7=IPT1.NUM(7,IEL) DO I1=1,IDIM XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN3-1)*(IDIM+1)+I1)) XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1)) XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1)) XCOOR((NN7-1)*(IDIM+1)+I1)= & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1))/3.0D0 ENDDO ENDDO ELSEIF(IPT1.ITYPEL .EQ. 11)THEN C C********** QUA9 C IF(LOQUA9)THEN C Elt already referred WRITE(IOIMP,*) 'Subroutine kdom2' WRITE(IOIMP,*) 'Mesh type not recognized' CALL ERREUR(5) ENDIF LOQUA9=.TRUE. C DO IEL=1,NBEL,1 NN1=IPT1.NUM(1,IEL) NN2=IPT1.NUM(2,IEL) NN3=IPT1.NUM(3,IEL) NN4=IPT1.NUM(4,IEL) NN5=IPT1.NUM(5,IEL) NN6=IPT1.NUM(6,IEL) NN7=IPT1.NUM(7,IEL) NN8=IPT1.NUM(8,IEL) NN9=IPT1.NUM(9,IEL) C C**************** Centre de l'elt C CALL KDOM7(NN1,NN3,NN5,NN7,XCEN) IF(IERR.NE.0)GOTO 9999 DO I1=1,IDIM XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN3-1)*(IDIM+1)+I1)) XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1)) XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN7-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1)) XCOOR((NN8-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN7-1)*(IDIM+1)+I1)+ & XCOOR((NN1-1)*(IDIM+1)+I1)) XCOOR((NN9-1)*(IDIM+1)+I1)=XCEN(I1) ENDDO ENDDO ELSEIF(IPT1.ITYPEL .EQ. 35)THEN C C********** TE15 C IF(LOTE15)THEN C Elt already referred WRITE(IOIMP,*) 'Subroutine kdom2' WRITE(IOIMP,*) 'Mesh type not recognized' CALL ERREUR(5) ENDIF LOTE15=.TRUE. C DO IEL=1,NBEL,1 NN1=IPT1.NUM(1,IEL) NN2=IPT1.NUM(2,IEL) NN3=IPT1.NUM(3,IEL) NN4=IPT1.NUM(4,IEL) NN5=IPT1.NUM(5,IEL) NN6=IPT1.NUM(6,IEL) NN7=IPT1.NUM(7,IEL) NN8=IPT1.NUM(8,IEL) NN9=IPT1.NUM(9,IEL) NN10=IPT1.NUM(10,IEL) NN11=IPT1.NUM(11,IEL) NN12=IPT1.NUM(12,IEL) NN13=IPT1.NUM(13,IEL) NN14=IPT1.NUM(14,IEL) NN15=IPT1.NUM(15,IEL) DO I1=1,IDIM XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN3-1)*(IDIM+1)+I1)) XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1)) XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN5-1)*(IDIM+1)+I1)+ & XCOOR((NN1-1)*(IDIM+1)+I1)) XCOOR((NN7-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN10-1)*(IDIM+1)+I1)+ & XCOOR((NN1-1)*(IDIM+1)+I1)) XCOOR((NN8-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN10-1)*(IDIM+1)+I1)+ & XCOOR((NN3-1)*(IDIM+1)+I1)) XCOOR((NN9-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN10-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1)) XCOOR((NN11-1)*(IDIM+1)+I1)= & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1))/3.0D0 XCOOR((NN12-1)*(IDIM+1)+I1)= & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN10-1)*(IDIM+1)+I1))/3.0D0 XCOOR((NN13-1)*(IDIM+1)+I1)= & (XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1)+ & XCOOR((NN10-1)*(IDIM+1)+I1))/3.0D0 XCOOR((NN14-1)*(IDIM+1)+I1)= & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1)+ & XCOOR((NN10-1)*(IDIM+1)+I1))/3.0D0 XCOOR((NN15-1)*(IDIM+1)+I1)= & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1)+ & XCOOR((NN10-1)*(IDIM+1)+I1))/4.0D0 ENDDO ENDDO ELSEIF(IPT1.ITYPEL .EQ. 36)THEN C C********** PY19 C IF(LOPY19)THEN C Elt already referred WRITE(IOIMP,*) 'Subroutine kdom2' WRITE(IOIMP,*) 'Mesh type not recognized' CALL ERREUR(5) ENDIF LOPY19=.TRUE. C DO IEL=1,NBEL,1 NN1=IPT1.NUM(1,IEL) NN2=IPT1.NUM(2,IEL) NN3=IPT1.NUM(3,IEL) NN4=IPT1.NUM(4,IEL) NN5=IPT1.NUM(5,IEL) NN6=IPT1.NUM(6,IEL) NN7=IPT1.NUM(7,IEL) NN8=IPT1.NUM(8,IEL) NN9=IPT1.NUM(9,IEL) NN10=IPT1.NUM(10,IEL) NN11=IPT1.NUM(11,IEL) NN12=IPT1.NUM(12,IEL) NN13=IPT1.NUM(13,IEL) NN14=IPT1.NUM(14,IEL) NN15=IPT1.NUM(15,IEL) NN16=IPT1.NUM(16,IEL) NN17=IPT1.NUM(17,IEL) NN18=IPT1.NUM(18,IEL) NN19=IPT1.NUM(19,IEL) C Computation of the center position CALL KDOM7(NN1,NN3,NN5,NN7,XCEN) IF(IERR.NE.0)GOTO 9999 DO I1=1,IDIM XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN3-1)*(IDIM+1)+I1)) XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1)) XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN5-1)*(IDIM+1)+I1)+ & XCOOR((NN7-1)*(IDIM+1)+I1)) XCOOR((NN8-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN7-1)*(IDIM+1)+I1)) XCOOR((NN9-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN13-1)*(IDIM+1)+I1)+ & XCOOR((NN1-1)*(IDIM+1)+I1)) XCOOR((NN10-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN13-1)*(IDIM+1)+I1)+ & XCOOR((NN3-1)*(IDIM+1)+I1)) XCOOR((NN11-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN13-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1)) XCOOR((NN12-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN13-1)*(IDIM+1)+I1)+ & XCOOR((NN7-1)*(IDIM+1)+I1)) XCOOR((NN14-1)*(IDIM+1)+I1)=XCEN(I1) XCOOR((NN15-1)*(IDIM+1)+I1)= & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN13-1)*(IDIM+1)+I1))/3.0D0 XCOOR((NN16-1)*(IDIM+1)+I1)= & (XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1)+ & XCOOR((NN13-1)*(IDIM+1)+I1))/3.0D0 XCOOR((NN17-1)*(IDIM+1)+I1)= & (XCOOR((NN5-1)*(IDIM+1)+I1)+ & XCOOR((NN7-1)*(IDIM+1)+I1)+ & XCOOR((NN13-1)*(IDIM+1)+I1))/3.0D0 XCOOR((NN18-1)*(IDIM+1)+I1)= & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN7-1)*(IDIM+1)+I1)+ & XCOOR((NN13-1)*(IDIM+1)+I1))/3.0D0 ENDDO C C************* To compute the center, we divide the pyramid into C 4 tetrahedra (see kdom3) C CALL KDOM3(NN1,NN3,NN5,NN7,NN14,NN13,XCEN,VOL) DO I1=1,3,1 XCOOR((NN19-1)*(IDIM+1)+I1)=XCEN(I1) ENDDO ENDDO ELSEIF(IPT1.ITYPEL .EQ. 34)THEN C C********** PR21 C IF(LOPR21)THEN C Elt already referred WRITE(IOIMP,*) 'Subroutine kdom2' WRITE(IOIMP,*) 'Mesh type not recognized' CALL ERREUR(5) ENDIF LOPR21=.TRUE. C DO IEL=1,NBEL,1 NN1=IPT1.NUM(1,IEL) NN2=IPT1.NUM(2,IEL) NN3=IPT1.NUM(3,IEL) NN4=IPT1.NUM(4,IEL) NN5=IPT1.NUM(5,IEL) NN6=IPT1.NUM(6,IEL) NN7=IPT1.NUM(7,IEL) NN8=IPT1.NUM(8,IEL) NN9=IPT1.NUM(9,IEL) NN10=IPT1.NUM(10,IEL) NN11=IPT1.NUM(11,IEL) NN12=IPT1.NUM(12,IEL) NN13=IPT1.NUM(13,IEL) NN14=IPT1.NUM(14,IEL) NN15=IPT1.NUM(15,IEL) NN16=IPT1.NUM(16,IEL) NN17=IPT1.NUM(17,IEL) NN18=IPT1.NUM(18,IEL) NN19=IPT1.NUM(19,IEL) NN20=IPT1.NUM(20,IEL) NN21=IPT1.NUM(21,IEL) DO I1=1,IDIM XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN3-1)*(IDIM+1)+I1)) XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1)) XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN5-1)*(IDIM+1)+I1)+ & XCOOR((NN1-1)*(IDIM+1)+I1)) XCOOR((NN7-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN10-1)*(IDIM+1)+I1)) XCOOR((NN8-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN12-1)*(IDIM+1)+I1)) XCOOR((NN9-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN5-1)*(IDIM+1)+I1)+ & XCOOR((NN14-1)*(IDIM+1)+I1)) XCOOR((NN11-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN10-1)*(IDIM+1)+I1)+ & XCOOR((NN12-1)*(IDIM+1)+I1)) XCOOR((NN13-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN12-1)*(IDIM+1)+I1)+ & XCOOR((NN14-1)*(IDIM+1)+I1)) XCOOR((NN15-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN10-1)*(IDIM+1)+I1)+ & XCOOR((NN14-1)*(IDIM+1)+I1)) XCOOR((NN19-1)*(IDIM+1)+I1)= & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1))/3.0D0 XCOOR((NN20-1)*(IDIM+1)+I1)= & (XCOOR((NN10-1)*(IDIM+1)+I1)+ & XCOOR((NN12-1)*(IDIM+1)+I1)+ & XCOOR((NN14-1)*(IDIM+1)+I1))/3.0D0 ENDDO C Computation of the QUA center positions CALL KDOM7(NN1,NN10,NN12,NN3,XCEN) IF(IERR.NE.0)GOTO 9999 DO I1=1,IDIM,1 XCOOR((NN16-1)*(IDIM+1)+I1)=XCEN(I1) ENDDO CALL KDOM7(NN5,NN3,NN12,NN14,XCEN) IF(IERR.NE.0)GOTO 9999 DO I1=1,IDIM,1 XCOOR((NN17-1)*(IDIM+1)+I1)=XCEN(I1) ENDDO CALL KDOM7(NN10,NN1,NN5,NN14,XCEN) IF(IERR.NE.0)GOTO 9999 DO I1=1,IDIM,1 XCOOR((NN18-1)*(IDIM+1)+I1)=XCEN(I1) ENDDO C C************* To compute the center, we divide the PR21 C into 2 tetrahedra and 3 pyramids C (1,3,5,PC), (10,14,12,PC), (1,10,12,3,16,PC), C (3,12,14,5,17,PC),(1,5,14,10,18,PC) C | C center of 3,12,14,5 C The results does not depend of the choice of C PC. C C C DO I1=1,3,1 XCOOR((NN21-1)*(IDIM+1)+I1)= & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1)+ & XCOOR((NN10-1)*(IDIM+1)+I1)+ & XCOOR((NN12-1)*(IDIM+1)+I1)+ & XCOOR((NN14-1)*(IDIM+1)+I1))/6.0D0 ENDDO C CALL KDOM3(NN10,NN12,NN3,NN1,NN16,NN21,X1,VOL1) CALL KDOM3(NN3,NN12,NN14,NN5,NN17,NN21,X2,VOL2) CALL KDOM3(NN1,NN5,NN14,NN10,NN18,NN21,X3,VOL3) CALL KDOM4(NN1,NN3,NN5,NN21,X4,VOL4) CALL KDOM4(NN10,NN14,NN12,NN21,X5,VOL5) VOL=VOL1+VOL2+VOL3+VOL4+VOL5 C DO I1=1,3,1 XCOOR((NN21-1)*(IDIM+1)+I1)=(VOL1*X1(I1)+VOL2*X2(I1) & +VOL3*X3(I1)+VOL4*X4(I1)+VOL5*X5(I1))/VOL ENDDO C C********** We check whether the center of gravity of the PR21 is C inside its domain C CALL KDOM8(NN1,NN2,NN3,NN4,NN5,NN6,NN7,NN8,NN9,NN10, & NN11,NN12,NN13,NN14,NN15,NN16,NN17,NN18,NN19,NN20, & NN21) IF(IERR.NE.0)GOTO 9999 ENDDO IF(NBSOUS .NE. 1) SEGDES IPT1 ELSEIF(IPT1.ITYPEL .EQ. 33)THEN C C********** CU27 C IF(LOCU27)THEN C Elt already referred WRITE(IOIMP,*) 'Subroutine kdom2' WRITE(IOIMP,*) 'Mesh type not recognized' CALL ERREUR(5) ENDIF LOCU27=.TRUE. C DO IEL=1,NBEL,1 NN1=IPT1.NUM(1,IEL) NN2=IPT1.NUM(2,IEL) NN3=IPT1.NUM(3,IEL) NN4=IPT1.NUM(4,IEL) NN5=IPT1.NUM(5,IEL) NN6=IPT1.NUM(6,IEL) NN7=IPT1.NUM(7,IEL) NN8=IPT1.NUM(8,IEL) NN9=IPT1.NUM(9,IEL) NN10=IPT1.NUM(10,IEL) NN11=IPT1.NUM(11,IEL) NN12=IPT1.NUM(12,IEL) NN13=IPT1.NUM(13,IEL) NN14=IPT1.NUM(14,IEL) NN15=IPT1.NUM(15,IEL) NN16=IPT1.NUM(16,IEL) NN17=IPT1.NUM(17,IEL) NN18=IPT1.NUM(18,IEL) NN19=IPT1.NUM(19,IEL) NN20=IPT1.NUM(20,IEL) NN21=IPT1.NUM(21,IEL) NN22=IPT1.NUM(22,IEL) NN23=IPT1.NUM(23,IEL) NN24=IPT1.NUM(24,IEL) NN25=IPT1.NUM(25,IEL) NN26=IPT1.NUM(26,IEL) NN27=IPT1.NUM(27,IEL) DO I1=1,IDIM XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN3-1)*(IDIM+1)+I1)) XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1)) XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN5-1)*(IDIM+1)+I1)+ & XCOOR((NN7-1)*(IDIM+1)+I1)) XCOOR((NN8-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN7-1)*(IDIM+1)+I1)) XCOOR((NN9-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN13-1)*(IDIM+1)+I1)) XCOOR((NN10-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN15-1)*(IDIM+1)+I1)) XCOOR((NN11-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN5-1)*(IDIM+1)+I1)+ & XCOOR((NN17-1)*(IDIM+1)+I1)) XCOOR((NN12-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN7-1)*(IDIM+1)+I1)+ & XCOOR((NN19-1)*(IDIM+1)+I1)) XCOOR((NN14-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN13-1)*(IDIM+1)+I1)+ & XCOOR((NN15-1)*(IDIM+1)+I1)) XCOOR((NN16-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN17-1)*(IDIM+1)+I1)+ & XCOOR((NN15-1)*(IDIM+1)+I1)) XCOOR((NN18-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN17-1)*(IDIM+1)+I1)+ & XCOOR((NN19-1)*(IDIM+1)+I1)) XCOOR((NN20-1)*(IDIM+1)+I1)=0.5D0* & (XCOOR((NN19-1)*(IDIM+1)+I1)+ & XCOOR((NN13-1)*(IDIM+1)+I1)) ENDDO C Computation of the QUA center positions CALL KDOM7(NN1,NN13,NN15,NN3,XCEN) IF(IERR.NE.0) GOTO 9999 DO I1=1,IDIM,1 XCOOR((NN21-1)*(IDIM+1)+I1)=XCEN(I1) ENDDO CALL KDOM7(NN3,NN15,NN17,NN5,XCEN) IF(IERR.NE.0)GOTO 9999 DO I1=1,IDIM,1 XCOOR((NN22-1)*(IDIM+1)+I1)=XCEN(I1) ENDDO CALL KDOM7(NN5,NN17,NN19,NN7,XCEN) IF(IERR.NE.0)GOTO 9999 DO I1=1,IDIM,1 XCOOR((NN23-1)*(IDIM+1)+I1)=XCEN(I1) ENDDO CALL KDOM7(NN7,NN19,NN13,NN1,XCEN) IF(IERR.NE.0)GOTO 9999 DO I1=1,IDIM,1 XCOOR((NN24-1)*(IDIM+1)+I1)=XCEN(I1) ENDDO CALL KDOM7(NN1,NN3,NN5,NN7,XCEN) IF(IERR.NE.0)GOTO 9999 DO I1=1,IDIM,1 XCOOR((NN25-1)*(IDIM+1)+I1)=XCEN(I1) ENDDO CALL KDOM7(NN13,NN19,NN17,NN15,XCEN) IF(IERR.NE.0)GOTO 9999 DO I1=1,IDIM,1 XCOOR((NN26-1)*(IDIM+1)+I1)=XCEN(I1) ENDDO C C************* To compute the center, we divide the cube C into 6 pyramids. C (1,3,5,7,25,PC) C (13,19,17,15,26,PC) C (1,13,15,3,21,PC) C (3,15,17,5,22,PC) C (5,17,19,7,23,PC) C (1,7,19,13,24,PC) C The results does not depend of the choice of C PC. C DO I1=1,3,1 XCOOR((NN27-1)*(IDIM+1)+I1)= & (XCOOR((NN1-1)*(IDIM+1)+I1)+ & XCOOR((NN3-1)*(IDIM+1)+I1)+ & XCOOR((NN5-1)*(IDIM+1)+I1)+ & XCOOR((NN7-1)*(IDIM+1)+I1)+ & XCOOR((NN13-1)*(IDIM+1)+I1)+ & XCOOR((NN15-1)*(IDIM+1)+I1)+ & XCOOR((NN17-1)*(IDIM+1)+I1)+ & XCOOR((NN19-1)*(IDIM+1)+I1))/8.0D0 ENDDO C CALL KDOM3(NN1,NN3,NN5,NN7,NN25,NN27,X1,VOL1) CALL KDOM3(NN13,NN19,NN17,NN15,NN26,NN27,X2,VOL2) CALL KDOM3(NN1,NN13,NN15,NN3,NN21,NN27,X3,VOL3) CALL KDOM3(NN3,NN15,NN17,NN5,NN22,NN27,X4,VOL4) CALL KDOM3(NN5,NN17,NN19,NN7,NN23,NN27,X5,VOL5) CALL KDOM3(NN1,NN7,NN19,NN13,NN24,NN27,X6,VOL6) VOL=VOL1+VOL2+VOL3+VOL4+VOL5+VOL6 DO I1=1,3,1 XCOOR((NN27-1)*(IDIM+1)+I1)=(VOL1*X1(I1)+VOL2*X2(I1) & +VOL3*X3(I1)+VOL4*X4(I1)+VOL5*X5(I1)+VOL6*X6(I1)) & /VOL ENDDO C C********** We check whether the center of gravity of the CU27 is C inside its domain C CALL KDOM9(NN1,NN2,NN3,NN4,NN5,NN6,NN7,NN8,NN9,NN10, & NN11,NN12,NN13,NN14,NN15,NN16,NN17,NN18,NN19,NN20, & NN21,NN22,NN23,NN24,NN25,NN26,NN27) IF(IERR.NE.0)GOTO 9999 ENDDO ELSE WRITE(IOIMP,*) 'Euler model' CALL ERREUR(21) ENDIF IF(NBSOUS .NE. 1) SEGDES IPT1 ENDDO C SEGDES MELEME C 9999 RETURN C END