Télécharger pro2.eso

Retour à la liste

Numérotation des lignes :

  1. C PRO2 SOURCE BP208322 16/11/18 21:20:15 9177
  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. do is= 1, nbsous
  216. if (ipt6.lisous(is).eq.mchel1.imache(iob)) goto 20
  217. enddo
  218. c verifie que les maillages sont disjoints : avertit
  219. do 22 is= 1, nbsous
  220. ip1 = ipt6.lisous(is)
  221. ip2 = mchel1.imache(iob)
  222. call interb(ip1,ip2,iret,iob1)
  223. if (iret.eq.1) goto 22
  224. write(ioimp,*) 'maillages non disjoints dans le MCHAML'
  225. * kich : on peut faire mieux, pardon.
  226. goto 23
  227. 22 continue
  228.  
  229. 23 nbsous = nbsous + 1
  230. IPT6.lisous(nbsous)=MCHEL1.IMACHE(IOB)
  231. 20 continue
  232. segadj ipt6
  233. C
  234. C ICI LE ZONAGE
  235. C============================================
  236. CALL ZONEG2(IPT6,IPT2,ICount,1)
  237. c write(6,*) ' dans pro2 sortie de zoneg2'
  238. C============================================
  239. C
  240. C-------------------------
  241. C TABLEAU des POINT PAS VUS
  242. C en previson de rectifications ulterieures
  243. NBOUB = NODES - IACCRO
  244. * write(6,*) ' nboubl ' , nboub
  245. SEGACT IPT2
  246.  
  247. * SEGINI IOUBLI
  248. * SEGACT IEXCLU
  249. IUOO = 0
  250. DO 21 I=1,NODES
  251. IP = IPT2.NUM(1,I)
  252. IF(IDEJVU(I).EQ.1) THEN
  253. go to 210
  254. ELSE
  255. IF(CKRIT(I).GT.CCRIT) THEN
  256. IDEJVU(I)=1
  257. IUOO= IUOO + 1
  258. * IPTOU(IUO,1)=IP
  259. ELSE
  260. * IP = IPT2.NUM(1,I)
  261. * write(6,*) ' ip ckrit (i)'
  262. * write(6,*) ip, ckrit(i), qqq(1,i),qqq(2,i),qqq(3,i)
  263. GO TO 21
  264. ENDIF
  265. endif
  266. 210 IPCA(IEINT(1,I))=1
  267. 21 CONTINUE
  268. iaccro=iaccro+iuoo
  269. * write(6,*) ' nodes iaccro',nodes,iaccro
  270. CC on donne le nombre de points non projetés
  271.  
  272. if (NBOUB.ne.iuoO) then
  273. INTERR(1)=NBOUB-IUOO
  274. * write(6,*) ' nombre de points non accrochés ',interr(1)
  275. CALL ERREUR(-322)
  276. endif
  277.  
  278. segact ipt6
  279. Ckich : on fabrique 1 chpoint par phase presente dans le MMODEL
  280. c stocke dans un listchpo, dans l ordre de MOTS (posibilite terme nul)
  281. if (ilphmo.lt.0) then
  282. jgn = 8
  283. jgm = 1
  284. segini mlmots
  285. mots(1) = liphas(1)
  286. else
  287. mlmots = ilphmo
  288. c n1 = nphas
  289. n1 = mots(/2)
  290. segini mlchpo
  291. endif
  292.  
  293. ipout = 0
  294. DO 700 IMP = 1,MOTS(/2)
  295. drphas = .false.
  296. do iphas = 1,NPHAS
  297. if(mots(imp).eq.liphas(iphas)) then
  298. drphas = .true.
  299. goto 703
  300. endif
  301. enddo
  302. 703 if(.not.drphas) then
  303. write(6,*) ' pas de projection phase ', mots(imp)
  304. goto 700
  305. endif
  306.  
  307. C** recherche du noms de composantes existantes dans le chpoint
  308. segini nomiii
  309. NA=0
  310. DO 220 i=1,nbsz
  311. ckich maintenant maillage disjoints if(ipca(i).eq.0) go to 220
  312. if (mchel1.conche(i)(17:24).ne.liphas(iphas)) goto 220
  313. mchaml = mchel1.ichaml(i)
  314. segact mchaml
  315. do 221 j=1,nomche(/2)
  316. IF( TYPCHE(J).NE.'REAL*8') GO TO 221
  317. IF(NA.NE.0) THEN
  318. do 222 k=1,nominc(/2)
  319. if( nominc(k).eq.nomche(j)(1:4)) go to 221
  320. 222 continue
  321. nominc(**)=nomche(j)
  322. ELSE
  323. nominc(**)=nomche(j)
  324. na=1
  325. ENDIF
  326. 221 continue
  327. 220 continue
  328. nnin=nominc(/2)
  329. nnnoe=iaccro
  330. * write(6,*) ' nnin nnoe ',nnin,nnnoe,iphas
  331. segini mtrav
  332. C
  333. C**** remplissage de mtrav en calculant les valeurs
  334. C
  335. do 311 iu = 1 , nominc(/2)
  336. inco(iu)=nominc(iu)
  337. 311 continue
  338. segsup nomiii
  339.  
  340. c recense les zones du mchel1 compatibles avec le maillage de ipt6 et la phase
  341. nliza6 = ipt6.lisous(/1)
  342. segini lliza
  343. do iza6 = 1,nliza6
  344. do ips = 1,nbsz
  345. IF (ipt6.lisous(iza6).eq.mchel1.imache(ips).AND.
  346. & liphas(iphas).eq.mchel1.conche(ips)(17:24)) THEN
  347. if (lliza(iza6).eq.0) then
  348. jg = 1
  349. segini mlenti
  350. LLIZA(IZA6) = mlenti
  351. lect(1) = ips
  352. else
  353. mlenti = lliza(iza6)
  354. jg = lect(/1)+1
  355. segadj mlenti
  356. lect(jg) = ips
  357. endif
  358. ENDIF
  359. enddo
  360. mlenti = lliza(iza6)
  361. enddo
  362.  
  363. segact ipt2
  364. ipla=0
  365. nbnn = nbnnma
  366. SEGINI ,WORK1
  367. do 224 I=1,nodes
  368. if(idejvu(i).eq.1) then
  369. ipla=ipla+1
  370. igeo(ipla)=ipt2.num(1,i)
  371. ip=ipt2.num(1,i)
  372. igeo(ipla)=ipt2.num(1,i)
  373. iza6=ieint(1,i)
  374. iel=ieint(2,i)
  375. *
  376. IF (LLIZA(iza6).EQ.0) goto 302
  377. mlenti = lliza(iza6)
  378. nza = lect(/1)
  379.  
  380. do 300 iz = 1,nza
  381. IZA=lect(iz)
  382.  
  383. 301 mchaml=mchel1.ichaml(iza)
  384. ipt1=mchel1.imache(iza)
  385. segact ipt1, mchaml
  386.  
  387.  
  388. if(ipmod.ne.0) then
  389. C recherche element fini pour fonction de forme dans qsijs
  390. do lk=1,nsous
  391. if(ipt1.eq.iiele(1,lk)) then
  392. kel=iiele(2,lk)
  393. goto 31
  394. endif
  395. enddo
  396. 31 continue
  397. else
  398. C par defaut on utilise les fonctions de transformation geometriques
  399. kel =IPT1.ITYPEL
  400. endif
  401. C
  402. C Coordonnées des sommets de l'élements IEL
  403. C
  404. nbnn=ipt1.num(/1)
  405. CALL DOXE(XCOOR,IDIM,NBNN,IPT1.NUM,IEL,CORDEL)
  406. C
  407. IDI2 = IDIM
  408. do ih=1,idim
  409. qsi(ih)=qqq(ih,i)
  410. enddo
  411. C
  412. C Coordonnées du point dans l'element de référence
  413. C
  414. CALL SHAPE(QSI(1),QSI(2),QSI(3),KEL,SHPP,IRET)
  415. * if ( qsi(1).lt.-1.00001.or.qsi(1).gt.100001) then
  416. * if ( qsi(2).lt.-1.00001.or.qsi(2).gt.100001) then
  417. * if ( qsi(3).lt.-1.00001.or.qsi(3).gt.100001) then
  418. * write(6,*) ip,iza,iel,kel,(qsi(ih),ih=1,idim)
  419. * write(6,*) ( shpp(1,ioup),ioup=1,nbnn)
  420. * endif
  421. * endif
  422. * endif
  423. C
  424. C Calcul des valeurs pour chaque type de composantes
  425. C
  426. do 200 il=1,ielval(/1)
  427. IF (TYPCHE(il).EQ.'REAL*8') THEN
  428. C
  429. do 225 kcol=1,nnin
  430. * write(6,2323)il,nomche(il)(1:4),inco(1)
  431. * 2323 format('IL ', I2, 'nomche', A4, ' inco' , A4)
  432. if(nomche(il)(1:4).eq.inco(kcol)) go to 226
  433. 225 continue
  434. 226 continue
  435. * if(kcol.ne.1) then
  436. * write(6,*) 'kcol nnin',kcol,nnin
  437. * write(6,*)' iza il ipla',iza,il,ipla
  438. * return
  439. * endif
  440. MELVAL = IELVAL(il)
  441. SEGACT MELVAL
  442. IEMN = MIN(IEL,VELCHE(/2))
  443. NHAR(kcol) = MCHEL1.INFCHE(iza,3)
  444. DO 150, NP = 1, NBNN
  445. IBMN = MIN(NP,velche(/1))
  446. bb(kcol,ipla)=SHPP(1,NP)*VELCHE(IBMN,IEMN)+bb(kcol,ipla)
  447. C
  448. 150 CONTINUE
  449. ibin(kcol,ipla)=ibin(kcol,ipla)+1
  450. * SEGDES MELVAL
  451. C
  452. ENDIF
  453. 200 CONTINUE
  454. 300 continue
  455. 302 continue
  456. endif
  457. 224 CONTINUE
  458.  
  459. do iza6 = 1,nliza6
  460. mlenti = lliza(iza6)
  461. segsup mlenti
  462. enddo
  463. segsup lliza
  464.  
  465. ITRAV=MTRAV
  466. DO 4370 IB=1,NNNOE
  467. DO 4370 IA=1,NNIN
  468. * write(6,fmt='('' ia ib bb ibin '',2i4,e12.5,i5)')ia,ib,bb(ia,ib),
  469. * * ibin(ia,ib)
  470. IF (IBIN(IA,IB).GT.1) THEN
  471. * write(6,*) ' ibin gt 1 ib ,ia' , ib ,ia, ibin(ia,ib)
  472. BB(IA,IB)=BB(IA,IB) / IBIN(IA,IB)
  473. IBIN(IA,IB)=1
  474. ENDIF
  475. 4370 CONTINUE
  476. * SEGDES MTRAV
  477. C
  478. C RECONSTUCTION DE LA PARTITION
  479. C
  480. CALL CRECHP(ITRAV,ICHPOI1)
  481. C
  482. IPOUT=ICHPOI1
  483. MCHPOI=ICHPOI1
  484. SEGACT MCHPOI*MOD
  485. * le champ est forcement diffus
  486. JATTRI(1)=1
  487. MTYPOI='PROJECTI'
  488. IFOPOI=IFOUR
  489. SEGDES MCHPOI
  490. SEGSUP MTRAV
  491. ckich dans le listchpoint
  492. * write(6,*) 'pro2',ilphmo,mchpoi,ipout
  493. if (ilphmo.gt.0) then
  494. * write(6,*) 'pro2',mchpoi,ipout,'iphas',liphas(iphas),iphas,imp
  495. ichpoi(imp) = MCHPOI
  496. endif
  497. 700 CONTINUE
  498. if (ipout.le.0) then
  499. call erreur(-79)
  500. return
  501. endif
  502. if (ilphmo.gt.0) ipout = mlchpo
  503. SEGSUP ICOUNT,MDATA
  504.  
  505. RETURN
  506. END
  507.  
  508.  
  509.  
  510.  
  511.  
  512.  

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