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  24/10/21    21:15:10     12042          
  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.       SEGMENT MPTVAL
  63.         INTEGER IPOS(NS),NSOF(NS)
  64.         INTEGER IVAL(NCOSOU)
  65.         CHARACTER*16 TYVAL(NCOSOU)
  66.       ENDSEGMENT
  67.  
  68.       SEGMENT MWRK1
  69.         REAL*8 DDHOOK(NSTRS,NSTRS),XDDL(LRE),XSTRS(NSTRS)
  70.         REAL*8 XE(3,NBBB),DDHOMU(NSTRS,NSTRS)
  71.         REAL*8 SHPWRK(6,NBNO),BGENE(LHOOK,LRE)
  72.         REAL*8 XE1(3,NBBB),XE2(3,NBBB),xstrs2(NSTRS)
  73.         REAL*8 xjac(3,3),valmat(20)
  74.       ENDSEGMENT
  75.  
  76.       SEGMENT MWRK2
  77.         REAL*8 TENS(9),tentra(9),xddls2(lre)
  78.         REAL*8 BGR(NGRA,LRE),BB(2,NGRA),gradi(ngra),R(ngra),u(ngra)
  79.       ENDSEGMENT
  80.  
  81.       SEGMENT MWRK3
  82.         REAL*8 BPSS(3,3),XEL(3,NBBB)
  83.         REAL*8 XNTH(LPP,LPP),XNTB(LPP,LPP),XNTT(LPP)
  84.       ENDSEGMENT
  85.  
  86.       SEGMENT MWRK5
  87.         REAL*8 XGENE(NSTN,LRN)
  88.       ENDSEGMENT
  89.  
  90.       SEGMENT MTRACE
  91.         REAL*8 TRACE(NBPTEL)
  92.       ENDSEGMENT
  93.  
  94.       CHARACTER*8 CMATE
  95.  
  96.       DIMENSION A(4,60),BBX(3,60),UDPGE(3)
  97.       DIMENSION IN(6),JN(6),ITAB(3,3),PP(4,4)
  98.       real*8 valcar(12),var(3)
  99.       real*8 cobma(lhook)
  100.  
  101.       DATA IN/1,2,3,1,1,2/
  102.       DATA JN/1,2,3,2,3,3/
  103.  
  104.       DATA ITAB(1,1),ITAB(1,2),ITAB(1,3)/1,4,5/
  105.       DATA ITAB(2,1),ITAB(2,2),ITAB(2,3)/4,2,6/
  106.       DATA ITAB(3,1),ITAB(3,2),ITAB(3,3)/5,6,3/
  107.  
  108.       real*8 s(2)
  109.  
  110.       s(1)=0.d0
  111.       s(2)=0.d0
  112.       kerr=0
  113.  
  114.       MWRK1  = 0
  115.       MWRK2  = 0
  116.       MWRK3  = 0
  117.       MWRK5  = 0
  118.       MTRACE = 0
  119.  
  120. C  Introduction du point autour duquel se fait le mouvement
  121. C  de la section en defo plane generalisee
  122. C  IIPDPG = numero du noeud/point support si defini pour le modele
  123. C  NDPGE > 0 si prise en compte du point support
  124.       IF (IIPDPG.GT.0) THEN
  125.         IF (IFOUR.EQ.-3) THEN
  126.           NDPGE=3
  127.           UDPGE(1)=UZDPG
  128.           UDPGE(2)=RYDPG
  129.           UDPGE(3)=RXDPG
  130.           IREF=(IIPDPG-1)*(IDIM+1)
  131.           XDPGE=XCOOR(IREF+1)
  132.           YDPGE=XCOOR(IREF+2)
  133.         ELSE IF (IFOUR.EQ.11) THEN
  134.           NDPGE=2
  135.           UDPGE(1)=UZDPG
  136.           UDPGE(2)=RXDPG
  137.           UDPGE(3)=XZero
  138.           XDPGE=XZero
  139.           YDPGE=XZero
  140.         ELSE IF (IFOUR.EQ. 7 .OR. IFOUR.EQ. 8 .OR. IFOUR.EQ. 9 .OR.
  141.      &           IFOUR.EQ.10 .OR. IFOUR.EQ.14) THEN
  142.           NDPGE=1
  143.           UDPGE(1)=UZDPG
  144.           UDPGE(2)=XZero
  145.           UDPGE(3)=XZero
  146.           XDPGE=XZero
  147.           YDPGE=XZero
  148.         else
  149.           write(ioimp,*) 'EPSI2 : ERREUR NDPGE'
  150.           call erreur(5)
  151.           return
  152.         ENDIF
  153.       ELSE
  154.         NDPGE=0
  155.         UDPGE(1)=UZDPG
  156.         UDPGE(2)=XZero
  157.         UDPGE(3)=XZero
  158.         XDPGE=XZero
  159.         YDPGE=XZero
  160.       ENDIF
  161.  
  162.       MELEME=IPMAIL
  163.       NBNN  =NUM(/1)
  164.       NBELEM=NUM(/2)
  165.  
  166.       NHRM=NIFOUR
  167.       MINTE=IPMINT
  168.       NBBB=NBNN
  169.  
  170. C Petite verification prealable (normalement inutile)
  171.       mptval = IVAEPS
  172.       if (NSTRS.ne.ival(/1)) then
  173.         write(ioimp,*) 'EPSI3 : incoherence NSTRS & IVAEPS'
  174.         call erreur(5)
  175.         return
  176.       endif
  177.       do icomp = 1, NSTRS
  178.         melval = IVAL(ICOMP)
  179.         if (melval.le.0) then
  180.           write(ioimp,*) 'EPSI3 : incoherence IVAEPS ival(',icomp,')=0'
  181.           call erreur(5)
  182.           return
  183.         endif
  184.         if (NBPTEL.NE.melval.velche(/1)) then
  185.           write(ioimp,*) 'EPSI3 : incoherence NSTRS & IVAEPS'
  186.           call erreur(5)
  187.           return
  188.         endif
  189.         if (NBELEM .NE. melval.velche(/2)) then
  190.           write(ioimp,*) 'EPSI3 : incoherence NSTRS & IVAEPS'
  191.           call erreur(5)
  192.           return
  193.         endif
  194.       enddo
  195.  
  196. C_______________________________________________________________________
  197. C
  198. C     numero des etiquettes      :
  199. C     etiquettes de 1 a 98 pour traitement specifique a l element
  200. C     dans la zone specifique a chaque element commencant par :
  201. C     5  continue
  202. C     element 5   etiquettes 1005 2005 3005 4005   ...
  203. C     44 continue
  204. C     element 44  etiquettes 1044 2044 3044 4044   ...
  205. C_______________________________________________________________________
  206. C
  207.       IF(MELE.GE.1.AND.MELE.LE.100) THEN
  208. C            CABL SEG2 SEG3 TRI3 TRI4 TRI6 TRI7 QUA4 QUA5 QUA8
  209.       GOTO (   99,  99,  99,   4,  99,   4,  99,   4,  99,   4
  210. C            QUA9 RAC2 RAC3 CUB8 CU20 PRI6 PR15 LIA3 LIA4 LIA6
  211.      1      ,  99,  99,  99,   4,   4,   4,   4,  99,  99,  99
  212. C            LIA8 MULT TET4 TE10 PYR5 PY13 COQ3  DKT POUT LISP
  213.      2      ,  99,  99,   4,   4,   4,   4,  99,  99,  99,  99
  214. C            FAC3 FAC4 FAC6 FAC8 LTR3 LQU4 LCU8 LPR6 LTE4 LPY5
  215.      3      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  216. C            COQ8 TUYA TUFI COQ2 POI1 BARR RACO LSU2 COQ4 LISM
  217.      4      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  218. C            COF3 RES2 LSU3 LSU4 LICO COQ6 CVS2 CVS3 CVT3 CVT6
  219.      5      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  220. C            CVQ4 CVQ8 THP5 TH13 THP6 TH15 THC8 TH20 ICT3 ICQ4
  221.      6      ,  99,  99,  99,  99,  99,  99,  99,  99,   4,   4
  222. C            ICT6 ICQ8 ICC8 ICT4 ICP6 IC20 IC10 IC15 TRIP QUAP
  223.      7      ,   4,   4,   4,   4,   4,   4,   4,   4,  79,  79
  224. C            CUBP TETP PRIP TIMO JOI2 JOI3 JOT3 JOI4 JOI6 JOI8
  225.      8      ,  79,  79,  79,  99,  99,  99,  99,  99,  99,  99
  226. C            LISC TRIH  DST LIC4 CERC TUYO LSE2 LITU HYT3 HYQ4
  227.      9      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99)
  228. c cccccc
  229.      .      ,MELE
  230.       ELSEIF(MELE.GE.101.AND.MELE.LE.200) THEN
  231. C            HYT4 HYP6 HYC8 TRIS QUAS POIS FOR3 JOP3 JOP6 JOP8
  232.       GOTO (   99,  99,  99,  99,  99,  99,  99,  80,  80,  80
  233. C            POL3 POL4 POL5 POL6 POL7 POL8 POL9 PO10 PO11 PO12
  234.      1      ,   4,   4,   4,   4,   4,   4,   4,   4,   4,   4
  235. C            PO13 PO14 BAR3 BAEX LIA2 QUAH CUBH ROT3 SEF2 TRF3
  236.      2      ,   4,   4,  99,  99,  99,  99,  99,  99,  99,  99
  237. C            QUF4 CUF8 PRF6 TEF4 PYF5 MSE3 MTR6 MQU9 MC27 MP18
  238.      3      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  239. C            MT10 MP14 SEF3 TRF7 QUF9 CF27 PF21 TF15 PF19 SEG6
  240.      4      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  241. C            TR21 QU36 C216 P126 TE56 PY91 TRH6 BSE2 BTR4 BQU5
  242.      5      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  243. C            BCU9 BPR7 BTE5 BPY6 FRO4 SEGS POJS JCT3 JCI4 JGI2
  244.      6      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  245. C            JGT3 JGI4 TRIQ QUAQ CUBQ TETQ PRIQ TRIR QUAR CUBR
  246.      7      ,  99,  99, 173, 173, 173, 173, 173, 173, 173, 173
  247. C            TETR PRIR Q4RI Q8RI JOQ3 JOQ6 JOQ8 JOR3 JOR6 JOR8
  248.      8      , 173, 173,   4,   4, 185, 185, 185, 185, 185, 185
  249. C            T1D2 T1D3 M1D2 M1D3 LC03 LC07 LC09 LC27 LC21 LC15
  250.      9      ,  99,  99,   4,   4,  99,  99,  99,  99,  99,  99)
  251. c cccccc
  252.      .      ,MELE-100
  253.       ELSEIF(MELE.GE.201.AND.MELE.LE.300) THEN
  254. C            LC19 LS03 LS07 LS09 LS27 LS21 LS15 LS19 BS03 BS07
  255.       GOTO (   99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  256. C            BS09 BS27 BS21 BS15 BS19 MC03 MC07 MC09 MC27 MC21
  257.      1      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  258. C            MC15 MC19 M103 M107 M109 M127 M121 M115 M119 MS03
  259.      2      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  260. C            MS07 MS09 MS27 MS21 MS15 MS19 QC03 QC07 QC09 QC27
  261.      3      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  262. C            QC21 QC15 QC19 Q103 Q107 Q109 Q127 Q121 Q115 Q119
  263.      4      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  264. C            QS03 QS07 QS09 QS27 QS21 QS15 QS19 CIFL SURE SHB8
  265.      5      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  266. C            CAF2 CAF3 XQ4R XC8R JOI1 ZCO2 ZCO3 ZCO4 TUY2 TUY3
  267.      6      ,  99,  99,  99,  99,  99,  99,  99,  99,  99,  99
  268. C            COS2 COA2 ICY5 IC13 CU27 PR21 TE15 PY19 C20R P15R
  269.      7      ,  99,  99,   4,   4,   4,   4,   4,   4,   4,   4)
  270. c cccccc
  271.      .      ,MELE-200
  272.       ENDIF
  273.       GOTO 99
  274. C
  275. C_______________________________________________________________________
  276. C
  277. C  elements massifs et elements incompressibles MECANIQUE
  278. C_______________________________________________________________________
  279. C
  280.   4   CONTINUE
  281.       IF (MFR.EQ.71 .OR. MFR.EQ.73) GOTO 97173
  282.  
  283. C IDERI <= 2 pour lineaire et quadratique et = 5 pour utilisateur
  284. C ===============================================================
  285.         IF ( IDERI.LE.2.OR.IDERI.EQ.5 ) THEN
  286.  
  287. C Elements massifs en FORMULATION 'MECANIQUE'
  288. C -------------------------------------------
  289.       NBNO=NBNN
  290.       NDDD=NDEP-NDPGE
  291. C
  292. C     Donnees liees a l'element de reference
  293. C
  294.       SEGINI,MWRK1
  295.       IF (Ideri.eq.2) SEGINI,MTRACE
  296. C
  297. C  boucle sur les elements
  298. C
  299.       DO 3004 IB=1,NBELEM
  300. C
  301. C     on cherche les deplacements
  302. C
  303.         MPTVAL=IVADEP
  304.         IE=1
  305.         DO IGAU=1,NBNN
  306.           DO ICOMP=1,NDDD
  307.             MELVAL=IVAL(ICOMP)
  308.             IGMN=MIN(IGAU,VELCHE(/1))
  309.             IBMN=MIN(IB  ,VELCHE(/2))
  310.             XDDL(IE)=VELCHE(IGMN,IBMN)
  311.             IE=IE+1
  312.           ENDDO
  313.         ENDDO
  314.         IF (NDPGE.GT.0) THEN
  315.           DO i=1,NDPGE
  316.             XDDL(IE)=UDPGE(i)
  317.             IE=IE+1
  318.           ENDDO
  319.         ENDIF
  320. C
  321. C     on cherche les coordonnees des noeuds de l element ib
  322. C
  323.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  324.  
  325. C on se met a mi-pas
  326.         if (ideri.eq.5) then
  327.           do iyu=1,nbnn
  328.             i_z = (iyu-1)*nddd
  329.             do i=1,idim
  330.               XE(i,iyu)= xe(i,iyu) + xddl( i + i_z )*0.5D0
  331.             enddo
  332.           enddo
  333.         endif
  334. C
  335. C     boucle sur les points de gauss
  336. C
  337.         ISDJC=0
  338. C
  339. C     Calcul des coeff de modification de b-barre (elts incompres)
  340. C= NOM  :   ICT3, ICQ4, ICT6, ICQ8, ICC8, ICT4, ICP6, IC20, IC10, IC15
  341. C= MELE :    69 ,  70 ,  71 ,  72 ,  73 ,  74 ,  75 ,  76 ,  77 ,  78
  342.         IF (MFR.EQ.31) THEN
  343.           CALL BBCALC(MELE,NBNN,MFR,IDIM,XE,
  344.      1                NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  345.      2                NSTRS,LRE,IFOUR,NHRM,A,BBX,SHPTOT,SHPWRK,
  346.      3                BGENE,XDPGE,YDPGE,PP)
  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)
  707.         ENDIF
  708. C
  709.         DO 594  IGAU=1,NBPTEL
  710.           CALL BMATST(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  711.      1                MELE,MFR,NBNN,LRE,IFOUR,NSTRS,NHRM,1.D0,XE,
  712.      2                SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE)
  713.  
  714.           IF (DJAC.EQ.0.D0) THEN
  715.             kerr=259
  716.             if (noer.eq.0) THEN
  717.               INTERR(1)=IB
  718.               CALL ERREUR(259)
  719.             endif
  720.             GOTO 9964
  721.           ENDIF
  722.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  723.  
  724. C En cas d'elements incompressibles : BGENE selon la methode B-BARRE
  725.           IF (MFR.EQ.31) THEN
  726.             CALL BBAR(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  727.      &                MELE,NBNN,LRE,IFOUR,NSTRS,XE,DJAC,A,BBX,BGENE)
  728.           ENDIF
  729. C
  730.           CALL BST(BGENE,XDDL,LRE,NSTRS,XSTRS)
  731. C dans xstrs on a les deformations II sur config mi pas
  732. C on va calculer grad u/2 puis decomposition polaire puis rtens
  733. C on retravaille sur config initiale
  734.          r_z=XZero
  735.          iipdpg=0
  736.          CALL BGRMAS(iGau,NOELE,NBNO,LRE,IFOUR,NGRA,NIFOUR,XE1,
  737.      .                r_z,SHPTOT,SHPWRK,BB,BGR,DJAC,IIPDPG)
  738.          do iou=1,lre
  739.            xddls2(iou)= 0.5D0 * xddl(iou)
  740.          enddo
  741.          CALL BGRDEP(BGR,NGRA,XDDLs2,LRE,GRADI)
  742. C on ajoute l'identite au gradient
  743.          if(idim.eQ.2) then
  744.            gradi(1)=gradi(1)+1.D0
  745.            gradi(4)=gradi(4)+1.D0
  746.          ELSE IF(IDIM.EQ.3) THEN
  747.            gradi(1)=gradi(1)+1.D0
  748.            gradi(5)=gradi(5)+1.D0
  749.            gradi(9)=gradi(9)+1.D0
  750.          ENDIF
  751.  
  752.          CALL POLA2(gradi,R,U,IDIM)
  753. C fait le rtens Rt.A.R on utilise u pour mettre Rt
  754. C et on met le tenseur dans le tableau tens
  755. C attention vu le stockage R est en fait Rt
  756.          if(idim.eq.2) then
  757.            U(1)=r(1)
  758.            u(2)=r(3)
  759.            U(3)=R(2)
  760.            u(4)=R(4)
  761.            tens(1)=xstrs(1)
  762.            tens(2)=xstrs(4)*0.5d0
  763.            tens(3)=xstrs(4)*0.5d0
  764.            tens(4)=xstrs(2)
  765.  
  766.          elseif(idim.eq.3) then
  767.            U(1)=r(1)
  768.            u(2)=r(4)
  769.            U(3)=R(7)
  770.            u(4)=R(2)
  771.            u(5)=r(5)
  772.            u(6)=r(8)
  773.            u(7)=r(3)
  774.            u(8)=r(6)
  775.            u(9)=r(9)
  776.            tens(1)=xstrs(1)
  777.            tens(2)=xstrs(4)*0.5D0
  778.            tens(3)=xstrs(5)*0.5D0
  779.            tens(4)=tens(2)
  780.            tens(5)=xstrs(2)
  781.            tens(6)=xstrs(6)*0.5D0
  782.            tens(7)=tens(3)
  783.            tens(8)=tens(6)
  784.            tens(9)=xstrs(3)
  785.          else
  786.            write(6,*)' idim est ni 2 ni 3 stop'
  787.            stop
  788.          endif
  789.  
  790.          CALL MULMAT(tentra,tens,U,IDIM,IDIM,IDIM)
  791.          CALL MULMAT(tens,R,Tentra,IDIM,IDIM,IDIM)
  792. C  tens contient le nouveau tenseur  on  va remplir xstrs
  793. C en 2 D epzz ne change pas
  794.          if(idim.eq.2) then
  795.            xstrs(1)=tens(1)
  796.            xstrs(2)=tens(4)
  797.            xstrs(4)=tens(2)*2.D0
  798.          else
  799.            xstrs(1)=tens(1)
  800.            xstrs(2)=tens(5)
  801.            xstrs(3)=tens(9)
  802.            xstrs(4)=tens(2)*2.D0
  803.            xstrs(5)=tens(3)*2.D0
  804.            xstrs(6)=tens(6)*2.D0
  805.         endif
  806. C
  807. C     remplissage du segment contenant les deformations
  808. C
  809.           MPTVAL=IVAEPS
  810.           DO ICOMP=1,NSTRS
  811.             MELVAL=IVAL(ICOMP)
  812.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  813.          ENDDO
  814. C
  815.  594   CONTINUE
  816. C
  817. C     fin de la boucle sur les points de gauss
  818. C
  819.         IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  820.           INTERR(1)=IB
  821.          if (noer.eq.1) then
  822.           kerr=195
  823.          else
  824.           CALL ERREUR(195)
  825.           GOTO 9964
  826.          endif
  827.         ENDIF
  828.  
  829.  394   CONTINUE
  830. C
  831. C    fin de la boucle sur les elements
  832. C
  833.  9964  CONTINUE
  834.       endif
  835. C
  836.       GOTO 510
  837.  
  838. C Elements massifs en FORMULATIONs 'ELECTROSTATIQUE' et 'DIFFUSION'
  839. C -----------------------------------------------------------------
  840. 97173 CONTINUE
  841.       NBNO = NBNN
  842.       NDDD = NDEP
  843.       SEGINI,MWRK1
  844. C-------------------------
  845. C Boucle sur les elements
  846. C-------------------------
  847.       DO IEL = 1, NBELEM
  848. C - Recuperation des coordonnees des noeuds de l element IEL
  849.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE)
  850. C - Recuperation des inconnues primales aux noeuds de l element IEL
  851.         MPTVAL = IVADEP
  852.         IE = 1
  853.         DO IGAU = 1, NBNN
  854.           DO ICOMP = 1, NDDD
  855.             MELVAL = IVAL(ICOMP)
  856.             IGMN = MIN(IGAU,VELCHE(/1))
  857.             IEMN = MIN(IEL ,VELCHE(/2))
  858.             XDDL(IE) = VELCHE(IGMN,IEMN)
  859.             IE = IE+1
  860.           ENDDO
  861.         ENDDO
  862. C--  --  --  --  --  --  --  --  --
  863. C - Boucle sur les points de Gauss
  864. C--  --  --  --  --  --  --  --  --
  865.         ISDJC=0
  866.         DO IGAU = 1, NBPTEL
  867. C -- Calcul de la matrice B et du jacobien au point de Gauss IGAU
  868.           IF (MFR.EQ.71) THEN
  869.             CALL BELEC(XE,SHPTOT(1,1,IGAU),NBNN,LHOOK,-1,
  870.      &                 SHPWRK,BGENE,DJAC)
  871.           ELSE IF (MFR.EQ.73) THEN
  872.             CALL BDIFF(XE,SHPTOT(1,1,IGAU),NBNN,LHOOK,-1,
  873.      &                 SHPWRK,BGENE,DJAC)
  874.           ENDIF
  875.           IF (DJAC.EQ.0.D0) THEN
  876.             kerr=259
  877.             if (noer.eq.0) THEN
  878.               INTERR(1)=IEL
  879.               CALL ERREUR(259)
  880.             endif
  881.             GOTO 98173
  882.           ENDIF
  883.           IF (DJAC.LT.0.D0) ISDJC = ISDJC+1
  884.           CALL BST(BGENE,XDDL,LRE,NSTRS,XSTRS)
  885. C -- Remplissage du segment contenant les "deformations" = -grad(.)
  886.           MPTVAL = IVAEPS
  887.           DO ICOMP = 1, NSTRS
  888.             MELVAL = IVAL(ICOMP)
  889.             VELCHE(IGAU,IEL) = XSTRS(ICOMP)
  890.           ENDDO
  891. C--  --  --  --  --  --  --  --  --
  892.         ENDDO
  893. C--  --  --  --  --  --  --  --  --
  894.         IF (ISDJC.NE.0 .AND. ISDJC.NE.NBPGAU) THEN
  895.           kerr=195
  896.           if (noer.eq.0) THEN
  897.             INTERR(1)=IEL
  898.             CALL ERREUR(195)
  899.             GOTO 98173
  900.           endif
  901.         ENDIF
  902. C-------------------------
  903.       ENDDO
  904. C-------------------------
  905. 98173 CONTINUE
  906.       GOTO 510
  907.  
  908. C_______________________________________________________________________
  909. C
  910. C     milieux poreux
  911. C_______________________________________________________________________
  912. C
  913.   79  CONTINUE
  914. C
  915. C    pour ces elements  nbbb = nombre de noeuds
  916. C                       nbno = nombre de fonctions de forme
  917. C
  918.       NBNO=IPORE
  919.       NSTN=1
  920.       LRN=NBNO-NBBB
  921.       LRB=LRE-LRN
  922. C
  923.       SEGINI,MWRK1,MWRK5
  924. C  Initialisation de MTRACE necessaire mais inutilise pour ces elements
  925.       IF (IREPS2.EQ.1) SEGINI MTRACE
  926. C
  927.       DO 3079  IB=1,NBELEM
  928. C
  929. C     on cherche  les coordonnees des noeuds de l element ib
  930. C
  931.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  932. C
  933. C     on cherche les deplacements
  934. C
  935.         MPTVAL=IVADEP
  936.         IE=1
  937.         DO IGAU=1,NBNN
  938.           DO ICOMP=1,NDEP-1
  939.             MELVAL=IVAL(ICOMP)
  940.             IGMN=MIN(IGAU,VELCHE(/1))
  941.             IBMN=MIN(IB  ,VELCHE(/2))
  942.             XDDL(IE)=VELCHE(IGMN,IBMN)
  943.             IE=IE+1
  944.           ENDDO
  945.         ENDDO
  946. C
  947. C     puis  les pressions
  948. C
  949.         MELVAL=IVAL(NDEP)
  950.         IBMN=MIN(IB,VELCHE(/2))
  951.         DO IGAU=1,LRN
  952.           IGAUSO=IBSOM(NSPOS(IELE)+IGAU-1)
  953.           IGMN=MIN(IGAUSO,VELCHE(/1))
  954.           XDDL(IE)=VELCHE(IGMN,IBMN)
  955.           IE=IE+1
  956.         ENDDO
  957. C
  958. C     boucle sur les points de gauss
  959. C
  960.         ISDJC=0
  961. C
  962.         DO 5079 IGAU=1,NBPTEL
  963. C
  964.           CALL BNPORE(IGAU,NBNO,NBBB,LRE,IFOUR,LHOOK,NSTN,NHRM,
  965.      &                1.D0,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,1)
  966. C
  967.           IF (DJAC.EQ.0.D0) THEN
  968.             INTERR(1)=IB
  969.             if (noer.eq.0) CALL ERREUR(259)
  970.             kerr=259
  971.             GOTO 9979
  972.           ENDIF
  973.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  974. C
  975.           CALL BST(BGENE,XDDL,LRE,LHOOK,XSTRS)
  976. C
  977. C     calcul des eps 2
  978. C
  979.           IF (IREPS2.EQ.1)
  980.      &      CALL BST2(SHPWRK,XDDL,XE,NBNN,IFOUR,XSTRS,TRACE,
  981.      &                IGAU,XDPGE,YDPGE,UDPGE,NHRM)
  982. C
  983. C     calcul de la pression
  984. C
  985.           XP=0.D0
  986.           DO ID=1,LRN
  987.             XP=XP+XGENE(1,ID)*XDDL(LRB+ID)
  988.           ENDDO
  989.           XSTRS(NSTRS)=XP
  990. C
  991. C     remplissage du segment contenant les deformations
  992. C
  993.           MPTVAL=IVAEPS
  994.           DO 7079 ICOMP=1,NSTRS
  995.             MELVAL=IVAL(ICOMP)
  996.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  997.  7079     CONTINUE
  998. C
  999.  5079   CONTINUE
  1000. C
  1001. C     fin de la boucle sur les points de gauss
  1002. C
  1003.         IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  1004.           INTERR(1)=IB
  1005.          if (noer.eq.1) then
  1006.           kerr=195
  1007.          else
  1008.           CALL ERREUR(195)
  1009.           GOTO 9979
  1010.          endif
  1011.         ENDIF
  1012. C
  1013.  3079 CONTINUE
  1014. C
  1015. C    fin de la boucle sur les elements
  1016. C
  1017.  9979 CONTINUE
  1018. C
  1019.       GOTO 510
  1020. C_______________________________________________________________________
  1021. C
  1022. C     milieux poreux - SUITE
  1023. C_______________________________________________________________________
  1024. C
  1025.  173  CONTINUE
  1026. C
  1027. C    pour ces elements  nbbb = nombre de noeuds
  1028. C                       nbno = nombre de fonctions de forme
  1029. C
  1030.       IF (MELE.GE.173.AND.MELE.LE.177) THEN
  1031.         IDECAP = 2
  1032.       ELSE IF (MELE.GE.178.AND.MELE.LE.182) THEN
  1033.         IDECAP = 3
  1034.       ENDIF
  1035. C
  1036.       NBNO=IPORE
  1037.       NSTN=IDECAP
  1038.       NSTB=4
  1039.       IF(IFOUR.EQ.1.OR.IFOUR.EQ.2) NSTB=6
  1040. C
  1041.       LPP=NBNO-NBBB
  1042.       LRN=IDECAP*LPP
  1043.       LRB=LRE-LRN
  1044. C
  1045.       SEGINI,MWRK1,MWRK5
  1046. C  Initialise de MTRACE necessaire mais inutilise pour cet element
  1047.       IF (IREPS2.EQ.1) SEGINI MTRACE
  1048. C
  1049.       DO 3173  IB=1,NBELEM
  1050. C
  1051. C     on cherche  les coordonnees des noeuds de l element ib
  1052. C
  1053.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1054. C
  1055. C     on cherche les deplacements
  1056. C
  1057.         MPTVAL=IVADEP
  1058.         IE=1
  1059.         DO IGAU=1,NBNN
  1060.           DO ICOMP=1,NDEP-IDECAP
  1061.             MELVAL=IVAL(ICOMP)
  1062.             IGMN=MIN(IGAU,VELCHE(/1))
  1063.             IBMN=MIN(IB  ,VELCHE(/2))
  1064.             XDDL(IE)=VELCHE(IGMN,IBMN)
  1065.             IE=IE+1
  1066.           ENDDO
  1067.         ENDDO
  1068. C
  1069. C     puis  les pressions
  1070. C
  1071.         DO IPR = 1,IDECAP
  1072.           MELVAL=IVAL(NDEP-IDECAP+IPR)
  1073.           DO IGAU=1,LPP
  1074.             IGAUSO=IBSOM(NSPOS(IELE)+IGAU-1)
  1075.             IGMN=MIN(IGAUSO,VELCHE(/1))
  1076.             IBMN=MIN(IB  ,VELCHE(/2))
  1077.             XDDL(IE)=VELCHE(IGMN,IBMN)
  1078.             IE=IE+1
  1079.           ENDDO
  1080.         ENDDO
  1081. C
  1082. C     boucle sur les points de gauss
  1083. C
  1084.         ISDJC=0
  1085. C
  1086.         DO 5173  IGAU=1,NBPTEL
  1087. C
  1088.           CALL BNQORE(IGAU,NBNO,NBBB,LRE,IFOUR,NSTB,NSTN,NHRM,
  1089.      &        1.D0,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,IDECAP,LHOOK,1)
  1090. C
  1091.           IF (DJAC.EQ.0.D0) THEN
  1092.             INTERR(1)=IB
  1093.             if (noer.eq.0) CALL ERREUR(259)
  1094.             kerr=259
  1095.             GOTO 9173
  1096.           ENDIF
  1097.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  1098. C
  1099.           CALL BST(BGENE,XDDL,LRE,LHOOK,XSTRS)
  1100. C
  1101. C     calcul des eps 2
  1102. C
  1103.           IF (IREPS2.EQ.1)
  1104.      &      CALL BST2(SHPWRK,XDDL,XE,NBNN,IFOUR,XSTRS,TRACE,
  1105.      &                IGAU,XDPGE,YDPGE,UDPGE,NHRM)
  1106. C
  1107. C     calcul des  pressions
  1108. C
  1109.           IE=LRB
  1110.           DO IPR=1,IDECAP
  1111.             XP=0.D0
  1112.             IPR1=(IPR-1)*LPP
  1113.             DO ID=1,LPP
  1114.               IE=IE+1
  1115.               XP=XP+XGENE(IPR,ID+IPR1)*XDDL(IE)
  1116.             ENDDO
  1117.             XSTRS(NSTRS-IDECAP+IPR)=XP
  1118.           ENDDO
  1119. C
  1120. C     remplissage du segment contenant les deformations
  1121. C
  1122.           MPTVAL=IVAEPS
  1123.           DO 7173 ICOMP=1,NSTRS
  1124.             MELVAL=IVAL(ICOMP)
  1125.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  1126.  7173     CONTINUE
  1127. C
  1128.  5173    CONTINUE
  1129. C
  1130. C     fin de la boucle sur les points de gauss
  1131. C
  1132.          IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  1133.            INTERR(1)=IB
  1134.          if (noer.eq.1) then
  1135.           kerr=195
  1136.          else
  1137.            CALL ERREUR(195)
  1138.            GOTO 9173
  1139.          endif
  1140.          ENDIF
  1141. C
  1142.  3173 CONTINUE
  1143. C
  1144. C    fin de la boucle sur les elements
  1145. C
  1146.  9173 CONTINUE
  1147. C
  1148.       GOTO 510
  1149.  
  1150. C_______________________________________________________________________
  1151. C
  1152. C     joints poreux
  1153. C_______________________________________________________________________
  1154. C
  1155.   80  CONTINUE
  1156. C
  1157. C    pour ces elements  nbbb = nombre de noeuds
  1158. C                       nbno = nombre de fonctions de forme
  1159. C
  1160.       NBNO=IPORE
  1161.       NSTN=1
  1162.       LRN=(NBNO-NBBB)*3/2
  1163.       LPP = LRN
  1164.       LRB=LRE-LRN
  1165.       NFAC=(3*NBBB-NBNO)/2
  1166. C
  1167.       SEGINI,MWRK1,MWRK3,MWRK5
  1168. C
  1169.       DO 3080  IB=1,NBELEM
  1170. C
  1171. C     on cherche d'abord les deplacements
  1172. C
  1173.         MPTVAL=IVADEP
  1174.         IE=1
  1175.         DO 4180 IGAU=1,NFAC
  1176.           DO 4280 ICOMP=1,NDEP-1
  1177.             MELVAL=IVAL(ICOMP)
  1178.             IGMN=MIN(IGAU,VELCHE(/1))
  1179.             IBMN=MIN(IB  ,VELCHE(/2))
  1180.             XDDL(IE)=VELCHE(IGMN,IBMN)
  1181.             IE=IE+1
  1182.  4280     CONTINUE
  1183.  4180   CONTINUE
  1184. C
  1185. C     puis les pressions
  1186. C
  1187.         MELVAL=IVAL(NDEP)
  1188.         DO 4080 IGAU=1,NBNN
  1189.           DO 4190 INSOM=1,NBSOM(IELE)
  1190.             IF (IGAU.EQ.IBSOM(NSPOS(IELE)+INSOM-1)) GO TO 4191
  1191. 4190      CONTINUE
  1192.           IF (IGAU.GT.NFAC) GO TO 4191
  1193.           GO TO 4080
  1194. 4191      CONTINUE
  1195.           IBMN=MIN(IB  ,VELCHE(/2))
  1196.           IGMN=MIN(IGAU,VELCHE(/1))
  1197.           XDDL(IE)=VELCHE(IGMN,IBMN)
  1198.           IE=IE+1
  1199.  4080   CONTINUE
  1200. C
  1201. C     on cherche  les coordonnees des noeuds de l element ib
  1202. C
  1203.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1204. C
  1205. C     calcul des exes locaux et des coordonnees locales
  1206. C
  1207.         CALL JOPLOC(XE,SHPTOT,NBBB,NBNO,IFOUR,XEL,BPSS)
  1208. C
  1209.         CALL INTDEL(XNTH,XNTB,XNTT,LRN,MELE)
  1210. C
  1211. C     boucle sur les points de gauss
  1212. C
  1213.         ISDJC=0
  1214. C
  1215.         DO 5080 IGAU=1,NBPTEL
  1216. C
  1217.           CALL BNPORJ(IGAU,NBNO,NBBB,LRE,IFOUR,LHOOK,NSTN,XE,XEL,
  1218.      &                SHPTOT,SHPWRK,BPSS,BGENE,XGENE,DJAC,1)
  1219. C
  1220.           IF (DJAC.EQ.0.D0) THEN
  1221.             INTERR(1)=IB
  1222.             if (noer.eq.0) CALL ERREUR(259)
  1223.             kerr=259
  1224.             GOTO 9980
  1225.           ENDIF
  1226.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  1227. C
  1228.           CALL BST(BGENE,XDDL,LRB,LHOOK,XSTRS)
  1229.  
  1230. C
  1231. C     calcul de la pression
  1232. C
  1233.           XP=0.D0
  1234.           DO 4480  ID=1,LRN
  1235.             XP=XP+XNTT(ID)*XGENE(1,ID)*XDDL(LRB+ID)
  1236.  4480     CONTINUE
  1237.           XSTRS(NSTRS)=XP
  1238. C
  1239. C     remplissage du segment contenant les deformations
  1240. C
  1241.           MPTVAL=IVAEPS
  1242.           DO 7080 ICOMP=1,NSTRS
  1243.             MELVAL=IVAL(ICOMP)
  1244.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  1245.  7080     CONTINUE
  1246. C
  1247.  5080   CONTINUE
  1248. C
  1249.         IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  1250.           INTERR(1)=IB
  1251.          if (noer.eq.1) then
  1252.           kerr=195
  1253.          else
  1254.           CALL ERREUR(195)
  1255.           GOTO 9980
  1256.          endif
  1257.         ENDIF
  1258. C
  1259.  3080 CONTINUE
  1260. C
  1261.  9980 CONTINUE
  1262.       GOTO 510
  1263.  
  1264. C_______________________________________________________________________
  1265. C
  1266. C     joints poreux  - SUITE
  1267. C_______________________________________________________________________
  1268. C
  1269.  185  CONTINUE
  1270. C
  1271. C    pour ces elements  nbbb = nombre de noeuds
  1272. C                       nbno = nombre de fonctions de forme
  1273. C
  1274.       IF (MELE.GE.185.AND.MELE.LE.187) THEN
  1275.         IDECAP = 2
  1276.       ELSE IF (MELE.GE.188.AND.MELE.LE.190) THEN
  1277.         IDECAP = 3
  1278.       ENDIF
  1279.  
  1280.       NBNO=IPORE
  1281.       NSTN=IDECAP
  1282.       NSTB=2
  1283.       IF(IFOUR.EQ.1.OR.IFOUR.EQ.2) NSTB=3
  1284.  
  1285.       LPP=(NBNO-NBBB)*3/2
  1286.       LRN=IDECAP*LPP
  1287.       LRB=LRE-LRN
  1288.  
  1289.       NFAC=(3*NBBB-NBNO)/2
  1290.  
  1291.       SEGINI,MWRK1,MWRK3,MWRK5
  1292.  
  1293.       DO 3185  IB=1,NBELEM
  1294. C
  1295. C     on cherche d'abord les deplacements
  1296. C
  1297.         MPTVAL=IVADEP
  1298.         IE=1
  1299.         DO 4185 IGAU=1,NFAC
  1300.           DO 4285 ICOMP=1,NDEP-IDECAP
  1301.             MELVAL=IVAL(ICOMP)
  1302.             IGMN=MIN(IGAU,VELCHE(/1))
  1303.             IBMN=MIN(IB  ,VELCHE(/2))
  1304.             XDDL(IE)=VELCHE(IGMN,IBMN)
  1305.             IE=IE+1
  1306.  4285     CONTINUE
  1307.  4185   CONTINUE
  1308. C
  1309. C     puis les pressions
  1310. C
  1311.         DO 4785 IPR=1,IDECAP
  1312.           MELVAL=IVAL(NDEP-IDECAP+IPR)
  1313.           DO 4085 IGAU=1,NBNN
  1314.             DO 4195 INSOM=1,NBSOM(IELE)
  1315.               IF (IGAU.EQ.IBSOM(NSPOS(IELE)+INSOM-1)) GO TO 4891
  1316. 4195        CONTINUE
  1317.             IF (IGAU.GT.NFAC) GO TO 4891
  1318.             GO TO 4085
  1319. 4891        CONTINUE
  1320.             IBMN=MIN(IB  ,VELCHE(/2))
  1321.             IGMN=MIN(IGAU,VELCHE(/1))
  1322.             XDDL(IE)=VELCHE(IGMN,IBMN)
  1323.             IE=IE+1
  1324.  4085     CONTINUE
  1325. 4785    CONTINUE
  1326. C
  1327. C     on cherche  les coordonnees des noeuds de l element ib
  1328. C
  1329.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE)
  1330. C
  1331. C     calcul des exes locaux et des coordonnees locales
  1332. C
  1333.         CALL JOPLOC(XE,SHPTOT,NBBB,NBNO,IFOUR,XEL,BPSS)
  1334. C
  1335.         CALL INTDEL(XNTH,XNTB,XNTT,LPP,MELE)
  1336. C
  1337. C     boucle sur les points de gauss
  1338. C
  1339.         ISDJC=0
  1340. C
  1341.         DO 5185  IGAU=1,NBPTEL
  1342. C
  1343.           CALL BNPQRJ(IGAU,NBNO,NBBB,LRE,IFOUR,LHOOK,NSTN,XE,XEL,
  1344.      &           SHPTOT,SHPWRK,BPSS,BGENE,XGENE,DJAC,IDECAP,NSTB,1)
  1345. C
  1346.           IF (DJAC.EQ.0.D0) THEN
  1347.             INTERR(1)=IB
  1348.             if (noer.eq.0) CALL ERREUR(259)
  1349.             kerr=259
  1350.             GOTO 9985
  1351.           ENDIF
  1352.           IF (DJAC.LT.0.D0) ISDJC=ISDJC+1
  1353. C
  1354.           CALL BST(BGENE,XDDL,LRB,LHOOK,XSTRS)
  1355. C
  1356. C     calcul de la pression
  1357. C
  1358.           IE=LRB
  1359.           DO 4985 IPR=1,IDECAP
  1360.             XP=0.D0
  1361.             IPR1=(IPR-1)*LPP
  1362.             DO 4485  ID=1,LPP
  1363.               IE=IE+1
  1364.               XP=XP+XNTT(ID)*XGENE(IPR,ID+IPR1)*XDDL(IE)
  1365.  4485       CONTINUE
  1366.             XSTRS(NSTRS-IDECAP+IPR)=XP
  1367.  4985     CONTINUE
  1368. C
  1369. C     remplissage du segment contenant les deformations
  1370. C
  1371.           MPTVAL=IVAEPS
  1372.           DO 7185 ICOMP=1,NSTRS
  1373.             MELVAL=IVAL(ICOMP)
  1374.             VELCHE(IGAU,IB)=XSTRS(ICOMP)
  1375.  7185     CONTINUE
  1376. C
  1377.  5185   CONTINUE
  1378. C
  1379.         IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN
  1380.           kerr=195
  1381.           INTERR(1)=IB
  1382.          if (noer.eq.0) then
  1383.           CALL ERREUR(195)
  1384.           GOTO 9985
  1385.          endif
  1386.         ENDIF
  1387. C
  1388.  3185 CONTINUE
  1389. C
  1390.  9985 CONTINUE
  1391. C
  1392.       GOTO 510
  1393. C____________________________________________________________________
  1394.  99   CONTINUE
  1395.       MOTERR(1:4)=NOMTP(MELE)
  1396.       MOTERR(9:12)='EPSI'
  1397.       CALL ERREUR(86)
  1398.  
  1399.  510  CONTINUE
  1400.       SEGSUP,MWRK1
  1401.       IF (MWRK2.NE.0)  SEGSUP,MWRK2
  1402.       IF (MWRK3.NE.0)  SEGSUP,MWRK3
  1403.       IF (MWRK5.NE.0)  SEGSUP,MWRK5
  1404.       IF (MTRACE.NE.0) SEGSUP MTRACE
  1405.  
  1406. c      RETURN
  1407.       END
  1408.  
  1409.  
  1410.  

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