Télécharger raff.eso

Retour à la liste

Numérotation des lignes :

raff
  1. C RAFF SOURCE MB234859 25/09/08 21:16:04 12358
  2. SUBROUTINE RAFF
  3. C----------------------------------------------------------------------C
  4. C Operateur RAFF C
  5. C----------------------------------------------------------------------C
  6. C
  7. C Syntaxe 1 : GEO2 = RAFF GEO1 CHPO1 ;
  8. C
  9. C Syntaxe 2 : LRE2 = RAFF LRE1 ENT1 ;
  10. C
  11. C----------------------------------------------------------------------C
  12. IMPLICIT REAL*8(A-H,O-Z)
  13. IMPLICIT INTEGER(I-N)
  14.  
  15. -INC PPARAM
  16. -INC CCOPTIO
  17. -INC CCGEOME
  18. c-INC CCREEL
  19.  
  20. -INC SMCOORD
  21. -INC SMELEME
  22. -INC SMCHPOI
  23. -INC SMMODEL
  24. -INC SMCHAML
  25.  
  26. INTEGER MCOL,NCOL
  27. C--------------------------------
  28. C B) VARIABLES RELATIVES A LA 2D !
  29. C--------------------------------
  30. SEGMENT KARETE(NBNDS,NCOL)
  31. SEGMENT KARET2(NBNDS,NCOL)
  32. SEGMENT KMILIE(NBNDS,NCOL)
  33. SEGMENT KARPOS(NBNDS)
  34. C
  35. C--------------------------------
  36. C C) VARIABLES RELATIVES A LA 3D !
  37. C--------------------------------
  38. SEGMENT IRELA(NBNDS)
  39. SEGMENT JPLANS(JPLA1,JPLA2)
  40. SEGMENT JPLAN2(JPLA1,JPLA2)
  41. SEGMENT JPLAN3(JPLA1,JPLA2)
  42. SEGMENT JPLCOM(JPLA1)
  43. SEGMENT JNOEFA(JNBFA,5)
  44. SEGMENT IWORK2(JNBSOM)
  45. SEGMENT JFARAF(JNBFA,LLLL)
  46. SEGMENT IWORK1(IJKLMN)
  47. C--------------------------
  48. C A) VARIABLES 'GENERALES' !
  49. C--------------------------
  50. SEGMENT ICPR(nbpts,2)
  51. SEGMENT LTYPE(NBLIS,2)
  52. SEGMENT LTYPE2(NBLIS2,4)
  53. SEGMENT LTYPE3(NBLIS3,2)
  54. SEGMENT LTYPE4(NBLIS3,2)
  55. SEGMENT IDENSI
  56. REAL*8 DENSI(NBNDS)
  57. ENDSEGMENT
  58. SEGMENT XDEN(nbpts)
  59. XPETI= 1E-30
  60. C
  61. C---------------------------------------------------
  62. C D) VARIABLES RELATIVES AUX MAILLAGES DE RELATIONS !
  63. C---------------------------------------------------
  64. POINTEUR IPT10.MELEME
  65. POINTEUR IPT11.MELEME
  66. POINTEUR IPT12.MELEME
  67. POINTEUR IPT13.MELEME
  68.  
  69. POINTEUR IPT14.MELEME
  70. POINTEUR IPT15.MELEME
  71. POINTEUR IPT16.MELEME
  72. SEGMENT MAIREL(NBMAIL)
  73. C=======================================================================
  74. C Lecture des donnees utilisateurs
  75. C=======================================================================
  76. C
  77. C SYNTAXE 2 : cas du LISTREEL
  78. C ---------------------------
  79. CALL LIROBJ('LISTREEL',ILRE1,0,IRETOU)
  80. IF (IERR.NE.0) RETURN
  81. C
  82. IF (IRETOU.NE.0) THEN
  83. ICAS = 0
  84. CALL LIRENT(KR1,0,IRETOU)
  85. IF (IERR.NE.0) RETURN
  86. IF (IRETOU.EQ.0) THEN
  87. CALL LIRREE(XR1,1,IRETOU)
  88. IF (IERR.NE.0) RETURN
  89. ICAS = 2
  90. ELSE
  91. ICAS = 1
  92. ENDIF
  93. CALL RAFLRE(ILRE1,KR1,XR1,ICAS,ILRE2)
  94. IF (IERR.NE.0) RETURN
  95. IF (ILRE2.NE.0) THEN
  96. CALL ECROBJ('LISTREEL',ILRE2)
  97. ELSE
  98. CALL ERREUR(5)
  99. ENDIF
  100. RETURN
  101. ENDIF
  102.  
  103. C SYNTAXE 1 : cas du MAILLAGE
  104. C ---------------------------
  105. C On lit le maillage a raffiner et le champoint de densite a affecter
  106. c ipt1 : maillage initial à raffiner
  107. c mchpo1 : chapoint de densité voulu(que devrai être définit sur ipt1)
  108. CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU)
  109. CALL ACTOBJ('MAILLAGE',IPT1,1)
  110. IF (IERR.NE.0) RETURN
  111. CALL LIROBJ('CHPOINT',MCHPO1,1,IRETOU)
  112. CALL ACTOBJ('CHPOINT',MCHPO1,1)
  113. IF (IERR.NE.0) RETURN
  114.  
  115. C=======================================================================
  116. C Initialisations - Declarations
  117. C=======================================================================
  118. NCOL=0
  119. MCOL=0
  120. NBOUCL=0
  121. NPTINI=nbpts
  122.  
  123. IJKLMN=4
  124. SEGINI IWORK1
  125. IWORK1(1)=4
  126. IWORK1(2)=8
  127. IWORK1(3)=6
  128. IWORK1(4)=10
  129.  
  130. SEGINI,XDEN
  131.  
  132. C=======================================================================
  133. C Tests sur les donnees
  134. C=======================================================================
  135. C On teste pour savoir si le maillage IPT1 (maillage 1) a un type
  136. C d'elements implemente.
  137. c TRI3 TRI6 QUA4 QUA8 CUB8 CUB26 PRY6 PR15 TET4 TET10 SURE
  138. C les pyramides 5 et 13 ont été débranchés
  139. c création du tableau
  140. NBSOU1 = IPT1.LISOUS(/1)
  141. NBLIS = MAX(1,NBSOU1)
  142. NSURE=0
  143. SEGINI LTYPE
  144. IPT2=IPT1
  145. DO 1 INB=1,NBLIS
  146. IF (NBSOU1.NE.0) IPT2=IPT1.LISOUS(INB)
  147. NT=IPT2.ITYPEL
  148. IF (NT.EQ.48) NSURE=NSURE+1
  149. LTYPE(INB,1)=INB
  150. LTYPE(INB,2)=NT
  151. IF (NT.EQ.16) THEN
  152. write (*,*) 'Attention le raffinement PRI6 a debrancher !'
  153. ENDIF
  154. IF (NT.EQ.23) THEN
  155. write (*,*) 'Attention le raffinement TET4 a debrancher !'
  156. ENDIF
  157. IF ((NT.NE.4).AND.(NT.NE.6).AND.(NT.NE.8).AND.(NT.NE.10).AND.
  158. + (NT.NE.14).AND.(NT.NE.15).AND.(NT.NE.16).AND.(NT.NE.17).AND.
  159. + (NT.NE.23).AND.(NT.NE.24).AND.(NT.NE.48)) THEN
  160. WRITE(*,*) 'Erreur, type d element non pris en compte'
  161. RETURN
  162. ENDIF
  163. c write(*,*) 'type d elem du maillage de depart', LTYPE(INB,1),
  164. c + LTYPE(INB,2)
  165.  
  166. C IF ((NT.NE.4).AND.(NT.NE.6).AND.(NT.NE.8).AND.(NT.NE.10).AND.
  167. C + (NT.NE.14).AND.(NT.NE.15).AND.(NT.NE.16).AND.(NT.NE.17).AND.
  168. C + (NT.NE.23).AND.(NT.NE.24).AND.(NT.NE.48).AND.(NT.NE.25).AND.
  169. C + (NT.NE.26)) RETURN
  170. IF ((NBSOU1.EQ.0).AND.(NT.EQ.48)) RETURN
  171. 1 CONTINUE
  172.  
  173. C
  174. C=======================================================================
  175. C Separation du maillage initial
  176. C=======================================================================
  177. C Si le maillage initial comporte des lisous de type 48, on le separe
  178. C en deux maillages dont un ne sera compose que des lisous de type 48.
  179. C Maillage d'elements de type 48 : IPT6
  180. C Maillage sans elements de type 48 : IPT5
  181. C IPT1 reste le maillage initial complet
  182. LCONF2=0
  183. IF (NSURE.EQ.0) GOTO 7
  184. LCONF2=1
  185. 7 CONTINUE
  186. C-------------------------------
  187. C A) Creation du maillage IPT6 !
  188. C-------------------------------
  189. C A.1/** Cas ou il n'y a qu'un lisous de type 48
  190. IF (NSURE.EQ.1) THEN
  191. DO 2 K=1,LTYPE(/1)
  192. IF (LTYPE(K,2).EQ.48) IPT6=IPT1.LISOUS(LTYPE(K,1))
  193. 2 CONTINUE
  194. ENDIF
  195. C
  196. C A.2/** Cas ou il y a plusieurs lisous de type 48
  197. IF (NSURE.GT.1) THEN
  198. NBREF=0
  199. NBSOUS=NSURE
  200. NBNN=0
  201. NBELEM=0
  202. SEGINI IPT6
  203. KLISOU=0
  204. DO 3 K=1,LTYPE(/1)
  205. IF (LTYPE(K,2).EQ.48) THEN
  206. KLISOU=KLISOU+1
  207. IPT6.LISOUS(KLISOU)=IPT1.LISOUS(LTYPE(K,1))
  208. ENDIF
  209. 3 CONTINUE
  210. ENDIF
  211. IF (NSURE.EQ.0) THEN
  212. NBREF=0
  213. NBSOUS=NSURE
  214. NBNN=0
  215. NBELEM=0
  216. SEGINI IPT6
  217. ENDIF
  218. C
  219. C------------------------------
  220. C B) Creation de IPT5 !
  221. C------------------------------
  222. C
  223. NBSOUS=LTYPE(/1)-NSURE
  224. NBREF=0
  225. NBNN=0
  226. NBELEM=0
  227. C B.0/** Cas ou il n'y a qu'un lisous en tout
  228. IF ((NBSOUS.EQ.1).AND.(NSURE.EQ.0)) THEN
  229. NBSOUS=0
  230. SEGINI IPT5
  231. IPT5=IPT1
  232. ENDIF
  233. C B.1/** Cas ou il n'y a qu'un lisous de type different de 48
  234. IF ((NBSOUS.EQ.1).AND.(NSURE.NE.0)) THEN
  235. SEGINI IPT5
  236. DO 4 K=1,LTYPE(/1)
  237. IF (LTYPE(K,2).NE.48) IPT5=IPT1.LISOUS(LTYPE(K,1))
  238. 4 CONTINUE
  239. ENDIF
  240. C
  241. C B.2/** Cas ou il y a plusieurs lisous de types differents de 48
  242. IF (NBSOUS.GT.1) THEN
  243. SEGINI IPT5
  244. KLISOU=0
  245. DO 5 K=1,LTYPE(/1)
  246. IF (LTYPE(K,2).NE.48) THEN
  247. KLISOU=KLISOU+1
  248. IPT5.LISOUS(KLISOU)=IPT1.LISOUS(LTYPE(K,1))
  249. ENDIF
  250. 5 CONTINUE
  251. ENDIF
  252.  
  253. c compteur de relaton globale
  254. LCONF3=LCONF2
  255. C=======================================================================
  256. C Debut de la boucle sur les itérations de Raffinement
  257. C=======================================================================
  258. 9 CONTINUE
  259.  
  260. c write(*,*) 'RAFF: ITERATION boucle 9',NBOUCL
  261. c write(*,*) 'ipt5 : ', ipt5
  262. c write(*,*) 'ipt6 : ', ipt6
  263. c write(*,*) 'lconf3 : ', lconf3
  264. C
  265. C=======================================================================
  266. C Renumerotation locale des noeuds du maillage
  267. C=======================================================================
  268. C ICPR(I,1)=J : Le Ieme noeud global est le Jeme local
  269. C ICPR(I,2)=J : Le Ieme noeud global appartient a J elements
  270.  
  271. SEGACT IPT5
  272.  
  273. NBOUCL=NBOUCL+1
  274. SEGINI ICPR
  275. NBNDS=0
  276. NCOL=0
  277. DO 11 INB=1,MAX(1,IPT5.LISOUS(/1))
  278. IF (IPT5.LISOUS(/1).EQ.0) THEN
  279. IPT2=IPT5
  280. ELSE
  281. IPT2=IPT5.LISOUS(INB)
  282. SEGACT IPT2
  283. ENDIF
  284. IF (IPT2.ITYPEL.EQ.48) GOTO 11
  285. DO 10 J=1,IPT2.NUM(/2)
  286. DO 110 I=1,IPT2.NUM(/1)
  287. ITMP=IPT2.NUM(I,J)
  288. IF (ICPR(ITMP,1).NE.0) THEN
  289. ICPR(ITMP,2)=ICPR(ITMP,2)+1
  290. NCOL=MAX(NCOL,ICPR(ITMP,2))
  291. ELSE
  292. NBNDS=NBNDS+1
  293. ICPR(ITMP,1)=NBNDS
  294. ICPR(ITMP,2)=1
  295. NCOL=MAX(NCOL,ICPR(ITMP,2))
  296. ENDIF
  297. 110 CONTINUE
  298. 10 CONTINUE
  299. 11 CONTINUE
  300. IF (IDIM.EQ.2) NCOL=NCOL*2
  301. IF (IDIM.EQ.3) NCOL=NCOL*4
  302.  
  303. C=======================================================================
  304. C Prise en compte des noeuds de relations
  305. C=======================================================================
  306. C IRELA(I,1) = 1 si le noeud I (en numerotation locale) est un noeud
  307. C de relation ("hanging node")
  308. C IRELA(I,1) = 0 sinon
  309. SEGINI IRELA
  310. IF (LCONF3.EQ.0) GOTO 1200
  311. SEGACT IPT6
  312. DO 1218 LAA=1,MAX(1,IPT6.LISOUS(/1))
  313. IF (IPT6.LISOUS(/1).EQ.0) THEN
  314. IPT2=IPT6
  315. ELSE
  316. IPT2=IPT6.LISOUS(LAA)
  317. ENDIF
  318. SEGACT IPT2
  319. IF (IPT2.ITYPEL.NE.48) GOTO 1218
  320. DO 1220 J=1,IPT2.NUM(/2)
  321. IRELA(ICPR(IPT2.NUM(1,J),1))=1
  322. 1220 CONTINUE
  323. 1218 CONTINUE
  324. 1200 CONTINUE
  325.  
  326. C
  327. C=======================================================================
  328. C Prise en compte du champoint de densite
  329. C=======================================================================
  330. C On affecte une densite a chacun des noeuds du maillage. On calcule
  331. C cette densite a partir des densites existantes dans le champoint, a
  332. C l'aide de la fonction RAFDEN.
  333. C Remarque : on ne prend pas en compte les densites qui pourraient
  334. C exister dans le tableau XCOOR.
  335. C
  336.  
  337. C------------------------------------------
  338. C A) Recuperation des valeurs du champoint !
  339. C------------------------------------------
  340. C XDEN(I,1) = densite affectee au noeud I (en numerotation globale)
  341. IF (NBOUCL.EQ.1) THEN
  342. SEGACT MCHPO1
  343. IF (MCHPO1.IPCHP(/1).NE.1) RETURN
  344. MSOUPO=MCHPO1.IPCHP(1)
  345. SEGACT MSOUPO
  346. IPT3=IGEOC
  347. MPOVA1=IPOVAL
  348. SEGACT IPT3
  349. SEGACT MPOVA1
  350. IF (MPOVA1.VPOCHA(/2).NE.1) RETURN
  351. NBPTCH=IPT3.NUM(/2)
  352. c SEGINI ICHPO
  353. DO 13 I=1,NBPTCH
  354. J=IPT3.NUM(1,I)
  355. c XR(I)=XCOOR((J-1)*(IDIM+1)+1)
  356. c YR(I)=XCOOR((J-1)*(IDIM+1)+2)
  357. c IF (IDIM.EQ.3) ZR(I)=XCOOR((J-1)*(IDIM+1)+3)
  358. c DEN(I)=MPOVA1.VPOCHA(I,1)
  359. c if(J.eq.115.or.J.eq.321.or.
  360. c . J.eq.6261.or.J.eq.5815) then
  361. C write(*,*) 'RAFF: #',J,'DEN(',I,'/',NBPTCH,')=',DEN(I)
  362. c endif
  363. XDEN(J)=MPOVA1.VPOCHA(I,1)
  364. 13 CONTINUE
  365. C Rem ;seulement pour les noeud "peres",
  366. C pour les "fils" raff2d remplit XDEN
  367. ENDIF
  368. C
  369. C--------------------------------------------------------------
  370. C B) Affectation d'une densite a chacun des noeuds du maillage !
  371. C--------------------------------------------------------------
  372. C IDENSI(I,1) = densite affectee au noeud I (en numerotation locale)
  373. SEGINI IDENSI
  374. DO 15 INB=1,MAX(1,IPT5.LISOUS(/1))
  375. IF (IPT5.LISOUS(/1).EQ.0) THEN
  376. IPT2=IPT5
  377. ELSE
  378. IPT2=IPT5.LISOUS(INB)
  379. SEGACT IPT2
  380. ENDIF
  381. c write(*,*) 'ipt2 : ', ipt2, 'type', IPT2.ITYPEL
  382. c write(*,*) IPT2.NUM(/2), IPT2.NUM(/1)
  383. c IF (IPT2.ITYPEL.EQ.48) GOTO 15
  384. DO 14 J=1,IPT2.NUM(/2)
  385. DO 114 I=1,IPT2.NUM(/1)
  386. NGLOB=IPT2.NUM(I,J)
  387. NLOC=ICPR(NGLOB,1)
  388. c XPT=XCOOR((NGLOB-1)*(IDIM+1)+1)
  389. c YPT=XCOOR((NGLOB-1)*(IDIM+1)+2)
  390. c ZPT=0
  391. c IF (IDIM.EQ.3) ZPT=XCOOR((NGLOB-1)*(IDIM+1)+3)
  392. c CALL RAFDEN(ICHPO,XPT,YPT,ZPT,IDIM,DENPT)
  393. c DENSI(NLOC)=DENPT
  394. DENSI(NLOC)=XDEN(NGLOB)
  395. c WRITE (*,*) 'densite au pt', NGLOB, XDEN(NGLOB), XPETI
  396. IF (XDEN(NGLOB).LT.XPETI) THEN
  397. WRITE (*,*) 'Champ de densité en entrée nul au point',
  398. + NGLOB, 'impossible de raffiner'
  399. RETURN
  400. ENDIF
  401. 114 CONTINUE
  402. 14 CONTINUE
  403. 15 CONTINUE
  404. C
  405. C=======================================================================
  406. C Remplissage du tableau KARETE
  407. C=======================================================================
  408. C On remplit le tableau KARETE avec les donnees provenant de l'analyse
  409. C du maillage.
  410. C KARETE(i,j)=k : le ieme noeud local forme une arete
  411. C avec le keme noeud local (k>i)
  412. C KARPOS(i)=k : compteur de position (k=nb d'aretes touchant i)
  413. C KARET2(i,j)=k : l'arete i-j appartient a k elements
  414. JPLA1=NBNDS
  415. JPLA2=NCOL
  416. SEGINI KARETE,KARPOS,KARET2
  417. NBARET=0
  418. NCOL=1
  419. DO 18 INB=1,MAX(1,IPT5.LISOUS(/1))
  420. IF (IPT5.LISOUS(/1).EQ.0) THEN
  421. IPT2=IPT5
  422. ELSE
  423. IPT2=IPT5.LISOUS(INB)
  424. SEGACT IPT2
  425. ENDIF
  426. c IF (IPT2.ITYPEL.EQ.48) GOTO 18
  427. ILPL=LPL(IPT2.ITYPEL)
  428. ILPT=LPT(IPT2.ITYPEL)
  429. DO 17 J=1,IPT2.NUM(/2)
  430. DO 117 K=1,ILPL*2-1,2
  431. NPTA=IPT2.NUM(KSEGM(ILPT+K-1),J)
  432. NPTB=IPT2.NUM(KSEGM(ILPT+K),J)
  433. NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
  434. NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
  435. NEXIST=0
  436. DO 16 I=1,MAX(1,KARPOS(NMIN))
  437. IF (KARETE(NMIN,I).EQ.NMAX) NEXIST=I
  438. 16 CONTINUE
  439. c création d'une nouvelle arrète
  440. IF (NEXIST.EQ.0) THEN
  441. IF ((KARPOS(NMIN)+1).GT.NCOL) THEN
  442. NCOL=KARPOS(NMIN)+1
  443. SEGADJ KARETE, KARET2
  444. ENDIF
  445. KARPOS(NMIN)=KARPOS(NMIN)+1
  446. KARETE(NMIN,KARPOS(NMIN))=NMAX
  447. KARET2(NMIN,KARPOS(NMIN))=1
  448. NBARET=NBARET+1
  449. c NCOL=MAX(NCOL,KARPOS(NMIN))
  450. ENDIF
  451. IF (NEXIST.NE.0) THEN
  452. KARET2(NMIN,NEXIST)=KARET2(NMIN,NEXIST)+1
  453. ENDIF
  454. 117 CONTINUE
  455. 17 CONTINUE
  456. 18 CONTINUE
  457. c SEGADJ KARETE,KARET2
  458. SEGDES KARETE,KARPOS,KARET2
  459. C
  460. C=======================================================================
  461. C Remplissage des tableaux JPLANS et JNOEFA
  462. C=======================================================================
  463. C On remplit les tableaux de description des faces dans le cas 3D.
  464. C JPLANS(i,j)=k : la jeme face dont le plus petit sommet en
  465. C numerotation locale est le sommet i est la face k
  466. C JNOEFA(i,1-3)=j-k-l : Les 3 plus petits sommets en numerotation
  467. C locale de la face i sont les noeuds j,k et l
  468. C JNOEFA(i,4-5)=j-k : La face identifiee par les 3 sommets precedents
  469. C appartient aux lisous de type j et k
  470. C JPLCOM(i)=k : compteur de position (k=nb de faces touchant i)
  471. IF (IDIM.NE.3) GOTO 40
  472. C----------------------------------------------------
  473. C A) Comptage du nombre d'elements total du maillage !
  474. C----------------------------------------------------
  475. NELTOT=0
  476. DO 50 INB=1,MAX(1,IPT5.LISOUS(/1))
  477. IF (IPT5.LISOUS(/1).EQ.0) THEN
  478. IPT2=IPT5
  479. ELSE
  480. IPT2=IPT5.LISOUS(INB)
  481. SEGACT IPT2
  482. ENDIF
  483. NELTOT=NELTOT+IPT2.NUM(/2)
  484. 50 CONTINUE
  485. C
  486. C-------------------------------------------------
  487. C B) Initialisation des tableaux JPLANS et JNOEFA !
  488. C-------------------------------------------------
  489. JNBFA=NELTOT*6
  490. JMAXPL=0
  491. JNUMFA=0
  492. SEGINI JPLANS,JPLCOM,JNOEFA
  493. C
  494. C-----------------------------------------------
  495. C C) Remplissage des tableaux JPLANS et JNOEFA !
  496. C-----------------------------------------------
  497. DO 51 INB=1,MAX(1,IPT5.LISOUS(/1))
  498. IF (IPT5.LISOUS(/1).EQ.0) THEN
  499. IPT2=IPT5
  500. ELSE
  501. IPT2=IPT5.LISOUS(INB)
  502. SEGACT IPT2
  503. ENDIF
  504. C IF (IPT2.ITYPEL.EQ.48) GOTO 51
  505. DO 52 KELEM=1,IPT2.NUM(/2)
  506. JTYPEL=IPT2.ITYPEL
  507. JLTEL1=LTEL(1,JTYPEL)
  508. JLTEL2=LTEL(2,JTYPEL)
  509. DO 53 IBB=1,JLTEL1
  510. JLDEL1=LDEL(1,JLTEL2-1+IBB)
  511. JTYFAC=IWORK1(JLDEL1)
  512. JLDEL2=LDEL(2,JLTEL2-1+IBB)
  513. JNBSOM=NBSOM(JTYFAC)
  514. JSPOS=NSPOS(JTYFAC)
  515. SEGINI IWORK2
  516. DO 54 IAA=1,JNBSOM
  517. NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),KELEM)
  518. IWORK2(IAA)=NGLOBA
  519. 54 CONTINUE
  520. C
  521. C----------------------------------------------------
  522. C D) Tri des sommets de la face par ordre croissant !
  523. C----------------------------------------------------
  524. C On va ranger les 3 plus petits sommets de la face (en numerotation
  525. C globale) dans les variables NPTA, NPTB, NPTC. On les classe par ordre
  526. C croissant.
  527. C
  528. C 1/ Classement des 3 sommets par ordre croissant (num globale)
  529. NPTA=nbpts+1
  530. NPTB=NPTA+1
  531. NPTC=NPTB+1
  532. DO 55 ICC=1,JNBSOM
  533. IF (IWORK2(ICC).LT.NPTA) THEN
  534. NPTC=NPTB
  535. NPTB=NPTA
  536. NPTA=IWORK2(ICC)
  537. ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  538. NPTC=NPTB
  539. NPTB=IWORK2(ICC)
  540. ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  541. NPTC=IWORK2(ICC)
  542. ENDIF
  543. 55 CONTINUE
  544. C
  545. C 2/ Classement des 3 sommets par ordre croissant (num locale)
  546. NPTA2=ICPR(NPTA,1)
  547. NPTB2=ICPR(NPTB,1)
  548. NPTC2=ICPR(NPTC,1)
  549. IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  550. NPTA=NPTA2
  551. NPTB=MIN(NPTB2,NPTC2)
  552. NPTC=MAX(NPTB2,NPTC2)
  553. ENDIF
  554. IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  555. NPTA=NPTB2
  556. NPTB=MIN(NPTA2,NPTC2)
  557. NPTC=MAX(NPTA2,NPTC2)
  558. ENDIF
  559. IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  560. NPTA=NPTC2
  561. NPTB=MIN(NPTA2,NPTB2)
  562. NPTC=MAX(NPTA2,NPTB2)
  563. ENDIF
  564. C
  565. C----------------------------------------
  566. C E) Renseignement des nouvelles faces !
  567. C----------------------------------------
  568. SEGACT LTYPE
  569. c IHERED=0
  570. c DO 1022 MOI=1,IPT2.NUM(/1)
  571. c IHERED=IHERED+IRELA(ICPR(IPT2.NUM(MOI,KELEM),1))
  572. c 1022 CONTINUE
  573.  
  574. C Si NPTA Ne correspond a aucune face : On cree une face
  575. IF (JPLCOM(NPTA).EQ.0) THEN
  576. JNUMFA=JNUMFA+1
  577. JPLCOM(NPTA)=JPLCOM(NPTA)+1
  578. JPLANS(NPTA,JPLCOM(NPTA))=JNUMFA
  579. JNOEFA(JNUMFA,1)=NPTA
  580. JNOEFA(JNUMFA,2)=NPTB
  581. JNOEFA(JNUMFA,3)=NPTC
  582. JNOEFA(JNUMFA,4)=LTYPE(INB,2)
  583. c IF (IHERED.NE.0) JNOEFA(JNUMFA,5)=LTYPE(INB,2)
  584. C On signale les sous face de faces supportant un sure
  585. C JNOEFA(JNUMFA,5)= 48
  586. IF (LCONF3.NE.0) THEN
  587. SEGACT IPT6
  588. DO 501 LAA1=1,MAX(1,IPT6.LISOUS(/1))
  589. IF (IPT6.LISOUS(/1).EQ.0) THEN
  590. IPT4=IPT6
  591. ELSE
  592. IPT4=IPT6.LISOUS(LAA1)
  593. ENDIF
  594. SEGACT IPT4
  595. IF (IPT4.ITYPEL.NE.48) GOTO 501
  596. IF (IPT4.ICOLOR(1).NE.3) GOTO 501
  597. DO 502,MAA1=1,IPT4.NUM(/2)
  598. NEXI1=0
  599. NEXI2=0
  600. DO 503 ICC1=1,JNBSOM
  601. IF (IWORK2(ICC1).EQ.IPT4.NUM(1, MAA1)) THEN
  602. NEXI1=NEXI1+1
  603. ENDIF
  604. 503 CONTINUE
  605. IF (NEXI1.GT.0) THEN
  606. DO 504 ICC2=1,JNBSOM
  607. DO 1504 INN2=2,IPT4.NUM(/1)
  608. IF (IWORK2(ICC2).EQ.IPT4.NUM(INN2,MAA1)) THEN
  609. NEXI2=NEXI2+1
  610. ENDIF
  611. 1504 CONTINUE
  612. 504 CONTINUE
  613. IF (NEXI2.GT.0) JNOEFA(JNUMFA,5)= 48
  614. ENDIF
  615. 502 CONTINUE
  616. 501 CONTINUE
  617. ENDIF
  618.  
  619.  
  620. C Si NPTA correspond a aumoins une face
  621. ELSE
  622. NEXIST=0
  623. DO 56 IDD=1,JPLCOM(NPTA)
  624. ITEMP1=JPLANS(NPTA,IDD)
  625. NAA=JNOEFA(ITEMP1,1)
  626. NBB=JNOEFA(ITEMP1,2)
  627. NCC=JNOEFA(ITEMP1,3)
  628. C Si la face que actuelle existe déja : on renseigne JNOEFA(NEXIST,5)
  629. IF (NAA.EQ.NPTA.AND.NBB.EQ.NPTB.AND.NCC.EQ.NPTC) THEN
  630. NEXIST=ITEMP1
  631. JNOEFA(NEXIST,5)=LTYPE(INB,2)
  632. ENDIF
  633.  
  634.  
  635. 56 CONTINUE
  636. C Si la face que actuelle n'existe pas : On cree une face
  637. IF (NEXIST.EQ.0) THEN
  638. JNUMFA=JNUMFA+1
  639. JPLCOM(NPTA)=JPLCOM(NPTA)+1
  640. JPLANS(NPTA,JPLCOM(NPTA))=JNUMFA
  641. JNOEFA(JNUMFA,1)=NPTA
  642. JNOEFA(JNUMFA,2)=NPTB
  643. JNOEFA(JNUMFA,3)=NPTC
  644. JNOEFA(JNUMFA,4)=LTYPE(INB,2)
  645. c IF (IHERED.NE.0) JNOEFA(JNUMFA,5)=LTYPE(INB,2)
  646.  
  647. C On signale les sous face de faces supportant un sure
  648. C JNOEFA(JNUMFA,5)= 48
  649. IF (LCONF3.NE.0) THEN
  650. SEGACT IPT6
  651. DO 505 LAA1=1,MAX(1,IPT6.LISOUS(/1))
  652. IF (IPT6.LISOUS(/1).EQ.0) THEN
  653. IPT4=IPT6
  654. ELSE
  655. IPT4=IPT6.LISOUS(LAA1)
  656. ENDIF
  657. SEGACT IPT4
  658. IF (IPT4.ITYPEL.NE.48) GOTO 505
  659. IF (IPT4.ICOLOR(1).NE.3) GOTO 505
  660. DO 506,MAA1=1,IPT4.NUM(/2)
  661. NEXI1=0
  662. NEXI2=0
  663. DO 507 ICC1=1,JNBSOM
  664. IF (IWORK2(ICC1).EQ.IPT4.NUM(1, MAA1)) THEN
  665. NEXI1=NEXI1+1
  666. ENDIF
  667. 507 CONTINUE
  668. IF (NEXI1.GT.0) THEN
  669. DO 508 ICC2=1,JNBSOM
  670. DO 1508 INN2=2,IPT4.NUM(/1)
  671. IF (IWORK2(ICC2).EQ.IPT4.NUM(INN2,MAA1)) THEN
  672. NEXI2=NEXI2+1
  673. ENDIF
  674. 1508 CONTINUE
  675. 508 CONTINUE
  676. IF (NEXI2.GT.0) JNOEFA(JNUMFA,5)= 48
  677. ENDIF
  678. 506 CONTINUE
  679. 505 CONTINUE
  680. ENDIF
  681. ENDIF
  682. ENDIF
  683.  
  684. JMAXPL=MAX(JMAXPL,JPLCOM(NPTA))
  685. 53 CONTINUE
  686. 52 CONTINUE
  687. 51 CONTINUE
  688. JPLA2=JMAXPL
  689. JNBFA=JNUMFA
  690. SEGADJ JPLANS,JPLCOM,JNOEFA
  691. SEGDES LTYPE
  692. C
  693. C=======================================================================
  694. C Pre-remplissage du tableau KMILIE
  695. C=======================================================================
  696. C On remplit le tableau lorsqu'il y a des maillages de relations de
  697. C type 48 (IPT6, LCONF3 non nul)
  698. C KARETE(i,j)=k ET KMILIE(i,j)=n : le noeud global n est situe
  699. C au milieu de l'arete [i,k]
  700. 40 CONTINUE
  701. SEGINI KMILIE
  702. SEGACT KARETE*MOD,KARPOS*MOD,KARET2*MOD
  703. IF (LCONF3.EQ.0) THEN
  704. GOTO 28
  705. ENDIF
  706. 22 CONTINUE
  707. NBRELA=0
  708. SEGACT IPT6
  709. DO 218 LAA=1,MAX(1,IPT6.LISOUS(/1))
  710. IF (IPT6.LISOUS(/1).EQ.0) THEN
  711. IPT2=IPT6
  712. ELSE
  713. IPT2=IPT6.LISOUS(LAA)
  714. ENDIF
  715. SEGACT IPT2
  716. IF (IPT2.ITYPEL.NE.48) GOTO 218
  717. IF ((IPT2.ICOLOR(1).NE.1).AND.(IPT2.ICOLOR(1).NE.2)) GOTO 218
  718. DO 220 J=1,IPT2.NUM(/2)
  719.  
  720. NPTA=IPT2.NUM(2,J)
  721. NPTB=IPT2.NUM(3,J)
  722. NPTC=IPT2.NUM(1,J)
  723. NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
  724. NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
  725. NEXIST=0
  726. DO 219 I=1,MAX(1,KARPOS(NMIN))
  727. IF (KARETE(NMIN,I).EQ.NMAX) NEXIST=I
  728. 219 CONTINUE
  729.  
  730. IF (NEXIST.NE.0) then
  731. KMILIE(NMIN,NEXIST)=NPTC
  732. c WRITE (*,*) 'REPORT ANCIENNE RELA AU NOEUD:', NPTC
  733. c write (*,*) 'karret2(nmin nexist)', KARET2(NMIN,NEXIST)
  734. KARET2(NMIN,NEXIST)=2
  735. ELSE
  736. IF (IDIM.EQ.2) goto 27
  737. IF ((KARPOS(NMIN)+1).GT.NCOL) THEN
  738. NCOL=KARPOS(NMIN)+1
  739. SEGADJ KARETE, KARET2, KMILIE
  740. ENDIF
  741. JPOSI = KARPOS(NMIN)
  742. KARPOS(NMIN)=JPOSI+1
  743. KARETE(NMIN,JPOSI+1)=NMAX
  744. KARET2(NMIN,JPOSI+1)=2
  745. KMILIE(NMIN,JPOSI+1)=NPTC
  746. c WRITE (*,*) 'création de l arrete:', NMIN, NMAX
  747. ENDIF
  748. C
  749. C-----------------------------------------------
  750. C Prise en compte des relations "hereditaires" !
  751. C-----------------------------------------------
  752. C Cette partie sert a traiter le cas ou un element serait raffine
  753. C plusieurs fois alors qu'a sa frontiere les autres elements ne le
  754. C seraient pas. Il y aurait donc plusieurs points sur une arete ayant
  755. C des relations de conformite et plus seulement le noeud milieu.
  756. C
  757. c dans l'exenple ci dessous où les noeuds C et D sont des hanging nodes
  758. c L'arette [C B] n'a pas encore été créée c'est ce qui est fait ici.
  759. c
  760. c | |
  761. c | |
  762. c |A C D |B
  763. c x---------X----X----X
  764. c | | | |
  765. c | | | |
  766. 27 CONTINUE
  767. IF (IDIM.EQ.2) THEN
  768. C 1/ Test sur la premiere arete : NMIN-NPTC
  769. IF (NEXIST.NE.0) THEN
  770. NEXIS2=0
  771. NMIN2=MIN(NMIN,ICPR(NPTC,1))
  772. NMAX2=MAX(NMIN,ICPR(NPTC,1))
  773. DO 24 I=1,MAX(1,KARPOS(NMIN2))
  774. IF (KARETE(NMIN2,I).EQ.NMAX2) NEXIS2=I
  775. 24 CONTINUE
  776. IF (NEXIS2.EQ.0) THEN
  777. NBRELA=NBRELA+1
  778. JNBCOL=KARETE(/2)
  779. JPOSI=KARPOS(NMIN2)
  780. IF (JNBCOL.GT.JPOSI) THEN
  781. KARPOS(NMIN2)=JPOSI+1
  782. KARETE(NMIN2,JPOSI+1)=NMAX2
  783. KARET2(NMIN2,JPOSI+1)=2
  784. ELSEIF (JNBCOL.EQ.JPOSI) THEN
  785. NCOL=JNBCOL+1
  786. SEGADJ KARETE,KMILIE,KARET2
  787. KARETE(NMIN2,JPOSI+1)=NMAX2
  788. KARPOS(NMIN2)=JPOSI+1
  789. KARET2(NMIN2,JPOSI+1)=2
  790. ENDIF
  791. c WRITE (*,*) 'création de l arrete:', NMIN2, NMAX2
  792. ELSE
  793. KARET2(NMIN2,NEXIS2)=2
  794. ENDIF
  795. C
  796. C 2/ Test sur la deuxieme arete : NMAX-NPTC
  797. NEXIS2=0
  798. NMIN2=MIN(NMAX,ICPR(NPTC,1))
  799. NMAX2=MAX(NMAX,ICPR(NPTC,1))
  800. DO 25 I=1,MAX(1,KARPOS(NMIN2))
  801. IF (KARETE(NMIN2,I).EQ.NMAX2) NEXIS2=I
  802. 25 CONTINUE
  803. IF (NEXIS2.EQ.0) THEN
  804. NBRELA=NBRELA+1
  805. JNBCOL=KARETE(/2)
  806. JPOSI=KARPOS(NMIN2)
  807. IF (JNBCOL.GT.JPOSI) THEN
  808. KARPOS(NMIN2)=JPOSI+1
  809. KARETE(NMIN2,JPOSI+1)=NMAX2
  810. KARET2(NMIN2,JPOSI+1)=2
  811. ELSEIF (JNBCOL.EQ.JPOSI) THEN
  812. NCOL=JNBCOL+1
  813. SEGADJ KARETE,KMILIE,KARET2
  814. KARETE(NMIN2,JPOSI+1)=NMAX2
  815. KARPOS(NMIN2)=JPOSI+1
  816. KARET2(NMIN2,JPOSI+1)=2
  817. ENDIF
  818. c WRITE (*,*) 'création de l arrete:', NMIN2, NMAX2
  819. ELSE
  820. KARET2(NMIN2,NEXIS2)=2
  821. ENDIF
  822. ENDIF
  823. ENDIF
  824. 220 CONTINUE
  825. 218 CONTINUE
  826. IF (NBRELA.NE.0) GOTO 22
  827. 28 CONTINUE
  828.  
  829. C
  830. C
  831. C=======================================================================
  832. C Pre-remplissage du tableau JFARAF
  833. C=======================================================================
  834. C On remplit le tableau JFARAF et JPLAN3 lorsqu'il y a des maillages de
  835. C relations de type 48 (IPT6, LCONF3 non nul).
  836. c
  837. C Si JPLANS(i,j)=k et JFARAF(k , (l-1)*2+1)= n : le noeud n (en
  838. C numérotation globale) est le hanging node au centre de la jieme face
  839. C dont le plus petit sommet en numerotation locale est le sommet i.
  840. C
  841. C Si de plus la face en question est un TRI6 et JFARAF(k , l*2)=m :
  842. C le noeud m (en numérotation globale) est le 4ieme noeud de l'element
  843. C 48 coul 4
  844. C
  845. C Si de plus la face en question est un QUA8-1 et JFARAF(k , l*2)=m :
  846. C le noeud m (en numérotation globale) est le 5eme noeud de l'element
  847. c 48 coul 4
  848. C
  849. C si JPLAN3(i,j)=k il y a k relation de compatibilité de faces associé a
  850. C la jieme face dont le plus petit sommet en numerotation locale est le
  851. C sommet i.
  852.  
  853. IF (IDIM.NE.3) GOTO 711
  854. JNBFA=JNOEFA(/1)
  855. LLLL=10
  856. SEGINI JFARAF
  857. SEGINI JPLAN2,JPLAN3
  858. IF (LCONF3.EQ.0) GOTO 711
  859. MPOS25=0
  860. JNBSOM=0
  861. 1287 CONTINUE
  862. NBRELA=0
  863. DO 670 MIE=1,MAX(1,IPT6.LISOUS(/1))
  864. IF (IPT6.LISOUS(/1).EQ.0) THEN
  865. IPT2=IPT6
  866. ELSE
  867. IPT2=IPT6.LISOUS(MIE)
  868. ENDIF
  869. SEGACT IPT2
  870. IF (IPT2.ITYPEL.NE.48) GOTO 670
  871. IF ((IPT2.ICOLOR(1).EQ.1).OR.(IPT2.ICOLOR(1).EQ.2)) GOTO 670
  872. LIGNUM=IPT2.NUM(/1)
  873. IF (IPT2.ICOLOR(1).EQ.4) MPOS25=4
  874. IF (IPT2.ICOLOR(1).EQ.5) MPOS25=5
  875. IF (IPT2.ICOLOR(1).EQ.3) JNBSOM=4
  876. IF (IPT2.ICOLOR(1).EQ.4) JNBSOM=3
  877. IF (IPT2.ICOLOR(1).EQ.5) JNBSOM=4
  878. IF (IPT2.ICOLOR(1).EQ.6) JNBSOM=4
  879. SEGINI IWORK2
  880. C
  881. C-------------------------
  882. C A) Recherche de la face !
  883. C-------------------------
  884. C 1/ Classement des 3 sommets par ordre croissant (num globale)
  885. DO 671,MAA=1,IPT2.NUM(/2)
  886. DO 673 MBB=1,JNBSOM
  887. IWORK2(MBB)=IPT2.NUM(LIGNUM-JNBSOM+MBB,MAA)
  888. 673 CONTINUE
  889. NPTA=nbpts+1
  890. NPTB=NPTA+1
  891. NPTC=NPTB+1
  892. DO 672 ICC=1,JNBSOM
  893. IF (IWORK2(ICC).LT.NPTA) THEN
  894. NPTC=NPTB
  895. NPTB=NPTA
  896. NPTA=IWORK2(ICC)
  897. ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  898. NPTC=NPTB
  899. NPTB=IWORK2(ICC)
  900. ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  901. NPTC=IWORK2(ICC)
  902. ENDIF
  903. 672 CONTINUE
  904. C
  905. C 2/ Classement des 3 sommets par ordre croissant (num locale)
  906. NPTA1=NPTA
  907. NPTB1=NPTB
  908. NPTC1=NPTC
  909. NPTA2=ICPR(NPTA,1)
  910. NPTB2=ICPR(NPTB,1)
  911. NPTC2=ICPR(NPTC,1)
  912. IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  913. NPTA=NPTA2
  914. NPTB=MIN(NPTB2,NPTC2)
  915. NPTC=MAX(NPTB2,NPTC2)
  916. ENDIF
  917. IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  918. NPTA=NPTB2
  919. NPTB=MIN(NPTA2,NPTC2)
  920. NPTC=MAX(NPTA2,NPTC2)
  921. ENDIF
  922. IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  923. NPTA=NPTC2
  924. NPTB=MIN(NPTA2,NPTB2)
  925. NPTC=MAX(NPTA2,NPTB2)
  926. ENDIF
  927. C
  928. C 3/ Recherche du numero de la face
  929. NEXIS2=0
  930. DO 674 IEE=1,JPLCOM(NPTA)
  931. MTMP=JPLANS(NPTA,IEE)
  932. JJ1=JNOEFA(MTMP,1)
  933. JJ2=JNOEFA(MTMP,2)
  934. JJ3=JNOEFA(MTMP,3)
  935. IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
  936. NEXIS2=IEE
  937. ENDIF
  938. 674 CONTINUE
  939. c 3.1/ si la face en question n'existe pas, dans le cas général on
  940. c s'arrète
  941. IF ((NEXIS2.EQ.0).AND.(IPT2.ICOLOR(1).NE.3)) THEN
  942. c WRITE (*,*) 'IPT2.ICOLOR(1)', IPT2.ICOLOR(1)
  943. WRITE (*,*) 'Raffinement impossible, le champ de densite
  944. + presente des variations trop brusques'
  945. WRITE (*,*) 'sure concerne', MAA
  946. WRITE (*,*) 'premier nd du sure concerne', IPT2.NUM(1,MAA)
  947. WRITE (*,*) 'numero du 1er nd de la face concernee', NPTA
  948. KARAF2=0
  949. c GOTO 679
  950. ENDIF
  951. c 3.2/ si la face en question n'existe pas dans le cas d'une relation
  952. c de face QUA4 il suffie de la créer
  953.  
  954. IF ((NEXIS2.EQ.0).AND.(IPT2.ICOLOR(1).EQ.3)) THEN
  955. JPLCOM(NPTA)= JPLCOM(NPTA)+1
  956. JNUMFA=JNUMFA+1
  957. JNUMFS=JNUMFA
  958. JNUMF1=JPLCOM(NPTA)
  959. JMAXPL=MAX(JMAXPL,JPLCOM(NPTA))
  960. c 3.2.1/ ajustement des segments
  961. JPLA2=JMAXPL
  962. JNBFA=JNUMFA
  963. SEGADJ JPLANS
  964. SEGADJ JNOEFA
  965. SEGADJ JPLAN3
  966. SEGADJ JFARAF
  967. C 3.2.1 Renseignement des iformation de la nouvelle face dans JPLANS et
  968. c JNOEFA, initialisation de JPLAN3
  969. JPLANS(NPTA,JNUMF1)=JNUMFA
  970. JPLAN3(NPTA,JNUMF1)=0
  971. JNOEFA(JNUMFA,1)= NPTA
  972. JNOEFA(JNUMFA,2)= NPTB
  973. JNOEFA(JNUMFA,3)= NPTC
  974. JNOEFA(JNUMFA,4)= 48
  975. JNOEFA(JNUMFA,5)= 48
  976. ENDIF
  977.  
  978. IF (NEXIS2.NE.0) THEN
  979. JNUMFS=JPLANS(NPTA,NEXIS2)
  980. ENDIF
  981. C
  982. C 4/ Recherche de la position de la face dans JPLANS
  983. NEXIS2=0
  984. DO 675 MAC=1,JPLCOM(NPTA)
  985. IF (JPLANS(NPTA,MAC).EQ.JNUMFS) NEXIS2=MAC
  986. 675 CONTINUE
  987. C
  988. C-------------------------------------
  989. C B) Renseignement de JFARAF et JPLAN3!
  990. C-------------------------------------
  991.  
  992.  
  993. JPLAN3(NPTA,NEXIS2)=JPLAN3(NPTA,NEXIS2)+1
  994. LTEMPO=(JPLAN3(NPTA,NEXIS2)-1)*2+1
  995. JFARAF(JNUMFS,LTEMPO)=IPT2.NUM(1,MAA)
  996. IF (MPOS25.NE.0) JFARAF(JNUMFS,LTEMPO+1)=IPT2.NUM(MPOS25,MAA)
  997. IF (JNOEFA(JNUMFS,5).EQ.0) JNOEFA(JNUMFS,5)=48
  998.  
  999.  
  1000. 671 CONTINUE
  1001. 670 CONTINUE
  1002. C IF (NBRELA.NE.0) GOTO 1287
  1003. 711 CONTINUE
  1004. C
  1005. C=======================================================================
  1006. C Calcul du critere et raffinement si necessaire
  1007. C=======================================================================
  1008. C Pour calculer le critere, on a besoin d'un modele et d'un chamelem
  1009. C constant egal a 1.
  1010. C
  1011. C------------------------------------------
  1012. C A) Creation d'un modele simple : MMODE1 !
  1013. C------------------------------------------
  1014. MN3 =1
  1015. NFOR=1
  1016. NMAT=2
  1017. NPARMO=0
  1018. NOBMOD=0
  1019. N1=1
  1020. SEGINI MMODE1
  1021. SEGINI IMODE1
  1022. MMODE1.KMODEL(1)=IMODE1
  1023. IMODE1.FORMOD(1)='MECANIQUE'
  1024. IMODE1.MATMOD(1)='ELASTIQUE'
  1025. IMODE1.MATMOD(2)='ISOTROPE'
  1026.  
  1027. KARAF2=0
  1028. C=======================================================================
  1029. C Boucle sur les eventuels sous-domaines de IPT5
  1030. C=======================================================================
  1031.  
  1032. c-----------------------------------------------------------------------
  1033. c initialisation de LTYPE3 qui va contenir les tableau LTYP2
  1034. c correspondant à chaque sous domaine de IPT5. Ces tableau LTYP2
  1035. c contiendront eux même les maillages raffinés issue des
  1036. c différents sous domaine de IPT5.
  1037. c LTYPE3(I,1)=LTYP2 pour le sous domaine I de IPT5
  1038. c LTYPE3(I,2)=IRR le permier indice d'un maillage 48 dans LTYP2
  1039. c-----------------------------------------------------------------------
  1040. NBLIS3=MAX(1,IPT5.LISOUS(/1))
  1041. SEGINI LTYPE3
  1042. MCONF=0
  1043. NSUR = 0
  1044. DO 33 INB=1,MAX(1,IPT5.LISOUS(/1))
  1045. NBNN=0
  1046. NBELEM=0
  1047. NBREF=0
  1048. NBSOUS=0
  1049. SEGINI IPT2
  1050. IF (IPT5.LISOUS(/1).EQ.0) THEN
  1051. IPT2=IPT5
  1052. ELSE
  1053. IPT2=IPT5.LISOUS(INB)
  1054. SEGACT IPT2
  1055. ENDIF
  1056.  
  1057. IMODE1.IMAMOD=IPT2
  1058. IMODE1.NEFMOD=IPT2.ITYPEL
  1059. CALL PRQUOI(IMODE1)
  1060.  
  1061. C-------------------------------------------------------
  1062. C B) Creation d'un chamelem constant egal a 1 : MCHEL1 !
  1063. C-------------------------------------------------------
  1064. CALL ZEROP(MMODE1,'NOEUD ',MCHEL1)
  1065. SEGACT MCHEL1*MOD
  1066. MCHEL1.IMACHE(1)=IPT2
  1067. IF (IDIM.EQ.2) MCHEL1.IFOCHE=-1
  1068. IF (IDIM.EQ.3) MCHEL1.IFOCHE=2
  1069. MCHAM1=MCHEL1.ICHAML(1)
  1070. SEGACT MCHAM1
  1071. MELVA1=MCHAM1.IELVAL(1)
  1072. SEGACT MELVA1*MOD
  1073. MELVA1.VELCHE(1,1)=1.D0
  1074. C
  1075. C-------------------------------------------------------------
  1076. C C) Calcul de la surface ou du volume des elements : MCHEL2 !
  1077. C-------------------------------------------------------------
  1078. CALL INTGCA(MMODE1,MCHEL1,0,1,IRET,XRET,MCHEL2)
  1079. IF (XRET.LT.XPETI) THEN
  1080. WRITE (*,*) 'Volume(resp surf) de la zone a raffiner nul'
  1081. RETURN
  1082. ENDIF
  1083. SEGACT IPT2
  1084.  
  1085. C---------------------------------------------------------------------
  1086. C D) Creation d'un chamelem indiquant les elements a raffiner MCHAM2 !
  1087. C---------------------------------------------------------------------
  1088. SEGACT KMILIE*MOD
  1089. SEGACT MCHEL2
  1090. MCHAM2=MCHEL2.ICHAML(1)
  1091. SEGACT MCHAM2
  1092. MELVA2=MCHAM2.IELVAL(1)
  1093. SEGACT MELVA2*MOD
  1094. NACREE=0
  1095. KARAF=0
  1096.  
  1097. SEGACT MELVA1
  1098. c boucles sur les éléments de ipt2
  1099. DO 32 J=1,IPT2.NUM(/2)
  1100.  
  1101. XMOY=0
  1102. DO 29 I=1,IPT2.NUM(/1)
  1103. NGLOB=IPT2.NUM(I,J)
  1104. NLOC=ICPR(NGLOB,1)
  1105. DENPT=DENSI(NLOC)
  1106. XMOY=XMOY+DENPT/(IPT2.NUM(/1))
  1107. 29 CONTINUE
  1108.  
  1109. XINT=MELVA2.VELCHE(1,J)
  1110. XDIM=IDIM*1.0
  1111. C Un element est a raffiner si la moyenne des densites affectees a
  1112. C l'ensemble de ses noeuds est inférieure a la racine carree de sa
  1113. C surface (en 2D) ou a la racine cubique de son volume (en 3D)
  1114. C WRITE (*,*) 'XMOY', XMOY, '(XINT**(1/XDIM))',(XINT**(1/XDIM))
  1115. IF ((XMOY-(XINT**(1/XDIM))).LT.0) THEN
  1116. c c y a t'il des +epsilon...?
  1117. c XDIFF = XMOY-(XINT**(1/XDIM))
  1118. c if(XDIFF.ge.0.and.XDIFF.lt.0.001*XMOY) then
  1119. c write(*,*) 'RAFF : XMOY,XDIFF=',XMOY,XDIFF,
  1120. c & 'pour l EF',J,(IPT2.NUM(iou,J),iou=1,IPT2.NUM(/1))
  1121. c endif
  1122. c IF (XDIFF.LT.0) THEN
  1123. MELVA2.VELCHE(1,J)=1
  1124. c KARAF nombre d'éléments à rafiner dans le sous domaine
  1125. KARAF=KARAF+1
  1126. ILPT=LPT(IPT2.ITYPEL)
  1127. ILPL=LPL(IPT2.ITYPEL)
  1128.  
  1129. c boucle sur les segments de l'élément à raffiner
  1130. DO 31 K=1,ILPL*2-1,2
  1131. NPTA=IPT2.NUM(KSEGM(ILPT+K-1),J)
  1132. NPTB=IPT2.NUM(KSEGM(ILPT+K),J)
  1133. NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
  1134. NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
  1135. NEXIST=0
  1136. DO 30 I=1,MAX(1,KARPOS(NMIN))
  1137. IF (KARETE(NMIN,I).EQ.NMAX) NEXIST=I
  1138. 30 CONTINUE
  1139. c On met dans kmilie les arêtes où il faudra mettre un nouveau noeud
  1140. c On leur associe la valeur -1
  1141. IF (KMILIE(NMIN,NEXIST).EQ.0) THEN
  1142. NACREE=NACREE+1
  1143. KMILIE(NMIN,NEXIST)=-1
  1144. ENDIF
  1145. 31 CONTINUE
  1146. C
  1147. C-------------------------------------------
  1148. C E) Comptage du nombre de faces a raffiner !
  1149. C-------------------------------------------
  1150. C 1/ Initialisations
  1151. IF (IDIM.NE.3) GOTO 32
  1152. JTYPEL=IPT2.ITYPEL
  1153. JLTEL1=LTEL(1,JTYPEL)
  1154. JLTEL2=LTEL(2,JTYPEL)
  1155. DO 73 IBB=1,JLTEL1
  1156. JLDEL1=LDEL(1,JLTEL2-1+IBB)
  1157. JTYFAC=IWORK1(JLDEL1)
  1158. JLDEL2=LDEL(2,JLTEL2-1+IBB)
  1159. C JNBPTS=NBNNE(JTYFAC)
  1160. JNBSOM=NBSOM(JTYFAC)
  1161. JSPOS=NSPOS(JTYFAC)
  1162. SEGINI IWORK2
  1163. DO 74 IAA=1,JNBSOM
  1164. NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),J)
  1165. IWORK2(IAA)=NGLOBA
  1166. 74 CONTINUE
  1167. C
  1168. C 2/ Classement des 3 sommets par ordre croissant (num globale)
  1169. NPTA=nbpts+1
  1170. NPTB=NPTA+1
  1171. NPTC=NPTB+1
  1172. DO 75 ICC=1,JNBSOM
  1173. IF (IWORK2(ICC).LT.NPTA) THEN
  1174. NPTC=NPTB
  1175. NPTB=NPTA
  1176. NPTA=IWORK2(ICC)
  1177. ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  1178. NPTC=NPTB
  1179. NPTB=IWORK2(ICC)
  1180. ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  1181. NPTC=IWORK2(ICC)
  1182. ENDIF
  1183. 75 CONTINUE
  1184. C
  1185. C 3/ Classement des 3 sommets par ordre croissant (num locale)
  1186. NPTA2=ICPR(NPTA,1)
  1187. NPTB2=ICPR(NPTB,1)
  1188. NPTC2=ICPR(NPTC,1)
  1189. IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  1190. NPTA=NPTA2
  1191. NPTB=MIN(NPTB2,NPTC2)
  1192. NPTC=MAX(NPTB2,NPTC2)
  1193. ENDIF
  1194. IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  1195. NPTA=NPTB2
  1196. NPTB=MIN(NPTA2,NPTC2)
  1197. NPTC=MAX(NPTA2,NPTC2)
  1198. ENDIF
  1199. IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  1200. NPTA=NPTC2
  1201. NPTB=MIN(NPTA2,NPTB2)
  1202. NPTC=MAX(NPTA2,NPTB2)
  1203. ENDIF
  1204. C
  1205. C 4/ Calcul du nombre de points a creer pour raffiner la face
  1206. NEXIS2=0
  1207. DO 76 IEE=1,JPLCOM(NPTA)
  1208. MTMP=JPLANS(NPTA,IEE)
  1209. JJ1=JNOEFA(MTMP,1)
  1210. JJ2=JNOEFA(MTMP,2)
  1211. JJ3=JNOEFA(MTMP,3)
  1212. IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
  1213. NEXIS2=IEE
  1214. ENDIF
  1215. 76 CONTINUE
  1216. IF (NEXIS2.NE.0) THEN
  1217. KFARAF=JPLAN2(NPTA,NEXIS2)
  1218. IF (KFARAF.EQ.0) THEN
  1219. NACREE=NACREE+NBINTE(JTYFAC)
  1220. JPLAN2(NPTA,NEXIS2)=JPLAN2(NPTA,NEXIS2)+1
  1221. ENDIF
  1222. ENDIF
  1223. 73 CONTINUE
  1224. ELSE
  1225.  
  1226. MELVA2.VELCHE(1,J)=0
  1227. ENDIF
  1228. 32 CONTINUE
  1229.  
  1230. c KARAF2 nombre d'éléments à rafiner sur tout les sous domaines
  1231. KARAF2=KARAF2+KARAF
  1232. c si aucun élément à raffiner on stoque directement IPT2 dans LTYPE3
  1233. c on stoque également les anciennes relations et on passe au
  1234. c sous-domaine suivant.
  1235.  
  1236. c ATTENTION SI NBOUCL=1 il est possible que des relations surabondantes
  1237.  
  1238. c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1239. c forcage : nb max itération de raff
  1240. c IF (NBOUCL.EQ.3) THEN
  1241. c KARAF2=0
  1242. c KARAF=0
  1243. c ENDIF
  1244. c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1245. IF (KARAF.EQ.0) THEN
  1246. IF (NBOUCL.EQ.1) THEN
  1247. IF (LCONF3.GT.0) THEN
  1248. SEGACT IPT6
  1249. NBLIS2=1+MAX(1,IPT6.LISOUS(/1))
  1250. SEGINI LTYPE2
  1251. LTYPE2(1,1)=1
  1252. LTYPE2(1,2)=IPT2.ITYPEL
  1253. LTYPE2(1,3)=IPT2
  1254. LTYPE2(1,4)=IPT2.ICOLOR(1)
  1255. LTYPE3(INB,2)=0
  1256. DO 676 MIE=1,MAX(1,IPT6.LISOUS(/1))
  1257. IF (IPT6.LISOUS(/1).EQ.0) THEN
  1258. IPT4=IPT6
  1259. ELSE
  1260. IPT4=IPT6.LISOUS(MIE)
  1261. ENDIF
  1262. NSUR=1
  1263. SEGACT IPT4
  1264. LTYPE2(1+MIE,1)=1+MIE
  1265. LTYPE2(1+MIE,2)=IPT4.ITYPEL
  1266. LTYPE2(1+MIE,3)=IPT4
  1267. LTYPE2(1+MIE,4)=IPT4.ICOLOR(1)
  1268. LTYPE3(INB,2)=2
  1269. 676 CONTINUE
  1270. ELSE
  1271. NBLIS2=1
  1272. SEGINI LTYPE2
  1273. LTYPE2(1,1)=1
  1274. LTYPE2(1,2)=IPT2.ITYPEL
  1275. LTYPE2(1,3)=IPT2
  1276. LTYPE2(1,4)=IPT2.ICOLOR(1)
  1277. LTYPE3(INB,2)=0
  1278. ENDIF
  1279. LTYPE3(INB,1)=LTYPE2
  1280. ELSE
  1281. SEGACT LTYPE4
  1282. LTYPE3(INB,1)=LTYPE4(INB,1)
  1283. LTYPE3(INB,2)=LTYPE4(INB,2)
  1284. ENDIF
  1285. GOTO 33
  1286. ENDIF
  1287. C
  1288. C----------------------------
  1289. C F) Raffinement du maillage !
  1290. C----------------------------
  1291. NML=2
  1292. C
  1293. C 1/ Raffinement d'un maillage 2D
  1294. C On recupere dans MELEME le maillage raffine sans ses relations de
  1295. C conformite et dans IPT8 les relations de ce maillage
  1296. C (IPT8 encore incomplet en sortie de RAFF2D)
  1297. c entrees de raff2D:
  1298. c - IPT2: maillage élémentaire à raffiner
  1299. c (ici sous domaine de ipt5)
  1300. c - ICPR: tableau de passage noeuds local/global
  1301. c (ici à partir de ipt5 interne à la boucle 9)
  1302. c - KARPOS: tableau du nb d'arretes par noeuds
  1303. c (ici à partir de ipt5, mis a jour avec ipt6)
  1304. c - KARETE: tableau des d'arretes KARETE(i,j)=k
  1305. c (ici à partir de ipt5, mis a jour avec ipt6)
  1306. c - KMILIE: tableau des hanging nodes associés au arrètes
  1307. c en entrée: à partir des relations de ipt6 et
  1308. c des nouveau noeuds à construires
  1309. c si relation issue de ipt6
  1310. c KMILIE(i,j)=n : le noeud global n est situe
  1311. c au milieu de l'arete [i,k]
  1312. c si nouvelle relation
  1313. c KMILIE(i,j)=-1 : il faut construire un noeud
  1314. c au milieu de l'arete [i,k]
  1315. c - MELVA2: champ par elem valait 1 pour les élem à raffiner
  1316. c et 0 pour les autres
  1317. c (ici à partir de ipt2 et de densi)
  1318. c - NACREE: nombre de noeuds à créer dans ipt2
  1319. c - KARAF: nombre d'élément à raffiner dans ipt2
  1320. c - IDCP : ??? variable non utilisé ???
  1321. c - XDEN : densité aux noeuds en notation globale
  1322. c (ici à partir de mchpo1)
  1323. c - KARET2: tableau du nombre d'éléments touchant les arretes
  1324. c (ici à partir de ipt5, mis a jour avec ipt6)
  1325. csorties de raff2D :
  1326. c - KMILIE: tableau des hanging nodes associés au arrètes
  1327. c en sortie: à partir des relations donnée en
  1328. c entrée et de celles crées
  1329. c si anciene relation à garder
  1330. c KMILIE(i,j)=n : le noeud global n est situe
  1331. c au milieu de l'arete [i,k]
  1332. c si anciene relation à suprimer
  1333. c KMILIE(i,j)=0 :
  1334. c si nouvelle relation
  1335. c KMILIE(i,j)=n : le noeud n à été construit au
  1336. c milieu de l'arete [i,k]
  1337. c - MELEME: Maillage élémentaire raffiné à partir de ipt2
  1338. c sans relations
  1339. c - IPT8 : Maillage de relation incomplet
  1340. c si J nouvel élément 48:
  1341. c IPT8.NUM(1,j) = nouveau hanging node
  1342. c IPT8.NUM(2:4,j) = autres noeuds de la relation
  1343. c IPT8.NUM(5,j) = 1 si arête a 2 noeuds
  1344. c 2 si arête a 3 noeuds
  1345. c si J nouvel élément 48:
  1346. c IPT8.NUM(1,j) = nouveau hanging node
  1347. c IPT8.NUM(2:5,j) = 0
  1348. c - XDEN : densité aux noeuds en notation globale
  1349. c interpolé aux nouveaux noeuds.
  1350. IF (IDIM.EQ.2) CALL RAFF2D(IPT2,ICPR,KARPOS,KARETE,KMILIE,
  1351. + MELVA2,NACREE,KARAF,MELEME,IDCP,IPT8,KARET2,XDEN)
  1352. C
  1353. C 2/ Raffinement d'un maillage 3D
  1354. IF (IDIM.EQ.3) THEN
  1355. NML=0
  1356. IF ((IPT2.ITYPEL.EQ.25).OR.(IPT2.ITYPEL.EQ.26)) NML=1
  1357. ENDIF
  1358. C a) Maillage 3D contenant des pyramides
  1359. IPT9=0
  1360. IF (NML.EQ.1) CALL RAFPYR(IPT2,ICPR,KARPOS,KARETE,KMILIE,
  1361. + MELVA2,NACREE,KARAF,MELEME,JPLANS,JPLAN3,JPLCOM,JNOEFA,
  1362. + IPT8,JFARAF,KARET2,IPT9,XDEN)
  1363. C b) Maillage 3D sans pyramides
  1364. IF (NML.EQ.0) CALL RAFF3D(IPT2,IPT6,ICPR,KARPOS,KARETE,KMILIE,
  1365. + MELVA2,NACREE,KARAF,MELEME,JPLANS,JPLAN3,JPLCOM,JNOEFA,
  1366. + IPT8,JFARAF,KARET2,XDEN)
  1367. C
  1368. C=======================================================================
  1369. C Creation du maillage de relations de conformite
  1370. C=======================================================================
  1371. C--------------------
  1372. C A) Initialisations !
  1373. C--------------------
  1374. IPT10=0
  1375. IPT11=0
  1376. IPT12=0
  1377. IPT13=0
  1378. IPT14=0
  1379. IPT15=0
  1380. IPT16=0
  1381. MCONF=0
  1382. c IF (KARAF2.EQ.0) GOTO 33
  1383. IF (IPT8.NE.0) THEN
  1384. C
  1385. C---------------------------------------
  1386. C B) Traitement des anciennes relations !
  1387. C---------------------------------------
  1388. C A ce stade, IPT8 contient les relations crees lors du dernier
  1389. C raffinement (noeud de relation + noeuds formant la relation) et les
  1390. C anciennes relations a reporter (noeud de relation seulement)
  1391.  
  1392. C 1/ Recherche de l'existence d'anciennes relations a reporter
  1393. NB10=0
  1394. MLONG=IPT8.NUM(/1)
  1395. DO 775 MM=1,IPT8.NUM(/2)
  1396. IF (IPT8.NUM(MLONG,MM).EQ.0) NB10=NB10+1
  1397. 775 CONTINUE
  1398.  
  1399. IF (NB10.NE.0) THEN
  1400. C
  1401. C 2/ Creation de IPT10
  1402. c MCOL=0
  1403. c NBNN=9
  1404. c NBELEM=NB10
  1405. c NBREF=0
  1406. c NBSOUS=0
  1407. c SEGINI IPT10
  1408. c IPT10.ITYPEL=48
  1409. C
  1410. C 3/ Recherche des anciennes relations a reporter
  1411. C Cette partie permet de compléter IPT8 avec les noeuds formant les
  1412. C relations, pour les anciennes relations à reporter
  1413. DO 628 MM11=1,IPT8.NUM(/2)
  1414. IF (IPT8.NUM(MLONG,MM11).NE.0) GOTO 628
  1415. INRELA=IPT8.NUM(1,MM11)
  1416. c On active pous le cas ou IPT6=IPT2.
  1417. SEGACT IPT6
  1418. DO 625 KOU=1,MAX(1,IPT6.LISOUS(/1))
  1419. IF (IPT6.LISOUS(/1).EQ.0) THEN
  1420. IPT2=IPT6
  1421. ELSE
  1422. IPT2=IPT6.LISOUS(KOU)
  1423. ENDIF
  1424. SEGACT IPT2
  1425. IF (IPT2.ITYPEL.NE.48) GOTO 625
  1426. DO 626 JKG=1,IPT2.NUM(/2)
  1427. INOEGL=IPT2.NUM(1,JKG)
  1428. IF (INOEGL.NE.INRELA) GOTO 626
  1429. MCOL=MCOL+1
  1430. DO 627 MTMP=1,IPT2.NUM(/1)
  1431. IPT8.NUM(MTMP,MM11)=IPT2.NUM(MTMP,JKG)
  1432. IPT8.NUM(MLONG,MM11)= IPT2.ICOLOR(JKG)
  1433. 627 CONTINUE
  1434. GOTO 628
  1435. 626 CONTINUE
  1436. 625 CONTINUE
  1437. 628 CONTINUE
  1438. ENDIF
  1439. c à ce stade le maillage de relation ipt8 est complet
  1440. C
  1441. C---------------------------------------
  1442. C C) Traitement des nouvelles relations !
  1443. C---------------------------------------
  1444. C Ici, IPT8 contient toutes les relations du nouveau maillage
  1445. C On les trie par type de relations
  1446.  
  1447. C 1/ Recherche de l'existence de nouvelles relations
  1448. NB11=0
  1449. NB12=0
  1450. NB13=0
  1451. NB14=0
  1452. NB15=0
  1453. NB16=0
  1454. MLONG=IPT8.NUM(/1)
  1455. DO 776 MM=1,IPT8.NUM(/2)
  1456. IF (IPT8.NUM(MLONG,MM).EQ.1) NB11=NB11+1
  1457. IF (IPT8.NUM(MLONG,MM).EQ.2) NB12=NB12+1
  1458. IF (IPT8.NUM(MLONG,MM).EQ.3) NB13=NB13+1
  1459. IF (IPT8.NUM(MLONG,MM).EQ.4) NB14=NB14+1
  1460. IF (IPT8.NUM(MLONG,MM).EQ.5) NB15=NB15+1
  1461. IF (IPT8.NUM(MLONG,MM).EQ.6) NB16=NB16+1
  1462. 776 CONTINUE
  1463. C
  1464. C 2/ Cas des relations 'aretes a 2 noeuds'
  1465. IF (NB11.NE.0) THEN
  1466. MCONF=MCONF+1
  1467. MCOL=0
  1468. NBNN=3
  1469. NBELEM=NB11
  1470. NBREF=0
  1471. NBSOUS=0
  1472. SEGINI IPT11
  1473. IPT11.ITYPEL=48
  1474. DO 778 MM11=1,IPT8.NUM(/2)
  1475. IF (IPT8.NUM(MLONG,MM11).NE.1) GOTO 778
  1476. MCOL=MCOL+1
  1477. DO 777 MTMP=1,3
  1478. IPT11.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
  1479. IPT11.ICOLOR(MCOL)=1
  1480. 777 CONTINUE
  1481. 778 CONTINUE
  1482. ENDIF
  1483. C
  1484. C 3/ Cas des relations 'aretes a 3 noeuds'
  1485. IF (NB12.NE.0) THEN
  1486. MCONF=MCONF+1
  1487. MCOL=0
  1488. NBNN=4
  1489. NBELEM=NB12
  1490. NBREF=0
  1491. NBSOUS=0
  1492. SEGINI IPT12
  1493. IPT12.ITYPEL=48
  1494. DO 668 MM11=1,IPT8.NUM(/2)
  1495. IF (IPT8.NUM(MLONG,MM11).NE.2) GOTO 668
  1496. MCOL=MCOL+1
  1497. DO 667 MTMP=1,4
  1498. IPT12.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
  1499. IPT12.ICOLOR(MCOL)=2
  1500. 667 CONTINUE
  1501. 668 CONTINUE
  1502. ENDIF
  1503. C
  1504. C 4/ Cas des relations 'face QUA4'
  1505. IF (NB13.NE.0) THEN
  1506. MCONF=MCONF+1
  1507. MCOL=0
  1508. NBNN=5
  1509. NBELEM=NB13
  1510. NBREF=0
  1511. NBSOUS=0
  1512. SEGINI IPT13
  1513. IPT13.ITYPEL=48
  1514. DO 658 MM11=1,IPT8.NUM(/2)
  1515. IF (IPT8.NUM(MLONG,MM11).NE.3) GOTO 658
  1516. MCOL=MCOL+1
  1517. DO 657 MTMP=1,5
  1518. IPT13.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
  1519. IPT13.ICOLOR(MCOL)=3
  1520. 657 CONTINUE
  1521. 658 CONTINUE
  1522. ENDIF
  1523. C
  1524. C 5/ Cas des relations 'face TRI6'
  1525. IF (NB14.NE.0) THEN
  1526. MCONF=MCONF+1
  1527. MCOL=0
  1528. NBNN=7
  1529. NBELEM=NB14
  1530. NBREF=0
  1531. NBSOUS=0
  1532. SEGINI IPT14
  1533. IPT14.ITYPEL=48
  1534. DO 648 MM11=1,IPT8.NUM(/2)
  1535. IF (IPT8.NUM(MLONG,MM11).NE.4) GOTO 648
  1536. MCOL=MCOL+1
  1537. DO 647 MTMP=1,7
  1538. IPT14.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
  1539. IPT14.ICOLOR(MCOL)=4
  1540. 647 CONTINUE
  1541. 648 CONTINUE
  1542. ENDIF
  1543. C
  1544. C 6/ Cas des relations 'face QUA8-1'
  1545. IF (NB15.NE.0) THEN
  1546. MCONF=MCONF+1
  1547. MCOL=0
  1548. NBNN=9
  1549. NBELEM=NB15
  1550. NBREF=0
  1551. NBSOUS=0
  1552. SEGINI IPT15
  1553. IPT15.ITYPEL=48
  1554. DO 638 MM11=1,IPT8.NUM(/2)
  1555. IF (IPT8.NUM(MLONG,MM11).NE.5) GOTO 638
  1556. MCOL=MCOL+1
  1557. DO 637 MTMP=1,9
  1558. IPT15.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
  1559. IPT15.ICOLOR(MCOL)=5
  1560. 637 CONTINUE
  1561. 638 CONTINUE
  1562. ENDIF
  1563. C
  1564. C 7/ Cas des relations 'face QUA8-2'
  1565. IF (NB16.NE.0) THEN
  1566. MCONF=MCONF+1
  1567. MCOL=0
  1568. NBNN=9
  1569. NBELEM=NB16
  1570. NBREF=0
  1571. NBSOUS=0
  1572. SEGINI IPT16
  1573. IPT16.ITYPEL=48
  1574. DO 678 MM11=1,IPT8.NUM(/2)
  1575. IF (IPT8.NUM(MLONG,MM11).NE.6) GOTO 678
  1576. MCOL=MCOL+1
  1577. DO 677 MTMP=1,9
  1578. IPT16.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
  1579. IPT16.ICOLOR(MCOL)=6
  1580. 677 CONTINUE
  1581. 678 CONTINUE
  1582. ENDIF
  1583. C
  1584. ENDIF
  1585. C
  1586. C--------------------------------------
  1587. C D) Creation du maillage de relations !
  1588. C--------------------------------------
  1589. C 1/ Initialisation de MAIREL
  1590. C MAIREL contient tous les maillages de relations, certains étant vides
  1591. C car il n'y a aucune relation du type considéré dans le nveau maillage
  1592. NBMAIL=6
  1593. SEGINI MAIREL
  1594. MAIREL(1)=IPT11
  1595. MAIREL(2)=IPT12
  1596. MAIREL(3)=IPT13
  1597. MAIREL(4)=IPT14
  1598. MAIREL(5)=IPT15
  1599. MAIREL(6)=IPT16
  1600. C
  1601. C 2/ Initialisation de LTYPE2
  1602. C LTYPE2 va contenir tous les sous-objets du nouveau maillage,
  1603. C relations et non relations
  1604. NBLIS2=MCONF+MAX(1,MELEME.LISOUS(/1))
  1605. SEGINI LTYPE2
  1606. C
  1607. C 3/ Renseignement des lisous non 48 dans LTYPE2
  1608. C S'il y a N sous-objets de non-relations dans le nouveau maillage,
  1609. C les N-ièmes premières colonnes de LTYPE2 contiennent ces sous-objets
  1610. DO 685 IZZ=1,MAX(1,MELEME.LISOUS(/1))
  1611. IF (MELEME.LISOUS(/1).EQ.0) THEN
  1612. IPT2=MELEME
  1613. ELSE
  1614. IPT2=MELEME.LISOUS(IZZ)
  1615. SEGACT IPT2
  1616. ENDIF
  1617. LTYPE2(IZZ,1)=IZZ
  1618. LTYPE2(IZZ,2)=IPT2.ITYPEL
  1619. LTYPE2(IZZ,3)=IPT2
  1620. LTYPE2(IZZ,4)=0
  1621. 685 CONTINUE
  1622. C
  1623. C 4/ Renseignement des lisous 48 dans LTYPE2
  1624. C Les colonnes suivantes contiennent les sous-objets de relations,
  1625. C chaque sous-objet correspondant a un type de relations de conformité
  1626. IZZ=NBLIS2-MCONF
  1627. DO 686 IXX=1,6
  1628. IF (MAIREL(IXX).EQ.0) GOTO 686
  1629. IPT2=MAIREL(IXX)
  1630. IZZ=IZZ+1
  1631. LTYPE2(IZZ,1)=IZZ
  1632. LTYPE2(IZZ,2)=IPT2.ITYPEL
  1633. LTYPE2(IZZ,3)=IPT2
  1634. LTYPE2(IZZ,4)=IPT2.ICOLOR(1)
  1635. 686 CONTINUE
  1636. C--------------------------------------
  1637. C E) Renseignement de LTYPE2 interne a
  1638. c la boucle 33 dans LTYPE3 externe a
  1639. c la boucle 33
  1640. C--------------------------------------C
  1641. C
  1642.  
  1643. LTYPE3(INB,1)=LTYPE2
  1644. IF (MCONF.EQ.0) THEN
  1645. LTYPE3(INB,2)=0
  1646. ELSE
  1647. LTYPE3(INB,2)=NBLIS2-MCONF+1
  1648. ENDIF
  1649. c fin de boucle sur les sous domaines d'ipt5
  1650. 33 CONTINUE
  1651. c SEGSUP IPT2
  1652.  
  1653.  
  1654. C Mise a jour des maillages en entrée sortie de la boucle 9
  1655. C IPT5, IPT6, IPT7
  1656. C=======================================================================
  1657. 679 CONTINUE
  1658. C--------------------------------------C
  1659. C A) initialisation
  1660. C--------------------------------------C
  1661.  
  1662. c 1/ initialisaton de IPT7 maillage stockant tout les sous maillage
  1663. c à la fin d'une itération de la boucle 9
  1664. NBS7=0
  1665. DO 331 I3=1,LTYPE3(/1)
  1666. LTYPE2=LTYPE3(I3, 1)
  1667. SEGACT LTYPE2
  1668. NBS7=NBS7+LTYPE2(/1)+1
  1669. 331 CONTINUE
  1670. IF (NBS7.EQ.1) NBS7=0
  1671. NBNN=0
  1672. NBELEM=0
  1673. NBREF=0
  1674. NBSOUS=NBS7
  1675. SEGINI IPT7
  1676.  
  1677. c 2/ initialisaton de IPT6 maillage stockant tout les sous maillage
  1678. c de relation à la fin d'une itération de la boucle 9
  1679. NBS6=0
  1680.  
  1681. DO 332 I3=1,LTYPE3(/1)
  1682. LTYPE2=LTYPE3(I3, 1)
  1683. SEGACT LTYPE2
  1684. IF (LTYPE3(I3,2).NE.0) NBS6=NBS6+(LTYPE2(/1)-LTYPE3(I3,2)+1)
  1685. 332 CONTINUE
  1686. IF (NBS6.EQ.1) NBS6=0
  1687. NBNN=0
  1688. NBELEM=0
  1689. NBREF=0
  1690. NBSOUS=NBS6
  1691. SEGINI IPT6
  1692. c 3/ initialisaton de IPT5 maillage stockant tout les sous maillage
  1693. c non relation à la fin d'une itération de la boucle 9
  1694. NBS5=0
  1695. DO 333 I3=1,LTYPE3(/1)
  1696. LTYPE2=LTYPE3(I3, 1)
  1697. SEGACT LTYPE2
  1698. IF (LTYPE3(I3,2).EQ.0) THEN
  1699. NBS5=NBS5+LTYPE2(/1)
  1700. ELSE
  1701. NBS5=NBS5+(LTYPE3(I3,2)-1)
  1702. ENDIF
  1703. 333 CONTINUE
  1704. IF (NBS5.EQ.1) NBS5=0
  1705. NBNN=0
  1706. NBELEM=0
  1707. NBREF=0
  1708. NBSOUS=NBS5
  1709. SEGINI IPT5
  1710. C--------------------------------------C
  1711. C B) assemblage des nouveaux maillage à
  1712. c partir des sous maillages stockés
  1713. c dans LTYP3
  1714. C--------------------------------------C
  1715. IND5=0
  1716. IND6=0
  1717. IND7=0
  1718. c 1/assemblage des nouveaux maillage non relation à partir des sous
  1719. c maillages stockés dans LTYP3
  1720.  
  1721. DO 35 I3=1,LTYPE3(/1)
  1722. c extraction d'un LTYPE2 de LTYPE3
  1723. LTYPE2=LTYPE3(I3, 1)
  1724. SEGACT LTYPE2
  1725. c NBNN=0
  1726. c NBELEM=0
  1727. c NBREF=0
  1728. c NBSOUS=0
  1729. DO 613 I2=1,MAX(1,(LTYPE3(I3,2)-1))
  1730. c extraction d'un MAILLAGE non relation de LTYPE2
  1731. IPT2=LTYPE2(I2,3)
  1732. SEGACT IPT2
  1733. IND5=IND5+1
  1734. IF (NBS5.EQ.0) THEN
  1735. IPT5=IPT2
  1736. c write (*,*) 'IPT5.' , IPT5
  1737. c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
  1738. ELSE
  1739. IPT5.LISOUS(IND5)=IPT2
  1740. c write (*,*) 'IPT5.LISOUS(',IND5,')' , IPT5.LISOUS(IND5)
  1741. c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
  1742. ENDIF
  1743. IND7=IND7+1
  1744. IF (NBS7.EQ.0) THEN
  1745. IPT7=IPT2
  1746. ELSE
  1747. IPT7.LISOUS(IND7)=IPT2
  1748. c write (*,*) 'IPT7.LISOUS(',IND7,')' , IPT7.LISOUS(IND7)
  1749. c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
  1750. ENDIF
  1751. 613 CONTINUE
  1752. 35 CONTINUE
  1753.  
  1754. c 2/assemblage des nouveaux maillage de relation à partir des derniers
  1755. c sous maillages stockés dans LTYPE3
  1756.  
  1757. DO 36 I3=1,LTYPE3(/1)
  1758. c extraction d'un LTYPE2 de LTYPE3
  1759. c write (*,*) 'I3', I3
  1760. c write (*,*) 'LTYPE3(I3,2)', LTYPE3(I3,2)
  1761. IF (LTYPE3(I3,2).EQ.0) GOTO 36
  1762. SEGACT LTYPE2
  1763. LTYPE2=LTYPE3(I3, 1)
  1764. NBNN=0
  1765. NBELEM=0
  1766. NBREF=0
  1767. NBSOUS=0
  1768. DO 614 I2=LTYPE3(I3,2), LTYPE2(/1)
  1769. c extraction d'un MAILLAGE de relation de LTYPE2
  1770. c write (*,*) 'I2', I2
  1771. SEGINI IPT2
  1772. IPT2=LTYPE2(I2,3)
  1773. SEGACT IPT2
  1774. c write (*,*) 'IPT2: Type', IPT2.ITYPEL,'nbelem', IPT2.NUM(/2)
  1775. c write (*,*) 'coul IPT2' , IPT2.ICOLOR(1)
  1776. IND6=IND6+1
  1777. IF (NBS6.EQ.0) THEN
  1778. IPT6=IPT2
  1779. c write (*,*) 'IPT6' , IPT6
  1780. c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
  1781. ELSE
  1782. IPT6.LISOUS(IND6)=IPT2
  1783.  
  1784. c write (*,*) 'IPT6.LISOUS(',IND6,')' , IPT6.LISOUS(IND6)
  1785. c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
  1786. ENDIF
  1787. 614 CONTINUE
  1788. 36 CONTINUE
  1789.  
  1790. c 3/ Concatenation des sous maillages de relation de même couleur
  1791.  
  1792. c GG 11 juin 2018 : attention il peut encore y avoir des problemes a la
  1793. c jonctions entre deux sous zones dans le cas ou on raffine plusieurs
  1794. c types d'element differents (passer cette etape de concatenation dans
  1795. c la boucle 33 ?
  1796.  
  1797. IF (NBS6.EQ.0) then
  1798. IND7=IND7+1
  1799. IPT7.LISOUS(IND7) = IPT6
  1800. c write (*,*) 'IPT7.LISOUS(',IND7,')' , IPT7.LISOUS(IND7)
  1801. c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
  1802.  
  1803. goto 370
  1804. ENDIF
  1805. NBS60 = 0
  1806. SEGACT IPT6
  1807.  
  1808. c boucle sur les couleurs
  1809. NBELEM = 0
  1810. NBNN= 0
  1811. NBREF= 0
  1812. NBSOUS=6
  1813. SEGINI IPT4
  1814. DO 37 ICOUL = 1, 6
  1815. c boucle sur les sous maillages de ipt6
  1816. NCOl = 0
  1817. DO 38 I6 = 1, NBS6
  1818. IPT3 = IPT6.LISOUS(I6)
  1819. SEGACT IPT3
  1820. IF ((IPT3.icolor(1)).EQ.ICOUL) then
  1821. NCOL = NCOL+1
  1822. IF (NCOL.EQ.1) THEN
  1823. SEGINI IPT2
  1824. IPT2 = IPT6.LISOUS(I6)
  1825.  
  1826. ENDIF
  1827. IF (NCOL.GE.2) THEN
  1828.  
  1829. c IPT3 = IPT6.LISOUS(I6)
  1830. c SEGACT IPT3
  1831. c Plusieurs sous maillages de même couleur on concatene IPT2 et IPT3
  1832. NBELEM = IPT2.NUM(/2)
  1833. NBNN= IPT2.NUM(/1)
  1834. NBREF=0
  1835. NBSOUS=0
  1836. NBELE0 = NBELEM
  1837. DO 381 IE3 = 1, IPT3.NUM(/2)
  1838. DO 382 IE2 = 1, NBELE0
  1839. IF ((IPT3.NUM(1,IE3)).EQ.(IPT2.NUM(1,IE2))) GOTO 381
  1840. 382 CONTINUE
  1841. NBELEM = NBELEM +1
  1842. SEGADJ IPT2
  1843. DO IN2 = 1, NBNN
  1844. IPT2.NUM(IN2 , NBELEM) = IPT3.NUM(IN2, IE3)
  1845. IPT2.icolor(NBELEM) = IPT3.icolor(IE3)
  1846. ENDDO
  1847. 381 CONTINUE
  1848. ENDIF
  1849. ENDIF
  1850. 38 CONTINUE
  1851. IF (NCOL.GE.1) THEN
  1852. IPT4.LISOUS(ICOUL) = IPT2
  1853. NBS60= NBS60+1
  1854. ENDIF
  1855.  
  1856. 37 CONTINUE
  1857.  
  1858. c mise a jour de ipt6 et ipt7
  1859. IF (NBS60.EQ.1) then
  1860. IPT6 = IPT2
  1861. segact ipt6
  1862.  
  1863. c write (*,*) 'IPT6' , IPT6
  1864. c write (*,*) 'type', IPT6.ITYPEL, 'Nbel' , IPT6.num(/2)
  1865. IND7=IND7+1
  1866. IPT7.LISOUS(IND7) = IPT2
  1867. c write (*,*) 'IPT7.LISOUS(',IND7,')' , IPT7.LISOUS(IND7)
  1868. c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
  1869. ELSE
  1870. SEGACT IPT4
  1871. NBELEM = 0
  1872. NBNN = 0
  1873. NBREF = 0
  1874. NBSOUS = NBS60
  1875. SEGINI IPT6
  1876. IS6=0
  1877. DO ICOUL = 1, 6
  1878. IF ((IPT4.LISOUS(ICOUL)).NE.0) THEN
  1879. IS6 = IS6+1
  1880. IPT6.LISOUS(IS6)= IPT4.LISOUS(ICOUL)
  1881. c write (*,*) 'IPT6.LISOUS(',IND6,')' , IPT6.LISOUS(IND6)
  1882. c write (*,*) 'type', IPT4.ITYPEL, 'Nbel' , IPT4.num(/2)
  1883.  
  1884. IND7=IND7+1
  1885. IPT7.LISOUS(IND7) = IPT4.LISOUS(ICOUL)
  1886. ENDIF
  1887. ENDDO
  1888. ENDIF
  1889. 370 CONTINUE
  1890. NBELEM = 0
  1891. NBNN = 0
  1892. NBREF = 0
  1893. NBSOUS = IND7
  1894.  
  1895. SEGADJ IPT7
  1896.  
  1897. SEGINI LTYPE4
  1898. LTYPE4 = LTYPE3
  1899.  
  1900. c mise a jour du nombre global de relations
  1901. LCONF3 = MCONF
  1902.  
  1903. C=======================================================================
  1904. C Test sur la progression du raffinement
  1905. C=======================================================================
  1906. C
  1907. C CB215821 Retrait du caractere *MOD le cas echeant
  1908. SEGACT IPT5,MCHPO1,MSOUPO,MPOVA1,IPT3,KMILIE
  1909. c SEGSUP ICPR,KARETE,KARPOS,KARET2,ICHPO,IDENSI
  1910. c SEGSUP ICPR,KARETE,KARPOS,KARET2,IDENSI
  1911. SEGSUP KARETE,KARPOS,KARET2
  1912. SEGSUP MCHAM1,MELVA1,MCHAM2,MELVA2,IMODE1
  1913. SEGSUP MMODE1,MCHEL1,MCHEL2
  1914.  
  1915. IF ((IDIM.EQ.3).AND.(NBOUCL.NE.1)) SEGSUP JPLANS,JPLAN2,JPLAN3,
  1916. +JPLCOM,JNOEFA,IWORK2,JFARAF, IPT9
  1917.  
  1918. C S'il y a eu a raffiner des elements lors de l'iteration courante,
  1919. C on passe a l'iteration suivante
  1920. C
  1921. c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1922. c forcage : nb max itération de raff
  1923. c IF (NBOUCL.EQ.3) THEN
  1924. c KARAF2=0
  1925. c KARAF=0
  1926. c ENDIF
  1927. c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1928. IF (KARAF2.NE.0) GOTO 9
  1929.  
  1930. C ### SORTIE TEMPORAIRE DU PROGRAMME ###
  1931. C Si on est a la premiere iteration de RAFF et qu'il n'y a pas
  1932. C d'elements a raffiner, on retourne en sortie le maillage de depart,
  1933. C donné en entree
  1934. IF ((KARAF2.EQ.0).AND.(NBOUCL.EQ.1)) THEN
  1935. WRITE(*,*) 'il n est pas necessaire de raffiner le maillage'
  1936. MELEME=IPT1
  1937. GOTO 998
  1938.  
  1939. ENDIF
  1940.  
  1941. C=======================================================================
  1942. C Mises a jour diverses
  1943. C=======================================================================
  1944. C ** On affecte une densite a chaque noeud que l'on a cree
  1945. c DO 37 I=NPTINI+1,XCOOR(/1)/(IDIM+1)
  1946. c XCOOR(I*(IDIM+1))=DENSI(ICPR(I,1))
  1947. c 37 CONTINUE
  1948.  
  1949. C=======================================================================
  1950. C Fin du programme
  1951. C=======================================================================
  1952.  
  1953. MELEME=IPT7
  1954.  
  1955. SEGSUP XDEN
  1956. 998 CONTINUE
  1957.  
  1958. CALL ACTOBJ('MAILLAGE',MELEME,1)
  1959. CALL ECROBJ('MAILLAGE',MELEME)
  1960.  
  1961. END
  1962.  
  1963.  
  1964.  
  1965.  

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