Télécharger la plaquette Cast3M


fr en

Télécharger supdep.eso

Retour à la liste

Numérotation des lignes : 

supdep

  1. C SUPDEP    SOURCE    PV090527  23/12/20    21:15:08     11813          
  2.       SUBROUTINE SUPDEP
  3. c=================================================================
  4. c cette procedure est appelée par SUPER
  5. c A partir du champ de deplacements interface et du champ d efforts
  6. c sur la sous structure, elle determine le champ de déplacement
  7. c sur toute la sous-structure
  8. c=================================================================
  9.       IMPLICIT INTEGER(I-N)
  10.       IMPLICIT REAL*8(A-H,O-Z)
  11. -INC SMSUPER
  12. -INC SMELEME
  13.  
  14. -INC PPARAM
  15. -INC CCOPTIO
  16. -INC SMMATRI
  17. -INC SMCHPOI
  18. -INC SMVECTD
  19. -INC SMRIGID
  20. -INC CCREEL
  21. c ????????????????????????????????????????????????????????????????
  22.       SEGMENT/BID/(BIDON(IIMAX+10)*D)
  23. c ????????????????????????????????????????????????????????????????
  24. c de la meme facon que dans MONDES , on va charger la maximum de
  25. c lignes en memoire pendant que l'on va gerer les efforts
  26. c IPLIG : première ligne en mémoire
  27. c IDLIG : dernière ligne en mémoire
  28. c----------------------------------------------------------------
  29.         PARAMETER (LPCRAY=10000)
  30.         segment ITRAA(NENS)
  31.         integer OOOVAL
  32.         logical leffort
  33.         character*4 mcle(1)
  34.         data mcle/'NOER'/
  35.  
  36.         dnorma=0.d0
  37.         dnormb=0.d0
  38. c creation d'une pile MRU
  39.         call OOOMRU(1)
  40.         call lirobj('SUPERELE',MSUPER,1,IRETOU)
  41.         if (ierr.ne.0) RETURN
  42. c lecture du déplacement interface
  43.         call lirobj('CHPOINT ',MDEPI,1,IRETOU)
  44.         if (ierr.ne.0) RETURN
  45. c lecture du  champ d'efforts si il existe
  46. c  27/11/2009  lecture obligatoire car risque d'erreur utilisateur
  47. c              due  a une mauvaise comprehension de l'operation
  48.         call lirobj('CHPOINT ',MFORC,1,IRETOU)
  49.         if (ierr.ne.0) RETURN
  50.         if (iretou .eq. 0) then
  51.             leffort = .false.
  52.         else
  53.             leffort = .true.
  54.         end if
  55.         noer=0
  56.         call lirmot(mcle,1,noer,0)
  57.  
  58.         segact MSUPER
  59. *       write(6,*) MRIGTO,MSUPEL,MSURAI,MBLOQU,MSUMAS,MCROUT
  60. c
  61. c recherche du nombre d'inconnues maitres
  62.         MRIGID = MSURAI
  63.         segact MRIGID
  64.         xMATRI = IRIGEL(4,1)
  65.         lagdua=msuper.islag
  66. *       write (6,*) ' msurai lagdua dans supdep ', mrigid,lagdua
  67.         segdes MRIGID
  68.         segact xMATRI
  69. *        XMATRI = IMATTT(1)
  70. *        segdes IMATRI
  71. *        segact XMATRI
  72.         nligra = RE(/1)
  73.         segdes XMATRI
  74. c
  75.         MCROUX = MCROUT
  76.         if(mcroux.eq.0) call erreur(1123)
  77.         if (ierr.ne.0) return
  78. * creation du meleme  des noeuds maitres
  79.        MRIGID=MRIGTO
  80.        SEGACT MRIGID
  81.  
  82.         meleme=msupel
  83.  
  84.         MRIGID = MSURAI
  85. c sauvegarde du deplacement interface
  86.         call reduir(MDEPI,MELEME,MDEPint)
  87.         call adchpo(mdepi,mdepint,mdext,1.d0,-1.d0)
  88. *       call ecchpo(mdepi,0)
  89. *       call ecchpo(mdepint,0)
  90. *        write (6,*) ' mdext '
  91. *       call ecchpo(mdext,0)
  92. c transformation du champ de déplacement en vecteur
  93. *       call ecmail(lagdua,0)
  94. **      call dbbch(mdepint,lagdua)
  95.         call chv1(MCROUX,MDEPInt,MVECTX,1)
  96. c normalisation
  97.         MMATRI = MCROUT
  98.         segact MMATRI
  99.         MDNOR = IDNORM
  100.         segact MDNOR
  101.         INC = DNOR(/1)
  102. c inbine : intger nombre inconnues esclaves
  103.  
  104.         inbine = INC - Nligra
  105. *       write (6,*) ' inbine inc nligra ',inbine,inc,nligra
  106.         MVECTD = MVECTX
  107.         segact MVECTD*mod
  108.         dnormb= 0.D0
  109.         do 1 ii1 = inbine+1,INC
  110.            VECTBB(ii1) = VECTBB(ii1) / DNOR(ii1)
  111. **         dnormb= max(dnormb,abs(VECTBB(ii1)))
  112.  1      continue
  113. **       dnormb = dnormb * xzprec*xzprec
  114.          dnormb = xpetit / xzprec
  115.  
  116.         MVECT1 = MVECTD
  117.         MILIGN = IILIGN
  118.         segact MILIGN
  119.         NNOE = ILIGN(/1)
  120.  
  121. c ????? idem que dans MONDES ??????
  122.         LTRK = MAX(1+LPCRAY,OOOVAL(1,4))
  123.         IIMAX = (((IJMAX+LTRK)/LTRK)+1)*LTRK
  124. c
  125. ***********************************************************************
  126. ***********************************************************************
  127. c pb à résoudre : L1t X +  L2t DEPI = D(-1) L1(-1) F
  128. c
  129. *  commencons par traiter le terme en effort
  130. c transformation du CHPOINT en vecteur
  131.         if (leffort) then
  132.           MCHPOI = MFORC
  133.           MRIGID = MRIGTO
  134.           segdes MRIGID
  135.           call copie2(mchpoi,mchpo1)
  136.           SEGACT MCHPO1
  137.         DO 432 I=1,mchpo1.IPCHP(/1)
  138.            MSOUPO=mchpo1.IPCHP(I)
  139.            SEGACT MSOUPO*MOD
  140.            IPT4=IGEOC
  141.            SEGINI,ipt5=ipt4
  142.            SEGDES ipt4
  143.            IGEOC=ipt5
  144.   432     CONTINUE
  145.           call dbbch(mchpo1,lagdua)
  146. *         call ecchpo(mchpo1 ,0)
  147.           call chv2(MCROUX,MCHPO1,MVECTX,1)
  148.           MVECTD = MVECTX
  149. c attention CHV2 desactive tout
  150.           segact MMATRI
  151.           segact MILIGN
  152.           segact MVECTD*MOD
  153.           MDIAG = IDIAG
  154.           segact MDIAG
  155. c
  156. c il faut normaliser ce vecteur
  157. c on va de plus recuperer l'indice du premier terme non nul
  158.           dnorma = 0.D0
  159.           inbi=vectbb(/1)
  160. *         write (6,*) ' inbine ',inbine
  161. *         write (6,*) ' vectbb-1 ',(vectbb(i),i=1,vectbb(/1))
  162. *         write (6,*) ' dnor     ',(dnor  (i),i=1,vectbb(/1))
  163.  
  164.  
  165.           do 2 ii1=1,inbine
  166.               VECTBB(ii1) = VECTBB(ii1) * DNOR(ii1)
  167. **            dnorma = max(dnorma,abs(VECTBB(ii1)))
  168.  2        continue
  169. *         write (6,*) ' vectbb-2 ',(vectbb(i),i=1,vectbb(/1))
  170. **        dnorma = dnorma * xzprec*xzprec
  171.           dnorma = xpetit/xzprec
  172. c
  173. c recherche du premier terme non nul
  174.           do 3 ii1=1,inbine
  175.              if (abs(VECTBB(ii1)).GE.dnorma) GOTO 4
  176.  3        continue
  177. *          write(6,*) ' attention vecteur force nulle'
  178.  4        continue
  179.           ipremf = ii1
  180.           if(dnorma.eq.0.d0) ipremf=max(1,inbine)
  181.  
  182. * effectuons maintenant le produit (L1)-1 EFFORT
  183. *-------------------------------------------------
  184. c ????????????
  185.           segini BID
  186.           ivalma = 0
  187. c ????????????
  188. c d'ou le bloc associé
  189.           INPREM = IPNO(ipremf)
  190. c et le bloc de la dernière ligne
  191.           INDERn=-1
  192.           if(inbine.ne.0)INDERN = IPNO(inbine)
  193.  
  194. c
  195.           IPLIG = INPREM
  196.           do 10 ii0=INPREM,INDERN
  197.              LIGN = ILIGN(ii0)
  198.              segact /ERR=20/ LIGN
  199.              IPRELL = IPREL
  200.              IVALMA = IVALMA+VAL(/1)
  201. *pvpv            if (IVALMA.GT.NGMAXY) GOTO 20
  202.              NA = IMMM(/1)
  203.              call supde2(IPPVV(1),IVPO(1),VAL(1),VECTBB(1),
  204.      >         na,inumli,inbine,ipremf,iprel,dnorma)
  205.  10       continue
  206.           segsup BID
  207.           IDLIG = INDERN
  208.           GOTO 60
  209. c on a eu des problèmes de memoire
  210.  20       continue
  211.           SEGSUP BID
  212.           IDLIG = ii0 - 1
  213.           itemp1 = ii0
  214. c c'est reparti
  215.           do 30 ii0=itemp1,INDERN
  216.             LIGN = ILIGN(ii0)
  217.             segact /ERR=22/ LIGN
  218.             NA = IMMM(/1)
  219.             IPRELL = IPREL
  220.             call supde2(IPPVV(1),IVPO(1),VAL(1),VECTBB(1),
  221.      >         na,inumli,inbine,ipremf,iprel,dnorma)
  222.             segdes LIGN*(NOMOD,MRU)
  223.             GOTO 30
  224. c il y a encore des problèmes
  225.  22         CONTINUE
  226. c il va falloir faire de la place ==> desactivation d'une ligne
  227. c reste-t-il des lignes à desactiver ?
  228.             if (IDLIG .lt. IPLIG) then
  229.               CALL ERREUR(5)
  230.             else
  231.               LIGN=ILIGN(IDLIG)
  232.               SEGDES LIGN*(NOMOD,MRU)
  233.               IDLIG=IDLIG-1
  234.             end if
  235.  30       continue
  236. c
  237. c muliplions maintenant par D-1
  238. c---------------------------------
  239.  60       continue
  240.           do 70 i=1,inbine
  241.              VECTBB(i) = VECTBB(i) * DIAG(i)
  242.  70       continue
  243.           segdes MDIAG
  244.         else
  245. c il n y a pas d'efforts
  246.            segini MVECTD
  247. c il faut preciser qu il n y a aucune ligne chargée
  248.            IDLIG = 0
  249.            IPLIG = NNOE
  250.         end if
  251. c
  252. c
  253. ***********************************************************************
  254. ***********************************************************************
  255. c Il faut maintenant effectuer le calcul L1t X1 + L2t X2 = Nouveau F
  256. c il ne s'agit que d'une remontée
  257. c
  258. c on va en profiter pour compter les mouvements d'ensemble
  259.         iaa = 0
  260.         segini itraa
  261. c
  262. c inumli : numero de la ligne en cours
  263.         inumli = INC
  264.         do 120 ii0=NNOE,1,-1
  265.            LIGN = ILIGN(ii0)
  266.  122       continue
  267.            if ((ii0.gt.IDLIG).or.(ii0.lt.IPLIG)) then
  268.              segact /ERR=124/ LIGN
  269.            else
  270.              IDLIG = IDLIG - 1
  271.            end if
  272.            NA = IMMM(/1)
  273.          IFIB=IVPO(/1)
  274.            call supde1(IPPVV(1),IVPO(1),VAL(1),VECTBB(1),
  275.      &      MVECT1.VECTBB(1),na,inumli,inbine,iprel,ifib,dnormb)
  276.            do ibb=1,NA
  277.             if (IMMM(ibb) .ne. 0) then
  278.              iaa = iaa+1
  279.              itraa(iaa)=iprel+ibb-1
  280.             end if
  281.            end do
  282.            segdes LIGN*(NOMOD,MRU)
  283.            GOTO 120
  284.  124       CONTINUE
  285. c encore des problèmes mémoires
  286.            if  (IDLIG .lt. IPLIG) then
  287.              CALL ERREUR(5)
  288.            else
  289.              LIG1=ILIGN(IDLIG)
  290.              SEGDES LIG1*(NOMOD,MRU)
  291.              IDLIG=IDLIG-1
  292.                  GOTO 122
  293.            end if
  294.  120    continue
  295.  
  296.         segsup MVECTD
  297.         MVECTD = MVECT1
  298. c
  299. *******************************************************************************
  300. *******************************************************************************
  301. c gestion des déplacements de corps rigide
  302. c meme chose que dans MONDES
  303. c                     ------
  304. c
  305.       if (nens .ne. 0) then
  306.         MINCPO = IINCPO
  307.         MIMIK = IIMIK
  308.         segact MINCPO,MIMIK
  309.         ipt2 = IGEOMA
  310.         segact ipt2
  311.         NNOE = INCPO(/2)
  312.         IINC1=INCPO(/1)
  313. c
  314.         XMA = xpetit
  315.         XMAL = xpetit
  316.         inan=0
  317.         do kk=1,INC
  318.           if (NOER.EQ.0.OR.(noer.eq.1.and.abs(vectbb(kk)).lt.xgrand))
  319.      &           goto 500
  320.           inan = inan + 1
  321.           vectbb(kk)=0.D0
  322.  500   continue
  323.          if (ittr(kk).eq.0) then
  324.            XMA=MAX(XMA,ABS(VECTBB(kk)))
  325.          else
  326.            XMAL=MAX(XMAL,ABS(VECTBB(kk)))
  327.          endif
  328.         end do
  329.         XMA = XMA * 1.d-10
  330.         XMAL = XMAL * 1.d-10
  331.         xmal = max(xma*1d-2,xmal)
  332. *       write (6,*) ' supdep cma cmal ',xma,xmal
  333. c boucle sur les mouvements d'ensemble
  334.         do 200 ia=1,NENS
  335.           i1=itraa(ia)
  336.           j=IPNO(i1)
  337.           do 210 k=1,iinc1
  338.             if(INCPO(K,J).eq.i1) goto 220
  339.  210      continue
  340.           call erreur(5)
  341.           return
  342.  220      continue
  343.           if ((ittr(i1).eq.0).and.(ABS(VECTBB(i1)).le.XMA)) GOTO 250
  344.           if ((ittr(i1).ne.0).and.(ABS(VECTBB(i1)).le.XMAL)) GOTO 250
  345.           MOTERR(1:4)=IMIK(k)
  346.           INTERR(1)=ipt2.NUM(1,J)
  347.           write (6,*) ' vectbb) ',vectbb(i1)
  348. **        CALL ERREUR(149)
  349. **        RETURN
  350.           call soucis(149)
  351.  251      continue
  352.  250      continue
  353.           jjk = ipt2.NUM(1,J)
  354.           IF(ITTR(I1).EQ.0) WRITE(IOIMP,280) JJK,IMIK(K)
  355.  
  356.           IF (IIMPI.NE.0 .AND. ITTR(I1).NE.0)
  357.      &                WRITE(IOIMP,290) JJK,IMIK(K)
  358.  
  359.  200    continue
  360.   280 FORMAT(' INDETERMINATION DETECTEE AU NOEUD ',I6,' INCONNUE ',
  361.      *   A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ',
  362.      * 'LA SUSDITE')
  363.   290 FORMAT(' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ',
  364.      * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ',
  365.      * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE ')
  366.  
  367.         segdes MINCPO,MIMIK
  368.           segdes ipt2
  369.       end if
  370. c
  371. *******************************************************************************
  372. *******************************************************************************
  373. * un petit coup de normalisation
  374.         do 300 ii1=1,INC
  375.             VECTBB(ii1) = VECTBB(ii1) * DNOR(ii1)
  376.  300    continue
  377.  
  378.         segdes MDNOR
  379.         segdes MVECTD
  380.         MVECTX = MVECTD
  381.         MMATRX = MMATRI
  382.         MRIGTX = MRIGTO
  383. *        write(6,*) ' mrigto mmatri mvectx ', mrigto,mmatri,mvectx
  384. c il faut transformer ce vecteur en chpoint
  385.         call VCH1(MMATRX,MVECTX,ISOLU,MRIGTX)
  386.         if (lagdua.ne.0) call dbbcf(isolu,lagdua)
  387. **      write (6,*) ' mdext '
  388. **      call ecchpo(mdext,0)
  389. **      write (6,*) ' isolu '
  390. **      call ecchpo(isolu,0)
  391.         call adchpo(isolu,mdext,mchpoi,1.d0,1.d0)
  392. *       mchpoi=isolu
  393.         segact MCHPOI*mod
  394.         JATTRI(1)=1
  395.         segdes MCHPOI
  396.  
  397.         segdes MMATRI
  398. c il faut maintenant rajouter les déplacements imposés éliminés (jrcond)
  399.         MRIGID = MRIGTO
  400. *        segact MRIGID
  401.         segdes MRIGID
  402. c
  403.         segdes MSUPER
  404. c il faut ajuster le type du champ
  405.  
  406. c on rajoute les multiplicateurs interfaces
  407.         call ecrobj('CHPOINT ',mchpoi)
  408.  
  409. c désactivation de la pile MRU
  410.         call OOOMRU(0)
  411.         return
  412.         end
  413.  
  414.  
  415.  
  416.  
  417.  
  418.  
  419.  
  420.  
  421.  
  422.  
  423.  
  424.  
  425.  
  426.  
  427.  
  428.  
  429.  
  430.  
  431.  
  432.  

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