Télécharger la plaquette Cast3M


fr en

Télécharger accro.eso

Retour à la liste

Numérotation des lignes : 

accro

  1. C ACCRO     SOURCE    CB215821  24/04/12    21:15:02     11897          
  2.  
  3. C========================================================================
  4. C   ACCROCHAGE D UN MAILLAGE  SUR UN MASSIF PAR L INTERMEDIAIRE
  5. C   DES FONCTIONS DE FORME
  6. C   SUBROUTINE APPELLEE PAR RELA2
  7. C   JMB      OCTOBRE  1996
  8. C========================================================================
  9.       SUBROUTINE ACCRO(igliss)
  10.  
  11.       IMPLICIT INTEGER(I-N)
  12.       IMPLICIT REAL*8 (A-H,O-Z)
  13.  
  14.  
  15. -INC PPARAM
  16. -INC CCOPTIO
  17. -INC CCREEL
  18. -INC CCHAMP
  19. -INC CCGEOME
  20.  
  21. -INC SMCOORD
  22. -INC SMELEME
  23. -INC SMRIGID
  24. -INC SMLMOTS
  25. -INC SMMODEL
  26.  
  27.       SEGMENT iexclu
  28.         INTEGER IDEJVU(NBPTI)
  29.       ENDSEGMENT
  30.  
  31.       SEGMENT icount
  32.         INTEGER NBPTOT
  33.         INTEGER MAXPZ(nbsz)
  34.         INTEGER IEINT(3,NODES)
  35.         REAL*8 QQQ(3,NODES)
  36.       ENDSEGMENT
  37.  
  38.       SEGMENT ioubli
  39.         INTEGER iptou(nboub)
  40.       ENDSEGMENT
  41.  
  42.       SEGMENT mtrav
  43.         INTEGER inomil(NBNOE1)
  44.         REAL*8 QSI(3)
  45.         REAL*8 SHPP(6,NBNOE1),XEL(3,NBNOE1)
  46.       ENDSEGMENT
  47.  
  48.       SEGMENT itange
  49.         INTEGER iextr(nbpti)
  50.         REAL*8 xperp1(3,nbpti), xperp2(3,nbpti)
  51.       ENDSEGMENT
  52.  
  53.       SEGMENT icpr(nbpti)
  54.       SEGMENT imul(idim,nbpt)
  55.       SEGMENT icomu(nbpt)
  56.  
  57.       DIMENSION ICOR(100),IMEL(100),imtt(100)
  58.       DIMENSION xpo(3),CCOE(100),Xboite(2,3)
  59.  
  60. C     Xboite : boite englobante de l'element pour accelerer un peu les choses
  61.  
  62.       LOGICAL ltelq
  63.  
  64.       idimp1 = idim + 1
  65.  
  66.       MAIL1  = 0
  67.       MAIL2  = 0
  68.       IFROCA = 0
  69.       MLMOTX = 0
  70.       CCRIT  = 1.D-5
  71.  
  72. C**** LECTURE des MAILLAGE et MMODEL ***********************************
  73.       CALL LIROBJ('MAILLAGE',MAIL1,1,iretou)
  74.       CALL ACTOBJ('MAILLAGE',MAIL1,1)
  75.       IF (IERR.NE.0) RETURN
  76.       CALL LIROBJ('MAILLAGE',MAIL2,0,iret2)
  77.       IF (IERR.NE.0) RETURN
  78.  
  79.       IF (iret2.NE.0) THEN        
  80.         CALL ACTOBJ('MAILLAGE',MAIL2,1)
  81.         MELEMI = MAIL2
  82.         IPP2   = MAIL1
  83.       ELSE
  84.         CALL LIROBJ('MMODEL  ',IPMODL,1,IFROCA)
  85.         CALL ACTOBJ('MMODEL  ',IPMODL,1)
  86.         IF (IERR.NE.0) RETURN
  87. c IPP2 sera extrait du modele IPMODL
  88.         MELEMI = MAIL1
  89.         IPP2   = 0
  90.       ENDIF
  91.  
  92. C**** LECTURE DES NOMS D'INCONNUE ET DU CRITERE ************************
  93.       CALL LIROBJ('LISTMOTS',MLMOTX,0,iretou)
  94.       IF(iretou .EQ. 1)THEN
  95.         mlmots = MLMOTX
  96.         SEGACT,mlmots
  97.       ENDIF
  98.       IF (IERR.NE.0) RETURN
  99.  
  100.       CALL LIRREE(CCRIT,0,iretou)
  101.       IF (IERR.NE.0) RETURN
  102.       CCRIT = -ABS(CCRIT)
  103.  
  104. C LISTE DES NOMS D'INCONNUES
  105. C  Options par defaut
  106.       IF (MLMOTX.EQ.0) THEN
  107.         NINC = IDIM
  108.         IF (IFOUR.EQ.1) NINC = NINC + 1
  109.         JGN = 4
  110.         JGM = NINC
  111.         SEGINI,mlmots
  112. C   2D axisymetrique
  113.         IF (IFOUR.EQ.0)  THEN
  114.           mlmots.MOTS(1) = 'UR  '
  115.           mlmots.MOTS(2) = 'UZ  '
  116. C   2D Fourier
  117.         ELSE IF (IFOUR.EQ.1)  THEN
  118.           mlmots.MOTS(1) = 'UR  '
  119.           mlmots.MOTS(2) = 'UZ  '
  120.           mlmots.MOTS(3) = 'UT  '
  121.         ELSE
  122.           mlmots.MOTS(1) = 'UX  '
  123.           mlmots.MOTS(2) = 'UY  '
  124.           IF (IDIM.EQ.3) mlmots.MOTS(3) = 'UZ  '
  125.         ENDIF
  126. C  On a lu les noms d'inconnues
  127.       ELSE
  128.         IF (igliss.EQ.1) THEN
  129.           CALL ERREUR(19)
  130.           RETURN
  131.         ENDIF
  132.         mlmots = MLMOTX
  133.         NINC = mlmots.MOTS(/2)
  134.         IF (NINC.LT.1 .OR. NINC.GT.100) THEN
  135.           CALL ERREUR(19)
  136.           RETURN
  137.         ENDIF
  138.       ENDIF
  139. C Identification des duals  des ddl  en question
  140. c bp 18.01.2011 : mieux vaut utiliser LNOMDD renseigné dans cchamp+bdata
  141.       DO i = 1, NINC
  142.         k = 0
  143.         CALL PLACE(NOMDD,LNOMDD,k,mlmots.MOTS(i))
  144.         IF (k.EQ.0) THEN
  145.           MOTERR(1:8) = mlmots.MOTS(i)
  146.           CALL ERREUR(931)
  147.           RETURN
  148.         ENDIF
  149.         ICOR(i) = k
  150.       ENDDO
  151.       SEGDES,mlmots
  152.       IF (MLMOTX.EQ.0) SEGSUP,mlmots
  153.  
  154. C**** Nombre de noeuds du segment MCOORD *******************************
  155.       SEGACT,mcoord*NOMOD
  156.       NBPTI = nbpts
  157.  
  158. C**** VERIFICATION DU MAILLAGE MELEMI **********************************
  159. C   CE MAILLAGE DOIT ETRE MASSIF
  160. C   TETRAEDRES ET PYRAMIDES SONT EXCLUS
  161.       IPT1 = MELEMI
  162.       NBSZ1 = IPT1.LISOUS(/1)
  163.       nbsz  = MAX(1,NBSZ1)
  164.       NBNOE1 = 0
  165.       meleme = IPT1
  166.       DO iob = 1, nbsz
  167.         IF (NBSZ1.NE.0) meleme = ipt1.LISOUS(iob)
  168.         ity  = meleme.ITYPEL
  169.         nbnn = meleme.NUM(/1)
  170.         NBNOE1 = MAX(NBNOE1,nbnn)
  171.         IF ((IDIM.EQ.2.AND.ity.LT.4) .OR.
  172.      &      (IDIM.EQ.3.AND.ity.LT.14)) THEN
  173.           CALL ERREUR(16)
  174.           RETURN
  175.         ENDIF
  176.       ENDDO
  177.  
  178. C**** CAS OU UN MODELE A ETE LU * EXTRATION DE IPP2 ********************
  179.       IF (IFROCA.NE.0) THEN
  180.         mmode1 = IPMODL
  181.         segini,icpr
  182.         nbpt = 0
  183. C creation des l'ensemble des points a accrocher et des multiplicateurs
  184. C correspondant. imul(i,j) est le ieme multiplicateur a accrocher au j eme
  185. C point local ( icpr(ia)=j)
  186. C --- puis Assemblage du maillage global (tous les câbles)
  187.         ltelq=.false.
  188.         do 500 isouf = 1, mmode1.kmodel(/1)
  189.           imode1 = mmode1.kmodel(isouf)
  190.           if (imode1.tymode(1).ne.'MMODEL') then
  191.             call erreur(975)
  192.             return
  193.           endif
  194.           ipt6 = imode1.imamod
  195.           if (ipt6.itypel.ne.22) then
  196.             call erreur(19)
  197.             return
  198.           endif
  199.           do iel = 1,ipt6.num(/2)
  200.             ia = ipt6.num(2,iel)
  201.             if (icpr(ia).eq.0) then
  202.               nbpt = nbpt+1
  203.               icpr(ia) = nbpt
  204.             endif
  205.           enddo
  206.           if (imode1.tymode(1).ne.'MMODEL') then
  207.             call erreur(975)
  208.             return
  209.           endif
  210.           mmode2 = imode1.ivamod(1)
  211.           imode2 = mmode2.kmodel(1)
  212.           ipt7 = imode2.imamod
  213.           if (isouf.eq.1) then
  214.             ipt5 = ipt7
  215.           else
  216.             call fuse(ipt5,ipt7,ipp4,ltelq)
  217.             ipt5 = ipp4
  218.           endif
  219.   500   continue
  220.         segini,imul,icomu
  221.         do 502 isouf = 1, mmode1.kmodel(/1)
  222.           imode1 = mmode1.kmodel(isouf)
  223.           ipt6 = imode1.imamod
  224.           do iel = 1, ipt6.num(/2)
  225.             ia = ipt6.num(2,iel)
  226.             ib = icpr(ia)
  227.             icomu(ib) = icomu(ib)+1
  228.             imul(icomu(ib),ib) = ipt6.num(1,iel)
  229.           enddo
  230.   502   continue
  231.         segsup,icomu
  232.         IPP2 = ipt5
  233.       ENDIF
  234. C**** FIN DU PRE-TRAITEMENT DANS LE CAS D'UN MODELE ********************
  235.  
  236. C**** TRAVAIL **********************************************************
  237.  
  238.       IPT2 = IPP2
  239. C--------------------------------------
  240. *---- cas de l'accrochage glissant ----
  241.       IF (igliss.EQ.1) THEN
  242. *       on va attacher a chaque point un vecteur tangent au cable
  243. *       si on est arrive avec un modele, on prend le maillage du modele
  244. *       en reference
  245.         segini itange
  246.         nbsz2 = ipt2.lisous(/1)
  247.         meleme = ipt2
  248.         do iob = 1, MAX(1,nbsz2)
  249.           if (nbsz2.ne.0) meleme = ipt2.lisous(i)
  250.           if (meleme.itypel.ne.2) then
  251.             call erreur(946)
  252.             return
  253.           endif
  254.           do j = 1, meleme.num(/2)
  255.             ia = num(1,j)
  256.             ib = num(2,j)
  257.             if (iextr(ia).gt.2.or.iextr(ib).gt.2) then
  258.               call erreur(947)
  259.               return
  260.             endif
  261.             ja = (ia-1)*idimp1
  262.             jb = (ib-1)*idimp1
  263.             xab = xcoor(jb+1) - xcoor(ja+1)
  264.             yab = xcoor(jb+2) - xcoor(ja+2)
  265.             zab = XZERO
  266.             if (idim.eq.3) zab = xcoor(jb+3) - xcoor(ja+3)
  267.             xlo = sqrt(xab * xab + yab * yab + zab * zab)
  268.             if (xlo.eq.XZERO) then
  269.               call erreur(945)
  270.               return
  271.             endif
  272.             xab = xab/xlo
  273.             yab = yab/xlo
  274.             zab = zab/xlo
  275.             if (iextr(ia).eq.0) then
  276.               xperp1(1,ia) = xab
  277.               xperp1(2,ia) = yab
  278.               xperp1(3,ia) = zab
  279.             else
  280.               xxx = 1.d0
  281.               if (iextr(ia).eq.1) xxx = -1.d0
  282.               zz1 = xab*xxx + xperp1(1,ia)
  283.               zz2 = yab*xxx + xperp1(2,ia)
  284.               zz3 = zab*xxx + xperp1(3,ia)
  285.               xlo = sqrt(zz1*zz1 + zz2*zz2 + zz3*zz3)
  286.               if (xlo.le.XZERO) then
  287.                 call erreur(945)
  288.                 return
  289.               endif
  290.               xperp1(1,ia) = zz1/xlo
  291.               xperp1(2,ia) = zz2/xlo
  292.               xperp1(3,ia) = zz3/xlo
  293.             endif
  294.             if (iextr(ib).eq.0) then
  295.               xperp1(1,ib) = xab
  296.               xperp1(2,ib) = yab
  297.               xperp1(3,ib) = zab
  298.             else
  299.               xxx = 1.d0
  300.               if (iextr(ia).eq.2) xxx = -1.d0
  301.               zz1 = xab*xxx + xperp1(1,ib)
  302.               zz2 = yab*xxx + xperp1(2,ib)
  303.               zz3 = zab*xxx + xperp1(3,ib)
  304.               xlo = sqrt(zz1*zz1 + zz2*zz2 + zz3*zz3)
  305.               if (xlo.le.XZERO) then
  306.                 call erreur(945)
  307.                 return
  308.               endif
  309.               xperp1(1,ib) = zz1/xlo
  310.               xperp1(2,ib) = zz2/xlo
  311.               xperp1(3,ib) = zz3/xlo
  312.             endif
  313.             iextr(ia) = iextr(ia) + 1
  314.             iextr(ib) = iextr(ib) + 2
  315.           enddo
  316.         enddo
  317.  
  318.         do 3659 j = 1, nbpti
  319.           if (iextr(j).eq.0) goto 3659
  320.           if (abs(xperp1(1,j)).gt.2.D-2 .or.
  321.      &        abs(xperp1(2,j)).gt.2.D-2) then
  322.             xperp2(1,j) = -xperp1(2,j)
  323.             xperp2(2,j) =  xperp1(1,j)
  324.           else
  325.             xperp2(1,j) = -xperp1(3,j)
  326.             xperp2(3,j) =  xperp1(1,j)
  327.           endif
  328.           if (idim.eq.3) then
  329.             xlo = sqrt (  xperp2(1,j)**2 + xperp2(2,j)**2
  330.      $                  + xperp2(3,j)**2 )
  331.             xperp2(1,j) = xperp2(1,j) / xlo
  332.             xperp2(2,j) = xperp2(2,j) / xlo
  333.             xperp2(3,j) = xperp2(3,j) / xlo
  334.             xaa = xperp1(3,j)*xperp2(2,j)-xperp1(2,j)*xperp2(3,j)
  335.             xbb = xperp1(1,j)*xperp2(3,j)-xperp1(3,j)*xperp2(1,j)
  336.             xcc = xperp1(2,j)*xperp2(1,j)-xperp1(1,j)*xperp2(2,j)
  337.             xperp1(1,j) = xaa
  338.             xperp1(2,j) = xbb
  339.             xperp1(3,j) = xcc
  340.           else
  341.             xperp1(1,j) = xperp2(1,j)
  342.             xperp1(2,j) = xperp2(2,j)
  343.           endif
  344.  3659   continue
  345.       ENDIF
  346. *---- fin du traitement specifique accrochage glissant ----
  347. C----------------------------------------------------------
  348.  
  349. C  ON VA CHANGER LE MAILLAGE IPP2 A ACCROCHER EN ELEMENTS POI1
  350. C  S'IL NE L EST PAS DEJA
  351.       IPT2 = IPP2
  352.       IF (IPT2.ITYPEL.NE.1) THEN
  353.         CALL CHANGE(IPP2,1)
  354.         IPT2 = IPP2
  355.         SEGACT,IPT2
  356.       ENDIF
  357. c On a maintenant : IPT2 = IPP2 qui est maillage de POI1 !
  358.       NODES = IPT2.NUM(/2)
  359.  
  360.       SEGINI,iexclu
  361.       IEXX = iexclu
  362.  
  363. c On a toujours : ipt1 = MELEMI
  364.       ipt1 = MELEMI
  365.  
  366. C On veut savoir s'il n'y a pas des noeuds de IPP2 qui ne sont pas
  367. C deja des noeuds d'un element de MELEMI, car pour RELA ACCRO, on n'a pas
  368. C besoin de les traiter (et de construire une matrice de rigidite).
  369. C Pour cela, pour tous les noeuds (ip) de MELEMI, idejvu(ip) = 2.
  370.       meleme = ipt1
  371.       DO iob = 1, NBSZ
  372.         IF (NBSZ1.NE.0) meleme = ipt1.LISOUS(iob)
  373.         NBNN1 = meleme.num(/1)
  374.         NBEL1 = meleme.num(/2)
  375.         DO iem = 1, NBEL1
  376.           DO ipt = 1, NBNN1
  377.             ip = meleme.NUM(ipt,iem)
  378.             iexclu.IDEJVU(ip) = 2
  379.           ENDDO
  380.         ENDDO
  381.       ENDDO
  382.  
  383.       SEGINI,icount
  384.       ICC = icount
  385.       icount.NBPTOT = 0
  386.  
  387.       SEGINI,mtrav
  388.  
  389. C idejvu(ip) = 1 si le point ip dans un element de MELEMI
  390. C---- BOUCLE SUR LES SOUS ZONES ----
  391.       meleme = ipt1
  392.       ITR=0
  393.       DO iob = 1, NBSZ
  394.         IF (NBSZ1.NE.0) meleme = ipt1.LISOUS(iob)
  395. C  ICI LE ZONAGE (Les segments sont actifs en E/S)
  396.         CALL ZONAG(iob,meleme,IPP2,ICC,IEXX,ITR)
  397.       ENDDO
  398. C---- FIN DE LA BOUCLE SUR LES SOUS ZONES ----
  399.  
  400.       IACCRO = icount.NBPTOT
  401.  
  402. C---- TABLEAU des POINT PAS VUS ----
  403. C-----------------------------------------------------------
  404.       NBOUB = NODES - IACCRO
  405. *     WRITE(ioimp,*) 'nombre de points exterieurs ' ,nboub
  406.       SEGINI,ioubli
  407.       iuo = 0
  408.       DO j = 1, NODES
  409.         ip = IPT2.NUM(1,j)
  410.         IF (iexclu.IDEJVU(ip).EQ.0) THEN
  411.           iuo = iuo + 1
  412.           ioubli.IPTOU(iuo) = ip
  413.         ENDIF
  414.       ENDDO
  415. *     WRITE(ioimp,*) 'nombre de points exterieurs ' ,nboub
  416.       IF (iuo.EQ.0) GOTO 9000
  417.  
  418. C-----------------------------------------------------------
  419. C---- RATTRAPAGE EVENTUEL DE POINTS SUR LA LIMITE ----
  420.       IPT1 = MELEMI
  421.       krat = 0
  422. C---- boucle sur les sous zones ----
  423.       meleme = IPT1
  424.       DO iob = 1, NBSZ
  425.         IF (NBSZ1.NE.0) meleme = IPT1.LISOUS(iob)
  426.         NBNN1 = meleme.num(/1)
  427.         NBEL1 = meleme.num(/2)
  428.         IELE = meleme.itypel
  429.         kd = kdegre(IELE)
  430. C element quadratiques
  431.         if (kd.eq.3) then
  432.           inomi1 = 0
  433.           nso = nbsom(IELE)
  434.           iad = nspos(IELE)-1
  435.           do 762 i = 1, NBNN1
  436.             do j = 1, nso
  437.               iso = ibsom(iad+j)
  438.               if (i.eq.iso) goto 762
  439.             enddo
  440.             inomi1 = inomi1 + 1
  441.             mtrav.inomil(inomi1) = i
  442.  762      continue
  443. C  elements  lineaires
  444.         else if (kd.eq.2) then
  445.           inomi1 = NBNN1
  446.           do i = 1, NBNN1
  447.             mtrav.inomil(i) = i
  448.           enddo
  449.         endif
  450.  
  451. C...... boucle sur les elements de la sz .......
  452.         DO 1121 iem = 1, NBEL1
  453.           IF (IERR .NE. 0) RETURN
  454.           CALL DOXE(xcoor,idim,NBNN1,meleme.num,iem,mtrav.xel)
  455.  
  456. C         Boite englobante de l'element
  457.           DO ii=1,idim
  458.             DO jj=1,2
  459.               IF(jj .EQ. 1)THEN
  460.                 Xboite(1,ii)= XGRAND
  461.               ELSE
  462.                 Xboite(2,ii)=-XGRAND
  463.               ENDIF
  464.             ENDDO
  465.           ENDDO
  466.           DO INOE=1,NBNN1
  467.             IPNO=IDIMP1*(meleme.num(INOE,iem)-1)
  468.             DO ii=1,idim
  469.               XCO=XCOOR(IPNO + ii)
  470.               Xboite(1,ii)=MIN(XCO,Xboite(1,ii))
  471.               Xboite(2,ii)=MAX(XCO,Xboite(2,ii))
  472.             ENDDO
  473.           ENDDO
  474.  
  475. C         On ajoute le critere de rattrapage (Attention CCRIT est negatif)
  476.           DO ii=1,idim
  477.             DO jj=1,2
  478.               IF(jj .EQ. 1)THEN
  479.                 Xboite(1,ii)= Xboite(1,ii) + CCRIT
  480.               ELSE
  481.                 Xboite(2,ii)= Xboite(2,ii) - CCRIT
  482.               ENDIF
  483.             ENDDO
  484.           ENDDO
  485.  
  486. C.......boucle sur les points candidats au rattrapage ......
  487.         DO 1120 IPT = 1, IUO
  488.  
  489.           ip = ioubli.IPTOU(ipt)
  490.           IF (ip.LT.0) GOTO 1120
  491. c***c*  Cas deja vu car correspond idejvu(ip) = 2 !
  492. c***          do jl = 1, NBNN1
  493. c***            if (ip.eq.meleme.num(jl,iem)) goto 1120
  494. c***          enddo
  495.  
  496. C         Verification de la boite englobante (Test rapide avant QSIJS)
  497.           IREFP = (IP-1)*IDIMP1    
  498.           DO ii=1,idim
  499.             XCO=XCOOR(IREFP+ii)
  500.             IF(XCO .LT. Xboite(1,ii))GOTO 1120
  501.             IF(XCO .GT. Xboite(2,ii))GOTO 1120
  502.             xpo(ii)=XCO
  503.           ENDDO
  504.  
  505. C         coordonnes intrinseques et fonctions de forme : qsijs
  506. C          WRITE(6,5777 ) ((mtrav.xel(ia,ib),ia=1,idim),ib=1,nbnn1)
  507. C 5777      format(6e12.5)
  508.  
  509.           CALL QSIJS(mtrav.xel,IELE,NBNN1,idim,xpo(1),
  510.      &               mtrav.SHPP,mtrav.qsi,iretou)
  511.           IF (iretou.NE.0) goto 1120
  512.           DO i = 1, inomi1
  513.             ilp = inomil(i)
  514.             IF (mtrav.shpp(1,ilp).lt.CCRIT) goto 1120
  515.           ENDDO
  516. C         ce point exterieur est rattrapable
  517.           krat = krat + 1
  518.           ioubli.iptou(ipt)   =-ip
  519.           NNP                 = icount.NBPTOT + 1
  520.           icount.IEINT(1,NNP) = ip
  521.           icount.IEINT(2,NNP) = iem
  522.           icount.IEINT(3,NNP) = iob
  523.           DO i = 1, idim
  524.             icount.QQQ(i,NNP) = mtrav.qsi(i)
  525.           ENDDO
  526.           icount.MAXPZ(iob) = icount.MAXPZ(iob) + 1
  527.           icount.NBPTOT = NNP
  528. C         WRITE(6,2375)ip,iem,(xpo(i1),i1=1,idim),(qsi(i2),i2=1,idim)
  529. C 2375    format(2i4,6e12.5)
  530. 1120    continue
  531. C.......fin de la boucle sur les points candidats au rattrapage ......
  532. 1121    CONTINUE
  533. C...... fin de la boucle sur les elements enveloppe de la sz .......
  534.       ENDDO
  535. C---- fin de la boucle sur les zones ----
  536.  
  537. C     comptage du nombre de points rattrapes
  538.       k = 0
  539.       do i = 1, IUO
  540.         IF (ioubli.IPTOU(I).LT.0) k = k+1
  541.       ENDDO
  542.       if (k.ne.krat) then
  543.       endif
  544.  
  545. *     WRITE(6,*) ' nombre de points rattrapes ' ,k
  546. C     on donne le nombre de points non projetés
  547. C     if (k.ne.iuo) then
  548. C       INTERR(1)=iuo-k
  549. C       CALL ERREUR(-322)
  550. C     endif
  551.       IACCRO = IACCRO + k
  552.  
  553. C---- FIN DE LA SEQUENCE DE RATTRAPAGE
  554. C-----------------------------------------------------------
  555.  9000 CONTINUE
  556.  
  557. c MESSAGE : Nombre de points accrochés %i1  sur  %i2 proposés
  558.       INTERR(1) = IACCRO
  559.       INTERR(2) = NODES
  560.       CALL ERREUR(-319)
  561. C
  562.       IPT1 = MELEMI
  563. C  NOMBRE DE SOUS ZONES TOUCHÉES  ET NOMBRE TOTAL DE NOEUDS CONCERNÉS
  564.       NBSZR = 0
  565.       ktot = 0
  566.       DO iob = 1, NBSZ
  567.         k = icount.MAXPZ(iob)
  568.         ktot = ktot + k
  569.         if (k.GT.0) NBSZR = NBSZR + 1
  570.       ENDDO
  571.       if (ktot.NE.icount.NBPTOT) then
  572.         write(ioimp,*) 'pb ktot nbptot',ktot,icount.NBPTOT
  573.         CALL ERREUR(5)
  574.       endif
  575.  
  576.       SEGACT,MCOORD*MOD
  577. C     dans le cas IFROCA ne 0 les multiplicateurs existent deja
  578.       IF (IFROCA.EQ.0) THEN
  579.         IF (igliss.EQ.1) NINC = NINC - 1
  580.         NBPTS = NBPTI + (icount.NBPTOT * NINC)
  581.         SEGADJ,MCOORD
  582.       ENDIF
  583.  
  584. C     INITIALISATION  DE MRIGID
  585.       NRIGEL = NBSZR * NINC
  586. C      NRIGE = 8
  587.       SEGINI,MRIGID
  588.       mrigid.IFORIG = IFOUR
  589.       mrigid.MTYMAT ='RIGIDITE'
  590.  
  591.       ipt1 = MELEMI
  592.       meleme = ipt1
  593.       ICZ = 0
  594.  
  595. C________________________________________________________________
  596. C____ BOUCLE SUR LES SOUS ZONES __________________________
  597.       DO 400 IOB = 1, NBSZ
  598.  
  599.         NBELEM = icount.MAXPZ(iob)
  600. *       WRITE(ioimp,*) 'SOUS ZONE ',IOB ,'NOMBRE DE POINTS ',NBELEM
  601.         IF (NBELEM.EQ.0) GOTO 395
  602.  
  603.         IF (NBSZ1.NE.0) meleme = ipt1.lisous(iob)
  604.         ICZ = ICZ + 1
  605.  
  606.         KEL   = meleme.ITYPEL
  607.         NBMA1 = meleme.NUM(/1)
  608.  
  609.         NHD = NBMA1+1
  610.  
  611. C     Dimensions des maillage (IPT3), xmatri, DESCR
  612.         NBNN   = NBMA1 + 2
  613.         NBSOUS = 0
  614.         NBREF  = 0
  615.  
  616.         NELRIG  = NBELEM
  617.         NLIGRD  = NBNN
  618.         IF (igliss.EQ.1) NLIGRD = 1 + IDIM*(NBMA1+1)
  619.         NLIGRP = NLIGRD
  620.  
  621. *       WRITE(ioimp,*) ' ninc' , ninc
  622.         DO iop = 1, NINC
  623.           SEGINI,IPT3
  624.           IPT3.ITYPEL = 22
  625.           SEGINI,XMATRI
  626.           XMATRI.SYMRE = 0
  627.           IMEL(iop) = IPT3
  628.           IMTT(iop) = XMATRI
  629.         ENDDO
  630.  
  631. C.... boucle sur les elements .......
  632.         xpo(3)= XZERO
  633.         ino = 0
  634.         DO 300 K = 1, NODES
  635.           if (icount.ieint(3,k).NE.iob) goto 300
  636.           ino = ino + 1
  637.           IP  = icount.IEINT(1,K)
  638.           IEL = icount.IEINT(2,K)
  639.           irefp = (ip-1) * idimp1
  640.           xpo(1) = xcoor(irefp+1)
  641.           xpo(2) = xcoor(irefp+2)
  642.           IF (IDIM.EQ.3) xpo(3) = xcoor(irefp+3)
  643.           qsi(1) = icount.QQQ(1,k)
  644.           qsi(2) = icount.QQQ(2,k)
  645.           qsi(3) = icount.QQQ(3,k)
  646.  
  647.           CALL SHAPE(QSI(1),QSI(2),QSI(3),KEL,SHPP,iret)
  648.  
  649.           CCOE(1) = 1.D0
  650.           NHD = NBMA1 + 1
  651.           DO KM = 1, NBMA1
  652.             CCOE(1+KM) = -SHPP(1,KM)
  653.           ENDDO
  654.  
  655.           DO 59 IK = 1, NINC
  656.  
  657.             IF (IFROCA.EQ.0) THEN
  658.               IDE = NBPTI * IDIMP1
  659.               XCOOR(IDE+1) = xpo(1)
  660.               XCOOR(IDE+2) = xpo(2)
  661.               XCOOR(IDE+3) = xpo(3)
  662.               NBPTI = NBPTI + 1
  663.             ELSE
  664.               ibp = icpr(ip)
  665.               NBPTI = imul(ik,ibp)
  666.             ENDIF
  667. C    REMPLISSAGE DU MELEME IPT3
  668.             ipt3 = imel(ik)
  669.             ipt3.NUM(1,ino) = NBPTI
  670.             ipt3.NUM(2,ino) = IP
  671.             DO i = 1, NBMA1
  672.               ipt3.NUM(2+i,ino) =  meleme.NUM(i,IEL)
  673.             ENDDO
  674.  
  675. C    REMPLISSAGE DU XMATRI
  676.           XMATRI = imtt(ik)
  677.  
  678. c    Cas glissant
  679.         IF (igliss.eq.1) THEN
  680.          xaa = xperp1(1,ip)
  681.          xbb = xperp1(2,ip)
  682.          if (idim.eq.3) xcc = xperp1(3,ip)
  683.          if (ik.eq.2) then
  684.            xaa = xperp2(1,ip)
  685.            xbb = xperp2(2,ip)
  686.            if (idim.eq.3) xcc = xperp2(3,ip)
  687.          endif
  688.          IF (IFROCA.NE.0.AND.ik.EQ.ninc) then
  689. C petit calcul du vecteur tangent
  690.            if( idim.eq.2) then
  691.              XAA = xperp1(2,ip)
  692.              XBB =-xperp1(1,ip)
  693.            else
  694.              XAA = xperp1(2,ip)*xperp2(3,ip)-xperp1(3,ip)*xperp2(2,ip)
  695.              XBB = xperp1(3,ip)*xperp2(1,ip)-xperp1(1,ip)*xperp2(3,ip)
  696.              XCC = xperp1(1,ip)*xperp2(2,ip)-xperp1(2,ip)*xperp2(1,ip)
  697.            endif
  698.          endif
  699.          ikm = 1
  700.          DO 337 km=1,nhd
  701.            CCC = -CCOE(KM)
  702.            ikm = ikm + 1
  703.  
  704.            RE(1,ikm,INO) = -CCC*XAA
  705.            RE(ikm,1,INO) = -CCC*XAA
  706.            ikm = ikm + 1
  707.            RE(1,ikm,INO) = -CCC*XBB
  708.            RE(ikm,1,INO) = -CCC*XBB
  709.            IF (IDIM.EQ.3) THEN
  710.              ikm = ikm + 1
  711.              RE(1,ikm,INO) = -CCC*XCC
  712.              RE(ikm,1,INO) = -CCC*XCC
  713.            ENDIF
  714.   337    CONTINUE
  715.  
  716. c       Cas non glissant
  717.         ELSE
  718.           DO KM = 1, NHD
  719.             CCC = CCOE(KM)
  720.             RE(1,KM+1,INO) = CCC
  721.             RE(1+KM,1,INO) = CCC
  722.           ENDDO
  723.         ENDIF
  724.   59  continue
  725. c     - fin de boucle sur les xmatri -
  726. 300     CONTINUE
  727. C.... fin de boucle sur les elements .......
  728.         if (ino.ne.nbelem) then
  729.           write(ioimp,*) 'Pb element zone ',iob,ino,nbelem
  730.           call erreur(5)
  731.         endif
  732.  
  733.         INRI = NINC*(ICZ-1) + 1
  734.  
  735.         DO iop = 1, NINC
  736.  
  737.           SEGINI DESCR
  738.           descr.LISINC(1) = 'LX'
  739.           descr.LISDUA(1) = 'FLX'
  740.           descr.NOELEP(1) = 1
  741.           descr.NOELED(1) = 1
  742. c    - Cas non glissant -
  743.           IF (igliss.eq.0) then
  744.             DO k = 2, NBNN
  745.               descr.LISINC(k) = NOMDD(ICOR(iop))
  746.               descr.LISDUA(k) = NOMDU(ICOR(iop))
  747.               descr.NOELEP(k) = k
  748.               descr.NOELED(k) = k
  749.             ENDDO
  750. c    - Cas glissant -
  751.           ELSE
  752.             DO km1 = 1, nhd
  753.               j = 1+(km1-1)*idim
  754.               descr.lisinc(j+1)= NOMDD(ICOR(1))
  755.               descr.lisinc(j+2)= NOMDD(ICOR(2))
  756.               descr.LISDUA(j+1)= NOMDU(ICOR(1))
  757.               descr.LISDUA(j+2)= NOMDU(ICOR(2))
  758.               descr.NOELEP(j+1)= KM1+1
  759.               descr.NOELEP(j+2)= KM1+1
  760.               descr.NOELED(j+1)= KM1+1
  761.               descr.NOELED(j+2)= KM1+1
  762.               if (idim.eq.3) then
  763.                 descr.lisinc(j +3)= NOMDD(ICOR(3))
  764.                 descr.LISDUA(j +3)= NOMDU(ICOR(3))
  765.                 descr.NOELEP(j +3)= KM1+1
  766.                 descr.NOELED(j +3)= KM1+1
  767.               endif
  768.             ENDDO
  769.           endif
  770.           SEGACT,DESCR
  771.           IPT3 = IMEL(iop)
  772.           SEGACT,IPT3
  773.           XMATRI = IMTT(iop)
  774.           SEGACT,xMATRI
  775.  
  776.           jop = INRI + iop - 1
  777.           IRIGEL(1,jop) = IPT3
  778. c          IRIGEL(2,jop) = 0
  779.           IRIGEL(3,jop) = DESCR
  780.           IRIGEL(4,jop) = XMATRI
  781.           IRIGEL(5,jop) = NIFOUR
  782. c          IRIGEL(6,jop) = 0
  783.           IF (IFROCA.NE.0.AND.iop.EQ.ninc) IRIGEL(6,jop) = 2
  784. c          IRIGEL(7,jop) = 0
  785. c          IRIGEL(8,jop) = 0
  786.           COERIG(jop) = 1.D0
  787.         ENDDO
  788.  
  789. 395     CONTINUE
  790.  
  791. 400   CONTINUE
  792. C____ FIN DE BOUCLE SUR LES SOUS ZONES __________________________
  793.  
  794.       if (icz.ne.nbszr) then
  795.         write(ioimp,*) 'Pb Rigi ',ICZ,NBSZR
  796.         call erreur(5)
  797.       endif
  798.  
  799. C**** MENAGE ***********************************************************
  800.       IF (IFROCA.EQ.0) SEGDES,MCOORD
  801.       SEGSUP,icount,iexclu,mtrav
  802.  
  803. *  elimination des termes trop petits dans les relations
  804.       call relasi(mrigid)
  805.       SEGACT,MRIGID
  806.  
  807. C**** ECRITURE DE LA RIGIDTE RESULTAT **********************************
  808.       CALL ECROBJ('RIGIDITE',MRIGID)
  809.  
  810.       END
  811.  
  812.  
  813.  
  814.  
  815.  
  816.  
  817.  
  818.  

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