Télécharger exincr.eso

Retour à la liste

Numérotation des lignes :

  1. C EXINCR SOURCE GOUNAND 11/06/14 21:15:22 7005
  2. SUBROUTINE EXINCR(RI1,MLMOT1,MLMOT2,RI2)
  3. IMPLICIT REAL*8 (A-H,O-Z)
  4. IMPLICIT INTEGER (I-N)
  5. C***********************************************************************
  6. C NOM : EXINCR
  7. C DESCRIPTION : Extrait d'une rigidité la sous-matrice
  8. C d'inconnues primales et duales celles données
  9. C en argument CH*4
  10. * Dans le cas d'une relation, il serait prudent
  11. * de créer un nouveau multiplicateur mais alors si
  12. * on recombine les matrices extraites, on ne reobtient
  13. * pas l'originale. D'autre part, si on ne crée pas
  14. * de nouveau multiplicateur reso buggera sur la
  15. * matrice recombinée (présence d'un même multiplicateur
  16. * dans deux matrices différentes) On choisit de ne pas
  17. * créer de multiplicateur ici.
  18. C
  19. C
  20. C LANGAGE : ESOPE
  21. C AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  22. C mél : gounand@semt2.smts.cea.fr
  23. C***********************************************************************
  24. C APPELES :
  25. C APPELES (E/S) :
  26. C APPELES (BLAS) :
  27. C APPELES (CALCUL) :
  28. C APPELE PAR :
  29. C***********************************************************************
  30. C SYNTAXE GIBIANE :
  31. C RIG2 = EXTR RIG1 LMOT1 LMOT2 ;
  32. C
  33. C ENTREES :
  34. C ENTREES/SORTIES :
  35. C SORTIES :
  36. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  37. C***********************************************************************
  38. C VERSION : v1, 01/06/2011, version initiale
  39. C HISTORIQUE : v1, 01/06/2011, création
  40. C HISTORIQUE :
  41. C HISTORIQUE :
  42. C***********************************************************************
  43. -INC CCOPTIO
  44. -INC SMRIGID
  45. -INC SMELEME
  46. -INC SMLMOTS
  47. *
  48. * Segment de correspondance : noms d'inconnues
  49. *
  50. segment trav1
  51. integer cor1p(nligp2),cor1d(nligd2)
  52. endsegment
  53. segment trav2
  54. integer cor1n(nbnn2),cor2n(nbnn)
  55. endsegment
  56. *
  57. * Segment extensibles pour ranger les rigidités
  58. * créés
  59. *
  60. segment imatno(0)
  61. segment idesno(0)
  62. segment igeono(0)
  63. segment coerno(0)
  64. segment irg8no(0)
  65. segment irg7no(0)
  66. segment irg6no(0)
  67. segment irg5no(0)
  68. segment irg2no(0)
  69. *
  70. logical lsym,lrela
  71. *
  72. * Executable statements
  73. *
  74. SEGACT ri1
  75. SEGACT MLMOT1,MLMOT2
  76. *
  77. segini imatno,idesno,igeono,coerno,irg7no,irg5no
  78. segini irg2no,irg6no,irg8no
  79. do 1000 irig=1,ri1.irigel(/2)
  80. des1=ri1.irigel(3,irig)
  81. segact des1
  82. nligrp=des1.lisinc(/2)
  83. nligrd=des1.lisdua(/2)
  84. nligp2=nligrp
  85. nligd2=nligrd
  86. segini trav1
  87. * On regarde le descripteur et les inconnues concernées
  88. nligp2=0
  89. do 410 iligrp=1,nligrp
  90. do 420 incp=1,mlmot1.mots(/2)
  91. if (mlmot1.mots(incp).eq.des1.lisinc(iligrp)) then
  92. nligp2=nligp2+1
  93. cor1p(nligp2)=iligrp
  94. goto 410
  95. endif
  96. 420 continue
  97. 410 continue
  98. if (nligp2.eq.0) goto 999
  99. nligd2=0
  100. do 510 iligrd=1,nligrd
  101. do 520 incd=1,mlmot2.mots(/2)
  102. if (mlmot2.mots(incd).eq.des1.lisdua(iligrd)) then
  103. nligd2=nligd2+1
  104. cor1d(nligd2)=iligrd
  105. goto 510
  106. endif
  107. 520 continue
  108. 510 continue
  109. * WRITE(IOIMP,*) 'irig=',irig
  110. * WRITE(IOIMP,*) 'cor1p='
  111. * WRITE (IOIMP,2020) (cor1p(I),I=1,cor1p(/1))
  112. * WRITE(IOIMP,*) 'cor1d='
  113. * WRITE (IOIMP,2020) (cor1d(I),I=1,cor1d(/1))
  114. * 2020 FORMAT (20(2X,I4) )
  115. if (nligd2.eq.0) goto 999
  116. * On fait un cas particulier pour les matrices de relations
  117. * Il faut vérifier qu'on a bien gardé la symétrie.
  118. * Il faut également que la relation ait gardé le multiplicateur
  119. * et porte au moins sur un ddl autre !
  120. ipt1=ri1.irigel(1,irig)
  121. segact ipt1
  122. lrela=.false.
  123. if (ipt1.itypel.eq.22) then
  124. if (nligp2.eq.nligd2) then
  125. do 700 ilig=1,nligp2
  126. if (cor1p(ilig).ne.cor1d(ilig)) goto 710
  127. 700 continue
  128. lrela=.true.
  129. 710 continue
  130. endif
  131. endif
  132. * l'ordre de ces deux lignes est important
  133. if (lrela.and.cor1p(1).ne.1) goto 998
  134. if (lrela.and.nligp2.eq.1) goto 998
  135. coerno(**)=ri1.coerig(irig)
  136. irg2no(**) =ri1.irigel(2,irig)
  137. irg5no(**) =ri1.irigel(5,irig)
  138. irg6no(**) =ri1.irigel(6,irig)
  139. irg8no(**) =ri1.irigel(8,irig)
  140. * Le cas particulier facile, où on garde toute la matrice
  141. * élémentaire : on recopie la matrice originale
  142. if (nligp2.eq.nligrp.and.nligd2.eq.nligrd) then
  143. igeono(**)=ri1.irigel(1,irig)
  144. idesno(**) =des1
  145. imatno(**) =ri1.irigel(4,irig)
  146. irg7no(**) =ri1.irigel(7,irig)
  147. goto 999
  148. endif
  149. * Sinon, il faut réduire la matrice et son descripteur, voire même
  150. * la géométrie
  151. * Réduction du descripteur
  152. nligrp=nligp2
  153. nligrd=nligd2
  154. segini des2
  155. do iligp2=1,nligp2
  156. iligp1=cor1p(iligp2)
  157. des2.lisinc(iligp2)=des1.lisinc(iligp1)
  158. des2.noelep(iligp2)=des1.noelep(iligp1)
  159. enddo
  160. do iligd2=1,nligd2
  161. iligd1=cor1d(iligd2)
  162. des2.lisdua(iligd2)=des1.lisdua(iligd1)
  163. des2.noeled(iligd2)=des1.noeled(iligd1)
  164. enddo
  165. idesno(**)=des2
  166. * Réduction de la matrice
  167. xmatr1=ri1.irigel(4,irig)
  168. segact xmatr1
  169. nelrig=xmatr1.re(/3)
  170. segini xmatr2
  171. do ielrig=1,nelrig
  172. do iligp2=1,nligp2
  173. iligp1=cor1p(iligp2)
  174. do iligd2=1,nligd2
  175. iligd1=cor1d(iligd2)
  176. xmatr2.re(iligd2,iligp2,ielrig)=
  177. $ xmatr1.re(iligd1,iligp1,ielrig)
  178. enddo
  179. enddo
  180. enddo
  181. segdes xmatr2
  182. segdes xmatr1
  183. imatno(**) =xmatr2
  184. * Réduction éventuelle de la géométrie si des noeuds
  185. * ne sont pas référencés
  186. nbnn=ipt1.num(/1)
  187. nbnn2=ipt1.num(/1)
  188. segini trav2
  189. nbnn2=0
  190. do in=1,des2.noelep(/1)
  191. ic1=cor2n(des2.noelep(in))
  192. if (ic1.eq.0) then
  193. nbnn2=nbnn2+1
  194. cor1n(nbnn2)=des2.noelep(in)
  195. cor2n(des2.noelep(in))=nbnn2
  196. endif
  197. enddo
  198. do in=1,des2.noeled(/1)
  199. ic1=cor2n(des2.noeled(in))
  200. if (ic1.eq.0) then
  201. nbnn2=nbnn2+1
  202. cor1n(nbnn2)=des2.noeled(in)
  203. cor2n(des2.noeled(in))=nbnn2
  204. endif
  205. enddo
  206. * Si tous les noeuds ont été vus, on peut peut-être
  207. * garder la géométrie
  208. if (nbnn2.eq.nbnn) then
  209. if (ipt1.itypel.eq.22.and.(.not.lrela)) then
  210. segini,ipt2=ipt1
  211. ipt2.itypel=28
  212. segdes ipt2
  213. else
  214. ipt2=ipt1
  215. endif
  216. else
  217. * Sinon, on ne garde que les noeuds vus dans le nouveau
  218. * maillage et on modifie le descripteur
  219. nbnn=nbnn2
  220. nbelem=ipt1.num(/1)
  221. nbsous=0
  222. nbref=0
  223. segini ipt2
  224. ipt2.itypel=28
  225. if (lrela) ipt2.itypel=22
  226. do iel=1,nbelem
  227. do ibnn2=1,nbnn2
  228. ibnn=cor1n(ibnn2)
  229. ipt2.num(ibnn2,iel)=ipt1.num(ibnn,iel)
  230. enddo
  231. enddo
  232. do iligrp=1,des2.noelep(/1)
  233. des2.noelep(iligrp)=cor2n(des2.noelep(iligrp))
  234. enddo
  235. do iligrd=1,des2.noeled(/1)
  236. des2.noeled(iligrd)=cor2n(des2.noeled(iligrd))
  237. enddo
  238. segdes ipt2
  239. endif
  240. segsup trav2
  241. igeono(**)=ipt2
  242. segdes des2
  243. * Garde-t-on la symétrie ou l'antisymétrie ? Par défaut non, mais si
  244. * les ddl primaux et duaux que l'on garde sont les mêmes, oui !
  245. lsym=.false.
  246. if (ri1.irigel(7,irig).ne.2) then
  247. if (nligp2.eq.nligd2) then
  248. do 600 ilig=1,nligp2
  249. if (cor1p(ilig).ne.cor1d(ilig)) goto 610
  250. 600 continue
  251. lsym=.true.
  252. 610 continue
  253. endif
  254. endif
  255. if (lsym) then
  256. irg7no(**) =ri1.irigel(7,irig)
  257. else
  258. irg7no(**) =2
  259. endif
  260. 998 continue
  261. segdes ipt1
  262. 999 continue
  263. segdes des1
  264. segsup trav1
  265. 1000 CONTINUE
  266. SEGDES MLMOT1,MLMOT2
  267. *
  268. * il ne reste plus qu'a ranger dans ri2 les nouvelles raideurs engendrées
  269. *
  270. NRIGEL=imatno(/1)
  271. SEGINI RI2
  272. RI2.MTYMAT=RI1.MTYMAT
  273. SEGDES RI1
  274. DO irig=1,NRIGEL
  275. RI2.COERIG(irig)=coerno(irig)
  276. RI2.IRIGEL(1,irig)=igeono(irig)
  277. RI2.IRIGEL(2,irig)=irg2no(irig)
  278. RI2.IRIGEL(3,irig)=idesno(irig)
  279. RI2.IRIGEL(4,irig)=imatno(irig)
  280. RI2.IRIGEL(5,irig)=irg5no(irig)
  281. RI2.IRIGEL(6,irig)=irg6no(irig)
  282. RI2.IRIGEL(7,irig)=irg7no(irig)
  283. RI2.IRIGEL(8,irig)=irg8no(irig)
  284. ENDDO
  285. SEGDES RI2
  286. segsup imatno,idesno,igeono,coerno,irg7no,irg5no
  287. segsup irg2no,irg6no,irg8no
  288. *
  289. * Normal termination
  290. *
  291. * IRET=0
  292. RETURN
  293. *
  294. * End of subroutine EXINCR
  295. *
  296. END
  297.  
  298.  
  299.  
  300.  
  301.  
  302.  
  303.  

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