Télécharger la plaquette Cast3M


fr en

Télécharger vercon.eso

Retour à la liste

Numérotation des lignes : 

vercon

  1. C VERCON    SOURCE    CB215821  16/12/05    22:28:37     9237           
  2.       SUBROUTINE VERCON(IPT1,ICONV)
  3. C-----------------------------------------------------------------------
  4. C
  5. C         VERIFICATION DE LA CONVEXITE D'UN MAILLAGE (2D, 3D)
  6. C
  7. C-----------------------------------------------------------------------
  8. C
  9. C     Appelé par VERMAI
  10. C
  11. C     Entrée :
  12. C       IPT1 : pointeur sur le maillage (MELEME) a tester.
  13. C              Ce maillage doit respecter certaines conditions:
  14. C              - en dimension 2, le contour doit etre fait de SEG2
  15. C              - en dimension 3, l'enveloppe doit etre faite de TRI3
  16. C
  17. C     Sortie :
  18. C       ICONV : entier indiquant si le maillage est convexe
  19. C               = 1 si le maillage est convexe
  20. C               = 0 sinon
  21. C
  22. C       Un message est emis si le contour, ou de l'enveloppe, ne
  23. C       respecte pas les conditions requises
  24. C
  25. C
  26. C     Note : en dimension 2, on ne supprime pas le MELEME IPT2 (contour)
  27. C            car lors de l'appel a PRCONT, le contour IPT2 est reference
  28. C            dans le maillage IPT1
  29.  
  30. C-----------------------------------------------------------------------
  31.       IMPLICIT INTEGER(I-N)
  32.       IMPLICIT REAL*8(A-H,O-Z)
  33. -INC SMCOORD
  34.  
  35. -INC PPARAM
  36. -INC CCOPTIO
  37. -INC SMELEME
  38. -INC SMLENTI
  39. C
  40.       SEGMENT,MADJACEN
  41.         INTEGER LEAC(NBL1,IDIM,2)
  42.       ENDSEGMENT
  43.       POINTEUR ILEA1.MADJACEN
  44. C
  45.       POINTEUR LNPT1.MLENTI
  46. C
  47. C     Tolerance pour la planeite des element adjacents
  48.       ZERO=1.E-14
  49. C
  50. C     1. CONSTRUCTION DU CONTOUR OU DE L'ENVELOPPE DU MAILLAGE
  51. C
  52. C     J'ecris le maillage dans la pile
  53.       CALL ECROBJ('MAILLAGE',IPT1)
  54.       IF (IERR.NE.0) RETURN
  55. C     puis je fais appel a l'operateur ad hoc pour calculer :
  56. C     - soit le contour, en dimension 2
  57. C     - soit l'enveloppe, en dimension 3
  58.       IF (IDIM.EQ.2) THEN
  59.         CALL PRCONT
  60.       ELSEIF (IDIM.EQ.3) THEN
  61.         CALL ENVVOL
  62.       ENDIF
  63. C     Je lis le contour/enveloppe dans la pile, je le designe par IPT2
  64.       CALL LIROBJ('MAILLAGE',IPT2,1,IRETOU)
  65.       IF (IERR.NE.0) RETURN
  66.       SEGACT,IPT2
  67. C
  68. C     2. TESTS SUR LE MAILLAGE DU CONTOUR/ENVELOPPE
  69. C
  70. C     Test si le type d'elements est bien celui attendu, a savoir:
  71. C     des TRI3 en dimension 3 ou bien des SEG2 en dimension 2
  72.       IF (IDIM.EQ.2) THEN
  73.         IF ((IPT2.ITYPEL).NE.2) THEN
  74.           ICONV=0
  75.           MOTERR(1:8) = 'utilise '
  76.           CALL ERREUR(926)
  77.           RETURN
  78.         ENDIF
  79.       ELSEIF (IDIM.EQ.3) THEN
  80.         IF ((IPT2.ITYPEL).NE.4) THEN
  81.           ICONV=0
  82.           MOTERR(1:8) = 'utilise '
  83.           CALL ERREUR(926)
  84.           RETURN
  85.         ENDIF
  86.       ENDIF
  87. C
  88. C     3. VERIFICATION DE LA CONVEXITE
  89. C
  90. C     3.1 EN DIMENSION 2
  91. C         On parcourt le contour de proche en proche et on verifie la
  92. C         convexite: le contour doit tourner toujours dans le meme sens
  93. C
  94.       IF (IDIM.EQ.2) THEN
  95.         I1=0
  96.         JG=IPT2.NUM(/2)
  97.         SEGINI,LNPT1
  98. C       On part du premier element du contour
  99.         NEL1=1
  100.         NNOA1=1
  101.         NNOB1=2
  102. C       Boucle sur les elements (indice I)
  103.         DO I=1,JG
  104.           LNPT1.LECT(I)=1
  105. C         Coordonnees des points A et B, extremites de l'element I
  106.           NREFA=IPT2.NUM(NNOA1,NEL1)
  107.           XA=XCOOR((NREFA-1)*(IDIM+1)+1)
  108.           YA=XCOOR((NREFA-1)*(IDIM+1)+2)
  109.           NREFB=IPT2.NUM(NNOB1,NEL1)
  110.           XB=XCOOR((NREFB-1)*(IDIM+1)+1)
  111.           YB=XCOOR((NREFB-1)*(IDIM+1)+2)
  112. C         Le voisin du dernier element est le premier
  113.           IF (I.EQ.JG) THEN
  114.             NEL2=1
  115.             NNOB2=1
  116.             NNOC2=2
  117.             GOTO 10
  118.           ENDIF
  119. C         Boucle sur les autres elements du contour
  120. C         on recherche le voisin de I partageant le point B
  121.           DO J=1,JG
  122. C           Si on a deja traite l'element J, alors on itere
  123.             IF (LNPT1.LECT(J).EQ.1) GOTO 20
  124. C           Boucle sur les noeuds de l'element J
  125.             DO NNOB2=1,IDIM
  126. C             On teste si le noeud NNOB2 correspond au noeud B
  127.               IF ((IPT2.NUM(NNOB2,J)).EQ.NREFB) THEN
  128.                 NEL2=J
  129.                 IF (NNOB2.EQ.1) NNOC2=2
  130.                 IF (NNOB2.EQ.2) NNOC2=1
  131.                 GOTO 10
  132.               ENDIF
  133.             ENDDO
  134. 20        CONTINUE
  135.           ENDDO
  136. 10        CONTINUE
  137. C         Coordonnees du point C
  138.           NREFC=IPT2.NUM(NNOC2,NEL2)
  139.           XC=XCOOR((NREFC-1)*(IDIM+1)+1)
  140.           YC=XCOOR((NREFC-1)*(IDIM+1)+2)
  141. C         Calcul de la composante Z du produit vectoriel AB^BC
  142. C         => on obtient ainsi le sens de rotation
  143.           ZPV=((XB-XA)*(YC-YB))-((YB-YA)*(XC-XB))
  144. C          CB215821 : Inutile de normer le resultat ==> SIGFPE si 2 pts
  145. C                     sont confondus... Le denominateur est toujours 
  146. C                     positif ou nul donc aucun interet de normer
  147. C          X0A=((XA**2)+(YA**2))**0.5
  148. C          X0B=((XB**2)+(YC**2))**0.5
  149. C          X0C=((XC**2)+(YC**2))**0.5
  150. C          ZPV1=ZPV/(X0A*X0B*X0C)
  151. C         Determination du sens de rotation initial
  152.           IF (I1.EQ.0) THEN
  153. C           Si le produit vectoriel est nul, on ne peut pas connaitre le
  154. C           sens initial, alors on passe a l'element suivant
  155.             IF (ABS(ZPV).LT.ZERO) THEN
  156.               GOTO 30
  157.             ELSE
  158.               SIG_0=SIGN(1D0,ZPV)
  159.               I1=1
  160.             ENDIF
  161.           ELSE
  162. C           Si le sens de rotation change, par rapport a l'initial,
  163. C           alors le maillage n'est pas convexe => on sort
  164.             IF (ABS(ZPV).LT.ZERO) THEN
  165.               GOTO 30
  166.             ELSE
  167.               SIG_I=SIGN(1D0,ZPV)
  168.               IF (SIG_I.NE.SIG_0) THEN
  169.                 ICONV=0
  170.                 SEGSUP,LNPT1
  171.                 RETURN
  172.               ENDIF
  173.             ENDIF
  174.           ENDIF
  175. 30        CONTINUE
  176. C         On passe a l'element suivant, B devient A
  177.           NEL1=NEL2
  178.           NNOA1=NNOB2
  179.           IF (NNOA1.EQ.1) NNOB1=2
  180.           IF (NNOA1.EQ.2) NNOB1=1
  181.         ENDDO
  182.         ICONV=1
  183.         SEGSUP,LNPT1
  184.         RETURN
  185. C
  186. C     3.2 EN DIMENSION 3
  187. C         On parcoure toutes les facettes de l'enveloppe et on recherche
  188. C         les 3 elements adjacents. Une fois trouves, on teste la
  189. C         convexite: les noeuds des elements adjacents doivent tous etre
  190. C         du meme cote du plan definit par la facette
  191. C
  192.       ELSEIF (IDIM.EQ.3) THEN
  193.         I1=0
  194.         NBL1=IPT2.NUM(/2)
  195.         SEGINI,ILEA1
  196. C       Boucle sur les elements (indice I)
  197.         DO I=1,NBL1
  198. C         Numeros des noeuds A, B et C de l'element I
  199.           NREFA=IPT2.NUM(1,I)
  200.           NREFB=IPT2.NUM(2,I)
  201.           NREFC=IPT2.NUM(3,I)
  202. C         Numeros des adjacents et des noeuds D, E et F
  203.           NELD=0
  204.           NELE=0
  205.           NELF=0
  206.           NREFD=0
  207.           NREFE=0
  208.           NREFF=0
  209. C         On part a la recherche des adjacents de l'element I
  210. C         Si des indices de la table de connectivite existent deja,
  211. C         alors  on va chercher les numeros des elements adjacents et
  212. C         ceux des noeuds D, E ou F
  213.           IF ((ILEA1.LEAC(I,1,1)).NE.0) THEN
  214.             NELD=ILEA1.LEAC(I,1,1)
  215.             NREFD=IPT2.NUM((ILEA1.LEAC(I,1,2)),NELD)
  216.           ENDIF
  217.           IF ((ILEA1.LEAC(I,2,1)).NE.0) THEN
  218.             NELE=ILEA1.LEAC(I,2,1)
  219.             NREFE=IPT2.NUM((ILEA1.LEAC(I,2,2)),NELE)
  220.           ENDIF
  221.           IF ((ILEA1.LEAC(I,3,1)).NE.0) THEN
  222.             NELF=ILEA1.LEAC(I,3,1)
  223.             NREFF=IPT2.NUM((ILEA1.LEAC(I,3,2)),NELF)
  224.           ENDIF
  225. C         Boucle sur les elements d'indice superieurs a I (indice II)
  226.           DO II=I+1,NBL1
  227. C           Si l'element II est un des adjacents connus, on itere
  228.             IF ((II.EQ.NELD).OR.(II.EQ.NELE).OR.(II.EQ.NELF)) GOTO 40
  229. C           Indicateurs si l'un des noeuds de II correspond a A, B ou C
  230.             IA=0
  231.             IB=0
  232.             IC=0
  233. C           Boucle sur les noeuds (indice JJ) de l'element II
  234.             DO JJ=1,IDIM
  235.               NREFX=IPT2.NUM(JJ,II)
  236. C             On teste si le noeud JJ correspond aux noeuds A, B ou C
  237.               IF (NREFX.EQ.NREFA) IA=JJ
  238.               IF (NREFX.EQ.NREFB) IB=JJ
  239.               IF (NREFX.EQ.NREFC) IC=JJ
  240. C             A t on trouve l'adjacent vu par A ?
  241.               IF ((IB.NE.0).AND.(IC.NE.0)) THEN
  242.                 IF ((IB+IC).EQ.3) NNOJJ=3
  243.                 IF ((IB+IC).EQ.4) NNOJJ=2
  244.                 IF ((IB+IC).EQ.5) NNOJJ=1
  245.                 NREFD=IPT2.NUM(NNOJJ,II)
  246.                 ILEA1.LEAC(I,1,1)=II
  247.                 ILEA1.LEAC(I,1,2)=NNOJJ
  248.                 ILEA1.LEAC(II,NNOJJ,1)=I
  249.                 ILEA1.LEAC(II,NNOJJ,2)=1
  250.               ENDIF
  251. C             A t on trouve l'adjacent vu par B ?
  252.               IF ((IA.NE.0).AND.(IC.NE.0)) THEN
  253.                 IF ((IA+IC).EQ.3) NNOJJ=3
  254.                 IF ((IA+IC).EQ.4) NNOJJ=2
  255.                 IF ((IA+IC).EQ.5) NNOJJ=1
  256.                 NREFE=IPT2.NUM(NNOJJ,II)
  257.                 ILEA1.LEAC(I,2,1)=II
  258.                 ILEA1.LEAC(I,2,2)=NNOJJ
  259.                 ILEA1.LEAC(II,NNOJJ,1)=I
  260.                 ILEA1.LEAC(II,NNOJJ,2)=2
  261.               ENDIF
  262. C             A t on trouve l'adjacent vu par C ?
  263.               IF ((IA.NE.0).AND.(IB.NE.0)) THEN
  264.                 IF ((IA+IB).EQ.3) NNOJJ=3
  265.                 IF ((IA+IB).EQ.4) NNOJJ=2
  266.                 IF ((IA+IB).EQ.5) NNOJJ=1
  267.                 NREFF=IPT2.NUM(NNOJJ,II)
  268.                 ILEA1.LEAC(I,3,1)=II
  269.                 ILEA1.LEAC(I,3,2)=NNOJJ
  270.                 ILEA1.LEAC(II,NNOJJ,1)=I
  271.                 ILEA1.LEAC(II,NNOJJ,2)=3
  272.               ENDIF
  273.             ENDDO
  274. 40          CONTINUE
  275. C           A t on trouve tous les adjacents ?
  276.             IF ((NREFD*NREFE*NREFF).NE.0) GOTO 50
  277.           ENDDO
  278. C         On a trouve les 3 voisins de l'element I
  279. 50        CONTINUE
  280. C         Coordonnees des noeuds A, B, C, D, E et F
  281.           XA=XCOOR((NREFA-1)*(IDIM+1)+1)
  282.           YA=XCOOR((NREFA-1)*(IDIM+1)+2)
  283.           ZA=XCOOR((NREFA-1)*(IDIM+1)+3)
  284.           XB=XCOOR((NREFB-1)*(IDIM+1)+1)
  285.           YB=XCOOR((NREFB-1)*(IDIM+1)+2)
  286.           ZB=XCOOR((NREFB-1)*(IDIM+1)+3)
  287.           XC=XCOOR((NREFC-1)*(IDIM+1)+1)
  288.           YC=XCOOR((NREFC-1)*(IDIM+1)+2)
  289.           ZC=XCOOR((NREFC-1)*(IDIM+1)+3)
  290.           XD=XCOOR((NREFD-1)*(IDIM+1)+1)
  291.           YD=XCOOR((NREFD-1)*(IDIM+1)+2)
  292.           ZD=XCOOR((NREFD-1)*(IDIM+1)+3)
  293.           XE=XCOOR((NREFE-1)*(IDIM+1)+1)
  294.           YE=XCOOR((NREFE-1)*(IDIM+1)+2)
  295.           ZE=XCOOR((NREFE-1)*(IDIM+1)+3)
  296.           XF=XCOOR((NREFF-1)*(IDIM+1)+1)
  297.           YF=XCOOR((NREFF-1)*(IDIM+1)+2)
  298.           ZF=XCOOR((NREFF-1)*(IDIM+1)+3)
  299. C         Calcul d'un vecteur normal au plan ABC, on prend le produit
  300. C         vectoriel n = AB^AC
  301.           XN=((YB-YA)*(ZC-ZA))-((ZB-ZA)*(YC-YA))
  302.           YN=((ZB-ZA)*(XC-XA))-((XB-XA)*(ZC-ZA))
  303.           ZN=((XB-XA)*(YC-YA))-((YB-YA)*(XC-XA))
  304. C         Calcul des produits scalaires n.AD n.AE et n.AF pour connaitre
  305. C         la position de D, E et F par rapport au plan ABC
  306.           PSD=(XN*(XD-XA))+(YN*(YD-YA))+(ZN*(ZD-ZA))
  307.           PSE=(XN*(XE-XA))+(YN*(YE-YA))+(ZN*(ZE-ZA))
  308.           PSF=(XN*(XF-XA))+(YN*(YF-YA))+(ZN*(ZF-ZA))
  309.           X0N=((XN**2)+(YN**2)+(ZN**2))**0.5
  310.           X0AD=(((XD-XA)**2)+((YD-YA)**2)+((ZD-ZA)**2))**0.5
  311.           X0AE=(((XE-XA)**2)+((YE-YA)**2)+((ZE-ZA)**2))**0.5
  312.           X0AF=(((XF-XA)**2)+((YF-YA)**2)+((ZF-ZA)**2))**0.5
  313.           PSD1=PSD/(X0N*X0AD)
  314.           PSE1=PSE/(X0N*X0AE)
  315.           PSF1=PSF/(X0N*X0AF)
  316. C         Cas particulier ou le point D est dans le plan ABC
  317.           IF (ABS(PSD1).LT.ZERO) THEN
  318.             SIGD=0.
  319.           ELSE
  320.             SIGD=SIGN(1D0,PSD)
  321.           ENDIF
  322. C         Cas particulier ou le point E est dans le plan ABC
  323.           IF (ABS(PSE1).LT.ZERO) THEN
  324.             SIGE=0.
  325.           ELSE
  326.             SIGE=SIGN(1D0,PSE)
  327.           ENDIF
  328. C         Cas particulier ou le point F est dans le plan ABC
  329.           IF (ABS(PSF1).LT.ZERO) THEN
  330.             SIGF=0.
  331.           ELSE
  332.             SIGF=SIGN(1D0,PSF)
  333.           ENDIF
  334. C         Si les 3 points ne sont pas tous situe du meme cote du plan
  335. C         ABC, alors le maillage n'est pas convexe => on sort
  336. C         - si D seulement est dans le plan ABC, on teste juste E et F
  337.           IF ((SIGD.EQ.0.).AND.(SIGE.NE.0.).AND.(SIGF.NE.0.)) THEN
  338.             IF (SIGE.NE.SIGF) THEN
  339.               ICONV=0
  340.               SEGSUP,ILEA1,IPT2
  341.               RETURN
  342.             ENDIF
  343.           ENDIF
  344. C         - si E seulement est dans le plan ABC, on teste juste D et F
  345.           IF ((SIGE.EQ.0.).AND.(SIGD.NE.0.).AND.(SIGF.NE.0.)) THEN
  346.             IF (SIGD.NE.SIGF) THEN
  347.               ICONV=0
  348.               SEGSUP,ILEA1,IPT2
  349.               RETURN
  350.             ENDIF
  351.           ENDIF
  352. C         - si F seulement est dans le plan ABC, on teste juste D et E
  353.           IF ((SIGF.EQ.0.).AND.(SIGD.NE.0.).AND.(SIGE.NE.0.)) THEN
  354.             IF (SIGD.NE.SIGE) THEN
  355.               ICONV=0
  356.               SEGSUP,ILEA1,IPT2
  357.               RETURN
  358.             ENDIF
  359.           ENDIF
  360. C         - cas general (D, E et F hors du plan ABC)
  361.           IF ((SIGD.NE.0.).AND.(SIGE.NE.0.).AND.(SIGF.NE.0.)) THEN
  362.             IF ((SIGD.NE.SIGE).OR.(SIGE.NE.SIGF)) THEN
  363.               ICONV=0
  364.               SEGSUP,ILEA1,IPT2
  365.               RETURN
  366.             ENDIF
  367.           ENDIF
  368. C         - dans tous les autres cas (2 ou 3 des points D, E, F sont
  369. C           dans le plan ABC), on itere
  370.         ENDDO
  371.         ICONV=1
  372.         SEGSUP,ILEA1,IPT2
  373.       ENDIF
  374.       RETURN
  375.       END
  376.  
  377.  
  378.  
  379.  

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