Télécharger frig3c.eso

Retour à la liste

Numérotation des lignes :

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

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