Télécharger la plaquette Cast3M


fr en

Télécharger kalpre.eso

Retour à la liste

Numérotation des lignes : 

kalpre

  1. C KALPRE    SOURCE    OF166741  24/10/03    21:15:23     12022          
  2.       SUBROUTINE KALPRE (MYMOD,INFOEL,SKFACE,XDEC,NELD)
  3.       IMPLICIT INTEGER(I-N)
  4.       IMPLICIT REAL*8 (A-H,O-Z)
  5. C----------------------------------------------------------------------
  6. C     Facteurs de forme OPTION CACHE 3D
  7. C     Calcul du decoupage des faces et initialisations
  8. C     SP appelle par facgen
  9. C   entree:
  10. C     MYMOD  : pointeur de l'objet modèle
  11. C     INFOEL : informations concernant le type des éléments des maillages
  12. C     SKFACE : pointeur sur l'objet 'faces' pour le calcul des FF
  13. C     SHC3D  : pointeur sur la surface de projection
  14. C     XDEC   : paramŠtre pour le decoupage des faces
  15. C     NELD   : nombre d'éléments virtuels
  16. C             (un élément COQ -> 2 facteurs de forme <=> 2 éléments virtuels)
  17. C   sortie:
  18. C     SKFACE : pointeur sur l'objet 'faces' pour le calcul des FF
  19. C     SHC3D  : pointeur sur la surface de projection
  20. C----------------------------------------------------------------------
  21.  
  22. -INC PPARAM
  23. -INC CCOPTIO
  24.  
  25. -INC SMCOORD
  26. -INC SMMODEL
  27. -INC SMELEME
  28. -INC TFFOR3D
  29.  
  30. C     -----------------------------------
  31.       SEGMENT SDECOU
  32.       INTEGER KDECOU(MFACE)
  33.       ENDSEGMENT
  34.       POINTEUR MYMOD.MMODEL
  35. C     -----------------------------------
  36. C     Stockage d'informations concernant le type des éléments des maillages
  37.       SEGMENT ,INFOEL
  38.         LOGICAL KCOQ(N1),KQUAD(N1)
  39.       ENDSEGMENT
  40. C     -----------------------------------
  41.       DIMENSION A(3,4),UU(4),IN(3)
  42.       DIMENSION A1(3,3),A2(3,3),U1(4),U2(4),G1(3),G2(3)
  43.       DIMENSION GG(3,400),SS(400)
  44.       LOGICAL ICOQ
  45. C
  46. C--------------------------------------------------------------------
  47. C
  48. C
  49.       DMIN=1.E-5
  50.       DMAX=1./DMIN
  51. C... nombre max de triangles de subdivision
  52.       NPM=20
  53. C... quadratiques
  54.       NSPA=1
  55.  
  56.       IF(INFOEL.EQ.0) THEN
  57.          ICOQ=.FALSE.
  58.       ELSE
  59.          ICOQ=.TRUE.
  60.          SEGACT INFOEL
  61.       ENDIF
  62. C
  63. C>>>  CREATION DE L'OBJET FACE
  64. C
  65.       NFACE = 0
  66.       NELD = 0
  67. C
  68.       MTYP = MYMOD.KMODEL(/1)
  69.       SEGACT INFOEL
  70. C
  71. C     On va compter le nombre d'éléments virtuels (égal au nombre de
  72. C     facteurs de forme par ligne) : NELD , et le nombre de faces : NFACE .
  73.       DO 10 ITYP=1,MTYP
  74.          IMODE1 = MYMOD.KMODEL(ITYP)
  75.          IPT1 = IMODE1.IMAMOD
  76.          NEL = IPT1.NUM(/2)
  77.          NPGEO  = IPT1.NUM(/1)
  78. C
  79. C        un élément COQ -> 2 facteurs de formes
  80.          IF (KCOQ(ITYP)) NEL=2*NEL
  81.          NELD = NELD + NEL
  82. C
  83. C        un élément carré -> 2 éléments triangles
  84.          IF(NPGEO.EQ.4.OR.NPGEO.EQ.8)  NEL = 2*NEL
  85.          NFACE=NFACE+NEL
  86.  10   CONTINUE
  87. C
  88.       IF (KIMP.GE.1) THEN
  89.       WRITE( 6,*) ' NOMBRE TOTAL D ELEMENTS ',NELD,'  DE FACES ',NFACE
  90.       ENDIF
  91. C
  92. C
  93.       MFACE = NFACE
  94.       MS = 3
  95.       SEGINI SKFACE
  96.       SEGINI SDECOU
  97. C
  98. C     KEL : INDICE ELEMENT DECOUPAGE (-> numéro global de la facette)
  99. C     IEL : INDICE ELEMENT D'ORIGINE (-> numero global de l'élément virtuel)
  100. C     Ces indices sont différents car on découpe les quadrangles en triangles
  101. C
  102.       KEL = 1
  103.       IEL = 0
  104. C
  105.       DO 100 ITYP = 1, MTYP
  106.          IMODE1 = MYMOD.KMODEL(ITYP)
  107.          IPT1 = IMODE1.IMAMOD
  108.  
  109.          NSGEO  = IPT1.NUM(/1)
  110.          IF(ICOQ.AND.KQUAD(ITYP)) THEN
  111.            NS = NSGEO/2
  112.          ELSE
  113.            NS=NSGEO
  114.          ENDIF
  115.          IF(ICOQ.AND.KQUAD(ITYP)) NSPA=2
  116.          NEL = IPT1.NUM(/2)
  117. C
  118. C...     TRI3 ou TRI6
  119.          IF (NS.EQ.3) THEN
  120. *////////
  121. *           WRITE (6,*) 'Le maillage ',IPT1,
  122. *    #      'est constitué d éléments triangles'
  123. *////////
  124.             DO 110 I = 1,NEL
  125.                IF (KCOQ(ITYP)) THEN
  126.                    IEL1 = IEL + (2*I-1)
  127.                ELSE
  128.                    IEL1 = IEL + I
  129.                ENDIF
  130.                KCORR(KEL) = IEL1
  131.  
  132.                DO 111 IS = 1,NSGEO,NSPA
  133.                   LS=IS
  134.                   IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2
  135. C**               NUK(IS,KEL) = IPT1.NUM(IS,I)
  136.                   NUK(LS,KEL) = IPT1.NUM(IS,I)
  137.                   IREF = (IDIM+1)*(NUK(LS,KEL)-1)
  138.                      DO 112 K = 1,KES
  139. C**                  A(K,LS) = XCOOR(IREF+K)
  140.                      A(K,LS) = XCOOR(IREF+K)
  141.  112                 CONTINUE
  142. *////////
  143. *       WRITE (6,*) 'noeud numero',IPT1.NUM(IS,I)
  144. *       WRITE (6,*) 'Coordonnées :',(A(K,IS),K=1,KES)
  145. *////////
  146.  111           CONTINUE
  147.                IF (KIMP.GE.3) THEN
  148.         WRITE (6,*) 'La facette ',KEL,'repose sur les points',
  149.      #                  (IPT1.NUM(IS,I),IS=1,NS)
  150.                ENDIF
  151.                CALL KNORM(KES,A,NS,UU,S(KEL))
  152.                DO 113 L = 1,4
  153.                   U(L,KEL) = UU(L)
  154.  113           CONTINUE
  155. C     WRITE(6,*) ' KEL  S U ',KEL,S(KEL),(U(K,NEL),K=1,4)
  156.  
  157.                KEL = KEL + 1
  158. *////////
  159.                IF (KCOQ(ITYP)) THEN
  160. C                  On remplit KCOOR , NUK , U et S pour l'élément inverse
  161.                    KCORR(KEL) = IEL1 + 1
  162.                    NUK(1,KEL) = NUK(3,KEL-1)
  163.                    NUK(2,KEL) = NUK(2,KEL-1)
  164.                    NUK(3,KEL) = NUK(1,KEL-1)
  165.                    DO 114 L = 1,4
  166.                       U(L,KEL) = - UU(L)
  167.  114               CONTINUE
  168.                    S(KEL) = S(KEL-1)
  169.                    KEL = KEL + 1
  170.                    ENDIF
  171. C
  172.  110        CONTINUE
  173. C
  174.             IF (KCOQ(ITYP)) THEN
  175.                 IEL = IEL + NEL*2
  176.             ELSE
  177.                 IEL = IEL + NEL
  178.             ENDIF
  179. C
  180. C...     QUA4 ou QUA8
  181.          ELSE
  182. *////////
  183. *           WRITE (6,*) 'Le maillage ',IPT1,
  184. *    #      'est constitué d éléments carrés'
  185. *////////
  186. C
  187.             DO 120 I = 1,NEL
  188. C
  189.             IF (KCOQ(ITYP)) THEN
  190.                 IEL1 = IEL + (2*I-1)
  191.             ELSE
  192.                 IEL1 = IEL + I
  193.             ENDIF
  194.                DO 121 IS = 1,NSGEO,NSPA
  195.                   LS=IS
  196.                   IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2
  197.                      IREF = (IDIM+1)*(IPT1.NUM(IS,I)-1)
  198. *////////
  199. *               WRITE (6,*) 'noeud numero',IPT1.NUM(IS,I)
  200. *               WRITE (6,*) 'IREF =',IREF
  201. *////////
  202.                      DO 122 K = 1,KES
  203. C**                  A(K,IS) = XCOOR(IREF+K)
  204.                      A(K,LS) = XCOOR(IREF+K)
  205.  122                 CONTINUE
  206. *////////
  207. *       WRITE (6,*) 'noeud numero',IPT1.NUM(IS,I)
  208. *       WRITE (6,*) 'Coordonnées :',(A(K,IS),K=1,KES)
  209. *////////
  210.  121           CONTINUE
  211. C
  212. C              On regarde comment on va découper le quadrangle
  213.                DIA1 = 0.
  214.                      DO 123 K = 1,KES
  215.                      DIA1 = DIA1 + (A(K,1)-A(K,3))**2
  216.  123                 CONTINUE
  217.                DIA2 = 0.
  218.                      DO 124 K = 1,KES
  219.                      DIA2 = DIA2 + (A(K,2)-A(K,4))**2
  220.  124                 CONTINUE
  221. *//////
  222. *     WRITE(6,*) ' DIA1 DIA2 ',DIA1,DIA2
  223. *//////
  224.  
  225. C          stockage des connectivités des traiangles
  226.  
  227.            IF(ICOQ.AND.KQUAD(ITYP)) THEN
  228.  
  229. * cas des qua8
  230.  
  231.                IF (DIA1.LE.DIA2) THEN
  232.                      NUK(1,KEL) = IPT1.NUM(1,I)
  233.                      NUK(2,KEL) = IPT1.NUM(3,I)
  234.                      NUK(3,KEL) = IPT1.NUM(5,I)
  235.                      KEL = KEL + 1
  236.                      IF (KCOQ(ITYP)) THEN
  237.                          NUK(1,KEL) = IPT1.NUM(5,I)
  238.                          NUK(2,KEL) = IPT1.NUM(3,I)
  239.                          NUK(3,KEL) = IPT1.NUM(1,I)
  240.                          KEL = KEL + 1
  241.                      ENDIF
  242.                      NUK(1,KEL) = IPT1.NUM(5,I)
  243.                      NUK(2,KEL) = IPT1.NUM(7,I)
  244.                      NUK(3,KEL) = IPT1.NUM(1,I)
  245.                      KEL = KEL + 1
  246.                      IF (KCOQ(ITYP)) THEN
  247.                          NUK(1,KEL) = IPT1.NUM(1,I)
  248.                          NUK(2,KEL) = IPT1.NUM(7,I)
  249.                          NUK(3,KEL) = IPT1.NUM(5,I)
  250.                          KEL = KEL + 1
  251.                      ENDIF
  252.                ELSE
  253.                      NUK(1,KEL) = IPT1.NUM(3,I)
  254.                      NUK(2,KEL) = IPT1.NUM(5,I)
  255.                      NUK(3,KEL) = IPT1.NUM(7,I)
  256.                      KEL = KEL + 1
  257.                      IF (KCOQ(ITYP)) THEN
  258.                          NUK(1,KEL) = IPT1.NUM(7,I)
  259.                          NUK(2,KEL) = IPT1.NUM(5,I)
  260.                          NUK(3,KEL) = IPT1.NUM(3,I)
  261.                          KEL = KEL + 1
  262.                      ENDIF
  263.                      NUK(1,KEL) = IPT1.NUM(7,I)
  264.                      NUK(2,KEL) = IPT1.NUM(1,I)
  265.                      NUK(3,KEL) = IPT1.NUM(3,I)
  266.                      KEL = KEL + 1
  267.                      IF (KCOQ(ITYP)) THEN
  268.                          NUK(1,KEL) = IPT1.NUM(3,I)
  269.                          NUK(2,KEL) = IPT1.NUM(1,I)
  270.                          NUK(3,KEL) = IPT1.NUM(7,I)
  271.                          KEL = KEL + 1
  272.                      ENDIF
  273.                ENDIF
  274.  
  275.            ELSE
  276.  
  277. * cas des qua4
  278.  
  279.                IF (DIA1.LE.DIA2) THEN
  280.                      NUK(1,KEL) = IPT1.NUM(1,I)
  281.                      NUK(2,KEL) = IPT1.NUM(2,I)
  282.                      NUK(3,KEL) = IPT1.NUM(3,I)
  283.                      KEL = KEL + 1
  284.                      IF (KCOQ(ITYP)) THEN
  285.                          NUK(1,KEL) = IPT1.NUM(3,I)
  286.                          NUK(2,KEL) = IPT1.NUM(2,I)
  287.                          NUK(3,KEL) = IPT1.NUM(1,I)
  288.                          KEL = KEL + 1
  289.                      ENDIF
  290.                      NUK(1,KEL) = IPT1.NUM(3,I)
  291.                      NUK(2,KEL) = IPT1.NUM(4,I)
  292.                      NUK(3,KEL) = IPT1.NUM(1,I)
  293.                      KEL = KEL + 1
  294.                      IF (KCOQ(ITYP)) THEN
  295.                          NUK(1,KEL) = IPT1.NUM(1,I)
  296.                          NUK(2,KEL) = IPT1.NUM(4,I)
  297.                          NUK(3,KEL) = IPT1.NUM(3,I)
  298.                          KEL = KEL + 1
  299.                      ENDIF
  300.                ELSE
  301.                      NUK(1,KEL) = IPT1.NUM(2,I)
  302.                      NUK(2,KEL) = IPT1.NUM(3,I)
  303.                      NUK(3,KEL) = IPT1.NUM(4,I)
  304.                      KEL = KEL + 1
  305.                      IF (KCOQ(ITYP)) THEN
  306.                          NUK(1,KEL) = IPT1.NUM(4,I)
  307.                          NUK(2,KEL) = IPT1.NUM(3,I)
  308.                          NUK(3,KEL) = IPT1.NUM(2,I)
  309.                          KEL = KEL + 1
  310.                      ENDIF
  311.                      NUK(1,KEL) = IPT1.NUM(4,I)
  312.                      NUK(2,KEL) = IPT1.NUM(1,I)
  313.                      NUK(3,KEL) = IPT1.NUM(2,I)
  314.                      KEL = KEL + 1
  315.                      IF (KCOQ(ITYP)) THEN
  316.                          NUK(1,KEL) = IPT1.NUM(2,I)
  317.                          NUK(2,KEL) = IPT1.NUM(1,I)
  318.                          NUK(3,KEL) = IPT1.NUM(4,I)
  319.                          KEL = KEL + 1
  320.                      ENDIF
  321.                ENDIF
  322.  
  323.            ENDIF
  324.  
  325. C
  326.                DO 125 IL = 1,2
  327. C
  328.                   IF (KCOQ(ITYP)) THEN
  329.                       INDICE = KEL - IL*2
  330.                   ELSE
  331.                       INDICE = KEL - IL
  332.                   ENDIF
  333.                   KCORR(INDICE) = IEL1
  334. C
  335.                   DO 126 IS = 1,3
  336.                     IREF = (IDIM+1)*(NUK(IS,INDICE)-1)
  337.                      DO 127 K = 1,KES
  338.                         A(K,IS) = XCOOR(IREF+K)
  339. C                       WRITE(6,*) 'A ',K,IS,A(K,IS)
  340.  127                 CONTINUE
  341.  126              CONTINUE
  342.                   CALL KNORM(KES,A,3,UU,S(INDICE))
  343.                   DO 128 L = 1,4
  344.                      U(L,INDICE) = UU(L)
  345.  128              CONTINUE
  346. CC    WRITE(6,*) ' KEL S U ',INDICE,S(INDICE),(U(K,INDICE),K=1,4)
  347.  
  348.  125            CONTINUE
  349. C
  350.                 IF (KCOQ(ITYP)) THEN
  351. C               On remplit KCOOR , U et S pour les éléments inverses
  352.                    DO 135 IL = 1,2
  353.                       INDICE = KEL - IL*2 +1
  354.                       KCORR(INDICE) = IEL1 +1
  355.                       DO 138 L = 1,4
  356.                          U(L,INDICE) = - U(L,INDICE-1)
  357.  138                  CONTINUE
  358.                       S(INDICE)=S(INDICE-1)
  359.  135               CONTINUE
  360.                 ENDIF
  361. C
  362.  120        CONTINUE
  363. C
  364.             IF (KCOQ(ITYP)) THEN
  365.                 IEL = IEL + NEL*2
  366.             ELSE
  367.                 IEL = IEL + NEL
  368.             ENDIF
  369. C
  370. C
  371.          ENDIF
  372. C
  373.  100  CONTINUE
  374. C
  375. C     On a pris toute l'information nécessaire concernant le maillage
  376.       SEGDES INFOEL
  377. *////////
  378.         IF(KIMP.GE.2) THEN
  379.         WRITE (6,*) 'Subdivision'
  380.         ENDIF
  381. *////////
  382. C-----------------------------------------------------------
  383.  
  384.       DO 200 K1=1,MFACE
  385.  
  386.           DO 202 K = 1,KES
  387.           G1(K)=0.
  388.          DO 201 ISS = 1,MS
  389.           IREF = (IDIM+1)*(NUK(ISS,K1)-1)
  390.           A1(K,ISS) = XCOOR(IREF+K)
  391.           G1(K) = G1(K) + A1(K,ISS)
  392.  201     CONTINUE
  393.           G1(K)=G1(K)/3.
  394.  202      CONTINUE
  395.           DO 203 K=1,KES+1
  396.           U1(K) = U(K,K1)
  397.  203      CONTINUE
  398. C
  399.           IF (XDEC.GE.0.01) THEN
  400.           DK1 = DMAX
  401.             DO 400 K2 = 1,MFACE
  402. C
  403.               DO 412 K=1,KES+1
  404.               U2(K) = U(K,K2)
  405.  412          CONTINUE
  406.  
  407.                DO 212 K = 1,KES
  408.                G2(K)=0.
  409.                DO 211 ISS = 1,MS
  410.                IREF = (IDIM+1)*(NUK(ISS,K2)-1)
  411.                A2(K,ISS) = XCOOR(IREF+K)
  412.                G2(K)=G2(K)+A2(K,ISS)
  413.  211           CONTINUE
  414.                G2(K)=G2(K)/3.
  415.  212           CONTINUE
  416.  
  417. C>>>  VISIBILITE A PRIORI
  418. C
  419.  
  420.            CALL KPRIOR(KES,MS,MS,A1,A2,U1,U2,KVU)
  421.            IF (KVU.NE.0) THEN
  422.              D1=G1(1)-G2(1)
  423.              D2=G1(2)-G2(2)
  424.              D3=G1(3)-G2(3)
  425.              DKK1 = SQRT(D1*D1+D2*D2+D3*D3)
  426. C            WRITE(6,*) ' DKK1 CK1 ',DKK1,CK1
  427.                IF(DKK1.GE.DMIN) THEN
  428.                CK1  = ABS(U1(1)*D1+ U1(2)*D2 + U1(3)*D3)/DKK1
  429.                  IF(CK1.GE.0.01) THEN
  430.                    DK1 =MIN(DK1,DKK1)
  431.                  ENDIF
  432.                ENDIF
  433.            ENDIF
  434.  400   CONTINUE
  435.              DR = DK1/XDEC
  436.        ELSE
  437.              DR = DMAX
  438.        ENDIF
  439.  
  440. C
  441.  901  CONTINUE
  442.            CALL KCREPA(DR,A1,S(K1),KES,MS,NPAT,GG,SS)
  443.  
  444.       IF(NPAT.GT.NPM) THEN
  445.        DR=DR*2.
  446.  
  447.          IF (KIMP.GE.2) THEN
  448.           WRITE(6,900) K1,NPAT
  449.  900      FORMAT(1X,' FACE ',I4,'   : SUBDIVISION EXCESSIVE ',I4)
  450.          ENDIF
  451.        GOTO 901
  452.       ENDIF
  453.            NPATCH = NPAT
  454.            KDECOU(K1)=NPAT
  455.            MSP = KES
  456.            SEGINI IPATCH
  457.  
  458.            DO 130 IP = 1,NPAT
  459.                SP(IP) = SS(IP)/S(K1)
  460.                DO 129 K=1,KES
  461.                GP(K,IP) = GG(K,IP)
  462.  129           CONTINUE
  463. C              WRITE(6,*) ' SP GP ',SP(IP),GP(1,IP),GP(2,IP),GP(3,IP)
  464.  130       CONTINUE
  465.  
  466.            KPATCH(K1) = IPATCH
  467.            SEGDES IPATCH
  468.  
  469.  200   CONTINUE
  470. C-----------------------------------------------------------
  471.  
  472.          IF(KIMP.GE.2) THEN
  473.            WRITE(6,*)
  474.            WRITE(6,*) 'NOMBRE DE SUBDIVISIONS PAR FACE '
  475.            WRITE(6,1000) (KDECOU(I),I=1,MFACE)
  476.  1000      FORMAT(1X,10(I4))
  477.          ENDIF
  478. C
  479.       SEGDES SKFACE
  480.       SEGSUP SDECOU
  481. C
  482.       RETURN
  483.       END
  484.  
  485.  
  486.  

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