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

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