Télécharger raff.eso

Retour à la liste

Numérotation des lignes :

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

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