Télécharger frig3c.eso

Retour à la liste

Numérotation des lignes :

  1. C FRIG3C SOURCE GF238795 18/02/01 21:15:41 9724
  2. SUBROUTINE FRIG3C (maifro,IPRIGI,IPCHJE,IPRIG2)
  3.  
  4. IMPLICIT INTEGER(I-N)
  5. IMPLICIT REAL*8(A-H,O-Z)
  6.  
  7. * Ce sous-programme calcule la raideur de frottement en 3D.
  8. * il a besoin pour cela du maillage de frottement et de la raideur
  9. * de contact (ou la raideur totale si c'est plus simple)
  10.  
  11. -INC CCOPTIO
  12. -INC CCREEL
  13. -INC CCGEOME
  14. -INC SMCHPOI
  15. -INC SMCOORD
  16. -INC SMELEME
  17. -INC SMRIGID
  18.  
  19. * icpr lx du contact ==> lx du frottement
  20. segment icpr1(xcoor(/1)/(idim+1))
  21. segment icpr2(xcoor(/1)/(idim+1))
  22. * xjeu champs de jeux initiaux
  23. segment xjeu(xcoor(/1)/(idim+1))
  24. *
  25. *
  26. * creation et remplissage de icpr
  27. *
  28. segact mcoord
  29. segini icpr1,icpr2
  30. nbp=0
  31. meleme=maifro
  32. segact meleme
  33. ipt1=meleme
  34. do is=1,max(1,lisous(/1))
  35. if (lisous(/1).ne.0) ipt1=lisous(is)
  36. segact ipt1
  37. if (ipt1.itypel.ne.22) then
  38. write (6,*) ' ipt1.itypel ',ipt1.itypel
  39. call erreur(16)
  40. endif
  41. if (ierr.ne.0) return
  42. do iel=1,ipt1.num(/2)
  43. il=ipt1.num(1,iel)
  44. if (icpr1(il).eq.0) then
  45. nbp=nbp+1
  46. icpr1(il)=ipt1.num(ipt1.num(/1)-1,iel)
  47. icpr2(il)=ipt1.num(ipt1.num(/1),iel)
  48. endif
  49. if(icpr1(il).ne.ipt1.num(ipt1.num(/1)-1,iel)) call erreur(5)
  50. enddo
  51. if (lisous(/1).ne.0) segdes ipt1
  52. enddo
  53. segdes meleme
  54.  
  55. * remplissage du champ de jeux (demi-frottement si jeu non nul)
  56.  
  57. segini xjeu
  58. mchpoi = IPCHJE
  59. segact mchpoi
  60. do 15 isoupo = 1, ipchp(/1)
  61. msoupo = ipchp(isoupo)
  62. segact msoupo
  63. mpoval=ipoval
  64. segact mpoval
  65. if (vpocha(/2).ne.1) call erreur(16)
  66. ipt8=igeoc
  67. segact ipt8
  68. if (ipt8.num(/1).ne.vpocha(/2)) call erreur(16)
  69. do 16 i=1,vpocha(/1)
  70. xjeu(ipt8.num(1,i))=vpocha(i,1)
  71. 16 continue
  72. segdes ipt8,mpoval,msoupo
  73. 15 continue
  74. segdes mchpoi
  75. IF (ierr.ne.0) return
  76.  
  77.  
  78.  
  79. *
  80. * boucle sur les raideurs de contact pour les transformer en frottement
  81. *
  82. mrigid=iprigi
  83. segact mrigid
  84. segini,ri1=mrigid
  85. do 10 ir=1,irigel(/2)
  86. ri1.irigel(1,ir)=0
  87. ri1.irigel(4,ir)=0
  88. meleme=irigel(1,ir)
  89. segact meleme
  90. ipt1=0
  91. if (itypel.ne.22) goto 10
  92. if (irigel(6,ir).eq.0) goto 10
  93. nbsous=0
  94. nbref=0
  95. nbnn=num(/1)
  96. nbelem=num(/2)*2
  97. segini,ipt1
  98. ipt1.itypel=22
  99. xmatri=irigel(4,ir)
  100. segact xmatri
  101. nligrd=re(/1)
  102. nligrp=re(/2)
  103. nelrig=re(/3)*2
  104. segini,xmatr1
  105. do iel=1,num(/2)
  106. * surface du premier triangle
  107. ip1=num(1,iel)
  108. ip2=num(2,iel)
  109. ip3=num(3,iel)
  110. xp1=xcoor((ip1-1)*4+1)
  111. yp1=xcoor((ip1-1)*4+2)
  112. zp1=xcoor((ip1-1)*4+3)
  113. xp2=xcoor((ip2-1)*4+1)
  114. yp2=xcoor((ip2-1)*4+2)
  115. zp2=xcoor((ip2-1)*4+3)
  116. xp3=xcoor((ip3-1)*4+1)
  117. yp3=xcoor((ip3-1)*4+2)
  118. zp3=xcoor((ip3-1)*4+3)
  119. xcr=max(abs(xp1-xp2),abs(xp2-xp3),abs(xp3-xp1),
  120. > abs(yp1-yp2),abs(yp2-yp3),abs(yp3-yp1),
  121. > abs(zp1-zp2),abs(zp2-zp3),abs(zp3-xp1))
  122. c xcr=sqrt(max(xcr1,xcr2,xcr3))
  123. xcr=0
  124. ** write (6,*) ' xcr xjeu ',xcr,xjeu(il)
  125. if (xjeu(il).gt.xcr) then
  126. xcof=0.5d0
  127. else
  128. xcof=1.d0
  129. endif
  130. * si mult de lagrange pas connu on a 0
  131. il=num(1,iel)
  132. * coefficient multiplicateur suivant le jeu par rapport a la taille de l'element
  133. * taille de l'element au carre
  134. iel1=2*iel-1
  135. if1=icpr1(il)
  136. ipt1.num(1,iel1)=if1
  137. ipt1.icolor(iel1)=icolor(iel)
  138. do in=2,ipt1.num(/1)
  139. ipt1.num(in,iel1)=num(in,iel)
  140. enddo
  141. sre=sqrt(re(1,nligrp-2,iel)**2+re(1,nligrp-1,iel)**2+
  142. > re(1,nligrp,iel)**2)
  143. if (abs(re(1,nligrp,iel)).lt.sre*0.5) then
  144. do ic=2,re(/2),3
  145. srev=sqrt(re(1,ic,iel)**2+re(1,ic+1,iel)**2+re(1,ic+2,iel)**2)
  146. srep=sqrt(re(1,ic,iel)**2+re(1,ic+1,iel)**2)
  147. * iel1 orthogonal a iel
  148. if (srep.lt.xpetit) srep=1.d0
  149. xmatr1.re(1,ic,iel1)=-re(1,ic+1,iel)*srev/srep
  150. xmatr1.re(1,ic+1,iel1)=re(1,ic,iel)*srev/srep
  151. xmatr1.re(1,ic+2,iel1)=0.d0
  152. enddo
  153. else
  154. do ic=2,re(/2),3
  155. srev=sqrt(re(1,ic,iel)**2+re(1,ic+1,iel)**2+re(1,ic+2,iel)**2)
  156. srep=sqrt(re(1,ic+2,iel)**2+re(1,ic,iel)**2)
  157. if (srep.lt.xpetit) srep=1.d0
  158. xmatr1.re(1,ic,iel1)=-re(1,ic+2,iel)*srev/srep
  159. xmatr1.re(1,ic+1,iel1)=0.d0
  160. xmatr1.re(1,ic+2,iel1)=re(1,ic,iel)*srev/srep
  161. enddo
  162. endif
  163. do ic=2,re(/1)
  164. xmatr1.re(ic,1,iel1)=xmatr1.re(1,ic,iel1)
  165. enddo
  166. iel2=2*iel
  167. if2=icpr2(il)
  168. ipt1.num(1,iel2)=if2
  169. ipt1.icolor(iel2)=icolor(iel)
  170. do in=2,ipt1.num(/1)
  171. ipt1.num(in,iel2)=num(in,iel)
  172. enddo
  173. sre=sqrt(re(1,nligrp-2,iel)**2+re(1,nligrp-1,iel)**2+
  174. > re(1,nligrp,iel)**2)
  175. do ic=2,re(/2),3
  176. if (sre.lt.xpetit) sre=1.d0
  177. * iel1 orthogonal a iel et iel1
  178. xmatr1.re(1,ic,iel2)=
  179. > (re(1,nligrp-1,iel)*xmatr1.re(1,ic+2,iel1)-
  180. > re(1,nligrp ,iel)*xmatr1.re(1,ic+1,iel1))/sre
  181. xmatr1.re(1,ic+1,iel2)=
  182. > (re(1,nligrp ,iel)*xmatr1.re(1,ic ,iel1)-
  183. > re(1,nligrp-2,iel)*xmatr1.re(1,ic+2,iel1))/sre
  184. xmatr1.re(1,ic+2,iel2)=
  185. > (re(1,nligrp-2,iel)*xmatr1.re(1,ic+1,iel1)-
  186. > re(1,nligrp-1,iel)*xmatr1.re(1,ic, iel1))/sre
  187. enddo
  188. do ic=2,re(/1)
  189. xmatr1.re(ic,1,iel2)=xmatr1.re(1,ic,iel2)
  190. enddo
  191. enddo
  192. segdes xmatri
  193. ri1.irigel(1,ir)=ipt1
  194. ri1.irigel(4,ir)=xmatr1
  195. ri1.irigel(6,ir)=2
  196.  
  197. 10 continue
  198. segdes mrigid
  199. *
  200. * boucle de compaction du resultat
  201. *
  202. mrigid=ri1
  203. irr=0
  204. do 100 ir=1,irigel(/2)
  205. meleme=irigel(1,ir)
  206. xmatri=irigel(4,ir)
  207. if (meleme.eq.0) goto 100
  208. ill=0
  209. do iel=1,num(/2)
  210. if (num(1,iel).ne.0) then
  211. ill=ill+1
  212. if (ill.ne.0) then
  213. do in=1,num(/1)
  214. num(in,ill)=num(in,iel)
  215. enddo
  216. icolor(ill)=icolor(iel)
  217. do ic=1,re(/1)
  218. re(1,ic,ill)=re(1,ic,iel)
  219. re(ic,1,ill)=re(ic,1,iel)
  220. enddo
  221. endif
  222. endif
  223. enddo
  224. if (ill.eq.0) goto 100
  225. if (ill.ne.num(/2)) then
  226. nbsous=0
  227. nbref=0
  228. nbnn=num(/1)
  229. nbelem=ill
  230. segadj meleme
  231. endif
  232. ** write (6,*) ' meleme sortie dans frig2c '
  233. ** call ecmail(meleme,0)
  234.  
  235.  
  236. irr=irr+1
  237. if (irr.ne.ir) then
  238. do ir1=1,irigel(/1)
  239. irigel(ir1,irr)=irigel(ir1,ir)
  240. enddo
  241. coerig(irr)=coerig(ir)
  242. endif
  243. 100 continue
  244. nrigel=irr
  245. if (irigel(/2).ne.irr) segadj mrigid
  246. iprig2=mrigid
  247. ** call prrigi(mrigid,1)
  248. segsup icpr1,icpr2,xjeu
  249. return
  250. end
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.  

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