Télécharger crer3d.eso

Retour à la liste

crer3d
C CRER3D    SOURCE    FD218221  24/02/07    21:15:08     11834          
       SUBROUTINE CRER3D(NBRINC,NTYPC,NINC,NDIMG,RTP,RTI,DELTA,
     # BETA,COHE,FRAC,GCH,bwPw,bgpg,SEFFR,SEFFI,SEFFO,EPSTR,EPSTI,
     # EPSTO,EPSLR,V33,V33T,V33S,V33ST,EPSTG,EPSPLG,SEFFG,FTHG,
     # NTMAX,FT,ACTIFT,DPFTDS,DGFTDS,LNUMCRT,NBRACT,PLASTT,Log_RTG,
     # NCMAX,FC,ACTIFC,DPFCDS,DGFCDS,LNUMCRC,NBRACC,PLASTC,RFP,RFI,
     # NLMAX,FL,ACTIFL,DPFLDS,DGFLDS,LNUMCRL,NBRACL,PLASTL,Log_RLG,
     # STOTR,STOTG,RESG,PLAST,RTG,RTLG,RFG,RFLG,reduc_resg,
     # PRECISION3D,ERR1,AFFICHE)
 
 
c      evaluation des criteres pour incl3d
c      la liste des citeres actif en traction (LNUMCRT ) et en compression (LNUMCRC)
c      est definie dans la repere generalise LNUMCRT(1) est le critere le plus
c      actif de traction-refermeture de valeur REST=FT(LNUMCRT(1)), 
 
c      il en est de même pour :
c      - les criteres de cisaillement aux interfaces des phase  et dans la matrice
c      - les critere de cisaillement, NUMCR est le numero de la phase plastifiee
c      si cette phase et la matrice elle peut inclure des interfaces dans sont
c      ecoulement
 
c      l ordre de rangement des criteres est le suivant (base generalisee des ecoulements)
c      idir=(iphase-1)*3+(itypt-1)*3+jdir 
c      avec idirg direction generalisee
c      iphase numero de la phanse 0 à ninc
c      ityp position zone d interet (1 inclusion 3,4,5 points d interet sur l interface)
c      jdir un des 3 direction principale des contraintes macroscopiques
 
 
c      (A.Sellier 2022/02/21)
c      maj 2022/02/22 ! pour le fun que des 2 ;)
c      maj 2022/27/02 amelioration par triage des criteres et traitement par groupes
c      maj 2022/03/30 criteres en forces thermodynamiques generalisees
c      2022/06/02 simplification pour version20I
c      2022/06/16 ecoulement par methode du gradient generalise
c      2023/06/26 correction methode gradient generalisee
 
       implicit real*8 (a-h,o-z)
       implicit integer (i-n)
 
c      ********variables externes ***************************************      
       integer NBRINC,NINC,IETAP,NTYPC
       logical AFFICHE
c      fractions volumiques des phases      
       real*8 FRAC(0:NBRINC)
c      indicateur de diffusion produits chimiques dans fissure
       real*8 GCH(0:NBRINC)       
c      contraintes totales au sens poro mecanique     
       real*8 SEFFR(-1:NBRINC,6,0:1)
       real*8 SEFFI(-1:NBRINC,6,0:1)
       real*8 SEFFO(-1:NBRINC,6,0:1)
c      contrainte hydrique
       real*8 SEW(0:NBRINC,0:1) 
c      resistance a la traction
       real*8 RTP(-1:NBRINC) 
c      resistance a la refermeture      
       real*8 RTI(-1:NBRINC) 
c      resistance a la refermeture de la phase
       real*8 RFP(-1:NBRINC)
c      resistance a la refermeture  de l interface      
       real*8 RFI(-1:NBRINC) 
c      proprite pour Drucker Pragger
       real*8 DELTA(0:NBRINC),BETA(0:NBRINC),COHE(0:NBRINC)           
c      deformations plastiques de traction      
       real*8 EPSTR(-1:NBRINC,6,0:1)
       real*8 EPSLR(-1:NBRINC,6,0:1)      
       real*8 EPSTI(-1:NBRINC,6,0:1)      
       real*8 EPSTO(-1:NBRINC,6,0:1)
       real*8 EPSRH(-1:NBRINC,6,0:1)
 
c      erreur
       integer err1       
c      precision
       real*8 precision3d
c      logique critere ACTIFT
       logical PLAST        
