Télécharger frig3c.eso

Retour à la liste

Numérotation des lignes :

  1. C FRIG3C SOURCE PV 17/03/16 21:15:01 9350
  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. xcr=sqrt(max(xcr1,xcr2,xcr3))
  123. ** write (6,*) ' xcr xjeu ',xcr,xjeu(il)
  124. if (xjeu(il).gt.xcr) then
  125. xcof=0.5d0
  126. else
  127. xcof=1.d0
  128. endif
  129. * si mult de lagrange pas connu on a 0
  130. il=num(1,iel)
  131. * coefficient multiplicateur suivant le jeu par rapport a la taille de l'element
  132. * taille de l'element au carre
  133. iel1=2*iel-1
  134. if1=icpr1(il)
  135. ipt1.num(1,iel1)=if1
  136. ipt1.icolor(iel1)=icolor(iel)
  137. do in=2,ipt1.num(/1)
  138. ipt1.num(in,iel1)=num(in,iel)
  139. enddo
  140. sre=sqrt(re(1,nligrp-2,iel)**2+re(1,nligrp-1,iel)**2+
  141. > re(1,nligrp,iel)**2)
  142. if (abs(re(1,nligrp,iel)).lt.sre*0.5) then
  143. do ic=2,re(/2),3
  144. srev=sqrt(re(1,ic,iel)**2+re(1,ic+1,iel)**2+re(1,ic+2,iel)**2)
  145. srep=sqrt(re(1,ic,iel)**2+re(1,ic+1,iel)**2)
  146. * iel1 orthogonal a iel
  147. if (srep.lt.xpetit) srep=1.d0
  148. xmatr1.re(1,ic,iel1)=-re(1,ic+1,iel)*srev/srep
  149. xmatr1.re(1,ic+1,iel1)=re(1,ic,iel)*srev/srep
  150. xmatr1.re(1,ic+2,iel1)=0.d0
  151. enddo
  152. else
  153. do ic=2,re(/2),3
  154. srev=sqrt(re(1,ic,iel)**2+re(1,ic+1,iel)**2+re(1,ic+2,iel)**2)
  155. srep=sqrt(re(1,ic+2,iel)**2+re(1,ic,iel)**2)
  156. if (srep.lt.xpetit) srep=1.d0
  157. xmatr1.re(1,ic,iel1)=-re(1,ic+2,iel)*srev/srep
  158. xmatr1.re(1,ic+1,iel1)=0.d0
  159. xmatr1.re(1,ic+2,iel1)=re(1,ic,iel)*srev/srep
  160. enddo
  161. endif
  162. do ic=2,re(/1)
  163. xmatr1.re(ic,1,iel1)=xmatr1.re(1,ic,iel1)
  164. enddo
  165. iel2=2*iel
  166. if2=icpr2(il)
  167. ipt1.num(1,iel2)=if2
  168. ipt1.icolor(iel2)=icolor(iel)
  169. do in=2,ipt1.num(/1)
  170. ipt1.num(in,iel2)=num(in,iel)
  171. enddo
  172. sre=sqrt(re(1,nligrp-2,iel)**2+re(1,nligrp-1,iel)**2+
  173. > re(1,nligrp,iel)**2)
  174. do ic=2,re(/2),3
  175. if (sre.lt.xpetit) sre=1.d0
  176. * iel1 orthogonal a iel et iel1
  177. xmatr1.re(1,ic,iel2)=
  178. > (re(1,nligrp-1,iel)*xmatr1.re(1,ic+2,iel1)-
  179. > re(1,nligrp ,iel)*xmatr1.re(1,ic+1,iel1))/sre
  180. xmatr1.re(1,ic+1,iel2)=
  181. > (re(1,nligrp ,iel)*xmatr1.re(1,ic ,iel1)-
  182. > re(1,nligrp-2,iel)*xmatr1.re(1,ic+2,iel1))/sre
  183. xmatr1.re(1,ic+2,iel2)=
  184. > (re(1,nligrp-2,iel)*xmatr1.re(1,ic+1,iel1)-
  185. > re(1,nligrp-1,iel)*xmatr1.re(1,ic, iel1))/sre
  186. enddo
  187. do ic=2,re(/1)
  188. xmatr1.re(ic,1,iel2)=xmatr1.re(1,ic,iel2)
  189. enddo
  190. enddo
  191. segdes xmatri
  192. ri1.irigel(1,ir)=ipt1
  193. ri1.irigel(4,ir)=xmatr1
  194. ri1.irigel(6,ir)=2
  195.  
  196. 10 continue
  197. segdes mrigid
  198. *
  199. * boucle de compaction du resultat
  200. *
  201. mrigid=ri1
  202. irr=0
  203. do 100 ir=1,irigel(/2)
  204. meleme=irigel(1,ir)
  205. xmatri=irigel(4,ir)
  206. if (meleme.eq.0) goto 100
  207. ill=0
  208. do iel=1,num(/2)
  209. if (num(1,iel).ne.0) then
  210. ill=ill+1
  211. if (ill.ne.0) then
  212. do in=1,num(/1)
  213. num(in,ill)=num(in,iel)
  214. enddo
  215. icolor(ill)=icolor(iel)
  216. do ic=1,re(/1)
  217. re(1,ic,ill)=re(1,ic,iel)
  218. re(ic,1,ill)=re(ic,1,iel)
  219. enddo
  220. endif
  221. endif
  222. enddo
  223. if (ill.eq.0) goto 100
  224. if (ill.ne.num(/2)) then
  225. nbsous=0
  226. nbref=0
  227. nbnn=num(/1)
  228. nbelem=ill
  229. segadj meleme
  230. endif
  231. ** write (6,*) ' meleme sortie dans frig2c '
  232. ** call ecmail(meleme,0)
  233.  
  234.  
  235. irr=irr+1
  236. if (irr.ne.ir) then
  237. do ir1=1,irigel(/1)
  238. irigel(ir1,irr)=irigel(ir1,ir)
  239. enddo
  240. coerig(irr)=coerig(ir)
  241. endif
  242. 100 continue
  243. nrigel=irr
  244. if (irigel(/2).ne.irr) segadj mrigid
  245. iprig2=mrigid
  246. ** call prrigi(mrigid,1)
  247. segsup icpr1,icpr2,xjeu
  248. return
  249. end
  250.  
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  

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