Télécharger raff.eso

Retour à la liste

Numérotation des lignes :

raff
  1. C RAFF SOURCE OF166741 24/10/07 21:15:44 12016
  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.  
  1060. C-------------------------------------------------------
  1061. C B) Creation d'un chamelem constant egal a 1 : MCHEL1 !
  1062. C-------------------------------------------------------
  1063. CALL ZEROP(MMODE1,'NOEUD ',MCHEL1)
  1064. SEGACT MCHEL1*MOD
  1065. MCHEL1.IMACHE(1)=IPT2
  1066. IF (IDIM.EQ.2) MCHEL1.IFOCHE=-1
  1067. IF (IDIM.EQ.3) MCHEL1.IFOCHE=2
  1068. MCHAM1=MCHEL1.ICHAML(1)
  1069. SEGACT MCHAM1
  1070. MELVA1=MCHAM1.IELVAL(1)
  1071. SEGACT MELVA1*MOD
  1072. MELVA1.VELCHE(1,1)=1.D0
  1073. C
  1074. C-------------------------------------------------------------
  1075. C C) Calcul de la surface ou du volume des elements : MCHEL2 !
  1076. C-------------------------------------------------------------
  1077. CALL INTGCA(MMODE1,MCHEL1,0,1,IRET,XRET,MCHEL2)
  1078. IF (XRET.LT.XPETI) THEN
  1079. WRITE (*,*) 'Volume(resp surf) de la zone a raffiner nul'
  1080. RETURN
  1081. ENDIF
  1082. SEGACT IPT2
  1083.  
  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. SEGACT MELVA1
  1097. c boucles sur les éléments de ipt2
  1098. DO 32 J=1,IPT2.NUM(/2)
  1099.  
  1100. XMOY=0
  1101. DO 29 I=1,IPT2.NUM(/1)
  1102. NGLOB=IPT2.NUM(I,J)
  1103. NLOC=ICPR(NGLOB,1)
  1104. DENPT=DENSI(NLOC)
  1105. XMOY=XMOY+DENPT/(IPT2.NUM(/1))
  1106. 29 CONTINUE
  1107.  
  1108. XINT=MELVA2.VELCHE(1,J)
  1109. XDIM=IDIM*1.0
  1110. C Un element est a raffiner si la moyenne des densites affectees a
  1111. C l'ensemble de ses noeuds est inférieure a la racine carree de sa
  1112. C surface (en 2D) ou a la racine cubique de son volume (en 3D)
  1113. C WRITE (*,*) 'XMOY', XMOY, '(XINT**(1/XDIM))',(XINT**(1/XDIM))
  1114. IF ((XMOY-(XINT**(1/XDIM))).LT.0) THEN
  1115. c c y a t'il des +epsilon...?
  1116. c XDIFF = XMOY-(XINT**(1/XDIM))
  1117. c if(XDIFF.ge.0.and.XDIFF.lt.0.001*XMOY) then
  1118. c write(*,*) 'RAFF : XMOY,XDIFF=',XMOY,XDIFF,
  1119. c & 'pour l EF',J,(IPT2.NUM(iou,J),iou=1,IPT2.NUM(/1))
  1120. c endif
  1121. c IF (XDIFF.LT.0) THEN
  1122. MELVA2.VELCHE(1,J)=1
  1123. c KARAF nombre d'éléments à rafiner dans le sous domaine
  1124. KARAF=KARAF+1
  1125. ILPT=LPT(IPT2.ITYPEL)
  1126. ILPL=LPL(IPT2.ITYPEL)
  1127.  
  1128. c boucle sur les segments de l'élément à raffiner
  1129. DO 31 K=1,ILPL*2-1,2
  1130. NPTA=IPT2.NUM(KSEGM(ILPT+K-1),J)
  1131. NPTB=IPT2.NUM(KSEGM(ILPT+K),J)
  1132. NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
  1133. NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
  1134. NEXIST=0
  1135. DO 30 I=1,MAX(1,KARPOS(NMIN))
  1136. IF (KARETE(NMIN,I).EQ.NMAX) NEXIST=I
  1137. 30 CONTINUE
  1138. c On met dans kmilie les arêtes où il faudra mettre un nouveau noeud
  1139. c On leur associe la valeur -1
  1140. IF (KMILIE(NMIN,NEXIST).EQ.0) THEN
  1141. NACREE=NACREE+1
  1142. KMILIE(NMIN,NEXIST)=-1
  1143. ENDIF
  1144. 31 CONTINUE
  1145. C
  1146. C-------------------------------------------
  1147. C E) Comptage du nombre de faces a raffiner !
  1148. C-------------------------------------------
  1149. C 1/ Initialisations
  1150. IF (IDIM.NE.3) GOTO 32
  1151. JTYPEL=IPT2.ITYPEL
  1152. JLTEL1=LTEL(1,JTYPEL)
  1153. JLTEL2=LTEL(2,JTYPEL)
  1154. DO 73 IBB=1,JLTEL1
  1155. JLDEL1=LDEL(1,JLTEL2-1+IBB)
  1156. JTYFAC=IWORK1(JLDEL1)
  1157. JLDEL2=LDEL(2,JLTEL2-1+IBB)
  1158. C JNBPTS=NBNNE(JTYFAC)
  1159. JNBSOM=NBSOM(JTYFAC)
  1160. JSPOS=NSPOS(JTYFAC)
  1161. SEGINI IWORK2
  1162. DO 74 IAA=1,JNBSOM
  1163. NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),J)
  1164. IWORK2(IAA)=NGLOBA
  1165. 74 CONTINUE
  1166. C
  1167. C 2/ Classement des 3 sommets par ordre croissant (num globale)
  1168. NPTA=nbpts+1
  1169. NPTB=NPTA+1
  1170. NPTC=NPTB+1
  1171. DO 75 ICC=1,JNBSOM
  1172. IF (IWORK2(ICC).LT.NPTA) THEN
  1173. NPTC=NPTB
  1174. NPTB=NPTA
  1175. NPTA=IWORK2(ICC)
  1176. ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  1177. NPTC=NPTB
  1178. NPTB=IWORK2(ICC)
  1179. ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  1180. NPTC=IWORK2(ICC)
  1181. ENDIF
  1182. 75 CONTINUE
  1183. C
  1184. C 3/ Classement des 3 sommets par ordre croissant (num locale)
  1185. NPTA2=ICPR(NPTA,1)
  1186. NPTB2=ICPR(NPTB,1)
  1187. NPTC2=ICPR(NPTC,1)
  1188. IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  1189. NPTA=NPTA2
  1190. NPTB=MIN(NPTB2,NPTC2)
  1191. NPTC=MAX(NPTB2,NPTC2)
  1192. ENDIF
  1193. IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  1194. NPTA=NPTB2
  1195. NPTB=MIN(NPTA2,NPTC2)
  1196. NPTC=MAX(NPTA2,NPTC2)
  1197. ENDIF
  1198. IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  1199. NPTA=NPTC2
  1200. NPTB=MIN(NPTA2,NPTB2)
  1201. NPTC=MAX(NPTA2,NPTB2)
  1202. ENDIF
  1203. C
  1204. C 4/ Calcul du nombre de points a creer pour raffiner la face
  1205. NEXIS2=0
  1206. DO 76 IEE=1,JPLCOM(NPTA)
  1207. MTMP=JPLANS(NPTA,IEE)
  1208. JJ1=JNOEFA(MTMP,1)
  1209. JJ2=JNOEFA(MTMP,2)
  1210. JJ3=JNOEFA(MTMP,3)
  1211. IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
  1212. NEXIS2=IEE
  1213. ENDIF
  1214. 76 CONTINUE
  1215. IF (NEXIS2.NE.0) THEN
  1216. KFARAF=JPLAN2(NPTA,NEXIS2)
  1217. IF (KFARAF.EQ.0) THEN
  1218. NACREE=NACREE+NBINTE(JTYFAC)
  1219. JPLAN2(NPTA,NEXIS2)=JPLAN2(NPTA,NEXIS2)+1
  1220. ENDIF
  1221. ENDIF
  1222. 73 CONTINUE
  1223. ELSE
  1224.  
  1225. MELVA2.VELCHE(1,J)=0
  1226. ENDIF
  1227. 32 CONTINUE
  1228.  
  1229. c KARAF2 nombre d'éléments à rafiner sur tout les sous domaines
  1230. KARAF2=KARAF2+KARAF
  1231. c si aucun élément à raffiner on stoque directement IPT2 dans LTYPE3
  1232. c on stoque également les anciennes relations et on passe au
  1233. c sous-domaine suivant.
  1234.  
  1235. c ATTENTION SI NBOUCL=1 il est possible que des relations surabondantes
  1236.  
  1237. c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1238. c forcage : nb max itération de raff
  1239. c IF (NBOUCL.EQ.3) THEN
  1240. c KARAF2=0
  1241. c KARAF=0
  1242. c ENDIF
  1243. c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1244. IF (KARAF.EQ.0) THEN
  1245. IF (NBOUCL.EQ.1) THEN
  1246. IF (LCONF3.GT.0) THEN
  1247. SEGACT IPT6
  1248. NBLIS2=1+MAX(1,IPT6.LISOUS(/1))
  1249. SEGINI LTYPE2
  1250. LTYPE2(1,1)=1
  1251. LTYPE2(1,2)=IPT2.ITYPEL
  1252. LTYPE2(1,3)=IPT2
  1253. LTYPE2(1,4)=IPT2.ICOLOR(1)
  1254. LTYPE3(INB,2)=0
  1255. DO 676 MIE=1,MAX(1,IPT6.LISOUS(/1))
  1256. IF (IPT6.LISOUS(/1).EQ.0) THEN
  1257. IPT4=IPT6
  1258. ELSE
  1259. IPT4=IPT6.LISOUS(MIE)
  1260. ENDIF
  1261. NSUR=1
  1262. SEGACT IPT4
  1263. LTYPE2(1+MIE,1)=1+MIE
  1264. LTYPE2(1+MIE,2)=IPT4.ITYPEL
  1265. LTYPE2(1+MIE,3)=IPT4
  1266. LTYPE2(1+MIE,4)=IPT4.ICOLOR(1)
  1267. LTYPE3(INB,2)=2
  1268. 676 CONTINUE
  1269. ELSE
  1270. NBLIS2=1
  1271. SEGINI LTYPE2
  1272. LTYPE2(1,1)=1
  1273. LTYPE2(1,2)=IPT2.ITYPEL
  1274. LTYPE2(1,3)=IPT2
  1275. LTYPE2(1,4)=IPT2.ICOLOR(1)
  1276. LTYPE3(INB,2)=0
  1277. ENDIF
  1278. LTYPE3(INB,1)=LTYPE2
  1279. ELSE
  1280. SEGACT LTYPE4
  1281. LTYPE3(INB,1)=LTYPE4(INB,1)
  1282. LTYPE3(INB,2)=LTYPE4(INB,2)
  1283. ENDIF
  1284. GOTO 33
  1285. ENDIF
  1286. C
  1287. C----------------------------
  1288. C F) Raffinement du maillage !
  1289. C----------------------------
  1290. NML=2
  1291. C
  1292. C 1/ Raffinement d'un maillage 2D
  1293. C On recupere dans MELEME le maillage raffine sans ses relations de
  1294. C conformite et dans IPT8 les relations de ce maillage
  1295. C (IPT8 encore incomplet en sortie de RAFF2D)
  1296. c entrees de raff2D:
  1297. c - IPT2: maillage élémentaire à raffiner
  1298. c (ici sous domaine de ipt5)
  1299. c - ICPR: tableau de passage noeuds local/global
  1300. c (ici à partir de ipt5 interne à la boucle 9)
  1301. c - KARPOS: tableau du nb d'arretes par noeuds
  1302. c (ici à partir de ipt5, mis a jour avec ipt6)
  1303. c - KARETE: tableau des d'arretes KARETE(i,j)=k
  1304. c (ici à partir de ipt5, mis a jour avec ipt6)
  1305. c - KMILIE: tableau des hanging nodes associés au arrètes
  1306. c en entrée: à partir des relations de ipt6 et
  1307. c des nouveau noeuds à construires
  1308. c si relation issue de ipt6
  1309. c KMILIE(i,j)=n : le noeud global n est situe
  1310. c au milieu de l'arete [i,k]
  1311. c si nouvelle relation
  1312. c KMILIE(i,j)=-1 : il faut construire un noeud
  1313. c au milieu de l'arete [i,k]
  1314. c - MELVA2: champ par elem valait 1 pour les élem à raffiner
  1315. c et 0 pour les autres
  1316. c (ici à partir de ipt2 et de densi)
  1317. c - NACREE: nombre de noeuds à créer dans ipt2
  1318. c - KARAF: nombre d'élément à raffiner dans ipt2
  1319. c - IDCP : ??? variable non utilisé ???
  1320. c - XDEN : densité aux noeuds en notation globale
  1321. c (ici à partir de mchpo1)
  1322. c - KARET2: tableau du nombre d'éléments touchant les arretes
  1323. c (ici à partir de ipt5, mis a jour avec ipt6)
  1324. csorties de raff2D :
  1325. c - KMILIE: tableau des hanging nodes associés au arrètes
  1326. c en sortie: à partir des relations donnée en
  1327. c entrée et de celles crées
  1328. c si anciene relation à garder
  1329. c KMILIE(i,j)=n : le noeud global n est situe
  1330. c au milieu de l'arete [i,k]
  1331. c si anciene relation à suprimer
  1332. c KMILIE(i,j)=0 :
  1333. c si nouvelle relation
  1334. c KMILIE(i,j)=n : le noeud n à été construit au
  1335. c milieu de l'arete [i,k]
  1336. c - MELEME: Maillage élémentaire raffiné à partir de ipt2
  1337. c sans relations
  1338. c - IPT8 : Maillage de relation incomplet
  1339. c si J nouvel élément 48:
  1340. c IPT8.NUM(1,j) = nouveau hanging node
  1341. c IPT8.NUM(2:4,j) = autres noeuds de la relation
  1342. c IPT8.NUM(5,j) = 1 si arête a 2 noeuds
  1343. c 2 si arête a 3 noeuds
  1344. c si J nouvel élément 48:
  1345. c IPT8.NUM(1,j) = nouveau hanging node
  1346. c IPT8.NUM(2:5,j) = 0
  1347. c - XDEN : densité aux noeuds en notation globale
  1348. c interpolé aux nouveaux noeuds.
  1349. IF (IDIM.EQ.2) CALL RAFF2D(IPT2,ICPR,KARPOS,KARETE,KMILIE,
  1350. + MELVA2,NACREE,KARAF,MELEME,IDCP,IPT8,KARET2,XDEN)
  1351. C
  1352. C 2/ Raffinement d'un maillage 3D
  1353. IF (IDIM.EQ.3) THEN
  1354. NML=0
  1355. IF ((IPT2.ITYPEL.EQ.25).OR.(IPT2.ITYPEL.EQ.26)) NML=1
  1356. ENDIF
  1357. C a) Maillage 3D contenant des pyramides
  1358. IPT9=0
  1359. IF (NML.EQ.1) CALL RAFPYR(IPT2,ICPR,KARPOS,KARETE,KMILIE,
  1360. + MELVA2,NACREE,KARAF,MELEME,JPLANS,JPLAN3,JPLCOM,JNOEFA,
  1361. + IPT8,JFARAF,KARET2,IPT9,XDEN)
  1362. C b) Maillage 3D sans pyramides
  1363. IF (NML.EQ.0) CALL RAFF3D(IPT2,IPT6,ICPR,KARPOS,KARETE,KMILIE,
  1364. + MELVA2,NACREE,KARAF,MELEME,JPLANS,JPLAN3,JPLCOM,JNOEFA,
  1365. + IPT8,JFARAF,KARET2,XDEN)
  1366. C
  1367. C=======================================================================
  1368. C Creation du maillage de relations de conformite
  1369. C=======================================================================
  1370. C--------------------
  1371. C A) Initialisations !
  1372. C--------------------
  1373. IPT10=0
  1374. IPT11=0
  1375. IPT12=0
  1376. IPT13=0
  1377. IPT14=0
  1378. IPT15=0
  1379. IPT16=0
  1380. MCONF=0
  1381. c IF (KARAF2.EQ.0) GOTO 33
  1382. IF (IPT8.NE.0) THEN
  1383. C
  1384. C---------------------------------------
  1385. C B) Traitement des anciennes relations !
  1386. C---------------------------------------
  1387. C A ce stade, IPT8 contient les relations crees lors du dernier
  1388. C raffinement (noeud de relation + noeuds formant la relation) et les
  1389. C anciennes relations a reporter (noeud de relation seulement)
  1390.  
  1391. C 1/ Recherche de l'existence d'anciennes relations a reporter
  1392. NB10=0
  1393. MLONG=IPT8.NUM(/1)
  1394. DO 775 MM=1,IPT8.NUM(/2)
  1395. IF (IPT8.NUM(MLONG,MM).EQ.0) NB10=NB10+1
  1396. 775 CONTINUE
  1397.  
  1398. IF (NB10.NE.0) THEN
  1399. C
  1400. C 2/ Creation de IPT10
  1401. c MCOL=0
  1402. c NBNN=9
  1403. c NBELEM=NB10
  1404. c NBREF=0
  1405. c NBSOUS=0
  1406. c SEGINI IPT10
  1407. c IPT10.ITYPEL=48
  1408. C
  1409. C 3/ Recherche des anciennes relations a reporter
  1410. C Cette partie permet de compléter IPT8 avec les noeuds formant les
  1411. C relations, pour les anciennes relations à reporter
  1412. DO 628 MM11=1,IPT8.NUM(/2)
  1413. IF (IPT8.NUM(MLONG,MM11).NE.0) GOTO 628
  1414. INRELA=IPT8.NUM(1,MM11)
  1415. c On active pous le cas ou IPT6=IPT2.
  1416. SEGACT IPT6
  1417. DO 625 KOU=1,MAX(1,IPT6.LISOUS(/1))
  1418. IF (IPT6.LISOUS(/1).EQ.0) THEN
  1419. IPT2=IPT6
  1420. ELSE
  1421. IPT2=IPT6.LISOUS(KOU)
  1422. ENDIF
  1423. SEGACT IPT2
  1424. IF (IPT2.ITYPEL.NE.48) GOTO 625
  1425. DO 626 JKG=1,IPT2.NUM(/2)
  1426. INOEGL=IPT2.NUM(1,JKG)
  1427. IF (INOEGL.NE.INRELA) GOTO 626
  1428. MCOL=MCOL+1
  1429. DO 627 MTMP=1,IPT2.NUM(/1)
  1430. IPT8.NUM(MTMP,MM11)=IPT2.NUM(MTMP,JKG)
  1431. IPT8.NUM(MLONG,MM11)= IPT2.ICOLOR(JKG)
  1432. 627 CONTINUE
  1433. GOTO 628
  1434. 626 CONTINUE
  1435. 625 CONTINUE
  1436. 628 CONTINUE
  1437. ENDIF
  1438. c à ce stade le maillage de relation ipt8 est complet
  1439. C
  1440. C---------------------------------------
  1441. C C) Traitement des nouvelles relations !
  1442. C---------------------------------------
  1443. C Ici, IPT8 contient toutes les relations du nouveau maillage
  1444. C On les trie par type de relations
  1445.  
  1446. C 1/ Recherche de l'existence de nouvelles relations
  1447. NB11=0
  1448. NB12=0
  1449. NB13=0
  1450. NB14=0
  1451. NB15=0
  1452. NB16=0
  1453. MLONG=IPT8.NUM(/1)
  1454. DO 776 MM=1,IPT8.NUM(/2)
  1455. IF (IPT8.NUM(MLONG,MM).EQ.1) NB11=NB11+1
  1456. IF (IPT8.NUM(MLONG,MM).EQ.2) NB12=NB12+1
  1457. IF (IPT8.NUM(MLONG,MM).EQ.3) NB13=NB13+1
  1458. IF (IPT8.NUM(MLONG,MM).EQ.4) NB14=NB14+1
  1459. IF (IPT8.NUM(MLONG,MM).EQ.5) NB15=NB15+1
  1460. IF (IPT8.NUM(MLONG,MM).EQ.6) NB16=NB16+1
  1461. 776 CONTINUE
  1462. C
  1463. C 2/ Cas des relations 'aretes a 2 noeuds'
  1464. IF (NB11.NE.0) THEN
  1465. MCONF=MCONF+1
  1466. MCOL=0
  1467. NBNN=3
  1468. NBELEM=NB11
  1469. NBREF=0
  1470. NBSOUS=0
  1471. SEGINI IPT11
  1472. IPT11.ITYPEL=48
  1473. DO 778 MM11=1,IPT8.NUM(/2)
  1474. IF (IPT8.NUM(MLONG,MM11).NE.1) GOTO 778
  1475. MCOL=MCOL+1
  1476. DO 777 MTMP=1,3
  1477. IPT11.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
  1478. IPT11.ICOLOR(MCOL)=1
  1479. 777 CONTINUE
  1480. 778 CONTINUE
  1481. ENDIF
  1482. C
  1483. C 3/ Cas des relations 'aretes a 3 noeuds'
  1484. IF (NB12.NE.0) THEN
  1485. MCONF=MCONF+1
  1486. MCOL=0
  1487. NBNN=4
  1488. NBELEM=NB12
  1489. NBREF=0
  1490. NBSOUS=0
  1491. SEGINI IPT12
  1492. IPT12.ITYPEL=48
  1493. DO 668 MM11=1,IPT8.NUM(/2)
  1494. IF (IPT8.NUM(MLONG,MM11).NE.2) GOTO 668
  1495. MCOL=MCOL+1
  1496. DO 667 MTMP=1,4
  1497. IPT12.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
  1498. IPT12.ICOLOR(MCOL)=2
  1499. 667 CONTINUE
  1500. 668 CONTINUE
  1501. ENDIF
  1502. C
  1503. C 4/ Cas des relations 'face QUA4'
  1504. IF (NB13.NE.0) THEN
  1505. MCONF=MCONF+1
  1506. MCOL=0
  1507. NBNN=5
  1508. NBELEM=NB13
  1509. NBREF=0
  1510. NBSOUS=0
  1511. SEGINI IPT13
  1512. IPT13.ITYPEL=48
  1513. DO 658 MM11=1,IPT8.NUM(/2)
  1514. IF (IPT8.NUM(MLONG,MM11).NE.3) GOTO 658
  1515. MCOL=MCOL+1
  1516. DO 657 MTMP=1,5
  1517. IPT13.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
  1518. IPT13.ICOLOR(MCOL)=3
  1519. 657 CONTINUE
  1520. 658 CONTINUE
  1521. ENDIF
  1522. C
  1523. C 5/ Cas des relations 'face TRI6'
  1524. IF (NB14.NE.0) THEN
  1525. MCONF=MCONF+1
  1526. MCOL=0
  1527. NBNN=7
  1528. NBELEM=NB14
  1529. NBREF=0
  1530. NBSOUS=0
  1531. SEGINI IPT14
  1532. IPT14.ITYPEL=48
  1533. DO 648 MM11=1,IPT8.NUM(/2)
  1534. IF (IPT8.NUM(MLONG,MM11).NE.4) GOTO 648
  1535. MCOL=MCOL+1
  1536. DO 647 MTMP=1,7
  1537. IPT14.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
  1538. IPT14.ICOLOR(MCOL)=4
  1539. 647 CONTINUE
  1540. 648 CONTINUE
  1541. ENDIF
  1542. C
  1543. C 6/ Cas des relations 'face QUA8-1'
  1544. IF (NB15.NE.0) THEN
  1545. MCONF=MCONF+1
  1546. MCOL=0
  1547. NBNN=9
  1548. NBELEM=NB15
  1549. NBREF=0
  1550. NBSOUS=0
  1551. SEGINI IPT15
  1552. IPT15.ITYPEL=48
  1553. DO 638 MM11=1,IPT8.NUM(/2)
  1554. IF (IPT8.NUM(MLONG,MM11).NE.5) GOTO 638
  1555. MCOL=MCOL+1
  1556. DO 637 MTMP=1,9
  1557. IPT15.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
  1558. IPT15.ICOLOR(MCOL)=5
  1559. 637 CONTINUE
  1560. 638 CONTINUE
  1561. ENDIF
  1562. C
  1563. C 7/ Cas des relations 'face QUA8-2'
  1564. IF (NB16.NE.0) THEN
  1565. MCONF=MCONF+1
  1566. MCOL=0
  1567. NBNN=9
  1568. NBELEM=NB16
  1569. NBREF=0
  1570. NBSOUS=0
  1571. SEGINI IPT16
  1572. IPT16.ITYPEL=48
  1573. DO 678 MM11=1,IPT8.NUM(/2)
  1574. IF (IPT8.NUM(MLONG,MM11).NE.6) GOTO 678
  1575. MCOL=MCOL+1
  1576. DO 677 MTMP=1,9
  1577. IPT16.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
  1578. IPT16.ICOLOR(MCOL)=6
  1579. 677 CONTINUE
  1580. 678 CONTINUE
  1581. ENDIF
  1582. C
  1583. ENDIF
  1584. C
  1585. C--------------------------------------
  1586. C D) Creation du maillage de relations !
  1587. C--------------------------------------
  1588. C 1/ Initialisation de MAIREL
  1589. C MAIREL contient tous les maillages de relations, certains étant vides
  1590. C car il n'y a aucune relation du type considéré dans le nveau maillage
  1591. NBMAIL=6
  1592. SEGINI MAIREL
  1593. MAIREL(1)=IPT11
  1594. MAIREL(2)=IPT12
  1595. MAIREL(3)=IPT13
  1596. MAIREL(4)=IPT14
  1597. MAIREL(5)=IPT15
  1598. MAIREL(6)=IPT16
  1599. C
  1600. C 2/ Initialisation de LTYPE2
  1601. C LTYPE2 va contenir tous les sous-objets du nouveau maillage,
  1602. C relations et non relations
  1603. NBLIS2=MCONF+MAX(1,MELEME.LISOUS(/1))
  1604. SEGINI LTYPE2
  1605. C
  1606. C 3/ Renseignement des lisous non 48 dans LTYPE2
  1607. C S'il y a N sous-objets de non-relations dans le nouveau maillage,
  1608. C les N-ièmes premières colonnes de LTYPE2 contiennent ces sous-objets
  1609. DO 685 IZZ=1,MAX(1,MELEME.LISOUS(/1))
  1610. IF (MELEME.LISOUS(/1).EQ.0) THEN
  1611. IPT2=MELEME
  1612. ELSE
  1613. IPT2=MELEME.LISOUS(IZZ)
  1614. SEGACT IPT2
  1615. ENDIF
  1616. LTYPE2(IZZ,1)=IZZ
  1617. LTYPE2(IZZ,2)=IPT2.ITYPEL
  1618. LTYPE2(IZZ,3)=IPT2
  1619. LTYPE2(IZZ,4)=0
  1620. 685 CONTINUE
  1621. C
  1622. C 4/ Renseignement des lisous 48 dans LTYPE2
  1623. C Les colonnes suivantes contiennent les sous-objets de relations,
  1624. C chaque sous-objet correspondant a un type de relations de conformité
  1625. IZZ=NBLIS2-MCONF
  1626. DO 686 IXX=1,6
  1627. IF (MAIREL(IXX).EQ.0) GOTO 686
  1628. IPT2=MAIREL(IXX)
  1629. IZZ=IZZ+1
  1630. LTYPE2(IZZ,1)=IZZ
  1631. LTYPE2(IZZ,2)=IPT2.ITYPEL
  1632. LTYPE2(IZZ,3)=IPT2
  1633. LTYPE2(IZZ,4)=IPT2.ICOLOR(1)
  1634. 686 CONTINUE
  1635. C--------------------------------------
  1636. C E) Renseignement de LTYPE2 interne a
  1637. c la boucle 33 dans LTYPE3 externe a
  1638. c la boucle 33
  1639. C--------------------------------------C
  1640. C
  1641.  
  1642. LTYPE3(INB,1)=LTYPE2
  1643. IF (MCONF.EQ.0) THEN
  1644. LTYPE3(INB,2)=0
  1645. ELSE
  1646. LTYPE3(INB,2)=NBLIS2-MCONF+1
  1647. ENDIF
  1648. c fin de boucle sur les sous domaines d'ipt5
  1649. 33 CONTINUE
  1650. c SEGSUP IPT2
  1651.  
  1652.  
  1653. C Mise a jour des maillages en entrée sortie de la boucle 9
  1654. C IPT5, IPT6, IPT7
  1655. C=======================================================================
  1656. 679 CONTINUE
  1657. C--------------------------------------C
  1658. C A) initialisation
  1659. C--------------------------------------C
  1660.  
  1661. c 1/ initialisaton de IPT7 maillage stockant tout les sous maillage
  1662. c à la fin d'une itération de la boucle 9
  1663. NBS7=0
  1664. DO 331 I3=1,LTYPE3(/1)
  1665. LTYPE2=LTYPE3(I3, 1)
  1666. SEGACT LTYPE2
  1667. NBS7=NBS7+LTYPE2(/1)+1
  1668. 331 CONTINUE
  1669. IF (NBS7.EQ.1) NBS7=0
  1670. NBNN=0
  1671. NBELEM=0
  1672. NBREF=0
  1673. NBSOUS=NBS7
  1674. SEGINI IPT7
  1675.  
  1676. c 2/ initialisaton de IPT6 maillage stockant tout les sous maillage
  1677. c de relation à la fin d'une itération de la boucle 9
  1678. NBS6=0
  1679.  
  1680. DO 332 I3=1,LTYPE3(/1)
  1681. LTYPE2=LTYPE3(I3, 1)
  1682. SEGACT LTYPE2
  1683. IF (LTYPE3(I3,2).NE.0) NBS6=NBS6+(LTYPE2(/1)-LTYPE3(I3,2)+1)
  1684. 332 CONTINUE
  1685. IF (NBS6.EQ.1) NBS6=0
  1686. NBNN=0
  1687. NBELEM=0
  1688. NBREF=0
  1689. NBSOUS=NBS6
  1690. SEGINI IPT6
  1691. c 3/ initialisaton de IPT5 maillage stockant tout les sous maillage
  1692. c non relation à la fin d'une itération de la boucle 9
  1693. NBS5=0
  1694. DO 333 I3=1,LTYPE3(/1)
  1695. LTYPE2=LTYPE3(I3, 1)
  1696. SEGACT LTYPE2
  1697. IF (LTYPE3(I3,2).EQ.0) THEN
  1698. NBS5=NBS5+LTYPE2(/1)
  1699. ELSE
  1700. NBS5=NBS5+(LTYPE3(I3,2)-1)
  1701. ENDIF
  1702. 333 CONTINUE
  1703. IF (NBS5.EQ.1) NBS5=0
  1704. NBNN=0
  1705. NBELEM=0
  1706. NBREF=0
  1707. NBSOUS=NBS5
  1708. SEGINI IPT5
  1709. C--------------------------------------C
  1710. C B) assemblage des nouveaux maillage à
  1711. c partir des sous maillages stockés
  1712. c dans LTYP3
  1713. C--------------------------------------C
  1714. IND5=0
  1715. IND6=0
  1716. IND7=0
  1717. c 1/assemblage des nouveaux maillage non relation à partir des sous
  1718. c maillages stockés dans LTYP3
  1719.  
  1720. DO 35 I3=1,LTYPE3(/1)
  1721. c extraction d'un LTYPE2 de LTYPE3
  1722. LTYPE2=LTYPE3(I3, 1)
  1723. SEGACT LTYPE2
  1724. c NBNN=0
  1725. c NBELEM=0
  1726. c NBREF=0
  1727. c NBSOUS=0
  1728. DO 613 I2=1,MAX(1,(LTYPE3(I3,2)-1))
  1729. c extraction d'un MAILLAGE non relation de LTYPE2
  1730. IPT2=LTYPE2(I2,3)
  1731. SEGACT IPT2
  1732. IND5=IND5+1
  1733. IF (NBS5.EQ.0) THEN
  1734. IPT5=IPT2
  1735. c write (*,*) 'IPT5.' , IPT5
  1736. c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
  1737. ELSE
  1738. IPT5.LISOUS(IND5)=IPT2
  1739. c write (*,*) 'IPT5.LISOUS(',IND5,')' , IPT5.LISOUS(IND5)
  1740. c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
  1741. ENDIF
  1742. IND7=IND7+1
  1743. IF (NBS7.EQ.0) THEN
  1744. IPT7=IPT2
  1745. ELSE
  1746. IPT7.LISOUS(IND7)=IPT2
  1747. c write (*,*) 'IPT7.LISOUS(',IND7,')' , IPT7.LISOUS(IND7)
  1748. c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
  1749. ENDIF
  1750. 613 CONTINUE
  1751. 35 CONTINUE
  1752.  
  1753. c 2/assemblage des nouveaux maillage de relation à partir des derniers
  1754. c sous maillages stockés dans LTYPE3
  1755.  
  1756. DO 36 I3=1,LTYPE3(/1)
  1757. c extraction d'un LTYPE2 de LTYPE3
  1758. c write (*,*) 'I3', I3
  1759. c write (*,*) 'LTYPE3(I3,2)', LTYPE3(I3,2)
  1760. IF (LTYPE3(I3,2).EQ.0) GOTO 36
  1761. SEGACT LTYPE2
  1762. LTYPE2=LTYPE3(I3, 1)
  1763. NBNN=0
  1764. NBELEM=0
  1765. NBREF=0
  1766. NBSOUS=0
  1767. DO 614 I2=LTYPE3(I3,2), LTYPE2(/1)
  1768. c extraction d'un MAILLAGE de relation de LTYPE2
  1769. c write (*,*) 'I2', I2
  1770. SEGINI IPT2
  1771. IPT2=LTYPE2(I2,3)
  1772. SEGACT IPT2
  1773. c write (*,*) 'IPT2: Type', IPT2.ITYPEL,'nbelem', IPT2.NUM(/2)
  1774. c write (*,*) 'coul IPT2' , IPT2.ICOLOR(1)
  1775. IND6=IND6+1
  1776. IF (NBS6.EQ.0) THEN
  1777. IPT6=IPT2
  1778. c write (*,*) 'IPT6' , IPT6
  1779. c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
  1780. ELSE
  1781. IPT6.LISOUS(IND6)=IPT2
  1782.  
  1783. c write (*,*) 'IPT6.LISOUS(',IND6,')' , IPT6.LISOUS(IND6)
  1784. c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
  1785. ENDIF
  1786. 614 CONTINUE
  1787. 36 CONTINUE
  1788.  
  1789. c 3/ Concatenation des sous maillages de relation de même couleur
  1790.  
  1791. c GG 11 juin 2018 : attention il peut encore y avoir des problemes a la
  1792. c jonctions entre deux sous zones dans le cas ou on raffine plusieurs
  1793. c types d'element differents (passer cette etape de concatenation dans
  1794. c la boucle 33 ?
  1795.  
  1796. IF (NBS6.EQ.0) then
  1797. IND7=IND7+1
  1798. IPT7.LISOUS(IND7) = IPT6
  1799. c write (*,*) 'IPT7.LISOUS(',IND7,')' , IPT7.LISOUS(IND7)
  1800. c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
  1801.  
  1802. goto 370
  1803. ENDIF
  1804. NBS60 = 0
  1805. SEGACT IPT6
  1806.  
  1807. c boucle sur les couleurs
  1808. NBELEM = 0
  1809. NBNN= 0
  1810. NBREF= 0
  1811. NBSOUS=6
  1812. SEGINI IPT4
  1813. DO 37 ICOUL = 1, 6
  1814. c boucle sur les sous maillages de ipt6
  1815. NCOl = 0
  1816. DO 38 I6 = 1, NBS6
  1817. IPT3 = IPT6.LISOUS(I6)
  1818. SEGACT IPT3
  1819. IF ((IPT3.icolor(1)).EQ.ICOUL) then
  1820. NCOL = NCOL+1
  1821. IF (NCOL.EQ.1) THEN
  1822. SEGINI IPT2
  1823. IPT2 = IPT6.LISOUS(I6)
  1824.  
  1825. ENDIF
  1826. IF (NCOL.GE.2) THEN
  1827.  
  1828. c IPT3 = IPT6.LISOUS(I6)
  1829. c SEGACT IPT3
  1830. c Plusieurs sous maillages de même couleur on concatene IPT2 et IPT3
  1831. NBELEM = IPT2.NUM(/2)
  1832. NBNN= IPT2.NUM(/1)
  1833. NBREF=0
  1834. NBSOUS=0
  1835. NBELE0 = NBELEM
  1836. DO 381 IE3 = 1, IPT3.NUM(/2)
  1837. DO 382 IE2 = 1, NBELE0
  1838. IF ((IPT3.NUM(1,IE3)).EQ.(IPT2.NUM(1,IE2))) GOTO 381
  1839. 382 CONTINUE
  1840. NBELEM = NBELEM +1
  1841. SEGADJ IPT2
  1842. DO IN2 = 1, NBNN
  1843. IPT2.NUM(IN2 , NBELEM) = IPT3.NUM(IN2, IE3)
  1844. IPT2.icolor(NBELEM) = IPT3.icolor(IE3)
  1845. ENDDO
  1846. 381 CONTINUE
  1847. ENDIF
  1848. ENDIF
  1849. 38 CONTINUE
  1850. IF (NCOL.GE.1) THEN
  1851. IPT4.LISOUS(ICOUL) = IPT2
  1852. NBS60= NBS60+1
  1853. ENDIF
  1854.  
  1855. 37 CONTINUE
  1856.  
  1857. c mise a jour de ipt6 et ipt7
  1858. IF (NBS60.EQ.1) then
  1859. IPT6 = IPT2
  1860. segact ipt6
  1861.  
  1862. c write (*,*) 'IPT6' , IPT6
  1863. c write (*,*) 'type', IPT6.ITYPEL, 'Nbel' , IPT6.num(/2)
  1864. IND7=IND7+1
  1865. IPT7.LISOUS(IND7) = IPT2
  1866. c write (*,*) 'IPT7.LISOUS(',IND7,')' , IPT7.LISOUS(IND7)
  1867. c write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
  1868. ELSE
  1869. SEGACT IPT4
  1870. NBELEM = 0
  1871. NBNN = 0
  1872. NBREF = 0
  1873. NBSOUS = NBS60
  1874. SEGINI IPT6
  1875. IS6=0
  1876. DO ICOUL = 1, 6
  1877. IF ((IPT4.LISOUS(ICOUL)).NE.0) THEN
  1878. IS6 = IS6+1
  1879. IPT6.LISOUS(IS6)= IPT4.LISOUS(ICOUL)
  1880. c write (*,*) 'IPT6.LISOUS(',IND6,')' , IPT6.LISOUS(IND6)
  1881. c write (*,*) 'type', IPT4.ITYPEL, 'Nbel' , IPT4.num(/2)
  1882.  
  1883. IND7=IND7+1
  1884. IPT7.LISOUS(IND7) = IPT4.LISOUS(ICOUL)
  1885. ENDIF
  1886. ENDDO
  1887. ENDIF
  1888. 370 CONTINUE
  1889. NBELEM = 0
  1890. NBNN = 0
  1891. NBREF = 0
  1892. NBSOUS = IND7
  1893.  
  1894. SEGADJ IPT7
  1895.  
  1896. SEGINI LTYPE4
  1897. LTYPE4 = LTYPE3
  1898.  
  1899. c mise a jour du nombre global de relations
  1900. LCONF3 = MCONF
  1901.  
  1902. C=======================================================================
  1903. C Test sur la progression du raffinement
  1904. C=======================================================================
  1905. C
  1906. C CB215821 Retrait du caractere *MOD le cas echeant
  1907. SEGACT IPT5,MCHPO1,MSOUPO,MPOVA1,IPT3,KMILIE
  1908. c SEGSUP ICPR,KARETE,KARPOS,KARET2,ICHPO,IDENSI
  1909. c SEGSUP ICPR,KARETE,KARPOS,KARET2,IDENSI
  1910. SEGSUP KARETE,KARPOS,KARET2
  1911. SEGSUP MCHAM1,MELVA1,MCHAM2,MELVA2,IMODE1
  1912. SEGSUP MMODE1,MCHEL1,MCHEL2
  1913.  
  1914. IF ((IDIM.EQ.3).AND.(NBOUCL.NE.1)) SEGSUP JPLANS,JPLAN2,JPLAN3,
  1915. +JPLCOM,JNOEFA,IWORK2,JFARAF, IPT9
  1916.  
  1917. C S'il y a eu a raffiner des elements lors de l'iteration courante,
  1918. C on passe a l'iteration suivante
  1919. C
  1920. c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1921. c forcage : nb max itération de raff
  1922. c IF (NBOUCL.EQ.3) THEN
  1923. c KARAF2=0
  1924. c KARAF=0
  1925. c ENDIF
  1926. c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1927. IF (KARAF2.NE.0) GOTO 9
  1928.  
  1929. C ### SORTIE TEMPORAIRE DU PROGRAMME ###
  1930. C Si on est a la premiere iteration de RAFF et qu'il n'y a pas
  1931. C d'elements a raffiner, on retourne en sortie le maillage de depart,
  1932. C donné en entree
  1933. IF ((KARAF2.EQ.0).AND.(NBOUCL.EQ.1)) THEN
  1934. WRITE(*,*) 'il n est pas necessaire de raffiner le maillage'
  1935. MELEME=IPT1
  1936. GOTO 998
  1937.  
  1938. ENDIF
  1939.  
  1940. C=======================================================================
  1941. C Mises a jour diverses
  1942. C=======================================================================
  1943. C ** On affecte une densite a chaque noeud que l'on a cree
  1944. c DO 37 I=NPTINI+1,XCOOR(/1)/(IDIM+1)
  1945. c XCOOR(I*(IDIM+1))=DENSI(ICPR(I,1))
  1946. c 37 CONTINUE
  1947.  
  1948. C=======================================================================
  1949. C Fin du programme
  1950. C=======================================================================
  1951.  
  1952. MELEME=IPT7
  1953.  
  1954. SEGSUP XDEN
  1955. 998 CONTINUE
  1956.  
  1957. CALL ACTOBJ('MAILLAGE',MELEME,1)
  1958. CALL ECROBJ('MAILLAGE',MELEME)
  1959.  
  1960. END
  1961.  
  1962.  
  1963.  

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