c      matrice de passage vers la base du critere le plus actif
       real*8 V33(3,3),V33T(3,3),V33D(3,3),V33DT(3,3)
       real*8 V33S(3,3),V33ST(3,3) 
c      donnes poro-mechaniques (eau, gel, dechrage)      
       real*8 bwPw(0:NBRINC),bgpg(0:NBRINC)
c      residu global pour la convergence
       real*8 RESG  
 
c      variables pour les decollement d interface en contrainte totale
       integer NLMAX
       real*8 FL(NLMAX)
       logical ACTIFL(NLMAX)
       real*8 DPFLDS(NLMAX,NDIMG),DGFLDS(NLMAX,NDIMG)
       integer LNUMCRL(NLMAX),NBRACL
       logical Log_RLG(NLMAX)
       logical PLASTL
       real*8 STOTR(-1:NBRINC,6,0:1)
       real*8 EPSLG(NDIMG),STOTG(NDIMG)
 
 
 
c      Variables internes definies dans l espace generalise
 
 
c      dimension de l espace generalise  
       integer NDIMG,NTMAX,NCMAX  
c      Res vecteur des residus des criteres
       real*8 Rest
c      indicateur de plasticite elementaire      
       logical PLASTT,PLASTC,Log_RTG(NTMAX)
c      contrainte en base generalisees et force thermo associee aux epsa
       real*8 SEFFG(NDIMG),FTHG(NDIMG)
c      deformation plastique de traction generalisee
       real*8 EPSTG(NDIMG),EPSPLG(NDIMG)       
c      derivee du critere par rapport à la contrainte
       real*8 DPFTDS(NTMAX,NDIMG),DPFCDS(NCMAX,NDIMG)  
c      derivee du critere par rapport à la contrainte
       real*8 DGFTDS(NTMAX,NDIMG),DGFCDS(NCMAX,NDIMG)    
c      resistances en base generalisee
       real*8 RTG(NDIMG),RTLG(NDIMG),RFG(NDIMG),RFLG(NDIMG)   
c      image des resistances dans la base des forces thermodynamiques
       real*8 FTRAC(NDIMG)
c      activite des criteres      
       logical ACTIFT(NTMAX),ACTIFC(NCMAX)
c      table du critere de rankine et dp  
       real*8 FT(NTMAX),FC(NCMAX)
c      numeror critre retenu pour les tensions refermeture      
       integer LNUMCRT(NTMAX),LNUMCRC(NCMAX)     
c      nombre de criteres actif, numero d ordre du dernier traitement 
       integer NBRACT
c      nombre de criteres actif cisaillement, numero d ordre du dernier traitement 
       integer NBRACC
c      coeff de viscosite sur l annulation des residus
       real*8 reduc_res       
 
c      pour DGFTDS, 1er indice = oredre, second indice = direction generalisee       
 
c     ********* variables locales **************************************
 
      real*8 x6(6),x3(3),x33(3,3)
      real*8 x06(6),x03(3),x16(6)
      real*8 vnorm,fmaxc,fmaxg,fmaxt,fminc,fming,fmint
      integer iphase,idir,idebut,jdir,idirg
 
      logical affiche_local
 
      real*8 press,y3(3),dilatance,frottement,cohesion,press_lim,taulim
      real*8 taueq,resc_lim,xpetit      
      integer posi,icr,ityp
      integer nbracc_i,ipt_i,icr_i,inc,idebut_inc
      logical cond_i
      real*8 SEFF3(3),FTH3(3),DPFCDS3(3),DGFCDS3(3),FC1,som
      logical ACTIFC1
 
c     il faut plastifier tous les pts d interface pour en 
c     plastifier un si cond_i estvrai
      cond_i=.false.
 
c     option d affichage local      
      affiche_local=affiche 
c      affiche_local=.true.
      if(affiche_local) then
        print*,'On est dans crer3d'
      end if
 
c     evaluation des criteres fin de pas
      posi=1     
 
c     ******************************************************************      
c     recherche des criteres actifs en traction 
c     ******************************************************************
 
c      lors de l evalauation en 0 il ne faut pas effacer 
c      les directions ni la liste des criteres actifs   
 
c      ------1er passage sans connaitre le critere le plus actif--------
 
c      variable de controle
       err1=0
c      petite valeur pour comparaison des criteres
       xpetit=0.d0      
 
