Télécharger frig3c.eso

Retour à la liste

Numérotation des lignes :

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

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