Télécharger psip3d.eso

Retour à la liste

Numérotation des lignes :

  1. C PSIP3D SOURCE BP208322 19/10/04 21:15:10 10326
  2. c
  3. SUBROUTINE PSIP3D(ideux,melmai,melfis,melfro,xcrit,
  4. & msoup1,msoup2,msoup3)
  5.  
  6. ***********************************************************************
  7. c OPERATEUR : PSIP
  8.  
  9. c APPELE PAR PSIPHI DANS LE CAS 3D
  10. c
  11. c CREATION : BP, 14/03/2012
  12. c MODIF : BP, 19/12/2013 plus d'appel au zonage
  13. c BP, 2018-06 correction d'une petite erreur + //isation
  14. c
  15. ***********************************************************************
  16. IMPLICIT REAL*8(A-H,O-Z)
  17. IMPLICIT INTEGER(I-N)
  18. -INC CCOPTIO
  19. -INC SMCOORD
  20. -INC SMELEME
  21. -INC SMCHPOI
  22. -INC CCREEL
  23. -INC CCASSIS
  24. C Declaration du COMMON pour le travail en parallele
  25. COMMON/psiabc/IPARAL
  26.  
  27. segment icpr(xcoor(/1)/(idim+1))
  28. segment mxicpr
  29. real*8 xicpr(xcoor(/1)/(idim+1))
  30. endsegment
  31. segment isup(npt)
  32. segment sfis
  33. integer ifis(nbelfi,2)
  34. endsegment
  35.  
  36. c MTRI3 : toutes les infos pour decrire un ensemble de TRI3
  37. c noeuds, coordonnees, vecteur unitaire des aretes, normale, aire
  38. SEGMENT MTRI3
  39. INTEGER KNODE(3,NTRI3)
  40. REAL*8 XNODE(3,NTRI3),YNODE(3,NTRI3),ZNODE(3,NTRI3)
  41. c REAL*8 XMIMA(NTRI3),YMIMA(NTRI3),ZMIMA(NTRI3)
  42. REAL*8 VARE12(4,NTRI3),VARE23(3,NTRI3),VARE31(3,NTRI3)
  43. REAL*8 VNORM(3,NTRI3),SURFAC(NTRI3)
  44. ENDSEGMENT
  45.  
  46. c SPARAL : pour la parallelisation
  47. C NBTHRD : nombre de threads demandes
  48. c KPT3 : ...
  49. SEGMENT SPARAL
  50. INTEGER NBTHRD
  51. INTEGER IERROR(NBTHR)
  52. INTEGER KPT3,KPOVA1,KPOVA2,KPOVA3
  53. INTEGER KTRI3,KXICPR,KSUP,KFIS
  54. INTEGER KCLE,KDEUX,KMIL5
  55. REAL*8 XCRIT1
  56. ENDSEGMENT
  57.  
  58. DATA EPSI/1.D-5/
  59. INTEGER ooperm(7,2)
  60. DATA ooperm/ 2, 4, 6, 1, 3, 5, 7,
  61. $ 4, 1, 5, 2, 6, 3, 7/
  62.  
  63. EXTERNAL PSI3Di
  64. LOGICAL BTHRD
  65.  
  66.  
  67. ***********************************************************************
  68. * Initialisation, Recup, Activation
  69. ***********************************************************************
  70.  
  71. idim1=idim+1
  72. c idbug=numero du noeud pour lequel on souhaite afficher le calcul
  73. idbug = 0
  74. c idbug = 16346
  75.  
  76. c des chpoints (les msoup sont deja actifs)
  77. mpova1=0
  78. mpova2=0
  79. mpova3=0
  80. if(msoup2.gt.0) mpova2=msoup2.ipoval
  81. if(msoup1.gt.0) mpova1=msoup1.ipoval
  82. if(msoup3.gt.0) mpova3=msoup3.ipoval
  83. c activation des maillages / desactivation
  84. c melmai : fait dans psiphi / ligne 1325 ou 1329
  85. c melfis : fait ligne 64 / ligne 1274
  86. c melfro : fait (si ideux>=1) ligne 78 / ligne 356
  87. ipt3=melmai
  88. c segact,ipt3
  89. npt=ipt3.num(/2)
  90.  
  91. c variables parfois testées
  92. imil5=0
  93.  
  94.  
  95. ***********************************************************************
  96. * Preparation du maillage de la fissure
  97. ***********************************************************************
  98. meleme=melfis
  99. segact meleme
  100. * transfo du front en TRI3
  101. if(itypel.ne.4) call change(melfis,4)
  102. meleme=melfis
  103. segact meleme
  104. c rem: inutile car change laisse melfis actif
  105. nseg=num(/2)
  106.  
  107. *------ Cas ou l'on calcule phi ET psi --------------------------------
  108. if (ideux.ge.1) then
  109. * on rearrange le melfis de facon à prevenir dans ifis(iel,1) quels
  110. * element touchent le front de fissure. Ceux ci sont ordonnées de facon
  111. * que le front de fissure soit du noeud 1 vers le noeud 2 de l'element
  112. ipt2=melfro
  113. segact ipt2
  114. segini icpr
  115. nbelem2 = ipt2.num(/2)
  116. nbnode2 = ipt2.num(/1)
  117. * petite numerotation locale
  118. * icpr(noeud)>0 si le noeud appartient au front de fissure
  119. * =numero du 1er element de melfro ou on a vu ce noeud
  120. * + on en profite pour calculer l'abscisse curviligne du front si ideux.eq.2
  121. * icpr(noeud)<0 si le noeud appartient au reste du contour de la fissure
  122. * =-numero du 1er element de melcon ou on a vu ce noeud
  123. if (ideux.ge.2) then
  124. segini,mxicpr
  125. xtau=0.d0
  126. ic=ipt2.num(1,1)
  127. x11= xcoor((ic-1)*idim1 +1)
  128. y11= xcoor((ic-1)*idim1 +2)
  129. z11= xcoor((ic-1)*idim1 +3)
  130. do ia=1,nbelem2
  131. do ib=1,nbnode2
  132. ic=ipt2.num(ib,ia)
  133. if (icpr(ic).eq.0) then
  134. icpr(ic)=ia
  135. x22= xcoor((ic-1)*idim1 +1)
  136. y22= xcoor((ic-1)*idim1 +2)
  137. z22= xcoor((ic-1)*idim1 +3)
  138. dxtau= sqrt( (x22-x11)*(x22-x11)+(y22-y11)*(y22-y11)
  139. $ +(z22-z11)*(z22-z11) )
  140. xtau=xtau+dxtau
  141. xicpr(ic)=xtau
  142. x11=x22
  143. y11=y22
  144. z11=z22
  145. c if(idebug.ge.1)
  146. c $ write(6,*) 'ia,ib,ic,xicpr(ic)=',ia,ib,ic,xicpr(ic)
  147. endif
  148. enddo
  149. enddo
  150. else
  151. do ia=1,nbelem2
  152. do ib=1,nbnode2
  153. ic=ipt2.num(ib,ia)
  154. if (icpr(ic).eq.0) then
  155. icpr(ic)=ia
  156. endif
  157. enddo
  158. enddo
  159. endif
  160. * bp (23/02/2012) : calcul du contour et fin de remplissage de icpr
  161. CALL ECROBJ('MAILLAGE',melfis)
  162. CALL PRCONT
  163. CALL LIROBJ('MAILLAGE',melcon,1,IRETOU)
  164. IF (IERR.NE.0) RETURN
  165. meleme = melfis
  166. IPT6 = melcon
  167. SEGACT,meleme,IPT6
  168. nbelCON = IPT6.num(/2)
  169. nbnoCON = IPT6.num(/1)
  170. do ia=1,nbelCON
  171. do ib=1,nbnoCON
  172. ic=IPT6.num(ib,ia)
  173. c noeud qui n'a pas ete vu dans melfro et qui appartient au contour
  174. if(icpr(ic).eq.0) icpr(ic)=-ia
  175. enddo
  176. enddo
  177. * nombre de segments appartenant au contour sans le front
  178. nbelem6 = nbelCON - nbelem2
  179. * remplissage effectif de ifis contenu dans le segment sfis
  180. i2firs = ipt2.num(1,1)
  181. i2last = ipt2.num(nbnode2,nbelem2)
  182. ivui2f = 0
  183. ivui2l = 0
  184. c write(6,*) 'i2firs,i2last=',i2firs,i2last
  185. nbelfi=nseg
  186. nbelem=nseg
  187. nbnn=3
  188. nbsous=0
  189. nbref=0
  190. c segini,ipt5=meleme
  191. * bp 20/12/2011 : on propose de remplir ipt5 avec :
  192. c + des elements du front (et du contour) d abord,
  193. c + des elements avec 1 point appartenant au front (et au contour) ensuite,
  194. c + les autres elements enfin.
  195. c bp 06/03/2012 : on propose de remplir ifis avec :
  196. c ifis(iseg,1) =2 si le segment [1-2] du iseg^ieme triangle appartient
  197. c au front de fissure, =3 si en plus c'est l extrémité du front,
  198. c =1 si l un des segments appartient au contour
  199. c =-1 si l'un des noeuds appartient au front
  200. c ifis(iseg,2) =+1+2+4 si les segments [1-2] [2-3] [3-1] du iseg^ieme triangle appartiennent au contour
  201. segini,ipt5
  202. ideb5 = 0
  203. imil5 = nbelem2
  204. ifin5 = nseg+1
  205. segini sfis
  206.  
  207. *--------boucle sur les elements de la fissure ----------------------
  208. do iel=1,nseg
  209. c write(*,*) 'boucle sur elements de fissure',iel,'/',nseg
  210.  
  211. * io est incrémenté de +1,+2,+4 si node 1,2,3 appartient au front
  212. * ioo est incrémenté de +1,+2,+4 si node 1,2,3 appartient au contour
  213. io=0
  214. ioo=0
  215. * ---boucle sur les noeuds des elements de la fissure ----------
  216. do ib=1,num(/1)
  217. ia=ib
  218. if(ib.eq.3) ia = 4
  219. ic = num(ib,iel)
  220. if(icpr(ic).gt.0) io=io+ia
  221. if(icpr(ic).lt.0) ioo=ioo+ia
  222. * cas particulier d un segment avec 1 point appartenant au contour et un autre
  223. if(ic.eq.i2firs.or.ic.eq.i2last) ioo=ioo+ia
  224. enddo
  225. * ---fin de boucle sur les noeuds des elements de la fissure ---
  226.  
  227. * remplissage de ipt5(i5,...) et ifis(i5,1)
  228. i5=0
  229. iperm=0
  230. * ---cas ou 1 segment appartient au front
  231. if (io.eq.3.or.io.gt.4) then
  232. ideb5=ideb5+1
  233. i5=ideb5
  234. * (eventuelle) permutation des noeuds pour avoir segment 1-2 = front
  235. if (io.eq.3) then
  236. ipt5.num(1,i5)=num(1,iel)
  237. ipt5.num(2,i5)=num(2,iel)
  238. ipt5.num(3,i5)=num(3,iel)
  239. elseif(io.eq.5) then
  240. ipt5.num(1,i5)=num(3,iel)
  241. ipt5.num(2,i5)=num(1,iel)
  242. ipt5.num(3,i5)=num(2,iel)
  243. iperm=1
  244. elseif(io.eq.6) then
  245. ipt5.num(1,i5)=num(2,iel)
  246. ipt5.num(2,i5)=num(3,iel)
  247. ipt5.num(3,i5)=num(1,iel)
  248. iperm=2
  249. elseif(io.eq.7) then
  250. write(ioimp,*) 'ERREUR: '
  251. $ ,'Plusieurs cotes d un meme triangle du maillage de '
  252. $ ,'la fissure converti en TRI3 appartiennent au front!'
  253. return
  254. endif
  255. * ifis( ,1)=3 si c'est une extremité
  256. if ((ipt5.num(1,i5)).eq.i2firs) then
  257. ifis(i5,1)=3
  258. ivui2f = ivui2f+1
  259. elseif((ipt5.num(2,i5)).eq.i2last) then
  260. ifis(i5,1)=3
  261. ivui2l = ivui2l+1
  262. else
  263. * ifis( ,1)=2 si c'est un tri3 du front avec des voisins
  264. ifis(i5,1)=2
  265. endif
  266. * ---cas ou 1 (ou 2) segment(s) appartient au contour
  267. elseif(ioo.eq.3.or.ioo.gt.4) then
  268. imil5=imil5+1
  269. i5=imil5
  270. * pas de permutation des noeuds
  271. ipt5.num(1,i5)=num(1,iel)
  272. ipt5.num(2,i5)=num(2,iel)
  273. ipt5.num(3,i5)=num(3,iel)
  274. * ifis( ,1)=1 si un segment appartient au contour
  275. ifis(i5,1)=1
  276. * ---cas ou seulement 1 point appartient au front
  277. elseif(io.eq.1.or.io.eq.2.or.io.eq.4) then
  278. ifin5=ifin5-1
  279. i5=ifin5
  280. * pas de permutation des noeuds
  281. ipt5.num(1,i5)=num(1,iel)
  282. ipt5.num(2,i5)=num(2,iel)
  283. ipt5.num(3,i5)=num(3,iel)
  284. * ifis( ,1)=-1 si 1 seul point du tri3 appartient au front
  285. ifis(i5,1)=-1
  286. * ---autres cas
  287. else
  288. ifin5=ifin5-1
  289. i5=ifin5
  290. * pas de permutation des noeuds
  291. ipt5.num(1,i5)=num(1,iel)
  292. ipt5.num(2,i5)=num(2,iel)
  293. ipt5.num(3,i5)=num(3,iel)
  294. endif
  295.  
  296. * remplissage de ifis(i5,2)
  297. * ---uniquement si 1 (ou 2) point(s) appartient au contour
  298. if(ioo.eq.3.or.ioo.gt.4) then
  299. * attention a l eventuelle permutation effectuee
  300. if(iperm.ne.0) ioo=ooperm(ioo,iperm)
  301. * -ioo=3 => noeuds 1-2 = segment 1 appartient au contour
  302. if (ioo.eq.3) then
  303. ifis(i5,2)=1
  304. * -ioo=5 => noeuds 3-1 = segment 3 appartient au contour
  305. elseif(ioo.eq.5) then
  306. ifis(i5,2)=4
  307. * -ioo=6 => noeuds 2-3 = segment 2 appartient au contour
  308. elseif(ioo.eq.6) then
  309. ifis(i5,2)=2
  310. * -cas ou 3 points appartiennent au contour => travail en +
  311. * pour determiner les 2 segments dans le contour
  312. elseif(ioo.eq.7) then
  313. c write(6,*)'tri3 avec 2 segments dans le contour!',iel,i5
  314. * rem: si les 3 points appartiennent au contour, on n'a pas de segments
  315. * appartenant au front
  316. ii1 = ipt5.num(1,i5)
  317. ii2 = ipt5.num(2,i5)
  318. ii3 = ipt5.num(3,i5)
  319. c il faut trouver quels sont les 2 segments appartenant au contour !
  320. c ia1 = -icpr(ii1)
  321. c ia2 = -icpr(ii2)
  322. c ia3 = -icpr(ii3)
  323. c dans de tres rares cas, on peut avoir des points appartenant aussi au front (i2firs et i2last)
  324. ia1 = abs(icpr(ii1))
  325. ia2 = abs(icpr(ii2))
  326. ia3 = abs(icpr(ii3))
  327. ia1min = min(ia1,ia2)
  328. ia1min = min(ia1min,ia3)
  329. ia1seg = 0
  330. do 11 ia1=ia1min,nbelCON
  331. ifis1 = 0
  332. do ib=1,nbnoCON
  333. if(ii1.eq.IPT6.num(ib,ia1)) ifis1=ifis1+1
  334. if(ii2.eq.IPT6.num(ib,ia1)) ifis1=ifis1+2
  335. if(ii3.eq.IPT6.num(ib,ia1)) ifis1=ifis1+4
  336. enddo
  337. c si on n'a trouve que 1 noeud dans le ia1 ieme SEG2 du contour, on continue a chercher
  338. if(ifis1.le.2) goto 11
  339. if(ifis1.eq.4) goto 11
  340. c on a trouve les 2 noeuds dans un SEG2 du contour
  341. if(ifis1.eq.3) ifis(i5,2)=ifis(i5,2)+1
  342. if(ifis1.eq.6) ifis(i5,2)=ifis(i5,2)+2
  343. if(ifis1.eq.5) ifis(i5,2)=ifis(i5,2)+4
  344. if(ifis1.ge.7) then
  345. write(ioimp,*) 'tri3 avec 3 segments dans le contour !'
  346. call erreur(26)
  347. endif
  348. if(ifis1.eq.3.or.ifis1.ge.5) ia1seg=ia1seg+1
  349. if(ia1seg.eq.2) goto 12
  350. 11 continue
  351. c on n a pas trouvé les 2 segments appartenant au contour
  352. write(ioimp,*)'on n a pas trouvé les 2 segments',
  353. $ ' du tri3 devant appartenant au contour !'
  354. call erreur(26)
  355. return
  356. c on a trouvé les 2 segments appartenant au contour
  357. 12 continue
  358. if(ifis(i5,2).lt.3.or.ifis(i5,2).ge.7) then
  359. write(ioimp,*) 'impossible de detecter les 2 segments ',
  360. $ 'appartenant au contour avec les noeuds ',ii1,ii2,ii3
  361. write(ioimp,*) 'ifis(ideb5,2)=',ifis(ideb5,2)
  362. write(ioimp,*) 'IPT6 pour les segments :'
  363. write(ioimp,*) ' ia1',ia1,(IPT6.num(iou,ia1),iou=1,2)
  364. write(ioimp,*) ' ia2',ia2,(IPT6.num(iou,ia2),iou=1,2)
  365. write(ioimp,*) ' ia3',ia3,(IPT6.num(iou,ia3),iou=1,2)
  366. call erreur(26)
  367. return
  368. endif
  369. endif
  370. c -fin du cas elseif(ioo.eq.7) then
  371. endif
  372. c --fin du cas if(ioo.eq.3.or.ioo.gt.4) then
  373.  
  374.  
  375.  
  376. enddo
  377. * fin de boucle sur les elements de la fissure -----------------
  378.  
  379. c if ((ifin5-ideb5).ne.1) then
  380. if (ideb5.ne.nbelem2.or.(ifin5-imil5).ne.1) then
  381. write(ioimp,*) 'pb avec le rearrangement de elements de fissure'
  382. write(ioimp,*) ideb5,nbelem2,imil5,ifin5
  383. call erreur(26)
  384. return
  385. endif
  386. * si le front n est pas ferme, il faut avoir detecte ses 2 extremites
  387. * sinon il est probable que la ligne du front melfro ne soit pas
  388. * orientée comme le bord du maillage melfis
  389. if ((ivui2f.ne.1).and.(ivui2l.ne.1)) then
  390. write(ioimp,*) 'ERREUR: Les extremites du front (noeuds',
  391. $ i2firs,i2last,') ont été detectees ',ivui2f, 'et ',ivui2l
  392. $ ,' fois au lieu de 1 fois attendue'
  393. write(ioimp,*) 'Verifiez la coherence avec l orientation'
  394. $ ,'du maillage de la fissure fourni'
  395. $ ,' et INVErser le sens du front le cas échéant.'
  396. c write(6,*) 'front de fissure='
  397. c CALL ECMAIL(ipt2,JENT2)
  398. c write(6,*) 'surface de fissure='
  399. c CALL ECMAIL(ipt5,JENT5)
  400. call erreur(26)
  401. return
  402. endif
  403. c segdes,meleme,IPT6
  404. segdes,meleme,IPT6,ipt2
  405. meleme=ipt5
  406. melfis=meleme
  407. segsup icpr
  408.  
  409. endif
  410. *------ fin du Cas ou l'on calcule phi ET psi --------------------------
  411.  
  412. * sont actifs: melmai(=struct changée en poi1), ipt5=meleme=melfis,
  413. * sont inactifs: ipt2=melfro
  414.  
  415.  
  416.  
  417. ***********************************************************************
  418. * debut du travail suivant valeur de icle
  419.  
  420. if (xcrit.eq.0.d0) then
  421. icle=0
  422. xcrit=1.D10
  423. else
  424. icle=1
  425. * ajout d'un epsilon 1.D-5 de rattrapage
  426. xcrit = 1.00001D0 * xcrit
  427. endif
  428.  
  429. melfis=meleme
  430. melpoi=melmai
  431.  
  432. ***********************************************************************
  433. * IL NE RESTE PLUS QU'A FAIRE LE TRAVAIL PROPREMENT DIT
  434. ***********************************************************************
  435.  
  436. * pour chaque segment de la fissure (melfis 2eme maillage lu)
  437. * on calcule le champ phi et psi des noeuds autour de l'element
  438. if(iimpi.ge.1) then
  439. c call gibtem(xkt)
  440. write(ioimp,*)'debut travail'
  441. endif
  442. *
  443. meleme=melfis
  444. c segact,meleme
  445. ipt3=melpoi
  446. c segact,ipt3
  447. npt=ipt3.num(/2)
  448. * attention on initialise isup avec les valeurs inverses,il faudra
  449. * le retourner
  450. c c segini isup,isupfi,xdis,xdmin,xdmin1,rrlim
  451. c segini isup,isupfi,xdis,xdmin,xdmin1
  452. segini isup
  453. c xdis,xdmin,xdmin1
  454. c if(ideux.ge.1) segini xdmin2,xdmin3
  455.  
  456. NTRI3=nseg
  457. SEGINI,MTRI3
  458. call cpu_time(start_time)
  459. *=====================================================================
  460. * Boucle sur les elements de la fissure =========================
  461.  
  462. do 444 iseg=1,nseg
  463.  
  464. c coordonnees du tri3
  465. ia=num(1,iseg)
  466. x11= xcoor((ia-1)*idim1 +1)
  467. y11= xcoor((ia-1)*idim1 +2)
  468. z11= xcoor((ia-1)*idim1 +3)
  469. ib=num(2,iseg)
  470. x22= xcoor((ib-1)*idim1 +1)
  471. y22= xcoor((ib-1)*idim1 +2)
  472. z22= xcoor((ib-1)*idim1 +3)
  473. ic=num(3,iseg)
  474. x33= xcoor((ic-1)*idim1 +1)
  475. y33= xcoor((ic-1)*idim1 +2)
  476. z33= xcoor((ic-1)*idim1 +3)
  477. c if (x11.ge.x22) then
  478. c xxmi=min(x22,x33) - xcrit
  479. c xxma=max(x11,x33) + xcrit
  480. c else
  481. c xxmi=min(x11,x33) - xcrit
  482. c xxma=max(x22,x33) + xcrit
  483. c endif
  484. c if (y11.ge.y22) then
  485. c yymi=min(y22,y33) - xcrit
  486. c yyma=max(y11,y33) + xcrit
  487. c else
  488. c yymi=min(y11,y33) - xcrit
  489. c yyma=max(y22,y33) + xcrit
  490. c endif
  491. c if (z11.ge.z22) then
  492. c zzmi=min(z22,z33) - xcrit
  493. c zzma=max(z11,z33) + xcrit
  494. c else
  495. c zzmi=min(z11,z33) - xcrit
  496. c zzma=max(z22,z33) + xcrit
  497. c endif
  498.  
  499. c vecteurs des aretes
  500. v12x=x22-x11
  501. v12y=y22-y11
  502. v12z=z22-z11
  503. v23x=x33-x22
  504. v23y=y33-y22
  505. v23z=z33-z22
  506. v31x=x11-x33
  507. v31y=y11-y33
  508. v31z=z11-z33
  509.  
  510. c vecteur normal
  511. vx=v12y*(z33-z11)-v12z*(y33-y11)
  512. vy=v12z*(x33-x11)-v12x*(z33-z11)
  513. vz=v12x*(y33-y11)-v12y*(x33-x11)
  514. bb=vx*vx + vy*vy+ vz*vz
  515. aa= sqrt(bb)
  516. if (aa.le.XPETIT) then
  517. write(ioimp,*) 'element fissure de longueur nulle ',aa
  518. write(ioimp,*) 'tri3.',iseg,'.',ia,'=',x11,y11,z11
  519. write(ioimp,*) 'tri3.',iseg,'.',ib,'=',x22,y22,z22
  520. write(ioimp,*) 'tri3.',iseg,'.',ic,'=',x33,y33,z33
  521. call erreur(26)
  522. return
  523. endif
  524. sur = bb/2.d0
  525. vx=vx/aa
  526. vy=vy/aa
  527. vz=vz/aa
  528. vxyz = max(abs(vx),abs(vy))
  529. vxyz = max(vxyz,abs(vz))
  530. if (vxyz.le.XPETIT) then
  531. write(ioimp,*)'element fissure de longueur nulle',aa,vx,vy,vz
  532. call erreur(26)
  533. return
  534. endif
  535. c xsur = 1.D-8*sur
  536.  
  537. c vecteurs des aretes (suite)
  538. v12lo= sqrt(v12x*v12x+v12y*v12y+v12z*v12z)
  539. v23lo= sqrt(v23x*v23x+v23y*v23y+v23z*v23z)
  540. v31lo= sqrt(v31x*v31x+v31y*v31y+v31z*v31z)
  541. v12x= v12x/v12lo
  542. v12y= v12y/v12lo
  543. v12z= v12z/v12lo
  544. v23x= v23x/v23lo
  545. v23y= v23y/v23lo
  546. v23z= v23z/v23lo
  547. v31x= v31x/v31lo
  548. v31y= v31y/v31lo
  549. v31z= v31z/v31lo
  550.  
  551. c REMPLISSAGE DU SEGMENT MTRI3
  552. MTRI3.KNODE(1,iseg)=ia
  553. MTRI3.KNODE(2,iseg)=ib
  554. MTRI3.KNODE(3,iseg)=ic
  555. MTRI3.XNODE(1,iseg)=x11
  556. MTRI3.XNODE(2,iseg)=x22
  557. MTRI3.XNODE(3,iseg)=x33
  558. MTRI3.YNODE(1,iseg)=y11
  559. MTRI3.YNODE(2,iseg)=y22
  560. MTRI3.YNODE(3,iseg)=y33
  561. MTRI3.ZNODE(1,iseg)=z11
  562. MTRI3.ZNODE(2,iseg)=z22
  563. MTRI3.ZNODE(3,iseg)=z33
  564. MTRI3.VARE12(1,iseg)=v12x
  565. MTRI3.VARE12(2,iseg)=v12y
  566. MTRI3.VARE12(3,iseg)=v12z
  567. MTRI3.VARE12(4,iseg)=v12lo
  568. MTRI3.VARE23(1,iseg)=v23x
  569. MTRI3.VARE23(2,iseg)=v23y
  570. MTRI3.VARE23(3,iseg)=v23z
  571. MTRI3.VARE31(1,iseg)=v31x
  572. MTRI3.VARE31(2,iseg)=v31y
  573. MTRI3.VARE31(3,iseg)=v31z
  574. MTRI3.VNORM(1,iseg)=vx
  575. MTRI3.VNORM(2,iseg)=vy
  576. MTRI3.VNORM(3,iseg)=vz
  577. MTRI3.SURFAC(iseg)=sur
  578.  
  579. 444 continue
  580. *=====================================================================
  581.  
  582.  
  583. *=====================================================================
  584. * 2EME PHASE DU TRAVAIL (DOUBLE BOUCLE)
  585. *=====================================================================
  586.  
  587. C FAUT-IL PASSER EN // ?
  588. c estimation grossiere
  589. c car on n'a trouve que des cas ou la // est benefique
  590. if (nseg*npt.lt.100) then
  591. NBTHR = 1
  592. BTHRD=.false.
  593. else
  594. NBTHR = NBTHRS
  595. BTHRD = .TRUE.
  596. CALL THREADII
  597. endif
  598.  
  599. C CREATION ET REMPLISSAGE DU SEGMENT POUR LA //iSATION
  600. SEGINI,SPARAL
  601. SPARAL.NBTHRD = NBTHR
  602. SPARAL.KPT3 = IPT3
  603. SPARAL.KTRI3 = MTRI3
  604. SPARAL.KCLE = icle
  605. SPARAL.KDEUX = ideux
  606. SPARAL.XCRIT1 = xcrit
  607. SPARAL.KSUP = isup
  608. SPARAL.KPOVA1 = mpova1
  609. SPARAL.KPOVA2 = mpova2
  610. SPARAL.KPOVA3 = mpova3
  611. SPARAL.KFIS = sfis
  612. SPARAL.KMIL5 = imil5
  613. SPARAL.Kxicpr = mxicpr
  614.  
  615. IF (BTHRD) THEN
  616. C Remplissage du 'COMMON/psiabc'
  617. IPARAL=SPARAL
  618. DO ith=2,NBTHR
  619. CALL THREADID(ith,PSI3Di)
  620. ENDDO
  621. CALL PSI3Di(1)
  622.  
  623. C Attente de la fin de tous les threads en cours de travail
  624. DO ith=2,NBTHR
  625. CALL THREADIF(ith)
  626. ENDDO
  627.  
  628. C On libere les Threads
  629. CALL THREADIS
  630.  
  631. C Verification de l'erreur (Apres liberation des THREADS)
  632. DO ith=1,NBTHR
  633. IRETOU=SPARAL.IERROR(ith)
  634. IF (IRETOU .GT. 0) THEN
  635. CALL ERREUR(IRETOU)
  636. RETURN
  637. ENDIF
  638. ENDDO
  639.  
  640. ELSE
  641. C Appel de la SUBROUTINE qui fait le travail
  642. CALL PSI3D0(1,SPARAL)
  643. C Verification de l'erreur
  644. IRETOU=SPARAL.IERROR(1)
  645. IF (IRETOU .GT. 0) THEN
  646. CALL ERREUR(IRETOU)
  647. RETURN
  648. ENDIF
  649.  
  650. ENDIF
  651.  
  652. ***********************************************************************
  653. * FIN DU TRAVAIL
  654. ***********************************************************************
  655.  
  656. * on inverse isup
  657. * AVANT : isup(i) =0 veut dire que le point na pas été traité
  658. * APRES : isup(i) =1 veut dire que le point peut etre enleve (=0 sinon)
  659. iaju=0
  660. do ia=1,isup(/1)
  661. if (isup(ia).eq.0) then
  662. iaju=iaju+1
  663. isup(ia)=1
  664. else
  665. isup(ia)=0
  666. endif
  667. enddo
  668. c if (iaju+igard.ne.isup(/1)) then
  669. c write(ioimp,*) 'psiphi : tout va mal '
  670. c $ ,iaju,'+',igard,'ne',isup(/1)
  671. c call erreur(26)
  672. c return
  673. c endif
  674. c bp : test ci dessus un peu chiant a garder en // et pas tre utile
  675. ipt1=ipt3
  676. ipt4=melfis
  677. segdes ipt4
  678. * iaju est le nombre de points a éliminer
  679. * isup(i) =1 veut dire que le point peut etre enleve (=0 sinon)
  680.  
  681. * faut-il ajuster les champs?
  682. if (iaju.ne.0) then
  683. * si oui on recree les bons chpoints
  684. nbelem=ipt1.num(/2)- iaju
  685. nbnn=1
  686. nbref=0
  687. nbsous=0
  688. segini ipt2
  689. ipt2.itypel=1
  690. nc=1
  691. n=nbelem
  692. segini mpova5
  693. *>>>bp&yt: ajout ligne suivante
  694. if(ideux.ge.1) segini mpova4
  695. if(ideux.ge.2) segini mpova6
  696. iel=0
  697. do i=1,ipt1.num(/2)
  698. if (isup(i).eq.0) then
  699. iel=iel+1
  700. ipt2.num(1,iel)=ipt1.num(1,i)
  701. mpova5.vpocha(iel,1)=mpova2.vpocha(i,1)
  702. *>>>bp&yt: ajout ligne suivante
  703. * on met le meme maillage pour les 2 level set (car toujours utilisées ensemble)
  704. if(ideux.ge.1) mpova4.vpocha(iel,1)=mpova1.vpocha(i,1)
  705. if(ideux.ge.2) mpova6.vpocha(iel,1)=mpova3.vpocha(i,1)
  706. endif
  707. enddo
  708. if (ideux.ge.2) then
  709. c segact msoup3*mod
  710. msoup3.ipoval=mpova6
  711. msoup3.igeoc=ipt2
  712. segdes msoup3,mpova6
  713. segsup,mpova3
  714. endif
  715. if (ideux.ge.1) then
  716. c segact msoup1*mod
  717. msoup1.ipoval=mpova4
  718. msoup1.igeoc=ipt2
  719. segdes msoup1,mpova4
  720. segsup,mpova1
  721. endif
  722. c segact msoup2*mod
  723. msoup2.ipoval=mpova5
  724. msoup2.igeoc=ipt2
  725. segdes msoup2,mpova5,ipt2
  726. segsup,mpova2
  727. segdes ipt1
  728. else
  729. if(ideux.ge.2) segdes,msoup3,mpova3
  730. if(ideux.ge.1) segdes,msoup1,mpova1
  731. segdes msoup2,mpova2,ipt1
  732. endif
  733. if (ideux.ge.1) then
  734. segsup,icpr,sfis
  735. if(ideux.ge.2) segsup,mxicpr
  736. endif
  737. if (icle.eq.1) then
  738. c segsup isup,isupfi,xdis,xdmin,xdmin1,rrlim
  739. c segsup isup,isupfi,xdis,xdmin,xdmin1
  740. segsup isup
  741. c if(ideux.ge.1) segsup xdmin2,xdmin3
  742. endif
  743. SEGSUP,MTRI3
  744.  
  745. RETURN
  746. END
  747.  
  748.  
  749.  
  750.  
  751.  
  752.  
  753.  
  754.  
  755.  
  756.  
  757.  

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