Télécharger la plaquette Cast3M


fr en

Télécharger pjmode.eso

Retour à la liste

Numérotation des lignes : 

pjmode

  1. C PJMODE    SOURCE    PV090527  26/04/30    21:15:59     12529          
  2.       SUBROUTINE PJMODE(ipmode)
  3.       IMPLICIT INTEGER(I-N)
  4.       IMPLICIT REAL*8 (A-H,O-Z)
  5. C=======================================================================
  6. C OPERATEUR PJBA :
  7. C           PROJECTION D'UN CHPOINT, D'UN CHARGEMENT OU D'UNE RIGIDITE
  8. C           SUR LES ELEMENTS D'UNE BASE MODALE B.
  9. C           LE RESULTAT EST DU MEME TYPE.
  10. C
  11. C SYNTAXE :
  12. C       *   FN  =  PJBA  B   OBJET           ;   SI BASE ELEMENTAIRE
  13. C       *   FN  =  PJBA  B  STR1  (N)  OBJET ;   SI BASE COMPLEXE
  14. C
  15. C         OBJET    POUVANT ETRE UNE FORCE OU UN CHARGEMENT,
  16. C                  OU UNE RIGIDITE DANS LE PREMIER CAS.
  17. C=======================================================================
  18. ***********************************************************
  19. *     PROJECTION D'UNE MATRICE SUR UNE BASE DE MODES      *
  20. * _______________________________________________________ *
  21. *                                                         *
  22. *     DATE : le 11 Avril 1995                             *
  23. *     AUTEUR : Nicolas BENECH                             *
  24. * _______________________________________________________ *
  25. *                                                         *
  26. *     MODULE(S) APPELANT(S) : PJBA                        *
  27. *                                                         *
  28. *     MODULE(S) APPELE(S) : ACCTAB, YTMX                  *
  29. * _______________________________________________________ *
  30. *                                                         *
  31. *    EN ENTREE :                                          *
  32. *        MRIGID : Matrice a projeter                      *
  33. *        MTAB1  : Base de modes, reels ou complexes       *
  34. *        'REEL' : indique que l'on utilise le produit     *
  35. *                 scalaire reel (pas de conjugaison)      *
  36. *                                                         *
  37. *    EN SORTIE :                                          *
  38. *        RI1 : Matrice projetee (partie reelle)           *
  39. *        RI2 : Matrice projetee (partie imaginaire)       *
  40. * _______________________________________________________ *
  41. *                                                         *
  42. *    REMARQUE :                                           *
  43. *         L'operation realisee est :                      *
  44. *                (MTAB1)t . MRIGID. MTAB1                 *
  45. *  Dans le cas complexe, la transposition est accompagnee *
  46. *   de la conjugaison (si REEL n'est pas mentionne).      *
  47. *
  48. * voir aussi PROJTA
  49. ***********************************************************
  50. *
  51. -INC SMCHPOI
  52. -INC SMCHARG
  53. -INC SMLCHPO
  54. -INC PPARAM
  55. -INC CCOPTIO
  56. -INC CCGEOME
  57. -INC CCREEL
  58. -INC SMELEME
  59. -INC CCHAMP
  60. -INC SMCHAML
  61. -INC SMMODEL
  62. -INC SMRIGID
  63. -INC SMCOORD
  64. -INC SMLMOTS
  65. -INC SMLENTI
  66.  
  67. C
  68. * Declarations
  69. *
  70.       PARAMETER(ZERO=0.D0)
  71.       REAL*8 XVAL, RMAX
  72.       CHARACTER*8 LETYPE
  73.       CHARACTER*8  TYPMOD,cmate
  74.       LOGICAL MODCOM,dedans,dchpo,l3,lr2,lirl
  75.       INTEGER I, J, NBMOD, POS, IREEL, IVALRE, IOBRE
  76.       REAL*8 XVALRE
  77.       LOGICAL LOGRE
  78.       segment plcf
  79.        integer lpref(ldepl),ldefo(ldepl),lmade(ldepl)
  80.        real*8  prmas(ldepl)
  81.       endsegment
  82.       segment pmod
  83.         integer kdefo(nbmod)
  84.       endsegment
  85.       segment prigmat
  86.         integer lrigmat(nrigmat,2+9)
  87.       endsegment
  88.       segment pmapmo
  89.          integer lmapmo(nmapmo),defpmo(nmapmo),dimpmo(nmapmo)
  90.          character*(LOCHPO) compmo(nmapmo)
  91.          real*8 coepmo(nmapmo)
  92.       endsegment
  93.       segment pcompo
  94.          character*4 mcol
  95.          real*8  valmod(nipmod)
  96.       endsegment
  97.       LOGICAL L0,L1,lcf
  98.       PARAMETER (ncod=8)
  99.       CHARACTER*(lochpo) IDDL,lcod(ncod),lcof(ncod),motinc
  100.       CHARACTER*8 TYPRET,CHARRE
  101.       data xlopre/1.d-11/
  102.       DATA KZERO/0/
  103.       data lcod/'UX','UY','UZ','RX','RY','RZ','UR','UT'/
  104.       data lcof/'FX','FY','FZ','MX','MY','MZ','FR','FT'/
  105.  
  106.       plcf = 0
  107.       jgn = lochpo
  108.       jgm = ncod
  109.       segini mlmot5
  110.       segini mlmot6
  111.        do io = 1,ncod
  112.         mlmot5.mots(io) = lcod(io)
  113.         mlmot6.mots(io) = lcof(io)
  114.        enddo
  115.  
  116.       modcom = .false.
  117.       dchpo = .false.
  118.       iriout = 0
  119.       iriout1 = 0
  120.       iriout2 = 0
  121.       mmodel = ipmode
  122.       n1 = kmodel(/1)
  123.       segini mmode1
  124.       jn = 0
  125.       do im = 1, n1
  126.        imodel = kmodel(im)
  127.        if (formod(1).eq.'MECANIQUE'.and.MATMOD(1).eq.'ELASTIQUE'
  128.      &.and.(MATMOD(2).eq.'MODAL'.or.MATMOD(2).eq.'STATIQUE')) then
  129.         jn = jn + 1
  130.         mmode1.kmodel(jn) = imodel
  131.        endif
  132.       enddo
  133.       if (jn.ne.0) then
  134.        n1 = jn
  135.        segadj mmode1
  136.        ipmode = mmode1
  137.       else
  138.        segsup mmode1
  139. * cas de projection non pr�vue
  140.        call erreur(5)
  141.        return
  142.       endif
  143.  
  144.       call lirobj('MCHAML  ',IPCAR1,1,iretou)
  145.       call actobj('MCHAML  ',IPCAR1,1)
  146.       if (ierr.ne.0) return
  147.  
  148.       ipchpo = 0
  149.       iprigi = 0
  150.       call lirobj('CHARGEME',IPCHAR,0,iretou)
  151.       if (iretou.eq.0) then
  152.         call lirobj('CHPOINT ',IPCHPO,0,iretou)
  153.         if(iretou .EQ. 1)call actobj('CHPOINT ',IPCHPO,1)
  154.       endif
  155.       if (iretou.eq.0) call lirobj('RIGIDITE',IPRIGI,0,iretou)
  156.  
  157.       if (iretou.eq.0) then
  158. * manque un op�rande
  159.         call erreur(5)
  160.         return
  161.       endif
  162.  
  163.       call reduaf (ipcar1,ipmode,IPCARA,1,iretr,kerre)
  164.       if (ierr.ne.0) return
  165.       if( iretr.ne.1) then
  166.         call erreur (kerre)
  167.         return
  168.       endif
  169.  
  170.        lcf = .false.
  171.        mmodel = ipmode
  172.        mchelm = ipcara
  173.       if (ipchar.ne.0) goto 100
  174.       if (iprigi.ne.0) goto 200
  175.       if (ipchpo.ne.0) then
  176.         n = 1
  177.         segini mcharg
  178.         ipchar = mcharg
  179.         segini icharg
  180.         kcharg(1) = icharg
  181.         ichpo1 = ipchpo
  182.         goto 100
  183.       endif
  184.  
  185.  
  186.  100   continue
  187.           MCHAR1=IPCHAR
  188.           SEGINI,MCHARG=MCHAR1
  189.           NBCHG=KCHARG(/1)
  190.           DO 10 INCHA=1,NBCHG
  191.             ICHAR1=KCHARG(INCHA)
  192.             SEGINI,ICHARG=ICHAR1
  193.             KCHARG(INCHA)=ICHARG
  194.             IP1=ICHPO1
  195. c
  196.       IRET = 0
  197. c
  198. c   deplacement impose => idepi=1
  199. c   force imposee      => idepi=0
  200. c
  201.       IDEPI =  0
  202. c     idepi = -1
  203.       KDEPI =  0
  204.       MCHPOI = IP1
  205.       CALL ACTOBJ('CHPOINT',IP1,1)
  206.       IF (MTYPOI.EQ.'FLX     ') IDEPI = 1
  207. c     if (idepi.lt.0) then
  208. c        moterr(1:8) = 'chpoint'
  209. c        call erreur(302)
  210. c        return
  211. c     endif
  212. c
  213.       NBNN   = 1
  214.       NBREF  = 0
  215.       NBSOUS = 0
  216. *
  217.       LDEPL = kmodel(/1)
  218.       if (.not.lcf) segini plcf
  219. c
  220. c
  221. c  **** on initialise le chpoint
  222. c
  223.       NSOUPO = 1
  224.       NAT=1
  225.       SEGINI,MCHPOI
  226.       IRET = MCHPOI
  227.       MTYPOI = '        '
  228.       MOCHDE=' J''AI ETE FABRIQUE PAR L''OPERATEUR PJBA'
  229.       IFOPOI = IFOUR
  230. *     champ de force nodal: nature discrete
  231.       JATTRI(1)=2
  232.       NC = 1
  233.       SEGINI,MSOUPO
  234.       IPCHP(1)  = MSOUPO
  235.       NOHARM(1) = NIFOUR
  236.       NOCOMP(1) = 'FALF    '
  237.  
  238.       do 101 inocomp=1,2
  239.  
  240.       N = LDEPL
  241.       SEGINI MPOVAL
  242.       IPOVAL = MPOVAL
  243. *
  244.       NBNN = 1
  245.       NBELEM = LDEPL
  246.       NBSOUS = 0
  247.       NBREF = 0
  248.       SEGINI MELEME
  249.       IGEOC = MELEME
  250.       ITYPEL = 1
  251.  
  252.         knum = 0
  253. c
  254. c  ****boucle sur les chpoints de depl
  255. c
  256.       DO 11 IM = 1,kmodel(/1)
  257.         imodel = kmodel(im)
  258.         nomid = lnomid(2)
  259.       if (.not.lcf) then
  260.         ipt1 = imamod
  261.         iptr = ipt1.num(1,1)
  262.         lpref(im) = iptr
  263.  
  264.       indc = 1
  265.  34   if (imache(indc).ne.imamod.or.conche(indc).ne.conmod) then
  266.           indc = indc + 1
  267.           if (indc.gt.imache(/1)) then
  268. * champ de caracteristiques incomplet
  269.            goto 99
  270.           endif
  271.           goto 34
  272.         endif
  273.  
  274.         mchaml = ichaml(indc)
  275.         do iij = 1, nomche(/2)
  276.           if (nomche(iij).eq.'DEFO') then
  277.             melval = ielval(iij)
  278.             ipp1 = ielche(1,1)
  279.            ldefo(im) = ipp1
  280.           endif
  281.           if (nomche(iij).eq.'MADE') then
  282.             melval = ielval(iij)
  283.             ipp2 = ielche(1,1)
  284.             lmade(im) = ipp2
  285.           endif
  286.           if (nomche(iij).eq.'MASS') then
  287.             melval = ielval(iij)
  288.             ymass = velche(1,1)
  289.             prmas(im) = ymass
  290.           endif
  291.          if(ldefo(im).gt.0.and.lmade(im).gt.0.and.
  292.      &prmas(im).gt.0) goto 35
  293.         enddo
  294.  35     continue
  295.         if (ldefo(im).eq.0) goto 99
  296.          if (prmas(im).le.0.and.cmatee(1:5).eq.'MODAL') goto 99
  297.          if (lmade(im).eq.0.and.cmatee(1:8).eq.'STATIQUE') goto 99
  298.  
  299.       endif
  300.  
  301.          if (NOCOMP(1).ne.lesobl(1)) goto 11
  302.          knum = knum + 1
  303.  
  304.       iptr = lpref(im)
  305.       ipp1 = ldefo(im)
  306.          NUM(1,knum) = IPTR
  307.          ICOLOR(knum) = IDCOUL
  308.          XRET = 0.D0
  309.          call xty1(ipp1,ip1,mlmot5,mlmot6,xret)
  310.         IF (IDEPI.NE.1) THEN
  311.         ELSE
  312. * ??
  313.           indn = 1
  314.  45      if (nomche(indn).ne.'FREQ') then
  315.          indn = indn + 1
  316.           if (indn.gt.nomche(/2)) then
  317. * pas la composante FREQ
  318.            goto 99
  319.           endif
  320.          goto 45
  321.          endif
  322.  
  323.          melval = ielval(indn)
  324.          x1 = velche(1,1)
  325.              OM = X1
  326.              OM = 2.D0 * XPI * OM
  327.              OM = OM * OM
  328.              XRET = -XRET / OM
  329.         ENDIF
  330.         IF (IFOUR .EQ. 1) THEN
  331.              IF (NIFOUR .NE. 0) THEN
  332.                 XRET = XRET*XPI
  333.              ELSE
  334.                 XRET = XRET*2.D0*XPI
  335.              ENDIF
  336.         ENDIF
  337.        VPOCHA(knum,1) = XRET
  338. *
  339.        if (cmatee(1:5).eq.'MODAL') then
  340.         ymass = prmas(im)
  341.        elseif (cmatee(1:8).eq.'STATIQUE') then
  342.         ipp2 = lmade(im)
  343.         call xty1(ipp1,ipp2,mlmot5,mlmot6,ymass)
  344.        else
  345.        endif
  346.       if (lmade(im).gt.0.and.ABS(XRET).gt.(1.d-10*ymass).and.
  347.      & ymass.gt.0.and.cmatee(1:5).eq.'MODAL') then
  348. * kich : on enleve la projection sur base modale - a creuser pour statique !
  349.        CALL ADCHPO(IP1,IPP2,IP2,1.d0,(XRET/ymass*(-1.d0)))
  350.        IP1 = IP2
  351.       endif
  352. *
  353.  11   CONTINUE
  354. *
  355.       lcf = .true.
  356. *
  357. *
  358.       if (knum.eq.kmodel(/1)) goto 102
  359.       if (inocomp.eq.1) then
  360.         if (knum.eq.0) then
  361.          NOCOMP(1) = 'FBET    '
  362.         else
  363.          N = knum
  364.          NBELEM = knum
  365.          segadj MPOVAL,MELEME
  366.          NSOUPO = 2
  367.          segadj MCHPOI
  368.          SEGINI,MSOUPO
  369.          IPCHP(2)  = MSOUPO
  370.          NOCOMP(1) = 'FBET    '
  371.         endif
  372.       endif
  373.  101  continue
  374.  
  375.  102  continue
  376.       N = knum
  377.       NBELEM = knum
  378.       segadj MPOVAL,MELEME
  379.  
  380.             IF(IERR.NE.0) RETURN
  381.             ICHPO1=IRET
  382.             SEGDES,ICHARG
  383.  10       CONTINUE
  384.       segsup mlmot5,mlmot6,plcf
  385.         if (ipchpo.gt.0) then
  386.           segsup icharg,mcharg
  387.           call actobj('CHPOINT ',iret,1)
  388.           call ecrobj('CHPOINT ',iret)
  389.           goto 999
  390.         endif
  391.           SEGDES,MCHARG
  392.           CALL ECROBJ('CHARGEME',MCHARG)
  393.  
  394.       goto 999
  395.  99   segsup mpoval,msoupo,mchpoi
  396.             call erreur(26)
  397.             return
  398.  
  399.  
  400.  200  continue
  401.       ipri1 = iprigi
  402.       call SEPA(ipri1,1)
  403.       ipri2 = iprigi
  404.       call SEPA(ipri2,2)
  405. *
  406. *
  407. *
  408. *
  409.       nmapmo = 100
  410.       kpmo = 0
  411.       segini pmapmo
  412.       do isous = 1,kmodel(/1)
  413.        imodel = kmodel(isous)
  414.        cmate = cmatee
  415.        meleme = imamod
  416.        if (itypel.ne.1) call erreur(5)
  417.        if (num(/1).ne.1) call erreur(5)
  418.        if (cmate.eq.'STATIQUE'.or.cmate.EQ.'MODAL') then
  419.          do ilp = 1,num(/2)
  420.            kpmo = kpmo + 1
  421.            if (kpmo.gt.nmapmo) then
  422.             nmapmo = nmapmo + 100
  423.             segadj pmapmo
  424.            endif
  425.            lmapmo(kpmo) = num(1,ilp)
  426.         if (cmate.eq.'STATIQUE') then
  427.            compmo(kpmo) = 'BETA    '
  428.         elseif (cmate.eq.'MODAL') then
  429.            compmo(kpmo) = 'ALFA    '
  430.         endif
  431.           do im = 1 , imache(/1)
  432.            if (imache(im).eq.imamod) then
  433.              if (conche(im).eq.conmod) then
  434.                mchaml = ichaml(im)
  435.                do iv = 1,ielval(/1)
  436.                  if (nomche(iv).eq.'DEFO') then
  437.                    melval = ielval(iv)
  438.                    ibmn = min(ilp,ielche(/2))
  439.                    defpmo(kpmo) = ielche(1,ibmn)
  440.                  endif
  441.                  if (nomche(iv).eq.'IDEF') then
  442.                    melval = ielval(iv)
  443.                    ibmn = min(ilp,ielche(/2))
  444.                    dimpmo(kpmo) = ielche(1,ibmn)
  445.                  endif
  446.                enddo
  447.              endif
  448.            endif
  449.           enddo
  450.  
  451.          enddo
  452.        endif
  453.       enddo
  454.  
  455.       nmapmo = kpmo
  456.       segadj pmapmo
  457.       nbmod = nmapmo
  458. *
  459.       N1 = NBMOD
  460.       nbcod = 8
  461.       segini pmod
  462.       SEGINI, MLCHP1
  463.       SEGINI, MLCHP2
  464.        jgm = 1
  465.        jgn = 4
  466.        segini mlmot4
  467. *
  468. * Constitution du maillage support et du segment descriptif
  469. *
  470.       NBNN = NBMOD
  471.       NBELEM = 1
  472.       NBSOUS = 0
  473.       NBREF = 0
  474.       SEGINI, MELEME
  475.         ITYPEL = 1
  476. *
  477.         NLIGRD=NBMOD
  478.         NLIGRP=NBMOD
  479.         SEGINI, DESCR
  480. *
  481.       mrigid = ipri1
  482.       segact mrigid
  483.       nrigel = coerig(/1)
  484.       if (nrigel.lt.1) goto 250
  485.       typmod = ' '
  486.       IREEL = -1
  487. C* POS ?      IF (POS.EQ.1) IREEL = 1
  488. *
  489.       LETYPE = ' '
  490.       DO 210  IM=1,NBMOD
  491.             IPT1 = 0
  492. *
  493.         imodel = kmodel(im)
  494.         ipt1 = imamod
  495.         iptr = ipt1.num(1,1)
  496. *  Cas reel ou cas complexe ?
  497. *
  498.           if (dimpmo(im).gt.0) TYPMOD = 'MODE_COM'
  499.  
  500.             IF (TYPMOD .EQ. 'MODE_COM') THEN
  501.               MODCOM=.TRUE.
  502.               mchpoi = defpmo(im)
  503.               MLCHP1.ICHPOI(IM) = MCHPOI
  504.               mchpoi = dimpmo(im)
  505.             MLCHP2.ICHPOI(IM) = MCHPOI
  506.             ELSE
  507.               MODCOM = .FALSE.
  508.               mchpoi = defpmo(im)
  509.               MLCHP1.ICHPOI(IM) = MCHPOI
  510.             ENDIF
  511. *
  512.             ipt1 = iptr
  513. *
  514.           MELEME.NUM(IM,1)=IPT1
  515. *
  516.          if (cmatee.eq.'MODAL') then
  517.             DESCR.LISINC(IM) = 'ALFA    '
  518.             DESCR.LISDUA(IM) = 'FALF    '
  519.          else if (cmatee.eq.'STATIQUE') then
  520.             DESCR.LISINC(IM) = 'BETA    '
  521.             DESCR.LISDUA(IM) = 'FBET    '
  522.          endif
  523.             DESCR.NOELEP(IM) = IM
  524.             DESCR.NOELED(IM) = IM
  525. *
  526.  210       CONTINUE
  527. *
  528. * Constitution des segments XMATRI
  529. *
  530.         NLIGRD=NBMOD
  531.         NLIGRP=NBMOD
  532.                 nelrig=1
  533. *
  534.         IF (LETYPE .EQ. 'ANNULE') THEN
  535.           rigrel=0
  536.           SEGINI, XMATR1
  537.           IF (MODCOM) THEN
  538.             rigrel=0
  539.             SEGINI, XMATR2
  540.             SEGDES, XMATR2
  541.           ENDIF
  542.           SEGDES, XMATR1
  543.           GOTO 55
  544.         ENDIF
  545. *
  546. * Cas reel
  547. *
  548.         rigrel=0
  549.         SEGINI, XMATR1
  550.           DO I=1, NBMOD
  551.             MCHPO1 = MLCHP1.ICHPOI(I)
  552.             DO J=1, NBMOD
  553.               MCHPO2 = MLCHP1.ICHPOI(J)
  554.               CALL YTMX (MCHPO2, MCHPO1, MRIGID, XVAL)
  555.               XMATR1.RE(I,J,1)=XVAL
  556.             ENDDO
  557.           ENDDO
  558.         SEGDES, XMATR1
  559. *
  560. * Cas complexe : calcul de termes complementaires
  561. *
  562.         IF (MODCOM) THEN
  563.           SEGACT, XMATR1*mod
  564.             DO I=1, NBMOD
  565.               MCHPO1 = MLCHP2.ICHPOI(I)
  566.               DO J=1, NBMOD
  567.                 MCHPO2 = MLCHP2.ICHPOI(J)
  568.                 CALL YTMX (MCHPO1, MCHPO2, MRIGID, XVAL)
  569.                 XMATR1.RE(I,J,1)=XMATR1.RE(I,J,1)-IREEL*XVAL
  570.               ENDDO
  571.             ENDDO
  572.           SEGDES, XMATR1
  573. *
  574.           rigrel=0
  575.           SEGINI, XMATR2
  576.             DO I=1, NBMOD
  577.               MCHPO1 = MLCHP1.ICHPOI(I)
  578.               DO J=1, NBMOD
  579.                 MCHPO2 = MLCHP2.ICHPOI(J)
  580.                 CALL YTMX (MCHPO1, MCHPO2, MRIGID, XVAL)
  581.                 XMATR2.RE(I,J,1)=XVAL
  582.               ENDDO
  583.             ENDDO
  584.             DO I=1, NBMOD
  585.               MCHPO1 = MLCHP2.ICHPOI(I)
  586.               DO J=1, NBMOD
  587.                 MCHPO2 = MLCHP1.ICHPOI(J)
  588.                 CALL YTMX (MCHPO1, MCHPO2, MRIGID, XVAL)
  589.                 XMATR2.RE(I,J,1)=XMATR2.RE(I,J,1)+IREEL*XVAL
  590.               ENDDO
  591.             ENDDO
  592.           SEGDES, XMATR2
  593.         ENDIF
  594. *
  595.       SEGACT, MRIGID
  596.         LETYPE = MRIGID.MTYMAT
  597.       SEGDES, MRIGID
  598. *
  599. * Creation des segments IMATRI
  600. *
  601.  55   NELRIG = 1
  602. *      SEGINI, IMATR1
  603. *        IMATR1.IMATTT(1) = XMATR1
  604.       SEGDES, xMATR1
  605.       IF (MODCOM) THEN
  606. *        SEGINI, IMATR2
  607. *          IMATR2.IMATTT(1) = XMATR2
  608.         SEGDES, xMATR2
  609.       ENDIF
  610. *
  611. * Creation des rigidites calculees
  612. *
  613.       NRIGE=7
  614.       NRIGEL=1
  615.       SEGINI, RI1
  616.         RI1.MTYMAT = LETYPE
  617.         RI1.IFORIG = IFOUR
  618.         RI1.IMGEO1 = 0
  619.         RI1.IMGEO2 = 0
  620.         RI1.COERIG = 1.D0
  621.         RI1.IRIGEL(1,1) = MELEME
  622.         RI1.IRIGEL(2,1) = 0
  623.         RI1.IRIGEL(3,1) = DESCR
  624.         RI1.IRIGEL(4,1) = xMATR1
  625.         RI1.IRIGEL(5,1) = NIFOUR
  626.         RI1.IRIGEL(6,1) = 0
  627.         RI1.IRIGEL(7,1) = 2
  628.         segact xmatr1*mod
  629.         xmatr1.symre=2
  630.         segdes xmatr1
  631.       SEGDES, RI1
  632.       IF (MODCOM) THEN
  633.         SEGINI, RI2 = RI1
  634.           RI2.IRIGEL(4,1) = xMATR2
  635.         SEGDES, RI2
  636.       ELSE
  637.         RI2 = 0
  638.         SEGSUP, MLCHP2
  639.       ENDIF
  640. *
  641.       iriout1 = ri1
  642.       iriout2 = ri2
  643.  
  644.  250  continue
  645.       mrigid = ipri2
  646.       segact mrigid
  647.       nrigel = coerig(/1)
  648.       if (nrigel.lt.1) goto 290
  649.       typmod = ' '
  650.  
  651.       nrigmat =100
  652.       kgmat = 0
  653.       segini prigmat
  654.  
  655.       KRIGEL = 0
  656.       nrigel = irigel(/2)
  657.       nrige = irigel(/1)
  658.       segini ri1
  659.       ri1.mtymat = mtymat
  660.       ri1.iforig = iforig
  661.       nrige0 = nrigel
  662.  
  663.       kige = 0
  664.       kige1 = 100
  665.       nrigel = kige1
  666.       segini ri2
  667.       ri2.mtymat = mtymat
  668.       ri2.iforig = iforig
  669.  
  670.       DO ire = 1,nrige0
  671.         meleme = irigel (1,ire)
  672.         segact meleme
  673.         if (itypel.ne.22) then
  674.            call erreur(977)
  675.            return
  676.         endif
  677.         nbelem = num(/2)
  678.         nbele0 = nbelem
  679.         descr = irigel(3,ire)
  680.         segact descr
  681.         nligrp0 = noelep(/1)
  682.         nligrd0 = noeled(/1)
  683.         nligrp = nligrp0 + nmapmo
  684.         nligrd = nligrd0 + nmapmo
  685.  
  686.         nbnn = num(/1)
  687.         nbsous = 0
  688.         nbref = 0
  689.         segini ipt2
  690.         ipt2.itypel = itypel
  691.         nbelem = 1
  692.         nbnn = nligrd
  693.         segini ipt1
  694.         ipt1.itypel = itypel
  695.         ri1.coerig(ire) = coerig(ire)
  696.         kele = 0
  697.  
  698.         xmatr1 = irigel(4,ire)
  699.         segact xmatr1
  700.         nelrig0 = xmatr1.re(/3)
  701.         nelrig = nelrig0 + nmapmo
  702.         rigrel=0
  703.         segini xmatr2
  704.         DO iele = 1,nbele0
  705.           ie2 = min(iele,nelrig0)
  706. *          xmatr1 = imatr1.imattt(ie2)
  707. *          segact xmatr1
  708.           nligrp = nligrp0 + nmapmo
  709.           nligrd = nligrd0 + nmapmo
  710.           nelrig=1
  711.          rigrel=0
  712.           segini des2,xmatri
  713.           des2.lisinc(1) = 'LX'
  714.           des2.lisdua(1) = 'FLX'
  715.           des2.noelep(1) = 1
  716.           des2.noeled(1) = 1
  717. * le premier point correspond aux multiplicateurs
  718.           CALL CREPO1 (ZERO, ZERO, ZERO, IPTS)
  719.           ipt1.num(1,1) = ipts
  720.           kgrp = 1
  721.           kirp = 1
  722.           do ipmo = 1,nmapmo
  723.             coepmo(ipmo) = 0.d0
  724.           enddo
  725.           do igrp = 2,nligrp0
  726.              jno = noelep(igrp)
  727.              motinc = lisinc(igrp)
  728.              IP1 = num(jno,iele)
  729. * recherche association noeud physique - points support déformée
  730.              do ilmat = 1,kgmat
  731.               if(lrigmat(ilmat,1).eq.ip1) goto 315
  732.              enddo
  733.  
  734.              kgmat = kgmat+1
  735.              ilmat = kgmat
  736.              if (kgmat.gt.nrigmat) then
  737.               nrigmat = nrigmat + 100
  738.               segadj prigmat
  739.              endif
  740.              kpb = 0
  741.              jg = 100
  742.              segini mlent3
  743.              lrigmat(kgmat,1) = ip1
  744.              do ikmo = 1, nmapmo
  745.                ichp1 = defpmo(ikmo)
  746.                call ecrcha('NOMU')
  747.                call ecrcha('MAIL')
  748.                call ecrobj('CHPOINT ',ichp1)
  749.                call extrai
  750.                call ecrobj('POINT   ',IP1)
  751.                call DANS
  752.                call lirlog(l3,1,iretou)
  753.                if(l3) then
  754.                  kpb = kpb + 1
  755.                  if (kpb.gt.jg) then
  756.                   jg = jg + 100
  757.                   segadj mlent3
  758.                  endif
  759.                  mlent3.lect(kpb) = ikmo
  760.                endif
  761.              enddo
  762.              jg = kpb
  763.              segadj mlent3
  764.              if (kpb.gt.0) then
  765.                   lrigmat(ilmat,2) = mlent3
  766.              else
  767.                   lrigmat(ilmat,2) = 0
  768.                 segsup mlent3
  769.              endif
  770.  
  771.  315  continue
  772.              ilr3 = lrigmat(ilmat,2)
  773.              if (ilr3.eq.0) goto 253
  774.              mlent3 = ilr3
  775.              segact mlent3
  776. * selection selon nom composante
  777.              mlmat = 0
  778.              do lmo = 1,9
  779.               if (motinc.eq.lcod(lmo)) mlmat = lmo+2
  780.              enddo
  781.              if (mlmat.eq.0) then
  782. *              WRITE(6,*) 'coefs pour cette composante non trouves'
  783.               goto 253
  784.              endif
  785.             if (lrigmat(ilmat,mlmat).ne.0) then
  786.               pcompo = lrigmat(ilmat,mlmat)
  787.               segact pcompo
  788.               nipmod = valmod(/1)
  789.               do ilg = 1,nipmod
  790.                 lkmo = mlent3.lect(ilg)
  791.                 coepmo(lkmo) = (valmod(ilg)* xmatr1.re(1,igrp,ie2))
  792.      &                     + coepmo(lkmo)
  793.               enddo
  794.             else
  795.              jg = mlent3.lect(/1)
  796.              nipmod = jg
  797.              segini pcompo
  798.              mcol = motinc
  799.              do ilg = 1,nipmod
  800.               lkmo = mlent3.lect(ilg)
  801.               ichp1 = defpmo(lkmo)
  802.               CALL EXTRA9(ICHP1,ip1,motinc,0,.false.,XFLOT,IRET)
  803.               coepmo(lkmo) = (xflot * xmatr1.re(1,igrp,ie2))
  804.      &                       + coepmo(lkmo)
  805.               valmod(ilg) = xflot
  806.              enddo
  807.              lrigmat(ilmat,mlmat) = pcompo
  808.             endif
  809.  
  810.   253    continue
  811.          enddo
  812.  
  813.          xmaut1 = 0.d0
  814.          do kpmo = 1,nmapmo
  815.           xmaut1 = max(xmaut1,ABS(coepmo(kpmo)))
  816.          enddo
  817.  
  818. * synthèse
  819.           do igrp = 2,nligrp0
  820.              jno = noelep(igrp)
  821.              motinc = lisinc(igrp)
  822.              IP1 = num(jno,iele)
  823.              lr2 = .false.
  824.             do jgmat = 1,kgmat
  825.               if(lrigmat(jgmat,1).eq.ip1) goto 325
  826.             enddo
  827. c            WRITE(6,*) 'bizarre, point dans l element non repertorie'
  828.             call erreur(5)
  829.             return
  830.  325        continue
  831.              mlmat = 0
  832.              do lmo = 1,9
  833.               if (motinc.eq.lcod(lmo)) mlmat = lmo+2
  834.              enddo
  835.              if (mlmat.eq.0) lr2 = .true.
  836.              if (lrigmat(jgmat,mlmat).eq.0) lr2 = .true.
  837.             if(lr2) then
  838.              jirp = 0
  839.              do iirp = 1,kgrp
  840.                if (ipt1.num(iirp,1).eq.ip1) jirp = iirp
  841.              enddo
  842. c recopie
  843.              kgrp = kgrp + 1
  844.                if (jirp.ne.0) then
  845.                 des2.noelep(kgrp) = des2.noelep(jirp)
  846.                 des2.noeled(kgrp) = des2.noeled(jirp)
  847.                else
  848.                 kirp = kirp + 1
  849.                 ipt1.num(kirp,1) = ip1
  850.                 des2.noelep(kgrp) = kirp
  851.                 des2.noeled(kgrp) = kirp
  852.                endif
  853.               des2.lisinc(kgrp) = lisinc(igrp)
  854.               des2.lisdua(kgrp) = lisdua(igrp)
  855.               re(1,kgrp,1) = xmatr1.re(1,igrp,ie2)
  856.               re(kgrp,1,1) = re(1,kgrp,1)
  857.             endif
  858. *
  859.           enddo
  860.  
  861.           do kpmo = 1,nmapmo
  862.             if (ABS(coepmo(kpmo)).gt.xlopre*xmaut1) then
  863.               kirp = kirp + 1
  864.               kgrp = kgrp + 1
  865.               ipt1.num(kirp,1) = lmapmo(kpmo)
  866.               des2.noelep(kgrp) = kirp
  867.               des2.noeled(kgrp) = kirp
  868.               motinc = compmo(kpmo)
  869.               des2.lisinc(kgrp) = motinc
  870.               if (motinc.eq.'ALFA    ') des2.lisdua(kgrp) = 'FALF    '
  871.               if (motinc.eq.'BETA    ') des2.lisdua(kgrp) = 'FBET    '
  872.               re(1,kgrp,1) = coepmo(kpmo)
  873.               re(kgrp,1,1) = re(1,kgrp,1)
  874.             endif
  875.           enddo
  876. *
  877.           lirl = .false.
  878.           if (kirp.ne.num(/1)) then
  879.             lirl = .true.
  880.           else
  881.             do io = 1,kirp
  882.               if (num(io,iele).ne.ipt1.num(io,1)) lirl=.true.
  883.             enddo
  884.           endif
  885. c creation d'un irigel
  886.           if (lirl) then
  887.            kige = kige + 1
  888.            if (kige.gt.kige1) then
  889.             nrigel = kige1 + 100
  890.             segadj ri2
  891.             kige1 = nrigel
  892.            endif
  893.            nbelem = 1
  894.            nbnn = kirp
  895.            segini ipt3
  896.            ipt3.itypel = itypel
  897.            do io =1,nbnn
  898.             ipt3.num(io,1) = ipt1.num(io,1)
  899.            enddo
  900.            nligrp = kgrp
  901.            nligrd = kgrp
  902.                    nelrig=1
  903.            rigrel=0
  904.            segadj xmatri,des2
  905. *           segini imatr3
  906. *           imatr3.imattt(1) = xmatri
  907.            segdes ipt3,des2,xmatri
  908.            RI2.IRIGEL(1,kige) = IPT3
  909.            RI2.IRIGEL(3,kige) = DES2
  910.            RI2.IRIGEL(4,kige) = xmatri
  911.            RI2.IRIGEL(2,kige) = 0
  912.            RI2.IRIGEL(5,kige) = irigel(5,ire)
  913.            RI2.IRIGEL(6,kige) = irigel(6,ire)
  914.            ri2.coerig(kige) = coerig(ire)
  915.           else
  916. * relation non modifiee pour cet element
  917.             kele = kele + 1
  918.             do ig = 1,nligrp0
  919.              ipt2.num(ig,kele) = ipt1.num(ig,1)
  920.             enddo
  921. *             imatr2.imattt(kele) = xmatr1
  922. * kich : a tester
  923.             do ju = 1,kgrp
  924.               xmatr2.re(1,ju,kele) = re(1,ju,1)
  925.               xmatr2.re(ju,1,kele) = re(ju,1,1)
  926.             enddo
  927.              segsup xmatri,des2
  928.           endif
  929.         ENDDO
  930.  
  931.         nbelem = kele
  932.         nelrig = kele
  933.         nligrd=xmatr2.re(/1)
  934.         nligrp=xmatr2.re(/2)
  935.         if (nbelem.gt.0) then
  936.           segadj ipt2
  937.           rigrel=0
  938.           segadj xmatr2
  939.           krigel = krigel + 1
  940.           RI1.IRIGEL(1,krigel) = IPT2
  941.           RI1.IRIGEL(3,krigel) = irigel(3,ire)
  942.           RI1.IRIGEL(4,krigel) = xmatr2
  943.           RI1.IRIGEL(2,krigel) = 0
  944.           RI1.IRIGEL(5,krigel) = irigel(5,ire)
  945.           RI1.IRIGEL(6,krigel) = irigel(6,ire)
  946.           segdes ipt2,xmatr2
  947.         else
  948.           segsup ipt2
  949.         endif
  950.         segsup ipt1
  951.        ENDDO
  952.  
  953.         iriout = 0
  954.         nrigel = krigel
  955.         segadj ri1
  956.         nrigel = kige
  957.         segadj ri2
  958.         segdes mrigid,ri1,ri2
  959.        if (kige.eq.0) segsup ri2
  960.        if (krigel.eq.0) segsup ri1
  961.        if (kige.gt.0.and.krigel.gt.0) then
  962. c        WRITE(6,*) 'fus', ri1,ri2,kige,krigel
  963.          call fusrig(ri1,ri2,iriout)
  964.          segsup ri1, ri2
  965.          return
  966.        endif
  967.        if (kige.gt.0) iriout = ri2
  968.        if (krigel.gt.0) iriout = ri1
  969.        if (iriout.eq.0) call erreur(-5)
  970. c       WRITE(6,*) 'iriout', iriout
  971.  
  972.  290  continue
  973.       if (iriout.ne.0) iriout3 = iriout
  974.       if (iriout1.ne.0) iriout3 = iriout1
  975.       if (iriout.ne.0.and.iriout1.ne.0) then
  976.          call fusrig(iriout, iriout1,iriout3)
  977.          ri1 = iriout
  978.          ri2 = iriout1
  979.          segsup ri1,ri2
  980.       endif
  981.  
  982.       call ecrobj('RIGIDITE',iriout3)
  983.       if (modcom)  call ecrobj('RIGIDITE',iriout2)
  984.  
  985.       goto 999
  986.  
  987.  199  continue
  988.       segsup  descr,meleme,mlchp1,mlchp2
  989.       call erreur(5)
  990.       return
  991.  
  992.  999  continue
  993.  
  994.       if (plcf.ne.0) segsup plcf
  995.  
  996.       END
  997.  
  998.  
  999.  
  1000.  
  1001.  
  1002.  
  1003.  
  1004.  

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