Télécharger impofu.eso

Retour à la liste

Numérotation des lignes :

  1. C IMPOFU SOURCE PV 17/04/12 21:15:00 9391
  2. * reunion des relations portant sur le meme multiplicateur de lagrange
  3. * si mchpoi est nul, on se contente de nettoyer les termes petits
  4. *
  5. subroutine impofu(mrigid,mchpoi)
  6. implicit real*8 (a-h,o-z)
  7. -INC SMRIGID
  8. -INC SMCHPOI
  9. -INC SMELEME
  10. -INC CCOPTIO
  11. -INC SMCOORD
  12. -INC CCGEOME
  13. -INC CCREEL
  14. *
  15. segment icpr(xcoor(/1)/(idim+1))
  16. segment val((nbpo-1)*idim,nblag)
  17. segment ielem(nbpo,nblag)
  18. segment nbnl(nblag)
  19. segment dnorm(nblag)
  20. segact mcoord
  21. segini icpr
  22. ** call prrigi(mrigid,0)
  23. ** call ecchpo(mchpoi,0)
  24. segact mrigid
  25. nblag=0
  26. * regroupement des raideurs dans un grand tableau dimensionne au max
  27. * on commence par reperer les multiplicateurs de lagrange et leur destination
  28. do 10 i=1,irigel(/2)
  29. meleme=irigel(1,i)
  30. segact meleme
  31. do 20 iel=1,num(/2)
  32. lag=num(1,iel)
  33. if (icpr(lag).eq.0) then
  34. nblag=nblag+1
  35. icpr(lag)=nblag
  36. * noter si formulation faible
  37. if (icolor(iel).eq.1) icpr(lag)=-nblag
  38. endif
  39. 20 continue
  40. 10 continue
  41. * sauvegarde d'un descr pour avoir les noms d'inconnues
  42. des1=irigel(3,1)
  43. segact des1
  44. * creation melme et xmatri a la taille max
  45. nbpo=7
  46. segini val,ielem,nbnl,dnorm
  47. * remplissage, d'abort le mult
  48. do i=1,icpr(/1)
  49. if (icpr(i).ne.0) then
  50. ielem(1,abs(icpr(i)))=i
  51. endif
  52. enddo
  53. do 100 i=1,irigel(/2)
  54. meleme=irigel(1,i)
  55. xmatri=irigel(4,i)
  56. segact xmatri,meleme
  57. imod=0
  58. do 110 iel=1,num(/2)
  59. lag=num(1,iel)
  60. kel=abs(icpr(lag))
  61. remax=-1.
  62. do iv=2,re(/1)
  63. remax=max(remax,abs(re(iv,1,iel)))
  64. enddo
  65. do 120 in=2,num(/1)
  66. ip=num(in,iel)
  67. * la contribution n'est pas significative on saute le noeud
  68. do iv=1,idim
  69. if (abs(re(1+(in-2)*idim+iv,1,iel)).gt.1d-6*remax) goto 160
  70. enddo
  71. * pour assurer que la somme des termes est nulle, on corrige le dernier point:
  72. if (imod.eq.0) segact xmatri*mod
  73. imod=1
  74. do iv=1,idim
  75. re(1+(num(/1)-2)*idim+iv,1,iel)=
  76. > re(1+(num(/1)-2)*idim+iv,1,iel)+
  77. > re(1+(in-2)*idim+iv,1,iel)
  78. enddo
  79. if (iimpi.ne.0) write (6,*) ' noeud elimine dans relation',
  80. > in,ip,re(1+(in-2)*idim+1 ,1,iel),remax
  81. goto 120
  82. 160 continue
  83. do 130 ir=2,nbpo
  84. if (ielem(ir,kel).eq.0) goto 150
  85. if (ielem(ir,kel).ne.ip) goto 130
  86. * le noeud est deja dans l'element, on ajoute les valeurs
  87. do iv=1,idim
  88. val((ir-2)*idim+iv,kel)=val((ir-2)*idim+iv,
  89. > kel)+re(1+(in-2)*idim+iv,1,iel)
  90. enddo
  91. goto 120
  92. 130 continue
  93. * le noeud n'est pas dans l'element, on ajoute le noeud et les valeurs
  94. 150 continue
  95. if (ir.gt.nbpo) then
  96. nbpo=ir
  97. segadj ielem,val
  98. endif
  99. ielem(ir,kel)=ip
  100. do iv=1,idim
  101. val((ir-2)*idim+iv,kel)=re(1+(in-2)*idim+iv,1,iel)
  102. enddo
  103. 120 continue
  104. 110 continue
  105. segsup meleme,xmatri
  106. descr=irigel(3,i)
  107. if (i.ne.1) segsup descr
  108. 100 continue
  109. segsup mrigid
  110. *
  111. * eclatement en paquets de meme nb de noeuds
  112. * et renormalisation
  113. *
  114.  
  115. * nb noeuds par element
  116. do iel=1,ielem(/2)
  117. do in=1,ielem(/1)
  118. if (ielem(in,iel).eq.0) goto 200
  119. enddo
  120. 200 continue
  121. nbnl(iel)=in-1
  122. enddo
  123. * elimination elem en double si point ligne
  124. if (idim.eq.3) then
  125. do 300 iel=1,ielem(/2)
  126. if (nbnl(iel).ne.4) goto 300
  127. do 301 jel=iel+1,ielem(/2)
  128. if (nbnl(jel).ne.4) goto 301
  129. if (ielem(4,iel).ne.ielem(4,jel)) goto 301
  130. if (ielem(2,iel)*ielem(3,iel).ne.ielem(2,jel)*ielem(3,jel))
  131. > goto 301
  132. if (ielem(2,iel)+ielem(3,iel).ne.ielem(2,jel)+ielem(3,jel))
  133. > goto 301
  134. if (impi.ne.0) write (6,*) 'elimination elem seg en double'
  135. nbnl(jel)=0
  136. 301 continue
  137. 300 continue
  138. endif
  139. * elimination elem en double si point point
  140. do 310 iel=1,ielem(/2)
  141. if (nbnl(iel).ne.3) goto 310
  142. do 311 jel=iel+1,ielem(/2)
  143. if (nbnl(jel).ne.3) goto 311
  144. if (ielem(3,iel).ne.ielem(3,jel)) goto 311
  145. if (ielem(2,iel).ne.ielem(2,jel)) goto 311
  146. if (iimpi.ne.0) write(6,*) 'elimination elem poin en double'
  147. nbnl(jel)=0
  148. 311 continue
  149. 310 continue
  150.  
  151. nrigel=0
  152. segini mrigid
  153. ichole = 0
  154. imgeo1 = 0
  155. imgeo2 = 0
  156. isupeq = 0
  157. iforig = ifomod
  158. mtymat='RIGIDITE'
  159. *
  160. nbsous=0
  161. nbref=0
  162. do 250 iel=1,ielem(/2)
  163. if (nbnl(iel).eq.0) goto 250
  164. nbnn=nbnl(iel)
  165. nbelem=0
  166. do 255 jel=1,ielem(/2)
  167. if (nbnl(jel).eq.nbnn) nbelem=nbelem+1
  168. 255 continue
  169. segini meleme
  170. itypel=22
  171. nligrd=(nbnn-1)*idim+1
  172. nligrp=nligrd
  173. nelrig=nbelem
  174. segini xmatri
  175. *
  176. segini descr
  177. lisinc(1)=des1.lisinc(1)
  178. lisdua(1)=des1.lisdua(1)
  179. noelep(1)=1
  180. noeled(1)=1
  181. do inc=2,nligrd
  182. lisinc(inc)=des1.lisinc(mod(inc-2,idim)+2)
  183. lisdua(inc)=des1.lisdua(mod(inc-2,idim)+2)
  184. noelep(inc)=(inc-2)/idim+2
  185. noeled(inc)=(inc-2)/idim+2
  186. enddo
  187.  
  188. nrigel=nrigel+1
  189. segadj mrigid
  190. irigel(1,nrigel)=meleme
  191. irigel(3,nrigel)=descr
  192. irigel(4,nrigel)=xmatri
  193. irigel(6,nrigel)=1
  194. coerig(nrigel)=1.d0
  195. kel=0
  196. do 260 jel=1,ielem(/2)
  197. if (nbnl(jel).ne.nbnn) goto 260
  198. kel=kel+1
  199. do 265 in=1,nbnn
  200. num(in,kel)=ielem(in,jel)
  201. 265 continue
  202. if (icpr(num(1,kel)).lt.0) icolor(kel)=1
  203. if (idim.eq.2.or.nbnn.eq.2) then
  204. xnorm2=(abs(val(1,jel))+abs(val(3,jel)))**2+
  205. > (abs(val(2,jel))+abs(val(4,jel)))**2
  206. else
  207. xnorm2=(abs(val(1,jel))+abs(val(4,jel))+abs(val(7,jel)))**2+
  208. > (abs(val(2,jel))+abs(val(5,jel))+abs(val(8,jel)))**2+
  209. > (abs(val(3,jel))+abs(val(6,jel))+abs(val(9,jel)))**2
  210. endif
  211. xnorm=sqrt(xnorm2)
  212. if (mchpoi.eq.0) xnorm=1.d0
  213. dnorm(abs(icpr(num(1,kel))))=xnorm
  214. do 270 in=1,idim*(nbnn-1)
  215. re(in+1,1,kel)=val(in,jel)/xnorm
  216. re(1,in+1,kel)=val(in,jel)/xnorm
  217. 270 continue
  218. nbnl(jel)=0
  219. 260 continue
  220. 250 continue
  221. segsup des1
  222. ** call prrigi(mrigid,0)
  223.  
  224. * et maintenant le second membre
  225. if (mchpoi.ne.0) then
  226. nc=1
  227. n=nblag
  228. segini mpoval
  229. nbnn=1
  230. nbelem=nblag
  231. segini meleme
  232. itypel=1
  233. segact mchpoi
  234. msoup1=ipchp(1)
  235. segact msoup1*mod
  236. mpova1=msoup1.ipoval
  237. segact mpova1
  238. ipt1=msoup1.igeoc
  239. segact ipt1
  240. do i=1,ipt1.num(/2)
  241. lag=ipt1.num(1,i)
  242. num(1,abs(icpr(lag)))=lag
  243. vpocha(abs(icpr(lag)),1)=mpova1.vpocha(i,1)/
  244. > dnorm(abs(icpr(lag)))+vpocha(abs(icpr(lag)),1)
  245. enddo
  246. segsup mpova1,ipt1
  247. msoup1.igeoc=meleme
  248. msoup1.ipoval=mpoval
  249. endif
  250. ***
  251. ** call ecchpo(mchpoi,0)
  252. segsup val,ielem,nbnl,dnorm
  253. return
  254. end
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.  
  262.  
  263.  
  264.  

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