Télécharger relasi.eso

Retour à la liste

Numérotation des lignes :

relasi
  1. C RELASI SOURCE PV 19/02/27 21:15:06 10128
  2. * simplification des relations
  3. * on elimine dans les relations les noeuds dont toutes les contributions sont < crit
  4. *
  5. subroutine relasi(mrigid)
  6. *
  7. implicit real*8 (a-h,o-z)
  8. logical travail
  9.  
  10. -INC PPARAM
  11. -INC CCOPTIO
  12. -INC SMRIGID
  13. -INC SMELEME
  14.  
  15. * segment de travail pour le regroupement des relations
  16. * garde indique en quel position du nouveau re on va
  17. * eleg indique la nouvelle position dans l'element d'un noeud garde
  18. segment trav
  19. integer garde(nligrd)
  20. integer eleg(nbnn)
  21. endsegment
  22. character*4 inco1,inco2
  23.  
  24. if (iimpi.eq.7) write(6,*) 'entree dans relasi'
  25. crit = 1d-10
  26. travail=.false.
  27. segact mrigid
  28. nrigel=0
  29. nrige1=0
  30. * ri1 va contenir les relations simplifiees
  31. * ri2 va contenir celles que l'on garde
  32. segini ri1
  33. segini,ri2=mrigid
  34. nbst=0
  35. do 200 ir=1,irigel(/2)
  36. ** write (6,*) ' boucle 200 ',ir
  37. meleme=irigel(1,ir)
  38. descr=irigel(3,ir)
  39. xmatri=irigel(4,ir)
  40. segact meleme,descr
  41. segact xmatri*mod
  42. if (itypel.ne.22) then
  43. goto 200
  44. endif
  45. ** fait au vol si besoin
  46. ipt2=meleme
  47. ** segini,ipt2=meleme
  48. nligrp=re(/1)
  49. nligrd=re(/2)
  50. nelrig=re(/3)
  51. xmatr2=xmatri
  52. ** segini,xmatr2
  53. ri2.irigel(1,ir)=ipt2
  54. ri2.irigel(4,ir)=xmatr2
  55. nbnn=num(/1)
  56. nbref=0
  57. nbsous=0
  58. nbsimp=0
  59. nbann=0
  60. nbnn=num(/1)
  61. segini trav
  62. do 100 iel=1,re(/3)
  63. * test deja traite
  64. if (ipt2.num(1,iel).eq.0) goto 100
  65. iok=1
  66. remax=0.d0
  67. * au cas ou la relation soit normalisee a autre chose que 1:
  68. do iv=2,re(/1)
  69. remax=max(remax,abs(re(iv,1,iel)))
  70. enddo
  71. critl = crit*remax
  72.  
  73. * test (et correction) qu'une inconnue apparait plusieurs fois dans la relation
  74. do iv1=1,nligrd
  75. ip1=num(noeled(iv1),iel)
  76. inco1=lisdua(iv1)
  77. do iv2=iv1+1,nligrd
  78. ip2=num(noeled(iv2),iel)
  79. inco2=lisdua(iv2)
  80. if (ip1.eq.ip2.and.inco1.eq.inco2) then
  81. re(1,iv1,iel)=re(1,iv1,iel)+re(1,iv2,iel)
  82. re(1,iv2,iel)=0.d0
  83. re(iv1,1,iel)=re(iv1,1,iel)+re(iv2,1,iel)
  84. re(iv2,1,iel)=0.d0
  85. endif
  86. enddo
  87. enddo
  88. * y a t il des simplifications dans cette relation?
  89. do iv=2,re(/1)
  90. if (abs(re(iv,1,iel)).gt.critl) iok=iok+1
  91. enddo
  92. if (iok.ne.re(/2)) then
  93. * il y a des termes a eliminer. on cree donc un nouveau paquet
  94. nligrp=iok
  95. nligrd=iok
  96. nelrig=re(/3)
  97. segini des1
  98. segini xmatr1
  99. nrige1=nrige1+1
  100. if (nrige1.gt.nrigel) then
  101. nrigel=nrige1+10000
  102. segadj ri1
  103. endif
  104. ri1.coerig(nrige1)=coerig(ir)
  105. ri1.irigel(2,nrige1)=irigel(2,ir)
  106. ri1.irigel(3,nrige1)=des1
  107. ri1.irigel(4,nrige1)=xmatr1
  108. do i=5,irigel(/1)
  109. ri1.irigel(i,nrige1)=irigel(i,ir)
  110. enddo
  111. * on recopie le lx
  112. ivc=1
  113. des1.lisinc(1)=lisinc(1)
  114. des1.lisdua(1)=lisdua(1)
  115. des1.noelep(1)=noelep(1)
  116. des1.noeled(1)=noeled(1)
  117. * recopie des valeurs en notant ou elles vont
  118. do iv=2,re(/1)
  119. if (abs(re(iv,1,iel)).gt.critl) then
  120. ivc=ivc+1
  121. garde(iv)=ivc
  122. xmatr1.re(1,ivc,1)=re(iv,1,iel)
  123. xmatr1.re(ivc,1,1)=re(iv,1,iel)
  124. des1.lisinc(ivc)=lisinc(iv)
  125. des1.lisdua(ivc)=lisdua(iv)
  126. des1.noelep(ivc)=noelep(iv)
  127. des1.noeled(ivc)=noeled(iv)
  128. else
  129. garde(iv)=0
  130. endif
  131. enddo
  132. * recreer le maillage support
  133. nbnn=num(/1)
  134. nbelem=re(/3)
  135. nbsous=0
  136. nbref=0
  137. segini ipt1
  138. ipt1.itypel=22
  139. ipt1.num(1,1)=num(1,iel)
  140. do in=2,nbnn
  141. ipt1.num(in,1)=-num(in,iel)
  142. enddo
  143. do ivc=2,des1.noelep(/1)
  144. ip=des1.noelep(ivc)
  145. ipt1.num(ip,1)=abs(ipt1.num(ip,1))
  146. enddo
  147. * dans ipt1 en negatif les noeud a supprimer
  148. inc=0
  149. * dans eleg la nouvelle position des noeuds
  150. do in=1,ipt1.num(/1)
  151. if (ipt1.num(in,1).gt.0) then
  152. inc = inc+1
  153. eleg(in)=inc
  154. ipt1.num(inc,1)=ipt1.num(in,1)
  155. do ivc=1,des1.noelep(/1)
  156. if (des1.noelep(ivc).eq.in) then
  157. des1.noelep(ivc)=inc
  158. des1.noeled(ivc)=inc
  159. endif
  160. enddo
  161. else
  162. eleg(in)=0
  163. endif
  164. enddo
  165. nbnn=inc
  166. segadj ipt1
  167. * maillage reconstitue et descr mis a jour
  168. ri1.irigel(1,nrige1)=ipt1
  169. * comme la relation a ete transferee dans un nouveau paquet on l'annule
  170. * on cree ipt2 et xmatr2 si ce n'a pas deja ete fait
  171. if (ipt2.eq.meleme) then
  172. segini,ipt2=meleme
  173. ri2.irigel(1,ir)=ipt2
  174. endif
  175. ipt2.num(1,iel)=0
  176. travail=.true.
  177. * maintenant, recherche des relations ayant les memes simplifications
  178. ** write (6,*) 'simplification',iel
  179. kel=1
  180. do 110 jel=iel+1,re(/3)
  181. remax=0.d0
  182. * au cas ou la relation soit normalisee a autre chose que 1:
  183. do iv=2,re(/1)
  184. remax=max(remax,abs(re(iv,1,jel)))
  185. enddo
  186. critl = crit*remax
  187. * deja traite?
  188. if (ipt2.num(1,jel).eq.0) goto 110
  189. do iv=2,re(/1)
  190. * meme simplification?
  191. if (abs(re(iv,1,jel)).gt.critl.and.garde(iv).eq.0)goto 110
  192. if (abs(re(iv,1,jel)).le.critl.and.garde(iv).ne.0)goto 110
  193. enddo
  194. ** write (6,*) 'meme simplification',jel
  195. * a rajouter dans xmatr1.re
  196. kel=kel+1
  197. do iv=2,re(/1)
  198. if (garde(iv).ne.0) then
  199. ivc=garde(iv)
  200. xmatr1.re(1,ivc,kel)=re(iv,1,jel)
  201. xmatr1.re(ivc,1,kel)=re(iv,1,jel)
  202. endif
  203. enddo
  204. do in=1,num(/1)
  205. if (eleg(in).ne.0) then
  206. ipt1.num(eleg(in),kel)=num(in,jel)
  207. endif
  208. enddo
  209. * annuler l'element
  210. ipt2.num(1,jel)=0
  211. 110 continue
  212. nelrig=kel
  213. segadj xmatr1
  214. nbelem=kel
  215. segadj ipt1
  216. nbsimp=nbsimp+kel
  217. endif
  218. 100 continue
  219. segsup trav
  220. * fin du do iel
  221. * maintenant, on retasse ri2.irigel(ir)
  222. if (ipt2.ne.meleme) then
  223. nligrp=re(/1)
  224. nligrd=re(/2)
  225. nelrig=re(/3)-nbsimp
  226. segini,xmatr2
  227. ri2.irigel(4,ir)=xmatr2
  228. nbelem=0
  229. nbnn=ipt2.num(/1)
  230. nbsous=0
  231. nbref=0
  232. do iel=1,ipt2.num(/2)
  233. if (ipt2.num(1,iel).ne.0) then
  234. nbelem=nbelem+1
  235. do in=1,ipt2.num(/1)
  236. ipt2.num(in,nbelem)=num(in,iel)
  237. enddo
  238. do ip=1,xmatr2.re(/2)
  239. do id=1,xmatr2.re(/1)
  240. xmatr2.re(id,ip,nbelem)=re(id,ip,iel)
  241. enddo
  242. enddo
  243. endif
  244. enddo
  245. if (nbelem.ne.nelrig) call erreur(5)
  246. segadj ipt2
  247. endif
  248. nbst = nbst + nbsimp
  249. * write(6,*) ' regroupement de relation ',nbsimp
  250. 200 continue
  251. * plus qu'a tout assembler dans le resultat.
  252. * au passage on elimine les matrices vides
  253. * mais avant, si on n'a rien eu a faire, on sort la matrice initiale qui est toujours dans mrigid
  254. if (.not.travail) return
  255.  
  256. nrige2=ri2.irigel(/2)
  257. if (iimpi.eq.7) write(6,*) ' nombre de relations simplifiees ',
  258. > nbst
  259. ** write(6,*) ' simplifications retenues annulees ',nrige1,nbann
  260. nrigel=nrige1+nrige2
  261. segadj ri2
  262. nrigel=0
  263. do ir=1,nrige2
  264. meleme=ri2.irigel(1,ir)
  265. segact meleme
  266. if (num(/2).ne.0) then
  267. nrigel=nrigel+1
  268. do is=1,ri2.irigel(/1)
  269. ri2.irigel(is,nrigel)=ri2.irigel(is,ir)
  270. enddo
  271. ri2.coerig(nrigel)=ri2.coerig(ir)
  272. endif
  273. enddo
  274. do ir=1,nrige1
  275. meleme=ri1.irigel(1,ir)
  276. segact meleme
  277. if (num(/2).ne.0) then
  278. nrigel=nrigel+1
  279. do is=1,ri2.irigel(/1)
  280. ri2.irigel(is,nrigel)=ri1.irigel(is,ir)
  281. enddo
  282. ri2.coerig(nrigel)=ri1.coerig(ir)
  283. endif
  284. enddo
  285. segsup ri1
  286. segadj ri2
  287.  
  288. mrigid = ri2
  289.  
  290. return
  291. end
  292.  
  293.  
  294.  
  295.  
  296.  
  297.  
  298.  
  299.  
  300.  
  301.  

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