Télécharger pro2.eso

Retour à la liste

Numérotation des lignes :

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

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