Télécharger raff.eso

Retour à la liste

Numérotation des lignes :

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

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