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    SP204843  25/07/03    21:15:05     12308          
  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,NOER)
  343.           IF (NOER.NE.0) THEN
  344.             CALL ERREUR(noer)
  345.             RETURN
  346.           ENDIF
  347.         ENDIF
  348.  
  349.         DO IGAU=1,NBPTEL
  350. C
  351.           CALL BMATST(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  352.      1                MELE,MFR,NBNN,LRE,IFOUR,NSTRS,NHRM,1.D0,XE,
  353.      2                SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE)
  354.  
  355.           IF (DJAC.EQ.0.D0) THEN
  356.             kerr=259
  357.             if (noer.eq.0) THEN
  358.               INTERR(1)=IB
  359.               CALL ERREUR(259)
  360.             endif
  361.             GOTO 9904
  362.           ENDIF
  363.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  364.  
  365. C En cas d'elements incompressibles : BGENE selon la methode B-BARRE
  366.           IF (MFR.EQ.31) THEN
  367.             CALL BBAR(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  368.      &                MELE,NBNN,LRE,IFOUR,NSTRS,XE,DJAC,A,BBX,BGENE)
  369.           ENDIF
  370. C
  371.           CALL BST(BGENE,XDDL,LRE,NSTRS,XSTRS)
  372. C
  373. C     calcul des eps 2
  374. C
  375.           IF (Ideri.eq.2)
  376.      &      CALL BST2(SHPWRK,XDDL,XE,NBNO,IFOUR,XSTRS,TRACE,
  377.      &                IGAU,XDPGE,YDPGE,UDPGE,NHRM)
  378. C
  379. C     remplissage du segment contenant les deformations
  380. C
  381.           MPTVAL=IVAEPS
  382.           DO ICOMP=1,NSTRS
  383.             MELVAL=IVAL(ICOMP)
  384.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  385.           ENDDO
  386. C
  387.         ENDDO
  388. C
  389. C     fin de la boucle sur les points de gauss
  390. C
  391. C**        IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  392.         IF (ISDJC.NE.0.AND.ISDJC.NE.NBPTEL) THEN
  393.           kerr=195
  394.           if (noer.eq.0) then
  395.             INTERR(1)=IB
  396.             CALL ERREUR(195)
  397.            endif
  398.           GOTO 9904
  399.         ENDIF
  400. C
  401. C   correction sur la partie quadratique de la deformation dans le cas
  402. C   des elements incompressibles
  403. C
  404.         IF (Ideri.eq.2) THEN
  405.           IF (MFR.EQ.31) THEN
  406.             CALL BBST2(TRACE,NBPTEL,IFOUR,MELE,POIGAU,QSIGAU,
  407.      &                 ETAGAU,DZEGAU,SHPTOT,NBNO,SHPWRK,XE,PP)
  408.             MPTVAL=IVAEPS
  409.             L=2
  410.             IF (IDIM.EQ.3 .OR. IFOUR.EQ.0) L=3
  411.             DO ICOMP=1,L
  412.               MELVAL=IVAL(ICOMP)
  413.               DO IGAU=1,NBPTEL
  414.                 VELCHE(IGAU,IB)=VELCHE(IGAU,IB)+TRACE(IGAU)
  415.               ENDDO
  416.             ENDDO
  417.           ENDIF
  418.         ENDIF
  419.  
  420.  3004 CONTINUE
  421. C
  422. C    fin de la boucle sur les elements
  423. C
  424.  9904 CONTINUE
  425.  
  426. C ===============================================================
  427. C Cas de la derivee de Truesdell
  428. C ===============================================================
  429.         ELSE IF (IDERI.EQ.3) THEN
  430.  
  431.       NBNO=NBNN
  432.       NDDD=NDEP-NDPGE
  433.       SEGINI,MWRK1
  434. C      IF (IREPS2.EQ.1) SEGINI,MTRACE
  435.  
  436. C       on cherche le max des variations des champs pour savoir s'il faut
  437. C       appeler hookis plusieurs fois
  438.       mptval=IVAMAT
  439.       nbgmat=0
  440.       nelmat=0
  441.       DO IM=1,NMATT
  442.         MELVAL=IVAL(IM)
  443.         IF (MELVAL.NE.0) THEN
  444.           nelmat=Max(nelmat,VELCHE(/2))
  445.           nbgmat=Max(nbgmat,VELCHE(/1))
  446.         ENDIF
  447.         VALMAT(IM) = 0.D0
  448.       ENDDO
  449. C
  450. C  boucle sur les elements
  451. C
  452.       DO 34  IB=1,NBELEM
  453. C
  454. C     on cherche les deplacements
  455. C
  456.         MPTVAL=IVADEP
  457.         IE=1
  458.         DO IGAU=1,NBNN
  459.           DO ICOMP=1,NDDD
  460.             MELVAL=IVAL(ICOMP)
  461.             IGMN=MIN(IGAU,VELCHE(/1))
  462.             IBMN=MIN(IB  ,VELCHE(/2))
  463.             XDDL(IE)=VELCHE(IGMN,IBMN)
  464.             IE=IE+1
  465.           ENDDO
  466.         ENDDO
  467.         IF (NDPGE.GT.0) THEN
  468.           DO i=1,NDPGE
  469.             XDDL(IE)=UDPGE(i)
  470.             IE=IE+1
  471.           ENDDO
  472.         ENDIF
  473. C
  474. C     on cherche  les coordonnees des noeuds de l element ib
  475. C
  476.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  477. C on ajoute aux coordonnees la moitie du deplacement pour faire
  478. C la configuration a mi-pas
  479.         do iyu=1,idim
  480.           i_z = (iyu-1)*nddd
  481.           do i=1,nbnn
  482.              XE(i,iyu)= XE(i,iyu) + xddl(i+i_z)*0.5D0
  483.           enddo
  484.         enddo
  485. C
  486. C     boucle sur les points de gauss
  487. C
  488.         ISDJC=0
  489. C
  490.         DO 54  IGAU=1,NBPTEL
  491.           CALL BMATST(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  492.      1                MELE,MFR,NBNN,LRE,IFOUR,NSTRS,NHRM,1.D0,XE,
  493.      2                SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE)
  494.  
  495.           IF (DJAC.EQ.0.D0) THEN
  496.             kerr=259
  497.             if (noer.eq.0) THEN
  498.               INTERR(1)=IB
  499.               CALL ERREUR(259)
  500.             endif
  501.             GOTO 994
  502.           ELSE IF (DJAC.LT.0.D0) THEN
  503.             ISDJC=ISDJC+1
  504.           ENDIF
  505. C
  506. C       on cherche les matrices de Hooke
  507. C
  508.           if(nbgmat+nelmat.gt.2 . or . ib+igau.eq.2) then
  509.             mptval=ivamat
  510.             DO IM=1,NMATT
  511.               MELVAL=IVAL(IM)
  512.               IF (MELVAL.NE.0) THEN
  513.                 IBMN=MIN(IB  ,VELCHE(/2))
  514.                 IGMN=MIN(IGAU,VELCHE(/1))
  515.                 VALMAT(IM)=VELCHE(IGMN,IBMN)
  516.               ELSE
  517.                 VALMAT(IM)=0.D0
  518.               ENDIF
  519.             ENDDO
  520.             kcas=2
  521.             CALL HOOKIS(VALMAT,VALCAR,VAR,MFR,IB,IGAU,EXCEN,EPAIST,
  522.      +                  INAT,MELE,NPINT,IFOUR,KCAS,NBGMAT,Nelmat,
  523.      +                  S,SECT,LHOOK,DDHOMU,DDHOOK,
  524.      +                  COBMA,XMOB,IRETOU)
  525.           endif
  526.           do i=1,nstrs
  527.             do iyu=1,nstrs
  528.              ddhomu(iyu,i)=ddhook(iyu,i)
  529.             enddo
  530.           enddo
  531.  
  532.           CALL DBST(BGENE,DDHomu,XDDL,LRE,NSTRS,XSTRS)
  533. C xstrs contient la contrainte on va faire pica xstrs zdep05
  534.           DO INO = 1, NBNN
  535.             i_z = (ino-1)*nddd
  536.             DO ID=1,IDIM
  537.               XE1(ID,INO)=XE(ID,INO)
  538.               XE2(ID,INO)=XE(ID,INO)-xddl( id + i_z )*0.5D0
  539.             ENDDO
  540.           ENDDO
  541.           DO IYU=1,3
  542.             DO i=1,3
  543.               XJAC(i,iyu)=0.D0
  544.             enddo
  545.           enddo
  546.           CALL HPRIME(XE1,NBNN,IDIM,SHPtot,IGAU,SHPWRK,DJAC)
  547. C
  548. C       CALCUL DE LA MATRICE  F
  549. C
  550.           DO IF=1,IDIM
  551.             DO IE=1,IDIM
  552.               R1 = 0.D0
  553.               DO ID=1,NBNN
  554.                 R1 = R1 + SHPWRK(IF+1,ID)*XE2(IE,ID)
  555.               ENDDO
  556.               XJAC(IE,IF) = R1
  557.             ENDDO
  558.           ENDDO
  559.         IF(IDIM.EQ.2) THEN
  560.            XJAC(3,3)=1.D0
  561.            IF(IFOUR.EQ.0) THEN
  562. C
  563. CC CAS AXISYMETRIQUE
  564. C
  565.              R1=0.D0
  566.              R2=0.D0
  567.              DO 150 ID=1,NBNN
  568.                R1=R1+SHPWRK(1,ID)*XE1(1,ID)
  569.                R2=R2+SHPWRK(1,ID)*XE2(1,ID)
  570.  150         CONTINUE
  571.              XJAC(3,3)=R2/(R1+1.D-20)
  572.            ENDIF
  573.         ENDIF
  574. CC  CALCUL DE DETERMINANT DE F
  575. C
  576.         IF(IDIM.EQ.3) THEN
  577.           DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3))
  578.           DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3))
  579.           DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2))
  580.         ELSE IF(IDIM.EQ.2) THEN
  581.           DETF=XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1)
  582.           DETF = DETF * XJAC (3,3)
  583.         ENDIF
  584.         DETF=1.D0/(DETF+1.D-20)
  585. C
  586. C     CALCUL DES CONTRAINTES DE CAUCHY
  587. C
  588.         DO ID=1,NSTRS
  589.           IND=IN(ID)
  590.           JND=JN(ID)
  591.           R1=0.D0
  592.           DO IE=1,IDIM
  593.             DO IF=1,IDIM
  594.               ICO=ITAB(IE,IF)
  595.               R1 = R1 + XSTRS(ICO)*XJAC(IND,IE)*XJAC(JND,IF)
  596.             ENDDO
  597.           ENDDO
  598.           XSTRS2(ID)= R1 * DETF
  599.         ENDDO
  600. C
  601. C PEGON  :  ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE
  602. C
  603.        IF(IDIM.EQ.2) THEN
  604.          xstrs2(3)=xstrs(3)*XJAC(3,3)*XJAC(3,3)*DETF
  605.        ENDIF
  606. C fin du calcul de capi dans dans xstrs2 la contrainte de kirchoff
  607. C on continu en calculant epsi sur config initiale
  608.        DO INO=1,NBNN
  609.          i_z = (ino-1) * nddd
  610.          DO ID=1,IDIM
  611.            XE(ID,INO)=XE2(ID,INO)+xddl( id + i_z )*0.5D0
  612.          ENDDO
  613.        ENDDO
  614. C      inversion loi de hooke
  615.        CALL INVALM(DDHOMU,LHOOK,LHOOK,KERRE,0.D0)
  616.        DO ID=1,LHOOK
  617.           R1 = 0.D0
  618.           DO J=1,LHOOK
  619.             R1 = R1 + DDHOMU(ID,J)*xstrs2(J)
  620.           ENDDO
  621.           xstrs(ID)= R1
  622.        ENDDO
  623. C
  624. C     remplissage du segment contenant les deformations
  625. C
  626.           MPTVAL=IVAEPS
  627.           DO ICOMP=1,NSTRS
  628.             MELVAL=IVAL(ICOMP)
  629.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  630.           ENDDO
  631.    54   continue
  632.  
  633.         IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  634.           kerr=195
  635.           if (noer.eq.0) then
  636.             INTERR(1)=IB
  637.             CALL ERREUR(195)
  638.             GOTO 994
  639.           endif
  640.         ENDIF
  641.    34 CONTINUE
  642.   994 CONTINUE
  643. C fin du truesdell
  644.  
  645. C ===============================================================
  646. C debut du jaumann
  647. C ===============================================================
  648.         ELSE IF (IDERI.EQ.4) THEN
  649.  
  650.       NBNO=NBNN
  651. C*    NDDD=NDEP
  652. C*    IF (IFOUR.EQ.-3) NDDD=NDEP-3
  653.       NDDD=NDEP-NDPGE
  654. C
  655.       SEGINI,MWRK1,MTRACE,MWRK2
  656.  
  657. C  boucle sur les elements
  658. C
  659.       DO 394  IB=1,NBELEM
  660. C
  661. C     on cherche les deplacements
  662. C
  663.         MPTVAL=IVADEP
  664.         IE=1
  665.         DO IGAU=1,NBNN
  666.           DO ICOMP=1,NDDD
  667.             MELVAL=IVAL(ICOMP)
  668.             IGMN=MIN(IGAU,VELCHE(/1))
  669.             IBMN=MIN(IB  ,VELCHE(/2))
  670.             XDDL(IE)=VELCHE(IGMN,IBMN)
  671.             IE=IE+1
  672.           ENDDO
  673.         ENDDO
  674.         IF (NDPGE.GT.0) THEN
  675.           DO i=1,NDPGE
  676.             XDDL(IE)=UDPGE(i)
  677.             IE=IE+1
  678.           ENDDO
  679.         ENDIF
  680. C
  681. C     on cherche  les coordonnees des noeuds de l element ib
  682. C
  683.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  684. C
  685. C   on se met sur la config a mi pas (XE) xe1 est la config initiale
  686. C
  687.         do iyu=1,nbnn
  688.           i_z = (iyu-1)*nddd
  689.           do iou=1,idim
  690.              XE1(iou,iyu)= xe(iou,iyu)
  691.              XE(iou,iyu)= xe(iou,iyu) + xddl( iou+ i_z )*0.5d0
  692.           enddo
  693.         enddo
  694. C
  695. C     boucle sur les points de gauss
  696. C
  697.         ISDJC=0
  698. C
  699. C     Calcul des coeff de modification de b-barre (elts incompres)
  700. C= NOM  :   ICT3, ICQ4, ICT6, ICQ8, ICC8, ICT4, ICP6, IC20, IC10, IC15
  701. C= MELE :    69 ,  70 ,  71 ,  72 ,  73 ,  74 ,  75 ,  76 ,  77 ,  78
  702.         IF (MFR.EQ.31) THEN
  703.           CALL BBCALC(MELE,NBNN,MFR,IDIM,XE,
  704.      1                NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  705.      2                NSTRS,LRE,IFOUR,NHRM,A,BBX,SHPTOT,SHPWRK,
  706.      3                BGENE,XDPGE,YDPGE,PP,NOER)
  707.           IF (NOER.NE.0) THEN
  708.             CALL ERREUR(noer)
  709.             RETURN
  710.           ENDIF
  711.         ENDIF
  712. C
  713.         DO 594  IGAU=1,NBPTEL
  714.           CALL BMATST(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  715.      1                MELE,MFR,NBNN,LRE,IFOUR,NSTRS,NHRM,1.D0,XE,
  716.      2                SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE)
  717.  
  718.           IF (DJAC.EQ.0.D0) THEN
  719.             kerr=259
  720.             if (noer.eq.0) THEN
  721.               INTERR(1)=IB
  722.               CALL ERREUR(259)
  723.             endif
  724.             GOTO 9964
  725.           ENDIF
  726.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  727.  
  728. C En cas d'elements incompressibles : BGENE selon la methode B-BARRE
  729.           IF (MFR.EQ.31) THEN
  730.             CALL BBAR(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  731.      &                MELE,NBNN,LRE,IFOUR,NSTRS,XE,DJAC,A,BBX,BGENE)
  732.           ENDIF
  733. C
  734.           CALL BST(BGENE,XDDL,LRE,NSTRS,XSTRS)
  735. C dans xstrs on a les deformations II sur config mi pas
  736. C on va calculer grad u/2 puis decomposition polaire puis rtens
  737. C on retravaille sur config initiale
  738.          r_z=XZero
  739.          iipdpg=0
  740.          CALL BGRMAS(iGau,NOELE,NBNO,LRE,IFOUR,NGRA,NIFOUR,XE1,
  741.      .                r_z,SHPTOT,SHPWRK,BB,BGR,DJAC,IIPDPG)
  742.          do iou=1,lre
  743.            xddls2(iou)= 0.5D0 * xddl(iou)
  744.          enddo
  745.          CALL BGRDEP(BGR,NGRA,XDDLs2,LRE,GRADI)
  746. C on ajoute l'identite au gradient
  747.          if(idim.eQ.2) then
  748.            gradi(1)=gradi(1)+1.D0
  749.            gradi(4)=gradi(4)+1.D0
  750.          ELSE IF(IDIM.EQ.3) THEN
  751.            gradi(1)=gradi(1)+1.D0
  752.            gradi(5)=gradi(5)+1.D0
  753.            gradi(9)=gradi(9)+1.D0
  754.          ENDIF
  755.  
  756.          CALL POLA2(gradi,R,U,IDIM)
  757. C fait le rtens Rt.A.R on utilise u pour mettre Rt
  758. C et on met le tenseur dans le tableau tens
  759. C attention vu le stockage R est en fait Rt
  760.          if(idim.eq.2) then
  761.            U(1)=r(1)
  762.            u(2)=r(3)
  763.            U(3)=R(2)
  764.            u(4)=R(4)
  765.            tens(1)=xstrs(1)
  766.            tens(2)=xstrs(4)*0.5d0
  767.            tens(3)=xstrs(4)*0.5d0
  768.            tens(4)=xstrs(2)
  769.  
  770.          elseif(idim.eq.3) then
  771.            U(1)=r(1)
  772.            u(2)=r(4)
  773.            U(3)=R(7)
  774.            u(4)=R(2)
  775.            u(5)=r(5)
  776.            u(6)=r(8)
  777.            u(7)=r(3)
  778.            u(8)=r(6)
  779.            u(9)=r(9)
  780.            tens(1)=xstrs(1)
  781.            tens(2)=xstrs(4)*0.5D0
  782.            tens(3)=xstrs(5)*0.5D0
  783.            tens(4)=tens(2)
  784.            tens(5)=xstrs(2)
  785.            tens(6)=xstrs(6)*0.5D0
  786.            tens(7)=tens(3)
  787.            tens(8)=tens(6)
  788.            tens(9)=xstrs(3)
  789.          else
  790.            write(6,*)' idim est ni 2 ni 3 stop'
  791.            stop
  792.          endif
  793.  
  794.          CALL MULMAT(tentra,tens,U,IDIM,IDIM,IDIM)
  795.          CALL MULMAT(tens,R,Tentra,IDIM,IDIM,IDIM)
  796. C  tens contient le nouveau tenseur  on  va remplir xstrs
  797. C en 2 D epzz ne change pas
  798.          if(idim.eq.2) then
  799.            xstrs(1)=tens(1)
  800.            xstrs(2)=tens(4)
  801.            xstrs(4)=tens(2)*2.D0
  802.          else
  803.            xstrs(1)=tens(1)
  804.            xstrs(2)=tens(5)
  805.            xstrs(3)=tens(9)
  806.            xstrs(4)=tens(2)*2.D0
  807.            xstrs(5)=tens(3)*2.D0
  808.            xstrs(6)=tens(6)*2.D0
  809.         endif
  810. C
  811. C     remplissage du segment contenant les deformations
  812. C
  813.           MPTVAL=IVAEPS
  814.           DO ICOMP=1,NSTRS
  815.             MELVAL=IVAL(ICOMP)
  816.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  817.          ENDDO
  818. C
  819.  594   CONTINUE
  820. C
  821. C     fin de la boucle sur les points de gauss
  822. C
  823.         IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  824.           INTERR(1)=IB
  825.          if (noer.eq.1) then
  826.           kerr=195
  827.          else
  828.           CALL ERREUR(195)
  829.           GOTO 9964
  830.          endif
  831.         ENDIF
  832.  
  833.  394   CONTINUE
  834. C
  835. C    fin de la boucle sur les elements
  836. C
  837.  9964  CONTINUE
  838.       endif
  839. C
  840.       GOTO 510
  841.  
  842. C Elements massifs en FORMULATIONs 'ELECTROSTATIQUE' et 'DIFFUSION'
  843. C -----------------------------------------------------------------
  844. 97173 CONTINUE
  845.       NBNO = NBNN
  846.       NDDD = NDEP
  847.       SEGINI,MWRK1
  848. C-------------------------
  849. C Boucle sur les elements
  850. C-------------------------
  851.       DO IEL = 1, NBELEM
  852. C - Recuperation des coordonnees des noeuds de l element IEL
  853.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE)
  854. C - Recuperation des inconnues primales aux noeuds de l element IEL
  855.         MPTVAL = IVADEP
  856.         IE = 1
  857.         DO IGAU = 1, NBNN
  858.           DO ICOMP = 1, NDDD
  859.             MELVAL = IVAL(ICOMP)
  860.             IGMN = MIN(IGAU,VELCHE(/1))
  861.             IEMN = MIN(IEL ,VELCHE(/2))
  862.             XDDL(IE) = VELCHE(IGMN,IEMN)
  863.             IE = IE+1
  864.           ENDDO
  865.         ENDDO
  866. C--  --  --  --  --  --  --  --  --
  867. C - Boucle sur les points de Gauss
  868. C--  --  --  --  --  --  --  --  --
  869.         ISDJC=0
  870.         DO IGAU = 1, NBPTEL
  871. C -- Calcul de la matrice B et du jacobien au point de Gauss IGAU
  872.           IF (MFR.EQ.71) THEN
  873.             CALL BELEC(XE,SHPTOT(1,1,IGAU),NBNN,LHOOK,-1,
  874.      &                 SHPWRK,BGENE,DJAC)
  875.           ELSE IF (MFR.EQ.73) THEN
  876.             CALL BDIFF(XE,SHPTOT(1,1,IGAU),NBNN,LHOOK,-1,
  877.      &                 SHPWRK,BGENE,DJAC)
  878.           ENDIF
  879.           IF (DJAC.EQ.0.D0) THEN
  880.             kerr=259
  881.             if (noer.eq.0) THEN
  882.               INTERR(1)=IEL
  883.               CALL ERREUR(259)
  884.             endif
  885.             GOTO 98173
  886.           ENDIF
  887.           IF (DJAC.LT.0.D0) ISDJC = ISDJC+1
  888.           CALL BST(BGENE,XDDL,LRE,NSTRS,XSTRS)
  889. C -- Remplissage du segment contenant les "deformations" = -grad(.)
  890.           MPTVAL = IVAEPS
  891.           DO ICOMP = 1, NSTRS
  892.             MELVAL = IVAL(ICOMP)
  893.             VELCHE(IGAU,IEL) = XSTRS(ICOMP)
  894.           ENDDO
  895. C--  --  --  --  --  --  --  --  --
  896.         ENDDO
  897. C--  --  --  --  --  --  --  --  --
  898.         IF (ISDJC.NE.0 .AND. ISDJC.NE.NBPGAU) THEN
  899.           kerr=195
  900.           if (noer.eq.0) THEN
  901.             INTERR(1)=IEL
  902.             CALL ERREUR(195)
  903.             GOTO 98173
  904.           endif
  905.         ENDIF
  906. C-------------------------
  907.       ENDDO
  908. C-------------------------
  909. 98173 CONTINUE
  910.       GOTO 510
  911.  
  912. C_______________________________________________________________________
  913. C
  914. C     milieux poreux
  915. C_______________________________________________________________________
  916. C
  917.   79  CONTINUE
  918. C
  919. C    pour ces elements  nbbb = nombre de noeuds
  920. C                       nbno = nombre de fonctions de forme
  921. C
  922.       NBNO=IPORE
  923.       NSTN=1
  924.       LRN=NBNO-NBBB
  925.       LRB=LRE-LRN
  926. C
  927.       SEGINI,MWRK1,MWRK5
  928. C  Initialisation de MTRACE necessaire mais inutilise pour ces elements
  929.       IF (IREPS2.EQ.1) SEGINI MTRACE
  930. C
  931.       DO 3079  IB=1,NBELEM
  932. C
  933. C     on cherche  les coordonnees des noeuds de l element ib
  934. C
  935.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  936. C
  937. C     on cherche les deplacements
  938. C
  939.         MPTVAL=IVADEP
  940.         IE=1
  941.         DO IGAU=1,NBNN
  942.           DO ICOMP=1,NDEP-1
  943.             MELVAL=IVAL(ICOMP)
  944.             IGMN=MIN(IGAU,VELCHE(/1))
  945.             IBMN=MIN(IB  ,VELCHE(/2))
  946.             XDDL(IE)=VELCHE(IGMN,IBMN)
  947.             IE=IE+1
  948.           ENDDO
  949.         ENDDO
  950. C
  951. C     puis  les pressions
  952. C
  953.         MELVAL=IVAL(NDEP)
  954.         IBMN=MIN(IB,VELCHE(/2))
  955.         DO IGAU=1,LRN
  956.           IGAUSO=IBSOM(NSPOS(IELE)+IGAU-1)
  957.           IGMN=MIN(IGAUSO,VELCHE(/1))
  958.           XDDL(IE)=VELCHE(IGMN,IBMN)
  959.           IE=IE+1
  960.         ENDDO
  961. C
  962. C     boucle sur les points de gauss
  963. C
  964.         ISDJC=0
  965. C
  966.         DO 5079 IGAU=1,NBPTEL
  967. C
  968.           CALL BNPORE(IGAU,NBNO,NBBB,LRE,IFOUR,LHOOK,NSTN,NHRM,
  969.      &                1.D0,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,1)
  970. C
  971.           IF (DJAC.EQ.0.D0) THEN
  972.             INTERR(1)=IB
  973.             if (noer.eq.0) CALL ERREUR(259)
  974.             kerr=259
  975.             GOTO 9979
  976.           ENDIF
  977.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  978. C
  979.           CALL BST(BGENE,XDDL,LRE,LHOOK,XSTRS)
  980. C
  981. C     calcul des eps 2
  982. C
  983.           IF (IREPS2.EQ.1)
  984.      &      CALL BST2(SHPWRK,XDDL,XE,NBNN,IFOUR,XSTRS,TRACE,
  985.      &                IGAU,XDPGE,YDPGE,UDPGE,NHRM)
  986. C
  987. C     calcul de la pression
  988. C
  989.           XP=0.D0
  990.           DO ID=1,LRN
  991.             XP=XP+XGENE(1,ID)*XDDL(LRB+ID)
  992.           ENDDO
  993.           XSTRS(NSTRS)=XP
  994. C
  995. C     remplissage du segment contenant les deformations
  996. C
  997.           MPTVAL=IVAEPS
  998.           DO 7079 ICOMP=1,NSTRS
  999.             MELVAL=IVAL(ICOMP)
  1000.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  1001.  7079     CONTINUE
  1002. C
  1003.  5079   CONTINUE
  1004. C
  1005. C     fin de la boucle sur les points de gauss
  1006. C
  1007.         IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  1008.           INTERR(1)=IB
  1009.          if (noer.eq.1) then
  1010.           kerr=195
  1011.          else
  1012.           CALL ERREUR(195)
  1013.           GOTO 9979
  1014.          endif
  1015.         ENDIF
  1016. C
  1017.  3079 CONTINUE
  1018. C
  1019. C    fin de la boucle sur les elements
  1020. C
  1021.  9979 CONTINUE
  1022. C
  1023.       GOTO 510
  1024. C_______________________________________________________________________
  1025. C
  1026. C     milieux poreux - SUITE
  1027. C_______________________________________________________________________
  1028. C
  1029.  173  CONTINUE
  1030. C
  1031. C    pour ces elements  nbbb = nombre de noeuds
  1032. C                       nbno = nombre de fonctions de forme
  1033. C
  1034.       IF (MELE.GE.173.AND.MELE.LE.177) THEN
  1035.         IDECAP = 2
  1036.       ELSE IF (MELE.GE.178.AND.MELE.LE.182) THEN
  1037.         IDECAP = 3
  1038.       ENDIF
  1039. C
  1040.       NBNO=IPORE
  1041.       NSTN=IDECAP
  1042.       NSTB=4
  1043.       IF(IFOUR.EQ.1.OR.IFOUR.EQ.2) NSTB=6
  1044. C
  1045.       LPP=NBNO-NBBB
  1046.       LRN=IDECAP*LPP
  1047.       LRB=LRE-LRN
  1048. C
  1049.       SEGINI,MWRK1,MWRK5
  1050. C  Initialise de MTRACE necessaire mais inutilise pour cet element
  1051.       IF (IREPS2.EQ.1) SEGINI MTRACE
  1052. C
  1053.       DO 3173  IB=1,NBELEM
  1054. C
  1055. C     on cherche  les coordonnees des noeuds de l element ib
  1056. C
  1057.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1058. C
  1059. C     on cherche les deplacements
  1060. C
  1061.         MPTVAL=IVADEP
  1062.         IE=1
  1063.         DO IGAU=1,NBNN
  1064.           DO ICOMP=1,NDEP-IDECAP
  1065.             MELVAL=IVAL(ICOMP)
  1066.             IGMN=MIN(IGAU,VELCHE(/1))
  1067.             IBMN=MIN(IB  ,VELCHE(/2))
  1068.             XDDL(IE)=VELCHE(IGMN,IBMN)
  1069.             IE=IE+1
  1070.           ENDDO
  1071.         ENDDO
  1072. C
  1073. C     puis  les pressions
  1074. C
  1075.         DO IPR = 1,IDECAP
  1076.           MELVAL=IVAL(NDEP-IDECAP+IPR)
  1077.           DO IGAU=1,LPP
  1078.             IGAUSO=IBSOM(NSPOS(IELE)+IGAU-1)
  1079.             IGMN=MIN(IGAUSO,VELCHE(/1))
  1080.             IBMN=MIN(IB  ,VELCHE(/2))
  1081.             XDDL(IE)=VELCHE(IGMN,IBMN)
  1082.             IE=IE+1
  1083.           ENDDO
  1084.         ENDDO
  1085. C
  1086. C     boucle sur les points de gauss
  1087. C
  1088.         ISDJC=0
  1089. C
  1090.         DO 5173  IGAU=1,NBPTEL
  1091. C
  1092.           CALL BNQORE(IGAU,NBNO,NBBB,LRE,IFOUR,NSTB,NSTN,NHRM,
  1093.      &        1.D0,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,IDECAP,LHOOK,1)
  1094. C
  1095.           IF (DJAC.EQ.0.D0) THEN
  1096.             INTERR(1)=IB
  1097.             if (noer.eq.0) CALL ERREUR(259)
  1098.             kerr=259
  1099.             GOTO 9173
  1100.           ENDIF
  1101.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  1102. C
  1103.           CALL BST(BGENE,XDDL,LRE,LHOOK,XSTRS)
  1104. C
  1105. C     calcul des eps 2
  1106. C
  1107.           IF (IREPS2.EQ.1)
  1108.      &      CALL BST2(SHPWRK,XDDL,XE,NBNN,IFOUR,XSTRS,TRACE,
  1109.      &                IGAU,XDPGE,YDPGE,UDPGE,NHRM)
  1110. C
  1111. C     calcul des  pressions
  1112. C
  1113.           IE=LRB
  1114.           DO IPR=1,IDECAP
  1115.             XP=0.D0
  1116.             IPR1=(IPR-1)*LPP
  1117.             DO ID=1,LPP
  1118.               IE=IE+1
  1119.               XP=XP+XGENE(IPR,ID+IPR1)*XDDL(IE)
  1120.             ENDDO
  1121.             XSTRS(NSTRS-IDECAP+IPR)=XP
  1122.           ENDDO
  1123. C
  1124. C     remplissage du segment contenant les deformations
  1125. C
  1126.           MPTVAL=IVAEPS
  1127.           DO 7173 ICOMP=1,NSTRS
  1128.             MELVAL=IVAL(ICOMP)
  1129.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  1130.  7173     CONTINUE
  1131. C
  1132.  5173    CONTINUE
  1133. C
  1134. C     fin de la boucle sur les points de gauss
  1135. C
  1136.          IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  1137.            INTERR(1)=IB
  1138.          if (noer.eq.1) then
  1139.           kerr=195
  1140.          else
  1141.            CALL ERREUR(195)
  1142.            GOTO 9173
  1143.          endif
  1144.          ENDIF
  1145. C
  1146.  3173 CONTINUE
  1147. C
  1148. C    fin de la boucle sur les elements
  1149. C
  1150.  9173 CONTINUE
  1151. C
  1152.       GOTO 510
  1153.  
  1154. C_______________________________________________________________________
  1155. C
  1156. C     joints poreux
  1157. C_______________________________________________________________________
  1158. C
  1159.   80  CONTINUE
  1160. C
  1161. C    pour ces elements  nbbb = nombre de noeuds
  1162. C                       nbno = nombre de fonctions de forme
  1163. C
  1164.       NBNO=IPORE
  1165.       NSTN=1
  1166.       LRN=(NBNO-NBBB)*3/2
  1167.       LPP = LRN
  1168.       LRB=LRE-LRN
  1169.       NFAC=(3*NBBB-NBNO)/2
  1170. C
  1171.       SEGINI,MWRK1,MWRK3,MWRK5
  1172. C
  1173.       DO 3080  IB=1,NBELEM
  1174. C
  1175. C     on cherche d'abord les deplacements
  1176. C
  1177.         MPTVAL=IVADEP
  1178.         IE=1
  1179.         DO 4180 IGAU=1,NFAC
  1180.           DO 4280 ICOMP=1,NDEP-1
  1181.             MELVAL=IVAL(ICOMP)
  1182.             IGMN=MIN(IGAU,VELCHE(/1))
  1183.             IBMN=MIN(IB  ,VELCHE(/2))
  1184.             XDDL(IE)=VELCHE(IGMN,IBMN)
  1185.             IE=IE+1
  1186.  4280     CONTINUE
  1187.  4180   CONTINUE
  1188. C
  1189. C     puis les pressions
  1190. C
  1191.         MELVAL=IVAL(NDEP)
  1192.         DO 4080 IGAU=1,NBNN
  1193.           DO 4190 INSOM=1,NBSOM(IELE)
  1194.             IF (IGAU.EQ.IBSOM(NSPOS(IELE)+INSOM-1)) GO TO 4191
  1195. 4190      CONTINUE
  1196.           IF (IGAU.GT.NFAC) GO TO 4191
  1197.           GO TO 4080
  1198. 4191      CONTINUE
  1199.           IBMN=MIN(IB  ,VELCHE(/2))
  1200.           IGMN=MIN(IGAU,VELCHE(/1))
  1201.           XDDL(IE)=VELCHE(IGMN,IBMN)
  1202.           IE=IE+1
  1203.  4080   CONTINUE
  1204. C
  1205. C     on cherche  les coordonnees des noeuds de l element ib
  1206. C
  1207.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1208. C
  1209. C     calcul des exes locaux et des coordonnees locales
  1210. C
  1211.         CALL JOPLOC(XE,SHPTOT,NBBB,NBNO,IFOUR,XEL,BPSS)
  1212. C
  1213.         CALL INTDEL(XNTH,XNTB,XNTT,LRN,MELE)
  1214. C
  1215. C     boucle sur les points de gauss
  1216. C
  1217.         ISDJC=0
  1218. C
  1219.         DO 5080 IGAU=1,NBPTEL
  1220. C
  1221.           CALL BNPORJ(IGAU,NBNO,NBBB,LRE,IFOUR,LHOOK,NSTN,XE,XEL,
  1222.      &                SHPTOT,SHPWRK,BPSS,BGENE,XGENE,DJAC,1)
  1223. C
  1224.           IF (DJAC.EQ.0.D0) THEN
  1225.             INTERR(1)=IB
  1226.             if (noer.eq.0) CALL ERREUR(259)
  1227.             kerr=259
  1228.             GOTO 9980
  1229.           ENDIF
  1230.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  1231. C
  1232.           CALL BST(BGENE,XDDL,LRB,LHOOK,XSTRS)
  1233.  
  1234. C
  1235. C     calcul de la pression
  1236. C
  1237.           XP=0.D0
  1238.           DO 4480  ID=1,LRN
  1239.             XP=XP+XNTT(ID)*XGENE(1,ID)*XDDL(LRB+ID)
  1240.  4480     CONTINUE
  1241.           XSTRS(NSTRS)=XP
  1242. C
  1243. C     remplissage du segment contenant les deformations
  1244. C
  1245.           MPTVAL=IVAEPS
  1246.           DO 7080 ICOMP=1,NSTRS
  1247.             MELVAL=IVAL(ICOMP)
  1248.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  1249.  7080     CONTINUE
  1250. C
  1251.  5080   CONTINUE
  1252. C
  1253.         IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  1254.           INTERR(1)=IB
  1255.          if (noer.eq.1) then
  1256.           kerr=195
  1257.          else
  1258.           CALL ERREUR(195)
  1259.           GOTO 9980
  1260.          endif
  1261.         ENDIF
  1262. C
  1263.  3080 CONTINUE
  1264. C
  1265.  9980 CONTINUE
  1266.       GOTO 510
  1267.  
  1268. C_______________________________________________________________________
  1269. C
  1270. C     joints poreux  - SUITE
  1271. C_______________________________________________________________________
  1272. C
  1273.  185  CONTINUE
  1274. C
  1275. C    pour ces elements  nbbb = nombre de noeuds
  1276. C                       nbno = nombre de fonctions de forme
  1277. C
  1278.       IF (MELE.GE.185.AND.MELE.LE.187) THEN
  1279.         IDECAP = 2
  1280.       ELSE IF (MELE.GE.188.AND.MELE.LE.190) THEN
  1281.         IDECAP = 3
  1282.       ENDIF
  1283.  
  1284.       NBNO=IPORE
  1285.       NSTN=IDECAP
  1286.       NSTB=2
  1287.       IF(IFOUR.EQ.1.OR.IFOUR.EQ.2) NSTB=3
  1288.  
  1289.       LPP=(NBNO-NBBB)*3/2
  1290.       LRN=IDECAP*LPP
  1291.       LRB=LRE-LRN
  1292.  
  1293.       NFAC=(3*NBBB-NBNO)/2
  1294.  
  1295.       SEGINI,MWRK1,MWRK3,MWRK5
  1296.  
  1297.       DO 3185  IB=1,NBELEM
  1298. C
  1299. C     on cherche d'abord les deplacements
  1300. C
  1301.         MPTVAL=IVADEP
  1302.         IE=1
  1303.         DO 4185 IGAU=1,NFAC
  1304.           DO 4285 ICOMP=1,NDEP-IDECAP
  1305.             MELVAL=IVAL(ICOMP)
  1306.             IGMN=MIN(IGAU,VELCHE(/1))
  1307.             IBMN=MIN(IB  ,VELCHE(/2))
  1308.             XDDL(IE)=VELCHE(IGMN,IBMN)
  1309.             IE=IE+1
  1310.  4285     CONTINUE
  1311.  4185   CONTINUE
  1312. C
  1313. C     puis les pressions
  1314. C
  1315.         DO 4785 IPR=1,IDECAP
  1316.           MELVAL=IVAL(NDEP-IDECAP+IPR)
  1317.           DO 4085 IGAU=1,NBNN
  1318.             DO 4195 INSOM=1,NBSOM(IELE)
  1319.               IF (IGAU.EQ.IBSOM(NSPOS(IELE)+INSOM-1)) GO TO 4891
  1320. 4195        CONTINUE
  1321.             IF (IGAU.GT.NFAC) GO TO 4891
  1322.             GO TO 4085
  1323. 4891        CONTINUE
  1324.             IBMN=MIN(IB  ,VELCHE(/2))
  1325.             IGMN=MIN(IGAU,VELCHE(/1))
  1326.             XDDL(IE)=VELCHE(IGMN,IBMN)
  1327.             IE=IE+1
  1328.  4085     CONTINUE
  1329. 4785    CONTINUE
  1330. C
  1331. C     on cherche  les coordonnees des noeuds de l element ib
  1332. C
  1333.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1334. C
  1335. C     calcul des exes locaux et des coordonnees locales
  1336. C
  1337.         CALL JOPLOC(XE,SHPTOT,NBBB,NBNO,IFOUR,XEL,BPSS)
  1338. C
  1339.         CALL INTDEL(XNTH,XNTB,XNTT,LPP,MELE)
  1340. C
  1341. C     boucle sur les points de gauss
  1342. C
  1343.         ISDJC=0
  1344. C
  1345.         DO 5185  IGAU=1,NBPTEL
  1346. C
  1347.           CALL BNPQRJ(IGAU,NBNO,NBBB,LRE,IFOUR,LHOOK,NSTN,XE,XEL,
  1348.      &           SHPTOT,SHPWRK,BPSS,BGENE,XGENE,DJAC,IDECAP,NSTB,1)
  1349. C
  1350.           IF (DJAC.EQ.0.D0) THEN
  1351.             INTERR(1)=IB
  1352.             if (noer.eq.0) CALL ERREUR(259)
  1353.             kerr=259
  1354.             GOTO 9985
  1355.           ENDIF
  1356.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  1357. C
  1358.           CALL BST(BGENE,XDDL,LRB,LHOOK,XSTRS)
  1359. C
  1360. C     calcul de la pression
  1361. C
  1362.           IE=LRB
  1363.           DO 4985 IPR=1,IDECAP
  1364.             XP=0.D0
  1365.             IPR1=(IPR-1)*LPP
  1366.             DO 4485  ID=1,LPP
  1367.               IE=IE+1
  1368.               XP=XP+XNTT(ID)*XGENE(IPR,ID+IPR1)*XDDL(IE)
  1369.  4485       CONTINUE
  1370.             XSTRS(NSTRS-IDECAP+IPR)=XP
  1371.  4985     CONTINUE
  1372. C
  1373. C     remplissage du segment contenant les deformations
  1374. C
  1375.           MPTVAL=IVAEPS
  1376.           DO 7185 ICOMP=1,NSTRS
  1377.             MELVAL=IVAL(ICOMP)
  1378.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  1379.  7185     CONTINUE
  1380. C
  1381.  5185   CONTINUE
  1382. C
  1383.         IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  1384.           kerr=195
  1385.           INTERR(1)=IB
  1386.          if (noer.eq.0) then
  1387.           CALL ERREUR(195)
  1388.           GOTO 9985
  1389.          endif
  1390.         ENDIF
  1391. C
  1392.  3185 CONTINUE
  1393. C
  1394.  9985 CONTINUE
  1395. C
  1396.       GOTO 510
  1397. C____________________________________________________________________
  1398.  99   CONTINUE
  1399.       MOTERR(1:4)=NOMTP(MELE)
  1400.       MOTERR(9:12)='EPSI'
  1401.       CALL ERREUR(86)
  1402.  
  1403.  510  CONTINUE
  1404.       SEGSUP,MWRK1
  1405.       IF (MWRK2.NE.0)  SEGSUP,MWRK2
  1406.       IF (MWRK3.NE.0)  SEGSUP,MWRK3
  1407.       IF (MWRK5.NE.0)  SEGSUP,MWRK5
  1408.       IF (MTRACE.NE.0) SEGSUP MTRACE
  1409.  
  1410. c      RETURN
  1411.       END
  1412.  
  1413.  
  1414.  
  1415.  
  1416.  

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