Télécharger impo31.eso

Retour à la liste

Numérotation des lignes :

  1. C IMPO31 SOURCE BP208322 16/11/18 21:17:39 9177
  2.  
  3. C Operateur 'IMPO' 'MAIL' ( | 'SYME' | ) Surf1 Surf2 ( Crit) ;
  4. C 3D | 'MESC' |
  5. C | 'FAIB' |
  6.  
  7. SUBROUTINE IMPO31
  8. *
  9. * modif octobre 2012 on ne cree pas les noeuds support des mult ici
  10. * on le fera dynamiquement dans impo32
  11. *
  12. * Et en decembre 2015 on affaiblit la formulation en regroupant les multiplicateurs de lagrange par element
  13. * ce qui rend inutile le developpement precedent
  14. *
  15.  
  16. IMPLICIT INTEGER(I-N)
  17. IMPLICIT REAL*8 (A-H,O-Z)
  18.  
  19. -INC CCOPTIO
  20.  
  21. -INC SMCOORD
  22. -INC SMELEME
  23. -INC CCGEOME
  24. segment lagran(xcoor(/1)/4)
  25.  
  26. SEGMENT IAUXTB(NBAUX)
  27. PARAMETER ( X1s4 = 0.25000000000000000000000000000000000000000D0 ,
  28. & X1s6 = 0.16666666666666666666666666666666666666667D0 )
  29. *
  30. LOGICAL faible
  31. PARAMETER (ncle = 3)
  32. CHARACTER*4 mcle(ncle)
  33. DATA mcle / 'MESC','FAIB','SYME' /
  34. *
  35. idimp1 = IDIM + 1
  36. lagran = 0
  37. *
  38. * Lecture d'un mot-cle eventuel : SYME(trique) ou MESC ou FAIB(le)
  39. *
  40. mesc = 0
  41. call lirmot(mcle,ncle,mesc,0)
  42. if (ierr.ne.0) return
  43. *
  44. * Lecture des deux surfaces en contact
  45. *
  46. call lirobj('MAILLAGE',ipt1,1,iretou)
  47. if (ierr.ne.0) return
  48. call lirobj('MAILLAGE',ipt2,1,iretou)
  49. if (ierr.ne.0) return
  50. *
  51. * Lecture eventuelle de la distance limite max(crit)
  52. *
  53. crit = 1.e+35
  54. call lirree(crit,0,iretou)
  55. if (ierr.ne.0) return
  56. cri2 = crit*crit
  57. * crit = abs(crit)
  58. *D write(ioimp,*) ' critere ',crit,cri2
  59. *
  60. * Lecture eventuelle de la distance limite min(critmin)
  61. *
  62. critmin = 0.d0
  63. call lirree(critmin,0,iretou)
  64. if (ierr.ne.0) return
  65. crimi2 = critmin*critmin
  66. * crit = abs(crit)
  67. *D write(ioimp,*) ' critere ',crit,cri2
  68. *
  69. * Activation et quelques verifications sur les maillages des surfaces
  70. *
  71. segact ipt1,ipt2
  72. if (ipt1.itypel.ne.4) call erreur(16)
  73. if (ipt2.itypel.ne.4) call erreur(16)
  74. if (ierr.ne.0) goto 900
  75. *
  76. nbpts0 = xcoor(/1) / idimp1
  77. *
  78. nel1 = ipt1.num(/2)
  79. nel2 = ipt2.num(/2)
  80. * noe1 = ipt1.num(/1)
  81. noe2 = ipt2.num(/1)
  82. * En fait : noe1 = noe2 = 3 car elements de type TRI3 !
  83. *
  84. * Increment sur le nombre d'elements de contact retenus pour agrandir
  85. * les longueurs des segments (maillage, mcoord) de maniere adequate
  86. incnel = 1000000
  87. *
  88. faible = mesc.eq.2
  89. if (mesc .ne. 1) mesc = 0
  90. * creation si il y a lieu du point support des multiplicateurs de lagrange non instanciées
  91. if (ilgni.eq.0) then
  92. nbpts = nbpts0 +1
  93. nbpts0 = nbpts
  94. segadj,mcoord
  95. ilgni = nbpts
  96. endif
  97. *
  98. *-----------------------------------------------------------------------
  99. *- FORMULATION "STANDARD" DU CONTACT -
  100. *-----------------------------------------------------------------------
  101. IF (faible) GOTO 500
  102. *
  103. * Changement decembre 2015 on ne créé qu'une seule relation par noeud
  104. * en sommant les contributions des éléments en visàvis
  105. *
  106. * Creation lagran qui contient le mult associe au noeud
  107. *
  108. if (lagran.eq.0) segini lagran
  109. *
  110. *
  111. * Nombre de points supports (multiplicateurs de Lagrange) associes
  112. * au contact (1 par noeud de contact)
  113. * Ajustement du segment mcoord
  114. *
  115. nbpts = nbpts0 + (incnel * 1)
  116. * segadj,mcoord
  117. *
  118. * tableau auxiliaire pour accelerer le test element deja cree
  119. *
  120. nbaux = nbpts0
  121. *
  122. * Creation du maillage de contact (meleme)
  123. *
  124. nbsous = 0
  125. nbref = 0
  126. nbnn = 5
  127. nbelem = incnel
  128. segini,meleme
  129. itypel = 22
  130. *
  131. nblag = nbpts0
  132. nbcoo = (nblag-1)*idimp1
  133. nbelc = 0
  134. nbpts = nbpts + nel1
  135. segadj,mcoord
  136. *
  137. ipas = 1
  138. *
  139. 5 continue
  140. segini iauxtb
  141. *
  142. * Boucle sur les elements de la 1ere surface de contact
  143. *
  144. do 10 iel=1,nel1
  145. nbeini = nbelc
  146. ip1 = ipt1.num(1,iel)
  147. ip2 = ipt1.num(2,iel)
  148. ip3 = ipt1.num(3,iel)
  149. ipv = (ip1-1)*idimp1
  150. xp1 = xcoor(ipv+1)
  151. yp1 = xcoor(ipv+2)
  152. zp1 = xcoor(ipv+3)
  153. ipv = (ip2-1)*idimp1
  154. xp2 = xcoor(ipv+1)
  155. yp2 = xcoor(ipv+2)
  156. zp2 = xcoor(ipv+3)
  157. ipv = (ip3-1)*idimp1
  158. xp3 = xcoor(ipv+1)
  159. yp3 = xcoor(ipv+2)
  160. zp3 = xcoor(ipv+3)
  161. xg1 = xp1+xp2+xp3
  162. yg1 = yp1+yp2+yp3
  163. zg1 = zp1+zp2+zp3
  164. *
  165. * Boucle sur les noeuds de la 2eme surface de contact
  166. *
  167. do 20 jel=1,nel2
  168. do 20 jnn=1,noe2
  169. jp=ipt2.num(jnn,jel)
  170. * verification que pas relation du point sur lui meme
  171. if (jp.eq.ip1) goto 20
  172. if (jp.eq.ip2) goto 20
  173. if (jp.eq.ip3) goto 20
  174. * verification que pas deja la relation
  175. if (iauxtb(jp).eq.iel) goto 20
  176. iauxtb(jp) = iel
  177. ipv = (jp-1)*idimp1
  178. xp=xcoor(ipv+1)
  179. yp=xcoor(ipv+2)
  180. zp=xcoor(ipv+3)
  181. * test distance
  182. d1=(xp1-xp)**2+(yp1-yp)**2+(zp1-zp)**2
  183. if (d1.le.cri2.and.d1.ge.crimi2) goto 22
  184. d1=(xp2-xp)**2+(yp2-yp)**2+(zp2-zp)**2
  185. if (d1.le.cri2.and.d1.ge.crimi2) goto 22
  186. d1=(xp3-xp)**2+(yp3-yp)**2+(zp3-zp)**2
  187. if (d1.le.cri2.and.d1.ge.crimi2) goto 22
  188. goto 20
  189. *
  190. 22 continue
  191. *
  192. * Mise a jour de la dimension des segments si necessaire
  193. *
  194. if (nbelc.eq.nbelem) then
  195. nbelem = nbelem + incnel
  196. write (6,*) ' ajustement impo31 ',nbelem
  197. segadj,meleme
  198. endif
  199. *
  200. * on incremente mcoord avec le pt support du multiplicateur si il y a lieu
  201. * on range dans le meleme l'element de contact associe
  202. *
  203. if (lagran(jp).eq.0) then
  204. nblag = nblag + 1
  205. nbcoo = nbcoo + idimp1
  206. xcoor(nbcoo+1) = xp
  207. xcoor(nbcoo+2) = yp
  208. xcoor(nbcoo+2) = zp
  209. xcoor(nbcoo+3) = 0.d0
  210. lagran(jp)=nblag
  211. endif
  212. *
  213. * on range dans le meleme
  214. *
  215. nbelc = nbelc + 1
  216. num(1,nbelc) = lagran(jp)
  217. num(2,nbelc) = ip1
  218. num(3,nbelc) = ip2
  219. num(4,nbelc) = ip3
  220. num(5,nbelc) = jp
  221. *
  222. 20 continue
  223. *
  224. 10 continue
  225. segsup iauxtb
  226. *
  227. * On recommence en inversant les deux surfaces en cas de traitement
  228. * symetrique du contact (option MESC non activee)
  229. *
  230. if (ipas.eq.1.and.mesc.ne.1) then
  231. ipas=ipas+1
  232. iaux=ipt1
  233. ipt1=ipt2
  234. ipt2=iaux
  235. nel1 = ipt1.num(/2)
  236. nel2 = ipt2.num(/2)
  237. noe2 = ipt2.num(/1)
  238. goto 5
  239. endif
  240. segsup lagran
  241. *
  242. GOTO 1000
  243. *
  244. *-----------------------------------------------------------------------
  245. *- FORMULATION "FAIBLE" DU CONTACT -
  246. *-----------------------------------------------------------------------
  247. 500 CONTINUE
  248. *
  249. * Le nb d'elements maximal a creer pour le maillage de contact FAIBLE :
  250. * 1 entre chaque element de chacune des surfaces
  251. *
  252. * Creation du maillage de contact (meleme)
  253. *
  254. nbnn = 7
  255. nbsous = 0
  256. nbref = 0
  257. nbelem = incnel
  258. segini,meleme
  259. itypel = 22
  260. *
  261. * Nombre de points supports (multiplicateurs de Lagrange) associes
  262. * au contact : 1 par element de contact
  263. * Ajustement du segment mcoord
  264. *
  265. nbpts = nbpts0 + nel1
  266. segadj,mcoord
  267. *
  268. nblag = nbpts0
  269. nbcoo = (nblag-1)*idimp1
  270. nbelc = 0
  271. *
  272. * Boucle sur les elements du maillage de la 1ere surface de contact
  273. *
  274. do 510 iel = 1, nel1
  275. ip1 = ipt1.num(1,iel)
  276. ip2 = ipt1.num(2,iel)
  277. ip3 = ipt1.num(3,iel)
  278. ipv = (ip1-1)*idimp1
  279. xp1 = xcoor(ipv+1)
  280. yp1 = xcoor(ipv+2)
  281. zp1 = xcoor(ipv+3)
  282. ipv = (ip2-1)*idimp1
  283. xp2 = xcoor(ipv+1)
  284. yp2 = xcoor(ipv+2)
  285. zp2 = xcoor(ipv+3)
  286. ipv = (ip3-1)*idimp1
  287. xp3 = xcoor(ipv+1)
  288. yp3 = xcoor(ipv+2)
  289. zp3 = xcoor(ipv+3)
  290. xg1 = xp1 + xp2 + xp3
  291. yg1 = yp1 + yp2 + yp3
  292. zg1 = zp1 + zp2 + zp3
  293. *
  294. * on incremente mcoord avec le pt support des multiplicateurs
  295. *
  296. nblag = nblag + 1
  297. nbcoo = nbcoo + idimp1
  298. *
  299. xcoor(nbcoo+1) = (xg1) /3
  300. xcoor(nbcoo+2) = (yg1) /3
  301. xcoor(nbcoo+3) = (zg1) /3
  302. xcoor(nbcoo+4) = 0.
  303. *
  304. * Boucle sur les elements du maillage de la 2e surface de contact
  305. *
  306. do 520 jel = 1, nel2
  307. jp1 = ipt2.num(1,jel)
  308. jp2 = ipt2.num(2,jel)
  309. jp3 = ipt2.num(3,jel)
  310. ipv = (jp1-1)*idimp1
  311. xp4 = xcoor(ipv+1)
  312. yp4 = xcoor(ipv+2)
  313. zp4 = xcoor(ipv+3)
  314. ipv = (jp2-1)*idimp1
  315. xp5 = xcoor(ipv+1)
  316. yp5 = xcoor(ipv+2)
  317. zp5 = xcoor(ipv+3)
  318. ipv = (jp3-1)*idimp1
  319. xp6 = xcoor(ipv+1)
  320. yp6 = xcoor(ipv+2)
  321. zp6 = xcoor(ipv+3)
  322.  
  323. * test distance
  324. d1 = (xp1-xp4)**2+(yp1-yp4)**2+(zp1-zp4)**2
  325. d2 = (xp2-xp4)**2+(yp2-yp4)**2+(zp2-zp4)**2
  326. d3 = (xp3-xp4)**2+(yp3-yp4)**2+(zp3-zp4)**2
  327. if ((d1.le.cri2.or.d2.le.cri2.or.d3.le.cri2).and.
  328. > (d1.ge.crimi2.or.d2.ge.crimi2.or.d3.ge.crimi2)) goto 550
  329.  
  330. d1 = (xp1-xp5)**2+(yp1-yp5)**2+(zp1-zp5)**2
  331. d2 = (xp2-xp5)**2+(yp2-yp5)**2+(zp2-zp5)**2
  332. d3 = (xp3-xp5)**2+(yp3-yp5)**2+(zp3-zp5)**2
  333. if ((d1.le.cri2.or.d2.le.cri2.or.d3.le.cri2).and.
  334. > (d1.ge.crimi2.or.d2.ge.crimi2.or.d3.ge.crimi2)) goto 550
  335.  
  336. d1 = (xp1-xp6)**2+(yp1-yp6)**2+(zp1-zp6)**2
  337. d2 = (xp2-xp6)**2+(yp2-yp6)**2+(zp2-zp6)**2
  338. d3 = (xp3-xp6)**2+(yp3-yp6)**2+(zp3-zp6)**2
  339. if ((d1.le.cri2.or.d2.le.cri2.or.d3.le.cri2).and.
  340. > (d1.ge.crimi2.or.d2.ge.crimi2.or.d3.ge.crimi2)) goto 550
  341.  
  342. goto 520
  343. *
  344. 550 continue
  345. *
  346. * Mise a jour de la dimension des segments si necessaire
  347. *
  348. nbelc=nbelc+1
  349. if (nbelc.eq.nbelem) then
  350. nbelem = nbelem + incnel
  351. segadj,meleme
  352. endif
  353. *
  354. * on range dans le meleme
  355. *
  356. num(1,nbelc) = nblag
  357. num(2,nbelc) = ip1
  358. num(3,nbelc) = ip2
  359. num(4,nbelc) = ip3
  360. num(5,nbelc) = jp1
  361. num(6,nbelc) = jp3
  362. num(7,nbelc) = jp2
  363. * on note la formulation dans icolor
  364. icolor(nbelc)=1
  365. *
  366. 520 continue
  367. *
  368. 510 continue
  369. *
  370. * GOTO 1000
  371. *
  372. *-----------------------------------------------------------------------
  373. * FIN :
  374. *-----------------------------------------------------------------------
  375. 1000 CONTINUE
  376. *D write(ioimp,*) ' nbpts nbpts0 nblag',nbpts,nbpts0,nblag
  377. *D write(ioimp,*) ' nbelem nbelc',nbelem,nbelc
  378. *
  379. * On renvoie le resultat mais au prealable on ajuste les longueurs
  380. *
  381. if (nbelem.lt.nbelc) then
  382. call erreur(5)
  383. segsup,meleme
  384. goto 900
  385. else if (nbelem.gt.nbelc) then
  386. nbelem = nbelc
  387. segadj,meleme
  388. nbpts = nblag
  389. segadj,mcoord
  390. endif
  391. segdes,meleme
  392. call ecrobj('MAILLAGE',meleme)
  393.  
  394. 900 CONTINUE
  395. segdes,ipt1,ipt2
  396.  
  397. return
  398. end
  399.  
  400.  
  401.  
  402.  
  403.  
  404.  
  405.  
  406.  

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