Sources Esope | Cast3M

mocu3d
C MOCU3D    SOURCE    CB215821  24/04/12    21:16:45     11897          
C MOCU3D    SOURCE    DC        98/06/03    21:45:21     3202
C MOCU3D     SOURCE    GOUNAND   97/10/13    22:21:06     2872
      SUBROUTINE MOCU3D(IPCURY,IPCURZ,IPAXIA,IPMODL,IPCAR,
     .   IPEPAX,IPMOMY,IPMOMZ,TOL,mtable,zveri)
**********************************************************************
*
*     OPERATEUR MOCU3D  (MOment/CUrvature)
*
* IPCURY :  Histoire en courbure
* IPCURZ :  Histoire en courbure
* IPAXIA :   Histoire en effort normal
* IPMODL : Modele
* IPCAR :    Materiau
* IPEPAX :  Histoire allongement axial
* IPMOMY : Histoire en moment
* IPMOMZ : Histoire en moment
* TOL : Tolerance
**********************************************************************
* D'APRES MOCUR DC 98
**********************************************************************
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C
      DIMENSION DEPSI(6),SIGMA(6),CRIGI(12),CMASS(12)
      logical bool,boolre,zveri
      character*8 char1,charre
C
 
-INC PPARAM
-INC CCOPTIO
-INC SMLREEL
-INC CCGEOME
-INC CCREEL
-INC SMMODEL
-INC SMCHAML
-INC SMTABLE
      POINTEUR IPCURY.MLREEL,IPCURZ.MLREEL,IPAXIA.MLREEL
      POINTEUR IPMOMY.MLREEL,IPMOMZ.MLREEL,IPEPAX.MLREEL
C
      NMAXIT = 50
      NSTRS2 = 6
      II     = 0
      XVA    = 0.D0
      BOOL   =.FALSE.
      IOB    = 0
      IRE    = 0
      XVARE  = 0.D0
 
C
*
      ICOUL1=IDCOUL
C
C------- NSTEPS --
C
      SEGACT,IPCURY,IPCURZ,IPAXIA
      NSTEPS=IPCURY.PROG(/1)
      IF(IPCURZ.PROG(/1).NE.NSTEPS)THEN
        SEGDES,IPCURY,IPCURZ,IPAXIA
        RETURN
      ENDIF
      IF(IPAXIA.PROG(/1).NE.NSTEPS)THEN
        SEGDES,IPCURY,IPCURZ,IPAXIA
        RETURN
      ENDIF
c
      if (zveri) then
         m=3
         segini mtable
         itab=mtable
         m=nsteps+1
         segini mtab1
         itab1=mtab1
c
         call ecctab (itab,'MOT     ',ii,xva,'CONTRAINTES',bool,iob,
     $         'TABLE   ',ire,xvare,char1,boolre,itab1)
         segini mtab2
         call ecctab (mtable,'MOT     ',ii,xva,'VARIABLES_INTERNES',
     $  bool,iob,
     $         'TABLE   ',ii,xva,char1,bool,mtab2)
         segini mtab3
         call ecctab (mtable,'MOT     ',ii,xva,'DEFORMATIONS',
     $  bool,iob,
     $         'TABLE   ',ii,xva,char1,bool,mtab3)
         segdes mtable
      endif
C
C-------- INITIALIZING ---------------
C
      CALL ZEROP(IPMODL,'CONTRAIN',IPSIG0)
      CALL ZEROP(IPMODL,'VARINTER',IPVAR0)
      DO IE1=1,6
        SIGMA(IE1)=0.D0
      END DO
      SZERO=0.0D0
      if (zveri) then
         call ecctab(mtab1,'ENTIER  ',0,xva,char1,bool,iob,
     $                  'MCHAML  ',ii,xva,char1,bool,ipsig0)
         call ecctab(mtab2,'ENTIER  ',0,xva,char1,bool,iob,
     $                  'MCHAML  ',ii,xva,char1,bool,ipvar0)
         call ecctab(mtab3,'ENTIER  ',0,xva,char1,bool,iob,
     $                  'FLOTTANT',ii,szero,char1,bool,iob)
      endif
      JG=NSTEPS
      SEGINI,IPMOMY,IPMOMZ,IPEPAX
C
C     INITIAL STIFNESS
C
      CALL FRIGIE(IPMODL,IPCAR,CRIGI,CMASS)
      XKFXEX=CRIGI(1)
      XKFXCY=CRIGI(3)
      XKFXCZ=-CRIGI(2)
