Télécharger la plaquette Cast3M


fr en

Télécharger precop.eso

Retour à la liste

Numérotation des lignes : 

precop

  1. C PRECOP    SOURCE    OF166741  25/02/21    21:18:10     12166          
  2.       SUBROUTINE PRECOP(IPMODL,IPCHA1,IPTAB,IPSTRS,IRAN,
  3.      &                  PS1,IPCHC1,IRET)
  4. C======================================================================C
  5. C                                                                      C
  6. C     ENTREES :                                                        C
  7. C                                                                      C
  8. C       IPMODL: POINTEUR SUR UN MMODEL                                 C
  9. C       IPCHA1: Pointeur sur le MCHAML de CARACTERISTIQUES             C
  10. C       PS1   : Force sous vérin
  11. C       IRAN  :  Pointeur sur maillage des points d application
  12. C       IPCHC1: POINTEUR SUR le MCHAML DE PRECONTRAINTE INITIALE
  13. C       IPTAB : pointeur sur la table des parametres frottement ....   C
  14. C                                                                      C
  15. C     SORTIES :                                                        C
  16. C                                                                      C
  17. C       IPSTRS:  MCHAML CONTENANT LES PRECONTAINTES ET LES FORCES      C
  18. C                DU CABLE SUR LE BETON                                 C
  19. C       IRET  :  1   cela c est bien passe
  20. C                0   probleme dans le traitement                       C
  21. C                                                                      C
  22. C======================================================================C
  23.       IMPLICIT INTEGER(I-N)
  24.       IMPLICIT REAL*8(A-H,O-Z)
  25.  
  26. -INC PPARAM
  27. -INC CCOPTIO
  28. -INC CCHAMP
  29.  
  30. -INC SMCHAML
  31. -INC SMELEME
  32. -INC SMINTE
  33. -INC SMMODEL
  34. -INC SMCOORD
  35. -INC SMTABLE
  36.  
  37. -INC TMPTVAL
  38.  
  39.       SEGMENT INFO
  40.          INTEGER INFELL(JG)
  41.       ENDSEGMENT
  42. * stockage  de pointeurs sur des segment sielc  crees dans splitag
  43.       segment siezo
  44.       integer iezon(*)
  45.       endsegment
  46.       segment sielc
  47.       integer ideb(*),ifin(*),nbcz,isens(2,*),idejvu(*)
  48.       endsegment
  49. C
  50.       SEGMENT ALTRAV
  51.          REAL*8 ANG(NAM+1),DANG(NAM),ACUR(NAM+1),DACUR(NAM)
  52.       ENDSEGMENT
  53. C
  54.       SEGMENT NOTYPE
  55.         CHARACTER*16 TYPE(NBTYPE)
  56.       ENDSEGMENT
  57. C
  58.       SEGMENT WRK3
  59.          REAL*8 X1(3),X2(3),X3(3),RL(3),RS(3)
  60.       ENDSEGMENT
  61. C
  62.       REAL*8 XVALIN,XVALRE
  63.       LOGICAL LOGRE,LOGIN
  64.       CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE
  65.       CHARACTER*4 MOPAR(6)
  66.  
  67.       DATA MOPAR/'F1  ','F2  ','GANC','RMU0','FPRG','RH10'/
  68. C
  69.       CHARACTER*8 CMATE,LNOM
  70.       CHARACTER*(NCONCH) CONM
  71.       PARAMETER ( NINF=3 )
  72.       PARAMETER ( ITEMAX=50 )
  73.       INTEGER INFOS(NINF)
  74.       DIMENSION FGG(ITEMAX),XLAM(ITEMAX)
  75.       logical ivers,lsupco
  76.  
  77.       ivers=.false.
  78.  
  79. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  80.  
  81.       if (ipchc1.ne.0) then
  82.       MCHELM  = ipchc1
  83.       segact MCHELM
  84. C   actuellement  un cable n a qu une zone de maillage
  85.       imail = imache(1)
  86.       endif
  87.  
  88.       IRET=0
  89. C     verification on doit avoir des elements soit BARR soit  CERC
  90.       MMODEL = ipmodl
  91.       segact mmodel
  92.       do i = 1,kmodel(/1)
  93.       imodel = kmodel(i)
  94.       segact imodel
  95.       if(ipchc1.ne.0.and.imail.ne.imamod) then
  96.       call erreur(21)
  97.       return
  98.       endif
  99.       meleme = imamod
  100.       segact meleme
  101.       if(itypel.ne.2) then
  102.         call erreur(16)
  103.         return
  104.       endif
  105.       if(ifomod.eq.0.and.itypel.ne.1) then
  106.       call erreur(21)
  107.       return
  108.       endif
  109.       enddo
  110.  
  111. C --- Vérification du support des MCHAML --------------
  112.  
  113.       if(ipchc1.gt.0) then
  114.          call quesup(ipmodl,ipchc1,5,0,iresig,iret2)
  115.          if(iresig.eq.9999) then
  116.             moterr(1:8)='CONTRAIN'
  117.             call erreur(124)
  118.             return
  119.          endif
  120.        endif
  121.  
  122.       call quesup(ipmodl,ipcha1,5,0,irecar,iret2)
  123.       if(irecar.eq.9999) then
  124.          moterr(1:8)='CARACTER'
  125.          call erreur(124)
  126.          return
  127.       endif
  128.       if(irecar.eq.1) then
  129.          ipcara = 0
  130.          iret = 0
  131.          call chasup(ipmodl,ipcha1,ipcara,iret,5)
  132.          if(iret.ne.0) return
  133.       else
  134.          ipcara = ipcha1
  135.       endif
  136. *---------------------------------------------------------
  137.       ipmod0 = ipmodl
  138. *---------------------------------------------------------
  139.       if(ifomod.ne.0.and.ifomod.ne.1) then
  140.       ipcar0 = ipcara
  141.       ipchc0 = ipchc1
  142.       call splitag(ipmodl,ipcara,ipchc1,ipmod2,ipcar2,
  143.      &      ipchc2,siezo)
  144.       ipmodl = ipmod2
  145.       ipcara = ipcar2
  146.       ipchc1 = ipchc2
  147.       else
  148. C  pas elegant mais on duplique  le modele
  149.       MMODEL = ipmodl
  150.       segini , mmode2=MMODEL
  151.       ipmodl = mmode2
  152.       endif
  153. C
  154.       ipmodt = ipmodl
  155.       MMODEL = IPMODL
  156. C
  157.       SEGACT MMODEL
  158.       NSOUS=KMODEL(/1)
  159. C
  160. C ----PREPARATION DU CHAMP DE SORTIE----------------
  161. C
  162.       mmode2 = ipmod0
  163.       segact mmode2
  164.       N1=mmode2.kmodel(/1)
  165.       L1=11
  166.       N3=6
  167.       segini mchel5
  168.       mchel5.titche='CONTRAINTES'
  169.       mchel5.ifoche=ifour
  170.  
  171.       do isous=1,n1
  172.       imodel=mmode2.kmodel(isous)
  173.       segact imodel
  174.       mchel5.imache(isous)=imamod
  175.       mchel5.conche(isous)=conmod
  176.       mele= nefmod
  177. C====
  178.       CALL ELQUOI(MELE,0,5,INFO,IMODEL)
  179.          IF (IERR.NE.0) THEN
  180.             SEGSUP MCHEL5
  181.             IRET=0
  182.             RETURN
  183.       ENDIF
  184.        MINTE= INFELL(11)
  185. C====
  186.          MCHEL5.INFCHE(ISOUS,1)=0
  187.          MCHEL5.INFCHE(ISOUS,2)=0
  188.          MCHEL5.INFCHE(ISOUS,3)=NIFOUR
  189.          MCHEL5.INFCHE(ISOUS,4)=MINTE
  190.          MCHEL5.INFCHE(ISOUS,5)=0
  191.          MCHEL5.INFCHE(ISOUS,6)=5
  192. C
  193.       N2=1
  194.       segini MCHAML
  195.       mchel5.ichaml(isous)=mchaml
  196.       nomche(1)='EFFX'
  197.       typche(1)='REAL*8'
  198.       meleme = imamod
  199.       segact meleme
  200.       N1PTEL=2
  201.       N1EL= NUM(/2)
  202.       N2PTEL=0
  203.       N2EL=0
  204.       segini melval
  205.       ielval(1)=melval
  206.       segsup info
  207.       enddo
  208.  
  209. C
  210. C ----FIN    PREPARATION DU CHAMP DE SORTIE----------------
  211. C
  212.  
  213. C   preparation de la recherche du sens d application des tensions
  214.        if(iran.ne.0)  then
  215.        ipt3=iran
  216.          segact ipt3
  217.          segini , ipt5=ipt3
  218.          ity=ipt5.itypel
  219.            if(ity.ne.1) then
  220.           call change(iran,1)
  221.           ipt5=iran
  222.           segact ipt5
  223.            endif
  224.        else
  225.        ipt5=0
  226.        endif
  227.  
  228. C
  229. C --- CREATION DU MCHELM  DE CALCUL  ( TEMPORAIRE )
  230. C
  231.       N1=NSOUS
  232. C
  233.       L1=11
  234.       N3=6
  235.       SEGINI MCHELM
  236.       TITCHE='CONTRAINTES'
  237.       IFOCHE=IFOUR
  238.       IPSTRS=MCHELM
  239. C
  240. C --- DEBUT DE LA BOUCLE SUR LES DIFFERENTES ZONES
  241. C
  242.       call oooeta(mcoord,ieta,imod)
  243.        if (ieta.ne.1) then
  244.            ieta=1
  245.            segact mcoord
  246.        endif
  247.       DO 500 ISOUS=1,NSOUS
  248. C
  249. C --- INITIALISATION
  250. C
  251.          MOCARA=0
  252.          MOMATR=0
  253.          MOSTRS=0
  254.          IVAMAT=0
  255.          IVACAR=0
  256.          IVASTR=0
  257.          ivasi0=0
  258. c
  259.          ipeffx = 0
  260.          ipyoun = 0
  261.          ipff   = 0
  262.          ipphif = 0
  263.          ipganc = 0
  264.          iprmu0 = 0
  265.          ipfprg = 0
  266.          iprh10 = 0
  267.          ipsect = 0
  268.          IMODEL=KMODEL(ISOUS)
  269.          SEGACT IMODEL
  270.          IPMAIL=IMAMOD
  271.          CONM  =CONMOD
  272.          IMACHE(ISOUS)=IPMAIL
  273.          CONCHE(ISOUS)=CONMOD
  274. C ...... TRAITEMENT DU MODELE ......................
  275.  
  276.          MELE=NEFMOD
  277.          MELEME=IPMAIL
  278. C
  279.          segact meleme
  280.          itpdel=itypel
  281.  
  282.          NFOR=FORMOD(/2)
  283.          NMAT=MATMOD(/2)
  284. C
  285. C    ... on ordonne le nouveau maillage (si on en a besoin) ...
  286.  
  287.          NBELEM = NUM(/2)
  288. C=====================================================
  289.          if(itypel.eq.1) then
  290.  
  291. c    ... dans le cas des POI1 on n'a rien a orienter ...
  292. C=====================================================
  293.          elseif(itypel.eq.2) then
  294.  
  295.            ipt2 = meleme
  296.             NOD1=NUM(1,1)
  297.             NOD2=NUM(2,num(/2))
  298. C           write(6,*)' noeuds extremites ' ,nod1,nod2
  299.          if(iran.eq.0)  then
  300.          ivers=.false.
  301.          else
  302.             do ino=1,ipt5.num(/2)
  303.             inodi=ipt5.num(1,ino)
  304.             if(inodi.eq.nod1) then
  305.             ivers=.false.
  306.             goto  364
  307.             elseif(inodi.eq.nod2) then
  308.             ivers=.true.
  309.             goto  364
  310.             endif
  311.             enddo
  312. C     aucune des extremites n'est  dans la liste
  313. C           on applique au point final
  314.             call erreur(833)
  315. C  attention aux desactivations
  316.             SEGSUP MCHELM
  317.             IRET=0
  318.             RETURN
  319. 364      continue
  320.          endif
  321. C=====================================================
  322. c    ... si ITYPEL différent de 1 et 2 ...
  323.          else
  324.  
  325. c       ... Type d''élément support inconnu !!!!!!!!!!! ............
  326.             SEGSUP MCHELM
  327.             IRET=0
  328.             RETURN
  329.             call erreur(19)
  330.             return
  331.          endif
  332. C=====================================================
  333. C
  334. C ------ NATURE DU MATERIAU
  335. C
  336.          CALL NOMATE(FORMOD,NFOR,MATMOD,NMAT,CMATE,MATE,INAT)
  337.          IF (CMATE.EQ.' ') THEN
  338.           CALL ERREUR(251)
  339.             SEGSUP MCHELM
  340.             IRET=0
  341.             RETURN
  342.           ENDIF
  343. C
  344. C ------ INFORMATION SUR L'ELEMENT FINI
  345. C
  346.          CALL ELQUOI(MELE,0,5,IPINF,IMODEL)
  347.          IF (IERR.NE.0) THEN
  348.             SEGSUP MCHELM
  349.             IRET=0
  350.             RETURN
  351.          ENDIF
  352. C
  353.          INFO=IPINF
  354.          MELE =INFELL(1)
  355.          MFR  =INFELL(13)
  356.          MINTE=INFELL(11)
  357.          NSTRS=INFELL(16)
  358.  
  359.          SEGINI WRK3
  360. C
  361. C ------ CREATION DU TABLEAU INFOS
  362. C
  363.          IPCHE1=0
  364.          IPCHE2=0
  365. c ... ATTENTION !!! A quoi ca sert ????????????? .......
  366.         CALL IDENT(IPMAIL,CONM,IPCHE1,IPCHE2,INFOS,IRTD)
  367.          IF (IRTD.EQ.0) THEN
  368.             SEGSUP MCHELM
  369.             SEGSUP INFO
  370.              IRET=0
  371.             RETURN
  372.          ENDIF
  373.  
  374. c    ... Remplissage du tableau INFCHE du MCHELM ...
  375.  
  376.          INFCHE(ISOUS,1)=0
  377.          INFCHE(ISOUS,2)=0
  378.          INFCHE(ISOUS,3)=NIFOUR
  379.          INFCHE(ISOUS,4)=MINTE
  380.          INFCHE(ISOUS,5)=0
  381.          INFCHE(ISOUS,6)=5
  382. C
  383. C ------ LECTURE de NBPGAU dans
  384. C
  385.          SEGACT MINTE
  386.          NBPGAU=POIGAU(/1)
  387. C
  388. C ------ RECHERCHE DES NOMS DE COMPOSANTES des contraintes ...
  389. C
  390.          if(lnomid(4).ne.0) then
  391.           nomid=lnomid(4)
  392.           segact nomid
  393.           mostrs=nomid
  394.           nstr=lesobl(/2)
  395.           nfac=lesfac(/2)
  396.           lsupco=.false.
  397.         else
  398.          lsupco=.true.
  399.          CALL IDCONT(IMODEL,IFOUR,MOSTRS,NSTR,NFAC)
  400.         endif
  401.          nomid = mostrs
  402.          segact,nomid
  403.          call place(lesobl,nstr,ipeffx,'EFFX')
  404. C
  405. C ------ VERIFICATION DE LEUR PRESENCE et ...
  406. C ------ RECUPERATION DES CONTRAINTES INITIALES (s'il le faut) ...
  407. C
  408.          if(IPCHC1.gt.0) then
  409.  
  410.             NBTYPE=INFELL(16)
  411.             SEGINI NOTYPE
  412.             MOTYPE=NOTYPE
  413.             TYPE(1)='REAL*8'
  414.  
  415.             call komcha(ipchc1,ipmail,conm,mostrs,motype,1,infos,3,
  416.      &                  ivasi0)
  417.  
  418.             IF (IERR.NE.0) THEN
  419.                NSTRS=0
  420.                GOTO 9990
  421.             ENDIF
  422.  
  423.             SEGSUP NOTYPE
  424.          endif
  425.  
  426. C
  427. C ------ CREATION DU MCHAML DE LA SOUS ZONE
  428. C
  429.          N2=NSTRS
  430.          SEGINI MCHAML
  431.          ICHAML(ISOUS)=MCHAML
  432.  
  433.          NSR=1
  434.          NCOSOR=NSTRS
  435.          SEGINI MPTVAL
  436.          IVASTR=MPTVAL
  437. C
  438. C ------ RECHERCHE DE LA TAILLE DES MELVAL A ALLOUER
  439. C
  440.          N1PTEL=nbpgau
  441.          N1EL=NBELEM
  442.          N2PTEL=0
  443.          N2EL=0
  444. C
  445.          NOMID=MOSTRS
  446.          SEGACT NOMID
  447. C    initialisation   des melvals  des precontaintes
  448.          DO 100 ICOMP=1,NSTRS
  449.             NOMCHE(ICOMP)=LESOBL(ICOMP)
  450.             TYPCHE(ICOMP)='REAL*8'
  451.             SEGINI MELVAL
  452.             IELVAL(ICOMP)=MELVAL
  453.             IVAL(ICOMP)=MELVAL
  454.  100     CONTINUE
  455. C            ipacara  preparé dans splitag
  456. C  recuperation du module de young  dans ivamat ( position 1 )
  457. C  recuperation de la section       dans ivacar ( position 2 )
  458.          mchel6=ipcara
  459.          segact mchel6
  460.          mcham6=mchel6.ichaml(isous)
  461.          segact mcham6
  462.          melval= mcham6.ielval(1)
  463.          segact melval
  464.          ivamat=melval
  465.          melval= mcham6.ielval(2)
  466.          segact melval
  467.          ivacar=melval
  468.  
  469. C
  470. C ------ CALCUL lui même --------------------------
  471. C
  472. C ... cas du CERC ...
  473. C
  474.  
  475.       IF(MELE.EQ.95) THEN
  476.          DO 2004 IB=1,NBELEM
  477.             MPTVAL=IVASTR
  478.             DO 1701 I=1,NSTRS
  479.                MELVAL=IVAL(I)
  480.                VELCHE(1,IB)=PS1
  481.  1701       CONTINUE
  482.          iret=1
  483.  2004    CONTINUE
  484.  
  485.        ELSE
  486. C  ----------------valeurs par defaut ----------------------
  487.  
  488.       f1 = 0.18d0
  489.       f2 = 0.002d0
  490.       ganc = 0.d0
  491.       rmu0 = 0.43d0
  492.       fprg = 1.7d9
  493.       rh10 = 2.5d0
  494.  
  495.       if(iptab.ne.0) then
  496.       MTABLE = iptab
  497.       segact mtable
  498.       nbind= mlotab
  499.           IVALIN = 0
  500.          XVALIN = 0.D0
  501.          LOGIN  = .TRUE.
  502.          IOBIN  = 0
  503.          TAPIND = 'MOT'
  504. C
  505.       TYPOBJ = '        '
  506.       Call acctab(MTABLE,'MOT',IVALIN,XVALIN,'FF',LOGIN,0,
  507.      &             TYPOBJ,IVELRE,XVALRE,CHARRE,LOGRE,IOBRE)
  508.       IF(TYPOBJ.EQ.'FLOTTANT') F1 = XVALRE
  509.       IF(TYPOBJ.EQ.'ENTIER  ') F1 = IVELRE
  510.       TYPOBJ = '        '
  511.       Call acctab(MTABLE,'MOT',IVALIN,XVALIN,'PHIF',LOGIN,0,
  512.      &             TYPOBJ,IVELRE,XVALRE,CHARRE,LOGRE,IOBRE)
  513.       IF(TYPOBJ.EQ.'FLOTTANT') F2 = XVALRE
  514.       IF(TYPOBJ.EQ.'ENTIER  ') F2 = IVELRE
  515.       TYPOBJ = '        '
  516.       Call acctab(MTABLE,'MOT',IVALIN,XVALIN,'GANC',LOGIN,0,
  517.      &             TYPOBJ,IVELRE,XVALRE,CHARRE,LOGRE,IOBRE)
  518.       IF(TYPOBJ.EQ.'FLOTTANT') GANC = XVALRE
  519.       IF(TYPOBJ.EQ.'ENTIER  ') GANC = IVELRE
  520.       TYPOBJ = '        '
  521.       Call acctab(MTABLE,'MOT',IVALIN,XVALIN,'RMU0',LOGIN,0,
  522.      &             TYPOBJ,IVELRE,XVALRE,CHARRE,LOGRE,IOBRE)
  523.       IF(TYPOBJ.EQ.'FLOTTANT') RMU0 = XVALRE
  524.       IF(TYPOBJ.EQ.'ENTIER  ') RMU0 = IVELRE
  525.       TYPOBJ = '        '
  526.       Call acctab(MTABLE,'MOT',IVALIN,XVALIN,'FPRG',LOGIN,0,
  527.      &             TYPOBJ,IVELRE,XVALRE,CHARRE,LOGRE,IOBRE)
  528.       IF(TYPOBJ.EQ.'FLOTTANT') FPRG = XVALRE
  529.       IF(TYPOBJ.EQ.'ENTIER  ') RPRG = IVELRE
  530.       TYPOBJ = '        '
  531.       Call acctab(MTABLE,'MOT',IVALIN,XVALIN,'RH10',LOGIN,0,
  532.      &             TYPOBJ,IVELRE,XVALRE,CHARRE,LOGRE,IOBRE)
  533.       IF(TYPOBJ.EQ.'FLOTTANT') RH10 = XVALRE
  534.       IF(TYPOBJ.EQ.'ENTIER  ') RH10 = IVELRE
  535.       segdes MTABLE
  536.  
  537. 4321  format(6E12.5)
  538.        endif
  539. C----------------------------------------------------------
  540. C       write(6,4321) f1,f2,ganc,rmu0,fprg,rh10
  541. C
  542. C ... Boucle sur les éléments dans lesquels on fait le calcul ...
  543. C
  544.  
  545.          SLON=0.0D0
  546.          FAI=0.0D0
  547.          ids= idim+1
  548.             segact,minte
  549.             if(nbpgau.ne.2) then
  550.                call erreur(5)
  551.                return
  552.             endif
  553. *
  554. *    NEW   1-ERE BOUCLE SUR LES ELEMENTS POUR TROUVER
  555. *     LA LONGUEUR D'INFLUENCE DU RECUL D'ANCRAGE
  556. *
  557.             NAM =NBELEM
  558.             SEGINI ALTRAV
  559.             IAM=0
  560.             SUM =0.D0
  561. *****       PRINT *,'IVERS=',IVERS
  562. *
  563.             DO 5005 IC=1,NBELEM
  564.  
  565.             if(ivers) then
  566. C    ordre inverse
  567.                JC= NBELEM+1-IC
  568.                NC3=ipt2.NUM(1,JC)
  569.                NC2=ipt2.NUM(2,JC)
  570.               IF(JC.EQ.NBELEM) THEN
  571.                  NC1=NC2
  572.               ELSE
  573.                  NC1=IPT2.NUM(2,JC+1)
  574.               ENDIF
  575.             else
  576. C    ordre normal
  577.                JC=IC
  578.                NC3=ipt2.NUM(2,JC)
  579.                NC2=ipt2.NUM(1,JC)
  580.                IF(JC.EQ.1) THEN
  581.                   NC1=NC2
  582.                ELSE
  583.                   NC1=IPT2.NUM(1,JC-1)
  584.                ENDIF
  585.             endif
  586. *****         PRINT *,'NC1=',NC1,' NC2=',NC2,' NC3=',NC3
  587.               JS1=(NC1-1)*IDS
  588.               JS2=(NC2-1)*IDS
  589.               JS3=(NC3-1)*IDS
  590.  
  591.                DO     IW=1,IDIM
  592.                   X1(IW)=XCOOR(JS1+IW)
  593.                   X2(IW)=XCOOR(JS2+IW)
  594.                   X3(IW)=XCOOR(JS3+IW)
  595.                enddo
  596. C          --- Distance entre points
  597.                DS1=0.0D0
  598.                DS2=0.0D0
  599.                DO     IW=1,IDIM
  600.                   RL(IW)=X3(IW)-X2(IW)
  601.                   RS(IW)=X2(IW)-X1(IW)
  602.                   DS1=DS1+RL(IW)*RL(IW)
  603.                   DS2=DS2+RS(IW)*RS(IW)
  604.                enddo
  605.                DS1=SQRT(DS1)
  606.                DS2=SQRT(DS2)
  607.                CS1=0.0D0
  608.                IF(DS2.NE.0.D0) THEN
  609.                   DO     IW=1,IDIM
  610.                      RL(IW)=RL(IW)/DS1
  611.                      RS(IW)=RS(IW)/DS2
  612.                      CS1=CS1+RL(IW)*RS(IW)
  613.                   enddo
  614.                ELSE
  615.                      CS1=1.D0
  616.                ENDIF
  617.                IF(CS1.GT.1.0) CS1=1.0D0
  618.                ALFA=ACOS(CS1)
  619.                IAM=IAM+1
  620.                ANG(IAM)=FAI
  621.                DANG(IAM)=ALFA
  622.                ACUR(IAM)=SLON
  623.                DACUR(IAM)=DS1
  624.                SLON=SLON+DS1
  625.                FAI=FAI+ALFA
  626.  
  627. *          PRINT *,'IAM=',IAM, ' ACUR =',ACUR(IAM),
  628. *     &    ' ANG=',ANG(IAM),' DACUR=',DACUR(IAM),' DANG=',DANG(IAM)
  629. *
  630.                IF(SLON.NE.0.D0) THEN
  631.                 SUM = SUM + FAI/SLON
  632.                ENDIF
  633. *
  634.  5005       CONTINUE
  635.             ANG(NAM+1)=FAI
  636.             ACUR(NAM+1)=SLON
  637. *
  638.             SLONT=SLON
  639.             FAIT=FAI
  640.             XLMBDA=0.D0
  641. *
  642.             IF(GANC.EQ.0.D0) GO TO 7999
  643. *
  644. *------------------------------------------------------------
  645. *    CALCUL DE LA LONGUEUR D'INFLUENCE DU RECUL D'ANCRAGE
  646. *              LE CAS ECHEANT  ( GANC > 0 )
  647. *------------------------------------------------------------
  648. *
  649. *      ON PREND LES MODULE ET SECTION DU 1-ER ELEMENT
  650. *      POUR LA 1-ERE APPROXIMATION DE XLAMBDA
  651. *
  652.             MELVAL =IVAMAT
  653.             EA = VELCHE(1,1)
  654.             MELVAL =IVACAR
  655.             SECT1 = VELCHE(1,1)
  656. *        write(6,*) ' f1, f2, slont, ps1, sect1, ea,ganc,sum,nbelem,fai'
  657. *       write(6,*)  f1, f2, slont, ps1, sect1, ea,ganc,sum,nbelem,fai
  658.             IF(F1.EQ.0.D0.OR.FAI.LE.1.D-4) THEN
  659.               IF(F2.EQ.0.D0) THEN
  660.                 XLMBDA=SLONT
  661.               ELSE
  662.                 aaa= F2*GANC*EA*SECT1/PS1
  663.                  bbb = aaa ** 0.5
  664.                  ccc = 1.-bbb
  665.                  if( ccc .le.0.) then
  666.                      call erreur(768)
  667.                      return
  668.                  else
  669. *               write(6,*) ' aaa,bbb,ccc',aaa,bbb,ccc
  670.                     XLMBDA = -(LOG(1.D0-
  671.      &                     (F2*GANC*EA*SECT1/PS1)**0.5))/F2
  672.                  endif
  673.               ENDIF
  674.             ELSE
  675.               ALSURL=SUM/FLOAT(NBELEM)
  676.               PP= F1*ALSURL + F2
  677.               IF(PP.EQ.0.D0) THEN
  678.                 XLMBDA=SLONT
  679.               ELSE
  680.                 XLMBDA = SQRT( (GANC*EA*SECT1)/(PS1*PP) )
  681.               ENDIF
  682. *          PRINT *,'ALSURL=',ALSURL
  683. *          PRINT *,'F2=',F2, 'PP=',PP
  684.             ENDIF
  685. *
  686. *    TEST SI XLAMBBDA > SLON
  687. *
  688.             IF(XLMBDA.GT.SLONT) XLMBDA=SLONT
  689.  
  690. *       PRINT *,'INITIALISATION DE XLMBDA=',XLMBDA
  691. *
  692.             IF(XLMBDA.LT.0.D0) THEN
  693. *****           PRINT *,' ATTENTION XLAMBDA INITIAL EST < 0 '
  694.                 CALL ERREUR(460)
  695.                 SEGSUP ALTRAV
  696.                 GO TO 9990
  697.             ENDIF
  698. *
  699.             ITER=0
  700.             LAST=0
  701.             IDICH=0
  702. *
  703. *         ITERATIONS
  704. *
  705.  8000       CONTINUE
  706.             ITER=ITER+1
  707. *       PRINT *,'ITER = ',ITER
  708.             IF(ITER.GT.ITEMAX) THEN
  709.                 CALL ERREUR(460)
  710.                 SEGSUP ALTRAV
  711.                 GO TO 9990
  712.             ENDIF
  713.             PSL=0.D0
  714.             SLON=0.0D0
  715.             FAI=0.0D0
  716.             IAM=0
  717. *
  718.             IF(LAST.EQ.0) THEN
  719. *         write(6,*) ' nam + 1',nam+1
  720. *         write(6,*) ' xlmbda acur(1)  acur(2) ang(1) ang(2) '
  721. *          write(6,*)  xlmbda,acur(1), acur(2) ,ang(1),ang(2)
  722. *             CALCUL DE L'ANGLE XALFA ASSOCIE A XLMBDA
  723.               CALL PRECN0(XLMBDA,ACUR,ANG,NAM+1,ANGBDA,ISUC)
  724.               IF(ISUC.EQ.NAM+1) ISUC=NAM
  725. *
  726. *****         PRINT *,'ISUC=',ISUC
  727.               IF(ISUC.LT.0) THEN
  728.                  CALL ERREUR(-ISUC)
  729.                  SEGSUP ALTRAV
  730.                  GO TO 9990
  731.               ENDIF
  732.  
  733.               SD=F1*ANGBDA + F2*XLMBDA
  734.               DSDDD=F1*DANG(ISUC)/DACUR(ISUC)+F2
  735. *
  736.             ELSE
  737.               SD=F1*FAIT + F2 * SLONT
  738.             ENDIF
  739. *
  740.             GG=0.D0
  741.             DGDD=0.D0
  742.             DO 8005 IC=1,NBELEM
  743.                if(ivers) then
  744.                   JC= NBELEM+1-IC
  745.                else
  746.                  JC=IC
  747.                endif
  748. *
  749.                IAM=IAM+1
  750.                ALFA = DANG(IAM)
  751.                DS1  = DACUR(IAM)
  752.  
  753.                do 8006  iptg=1,nbpgau
  754.  
  755.                if(ivers) then
  756.                ig= nbpgau+1-iptg
  757.                else
  758.                ig = iptg
  759.                endif
  760.  
  761. *         ---- Récupération des caractéristiques matérielles et géométriques ----
  762.  
  763.                melval =ivamat
  764.                ea = velche(min(ig,velche(/1)),
  765.      &                       min(JC,velche(/2)))
  766.  
  767.                melval =ivacar
  768.                sect1 = velche(min(ig,velche(/1)),
  769.      &                        min(JC,velche(/2)))
  770.  
  771.                SABS=SLON + DS1*(1+QSIGAU(IPTG))/2
  772.                IF(SABS.GT.XLMBDA) GO TO 8007
  773.  
  774.                SP=F1*FAI + F2*SABS
  775.                IF(LAST.EQ.0) THEN
  776.                  GG = GG +
  777.      &           PS1 * (EXP(-SP) - EXP(SP-2.D0*SD))*DS1*0.5D0
  778.      &                                   /(EA*SECT1)
  779.  
  780.                  DGDD = DGDD +
  781.      &           PS1 * (EXP(SP-2.D0*SD))*DSDDD*DS1/(EA*SECT1)
  782.                ELSE
  783.                  GG = GG + PS1 * (EXP(-SP))*DS1*0.5D0
  784.      &                                   /(EA*SECT1)
  785.                  DGDD = DGDD + (EXP(SP-SD))*DS1*0.5D0
  786.      &                                   /(EA*SECT1)
  787.                ENDIF
  788.  
  789.  8006          continue
  790.                SLON=SLON+DS1
  791.                FAI=FAI+ALFA
  792.  8005       CONTINUE
  793. *
  794. *        TEST DE CONVERGENCE
  795. *
  796.  8007       CONTINUE
  797. *****       PRINT *,'GG=',GG, ' GANC=',GANC,' FG =',GG-GANC
  798. *****       PRINT *,'DGDD=', DGDD
  799.  
  800.             FGG(ITER)=GG-GANC
  801.             XLAM(ITER)=XLMBDA
  802.  
  803. *
  804. *        CAS OU L'EFFET DE L'ANCRAGE DEPASSE LA LONGUEUR
  805. *        DU CABLE
  806. *
  807. *****    PRINT *,'LAST=',LAST
  808.          IF(LAST.EQ.1) THEN
  809.  
  810.              IF(ABS(GG-GANC)/GANC.LT.1.D-3) THEN
  811.                 PSL=0.D0
  812.              ELSE
  813.                 PSL=(GG-GANC)/DGDD
  814.              ENDIF
  815.  
  816. *****        PRINT *,'PSL=',PSL
  817.              IF(PSL.GT.PS1*EXP(-SD)) THEN
  818. *****           PRINT *,' ATTENTION PSL EST TROP GRAND '
  819.              ENDIF
  820.              GO TO 7999
  821.          ENDIF
  822.  
  823.  
  824.          IF(ABS(GG-GANC)/GANC.LT.1.D-3) GO TO 7999
  825. *
  826.          IF(ITER.GE.6 .AND.IDICH.EQ.0.AND.
  827.      &                     FGG(ITER)*FGG(ITER-1).LT.0.D0) THEN
  828.             IDICH=1
  829.             XLAM1=XLAM(ITER-1)
  830.             XLAM2=XLAM(ITER)
  831.             FG1=FGG(ITER-1)
  832.             FG2=FGG(ITER)
  833.          ENDIF
  834. *****    PRINT *,'IDICH=',IDICH
  835. *
  836.          IF(IDICH.EQ.0) THEN
  837.             XLMBDA=XLMBDA - (GG-GANC)/DGDD
  838.          ELSE
  839. *****       PRINT *,'XLAM(ITER)=',XLAM(ITER)
  840. *****       PRINT *, '  FGG(ITER)=',FGG(ITER)
  841.  
  842.             IF(FG1*FGG(ITER).LT.0.D0) THEN
  843.                XLAM2=XLAM(ITER)
  844.                FG2=FGG(ITER)
  845.             ELSE IF (FG1*FGG(ITER).GT.0.D0) THEN
  846.                XLAM1=XLAM(ITER)
  847.                FG1=FGG(ITER)
  848.             ENDIF
  849. *
  850. *****       PRINT *,'XLAM1=',XLAM1,'  FG1=',FG1
  851. *****       PRINT *,'XLAM2=',XLAM2,'  FG2=',FG2
  852.             DLAM=XLAM2-XLAM1
  853.             PROF=FG1/(FG1-FG2)
  854.             XLMBDA=XLAM1+PROF*DLAM
  855.          ENDIF
  856. *
  857. *****      PRINT *,'NOUVELLE XLMBDA = ', XLMBDA
  858.  
  859.            IF(XLMBDA.GE.SLONT) THEN
  860.               IF(ITER.EQ.ITEMAX-1) THEN
  861.                  LAST=1
  862.               ELSE
  863.                  XLMBDA=SLONT
  864.               ENDIF
  865.            ENDIF
  866. *
  867. *          si on trouve XLMBDA < 0 , on reinitialise
  868. *
  869.            IF(XLMBDA.LT.0.D0) THEN
  870.                 IF(ITER.GE.2) THEN
  871.                    XLMBDA=(XLAM(ITER)+XLAM(1))/ITER
  872.                 ELSE
  873.                    XLMBDA=XLAM(1)/2.D0
  874.                 ENDIF
  875. *****           PRINT *,'XLMBDA  REINITIALISEE = ', XLMBDA
  876.            ENDIF
  877. *
  878. *          SORTIE SI IDICH=1 ET XLMBDA CONVERGE
  879. *
  880.            IF(IDICH.EQ.1.AND.ITER.EQ.ITEMAX) THEN
  881.               IF(XLMBDA.NE.0.D0) THEN
  882.                  IF(ABS(XLMBDA-XLAM(ITER))/XLMBDA.
  883.      &                 LT.1.D-3) GO TO 7999
  884.               ENDIF
  885.            ENDIF
  886. *
  887.            GO TO 8000
  888. *
  889. *    2-EME BOUCLE SUR LES ELEMENTS
  890. *
  891.  7999       CONTINUE
  892.             SLON=0.0D0
  893.             FAI=0.0D0
  894.             IAM=0
  895.             DO 4005 IC=1,NBELEM
  896.  
  897.             if(ivers) then
  898. C    ordre inverse
  899.                JC= NBELEM+1-IC
  900.             else
  901. C    ordre normal
  902.               JC=IC
  903.             endif
  904.  
  905. C     ---- ON RECUPERE LES GRANDEURS GEOMETRIQUES -----------
  906.  
  907.                IAM=IAM+1
  908.                ALFA = DANG(IAM)
  909.                DS1  = DACUR(IAM)
  910.  
  911.                do 12345 iptg=1,nbpgau
  912.  
  913.                if(ivers) then
  914.                ig= nbpgau+1-iptg
  915.                else
  916.                ig = iptg
  917.                endif
  918.  
  919. C         ---- Récupération des EFFX deja présents -----------
  920.  
  921.                IF (IPCHC1.EQ.0) THEN
  922.                   SINI = 0.D0
  923.                ELSE
  924.                   mptval = ivasi0
  925.                   MELVAL = IVAL(ipeffx)
  926.                   SINI =velche(min(ig,velche(/1)),
  927.      &                         min(JC,velche(/2)))
  928.                ENDIF
  929.  
  930. C         ---- Récupération des caractéristiques matérielles et géométriques ----
  931.  
  932.                melval =ivamat
  933.                ea = velche(min(ig,velche(/1)),
  934.      &                       min(JC,velche(/2)))
  935.  
  936.                melval =ivacar
  937.                sect1 = velche(min(ig,velche(/1)),
  938.      &                        min(JC,velche(/2)))
  939.  
  940. C          --- Calcul des pertes de précontrainte -------------------
  941.  
  942.                ICOMP = 0
  943.                call precn1(PS1,ea,f1,f2,ganc,RMU0,FPRG,RH10,
  944.      &                     SLON+DS1*(1+qsigau(iptg))/2,FAI,SINI,
  945.      &                     SECT1,XLMBDA,SD,LAST,PSL,PSOUT,ICOMP)
  946.                IF (ICOMP.EQ.1) THEN
  947.                   CALL ERREUR(768)
  948.                   goto 9990
  949.                   RETURN
  950.                ENDIF
  951.  
  952. C          --- On range les résultats
  953.  
  954.                iret=1
  955.                MPTVAL=IVASTR
  956.                DO 1700 I=1,NSTRS
  957.                   MELVAL=IVAL(I)
  958.                   VELCHE(ig,JC)=PSOUT
  959.  1700          CONTINUE
  960.  
  961. 12345          continue
  962.                SLON=SLON+DS1
  963.                FAI=FAI+ALFA
  964. C            write(6,*) 'slon et fai ',slon,fai
  965.  
  966.  4005       CONTINUE
  967.             SEGSUP ALTRAV
  968. C----------------------    endif  sur  CERC ou  ELEMENTS LINEAIRE
  969.       END IF
  970.  
  971. C
  972.  9990 CONTINUE
  973.  
  974.       SEGSUP WRK3
  975.       IF(IERR.NE.0)THEN
  976.          SEGSUP MCHAML
  977.       ENDIF
  978. c
  979. C      CALL DTMVAL(IVACAR,1)
  980.  
  981.       if(IPCHC1.gt.0) call dtmval(ivasi0,1)
  982. C
  983.       IF(IERR.NE.0)THEN
  984.           CALL DTMVAL(IVASTR,3)
  985.       ELSE
  986.           CALL DTMVAL(IVASTR,1)
  987.       ENDIF
  988. C
  989.       IF(MOMATR.NE.0)THEN
  990.          NOMID=MOMATR
  991.          SEGSUP NOMID
  992.       ENDIF
  993.       IF(MOCARA.NE.0)THEN
  994.          NOMID=MOCARA
  995.          SEGSUP NOMID
  996.       ENDIF
  997.       IF(MOSTRS.NE.0)THEN
  998.          NOMID=MOSTRS
  999.          if(lsupco)SEGSUP NOMID
  1000.       ENDIF
  1001. C
  1002.       SEGSUP INFO
  1003.       if(iret.eq.0) return
  1004. C
  1005.   500 CONTINUE
  1006. c-------- Fin de la boucle sur les zones du MODELE  TEMPORAIRE
  1007. C
  1008. C =========== rangement dans le MCHELM  final   MCHEL5
  1009.       if(ifomod.ne.0.and.ifomod.ne.1) then
  1010. C           les cables standards
  1011.       segact siezo
  1012. C      segact mchel5
  1013.       MCHELM = IPSTRS
  1014.       segact MCHELM
  1015. C
  1016. C  boucle sur les zones du modele ou champ initial
  1017.  
  1018. C
  1019.         ibc= 0
  1020.       do 3010 isous=1,mchel5.imache(/1)
  1021.       sielc = iezon(isous)
  1022.       segact sielc
  1023. C
  1024.         mcham5= mchel5.ichaml(isous)
  1025.         segact mcham5
  1026.         melva5 =mcham5.ielval(1)
  1027.         segact melva5*mod
  1028. C  boucle sur les partitions provisoires de la sz
  1029. C
  1030.       inel = 0
  1031.         do 3020 ik=1,nbcz
  1032.         ibc=ibc+1
  1033.         mchaml=ichaml(ibc)
  1034.         segact mchaml
  1035.         melval=ielval(1)
  1036.         segact melval
  1037. C
  1038.       do iii=1,velche(/2)
  1039.       inel = inel+1
  1040.       iel=isens(1,inel)
  1041.         if(isens(2,inel).eq.1) then
  1042.       do ip=1,velche(/1)
  1043.       melva5.velche(ip,iel)=velche(ip,iii)
  1044.       enddo
  1045.         else
  1046.       do ip=1,velche(/1)
  1047.       melva5.velche(3-ip,iel)=velche(ip,iii)
  1048.       enddo
  1049.         endif
  1050.       enddo
  1051. 3020  continue
  1052.       segsup sielc
  1053. 3010  continue
  1054.       call DTCHAM(ipcar2)
  1055.       if(ipchc2.ne.0) call DTCHAM(ipchc2)
  1056.       segsup ipt5
  1057.       else
  1058. C          les cerc
  1059.       icha1 = MCHEL5
  1060.       MCHEL5 = MCHELM
  1061.       MCHELM = icha1
  1062.       endif
  1063. 333   format(i4,2E12.5)
  1064. C  maintenant les destructions d objets temporaires
  1065.       call DTCHAM(MCHELM)
  1066. C      call DTMODL(ipmodt)
  1067.       mmodel = ipmodt
  1068.       segsup mmodel
  1069. C
  1070.       IPSTRS =MCHEL5
  1071.       iret = 1
  1072.  
  1073.       RETURN
  1074.       END
  1075.  
  1076.  
  1077.  

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