Télécharger impos1.eso

Retour à la liste

Numérotation des lignes :

  1. C IMPOS1 SOURCE PV 15/12/02 21:15:13 8720
  2.  
  3. C Operateur 'IMPO' 'MAIL' ( | 'MESC' | ) Lig1 Lig2 ( Crit) ;
  4. C 2D | 'FAIB' |
  5. C | 'SYME' |
  6.  
  7. SUBROUTINE IMPOS1
  8.  
  9. IMPLICIT INTEGER(I-N)
  10. IMPLICIT REAL*8 (A-H,O-Z)
  11.  
  12. -INC CCOPTIO
  13. -INC SMCOORD
  14. -INC SMELEME
  15. segment lagran(xcoor(/1)/3)
  16.  
  17. PARAMETER ( X1s3 = 0.33333333333333333333333333333333333333333D0 ,
  18. & X1s4 = 0.25000000000000000000000000000000000000000D0 )
  19.  
  20. logical faible
  21. PARAMETER (ncle = 3)
  22. character*4 mcle(ncle)
  23. data mcle / 'MESC','FAIB','SYME' /
  24. *
  25. idimp1 = IDIM + 1
  26. lagran=0
  27. *
  28. * Lecture du mot eventuel SYME(trique) ou MESC ou FAIB(le)
  29. *
  30. mesc = 0
  31. call lirmot(mcle,ncle,mesc,0)
  32. if (ierr.ne.0) return
  33. *
  34. * Lecture des deux lignes en contact
  35. *
  36. call lirobj('MAILLAGE',ipt1,1,iretou)
  37. if (ierr.ne.0) return
  38. call lirobj('MAILLAGE',ipt2,1,iretou)
  39. if (ierr.ne.0) return
  40. *
  41. * Lecture eventuelle de la distance limite (crit)
  42. *
  43. crit=1.e+35
  44. call lirree(crit,0,iretou)
  45. if (ierr.ne.0) return
  46. cri2 = crit*crit
  47. * crit = abs(crit)
  48. *DBG write(ioimp,*) ' critere ',crit,cri2
  49. *
  50. * Petites verifications sur les lignes
  51. *
  52. segact,ipt1,ipt2
  53. if (ipt1.itypel.ne.2) call erreur(16)
  54. if (ipt2.itypel.ne.2) call erreur(16)
  55. if (ierr.ne.0) goto 900
  56. *
  57. * Quelques initialisations
  58. *
  59. segact mcoord
  60. nbpts0 = xcoor(/1) / idimp1
  61. *
  62. nel1 = ipt1.num(/2)
  63. nel2 = ipt2.num(/2)
  64. * noe1 = ipt1.num(/1)
  65. noe2 = ipt2.num(/1)
  66. * En fait : noe2 = noe1 = 2 car elements de type SEG2 !
  67. *
  68. * Increment sur le nombre d'elements de contact retenus pour agrandir
  69. * les longueurs des segments (maillage, mcoord) de maniere adequate
  70. incnel = 10000
  71. *
  72. faible = mesc.eq.2
  73. if (mesc.ne.1) mesc = 0
  74.  
  75. *=======================================================================
  76. * Formulation forte ("standard") du contact
  77. *=======================================================================
  78. if (faible) goto 200
  79. *
  80. * Changement decembre 2015 on ne créé qu'une seule relation par noeud
  81. * en sommant les contributions des éléments en visàvis
  82. *
  83. * Creation lagran qui contient le mult associe au noeud
  84. *
  85. if (lagran.eq.0) segini lagran
  86. *
  87. * Creation du maillage de contact (meleme)
  88. *
  89. nbsous = 0
  90. nbref = 0
  91. nbnn = 4
  92. nbelem = incnel
  93. segini,meleme
  94. itypel = 22
  95. *
  96. * Nombre de points supports (multiplicateurs de Lagrange) associes
  97. * au contact : 1 par noeud de la deuxieme ligne au max
  98. * Ajustement du segment mcoord
  99. *
  100. nbpts = nbpts0 + (incnel * 1)
  101. segadj,mcoord
  102. *
  103. nblag = nbpts0
  104. nbcoo = (nblag-1)*idimp1
  105. nbelc = 0
  106. *
  107. ipas = 1
  108. *
  109. 5 continue
  110. *
  111. * On boucle sur les elements de la 1ere ligne de contact
  112. *
  113. do 10 iel = 1, nel1
  114. nbeini = nbelc
  115. ip1=ipt1.num(1,iel)
  116. ip2=ipt1.num(2,iel)
  117. xp1=xcoor((ip1-1)*idimp1+1)
  118. yp1=xcoor((ip1-1)*idimp1+2)
  119. xp2=xcoor((ip2-1)*idimp1+1)
  120. yp2=xcoor((ip2-1)*idimp1+2)
  121. xg1 = xp1+xp2
  122. yg1 = yp1+yp2
  123. *
  124. * Boucle sur les noeuds de la 2eme ligne de contact
  125. *
  126. do 20 jel=1,nel2
  127. do 20 jnn=1,noe2
  128. jp=ipt2.num(jnn,jel)
  129. * verification que pas relation du point sur lui meme
  130. if (jp.eq.ip1) goto 20
  131. if (jp.eq.ip2) goto 20
  132. * verification que pas deja la relation
  133. do 21 irela = nbeini+1, nbelc
  134. if (jp.eq.num(4,irela)) goto 20
  135. 21 continue
  136. xp=xcoor((jp-1)*idimp1+1)
  137. yp=xcoor((jp-1)*idimp1+2)
  138. * test distance
  139. d1 = (xp1-xp)**2+(yp1-yp)**2
  140. d2 = (xp2-xp)**2+(yp2-yp)**2
  141. if (cri2.lt.d1.and.cri2.lt.d2) goto 20
  142. *
  143. * Mise a jour de la dimension des segments si necessaire
  144. *
  145. if (nbelc.eq.nbelem) then
  146. nbelem = nbelem + incnel
  147. segadj,meleme
  148. nbpts = nbpts + incnel
  149. segadj,mcoord
  150. endif
  151. *
  152. * on incremente mcoord avec le pt support du multiplicateur si il y a lieu
  153. * on range dans le meleme l'element de contact associe
  154. *
  155. if (lagran(jp).eq.0) then
  156. nblag = nblag + 1
  157. nbcoo = nbcoo + idimp1
  158. xcoor(nbcoo+1) = xp
  159. xcoor(nbcoo+2) = yp
  160. xcoor(nbcoo+3) = 0.d0
  161. lagran(jp)=nblag
  162. endif
  163. *
  164. nbelc = nbelc + 1
  165. num(1,nbelc) = lagran(jp)
  166. num(2,nbelc) = ip1
  167. num(3,nbelc) = ip2
  168. num(4,nbelc) = jp
  169. *
  170. 20 continue
  171. *
  172. 10 continue
  173. *
  174. * On recommence en inversant les deux lignes en cas de traitement
  175. * symetrique du contact (option MESC non activee)
  176. *
  177. if (ipas.eq.1.and.mesc.ne.1) then
  178. ipas=ipas+1
  179. iaux=ipt1
  180. ipt1=ipt2
  181. ipt2=iaux
  182. nel1 = ipt1.num(/2)
  183. nel2 = ipt2.num(/2)
  184. noe2 = ipt2.num(/1)
  185. goto 5
  186. endif
  187. *
  188. segsup lagran
  189. goto 1000
  190.  
  191. *=======================================================================
  192. * Formulation "faible" du contact
  193. *=======================================================================
  194. 200 continue
  195. *
  196. * Le nb d'elements maximal a creer pour le maillage de contact FAIBLE :
  197. * 1 entre chaque element de chacune des lignes
  198. *
  199. * Creation du maillage de contact (meleme)
  200. *
  201. nbsous= 0
  202. nbref = 0
  203. nbnn = 5
  204. nbelem = incnel
  205. segini,meleme
  206. itypel = 22
  207. *
  208. * Nomre de points supports (multiplicateurs de Lagrange) associes
  209. * au contact : 1 par element de contact
  210. * Ajustement du segment mcoord
  211. *
  212. nbpts = nbpts0 + nel1
  213. segadj,mcoord
  214. *
  215. nblag = nbpts0
  216. nbcoo = (nblag-1)*idimp1
  217. nbelc = 0
  218. *
  219. * Boucle sur les elements de la 1ere ligne de contact
  220. *
  221. do 210 iel = 1, nel1
  222. ip1 = ipt1.num(1,iel)
  223. ip2 = ipt1.num(2,iel)
  224. ipv = (ip1-1)*idimp1
  225. xp1 = xcoor(ipv+1)
  226. yp1 = xcoor(ipv+2)
  227. ipv = (ip2-1)*idimp1
  228. xp2 = xcoor(ipv+1)
  229. yp2 = xcoor(ipv+2)
  230. xg1 = xp1 + xp2
  231. yg1 = yp1 + yp2
  232. *
  233. * Boucle sur les elements du maillage de la 2e ligne de contact
  234. *
  235. nblag = nblag + 1
  236. *
  237. * on incremente mcoord avec les pts support des multiplicateurs
  238. * on range dans le meleme l'element de contact correspondant
  239. *
  240. nbcoo = nbcoo + idimp1
  241. *
  242. xcoor(nbcoo+1) = xg1/2.d0
  243. xcoor(nbcoo+2) = yg1/2.d0
  244. xcoor(nbcoo+3) = 0.d0
  245. *
  246. do 220 jel=1,nel2
  247. jp1 = ipt2.num(1,jel)
  248. jp2 = ipt2.num(2,jel)
  249. ipv = (jp1-1)*idimp1
  250. xp3 = xcoor(ipv+1)
  251. yp3 = xcoor(ipv+2)
  252. ipv = (jp2-1)*idimp1
  253. xp4 = xcoor(ipv+1)
  254. yp4 = xcoor(ipv+2)
  255.  
  256. * test distance
  257. d1 = (xp1-xp3)**2+(yp1-yp3)**2
  258. d2 = (xp2-xp3)**2+(yp2-yp3)**2
  259. if (d1.le.cri2.or.d2.le.cri2) goto 250
  260. d1 = (xp1-xp4)**2+(yp1-yp4)**2
  261. d2 = (xp2-xp4)**2+(yp2-yp4)**2
  262. if (d1.le.cri2.or.d2.le.cri2) goto 250
  263. goto 220
  264. *
  265. 250 continue
  266. *
  267. * Mise a jour de la dimension des segments si necessaire
  268. *
  269. nbelc = nbelc + 1
  270. if (nbelc.eq.nbelem) then
  271. nbelem = nbelem + incnel
  272. segadj,meleme
  273. endif
  274. *
  275. num(1,nbelc) = nblag
  276. num(2,nbelc) = ip1
  277. num(3,nbelc) = ip2
  278. num(4,nbelc) = jp1
  279. num(5,nbelc) = jp2
  280. icolor(nbelc)=1
  281. *
  282. 220 continue
  283. *
  284. 210 continue
  285. *
  286. * goto 1000
  287. *
  288. *=======================================================================
  289. * Fin
  290. *=======================================================================
  291. 1000 continue
  292. *DBG write(ioimp,*) ' nbpts nbpts0 nblag',nbpts,nbpts0,nblag
  293. *DBG write(ioimp,*) ' nbelem nbelc',nbelem,nbelc
  294. *
  295. * On renvoie le resultat mais au prealable on ajuste les longueurs
  296. *
  297. if (nbelem.lt.nbelc) then
  298. call erreur(5)
  299. segsup,meleme
  300. goto 900
  301. else if (nbelem.gt.nbelc) then
  302. nbelem = nbelc
  303. segadj,meleme
  304. nbpts = nblag
  305. segadj,mcoord
  306. endif
  307. segdes,meleme
  308. call ecrobj('MAILLAGE',meleme)
  309.  
  310. 900 continue
  311. segdes,ipt1,ipt2
  312.  
  313. return
  314. end
  315.  
  316.  
  317.  
  318.  
  319.  

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