Télécharger la plaquette Cast3M


fr en

Télécharger pavpoi.eso

Retour à la liste

Numérotation des lignes : 

pavpoi

  1. C PAVPOI    SOURCE    BP208322  15/06/22    21:21:10     8543
  2.         SUBROUTINE PAVPOI(IWRK,IELE,NDEC,NNO)
  3. ****************************************************************************
  4. *      appelee par prom  projection par minimisation  d une integrale
  5. *      FONCTION :
  6. *      distribue une collection de point d integration dans
  7. *      les elements de reference et les rend a prom
  8. *  ENTREES  :
  9. *       NDEC nomdre de points dans une direction ( isotrope pour l instant *
  10. *       iele    numero du support geomtrique de l element
  11. *       NNO  nombre de points de l element ( utile pour polygone uniquement)
  12. *  SORTIES  :
  13. *       Tableau XX dans le segment iwrk  et NDD
  14. *****************************************************************************
  15.       IMPLICIT INTEGER(I-N)
  16.         IMPLICIT REAL*8 (A-H,O-Z)
  17.  
  18. -INC PPARAM
  19. -INC CCOPTIO
  20. -INC CCHAMP
  21. -INC CCREEL
  22. *-
  23.         PARAMETER ( O1=1.D0,O2=2.D0)
  24.         SEGMENT MTR
  25.         REAL*8 G(3,6),ppd(6)
  26.         ENDSEGMENT
  27.         SEGMENT WRK5
  28.         REAL*8 W(200),A(NBN1,NBN1),B(NBN1,NBN1),XPU(3),RHS(NBN1)
  29.         REAL*8 SHPXXX(6,NBNN,200),XX(3,200),VAL(NBN1),POIDS(200)
  30.         INTEGER NDD
  31.         ENDSEGMENT
  32. C
  33.         WRK5=IWRK
  34.         ip = 0
  35. C       write(6,*)' xx(/1) xx(/2)',xx(/1) ,xx(/2)
  36. C
  37.         goto (99,2, 2,4,4,4,4,8,8,8,8,100,100,8,8,4,4,99,99,
  38.      &  99,99,99,23,23,23,23,99,99,99,99,99,32,99,99,99,99,99,
  39.      &  99),iele
  40.  
  41. C  ELEMENT DE REFERENCE : SEG2 - SEG3
  42. C ------------------------
  43.  2    CONTINUE
  44.       delx=O2/ndec
  45.       x=-O1-delx/O2
  46.       DO i=1,ndec
  47.         xx(1,i)=x+delx*i
  48.         poids(i)=delx
  49.       ENDDO
  50.       ip=ndec
  51.       GOTO 100
  52.  
  53. 4      CONTINUE
  54. C       ELEMENT DE REFERENCE TRIANGLES +  BASE PRISME
  55.        ppp = O1/(O2*ndec*ndec)
  56.        r= O1/(3.D0*ndec)
  57.        p=O1/ndec
  58.        x0 =-p+r
  59.        y0 =x0
  60.        x1=x0+r
  61.        y1=x1
  62.        ppp = O1/(ndec*ndec*O2)
  63. C
  64.        do k=1,ndec
  65.        do i=1,ndec-k+1
  66.        ip=ip+1
  67.        xx(1,ip) = x0+p*i
  68.        xx(2,ip) = y0+p*k
  69.        poids(ip)= ppp
  70.        if(i.eq.ndec-k+1) goto 50
  71.        ip=ip+1
  72.        xx(1,ip) = x1+p*i
  73.        xx(2,ip) = y1+p*k
  74.        poids(ip)= ppp
  75. 50    continue
  76.        enddo
  77.        enddo
  78.         if(iele.eq.16.or.iele.eq.17) goto 16
  79.         GOTO 100
  80.  
  81. 8      CONTINUE
  82. C       ELEMENTS DE REFERENCE  QUADRANGLE  +   BASE CUBE
  83.       delx = O2/ndec
  84.       dely=delx
  85.       ppp = O2*O2/(ndec*ndec)
  86.       x= -O1 - delx/O2
  87.       y= -O1 + dely/O2
  88.       do i=1,ndec
  89.       do j=1,ndec
  90.       ip=ip+1
  91.       xx(1,ip)=x+delx*j
  92.       xx(2,ip)=y
  93.       poids(ip) = ppp
  94.       enddo
  95.       y = y + dely
  96.       enddo
  97.       if(iele.eq.14.or.iele.eq.15) goto 14
  98.       GOTO 100
  99. C
  100. 14     CONTINUE
  101. C      ELEMENTS  DE REFERENCE CUBES
  102.        delz = delx
  103.        ppp = ppp*delz
  104.        z0= -O1+delz/O2
  105. C
  106.        do i=1,ndec*ndec
  107.        xx(3,i)=z0
  108.        poids(i)=ppp
  109.        enddo
  110. C
  111.        ip= ndec*ndec
  112.         do k=1,ndec-1
  113.         do i=1,ndec*ndec
  114.        ip=ip+1
  115.        xx(1,ip)=xx(1,i)
  116.        xx(2,ip)=xx(2,i)
  117.        xx(3,ip)=xx(3,i)+delz*k
  118.        poids(ip)=ppp
  119.         enddo
  120.         enddo
  121.         GOTO 100
  122. C
  123. 16     CONTINUE
  124. C      ELEMENTS  DE REFRENCE PRISMES
  125.        delz = O2/ndec
  126.        z0= -O1+delz/O2
  127.        ppp = ppp*delz
  128. C
  129.        do i=1,ndec*ndec
  130.        xx(3,i)=z0
  131.        poids(i)=ppp
  132.        enddo
  133. C
  134.        ip= ndec*ndec
  135.         do k=1,ndec-1
  136.         do i=1,ndec*ndec
  137.        ip=ip+1
  138.        xx(1,ip)=xx(1,i)
  139.        xx(2,ip)=xx(2,i)
  140.        xx(3,ip)=xx(3,i)+delz*k
  141.        poids(ip)=ppp
  142.         enddo
  143.         enddo
  144.        GOTO 100
  145. 23     CONTINUE
  146. C      ELEMENTS  DE REFERENCE TETRAEDRES
  147.        p =O1/ndec
  148.        unq= O1/(O2*O2)
  149.        und=O1/O2
  150.        ppp = O1/(6.D0*ndec*ndec*ndec)
  151.        segini mtr
  152.        g(1,1) = unq
  153.        g(2,1) = unq
  154.        g(3,1) = unq
  155.        g(1,2) = O1-unq
  156.        g(2,2) = unq
  157.        g(3,2) = und
  158.        g(1,3) = unq
  159.        g(2,3) = O1-unq
  160.        g(3,3) = und
  161.        g(1,4) = und
  162.        g(2,4) = und
  163.        g(3,4) = O1-unq
  164.        g(1,5) = und
  165.        g(2,5) = und
  166.        g(3,5) = unq
  167.        g(1,6) = O1-unq
  168.        g(2,6) = O1-unq
  169.        g(3,6) = O1-unq
  170.         do i=1,6
  171.         do j=1,3
  172.         g(j,i)=g(j,i)/ndec-p
  173.         enddo
  174.         enddo
  175.  
  176.         do k=1,ndec
  177.         tz = p*k
  178.         do j=1,ndec-k+1
  179.         ty = p*j
  180.         do i=1,ndec-j+1-k+1
  181.         tx = p*i
  182. C
  183.         do 500 l=1,6
  184.         x=g(1,l)+tx
  185.         y=g(2,l)+ty
  186.         z=g(3,l)+tz
  187.         if((x+y+z).GT.O1) goto 500
  188.         ip = ip+1
  189.         xx(1,ip)=x
  190.         xx(2,ip)=y
  191.         xx(3,ip)=z
  192.         poids(ip) = ppp
  193. 500     continue
  194.         enddo
  195.         enddo
  196.         enddo
  197.         segsup mtr
  198.        IF (iele.eq.25.or.iele.eq.26 )  goto 25
  199.        GOTO 100
  200. 25     CONTINUE
  201. C      ELEMENTS  DE REFERENCE PYRAMIDES ( decoupee en 4 tetraedres)
  202.        ipo=ip
  203.         do ik=1,ipo
  204.        xx(1,ip+ik) = -xx(1,ik)
  205.        xx(2,ip+ik) =  xx(2,ik)
  206.        xx(3,ip+ik) =  xx(3,ik)
  207.        enddo
  208.        ipo=ipo*2
  209.        ip=ipo
  210.        do ik=1,ipo
  211.        xx(1,ip+ik) =  xx(1,ik)
  212.        xx(2,ip+ik) = -xx(2,ik)
  213.        xx(3,ip+ik) =  xx(3,ik)
  214.        enddo
  215.        ip = ipo*2
  216. C                    les poids
  217.        do ik=1,ip
  218.        poids(ik)=ppp
  219.        enddo
  220.        goto 100
  221. 32     continue
  222. C    -----------   les element polygonaux
  223.        alpa2= xpi/nno
  224.        alpa = alpa2*O2
  225.        px=O1/ndec
  226.        py =px * sin(alpa)
  227.        ppp = sin(alpa)/(ndec*ndec)/o2
  228.        rx= O2/3.D0*cos(alpa2)*cos(alpa2)/ndec
  229.        ry =O2/3.D0*cos(alpa2)*sin(alpa2)/ndec
  230. C
  231.        x0 =-px+rx
  232.        y0 =-py+ry
  233.        x1=x0+rx
  234.        y1=y0+ry
  235.  
  236. C      premier triangle
  237.        do k=1,ndec
  238.        do i=1,ndec-k+1
  239.        ip=ip+1
  240.        xx(1,ip) = x0+px*i
  241.        xx(2,ip) = y0+py*k
  242.        poids(ip)= ppp
  243.        if(i.eq.ndec-k+1) goto 57
  244.        ip=ip+1
  245.        xx(1,ip) = x1+px*i
  246.        xx(2,ip) = y1+py*k
  247.        poids(ip)= ppp
  248. 57    continue
  249.        enddo
  250.        x0 = x0 + px*cos(alpa)
  251.        enddo
  252.        ipo=ip
  253. C    tes triangles complementaires
  254.        do ik=1,nno-1
  255.        alpha= O2*xpi/nno*ik
  256.        cc = cos(alpha)
  257.        ss = sin(alpha)
  258.                 do j=1,ipo
  259.                 ip =ip+1
  260.                 xx(1,ip)=cc*xx(1,j)-ss*xx(2,j)
  261.                 xx(2,ip)=ss*xx(1,j)+cc*xx(2,j)
  262.                 poids(ip)=ppp
  263.                 enddo
  264.        enddo
  265. C       write(6,3333) ((xx(i,j),i=1,2),j=1,ip)
  266. 3333   format(2E12.5)
  267. C
  268.        goto 100
  269. 99     CONTINUE
  270.        MOTERR(1:4)=NOMTP(IELE)
  271.        MOTERR(9:14)='PAVPOI'
  272.        CALL ERREUR(86)
  273. C      ELEMENTS  NON IMPLEMENTES
  274.  
  275. 100    continue
  276.        ndd =ip
  277. C       write(6,*) 'pavpoi ip poids  som poids ' ,ip,ppp,(ip*ppp)
  278.        return
  279.  
  280.  
  281.        END
  282.  
  283.  
  284.  
  285.  
  286.  
  287.  
  288.  
  289.  
  290.  
  291.  
  292.  

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