Télécharger pro2.eso

Retour à la liste

Numérotation des lignes :

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

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