MAKE_INVARIANTS_AND_MAIN_STRAIN.f90

Go to the documentation of this file.

!    Copyright (C) 2012 The SPEED FOUNDATION
!    Author: Ilario Mazzieri
!
!    This file is part of SPEED.
!
!    SPEED is free software; you can redistribute it and/or modify it
!    under the terms of the GNU Affero General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    SPEED is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!    Affero General Public License for more details.
!
!    You should have received a copy of the GNU Affero General Public License
!    along with SPEED.  If not, see <http://www.gnu.org/licenses/>.
 
 
 
 
      subroutine make_invariants_and_main_strain(nn,&                                         
                                                duxdx,duydx,duzdx,&         
                                                duxdy,duydy,duzdy,&         
                                                duxdz,duydz,duzdz,&         
                                                R,yon)
                                            
      implicit none
 
 
      integer*4 :: nn
      integer*4 :: i,j,k
 
      integer*4, dimension(nn,nn,nn) :: yon        
 
      real*8 :: exx,eyy,ezz,exy,eyz,ezx
      real*8 :: r1,r2,r3                        
      real*8 :: invar1,invar2,invar3
      real*8 :: strain_i,strain_ii,strain_iii
      real*8 :: r_average_elem
 
      real*8, dimension(nn,nn,nn) :: duxdx,duxdy,duxdz,duydx,duydy,duydz,duzdx,duzdy,duzdz
      real*8, dimension(nn,nn,nn) :: r
 
 
      r_average_elem = 0.0d0
      
      do i = 1,nn
         do j = 1,nn
            do k = 1,nn
 
                if (yon(i,j,k).eq.1) then
                
                        !+-----------------------------------------
                        !! Invariants of the strain tensor sigma_ij
 
                        exx = duxdx(i,j,k)
                        eyy = duydy(i,j,k)
                        ezz = duzdz(i,j,k)
 
                        exy = .5 * (duxdy(i,j,k) + duydx(i,j,k))
                        eyz = .5 * (duydz(i,j,k) + duzdy(i,j,k))
                        ezx = .5 * (duzdx(i,j,k) + duxdz(i,j,k))
 
                        invar1 = exx + eyy + ezz
 
                        invar2 = exx * eyy + &
                                 eyy * ezz + &
                                 ezz* exx - &
                                 exy**2 - &
                                 eyz**2 - &
                                 ezx**2
                                                                  
                        invar3 = exx * eyy * ezz + &
                                 2 * exy * eyz * (-ezx) - &  
                                 exx * eyz**2 - &
                                 eyy * ezx**2 - &
                                 ezz * exy**2 
 
                        !+-----------------------------------------
                        !! Principal strain computations
 
                        call cubic(-invar1,invar2,-invar3,&
                                  strain_i,strain_ii,strain_iii)
 
                        !call ZerosPolyO3(-invar1,invar2,-invar3,&
                        !          strain_I,strain_II,strain_III)
 
 
                        !+-----------------------------------------
                        !! Max shear strain computation
 
                        r(i,j,k) = (max(strain_i,strain_ii,strain_iii) - &
                                        min(strain_i,strain_ii,strain_iii) ) / 2
 
                        r_average_elem = r_average_elem + r(i,j,k)
                                        
                        !! 
                        !+-----------------------------------------
 
                endif !if (yon(i,j,k).eq.1) then
                                          
            enddo
         enddo
      enddo
 
      r_average_elem = r_average_elem / (nn*nn*nn)
 
      do i = 1,nn
        do j = 1,nn
           do k = 1,nn
 
                r(i,j,k) = r_average_elem
 
           enddo
        enddo
      enddo
 
      return
      
      end subroutine make_invariants_and_main_strain