c     ------debut de la recherche d activite des criteres--------------- 
 
 
       PLAST=.false.
       PLASTT=.false.
c      initialisation des indicateur directionnels en base generalisee
       if(affiche_local) then
         print*,'Dans crer3d on initialise les indicateur de Trac'
       end if
c      sur toutes les directions generalisee potentiellement actives
 
       do idir=1,NDIMG
         FTHG(idir)=0.d0
         SEFFG(idir)=0.d0
         STOTG(idir)=0.d0
       end do
 
 
c     ------- Base d evaluation des criteres ---------------------------
 
 
c       les retours se font en base principale des contraintes 
c       effectives moyennes
        iphase=-1
        posi=1            
        do idir=1,6
            x6(idir)=SEFFR(iphase,idir,posi)
        end do      
        call x6x33(x6,x33)    
        call b3_v33(x33,x3,V33S)
        call traps1(V33ST,V33S,3)
        do i=1,3
            do j=1,3
                V33(i,j)=v33S(i,j)
                V33T(i,j)=V33ST(i,j) 
            end do
        end do
 
 
c     ----passage dans la base des contraintes generalisee -------------
 
c     ----pour les composantes moyennes des phases ---------------------
 
      do iphase=0,NINC
        if(affiche_local) then
            print*,'Dans crer3d traction radiale phase:', iphase
        end if
c       (ie des deformation elastiques)     
        if(iphase.eq.0) then
            idebut=9*NINC        
        else
            idebut=9*(iphase-1)
        end if
c       recup contraintes totales et changement de base
        do idir=1,6
            x6(idir)=SEFFR(iphase,idir,posi)
        end do
c       projection sur base de l increment de charge       
        call prjc3d(x6,v33,v33t,x3,x16,.false.)
        do idir=1,3
            SEFFG(idebut+idir)=x3(idir)
        end do
c       passage de la deformation plastique en base principales
        do idir=1,6
            x06(idir)=EPSTR(iphase,idir,posi)
        end do                
        call prjc3d(x06,v33,v33t,x03,x16,.false.)
        do idir=1,3
            EPSTG(idebut+idir)=x03(idir)
        end do
      end do
 
c     ---- composantes d interface radiale ----------------------------- 
 
      do iphase=1,NINC
        if(affiche_local) then
            print*,'Criter1inc3d traction radiale interface :',iphase     
        end if       
c       position dans vecteur des contraintes     
        idebut=9*(iphase-1)+3
c       recup contraintes totales et changement de base
        do idir=1,6
            x6(idir)=SEFFI(iphase,idir,posi)
        end do
c       projection sur base de l increment de charge        
        call prjc3d(x6,v33,v33t,x3,x16,.false.)
        do idir=1,3
            SEFFG(idebut+idir)=x3(idir)
        end do
c       passage de la deformation plastique en base principales
        do idir=1,6
            x06(idir)=EPSTI(iphase,idir,posi)
        end do                
        call prjc3d(x06,v33,v33t,x03,x16,.false.)
        do idir=1,3
            EPSTG(idebut+idir)=x03(idir)
        end do       
      end do
 
c     ---- composantes de contrainte totale  dans l inclusion ---------- 
 
      do iphase=1,NINC
        if(affiche_local) then
            print*,'Criter1inc3d traction radiale interface :',iphase     
        end if       
c       position dans vecteur des contraintes     
        idebut=9*(iphase-1)
c       recup contraintes totales et changement de base
        do idir=1,6
            x6(idir)=STOTR(iphase,idir,posi)
        end do
c       projection sur base de l increment de charge        
        call prjc3d(x6,v33,v33t,x3,x16,.false.)
        do idir=1,3
            STOTG(idebut+idir)=x3(idir)
        end do
c       passage de la deformation plastique en base principales
        do idir=1,6
            x06(idir)=EPSLR(iphase,idir,posi)
        end do                
        call prjc3d(x06,v33,v33t,x03,x16,.false.)
        do idir=1,3
            EPSLG(idebut+idir)=x03(idir)
        end do       
      end do
 
c     ---- composantes orthoradiales dans les interfaces ---------------
 
      do iphase=1,NINC      
        if(affiche_local) then
            print*,'Criter1inc3d traction orthoradiale phase:',iphase
        end if        
c       position dans vecteur des contraintes     
        idebut=9*(iphase-1)+6
