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  24/07/25    21:15:02     11950          
  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.  
  493.           CALL QSIJS(mtrav.xel,IELE,NBNN1,IDIM,xpo(1),
  494.      &               mtrav.SHPP,mtrav.qsi,iretou)
  495.           IF (iretou.NE.0) goto 1120
  496.           DO i = 1, inomi1
  497.             ilp = inomil(i)
  498.             IF (mtrav.shpp(1,ilp).lt.CCRIT) goto 1120
  499.           ENDDO
  500. C         ce point exterieur est rattrapable
  501.           krat = krat + 1
  502.           ioubli.iptou(ipt)   =-ip
  503.           NNP                 = icount.NBPTOT + 1
  504.           icount.IEINT(1,NNP) = ip
  505.           icount.IEINT(2,NNP) = iem
  506.           icount.IEINT(3,NNP) = iob
  507.           DO i = 1, IDIM
  508.             icount.QQQ(i,NNP) = mtrav.qsi(i)
  509.           ENDDO
  510.           icount.MAXPZ(iob) = icount.MAXPZ(iob) + 1
  511.           icount.NBPTOT = NNP
  512. C         WRITE(6,2375)ip,iem,(xpo(i1),i1=1,idim),(qsi(i2),i2=1,idim)
  513. C 2375    format(2i4,6e12.5)
  514. 1120    continue
  515. C.......fin de la boucle sur les points candidats au rattrapage ......
  516. 1121    CONTINUE
  517. C...... fin de la boucle sur les elements enveloppe de la sz .......
  518.       ENDDO
  519. C---- fin de la boucle sur les zones ----
  520.  
  521. C     comptage du nombre de points rattrapes
  522.       k = 0
  523.       do i = 1, IUO
  524.         IF (ioubli.IPTOU(I).LT.0) k = k+1
  525.       ENDDO
  526.       if (k.ne.krat) then
  527.       endif
  528.  
  529. *     WRITE(6,*) ' nombre de points rattrapes ' ,k
  530. C     on donne le nombre de points non projetés
  531. C     if (k.ne.iuo) then
  532. C       INTERR(1)=iuo-k
  533. C       CALL ERREUR(-322)
  534. C     endif
  535.       IACCRO = IACCRO + k
  536.  
  537. C---- FIN DE LA SEQUENCE DE RATTRAPAGE
  538. C-----------------------------------------------------------
  539.  9000 CONTINUE
  540.  
  541. c MESSAGE : Nombre de points accrochés %i1  sur  %i2 proposés
  542.       INTERR(1) = IACCRO
  543.       INTERR(2) = NODES
  544.       CALL ERREUR(-319)
  545.  
  546.       IPT1 = MELEMI
  547. C  NOMBRE DE SOUS-ZONES TOUCHEES ET NOMBRE TOTAL DE NOEUDS CONCERNES
  548.       NBSZR = 0
  549.       ktot = 0
  550.       DO iob = 1, NBSZ
  551.         k = icount.MAXPZ(iob)
  552.         ktot = ktot + k
  553.         if (k.GT.0) NBSZR = NBSZR + 1
  554.       ENDDO
  555.       if (ktot.NE.icount.NBPTOT) then
  556.         write(ioimp,*) 'pb ktot nbptot',ktot,icount.NBPTOT
  557.         CALL ERREUR(5)
  558.       endif
  559.  
  560.       SEGACT,MCOORD*MOD
  561. C     dans le cas IFROCA ne 0 les multiplicateurs existent deja
  562.       IF (IFROCA.EQ.0) THEN
  563.         IF (igliss.EQ.1) NINC = NINC - 1
  564.         NBPTS = NBPTI + (icount.NBPTOT * NINC)
  565.         SEGADJ,MCOORD
  566.       ENDIF
  567.  
  568. C     INITIALISATION  DE MRIGID
  569.       NRIGEL = NBSZR * NINC
  570. C      NRIGE = 8
  571.       SEGINI,MRIGID
  572.       mrigid.IFORIG = IFOUR
  573.       mrigid.MTYMAT ='RIGIDITE'
  574.  
  575.       ipt1 = MELEMI
  576.       meleme = ipt1
  577.       ICZ = 0
  578.  
  579. C________________________________________________________________
  580. C____ BOUCLE SUR LES SOUS ZONES __________________________
  581.       DO 400 IOB = 1, NBSZ
  582.  
  583.         NBELEM = icount.MAXPZ(iob)
  584. c*        WRITE(ioimp,*) 'SOUS ZONE ',IOB ,'NOMBRE DE POINTS ',NBELEM
  585.         IF (NBELEM.EQ.0) GOTO 395
  586.  
  587.         IF (NBSZ1.NE.0) meleme = ipt1.lisous(iob)
  588.         ICZ = ICZ + 1
  589.  
  590.         KEL   = meleme.ITYPEL
  591.         NBMA1 = meleme.NUM(/1)
  592.  
  593. * Cas particulier des POLYGONES : meleme.itypel = 32 et mele = 108 + nbnoe
  594.         IF (KEL .EQ. 32) THEN
  595.           MELE = 108 + NBMA1
  596.         ENDIF
  597.  
  598.         NHD = NBMA1+1
  599.  
  600. C     Dimensions des maillage (IPT3), xmatri, DESCR
  601.         NBNN   = NBMA1 + 2
  602.         NBSOUS = 0
  603.         NBREF  = 0
  604.  
  605.         NELRIG  = NBELEM
  606.         NLIGRD  = NBNN
  607.         IF (igliss.EQ.1) NLIGRD = 1 + IDIM*(NBMA1+1)
  608.         NLIGRP = NLIGRD
  609.  
  610. c*       WRITE(ioimp,*) ' ninc' , ninc
  611.         DO iop = 1, NINC
  612.           SEGINI,IPT3
  613.           IPT3.ITYPEL = 22
  614.           SEGINI,XMATRI
  615.           XMATRI.SYMRE = 0
  616.           IMEL(iop) = IPT3
  617.           IMTT(iop) = XMATRI
  618.         ENDDO
  619.  
  620. C.... boucle sur les elements .......
  621.         xpo(3)= XZERO
  622.         ino = 0
  623.         DO 300 K = 1, NODES
  624.           if (icount.ieint(3,k).NE.iob) goto 300
  625.           ino = ino + 1
  626.           IP  = icount.IEINT(1,K)
  627.           IEL = icount.IEINT(2,K)
  628.           irefp = (ip-1) * IDIMP1
  629.           xpo(1) = xcoor(irefp+1)
  630.           xpo(2) = xcoor(irefp+2)
  631.           IF (IDIM.EQ.3) xpo(3) = xcoor(irefp+3)
  632.           qsi(1) = icount.QQQ(1,k)
  633.           qsi(2) = icount.QQQ(2,k)
  634.           qsi(3) = icount.QQQ(3,k)
  635.  
  636. * Cas particulier des POLYGONES :
  637.           IF (KEL .EQ. 32) THEN
  638.             CALL SHPOLY(QSI(1),QSI(2),QSI(3),MELE,SHPP,iret)
  639.           ELSE
  640.             CALL SHAPE(QSI(1),QSI(2),QSI(3),KEL,SHPP,iret)
  641.           ENDIF
  642.           if (iret.eq.0) then
  643.             moterr(1:4)=NOMS(KEL)
  644.             call erreur(68)
  645.             return
  646.           endif
  647.  
  648.           CCOE(1) = 1.D0
  649.           NHD = NBMA1 + 1
  650.           DO KM = 1, NBMA1
  651.             CCOE(1+KM) = -SHPP(1,KM)
  652.           ENDDO
  653.  
  654.           DO 59 IK = 1, NINC
  655.  
  656.             IF (IFROCA.EQ.0) THEN
  657.               IDE = NBPTI * IDIMP1
  658.               XCOOR(IDE+1) = xpo(1)
  659.               XCOOR(IDE+2) = xpo(2)
  660.               XCOOR(IDE+3) = xpo(3)
  661.               NBPTI = NBPTI + 1
  662.             ELSE
  663.               ibp = icpr(ip)
  664.               NBPTI = imul(ik,ibp)
  665.             ENDIF
  666. C    REMPLISSAGE DU MELEME IPT3
  667.             ipt3 = imel(ik)
  668.             ipt3.NUM(1,ino) = NBPTI
  669.             ipt3.NUM(2,ino) = IP
  670.             DO i = 1, NBMA1
  671.               ipt3.NUM(2+i,ino) =  meleme.NUM(i,IEL)
  672.             ENDDO
  673.  
  674. C    REMPLISSAGE DU XMATRI
  675.           XMATRI = imtt(ik)
  676.  
  677. c    Cas glissant
  678.         IF (igliss.eq.1) THEN
  679.          xaa = xperp1(1,ip)
  680.          xbb = xperp1(2,ip)
  681.          if (IDIM.eq.3) xcc = xperp1(3,ip)
  682.          if (ik.eq.2) then
  683.            xaa = xperp2(1,ip)
  684.            xbb = xperp2(2,ip)
  685.            if (IDIM.eq.3) xcc = xperp2(3,ip)
  686.          endif
  687.          IF (IFROCA.NE.0.AND.ik.EQ.ninc) then
  688. C petit calcul du vecteur tangent
  689.            if( IDIM.eq.2) then
  690.              XAA = xperp1(2,ip)
  691.              XBB =-xperp1(1,ip)
  692.            else
  693.              XAA = xperp1(2,ip)*xperp2(3,ip)-xperp1(3,ip)*xperp2(2,ip)
  694.              XBB = xperp1(3,ip)*xperp2(1,ip)-xperp1(1,ip)*xperp2(3,ip)
  695.              XCC = xperp1(1,ip)*xperp2(2,ip)-xperp1(2,ip)*xperp2(1,ip)
  696.            endif
  697.          endif
  698.          ikm = 1
  699.          DO km = 1, nhd
  700.            CCC = -CCOE(KM)
  701.            ikm = ikm + 1
  702.  
  703.            RE(1,ikm,INO) = -CCC*XAA
  704.            RE(ikm,1,INO) = -CCC*XAA
  705.            ikm = ikm + 1
  706.            RE(1,ikm,INO) = -CCC*XBB
  707.            RE(ikm,1,INO) = -CCC*XBB
  708.            IF (IDIM.EQ.3) THEN
  709.              ikm = ikm + 1
  710.              RE(1,ikm,INO) = -CCC*XCC
  711.              RE(ikm,1,INO) = -CCC*XCC
  712.            ENDIF
  713.          ENDDO
  714.  
  715. c       Cas non glissant
  716.         ELSE
  717.           DO KM = 1, NHD
  718.             CCC = CCOE(KM)
  719.             RE(1,KM+1,INO) = CCC
  720.             RE(1+KM,1,INO) = CCC
  721.           ENDDO
  722.         ENDIF
  723.   59  continue
  724. c     - fin de boucle sur les xmatri -
  725. 300     CONTINUE
  726. C.... fin de boucle sur les elements .......
  727.         if (ino.ne.nbelem) then
  728.           write(ioimp,*) 'Pb element zone ',iob,ino,nbelem
  729.           call erreur(5)
  730.         endif
  731.  
  732.         INRI = NINC*(ICZ-1) + 1
  733.  
  734.         DO iop = 1, NINC
  735.  
  736.           SEGINI DESCR
  737.           descr.LISINC(1) = 'LX  '
  738.           descr.LISDUA(1) = 'FLX '
  739.           descr.NOELEP(1) = 1
  740.           descr.NOELED(1) = 1
  741. c    - Cas non glissant -
  742.           IF (igliss.eq.0) then
  743.             DO k = 2, NBNN
  744.               descr.LISINC(k) = NOMDD(ICOR(iop))
  745.               descr.LISDUA(k) = NOMDU(ICOR(iop))
  746.               descr.NOELEP(k) = k
  747.               descr.NOELED(k) = k
  748.             ENDDO
  749. c    - Cas glissant -
  750.           ELSE
  751.             DO km1 = 1, nhd
  752.               j = 1+(km1-1)*IDIM
  753.               descr.lisinc(j+1)= NOMDD(ICOR(1))
  754.               descr.lisinc(j+2)= NOMDD(ICOR(2))
  755.               descr.LISDUA(j+1)= NOMDU(ICOR(1))
  756.               descr.LISDUA(j+2)= NOMDU(ICOR(2))
  757.               descr.NOELEP(j+1)= KM1+1
  758.               descr.NOELEP(j+2)= KM1+1
  759.               descr.NOELED(j+1)= KM1+1
  760.               descr.NOELED(j+2)= KM1+1
  761.               if (IDIM.eq.3) then
  762.                 descr.lisinc(j +3)= NOMDD(ICOR(3))
  763.                 descr.LISDUA(j +3)= NOMDU(ICOR(3))
  764.                 descr.NOELEP(j +3)= KM1+1
  765.                 descr.NOELED(j +3)= KM1+1
  766.               endif
  767.             ENDDO
  768.           endif
  769.           SEGACT,DESCR
  770.           IPT3 = IMEL(iop)
  771.           SEGACT,IPT3
  772.           XMATRI = IMTT(iop)
  773.           SEGACT,xMATRI
  774.  
  775.           jop = INRI + iop - 1
  776.           IRIGEL(1,jop) = IPT3
  777. c          IRIGEL(2,jop) = 0
  778.           IRIGEL(3,jop) = DESCR
  779.           IRIGEL(4,jop) = XMATRI
  780.           IRIGEL(5,jop) = NIFOUR
  781. c          IRIGEL(6,jop) = 0
  782.           IF (IFROCA.NE.0.AND.iop.EQ.ninc) IRIGEL(6,jop) = 2
  783. c          IRIGEL(7,jop) = 0
  784. c          IRIGEL(8,jop) = 0
  785.           COERIG(jop) = 1.D0
  786.         ENDDO
  787.  
  788. 395     CONTINUE
  789.  
  790. 400   CONTINUE
  791. C____ FIN DE BOUCLE SUR LES SOUS ZONES __________________________
  792.  
  793.       if (icz.ne.nbszr) then
  794.         write(ioimp,*) 'Pb Rigi ',ICZ,NBSZR
  795.         call erreur(5)
  796.       endif
  797.  
  798. C**** MENAGE ***********************************************************
  799.       IF (IFROCA.EQ.0) SEGDES,MCOORD
  800.       SEGSUP,icount,iexclu,mtrav
  801.  
  802. *  elimination des termes trop petits dans les relations
  803.       call relasi(mrigid)
  804.       SEGACT,MRIGID
  805.  
  806. C**** ECRITURE DE LA RIGIDTE RESULTAT **********************************
  807.       CALL ECROBJ('RIGIDITE',MRIGID)
  808.  
  809. c      return
  810.       END
  811.  
  812.  
  813.  

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