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    CB215821  25/04/23    21:15:32     12247          
  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.           SEGINI, XMATR1
  536.           IF (MODCOM) THEN
  537.             SEGINI, XMATR2
  538.             SEGDES, XMATR2
  539.           ENDIF
  540.           SEGDES, XMATR1
  541.           GOTO 55
  542.         ENDIF
  543. *
  544. * Cas reel
  545. *
  546.         SEGINI, XMATR1
  547.           DO I=1, NBMOD
  548.             MCHPO1 = MLCHP1.ICHPOI(I)
  549.             DO J=1, NBMOD
  550.               MCHPO2 = MLCHP1.ICHPOI(J)
  551.               CALL YTMX (MCHPO2, MCHPO1, MRIGID, XVAL)
  552.               XMATR1.RE(I,J,1)=XVAL
  553.             ENDDO
  554.           ENDDO
  555.         SEGDES, XMATR1
  556. *
  557. * Cas complexe : calcul de termes complementaires
  558. *
  559.         IF (MODCOM) THEN
  560.           SEGACT, XMATR1*mod
  561.             DO I=1, NBMOD
  562.               MCHPO1 = MLCHP2.ICHPOI(I)
  563.               DO J=1, NBMOD
  564.                 MCHPO2 = MLCHP2.ICHPOI(J)
  565.                 CALL YTMX (MCHPO1, MCHPO2, MRIGID, XVAL)
  566.                 XMATR1.RE(I,J,1)=XMATR1.RE(I,J,1)-IREEL*XVAL
  567.               ENDDO
  568.             ENDDO
  569.           SEGDES, XMATR1
  570. *
  571.           SEGINI, XMATR2
  572.             DO I=1, NBMOD
  573.               MCHPO1 = MLCHP1.ICHPOI(I)
  574.               DO J=1, NBMOD
  575.                 MCHPO2 = MLCHP2.ICHPOI(J)
  576.                 CALL YTMX (MCHPO1, MCHPO2, MRIGID, XVAL)
  577.                 XMATR2.RE(I,J,1)=XVAL
  578.               ENDDO
  579.             ENDDO
  580.             DO I=1, NBMOD
  581.               MCHPO1 = MLCHP2.ICHPOI(I)
  582.               DO J=1, NBMOD
  583.                 MCHPO2 = MLCHP1.ICHPOI(J)
  584.                 CALL YTMX (MCHPO1, MCHPO2, MRIGID, XVAL)
  585.                 XMATR2.RE(I,J,1)=XMATR2.RE(I,J,1)+IREEL*XVAL
  586.               ENDDO
  587.             ENDDO
  588.           SEGDES, XMATR2
  589.         ENDIF
  590. *
  591.       SEGACT, MRIGID
  592.         LETYPE = MRIGID.MTYMAT
  593.       SEGDES, MRIGID
  594. *
  595. * Creation des segments IMATRI
  596. *
  597.  55   NELRIG = 1
  598. *      SEGINI, IMATR1
  599. *        IMATR1.IMATTT(1) = XMATR1
  600.       SEGDES, xMATR1
  601.       IF (MODCOM) THEN
  602. *        SEGINI, IMATR2
  603. *          IMATR2.IMATTT(1) = XMATR2
  604.         SEGDES, xMATR2
  605.       ENDIF
  606. *
  607. * Creation des rigidites calculees
  608. *
  609.       NRIGE=7
  610.       NRIGEL=1
  611.       SEGINI, RI1
  612.         RI1.MTYMAT = LETYPE
  613.         RI1.IFORIG = IFOUR
  614.         RI1.IMGEO1 = 0
  615.         RI1.IMGEO2 = 0
  616.         RI1.COERIG = 1.D0
  617.         RI1.IRIGEL(1,1) = MELEME
  618.         RI1.IRIGEL(2,1) = 0
  619.         RI1.IRIGEL(3,1) = DESCR
  620.         RI1.IRIGEL(4,1) = xMATR1
  621.         RI1.IRIGEL(5,1) = NIFOUR
  622.         RI1.IRIGEL(6,1) = 0
  623.         RI1.IRIGEL(7,1) = 2
  624.         segact xmatr1*mod
  625.         xmatr1.symre=2
  626.         segdes xmatr1
  627.       SEGDES, RI1
  628.       IF (MODCOM) THEN
  629.         SEGINI, RI2 = RI1
  630.           RI2.IRIGEL(4,1) = xMATR2
  631.         SEGDES, RI2
  632.       ELSE
  633.         RI2 = 0
  634.         SEGSUP, MLCHP2
  635.       ENDIF
  636. *
  637.       iriout1 = ri1
  638.       iriout2 = ri2
  639.  
  640.  250  continue
  641.       mrigid = ipri2
  642.       segact mrigid
  643.       nrigel = coerig(/1)
  644.       if (nrigel.lt.1) goto 290
  645.       typmod = ' '
  646.  
  647.       nrigmat =100
  648.       kgmat = 0
  649.       segini prigmat
  650.  
  651.       KRIGEL = 0
  652.       nrigel = irigel(/2)
  653.       nrige = irigel(/1)
  654.       segini ri1
  655.       ri1.mtymat = mtymat
  656.       ri1.iforig = iforig
  657.       nrige0 = nrigel
  658.  
  659.       kige = 0
  660.       kige1 = 100
  661.       nrigel = kige1
  662.       segini ri2
  663.       ri2.mtymat = mtymat
  664.       ri2.iforig = iforig
  665.  
  666.       DO ire = 1,nrige0
  667.         meleme = irigel (1,ire)
  668.         segact meleme
  669.         if (itypel.ne.22) then
  670.            call erreur(977)
  671.            return
  672.         endif
  673.         nbelem = num(/2)
  674.         nbele0 = nbelem
  675.         descr = irigel(3,ire)
  676.         segact descr
  677.         nligrp0 = noelep(/1)
  678.         nligrd0 = noeled(/1)
  679.         nligrp = nligrp0 + nmapmo
  680.         nligrd = nligrd0 + nmapmo
  681.  
  682.         nbnn = num(/1)
  683.         nbsous = 0
  684.         nbref = 0
  685.         segini ipt2
  686.         ipt2.itypel = itypel
  687.         nbelem = 1
  688.         nbnn = nligrd
  689.         segini ipt1
  690.         ipt1.itypel = itypel
  691.         ri1.coerig(ire) = coerig(ire)
  692.         kele = 0
  693.  
  694.         xmatr1 = irigel(4,ire)
  695.         segact xmatr1
  696.         nelrig0 = xmatr1.re(/3)
  697.         nelrig = nelrig0 + nmapmo
  698.         segini xmatr2
  699.         DO iele = 1,nbele0
  700.           ie2 = min(iele,nelrig0)
  701. *          xmatr1 = imatr1.imattt(ie2)
  702. *          segact xmatr1
  703.           nligrp = nligrp0 + nmapmo
  704.           nligrd = nligrd0 + nmapmo
  705.           nelrig=1
  706.           segini des2,xmatri
  707.           des2.lisinc(1) = 'LX'
  708.           des2.lisdua(1) = 'FLX'
  709.           des2.noelep(1) = 1
  710.           des2.noeled(1) = 1
  711. * le premier point correspond aux multiplicateurs
  712.           CALL CREPO1 (ZERO, ZERO, ZERO, IPTS)
  713.           ipt1.num(1,1) = ipts
  714.           kgrp = 1
  715.           kirp = 1
  716.           do ipmo = 1,nmapmo
  717.             coepmo(ipmo) = 0.d0
  718.           enddo
  719.           do igrp = 2,nligrp0
  720.              jno = noelep(igrp)
  721.              motinc = lisinc(igrp)
  722.              IP1 = num(jno,iele)
  723. * recherche association noeud physique - points support déformée
  724.              do ilmat = 1,kgmat
  725.               if(lrigmat(ilmat,1).eq.ip1) goto 315
  726.              enddo
  727.  
  728.              kgmat = kgmat+1
  729.              ilmat = kgmat
  730.              if (kgmat.gt.nrigmat) then
  731.               nrigmat = nrigmat + 100
  732.               segadj prigmat
  733.              endif
  734.              kpb = 0
  735.              jg = 100
  736.              segini mlent3
  737.              lrigmat(kgmat,1) = ip1
  738.              do ikmo = 1, nmapmo
  739.                ichp1 = defpmo(ikmo)
  740.                call ecrcha('NOMU')
  741.                call ecrcha('MAIL')
  742.                call ecrobj('CHPOINT ',ichp1)
  743.                call extrai
  744.                call ecrobj('POINT   ',IP1)
  745.                call DANS
  746.                call lirlog(l3,1,iretou)
  747.                if(l3) then
  748.                  kpb = kpb + 1
  749.                  if (kpb.gt.jg) then
  750.                   jg = jg + 100
  751.                   segadj mlent3
  752.                  endif
  753.                  mlent3.lect(kpb) = ikmo
  754.                endif
  755.              enddo
  756.              jg = kpb
  757.              segadj mlent3
  758.              if (kpb.gt.0) then
  759.                   lrigmat(ilmat,2) = mlent3
  760.              else
  761.                   lrigmat(ilmat,2) = 0
  762.                 segsup mlent3
  763.              endif
  764.  
  765.  315  continue
  766.              ilr3 = lrigmat(ilmat,2)
  767.              if (ilr3.eq.0) goto 253
  768.              mlent3 = ilr3
  769.              segact mlent3
  770. * selection selon nom composante
  771.              mlmat = 0
  772.              do lmo = 1,9
  773.               if (motinc.eq.lcod(lmo)) mlmat = lmo+2
  774.              enddo
  775.              if (mlmat.eq.0) then
  776. *              WRITE(6,*) 'coefs pour cette composante non trouves'
  777.               goto 253
  778.              endif
  779.             if (lrigmat(ilmat,mlmat).ne.0) then
  780.               pcompo = lrigmat(ilmat,mlmat)
  781.               segact pcompo
  782.               nipmod = valmod(/1)
  783.               do ilg = 1,nipmod
  784.                 lkmo = mlent3.lect(ilg)
  785.                 coepmo(lkmo) = (valmod(ilg)* xmatr1.re(1,igrp,ie2))
  786.      &                     + coepmo(lkmo)
  787.               enddo
  788.             else
  789.              jg = mlent3.lect(/1)
  790.              nipmod = jg
  791.              segini pcompo
  792.              mcol = motinc
  793.              do ilg = 1,nipmod
  794.               lkmo = mlent3.lect(ilg)
  795.               ichp1 = defpmo(lkmo)
  796.               CALL EXTRA9(ICHP1,ip1,motinc,0,.false.,XFLOT,IRET)
  797.               coepmo(lkmo) = (xflot * xmatr1.re(1,igrp,ie2))
  798.      &                       + coepmo(lkmo)
  799.               valmod(ilg) = xflot
  800.              enddo
  801.              lrigmat(ilmat,mlmat) = pcompo
  802.             endif
  803.  
  804.   253    continue
  805.          enddo
  806.  
  807.          xmaut1 = 0.d0
  808.          do kpmo = 1,nmapmo
  809.           xmaut1 = max(xmaut1,ABS(coepmo(kpmo)))
  810.          enddo
  811.  
  812. * synthèse
  813.           do igrp = 2,nligrp0
  814.              jno = noelep(igrp)
  815.              motinc = lisinc(igrp)
  816.              IP1 = num(jno,iele)
  817.              lr2 = .false.
  818.             do jgmat = 1,kgmat
  819.               if(lrigmat(jgmat,1).eq.ip1) goto 325
  820.             enddo
  821. c            WRITE(6,*) 'bizarre, point dans l element non repertorie'
  822.             call erreur(5)
  823.             return
  824.  325        continue
  825.              mlmat = 0
  826.              do lmo = 1,9
  827.               if (motinc.eq.lcod(lmo)) mlmat = lmo+2
  828.              enddo
  829.              if (mlmat.eq.0) lr2 = .true.
  830.              if (lrigmat(jgmat,mlmat).eq.0) lr2 = .true.
  831.             if(lr2) then
  832.              jirp = 0
  833.              do iirp = 1,kgrp
  834.                if (ipt1.num(iirp,1).eq.ip1) jirp = iirp
  835.              enddo
  836. c recopie
  837.              kgrp = kgrp + 1
  838.                if (jirp.ne.0) then
  839.                 des2.noelep(kgrp) = des2.noelep(jirp)
  840.                 des2.noeled(kgrp) = des2.noeled(jirp)
  841.                else
  842.                 kirp = kirp + 1
  843.                 ipt1.num(kirp,1) = ip1
  844.                 des2.noelep(kgrp) = kirp
  845.                 des2.noeled(kgrp) = kirp
  846.                endif
  847.               des2.lisinc(kgrp) = lisinc(igrp)
  848.               des2.lisdua(kgrp) = lisdua(igrp)
  849.               re(1,kgrp,1) = xmatr1.re(1,igrp,ie2)
  850.               re(kgrp,1,1) = re(1,kgrp,1)
  851.             endif
  852. *
  853.           enddo
  854.  
  855.           do kpmo = 1,nmapmo
  856.             if (ABS(coepmo(kpmo)).gt.xlopre*xmaut1) then
  857.               kirp = kirp + 1
  858.               kgrp = kgrp + 1
  859.               ipt1.num(kirp,1) = lmapmo(kpmo)
  860.               des2.noelep(kgrp) = kirp
  861.               des2.noeled(kgrp) = kirp
  862.               motinc = compmo(kpmo)
  863.               des2.lisinc(kgrp) = motinc
  864.               if (motinc.eq.'ALFA    ') des2.lisdua(kgrp) = 'FALF    '
  865.               if (motinc.eq.'BETA    ') des2.lisdua(kgrp) = 'FBET    '
  866.               re(1,kgrp,1) = coepmo(kpmo)
  867.               re(kgrp,1,1) = re(1,kgrp,1)
  868.             endif
  869.           enddo
  870. *
  871.           lirl = .false.
  872.           if (kirp.ne.num(/1)) then
  873.             lirl = .true.
  874.           else
  875.             do io = 1,kirp
  876.               if (num(io,iele).ne.ipt1.num(io,1)) lirl=.true.
  877.             enddo
  878.           endif
  879. c creation d'un irigel
  880.           if (lirl) then
  881.            kige = kige + 1
  882.            if (kige.gt.kige1) then
  883.             nrigel = kige1 + 100
  884.             segadj ri2
  885.             kige1 = nrigel
  886.            endif
  887.            nbelem = 1
  888.            nbnn = kirp
  889.            segini ipt3
  890.            ipt3.itypel = itypel
  891.            do io =1,nbnn
  892.             ipt3.num(io,1) = ipt1.num(io,1)
  893.            enddo
  894.            nligrp = kgrp
  895.            nligrd = kgrp
  896.                    nelrig=1
  897.            segadj xmatri,des2
  898. *           segini imatr3
  899. *           imatr3.imattt(1) = xmatri
  900.            segdes ipt3,des2,xmatri
  901.            RI2.IRIGEL(1,kige) = IPT3
  902.            RI2.IRIGEL(3,kige) = DES2
  903.            RI2.IRIGEL(4,kige) = xmatri
  904.            RI2.IRIGEL(2,kige) = 0
  905.            RI2.IRIGEL(5,kige) = irigel(5,ire)
  906.            RI2.IRIGEL(6,kige) = irigel(6,ire)
  907.            ri2.coerig(kige) = coerig(ire)
  908.           else
  909. * relation non modifiee pour cet element
  910.             kele = kele + 1
  911.             do ig = 1,nligrp0
  912.              ipt2.num(ig,kele) = ipt1.num(ig,1)
  913.             enddo
  914. *             imatr2.imattt(kele) = xmatr1
  915. * kich : a tester
  916.             do ju = 1,kgrp
  917.               xmatr2.re(1,ju,kele) = re(1,ju,1)
  918.               xmatr2.re(ju,1,kele) = re(ju,1,1)
  919.             enddo
  920.              segsup xmatri,des2
  921.           endif
  922.         ENDDO
  923.  
  924.         nbelem = kele
  925.         nelrig = kele
  926.         nligrd=xmatr2.re(/1)
  927.         nligrp=xmatr2.re(/2)
  928.         if (nbelem.gt.0) then
  929.           segadj ipt2
  930.           segadj xmatr2
  931.           krigel = krigel + 1
  932.           RI1.IRIGEL(1,krigel) = IPT2
  933.           RI1.IRIGEL(3,krigel) = irigel(3,ire)
  934.           RI1.IRIGEL(4,krigel) = xmatr2
  935.           RI1.IRIGEL(2,krigel) = 0
  936.           RI1.IRIGEL(5,krigel) = irigel(5,ire)
  937.           RI1.IRIGEL(6,krigel) = irigel(6,ire)
  938.           segdes ipt2,xmatr2
  939.         else
  940.           segsup ipt2
  941.         endif
  942.         segsup ipt1
  943.        ENDDO
  944.  
  945.         iriout = 0
  946.         nrigel = krigel
  947.         segadj ri1
  948.         nrigel = kige
  949.         segadj ri2
  950.         segdes mrigid,ri1,ri2
  951.        if (kige.eq.0) segsup ri2
  952.        if (krigel.eq.0) segsup ri1
  953.        if (kige.gt.0.and.krigel.gt.0) then
  954. c        WRITE(6,*) 'fus', ri1,ri2,kige,krigel
  955.          call fusrig(ri1,ri2,iriout)
  956.          segsup ri1, ri2
  957.          return
  958.        endif
  959.        if (kige.gt.0) iriout = ri2
  960.        if (krigel.gt.0) iriout = ri1
  961.        if (iriout.eq.0) call erreur(-5)
  962. c       WRITE(6,*) 'iriout', iriout
  963.  
  964.  290  continue
  965.       if (iriout.ne.0) iriout3 = iriout
  966.       if (iriout1.ne.0) iriout3 = iriout1
  967.       if (iriout.ne.0.and.iriout1.ne.0) then
  968.          call fusrig(iriout, iriout1,iriout3)
  969.          ri1 = iriout
  970.          ri2 = iriout1
  971.          segsup ri1,ri2
  972.       endif
  973.  
  974.       call ecrobj('RIGIDITE',iriout3)
  975.       if (modcom)  call ecrobj('RIGIDITE',iriout2)
  976.  
  977.       goto 999
  978.  
  979.  199  continue
  980.       segsup  descr,meleme,mlchp1,mlchp2
  981.       call erreur(5)
  982.       return
  983.  
  984.  999  continue
  985.  
  986.       if (plcf.ne.0) segsup plcf
  987.  
  988.       END
  989.  
  990.  
  991.  
  992.  
  993.  
  994.  

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