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

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