C
C     NULL DEFORMATION INCREMENT
C
      DEPSI(2)=0.D0
      DEPSI(3)=0.D0
      DEPSI(4)=0.D0
C
      YAUX0=0.D0
      ZAUX0=0.D0
      ZAXI0=0.D0
      DEFOA=0.D0
C
C-------- LOOP ---------------
C
      DO 2000 I=1,NSTEPS
C
        YAUX1=IPCURY.PROG(I)
        ZAUX1=IPCURZ.PROG(I)
        ZAXI1=IPAXIA.PROG(I)
 
        DYCURV=YAUX1-YAUX0
        DZCURV=ZAUX1-ZAUX0
        DFAXIA=ZAXI1-ZAXI0
C
        SIGPRE=ZAXI0
        XKFXEX=CRIGI(1)
        DSITER=(DFAXIA-XKFXCY*DYCURV-XKFXCZ*DZCURV)/XKFXEX
        SZERO =DSITER
C
        DEPSI(5)=DYCURV
        DEPSI(6)=DZCURV
C
        DO 5000 NITER=1,NMAXIT
          DEPSI(1)=SZERO
C
          CALL FCOUL1(DEPSI,IPMODL,IPSIG0,IPVAR0,IPCAR,xzero,xzero,
     $                SIGMA,IPSIGF,IPVARF,IRETO,NSTRS2)
          IF(IRETO.EQ.0) THEN
            SEGSUP,IPMOMY,IPMOMZ,IPEPAX
            RETURN
          ENDIF
C
          DSITET=(ZAXI1-SIGMA(1))/XKFXEX
C
          IF(IIMPI.EQ.1789)THEN
            write(IOIMP,999)I,NITER,XKFXEX,DSITET
 999        format(1x,'I=',I4,' NITER=',I2,' XKFXEX=',
     $               1pd12.5,' DSITET=',1pd12.5)
          ENDIF
C
          IF(ABS(DSITET).LT.TOL)GOTO 5010
C
C  on tente une approximation quasi-newton apres la prediction ...
C
          IF(NITER.GT.1) XKFXEX=(SIGMA(1)-SIGPRE)/DSITER
          SIGPRE=SIGMA(1)
          DSITER=(ZAXI1-SIGMA(1))/XKFXEX
          SZERO=SZERO+DSITER
C
          CALL DTCHAM(IPSIGF)
          CALL DTCHAM(IPVARF)
C
5000    CONTINUE
C
        INTERR(1)=NMAXIT
        CALL ERREUR(638)
        CALL DTCHAM(IPSIG0)
        CALL DTCHAM(IPVAR0)
        SEGSUP,IPMOMY,IPMOMZ,IPEPAX
        SEGDES,IPCURY,IPCURZ,IPAXIA
        RETURN
C
5010    ZAXI0=SIGMA(1)
        ZAUX0=ZAUX1
        YAUX0=YAUX1
        DEFOA=DEFOA+SZERO
C
        if (.not.zveri) then
           CALL DTCHAM(IPSIG0)
           CALL DTCHAM(IPVAR0)
        endif
        IPSIG0=IPSIGF
        IPVAR0=IPVARF
c
        if (zveri) then
           call ecctab(mtab1,'ENTIER',I,xva,char1,bool,iob,
     $                  'MCHAML',ii,xva,char1,bool,ipsig0)
           call ecctab(mtab2,'ENTIER',I,xva,char1,bool,iob,
     $                  'MCHAML',ii,xva,char1,bool,ipvar0)
           call ecctab(mtab3,'ENTIER  ',I,xva,char1,bool,iob,
     $                  'FLOTTANT  ',ii,DEFOA,char1,bool,iob)
       endif
C
 
 
        IPMOMY.PROG(I)=SIGMA(5)
        IPMOMZ.PROG(I)=SIGMA(6)
        IPEPAX.PROG(I)=DEFOA
C
C
2000  CONTINUE
C
      if (zveri) then
         segdes,mtab1,mtab2,mtab3
      else
         CALL DTCHAM(IPSIG0)
         CALL DTCHAM(IPVAR0)
      endif
      SEGDES,IPCURY,IPCURZ,IPAXIA
C
      SEGDES,IPMOMY,IPMOMZ,IPEPAX
C
      RETURN
      END