Télécharger raff.eso

Retour à la liste

Numérotation des lignes :

  1. C RAFF SOURCE CB215821 19/08/20 21:21:10 10287
  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. -INC CCOPTIO
  15. -INC SMCOORD
  16. -INC CCGEOME
  17. -INC SMELEME
  18. -INC SMCHPOI
  19. -INC SMMODEL
  20. -INC SMCHAML
  21. c-INC CCREEL
  22. C
  23. INTEGER MCOL,NCOL
  24. C--------------------------------
  25. C B) VARIABLES RELATIVES A LA 2D !
  26. C--------------------------------
  27. SEGMENT KARETE(NBNDS,NCOL)
  28. SEGMENT KARET2(NBNDS,NCOL)
  29. SEGMENT KMILIE(NBNDS,NCOL)
  30. SEGMENT KARPOS(NBNDS)
  31. C
  32. C--------------------------------
  33. C C) VARIABLES RELATIVES A LA 3D !
  34. C--------------------------------
  35. SEGMENT IRELA(NBNDS)
  36. SEGMENT JPLANS(JPLA1,JPLA2)
  37. SEGMENT JPLAN2(JPLA1,JPLA2)
  38. SEGMENT JPLAN3(JPLA1,JPLA2)
  39. SEGMENT JPLCOM(JPLA1)
  40. SEGMENT JNOEFA(JNBFA,5)
  41. SEGMENT IWORK2(JNBSOM)
  42. SEGMENT JFARAF(JNBFA,LLLL)
  43. SEGMENT IWORK1(IJKLMN)
  44. C--------------------------
  45. C A) VARIABLES 'GENERALES' !
  46. C--------------------------
  47. SEGMENT ICPR((XCOOR(/1)/(IDIM+1)),2)
  48. SEGMENT LTYPE(NBLIS,2)
  49. SEGMENT LTYPE2(NBLIS2,4)
  50. SEGMENT LTYPE3(NBLIS3,2)
  51. SEGMENT LTYPE4(NBLIS3,2)
  52. c SEGMENT ICHPO
  53. c REAL*8 XR(NBPTCH),YR(NBPTCH),ZR(NBPTCH),DEN(NBPTCH)
  54. c ENDSEGMENT
  55. SEGMENT IDENSI
  56. REAL*8 DENSI(NBNDS)
  57. ENDSEGMENT
  58. SEGMENT XDEN((XCOOR(/1)/(IDIM+1)))
  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. CALL LIRENT(KR1,1,IRETOU)
  84. IF (IERR.NE.0) RETURN
  85. CALL RAFLRE(ILRE1,KR1,ILRE2)
  86. IF (IERR.NE.0) RETURN
  87. IF (ILRE2.NE.0) THEN
  88. CALL ECROBJ('LISTREEL',ILRE2)
  89. ELSE
  90. CALL ERREUR(5)
  91. ENDIF
  92. RETURN
  93. ENDIF
  94. C
  95. C SYNTAXE 1 : cas du MAILLAGE
  96. C ---------------------------
  97. C On lit le maillage a raffiner et le champoint de densite a affecter
  98. c ipt1 : maillage initial à raffiner
  99. c mchpo1 : chapoint de densité voulu(que devrai être définit sur ipt1)
  100. CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU)
  101. CALL ACTOBJ('MAILLAGE',IPT1,1)
  102. IF (IERR.NE.0) RETURN
  103. CALL LIROBJ('CHPOINT',MCHPO1,1,IRETOU)
  104. CALL ACTOBJ('CHPOINT',MCHPO1,1)
  105. IF (IERR.NE.0) RETURN
  106. C
  107. C=======================================================================
  108. C Initialisations - Declarations
  109. C=======================================================================
  110. NCOL=0
  111. MCOL=0
  112. NBOUCL=0
  113. NPTINI=(XCOOR(/1)/(IDIM+1))
  114. SEGACT IPT1
  115.  
  116. IJKLMN=4
  117. SEGINI IWORK1
  118. IWORK1(1)=4
  119. IWORK1(2)=8
  120. IWORK1(3)=6
  121. IWORK1(4)=10
  122.  
  123. SEGINI,XDEN
  124.  
  125. C
  126. C=======================================================================
  127. C Tests sur les donnees
  128. C=======================================================================
  129. C On teste pour savoir si le maillage IPT1 (maillage 1) a un type
  130. C d'elements implemante.
  131. c TRI3 TRI6 QUA4 QUA8 CUB8 CUB26 PRY6 PR15 TET4 TET10 SURE
  132. C les pyramides 5 et 13 ont été débranchés
  133. c création du tableau
  134. NBLIS=MAX(1,IPT1.LISOUS(/1))
  135. NSURE=0
  136. SEGINI LTYPE
  137. DO 1 INB=1,MAX(1,IPT1.LISOUS(/1))
  138. IF (IPT1.LISOUS(/1).EQ.0) THEN
  139. IPT2=IPT1
  140. ELSE
  141. IPT2=IPT1.LISOUS(INB)
  142. SEGACT IPT2
  143. ENDIF
  144. NT=IPT2.ITYPEL
  145. IF (NT.EQ.48) NSURE=NSURE+1
  146. LTYPE(INB,1)=INB
  147. LTYPE(INB,2)=NT
  148. IF (NT.EQ.16) THEN
  149. write (*,*) 'Attention le raffinement PRI6 a debrancher !'
  150. ENDIF
  151. IF (NT.EQ.23) THEN
  152. write (*,*) 'Attention le raffinement TET4 a debrancher !'
  153. ENDIF
  154. IF ((NT.NE.4).AND.(NT.NE.6).AND.(NT.NE.8).AND.(NT.NE.10).AND.
  155. + (NT.NE.14).AND.(NT.NE.15).AND.(NT.NE.16).AND.(NT.NE.17).AND.
  156. + (NT.NE.23).AND.(NT.NE.24).AND.(NT.NE.48)) THEN
  157. WRITE(*,*) 'Erreur, type d element non pris en compte'
  158. RETURN
  159. ENDIF
  160. c write(*,*) 'type d elem du maillage de depart', LTYPE(INB,1),
  161. c + LTYPE(INB,2)
  162.  
  163. C IF ((NT.NE.4).AND.(NT.NE.6).AND.(NT.NE.8).AND.(NT.NE.10).AND.
  164. C + (NT.NE.14).AND.(NT.NE.15).AND.(NT.NE.16).AND.(NT.NE.17).AND.
  165. C + (NT.NE.23).AND.(NT.NE.24).AND.(NT.NE.48).AND.(NT.NE.25).AND.
  166. C + (NT.NE.26)) RETURN
  167. IF ((IPT1.LISOUS(/1).EQ.0).AND.(NT.EQ.48)) RETURN
  168. 1 CONTINUE
  169.  
  170. C
  171. C=======================================================================
  172. C Separation du maillage initial
  173. C=======================================================================
  174. C Si le maillage initial comporte des lisous de type 48, on le separe
  175. C en deux maillages dont un ne sera compose que des lisous de type 48.
  176. C Maillage d'elements de type 48 : IPT6
  177. C Maillage sans elements de type 48 : IPT5
  178. C IPT1 reste le maillage initial complet
  179. LCONF2=0
  180. IF (NSURE.EQ.0) GOTO 7
  181. LCONF2=1
  182. 7 CONTINUE
  183. C-------------------------------
  184. C A) Creation du maillage IPT6 !
  185. C-------------------------------
  186. C A.1/** Cas ou il n'y a qu'un lisous de type 48
  187. IF (NSURE.EQ.1) THEN
  188. DO 2 K=1,LTYPE(/1)
  189. IF (LTYPE(K,2).EQ.48) IPT6=IPT1.LISOUS(LTYPE(K,1))
  190. 2 CONTINUE
  191. ENDIF
  192. C
  193. C A.2/** Cas ou il y a plusieurs lisous de type 48
  194. IF (NSURE.GT.1) THEN
  195. NBREF=0
  196. NBSOUS=NSURE
  197. NBNN=0
  198. NBELEM=0
  199. SEGINI IPT6
  200. KLISOU=0
  201. DO 3 K=1,LTYPE(/1)
  202. IF (LTYPE(K,2).EQ.48) THEN
  203. KLISOU=KLISOU+1
  204. IPT6.LISOUS(KLISOU)=IPT1.LISOUS(LTYPE(K,1))
  205. ENDIF
  206. 3 CONTINUE
  207. ENDIF
  208. IF (NSURE.EQ.0) THEN
  209. NBREF=0
  210. NBSOUS=NSURE
  211. NBNN=0
  212. NBELEM=0
  213. SEGINI IPT6
  214. ENDIF
  215. C
  216. C------------------------------
  217. C B) Creation de IPT5 !
  218. C------------------------------
  219. C
  220. NBSOUS=LTYPE(/1)-NSURE
  221. NBREF=0
  222. NBNN=0
  223. NBELEM=0
  224. C B.0/** Cas ou il n'y a qu'un lisous en tout
  225. IF ((NBSOUS.EQ.1).AND.(NSURE.EQ.0)) THEN
  226. NBSOUS=0
  227. SEGINI IPT5
  228. IPT5=IPT1
  229. ENDIF
  230. C B.1/** Cas ou il n'y a qu'un lisous de type different de 48
  231. IF ((NBSOUS.EQ.1).AND.(NSURE.NE.0)) THEN
  232. SEGINI IPT5
  233. DO 4 K=1,LTYPE(/1)
  234. IF (LTYPE(K,2).NE.48) IPT5=IPT1.LISOUS(LTYPE(K,1))
  235. 4 CONTINUE
  236. ENDIF
  237. C
  238. C B.2/** Cas ou il y a plusieurs lisous de types differents de 48
  239. IF (NBSOUS.GT.1) THEN
  240. SEGINI IPT5
  241. KLISOU=0
  242. DO 5 K=1,LTYPE(/1)
  243. IF (LTYPE(K,2).NE.48) THEN
  244. KLISOU=KLISOU+1
  245. IPT5.LISOUS(KLISOU)=IPT1.LISOUS(LTYPE(K,1))
  246. ENDIF
  247. 5 CONTINUE
  248. ENDIF
  249.  
  250. c compteur de relaton globale
  251. LCONF3=LCONF2
  252. C=======================================================================
  253. C Debut de la boucle sur les itérations de Raffinement
  254. C=======================================================================
  255. 9 CONTINUE
  256.  
  257. c write(*,*) 'RAFF: ITERATION boucle 9',NBOUCL
  258. c write(*,*) 'ipt5 : ', ipt5
  259. c write(*,*) 'ipt6 : ', ipt6
  260. c write(*,*) 'lconf3 : ', lconf3
  261. C
  262. C=======================================================================
  263. C Renumerotation locale des noeuds du maillage
  264. C=======================================================================
  265. C ICPR(I,1)=J : Le Ieme noeud global est le Jeme local
  266. C ICPR(I,2)=J : Le Ieme noeud global appartient a J elements
  267.  
  268. SEGACT IPT5
  269.  
  270. NBOUCL=NBOUCL+1
  271. SEGINI ICPR
  272. NBNDS=0
  273. NCOL=0
  274. DO 11 INB=1,MAX(1,IPT5.LISOUS(/1))
  275. IF (IPT5.LISOUS(/1).EQ.0) THEN
  276. IPT2=IPT5
  277. ELSE
  278. IPT2=IPT5.LISOUS(INB)
  279. SEGACT IPT2
  280. ENDIF
  281. IF (IPT2.ITYPEL.EQ.48) GOTO 11
  282. DO 10 J=1,IPT2.NUM(/2)
  283. DO 10 I=1,IPT2.NUM(/1)
  284. ITMP=IPT2.NUM(I,J)
  285. IF (ICPR(ITMP,1).NE.0) THEN
  286. ICPR(ITMP,2)=ICPR(ITMP,2)+1
  287. NCOL=MAX(NCOL,ICPR(ITMP,2))
  288. ELSE
  289. NBNDS=NBNDS+1
  290. ICPR(ITMP,1)=NBNDS
  291. ICPR(ITMP,2)=1
  292. NCOL=MAX(NCOL,ICPR(ITMP,2))
  293. ENDIF
  294. 10 CONTINUE
  295. 11 CONTINUE
  296. IF (IDIM.EQ.2) NCOL=NCOL*2
  297. IF (IDIM.EQ.3) NCOL=NCOL*4
  298. C
  299. C
  300. C
  301. C
  302. C=======================================================================
  303. C Prise en compte des noeuds de relations
  304. C=======================================================================
  305. C IRELA(I,1) = 1 si le noeud I (en numerotation locale) est un noeud
  306. C de relation ("hanging node")
  307. C IRELA(I,1) = 0 sinon
  308. SEGINI IRELA
  309. IF (LCONF3.EQ.0) GOTO 1200
  310. SEGACT IPT6
  311. DO 1218 LAA=1,MAX(1,IPT6.LISOUS(/1))
  312. IF (IPT6.LISOUS(/1).EQ.0) THEN
  313. IPT2=IPT6
  314. ELSE
  315. IPT2=IPT6.LISOUS(LAA)
  316. ENDIF
  317. SEGACT IPT2
  318. IF (IPT2.ITYPEL.NE.48) GOTO 1218
  319. DO 1220 J=1,IPT2.NUM(/2)
  320. IRELA(ICPR(IPT2.NUM(1,J),1))=1
  321. 1220 CONTINUE
  322. 1218 CONTINUE
  323. 1200 CONTINUE
  324.  
  325. C
  326. C=======================================================================
  327. C Prise en compte du champoint de densite
  328. C=======================================================================
  329. C On affecte une densite a chacun des noeuds du maillage. On calcule
  330. C cette densite a partir des densites existantes dans le champoint, a
  331. C l'aide de la fonction RAFDEN.
  332. C Remarque : on ne prend pas en compte les densites qui pourraient
  333. C exister dans le tableau XCOOR.
  334. C
  335.  
  336. C------------------------------------------
  337. C A) Recuperation des valeurs du champoint !
  338. C------------------------------------------
  339. C XDEN(I,1) = densite affectee au noeud I (en numerotation globale)
  340. IF (NBOUCL.EQ.1) THEN
  341. SEGACT MCHPO1
  342. IF (MCHPO1.IPCHP(/1).NE.1) RETURN
  343. MSOUPO=MCHPO1.IPCHP(1)
  344. SEGACT MSOUPO
  345. IPT3=IGEOC
  346. MPOVA1=IPOVAL
  347. SEGACT IPT3
  348. SEGACT MPOVA1
  349. IF (MPOVA1.VPOCHA(/2).NE.1) RETURN
  350. NBPTCH=IPT3.NUM(/2)
  351. c SEGINI ICHPO
  352. DO 13 I=1,NBPTCH
  353. J=IPT3.NUM(1,I)
  354. c XR(I)=XCOOR((J-1)*(IDIM+1)+1)
  355. c YR(I)=XCOOR((J-1)*(IDIM+1)+2)
  356. c IF (IDIM.EQ.3) ZR(I)=XCOOR((J-1)*(IDIM+1)+3)
  357. c DEN(I)=MPOVA1.VPOCHA(I,1)
  358. c if(J.eq.115.or.J.eq.321.or.
  359. c . J.eq.6261.or.J.eq.5815) then
  360. C write(*,*) 'RAFF: #',J,'DEN(',I,'/',NBPTCH,')=',DEN(I)
  361. c endif
  362. XDEN(J)=MPOVA1.VPOCHA(I,1)
  363. 13 CONTINUE
  364. C Rem ;seulement pour les noeud "peres",
  365. C pour les "fils" raff2d remplit XDEN
  366. ENDIF
  367. C
  368. C--------------------------------------------------------------
  369. C B) Affectation d'une densite a chacun des noeuds du maillage !
  370. C--------------------------------------------------------------
  371. C IDENSI(I,1) = densite affectee au noeud I (en numerotation locale)
  372. SEGINI IDENSI
  373. DO 15 INB=1,MAX(1,IPT5.LISOUS(/1))
  374. IF (IPT5.LISOUS(/1).EQ.0) THEN
  375. IPT2=IPT5
  376. ELSE
  377. IPT2=IPT5.LISOUS(INB)
  378. SEGACT IPT2
  379. ENDIF
  380. c write(*,*) 'ipt2 : ', ipt2, 'type', IPT2.ITYPEL
  381. c write(*,*) IPT2.NUM(/2), IPT2.NUM(/1)
  382. c IF (IPT2.ITYPEL.EQ.48) GOTO 15
  383. DO 14 J=1,IPT2.NUM(/2)
  384. DO 14 I=1,IPT2.NUM(/1)
  385. NGLOB=IPT2.NUM(I,J)
  386. NLOC=ICPR(NGLOB,1)
  387. c XPT=XCOOR((NGLOB-1)*(IDIM+1)+1)
  388. c YPT=XCOOR((NGLOB-1)*(IDIM+1)+2)
  389. c ZPT=0
  390. c IF (IDIM.EQ.3) ZPT=XCOOR((NGLOB-1)*(IDIM+1)+3)
  391. c CALL RAFDEN(ICHPO,XPT,YPT,ZPT,IDIM,DENPT)
  392. c DENSI(NLOC)=DENPT
  393. DENSI(NLOC)=XDEN(NGLOB)
  394. c WRITE (*,*) 'densite au pt', NGLOB, XDEN(NGLOB), XPETI
  395. IF (XDEN(NGLOB).LT.XPETI) THEN
  396. WRITE (*,*) 'Champ de densité en entrée nul au point',
  397. + NGLOB, 'impossible de raffiner'
  398. RETURN
  399. ENDIF
  400. 14 CONTINUE
  401. 15 CONTINUE
  402. C
  403. C=======================================================================
  404. C Remplissage du tableau KARETE
  405. C=======================================================================
  406. C On remplit le tableau KARETE avec les donnees provenant de l'analyse
  407. C du maillage.
  408. C KARETE(i,j)=k : le ieme noeud local forme une arete
  409. C avec le keme noeud local (k>i)
  410. C KARPOS(i)=k : compteur de position (k=nb d'aretes touchant i)
  411. C KARET2(i,j)=k : l'arete i-j appartient a k elements
  412. JPLA1=NBNDS
  413. JPLA2=NCOL
  414. SEGINI KARETE,KARPOS,KARET2
  415. NBARET=0
  416. NCOL=1
  417. DO 18 INB=1,MAX(1,IPT5.LISOUS(/1))
  418. IF (IPT5.LISOUS(/1).EQ.0) THEN
  419. IPT2=IPT5
  420. ELSE
  421. IPT2=IPT5.LISOUS(INB)
  422. SEGACT IPT2
  423. ENDIF
  424. c IF (IPT2.ITYPEL.EQ.48) GOTO 18
  425. ILPL=LPL(IPT2.ITYPEL)
  426. ILPT=LPT(IPT2.ITYPEL)
  427. DO 17 J=1,IPT2.NUM(/2)
  428. DO 17 K=1,ILPL*2-1,2
  429. NPTA=IPT2.NUM(KSEGM(ILPT+K-1),J)
  430. NPTB=IPT2.NUM(KSEGM(ILPT+K),J)
  431. NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
  432. NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
  433. NEXIST=0
  434. DO 16 I=1,MAX(1,KARPOS(NMIN))
  435. IF (KARETE(NMIN,I).EQ.NMAX) NEXIST=I
  436. 16 CONTINUE
  437. c création d'une nouvelle arrète
  438. IF (NEXIST.EQ.0) THEN
  439. IF ((KARPOS(NMIN)+1).GT.NCOL) THEN
  440. NCOL=KARPOS(NMIN)+1
  441. SEGADJ KARETE, KARET2
  442. ENDIF
  443. KARPOS(NMIN)=KARPOS(NMIN)+1
  444. KARETE(NMIN,KARPOS(NMIN))=NMAX
  445. KARET2(NMIN,KARPOS(NMIN))=1
  446. NBARET=NBARET+1
  447. c NCOL=MAX(NCOL,KARPOS(NMIN))
  448. ENDIF
  449. IF (NEXIST.NE.0) THEN
  450. KARET2(NMIN,NEXIST)=KARET2(NMIN,NEXIST)+1
  451. ENDIF
  452. 17 CONTINUE
  453. 18 CONTINUE
  454. c SEGADJ KARETE,KARET2
  455. SEGDES KARETE,KARPOS,KARET2
  456. C
  457. C=======================================================================
  458. C Remplissage des tableaux JPLANS et JNOEFA
  459. C=======================================================================
  460. C On remplit les tableaux de description des faces dans le cas 3D.
  461. C JPLANS(i,j)=k : la jeme face dont le plus petit sommet en
  462. C numerotation locale est le sommet i est la face k
  463. C JNOEFA(i,1-3)=j-k-l : Les 3 plus petits sommets en numerotation
  464. C locale de la face i sont les noeuds j,k et l
  465. C JNOEFA(i,4-5)=j-k : La face identifiee par les 3 sommets precedents
  466. C appartient aux lisous de type j et k
  467. C JPLCOM(i)=k : compteur de position (k=nb de faces touchant i)
  468. IF (IDIM.NE.3) GOTO 40
  469. C----------------------------------------------------
  470. C A) Comptage du nombre d'elements total du maillage !
  471. C----------------------------------------------------
  472. NELTOT=0
  473. DO 50 INB=1,MAX(1,IPT5.LISOUS(/1))
  474. IF (IPT5.LISOUS(/1).EQ.0) THEN
  475. IPT2=IPT5
  476. ELSE
  477. IPT2=IPT5.LISOUS(INB)
  478. SEGACT IPT2
  479. ENDIF
  480. NELTOT=NELTOT+IPT2.NUM(/2)
  481. 50 CONTINUE
  482. C
  483. C-------------------------------------------------
  484. C B) Initialisation des tableaux JPLANS et JNOEFA !
  485. C-------------------------------------------------
  486. JNBFA=NELTOT*6
  487. JMAXPL=0
  488. JNUMFA=0
  489. SEGINI JPLANS,JPLCOM,JNOEFA
  490. C
  491. C-----------------------------------------------
  492. C C) Remplissage des tableaux JPLANS et JNOEFA !
  493. C-----------------------------------------------
  494. DO 51 INB=1,MAX(1,IPT5.LISOUS(/1))
  495. IF (IPT5.LISOUS(/1).EQ.0) THEN
  496. IPT2=IPT5
  497. ELSE
  498. IPT2=IPT5.LISOUS(INB)
  499. SEGACT IPT2
  500. ENDIF
  501. C IF (IPT2.ITYPEL.EQ.48) GOTO 51
  502. DO 52 KELEM=1,IPT2.NUM(/2)
  503. JTYPEL=IPT2.ITYPEL
  504. JLTEL1=LTEL(1,JTYPEL)
  505. JLTEL2=LTEL(2,JTYPEL)
  506. DO 53 IBB=1,JLTEL1
  507. JLDEL1=LDEL(1,JLTEL2-1+IBB)
  508. JTYFAC=IWORK1(JLDEL1)
  509. JLDEL2=LDEL(2,JLTEL2-1+IBB)
  510. JNBSOM=NBSOM(JTYFAC)
  511. JSPOS=NSPOS(JTYFAC)
  512. SEGINI IWORK2
  513. DO 54 IAA=1,JNBSOM
  514. NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),KELEM)
  515. IWORK2(IAA)=NGLOBA
  516. 54 CONTINUE
  517. C
  518. C----------------------------------------------------
  519. C D) Tri des sommets de la face par ordre croissant !
  520. C----------------------------------------------------
  521. C On va ranger les 3 plus petits sommets de la face (en numerotation
  522. C globale) dans les variables NPTA, NPTB, NPTC. On les classe par ordre
  523. C croissant.
  524. C
  525. C 1/ Classement des 3 sommets par ordre croissant (num globale)
  526. NPTA=(XCOOR(/1)/(IDIM+1))+1
  527. NPTB=NPTA+1
  528. NPTC=NPTB+1
  529. DO 55 ICC=1,JNBSOM
  530. IF (IWORK2(ICC).LT.NPTA) THEN
  531. NPTC=NPTB
  532. NPTB=NPTA
  533. NPTA=IWORK2(ICC)
  534. ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  535. NPTC=NPTB
  536. NPTB=IWORK2(ICC)
  537. ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  538. NPTC=IWORK2(ICC)
  539. ENDIF
  540. 55 CONTINUE
  541. C
  542. C 2/ Classement des 3 sommets par ordre croissant (num locale)
  543. NPTA2=ICPR(NPTA,1)
  544. NPTB2=ICPR(NPTB,1)
  545. NPTC2=ICPR(NPTC,1)
  546. IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  547. NPTA=NPTA2
  548. NPTB=MIN(NPTB2,NPTC2)
  549. NPTC=MAX(NPTB2,NPTC2)
  550. ENDIF
  551. IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  552. NPTA=NPTB2
  553. NPTB=MIN(NPTA2,NPTC2)
  554. NPTC=MAX(NPTA2,NPTC2)
  555. ENDIF
  556. IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  557. NPTA=NPTC2
  558. NPTB=MIN(NPTA2,NPTB2)
  559. NPTC=MAX(NPTA2,NPTB2)
  560. ENDIF
  561. C
  562. C----------------------------------------
  563. C E) Renseignement des nouvelles faces !
  564. C----------------------------------------
  565. SEGACT LTYPE
  566. c IHERED=0
  567. c DO 1022 MOI=1,IPT2.NUM(/1)
  568. c IHERED=IHERED+IRELA(ICPR(IPT2.NUM(MOI,KELEM),1))
  569. c 1022 CONTINUE
  570.  
  571. C Si NPTA Ne correspond a aucune face : On cree une face
  572. IF (JPLCOM(NPTA).EQ.0) THEN
  573. JNUMFA=JNUMFA+1
  574. JPLCOM(NPTA)=JPLCOM(NPTA)+1
  575. JPLANS(NPTA,JPLCOM(NPTA))=JNUMFA
  576. JNOEFA(JNUMFA,1)=NPTA
  577. JNOEFA(JNUMFA,2)=NPTB
  578. JNOEFA(JNUMFA,3)=NPTC
  579. JNOEFA(JNUMFA,4)=LTYPE(INB,2)
  580. c IF (IHERED.NE.0) JNOEFA(JNUMFA,5)=LTYPE(INB,2)
  581. C On signale les sous face de faces supportant un sure
  582. C JNOEFA(JNUMFA,5)= 48
  583. IF (LCONF3.NE.0) THEN
  584. SEGACT IPT6
  585. DO 501 LAA1=1,MAX(1,IPT6.LISOUS(/1))
  586. IF (IPT6.LISOUS(/1).EQ.0) THEN
  587. IPT4=IPT6
  588. ELSE
  589. IPT4=IPT6.LISOUS(LAA1)
  590. ENDIF
  591. SEGACT IPT4
  592. IF (IPT4.ITYPEL.NE.48) GOTO 501
  593. IF (IPT4.ICOLOR(1).NE.3) GOTO 501
  594. DO 502,MAA1=1,IPT4.NUM(/2)
  595. NEXI1=0
  596. NEXI2=0
  597. DO 503 ICC1=1,JNBSOM
  598. IF (IWORK2(ICC1).EQ.IPT4.NUM(1, MAA1)) THEN
  599. NEXI1=NEXI1+1
  600. ENDIF
  601. 503 CONTINUE
  602. IF (NEXI1.GT.0) THEN
  603. DO 504 ICC2=1,JNBSOM
  604. DO 504 INN2=2,IPT4.NUM(/1)
  605. IF (IWORK2(ICC2).EQ.IPT4.NUM(INN2,MAA1)) THEN
  606. NEXI2=NEXI2+1
  607. ENDIF
  608. 504 CONTINUE
  609. IF (NEXI2.GT.0) JNOEFA(JNUMFA,5)= 48
  610. ENDIF
  611. 502 CONTINUE
  612. 501 CONTINUE
  613. ENDIF
  614.  
  615.  
  616. C Si NPTA correspond a aumoins une face
  617. ELSE
  618. NEXIST=0
  619. DO 56 IDD=1,JPLCOM(NPTA)
  620. ITEMP1=JPLANS(NPTA,IDD)
  621. NAA=JNOEFA(ITEMP1,1)
  622. NBB=JNOEFA(ITEMP1,2)
  623. NCC=JNOEFA(ITEMP1,3)
  624. C Si la face que actuelle existe déja : on renseigne JNOEFA(NEXIST,5)
  625. IF (NAA.EQ.NPTA.AND.NBB.EQ.NPTB.AND.NCC.EQ.NPTC) THEN
  626. NEXIST=ITEMP1
  627. JNOEFA(NEXIST,5)=LTYPE(INB,2)
  628. ENDIF
  629.  
  630.  
  631. 56 CONTINUE
  632. C Si la face que actuelle n'existe pas : On cree une face
  633. IF (NEXIST.EQ.0) THEN
  634. JNUMFA=JNUMFA+1
  635. JPLCOM(NPTA)=JPLCOM(NPTA)+1
  636. JPLANS(NPTA,JPLCOM(NPTA))=JNUMFA
  637. JNOEFA(JNUMFA,1)=NPTA
  638. JNOEFA(JNUMFA,2)=NPTB
  639. JNOEFA(JNUMFA,3)=NPTC
  640. JNOEFA(JNUMFA,4)=LTYPE(INB,2)
  641. c IF (IHERED.NE.0) JNOEFA(JNUMFA,5)=LTYPE(INB,2)
  642.  
  643. C On signale les sous face de faces supportant un sure
  644. C JNOEFA(JNUMFA,5)= 48
  645. IF (LCONF3.NE.0) THEN
  646. SEGACT IPT6
  647. DO 505 LAA1=1,MAX(1,IPT6.LISOUS(/1))
  648. IF (IPT6.LISOUS(/1).EQ.0) THEN
  649. IPT4=IPT6
  650. ELSE
  651. IPT4=IPT6.LISOUS(LAA1)
  652. ENDIF
  653. SEGACT IPT4
  654. IF (IPT4.ITYPEL.NE.48) GOTO 505
  655. IF (IPT4.ICOLOR(1).NE.3) GOTO 505
  656. DO 506,MAA1=1,IPT4.NUM(/2)
  657. NEXI1=0
  658. NEXI2=0
  659. DO 507 ICC1=1,JNBSOM
  660. IF (IWORK2(ICC1).EQ.IPT4.NUM(1, MAA1)) THEN
  661. NEXI1=NEXI1+1
  662. ENDIF
  663. 507 CONTINUE
  664. IF (NEXI1.GT.0) THEN
  665. DO 508 ICC2=1,JNBSOM
  666. DO 508 INN2=2,IPT4.NUM(/1)
  667. IF (IWORK2(ICC2).EQ.IPT4.NUM(INN2,MAA1)) THEN
  668. NEXI2=NEXI2+1
  669. ENDIF
  670. 508 CONTINUE
  671. IF (NEXI2.GT.0) JNOEFA(JNUMFA,5)= 48
  672. ENDIF
  673. 506 CONTINUE
  674. 505 CONTINUE
  675. ENDIF
  676. ENDIF
  677. ENDIF
  678.  
  679. JMAXPL=MAX(JMAXPL,JPLCOM(NPTA))
  680. 53 CONTINUE
  681. 52 CONTINUE
  682. 51 CONTINUE
  683. JPLA2=JMAXPL
  684. JNBFA=JNUMFA
  685. SEGADJ JPLANS,JPLCOM,JNOEFA
  686. SEGDES LTYPE
  687. C
  688. C=======================================================================
  689. C Pre-remplissage du tableau KMILIE
  690. C=======================================================================
  691. C On remplit le tableau lorsqu'il y a des maillages de relations de
  692. C type 48 (IPT6, LCONF3 non nul)
  693. C KARETE(i,j)=k ET KMILIE(i,j)=n : le noeud global n est situe
  694. C au milieu de l'arete [i,k]
  695. 40 CONTINUE
  696. SEGINI KMILIE
  697. SEGACT KARETE*MOD,KARPOS*MOD,KARET2*MOD
  698. IF (LCONF3.EQ.0) THEN
  699. GOTO 28
  700. ENDIF
  701. 22 CONTINUE
  702. NBRELA=0
  703. SEGACT IPT6
  704. DO 218 LAA=1,MAX(1,IPT6.LISOUS(/1))
  705. IF (IPT6.LISOUS(/1).EQ.0) THEN
  706. IPT2=IPT6
  707. ELSE
  708. IPT2=IPT6.LISOUS(LAA)
  709. ENDIF
  710. SEGACT IPT2
  711. IF (IPT2.ITYPEL.NE.48) GOTO 218
  712. IF ((IPT2.ICOLOR(1).NE.1).AND.(IPT2.ICOLOR(1).NE.2)) GOTO 218
  713. DO 220 J=1,IPT2.NUM(/2)
  714.  
  715. NPTA=IPT2.NUM(2,J)
  716. NPTB=IPT2.NUM(3,J)
  717. NPTC=IPT2.NUM(1,J)
  718. NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
  719. NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
  720. NEXIST=0
  721. DO 219 I=1,MAX(1,KARPOS(NMIN))
  722. IF (KARETE(NMIN,I).EQ.NMAX) NEXIST=I
  723. 219 CONTINUE
  724.  
  725. IF (NEXIST.NE.0) then
  726. KMILIE(NMIN,NEXIST)=NPTC
  727. c WRITE (*,*) 'REPORT ANCIENNE RELA AU NOEUD:', NPTC
  728. c write (*,*) 'karret2(nmin nexist)', KARET2(NMIN,NEXIST)
  729. KARET2(NMIN,NEXIST)=2
  730. ELSE
  731. IF (IDIM.EQ.2) goto 27
  732. IF ((KARPOS(NMIN)+1).GT.NCOL) THEN
  733. NCOL=KARPOS(NMIN)+1
  734. SEGADJ KARETE, KARET2, KMILIE
  735. ENDIF
  736. JPOSI = KARPOS(NMIN)
  737. KARPOS(NMIN)=JPOSI+1
  738. KARETE(NMIN,JPOSI+1)=NMAX
  739. KARET2(NMIN,JPOSI+1)=2
  740. KMILIE(NMIN,JPOSI+1)=NPTC
  741. c WRITE (*,*) 'création de l arrete:', NMIN, NMAX
  742. ENDIF
  743. C
  744. C-----------------------------------------------
  745. C Prise en compte des relations "hereditaires" !
  746. C-----------------------------------------------
  747. C Cette partie sert a traiter le cas ou un element serait raffine
  748. C plusieurs fois alors qu'a sa frontiere les autres elements ne le
  749. C seraient pas. Il y aurait donc plusieurs points sur une arete ayant
  750. C des relations de conformite et plus seulement le noeud milieu.
  751. C
  752. c dans l'exenple ci dessous où les noeuds C et D sont des hanging nodes
  753. c L'arette [C B] n'a pas encore été créée c'est ce qui est fait ici.
  754. c
  755. c | |
  756. c | |
  757. c |A C D |B
  758. c x---------X----X----X
  759. c | | | |
  760. c | | | |
  761. 27 CONTINUE
  762. IF (IDIM.EQ.2) THEN
  763. C 1/ Test sur la premiere arete : NMIN-NPTC
  764. IF (NEXIST.NE.0) THEN
  765. NEXIS2=0
  766. NMIN2=MIN(NMIN,ICPR(NPTC,1))
  767. NMAX2=MAX(NMIN,ICPR(NPTC,1))
  768. DO 24 I=1,MAX(1,KARPOS(NMIN2))
  769. IF (KARETE(NMIN2,I).EQ.NMAX2) NEXIS2=I
  770. 24 CONTINUE
  771. IF (NEXIS2.EQ.0) THEN
  772. NBRELA=NBRELA+1
  773. JNBCOL=KARETE(/2)
  774. JPOSI=KARPOS(NMIN2)
  775. IF (JNBCOL.GT.JPOSI) THEN
  776. KARPOS(NMIN2)=JPOSI+1
  777. KARETE(NMIN2,JPOSI+1)=NMAX2
  778. KARET2(NMIN2,JPOSI+1)=2
  779. ELSEIF (JNBCOL.EQ.JPOSI) THEN
  780. NCOL=JNBCOL+1
  781. SEGADJ KARETE,KMILIE,KARET2
  782. KARETE(NMIN2,JPOSI+1)=NMAX2
  783. KARPOS(NMIN2)=JPOSI+1
  784. KARET2(NMIN2,JPOSI+1)=2
  785. ENDIF
  786. c WRITE (*,*) 'création de l arrete:', NMIN2, NMAX2
  787. ELSE
  788. KARET2(NMIN2,NEXIS2)=2
  789. ENDIF
  790. C
  791. C 2/ Test sur la deuxieme arete : NMAX-NPTC
  792. NEXIS2=0
  793. NMIN2=MIN(NMAX,ICPR(NPTC,1))
  794. NMAX2=MAX(NMAX,ICPR(NPTC,1))
  795. DO 25 I=1,MAX(1,KARPOS(NMIN2))
  796. IF (KARETE(NMIN2,I).EQ.NMAX2) NEXIS2=I
  797. 25 CONTINUE
  798. IF (NEXIS2.EQ.0) THEN
  799. NBRELA=NBRELA+1
  800. JNBCOL=KARETE(/2)
  801. JPOSI=KARPOS(NMIN2)
  802. IF (JNBCOL.GT.JPOSI) THEN
  803. KARPOS(NMIN2)=JPOSI+1
  804. KARETE(NMIN2,JPOSI+1)=NMAX2
  805. KARET2(NMIN2,JPOSI+1)=2
  806. ELSEIF (JNBCOL.EQ.JPOSI) THEN
  807. NCOL=JNBCOL+1
  808. SEGADJ KARETE,KMILIE,KARET2
  809. KARETE(NMIN2,JPOSI+1)=NMAX2
  810. KARPOS(NMIN2)=JPOSI+1
  811. KARET2(NMIN2,JPOSI+1)=2
  812. ENDIF
  813. c WRITE (*,*) 'création de l arrete:', NMIN2, NMAX2
  814. ELSE
  815. KARET2(NMIN2,NEXIS2)=2
  816. ENDIF
  817. ENDIF
  818. ENDIF
  819. 220 CONTINUE
  820. 218 CONTINUE
  821. IF (NBRELA.NE.0) GOTO 22
  822. 28 CONTINUE
  823.  
  824. C
  825. C
  826. C=======================================================================
  827. C Pre-remplissage du tableau JFARAF
  828. C=======================================================================
  829. C On remplit le tableau JFARAF et JPLAN3 lorsqu'il y a des maillages de
  830. C relations de type 48 (IPT6, LCONF3 non nul).
  831. c
  832. C Si JPLANS(i,j)=k et JFARAF(k , (l-1)*2+1)= n : le noeud n (en
  833. C numérotation globale) est le hanging node au centre de la jieme face
  834. C dont le plus petit sommet en numerotation locale est le sommet i.
  835. C
  836. C Si de plus la face en question est un TRI6 et JFARAF(k , l*2)=m :
  837. C le noeud m (en numérotation globale) est le 4ieme noeud de l'element
  838. C 48 coul 4
  839. C
  840. C Si de plus la face en question est un QUA8-1 et JFARAF(k , l*2)=m :
  841. C le noeud m (en numérotation globale) est le 5eme noeud de l'element
  842. c 48 coul 4
  843. C
  844. C si JPLAN3(i,j)=k il y a k relation de compatibilité de faces associé a
  845. C la jieme face dont le plus petit sommet en numerotation locale est le
  846. C sommet i.
  847.  
  848. IF (IDIM.NE.3) GOTO 711
  849. JNBFA=JNOEFA(/1)
  850. LLLL=10
  851. SEGINI JFARAF
  852. SEGINI JPLAN2,JPLAN3
  853. IF (LCONF3.EQ.0) GOTO 711
  854. MPOS25=0
  855. JNBSOM=0
  856. 1287 CONTINUE
  857. NBRELA=0
  858. DO 670 MIE=1,MAX(1,IPT6.LISOUS(/1))
  859. IF (IPT6.LISOUS(/1).EQ.0) THEN
  860. IPT2=IPT6
  861. ELSE
  862. IPT2=IPT6.LISOUS(MIE)
  863. ENDIF
  864. SEGACT IPT2
  865. IF (IPT2.ITYPEL.NE.48) GOTO 670
  866. IF ((IPT2.ICOLOR(1).EQ.1).OR.(IPT2.ICOLOR(1).EQ.2)) GOTO 670
  867. LIGNUM=IPT2.NUM(/1)
  868. IF (IPT2.ICOLOR(1).EQ.4) MPOS25=4
  869. IF (IPT2.ICOLOR(1).EQ.5) MPOS25=5
  870. IF (IPT2.ICOLOR(1).EQ.3) JNBSOM=4
  871. IF (IPT2.ICOLOR(1).EQ.4) JNBSOM=3
  872. IF (IPT2.ICOLOR(1).EQ.5) JNBSOM=4
  873. IF (IPT2.ICOLOR(1).EQ.6) JNBSOM=4
  874. SEGINI IWORK2
  875. C
  876. C-------------------------
  877. C A) Recherche de la face !
  878. C-------------------------
  879. C 1/ Classement des 3 sommets par ordre croissant (num globale)
  880. DO 671,MAA=1,IPT2.NUM(/2)
  881. DO 673 MBB=1,JNBSOM
  882. IWORK2(MBB)=IPT2.NUM(LIGNUM-JNBSOM+MBB,MAA)
  883. 673 CONTINUE
  884. NPTA=(XCOOR(/1)/(IDIM+1))+1
  885. NPTB=NPTA+1
  886. NPTC=NPTB+1
  887. DO 672 ICC=1,JNBSOM
  888. IF (IWORK2(ICC).LT.NPTA) THEN
  889. NPTC=NPTB
  890. NPTB=NPTA
  891. NPTA=IWORK2(ICC)
  892. ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  893. NPTC=NPTB
  894. NPTB=IWORK2(ICC)
  895. ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  896. NPTC=IWORK2(ICC)
  897. ENDIF
  898. 672 CONTINUE
  899. C
  900. C 2/ Classement des 3 sommets par ordre croissant (num locale)
  901. NPTA1=NPTA
  902. NPTB1=NPTB
  903. NPTC1=NPTC
  904. NPTA2=ICPR(NPTA,1)
  905. NPTB2=ICPR(NPTB,1)
  906. NPTC2=ICPR(NPTC,1)
  907. IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  908. NPTA=NPTA2
  909. NPTB=MIN(NPTB2,NPTC2)
  910. NPTC=MAX(NPTB2,NPTC2)
  911. ENDIF
  912. IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  913. NPTA=NPTB2
  914. NPTB=MIN(NPTA2,NPTC2)
  915. NPTC=MAX(NPTA2,NPTC2)
  916. ENDIF
  917. IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  918. NPTA=NPTC2
  919. NPTB=MIN(NPTA2,NPTB2)
  920. NPTC=MAX(NPTA2,NPTB2)
  921. ENDIF
  922. C
  923. C 3/ Recherche du numero de la face
  924. NEXIS2=0
  925. DO 674 IEE=1,JPLCOM(NPTA)
  926. MTMP=JPLANS(NPTA,IEE)
  927. JJ1=JNOEFA(MTMP,1)
  928. JJ2=JNOEFA(MTMP,2)
  929. JJ3=JNOEFA(MTMP,3)
  930. IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
  931. NEXIS2=IEE
  932. ENDIF
  933. 674 CONTINUE
  934. c 3.1/ si la face en question n'existe pas, dans le cas général on
  935. c s'arrète
  936. IF ((NEXIS2.EQ.0).AND.(IPT2.ICOLOR(1).NE.3)) THEN
  937. c WRITE (*,*) 'IPT2.ICOLOR(1)', IPT2.ICOLOR(1)
  938. WRITE (*,*) 'Raffinement impossible, le champ de densite
  939. + presente des variations trop brusques'
  940. WRITE (*,*) 'sure concerne', MAA
  941. WRITE (*,*) 'premier nd du sure concerne', IPT2.NUM(1,MAA)
  942. WRITE (*,*) 'numero du 1er nd de la face concernee', NPTA
  943. KARAF2=0
  944. c GOTO 679
  945. ENDIF
  946. c 3.2/ si la face en question n'existe pas dans le cas d'une relation
  947. c de face QUA4 il suffie de la créer
  948.  
  949. IF ((NEXIS2.EQ.0).AND.(IPT2.ICOLOR(1).EQ.3)) THEN
  950. JPLCOM(NPTA)= JPLCOM(NPTA)+1
  951. JNUMFA=JNUMFA+1
  952. JNUMFS=JNUMFA
  953. JNUMF1=JPLCOM(NPTA)
  954. JMAXPL=MAX(JMAXPL,JPLCOM(NPTA))
  955. c 3.2.1/ ajustement des segments
  956. JPLA2=JMAXPL
  957. JNBFA=JNUMFA
  958. SEGADJ JPLANS
  959. SEGADJ JNOEFA
  960. SEGADJ JPLAN3
  961. SEGADJ JFARAF
  962. C 3.2.1 Renseignement des iformation de la nouvelle face dans JPLANS et
  963. c JNOEFA, initialisation de JPLAN3
  964. JPLANS(NPTA,JNUMF1)=JNUMFA
  965. JPLAN3(NPTA,JNUMF1)=0
  966. JNOEFA(JNUMFA,1)= NPTA
  967. JNOEFA(JNUMFA,2)= NPTB
  968. JNOEFA(JNUMFA,3)= NPTC
  969. JNOEFA(JNUMFA,4)= 48
  970. JNOEFA(JNUMFA,5)= 48
  971. ENDIF
  972.  
  973. IF (NEXIS2.NE.0) THEN
  974. JNUMFS=JPLANS(NPTA,NEXIS2)
  975. ENDIF
  976. C
  977. C 4/ Recherche de la position de la face dans JPLANS
  978. NEXIS2=0
  979. DO 675 MAC=1,JPLCOM(NPTA)
  980. IF (JPLANS(NPTA,MAC).EQ.JNUMFS) NEXIS2=MAC
  981. 675 CONTINUE
  982. C
  983. C-------------------------------------
  984. C B) Renseignement de JFARAF et JPLAN3!
  985. C-------------------------------------
  986.  
  987.  
  988. JPLAN3(NPTA,NEXIS2)=JPLAN3(NPTA,NEXIS2)+1
  989. LTEMPO=(JPLAN3(NPTA,NEXIS2)-1)*2+1
  990. JFARAF(JNUMFS,LTEMPO)=IPT2.NUM(1,MAA)
  991. IF (MPOS25.NE.0) JFARAF(JNUMFS,LTEMPO+1)=IPT2.NUM(MPOS25,MAA)
  992. IF (JNOEFA(JNUMFS,5).EQ.0) JNOEFA(JNUMFS,5)=48
  993.  
  994.  
  995. 671 CONTINUE
  996. 670 CONTINUE
  997. C IF (NBRELA.NE.0) GOTO 1287
  998. 711 CONTINUE
  999. C
  1000. C=======================================================================
  1001. C Calcul du critere et raffinement si necessaire
  1002. C=======================================================================
  1003. C Pour calculer le critere, on a besoin d'un modele et d'un chamelem
  1004. C constant egal a 1.
  1005. C
  1006. C------------------------------------------
  1007. C A) Creation d'un modele simple : MMODE1 !
  1008. C------------------------------------------
  1009. MN3=0
  1010. NFOR=1
  1011. NMAT=2
  1012. NPARMO=0
  1013. NOBMOD=0
  1014. N1=1
  1015. SEGINI MMODE1
  1016. SEGINI IMODE1
  1017. MMODE1.KMODEL(1)=IMODE1
  1018. IMODE1.FORMOD(1)='MECANIQUE'
  1019. IMODE1.MATMOD(1)='ELASTIQUE'
  1020. IMODE1.MATMOD(2)='ISOTROPE'
  1021.  
  1022. C
  1023. KARAF2=0
  1024. C=======================================================================
  1025. C Boucle sur les eventuels sous-domaines de IPT5
  1026. C=======================================================================
  1027.  
  1028. c-----------------------------------------------------------------------
  1029. c initialisation de LTYPE3 qui va contenir les tableau LTYP2
  1030. c correspondant à chaque sous domaine de IPT5. Ces tableau LTYP2
  1031. c contiendront eux même les maillages raffinés issue des
  1032. c différents sous domaine de IPT5.
  1033. c LTYPE3(I,1)=LTYP2 pour le sous domaine I de IPT5
  1034. c LTYPE3(I,2)=IRR le permier indice d'un maillage 48 dans LTYP2
  1035. c-----------------------------------------------------------------------
  1036. NBLIS3=MAX(1,IPT5.LISOUS(/1))
  1037. SEGINI LTYPE3
  1038. MCONF=0
  1039. NSUR = 0
  1040. DO 33 INB=1,MAX(1,IPT5.LISOUS(/1))
  1041. NBNN=0
  1042. NBELEM=0
  1043. NBREF=0
  1044. NBSOUS=0
  1045. SEGINI IPT2
  1046. IF (IPT5.LISOUS(/1).EQ.0) THEN
  1047. IPT2=IPT5
  1048. ELSE
  1049. IPT2=IPT5.LISOUS(INB)
  1050. SEGACT IPT2
  1051. ENDIF
  1052.  
  1053. IMODE1.IMAMOD=IPT2
  1054. IMODE1.NEFMOD=IPT2.ITYPEL
  1055. C
  1056. C-------------------------------------------------------
  1057. C B) Creation d'un chamelem constant egal a 1 : MCHEL1 !
  1058. C-------------------------------------------------------
  1059. CALL ZEROP(MMODE1,'NOEUD ',MCHEL1)
  1060. SEGACT MCHEL1*MOD
  1061. MCHEL1.IMACHE(1)=IPT2
  1062. IF (IDIM.EQ.2) MCHEL1.IFOCHE=-1
  1063. IF (IDIM.EQ.3) MCHEL1.IFOCHE=2
  1064. MCHAM1=MCHEL1.ICHAML(1)
  1065. SEGACT MCHAM1
  1066. MELVA1=MCHAM1.IELVAL(1)
  1067. SEGACT MELVA1*MOD
  1068. c DO 399 J=1,IPT2.NUM(/2)
  1069. c MELVA1.VELCHE(1,J)=1.0
  1070. c 399 CONTINUE
  1071. MELVA1.VELCHE(1,1)=1.D0
  1072. C
  1073. C-------------------------------------------------------------
  1074. C C) Calcul de la surface ou du volume des elements : MCHEL2 !
  1075. C-------------------------------------------------------------
  1076. CALL INTGCA(MMODE1,MCHEL1,0,1,IRET,XRET,MCHEL2)
  1077. IF (XRET.LT.XPETI) THEN
  1078. WRITE (*,*) 'Volume(resp surf) de la zone a raffiner nul'
  1079. RETURN
  1080. ENDIF
  1081. SEGACT IPT2
  1082.  
  1083. C
  1084. C---------------------------------------------------------------------
  1085. C D) Creation d'un chamelem indiquant les elements a raffiner MCHAM2 !
  1086. C---------------------------------------------------------------------
  1087. SEGACT KMILIE*MOD
  1088. SEGACT MCHEL2
  1089. MCHAM2=MCHEL2.ICHAML(1)
  1090. SEGACT MCHAM2
  1091. MELVA2=MCHAM2.IELVAL(1)
  1092. SEGACT MELVA2*MOD
  1093. NACREE=0
  1094. KARAF=0
  1095.  
  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=(XCOOR(/1)/(IDIM+1))+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)
  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 probleme a la
  1793. c jonctions entre deux sous zones dans le cas ou on raffine plusieurs
  1794. c types d'elments différents (passer cette etape de concaténation 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.  
  1904. C
  1905. C=======================================================================
  1906. C Test sur la progression du raffinement
  1907. C=======================================================================
  1908. C
  1909. C
  1910. C CB215821 Retrait du caractere *MOD le cas echeant
  1911. SEGACT IPT5,MCHPO1,MSOUPO,MPOVA1,IPT3,KMILIE
  1912. c SEGSUP ICPR,KARETE,KARPOS,KARET2,ICHPO,IDENSI
  1913. c SEGSUP ICPR,KARETE,KARPOS,KARET2,IDENSI
  1914. SEGSUP KARETE,KARPOS,KARET2
  1915. SEGSUP MCHAM1,MELVA1,MCHAM2,MELVA2,IMODE1
  1916. SEGSUP MMODE1,MCHEL1,MCHEL2
  1917.  
  1918. IF ((IDIM.EQ.3).AND.(NBOUCL.NE.1)) SEGSUP JPLANS,JPLAN2,JPLAN3,
  1919. +JPLCOM,JNOEFA,IWORK2,JFARAF, IPT9
  1920.  
  1921.  
  1922. C S'il y a eu a raffiner des elements lors de l'iteration courante,
  1923. C on passe a l'iteration suivante
  1924. C
  1925. c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1926. c forcage : nb max itération de raff
  1927. c IF (NBOUCL.EQ.3) THEN
  1928. c KARAF2=0
  1929. c KARAF=0
  1930. c ENDIF
  1931. c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1932. IF (KARAF2.NE.0) GOTO 9
  1933.  
  1934.  
  1935.  
  1936. C
  1937. C ### SORTIE TEMPORAIRE DU PROGRAMME ###
  1938. C Si on est a la premiere iteration de RAFF et qu'il n'y a pas
  1939. C d'elements a raffiner, on retourne en sortie le maillage de depart,
  1940. C donné en entree
  1941. IF ((KARAF2.EQ.0).AND.(NBOUCL.EQ.1)) THEN
  1942. WRITE(*,*) 'il n est pas necessaire de raffiner le maillage'
  1943. MELEME=IPT1
  1944. GOTO 998
  1945.  
  1946. ENDIF
  1947.  
  1948.  
  1949. C
  1950. C
  1951. C=======================================================================
  1952. C Mises a jour diverses
  1953. C=======================================================================
  1954. C ** On affecte une densite a chaque noeud que l'on a cree
  1955. c DO 37 I=NPTINI+1,XCOOR(/1)/(IDIM+1)
  1956. c XCOOR(I*(IDIM+1))=DENSI(ICPR(I,1))
  1957. c 37 CONTINUE
  1958.  
  1959. C
  1960. C=======================================================================
  1961. C Fin du programme
  1962. C=======================================================================
  1963.  
  1964. MELEME=IPT7
  1965.  
  1966. SEGSUP XDEN
  1967. 998 CONTINUE
  1968.  
  1969. CALL ACTOBJ('MAILLAGE',MELEME,1)
  1970. CALL ECROBJ('MAILLAGE',MELEME)
  1971.  
  1972. END
  1973.  
  1974.  
  1975.  

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