Télécharger pro2.eso

Retour à la liste

Numérotation des lignes :

  1. C PRO2 SOURCE PV 20/03/30 21:22:44 10567
  2. SUBROUTINE PRO2(IPGEA,IPCHEL,isort,IPOUT,ilphmo)
  3. C=====================================================================
  4. C
  5. C SOUS-PROGRAMME DE L'OPERATEUR PROIET
  6. C
  7. C ENTREES :
  8. C
  9. C IPGEA = POINTEUR DE UN OBJET MAILLAGE.
  10. C IPCHEL = POINTEUR DE UN OBJET CHAMELEM
  11. C isort 0 on veut un chpo brutal a 1 sous zone non repartitionné
  12. C en sortie pointuer su un segment ipca
  13. C 1 on veut un champ conforme partitionne
  14. c ilphmo pointeur sur un listmot phases du modele cible
  15. C SORTIES :
  16. C
  17. C IPOUT = POINTEUR SUR UN OBJET DE TYPE CHPOINT
  18. C
  19. C=====================================================================
  20. C
  21. C TRES INSPIRE DE LA FONCTION ACCRO DEVELOPPEE PAR JMB
  22. C
  23. C=====================================================================
  24.  
  25. IMPLICIT INTEGER(I-N)
  26. IMPLICIT REAL*8 (A-H,O-Z)
  27. *-
  28. -INC CCREEL
  29.  
  30. -INC PPARAM
  31. -INC CCOPTIO
  32. -INC SMCOORD
  33. -INC SMELEME
  34. -INC SMCHPOI
  35. -INC SMLCHPO
  36. -INC SMLMOTS
  37. -INC SMCHAML
  38. -INC TMTRAV
  39. -INC CCGEOME
  40. -INC SMMODEL
  41. -INC SMLENTI
  42. SEGMENT nomiii
  43. character*4 NOMINC(0)
  44. endsegment
  45. SEGMENT MDATA
  46. INTEGER IPCA(NBSZ)
  47. ENDSEGMENT
  48. SEGMENT IPEX(NODES)
  49. SEGMENT ISOM(NBSOM(IPT3.ITYPEL))
  50. SEGMENT MTR1
  51. CHARACTER*4 IPCOM(0)
  52. ENDSEGMENT
  53. SEGMENT SEPHAS
  54. CHARACTER*8 LIPHAS(NPHAS)
  55. ENDSEGMENT
  56. SEGMENT LLIZA(NLIZA6)
  57. SEGMENT/MTR4/(IPHAR(0))
  58. segment inomil(0)
  59. SEGMENT/MTRA/(NOPOIN(nbpts))
  60. *
  61. LOGICAL drphas
  62. C
  63. segment info
  64. integer infell(JG)
  65. endsegment
  66. C
  67. SEGMENT ICOUNT
  68. INTEGER IEINT(2,NODES),IDEJVU(NODES),IACCRO
  69. * ient (1,i)= j, ient(2,i)=k veut dire que le ieme point de ipt2 est dans
  70. * le k ieme element de la jeme sous zone.
  71. * idejvu(i)=1 veut dire que le ieme noeud de ipt2 est deja traité
  72. REAL*8 QQQ(3,NODES),CKRIT(NODES)
  73. * CKRIT est le meilleur des criteres déja rencontré
  74. ENDSEGMENT
  75. SEGMENT ICOSOM
  76. INTEGER IPTT(NBSZ),ICENV(NBSZ),IVOL(NBSZ)
  77. ENDSEGMENT
  78. SEGMENT WORK1
  79. REAL*8 CORDEL(3,NBNN),QSI(3),XPO(3),SHPP(6,NBNN),V1(3),V2(3)
  80. REAL*8 XEL(3,NBNN),BPSS(3,3),XPT(3),cored(3,NBNN),V3(3)
  81. REAL*8 R(3),U(3),W(3,3)
  82. INTEGER II(3)
  83. ENDSEGMENT
  84. segment siele
  85. integer iiele(2,nsous)
  86. endsegment
  87. external shape
  88.  
  89. C lecture optionnelle du modele associe au MCHAML origine
  90. IPMOD=0
  91. SIELE= 0
  92. call LIROBJ('MMODEL ',IPMOD,0,iret)
  93. * write(6,*) ' dans pro2 lecture de ipmod ',ipmod
  94. if(iret.eq.1) then
  95. call ACTOBJ('MMODEL ',IPMOD,1)
  96. mmodel = ipmod
  97. NSOUS=KMODEL(/1)
  98. segini siele
  99. DO 30 ISOUS = 1, NSOUS
  100. IMODEL = KMODEL(ISOUS)
  101. MELEME = IMAMOD
  102. MELE = NEFMOD
  103. C COQUE INTEGREE OU PAS ?
  104. IF(INFMOD(/1).NE.0)THEN
  105. NPINT=INFMOD(1)
  106. ELSE
  107. NPINT=0
  108. ENDIF
  109. C
  110. call elquoi(mele,npint,6,ipinf,imodel)
  111. if (ierr.ne.0) then
  112. C#MC segsup mchelm
  113. iret=0
  114. return
  115. endif
  116. C
  117. info=ipinf
  118. iiele(1,isous)=IMAMOD
  119. iiele(2,isous)= infelL(14)
  120. segsup info
  121. 30 CONTINUE
  122. endif
  123. C lecture optionnelle du critere de rattrapage de points exterieurs
  124. ccrit=1.D-5
  125. call lirree(ccrit,0,iretou)
  126. ccrit = -ABS(ccrit)
  127. * write(6,*) ' dans pro2 lecture de ccrit', ccrit
  128. C
  129. C ON VA CHANGER LE MAILLAGE RECEPTEUR EN ELEMENTS POI1
  130. C SI IL NE L EST PAS DEJA
  131.  
  132. IPT2 = IPGEA
  133. SEGACT IPT2
  134. IF(IPT2.ITYPEL.NE.1) CALL CHANGE(IPGEA,1)
  135. IF(IERR.NE.0) RETURN
  136. C
  137. IPT2 = IPGEA
  138. SEGACT IPT2
  139. NODES = IPT2.NUM(/2)
  140. SEGINI ICOUNT
  141. IACCRO = 0
  142. do 3259 iu=1,nodes
  143. ckrit(iu)=-1.d10
  144. 3259 continue
  145. C
  146. C LE MAILLAGE SUPPORT DU CHAMELEM DOIT ETRE MASSIF
  147. C
  148. MCHEL1 = IPCHEL
  149. SEGACT MCHEL1
  150. NBSZ = MCHEL1.IMACHE(/1)
  151. * write(6,*) ' nbsz',nbsz
  152. C
  153. nbnnma=0
  154. DO 10 IOB=1, NBSZ
  155. C
  156. IPT1=MCHEL1.IMACHE(IOB)
  157. SEGACT IPT1
  158. C
  159. ITY = IPT1.ITYPEL
  160. C En dimension 1 : EF massifs s'appuient sur SEG2 ou SEG3
  161. IF (IDIM.EQ.1) THEN
  162. IF (ITY.NE.2.AND.ITY.NE.3) THEN
  163. CALL ERREUR(16)
  164. RETURN
  165. ENDIF
  166. ELSE
  167. C En dimension 2,3 :
  168. cbp IF (ITY.LT.4) THEN
  169. IF (ITY.LT.4.OR.(IDIM.eq.3.and.ITY.LT.14)) THEN
  170. write(ioimp,*) 'IDIM,ITYPEL=',IDIM,ITY
  171. CALL ERREUR(16)
  172. RETURN
  173. ENDIF
  174. ENDIF
  175. nbnnma=max(nbnnma,ipt1.num(/1))
  176. C
  177. 10 CONTINUE
  178. C SEGACT MCOORD
  179.  
  180. SEGINI MDATA
  181. Ckich on regarde les phases du MCHAML a projeter
  182. nphas = nbsz
  183. segini sephas
  184. liphas(1) = mchel1.conche(1)(17:24)
  185. nphas = 1
  186. if ( mchel1.conche(/2).gt.1) then
  187. do 11 jpha =2, mchel1.conche(/2)
  188. do k = 1,nphas
  189. if (mchel1.conche(jpha)(17:24).eq.liphas(k)) goto 11
  190. enddo
  191. nphas = nphas + 1
  192. liphas(nphas) = mchel1.conche(jpha)(17:24)
  193. 11 continue
  194. segadj sephas
  195. endif
  196. Ckich on fabrique un maillage contenant tous les maillages du champ
  197. c de pointeurs distincts
  198. nbsous=nbsz
  199. nbelem=0
  200. nbnn=0
  201. nbref=0
  202. segini ipt6
  203.  
  204. nbsous = 1
  205. ipt6.lisous(1) = mchel1.imache(1)
  206. DO 20 IOB=2,NBSZ
  207.  
  208. ipt5 = mchel1.imache(iob)
  209. segact ipt5
  210. if (ipt5.itypel.eq.48) goto 20
  211. do is= 1, nbsous
  212. if (ipt6.lisous(is).eq.mchel1.imache(iob)) goto 20
  213. enddo
  214. c verifie que les maillages sont disjoints : avertit
  215. do 22 is= 1, nbsous
  216. ip1 = ipt6.lisous(is)
  217. ip2 = mchel1.imache(iob)
  218. call interb(ip1,ip2,iret,iob1)
  219. if (iret.eq.1) goto 22
  220. write(ioimp,*) 'maillages non disjoints dans le MCHAML'
  221. * kich : on peut faire mieux, pardon.
  222. goto 23
  223. 22 continue
  224.  
  225. 23 nbsous = nbsous + 1
  226. IPT6.lisous(nbsous)=MCHEL1.IMACHE(IOB)
  227. 20 continue
  228. segadj ipt6
  229. C
  230. C ICI LE ZONAGE
  231. C============================================
  232. CALL ZONEG2(IPT6,IPT2,ICount,1)
  233. c write(6,*) ' dans pro2 sortie de zoneg2'
  234. C============================================
  235. C
  236. C-------------------------
  237. C TABLEAU des POINT PAS VUS
  238. C en previson de rectifications ulterieures
  239. NBOUB = NODES - IACCRO
  240. * write(6,*) ' nboubl ' , nboub
  241. SEGACT IPT2
  242.  
  243. * SEGINI IOUBLI
  244. * SEGACT IEXCLU
  245. IUOO = 0
  246. DO 21 I=1,NODES
  247. * IP = IPT2.NUM(1,I)
  248. IF(IDEJVU(I).EQ.1) THEN
  249. go to 210
  250. ELSE
  251. IF(CKRIT(I).GT.CCRIT) THEN
  252. IDEJVU(I)=1
  253. IUOO= IUOO + 1
  254. * IPTOU(IUO,1)=IP
  255. ELSE
  256. * IP = IPT2.NUM(1,I)
  257. * write(6,*) ' ip ckrit (i)'
  258. * write(6,*) ip, ckrit(i), qqq(1,i),qqq(2,i),qqq(3,i)
  259. GO TO 21
  260. ENDIF
  261. endif
  262. 210 IPCA(IEINT(1,I))=1
  263. 21 CONTINUE
  264. iaccro=iaccro+iuoo
  265. * write(6,*) ' nodes iaccro',nodes,iaccro
  266. CC on donne le nombre de points non projetés
  267.  
  268. if (NBOUB.ne.iuoO) then
  269. inonac=NBOUB-IUOO
  270. INTERR(1)=inonac
  271. * write(6,*) ' nombre de points non accrochés ',interr(1)
  272. CALL ERREUR(-322)
  273. CALL SOUCIS(inonac)
  274. endif
  275.  
  276. segact ipt6
  277. Ckich : on fabrique 1 chpoint par phase presente dans le MMODEL
  278. c stocke dans un listchpo, dans l ordre de MOTS (posibilite terme nul)
  279. if (ilphmo.lt.0) then
  280. jgn = 8
  281. jgm = 1
  282. segini mlmots
  283. mots(1) = liphas(1)
  284. else
  285. mlmots = ilphmo
  286. c n1 = nphas
  287. n1 = mots(/2)
  288. segini mlchpo
  289. endif
  290.  
  291. ipout = 0
  292. DO 700 IMP = 1,MOTS(/2)
  293. drphas = .false.
  294. do iphas = 1,NPHAS
  295. if(mots(imp).eq.liphas(iphas)) then
  296. drphas = .true.
  297. goto 703
  298. endif
  299. enddo
  300. 703 if(.not.drphas) then
  301. write(6,*) ' pas de projection phase ', mots(imp)
  302. goto 700
  303. endif
  304.  
  305. C** recherche du noms de composantes existantes dans le chpoint
  306. segini nomiii
  307. NA=0
  308. DO 220 i=1,nbsz
  309. ckich maintenant maillage disjoints if(ipca(i).eq.0) go to 220
  310. if (mchel1.conche(i)(17:24).ne.liphas(iphas)) goto 220
  311. mchaml = mchel1.ichaml(i)
  312. segact mchaml
  313. do 221 j=1,nomche(/2)
  314. IF( TYPCHE(J).NE.'REAL*8') GO TO 221
  315. IF(NA.NE.0) THEN
  316. do 222 k=1,nominc(/2)
  317. if( nominc(k).eq.nomche(j)(1:4)) go to 221
  318. 222 continue
  319. nominc(**)=nomche(j)
  320. ELSE
  321. nominc(**)=nomche(j)
  322. na=1
  323. ENDIF
  324. 221 continue
  325. 220 continue
  326. nnin=nominc(/2)
  327. nnnoe=iaccro
  328. * write(6,*) ' nnin nnoe ',nnin,nnnoe,iphas
  329. segini mtrav
  330. C
  331. C**** remplissage de mtrav en calculant les valeurs
  332. C
  333. do 311 iu = 1 , nominc(/2)
  334. inco(iu)=nominc(iu)
  335. 311 continue
  336. segsup nomiii
  337.  
  338. c recense les zones du mchel1 compatibles avec le maillage de ipt6 et la phase
  339. nliza6 = ipt6.lisous(/1)
  340. segini lliza
  341. do iza6 = 1,nliza6
  342. do ips = 1,nbsz
  343. IF (ipt6.lisous(iza6).eq.mchel1.imache(ips).AND.
  344. & liphas(iphas).eq.mchel1.conche(ips)(17:24)) THEN
  345. if (lliza(iza6).eq.0) then
  346. jg = 1
  347. segini mlenti
  348. LLIZA(IZA6) = mlenti
  349. lect(1) = ips
  350. else
  351. mlenti = lliza(iza6)
  352. jg = lect(/1)+1
  353. segadj mlenti
  354. lect(jg) = ips
  355. endif
  356. ENDIF
  357. enddo
  358. mlenti = lliza(iza6)
  359. enddo
  360.  
  361. segact ipt2
  362. ipla=0
  363. nbnn = nbnnma
  364. SEGINI ,WORK1
  365. do 224 I=1,nodes
  366. if(idejvu(i).eq.1) then
  367. ipla=ipla+1
  368. igeo(ipla)=ipt2.num(1,i)
  369. ip=ipt2.num(1,i)
  370. igeo(ipla)=ipt2.num(1,i)
  371. iza6=ieint(1,i)
  372. iel=ieint(2,i)
  373. *
  374. IF (LLIZA(iza6).EQ.0) goto 302
  375. mlenti = lliza(iza6)
  376. nza = lect(/1)
  377.  
  378. do 300 iz = 1,nza
  379. IZA=lect(iz)
  380.  
  381. mchaml=mchel1.ichaml(iza)
  382. ipt1=mchel1.imache(iza)
  383. segact ipt1, mchaml
  384.  
  385.  
  386. if(ipmod.ne.0) then
  387. C recherche element fini pour fonction de forme dans qsijs
  388. do lk=1,nsous
  389. if(ipt1.eq.iiele(1,lk)) then
  390. kel=iiele(2,lk)
  391. goto 31
  392. endif
  393. enddo
  394. 31 continue
  395. else
  396. C par defaut on utilise les fonctions de transformation geometriques
  397. kel =IPT1.ITYPEL
  398. endif
  399. C
  400. C Coordonnées des sommets de l'élements IEL
  401. C
  402. nbnn=ipt1.num(/1)
  403. CALL DOXE(XCOOR,IDIM,NBNN,IPT1.NUM,IEL,CORDEL)
  404. C
  405. IDI2 = IDIM
  406. do ih=1,idim
  407. qsi(ih)=qqq(ih,i)
  408. enddo
  409. C
  410. C Coordonnées du point dans l'element de référence
  411. C
  412. CALL SHAPE(QSI(1),QSI(2),QSI(3),KEL,SHPP,IRET)
  413. * if ( qsi(1).lt.-1.00001.or.qsi(1).gt.1.00001) then
  414. * if ( qsi(2).lt.-1.00001.or.qsi(2).gt.1.00001) then
  415. * if ( qsi(3).lt.-1.00001.or.qsi(3).gt.1.00001) then
  416. * write(6,*) ip,iza,iel,kel,(qsi(ih),ih=1,idim)
  417. * write(6,*) ( shpp(1,ioup),ioup=1,nbnn)
  418. * endif
  419. * endif
  420. * endif
  421. C
  422. C Calcul des valeurs pour chaque type de composantes
  423. C
  424. do 200 il=1,ielval(/1)
  425. IF (TYPCHE(il).EQ.'REAL*8') THEN
  426. C
  427. do 225 kcol=1,nnin
  428. * write(6,2323)il,nomche(il)(1:4),inco(1)
  429. * 2323 format('IL ', I2, 'nomche', A4, ' inco' , A4)
  430. if(nomche(il)(1:4).eq.inco(kcol)) go to 226
  431. 225 continue
  432. 226 continue
  433. * if(kcol.ne.1) then
  434. * write(6,*) 'kcol nnin',kcol,nnin
  435. * write(6,*)' iza il ipla',iza,il,ipla
  436. * return
  437. * endif
  438. MELVAL = IELVAL(il)
  439. SEGACT MELVAL
  440. IEMN = MIN(IEL,VELCHE(/2))
  441. NHAR(kcol) = MCHEL1.INFCHE(iza,3)
  442. DO 150, NP = 1, NBNN
  443. IBMN = MIN(NP,velche(/1))
  444. bb(kcol,ipla)=SHPP(1,NP)*VELCHE(IBMN,IEMN)
  445. $ +bb(kcol,ipla)
  446. C
  447. 150 CONTINUE
  448. ibin(kcol,ipla)=ibin(kcol,ipla)+1
  449. C
  450. ENDIF
  451. 200 CONTINUE
  452. 300 continue
  453. 302 continue
  454. endif
  455. 224 CONTINUE
  456.  
  457. do iza6 = 1,nliza6
  458. mlenti = lliza(iza6)
  459. segsup mlenti
  460. enddo
  461. segsup lliza
  462.  
  463. ITRAV=MTRAV
  464. DO 4370 IB=1,NNNOE
  465. DO 43701 IA=1,NNIN
  466. * write(6,fmt='('' ia ib bb ibin '',2i4,e12.5,i5)')ia,ib,bb(ia,ib),
  467. * * ibin(ia,ib)
  468. IF (IBIN(IA,IB).GT.1) THEN
  469. * write(6,*) ' ibin gt 1 ib ,ia' , ib ,ia, ibin(ia,ib)
  470. BB(IA,IB)=BB(IA,IB) / IBIN(IA,IB)
  471. IBIN(IA,IB)=1
  472. ENDIF
  473. 43701 CONTINUE
  474. 4370 CONTINUE
  475. C
  476. C RECONSTUCTION DE LA PARTITION
  477. C
  478. CALL CRECHP(ITRAV,ICHPOI1)
  479. C
  480. IPOUT=ICHPOI1
  481. MCHPOI=ICHPOI1
  482. * le champ est forcement diffus
  483. JATTRI(1)=1
  484. MTYPOI='PROJECTI'
  485. IFOPOI=IFOUR
  486. SEGSUP MTRAV
  487. ckich dans le listchpoint
  488. * write(6,*) 'pro2',ilphmo,mchpoi,ipout
  489. if (ilphmo.gt.0) then
  490. * write(6,*) 'pro2',mchpoi,ipout,'iphas',liphas(iphas),iphas,imp
  491. ichpoi(imp) = MCHPOI
  492. endif
  493. 700 CONTINUE
  494. if (ipout.le.0) then
  495. call erreur(-79)
  496. return
  497. endif
  498. if (ilphmo.gt.0) ipout = mlchpo
  499. SEGSUP ICOUNT,MDATA
  500.  
  501. END
  502.  
  503.  
  504.  
  505.  

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