c       recup contrainte totale et digonalisation
        do idir=1,6
            x6(idir)=SEFFO(iphase,idir,posi)
        end do      
c       projection sur base de l increment de charge         
        call prjc3d(x6,v33,v33t,x3,x16,.false.)
        do idir=1,3
            SEFFG(idebut+idir)=x3(idir)
        end do
c       passage de la deformation plastique en base principales
        do idir=1,6
            x06(idir)=EPSTO(iphase,idir,posi)
        end do                
        call prjc3d(x06,v33,v33t,x03,x16,.false.)
        do idir=1,3
            EPSTG(idebut+idir)=x03(idir)
        end do 
      end do
 
 
c***********************************************************************       
c     Forces thermodynamique 
c***********************************************************************
 
      call ftnc3d(NBRINC,NINC,NDIMG,SEFFG,FRAC,GCH,
     # bwPw,bgPg,FTHG,AFFICHE_local,ERR1)
 
 
c***********************************************************************
c     criteres de traction diffuse
c***********************************************************************
 
       do idir=1,NTMAX
         Log_RTG(idir)=.false.
         ACTIFT(idir)=.false. 
         FT(idir)=0.d0
         do jdir=1,NDIMG
            DPFTDS(idir,jdir)=0.d0
            DGFTDS(idir,jdir)=0.d0
         end do            
       end do
 
      do iphase=0,ninc
        if(iphase.eq.0) then        
            idebut=9*ninc
            do icr=idebut+1,idebut+3
                idir=icr            
                RTG(idir)=RTP(0)
                RFG(idir)=RFP(0)
c               a l interieur des phases on a FTH=SEFF
                call RanKinc3d(NDIMG,NTMAX,icr,idir,
     #          SEFFG,EPSTG,RTG,RFG,SEFFG,FT,DPFTDS,DGFTDS,ACTIFT,
     #          Log_RTG,precision3d,affiche,err1)                
            end do 
        else
c           boucle sur les positions des zones d interet            
            do ityp=1,3
                idebut=(iphase-1)*9+(ityp-1)*3
                do idir=idebut+1,idebut+3
                    icr=idir
                    if (ityp.eq.1) then
c                       ityp=1 radial ds la phase    
                        RTG(idir)=RTP(iphase)
                        RFG(idir)=RFP(iphase)
c                       a l interieur des phases on a FTH=SEFF
                        call RanKinc3d(NDIMG,NTMAX,icr,idir,
     #                  SEFFG,EPSTG,RTG,RFG,SEFFG,FT,DPFTDS,DGFTDS,
     #                  ACTIFT,Log_RTG,precision3d,affiche,err1)
                    else
c                       ityp=2 radial dans l interface cote matrice (non actif)                  
c                       ityp=3 ortho radial ds la matrice (activable)                   
                        RTG(idir)=RTP(0)
                        RFG(idir)=RFP(0)
                        if(ityp.eq.3) then
c                         on evalue le critere orthoradial                          
                          call RanKinc3d(NDIMG,NTMAX,icr,idir,
     #                    SEFFG,EPSTG,RTG,RFG,SEFFG,FT,DPFTDS,DGFTDS,
     #                    ACTIFT,Log_RTG,precision3d,affiche,err1)
                        else
c                         le critere radial de traction cote matrice estvrai
c                         desactive au profit du critere de decollement 
c                         qui est traite dans les citere d'interface localise (eps pl l)
                          ACTIFT(icr)=.false.
                          FT(icr)=0.d0
                          DPFTDS(icr,idir)=0.d0
                          DGFTDS(icr,idir)=0.d0
                        end if
                    end if
                end do
            end do
        end if
      end do
 
 
c     activite de la plasticite de traction diffuse
      PLASTT=.false.
      do icr=1,NTMAX 
        PLASTT=PLASTT.or.ACTIFT(icr)
      end do     
 
 
c     selection des criteres de traction diffuse    
      if(PLASTT) then  
 
c       triage des criteres du plus grand au plus petit ds LNUMCRT 
        call ordr3d(NTMAX,LNUMCRT,FT,affiche_local,ACTIFT,NBRACT)
 
c       affichage des resultats de l analyse des criteres actifs         
        if((NBRACT.ne.0).and.(affiche_local)) then
             print*,'Dans crer3d TRACTION EFFECTIVE PHASES'
             write(*,'(A9,I2,1X,/)')
     #       'NBRACT=',NBRACT
              do icr=1,NTMAX
                 write(*,'(A10,I2,1X,5(A7,E10.3,1X),A5,1X,L1)')
     #           'Direction:',LNUMCRT(icr),
     #           'DGFTDS=',DGFTDS(LNUMCRT(icr),LNUMCRT(icr)),
     #           'DPFTDS=',DPFTDS(LNUMCRT(icr),LNUMCRT(icr)),
     #           'FT=',FT(LNUMCRT(icr)),
     #           'FTHG=',FTHG(LNUMCRT(icr)),
     #           'SEFFG=',SEFFG(LNUMCRT(icr)),
     #           'Activ:',ACTIFT(LNUMCRT(icr))
              end do
c             read*
        end if 
 
       else
 
          NBRACT=0
 
       end if
 
c***********************************************************************
c     criteres de traction par decollement  localise aux interfaces 
c***********************************************************************
c      goto 10
      do idir=1,NLMAX
         Log_RLG(idir)=.false.
         ACTIFL(idir)=.false. 
         do jdir=1,NDIMG
             DPFLDS(idir,jdir)=0.d0
             DGFLDS(idir,jdir)=0.d0
         end do             
         FL(idir)=0.d0
      end do
 
      do iphase=1,ninc 
c         print*,'critere interface phase',iphase 
         idebut=(iphase-1)*9
         do idirg=idebut+1,idebut+3
c             print*,'direction',idirg
             icr=idirg
             RTLG(idirg)=RTI(iphase)
             RFLG(idirg)=RFI(iphase)
c            aux interfaces FTH=STOTG
            call RanKinc3d(NDIMG,NLMAX,icr,idirg,
     #      STOTG,EPSLG,RTLG,RFLG,STOTG,FL,DPFLDS,DGFLDS,ACTIFL,
     #      Log_RLG,precision3d,affiche,err1)
         end do
      end do
 
c     activite de la plasticite d interface localisee 
      PLASTL=.false.
      do icr=1,NLMAX 
        PLASTL=PLASTL.or.ACTIFL(icr)
      end do     
 
c     selection des criteres a traiter      
      if(PLASTL) then  
 
c       triage des criteres du plus grand au plus petit ds LNUMCRT 
        call ordr3d(NLMAX,LNUMCRL,FL,affiche_local,ACTIFL,NBRACL)
 
 
c       affichage des resultats de l analyse des criteres actifs         
        if((NBRACL.ne.0).and.(affiche_local)) then
             print*,'Dans crer3d TRACTION TOTALE INTERFACES'
             write(*,'(A9,I2,1X,/)')
     #       'NBRACL=',NBRACL
              do icr=1,NLMAX
                 write(*,'(A10,I2,1X,5(A7,E10.3,1X),A5,1X,L1)')
     #           'Direction:',LNUMCRL(icr),
     #           'DGFTDS=',DGFLDS(LNUMCRL(icr),LNUMCRL(icr)),
     #           'DPFTDS=',DPFLDS(LNUMCRL(icr),LNUMCRL(icr)),
     #           'FT=',FL(LNUMCRL(icr)),
     #           'FTHG=',FTHG(LNUMCRL(icr)),
     #           'STOTG=',STOTG(LNUMCRL(icr)),
     #           'Activ:',ACTIFL(LNUMCRL(icr))
              end do
c             read*
        end if 
 
       else
 
         NBRACL=0
 
       end if
 
 
c***********************************************************************
c     criteres de cisaillement 
c***********************************************************************
 
c     ---- initialisations 
 
c       lors de l evalauation en 0 il ne faut pas effacer 
c       les directions, ni la liste des criteres actifs 
 
10      PLASTC=.false.
        TAUEQ=0.d0
        PRESS=0.d0
        NBRACC=0
c       direction pour le critere le plus actif en cisaillement
        do icr=1,NCMAX
            FC(icr)=0.d0
            ACTIFC(icr)=.false.        
            do idir=1,NDIMG
                DGFCDS(icr,idir)=0.d0
                DPFCDS(icr,idir)=0.d0                
            end do 
        end do  
 
c     ---- desactivation possible de DP --------------------------------
 
      if((NBRACT.gt.0).or.(NBRACL.gt.0)) then
c       on n examine pas les cisaillement tant que des 
c       tractions sont actifs      
        goto 40
      end if
 
c     ------------------------------------------------------------------ 
 
      if(affiche_local) then
            print*,'Dans crer3d on demarre le cisaillement'
      end if 
c     boucle sur les phases 
 
      do iphase=0,ninc
 
c       affectation d un numero au critere
        if(iphase.eq.0) then
 
c           -- matrice  en cisaillement --------------------------------
 
            IDEBUT=9*ninc
c           type de critere (1 dans la phase, 2 dans l interface)
            ityp=1            
c           numero du critere            
            icr=ntypc*ninc+ityp
 
c          recuperation des contraintes         
            do idir=1,3
                SEFF3(idir)=SEFFG(IDEBUT+IDIR)
                FTH3(idir)=FTHG(IDEBUT+IDIR)
            end do
c           initialisation de l activite du critere et de sa valeur            
            actifc1=.false.
            FC1=0.d0
c           calcul du critere de cisaillement                
            call dpin3d(NBRINC,NCMAX,NDIMG,IPHASE,
     #      BETA,DELTA,COHE,RTP,SEFF3,FC1,DPFCDS3,DGFCDS3,
     #      ACTIFC1,FTH3,PRECISION3D,AFFICHE,ERR1)
c           stockage ds la variables generalisees     
            if(ACTIFC1) then
                ACTIFC(icr)=.true.
                NBRACC=NBRACC+1
                LNUMCRC(NBRACC)=icr
                FC(icr)=FC1
                do idir=1,3
                    DPFCDS(icr,idebut+idir)=DPFCDS3(idir)
                    DGFCDS(icr,idebut+idir)=DGFCDS3(idir)                      
                end do   
            else
                ACTIFC(icr)=.false.            
            end if
 
        else 
 
c           ---- inclusion en cisaillement -----------------------------
 
            IDEBUT=9*(iphase-1)
c           type de critere (1 dans la phase, 2 dans l interface)
            ityp=1            
c           numero du critere            
            icr=ntypc*(iphase-1)+ityp
 
c          recuperation des contraintes         
            do idir=1,3
                SEFF3(idir)=SEFFG(IDEBUT+IDIR)
                FTH3(idir)=FTHG(IDEBUT+IDIR)
            end do
c           initialisation de l activite du critere et de sa valeur            
            actifc1=.false.
            FC1=0.d0
c           calcul du critere de cisaillement                
            call dpin3d(NBRINC,NCMAX,NDIMG,IPHASE,
     #      BETA,DELTA,COHE,RTP,SEFF3,FC1,DPFCDS3,DGFCDS3,
     #      ACTIFC1,FTH3,PRECISION3D,AFFICHE,ERR1)
c           stockage ds la variables generalisees     
            if(ACTIFC1) then
                ACTIFC(icr)=.true.
                NBRACC=NBRACC+1
                LNUMCRC(NBRACC)=icr
                FC(icr)=FC1
                do idir=1,3
                    DPFCDS(icr,idebut+idir)=DPFCDS3(idir)
                    DGFCDS(icr,idebut+idir)=DGFCDS3(idir)
                end do
            else
                ACTIFC(icr)=.false.
            end if
 
c           --- interface de l inclusion en cisaillement ---------------            
 
c           type de critere (1 dans la phase, 2 dans l interface)
            ityp=2            
c           numero du critere            
            icr=ntypc*(iphase-1)+ityp
 
c           recuperation des contraintes au sens thermodynamique         
            do idir=1,3
                SEFF3(idir)=
     #          (SEFFG(IDEBUT+3+IDIR)+2.d0*SEFFG(IDEBUT+6+IDIR))/3.d0
                FTH3(idir)=
     #          (FTHG(IDEBUT+3+IDIR)+2.d0*FTHG(IDEBUT+6+IDIR))/3.d0
            end do
c           initialisation de l activite du critere et de sa valeur            
            actifc1=.false.
            FC1=0.d0
c           calcul du critere de cisaillement                
            call dpin3d(NBRINC,NCMAX,NDIMG,IPHASE,
     #      BETA,DELTA,COHE,RTP,SEFF3,FC1,DPFCDS3,DGFCDS3,
     #      ACTIFC1,FTH3,PRECISION3D,AFFICHE,ERR1)
c           stockage ds la variables generalisees     
            if(ACTIFC1) then
                do idir=1,3
                    DPFCDS(icr,idebut+3+idir)=DPFCDS3(idir)/3.d0
                    DPFCDS(icr,idebut+6+idir)=DPFCDS3(idir)*2.d0/3.d0
                    DGFCDS(icr,idebut+3+idir)=DGFCDS3(idir)/3.d0
                    DGFCDS(icr,idebut+6+idir)=DGFCDS3(idir)*2.d0/3.d0
                end do
                ACTIFC(icr)=.true.
                NBRACC=NBRACC+1
                LNUMCRC(NBRACC)=icr
                FC(icr)=FC1 
            else
                ACTIFC(icr)=.false.               
            end if
 
        end if       
 
        if(err1.eq.1) then
            print*,'Erreur dans le calcul du cisaillement'
            print*,'dans crer3d'
            return
        end if
 
      end do
 
c     bilan des activites de cisaillement
      if(NBRACC.ne.0) then      
         PLASTC=.true.
      end if      
 
      if(PLASTC) then
 
c       -- triage des criteres du plus grand au plus petit ds LNUMCRC --
 
        call ordr3d(NCMAX,LNUMCRC,FC,affiche_local,ACTIFC,NBRACC)
 
 
c       --- affichage des resultats de l analyse des criteres actifs ---
 
        if((NBRACC.ne.0).and.(affiche_local)) then
             print*,'Dans crer3d'
             write(*,'(A9,I2)') 'NBRACC=',NBRACC      
             do icr=1,NBRACC
                 write(*,'(A10,I2,1X,A3,E10.3,1X,A6,1X,L1)')
     #           'Citere:',LNUMCRC(icr),
     #           'FC=',FT(LNUMCRC(icr)),
     #           'Activ:',ACTIFC(LNUMCRC(icr))
                 do idir=1,NDIMG
                     write(*,'(A9,I3,1X,2(A7,E10.3,1X))')
     #               'Direction',idir,
     #               'DGFCDS=',DGFCDS(LNUMCRC(icr),idir),
     #               'DPFCDS=',DPFCDS(LNUMCRC(icr),idir)
                 end do
             end do
c            read*
        end if 
 
       else
 
         NBRACC=0
 
       end if       
 
c     ------- fin zone critere de cisaillement -------------------------
 
 
c***********************************************************************
c     bilan de franchissement des criteres
c***********************************************************************
 
40    plast=plastt.or.plastc.or.plastl
 
c     initialisation critere global
      resg=0.d0 
 
c     evaluation du critere global
      if(plast) then
 
        if(plastt) then
            do i=1,nbract
                j=lnumcrt(i)
                RESG=RESG+FT(j)**2
            end do
C             fmaxt=ft(lnumcrt(1))
C             fmint=ft(lnumcrt(nbract))
C             coefft=fmint/fmaxt
C             fmaxg=fmaxt
C             fming=fmint
C             coeffg=coefft
        end if
 
        if(plastl) then
            do i=1,nbracl
                j=lnumcrl(i)
                RESG=RESG+FL(j)**2
            end do
C             fmaxl=fl(lnumcrl(1))
C             fminl=fl(lnumcrl(nbract))
C             coeffl=fminl/fmaxl
C             fmaxg=fmaxl
C             fming=fminl
C             coeffg=coeffl
        end if
 
        if(plastc) then          
            do i=1,nbracc
                j=lnumcrc(i)
                RESG=RESG+FC(j)**2
            end do
C             fmaxc=fc(lnumcrc(1))
C             fminc=fc(lnumcrc(nbracc))
C             coeffc=fminc/fmaxc   
C             fmaxg=fmaxc
C             fming=fminc 
C             coeffg=coeffc
        end if
 
        RESG=sqrt(RESG)
 
C         if(plastc.and.plastt) then
C             coeffg=min(coefft,coeffc)
 
C         else if(plastc.and.plastl) then
C             coeffg=min(coeffl,coeffc)
 
C         else if(plastt.and.plastl) then
C             coeffg=min(coeffl,coefft)
 
C         else if(plastt.and.plastl.and.plastc) then
C             coeffg=min(coeffl,coefft,coeffc)
C         end if
 
c        reduc_resg=coeffg
        reduc_resg=1.d0
c        print*,'Coefficient de relaxation du residu:',reduc_resg
      end if        
 
      if(affiche_local) then     
        print*,'Fin de crer3d'
*       read*
      end if        
 
      return
 
      end
 
 
 
 
 
 
 
 

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