Esope source files

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