Télécharger pavpoi.eso

Retour à la liste

Numérotation des lignes :

  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. -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'
  270. CALL ERREUR(86)
  271. C ELEMENTS NON IMPLEMENTES
  272.  
  273. 100 continue
  274. ndd =ip
  275. C write(6,*) 'pavpoi ip poids som poids ' ,ip,ppp,(ip*ppp)
  276. return
  277.  
  278.  
  279. END
  280.  
  281.  
  282.  
  283.  
  284.  
  285.  
  286.  
  287.  
  288.  
  289.  
  290.  

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