Sources Esope | Cast3M

mocu2d
C MOCU2D    SOURCE    CB215821  24/04/12    21:16:44     11897          
      SUBROUTINE MOCU2D(IPCURZ,IPAXIA,IPMODL,IPCAR,IPEPAX,IPMOMZ
     . ,TOL,mtable,zveri)
**********************************************************************
*
*     OPERATEUR MOCU2D  (MOment/CUrvature)
*
* IPCURZ :  Histoire en courbure
* IPAXIA :   Histoire en effort normal
* IPMODL : Modele
* IPCAR :    Materiau
* IPEPAX :  Histoire allongement axial
* 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(3),SIGMA(3),CRIGI(12),CMASS(12)
      logical bool,zveri
      character*1 char1
 
C
-INC PPARAM
-INC CCOPTIO
-INC SMLREEL
-INC CCGEOME
-INC CCREEL
-INC SMMODEL
-INC SMCHAML
-INC SMTABLE
      POINTEUR IPCURZ.MLREEL,IPAXIA.MLREEL
      POINTEUR IPMOMZ.MLREEL,IPEPAX.MLREEL
C
      NMAXIT = 50
      NSTRS2 = 3
      II     = 0
      XVA    = 0.D0
      BOOL   =.FALSE.
      IOB    = 0
C
*
      ICOUL1=IDCOUL
C
C------- NSTEPS --
C
      SEGACT,IPCURZ,IPAXIA
      NSTEPS=IPCURZ.PROG(/1)
      IF(IPAXIA.PROG(/1).NE.NSTEPS)THEN
        SEGDES,IPCURZ,IPAXIA
        RETURN
      ENDIF
C
      if (zveri) then
         m=3
         segini mtable
         m=nsteps + 1
         segini mtab1
         call ecctab (mtable,'MOT     ',ii,xva,'CONTRAINTES',bool,iob,
     $         'TABLE   ',ii,xva,char1,bool,mtab1)
         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,3
        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
C
      JG=NSTEPS
      SEGINI,IPMOMZ,IPEPAX
C
C     INITIAL STIFNESS
C
      CALL FRIGIE(IPMODL,IPCAR,CRIGI,CMASS)
      XKFXEX=CRIGI(1)
      XKFXCZ=-CRIGI(2)
C
C     NULL DEFORMATION INCREMENT
C
      DEPSI(2)=0.D0
C
      ZAUX0=0.D0
      ZAXI0=0.D0
      DEFOA=0.D0
C
C-------- LOOP ---------------
C
      DO 2000 I=1,NSTEPS
C
        ZAUX1=IPCURZ.PROG(I)
        ZAXI1=IPAXIA.PROG(I)
        DZCURV=ZAUX1-ZAUX0
        DFAXIA=ZAXI1-ZAXI0
C
        SIGPRE=ZAXI0
        XKFXEX=CRIGI(1)
        DSITER=(DFAXIA-XKFXCZ*DZCURV)/XKFXEX
        SZERO =DSITER
C
        DEPSI(3)=DZCURV
C
        DO 5000 NITER=1,NMAXIT
          DEPSI(1)=SZERO
C
**        write(6,*) 'MOCU2D definition incorrecte de time0 et timef'
          time0=0.D0
          timef=0.D0
          CALL FCOUL1(DEPSI,IPMODL,IPSIG0,IPVAR0,IPCAR,time0,timef,
     $                SIGMA,IPSIGF,IPVARF,IRETO,NSTRS2)
          IF(IRETO.EQ.0) THEN
            SEGSUP,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,IPMOMZ,IPEPAX
        SEGDES,IPCURZ,IPAXIA
        RETURN
C
5010    ZAXI0=SIGMA(1)
        ZAUX0=ZAUX1
        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',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  ',I,xva,char1,bool,iob,
     $                  'FLOTTANT  ',ii,DEFOA,char1,bool,iob)
       endif
C
        IPMOMZ.PROG(I)=SIGMA(3)
        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,IPCURZ,IPAXIA
C
      SEGDES,IPMOMZ,IPEPAX
C
      RETURN
      END