Télécharger numop2.eso

Retour à la liste

Numérotation des lignes :

  1. C NUMOP2 SOURCE PV 19/04/10 21:15:02 10190
  2. C RACINE DE LA NUMEROTATION POUR LA SORTIE SUR FAC
  3. C
  4. C methode utilisee: NESTED DISSECTION.
  5.  
  6. C ce programme refait l'indicage de la matrice afin de minimiser
  7. C le profil.
  8.  
  9.  
  10. SUBROUTINE NUMOP2(MELEME,ICPR,NODES)
  11. IMPLICIT INTEGER(I-N)
  12. -INC SMELEME
  13. -INC SMCOORD
  14. -INC CCOPTIO
  15. -INC CCASSIS
  16.  
  17. SEGMENT JMEM(NODES+1),JMEMN(NODES+1)
  18. C JMEM et JMEMN contiennent le nombre d'element auquel appartient un noeud
  19.  
  20. SEGMENT JNT(NODES)
  21. C JNT contient la nouvelle numerotation
  22.  
  23. SEGMENT ICPR(XCOOR(/1)/(IDIM+1))
  24. C ICPR au debut contient l'ancienne numerotation ,
  25. C a la fin la nouvelle.
  26.  
  27. SEGMENT IADJ(NODES+1)
  28. SEGMENT JADJC(0)
  29. C IADJ(i) pointe sur JADJC qui contient les voisins de i entre
  30. C IADJ(i) et IADJ(i+1)-1
  31.  
  32. SEGMENT LAGRAN(NB)
  33. C contient les noeud de lagrange et les noeuds les suivant directement
  34. C cf element de type 22
  35. SEGMENT INVLAG(NB1)
  36.  
  37. SEGMENT BOOLEEN
  38. LOGICAL BOOL(NODES)
  39. ENDSEGMENT
  40. C BOOL(i) = true si le noeud i a ete deja mentionne dans la liste
  41. C des voisins JADJC.
  42.  
  43. SEGMENT IMEMOIR(NBV),LMEMOIR(NBV)
  44. C contient les elements appartenant a chaque noeud,sous forme de liste.
  45.  
  46. INTEGER ELEM
  47. C nom d'un element
  48.  
  49. INTEGER N
  50.  
  51.  
  52. SEGMENT MASQUE
  53. LOGICAL MASQ(NODES)
  54. ENDSEGMENT
  55. C MASQ(X)=.TRUE. si le noeud X n'a pas ete renumerote;
  56. C .FALSE. si il l'a ete.
  57.  
  58. INTEGER DIM,DIMSEP
  59. C DIM= nombre de noeuds renumerotes.
  60.  
  61. INTEGER PIVOT
  62. C PIVOT est le noeud utile a la division du domaine.
  63.  
  64. SEGMENT IPOS(NODES*3)
  65. C est le vecteur contenant la numerotation dans les deux sens,de 1 a NODES
  66. C puis de NODES+1 a 2* NODES, cf la subroutine SEPAR
  67. C
  68. C segments utilisés dans sepa2
  69. C
  70. SEGMENT NRELONG(NODES*nbthr)
  71. C NRELONG contient pour chaque noeud sa profondeur.
  72.  
  73. SEGMENT NOELON(NODES*nbthr)
  74. SEGMENT NOEL2(NODES)
  75. SEGMENT LONDIM(NODES*nbthr)
  76. C NOELON contient les noeuds de profondeur LONG.
  77. C DIMLON= dimension de NOELON.
  78.  
  79.  
  80. C**********************************
  81.  
  82. C debut du program
  83.  
  84. C**********************************
  85.  
  86.  
  87. * pour que izero soit de la bonne longueur
  88. izero=0
  89. C initialisation
  90. C*******************************
  91.  
  92.  
  93. C norme d'erreur
  94. SEGACT ICPR*MOD
  95. NODES=ICPR(/1)
  96. SEGACT MELEME
  97. C icpr: numero des noeuds.
  98. C meleme: objet de maillage (cf assem2.eso)
  99.  
  100. DO 10 I=1,ICPR(/1)
  101. ICPR(I)=0
  102. 10 CONTINUE
  103.  
  104. IPT1=MELEME
  105. IKOU=0
  106. NBV=0
  107. NB1=0
  108. NB2=0
  109.  
  110. DO 100 IO=1,MAX(1,LISOUS(/1))
  111. IF (LISOUS(/1).GT.0) THEN
  112. IPT1=LISOUS(IO)
  113. SEGACT IPT1
  114. ENDIF
  115. C on cree la numerotation des noeuds.
  116. C 'nb noeuds/element'=IPT1.NUM(/1)
  117. C 'nb element'=IPT1.NUM(/2)
  118. IF(abs(IPT1.ITYPEL).EQ.22) THEN
  119. NB1=NB1+IPT1.NUM(/2)
  120. NB2=MAX(NB2,IPT1.NUM(/1))
  121. C NB1= nbre d'éléments de type 22.
  122. C NB2=nbre de noeuds/élément maximum parmi
  123. C les éléments de type 22.
  124. ENDIF
  125. DO 150 J=1,IPT1.NUM(/2)
  126. DO 150 I=1,IPT1.NUM(/1)
  127. IJ=IPT1.NUM(I,J)
  128. C IJ est le Ième noeud du Jème élément.
  129. IF (ICPR(IJ).EQ.0) THEN
  130. C s'il est déjà numéroté, on ne fait rien.
  131. IKOU=IKOU+1
  132. ICPR(IJ)=IKOU
  133. ENDIF
  134. 150 CONTINUE
  135. 100 CONTINUE
  136.  
  137. NODES=IKOU
  138. NB=NB2*NB1
  139.  
  140. C***** initalisation des segments*********
  141.  
  142. SEGINI IADJ,JADJC,JMEM,JMEMN,LAGRAN
  143. SEGINI BOOLEEN,JNT
  144.  
  145. DO 20 I=1,NODES+1
  146. IADJ(I)=0
  147. JMEM(I)=0
  148. JMEMN(I)=0
  149. 20 CONTINUE
  150.  
  151. C******************************************
  152.  
  153. IPT1=MELEME
  154. IADJ(1)=1
  155. INC=0
  156. DO 200 IO=1,MAX(1,LISOUS(/1))
  157.  
  158. IF (LISOUS(/1).NE.0) THEN
  159. IPT1=LISOUS(IO)
  160. ENDIF
  161.  
  162. DO 210 J=1,IPT1.NUM(/2)
  163. IF(abs(IPT1.ITYPEL).EQ.22) THEN
  164. is=sign(1,ipt1.itypel)
  165. DO 220 I=1,IPT1.NUM(/1)
  166. C chaque element de type 22 a au plus NB2 noeuds.
  167. LAGRAN(INC*NB2+I)=ICPR(IPT1.NUM(I,J))*is
  168. C les noeuds de l'elements de type 22
  169. C sont ranges dans le vecteur LAGRAN.
  170. 220 CONTINUE
  171. DO 225 I=IPT1.NUM(/1)+1,NB2
  172. LAGRAN(INC*NB2+I)=0
  173. C comme on a alloue la place memoire maximale,
  174. C on remplit les cases restantes avec des 0.
  175. 225 CONTINUE
  176. INC=INC+1
  177. C INC=le nbre d'elements de type 22.
  178. ENDIF
  179.  
  180. DO 230 I=1,IPT1.NUM(/1)
  181. IJ=ICPR(IPT1.NUM(I,J))+1
  182. JMEM(IJ)=JMEM(IJ)+1
  183. C JMEM(I+1): nb elements auquel le noeud I appartient
  184. 230 CONTINUE
  185. 210 CONTINUE
  186.  
  187. 200 CONTINUE
  188.  
  189.  
  190. JMEM(1)=1
  191. DO 30 I=1,NODES
  192. JMEM(I+1)=JMEM(I)+JMEM(I+1)
  193. C JMEM(I+1)=indice de depart des elements
  194. C auxquels le noeud I appartient.
  195. 30 CONTINUE
  196. NBV=JMEM(NODES+1)
  197. C NBV= dimension de IMEMOIR.
  198. SEGINI IMEMOIR,LMEMOIR
  199.  
  200.  
  201.  
  202. IPT1=MELEME
  203.  
  204. DO 300 IO=MAX(1,LISOUS(/1)),1,-1
  205. IF (LISOUS(/1).NE.0) THEN
  206. IPT1=LISOUS(IO)
  207. ENDIF
  208. DO 350 J=1,IPT1.NUM(/2)
  209. DO 350 I=1,IPT1.NUM(/1)
  210. IJ=ICPR(IPT1.NUM(I,J))
  211. JMEMN(IJ+1)=JMEMN(IJ+1)+1
  212. IMEMOIR(JMEM(IJ)+JMEMN(IJ+1)-1)=J
  213. LMEMOIR(JMEM(IJ)+JMEMN(IJ+1)-1)=IO
  214. C on range dans IMEMOIR tous les elements des sous-objets
  215. C IO auxquels appartient le noeud ICPR(IPT1.NUM(I,J)).
  216. C On connait pour chaque element, le sous-objet auquel
  217. C il appartient grace a LMEMOIR
  218. 350 CONTINUE
  219. 300 CONTINUE
  220.  
  221. DO 410 J=1,NODES
  222. BOOL(J)=.FALSE.
  223. 410 CONTINUE
  224. DO 400 I=1,NODES
  225. IADJ(I+1)=IADJ(I)
  226. DO 420 J=JMEM(I),JMEM(I+1)-1
  227. ELEM=IMEMOIR(J)
  228. C ELEM=element auquel appartient le noeud I.
  229.  
  230. IPT1=MELEME
  231. IF (LISOUS(/1).NE.0) IPT1=LISOUS(LMEMOIR(J))
  232. DO 430 K=1,IPT1.NUM(/1)
  233. C k representatif du nb de noeuds par elements.
  234. IK=ICPR(IPT1.NUM(K,ELEM))
  235. IF ((I.NE.IK).AND.
  236. & (.NOT.(BOOL(IK)))) THEN
  237. C si i n'est pas egal a un des nouveaux numeros des noeuds
  238. C de l'element ELEM et si il n'appartient pas deja a l'ens des
  239. C voisins du noeud i(jadjc(i)),alors on le rajoute.
  240. C JADJC(IADJ(I+1))=IK
  241. JADJC(**)=IK
  242. IADJ(I+1)=IADJ(I+1)+1
  243. BOOL(IK)=.TRUE.
  244. ENDIF
  245. 430 CONTINUE
  246. 420 CONTINUE
  247. * remise a faux de bool
  248. DO 412 J=IADJ(I),IADJ(I+1)-1
  249. IK=JADJC(J)
  250. BOOL(IK)=.FALSE.
  251. 412 CONTINUE
  252. 400 CONTINUE
  253.  
  254. SEGSUP JMEM,JMEMN,IMEMOIR,LMEMOIR,BOOLEEN
  255.  
  256.  
  257.  
  258. C**************************************************************************
  259.  
  260.  
  261. C affectation
  262. C************************
  263.  
  264.  
  265. if (nbthrs.gt.1) call threadii
  266. SEGINI IPOS,MASQUE
  267. IPOSMAX=0
  268. ** write (6,*) ' nodes ',nodes
  269. DO 50 I=1,NODES
  270. MASQ(I)=.TRUE.
  271. IPOS(I)=0
  272. IPOS(NODES+I)=0
  273. IPOS(2*NODES+I)=0
  274. 50 CONTINUE
  275. C initialement, les noeuds ne sont pas masques,ont donc
  276. C une position nulle.
  277.  
  278. DIM=0
  279. C le nombre de noeuds renumerotes DIM est initialement egal a zero.
  280. C on initialise un premier separateur avec lisous(1) et lisous(2)
  281. C qui contiennent les noeuds maitres du super element (si appele par assem4)
  282. mdomn=0
  283. iposv=ipos(mdomn+1)
  284. iposmax=iposmax+3
  285. ipos(iposmax+1-2)=iposv+1
  286. ipos(iposmax+1-1)=iposv+2
  287. ipos(iposmax+1-0)=iposv+3
  288. dimsep=0
  289. do io=1,2
  290. ipt1=lisous(io)
  291. if (ipt1.itypel.eq.0 ) then
  292. do j=1,ipt1.num(/2)
  293. do i=1,ipt1.num(/1)
  294. ip=icpr(ipt1.num(i,j))
  295. if (masq(ip)) then
  296. masq(ip)=.false.
  297. ipos(ip+nodes)=iposmax-2
  298. dimsep=dimsep+1
  299. endif
  300. enddo
  301. enddo
  302. endif
  303. enddo
  304. do i=1,nodes
  305. if (masq(i)) then
  306. ipos(i+nodes)=iposmax
  307. ipos(i+2*nodes)=ipos(i+2*nodes)+1
  308. endif
  309. enddo
  310. ** write(6,*) ' dimsep dans numop2 ',dimsep
  311. dim=dim+dimsep
  312. C ****************************************
  313. C boucle principale
  314. NS=NODES
  315. ns=ns-dimsep
  316.  
  317. nbthr=min(64,nbthrs)
  318. SEGINI NRELONG,NOELON,noel2,londim
  319. ** write (6,*) ' avant appel sepa2 '
  320. DO 500 I=1,NODES
  321. 550 IF(.NOT.MASQ(I)) GOTO 500
  322. C si le noeud est masque alors ne rien faire: il est deja
  323. C renumerote. On passe au noeud suivant.
  324.  
  325. PIVOT=I
  326.  
  327. CALL SEPA2(IADJ,JADJC,PIVOT,MASQUE,DIMSEP,NS,
  328. > IPOS,NODES,IPOSMAX,nrelong,noelon,noel2,
  329. > londim,nbthr,izero)
  330. if (ierr.ne.0) return
  331. C separe le domaine d'etude en 2 parties.
  332. C on decrit le domaine d'etude a partir du pivot et on cherche la
  333. C longueur maximale en decrivant les voisins de pivot, et leurs
  334. C voisins... jusqu'a rencontrer un voisin masque. On cree alors
  335. C une nouvelle separation.
  336. C les noeuds masques delimitent la separation du domaine.
  337.  
  338.  
  339. DIM=DIM+DIMSEP
  340. NS=NS-DIMSEP
  341. C la dimension de noeuds renumerotes est augmente de DIMSEP.
  342. C Celle de noeuds a renumeroter est diminue de DIMSEP.
  343.  
  344. * IF (DIM.GE.NODES) GOTO 600
  345. C si tous les noeuds ont ete renumerotes, on arrete.
  346.  
  347. GOTO 550
  348.  
  349. 500 CONTINUE
  350. ** write (6,*) ' apres appel sepa2 '
  351.  
  352. SEGSUP NRELONG,NOELON,noel2,londim
  353.  
  354. 600 CONTINUE
  355. if (nbthrs.gt.1) call threadis
  356. *
  357. * tri dans chaque zone
  358. do 610 lpoint=1,nodes
  359. mdomn=ipos(lpoint+nodes)
  360. iposv=ipos(mdomn+1)
  361. iposi=3*nodes
  362. do 620 kk=iadj(lpoint),iadj(lpoint+1)-1
  363. k=jadjc(kk)
  364. iposk=ipos(ipos(k+nodes)+1)
  365. if (iposk.ne.iposv) then
  366. iposi=min(iposi,iposk)
  367. endif
  368. 620 continue
  369. ipos(lpoint+2*nodes)=iposi*(2*(nb1+2))
  370. 610 continue
  371.  
  372. *
  373. * mise a la bonne place des multiplicateurs de Lagrange
  374. do 700 J=0,NB1-1
  375. jp1=j+1
  376. iposvs=2**30
  377. iposvr=-2**30
  378. * ips est le premier noeud concerne par la relation, ipr le dernier
  379. ipr=1
  380. ips=1
  381. mdomnr=0
  382. mdomns=0
  383. if (nb2.eq.2) goto 710
  384. * write (6,*) 'numop2 ',(lagran(J*NB2+il),il=1,nb2)
  385. do 800 il=3,nb2
  386. ip = abs(LAGRAN(J*NB2+il))
  387. * if (ip.eq.0) write (6,*) ' prob numop2 '
  388. if (ip.eq.0) goto 800
  389. mdomn=ipos(ip+nodes)
  390. iposv=ipos(mdomn+1)
  391. if (iposv.gt.iposvr) then
  392. ipr=ip
  393. iposvr=iposv
  394. mdomnr=mdomn
  395. ipaur=ipos(ip+2*nodes)
  396. elseif (iposv.eq.iposvr) then
  397. if (ipos(ip+2*nodes).gt.ipaur) then
  398. ipaur=ipos(ip+2*nodes)
  399. ipr=ip
  400. endif
  401. endif
  402. if (iposv.lt.iposvs) then
  403. ips=ip
  404. iposvs=iposv
  405. mdomns=mdomn
  406. ipaus=ipos(ip+2*nodes)
  407. elseif (iposv.eq.iposvs) then
  408. if (ipos(ip+2*nodes).lt.ipaus) then
  409. ipaus=ipos(ip+2*nodes)
  410. ips=ip
  411. endif
  412. endif
  413. 800 continue
  414. 710 continue
  415. *pv if (ipr.eq.0) write (6,*) ' probleme ipr numop2'
  416. *pv if (ips.eq.0) write (6,*) ' probleme ips numop2'
  417. *
  418. * le premier mult avant le premier noeud
  419. ip=abs(LAGRAN(J*NB2+1))
  420. IPOS(IP+2*NODES)= ipaus-jp1-1
  421. * numeroter le frottement apres le contact
  422. if (LAGRAN(J*NB2+1).lt.0) IPOS(IP+2*NODES)=ipaus-jp1
  423. IPOS(IP+nodes)=mdomns
  424. ** write (6,*) 'premier mult ',ip,ipos(ip+nodes),ipos(ip+2*nodes)
  425. *
  426. * le deuxieme mult apres le dernier noeud
  427. ip=abs(LAGRAN(J*NB2+2))
  428. IPOS(IP+2*NODES)= ipaur+jp1+1
  429. * numeroter le frottement apres le contact
  430. if (LAGRAN(J*NB2+2).lt.0) IPOS(IP+2*NODES)=ipaur+jp1
  431. IPOS(IP+nodes)=mdomnr
  432. ** write (6,*) 'deuxieme mult ',ip,ipos(ip+nodes),ipos(ip+2*nodes)
  433. *
  434. 700 continue
  435. SEGSUP IADJ,JADJC,MASQUE
  436. * ok maintenant on trie
  437. CALL SORTI2(IPOS,JNT,NODES)
  438.  
  439. C***************************************************************************
  440.  
  441. DO 860 I=1,ICPR(/1)
  442. IF(ICPR(I).NE.0) ICPR(I)=JNT(ICPR(I))
  443. C numerotation finale.
  444. 860 CONTINUE
  445.  
  446.  
  447. SEGSUP JNT,IPOS,LAGRAN
  448.  
  449.  
  450. RETURN
  451. END
  452.  
  453.  
  454.  
  455.  
  456.  
  457.  
  458.  
  459.  
  460.  
  461.  
  462.  
  463.  
  464.  
  465.  
  466.  
  467.  
  468.  
  469.  
  470.  
  471.  
  472.  
  473.  
  474.  
  475.  
  476.  
  477.  
  478.  
  479.  
  480.  

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