Télécharger la plaquette Cast3M


fr en

Télécharger impos1.eso

Retour à la liste

Numérotation des lignes : 

impos1

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

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