Télécharger fuchpo.eso

Retour à la liste

Numérotation des lignes :

  1. C FUCHPO SOURCE CB215821 18/09/27 21:15:20 9936
  2. SUBROUTINE FUCHPO(IP1,IP2,IRET)
  3. C======================================================================
  4. C fonction:
  5. C sous routine pour fusionner deux champs par points diffus
  6. C
  7. C arguments:
  8. C ip1 (E) pointeur sur le premier des deux champ par point
  9. C ip2 (E) pointeur sur le second des deux champ par point
  10. C iret (S) pointeur sur le champ par point resultat
  11. C
  12. C variables:
  13. C
  14. C * mtrav : - bb(i,j) est la valeur de la ieme inconnue de champ pour
  15. C le jieme noeud du tableau igeo .
  16. C - inco(nnin) contient le nom des nnin inconnues differentes
  17. C - ibin(i,j)=1 ou 0 indique que la ieme inconnue du champ
  18. C existe pour le jieme noeud du tableau igeo .
  19. C - igeo(i) est le numero a mettre dans un objet meleme pour
  20. C referencer le ieme noeud .
  21. C
  22. C * mtra : - nopoin(i) adresse de colonne dans bb et ibin des valeurs
  23. C correspondant au noeud i .
  24. C
  25. C * mtr1 : - ipcom liste des noms des inconnues permet de creer inco .
  26. C
  27. C * mtr2 : - iicom adresse dans ipcom des inconnues correspondant au
  28. C 2ieme ch point .
  29. C
  30. C * mtr3 : - index tableau de correspondance entre les supports geome-
  31. C triques du 1er chpoint et du 2ieme chpoint .
  32. C
  33. C auteur: A de Gayffier 13/06/94
  34. C======================================================================
  35. IMPLICIT INTEGER(I-N)
  36. IMPLICIT REAL*8(A-H,O-Z)
  37. -INC SMCHPOI
  38. -INC CCOPTIO
  39. -INC SMELEME
  40. -INC SMCOORD
  41. -INC TMTRAV
  42. SEGMENT/MTRA/(NOPOIN(XCOOR(/1)/(IDIM+1)))
  43. SEGMENT MTR1
  44. CHARACTER*4 IPCOM(0)
  45. ENDSEGMENT
  46. SEGMENT/MTR2/(IICOM(0))
  47. SEGMENT/MTR3/(INDEX(0))
  48. SEGMENT/MTR4/(IPHAR(0))
  49. C ordre de grandeur des composantes
  50. SEGMENT/MTR5/(DMOY(NNIN))
  51. C
  52. DIMENSION IPO(2)
  53. CHARACTER*8 MOT
  54. DATA UN,ZERO/1.D0,0.D0/
  55. character*4 mcle(1)
  56. data mcle/'NOER'/
  57. C
  58. noer=0
  59. call lirmot(mcle,1,noer,0)
  60. segact mcoord
  61. IRET = 0
  62. MCHPO1 = IP1
  63. MCHPO2 = IP2
  64. SEGACT MCHPO1,MCHPO2
  65. C
  66. C verification de la compatibilité des natures
  67. C
  68. NAT1 = MCHPO1.JATTRI(1)
  69. NAT2 = MCHPO2.JATTRI(1)
  70. C
  71. IF ( (NAT1*NAT2) .EQ. 0) THEN
  72. C une des deux nature est indeterminée
  73. CALL ERREUR(650)
  74. SEGDES MCHPO1,MCHPO2
  75. RETURN
  76. ELSE
  77. NANA = (NAT1 + NAT2) / NAT1 / NAT2
  78. IF ((NAT1 .EQ. 2) .AND. (NAT2 .EQ. 2)) THEN
  79. C les deux champ sont discrets: on somme les composantes communes
  80. CALL ADCHPO(IP1,IP2,IRET,1D0,1D0)
  81. RETURN
  82. ELSE IF ((NAT1 .NE. 1) .OR. (NAT2 .NE. 1)) THEN
  83. C les natures ne sont pas compatibles
  84. CALL ERREUR(649)
  85. RETURN
  86. ENDIF
  87. ENDIF
  88. C
  89. C les deux champ sont de nature diffuse
  90. C on moyenne les composantes communes
  91. C
  92. IF ( IP1 .NE. IP2) GOTO 60
  93. C
  94. C *** cas ou les 2 pointeurs ipo1 et ipo2 sont egaux
  95. C
  96. MCHPO1=IP1
  97. SEGACT MCHPO1
  98. NSOUPO=MCHPO1.IPCHP(/1)
  99. NAT=MCHPO1.JATTRI(/1)
  100. SEGINI MCHPOI
  101. DO 10 I=1,NAT
  102. JATTRI(I)=MCHPO1.JATTRI(I)
  103. 10 CONTINUE
  104. IRET=MCHPOI
  105. MOCHDE=MCHPO1.MOCHDE
  106. MTYPOI=MCHPO1.MTYPOI
  107. IFOPOI=IFOUR
  108. DO 50 IA=1,NSOUPO
  109. MSOUP1=MCHPO1.IPCHP(IA)
  110. SEGACT MSOUP1
  111. NC=MSOUP1.NOCOMP(/2)
  112. SEGINI MSOUPO
  113. IPCHP(IA)=MSOUPO
  114. DO 20 IB=1,NC
  115. NOCOMP(IB)=MSOUP1.NOCOMP(IB)
  116. NOHARM(IB)=MSOUP1.NOHARM(IB)
  117. 20 CONTINUE
  118. IGEOC=MSOUP1.IGEOC
  119. MPOVA1=MSOUP1.IPOVAL
  120. SEGACT MPOVA1
  121. N =MPOVA1.VPOCHA(/1)
  122. NC1=MPOVA1.VPOCHA(/2)
  123. IF(NC1.EQ.NC) GOTO 30
  124. C
  125. C erreur pb dimension tableau voir routine adchpo
  126. C
  127. IRET=0
  128. SEGDES MPOVA1,MSOUP1,MCHPO1
  129. SEGSUP MSOUPO,MCHPOI
  130. CALL ERREUR(114)
  131. C on sort de la sous routine
  132. GOTO 666
  133. C
  134. 30 CONTINUE
  135. SEGINI MPOVAL
  136. IPOVAL=MPOVAL
  137. DO 40 IC=1,NC
  138. DO 40 IB=1,N
  139. VPOCHA(IB,IC)=MPOVA1.VPOCHA(IB,IC)
  140. 40 CONTINUE
  141. SEGDES MPOVAL,MPOVA1,MSOUPO,MSOUP1
  142. 50 CONTINUE
  143. C
  144. SEGDES MCHPOI,MCHPO1
  145. C on sort de la sous routine
  146. GOTO 666
  147. C
  148. C *** cas ou les pointeurs ip1 et ip2 sont differents
  149. C
  150. 60 CONTINUE
  151. MCHPO1=IP1
  152. MCHPO2=IP2
  153. IPO(1)=IP1
  154. IPO(2)=IP2
  155. SEGACT MCHPO1,MCHPO2
  156. MOT=MCHPO1.MTYPOI
  157. IF(MOT.NE.MCHPO2.MTYPOI) THEN
  158. MOT='ET OU +'
  159. ENDIF
  160. SEGDES MCHPO1,MCHPO2
  161. C
  162. C on verifie que les nbres de sous paquets sont egaux
  163. C
  164. MCHPO1=IP1
  165. MCHPO2=IP2
  166. SEGACT MCHPO1,MCHPO2
  167. NSOUP1=MCHPO1.IPCHP(/1)
  168. NSOUP2=MCHPO2.IPCHP(/1)
  169. NAT1 = MCHPO1.JATTRI(/1)
  170. NAT2 = MCHPO2.JATTRI(/1)
  171. IF(NSOUP1.EQ.NSOUP2) GO TO 75
  172. SEGDES MCHPO1,MCHPO2
  173. C traitement par la methode générale
  174. GO TO 300
  175. C
  176. C on regarde si les supports geometriques sont identiques a une
  177. C permutation pres
  178. C
  179. 75 SEGINI MTR3
  180. DO 100 IA=1,NSOUP1
  181. MSOUP1=MCHPO1.IPCHP(IA)
  182. SEGACT MSOUP1
  183. DO 80 IB=1,NSOUP2
  184. MSOUP2=MCHPO2.IPCHP(IB)
  185. SEGACT MSOUP2
  186. IF(MSOUP1.IGEOC.EQ.MSOUP2.IGEOC) GO TO 90
  187. SEGDES MSOUP2
  188. 80 CONTINUE
  189. C
  190. C il n y a pas egalite des supports geometriques a une permutation
  191. C pres
  192. C
  193. SEGDES MSOUP1,MCHPO1,MCHPO2
  194. SEGSUP MTR3
  195. C traitement par la methode générale
  196. GO TO 300
  197. C
  198. 90 CONTINUE
  199. C la permutation est rangée dans index
  200. INDEX(**)=IB
  201. SEGDES MSOUP1,MSOUP2
  202. 100 CONTINUE
  203. C
  204. C *** cas ou il y a egalite des supports geometriques a une permutation
  205. C pres
  206. C
  207. NSOUPO=NSOUP1
  208. NAT=MAX(NAT1,NAT2,1)
  209. SEGINI MCHPOI
  210. JATTRI(1) = 1
  211. IRET=MCHPOI
  212. MTYPOI=MOT
  213. MOCHDE=MCHPO1.MOCHDE
  214. IFOPOI=IFOUR
  215. C
  216. DO 250 IA=1,NSOUP1
  217. MSOUP1=MCHPO1.IPCHP(IA)
  218. MSOUP2=MCHPO2.IPCHP(INDEX(IA))
  219. SEGACT MSOUP1,MSOUP2
  220. C
  221. C comparaison des noms des composantes
  222. C
  223. SEGINI MTR1,MTR4
  224. NC1=MSOUP1.NOCOMP(/2)
  225. NC2=MSOUP2.NOCOMP(/2)
  226. DO 130 IB=1,NC1
  227. IPCOM(**)=MSOUP1.NOCOMP(IB)
  228. IPHAR(**)=MSOUP1.NOHARM(IB)
  229. 130 CONTINUE
  230. SEGINI MTR2
  231. DO 160 IB=1,NC2
  232. DO 140 IC=1,NC1
  233. IF(MSOUP2.NOCOMP(IB).NE.MSOUP1.NOCOMP(IC)) GOTO 140
  234. IF(MSOUP2.NOHARM(IB).EQ.MSOUP1.NOHARM(IC)) GOTO 150
  235. 140 CONTINUE
  236. C la composante du IB n'est pas commune
  237. IPCOM(**)=MSOUP2.NOCOMP(IB)
  238. IPHAR(**)=MSOUP2.NOHARM(IB)
  239. IICOM(**)=IPCOM(/2)
  240. GO TO 160
  241. 150 CONTINUE
  242. C la composante est commune
  243. IICOM(**)=IC
  244. 160 CONTINUE
  245. C
  246. MPOVA1=MSOUP1.IPOVAL
  247. MPOVA2=MSOUP2.IPOVAL
  248. SEGACT MPOVA1,MPOVA2
  249. N1=MPOVA1.VPOCHA(/1)
  250. N2=MPOVA2.VPOCHA(/1)
  251. NCX1=MPOVA1.VPOCHA(/2)
  252. NCX2=MPOVA2.VPOCHA(/2)
  253. IF(NCX1.NE.NC1) GOTO 170
  254. IF(NCX2.NE.NC2) GOTO 170
  255. IF(N1.NE.N2) GOTO 170
  256. GOTO 180
  257. 170 CONTINUE
  258. SEGDES MSOUP1,MSOUP2,MPOVA1,MPOVA2,MCHPO1,MCHPO2
  259. SEGSUP MTR1,MTR2,MTR3,MCHPOI,MTR4
  260. C
  261. C pb avec les dimensions des chpoints
  262. C
  263. CALL ERREUR(114)
  264. IRET=0
  265. C on sort de la sous routine
  266. GOTO 666
  267. C
  268. 180 CONTINUE
  269. N=N1
  270. NC=IPCOM(/2)
  271. SEGINI MPOVAL
  272. NNIN = NC
  273. SEGINI MTR5
  274. C
  275. C mise a 0 de vpocha
  276. C
  277. * DO 190 IB=1,N
  278. * DO 190 IC=1,NC
  279. * VPOCHA(IB,IC)=ZERO
  280. * 190 CONTINUE
  281. C
  282. C addition des chpoints
  283. C
  284. C on place les valeurs de MCHPO1
  285. DO 210 IC=1,NC1
  286. DO 200 IB=1,N
  287. VPOCHA(IB,IC) = MPOVA1.VPOCHA(IB,IC)+VPOCHA(IB,IC)
  288. DMOY(IC) = DMOY(IC) + ABS(VPOCHA(IB,IC)/N)
  289. 200 CONTINUE
  290. IF (IIMPI.EQ.123)
  291. & write (IOIMP,*) ' ic dmoy(ic) ',ic,dmoy(ic)
  292. 210 CONTINUE
  293. C
  294. DO 230 IC=1,NC2
  295. IIC=IICOM(IC)
  296. DO 220 IB=1,N
  297. IF ( IIC .LE. NC1 ) THEN
  298. C il s'agit d'ne composante commune
  299. XX1 = MPOVA2.VPOCHA(IB,IC)
  300. XX2 = VPOCHA(IB,IIC)
  301. DXX = ABS ( XX2 - XX1)
  302. * SXX = MIN(ABS ( XX1 + XX2 ) / 2.D0,1.D-50)
  303. SXX = DMOY(IIC)
  304. IF ( (DXX / SXX) .LT. 1.D-3 .or.noer.eq.1) THEN
  305. VPOCHA(IB,IIC)= ( XX1 + XX2 ) / 2.D0
  306. ELSE
  307. C les valeurs des champ diffus au meme point sont différentes
  308. IF (IIMPI.EQ.123)
  309. & write (IOIMP,*) xx1,xx2,SXX
  310. CALL ERREUR(651)
  311. SEGDES MSOUP1,MSOUP2,MPOVA1,MPOVA2,MCHPO1,MCHPO2
  312. SEGSUP MTR1,MTR2,MTR3,MCHPOI,MTR4
  313. C on sort
  314. GOTO 666
  315. ENDIF
  316. ELSE
  317. C composantes non communes
  318. VPOCHA(IB,IIC) = MPOVA2.VPOCHA(IB,IC)+VPOCHA(IB,IIC)
  319. ENDIF
  320. 220 CONTINUE
  321. 230 CONTINUE
  322. SEGDES MPOVA2
  323. SEGDES MPOVA1
  324. C
  325. SEGINI MSOUPO
  326. DO 240 IB=1,NC
  327. NOCOMP(IB)=IPCOM(IB)
  328. NOHARM(IB)=IPHAR(IB)
  329. 240 CONTINUE
  330. SEGSUP MTR1,MTR2,MTR4
  331. IPOVAL=MPOVAL
  332. SEGDES MPOVAL
  333. IPT2=MSOUP1.IGEOC
  334. **** SEGINI,IPT1=IPT2
  335. IPT1 = IPT2
  336. IGEOC=IPT1
  337. SEGDES IPT1
  338. IPCHP(IA)=MSOUPO
  339. SEGDES MSOUPO,MSOUP1,MSOUP2
  340. SEGSUP MTR5
  341. 250 CONTINUE
  342. C
  343. SEGSUP MTR3
  344. SEGDES MCHPOI,MCHPO1,MCHPO2
  345. C on sort
  346. GOTO 666
  347. C
  348. C *** cas ou les supports geometriques ne sont pas identiques
  349. C a une permutation pres
  350. C
  351. 300 CONTINUE
  352. C
  353. C **** a-t-on affaires a des champoints vides?
  354. C
  355. MCHPOI=IP1
  356. SEGACT MCHPOI
  357. NS1=IPCHP(/1)
  358. MCHPO2=IP2
  359. SEGACT MCHPO2
  360. NS2=MCHPO2.IPCHP(/1)
  361. IF(NS1*NS2.NE.0) GO TO 3001
  362. SEGDES MCHPOI,MCHPO2
  363. IF(NS1+NS2.NE.0) THEN
  364. C un seul des chpoints est vide
  365. IF(NS1.EQ.0) IP1=IP2
  366. CALL ECRCHA('GEOM')
  367. CALL ECROBJ('CHPOINT ',IP1)
  368. CALL COPIER
  369. CALL LIROBJ('CHPOINT',IRET,1,IRETOU)
  370. ELSE
  371. C les deux chpoints sont vides
  372. NSOUPO=0
  373. NAT=1
  374. SEGINI MCHPOI
  375. IRET = MCHPOI
  376. ENDIF
  377. RETURN
  378. 3001 CONTINUE
  379. C
  380. SEGINI MTRA,MTR1,MTR4
  381. C
  382. C mise a zero de nopoin
  383. C
  384. * DO 305 IA=1,XCOOR(/1)/(IDIM+1)
  385. * NOPOIN(IA)=0
  386. * 305 CONTINUE
  387. C
  388. MCHPOI=IP1
  389. SEGACT MCHPOI
  390. MSOUPO=IPCHP(1)
  391. SEGACT MSOUPO
  392. IPCOM(**)=NOCOMP(1)
  393. IPHAR(**)=NOHARM(1)
  394. SEGDES MSOUPO,MCHPOI
  395. NC=1
  396. IK=0
  397. C
  398. C creation de nopoin et de ipcom
  399. C
  400. DO 360 ICH=1,2
  401. MCHPOI=IPO(ICH)
  402. SEGACT MCHPOI
  403. NSOUPO=IPCHP(/1)
  404. C
  405. C boucle sur les sous references d 1 chpoint
  406. C
  407. DO 350 IA=1,NSOUPO
  408. MSOUPO=IPCHP(IA)
  409. SEGACT MSOUPO
  410. MELEME=IGEOC
  411. SEGACT MELEME
  412. NBNN =NUM(/1)
  413. NBELEM=NUM(/2)
  414. DO 310 IB=1,NBELEM
  415. K=NUM(1,IB)
  416. IF(NOPOIN(K).NE.0) GO TO 310
  417. IK=IK+1
  418. NOPOIN(K)=IK
  419. 310 CONTINUE
  420. SEGDES MELEME
  421. NCBBB=NOCOMP(/2)
  422. DO 330 IB=1,NCBBB
  423. NC=IPCOM(/2)
  424. DO 320 IC=1,NC
  425. IF(IPCOM(IC).NE.NOCOMP(IB)) GO TO 320
  426. IF(IPHAR(IC).EQ.NOHARM(IB)) GO TO 330
  427. 320 CONTINUE
  428. IPCOM(**)=NOCOMP(IB)
  429. IPHAR(**)=NOHARM(IB)
  430. NC=NC+1
  431. 330 CONTINUE
  432. SEGDES MSOUPO
  433. 350 CONTINUE
  434. SEGDES MCHPOI
  435. 360 CONTINUE
  436. C
  437. NNIN=NC
  438. NNNOE=IK
  439. SEGINI MTRAV
  440. C
  441. C initialisation a zero des tableaux
  442. C
  443. SEGINI MTR5
  444. C
  445. * DO 370 IB=1,NNNOE
  446. * DO 370 IA=1,NNIN
  447. * BB(IA,IB)=ZERO
  448. * IBIN(IA,IB)=0
  449. * IMOY(IA,IB) = 0
  450. * 370 CONTINUE
  451. C
  452. C creation de inco
  453. C
  454. DO 380 IA=1,NNIN
  455. INCO(IA)=IPCOM(IA)
  456. NHAR(IA)=IPHAR(IA)
  457. 380 CONTINUE
  458. C
  459. C creation de bb,ibin,igeo
  460. C
  461. DO 450 ICH=1,2
  462. MCHPOI=IPO(ICH)
  463. SEGACT MCHPOI
  464. NSOUPO=IPCHP(/1)
  465. DO 430 IA=1,NSOUPO
  466. MSOUPO=IPCHP(IA)
  467. SEGACT MSOUPO
  468. MELEME=IGEOC
  469. SEGACT MELEME
  470. MPOVAL=IPOVAL
  471. SEGACT MPOVAL
  472. NBELEM=NUM(/2)
  473. N=VPOCHA(/1)
  474. NC=VPOCHA(/2)
  475. NC1=NOCOMP(/2)
  476. C
  477. DO 420 IB=1,NC1
  478. DO 390 IC=1,NNIN
  479. IF(NOCOMP(IB).NE.IPCOM(IC)) GO TO 390
  480. IF(NOHARM(IB).EQ.IPHAR(IC)) GO TO 400
  481. 390 CONTINUE
  482. 400 CONTINUE
  483. DO 411 ID=1,NBELEM
  484. DMOY(IB)=DMOY(IB)+ABS(VPOCHA(ID,IB)/NBELEM)
  485. 411 CONTINUE
  486. DO 410 ID=1,NBELEM
  487. KI=NOPOIN(NUM(1,ID))
  488. IGEO(KI)=NUM(1,ID)
  489. IF ( IBIN(IC,KI) .EQ. 1) THEN
  490. C la valeur au point est defini dans les deux champs
  491. XX1 = BB(IC,KI)
  492. XX2 = VPOCHA(ID,IB)
  493. DXX = ABS ( XX2 - XX1 )
  494. SXX = DMOY(IB)
  495. IF ( DXX .LE. (1.D-4*SXX).or.noer.eq.1) THEN
  496. BB(IC,KI) = ( XX1 + XX2 ) / 2.D0
  497. ELSE
  498. C les valeurs des champs au meme point sont différentes
  499. IF (IIMPI.EQ.123)
  500. & write (IOIMP,*) xx1,xx2,sxx
  501. CALL ERREUR (651)
  502. SEGSUP MTRAV,MTRA,MTR1,MTR4,MTR5
  503. GOTO 666
  504. ENDIF
  505. ELSE
  506. BB(IC,KI)=BB(IC,KI)+VPOCHA(ID,IB)
  507. * DMOY(IC) = DMOY(IC) +ABS(BB(IC,KI)/NNNOE)
  508. ENDIF
  509. IBIN(IC,KI)=1
  510. 410 CONTINUE
  511. 420 CONTINUE
  512. SEGDES MELEME,MSOUPO,MPOVAL
  513. 430 CONTINUE
  514. SEGDES MCHPOI
  515. 450 CONTINUE
  516. ITRAV=MTRAV
  517. SEGDES MTRAV
  518. C
  519. C reconstuction de la partition
  520. C
  521. CALL CRECHP(ITRAV,ICHPOI)
  522. C
  523. SEGSUP MTRAV,MTRA,MTR1,MTR4,MTR5
  524. IRET=ICHPOI
  525. MCHPOI=ICHPOI
  526. MCHPO1 = IP1
  527. MCHPO2 = IP2
  528. SEGACT MCHPO1,MCHPO2
  529. NAT1 = MCHPO1.JATTRI(/1)
  530. NAT2 = MCHPO2.JATTRI(/1)
  531. NAT=MAX(NAT1,NAT2,1)
  532. NSOUPO = IPCHP(/1)
  533. SEGADJ MCHPOI
  534. IRET=MCHPOI
  535. MTYPOI=MOT
  536. JATTRI(1) = 1
  537. C SEGACT MCHPOI*NOMOD
  538. C SEGDES MCHPO1,MCHPO2
  539. C
  540. 666 CONTINUE
  541. RETURN
  542. END
  543.  
  544.  
  545.  

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