Télécharger la plaquette Cast3M


fr en

Télécharger bsigm3.eso

Retour à la liste

Numérotation des lignes : 

bsigm3

  1. C BSIGM3    SOURCE    OF166741  25/02/21    21:15:13     12166          
  2.       SUBROUTINE BSIGM3(IPMAIL,LRE,NSTRS,LW,IVACAR,NCARR,IVECT,
  3.      &  MELE,CMATE,IVASTR,ISOUS,NBPGAU,NBPTEL,IPMINT,NFOR,IVAFOR
  4.      &   ,ADPG,BDPG,CDPG,IIPDPG,IVAMAT,NMATT,MFR,dcmate)
  5. *----------------------------------------------------------------------
  6. *        _______________________________                              *
  7. *        |                              |                             *
  8. *        |  calcul des forces aux noeuds|                             *
  9. *        |______________________________|                             *
  10. *                                                                     *
  11. *      poutre,tuyau,linespring,tuyau fissure,barre,cerce,tuyo         *
  12. *      poutre de timoschenko,shb8,joi1,zone_cohesive
  13. *                                                                     *
  14. *---------------------------------------------------------------------*
  15. *                                                                     *
  16. *   entrees :                                                         *
  17. *   ________                                                          *
  18. *                                                                     *
  19. *        ipmail   pointeur sur un segment  meleme                     *
  20. *        lre      nombre de ddl dans la matrice de rigidite           *
  21. *        nstrs    nombre de composante de contraintes/deformations    *
  22. *        lw       dimension du tableau de travail de l'element        *
  23. *        ivacar   pointeur sur les chamelems de caracteristiques géo. *
  24. *        ncarr    nombre de caracteristiques geometriques             *
  25. *        ivamat   pointeur sur les chamelems de caracteristiques mat. *
  26. *        nmatt    nombre de caracteristiques matériau                 *
  27. *        ivect    flag indiquant si on a entree les axes locaux       *
  28. *        mele     numero de l'element fini                            *
  29. *        ivastr   pointeur sur un segment mptval contenant les        *
  30. *                 les melvals de contraints                           *
  31. *        isous     numero de la sous-zone                             *
  32. *        nbpgau   nombre de points d'integration pour les contraintes *
  33. *        nbptel   nombre de points par element                        *
  34. *        ipmint   pointeur sur un segment minte                       *
  35. *        nfor     nombre de composantes de forces                     *
  36. *                                                                     *
  37. *   sorties :                                                         *
  38. *   ________                                                          *
  39. *                                                                     *
  40. *        ivafor   pointeur sur un segment mptval contenant les        *
  41. *                 les melvals de forces                               *
  42.  
  43. *                                                                     *
  44. *---------------------------------------------------------------------*
  45.       IMPLICIT INTEGER(I-N)
  46.       IMPLICIT REAL*8(A-H,O-Z)
  47.  
  48. -INC PPARAM
  49. -INC CCOPTIO
  50. -INC CCHAMP
  51. -INC CCREEL
  52.  
  53. -INC SMCHAML
  54. -INC SMCHPOI
  55. -INC SMELEME
  56. -INC SMCOORD
  57. -INC SMMODEL
  58. -INC SMINTE
  59. -INC SMLMOTS
  60.  
  61. -INC TMPTVAL
  62.  
  63.       SEGMENT WRK1
  64.        REAL*8 XFORC(LRE), XSTRS(NSTRS), XE(3,NBBB)
  65.       ENDSEGMENT
  66. *
  67.       SEGMENT WRK2
  68.        REAL*8 SHPWRK(6,NBNN), BGENE(NSTRS,LRE)
  69.       ENDSEGMENT
  70. *
  71.       SEGMENT WRK3
  72.        REAL*8 WORK(LW)
  73.       ENDSEGMENT
  74. *
  75.       SEGMENT WRK4
  76.        REAL*8 BPSS(3,3), XEL(3,NBBB), XFOLO(LRE)
  77.       ENDSEGMENT
  78. *
  79.       SEGMENT WRK5
  80.        REAL*8 XGENE(NSTN,LRN)
  81.       ENDSEGMENT
  82. *  segment pour shb8
  83.        SEGMENT WRK7
  84.         REAL*8 PROPEL(24),D(3)
  85.         REAL*8 rel(24,24),work1(30)
  86.        ENDSEGMENT
  87. *
  88.       CHARACTER*4 lesinc(7),lesdua(7)
  89.       DATA lesinc/'UX','UY','UZ','RX','RY','RZ','UR'/
  90.       DATA lesdua/'FX','FY','FZ','MX','MY','MZ','FR'/
  91.       CHARACTER*8 CMATE
  92.       logical dcmat2,dcmate
  93.       dcmat2 = .false.
  94.  
  95.       MELEME=IPMAIL
  96. c*      SEGACT MELEME
  97.       NBNN=NUM(/1)
  98.       NBELEM=NUM(/2)
  99. *
  100. *    introduction des coord du point autour duquel se fait le
  101. *    mouvement de la section en defo plane generalisee
  102. *    et initialisation des forces au noeud support de la defo
  103. *    plane generalisee
  104. *
  105.         ADPG=0.D0
  106.         BDPG=0.D0
  107.         CDPG=0.D0
  108.       IF (IFOUR.EQ.-3)THEN
  109.         IP=IIPDPG
  110. c*        SEGACT MCOORD
  111.         IREF=(IP-1)*(IDIM+1)
  112.         XDPGE=XCOOR(IREF+1)
  113.         YDPGE=XCOOR(IREF+2)
  114.       ELSE
  115.         XDPGE=0.D0
  116.         YDPGE=0.D0
  117.       ENDIF
  118. *
  119.       NHRM=NIFOUR
  120.       MINTE=IPMINT
  121. *
  122. c_______________________________________________________________________
  123. c_______________________________________________________________________
  124. c
  125. c     numero des etiquettes      :
  126. c     etiquettes de 1 a 98 pour traitement specifique a l element
  127. c     dans la zone specifique a chaque element commencant par :
  128. c     5  continue
  129. c     element 5   etiquettes 1005 2005 3005 4005   ...
  130. c     44 continue
  131. c     element 44  etiquettes 1044 2044 3044 4044   ...
  132. c_______________________________________________________________________
  133. c
  134.       IF (MELE.LE.100)
  135.      &GOTO (99,2,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  136.      1      99,99,99,99,99,99,99,99,29,30,99,99,99,99,99,99,99,99,99,99,
  137.      2      99,29,43,99,45,46,99,99,99,30,99,99,99,99,99,99,99,99,99,99,
  138.      3      99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  139.      4      99,99,99,29,99,99,99,99,99,99,99,99,99,99,46,96,99,99,99,99
  140.      5      ),MELE
  141.       IF (MELE.LE.200)
  142.      &GOTO (99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  143.      1      99,99,46,124,125,34,34,34,34,34,34,34,34,34,34,34,34,34,34,
  144.      2      34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,
  145.      3      34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,
  146.      4      34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,
  147.      5      34),MELE-100
  148.       IF (MELE.LE.300)
  149.      &GOTO (34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,
  150.      1      34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,
  151.      2      34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,34,
  152.      3      260,34,34,34,34,265,266,266,266,99,99,271,272),MELE-200
  153. c
  154.  34   CONTINUE
  155.       GOTO 99
  156. C____________________________________________________________________
  157. C
  158. C     ELEMENT SEG2 (pour IMPEDANCE)
  159. C____________________________________________________________________
  160. C
  161.   2   CONTINUE
  162.  
  163.       MPTVAL = IVAFOR
  164.  
  165.       MPTVA1 = IVASTR
  166.       numstr = mptva1.ival(/1)
  167.  
  168.       IF (CMATE.EQ.'IMPELAST'.OR.CMATE.EQ.'IMPVOIGT'.OR.
  169.      &    CMATE.EQ.'IMPREUSS'.OR.CMATE.EQ.'IMPCOMPL') THEN
  170. * detection impedance "sur mesure"
  171.        if (mptval.ival(/1).EQ.numstr*2) dcmat2 = .true.
  172.       ENDIF
  173.  
  174.       DO 310 IB=1,NBELEM
  175. C
  176. C  RANGEMENT DANS MELVAL
  177. C
  178.       DO     IGAU=1,NBNN
  179.       DO     ICOMP=1,NFOR
  180.         if (dcmat2) then
  181.           if(icomp.le.numstr) then
  182.          melva1 = mptva1.ival(icomp)
  183.           else
  184.          melva1 = mptva1.ival(icomp - numstr)
  185.           endif
  186.           if (igau.lt.2) then
  187.          melval = ival(icomp)
  188.           else
  189.          melval = ival(icomp + nfor)
  190.           endif
  191.         else
  192.          MELVA1=MPTVA1.IVAL(ICOMP)
  193.          MELVAL=mptval.IVAL(ICOMP)
  194.         endif
  195.          IBMN=MIN(IB  ,VELCHE(/2))
  196.          IBM1=MIN(IB  ,MELVA1.VELCHE(/2))
  197.          IGM1=MIN(IGAU,MELVA1.VELCHE(/1))
  198.          VELCHE(IGAU,IBMN)= MELVA1.VELCHE(IGM1,IBM1)
  199.       enddo    
  200.       enddo    
  201.  310  CONTINUE
  202.       GO TO 510
  203. c_______________________________________________________________________
  204. c_______________________________________________________________________
  205. c
  206. c  elements poutre et poutre de timoschenko
  207. c_______________________________________________________________________
  208. c
  209.  29   CONTINUE
  210.       if (dcmate) goto 2
  211.       NBBB=NBNN
  212.       SEGINI WRK1,WRK3
  213. c
  214.       NCARR1=NCARR
  215. c
  216.       DO 3029 IB=1,NBELEM
  217. c
  218. c  on cherche les coordonnees des noeuds de l elementib
  219. c
  220.       CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  221. c
  222. c  il faudrait aussi modifier le vecteur local de la poutre
  223. c
  224. c     rangement des caracteristiques dans work
  225. C
  226. C
  227.       MPTVAL=IVACAR
  228.       if(mptval.ne.0) then
  229.       DO 6029 IC=1,NCARR1
  230.         WORK(IC)=0.D0
  231.         MELVAL=IVAL(IC)
  232.         IF (MELVAL.NE.0) THEN
  233.           IBMN=MIN(IB,VELCHE(/2))
  234.           DO 4029 IGAU=1,NBNN
  235.             IGMN=MIN(IGAU,VELCHE(/1))
  236.             WORK(IC)=WORK(IC)+VELCHE(IGMN,IBMN)
  237.  4029     CONTINUE
  238.           WORK(IC)=WORK(IC)/NBNN
  239.         ENDIF
  240.  6029 CONTINUE
  241.       endif
  242. c
  243. c  cas des tuyaux - on calcule les caracteristiques de la poutre
  244. c                   equivalente
  245. c
  246.       IF(MELE.EQ.42) THEN
  247.         CISA=WORK(4)
  248.         VX  =WORK(5)
  249.         VY  =WORK(6)
  250.         VZ  =WORK(7)
  251.         CALL TUYCAR(WORK,CISA,VX,VY,VZ,KERRE,0)
  252.         IF (KERRE.EQ.77) THEN
  253.           CALL ERREUR(77)
  254.           GOTO 510
  255.         ENDIF
  256.       ENDIF
  257. c
  258. c  on cherche les contraintes -
  259. c
  260.       MPTVAL=IVASTR
  261.       IE=9
  262.       DO      IGAU=1,2
  263.       DO      ICOMP=1,NSTRS
  264.         IE=IE+1
  265.         MELVAL=IVAL(ICOMP)
  266.         IBMN=MIN(IB  ,VELCHE(/2))
  267.         IGMN=MIN(IGAU,VELCHE(/1))
  268.         WORK(IE)=VELCHE(IGMN,IBMN)
  269.       enddo    
  270.       enddo    
  271. c
  272. c  on calcule les forces internes
  273. c
  274.       IF(MELE.EQ.84) THEN
  275.         IF(CMATE.NE.'SECTION') THEN
  276.           IF (IFOUR.EQ.-2.OR.IFOUR.EQ.-1.OR.IFOUR.EQ.-3) THEN
  277.            CALL TIMBS2(XFORC,XE,WORK(10),WORK(22),KERRE)
  278.           ELSE
  279.            CALL TIMBSG(XFORC,WORK(7),XE,WORK(10),WORK(22),KERRE)
  280.           ENDIF
  281.         ELSE
  282.           IF (IFOUR.EQ.-2.OR.IFOUR.EQ.-1.OR.IFOUR.EQ.-3) THEN
  283.            CALL TIMBS2(XFORC,XE,WORK(10),WORK(22),KERRE)
  284.           ELSE
  285.             CALL TIMBSG(XFORC,WORK(1),XE,WORK(10),WORK(22),KERRE)
  286.           ENDIF
  287.         ENDIF
  288.       ELSE
  289.           IF (IFOUR.EQ.-2.OR.IFOUR.EQ.-1.OR.IFOUR.EQ.-3) THEN
  290.             CALL POUBS2(XFORC,WORK,XE,WORK(10),WORK(22),KERRE)
  291.           ELSE
  292.             CALL POUBSG(XFORC,WORK,XE,WORK(10),WORK(22),KERRE)
  293.           ENDIF
  294.       ENDIF
  295.       IF(KERRE.NE.0) THEN
  296.         INTERR(1)=ISOUS
  297.         INTERR(2)=IB
  298.         SEGSUP WRK1,WRK3
  299.         CALL ERREUR(128)
  300.         GO TO 510
  301.       ENDIF
  302. c
  303. c  rangement dans melval
  304. c
  305.       IE=0
  306.       MPTVAL=IVAFOR
  307.       DO      IGAU=1,NBNN
  308.       DO      ICOMP=1,NFOR
  309.         IE=IE+1
  310.         MELVAL=IVAL(ICOMP)
  311.         IBMN=MIN(IB  ,VELCHE(/2))
  312.         VELCHE(IGAU,IBMN)=XFORC(IE)
  313.       enddo    
  314.       enddo    
  315. 3029  CONTINUE
  316.       SEGSUP WRK1,WRK3
  317.       GO TO 510
  318. c_______________________________________________________________________
  319. c
  320. c     elements lisp et lism
  321. c_______________________________________________________________________
  322. c
  323.   30  CONTINUE
  324.       NBBB=NBNN
  325.       NBNO=SHPTOT(/2)
  326.       SEGINI WRK1,WRK3,WRK4
  327. c
  328.       DO 3030  IB=1,NBELEM
  329. c
  330. c     on cherche les coordonnees des noeuds de l element ib
  331. c
  332.       CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  333. c
  334. c
  335. c     on cherche les contraintes
  336. c
  337.       IE=0
  338.       MPTVAL=IVASTR
  339.       DO      IGAU=1,NBPGAU
  340.       DO      ICOMP=1,NSTRS
  341.         IE=IE+1
  342.         MELVAL=IVAL(ICOMP)
  343.         IGMN=MIN(IGAU,VELCHE(/1))
  344.         IBMN=MIN(IB  ,VELCHE(/2))
  345.         WORK(IE)=VELCHE(IGMN,IBMN)
  346.       enddo     
  347.       enddo     
  348. c
  349. c     on cherche les caracteristiques
  350. c
  351.       MPTVAL=IVACAR
  352.       DO      IGAU=1,NBPGAU
  353.       DO      ICOMP=1,NCARR
  354.         IE=IE+1
  355.         MELVAL=IVAL(ICOMP)
  356.         IF (MELVAL.NE.0) THEN
  357.           IGMN=MIN(IGAU,VELCHE(/1))
  358.           IBMN=MIN(IB  ,VELCHE(/2))
  359.           WORK(IE)=VELCHE(IGMN,IBMN)
  360.         ELSE
  361.           WORK(IE)=0.D0
  362.         ENDIF
  363.       enddo     
  364.       enddo     
  365. c
  366. c     on calcule b*sigma
  367. c
  368.       ICNT=NBPGAU*NSTRS+1
  369.       CALL LISPBS(WORK(1),WORK(ICNT),POIGAU,SHPTOT,
  370.      1    NBPGAU,NBNO,XE,XFOLO,BPSS,XFORC)
  371. c
  372. c    rangement dans melval
  373. c
  374.       IE=0
  375.       MPTVAL=IVAFOR
  376.       DO      IGAU=1,NBNN
  377.       DO      ICOMP=1,6
  378.         IE=IE+1
  379.         MELVAL=IVAL(ICOMP)
  380.         IBMN=MIN(IB  ,VELCHE(/2))
  381.         VELCHE(IGAU,IBMN)=XFORC(IE)
  382.       enddo     
  383.       enddo     
  384.  3030 CONTINUE
  385.       SEGSUP WRK1,WRK3,WRK4
  386.       GOTO 510
  387. c_______________________________________________________________________
  388. c
  389. c    element tuyau fissure
  390. c_______________________________________________________________________
  391. c
  392. 43    CONTINUE
  393.       NBBB=NBNN
  394.       SEGINI WRK1,WRK3
  395. c
  396.       DO 3043 IB=1,NBELEM
  397. c
  398. c   on cherche les contraintes
  399. c
  400.       IE=0
  401.       MPTVAL=IVASTR
  402.       DO      IGAU=1,NBPTEL
  403.       DO      ICOMP=1,NSTRS
  404.         IE=IE+1
  405.         MELVAL=IVAL(ICOMP)
  406.         IGMN=MIN(IGAU,VELCHE(/1))
  407.         IBMN=MIN(IB  ,VELCHE(/2))
  408.         WORK(IE)=VELCHE(IGMN,IBMN)
  409.       enddo    
  410.       enddo    
  411. c
  412. c     on cherche les caracteristiques
  413. c
  414.       MPTVAL=IVACAR
  415.       DO 5043 ICOMP=1,NCARR
  416.         MELVAL=IVAL(ICOMP)
  417.         IF (MELVAL.NE.0) THEN
  418.           IBMN=MIN(IB  ,VELCHE(/2))
  419.           WORK(ICOMP+8)=VELCHE(1,IBMN)
  420.         ELSE
  421.           WORK(ICOMP+8)=0.D0
  422.         ENDIF
  423. 5043  CONTINUE
  424. c
  425. c     on calcule les forces internes
  426. c
  427.       CALL TUFIBS(XFORC,WORK,WORK(9),WORK(18),KERRE)
  428. c
  429. c    rangement dans melval
  430. c
  431.       IE=0
  432.       MPTVAL=IVAFOR
  433.       DO      IGAU=1,NBNN
  434.       DO      ICOMP=1,NFOR
  435.         IE=IE+1
  436.         MELVAL=IVAL(ICOMP)
  437.         IBMN=MIN(IB  ,VELCHE(/2))
  438.         VELCHE(IGAU,IBMN)=XFORC(IE)
  439.       enddo     
  440.       enddo     
  441. 3043  CONTINUE
  442.       SEGSUP WRK1,WRK3
  443.       GO TO 510
  444. c_______________________________________________________________________
  445. c
  446. c  element point (poi1) defos planes generalisees/materiau IMPEDANCE
  447. c_______________________________________________________________________
  448. c
  449.  45   CONTINUE
  450.       NBBB=NBNN
  451. *
  452.       IF ((CMATE.EQ.'IMPELAST').OR.(CMATE.EQ.'IMPVOIGT').OR.
  453.      &    (CMATE.EQ.'IMPREUSS').OR.(cmate.eq.'IMPCOMPL').or.
  454.      &    (cmate.eq.'MODAL   ').or.(cmate.eq.'STATIQUE')) THEN
  455.        mptva1 = ivastr
  456.        mptval = ivafor
  457.        segact mptva1,mptval
  458.        do iv = 1,ival(/1)
  459.         melva1 = mptva1.ival(iv)
  460.         melval = ival(iv)
  461.         segact melva1,melval*mod
  462.         DO  IB=1,NBELEM
  463.  
  464.          DO 4145 IGAU=1,NBNN
  465.             IGMN=MIN(IGAU,MELVA1.VELCHE(/1))
  466.             IBMN=MIN(IB  ,MELVA1.VELCHE(/2))
  467.             valstr = MELVA1.VELCHE(IGMN,IBMN)
  468.             VELCHE(IGMN,IBMN) = valstr
  469.  4145    CONTINUE
  470.         ENDDO
  471.  
  472.        enddo
  473.  
  474.        GOTO 510
  475.       ENDIF
  476. *
  477.       IF(MELE.EQ.45.AND.IFOUR.NE.-3) THEN
  478.         GO TO 99
  479.       ENDIF
  480. *
  481.       SEGINI WRK1,WRK3
  482. c
  483.       DO 3045 IB=1,NBELEM
  484. c
  485. c  on cherche les coordonnees des noeuds de l element ib
  486. c
  487.       CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  488. c
  489. c  mise a zero des forces internes
  490. c
  491.       CALL ZERO(XFORC,1,LRE)
  492. c
  493. c  on cherche l'effort
  494. c
  495.       MPTVAL=IVASTR
  496.       MELVAL=IVAL(1)
  497.       IBMN=MIN(IB  ,VELCHE(/2))
  498.       EFFORT=VELCHE(1,IBMN)
  499. c
  500. c  on calcule les forces internes
  501. c
  502.        CALL PO1BSG(XE,LRE,XDPGE,YDPGE,EFFORT,XFORC)
  503.        ADPG=ADPG+XFORC(3)
  504.        BDPG=BDPG+XFORC(4)
  505.        CDPG=CDPG+XFORC(5)
  506. c
  507. c  rangement dans melval
  508. c
  509.       IE=0
  510.       MPTVAL=IVAFOR
  511.       DO      IGAU=1,NBNN
  512.       DO      ICOMP=1,NFOR-3
  513.         IE=IE+1
  514.         MELVAL=IVAL(ICOMP)
  515.         IBMN=MIN(IB  ,VELCHE(/2))
  516.         VELCHE(IGAU,IBMN)=XFORC(IE)
  517.       enddo    
  518.       enddo    
  519. 3045  CONTINUE
  520. c
  521.       SEGSUP WRK1,WRK3
  522.       GO TO 510
  523. c_______________________________________________________________________
  524. c
  525. c  elements barre  et cerce
  526. c_______________________________________________________________________
  527. c
  528.  46   CONTINUE
  529.       NBBB=NBNN
  530. *
  531.       IF(MELE.EQ.95.AND.IFOUR.NE.0.AND.IFOUR.NE.1) THEN
  532.         GO TO 99
  533.       ENDIF
  534. *
  535.       SEGINI WRK1,WRK3
  536. c
  537.       DO 3046 IB=1,NBELEM
  538. c
  539. c  on cherche les coordonnees des noeuds de l element ib
  540. c
  541.       CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  542. C
  543.       IF(MELE.EQ.123) THEN
  544. c
  545. c  on cherche l'effort
  546. c
  547.         IE=0
  548.         MPTVAL=IVASTR
  549.         DO      IGAU=1,NBPTEL
  550.         DO      ICOMP=1,NSTRS
  551.            IE=IE+1
  552.            MELVAL=IVAL(ICOMP)
  553.            IGMN=MIN(IGAU,VELCHE(/1))
  554.            IBMN=MIN(IB  ,VELCHE(/2))
  555.            WORK(IE)=VELCHE(IGMN,IBMN)
  556.         enddo    
  557.         enddo    
  558. c
  559. c  on calcule les forces internes
  560. c
  561.         CALL BARBS3(XFORC,XE,WORK,KERRE,QSIGAU,POIGAU,NBPGAU,IB)
  562. c
  563.       ELSE
  564. c
  565. c    on cherche l'effort
  566. c
  567.         MPTVAL=IVASTR
  568.         MELVAL=IVAL(1)
  569.         NPPTEL=VELCHE(/1)
  570.         IBMN=MIN(IB  ,VELCHE(/2))
  571.         IF(NPPTEL.EQ.1) THEN
  572.           EFFORT=VELCHE(1,IBMN)
  573.         ELSE IF(NPPTEL.EQ.2) THEN
  574.           EFFOR1=VELCHE(1,IBMN)
  575.           EFFOR2=VELCHE(2,IBMN)
  576.           EFFORT=0.5D0*(EFFOR1+EFFOR2)
  577.         ENDIF
  578. c
  579. c  on calcule les forces internes
  580. c
  581.         IF(MELE.EQ.46) CALL BARBSG(XFORC,XE,EFFORT,KERRE)
  582.         IF(MELE.EQ.95) CALL CERBSG(XFORC,XE,EFFORT,KERRE)
  583.       ENDIF
  584. C
  585.       IF(KERRE.NE.0) THEN
  586.         INTERR(1)=ISOUS
  587.         INTERR(2)=IB
  588.         SEGSUP WRK1,WRK3
  589.         IF(MELE.EQ.46) CALL ERREUR(128)
  590.         IF(MELE.EQ.95) CALL ERREUR(128)
  591.         GO TO 510
  592.       ENDIF
  593. c
  594. c  rangement dans melval
  595. c
  596.       NFOD=NFOR
  597.       IF (IFOUR.EQ.-3.AND.MELE.EQ.46) NFOD=NFOR-3
  598.       IE=0
  599.       MPTVAL=IVAFOR
  600.       DO      IGAU=1,NBNN
  601.       DO      ICOMP=1,NFOD
  602.         IE=IE+1
  603.         MELVAL=IVAL(ICOMP)
  604.         IBMN=MIN(IB  ,VELCHE(/2))
  605.         VELCHE(IGAU,IBMN)=XFORC(IE)
  606.       enddo    
  607.       enddo    
  608. 3046  CONTINUE
  609.       SEGSUP WRK1,WRK3
  610.       GO TO 510
  611. C_______________________________________________________________________
  612. C
  613. C     ELEMENT BARRE 3D EXCENTRE (BAEX)
  614. C_______________________________________________________________________
  615. C
  616.  124  CONTINUE
  617.       NBBB=NBNN
  618.       NSTN=NBNN
  619.       LRN =LRE
  620.       SEGINI WRK1,WRK3,WRK5
  621. C
  622. C     BOUCLE DE CALCUL POUR LES DIFFERENTS ELEMENTS
  623. C
  624.       KERRE=0
  625.       DO 3124 IB=1,NBELEM
  626. C
  627. C  ON RECUPERE LA SECTION DE L'ELEMENT, SES EXCENTREMENTS ET SON
  628. C  ORIENTATION. LES CARACTERISTIQUES SONT RANGEES DANS WORK
  629. C  SELON L'ORDRE SUIVANT: SECT EXCZ EXCY VX VY VZ
  630. C
  631.       MPTVAL=IVACAR
  632.       DO IC=1,NCARR
  633.         IF(IVAL(IC).NE.0) THEN
  634.           MELVAL=IVAL(IC)
  635.           IBMN=MIN(IB,VELCHE(/2))
  636.           WORK(IC)=VELCHE(1,IBMN)
  637.         ELSE
  638.           WORK(IC)=0.D0
  639.         ENDIF
  640.       END DO
  641.       SECT=WORK(1)
  642. C
  643. C   XGENE  STOCKE LA MATRICE DE PASSAGE DE L'ELEMENT EXCENTRE
  644. C
  645.       CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  646.       CALL MAPAEX(XE,NBNN,WORK,AL,XGENE,LRE,KERRE)
  647.       IF(KERRE.NE.0) THEN
  648.         INTERR(1)=ISOUS
  649.         INTERR(2)=IB
  650.         IF(KERRE.EQ.1) CALL ERREUR(128)
  651.       ENDIF
  652. C
  653. C  MISE A ZERO DES FORCES INTERNES
  654. C
  655.       CALL ZERO(XFORC,1,LRE)
  656. C
  657. C  ON CHERCHE L'EFFORT
  658. C
  659.       MPTVAL=IVASTR
  660.       MELVAL=IVAL(1)
  661.       IBMN=MIN(IB  ,VELCHE(/2))
  662.       NPPTEL=VELCHE(/1)
  663.       IF(NPPTEL.EQ.1) THEN
  664.         EFFORT=VELCHE(1,IBMN)
  665.       ELSE IF(NPPTEL.EQ.2) THEN
  666.         EFFOR1=VELCHE(1,IBMN)
  667.         EFFOR2=VELCHE(2,IBMN)
  668.         EFFORT=0.5D0*(EFFOR1+EFFOR2)
  669.       ENDIF
  670. CC      EFFORT=SECT*EFFORT
  671. C
  672. C  ON CALCULE LES FORCES INTERNES
  673. C
  674.       CALL BARINT(XFORC,XGENE,EFFORT,LRE)
  675. C
  676. C  RANGEMENT DANS MELVAL
  677. C
  678.       IE=0
  679.       MPTVAL=IVAFOR
  680.       DO      IGAU=1,NBNN
  681.       DO      ICOMP=1,NFOR
  682.         IE=IE+1
  683.         MELVAL=IVAL(ICOMP)
  684.         IBMN=MIN(IB  ,VELCHE(/2))
  685.         VELCHE(IGAU,IBMN)=XFORC(IE)
  686.       enddo    
  687.       enddo    
  688. 3124  CONTINUE
  689.       SEGSUP WRK1,WRK3,WRK5
  690.       GO TO 510
  691. c cccccccccccccccccccccccccccccccccccccccccccccccccccccc c
  692. c    element coaxial COS2 (3D pour liaison acier-beton)  c
  693. c cccccccccccccccccccccccccccccccccccccccccccccccccccccc c
  694. 271   continue
  695.       NBBB=NBNN
  696.       NSTN=NBNN
  697.       LRN =LRE
  698.       SEGINI WRK1,WRK3,WRK4,WRK5
  699.       do 2713 ib=1,nbelem
  700. C
  701.       CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  702.       vx1= xe(1,2) - xe(1,1)
  703.       vy1= xe(2,2) - xe(2,1)
  704.       vz1=0.d0
  705.       if(idim.eq.3) vz1= xe(3,2)-xe(3,1)
  706.       xl= sqrt( vx1*vx1 + vy1*vy1 + vz1*vz1)
  707.       vx1= vx1 / xl
  708.       vy1= vy1/ xl
  709.       if(idim.eq.3) THEN
  710.         vz1=vz1 / xl
  711.       ENDIF
  712.  
  713.       IE=0
  714.       MPTVAL=IVASTR
  715.       DO      IGAU=1,2
  716.       DO      ICOMP=1,NSTRS
  717.         IE=IE+1
  718.         MELVAL=IVAL(ICOMP)
  719.         IGMN=MIN(IGAU,VELCHE(/1))
  720.         IBMN=MIN(IB  ,VELCHE(/2))
  721.         WORK(IE)=VELCHE(IGMN,IBMN)
  722.       enddo    
  723.       enddo    
  724. C
  725.       MPTVAL=IVAcar
  726.       DO 2712 ICOMP=1,NCARR
  727.           MELVAL=IVAL(ICOMP)
  728.           IGMN = VELCHE(/1)
  729.           IBMN=MIN(IB,VELCHE(/2))
  730.           SECA =VELCHE(IGMN,IBMN)
  731. 2712  CONTINUE
  732.       diam = sqrt(4.d0*SECA/xpi)
  733. C
  734. C  MISE A ZERO DES FORCES INTERNES
  735. C
  736.       CALL ZERO(XFORC,1,LRE)
  737. C
  738. C
  739. C FORCES DANS LA DIRECTION TANGENTIELLE
  740. C
  741.       ag = 1.d0-0.5773502691896257645d0
  742.       agg = ag/(2.d0*(1-ag))
  743.       t11 = work(1) + (work(1) -work(1 + idim))*agg
  744.       t21 = work(1 + idim) + (work(1 + idim) -work(1))*agg
  745.       FO11 = xpi*diam*xl*(t11/2.d0 + (t21 - t11)/8.d0)
  746.       FO21 = xpi*diam*xl*(t21/2.d0 + (t11 - t21)/8.d0)
  747. C
  748. C FORCES DANS LA DIRECTION NORMALE
  749. C
  750.       t12 = WORK(2) + (work(2) -work(2 + idim))*agg
  751.       t22 = WORK(2 + IDIM)+ (work(2 + idim) -work(2))*agg
  752.       FO12 = -1.d0*diam*xl*(t12/2.d0 + (t22 - t12)/8.d0)
  753.       FO22 = -1.d0*diam*xl*(t22/2.d0 + (t12 - t22)/8.d0)
  754. c      write(6,*) 'xstrs 1',ib,t11,t12,t12
  755. c      write(6,*) 'xstrs 2',ib,t21,t22,t22
  756.       IF  (IDIM.EQ.2) THEN
  757.            XFORC(1) = (FO11*VX1 + FO12*VY1) + XFORC(1)
  758.            XFORC(3) = (FO21*VX1 + FO22*VY1) + XFORC(3)
  759.            XFORC(2) = (FO11*VY1 + FO12*VX1) + XFORC(2)
  760.            XFORC(4) = (FO21*VY1 + FO22*VX1) + XFORC(4)
  761.            XFORC(5) = (-1.d0*(FO11*VX1 + FO12*VY1)) + XFORC(5)
  762.            XFORC(7) = (-1.d0*(FO21*VX1 + FO22*VY1)) + XFORC(7)
  763.            XFORC(6) = (-1.d0*(FO11*VY1 + FO12*VX1)) + XFORC(6)
  764.            XFORC(8) = (-1.d0*(FO21*VY1 + FO22*VX1)) + XFORC(8)
  765.       ELSE IF (IDIM.EQ.3) THEN
  766.           IF (vy1.EQ.0.0D0.AND.vz1.EQ.0.0D0) THEN
  767.             vx22 = 0.0D0
  768.             vy22 = 1.0D0
  769.             vz22 = 0.0D0
  770.           ELSE IF ((vx1.EQ.0.0D0.AND.vy1.EQ.0.0D0).OR.
  771.      .            (vx1.EQ.0.0D0.AND.vz1.EQ.0.0D0)) THEN
  772.             vx22 = 1.0D0
  773.             vy22 = 0.0D0
  774.             vz22 = 0.0D0
  775.           ELSE IF (vy1.NE.0.0D0.AND.vz1.NE.0.0D0) THEN
  776.             Vx22 = 0.0D0
  777.             Vy22 = -vz1
  778.             Vz22 = vy1
  779.             LLL = 1
  780.           ELSE IF (vx1.NE.0.0D0.AND.vz1.NE.0.0D0) THEN
  781.             Vx22 = -vz1
  782.             Vy22 = 0.0D0
  783.             Vz22 = vx1
  784.             LLL = 1
  785.           ELSE IF (vy1.NE.0.0D0.AND.vx1.NE.0.0D0) THEN
  786.             Vx22 = -vy1
  787.             Vy22 = vx1
  788.             Vz22 = 0.0D0
  789.             LLL = 1
  790.           END IF
  791.           xl22 = sqrt((vx22*vx22) + (vy22*vy22)+ (vz22*vz22))
  792.           vx22 = vx22/xl22
  793.           vy22 = vy22/xl22
  794.           vz22 = vz22/xl22
  795.           vx3 = (vy1*Vz22) - (vz1*Vy22)
  796.           vy3 = (vz1*Vx22) - (vx1*vz22)
  797.           vz3 = (vx1*vy22) - (vy1*Vx22)
  798.           xl3 = sqrt((vx3*vx3) + (vy3*vy3)+ (vz3*vz3))
  799.           vx3 = vx3/xl3
  800.           vy3 = vy3/xl3
  801.           vz3 = vz3/xl3
  802.           vx2 = (vy3*vz1) - (vz3*vy1)
  803.           vy2 = (vz3*vx1) - (vx3*vz1)
  804.           vz2 = (vx3*vy1) - (vy3*vx1)
  805.           xl2 = sqrt((vx2*vx2) + (vy2*vy2)+ (vz2*vz2))
  806.           vx2 = vx2/xl2
  807.           vy2 = vy2/xl2
  808.           vz2 = vz2/xl2
  809.           F1 = FO11*vx1 + FO12*vx2 + FO12*vx3
  810.           F2 = FO11*vy1 + FO12*vy2 + FO12*vy3
  811.           F3 = FO11*vz1 + FO12*vz2 + FO12*vz3
  812.           F4 = FO21*vx1 + FO22*vx2 + FO22*vx3
  813.           F5 = FO21*vy1 + FO22*vy2 + FO22*vy3
  814.           F6 = FO21*vz1 + FO22*vz2 + FO22*vz3
  815.           XFORC(1) = F1 + XFORC(1)
  816.           XFORC(2) = F2 + XFORC(2)
  817.           XFORC(3) = F3 + XFORC(3)
  818.           XFORC(4) = F4 + XFORC(4)
  819.           XFORC(5) = F5 + XFORC(5)
  820.           XFORC(6) = F6 + XFORC(6)
  821.           XFORC(7) =  -1.d0*F4 + XFORC(7)
  822.           XFORC(8) =  -1.d0*F5 + XFORC(8)
  823.           XFORC(9) =  -1.d0*F6 + XFORC(9)
  824.           XFORC(10) = -1.d0*F1 + XFORC(10)
  825.           XFORC(11) = -1.d0*F2 + XFORC(11)
  826.           XFORC(12) = -1.d0*F3 + XFORC(12)
  827.        ENDIF
  828. c       write(6,*) 'xforc cos2', (xforc(io),io = 1,lre)
  829. C
  830. C   RANGEMENT DANS MELVAL
  831. C
  832.       IE=0
  833.       MPTVAL=IVAFOR
  834. C
  835. C         NODE=4= NOMBRE DE NOEUDS
  836. C         ICOMP=NSTRS= NOMBRE DE COMPOSANTES DE CONTRAINTES PAR NOEUD
  837. C
  838.       DO      NODE=1,4
  839.       DO      ICOMP=1,NSTRS
  840.         IE=IE+1
  841.         MELVAL=IVAL(ICOMP)
  842.         IBMN=MIN(IB  ,VELCHE(/2))
  843.         VELCHE(NODE,IBMN)=XFORC(IE)
  844.       enddo     
  845.       enddo     
  846.  
  847.  2713 continue
  848.  
  849.         SEGSUP WRK1,WRK3,WRK4,WRK5
  850.          GO TO 510
  851.  
  852. C_______________________________________________________________________
  853. C
  854. C  ELEMENT COAXIAL (COA2)
  855. C_______________________________________________________________________
  856. C
  857. 272   continue
  858.       NBNO=NBNN
  859.       NBBB=NBNN
  860.       SEGINI WRK1,WRK2,WRK3,WRK4
  861.       LW = 24
  862. C
  863.       DO 2721 IB=1,NBELEM
  864. C
  865. C  ON CHERCHE LES COORDONNEES DES NOEUDS DE L ELEMENT IB
  866. C
  867.       CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  868. C
  869. C  MISE A ZERO DES FORCES INTERNES
  870. C
  871.       CALL ZERO(XFORC,1,LRE)
  872. C
  873. C     REPERE LOCAL
  874. C
  875.       CALL CO2LOC(XE,SHPTOT,NBNN,XEL,BPSS,NOQUAL,IDIM)
  876. C
  877. C  BOUCLE SUR LES POINTS DE GAUSS
  878. C
  879.       DO 2722 IGAU=1,NBPGAU
  880. C
  881. C     CALCUL DE LA MATRICE B ET DU JACOBIEN EN IGAU
  882. C
  883.       CALL BCO2(IGAU,MFR,IFOUR,NIFOUR,XEL,BPSS,SHPTOT,SHPWRK,
  884.      .               BGENE,DJAC,IRRT,IDIM,NBNN,NSTRS,LRE)
  885.       IF(IRRT.NE.0) THEN
  886.              INTERR(1)=IB
  887.              CALL ERREUR(764)
  888.              GOTO 9985
  889.       END IF
  890. C
  891. C  ON CHERCHE LES CONTRAINTES -
  892. C
  893.       MPTVAL=IVASTR
  894.       DO 2723 ICOMP=1,NSTRS
  895.         MELVAL=IVAL(ICOMP)
  896.         IGMN=MIN(IGAU,VELCHE(/1))
  897.         IBMN=MIN(IB  ,VELCHE(/2))
  898.         XSTRS(ICOMP)=VELCHE(IGMN,IBMN)
  899.  2723 CONTINUE
  900. C
  901. C
  902. C     ON CHERCHE LES COORDONNEES DES NOEUDS DE L ELEMENT IB
  903. C
  904.       CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  905.       vx1= xe(1,2) - xe(1,1)
  906.       vy1= xe(2,2) - xe(2,1)
  907.       vz1=0.d0
  908.       if(idim.eq.3) vz1= xe(3,2)-xe(3,1)
  909.       xl= sqrt( vx1*vx1 + vy1*vy1 + vz1*vz1)
  910.       vx1= vx1 / xl
  911.       vy1= vy1/ xl
  912.       if(idim.eq.3) THEN
  913.         vz1=vz1 / xl
  914.       ENDIF
  915. C
  916. C recuperation de la section et calcul du diamètre
  917. C
  918.       MPTVAL=IVACAR
  919.       DO 2729 ICOMP=1,NCARR
  920.           MELVAL=IVAL(ICOMP)
  921.           IGMN = VELCHE(/1)
  922.           IBMN=MIN(IB,VELCHE(/2))
  923.           SECA =VELCHE(IGMN,IBMN)
  924. 2729  CONTINUE
  925.       diam = sqrt(4.d0*SECA/xpi)
  926. C
  927. C  ON CALCULE B*EFFORTS
  928. C
  929.       DJAC=DJAC*POIGAU(IGAU)
  930.       DO i=1,LRE
  931.             cc = 0.D0
  932.             DO j=1,NSTRS
  933.                cc = cc + (BGENE(j,i) * XSTRS(j))
  934.             ENDDO
  935.             WORK(i) = xl * cc * DJAC * diam
  936.       ENDDO
  937.  
  938.       IF (IDIM.EQ.2) THEN
  939.           WORK(1) = WORK(1) * xpi
  940.           WORK(3) = WORK(3) * xpi
  941.           WORK(5) = WORK(5) * xpi
  942.           WORK(7) = WORK(7) * xpi
  943.           WORK(1+4*idim) = WORK(1)*vx1 - WORK(2)*vy1
  944.           WORK(2+4*idim) = WORK(1)*vy1 + WORK(2)*vx1
  945.           WORK(3+4*idim) = WORK(3)*vx1 - WORK(4)*vy1
  946.           WORK(4+4*idim) = WORK(3)*vy1 + WORK(4)*vx1
  947.           WORK(5+4*idim) = WORK(5)*vx1 - WORK(6)*vy1
  948.           WORK(6+4*idim) = WORK(5)*vy1 + WORK(6)*vx1
  949.           WORK(7+4*idim) = WORK(7)*vx1 - WORK(8)*vy1
  950.           WORK(8+4*idim) = WORK(7)*vy1 + WORK(8)*vx1
  951.       ELSE IF (IDIM.EQ.3) THEN
  952.           WORK(1) = WORK(1) * xpi
  953.           WORK(4) = WORK(4) * xpi
  954.           WORK(7) = WORK(7) * xpi
  955.           WORK(10) = WORK(10) * xpi
  956.           IF (vy1.EQ.0.0D0.AND.vz1.EQ.0.0D0) THEN
  957.             vx22 = 0.0D0
  958.             vy22 = 1.0D0
  959.             vz22 = 0.0D0
  960.           ELSE IF ((vx1.EQ.0.0D0.AND.vy1.EQ.0.0D0).OR.
  961.      .            (vx1.EQ.0.0D0.AND.vz1.EQ.0.0D0)) THEN
  962.             vx22 = 1.0D0
  963.             vy22 = 0.0D0
  964.             vz22 = 0.0D0
  965.           ELSE IF (vy1.NE.0.0D0.AND.vz1.NE.0.0D0) THEN
  966.             Vx22 = 0.0D0
  967.             Vy22 = -vz1
  968.             Vz22 = vy1
  969.             LLL = 1
  970.           ELSE IF (vx1.NE.0.0D0.AND.vz1.NE.0.0D0) THEN
  971.             Vx22 = -vz1
  972.             Vy22 = 0.0D0
  973.             Vz22 = vx1
  974.             LLL = 1
  975.           ELSE IF (vy1.NE.0.0D0.AND.vx1.NE.0.0D0) THEN
  976.             Vx22 = -vy1
  977.             Vy22 = vx1
  978.             Vz22 = 0.0D0
  979.             LLL = 1
  980.           END IF
  981.           xl22 = sqrt((vx22*vx22) + (vy22*vy22)+ (vz22*vz22))
  982.           vx22 = vx22/xl22
  983.           vy22 = vy22/xl22
  984.           vz22 = vz22/xl22
  985.           vx3 = (vy1*Vz22) - (vz1*Vy22)
  986.           vy3 = (vz1*Vx22) - (vx1*vz22)
  987.           vz3 = (vx1*vy22) - (vy1*Vx22)
  988.           xl3 = sqrt((vx3*vx3) + (vy3*vy3)+ (vz3*vz3))
  989.           vx3 = vx3/xl3
  990.           vy3 = vy3/xl3
  991.           vz3 = vz3/xl3
  992.           vx2 = (vy3*vz1) - (vz3*vy1)
  993.           vy2 = (vz3*vx1) - (vx3*vz1)
  994.           vz2 = (vx3*vy1) - (vy3*vx1)
  995.           xl2 = sqrt((vx2*vx2) + (vy2*vy2)+ (vz2*vz2))
  996.           vx2 = vx2/xl2
  997.           vy2 = vy2/xl2
  998.           vz2 = vz2/xl2
  999.           WORK(1+4*idim) = WORK(1)*vx1 + WORK(2)*vx2 + WORK(3)*vx3
  1000.           WORK(2+4*idim) = WORK(1)*vy1 + WORK(2)*vy2 + WORK(3)*vy3
  1001.           WORK(3+4*idim) = WORK(1)*vz1 + WORK(2)*vz2 + WORK(3)*vz3
  1002.           WORK(4+4*idim) = WORK(4)*vx1 + WORK(5)*vx2 + WORK(6)*vx3
  1003.           WORK(5+4*idim) = WORK(4)*vy1 + WORK(5)*vy2 + WORK(6)*vy3
  1004.           WORK(6+4*idim) = WORK(4)*vz1 + WORK(5)*vz2 + WORK(6)*vz3
  1005.           WORK(7+4*idim) = WORK(7)*vx1 + WORK(8)*vx2 + WORK(9)*vx3
  1006.           WORK(8+4*idim) = WORK(7)*vy1 + WORK(8)*vy2 + WORK(9)*vy3
  1007.           WORK(9+4*idim) = WORK(7)*vz1 + WORK(8)*vz2 + WORK(9)*vz3
  1008.           WORK(10+4*idim) = WORK(10)*vx1 + WORK(11)*vx2 + WORK(12)*vx3
  1009.           WORK(11+4*idim) = WORK(10)*vy1 + WORK(11)*vy2 + WORK(12)*vy3
  1010.           WORK(12+4*idim) = WORK(10)*vz1 + WORK(11)*vz2 + WORK(12)*vz3
  1011.           DO i=1,LRE
  1012.             XFORC(i) = WORK(i+4*idim) + XFORC(i)
  1013.           ENDDO
  1014.       END IF
  1015.  2722 CONTINUE
  1016. c      write(6,*) 'xforc coa2', (xforc(io),io = 1,lre)
  1017. C
  1018. C   RANGEMENT DANS MELVAL
  1019. C
  1020.       IE=0
  1021.       MPTVAL=IVAFOR
  1022. C
  1023.       DO      NODE=1,NBNN
  1024. **  AM  30/04/19
  1025. **    DO 2724 ICOMP=1,NSTRS
  1026.       DO      ICOMP=1,IDIM
  1027.         IE=IE+1
  1028.         MELVAL=IVAL(ICOMP)
  1029.         IBMN=MIN(IB  ,VELCHE(/2))
  1030.         VELCHE(NODE,IBMN)=XFORC(IE)
  1031.       enddo     
  1032.       enddo     
  1033.  2721 CONTINUE
  1034.  9985 CONTINUE
  1035.         SEGSUP WRK1,WRK2,WRK3,WRK4
  1036.          GO TO 510
  1037.  
  1038. *-----------------------------------------------
  1039. *                     element shb8
  1040. *------------------------------------------------
  1041.  260  CONTINUE
  1042.       NBBB=NBNN
  1043.       SEGINI WRK1,WRK7
  1044.       DO 3260  IB=1,NBELEM
  1045. c
  1046. c     on cherche les coordonnees des noeuds de l element ib
  1047. c
  1048.       CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1049. c
  1050. c     on cherche les contraintes
  1051. c
  1052.       IE=0
  1053.       MPTVAL=IVASTR
  1054.       DO      IGAU=1,NBPGAU
  1055.       DO      ICOMP=1,NSTRS
  1056.         IE=IE+1
  1057.         MELVAL=IVAL(ICOMP)
  1058.         IGMN=MIN(IGAU,VELCHE(/1))
  1059.         IBMN=MIN(IB  ,VELCHE(/2))
  1060.         WORK1(IE)=VELCHE(IGMN,IBMN)
  1061.       enddo    
  1062.       enddo    
  1063. c
  1064. c     on cherche les materiaux
  1065. c
  1066.       yyo=0.d0
  1067.       ynu=0.d0
  1068.       MPTVAL=IVAMAT
  1069.       segact mptval
  1070.       DO      IGAU=1,NBPGAU
  1071.       DO      ICOMP=1,2
  1072.         MELVAL=IVAL(ICOMP)
  1073.         IF (MELVAL.NE.0) THEN
  1074.           IGMN=MIN(IGAU,VELCHE(/1))
  1075.           IBMN=MIN(IB  ,VELCHE(/2))
  1076.           if(icomp.eq.1) yyo=yyo+VELCHE(IGMN,IBMN)/nbpgau
  1077.           if(icomp.eq.2)ynu=ynu+VELCHE(IGMN,IBMN)/nbpgau
  1078.         ENDIF
  1079.       enddo    
  1080.       enddo    
  1081. c
  1082. c     on calcule b*sigma
  1083. c
  1084.       d(1)=yyo
  1085.       d(2)=ynu
  1086.       d(3)=0
  1087.       CALL shb8(8,XE,D,propel,work1,rel,xforc)
  1088. c
  1089. c    rangement dans melval
  1090. c
  1091.       IE=0
  1092.       MPTVAL=IVAFOR
  1093.       DO      IGAU=1,Nbnn
  1094.       DO      ICOMP=1,3
  1095.         IE=IE+1
  1096.         MELVAL=IVAL(ICOMP)
  1097.         IBMN=MIN(IB  ,VELCHE(/2))
  1098.         VELCHE(IGAU,IBMN)=XFORC(IE)
  1099.       enddo    
  1100.       enddo    
  1101.  3260 CONTINUE
  1102.       SEGSUP WRK1,WRk7
  1103.       GO TO 510
  1104. C_______________________________________________________________________
  1105. C
  1106. C     LIA2 : element de liaison a 2 noeuds (6 ddl par
  1107. C noeuds)
  1108. C_______________________________________________________________________
  1109. C
  1110. 125   CONTINUE
  1111.       NBBB=NBNN
  1112.       NSTN=3
  1113.       LRN =3
  1114.       SEGINI WRK1,WRK3,WRK5
  1115. C
  1116.       KERRE=0
  1117.       DO 3109 IB=1,NBELEM
  1118. C
  1119. C  MISE A ZERO DES FORCES INTERNES
  1120. C
  1121.       CALL ZERO(XFORC,1,LRE)
  1122. C
  1123.       MPTVAL=IVACAR
  1124.       DO IC=1,NCARR
  1125.         IF(IVAL(IC).NE.0) THEN
  1126.           MELVAL=IVAL(IC)
  1127.           IBMN=MIN(IB,VELCHE(/2))
  1128.           WORK(IC)=VELCHE(1,IBMN)
  1129.         ELSE
  1130.           WORK(IC)=0.D0
  1131.         ENDIF
  1132.       END DO
  1133. C
  1134.       CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1135.       CALL MAPALI(XE,NBNN,WORK,XGENE,KERRE)
  1136.       IF(KERRE.NE.0) THEN
  1137.         INTERR(1)=ISOUS
  1138.         INTERR(2)=IB
  1139.         IF(KERRE.EQ.1) CALL ERREUR(128)
  1140.       ENDIF
  1141. C
  1142. C  ON CHERCHE LES CONTRAINTES -
  1143. C
  1144.       MPTVAL=IVASTR
  1145.       IE=0
  1146.       DO      IGAU=1,2
  1147.       DO      ICOMP=1,NSTRS
  1148.         IE=IE+1
  1149.         MELVAL=IVAL(ICOMP)
  1150.         IBMN=MIN(IB  ,VELCHE(/2))
  1151.         IGMN=MIN(IGAU,VELCHE(/1))
  1152.         XSTRS(IE)=VELCHE(IGMN,IBMN)
  1153.       enddo     
  1154.       enddo     
  1155. C
  1156. C  ON CALCULE LES FORCES INTERNES   !!!!  a  ameliorer
  1157. C
  1158.       CALL LIAINT(XSTRS,XGENE,XFORC,LRE,NSTRS)
  1159. C
  1160. C  RANGEMENT DANS MELVAL
  1161. C
  1162.       IE=0
  1163.       MPTVAL=IVAFOR
  1164.       DO      IGAU=1,NBNN
  1165.       DO      ICOMP=1,NFOR
  1166.         IE=IE+1
  1167.         MELVAL=IVAL(ICOMP)
  1168.         IBMN=MIN(IB  ,VELCHE(/2))
  1169.         VELCHE(IGAU,IBMN)=XFORC(IE)
  1170.       enddo     
  1171.       enddo     
  1172. 3109  CONTINUE
  1173.       SEGSUP WRK1,WRK3,WRK5
  1174.       GO TO 510
  1175. C_______________________________________________________________________
  1176. C
  1177. C     JOI1 : element de liaison a 2 noeuds (6 ddl par
  1178. C noeuds)
  1179. C_______________________________________________________________________
  1180. C
  1181. 265   CONTINUE
  1182.       NBBB=NBNN
  1183.       NSTN=3
  1184.       LRN =3
  1185.       SEGINI WRK1,WRK3,WRK4
  1186. C
  1187.       KERRE=0
  1188.       DO 3110 IB=1,NBELEM
  1189. C
  1190. C  MISE A ZERO DES FORCES INTERNES
  1191. C
  1192.       CALL ZERO(XFORC,1,LRE)
  1193. C
  1194.       MPTVAL=IVAMAT
  1195.       DO IC=1,NMATT
  1196.         IF(IVAL(IC).NE.0) THEN
  1197.           MELVAL=IVAL(IC)
  1198.           IBMN=MIN(IB,VELCHE(/2))
  1199.           WORK(IC)=VELCHE(1,IBMN)
  1200.         ELSE
  1201.           WORK(IC)=0.D0
  1202.         ENDIF
  1203.       END DO
  1204. C
  1205.       CALL MAPALU(NMATT,WORK,BPSS,IDIM)
  1206. C
  1207. C  ON CHERCHE LES CONTRAINTES -
  1208. C
  1209.       MPTVAL=IVASTR
  1210.  
  1211.       IE=0
  1212.       DO 7110 ICOMP=1,NSTRS
  1213.         IE=IE+1
  1214.         MELVAL=IVAL(ICOMP)
  1215.         IBMN=MIN(IB  ,VELCHE(/2))
  1216.         XSTRS(IE)=VELCHE(1,IBMN)
  1217. 7110  CONTINUE
  1218. C
  1219. C  ON CALCULE LES FORCES INTERNES
  1220. C
  1221.       CALL INTJOI(XSTRS,XFORC,NSTRS)
  1222. C
  1223. C  ON PASSE LES CONTRAINTES DANS LE REPERE GLOBAL
  1224. C
  1225.       IAW1 = 101
  1226.       IAW2=IAW1+LRE
  1227.       CALL JOIGLV(XFORC,BPSS,WORK(IAW1),WORK(IAW2),LRE,IDIM)
  1228. C
  1229. C  RANGEMENT DANS MELVAL
  1230. C
  1231.       IE=0
  1232.       MPTVAL=IVAFOR
  1233.       DO      IGAU=1,NBNN
  1234.       DO      ICOMP=1,NFOR
  1235.         IE=IE+1
  1236.         MELVAL=IVAL(ICOMP)
  1237.         IBMN=MIN(IB  ,VELCHE(/2))
  1238.         VELCHE(IGAU,IBMN)=XFORC(IE)
  1239.       enddo     
  1240.       enddo     
  1241. 3110  CONTINUE
  1242.       SEGSUP WRK1,WRK3,WRK4
  1243.       GO TO 510
  1244. c_______________________________________________________________________
  1245. c
  1246. c  element tuyo
  1247. c_______________________________________________________________________
  1248. c
  1249.  96   CONTINUE
  1250.       NBBB=NBNN
  1251.       SEGINI WRK1,WRK3
  1252. c
  1253.       DO 3096 IB=1,NBELEM
  1254. c
  1255. c  on cherche les coordonnees des noeuds de l elementib
  1256. c
  1257.       CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1258. c
  1259. c  il faudrait aussi modifier le vecteur local de la poutre
  1260. c
  1261. c  mise a zero des forces internes
  1262. c
  1263.       CALL ZERO(XFORC,1,LRE)
  1264. c
  1265. c     rangement des caracteristiques dans work
  1266. c
  1267.       MPTVAL=IVACAR
  1268.       DO 6096 IC=1,NCARR
  1269.         IF (IVAL(IC).NE.0) THEN
  1270.           MELVAL=IVAL(IC)
  1271.           IBMN=MIN(IB,VELCHE(/2))
  1272.           WORK(IC)=VELCHE(1,IBMN)
  1273.         ELSE
  1274.           WORK(IC)=0.D0
  1275.         ENDIF
  1276.  6096 CONTINUE
  1277. c
  1278. c  cas ou on a lu le mot vecteur
  1279. c
  1280. c
  1281. c
  1282. c  cas des tuyaux - on calcule les caracteristiques de la poutre
  1283. c                   equivalente
  1284. c
  1285. *Bizarre ici MELE = 96 ???
  1286.       IF(MELE.EQ.42) THEN
  1287.         CISA=WORK(4)
  1288.         VX  =WORK(5)
  1289.         VY  =WORK(6)
  1290.         VZ  =WORK(7)
  1291.         CALL TUYCAR(WORK,CISA,VX,VY,VZ,KERRE,0)
  1292.         IF (KERRE.EQ.77) THEN
  1293.           CALL ERREUR(77)
  1294.           GOTO 510
  1295.         ENDIF
  1296.       ENDIF
  1297. c
  1298. c  on cherche les contraintes -
  1299. c
  1300.       MPTVAL=IVASTR
  1301.       IE=9
  1302.       DO      IGAU=1,2
  1303.       DO      ICOMP=1,NSTRS
  1304.         IE=IE+1
  1305.         MELVAL=IVAL(ICOMP)
  1306.         IBMN=MIN(IB  ,VELCHE(/2))
  1307.         IGMN=MIN(IGAU,VELCHE(/1))
  1308.         WORK(IE)=VELCHE(IGMN,IBMN)
  1309.       enddo    
  1310.       enddo    
  1311. c
  1312. c  on calcule les forces internes
  1313. c
  1314.       CALL POUBSG(XFORC,WORK,XE,WORK(10),WORK(22),KERRE)
  1315.       IF(KERRE.EQ.0) GO TO 5096
  1316.       INTERR(1)=ISOUS
  1317.       INTERR(2)=IB
  1318.       SEGSUP WRK1,WRK3
  1319.       CALL ERREUR(128)
  1320.       GO TO 510
  1321. 5096  CONTINUE
  1322. c
  1323. c  rangement dans melval
  1324. c
  1325.       IE=0
  1326.       MPTVAL=IVAFOR
  1327.       DO      IGAU=1,NBNN
  1328.       DO      ICOMP=1,NFOR
  1329.         IE=IE+1
  1330.         MELVAL=IVAL(ICOMP)
  1331.         IBMN=MIN(IB  ,VELCHE(/2))
  1332.         VELCHE(IGAU,IBMN)=XFORC(IE)
  1333.       enddo     
  1334.       enddo     
  1335. 3096  CONTINUE
  1336.       SEGSUP WRK1,WRK3
  1337.       GO TO 510
  1338.  
  1339. C_______________________________________________________________________
  1340. C
  1341. C     ELEMENTS ZONE_COHESIVE ZOC2,ZOC3,ZOC4
  1342. C_______________________________________________________________________
  1343. C
  1344.   266 CONTINUE
  1345.  
  1346.       NDIM = 2
  1347.       IF(IFOUR.GT.0) NDIM = 3
  1348.       NBBB=NBNN
  1349.       SEGINI WRK1,WRK2,WRK4
  1350. C
  1351.       DO 3266 IB=1,NBELEM
  1352. C
  1353. C  ON CHERCHE LES COORDONNEES DES NOEUDS DE L ELEMENT IB
  1354. C
  1355.       CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1356. C
  1357. C  MISE A ZERO DES FORCES INTERNES
  1358. C
  1359.       CALL ZERO(XFORC,1,LRE)
  1360. C
  1361. C  BOUCLE SUR LES POINTS DE GAUSS
  1362. C
  1363.       DO 6266 IGAU=1,NBPGAU
  1364. C
  1365.        CALL ZCOLOC(XE,SHPTOT,NBNN,MELE,IFOUR,IGAU,BPSS)
  1366. C
  1367.        CALL BZCO(IGAU,MFR,IFOUR,NIFOUR,XE,BPSS,SHPTOT,
  1368.      .          NSTRS,NBNN,LRE,MELE,SHPWRK,BGENE,DJAC,IERT)
  1369.        IF (IERT.NE.0) THEN
  1370.          INTERR(1)=IB
  1371.          CALL ERREUR(612)
  1372.          GOTO 99266
  1373.        ENDIF
  1374. C
  1375. C  ON CHERCHE LES CONTRAINTES -
  1376. C
  1377.       MPTVAL=IVASTR
  1378.       DO 7266 ICOMP=1,NSTRS
  1379.         MELVAL=IVAL(ICOMP)
  1380.         IGMN=MIN(IGAU,VELCHE(/1))
  1381.         IBMN=MIN(IB  ,VELCHE(/2))
  1382.         XSTRS(ICOMP)=VELCHE(IGMN,IBMN)
  1383.  7266 CONTINUE
  1384. C
  1385. C  ON CALCULE B*EFFORTS
  1386. C
  1387.       DJAC=DJAC*POIGAU(IGAU)
  1388.       CALL BSIG(BGENE,XSTRS,NSTRS,LRE,DJAC,XFORC)
  1389.  6266 CONTINUE
  1390. C
  1391. C  TRAITEMENT DE XFORC ET RANGEMENT DANS MELVAL
  1392. C
  1393.       IE=0
  1394.       MPTVAL=IVAFOR
  1395. C
  1396.       DO      NODE=1,NBNN
  1397.       DO      ICOMP=1,NDIM
  1398.         IE=IE+1
  1399.         MELVAL=IVAL(ICOMP)
  1400.         IBMN=MIN(IB  ,VELCHE(/2))
  1401.         VELCHE(NODE,IBMN)=XFORC(IE)
  1402.       enddo    
  1403.       enddo    
  1404.  3266 CONTINUE
  1405.  
  1406. 99266 CONTINUE
  1407.       SEGSUP WRK1,WRK2,WRK4
  1408.       GOTO 510
  1409. c_______________________________________________________________________
  1410.    99 CONTINUE
  1411.       MOTERR(1:4)=NOMTP(MELE)
  1412.       MOTERR(5:12)='BSIGMA'
  1413.       CALL ERREUR(86)
  1414. *
  1415.   510 CONTINUE
  1416.       RETURN
  1417.       END
  1418.  
  1419.  
  1420.  

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