Télécharger numop2.eso

Retour à la liste

Numérotation des lignes :

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

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