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

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