Télécharger la plaquette Cast3M


fr en

Télécharger impo31.eso

Retour à la liste

Numérotation des lignes : 

impo31

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

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