Télécharger la plaquette Cast3M


fr en

Télécharger inclus.eso

Retour à la liste

Numérotation des lignes : 

inclus

  1. C INCLUS    SOURCE    PV        20/03/30    21:20:00     10567          
  2. C
  3. C   CE SOUS PROGRAMME EXTRAIT DU PREMIER OBJET L'ENSEMBLE DES ELEMENTS
  4. C   STRICTEMENT INCLUS DANS LE DEUXIEME OBJET
  5. C
  6.       SUBROUTINE INCLUS
  7.       IMPLICIT INTEGER(I-N)
  8.       IMPLICIT REAL*8 (A-H,O-Z)
  9. -INC CCREEL
  10. *-
  11. -INC SMELEME
  12.  
  13. -INC PPARAM
  14. -INC CCOPTIO
  15. -INC SMCOORD
  16. -INC CCGEOME
  17.       CHARACTER*4 LEMOT(5)
  18.       DATA LEMOT/'BARY','STRI','LARG','VOLU','NOID'/
  19.       SEGMENT /FER/(NFI(ITT),MAI(IPP),ITOUR)
  20.       SEGMENT ICPR(nbpts)
  21.  
  22.       CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU)
  23.       CALL ACTOBJ('MAILLAGE',IPT1,1)
  24.       CALL LIROBJ('MAILLAGE',IPT2,1,IRETOU)
  25.       CALL ACTOBJ('MAILLAGE',IPT2,1)
  26.       IF (IERR.NE.0) RETURN
  27.  
  28. *     On fixe les options
  29. *     Par défaut : non barycentrique, appartenance large, et erreur si
  30. *     maillage vide
  31.       IBARY=0
  32.       IFRON=0
  33.       IVERI=0
  34.       DO I=1,3
  35.          CALL LIRMOT(LEMOT,5,IRET,0)
  36.          IF (IRET.EQ.1) IBARY=1
  37.          IF (IRET.EQ.2) IFRON=1
  38.          IF (IRET.EQ.3) IFRON=0
  39.          IF (IRET.EQ.4) THEN
  40. *           Cas de maillage volumique
  41.             CALL INCLU3(IPT1,IPT2)
  42.             RETURN
  43.          ENDIF
  44.          IF (IRET.EQ.5) IVERI=1
  45.          IF (IRET.EQ.0) GOTO 20
  46.       ENDDO
  47.  
  48.  20   CONTINUE
  49.  
  50. C CRITERE D'INCLUSION PAS UTILISE EN 2D...
  51.       CALL LIRREE(XCRITT,0,IRET)
  52.       XCRITT=0.D0
  53.  
  54. *     Cas d'une surface plane.
  55.       SEGACT IPT2
  56.  
  57. C     Appel éventuel à CONTOUR
  58.       IF ((KSURF(IPT2.ITYPEL).EQ.0).AND.(IPT2.LISOUS(/1).EQ.0)) GOTO 25
  59.       CALL ECROBJ('MAILLAGE',IPT2)
  60.       CALL PRCONT
  61.       CALL LIROBJ('MAILLAGE',IPT2,1,IRETOU)
  62.       IF (IERR.NE.0) RETURN
  63.       SEGACT IPT2
  64.  
  65.  25   CONTINUE
  66. *     Transformation éventuelle en SEG2 (SEG3 ou maillage complexe)
  67.       IF(IPT2.ITYPEL.EQ.3.OR.IPT2.LISOUS(/1).NE.0)
  68.      $     CALL CHANGE (IPT2,2)
  69.  
  70. *     Remplissage segment FER
  71.       ICLE=0
  72.       CALL AVTRSF(IPT2,FER,ICLE)
  73.       IF (IERR.NE.0) RETURN
  74.  
  75. *     Détermination du plan par son centre de gravité et son vecteur
  76. *     normal unitaire (pour changement de variable)
  77.       SEGACT MCOORD
  78.       SEGINI ICPR
  79.       DO 1 I=1,ICPR(/1)
  80.          ICPR(I)=0
  81.  1    CONTINUE
  82.       XGRAV=0.
  83.       YGRAV=0.
  84.       ZGRAV=0.
  85.       DO 30 I=MAI(1)+1,MAI(ITOUR+1)
  86.          IREF=NFI(I)*(IDIM+1)-IDIM
  87.          XGRAV=XGRAV+XCOOR(IREF)
  88.          YGRAV=YGRAV+XCOOR(IREF+1)
  89.          ZGRAV=ZGRAV+XCOOR(IREF+2)
  90.  30   CONTINUE
  91.       XGRAV=XGRAV/FLOAT(MAI(ITOUR+1)-MAI(1))
  92.       YGRAV=YGRAV/FLOAT(MAI(ITOUR+1)-MAI(1))
  93.       ZGRAV=ZGRAV/FLOAT(MAI(ITOUR+1)-MAI(1))
  94.       IF (IDIM.EQ.2) ZGRAV=0.
  95.       XNORM=0.D0
  96.       YNORM=0.D0
  97.       ZNORM=0.D0
  98.       DO 42 IT=1,ITOUR
  99.          IREF=NFI(MAI(IT+1))*(IDIM+1)-IDIM
  100.          XV1=XCOOR(IREF)-XGRAV
  101.          YV1=XCOOR(IREF+1)-YGRAV
  102.          ZV1=XCOOR(IREF+2)-ZGRAV
  103.          IF (IDIM.EQ.2) ZV1=0.
  104.          DO 421 I=MAI(IT-1+1)+1,MAI(IT+1)
  105.             IREF=(NFI(I))*(IDIM+1)-IDIM
  106.             XV2=XCOOR(IREF)-XGRAV
  107.             YV2=XCOOR(IREF+1)-YGRAV
  108.             ZV2=0.
  109.             IF (IDIM.GE.3) ZV2=XCOOR(IREF+2)-ZGRAV
  110.             XNORM=XNORM+YV1*ZV2-ZV1*YV2
  111.             YNORM=YNORM+ZV1*XV2-XV1*ZV2
  112.             ZNORM=ZNORM+XV1*YV2-YV1*XV2
  113.             XV1=XV2
  114.             YV1=YV2
  115.             ZV1=ZV2
  116.  421     CONTINUE
  117.  42   CONTINUE
  118.       DNORM=SQRT(XNORM**2+YNORM**2+ZNORM**2)
  119.       XNORM=XNORM/DNORM
  120.       YNORM=YNORM/DNORM
  121.       ZNORM=ZNORM/DNORM
  122.  
  123.       MELEME=IPT1
  124.       SEGACT MELEME
  125.       NBS   =LISOUS(/1)
  126.       NBSOUS=NBS
  127.       IPT8=0
  128.       IF (NBSOUS.NE.0) THEN
  129.          NBREF=0
  130.          NBNN=0
  131.          NBELEM=0
  132.          SEGINI IPT8
  133.          IOB1=0
  134.       ENDIF
  135.  
  136.       IF(IBARY.EQ.0) THEN
  137. *        Eligibilité d'après chaque point de l'élément
  138.          DO 7001 IOB=1,MAX(1,NBS)
  139.             IF (NBS.NE.0) THEN
  140.                IPT1=LISOUS(IOB)
  141.                SEGACT IPT1
  142.             ENDIF
  143.             NBELEM=IPT1.NUM(/2)
  144.             NBNN=IPT1.NUM(/1)
  145.             NBREF=0
  146.             NBSOUS=0
  147.             SEGINI IPT4
  148.             IPT4.ITYPEL=IPT1.ITYPEL
  149.             DO 2 I=1,NBNN
  150.                DO 2222 J=1,NBELEM
  151.                   IPT4.NUM(I,J)=0
  152.  2222          CONTINUE
  153.  2          CONTINUE
  154.  
  155.             ICO=0
  156.             DO 10 J=1,NBELEM
  157.                DO 11 I=1,NBNN
  158.                   IP=IPT1.NUM(I,J)
  159.                   IF (ICPR(IP).EQ.-1) GOTO 10
  160.                   IF (ICPR(IP).EQ.1) GOTO 11
  161.                   IREF=IP*(IDIM+1)-IDIM
  162.                   XP=XCOOR(IREF)
  163.                   YP=XCOOR(IREF+1)
  164.                   ZP=XCOOR(IREF+2)
  165.                   IF (IDIM.EQ.2) ZP=0.
  166.                   ANG=0.D0
  167.                   DO 12 IT=1,ITOUR
  168.                      IF (IP.EQ.NFI(MAI(IT+1)))THEN
  169.                         IF(IFRON.EQ.0) GO TO 13
  170.                         GO TO 18
  171.                      ENDIF
  172.                      IREF=NFI(MAI(IT+1))*(IDIM+1)-IDIM
  173.                      XV1=XCOOR(IREF)-XP
  174.                      YV1=XCOOR(IREF+1)-YP
  175.                      ZV1=XCOOR(IREF+2)-ZP
  176.                      IF (IDIM.EQ.2) ZV1=0.
  177.                      DO 120 M=MAI(IT-1+1)+1,MAI(IT+1)
  178.                         IF (IP.EQ.NFI(M))THEN
  179.                            IF(IFRON.EQ.0) GO TO 13
  180.                            GO TO 18
  181.                         ENDIF
  182.                         IREF=NFI(M)*(IDIM+1)-IDIM
  183.                         XV2=XCOOR(IREF)-XP
  184.                         YV2=XCOOR(IREF+1)-YP
  185.                         ZV2=XCOOR(IREF+2)-ZP
  186.                         IF (IDIM.EQ.2) ZV2=0.
  187.                         V1=XNORM*(YV1*ZV2-ZV1*YV2)+YNORM*(ZV1*XV2-XV1
  188.      $                       *ZV2)+ZNORM*(XV1*YV2-YV1*XV2)
  189.                         V2=XV1*XV2+YV1*YV2+ZV1*ZV2
  190.                         IF( IFRON.EQ.1) THEN
  191.                            V6 = ( XV1 -XV2)*( XV1 -XV2) + ( YV1 -YV2)* (
  192.      $                          YV1 -YV2)+ ( ZV1-ZV2)* ( ZV1-ZV2)
  193.                            IF( V2.LE.0. . OR. ( V2 . LT . V6* 1.D-4))
  194.      $                          THEN
  195.                               V4 =SQRT (XV1* XV1 + YV1 * YV1 + ZV1 * ZV1
  196.      $                             )
  197.                               V5 =SQRT  (XV2* XV2 + YV2 * YV2 + ZV2
  198.      $                             * ZV2  )
  199.                               V3 = ((YV1*ZV2-ZV1*YV2)*(YV1*ZV2-ZV1*YV2)
  200.      $                             +(ZV1*XV2-XV1*ZV2)*(ZV1*XV2-XV1*ZV2)
  201.      $                             +(XV1*YV2-YV1*XV2) *(XV1*YV2-YV1*XV2)
  202.      $                             ) **0.5D0
  203.                               IF( V4* 1.D4. LT. V5 . OR . V5*1.D4 .LT.V4
  204.      $                             ) GO TO 18
  205.                               IF(  V3*1.D4   . LT . ( ABS ( V2) )  )
  206.      #                             go to 18
  207.                            ENDIF
  208.                         ENDIF
  209.                         IF (V1.EQ.0.D0.AND.V2.EQ.0.D0) GOTO 121
  210.                         ANG=ANG+ATAN2(V1,V2)
  211.  121                    CONTINUE
  212.                         XV1=XV2
  213.                         YV1=YV2
  214.                         ZV1=ZV2
  215.  120                 CONTINUE
  216.  12               CONTINUE
  217.                   IF (ABS(ANG).GT.XPI) GOTO 13
  218.  18               ICPR(IP)=-1
  219.                   GOTO 10
  220.  13               ICPR(IP)=1
  221.  11            CONTINUE
  222. *              Cet élément est ok, on le met de côté.
  223.                DO 14 L=1,NBNN
  224.                   IPT4.NUM(L,J)=IPT1.NUM(L,J)
  225.  14            CONTINUE
  226.                IPT4.ICOLOR(J)=IPT1.ICOLOR(J)
  227.                ICO=ICO+1
  228.  10         CONTINUE
  229.  
  230.             NBELEM=ICO
  231.             IF (NBELEM.EQ.0) GOTO 7001
  232.  
  233. *           sauvegarde partition ok
  234.             SEGINI IPT3
  235.             IPT3.ITYPEL=IPT1.ITYPEL
  236.             JJ=0
  237.             DO 21 J=1,IPT4.NUM(/2)
  238.                IF (IPT4.NUM(1,J).EQ.0) GOTO 21
  239.                JJ=JJ+1
  240.                IPT3.ICOLOR(JJ)=IPT4.ICOLOR(J)
  241.                DO 22 I=1,NBNN
  242.                   IPT3.NUM(I,JJ)=IPT4.NUM(I,J)
  243.  22            CONTINUE
  244.  21         CONTINUE
  245.             SEGSUP IPT4
  246.  
  247.             IF (NBS.NE.0) THEN
  248.                IOB1=IOB1+1
  249.                IPT8.LISOUS(IOB1)=IPT3
  250.             ENDIF
  251.  7001    CONTINUE
  252.       ELSE
  253. *        Eligibilité d'après la position du barycentre
  254.          DO 7051 IOB=1,MAX(1,NBS)
  255.             IF (NBS.NE.0) THEN
  256.                IPT1=LISOUS(IOB)
  257.                SEGACT IPT1
  258.             ENDIF
  259.             NBELEM=IPT1.NUM(/2)
  260.             NBNN=IPT1.NUM(/1)
  261.             NBREF=0
  262.             NBSOUS=0
  263.             SEGINI IPT4
  264.             IPT4.ITYPEL=IPT1.ITYPEL
  265.             DO 52 I=1,NBNN
  266.                DO 521 J=1,NBELEM
  267.                   IPT4.NUM(I,J)=0
  268.  521           CONTINUE
  269.  52         CONTINUE
  270.             ICO=0
  271.             DO 60 J=1,NBELEM
  272.                XP=0.D0
  273.                YP=0.D0
  274.                ZP=0.D0
  275.                DO 61 I=1,NBNN
  276.                   IP=IPT1.NUM(I,J)
  277.                   IREF=IP*(IDIM+1)-IDIM
  278.                   XP=XCOOR(IREF) +XP
  279.                   YP=XCOOR(IREF+1)+YP
  280.                   ZP=XCOOR(IREF+2)  +ZP
  281.  61            CONTINUE
  282.                IF (IDIM.EQ.2) ZP=0.D0
  283.                XP = XP / FLOAT(NBNN)
  284.                YP = YP / FLOAT(NBNN)
  285.                ZP = ZP / FLOAT(NBNN)
  286.                ANG=0.D0
  287.                DO 62 IT=1,ITOUR
  288.                   IREF=NFI(MAI(IT+1))*(IDIM+1)-IDIM
  289.                   XV1=XCOOR(IREF)-XP
  290.                   YV1=XCOOR(IREF+1)-YP
  291.                   ZV1=XCOOR(IREF+2)-ZP
  292.                   IF (IDIM.EQ.2) ZV1=0.
  293.                   DO 621 M=MAI(IT-1+1)+1,MAI(IT+1)
  294.                      IREF=NFI(M)*(IDIM+1)-IDIM
  295.                      XV2=XCOOR(IREF)-XP
  296.                      YV2=XCOOR(IREF+1)-YP
  297.                      ZV2=XCOOR(IREF+2)-ZP
  298.                      IF (IDIM.EQ.2) ZV2=0.
  299.                      V1=XNORM*(YV1*ZV2-ZV1*YV2)+YNORM*(ZV1*XV2-XV1*ZV2)
  300.      $                    +ZNORM*(XV1*YV2-YV1*XV2)
  301.                      V2=XV1*XV2+YV1*YV2+ZV1*ZV2
  302.                      IF (V1.EQ.0.D0.AND.V2.EQ.0.D0) GOTO 171
  303.                      ANG=ANG+ATAN2(V1,V2)
  304.  171                 CONTINUE
  305.                      XV1=XV2
  306.                      YV1=YV2
  307.                      ZV1=ZV2
  308.  621              CONTINUE
  309.  62            CONTINUE
  310.                IF (ABS(ANG).GT.XPI) GOTO 63
  311.                GOTO 60
  312.  63            CONTINUE
  313. *              Cet élément est ok, on le met de côté.
  314.                DO 64 L=1,NBNN
  315.                   IPT4.NUM(L,J)=IPT1.NUM(L,J)
  316.  64            CONTINUE
  317.                IPT4.ICOLOR(J)=IPT1.ICOLOR(J)
  318.                ICO=ICO+1
  319.  60         CONTINUE
  320.  
  321.             NBELEM=ICO
  322.             IF (NBELEM.EQ.0) GOTO 7051
  323.  
  324. *           sauvegarde partition ok
  325.             SEGINI IPT3
  326.             IPT3.ITYPEL=IPT1.ITYPEL
  327.             JJ=0
  328.             DO 71 J=1,IPT4.NUM(/2)
  329.                IF (IPT4.NUM(1,J).EQ.0) GOTO 71
  330.                JJ=JJ+1
  331.                IPT3.ICOLOR(JJ)=IPT4.ICOLOR(J)
  332.                DO 72 I=1,NBNN
  333.                   IPT3.NUM(I,JJ)=IPT4.NUM(I,J)
  334.  72            CONTINUE
  335.  71         CONTINUE
  336.             SEGSUP IPT4
  337.             IF (NBS.NE.0) THEN
  338.                IOB1=IOB1+1
  339.                IPT8.LISOUS(IOB1)=IPT3
  340.             ENDIF
  341.  7051    CONTINUE
  342.       ENDIF
  343.  
  344. *     Ecriture du maillage résultat et sortie
  345. *     ---------------------------------------
  346.       IF (((NBS.EQ.0).AND.(NBELEM.EQ.0)).OR.
  347.      &    ((NBS.NE.0).AND.(IOB1.EQ.0))) THEN
  348. *        on n'a rien trouvé
  349.          IF (IVERI.EQ.0) THEN
  350. *           Tache impossible. Probablement données erronées
  351.             CALL ERREUR(26)
  352.             RETURN
  353.          ELSE
  354. *           on renvoie un maillage vide
  355.             GOTO 7006
  356.          ENDIF
  357.       ENDIF
  358.       IF (NBS.EQ.0) THEN
  359. *        Si maillage simple en entrée
  360. *        On écrit directement IPT3
  361.          GOTO 7002
  362.       ELSE
  363. *        Si maillage complexe en entrée
  364. *        IOB1 : nombre de partitions retenues
  365.          IF     (IOB1.EQ.NBS) THEN
  366.             GOTO 7003
  367.          ELSEIF (IOB1.EQ.1)   THEN
  368.             GOTO 7005
  369.          ELSE
  370.             GOTO 7004
  371.          ENDIF
  372.       ENDIF
  373.  
  374. *     Maillage complexe, avec moins de partitions que l'objet entré
  375.  7004 CONTINUE
  376.       NBSOUS=IOB1
  377.       NBREF=0
  378.       NBNN=0
  379.       NBELEM=0
  380.       SEGINI IPT3
  381.       DO IOB=1,NBSOUS
  382.          IPT3.LISOUS(IOB)=IPT8.LISOUS(IOB)
  383.       ENDDO
  384.       SEGSUP IPT8
  385.       GOTO 7002
  386.  
  387. *     Maillage simple, avec moins de partitions que l'objet entré
  388.  7005 CONTINUE
  389.       IPT3=IPT8.LISOUS(1)
  390.       SEGSUP IPT8
  391.       GOTO 7002
  392.  
  393. *     Maillage avec autant de partitions que l'objet entré
  394. *     (inutile de recréer un nouveau maillage)
  395.  7003 CONTINUE
  396.       IPT3=IPT8
  397.       GOTO 7002
  398.  
  399. *     Maillage vide
  400.  7006 CONTINUE
  401.       NBSOUS=0
  402.       NBREF=0
  403.       NBNN=0
  404.       NBELEM=0
  405.       SEGINI IPT3
  406.       IF (NBS.NE.0) SEGSUP IPT8
  407.  
  408.  7002 CONTINUE
  409.       SEGSUP ICPR,FER
  410.       CALL ACTOBJ('MAILLAGE',IPT3,1)
  411.       CALL ECROBJ('MAILLAGE',IPT3)
  412.  
  413.       END
  414.  
  415.  
  416.  
  417.  

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