Télécharger raff.eso

Retour à la liste

Numérotation des lignes :

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

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