Télécharger la plaquette Cast3M


fr en

Télécharger kdom7.eso

Retour à la liste

Numérotation des lignes : 

kdom7

  1. C KDOM7     SOURCE    KK2000    14/04/10    21:15:14     8032
  2.       SUBROUTINE KDOM7(IP1,IP2,IP3,IP4,XCEN)
  3. C
  4. C************************************************************************
  5. C
  6. C PROJET            :  CASTEM 2000
  7. C
  8. C NOM               :  KDOM7
  9. C
  10. C DESCRIPTION       :  Subroutine called by LEKMA0
  11. C                      Given a QUA4 (IP1,IP2,IP3,IP4), this subroutine
  12. C                      compute its center by averaging the centers of
  13. C                      the two possible configuration it can be
  14. C                      divided in
  15. C
  16. C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  17. C
  18. C AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/LTMF
  19. C
  20. C************************************************************************
  21. C
  22. C
  23.  
  24.       IMPLICIT INTEGER(I-N)
  25.       IMPLICIT REAL*8(A-H,O-Z)
  26.  
  27.       INTEGER IP1,IP2,IP3, IP4, I1,I2, IPCEN, IPL, IPR
  28.       REAL*8 VOL, XCEN(3),VOLT,DX0,DY0,DZ0,DX1,DY1,DZ1,CEL
  29.      &     ,RNOX0,RNOY0,RNOZ0,RNOX1,RNOY1,RNOZ1,XCEN1(3),VOLT1,PSCA
  30.      &     ,RNOT
  31.       LOGICAL LOGREG
  32.  
  33. -INC PPARAM
  34. -INC CCOPTIO
  35. -INC SMCOORD
  36. C
  37.       DZ0=0.0D0
  38.       DZ1=0.0D0
  39.       RNOT=0.0D0
  40. C
  41. C**** We chose as positive direction the one which
  42. C     has the maximum modulus IP1IP(1+1) x IP1IP(1+3)
  43. C
  44.       DO I1=1,4
  45.          IF(I1 .EQ. 1)THEN
  46. C           We consider the triangles IP1, IP2, IP4
  47.             IPCEN=IP1
  48.             IPR=IP2
  49.             IPL=IP4
  50.          ELSEIF(I1 .EQ. 2)THEN
  51.             IPCEN=IP2
  52.             IPR=IP3
  53.             IPL=IP1
  54.          ELSEIF(I1 .EQ. 3)THEN
  55.             IPCEN=IP3
  56.             IPR=IP4
  57.             IPL=IP2
  58.          ELSEIF(I1 .EQ. 4)THEN
  59.             IPCEN=IP4
  60.             IPR=IP1
  61.             IPL=IP3
  62.          ENDIF
  63.          DX0=XCOOR((IPR-1)*(IDIM+1)+1)-
  64.      &        XCOOR((IPCEN-1)*(IDIM+1)+1)
  65.          DY0=XCOOR((IPR-1)*(IDIM+1)+2)-
  66.      &        XCOOR((IPCEN-1)*(IDIM+1)+2)
  67.          IF(IDIM. EQ. 3) DZ0=XCOOR((IPR-1)*(IDIM+1)+3)-
  68.      &        XCOOR((IPCEN-1)*(IDIM+1)+3)
  69.          DX1=XCOOR((IPL-1)*(IDIM+1)+1)-
  70.      &        XCOOR((IPCEN-1)*(IDIM+1)+1)
  71.          DY1=XCOOR((IPL-1)*(IDIM+1)+2)-
  72.      &        XCOOR((IPCEN-1)*(IDIM+1)+2)
  73.          IF(IDIM. EQ. 3) DZ1=XCOOR((IPL-1)*(IDIM+1)+3)-
  74.      &        XCOOR((IPCEN-1)*(IDIM+1)+3)
  75.          RNOX1=(DY0*DZ1-DY1*DZ0)
  76.          RNOY1=(DZ0*DX1-DZ1*DX0)
  77.          RNOZ1=(DX0*DY1-DX1*DY0)
  78.          VOL=((RNOX1**2)+(RNOY1**2)+(RNOZ1**2))**0.5D0
  79.          VOL=VOL/0.5D0
  80.          IF(VOL .GT. RNOT)THEN
  81.             RNOT=VOL
  82.             RNOX0=RNOX1
  83.             RNOY0=RNOY1
  84.             RNOZ0=RNOZ1
  85.          ENDIF
  86.       ENDDO
  87. C
  88. C**** First division in triangles
  89. C     We consider the triangles IP3, IP1, IP2 and
  90. C     IP4, IP1, IP3
  91. C
  92.       XCEN1(1)=0.0D0
  93.       XCEN1(2)=0.0D0
  94.       XCEN1(3)=0.0D0
  95.       XCEN(1)=0.0D0
  96.       XCEN(2)=0.0D0
  97.       XCEN(3)=0.0D0
  98.       VOLT=0.0D0
  99.       VOLT1=0.0D0
  100.       LOGREG=.TRUE.
  101. C
  102.       DO I1=1,2
  103.          IF(I1 .EQ. 1)THEN
  104. C           We consider the triangles IP3, IP1, IP2
  105.             IPCEN=IP1
  106.             IPR=IP2
  107.             IPL=IP3
  108.          ELSEIF(I1 .EQ. 2)THEN
  109.             IPCEN=IP1
  110.             IPR=IP3
  111.             IPL=IP4
  112.          ENDIF
  113.          DX0=XCOOR((IPR-1)*(IDIM+1)+1)-
  114.      &        XCOOR((IPCEN-1)*(IDIM+1)+1)
  115.          DY0=XCOOR((IPR-1)*(IDIM+1)+2)-
  116.      &        XCOOR((IPCEN-1)*(IDIM+1)+2)
  117.          IF(IDIM. EQ. 3) DZ0=XCOOR((IPR-1)*(IDIM+1)+3)-
  118.      &        XCOOR((IPCEN-1)*(IDIM+1)+3)
  119.          DX1=XCOOR((IPL-1)*(IDIM+1)+1)-
  120.      &        XCOOR((IPCEN-1)*(IDIM+1)+1)
  121.          DY1=XCOOR((IPL-1)*(IDIM+1)+2)-
  122.      &        XCOOR((IPCEN-1)*(IDIM+1)+2)
  123.          IF(IDIM. EQ. 3) DZ1=XCOOR((IPL-1)*(IDIM+1)+3)-
  124.      &        XCOOR((IPCEN-1)*(IDIM+1)+3)
  125.          RNOX1=(DY0*DZ1-DY1*DZ0)
  126.          RNOY1=(DZ0*DX1-DZ1*DX0)
  127.          RNOZ1=(DX0*DY1-DX1*DY0)
  128.          VOL=((RNOX1**2)+(RNOY1**2)+(RNOZ1**2))**0.5D0
  129.          VOL=VOL*0.5D0
  130.          PSCA=(RNOX0*RNOX1)+(RNOY0*RNOY1)+(RNOZ0*RNOZ1)
  131.          IF((VOL .LE. (RNOT*1.0D-6)) .OR. (PSCA.LE.0))THEN
  132.             LOGREG=.FALSE.
  133.          ENDIF
  134.          VOLT1=VOLT1+VOL
  135.          DO I2=1,IDIM,1
  136.             CEL=(XCOOR((IPCEN-1)*(IDIM+1)+I2)+
  137.      &           XCOOR((IPL-1)*(IDIM+1)+I2)+XCOOR((IPR-1)*(IDIM+1)+I2))
  138.      &           /3.0D0
  139.             XCEN1(I2)=XCEN1(I2)+(VOL*CEL)
  140.          ENDDO
  141.       ENDDO
  142. C
  143. C**** If one of the area is negative or almost zero the division is not
  144. C     possible (LOGREG=.FALSE.)
  145. C
  146.       IF(LOGREG)THEN
  147.          VOLT=VOLT1
  148.          DO I2=1,IDIM,1
  149.             XCEN(I2)=XCEN1(I2)
  150.          ENDDO
  151.       ENDIF
  152. C
  153. C**** Second division in triangles
  154. C     We consider the triangles IP3, IP2, IP4 and
  155. C     IP4, IP2, IP1
  156. C
  157. C
  158.       VOLT1=0.0D0
  159.       LOGREG=.TRUE.
  160.       XCEN1(1)=0.0D0
  161.       XCEN1(2)=0.0D0
  162.       XCEN1(3)=0.0D0
  163. C
  164.       DO I1=1,2
  165.          IF(I1 .EQ. 1)THEN
  166. C           We consider the triangles
  167.             IPCEN=IP2
  168.             IPR=IP3
  169.             IPL=IP4
  170.          ELSEIF(I1 .EQ. 2)THEN
  171.             IPCEN=IP2
  172.             IPR=IP4
  173.             IPL=IP1
  174.          ENDIF
  175.          DX0=XCOOR((IPR-1)*(IDIM+1)+1)-
  176.      &        XCOOR((IPCEN-1)*(IDIM+1)+1)
  177.          DY0=XCOOR((IPR-1)*(IDIM+1)+2)-
  178.      &        XCOOR((IPCEN-1)*(IDIM+1)+2)
  179.          IF(IDIM. EQ. 3) DZ0=XCOOR((IPR-1)*(IDIM+1)+3)-
  180.      &        XCOOR((IPCEN-1)*(IDIM+1)+3)
  181.          DX1=XCOOR((IPL-1)*(IDIM+1)+1)-
  182.      &        XCOOR((IPCEN-1)*(IDIM+1)+1)
  183.          DY1=XCOOR((IPL-1)*(IDIM+1)+2)-
  184.      &        XCOOR((IPCEN-1)*(IDIM+1)+2)
  185.          IF(IDIM. EQ. 3) DZ1=XCOOR((IPL-1)*(IDIM+1)+3)-
  186.      &        XCOOR((IPCEN-1)*(IDIM+1)+3)
  187.          RNOX1=(DY0*DZ1-DY1*DZ0)
  188.          RNOY1=(DZ0*DX1-DZ1*DX0)
  189.          RNOZ1=(DX0*DY1-DX1*DY0)
  190.          VOL=((RNOX1**2)+(RNOY1**2)+(RNOZ1**2))**0.5D0
  191.          VOL=VOL*0.5D0
  192.          PSCA=(RNOX0*RNOX1)+(RNOY0*RNOY1)+(RNOZ0*RNOZ1)
  193.          IF((VOL .LE. (RNOT*1.0D-6)) .OR. (PSCA.LE.0))THEN
  194.             LOGREG=.FALSE.
  195.          ENDIF
  196.          VOLT1=VOLT1+VOL
  197.          DO I2=1,IDIM,1
  198.             CEL=(XCOOR((IPCEN-1)*(IDIM+1)+I2)+
  199.      &           XCOOR((IPL-1)*(IDIM+1)+I2)+XCOOR((IPR-1)*(IDIM+1)+I2))
  200.      &           /3.0D0
  201.             XCEN1(I2)=XCEN1(I2)+(VOL*CEL)
  202.          ENDDO
  203.       ENDDO
  204. C
  205.       IF(LOGREG)THEN
  206.          VOLT=VOLT+VOLT1
  207.          DO I2=1,IDIM,1
  208.             XCEN(I2)=XCEN(I2)+XCEN1(I2)
  209.          ENDDO
  210.       ENDIF
  211. C
  212.       DO I2=1,IDIM,1
  213.          XCEN(I2)=XCEN(I2)/VOLT
  214.       ENDDO
  215. C
  216. C**** Finally, we check whether the center of gravity is "inside"
  217. C     the domain, i.e. whether the normals CENIPi x CENIP(i+1)
  218. C     have the same direction of (RNOX0,RNOY0,RNOZ0)
  219. C
  220.       LOGREG=.FALSE.
  221.       DO I1=1,4
  222.          IF(I1 .EQ. 1)THEN
  223. C           We consider the triangles CEN, IP1, IP2
  224.             IPR=IP1
  225.             IPL=IP2
  226.          ELSEIF(I1 .EQ. 2)THEN
  227.             IPR=IP2
  228.             IPL=IP3
  229.          ELSEIF(I1 .EQ. 3)THEN
  230.             IPR=IP3
  231.             IPL=IP4
  232.          ELSEIF(I1 .EQ. 4)THEN
  233.             IPR=IP4
  234.             IPL=IP1
  235.          ENDIF
  236.          DX0=XCOOR((IPR-1)*(IDIM+1)+1)-
  237.      &        XCEN(1)
  238.          DY0=XCOOR((IPR-1)*(IDIM+1)+2)-
  239.      &        XCEN(2)
  240.          IF(IDIM. EQ. 3) DZ0=XCOOR((IPR-1)*(IDIM+1)+3)-
  241.      &        XCEN(3)
  242.          DX1=XCOOR((IPL-1)*(IDIM+1)+1)-
  243.      &        XCEN(1)
  244.          DY1=XCOOR((IPL-1)*(IDIM+1)+2)-
  245.      &        XCEN(2)
  246.          IF(IDIM. EQ. 3) DZ1=XCOOR((IPL-1)*(IDIM+1)+3)-
  247.      &        XCEN(3)
  248.          RNOX1=(DY0*DZ1-DY1*DZ0)
  249.          RNOY1=(DZ0*DX1-DZ1*DX0)
  250.          RNOZ1=(DX0*DY1-DX1*DY0)
  251.          PSCA=(RNOX0*RNOX1)+(RNOY0*RNOY1)+(RNOZ0*RNOZ1)
  252.          IF(PSCA .LE. 0.0D0)THEN
  253.             LOGREG=.TRUE.
  254.          ENDIF
  255.       ENDDO
  256. C
  257.       IF(LOGREG)THEN
  258. C
  259.          WRITE(IOIMP,*) 'Subroutine kdom7.eso'
  260.          WRITE(IOIMP,*) 'QUA4'
  261. C
  262. C        Erreur -107 0
  263. C        Liste des noeuds de l'élément :
  264. C
  265.          CALL ERREUR(-107)
  266.          WRITE(IOIMP,*) IP1, IP2, IP3, IP4
  267. C
  268. C        Erreur 845 2
  269. C        Maillage donné incorrect ?!!!
  270. C
  271.          CALL ERREUR(845)
  272. C
  273.       ENDIF
  274. C
  275.       RETURN
  276.       END
  277.  
  278.  
  279.  
  280.  

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