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

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