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

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