Télécharger numopt.eso

Retour à la liste

Numérotation des lignes :

  1. C NUMOPT SOURCE PV 19/04/28 14:42:59 10210
  2. C RACINE DE LA NUMEROTATION POUR LA SORTIE SUR FAC
  3. C
  4. SUBROUTINE NUMOPT(MELEME,ICPR,NODES)
  5. IMPLICIT INTEGER(I-N)
  6. -INC SMELEME
  7. -INC SMCOORD
  8. -INC CCOPTIO
  9. -INC CCPRECO
  10. CHARACTER*4 MVAL
  11. DATA MVAL/'NOOP'/
  12. SEGMENT MEMJT1(NKON)
  13. SEGMENT MEMJT2(NKON)
  14. SEGMENT IPOME(NODES+1)
  15. SEGMENT JNT(NODES)
  16. SEGMENT JMEM(NODES)
  17. SEGMENT ICPR(XCOOR(/1)/(IDIM+1))
  18. SEGMENT ICPREN(XCOOR(/1)/(IDIM+1))
  19. SEGMENT IDCP(XCOOR(/1)/(IDIM+1))
  20. SEGMENT JOINT(NODES)
  21. SEGMENT NEWJT(NODES)
  22. icpren=0
  23. segact mcoord
  24. segact meleme
  25. * cas trivial ou le maillage est fait de points isoles
  26. CALL LIRMOT(MVAL,1,IRET,0)
  27. if (itypel.eq.1.or.iret.eq.1.or.nucrou.eq.-1) then
  28. segact icpr*MOD
  29. ia=0
  30. ipt1=meleme
  31. segact meleme
  32. do k=1,max(1,lisous(/1))
  33. if (lisous(/1).ne.0) ipt1=lisous(k)
  34. segact ipt1
  35. do j= ipt1.num(/2),1,-1
  36. do i=1,ipt1.num(/1)
  37. if (icpr(ipt1.num(i,j)).eq.0) then
  38. ia=ia+1
  39. icpr(ipt1.num(i,j))=ia
  40. endif
  41. enddo
  42. enddo
  43. enddo
  44. nodes=ia
  45. return
  46. endif
  47. nucro=0
  48. * si plus de 4000 noeuds dans le meleme numérotation NS
  49. segact icpr*mod
  50. do i=1,icpr(/1)
  51. icpr(i)=0
  52. enddo
  53. segact meleme
  54. ia=0
  55. ib=0
  56. ipd=0
  57. ipt1=meleme
  58. ifrot=0
  59. do 60 isous=1,max(1,lisous(/1))
  60. if (lisous(/1).ne.0) ipt1=lisous(isous)
  61. segact ipt1
  62. if (ipt1.itypel.eq.-22) then
  63. * write (6,*) ' frottement '
  64. ifrot=1
  65. endif
  66. do 50 i=1,ipt1.num(/2)
  67. do 50 j=1,ipt1.num(/1)
  68. if (icpr(ipt1.num(j,i)).eq.0) then
  69. ia=ia+1
  70. icpr(ipt1.num(j,i))=ia
  71. endif
  72. 50 continue
  73. if (ipt1.num(/1).gt.ib) then
  74. ib=ipt1.num(/1)
  75. ipd=icpr(ipt1.num(1,1))
  76. endif
  77. 60 continue
  78. if (ia.lt.ib*2 .and.ia.gt.100) then
  79. ** write (6,*) ' superelement detecte ',ia,ib
  80. ia=1
  81. else
  82. ipd=0
  83. endif
  84. do i=1,icpr(/1)
  85. icpr(i)=0
  86. enddo
  87. * choix de la numerotation en fonction du nombre de noeuds
  88. if (ia.gt.4000) nucro=2
  89. if (ia.le.4000) nucro=0
  90. * Reverse CM en iteratif
  91. IF(NUCROU.EQ.1) nucro=0
  92. IF(NUCROU.EQ.2) nucro=2
  93. * en cas de frottement forcément nested dissection pour placer correctement les mult de frot
  94. if (ifrot.eq.1) nucro=2
  95. * ajuster le nombre de zero consecutif
  96. if (nucro.eq.0) izrosf=48
  97. if (nucro.eq.2) izrosf=16
  98. if(ipd.ne.0) izrosf=16
  99. CALL LIRMOT(MVAL,1,IRET,0)
  100. if (iret.eq.0) then
  101. * a t'on deja fait l'operation ???
  102. segini icpren,idcp
  103. incren=0
  104. do isous=2,lisous(/1)
  105. ipt3=lisous(isous)
  106. segact ipt3
  107. * SG 2011/10/07 on a échangé l'ordre des boucles pour le remettre dans
  108. * l'ordre "usuel"
  109. do j=1,ipt3.num(/2)
  110. do i=1,ipt3.num(/1)
  111. ip=ipt3.num(i,j)
  112. if (icpren(ip).eq.0) then
  113. incren=incren+1
  114. icpren(ip)=incren
  115. endif
  116. enddo
  117. enddo
  118. enddo
  119. do i=1,icpren(/1)
  120. if (icpren(i).ne.0) idcp(icpren(i))=i
  121. enddo
  122.  
  123. do 150 ip=1,30
  124. if (prenum(ip).ne.0) then
  125. ipt1=prenum(ip)
  126. segact ipt1
  127. * write (6,*) ' prenum ',ip,lisous(/1),ipt1.lisous(/1)
  128. if (lisous(/1)+1.eq.ipt1.lisous(/1)) then
  129. *SG Attention ! On a osé stocker le type de renumérotation qui avait été
  130. *SG calculée dans le itypel
  131. ipt2=ipt1.lisous(ipt1.lisous(/1))
  132. segact ipt2
  133. nucroa=ipt2.itypel
  134. segdes ipt2
  135. if (nucro.ne.nucroa) goto 11
  136. * verif maillage identique
  137. * write (6,*) ' maillage 1 ',lisous(/1)
  138. * write (6,*) ' maillage 2 ',ipt1.lisous(/1)
  139. do isous=2,lisous(/1)
  140. ipt3=lisous(isous)
  141. ipt4=ipt1.lisous(isous)
  142. segact ipt3,ipt4
  143. * write (6,*) ipt3.num(/1),ipt4.num(/1),ipt3.num(/2),ipt4.num(/2)
  144. * write (6,*) ' itypel ',ipt3.itypel,ipt4.itypel
  145. if (ipt3.itypel.ne.ipt4.itypel) goto 10
  146. if (ipt3.num(/1).ne.ipt4.num(/1)) goto 10
  147. if (ipt3.num(/2).ne.ipt4.num(/2)) goto 10
  148. * SG 2011/10/07 on a échangé l'ordre des boucles pour le remettre dans
  149. * l'ordre "usuel"
  150. do j=1,ipt3.num(/2)
  151. do i=1,ipt3.num(/1)
  152. * write (6,*) i,j,ipt3.num(i,j),
  153. * > icpren(ipt3.num(i,j)),ipt4.num(i,j)
  154. if (icpren(ipt3.num(i,j)).ne.ipt4.num(i,j)) goto 10
  155. enddo
  156. enddo
  157. segdes ipt3,ipt4
  158. enddo
  159. if (iimpi.ne.0) write (6,*)
  160. > ' preconditionnement de la numerotation ',ip,
  161. > (prenum(i),i=1,30)
  162. ipt2=ipt1.lisous(ipt1.lisous(/1))
  163. segact ipt2
  164. segact icpr*MOD
  165. nodes=ipt2.num(/2)
  166. * write (6,*) ' icpr(/1),nodes ',icpr(/1),nodes
  167. do i=1,nodes
  168. if (idcp(ipt2.num(1,i)).ne.0) icpr(idcp(ipt2.num(1,i)))=i
  169. enddo
  170. ia=nodes
  171. segdes ipt1,ipt2
  172. * write (6,*) ' apres precond dans numopt '
  173. * write (6,*) (icpr(i),i=1,icpr(/1))
  174. return
  175.  
  176. 10 continue
  177. segdes ipt3,ipt4
  178. 11 continue
  179. endif
  180. endif
  181. 150 continue
  182. if (nucro.eq.1) then
  183. call numop1(MELEME,ICPR,NODES)
  184. goto 300
  185. endif
  186. if (nucro.eq.2) then
  187. call numop2(MELEME,ICPR,NODES)
  188. segact icpr*mod
  189. goto 300
  190. endif
  191. endif
  192. if (iret.eq.1) call refus
  193. IENORM=2000000000
  194. NODES=XCOOR(/1)/(IDIM+1)
  195. SEGACT ICPR*MOD
  196. SEGACT MELEME
  197. DO 200 I=1,ICPR(/1)
  198. 200 ICPR(I)=0
  199. IPT1=MELEME
  200. IKOU=0
  201. DO 202 IO=1,MAX(1,LISOUS(/1))
  202. IF (LISOUS(/1).NE.0) THEN
  203. IPT1=LISOUS(IO)
  204. SEGACT IPT1
  205. ENDIF
  206. DO 203 J=1,IPT1.NUM(/2)
  207. * inversion car on reinerse a la fin
  208. DO 203 I=1,IPT1.NUM(/1)
  209. IJ=IPT1.NUM(I,J)
  210. IF (ICPR(IJ).NE.0) GOTO 203
  211. IKOU=IKOU+1
  212. ICPR(IJ)=IKOU
  213. 203 CONTINUE
  214. 202 CONTINUE
  215. NODES=IKOU
  216. SEGINI JNT,JMEM,NEWJT,JOINT,IPOME
  217. DO 1 I=1,NODES
  218. 1 JMEM(I)=0
  219. IPT1=MELEME
  220. DO 3 IO=1,MAX(1,LISOUS(/1))
  221. IF (LISOUS(/1).NE.0) IPT1=LISOUS(IO)
  222. DO 4 J=1,IPT1.NUM(/2)
  223. DO 4 I=1,IPT1.NUM(/1)
  224. JMEM(ICPR(IPT1.NUM(I,J)))=JMEM(ICPR(IPT1.NUM(I,J)))+1
  225. 4 CONTINUE
  226. 3 CONTINUE
  227. IPOME(1)=0
  228. DO 6 I=1,NODES
  229. IPOME(I+1)=IPOME (I) + JMEM(I)
  230. 6 CONTINUE
  231. NKON=IPOME(NODES+1)
  232. SEGINI MEMJT1,memjt2
  233. DO 90 J=1,NODES
  234. 90 JMEM(J)=0
  235. IPT1=MELEME
  236. DO 101 IO=1,MAX(1,LISOUS(/1))
  237. IF (LISOUS(/1).NE.0) IPT1=LISOUS(IO)
  238. DO 100 I=1,IPT1.NUM(/2)
  239. DO 100 J=1,IPT1.NUM(/1)
  240. IND=ICPR(IPT1.NUM(J,I))
  241. JMEM(IND)=JMEM(IND)+1
  242. MEMJT1(IPOME(IND)+JMEM(IND))=I
  243. MEMJT2(IPOME(IND)+JMEM(IND))=IO
  244. 100 CONTINUE
  245. 101 CONTINUE
  246. IDIFF=1
  247. CALL OPTNUM(MELEME,MEMJT1,memjt2,JMEM,JNT,NEWJT,JOINT,IDIFF,IPOME,
  248. # ICPR,NODES,NOOPTI,ipd)
  249. * ne pas inverse si superelement pour conserver sa numerotation
  250. if (ipd.eq.0) then
  251. DO 110 I=1,NODES
  252. JNT(I)=NODES-JNT(I)+1
  253. 110 CONTINUE
  254. endif
  255. C PERMUTER LES COORDONNEES ET CORRIGER IPAD
  256. DO 111 I=1,icpr(/1)
  257. IF (ICPR(I).EQ.0) GOTO 111
  258. ICPR(I)=JNT(ICPR(I))
  259. 111 CONTINUE
  260. IF(NOOPTI .EQ.0) GO TO 116
  261. IA=0
  262. DO 115 I=1,icpr(/1)
  263. IF(ICPR(I).EQ.0) GO TO 115
  264. IA=IA+1
  265. ICPR(I)=IA
  266. 115 CONTINUE
  267. 116 CONTINUE
  268. SEGSUP MEMJT1,MEMJT2,JMEM,NEWJT,JOINT,IPOME
  269. SEGSUP JNT
  270. * Mise en commentaire du retuen par SG car comme cela on préconditionne
  271. * aussi le Cuthill-McKee
  272. * RETURN
  273. 300 continue
  274. * sauver le resultat au cas ou on veuille recommencer
  275. segini,ipt1=meleme
  276. nbsous=ipt1.lisous(/1)+1
  277. nbref=ipt1.lisref(/1)
  278. nbnn=ipt1.num(/1)
  279. nbelem=ipt1.num(/2)
  280. segadj ipt1
  281. nbnn=1
  282. nbsous=0
  283. nbref=0
  284. nbelem=nodes
  285. segini ipt2
  286. *SG Attention ! On ose stocker le type de renumérotation qui a été
  287. *SG calculée dans le itypel
  288. ipt2.itypel=nucro
  289. ia=0
  290. do 120 i=1,icpr(/1)
  291. if (icpr(i).gt.nodes.or.icpr(i).eq.0) goto 120
  292. ia=ia+1
  293. ipt2.num(1,icpr(i))=icpren(i)
  294. 120 continue
  295. * call ecmail(ipt2)
  296. ipt1.lisous(ipt1.lisous(/1))=ipt2
  297. segdes ipt2
  298. do 125 is=1,ipt1.lisous(/1)-1
  299. ipt2=ipt1.lisous(is)
  300. segini,ipt3=ipt2
  301. segdes ipt2
  302. * SG 2011/10/07 on a échangé l'ordre des boucles pour le remettre dans
  303. * l'ordre "usuel"
  304. do il=1,ipt3.num(/2)
  305. do ip=1,ipt3.num(/1)
  306. ipt3.num(ip,il)=icpren(ipt3.num(ip,il))
  307. enddo
  308. enddo
  309. segdes ipt3
  310. ipt1.lisous(is)=ipt3
  311. 125 continue
  312. if (ipt1.lisous(/1).eq.1) goto 131
  313. segdes ipt1
  314. do 130 ip=29,1,-1
  315. prenum(ip+1)=prenum(ip)
  316. 130 continue
  317. prenum(1)=ipt1
  318. 131 continue
  319. if (icpren.ne.0) segsup icpren,idcp
  320. segact icpr*mod
  321. * write (6,*) ' fin de numopt '
  322. * write (6,*) (icpr(i),i=1,icpr(/1))
  323. *
  324. RETURN
  325. END
  326.  
  327.  
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  
  334.  
  335.  
  336.  
  337.  
  338.  
  339.  
  340.  
  341.  
  342.  
  343.  
  344.  

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