Télécharger pro2.eso

Retour à la liste

Numérotation des lignes :

pro2
  1. C PRO2 SOURCE CB215821 22/07/20 15:39:44 11411
  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*(LOCOMP) 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.  
  152. * Si MCHAML vide : erreur
  153. if (NBSZ.eq.0) then
  154. moterr(1:8) = 'MCHAML '
  155. call erreur(1027)
  156. return
  157. endif
  158.  
  159. * write(6,*) ' nbsz',nbsz
  160. C
  161. nbnnma=0
  162. DO 10 IOB=1, NBSZ
  163. C
  164. IPT1=MCHEL1.IMACHE(IOB)
  165. SEGACT IPT1
  166. C
  167. ITY = IPT1.ITYPEL
  168. C En dimension 1 : EF massifs s'appuient sur SEG2 ou SEG3
  169. IF (IDIM.EQ.1) THEN
  170. IF (ITY.NE.2.AND.ITY.NE.3) THEN
  171. CALL ERREUR(16)
  172. RETURN
  173. ENDIF
  174. ELSE
  175. C En dimension 2,3 :
  176. cbp IF (ITY.LT.4) THEN
  177. IF (ITY.LT.4.OR.(IDIM.eq.3.and.ITY.LT.14)) THEN
  178. write(ioimp,*) 'IDIM,ITYPEL=',IDIM,ITY
  179. CALL ERREUR(16)
  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. inonac=NBOUB-IUOO
  278. INTERR(1)=inonac
  279. * write(6,*) ' nombre de points non accrochés ',interr(1)
  280. CALL ERREUR(-322)
  281. CALL SOUCIS(inonac)
  282. endif
  283.  
  284. segact ipt6
  285. Ckich : on fabrique 1 chpoint par phase presente dans le MMODEL
  286. c stocke dans un listchpo, dans l ordre de MOTS (posibilite terme nul)
  287. if (ilphmo.lt.0) then
  288. jgn = 8
  289. jgm = 1
  290. segini mlmots
  291. mots(1) = liphas(1)
  292. else
  293. mlmots = ilphmo
  294. c n1 = nphas
  295. n1 = mots(/2)
  296. segini mlchpo
  297. endif
  298.  
  299. ipout = 0
  300. DO 700 IMP = 1,MOTS(/2)
  301. drphas = .false.
  302. do iphas = 1,NPHAS
  303. if(mots(imp).eq.liphas(iphas)) then
  304. drphas = .true.
  305. goto 703
  306. endif
  307. enddo
  308. 703 if(.not.drphas) then
  309. write(6,*) ' pas de projection phase ', mots(imp)
  310. goto 700
  311. endif
  312.  
  313. C** recherche du noms de composantes existantes dans le chpoint
  314. segini nomiii
  315. NA=0
  316. DO 220 i=1,nbsz
  317. ckich maintenant maillage disjoints if(ipca(i).eq.0) go to 220
  318. if (mchel1.conche(i)(17:24).ne.liphas(iphas)) goto 220
  319. mchaml = mchel1.ichaml(i)
  320. segact mchaml
  321. do 221 j=1,nomche(/2)
  322. IF( TYPCHE(J).NE.'REAL*8') GO TO 221
  323. IF(NA.NE.0) THEN
  324. do 222 k=1,nominc(/2)
  325. if( nominc(k).eq.nomche(j)) go to 221
  326. 222 continue
  327. nominc(**)=nomche(j)
  328. ELSE
  329. nominc(**)=nomche(j)
  330. na=1
  331. ENDIF
  332. 221 continue
  333. 220 continue
  334. nnin=nominc(/2)
  335. nnnoe=iaccro
  336. * write(6,*) ' nnin nnoe ',nnin,nnnoe,iphas
  337. segini mtrav
  338. C
  339. C**** remplissage de mtrav en calculant les valeurs
  340. C
  341. do 311 iu = 1 , nominc(/2)
  342. inco(iu)=nominc(iu)
  343. 311 continue
  344. segsup nomiii
  345.  
  346. c recense les zones du mchel1 compatibles avec le maillage de ipt6 et la phase
  347. nliza6 = ipt6.lisous(/1)
  348. segini lliza
  349. do iza6 = 1,nliza6
  350. do ips = 1,nbsz
  351. IF (ipt6.lisous(iza6).eq.mchel1.imache(ips).AND.
  352. & liphas(iphas).eq.mchel1.conche(ips)(17:24)) THEN
  353. if (lliza(iza6).eq.0) then
  354. jg = 1
  355. segini mlenti
  356. LLIZA(IZA6) = mlenti
  357. lect(1) = ips
  358. else
  359. mlenti = lliza(iza6)
  360. jg = lect(/1)+1
  361. segadj mlenti
  362. lect(jg) = ips
  363. endif
  364. ENDIF
  365. enddo
  366. mlenti = lliza(iza6)
  367. enddo
  368.  
  369. segact ipt2
  370. ipla=0
  371. nbnn = nbnnma
  372. SEGINI ,WORK1
  373. do 224 I=1,nodes
  374. if(idejvu(i).eq.1) then
  375. ipla=ipla+1
  376. igeo(ipla)=ipt2.num(1,i)
  377. ip=ipt2.num(1,i)
  378. igeo(ipla)=ipt2.num(1,i)
  379. iza6=ieint(1,i)
  380. iel=ieint(2,i)
  381. *
  382. IF (LLIZA(iza6).EQ.0) goto 302
  383. mlenti = lliza(iza6)
  384. nza = lect(/1)
  385.  
  386. do 300 iz = 1,nza
  387. IZA=lect(iz)
  388.  
  389. mchaml=mchel1.ichaml(iza)
  390. ipt1=mchel1.imache(iza)
  391. segact ipt1, mchaml
  392.  
  393.  
  394. if(ipmod.ne.0) then
  395. C recherche element fini pour fonction de forme dans qsijs
  396. do lk=1,nsous
  397. if(ipt1.eq.iiele(1,lk)) then
  398. kel=iiele(2,lk)
  399. goto 31
  400. endif
  401. enddo
  402. 31 continue
  403. else
  404. C par defaut on utilise les fonctions de transformation geometriques
  405. kel =IPT1.ITYPEL
  406. endif
  407. C
  408. C Coordonnées des sommets de l'élements IEL
  409. C
  410. nbnn=ipt1.num(/1)
  411. CALL DOXE(XCOOR,IDIM,NBNN,IPT1.NUM,IEL,CORDEL)
  412. C
  413. IDI2 = IDIM
  414. do ih=1,idim
  415. qsi(ih)=qqq(ih,i)
  416. enddo
  417. C
  418. C Coordonnées du point dans l'element de référence
  419. C
  420. CALL SHAPE(QSI(1),QSI(2),QSI(3),KEL,SHPP,IRET)
  421. * if ( qsi(1).lt.-1.00001.or.qsi(1).gt.1.00001) then
  422. * if ( qsi(2).lt.-1.00001.or.qsi(2).gt.1.00001) then
  423. * if ( qsi(3).lt.-1.00001.or.qsi(3).gt.1.00001) then
  424. * write(6,*) ip,iza,iel,kel,(qsi(ih),ih=1,idim)
  425. * write(6,*) ( shpp(1,ioup),ioup=1,nbnn)
  426. * endif
  427. * endif
  428. * endif
  429. C
  430. C Calcul des valeurs pour chaque type de composantes
  431. C
  432. do 200 il=1,ielval(/1)
  433. IF (TYPCHE(il).EQ.'REAL*8') THEN
  434. C
  435. do 225 kcol=1,nnin
  436. * write(6,2323)il,nomche(il)(1:4),inco(1)
  437. * 2323 format('IL ', I2, 'nomche', A4, ' inco' , A4)
  438. if(nomche(il).eq.inco(kcol)) go to 226
  439. 225 continue
  440. 226 continue
  441. * if(kcol.ne.1) then
  442. * write(6,*) 'kcol nnin',kcol,nnin
  443. * write(6,*)' iza il ipla',iza,il,ipla
  444. * return
  445. * endif
  446. MELVAL = IELVAL(il)
  447. SEGACT MELVAL
  448. IEMN = MIN(IEL,VELCHE(/2))
  449. NHAR(kcol) = MCHEL1.INFCHE(iza,3)
  450. DO 150, NP = 1, NBNN
  451. IBMN = MIN(NP,velche(/1))
  452. bb(kcol,ipla)=SHPP(1,NP)*VELCHE(IBMN,IEMN)
  453. $ +bb(kcol,ipla)
  454. C
  455. 150 CONTINUE
  456. ibin(kcol,ipla)=ibin(kcol,ipla)+1
  457. C
  458. ENDIF
  459. 200 CONTINUE
  460. 300 continue
  461. 302 continue
  462. endif
  463. 224 CONTINUE
  464.  
  465. do iza6 = 1,nliza6
  466. mlenti = lliza(iza6)
  467. segsup mlenti
  468. enddo
  469. segsup lliza
  470.  
  471. ITRAV=MTRAV
  472. DO 4370 IB=1,NNNOE
  473. DO 43701 IA=1,NNIN
  474. * write(6,fmt='('' ia ib bb ibin '',2i4,e12.5,i5)')ia,ib,bb(ia,ib),
  475. * * ibin(ia,ib)
  476. IF (IBIN(IA,IB).GT.1) THEN
  477. * write(6,*) ' ibin gt 1 ib ,ia' , ib ,ia, ibin(ia,ib)
  478. BB(IA,IB)=BB(IA,IB) / IBIN(IA,IB)
  479. IBIN(IA,IB)=1
  480. ENDIF
  481. 43701 CONTINUE
  482. 4370 CONTINUE
  483. C
  484. C RECONSTUCTION DE LA PARTITION
  485. C
  486. CALL CRECHP(ITRAV,ICHPOI1)
  487. C
  488. IPOUT=ICHPOI1
  489. MCHPOI=ICHPOI1
  490. * le champ est forcement diffus
  491. JATTRI(1)=1
  492. MTYPOI='PROJECTI'
  493. IFOPOI=IFOUR
  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. END
  510.  
  511.  
  512.  
  513.  
  514.  
  515.  
  516.  
  517.  

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