Télécharger reduaf.eso

Retour à la liste

Numérotation des lignes :

  1. C REDUAF SOURCE CB215821 16/12/05 22:04:30 9237
  2.  
  3. * Reduction du champ par element jchelm sur le modele mmodtm
  4. * Le resultat est le champ par element mchel2 pour iret = 1 (KERRE=0),
  5. * sinon en cas d'erreur mchel2 = 0 pour iret = 0 (KERRE = num. erreur)
  6. * En sortie le champ mchel2 est un champ entierement actif.
  7.  
  8. SUBROUTINE REDUAF (jchelm,mmodtm,mchel2,istri,iret,KERRE)
  9.  
  10. IMPLICIT REAL*8(A-H,O-Z)
  11. IMPLICIT INTEGER (I-N)
  12.  
  13. -INC CCOPTIO
  14.  
  15. -INC SMCHAML
  16. -INC SMMODEL
  17.  
  18. -INC SMCOORD
  19. -INC SMELEME
  20. -INC SMLENTI
  21. -INC CCPRECO
  22.  
  23. segment izone(NZ,NSMOD)
  24. segment ismel(NZ,NSMOD)
  25.  
  26. segment icpr(nbpts)
  27. segment inde(jg)
  28.  
  29. CHARACTER*(NCONCH) conloc
  30. CHARACTER*(8) nomloc
  31. CHARACTER*(16) typloc
  32.  
  33. LOGICAL BDARCY
  34.  
  35. *G if (iimpi.eq.7203) then
  36. *G write(ioimp,*) 'Entree dans reduaf',mmodtm,jchelm
  37. *G call zpchel(jchelm,1)
  38. *G endif
  39.  
  40. iret = 1
  41. KERRE = 0
  42. mchel2 = 0
  43. BDARCY = .FALSE.
  44.  
  45. * -----------------------------------
  46. * Activation de tous les sous-modeles
  47. * -----------------------------------
  48. mmodel = mmodtm
  49.  
  50. SEGACT,mmodel
  51.  
  52. NSMOD = mmodel.kmodel(/1)
  53. DO is = 1, NSMOD
  54. imodel = mmodel.kmodel(is)
  55. SEGACT,imodel
  56. C Verification si on a un modele de DARCY actuellement incompatible
  57. C Car il se servent du MAILLAGE dans la TABLE DOMAINE et pas celui
  58. C contenu dans le MMODEL
  59. IDARC=0
  60. CALL PLACE(imodel.FORMOD,imodel.FORMOD(/2),IDARC,'DARCY')
  61. IF (IDARC .NE. 0) BDARCY = .TRUE.
  62. ENDDO
  63.  
  64. IF (BDARCY) THEN
  65. mchel2=jchelm
  66. C SEGACT complet de mchel2
  67. mchelm=mchel2
  68. segact,mchelm
  69. do j=1,ichaml(/1)
  70. mchaml=ichaml(j)
  71. segact,mchaml
  72. do k=1,ielval(/1)
  73. melval=ielval(k)
  74. segact,melval
  75. enddo
  76. enddo
  77. RETURN
  78. ENDIF
  79.  
  80.  
  81. C ---------------------------------------------------------------------
  82. C Verification que le MCHAML n'est pas deja dans le CCPRECO
  83. C ---------------------------------------------------------------------
  84. ith = 0
  85. CALL ooonth(ith)
  86. ith1 = ith + 1
  87.  
  88. ITAILL = NBPRRE(ith1)
  89. DO 201 IPREC1 = 1, ITAILL
  90. IF (PRECMO(IPREC1,ith1) .NE. mmodtm) GOTO 201
  91. IF ((PRECM1(IPREC1,ith1) .EQ. jchelm) .OR.
  92. & (PRECM2(IPREC1,ith1) .EQ. jchelm)) THEN
  93. mchel2 = PRECM2(IPREC1,ith1)
  94. C IF ((PRECM1(IPREC1,ith1) .EQ. jchelm)) THEN
  95. C PRINT *,'REDUAF',ith,'DEJA en MEMOIRE 1',IPREC1
  96. C ELSE
  97. C PRINT *,'REDUAF',ith,'DEJA en MEMOIRE 2',IPREC1
  98. C ENDIF
  99. C IF (IPREC1 .EQ. NPREDU) THEN
  100. C PRINT *,' CCPRECO trop petit :',IPREC1
  101. C CALL ERREUR(5)
  102. C ENDIF
  103. C CALL TRBAC
  104. C Mise a jour du preconditionnement dans CCPRECO : Deplacement en position 1
  105. IF (IPREC1 .EQ. 1) RETURN
  106. DO IPREC2 = IPREC1,2,-1
  107. PRECMO(IPREC2,ith1) = PRECMO(IPREC2 - 1,ith1)
  108. PRECM1(IPREC2,ith1) = PRECM1(IPREC2 - 1,ith1)
  109. PRECM2(IPREC2,ith1) = PRECM2(IPREC2 - 1,ith1)
  110. ENDDO
  111. PRECMO(1,ith1) = mmodtm
  112. PRECM1(1,ith1) = jchelm
  113. PRECM2(1,ith1) = mchel2
  114.  
  115. C SEGACT complet de mchel2
  116. mchelm=mchel2
  117. C PRINT *,mchelm
  118. segact,mchelm
  119. do j=1,ichaml(/1)
  120. mchaml=ichaml(j)
  121. segact,mchaml
  122. do k=1,ielval(/1)
  123. melval=ielval(k)
  124. segact,melval
  125. enddo
  126. enddo
  127. RETURN
  128. ENDIF
  129. 201 CONTINUE
  130.  
  131. 1 CONTINUE
  132.  
  133. C Mise a jour du preconditionnement dans CCPRECO
  134. C Glissement des valeurs vers le bas
  135. ITAILL = MIN(ITAILL + 1, NPREDU)
  136. NBPRRE(ith1) = ITAILL
  137. C PRINT *,'REDUAF : On a sauve ', ITH, ITAILL,' TRIPLETS'
  138. DO IPRECO = ITAILL,2,-1
  139. PRECMO(IPRECO,ith1) = PRECMO(IPRECO - 1,ith1)
  140. PRECM1(IPRECO,ith1) = PRECM1(IPRECO - 1,ith1)
  141. PRECM2(IPRECO,ith1) = PRECM2(IPRECO - 1,ith1)
  142. ENDDO
  143. PRECMO(1,ith1) = mmodtm
  144. PRECM1(1,ith1) = jchelm
  145.  
  146. mchelm = jchelm
  147. SEGACT,mchelm
  148. NZ = mchelm.imache(/1)
  149. L1 = mchelm.titche(/1)
  150. N3 = mchelm.infche(/2)
  151.  
  152. * -----------------------------------------
  153. * Cas tres particulier de MCHELM resultat :
  154. * -----------------------------------------
  155. IF (NZ.EQ.0) THEN
  156. *G if (iimpi.eq.7203) write(ioimp,*) 'CAS PARTICULIER NZ = 0'
  157. SEGINI,mchel2=mchelm
  158. SEGDES,mchelm
  159. *// SEGDES,mchel2
  160.  
  161. C Mise a jour du preconditionnement dans CCPRECO (Nouveau champ mchel2)
  162. PRECM2(1,ith1) = mchel2
  163. C PRINT *,'REDUAF',ith,mmodtm,jchelm,'PAS en MEMOIRE',mchel2
  164. RETURN
  165. ENDIF
  166.  
  167. * Cas formulation melange: on ajoute les modeles de phases en les activant
  168. nvim = 0
  169. IF (istri.NE.1) THEN
  170. DO is = 1, NSMOD
  171. imodel = kmodel(is)
  172. if (formod(1).EQ.'MELANGE') nvim = nvim + ivamod(/1)
  173. ENDDO
  174. IF (nvim.ne.0) THEN
  175. n1 = NSMOD + nvim
  176. SEGINI,mmode1
  177. nc = 0
  178. DO is = 1, NSMOD
  179. nc = nc + 1
  180. imodel = kmodel(is)
  181. mmode1.kmodel(nc) = imodel
  182. IF (formod(1).EQ.'MELANGE') then
  183. DO j = 1, ivamod(/1)
  184. IF (tymode(j).NE.'IMODEL') THEN
  185. KERRE = 21
  186. ELSE
  187. imode1 = ivamod(j)
  188. SEGACT,imode1
  189. nc = nc + 1
  190. mmode1.kmodel(nc) = imode1
  191. ENDIF
  192. ENDDO
  193. ENDIF
  194. ENDDO
  195. NSMOD = nc
  196. mmodel = mmode1
  197. IF (KERRE.NE.0) GOTO 9010
  198. ENDIF
  199. ENDIF
  200.  
  201. * Quelques initialisations :
  202. SEGINI,izone,ismel
  203. * mlent2 contient le nombre d'elements du maillage de chaque sous-modele.
  204. jg = NSMOD
  205. SEGINI,mlent2
  206. * mlent3 contient les intersections entre les maillages determinees :
  207. * mlent3.lect(i3) avec ismel(iz,is) = i3 correspond a l'intersection
  208. * entre le maillage du sous-modele is et la sous-zone iz du champ si
  209. * la valeur de i3 n'est pas nulle !
  210. jg = NSMOD * NZ
  211. SEGINI,mlent3
  212. NL3 = 0
  213. ISOZM = 0
  214.  
  215. *? SEGACT,mcoord
  216. nbpts = mcoord.xcoor(/1) / (idim+1) + 1
  217. np1 = nbpts - 1
  218. icpr = 0
  219. inde = 0
  220. *
  221. * Regroupement des zones directement appariees avec un sous-modele
  222. * Recherche des zones pouvant intersecter le maillage d'un sous-modele
  223. DO 100 is = 1, NSMOD
  224. imodel = mmodel.kmodel(is)
  225. IF (imodel.nefmod.EQ.22) GOTO 100
  226. meleme = imodel.imamod
  227. SEGACT,meleme
  228. itypm = meleme.itypel
  229. mlent2.lect(is) = meleme.num(/2)
  230. * On parcourt tous les NZ chamelem elementaires.
  231. DO 101 iz = 1, NZ
  232. conloc = mchelm.conche(iz)
  233. IF (conloc.NE.' ' .AND. conloc.NE.imodel.conmod) GOTO 101
  234. ixx = 0
  235. ipt1 = mchelm.imache(iz)
  236. * Correspondance maillage sous-zone et sous-modele
  237. IF (ipt1.EQ.meleme) THEN
  238. ixx = 1
  239. * Pas de correspondance directe, recherche intersection potentielle
  240. ELSE
  241. SEGACT,ipt1
  242. IF (ipt1.itypel.NE.itypm) GOTO 102
  243. * On va regarder si on n a pas deja evalue l'intersection :
  244. * (meme sous-modele is et sous-zone precedente ia<iz)
  245. DO ia = 1, iz-1
  246. IF (ipt1.EQ.mchelm.imache(ia)) THEN
  247. IF (ismel(ia,is).GT.0) THEN
  248. ixx = -2
  249. ismel(iz,is) = ismel(ia,is)
  250. GOTO 102
  251. ENDIF
  252. ENDIF
  253. ENDDO
  254. * (meme sous-zone iz et sous-modele ia<is)
  255. DO 103 ia = 1, is-1
  256. imode2 = mmodel.kmodel(ia)
  257. IF (imode2.nefmod.EQ.22) GOTO 103
  258. ipt2 = imode2.imamod
  259. IF (ipt2.EQ.meleme) THEN
  260. IF (ismel(iz,ia).GT.0) THEN
  261. ixx = -3
  262. ismel(iz,is) = ismel(iz,ia)
  263. GOTO 102
  264. ENDIF
  265. ENDIF
  266. 103 CONTINUE
  267. * On va determiner l'intersection de ipt1 et meleme :
  268. * On va creer un tableau (liste enti) de correspondance des elements
  269. * de ipt1 qui sont dans meleme
  270. nbno1 = ipt1.num(/1)
  271. nbel1 = ipt1.num(/2)
  272. IF (icpr.EQ.0) THEN
  273. SEGINI,icpr
  274. ELSE
  275. DO j = 1, nbpts
  276. icpr(j) = 0
  277. ENDDO
  278. ENDIF
  279. DO j = 1, nbel1
  280. DO m = 1, nbno1
  281. ib = ipt1.num(m,j)
  282. icpr(ib) = icpr(ib) + 1
  283. ENDDO
  284. ENDDO
  285. iprec = icpr(1)
  286. DO j = 2, np1
  287. iprec = iprec + icpr(j)
  288. icpr(j) = iprec
  289. ENDDO
  290. jg = icpr(np1)
  291. icpr(nbpts) = jg
  292. IF (inde.EQ.0) THEN
  293. SEGINI,inde
  294. ELSE
  295. IF (jg.GT.inde(/1)) THEN
  296. SEGADJ,inde
  297. ENDIF
  298. DO j = 1, jg
  299. inde(j) = 0
  300. ENDDO
  301. ENDIF
  302. DO j = 1, nbel1
  303. DO m = 1, nbno1
  304. ib = ipt1.num(m,j)
  305. ia = icpr(ib)
  306. inde(ia) = j
  307. icpr(ib) = ia - 1
  308. ENDDO
  309. ENDDO
  310. * Fin du travail preparatoire pour le maillage ipt1
  311. ipt2 = imodel.imamod
  312. nbno2 = ipt2.num(/1)
  313. nbel2 = ipt2.num(/2)
  314. c* ipt2 = imodel.imamod = meleme
  315. c* nbno2 = ipt2.num(/1) = nbno1
  316. c* nbel2 = ipt2.num(/2) = mlent2.lect(is)
  317. C on fabrique le mlenti de correspondance
  318. * on dimensionne au nombre d elements de ipt2 = sous-modele is
  319. jg = nbel2
  320. SEGINI,mlenti
  321. ibon = 0
  322. DO 110 iel2 = 1, nbel2
  323. ia = ipt2.num(1,iel2)
  324. ideb = icpr(ia)+1
  325. ifin = icpr(ia+1)
  326. IF (ifin.LT.ideb) GOTO 110
  327. DO 111 ib = ideb, ifin
  328. iel1 = inde(ib)
  329. DO j = 1, nbno1
  330. IF (ipt2.num(j,iel2).NE.ipt1.num(j,iel1)) GOTO 111
  331. ENDDO
  332. ibon = ibon + 1
  333. mlenti.lect(iel2) = iel1
  334. GOTO 110
  335. 111 CONTINUE
  336. 110 CONTINUE
  337. * On a determine s'il y avait bien intersection ou non
  338. IF (ibon.EQ.0) THEN
  339. ixx = 0
  340. ismel(iz,is) = 0
  341. SEGSUP,mlenti
  342. ELSE
  343. * Si on a trouve plus d'elements dans l'intersection que dans ipt1 !
  344. IF (ibon.GT.nbel1) THEN
  345. write(ioimp,*) 'REDUAF : Etiquette 11x intersection ?'
  346. ENDIF
  347. NL3 = NL3 + 1
  348. mlent3.lect(NL3) = mlenti
  349. ixx = -1
  350. ismel(iz,is) = NL3
  351. * SEGDES,mlenti
  352. ENDIF
  353. 102 CONTINUE
  354. SEGDES,ipt1
  355. ENDIF
  356. *G write(*,*) ' -',iz,is,ixx,ismel(iz,is)
  357. * Sous-zone du mchelm a traiter
  358. IF (ixx.NE.0) THEN
  359. DO 105 ia = 1, iz-1
  360. ib = izone(ia,is)
  361. IF (ib.EQ.0) GOTO 105
  362. IF (conche(ia).NE.conloc) GOTO 105
  363. DO k = 1, N3
  364. IF (k.NE.4) THEN
  365. IF (infche(ia,k).NE.infche(iz,k)) GOTO 105
  366. ENDIF
  367. ENDDO
  368. izone(iz,is) = ib
  369. GOTO 106
  370. 105 CONTINUE
  371. ISOZM = ISOZM + 1
  372. izone(iz,is) = ISOZM
  373. 106 CONTINUE
  374. ENDIF
  375. *G write(*,*) ' -',iz,is,ixx,izone(iz,is)
  376. 101 CONTINUE
  377. SEGDES,meleme
  378. 100 CONTINUE
  379.  
  380. IF (icpr.NE.0) SEGSUP,icpr
  381. IF (inde.NE.0) SEGSUP,inde
  382. *
  383. * ---------------------------------
  384. * Construction du MCHELM resultat :
  385. * ---------------------------------
  386. * Grace au traitement ci-dessus (boucle 105), ISOZM correspond a N1 :
  387. N1 = ISOZM
  388. L1 = mchelm.titche(/1)
  389. N3 = mchelm.infche(/2)
  390. SEGINI,mchel2
  391. mchel2.titche = mchelm.titche
  392. mchel2.ifoche = mchelm.ifoche
  393.  
  394. * Pour chaque sous-modele is, on regroupe les sous-zones du mchelm iz
  395. * associees (izone(iz,is) > 0) :
  396. DO 200 is = 1, NSMOD
  397. imodel = kmodel(is)
  398. IF (imodel.nefmod.EQ.22) GOTO 200
  399. ipt2 = imodel.imamod
  400. nbel2 = mlent2.lect(is)
  401.  
  402. DO 210 iz = 1, NZ
  403. in1 = izone(iz,is)
  404. IF (in1.LE.0) GOTO 210
  405. mchaml = mchelm.ichaml(iz)
  406. SEGACT,mchaml
  407. * Cas particulier du mchaml sans composante (on ne fait rien) :
  408. n21 = mchaml.ielval(/1)
  409. IF (n21.EQ.0) GOTO 211
  410. *
  411. IF (mchel2.imache(in1).EQ.0) THEN
  412. *G write(ioimp,*) ' Cas 1 :',mchel2.imache(in1)
  413. mchel2.conche(in1) = mchelm.conche(iz)
  414. mchel2.imache(in1) = ipt2
  415. DO k = 1, N3
  416. mchel2.infche(in1,k) = mchelm.infche(iz,k)
  417. ENDDO
  418. n22 = 0
  419. n2 = n22 + n21
  420. SEGINI,mcham2
  421. mchel2.ichaml(in1) = mcham2
  422. ELSE
  423. *G write(ioimp,*) ' Cas 2 :',mchel2.imache(in1)
  424. mcham2 = mchel2.ichaml(in1)
  425. n22 = mcham2.ielval(/1)
  426. n2 = n22 + n21
  427. SEGADJ,mcham2
  428. ENDIF
  429. mlenti = ismel(iz,is)
  430. IF (mlenti.GT.0) mlenti = mlent3.lect(mlenti)
  431. *G write(ioimp,*) ' :',iz,is,mlenti,n22,n21,n2
  432. DO i = 1, n21
  433. nomloc = mchaml.nomche(i)
  434. iplac = 0
  435. IF (n22.NE.0) THEN
  436. CALL PLACE(mcham2.nomche(1),n22,iplac,nomloc)
  437. ENDIF
  438. typloc = mchaml.typche(i)
  439. melval = mchaml.ielval(i)
  440. SEGACT,melval
  441. nbpi1 = MAX(melval.velche(/1),melval.ielche(/1))
  442. nbel1 = MAX(melval.velche(/2),melval.ielche(/2))
  443. IF (nbel1.GT.1) nbel1 = nbel2
  444. IF (iplac.EQ.0) THEN
  445. iplac = n22 + i
  446. mcham2.nomche(iplac) = nomloc
  447. mcham2.typche(iplac) = typloc
  448. IF (typloc.EQ.'REAL*8 ') THEN
  449. n1ptel = nbpi1
  450. n1el = nbel2
  451. n2ptel = 0
  452. n2el = 0
  453. ELSE
  454. n1ptel = 0
  455. n1el = 0
  456. n2ptel = nbpi1
  457. n2el = nbel2
  458. ENDIF
  459. SEGINI,melva2
  460. mcham2.ielval(iplac) = melva2
  461. ELSE
  462. * incompatibilite du type de composante entre champs
  463. IF (mcham2.typche(iplac).NE.typloc) THEN
  464. KERRE = 917
  465. MOTERR(1:4) = nomloc
  466. MOTERR(5:21) = typloc
  467. MOTERR(22:38) = mcham2.typche(iplac)
  468. GOTO 9000
  469. ENDIF
  470. melva2 = mcham2.ielval(iplac)
  471. SEGACT,melva2*MOD
  472. ENDIF
  473. * On ajoute melval a melva2 en tenant compte de l'intersection entre
  474. * les maillages (mlenti = 0 si maillage identique, >0 sinon)
  475. * "Extension" de melva2 si besoin par rapport a melval (appel a MELEXT)
  476. * sera effectuee en prealable de l'addition des valeurs dans MELADD.
  477. CALL MELADD(melva2,melval,typloc,mlenti,KERRE)
  478. SEGDES,melval
  479. IF (KERRE.NE.0) GOTO 9000
  480. ENDDO
  481. *
  482. 211 CONTINUE
  483. SEGDES,mchaml
  484. 210 CONTINUE
  485. 200 CONTINUE
  486.  
  487. C Compactage du champ resultat :
  488. C ------------------------------
  489. n1max = n1
  490. n1 = 0
  491. DO 310 i = 1, n1max
  492. mcham2 = mchel2.ichaml(i)
  493. IF (mcham2.EQ.0) GOTO 310
  494. * on compacte les composantes (s'il y en a bien sur !)
  495. n22 = mcham2.ielval(/1)
  496. IF (n22.EQ.0) GOTO 312
  497. n2 = 0
  498. DO 311 j = 1, n22
  499. melva2 = mcham2.ielval(j)
  500. IF (melva2.EQ.0) GOTO 311
  501. CALL COMRED(melva2)
  502. n2 = n2 + 1
  503. mcham2.nomche(n2) = mcham2.nomche(j)
  504. mcham2.typche(n2) = mcham2.typche(j)
  505. mcham2.ielval(n2) = melva2
  506. *// segdes,melva2
  507. 311 CONTINUE
  508. IF (n2.EQ.0) GOTO 310
  509. IF (n2.NE.n22) SEGADJ,mcham2
  510. 312 CONTINUE
  511. *// segdes,mcham2
  512. n1 = n1 + 1
  513. mchel2.conche(n1) = mchel2.conche(i)
  514. mchel2.imache(n1) = mchel2.imache(i)
  515. mchel2.ichaml(n1) = mcham2
  516. DO j = 1, N3
  517. mchel2.infche(n1,j) = mchel2.infche(i,j)
  518. ENDDO
  519. 310 CONTINUE
  520. IF (n1.NE.n1max) SEGADJ,mchel2
  521. *// segdes,mchel2
  522.  
  523. * On sort un champ vide s'il n'y a pas de zone commune :
  524. c* IF (n1.EQ.0) THEN
  525. c**G if (iimpi.eq.7203) write(ioimp,*) 'N1 = 0 apres traitement'
  526. c* KERRE = 21
  527. c* ENDIF
  528.  
  529. 9000 CONTINUE
  530. * Destruction des segments de travail devenus inutiles :
  531. DO i = 1, NL3
  532. mlenti = mlent3.lect(i)
  533. SEGSUP,mlenti
  534. ENDDO
  535. SEGSUP,izone,ismel,mlent3,mlent2
  536.  
  537. 9010 CONTINUE
  538. * Desactivation des entrees : modele et chamelem
  539. C CB215821 : On laisse ouvert le MMODEL et les IMODEL en sortie
  540. C DO i = 1, NSMOD
  541. C imodel = kmodel(i)
  542. C SEGDES,imodel
  543. C ENDDO
  544. IF (nvim.NE.0) SEGSUP,mmodel
  545. C mmodel = mmodtm
  546. C SEGDES,mmodel
  547.  
  548. SEGDES,mchelm
  549.  
  550. IF (KERRE.NE.0) THEN
  551. iret = 0
  552. mchel2 = 0
  553. ENDIF
  554.  
  555. *G if (iimpi.eq.7203) then
  556. *G write(ioimp,*) 'Sortie de reduaf',mchel2,kerre
  557. *G if (kerre.eq.0) call zpchel(mchel2,1)
  558. *G endif
  559.  
  560. C Mise a jour du preconditionnement dans CCPRECO (Nouveau champ mchel2)
  561. PRECM2(1,ith1) = mchel2
  562.  
  563. C PRINT *,'REDUAF',ith,mmodtm,jchelm,'PAS en MEMOIRE',mchel2
  564. RETURN
  565. END
  566.  
  567.  
  568.  

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