Télécharger raff3d.eso

Retour à la liste

Numérotation des lignes :

  1. C RAFF3D SOURCE GG250959 17/03/15 21:15:10 9344
  2. SUBROUTINE RAFF3D(IPT2,IPT3,ICPR,KARPOS,KARETE,KMILIE,MELVA2,
  3. .NACREE,KARAF,IPT4,JPLANS,JPLAN3,JPLCOM,JNOEFA,IPT7,JFARAF,KARET2,
  4. .XDEN)
  5. c entrees de raff3D:
  6. c - IPT2: maillage élémentaire à raffiner
  7. c - IPT3: maillage de SURE précédant
  8. c - ICPR: tableau de passage noeuds local/global
  9. c - KARPOS: tableau du nb d'arretes par noeuds
  10. c - KARETE: tableau des d'arretes KARETE(i,j)=k
  11. c - KMILIE: tableau des hanging nodes associés au arêtes
  12. c si ancienne relation
  13. c KMILIE(i,j)=n : le noeud global n est situe
  14. c au milieu de l'arete [i,k]
  15. c si nouvelle relation
  16. c KMILIE(i,j)=-1 : il faut construire un noeud
  17. c au milieu de l'arete [i,k]
  18. c - MELVA2: champ par elem valait 1 pour les élem à raffiner
  19. c et 0 pour les autres
  20. c - NACREE: nombre de noeuds à créer dans ipt2
  21. c - KARAF: nombre d'élément à raffiner dans ipt2
  22. c - JPLANS : tableau des faces JPLAN(i,j)=k
  23. c - JPLAN3 : tableau du nombre de relation de compat par face
  24. c - JPLCOM : tableau du nb de faces par noeuds
  25. c - JNOEFA : tableau des noeuds et éléments associé au faces
  26. c - JFARAF : tableau des hanging nodes associés au faces pour
  27. c les anciennes relations.
  28. c - XDEN : densité aux noeuds en notation globale
  29. c - KARET2: tableau du nombre d'éléments touchant les arêtes
  30. c sorties:
  31. c - KMILIE: tableau des hanging nodes associés au arêtes
  32. c (en sortie: à partir des relations donnée en
  33. c entrée et de celles crées
  34. c si anciene relation à garder
  35. c KMILIE(i,j)=n : le noeud global n est situe
  36. c au milieu de l'arete [i,k]
  37. c si anciene relation à supprimer
  38. c KMILIE(i,j)=0 :
  39. c si nouvelle relation
  40. c KMILIE(i,j)=n : le noeud n à été construit au
  41. c milieu de l'arete [i,k]
  42. c - JPLAN3 : tableau du nombre de relation de compat par face
  43. c mis a jour
  44. c - JFARAF : tableau des hanging nodes associés au faces pour
  45. c mis a jour.
  46. c - MELEME: Maillage élémentaire raffiné à partir de ipt2
  47. c sans relations
  48. c - IPT7 : Maillage de relations créé incomplet
  49. c - XDEN : densité aux noeuds en notation globale
  50. c interpolé aux nouveaux noeuds.
  51. IMPLICIT REAL*8(A-H,O-Z)
  52. IMPLICIT INTEGER(I-N)
  53. -INC CCOPTIO
  54. -INC SMCOORD
  55. -INC CCGEOME
  56. -INC SMELEME
  57. -INC SMCHPOI
  58. -INC SMMODEL
  59. -INC SMCHAML
  60. C
  61. C======================================================================
  62. C Declarations
  63. C======================================================================
  64. SEGMENT ICPR((XCOOR(/1)/(IDIM+1)),2)
  65. SEGMENT XDEN((XCOOR(/1)/(IDIM+1)))
  66. SEGMENT KARETE(NBNDS,NCOL)
  67. SEGMENT KARET2(NBNDS,NCOL)
  68. SEGMENT KMILIE(NBNDS,NCOL)
  69. SEGMENT KARPOS(NBNDS)
  70. SEGMENT JPLANS(JPLA1,JPLA2)
  71. SEGMENT JPLAN3(JPLA1,JPLA2)
  72. SEGMENT JPLCOM(JPLA1)
  73. SEGMENT JNOEFA(JNBFA,5)
  74. SEGMENT JFARAF(JNBFA,LLLL)
  75. SEGMENT NUMNOE(INUMNO)
  76. SEGMENT IWORK2(JNBSOM)
  77. SEGMENT IWORK1(IJKLMN)
  78. SEGMENT IWORK4(JNBSOM)
  79. SEGMENT IWORK3(IJKLMN)
  80.  
  81. REAL*8 LON68,LON510,LON47,LON29
  82. C
  83. C======================================================================
  84. C Initialisations
  85. C======================================================================
  86. SEGACT JPLANS,JPLCOM,JNOEFA*MOD,JPLAN3*MOD, IPT3
  87. SEGACT IPT2,ICPR,KARPOS,KARETE,KMILIE*MOD,MELVA2,JFARAF*MOD
  88. IJKLMN=4
  89. ipt7=0
  90. SEGINI IWORK1
  91. IWORK1(1)=4
  92. IWORK1(2)=8
  93. IWORK1(3)=6
  94. IWORK1(4)=10
  95. C JNBFA=JNOEFA(/1)
  96. C LLLL=2
  97. C IF (IPT2.ITYPEL.EQ.15) LLLL=10
  98. C IF ((IPT2.ITYPEL.EQ.17).OR.(IPT2.ITYPEL.EQ.24)) LLLL=6
  99. C LLLL=10
  100. C SEGINI JFARAF
  101. NCOMPL=LNELM(1,(IPT2.ITYPEL-1)*2+2)
  102. IF (NCOMPL.NE.0) GOTO 999
  103. NBNN=NBNNE(LNELM(2,(IPT2.ITYPEL-1)*2+1))
  104. INELM=LNELM(1,(IPT2.ITYPEL-1)*2+1)
  105. C
  106. C Creation du squelette du maillage resultat
  107. NBELEM=IPT2.NUM(/2)-KARAF+INELM*KARAF
  108. NBSOUS=0
  109. NBREF=0
  110. SEGINI IPT4
  111. IPT4.ITYPEL=LNELM(2,(IPT2.ITYPEL-1)*2+1)
  112. LPO2=LPOS2(IPT2.ITYPEL)
  113. NBPT0=XCOOR(/1)/(IDIM+1)
  114. NBPT1=XCOOR(/1)/(IDIM+1)
  115. NBPTS=NBPT1+NACREE+KARAF*NBINTE(IPT2.ITYPEL)
  116. SEGADJ MCOORD,XDEN
  117. INUMNO=NBRAF(IPT2.ITYPEL)
  118. SEGINI NUMNOE
  119. LCOMP=1
  120. C
  121. C======================================================================
  122. C Phase de raffinement 3D
  123. C======================================================================
  124. C=======================================
  125. C A) Boucle sur les elements a raffiner !
  126. C=======================================
  127. DO 6 IARAF=1,MELVA2.VELCHE(/2)
  128. IF (MELVA2.VELCHE(1,IARAF).NE.1) THEN
  129. DO 1 IJKL=1,IPT2.NUM(/1)
  130. IPT4.NUM(IJKL,NBELEM)=IPT2.NUM(IJKL,IARAF)
  131. 1 CONTINUE
  132. NBELEM=NBELEM-1
  133. GOTO 6
  134. ENDIF
  135. JCOMPT=0
  136. JPOS5=LPOS5(IPT2.ITYPEL)
  137. C
  138. C==========================================================
  139. C B) Boucle sur les noeuds a creer pour raffiner l'element !
  140. C==========================================================
  141. DO 4 I=1,NBRAF(IPT2.ITYPEL)
  142. JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1)
  143. JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1)
  144. JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT
  145. LTYPNO=JTYPNO(JPOS5-1+I)
  146. c write(*,*) 'LTYPNO', LTYPNO
  147. C
  148. C-------------------------------
  149. C ** B.1 / On est sur une arete !
  150. C-------------------------------
  151. IF (LTYPNO.EQ.0) THEN
  152. NPTA=IPT2.NUM(LISNOE(JLISN),IARAF)
  153. NPTB=IPT2.NUM(LISNOE(JLISN+1),IARAF)
  154. NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
  155. NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
  156. DO 2 K=1,MAX(1,KARPOS(NMIN))
  157. IF (KARETE(NMIN,K).EQ.NMAX) NEXIST=K
  158. 2 CONTINUE
  159. IF (KMILIE(NMIN,NEXIST).GT.0) THEN
  160. c si le noeud existe deja on rempli numnoe avec Kmilie
  161. NUMNOE(I)=KMILIE(NMIN,NEXIST)
  162. JCOMPT=JCOMPT+JLONG
  163. GOTO 4
  164. ELSE
  165. c sinon on ajoute un noeud a la fin de numnoe et on le met dans kmilie
  166. NBPT1=NBPT1+1
  167. NUMNOE(I)=NBPT1
  168. KMILIE(NMIN,NEXIST)=NBPT1
  169. ENDIF
  170.  
  171. c IF ((NBPT1.eq.10330).or.(NBPT1.eq.10341)) THEN
  172. c WRITE(*,*) 'NUMNOE(I)', NUMNOE(I)
  173. c WRITE(*,*) 'NPTA', NPTA, 'NPTB', NPTB
  174. c WRITE(*,*) 'NMIN', NMIN, 'NMAX', NMAX ,'NEXI', NEXIST
  175. c WRITE(*,*) 'KMILIE(NMIN,NEXIST)', KMILIE(NMIN,NEXIST)
  176. c WRITE(*,*) 'KARETE(NMIN,NEXIST)', KARETE(NMIN,NEXIST)
  177. c endif
  178. C
  179. C------------------------------
  180. C ** B.2 / On est sur une face !
  181. C------------------------------
  182. C Il faut identifier cette face (numero, sommets...)
  183. ELSEIF ((LTYPNO.GT.0).AND.(LTYPNO.LT.7)) THEN
  184. C => a) On initialise toutes les donnees relatives a cette face
  185. JTYPEL=IPT2.ITYPEL
  186. JLTEL2=LTEL(1,JTYPEL)-1+LTYPNO
  187. JLTEL2=LTEL(2,JTYPEL)-1+LTYPNO
  188. JLDEL1=LDEL(1,JLTEL2)
  189. JTYFAC=IWORK1(JLDEL1)
  190. JLDEL2=LDEL(2,JLTEL2)
  191. JNBSOM=NBSOM(JTYFAC)
  192. JSPOS=NSPOS(JTYFAC)
  193. C => b) On identifie les sommets de la face (n° global)
  194. SEGINI IWORK2
  195. DO 10 IAA=1,JNBSOM
  196. NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),IARAF)
  197. IWORK2(IAA)=NGLOBA
  198. 10 CONTINUE
  199. C => c) On classe ces sommets par ordre croissant (NPTA < NPTB < NPTC)
  200. NPTA=(XCOOR(/1)/(IDIM+1))+1
  201. NPTB=NPTA+1
  202. NPTC=NPTB+1
  203. DO 11 ICC=1,JNBSOM
  204. IF (IWORK2(ICC).LT.NPTA) THEN
  205. NPTC=NPTB
  206. NPTB=NPTA
  207. NPTA=IWORK2(ICC)
  208. ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  209. NPTC=NPTB
  210. NPTB=IWORK2(ICC)
  211. ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  212. NPTC=IWORK2(ICC)
  213. ENDIF
  214. 11 CONTINUE
  215. C => d) On passe ces sommets en n° locale
  216. NPTA1=NPTA
  217. NPTB1=NPTB
  218. NPTC1=NPTC
  219.  
  220. NPTA2=ICPR(NPTA,1)
  221. NPTB2=ICPR(NPTB,1)
  222. NPTC2=ICPR(NPTC,1)
  223. IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  224. NPTA=NPTA2
  225. NPTB=MIN(NPTB2,NPTC2)
  226. NPTC=MAX(NPTB2,NPTC2)
  227. ENDIF
  228. IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  229. NPTA=NPTB2
  230. NPTB=MIN(NPTA2,NPTC2)
  231. NPTC=MAX(NPTA2,NPTC2)
  232. ENDIF
  233. IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  234. NPTA=NPTC2
  235. NPTB=MIN(NPTA2,NPTB2)
  236. NPTC=MAX(NPTA2,NPTB2)
  237. ENDIF
  238. C => e) On cherche le numero de la face
  239. NEXIS2=0
  240. DO 12 IEE=1,JPLCOM(NPTA)
  241. MTMP=JPLANS(NPTA,IEE)
  242. JJ1=JNOEFA(MTMP,1)
  243. JJ2=JNOEFA(MTMP,2)
  244. JJ3=JNOEFA(MTMP,3)
  245. IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
  246. NEXIS2=IEE
  247. ENDIF
  248. 12 CONTINUE
  249. JNUMFA=JPLANS(NPTA,NEXIS2)
  250. C
  251. C---------------------------------
  252. C ** B.3 / Raffinement de la face !
  253. C---------------------------------
  254. KTEST1=JPLAN3(NPTA,NEXIS2)
  255. KTEST2=NBINTE(JTYFAC)
  256. C => a) Si la face est de type QUA4 (un seul noeud a creer, au milieu)
  257. IF (JTYFAC.EQ.8) THEN
  258. IF (KTEST1.LT.KTEST2) THEN
  259. c a.1) si le noeud milieu n'existe pas on l'ajoute a la fin de numnoe
  260. c et on rempli JFARAF et JPLAN3
  261. NBPT1=NBPT1+1
  262. NUMNOE(I)=NBPT1
  263. JPLAN3(NPTA,NEXIS2)=JPLAN3(NPTA,NEXIS2)+1
  264. JFARAF(JNUMFA,1)=NBPT1
  265. ELSEIF (KTEST1.EQ.KTEST2) THEN
  266. c a.2) si le noeud milieu existe deja on rempli numnoe grace à JFARAF
  267. c et on met celui-ci a zéro
  268. NUMNOE(I)=JFARAF(JNUMFA,1)
  269. JFARAF(JNUMFA,1)=0
  270. JCOMPT=JCOMPT+JLONG
  271. GOTO 4
  272. ELSE
  273. WRITE(*,*) 'ERREUR'
  274. ENDIF
  275. c IF ((NBPT1.eq.10248).or.(NBPT1.eq.10276)) THEN
  276. c IF ((NPTA.eq.235).or.(NPTA.eq.236).or.(NPTA.eq.237)) THEN
  277. c WRITE(*,*) 'NUMNOE(I)', NUMNOE(I)
  278. c WRITE(*,*) 'NPTA', NPTA1, 'NPTB', NPTB1, 'NPTC', NPTC1
  279. c WRITE(*,*) 'NPTA', NPTA, 'NPTB', NPTB, 'NPTC', NPTC
  280. c WRITE(*,*) 'JNUMFA', JNUMFA, 'NEXI', NEXIS2
  281. c WRITE(*,*) 'JPLAN3(NPTA,NEXIS2)', JPLAN3(NPTA,NEXIS2)
  282. c WRITE(*,*) 'JFARAF(JNUMFA,1)', JFARAF(JNUMFA,1)
  283. c WRITE(*,*) 'JNOEFA(JNUMFA,1)', JNOEFA(JNUMFA,1)
  284. c WRITE(*,*) 'JNOEFA(JNUMFA,2)', JNOEFA(JNUMFA,2)
  285. c WRITE(*,*) 'JNOEFA(JNUMFA,3)', JNOEFA(JNUMFA,3)
  286. c WRITE(*,*) 'JNOEFA(JNUMFA,4)', JNOEFA(JNUMFA,4)
  287. c WRITE(*,*) 'JNOEFA(JNUMFA,5)', JNOEFA(JNUMFA,5)
  288. c endif
  289. ENDIF
  290. C => b) Si la face n'est pas de type QUA4 (donc de type TRI6 ou QUA8)
  291. IF (JTYFAC.NE.8) THEN
  292. C b.1) Si il manque des noeuds dans la face
  293. C On ajoute un noeud à la fin de numnoe
  294. IF (KTEST1.LT.KTEST2) THEN
  295. NBPT1=NBPT1+1
  296. NUMNOE(I)=NBPT1
  297. JPLAN3(NPTA,NEXIS2)=JPLAN3(NPTA,NEXIS2)+1
  298. NEXIS3=0
  299. XCO2=0.25
  300. DO 13 KBB=1,JLONG
  301. XCO1=XCOEFF(JPOS1-1+KBB)
  302. IF (XCO1.EQ.XCO2) NEXIS3=KBB
  303. 13 CONTINUE
  304. JFARAF(JNUMFA,2*(KTEST1)+1)=NBPT1
  305. IF (NEXIS3.NE.0) THEN
  306. C b.1.1) si le point est sur une arête du triangle interieur
  307. C (pour un TRI6)
  308. C si le point est sur une arête de la coie interieure
  309. C (pour un QUA8)
  310.  
  311. C on renseigne la partie paire de JFARAF qui doit déja
  312. C exister ????
  313. NGLOB=IPT2.NUM(LISNOE(JLISN-1+NEXIS3),IARAF)
  314. JFARAF(JNUMFA,2*(KTEST1)+2)=NGLOB
  315. ELSE
  316. C b.1.2) si le point est le milieu d'un QUA8
  317.  
  318. C on met la partie paire de JFARAF à zéro
  319. JFARAF(JNUMFA,2*(KTEST1)+2)=0
  320. ENDIF
  321. ELSEIF (KTEST1.EQ.KTEST2) THEN
  322. C b.2) Si il ne manque pas de noeuds dans la face
  323. NEXIS3=0
  324. NEXIS4=0
  325. XCO2=0.25
  326. DO 14 KBB=1,JLONG
  327. XCO1=XCOEFF(JPOS1-1+KBB)
  328. IF (XCO1.EQ.XCO2) NEXIS3=KBB
  329. 14 CONTINUE
  330. IF (NEXIS3.NE.0) THEN
  331. C b.2.1) si le point est sur une arête du triangle interieur
  332. C (pour un TRI6)
  333. C si le point est sur une arête de la croie interieure
  334. C (pour un QUA8)
  335. NGLOB=IPT2.NUM(LISNOE(JLISN-1+NEXIS3),IARAF)
  336. DO 15 KAA=2,JFARAF(/2),2
  337. IF (JFARAF(JNUMFA,KAA).EQ.NGLOB) NEXIS4=KAA
  338. 15 CONTINUE
  339. c on rempli numnoe grace à JFARAF et on met la partie
  340. c associée au noeud en question de celui-ci a 0
  341. NUMNOE(I)=JFARAF(JNUMFA,NEXIS4-1)
  342. JFARAF(JNUMFA,NEXIS4-1)=0
  343. ELSE
  344. C b.2.2) si le point est le milieu d'un QUA8
  345. DO 16 KAA=2,JFARAF(/2),2
  346. IF (JFARAF(JNUMFA,KAA).EQ.0) NEXIS4=KAA
  347. 16 CONTINUE
  348. c on rempli numnoe grace à JFARAF et on met la partie
  349. c associée au noeud en question de celui-ci a 0
  350. NUMNOE(I)=JFARAF(JNUMFA,NEXIS4-1)
  351. JFARAF(JNUMFA,NEXIS4-1)=0
  352. ENDIF
  353. JCOMPT=JCOMPT+JLONG
  354. GOTO 4
  355. ENDIF
  356. ENDIF
  357. C
  358. C------------------------------------------------------
  359. C ** B.4 / On est a l'interieur du volume de l'element !
  360. C------------------------------------------------------
  361. ELSEIF (LTYPNO.EQ.7) THEN
  362. NBPT1=NBPT1+1
  363. NUMNOE(I)=NBPT1
  364. ENDIF
  365. C allocation de plus de mémoire dans xcoord et xden si nécéssaire
  366. IF (NBPT1.EQ.NBPTS) THEN
  367. NBPTS=NBPTS+200
  368. SEGADJ MCOORD,XDEN
  369. ENDIF
  370. C
  371. C==============================
  372. C C) Creation du nouveau point !
  373. C==============================
  374. C On continue ici que lorsque l'on doit creer un nouveau point
  375. XPT=0.
  376. YPT=0.
  377. ZPT=0.
  378. XDEN1=0.D0
  379. XDENMIN = 1000.
  380. DO 3 J=1,JLONG
  381. NGLOB=IPT2.NUM(LISNOE(JLISN-1+J),IARAF)
  382. XINI=XCOOR((NGLOB-1)*(IDIM+1)+1)
  383. YINI=XCOOR((NGLOB-1)*(IDIM+1)+2)
  384. ZINI=XCOOR((NGLOB-1)*(IDIM+1)+3)
  385. XPT=XPT+XINI*XCOEFF(JPOS1-1+J)
  386. YPT=YPT+YINI*XCOEFF(JPOS1-1+J)
  387. ZPT=ZPT+ZINI*XCOEFF(JPOS1-1+J)
  388.  
  389. XDEN1=XDEN1+XDEN(NGLOB)*XCOEFF(JPOS1-1+J)
  390. IF (XDEN(NGLOB).LT.XDENMIN) XDENMIN = XDEN(NGLOB)
  391. c IF (NBPT1.eq.10219) THEN
  392. c write (*,*) 'nglob', nglob ,'nbpt1', nbpt1
  393. c write (*,*) 'xpt ', xpt
  394. c write (*,*) 'ypt ', ypt
  395. c write (*,*) 'zpt ', zpt
  396. c write (*,*) 'xden1 ', xden1
  397. c write (*,*) 'xcoeff1 ', XCOEFF(JPOS1-1+J)
  398. c endif
  399. 3 CONTINUE
  400.  
  401. c seuil de densité a xdenmin
  402. IF (XDEN1.LT.XDENMIN) XDEN1 = XDENMIN
  403.  
  404. XCOOR((NBPT1-1)*(IDIM+1)+1)=XPT
  405. XCOOR((NBPT1-1)*(IDIM+1)+2)=YPT
  406. XCOOR((NBPT1-1)*(IDIM+1)+3)=ZPT
  407. XCOOR((NBPT1-1)*(IDIM+1)+4)=XDEN1
  408. XDEN(NBPT1)=XDEN1
  409. JCOMPT=JCOMPT+JLONG
  410.  
  411. C======================================================================
  412. 4 CONTINUE
  413. C======================================================================
  414.  
  415. JPOS4=LPOS4(IPT2.ITYPEL)
  416. C
  417. C Cas des tetraedres
  418. LONDIA=0
  419. JTYPEL=IPT2.ITYPEL
  420. IF (JTYPEL.EQ.23) THEN
  421. LPT6=NUMNOE(6)
  422. LPT5=NUMNOE(5)
  423. LPT8=NUMNOE(8)
  424. LPT10=NUMNOE(10)
  425. XPT6=XCOOR((LPT6-1)*(IDIM+1)+1)
  426. XPT5=XCOOR((LPT5-1)*(IDIM+1)+1)
  427. XPT8=XCOOR((LPT8-1)*(IDIM+1)+1)
  428. XPT10=XCOOR((LPT10-1)*(IDIM+1)+1)
  429. YPT6=XCOOR((LPT6-1)*(IDIM+1)+2)
  430. YPT5=XCOOR((LPT5-1)*(IDIM+1)+2)
  431. YPT8=XCOOR((LPT8-1)*(IDIM+1)+2)
  432. YPT10=XCOOR((LPT10-1)*(IDIM+1)+2)
  433. ZPT6=XCOOR((LPT6-1)*(IDIM+1)+3)
  434. ZPT5=XCOOR((LPT5-1)*(IDIM+1)+3)
  435. ZPT8=XCOOR((LPT8-1)*(IDIM+1)+3)
  436. ZPT10=XCOOR((LPT10-1)*(IDIM+1)+3)
  437. LON68=SQRT((XPT8-XPT6)**2+(YPT8-YPT6)**2+(ZPT8-ZPT6)**2)
  438. LON510=SQRT((XPT10-XPT5)**2+(YPT10-YPT5)**2+(ZPT10-ZPT5)**2)
  439. IF (LON510.LT.LON68) LONDIA=4*4
  440. ENDIF
  441. IF (JTYPEL.EQ.24) THEN
  442. LPT2=NUMNOE(2)
  443. LPT4=NUMNOE(4)
  444. LPT7=NUMNOE(7)
  445. LPT9=NUMNOE(9)
  446. XPT2=XCOOR((LPT2-1)*(IDIM+1)+1)
  447. XPT4=XCOOR((LPT4-1)*(IDIM+1)+1)
  448. XPT7=XCOOR((LPT7-1)*(IDIM+1)+1)
  449. XPT9=XCOOR((LPT9-1)*(IDIM+1)+1)
  450. YPT2=XCOOR((LPT2-1)*(IDIM+1)+2)
  451. YPT4=XCOOR((LPT4-1)*(IDIM+1)+2)
  452. YPT7=XCOOR((LPT7-1)*(IDIM+1)+2)
  453. YPT9=XCOOR((LPT9-1)*(IDIM+1)+2)
  454. ZPT2=XCOOR((LPT2-1)*(IDIM+1)+3)
  455. ZPT4=XCOOR((LPT4-1)*(IDIM+1)+3)
  456. ZPT7=XCOOR((LPT7-1)*(IDIM+1)+3)
  457. ZPT9=XCOOR((LPT9-1)*(IDIM+1)+3)
  458. LON47=SQRT((XPT7-XPT4)**2+(YPT7-YPT4)**2+(ZPT7-ZPT4)**2)
  459. LON29=SQRT((XPT9-XPT2)**2+(YPT9-YPT2)**2+(ZPT9-ZPT2)**2)
  460. IF (LON29.LT.LON47) LONDIA=4*10
  461. ENDIF
  462. C
  463. C===================================
  464. C D) Creation des nouveaux elements !
  465. C===================================
  466. C On remplit la portion de IPT4 relative aux elements crees a partir
  467. C de la division de l'element IARAF (indice de boucle 1).
  468. C Cette portion de IPT4 contient les colonnes dont la valeur s'etend
  469. C de INELM*(LCOMP-1)+1 a INELM*LCOMP.
  470. DO 5 J=1,INELM
  471. DO 5 I=1,IPT4.NUM(/1)
  472. NTEMP=LIELM(JPOS4-1+NBNN*(J-1)+I)
  473. IF (((JTYPEL.EQ.23).OR.(JTYPEL.EQ.24)).AND.(J.GT.4)) THEN
  474. NTEMP=LIELM(JPOS4-1+NBNN*(J-1)+I+LONDIA)
  475. ENDIF
  476. IF (NTEMP.GT.NBNN) THEN
  477. IPT4.NUM(I,INELM*(LCOMP-1)+J)=NUMNOE(NTEMP-NBNN)
  478. ELSE
  479. IPT4.NUM(I,INELM*(LCOMP-1)+J)=IPT2.NUM(NTEMP,IARAF)
  480. ENDIF
  481. 5 CONTINUE
  482. LCOMP=LCOMP+1
  483. 6 CONTINUE
  484. C
  485. NBPTS=NBPT1
  486. C SEGADJ MCOORD,XDEN
  487. C
  488. C=======================================================================
  489. C Preparation du maillage de relations
  490. C=======================================================================
  491. SEGSUP KARET2
  492. NBNDS=KARETE(/1)
  493. NCOL=KARETE(/2)
  494. SEGINI KARET2
  495. C==================================
  496. C A) Relations dues aux faces (3D) !
  497. C==================================
  498. C Tous les noeuds qui restent dans JFARAF sont a creer en tant que
  499. C relations de conformite ou des bords du domaine
  500.  
  501. C-----------------------------------------------------------------------
  502. C ** A.1 / Prise en compte des relation
  503. C d'arretes pour les arretes de chaque faces
  504. C-----------------------------------------------------------------------
  505.  
  506. C-----------------------------------------------------------------------
  507. C** A.1.1/ faces raffinée
  508. C-----------------------------------------------------------------------
  509.  
  510. c write (*,*) 'relations dues au faces'
  511. JTYPEL=IPT2.ITYPEL
  512. JLTEL1=LTEL(1,JTYPEL)
  513. JLTEL2=LTEL(2,JTYPEL)
  514. IF (JTYPEL.NE.14) GOTO 141
  515. c Boucle sur les elements de l'ancien maillage
  516. DO 121 KELEM1=1,IPT2.NUM(/2)
  517. c Boucle sur les faces de l'ancien maillage
  518. DO 122 IBB1=1,JLTEL1
  519. JLDEL1=LDEL(1,JLTEL2-1+IBB1)
  520. JTYFAC=IWORK1(JLDEL1)
  521. JLDEL2=LDEL(2,JLTEL2-1+IBB1)
  522. JNBSOM=NBSOM(JTYFAC)
  523. JSPOS=NSPOS(JTYFAC)
  524. SEGINI IWORK2
  525. C
  526. C 1/ trouver le numéro de la face
  527. DO 123 IAA=1,JNBSOM
  528. NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),KELEM1)
  529. IWORK2(IAA)=NGLOBA
  530. 123 CONTINUE
  531. C
  532. C 1.1/ Classement des 3 sommets par ordre croissant de n° globale
  533. NPTA=(XCOOR(/1)/(IDIM+1))+1
  534. NPTB=NPTA+1
  535. NPTC=NPTB+1
  536. DO 124 ICC=1,JNBSOM
  537. IF (IWORK2(ICC).LT.NPTA) THEN
  538. NPTC=NPTB
  539. NPTB=NPTA
  540. NPTA=IWORK2(ICC)
  541. ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  542. NPTC=NPTB
  543. NPTB=IWORK2(ICC)
  544. ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  545. NPTC=IWORK2(ICC)
  546. ENDIF
  547. 124 CONTINUE
  548. C
  549. C1.2/ Classement des 3 sommets par ordre croissant de n° locale
  550. NPTA2=ICPR(NPTA,1)
  551. NPTB2=ICPR(NPTB,1)
  552. NPTC2=ICPR(NPTC,1)
  553. IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  554. NPTA=NPTA2
  555. NPTB=MIN(NPTB2,NPTC2)
  556. NPTC=MAX(NPTB2,NPTC2)
  557. ENDIF
  558. IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  559. NPTA=NPTB2
  560. NPTB=MIN(NPTA2,NPTC2)
  561. NPTC=MAX(NPTA2,NPTC2)
  562. ENDIF
  563. IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  564. NPTA=NPTC2
  565. NPTB=MIN(NPTA2,NPTB2)
  566. NPTC=MAX(NPTA2,NPTB2)
  567. ENDIF
  568. C
  569. C 1.3/ Recherche du numero de la face
  570. NEXIS2=0
  571. DO 125 IEE=1,JPLCOM(NPTA)
  572. MTMP=JPLANS(NPTA,IEE)
  573. JJ1=JNOEFA(MTMP,1)
  574. JJ2=JNOEFA(MTMP,2)
  575. JJ3=JNOEFA(MTMP,3)
  576. IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
  577. NEXIS2=IEE
  578. ENDIF
  579. 125 CONTINUE
  580.  
  581. JNUMF1=JPLANS(NPTA,NEXIS2)
  582.  
  583. C
  584. C 2/ Exclusion des faces non voulues
  585.  
  586. C on exlu les faces où un noeud existait déja
  587. IF (JFARAF(JNUMF1,1).EQ.0) GOTO 122
  588. C on exlu les faces de bords
  589. IF (JNOEFA(JNUMF1,5).EQ.0) GOTO 122
  590.  
  591. C 3/ Recherche des arretes de la faces pour y mettre des relations
  592.  
  593. IF (JTYPEL.NE.14) GOTO 137
  594. NBAR1=0
  595. DO 138 ISS=1, (JNBSOM-1)
  596. DO 138 JSS=(ISS+1), JNBSOM
  597.  
  598. NMIN=MIN(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1))
  599. NMAX=MAX(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1))
  600.  
  601. DO 139 IRR=1,MAX(1,KARPOS(NMIN))
  602. IF (KARETE(NMIN,IRR).EQ.NMAX) THEN
  603. NBAR1=NBAR1+1
  604. KARET2(NMIN,IRR)=KARET2(NMIN,IRR)+1
  605. c write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR)
  606. c IF (JFARAF(JNUMF1,1).EQ.15748) THEN
  607. c write (*,*) ISS,JSS
  608. c write (*,*) '!!! NMIN', NMIN, 'NMAX', NMAX, '!!!'
  609. c write (*,*) IWORK2(ISS), IWORK2(JSS)
  610. c write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR)
  611. c write (*,*) 'KARET2(NMIN,IRR)',KARET2(NMIN,IRR)
  612. c ENDIF
  613. ENDIF
  614. 139 CONTINUE
  615. 138 CONTINUE
  616. c WRITE (*,*) 'Nbarrete', NBAR1
  617. c IF (NBAR1.NE.4) THEN
  618. c write (*,*) IWORK2(1), IWORK2(2)
  619. c write (*,*) IWORK2(3), IWORK2(4)
  620. c ENDIF
  621. 137 CONTINUE
  622. 122 CONTINUE
  623. 121 CONTINUE
  624.  
  625. 141 CONTINUE
  626.  
  627. C-----------------------------------------------------------------------
  628. C** A.1.2/ relations de faces conservées dans le maillage précédent
  629. C-----------------------------------------------------------------------
  630.  
  631. SEGINI IPT5
  632. DO 670 MIE=1,MAX(1,IPT3.LISOUS(/1))
  633. IF (IPT3.LISOUS(/1).EQ.0) THEN
  634. IPT5=IPT3
  635. ELSE
  636. IPT5=IPT3.LISOUS(MIE)
  637. ENDIF
  638. SEGACT IPT5
  639. c write (*,*) 'ipt3 raff 3D' , ipt3, IPT3.ITYPEL, IPT3.ICOLOR(1)
  640. c write (*,*) 'ipt5 raff 3D' , ipt5, IPT5.ITYPEL, IPT5.ICOLOR(1)
  641. IF (IPT5.ITYPEL.NE.48) GOTO 670
  642. IF ((IPT5.ICOLOR(1).NE.3)) GOTO 670
  643. LIGNUM=IPT2.NUM(/1)
  644. JNBSOM=4
  645. SEGINI IWORK2
  646.  
  647. C
  648. C 1/ trouver le numéro de la face
  649. C 1.1/ Classement des 3 sommets par ordre croissant (num globale)
  650. DO 671,MAA=1,IPT5.NUM(/2)
  651. DO 673 MBB=1,JNBSOM
  652. c write (*,*) 'MAA' , MAA, 'MBB', MBB
  653. IWORK2(MBB)=IPT5.NUM(MBB+1,MAA)
  654.  
  655. 673 CONTINUE
  656. NPTA=(XCOOR(/1)/(IDIM+1))+1
  657. NPTB=NPTA+1
  658. NPTC=NPTB+1
  659. DO 672 ICC=1,JNBSOM
  660. IF (IWORK2(ICC).LT.NPTA) THEN
  661. NPTC=NPTB
  662. NPTB=NPTA
  663. NPTA=IWORK2(ICC)
  664. ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  665. NPTC=NPTB
  666. NPTB=IWORK2(ICC)
  667. ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  668. NPTC=IWORK2(ICC)
  669. ENDIF
  670. 672 CONTINUE
  671. C
  672. C 1.2/ Classement des 3 sommets par ordre croissant (num locale)
  673. NPTA1=NPTA
  674. NPTB1=NPTB
  675. NPTC1=NPTC
  676. NPTA2=ICPR(NPTA,1)
  677. NPTB2=ICPR(NPTB,1)
  678. NPTC2=ICPR(NPTC,1)
  679. IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  680. NPTA=NPTA2
  681. NPTB=MIN(NPTB2,NPTC2)
  682. NPTC=MAX(NPTB2,NPTC2)
  683. ENDIF
  684. IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  685. NPTA=NPTB2
  686. NPTB=MIN(NPTA2,NPTC2)
  687. NPTC=MAX(NPTA2,NPTC2)
  688. ENDIF
  689. IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  690. NPTA=NPTC2
  691. NPTB=MIN(NPTA2,NPTB2)
  692. NPTC=MAX(NPTA2,NPTB2)
  693. ENDIF
  694. C
  695. C 1.3/ Recherche du numero de la face
  696. NEXIS2=0
  697. DO 674 IEE=1,JPLCOM(NPTA)
  698. MTMP=JPLANS(NPTA,IEE)
  699. JJ1=JNOEFA(MTMP,1)
  700. JJ2=JNOEFA(MTMP,2)
  701. JJ3=JNOEFA(MTMP,3)
  702. IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
  703. NEXIS2=IEE
  704. ENDIF
  705. 674 CONTINUE
  706. JNUMF1=JPLANS(NPTA,NEXIS2)
  707. c IF (MAA.EQ.5) THEN
  708. c write (*,*) 'MAA', MAA, 'nexis2', nexis2
  709. c write (*,*) IWORK2(1), IWORK2(2), IWORK2(3), IWORK2(4)
  710. c write (*,*) 'JFARAF(JNUMF1,1)', JFARAF(JNUMF1,1)
  711. c write (*,*) 'JNOEFA(JNUMF1,5)',JNOEFA(JNUMF1,5)
  712. c ENDIF
  713. C
  714. C 2/ Exclusion des faces non voulues
  715. IF (NEXIS2.EQ.0) THEN
  716. c WRITE (*,*)'FACE NON TROUVE SURE:', MAA
  717. GOTO 671
  718. ENDIF
  719. C on exlu les faces où un noeud existait déja
  720. IF (JFARAF(JNUMF1,1).EQ.0) GOTO 671
  721.  
  722. C on exlu les faces de bords
  723. IF (JNOEFA(JNUMF1,5).EQ.0) GOTO 671
  724. NBAR1=0
  725. C
  726. C 3/ Recherche des arretes de la faces pour y mettre des relations
  727.  
  728. DO 135 ISS=1, (JNBSOM-1)
  729. DO 135 JSS=(ISS+1), JNBSOM
  730.  
  731. NMIN=MIN(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1))
  732. NMAX=MAX(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1))
  733.  
  734. DO 136 IRR=1,MAX(1,KARPOS(NMIN))
  735. IF (KARETE(NMIN,IRR).EQ.NMAX) THEN
  736. NBAR1=NBAR1+1
  737. KARET2(NMIN,IRR)=KARET2(NMIN,IRR)+1
  738. c write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR)
  739. c IF (MAA.EQ.5) THEN
  740. c write (*,*) ISS,JSS
  741. c write (*,*) '!!! NMIN', NMIN, 'NMAX', NMAX, '!!!'
  742. c write (*,*) IWORK2(ISS), IWORK2(JSS)
  743. c write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR)
  744. c write (*,*) 'KARET2(NMIN,IRR)',KARET2(NMIN,IRR)
  745. c ENDIF
  746. ENDIF
  747. 136 CONTINUE
  748. 135 CONTINUE
  749.  
  750. 671 CONTINUE
  751.  
  752. 670 CONTINUE
  753. SEGDES IPT5
  754. C---------------------------------
  755. C ** A.2 / Initialisation de IPT5 !
  756. C---------------------------------
  757.  
  758.  
  759.  
  760. C 1/ Comptage du nombre de noeuds soumis a des relations
  761. NBRELA=0
  762. DO 114 IHF=1,JFARAF(/1)
  763. C on exlu les faces de bords
  764. IF (JNOEFA(IHF,5).EQ.0) GOTO 114
  765. C on exlu les faces où un noeud existait déja
  766. IF (JFARAF(IHF,1).EQ.0) GOTO 114
  767. DO 113 JF=1,JFARAF(/2),2
  768. IF (JFARAF(IHF,JF).NE.0) NBRELA=NBRELA+1
  769. 113 CONTINUE
  770. 114 CONTINUE
  771. C
  772. C 2/ Creation de IPT5
  773. IF (NBRELA.EQ.0) THEN
  774. IPT5=0
  775. GOTO 210
  776. ENDIF
  777. C IF (IPT2.ITYPEL.EQ.14) NBNN=4
  778. C IF (IPT2.ITYPEL.EQ.15) NBNN=8
  779. C IF (IPT2.ITYPEL.EQ.17) NBNN=5
  780. C IF (IPT2.ITYPEL.EQ.24) NBNN=5
  781. NBNN=10
  782. NBELEM=NBRELA
  783. NBSOUS=0
  784. NBREF=0
  785. SEGINI IPT5
  786. IPT5.ITYPEL=48
  787. C
  788. C 3/ Renseignement des noeuds support des relations
  789. NBRELA=0
  790. DO 116 IPF=1,JFARAF(/1)
  791. IF (JNOEFA(IPF,5).EQ.0) GOTO 116
  792. IF (JFARAF(IPF,1).EQ.0) GOTO 116
  793. DO 115 JF=1,JFARAF(/2),2
  794. IF (JFARAF(IPF,JF).EQ.0) GOTO 115
  795. NBRELA=NBRELA+1
  796. IPT5.NUM(1,NBRELA)=JFARAF(IPF,JF)
  797. 115 CONTINUE
  798. 116 CONTINUE
  799. c write (*,*) 'NBRELA FACE ', NBRELA
  800. C
  801. C-----------------------------------------------------
  802. C ** A.3 / Recherche des noeuds formant les relations !
  803. C-----------------------------------------------------
  804. C 1/ Boucle sur l'ensemble des noeuds a creer
  805. DO 200 IARAF=1,MELVA2.VELCHE(/2)
  806. IF (MELVA2.VELCHE(1,IARAF).NE.1) GOTO 200
  807. JCOMPT=0
  808. JPOS5=LPOS5(IPT2.ITYPEL)
  809. DO 190 I=1,NBRAF(IPT2.ITYPEL)
  810. JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1)
  811. JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1)
  812. JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT
  813. LTYPNO=JTYPNO(JPOS5-1+I)
  814. IF ((LTYPNO.EQ.0).OR.(LTYPNO.EQ.7)) GOTO 189
  815. C
  816. C 2/ Preparation pour trouver le noeud et la face en question
  817. JTYPEL=IPT2.ITYPEL
  818. JLTEL2=LTEL(2,JTYPEL)-1+LTYPNO
  819. JLDEL1=LDEL(1,JLTEL2)
  820. JTYFAC=IWORK1(JLDEL1)
  821. JLDEL2=LDEL(2,JLTEL2)
  822. JNBSOM=NBSOM(JTYFAC)
  823. JSPOS=NSPOS(JTYFAC)
  824. SEGINI IWORK2
  825. C
  826. C 3/ Classement des 3 sommets par ordre croissant de n° globale
  827. DO 100 IAA=1,JNBSOM
  828. NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),IARAF)
  829. IWORK2(IAA)=NGLOBA
  830. 100 CONTINUE
  831. NPTA=(XCOOR(/1)/(IDIM+1))+1
  832. NPTB=NPTA+1
  833. NPTC=NPTB+1
  834. DO 110 ICC=1,JNBSOM
  835. IF (IWORK2(ICC).LT.NPTA) THEN
  836. NPTC=NPTB
  837. NPTB=NPTA
  838. NPTA=IWORK2(ICC)
  839. ELSEIF (IWORK2(ICC).LT.NPTB) THEN
  840. NPTC=NPTB
  841. NPTB=IWORK2(ICC)
  842. ELSEIF (IWORK2(ICC).LT.NPTC) THEN
  843. NPTC=IWORK2(ICC)
  844. ENDIF
  845. 110 CONTINUE
  846. C
  847. C 4/ Classement des 3 sommets par ordre croissant de n° locale
  848. NPTA2=ICPR(NPTA,1)
  849. NPTB2=ICPR(NPTB,1)
  850. NPTC2=ICPR(NPTC,1)
  851. IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
  852. NPTA=NPTA2
  853. NPTB=MIN(NPTB2,NPTC2)
  854. NPTC=MAX(NPTB2,NPTC2)
  855. ENDIF
  856. IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
  857. NPTA=NPTB2
  858. NPTB=MIN(NPTA2,NPTC2)
  859. NPTC=MAX(NPTA2,NPTC2)
  860. ENDIF
  861. IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
  862. NPTA=NPTC2
  863. NPTB=MIN(NPTA2,NPTB2)
  864. NPTC=MAX(NPTA2,NPTB2)
  865. ENDIF
  866. C
  867. C 5/ Recherche du numero de la face
  868. NEXIS2=0
  869. DO 120 IEE=1,JPLCOM(NPTA)
  870. MTMP=JPLANS(NPTA,IEE)
  871. JJ1=JNOEFA(MTMP,1)
  872. JJ2=JNOEFA(MTMP,2)
  873. JJ3=JNOEFA(MTMP,3)
  874. IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
  875. NEXIS2=IEE
  876. ENDIF
  877. 120 CONTINUE
  878. JNUMFA=JPLANS(NPTA,NEXIS2)
  879.  
  880. C
  881. C 6/ Recherche du numero global du point
  882. IF (JNOEFA(JNUMFA,5).EQ.0) GOTO 189
  883. IF (JTYFAC.EQ.8) INOEGL=JFARAF(JNUMFA,1)
  884. IF (JTYFAC.NE.8) THEN
  885. NEXIS3=0
  886. NEXIS4=0
  887. XCO2=0.25
  888. DO 140 KBB=1,JLONG
  889. XCO1=XCOEFF(JPOS1-1+KBB)
  890. IF (XCO1.EQ.XCO2) NEXIS3=KBB
  891. 140 CONTINUE
  892. IF (NEXIS3.NE.0) THEN
  893. NGLOB=IPT2.NUM(LISNOE(JLISN-1+NEXIS3),IARAF)
  894. DO 150 KAA=2,JFARAF(/2),2
  895. IF (JFARAF(JNUMFA,KAA).EQ.NGLOB) NEXIS4=KAA
  896. 150 CONTINUE
  897. INOEGL=JFARAF(JNUMFA,NEXIS4-1)
  898. ELSE
  899. DO 160 KAA=2,JFARAF(/2),2
  900. IF (JFARAF(JNUMFA,KAA).EQ.0) NEXIS4=KAA
  901. 160 CONTINUE
  902. INOEGL=JFARAF(JNUMFA,NEXIS4-1)
  903. ENDIF
  904. ENDIF
  905.  
  906. C
  907. C------------------------------
  908. C ** A.4 / Remplissage de IPT5 !
  909. C------------------------------
  910. C 1/ Recherche de la position du point dans IPT5
  911. NEXIS5=0
  912. DO 170 IGG=1,IPT5.NUM(/2)
  913. IF (INOEGL.EQ.IPT5.NUM(1,IGG)) NEXIS5=IGG
  914. 170 CONTINUE
  915. IF (NEXIS5.EQ.0) GOTO 189
  916. IF (INOEGL.EQ.139) THEN
  917. c write (*,*) INOEGL, IWORK2(1), IWORK2(2), IWORK2(3)
  918. ENDIF
  919.  
  920. C
  921. C 2/ Renseignement des points formant les relations
  922. DO 180 IHH=1,JLONG
  923. IPT5.NUM(1+IHH,NEXIS5)=IPT2.NUM(LISNOE(JLISN-1+IHH),IARAF)
  924. 180 CONTINUE
  925. IF (JLONG.EQ.4) IPT5.NUM(10,NEXIS5)=3
  926. IF (JLONG.EQ.5) IPT5.NUM(10,NEXIS5)=4
  927. IF (JLONG.EQ.8) THEN
  928. IF (JPOS1.EQ.16) IPT5.NUM(10,NEXIS5)=5
  929. IF (JPOS1.EQ.24) IPT5.NUM(10,NEXIS5)=6
  930. ENDIF
  931. 189 CONTINUE
  932. JCOMPT=JCOMPT+JLONG
  933. 190 CONTINUE
  934. 200 CONTINUE
  935. 210 CONTINUE
  936. C
  937. C===================================
  938. C B) Relations dues aux aretes (2D) !
  939. C===================================
  940. C On cree un maillage IPT6 contenant tous les noeuds soumis a des
  941. C relations
  942. C
  943.  
  944. c ILPL=LPL(IPT4.ITYPEL)
  945. c ILPT=LPT(IPT4.ITYPEL)
  946. c DO 317 J=1,IPT4.NUM(/2)
  947. c DO 317 K=1,ILPL*2-1,2
  948. c NPTA=IPT4.NUM(KSEGM(ILPT+K-1),J)
  949. c NPTB=IPT4.NUM(KSEGM(ILPT+K),J)
  950. c
  951. c if ((NPTB.eq.129).and.(NPTA.eq.9087)) then
  952. c write(*,*) '!ici!',NBPT0
  953. c endif
  954. c IF((NPTA.GT.NBPT0).OR.(NPTB.GT.NBPT0)) THEN
  955. c GOTO 317
  956. c ENDIF
  957. c NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
  958. c NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
  959. c if ((NMIN.eq.235).and.(NMAX.eq.236)) then
  960. c write(*,*) 'nmin', nmin, 'nmax', nmax
  961. c write(*,*) 'NPTA',NPTA, 'NPTB', NPTB
  962. c endif
  963. c NEXIST=0
  964. c DO 316 I=1,MAX(1,KARPOS(NMIN))
  965. c IF (KARETE(NMIN,I).EQ.NMAX) THEN
  966. c KARET2(NMIN,I)=KARET2(NMIN,I)+1
  967. c
  968. c IF (KMILIE(NMIN,I).gt.0) then
  969. c prise en compte des arretes multiples sur le nouveau maillage
  970. c dans l'exemple ci dessous où les noeuds C et D sont des hanging nodes
  971. c L'arette [C B] n'a pas encore été prise en compte dans karet2 c'est
  972. c ce qui est fait ici.
  973. c
  974. c | |
  975. c | |
  976. c |A C D |B
  977. c x---------X----X----X
  978. c | | | |
  979. c | | | |
  980.  
  981. c1/ Test sur la premiere arete : NMIN-KMILIE(NMIN,I)
  982. c NMILI = ICPR(KMILIE(NMIN,I),1)
  983. c NEXIS2=0
  984. c NMIN2=MIN(NMIN,NMILI)
  985. c NMAX2=MAX(NMIN,NMILI)
  986. c DO 318 I2=1,KMILIE(/2)
  987. c IF (KMILIE(NMIN2,I2).GT.0) THEN
  988. c KARET2(NMIN2,I2)=KARET2(NMIN2,I2)+1
  989. c ENDIF
  990. c 318 CONTINUE
  991. C
  992. C 2/ Test sur la deuxieme arete : NMAX-NPTC
  993. c NEXIS2=0
  994. c NMIN2=MIN(NMAX,NMILI)
  995. c NMAX2=MAX(NMAX,NMILI)
  996. c DO 319 I2=1,KMILIE(/2)
  997. c IF (KMILIE(NMIN2,I2).GT.0) THEN
  998. c KARET2(NMIN2,I2)=KARET2(NMIN2,I2)+1
  999. c ENDIF
  1000. c 319 CONTINUE
  1001. c ENDIF
  1002. c ENDIF
  1003. c 316 CONTINUE
  1004. c 317 CONTINUE
  1005.  
  1006. C 1/ Comptage du nombre de noeuds soumis a des relations
  1007. NBELEM=0
  1008. DO 27 J=1,KMILIE(/2)
  1009. DO 27 I=1,KMILIE(/1)
  1010.  
  1011. c IF (KMILIE(I,J).eq.10330) write(*,*) '!kmilie', I, J, '=10330!'
  1012. c IF (KMILIE(I,J).eq.10330) write(*,*) 'KARET2(I,J)', KARET2(I,J)
  1013. c IF (KMILIE(I,J).eq.10341) write(*,*) '!kmilie', I, J, '=10341!'
  1014. c IF (KMILIE(I,J).eq.10341) write(*,*) 'KARET2(I,J)', KARET2(I,J)
  1015. IF (KARET2(I,J).EQ.0) GOTO 27
  1016. IF (KMILIE(I,J).GT.0) NBELEM=NBELEM+1
  1017.  
  1018. 27 CONTINUE
  1019. c write (*,*) 'NBELEM', NBELEM
  1020. C
  1021. C 2/ Creation de IPT6
  1022. IPT6=0
  1023. IF (NBELEM.EQ.0) GOTO 999
  1024. NBNN=5
  1025. NBREF=0
  1026. NBSOUS=0
  1027. SEGINI IPT6
  1028. IPT6.ITYPEL=48
  1029. C
  1030. C 3/ Renseignement des noeuds support des relations
  1031. DO 28 J=1,KMILIE(/2)
  1032. DO 28 I=1,KMILIE(/1)
  1033. IF (KARET2(I,J).EQ.0) GOTO 28
  1034. IF (KMILIE(I,J).GT.0) THEN
  1035. NBREF=NBREF+1
  1036. IPT6.NUM(1,NBREF)=KMILIE(I,J)
  1037. ENDIF
  1038. 28 CONTINUE
  1039. C WRITE (*,*) 'nombre de noeuds supports de rela seg', NBREF
  1040. C
  1041. C 4/ Recherche des noeuds formant les relations
  1042. DO 24 IARAF=1,MELVA2.VELCHE(/2)
  1043. IF (MELVA2.VELCHE(1,IARAF).NE.1) GOTO 24
  1044. JCOMPT=0
  1045. JPOS5=LPOS5(IPT2.ITYPEL)
  1046. DO 23 I=1,NBRAF(IPT2.ITYPEL)
  1047. JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1)
  1048. JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1)
  1049. JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT
  1050. LTYPNO=JTYPNO(JPOS5-1+I)
  1051. IF (LTYPNO.NE.0) GOTO 22
  1052. NPTA=IPT2.NUM(LISNOE(JLISN),IARAF)
  1053. NPTB=IPT2.NUM(LISNOE(JLISN+1),IARAF)
  1054. NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
  1055. NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
  1056. DO 29 K=1,MAX(1,KARPOS(NMIN))
  1057. IF (KARETE(NMIN,K).EQ.NMAX) NEXIST=K
  1058. 29 CONTINUE
  1059. IF (KARET2(NMIN,NEXIST).EQ.0) GOTO 22
  1060. IF (KMILIE(NMIN,NEXIST).EQ.0) GOTO 22
  1061. NEXIS5=0
  1062. DO 20 MM=1,IPT6.NUM(/2)
  1063. INOEGL=KMILIE(NMIN,NEXIST)
  1064. INRELA=IPT6.NUM(1,MM)
  1065. IF (INOEGL.EQ.INRELA) NEXIS5=MM
  1066. 20 CONTINUE
  1067. C
  1068. C 5/ Renseignement des noeuds formant les relations
  1069. DO 21 IHH=1,JLONG
  1070. IPT6.NUM(1+IHH,NEXIS5)=IPT2.NUM(LISNOE(JLISN-1+IHH),IARAF)
  1071. 21 CONTINUE
  1072. IF (JLONG.EQ.2) IPT6.NUM(5,NEXIS5)=1
  1073. IF (JLONG.EQ.3) IPT6.NUM(5,NEXIS5)=2
  1074. 22 CONTINUE
  1075. JCOMPT=JCOMPT+JLONG
  1076. 23 CONTINUE
  1077. 24 CONTINUE
  1078.  
  1079. 444 CONTINUE
  1080.  
  1081. C
  1082. C============================================
  1083. C C) Creation du maillage de relations final !
  1084. C============================================
  1085. IF (IPT5.EQ.0) THEN
  1086. IPT7=IPT6
  1087. GOTO 999
  1088. ENDIF
  1089. NBELEM=IPT5.NUM(/2)+IPT6.NUM(/2)
  1090. C NBNN=MAX(IPT5.NUM(/1),IPT6.NUM(/1))
  1091. NBNN=10
  1092. NBREF=0
  1093. NBSOUS=0
  1094. SEGINI IPT7
  1095. IPT7.ITYPEL=48
  1096. DO 42 NEO=1,IPT5.NUM(/2)
  1097. DO 42 MOR=1,10
  1098. IPT7.NUM(MOR,NEO)=IPT5.NUM(MOR,NEO)
  1099. IPT7.ICOLOR(NEO)=IPT5.ICOLOR(NEO)
  1100. 42 CONTINUE
  1101. NN5 = IPT5.NUM(/2)
  1102. DO 43 NEO=IPT5.NUM(/2)+1,IPT6.NUM(/2)+IPT5.NUM(/2)
  1103. DO 43 MOR=1,IPT6.NUM(/1)
  1104. IF (MOR.LT.IPT6.NUM(/1)) THEN
  1105. IPT7.NUM(MOR,NEO)=IPT6.NUM(MOR,NEO-NN5)
  1106. ELSE
  1107. IPT7.NUM(10,NEO)=IPT6.NUM(MOR,NEO-NN5)
  1108. ENDIF
  1109. IPT7.ICOLOR(NEO)=IPT6.ICOLOR(NEO-NN5)
  1110. 43 CONTINUE
  1111. SEGDES IPT5, IPT6
  1112.  
  1113. C=======================================================================
  1114. C Fin du programme
  1115. C=======================================================================
  1116. 999 CONTINUE
  1117. c write (*,*) 'IPT2', IPT2
  1118. c write (*,*) 'IPT3', IPT3
  1119. c write (*,*) 'IPT4', IPT4
  1120. c write (*,*) 'IPT5', IPT5
  1121. c write (*,*) 'IPT6', IPT6
  1122. c write (*,*) 'IPT7', IPT7
  1123. RETURN
  1124. END
  1125.  
  1126.  
  1127.  
  1128.  
  1129.  
  1130.  
  1131.  
  1132.  
  1133.  
  1134.  
  1135.  
  1136.  

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