Télécharger exincr.eso

Retour à la liste

Numérotation des lignes :

  1. C EXINCR SOURCE JC220346 19/11/08 21:15:01 10366
  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. segadj,trav1
  117. * On fait un cas particulier pour les matrices de relations
  118. * Il faut vérifier qu'on a bien gardé la symétrie.
  119. * Il faut également que la relation ait gardé le multiplicateur
  120. * et porte au moins sur un ddl autre !
  121. ipt1=ri1.irigel(1,irig)
  122. segact ipt1
  123. lrela=.false.
  124. if (ipt1.itypel.eq.22) then
  125. if (nligp2.eq.nligd2) then
  126. do 700 ilig=1,nligp2
  127. if (cor1p(ilig).ne.cor1d(ilig)) goto 710
  128. 700 continue
  129. lrela=.true.
  130. 710 continue
  131. endif
  132. endif
  133. * l'ordre de ces deux lignes est important
  134. if (lrela.and.cor1p(1).ne.1) goto 998
  135. if (lrela.and.nligp2.eq.1) goto 998
  136. coerno(**)=ri1.coerig(irig)
  137. irg2no(**) =ri1.irigel(2,irig)
  138. irg5no(**) =ri1.irigel(5,irig)
  139. irg6no(**) =ri1.irigel(6,irig)
  140. irg8no(**) =ri1.irigel(8,irig)
  141. * Le cas particulier facile, où on garde toute la matrice
  142. * élémentaire : on recopie la matrice originale
  143. if (nligp2.eq.nligrp.and.nligd2.eq.nligrd) then
  144. igeono(**)=ri1.irigel(1,irig)
  145. idesno(**) =des1
  146. imatno(**) =ri1.irigel(4,irig)
  147. irg7no(**) =ri1.irigel(7,irig)
  148. goto 999
  149. endif
  150. * Sinon, il faut réduire la matrice et son descripteur, voire même
  151. * la géométrie
  152. * Réduction du descripteur
  153. nligrp=nligp2
  154. nligrd=nligd2
  155. segini des2
  156. do iligp2=1,nligp2
  157. iligp1=cor1p(iligp2)
  158. des2.lisinc(iligp2)=des1.lisinc(iligp1)
  159. des2.noelep(iligp2)=des1.noelep(iligp1)
  160. enddo
  161. do iligd2=1,nligd2
  162. iligd1=cor1d(iligd2)
  163. des2.lisdua(iligd2)=des1.lisdua(iligd1)
  164. des2.noeled(iligd2)=des1.noeled(iligd1)
  165. enddo
  166. idesno(**)=des2
  167. * Réduction de la matrice
  168. xmatr1=ri1.irigel(4,irig)
  169. segact xmatr1
  170. nelrig=xmatr1.re(/3)
  171. segini xmatr2
  172. do ielrig=1,nelrig
  173. do iligp2=1,nligp2
  174. iligp1=cor1p(iligp2)
  175. do iligd2=1,nligd2
  176. iligd1=cor1d(iligd2)
  177. xmatr2.re(iligd2,iligp2,ielrig)=
  178. & xmatr1.re(iligd1,iligp1,ielrig)
  179. enddo
  180. enddo
  181. enddo
  182. segdes xmatr2
  183. segdes xmatr1
  184. imatno(**) =xmatr2
  185. * Réduction éventuelle de la géométrie si des noeuds
  186. * ne sont pas référencés
  187. nbnn=ipt1.num(/1)
  188. nbnn2=ipt1.num(/1)
  189. segini trav2
  190. nbnn2=0
  191. do in=1,des2.noelep(/1)
  192. ic1=cor2n(des2.noelep(in))
  193. if (ic1.eq.0) then
  194. nbnn2=nbnn2+1
  195. cor1n(nbnn2)=des2.noelep(in)
  196. cor2n(des2.noelep(in))=nbnn2
  197. endif
  198. enddo
  199. do in=1,des2.noeled(/1)
  200. ic1=cor2n(des2.noeled(in))
  201. if (ic1.eq.0) then
  202. nbnn2=nbnn2+1
  203. cor1n(nbnn2)=des2.noeled(in)
  204. cor2n(des2.noeled(in))=nbnn2
  205. endif
  206. enddo
  207. * Si tous les noeuds ont été vus, on peut peut-être
  208. * garder la géométrie
  209. if (nbnn2.eq.nbnn) then
  210. if (ipt1.itypel.eq.22.and.(.not.lrela)) then
  211. segini,ipt2=ipt1
  212. ipt2.itypel=28
  213. segdes ipt2
  214. else
  215. ipt2=ipt1
  216. endif
  217. else
  218. * Sinon, on ne garde que les noeuds vus dans le nouveau
  219. * maillage et on modifie le descripteur
  220. nbnn=nbnn2
  221. nbelem=ipt1.num(/2)
  222. nbsous=0
  223. nbref=0
  224. segini ipt2
  225. ipt2.itypel=28
  226. if (lrela) ipt2.itypel=22
  227. do iel=1,nbelem
  228. do ibnn2=1,nbnn2
  229. ibnn=cor1n(ibnn2)
  230. ipt2.num(ibnn2,iel)=ipt1.num(ibnn,iel)
  231. enddo
  232. enddo
  233. do iligrp=1,des2.noelep(/1)
  234. des2.noelep(iligrp)=cor2n(des2.noelep(iligrp))
  235. enddo
  236. do iligrd=1,des2.noeled(/1)
  237. des2.noeled(iligrd)=cor2n(des2.noeled(iligrd))
  238. enddo
  239. segdes ipt2
  240. endif
  241. segsup trav2
  242. igeono(**)=ipt2
  243. segdes des2
  244. * Garde-t-on la symétrie ou l'antisymétrie ? Par défaut non, mais si
  245. * les ddl primaux et duaux que l'on garde sont les mêmes, oui !
  246. lsym=.false.
  247. if (ri1.irigel(7,irig).ne.2) then
  248. if (nligp2.eq.nligd2) then
  249. do 600 ilig=1,nligp2
  250. if (cor1p(ilig).ne.cor1d(ilig)) goto 610
  251. 600 continue
  252. lsym=.true.
  253. 610 continue
  254. endif
  255. endif
  256. if (lsym) then
  257. irg7no(**) =ri1.irigel(7,irig)
  258. else
  259. irg7no(**) =2
  260. endif
  261. 998 continue
  262. segdes ipt1
  263. 999 continue
  264. segdes des1
  265. segsup trav1
  266. 1000 CONTINUE
  267. SEGDES MLMOT1,MLMOT2
  268. *
  269. * il ne reste plus qu'a ranger dans ri2 les nouvelles raideurs engendrées
  270. *
  271. NRIGEL=imatno(/1)
  272. SEGINI RI2
  273. RI2.MTYMAT=RI1.MTYMAT
  274. SEGDES RI1
  275. DO irig=1,NRIGEL
  276. RI2.COERIG(irig)=coerno(irig)
  277. RI2.IRIGEL(1,irig)=igeono(irig)
  278. RI2.IRIGEL(2,irig)=irg2no(irig)
  279. RI2.IRIGEL(3,irig)=idesno(irig)
  280. RI2.IRIGEL(4,irig)=imatno(irig)
  281. RI2.IRIGEL(5,irig)=irg5no(irig)
  282. RI2.IRIGEL(6,irig)=irg6no(irig)
  283. RI2.IRIGEL(7,irig)=irg7no(irig)
  284. xmatri=ri2.irigel(4,irig)
  285. segact xmatri*mod
  286. xmatri.symre=ri2.irigel(7,irig)
  287. segdes xmatri
  288. RI2.IRIGEL(8,irig)=irg8no(irig)
  289. ENDDO
  290. SEGDES RI2
  291. segsup imatno,idesno,igeono,coerno,irg7no,irg5no
  292. segsup irg2no,irg6no,irg8no
  293. *
  294. * Normal termination
  295. *
  296. * IRET=0
  297. RETURN
  298. *
  299. * End of subroutine EXINCR
  300. *
  301. END
  302.  
  303.  
  304.  
  305.  
  306.  
  307.  
  308.  
  309.  
  310.  

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