Télécharger la plaquette Cast3M


fr en

Télécharger epsi2.eso

Retour à la liste

Numérotation des lignes : 

epsi2

  1. C EPSI2     SOURCE    OF166741  25/02/21    21:16:12     12166          
  2.  
  3.       SUBROUTINE EPSI2(IPMAIL,IPMINT,MELE,IELE,
  4.      &  IVADEP,NBPTEL,LRE,NSTRS,LHOOK,
  5.      &  MFR,NDEP,IPORE,IREPS2,NBPGAU,IVAEPS,UZDPG,RYDPG,RXDPG,IIPDPG,
  6.      &  IDERI,ivamat,ivade2,mate,nmatt,cmate,ngra,noer,kerr)
  7.  
  8. C---------------------------------------------------------------------*
  9. C
  10. C            calcul des deformations
  11. C
  12. C             massif, poreux, joints poreux, incompressibles
  13. C---------------------------------------------------------------------*
  14. C                                                                     *
  15. C   entrees :                                                         *
  16. C   ________                                                          *
  17. C                                                                     *
  18. C        ipmail   pointeur sur un segment  meleme                     *
  19. C        ipmint   pointeur sur un segment minte                       *
  20. C        mele     numero de l'element fini                            *
  21. C        iele     numero geometrique de l'element                     *
  22. C        nbpgau   nombre de point d'integration pour la rigidite      *
  23. C        ivadep   pointeur sur le chamelem de deplacements            *
  24. C        nbptel   nombre de points par element                        *
  25. C        lre      nombre de ddl dans la matrice de rigidite           *
  26. C        nstrs    nombre de composante de contraintes/deformations    *
  27. C                 pour une matrice de hooke                           *
  28. C        lhook    dimension de la matrice de hooke                    *
  29. C        mfr      numero de la formulation de l'element fini          *
  30. C        ndep     nombre de composantes de deplacements               *
  31. C        ipore    nombre de fonctions de forme                        *
  32. C        iresp2   flag pour indiquer si on veut les contraintes       *
  33. C                  de piola-kirchhoff                                 *
  34. C        uzdpg    = deformation au point nsdpge support de la         *
  35. C        rydpf    = deformation plane generalisee                     *
  36. C        rxdpg    =                                                   *
  37. C                                                                     *
  38. C   sorties :                                                         *
  39. C   ________                                                          *
  40. C                                                                     *
  41. C        ivaeps   pointeur sur un segment mptval contenant les        *
  42. C                 les melvals de deformations                         *
  43. C---------------------------------------------------------------------*
  44. C Pour MEMOIRE : Si MELE element incompressible alors MFR = 31
  45. C---------------------------------------------------------------------*
  46.  
  47.       IMPLICIT INTEGER(I-N)
  48.       IMPLICIT REAL*8(A-H,O-Z)
  49.  
  50. -INC PPARAM
  51. -INC CCOPTIO
  52. -INC CCREEL
  53. -INC CCHAMP
  54. -INC CCGEOME
  55.  
  56. -INC SMCOORD
  57. -INC SMCHAML
  58. -INC SMCHPOI
  59. -INC SMELEME
  60. -INC SMINTE
  61.  
  62. -INC TMPTVAL
  63.  
  64.       SEGMENT MWRK1
  65.         REAL*8 DDHOOK(NSTRS,NSTRS),XDDL(LRE),XSTRS(NSTRS)
  66.         REAL*8 XE(3,NBBB),DDHOMU(NSTRS,NSTRS)
  67.         REAL*8 SHPWRK(6,NBNO),BGENE(LHOOK,LRE)
  68.         REAL*8 XE1(3,NBBB),XE2(3,NBBB),xstrs2(NSTRS)
  69.         REAL*8 xjac(3,3),valmat(20)
  70.       ENDSEGMENT
  71.  
  72.       SEGMENT MWRK2
  73.         REAL*8 TENS(9),tentra(9),xddls2(lre)
  74.         REAL*8 BGR(NGRA,LRE),BB(2,NGRA),gradi(ngra),R(ngra),u(ngra)
  75.       ENDSEGMENT
  76.  
  77.       SEGMENT MWRK3
  78.         REAL*8 BPSS(3,3),XEL(3,NBBB)
  79.         REAL*8 XNTH(LPP,LPP),XNTB(LPP,LPP),XNTT(LPP)
  80.       ENDSEGMENT
  81.  
  82.       SEGMENT MWRK5
  83.         REAL*8 XGENE(NSTN,LRN)
  84.       ENDSEGMENT
  85.  
  86.       SEGMENT MTRACE
  87.         REAL*8 TRACE(NBPTEL)
  88.       ENDSEGMENT
  89.  
  90.       CHARACTER*8 CMATE
  91.  
  92.       DIMENSION A(4,60),BBX(3,60),UDPGE(3)
  93.       DIMENSION IN(6),JN(6),ITAB(3,3),PP(4,4)
  94.       real*8 valcar(12),var(3)
  95.       real*8 cobma(lhook)
  96.  
  97.       DATA IN/1,2,3,1,1,2/
  98.       DATA JN/1,2,3,2,3,3/
  99.  
  100.       DATA ITAB(1,1),ITAB(1,2),ITAB(1,3)/1,4,5/
  101.       DATA ITAB(2,1),ITAB(2,2),ITAB(2,3)/4,2,6/
  102.       DATA ITAB(3,1),ITAB(3,2),ITAB(3,3)/5,6,3/
  103.  
  104.       real*8 s(2)
  105.  
  106.       s(1)=0.d0
  107.       s(2)=0.d0
  108.       kerr=0
  109.  
  110.       MWRK1  = 0
  111.       MWRK2  = 0
  112.       MWRK3  = 0
  113.       MWRK5  = 0
  114.       MTRACE = 0
  115.  
  116. C  Introduction du point autour duquel se fait le mouvement
  117. C  de la section en defo plane generalisee
  118. C  IIPDPG = numero du noeud/point support si defini pour le modele
  119. C  NDPGE > 0 si prise en compte du point support
  120.       IF (IIPDPG.GT.0) THEN
  121.         IF (IFOUR.EQ.-3) THEN
  122.           NDPGE=3
  123.           UDPGE(1)=UZDPG
  124.           UDPGE(2)=RYDPG
  125.           UDPGE(3)=RXDPG
  126.           IREF=(IIPDPG-1)*(IDIM+1)
  127.           XDPGE=XCOOR(IREF+1)
  128.           YDPGE=XCOOR(IREF+2)
  129.         ELSE IF (IFOUR.EQ.11) THEN
  130.           NDPGE=2
  131.           UDPGE(1)=UZDPG
  132.           UDPGE(2)=RXDPG
  133.           UDPGE(3)=XZero
  134.           XDPGE=XZero
  135.           YDPGE=XZero
  136.         ELSE IF (IFOUR.EQ. 7 .OR. IFOUR.EQ. 8 .OR. IFOUR.EQ. 9 .OR.
  137.      &           IFOUR.EQ.10 .OR. IFOUR.EQ.14) THEN
  138.           NDPGE=1
  139.           UDPGE(1)=UZDPG
  140.           UDPGE(2)=XZero
  141.           UDPGE(3)=XZero
  142.           XDPGE=XZero
  143.           YDPGE=XZero
  144.         else
  145.           write(ioimp,*) 'EPSI2 : ERREUR NDPGE'
  146.           call erreur(5)
  147.           return
  148.         ENDIF
  149.       ELSE
  150.         NDPGE=0
  151.         UDPGE(1)=UZDPG
  152.         UDPGE(2)=XZero
  153.         UDPGE(3)=XZero
  154.         XDPGE=XZero
  155.         YDPGE=XZero
  156.       ENDIF
  157.  
  158.       MELEME=IPMAIL
  159.       NBNN  =NUM(/1)
  160.       NBELEM=NUM(/2)
  161.  
  162.       NHRM=NIFOUR
  163.       MINTE=IPMINT
  164.       NBBB=NBNN
  165.  
  166. C Petite verification prealable (normalement inutile)
  167.       mptval = IVAEPS
  168.       if (NSTRS.ne.ival(/1)) then
  169.         write(ioimp,*) 'EPSI3 : incoherence NSTRS & IVAEPS'
  170.         call erreur(5)
  171.         return
  172.       endif
  173.       do icomp = 1, NSTRS
  174.         melval = IVAL(ICOMP)
  175.         if (melval.le.0) then
  176.           write(ioimp,*) 'EPSI3 : incoherence IVAEPS ival(',icomp,')=0'
  177.           call erreur(5)
  178.           return
  179.         endif
  180.         if (NBPTEL.NE.melval.velche(/1)) then
  181.           write(ioimp,*) 'EPSI3 : incoherence NSTRS & IVAEPS'
  182.           call erreur(5)
  183.           return
  184.         endif
  185.         if (NBELEM .NE. melval.velche(/2)) then
  186.           write(ioimp,*) 'EPSI3 : incoherence NSTRS & IVAEPS'
  187.           call erreur(5)
  188.           return
  189.         endif
  190.       enddo
  191.  
  192. C_______________________________________________________________________
  193. C
  194. C     numero des etiquettes      :
  195. C     etiquettes de 1 a 98 pour traitement specifique a l element
  196. C     dans la zone specifique a chaque element commencant par :
  197. C     5  continue
  198. C     element 5   etiquettes 1005 2005 3005 4005   ...
  199. C     44 continue
  200. C     element 44  etiquettes 1044 2044 3044 4044   ...
  201. C_______________________________________________________________________
  202. C
  203.       IF(MELE.GE.1.AND.MELE.LE.100) THEN
  204. C            CABL SEG2 SEG3 TRI3 TRI4 TRI6 TRI7 QUA4 QUA5 QUA8
  205.       GOTO (   99,  99,  99,   4,  99,   4,  99,   4,  99,   4
  206. C            QUA9 RAC2 RAC3 CUB8 CU20 PRI6 PR15 LIA3 LIA4 LIA6
  207.      1      ,  99,  99,  99,   4,   4,   4,   4,  99,  99,  99
  208. C            LIA8 MULT TET4 TE10 PYR5 PY13 COQ3  DKT POUT LISP
  209.      2      ,  99,  99,   4,   4,   4,   4,  99,  99,  99,  99
  210. C            FAC3 FAC4 FAC6 FAC8 LTR3 LQU4 LCU8 LPR6 LTE4 LPY5
  211.      3      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  212. C            COQ8 TUYA TUFI COQ2 POI1 BARR RACO LSU2 COQ4 LISM
  213.      4      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  214. C            COF3 RES2 LSU3 LSU4 LICO COQ6 CVS2 CVS3 CVT3 CVT6
  215.      5      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  216. C            CVQ4 CVQ8 THP5 TH13 THP6 TH15 THC8 TH20 ICT3 ICQ4
  217.      6      ,  99,  99,  99,  99,  99,  99,  99,  99,   4,   4
  218. C            ICT6 ICQ8 ICC8 ICT4 ICP6 IC20 IC10 IC15 TRIP QUAP
  219.      7      ,   4,   4,   4,   4,   4,   4,   4,   4,  79,  79
  220. C            CUBP TETP PRIP TIMO JOI2 JOI3 JOT3 JOI4 JOI6 JOI8
  221.      8      ,  79,  79,  79,  99,  99,  99,  99,  99,  99,  99
  222. C            LISC TRIH  DST LIC4 CERC TUYO LSE2 LITU HYT3 HYQ4
  223.      9      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99)
  224. c cccccc
  225.      .      ,MELE
  226.       ELSEIF(MELE.GE.101.AND.MELE.LE.200) THEN
  227. C            HYT4 HYP6 HYC8 TRIS QUAS POIS FOR3 JOP3 JOP6 JOP8
  228.       GOTO (   99,  99,  99,  99,  99,  99,  99,  80,  80,  80
  229. C            POL3 POL4 POL5 POL6 POL7 POL8 POL9 PO10 PO11 PO12
  230.      1      ,   4,   4,   4,   4,   4,   4,   4,   4,   4,   4
  231. C            PO13 PO14 BAR3 BAEX LIA2 QUAH CUBH ROT3 SEF2 TRF3
  232.      2      ,   4,   4,  99,  99,  99,  99,  99,  99,  99,  99
  233. C            QUF4 CUF8 PRF6 TEF4 PYF5 MSE3 MTR6 MQU9 MC27 MP18
  234.      3      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  235. C            MT10 MP14 SEF3 TRF7 QUF9 CF27 PF21 TF15 PF19 SEG6
  236.      4      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  237. C            TR21 QU36 C216 P126 TE56 PY91 TRH6 BSE2 BTR4 BQU5
  238.      5      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  239. C            BCU9 BPR7 BTE5 BPY6 FRO4 SEGS POJS JCT3 JCI4 JGI2
  240.      6      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  241. C            JGT3 JGI4 TRIQ QUAQ CUBQ TETQ PRIQ TRIR QUAR CUBR
  242.      7      ,  99,  99, 173, 173, 173, 173, 173, 173, 173, 173
  243. C            TETR PRIR Q4RI Q8RI JOQ3 JOQ6 JOQ8 JOR3 JOR6 JOR8
  244.      8      , 173, 173,   4,   4, 185, 185, 185, 185, 185, 185
  245. C            T1D2 T1D3 M1D2 M1D3 LC03 LC07 LC09 LC27 LC21 LC15
  246.      9      ,  99,  99,   4,   4,  99,  99,  99,  99,  99,  99)
  247. c cccccc
  248.      .      ,MELE-100
  249.       ELSEIF(MELE.GE.201.AND.MELE.LE.300) THEN
  250. C            LC19 LS03 LS07 LS09 LS27 LS21 LS15 LS19 BS03 BS07
  251.       GOTO (   99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  252. C            BS09 BS27 BS21 BS15 BS19 MC03 MC07 MC09 MC27 MC21
  253.      1      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  254. C            MC15 MC19 M103 M107 M109 M127 M121 M115 M119 MS03
  255.      2      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  256. C            MS07 MS09 MS27 MS21 MS15 MS19 QC03 QC07 QC09 QC27
  257.      3      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  258. C            QC21 QC15 QC19 Q103 Q107 Q109 Q127 Q121 Q115 Q119
  259.      4      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  260. C            QS03 QS07 QS09 QS27 QS21 QS15 QS19 CIFL SURE SHB8
  261.      5      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  262. C            CAF2 CAF3 XQ4R XC8R JOI1 ZCO2 ZCO3 ZCO4 TUY2 TUY3
  263.      6      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  264. C            COS2 COA2 ICY5 IC13 CU27 PR21 TE15 PY19 C20R P15R
  265.      7      ,  99,  99,   4,   4,   4,   4,   4,   4,   4,   4)
  266. c cccccc
  267.      .      ,MELE-200
  268.       ENDIF
  269.       GOTO 99
  270. C
  271. C_______________________________________________________________________
  272. C
  273. C  elements massifs et elements incompressibles MECANIQUE
  274. C_______________________________________________________________________
  275. C
  276.   4   CONTINUE
  277.       IF (MFR.EQ.71 .OR. MFR.EQ.73) GOTO 97173
  278.  
  279. C IDERI <= 2 pour lineaire et quadratique et = 5 pour utilisateur
  280. C ===============================================================
  281.         IF ( IDERI.LE.2.OR.IDERI.EQ.5 ) THEN
  282.  
  283. C Elements massifs en FORMULATION 'MECANIQUE'
  284. C -------------------------------------------
  285.       NBNO=NBNN
  286.       NDDD=NDEP-NDPGE
  287. C
  288. C     Donnees liees a l'element de reference
  289. C
  290.       SEGINI,MWRK1
  291.       IF (Ideri.eq.2) SEGINI,MTRACE
  292. C
  293. C  boucle sur les elements
  294. C
  295.       DO 3004 IB=1,NBELEM
  296. C
  297. C     on cherche les deplacements
  298. C
  299.         MPTVAL=IVADEP
  300.         IE=1
  301.         DO IGAU=1,NBNN
  302.           DO ICOMP=1,NDDD
  303.             MELVAL=IVAL(ICOMP)
  304.             IGMN=MIN(IGAU,VELCHE(/1))
  305.             IBMN=MIN(IB  ,VELCHE(/2))
  306.             XDDL(IE)=VELCHE(IGMN,IBMN)
  307.             IE=IE+1
  308.           ENDDO
  309.         ENDDO
  310.         IF (NDPGE.GT.0) THEN
  311.           DO i=1,NDPGE
  312.             XDDL(IE)=UDPGE(i)
  313.             IE=IE+1
  314.           ENDDO
  315.         ENDIF
  316. C
  317. C     on cherche les coordonnees des noeuds de l element ib
  318. C
  319.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  320.  
  321. C on se met a mi-pas
  322.         if (ideri.eq.5) then
  323.           do iyu=1,nbnn
  324.             i_z = (iyu-1)*nddd
  325.             do i=1,idim
  326.               XE(i,iyu)= xe(i,iyu) + xddl( i + i_z )*0.5D0
  327.             enddo
  328.           enddo
  329.         endif
  330. C
  331. C     boucle sur les points de gauss
  332. C
  333.         ISDJC=0
  334. C
  335. C     Calcul des coeff de modification de b-barre (elts incompres)
  336. C= NOM  :   ICT3, ICQ4, ICT6, ICQ8, ICC8, ICT4, ICP6, IC20, IC10, IC15
  337. C= MELE :    69 ,  70 ,  71 ,  72 ,  73 ,  74 ,  75 ,  76 ,  77 ,  78
  338.         IF (MFR.EQ.31) THEN
  339.           CALL BBCALC(MELE,NBNN,MFR,IDIM,XE,
  340.      1                NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  341.      2                NSTRS,LRE,IFOUR,NHRM,A,BBX,SHPTOT,SHPWRK,
  342.      3                BGENE,XDPGE,YDPGE,PP)
  343.         ENDIF
  344.  
  345.         DO IGAU=1,NBPTEL
  346. C
  347.           CALL BMATST(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  348.      1                MELE,MFR,NBNN,LRE,IFOUR,NSTRS,NHRM,1.D0,XE,
  349.      2                SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE)
  350.  
  351.           IF (DJAC.EQ.0.D0) THEN
  352.             kerr=259
  353.             if (noer.eq.0) THEN
  354.               INTERR(1)=IB
  355.               CALL ERREUR(259)
  356.             endif
  357.             GOTO 9904
  358.           ENDIF
  359.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  360.  
  361. C En cas d'elements incompressibles : BGENE selon la methode B-BARRE
  362.           IF (MFR.EQ.31) THEN
  363.             CALL BBAR(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  364.      &                MELE,NBNN,LRE,IFOUR,NSTRS,XE,DJAC,A,BBX,BGENE)
  365.           ENDIF
  366. C
  367.           CALL BST(BGENE,XDDL,LRE,NSTRS,XSTRS)
  368. C
  369. C     calcul des eps 2
  370. C
  371.           IF (Ideri.eq.2)
  372.      &      CALL BST2(SHPWRK,XDDL,XE,NBNO,IFOUR,XSTRS,TRACE,
  373.      &                IGAU,XDPGE,YDPGE,UDPGE,NHRM)
  374. C
  375. C     remplissage du segment contenant les deformations
  376. C
  377.           MPTVAL=IVAEPS
  378.           DO ICOMP=1,NSTRS
  379.             MELVAL=IVAL(ICOMP)
  380.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  381.           ENDDO
  382. C
  383.         ENDDO
  384. C
  385. C     fin de la boucle sur les points de gauss
  386. C
  387. C**        IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  388.         IF (ISDJC.NE.0.AND.ISDJC.NE.NBPTEL) THEN
  389.           kerr=195
  390.           if (noer.eq.0) then
  391.             INTERR(1)=IB
  392.             CALL ERREUR(195)
  393.            endif
  394.           GOTO 9904
  395.         ENDIF
  396. C
  397. C   correction sur la partie quadratique de la deformation dans le cas
  398. C   des elements incompressibles
  399. C
  400.         IF (Ideri.eq.2) THEN
  401.           IF (MFR.EQ.31) THEN
  402.             CALL BBST2(TRACE,NBPTEL,IFOUR,MELE,POIGAU,QSIGAU,
  403.      &                 ETAGAU,DZEGAU,SHPTOT,NBNO,SHPWRK,XE,PP)
  404.             MPTVAL=IVAEPS
  405.             L=2
  406.             IF (IDIM.EQ.3 .OR. IFOUR.EQ.0) L=3
  407.             DO ICOMP=1,L
  408.               MELVAL=IVAL(ICOMP)
  409.               DO IGAU=1,NBPTEL
  410.                 VELCHE(IGAU,IB)=VELCHE(IGAU,IB)+TRACE(IGAU)
  411.               ENDDO
  412.             ENDDO
  413.           ENDIF
  414.         ENDIF
  415.  
  416.  3004 CONTINUE
  417. C
  418. C    fin de la boucle sur les elements
  419. C
  420.  9904 CONTINUE
  421.  
  422. C ===============================================================
  423. C Cas de la derivee de Truesdell
  424. C ===============================================================
  425.         ELSE IF (IDERI.EQ.3) THEN
  426.  
  427.       NBNO=NBNN
  428.       NDDD=NDEP-NDPGE
  429.       SEGINI,MWRK1
  430. C      IF (IREPS2.EQ.1) SEGINI,MTRACE
  431.  
  432. C       on cherche le max des variations des champs pour savoir s'il faut
  433. C       appeler hookis plusieurs fois
  434.       mptval=IVAMAT
  435.       nbgmat=0
  436.       nelmat=0
  437.       DO IM=1,NMATT
  438.         MELVAL=IVAL(IM)
  439.         IF (MELVAL.NE.0) THEN
  440.           nelmat=Max(nelmat,VELCHE(/2))
  441.           nbgmat=Max(nbgmat,VELCHE(/1))
  442.         ENDIF
  443.         VALMAT(IM) = 0.D0
  444.       ENDDO
  445. C
  446. C  boucle sur les elements
  447. C
  448.       DO 34  IB=1,NBELEM
  449. C
  450. C     on cherche les deplacements
  451. C
  452.         MPTVAL=IVADEP
  453.         IE=1
  454.         DO IGAU=1,NBNN
  455.           DO ICOMP=1,NDDD
  456.             MELVAL=IVAL(ICOMP)
  457.             IGMN=MIN(IGAU,VELCHE(/1))
  458.             IBMN=MIN(IB  ,VELCHE(/2))
  459.             XDDL(IE)=VELCHE(IGMN,IBMN)
  460.             IE=IE+1
  461.           ENDDO
  462.         ENDDO
  463.         IF (NDPGE.GT.0) THEN
  464.           DO i=1,NDPGE
  465.             XDDL(IE)=UDPGE(i)
  466.             IE=IE+1
  467.           ENDDO
  468.         ENDIF
  469. C
  470. C     on cherche  les coordonnees des noeuds de l element ib
  471. C
  472.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  473. C on ajoute aux coordonnees la moitie du deplacement pour faire
  474. C la configuration a mi-pas
  475.         do iyu=1,idim
  476.           i_z = (iyu-1)*nddd
  477.           do i=1,nbnn
  478.              XE(i,iyu)= XE(i,iyu) + xddl(i+i_z)*0.5D0
  479.           enddo
  480.         enddo
  481. C
  482. C     boucle sur les points de gauss
  483. C
  484.         ISDJC=0
  485. C
  486.         DO 54  IGAU=1,NBPTEL
  487.           CALL BMATST(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  488.      1                MELE,MFR,NBNN,LRE,IFOUR,NSTRS,NHRM,1.D0,XE,
  489.      2                SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE)
  490.  
  491.           IF (DJAC.EQ.0.D0) THEN
  492.             kerr=259
  493.             if (noer.eq.0) THEN
  494.               INTERR(1)=IB
  495.               CALL ERREUR(259)
  496.             endif
  497.             GOTO 994
  498.           ELSE IF (DJAC.LT.0.D0) THEN
  499.             ISDJC=ISDJC+1
  500.           ENDIF
  501. C
  502. C       on cherche les matrices de Hooke
  503. C
  504.           if(nbgmat+nelmat.gt.2 . or . ib+igau.eq.2) then
  505.             mptval=ivamat
  506.             DO IM=1,NMATT
  507.               MELVAL=IVAL(IM)
  508.               IF (MELVAL.NE.0) THEN
  509.                 IBMN=MIN(IB  ,VELCHE(/2))
  510.                 IGMN=MIN(IGAU,VELCHE(/1))
  511.                 VALMAT(IM)=VELCHE(IGMN,IBMN)
  512.               ELSE
  513.                 VALMAT(IM)=0.D0
  514.               ENDIF
  515.             ENDDO
  516.             kcas=2
  517.             CALL HOOKIS(VALMAT,VALCAR,VAR,MFR,IB,IGAU,EXCEN,EPAIST,
  518.      +                  INAT,MELE,NPINT,IFOUR,KCAS,NBGMAT,Nelmat,
  519.      +                  S,SECT,LHOOK,DDHOMU,DDHOOK,
  520.      +                  COBMA,XMOB,IRETOU)
  521.           endif
  522.           do i=1,nstrs
  523.             do iyu=1,nstrs
  524.              ddhomu(iyu,i)=ddhook(iyu,i)
  525.             enddo
  526.           enddo
  527.  
  528.           CALL DBST(BGENE,DDHomu,XDDL,LRE,NSTRS,XSTRS)
  529. C xstrs contient la contrainte on va faire pica xstrs zdep05
  530.           DO INO = 1, NBNN
  531.             i_z = (ino-1)*nddd
  532.             DO ID=1,IDIM
  533.               XE1(ID,INO)=XE(ID,INO)
  534.               XE2(ID,INO)=XE(ID,INO)-xddl( id + i_z )*0.5D0
  535.             ENDDO
  536.           ENDDO
  537.           DO IYU=1,3
  538.             DO i=1,3
  539.               XJAC(i,iyu)=0.D0
  540.             enddo
  541.           enddo
  542.           CALL HPRIME(XE1,NBNN,IDIM,SHPtot,IGAU,SHPWRK,DJAC)
  543. C
  544. C       CALCUL DE LA MATRICE  F
  545. C
  546.           DO IF=1,IDIM
  547.             DO IE=1,IDIM
  548.               R1 = 0.D0
  549.               DO ID=1,NBNN
  550.                 R1 = R1 + SHPWRK(IF+1,ID)*XE2(IE,ID)
  551.               ENDDO
  552.               XJAC(IE,IF) = R1
  553.             ENDDO
  554.           ENDDO
  555.         IF(IDIM.EQ.2) THEN
  556.            XJAC(3,3)=1.D0
  557.            IF(IFOUR.EQ.0) THEN
  558. C
  559. CC CAS AXISYMETRIQUE
  560. C
  561.              R1=0.D0
  562.              R2=0.D0
  563.              DO 150 ID=1,NBNN
  564.                R1=R1+SHPWRK(1,ID)*XE1(1,ID)
  565.                R2=R2+SHPWRK(1,ID)*XE2(1,ID)
  566.  150         CONTINUE
  567.              XJAC(3,3)=R2/(R1+1.D-20)
  568.            ENDIF
  569.         ENDIF
  570. CC  CALCUL DE DETERMINANT DE F
  571. C
  572.         IF(IDIM.EQ.3) THEN
  573.           DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3))
  574.           DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3))
  575.           DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2))
  576.         ELSE IF(IDIM.EQ.2) THEN
  577.           DETF=XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1)
  578.           DETF = DETF * XJAC (3,3)
  579.         ENDIF
  580.         DETF=1.D0/(DETF+1.D-20)
  581. C
  582. C     CALCUL DES CONTRAINTES DE CAUCHY
  583. C
  584.         DO ID=1,NSTRS
  585.           IND=IN(ID)
  586.           JND=JN(ID)
  587.           R1=0.D0
  588.           DO IE=1,IDIM
  589.             DO IF=1,IDIM
  590.               ICO=ITAB(IE,IF)
  591.               R1 = R1 + XSTRS(ICO)*XJAC(IND,IE)*XJAC(JND,IF)
  592.             ENDDO
  593.           ENDDO
  594.           XSTRS2(ID)= R1 * DETF
  595.         ENDDO
  596. C
  597. C PEGON  :  ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  598. C
  599.        IF(IDIM.EQ.2) THEN
  600.          xstrs2(3)=xstrs(3)*XJAC(3,3)*XJAC(3,3)*DETF
  601.        ENDIF
  602. C fin du calcul de capi dans dans xstrs2 la contrainte de kirchoff
  603. C on continu en calculant epsi sur config initiale
  604.        DO INO=1,NBNN
  605.          i_z = (ino-1) * nddd
  606.          DO ID=1,IDIM
  607.            XE(ID,INO)=XE2(ID,INO)+xddl( id + i_z )*0.5D0
  608.          ENDDO
  609.        ENDDO
  610. C      inversion loi de hooke
  611.        CALL INVALM(DDHOMU,LHOOK,LHOOK,KERRE,0.D0)
  612.        DO ID=1,LHOOK
  613.           R1 = 0.D0
  614.           DO J=1,LHOOK
  615.             R1 = R1 + DDHOMU(ID,J)*xstrs2(J)
  616.           ENDDO
  617.           xstrs(ID)= R1
  618.        ENDDO
  619. C
  620. C     remplissage du segment contenant les deformations
  621. C
  622.           MPTVAL=IVAEPS
  623.           DO ICOMP=1,NSTRS
  624.             MELVAL=IVAL(ICOMP)
  625.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  626.           ENDDO
  627.    54   continue
  628.  
  629.         IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  630.           kerr=195
  631.           if (noer.eq.0) then
  632.             INTERR(1)=IB
  633.             CALL ERREUR(195)
  634.             GOTO 994
  635.           endif
  636.         ENDIF
  637.    34 CONTINUE
  638.   994 CONTINUE
  639. C fin du truesdell
  640.  
  641. C ===============================================================
  642. C debut du jaumann
  643. C ===============================================================
  644.         ELSE IF (IDERI.EQ.4) THEN
  645.  
  646.       NBNO=NBNN
  647. C*    NDDD=NDEP
  648. C*    IF (IFOUR.EQ.-3) NDDD=NDEP-3
  649.       NDDD=NDEP-NDPGE
  650. C
  651.       SEGINI,MWRK1,MTRACE,MWRK2
  652.  
  653. C  boucle sur les elements
  654. C
  655.       DO 394  IB=1,NBELEM
  656. C
  657. C     on cherche les deplacements
  658. C
  659.         MPTVAL=IVADEP
  660.         IE=1
  661.         DO IGAU=1,NBNN
  662.           DO ICOMP=1,NDDD
  663.             MELVAL=IVAL(ICOMP)
  664.             IGMN=MIN(IGAU,VELCHE(/1))
  665.             IBMN=MIN(IB  ,VELCHE(/2))
  666.             XDDL(IE)=VELCHE(IGMN,IBMN)
  667.             IE=IE+1
  668.           ENDDO
  669.         ENDDO
  670.         IF (NDPGE.GT.0) THEN
  671.           DO i=1,NDPGE
  672.             XDDL(IE)=UDPGE(i)
  673.             IE=IE+1
  674.           ENDDO
  675.         ENDIF
  676. C
  677. C     on cherche  les coordonnees des noeuds de l element ib
  678. C
  679.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  680. C
  681. C   on se met sur la config a mi pas (XE) xe1 est la config initiale
  682. C
  683.         do iyu=1,nbnn
  684.           i_z = (iyu-1)*nddd
  685.           do iou=1,idim
  686.              XE1(iou,iyu)= xe(iou,iyu)
  687.              XE(iou,iyu)= xe(iou,iyu) + xddl( iou+ i_z )*0.5d0
  688.           enddo
  689.         enddo
  690. C
  691. C     boucle sur les points de gauss
  692. C
  693.         ISDJC=0
  694. C
  695. C     Calcul des coeff de modification de b-barre (elts incompres)
  696. C= NOM  :   ICT3, ICQ4, ICT6, ICQ8, ICC8, ICT4, ICP6, IC20, IC10, IC15
  697. C= MELE :    69 ,  70 ,  71 ,  72 ,  73 ,  74 ,  75 ,  76 ,  77 ,  78
  698.         IF (MFR.EQ.31) THEN
  699.           CALL BBCALC(MELE,NBNN,MFR,IDIM,XE,
  700.      1                NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  701.      2                NSTRS,LRE,IFOUR,NHRM,A,BBX,SHPTOT,SHPWRK,
  702.      3                BGENE,XDPGE,YDPGE,PP)
  703.         ENDIF
  704. C
  705.         DO 594  IGAU=1,NBPTEL
  706.           CALL BMATST(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  707.      1                MELE,MFR,NBNN,LRE,IFOUR,NSTRS,NHRM,1.D0,XE,
  708.      2                SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE)
  709.  
  710.           IF (DJAC.EQ.0.D0) THEN
  711.             kerr=259
  712.             if (noer.eq.0) THEN
  713.               INTERR(1)=IB
  714.               CALL ERREUR(259)
  715.             endif
  716.             GOTO 9964
  717.           ENDIF
  718.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  719.  
  720. C En cas d'elements incompressibles : BGENE selon la methode B-BARRE
  721.           IF (MFR.EQ.31) THEN
  722.             CALL BBAR(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  723.      &                MELE,NBNN,LRE,IFOUR,NSTRS,XE,DJAC,A,BBX,BGENE)
  724.           ENDIF
  725. C
  726.           CALL BST(BGENE,XDDL,LRE,NSTRS,XSTRS)
  727. C dans xstrs on a les deformations II sur config mi pas
  728. C on va calculer grad u/2 puis decomposition polaire puis rtens
  729. C on retravaille sur config initiale
  730.          r_z=XZero
  731.          iipdpg=0
  732.          CALL BGRMAS(iGau,NOELE,NBNO,LRE,IFOUR,NGRA,NIFOUR,XE1,
  733.      .                r_z,SHPTOT,SHPWRK,BB,BGR,DJAC,IIPDPG)
  734.          do iou=1,lre
  735.            xddls2(iou)= 0.5D0 * xddl(iou)
  736.          enddo
  737.          CALL BGRDEP(BGR,NGRA,XDDLs2,LRE,GRADI)
  738. C on ajoute l'identite au gradient
  739.          if(idim.eQ.2) then
  740.            gradi(1)=gradi(1)+1.D0
  741.            gradi(4)=gradi(4)+1.D0
  742.          ELSE IF(IDIM.EQ.3) THEN
  743.            gradi(1)=gradi(1)+1.D0
  744.            gradi(5)=gradi(5)+1.D0
  745.            gradi(9)=gradi(9)+1.D0
  746.          ENDIF
  747.  
  748.          CALL POLA2(gradi,R,U,IDIM)
  749. C fait le rtens Rt.A.R on utilise u pour mettre Rt
  750. C et on met le tenseur dans le tableau tens
  751. C attention vu le stockage R est en fait Rt
  752.          if(idim.eq.2) then
  753.            U(1)=r(1)
  754.            u(2)=r(3)
  755.            U(3)=R(2)
  756.            u(4)=R(4)
  757.            tens(1)=xstrs(1)
  758.            tens(2)=xstrs(4)*0.5d0
  759.            tens(3)=xstrs(4)*0.5d0
  760.            tens(4)=xstrs(2)
  761.  
  762.          elseif(idim.eq.3) then
  763.            U(1)=r(1)
  764.            u(2)=r(4)
  765.            U(3)=R(7)
  766.            u(4)=R(2)
  767.            u(5)=r(5)
  768.            u(6)=r(8)
  769.            u(7)=r(3)
  770.            u(8)=r(6)
  771.            u(9)=r(9)
  772.            tens(1)=xstrs(1)
  773.            tens(2)=xstrs(4)*0.5D0
  774.            tens(3)=xstrs(5)*0.5D0
  775.            tens(4)=tens(2)
  776.            tens(5)=xstrs(2)
  777.            tens(6)=xstrs(6)*0.5D0
  778.            tens(7)=tens(3)
  779.            tens(8)=tens(6)
  780.            tens(9)=xstrs(3)
  781.          else
  782.            write(6,*)' idim est ni 2 ni 3 stop'
  783.            stop
  784.          endif
  785.  
  786.          CALL MULMAT(tentra,tens,U,IDIM,IDIM,IDIM)
  787.          CALL MULMAT(tens,R,Tentra,IDIM,IDIM,IDIM)
  788. C  tens contient le nouveau tenseur  on  va remplir xstrs
  789. C en 2 D epzz ne change pas
  790.          if(idim.eq.2) then
  791.            xstrs(1)=tens(1)
  792.            xstrs(2)=tens(4)
  793.            xstrs(4)=tens(2)*2.D0
  794.          else
  795.            xstrs(1)=tens(1)
  796.            xstrs(2)=tens(5)
  797.            xstrs(3)=tens(9)
  798.            xstrs(4)=tens(2)*2.D0
  799.            xstrs(5)=tens(3)*2.D0
  800.            xstrs(6)=tens(6)*2.D0
  801.         endif
  802. C
  803. C     remplissage du segment contenant les deformations
  804. C
  805.           MPTVAL=IVAEPS
  806.           DO ICOMP=1,NSTRS
  807.             MELVAL=IVAL(ICOMP)
  808.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  809.          ENDDO
  810. C
  811.  594   CONTINUE
  812. C
  813. C     fin de la boucle sur les points de gauss
  814. C
  815.         IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  816.           INTERR(1)=IB
  817.          if (noer.eq.1) then
  818.           kerr=195
  819.          else
  820.           CALL ERREUR(195)
  821.           GOTO 9964
  822.          endif
  823.         ENDIF
  824.  
  825.  394   CONTINUE
  826. C
  827. C    fin de la boucle sur les elements
  828. C
  829.  9964  CONTINUE
  830.       endif
  831. C
  832.       GOTO 510
  833.  
  834. C Elements massifs en FORMULATIONs 'ELECTROSTATIQUE' et 'DIFFUSION'
  835. C -----------------------------------------------------------------
  836. 97173 CONTINUE
  837.       NBNO = NBNN
  838.       NDDD = NDEP
  839.       SEGINI,MWRK1
  840. C-------------------------
  841. C Boucle sur les elements
  842. C-------------------------
  843.       DO IEL = 1, NBELEM
  844. C - Recuperation des coordonnees des noeuds de l element IEL
  845.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE)
  846. C - Recuperation des inconnues primales aux noeuds de l element IEL
  847.         MPTVAL = IVADEP
  848.         IE = 1
  849.         DO IGAU = 1, NBNN
  850.           DO ICOMP = 1, NDDD
  851.             MELVAL = IVAL(ICOMP)
  852.             IGMN = MIN(IGAU,VELCHE(/1))
  853.             IEMN = MIN(IEL ,VELCHE(/2))
  854.             XDDL(IE) = VELCHE(IGMN,IEMN)
  855.             IE = IE+1
  856.           ENDDO
  857.         ENDDO
  858. C--  --  --  --  --  --  --  --  --
  859. C - Boucle sur les points de Gauss
  860. C--  --  --  --  --  --  --  --  --
  861.         ISDJC=0
  862.         DO IGAU = 1, NBPTEL
  863. C -- Calcul de la matrice B et du jacobien au point de Gauss IGAU
  864.           IF (MFR.EQ.71) THEN
  865.             CALL BELEC(XE,SHPTOT(1,1,IGAU),NBNN,LHOOK,-1,
  866.      &                 SHPWRK,BGENE,DJAC)
  867.           ELSE IF (MFR.EQ.73) THEN
  868.             CALL BDIFF(XE,SHPTOT(1,1,IGAU),NBNN,LHOOK,-1,
  869.      &                 SHPWRK,BGENE,DJAC)
  870.           ENDIF
  871.           IF (DJAC.EQ.0.D0) THEN
  872.             kerr=259
  873.             if (noer.eq.0) THEN
  874.               INTERR(1)=IEL
  875.               CALL ERREUR(259)
  876.             endif
  877.             GOTO 98173
  878.           ENDIF
  879.           IF (DJAC.LT.0.D0) ISDJC = ISDJC+1
  880.           CALL BST(BGENE,XDDL,LRE,NSTRS,XSTRS)
  881. C -- Remplissage du segment contenant les "deformations" = -grad(.)
  882.           MPTVAL = IVAEPS
  883.           DO ICOMP = 1, NSTRS
  884.             MELVAL = IVAL(ICOMP)
  885.             VELCHE(IGAU,IEL) = XSTRS(ICOMP)
  886.           ENDDO
  887. C--  --  --  --  --  --  --  --  --
  888.         ENDDO
  889. C--  --  --  --  --  --  --  --  --
  890.         IF (ISDJC.NE.0 .AND. ISDJC.NE.NBPGAU) THEN
  891.           kerr=195
  892.           if (noer.eq.0) THEN
  893.             INTERR(1)=IEL
  894.             CALL ERREUR(195)
  895.             GOTO 98173
  896.           endif
  897.         ENDIF
  898. C-------------------------
  899.       ENDDO
  900. C-------------------------
  901. 98173 CONTINUE
  902.       GOTO 510
  903.  
  904. C_______________________________________________________________________
  905. C
  906. C     milieux poreux
  907. C_______________________________________________________________________
  908. C
  909.   79  CONTINUE
  910. C
  911. C    pour ces elements  nbbb = nombre de noeuds
  912. C                       nbno = nombre de fonctions de forme
  913. C
  914.       NBNO=IPORE
  915.       NSTN=1
  916.       LRN=NBNO-NBBB
  917.       LRB=LRE-LRN
  918. C
  919.       SEGINI,MWRK1,MWRK5
  920. C  Initialisation de MTRACE necessaire mais inutilise pour ces elements
  921.       IF (IREPS2.EQ.1) SEGINI MTRACE
  922. C
  923.       DO 3079  IB=1,NBELEM
  924. C
  925. C     on cherche  les coordonnees des noeuds de l element ib
  926. C
  927.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  928. C
  929. C     on cherche les deplacements
  930. C
  931.         MPTVAL=IVADEP
  932.         IE=1
  933.         DO IGAU=1,NBNN
  934.           DO ICOMP=1,NDEP-1
  935.             MELVAL=IVAL(ICOMP)
  936.             IGMN=MIN(IGAU,VELCHE(/1))
  937.             IBMN=MIN(IB  ,VELCHE(/2))
  938.             XDDL(IE)=VELCHE(IGMN,IBMN)
  939.             IE=IE+1
  940.           ENDDO
  941.         ENDDO
  942. C
  943. C     puis  les pressions
  944. C
  945.         MELVAL=IVAL(NDEP)
  946.         IBMN=MIN(IB,VELCHE(/2))
  947.         DO IGAU=1,LRN
  948.           IGAUSO=IBSOM(NSPOS(IELE)+IGAU-1)
  949.           IGMN=MIN(IGAUSO,VELCHE(/1))
  950.           XDDL(IE)=VELCHE(IGMN,IBMN)
  951.           IE=IE+1
  952.         ENDDO
  953. C
  954. C     boucle sur les points de gauss
  955. C
  956.         ISDJC=0
  957. C
  958.         DO 5079 IGAU=1,NBPTEL
  959. C
  960.           CALL BNPORE(IGAU,NBNO,NBBB,LRE,IFOUR,LHOOK,NSTN,NHRM,
  961.      &                1.D0,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,1)
  962. C
  963.           IF (DJAC.EQ.0.D0) THEN
  964.             INTERR(1)=IB
  965.             if (noer.eq.0) CALL ERREUR(259)
  966.             kerr=259
  967.             GOTO 9979
  968.           ENDIF
  969.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  970. C
  971.           CALL BST(BGENE,XDDL,LRE,LHOOK,XSTRS)
  972. C
  973. C     calcul des eps 2
  974. C
  975.           IF (IREPS2.EQ.1)
  976.      &      CALL BST2(SHPWRK,XDDL,XE,NBNN,IFOUR,XSTRS,TRACE,
  977.      &                IGAU,XDPGE,YDPGE,UDPGE,NHRM)
  978. C
  979. C     calcul de la pression
  980. C
  981.           XP=0.D0
  982.           DO ID=1,LRN
  983.             XP=XP+XGENE(1,ID)*XDDL(LRB+ID)
  984.           ENDDO
  985.           XSTRS(NSTRS)=XP
  986. C
  987. C     remplissage du segment contenant les deformations
  988. C
  989.           MPTVAL=IVAEPS
  990.           DO 7079 ICOMP=1,NSTRS
  991.             MELVAL=IVAL(ICOMP)
  992.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  993.  7079     CONTINUE
  994. C
  995.  5079   CONTINUE
  996. C
  997. C     fin de la boucle sur les points de gauss
  998. C
  999.         IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  1000.           INTERR(1)=IB
  1001.          if (noer.eq.1) then
  1002.           kerr=195
  1003.          else
  1004.           CALL ERREUR(195)
  1005.           GOTO 9979
  1006.          endif
  1007.         ENDIF
  1008. C
  1009.  3079 CONTINUE
  1010. C
  1011. C    fin de la boucle sur les elements
  1012. C
  1013.  9979 CONTINUE
  1014. C
  1015.       GOTO 510
  1016. C_______________________________________________________________________
  1017. C
  1018. C     milieux poreux - SUITE
  1019. C_______________________________________________________________________
  1020. C
  1021.  173  CONTINUE
  1022. C
  1023. C    pour ces elements  nbbb = nombre de noeuds
  1024. C                       nbno = nombre de fonctions de forme
  1025. C
  1026.       IF (MELE.GE.173.AND.MELE.LE.177) THEN
  1027.         IDECAP = 2
  1028.       ELSE IF (MELE.GE.178.AND.MELE.LE.182) THEN
  1029.         IDECAP = 3
  1030.       ENDIF
  1031. C
  1032.       NBNO=IPORE
  1033.       NSTN=IDECAP
  1034.       NSTB=4
  1035.       IF(IFOUR.EQ.1.OR.IFOUR.EQ.2) NSTB=6
  1036. C
  1037.       LPP=NBNO-NBBB
  1038.       LRN=IDECAP*LPP
  1039.       LRB=LRE-LRN
  1040. C
  1041.       SEGINI,MWRK1,MWRK5
  1042. C  Initialise de MTRACE necessaire mais inutilise pour cet element
  1043.       IF (IREPS2.EQ.1) SEGINI MTRACE
  1044. C
  1045.       DO 3173  IB=1,NBELEM
  1046. C
  1047. C     on cherche  les coordonnees des noeuds de l element ib
  1048. C
  1049.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1050. C
  1051. C     on cherche les deplacements
  1052. C
  1053.         MPTVAL=IVADEP
  1054.         IE=1
  1055.         DO IGAU=1,NBNN
  1056.           DO ICOMP=1,NDEP-IDECAP
  1057.             MELVAL=IVAL(ICOMP)
  1058.             IGMN=MIN(IGAU,VELCHE(/1))
  1059.             IBMN=MIN(IB  ,VELCHE(/2))
  1060.             XDDL(IE)=VELCHE(IGMN,IBMN)
  1061.             IE=IE+1
  1062.           ENDDO
  1063.         ENDDO
  1064. C
  1065. C     puis  les pressions
  1066. C
  1067.         DO IPR = 1,IDECAP
  1068.           MELVAL=IVAL(NDEP-IDECAP+IPR)
  1069.           DO IGAU=1,LPP
  1070.             IGAUSO=IBSOM(NSPOS(IELE)+IGAU-1)
  1071.             IGMN=MIN(IGAUSO,VELCHE(/1))
  1072.             IBMN=MIN(IB  ,VELCHE(/2))
  1073.             XDDL(IE)=VELCHE(IGMN,IBMN)
  1074.             IE=IE+1
  1075.           ENDDO
  1076.         ENDDO
  1077. C
  1078. C     boucle sur les points de gauss
  1079. C
  1080.         ISDJC=0
  1081. C
  1082.         DO 5173  IGAU=1,NBPTEL
  1083. C
  1084.           CALL BNQORE(IGAU,NBNO,NBBB,LRE,IFOUR,NSTB,NSTN,NHRM,
  1085.      &        1.D0,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,IDECAP,LHOOK,1)
  1086. C
  1087.           IF (DJAC.EQ.0.D0) THEN
  1088.             INTERR(1)=IB
  1089.             if (noer.eq.0) CALL ERREUR(259)
  1090.             kerr=259
  1091.             GOTO 9173
  1092.           ENDIF
  1093.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  1094. C
  1095.           CALL BST(BGENE,XDDL,LRE,LHOOK,XSTRS)
  1096. C
  1097. C     calcul des eps 2
  1098. C
  1099.           IF (IREPS2.EQ.1)
  1100.      &      CALL BST2(SHPWRK,XDDL,XE,NBNN,IFOUR,XSTRS,TRACE,
  1101.      &                IGAU,XDPGE,YDPGE,UDPGE,NHRM)
  1102. C
  1103. C     calcul des  pressions
  1104. C
  1105.           IE=LRB
  1106.           DO IPR=1,IDECAP
  1107.             XP=0.D0
  1108.             IPR1=(IPR-1)*LPP
  1109.             DO ID=1,LPP
  1110.               IE=IE+1
  1111.               XP=XP+XGENE(IPR,ID+IPR1)*XDDL(IE)
  1112.             ENDDO
  1113.             XSTRS(NSTRS-IDECAP+IPR)=XP
  1114.           ENDDO
  1115. C
  1116. C     remplissage du segment contenant les deformations
  1117. C
  1118.           MPTVAL=IVAEPS
  1119.           DO 7173 ICOMP=1,NSTRS
  1120.             MELVAL=IVAL(ICOMP)
  1121.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  1122.  7173     CONTINUE
  1123. C
  1124.  5173    CONTINUE
  1125. C
  1126. C     fin de la boucle sur les points de gauss
  1127. C
  1128.          IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  1129.            INTERR(1)=IB
  1130.          if (noer.eq.1) then
  1131.           kerr=195
  1132.          else
  1133.            CALL ERREUR(195)
  1134.            GOTO 9173
  1135.          endif
  1136.          ENDIF
  1137. C
  1138.  3173 CONTINUE
  1139. C
  1140. C    fin de la boucle sur les elements
  1141. C
  1142.  9173 CONTINUE
  1143. C
  1144.       GOTO 510
  1145.  
  1146. C_______________________________________________________________________
  1147. C
  1148. C     joints poreux
  1149. C_______________________________________________________________________
  1150. C
  1151.   80  CONTINUE
  1152. C
  1153. C    pour ces elements  nbbb = nombre de noeuds
  1154. C                       nbno = nombre de fonctions de forme
  1155. C
  1156.       NBNO=IPORE
  1157.       NSTN=1
  1158.       LRN=(NBNO-NBBB)*3/2
  1159.       LPP = LRN
  1160.       LRB=LRE-LRN
  1161.       NFAC=(3*NBBB-NBNO)/2
  1162. C
  1163.       SEGINI,MWRK1,MWRK3,MWRK5
  1164. C
  1165.       DO 3080  IB=1,NBELEM
  1166. C
  1167. C     on cherche d'abord les deplacements
  1168. C
  1169.         MPTVAL=IVADEP
  1170.         IE=1
  1171.         DO 4180 IGAU=1,NFAC
  1172.           DO 4280 ICOMP=1,NDEP-1
  1173.             MELVAL=IVAL(ICOMP)
  1174.             IGMN=MIN(IGAU,VELCHE(/1))
  1175.             IBMN=MIN(IB  ,VELCHE(/2))
  1176.             XDDL(IE)=VELCHE(IGMN,IBMN)
  1177.             IE=IE+1
  1178.  4280     CONTINUE
  1179.  4180   CONTINUE
  1180. C
  1181. C     puis les pressions
  1182. C
  1183.         MELVAL=IVAL(NDEP)
  1184.         DO 4080 IGAU=1,NBNN
  1185.           DO 4190 INSOM=1,NBSOM(IELE)
  1186.             IF (IGAU.EQ.IBSOM(NSPOS(IELE)+INSOM-1)) GO TO 4191
  1187. 4190      CONTINUE
  1188.           IF (IGAU.GT.NFAC) GO TO 4191
  1189.           GO TO 4080
  1190. 4191      CONTINUE
  1191.           IBMN=MIN(IB  ,VELCHE(/2))
  1192.           IGMN=MIN(IGAU,VELCHE(/1))
  1193.           XDDL(IE)=VELCHE(IGMN,IBMN)
  1194.           IE=IE+1
  1195.  4080   CONTINUE
  1196. C
  1197. C     on cherche  les coordonnees des noeuds de l element ib
  1198. C
  1199.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1200. C
  1201. C     calcul des exes locaux et des coordonnees locales
  1202. C
  1203.         CALL JOPLOC(XE,SHPTOT,NBBB,NBNO,IFOUR,XEL,BPSS)
  1204. C
  1205.         CALL INTDEL(XNTH,XNTB,XNTT,LRN,MELE)
  1206. C
  1207. C     boucle sur les points de gauss
  1208. C
  1209.         ISDJC=0
  1210. C
  1211.         DO 5080 IGAU=1,NBPTEL
  1212. C
  1213.           CALL BNPORJ(IGAU,NBNO,NBBB,LRE,IFOUR,LHOOK,NSTN,XE,XEL,
  1214.      &                SHPTOT,SHPWRK,BPSS,BGENE,XGENE,DJAC,1)
  1215. C
  1216.           IF (DJAC.EQ.0.D0) THEN
  1217.             INTERR(1)=IB
  1218.             if (noer.eq.0) CALL ERREUR(259)
  1219.             kerr=259
  1220.             GOTO 9980
  1221.           ENDIF
  1222.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  1223. C
  1224.           CALL BST(BGENE,XDDL,LRB,LHOOK,XSTRS)
  1225.  
  1226. C
  1227. C     calcul de la pression
  1228. C
  1229.           XP=0.D0
  1230.           DO 4480  ID=1,LRN
  1231.             XP=XP+XNTT(ID)*XGENE(1,ID)*XDDL(LRB+ID)
  1232.  4480     CONTINUE
  1233.           XSTRS(NSTRS)=XP
  1234. C
  1235. C     remplissage du segment contenant les deformations
  1236. C
  1237.           MPTVAL=IVAEPS
  1238.           DO 7080 ICOMP=1,NSTRS
  1239.             MELVAL=IVAL(ICOMP)
  1240.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  1241.  7080     CONTINUE
  1242. C
  1243.  5080   CONTINUE
  1244. C
  1245.         IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  1246.           INTERR(1)=IB
  1247.          if (noer.eq.1) then
  1248.           kerr=195
  1249.          else
  1250.           CALL ERREUR(195)
  1251.           GOTO 9980
  1252.          endif
  1253.         ENDIF
  1254. C
  1255.  3080 CONTINUE
  1256. C
  1257.  9980 CONTINUE
  1258.       GOTO 510
  1259.  
  1260. C_______________________________________________________________________
  1261. C
  1262. C     joints poreux  - SUITE
  1263. C_______________________________________________________________________
  1264. C
  1265.  185  CONTINUE
  1266. C
  1267. C    pour ces elements  nbbb = nombre de noeuds
  1268. C                       nbno = nombre de fonctions de forme
  1269. C
  1270.       IF (MELE.GE.185.AND.MELE.LE.187) THEN
  1271.         IDECAP = 2
  1272.       ELSE IF (MELE.GE.188.AND.MELE.LE.190) THEN
  1273.         IDECAP = 3
  1274.       ENDIF
  1275.  
  1276.       NBNO=IPORE
  1277.       NSTN=IDECAP
  1278.       NSTB=2
  1279.       IF(IFOUR.EQ.1.OR.IFOUR.EQ.2) NSTB=3
  1280.  
  1281.       LPP=(NBNO-NBBB)*3/2
  1282.       LRN=IDECAP*LPP
  1283.       LRB=LRE-LRN
  1284.  
  1285.       NFAC=(3*NBBB-NBNO)/2
  1286.  
  1287.       SEGINI,MWRK1,MWRK3,MWRK5
  1288.  
  1289.       DO 3185  IB=1,NBELEM
  1290. C
  1291. C     on cherche d'abord les deplacements
  1292. C
  1293.         MPTVAL=IVADEP
  1294.         IE=1
  1295.         DO 4185 IGAU=1,NFAC
  1296.           DO 4285 ICOMP=1,NDEP-IDECAP
  1297.             MELVAL=IVAL(ICOMP)
  1298.             IGMN=MIN(IGAU,VELCHE(/1))
  1299.             IBMN=MIN(IB  ,VELCHE(/2))
  1300.             XDDL(IE)=VELCHE(IGMN,IBMN)
  1301.             IE=IE+1
  1302.  4285     CONTINUE
  1303.  4185   CONTINUE
  1304. C
  1305. C     puis les pressions
  1306. C
  1307.         DO 4785 IPR=1,IDECAP
  1308.           MELVAL=IVAL(NDEP-IDECAP+IPR)
  1309.           DO 4085 IGAU=1,NBNN
  1310.             DO 4195 INSOM=1,NBSOM(IELE)
  1311.               IF (IGAU.EQ.IBSOM(NSPOS(IELE)+INSOM-1)) GO TO 4891
  1312. 4195        CONTINUE
  1313.             IF (IGAU.GT.NFAC) GO TO 4891
  1314.             GO TO 4085
  1315. 4891        CONTINUE
  1316.             IBMN=MIN(IB  ,VELCHE(/2))
  1317.             IGMN=MIN(IGAU,VELCHE(/1))
  1318.             XDDL(IE)=VELCHE(IGMN,IBMN)
  1319.             IE=IE+1
  1320.  4085     CONTINUE
  1321. 4785    CONTINUE
  1322. C
  1323. C     on cherche  les coordonnees des noeuds de l element ib
  1324. C
  1325.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1326. C
  1327. C     calcul des exes locaux et des coordonnees locales
  1328. C
  1329.         CALL JOPLOC(XE,SHPTOT,NBBB,NBNO,IFOUR,XEL,BPSS)
  1330. C
  1331.         CALL INTDEL(XNTH,XNTB,XNTT,LPP,MELE)
  1332. C
  1333. C     boucle sur les points de gauss
  1334. C
  1335.         ISDJC=0
  1336. C
  1337.         DO 5185  IGAU=1,NBPTEL
  1338. C
  1339.           CALL BNPQRJ(IGAU,NBNO,NBBB,LRE,IFOUR,LHOOK,NSTN,XE,XEL,
  1340.      &           SHPTOT,SHPWRK,BPSS,BGENE,XGENE,DJAC,IDECAP,NSTB,1)
  1341. C
  1342.           IF (DJAC.EQ.0.D0) THEN
  1343.             INTERR(1)=IB
  1344.             if (noer.eq.0) CALL ERREUR(259)
  1345.             kerr=259
  1346.             GOTO 9985
  1347.           ENDIF
  1348.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  1349. C
  1350.           CALL BST(BGENE,XDDL,LRB,LHOOK,XSTRS)
  1351. C
  1352. C     calcul de la pression
  1353. C
  1354.           IE=LRB
  1355.           DO 4985 IPR=1,IDECAP
  1356.             XP=0.D0
  1357.             IPR1=(IPR-1)*LPP
  1358.             DO 4485  ID=1,LPP
  1359.               IE=IE+1
  1360.               XP=XP+XNTT(ID)*XGENE(IPR,ID+IPR1)*XDDL(IE)
  1361.  4485       CONTINUE
  1362.             XSTRS(NSTRS-IDECAP+IPR)=XP
  1363.  4985     CONTINUE
  1364. C
  1365. C     remplissage du segment contenant les deformations
  1366. C
  1367.           MPTVAL=IVAEPS
  1368.           DO 7185 ICOMP=1,NSTRS
  1369.             MELVAL=IVAL(ICOMP)
  1370.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  1371.  7185     CONTINUE
  1372. C
  1373.  5185   CONTINUE
  1374. C
  1375.         IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  1376.           kerr=195
  1377.           INTERR(1)=IB
  1378.          if (noer.eq.0) then
  1379.           CALL ERREUR(195)
  1380.           GOTO 9985
  1381.          endif
  1382.         ENDIF
  1383. C
  1384.  3185 CONTINUE
  1385. C
  1386.  9985 CONTINUE
  1387. C
  1388.       GOTO 510
  1389. C____________________________________________________________________
  1390.  99   CONTINUE
  1391.       MOTERR(1:4)=NOMTP(MELE)
  1392.       MOTERR(9:12)='EPSI'
  1393.       CALL ERREUR(86)
  1394.  
  1395.  510  CONTINUE
  1396.       SEGSUP,MWRK1
  1397.       IF (MWRK2.NE.0)  SEGSUP,MWRK2
  1398.       IF (MWRK3.NE.0)  SEGSUP,MWRK3
  1399.       IF (MWRK5.NE.0)  SEGSUP,MWRK5
  1400.       IF (MTRACE.NE.0) SEGSUP MTRACE
  1401.  
  1402. c      RETURN
  1403.       END
  1404.  
  1405.  
  1406.  
  1407.  

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