Télécharger relasi.eso

Retour à la liste

Numérotation des lignes :

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

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