constitutitions of "j2" and "dislobased" are referenced in constitutive and included to the Marc interface

trunk/constitutive.f90

@@ -81,6 +81,14 @@ subroutine constitutive_init()
            constitutive_sizeDotState(g,i,e) =    constitutive_j2_sizeDotState(myInstance)
            constitutive_sizeState(g,i,e) =       constitutive_j2_sizeState(myInstance)
            constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(myInstance)
+         case (constitutive_dislobased_label)
+           allocate(constitutive_state_old(g,i,e)%p(constitutive_dislobased_sizeState(myInstance)))
+           allocate(constitutive_state_new(g,i,e)%p(constitutive_dislobased_sizeState(myInstance)))
+           constitutive_state_new(g,i,e)%p = constitutive_dislobased_stateInit(g,i,e)
+           constitutive_state_old(g,i,e)%p = constitutive_dislobased_stateInit(g,i,e)
+           constitutive_sizeDotState(g,i,e) =    constitutive_dislobased_sizeDotState(myInstance)
+           constitutive_sizeState(g,i,e) =       constitutive_dislobased_sizeState(myInstance)
+           constitutive_sizePostResults(g,i,e) = constitutive_dislobased_sizePostResults(myInstance)
          case default
            call IO_error(200,material_phase(g,i,e))      ! unknown constitution
        end select
@@ -237,7 +245,7 @@ function constitutive_dotState(Tstar_v,Temperature,ipc,ip,el)
    case (constitutive_j2_label)
      constitutive_dotState = constitutive_j2_dotState(Tstar_v,Temperature,constitutive_state_new,ipc,ip,el)
    case (constitutive_dislobased_label)
-     call constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state_new,ipc,ip,el)
+     constitutive_dotState = constitutive_dislobased_dotState(Tstar_v,Temperature,constitutive_state_new,ipc,ip,el)
 
  end select
  return
@@ -274,7 +282,7 @@ pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el)
    case (constitutive_j2_label)
      constitutive_postResults = constitutive_j2_postResults(Tstar_v,Temperature,dt,constitutive_state_new,ipc,ip,el)
    case (constitutive_dislobased_label)
-     call constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state_new,ipc,ip,el)
+     constitutive_postResults = constitutive_dislobased_postResults(Tstar_v,Temperature,dt,constitutive_state_new,ipc,ip,el)
 
  end select
 

trunk/constitutive_j2.f90

@@ -0,0 +1,391 @@
+
+!*****************************************************
+!*      Module: CONSTITUTIVE_J2        				 *
+!*****************************************************
+!* contains:                                         *
+!* - constitutive equations                          *
+!* - parameters definition                           *
+!*****************************************************
+
+!	[Alu]
+!	constitution            j2
+!	(output)                flowstress
+!	(output)                strainrate
+!	c11                     110.9e9    # (3 C11 + 2 C12 + 2 C44) / 5
+!	c12                     58.34e9    # (1 C11 + 4 C12 - 1 C44) / 5
+!	taylorfactor            3
+!	s0                      31e6
+!	gdot0                   0.001
+!	n                       20
+!	h0                      75e6
+!	s_sat                   63e6
+!	w0                      2.25
+
+MODULE constitutive_j2
+
+!*** Include other modules ***
+ use prec, only: pReal,pInt
+ implicit none
+
+ character (len=*), parameter :: constitutive_j2_label = 'j2'
+ 
+ integer(pInt),   dimension(:),     allocatable :: constitutive_j2_sizeDotState, &
+                                                   constitutive_j2_sizeState, &
+                                                   constitutive_j2_sizePostResults
+ character(len=64), dimension(:,:), allocatable :: constitutive_j2_output
+ real(pReal), dimension(:),     allocatable :: constitutive_j2_C11
+ real(pReal), dimension(:),     allocatable :: constitutive_j2_C12
+ real(pReal), dimension(:,:,:), allocatable :: constitutive_j2_Cslip_66
+!* Visco-plastic constitutive_j2 parameters
+ real(pReal), dimension(:),     allocatable :: constitutive_j2_fTaylor
+ real(pReal), dimension(:),     allocatable :: constitutive_j2_s0
+ real(pReal), dimension(:),     allocatable :: constitutive_j2_gdot0
+ real(pReal), dimension(:),     allocatable :: constitutive_j2_n
+ real(pReal), dimension(:),     allocatable :: constitutive_j2_h0
+ real(pReal), dimension(:),     allocatable :: constitutive_j2_s_sat
+ real(pReal), dimension(:),     allocatable :: constitutive_j2_w0
+
+
+CONTAINS
+!****************************************
+!* - constitutive_init
+!* - constitutive_homogenizedC
+!* - constitutive_microstructure
+!* - constitutive_LpAndItsTangent
+!* - consistutive_dotState
+!* - consistutive_postResults
+!****************************************
+
+
+subroutine constitutive_j2_init(file)
+!**************************************
+!*      Module initialization         *
+!**************************************
+ use prec, only: pInt, pReal
+ use math, only: math_Mandel3333to66, math_Voigt66to3333
+ use IO
+ use material
+ integer(pInt), intent(in) :: file
+ integer(pInt), parameter :: maxNchunks = 7
+ integer(pInt), dimension(1+2*maxNchunks) :: positions
+ integer(pInt) section, maxNinstance, i,j,k,l, output
+ character(len=64) tag
+ character(len=1024) line
+ 
+ maxNinstance = count(phase_constitution == constitutive_j2_label)
+ if (maxNinstance == 0) return
+
+ allocate(constitutive_j2_sizeDotState(maxNinstance)) ;   constitutive_j2_sizeDotState = 0_pInt
+ allocate(constitutive_j2_sizeState(maxNinstance)) ;      constitutive_j2_sizeState = 0_pInt
+ allocate(constitutive_j2_sizePostResults(maxNinstance)); constitutive_j2_sizePostResults = 0_pInt
+ allocate(constitutive_j2_output(maxval(phase_Noutput), &
+                                               maxNinstance)) ;         constitutive_j2_output = ''
+ allocate(constitutive_j2_C11(maxNinstance)) ;            constitutive_j2_C11 = 0.0_pReal
+ allocate(constitutive_j2_C12(maxNinstance)) ;            constitutive_j2_C12 = 0.0_pReal
+ allocate(constitutive_j2_Cslip_66(6,6,maxNinstance)) ;   constitutive_j2_Cslip_66 = 0.0_pReal
+ allocate(constitutive_j2_fTaylor(maxNinstance)) ;        constitutive_j2_fTaylor = 0.0_pReal
+ allocate(constitutive_j2_s0(maxNinstance)) ;             constitutive_j2_s0 = 0.0_pReal
+ allocate(constitutive_j2_gdot0(maxNinstance)) ;          constitutive_j2_gdot0 = 0.0_pReal
+ allocate(constitutive_j2_n(maxNinstance)) ;              constitutive_j2_n = 0.0_pReal
+ allocate(constitutive_j2_h0(maxNinstance)) ;             constitutive_j2_h0 = 0.0_pReal
+ allocate(constitutive_j2_s_sat(maxNinstance)) ;          constitutive_j2_s_sat = 0.0_pReal
+ allocate(constitutive_j2_w0(maxNinstance)) ;             constitutive_j2_w0 = 0.0_pReal
+ 
+ rewind(file)
+ line = ''
+ section = 0
+ 
+ do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase')     ! wind forward to <phase>
+   read(file,'(a1024)',END=100) line
+ enddo
+
+ do                                                       ! read thru sections of phase part
+   read(file,'(a1024)',END=100) line
+   if (IO_isBlank(line)) cycle                            ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') exit                ! stop at next part
+   if (IO_getTag(line,'[',']') /= '') then                ! next section
+     section = section + 1
+     output = 0                                           ! reset output counter
+   endif
+   if (section > 0 .and. phase_constitution(section) == constitutive_j2_label) then  ! one of my sections
+     i = phase_constitutionInstance(section)     ! which instance of my constitution is present phase
+     positions = IO_stringPos(line,maxNchunks)
+     tag = IO_lc(IO_stringValue(line,positions,1))        ! extract key
+     select case(tag)
+       case ('(output)')
+         output = output + 1
+         constitutive_j2_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
+       case ('c11')
+              constitutive_j2_C11(i) = IO_floatValue(line,positions,2)
+       case ('c12')
+              constitutive_j2_C12(i) = IO_floatValue(line,positions,2)
+       case ('s0')
+              constitutive_j2_s0(i) = IO_floatValue(line,positions,2)
+       case ('gdot0')
+              constitutive_j2_gdot0(i) = IO_floatValue(line,positions,2)
+       case ('n')
+              constitutive_j2_n(i) = IO_floatValue(line,positions,2)
+       case ('h0')
+              constitutive_j2_h0(i) = IO_floatValue(line,positions,2)
+       case ('s_sat')
+              constitutive_j2_s_sat(i) = IO_floatValue(line,positions,2)
+       case ('w0')
+              constitutive_j2_w0(i) = IO_floatValue(line,positions,2)
+       case ('taylorfactor')
+              constitutive_j2_fTaylor(i) = IO_floatValue(line,positions,2)
+     end select
+   endif
+ enddo
+
+100 do i = 1,maxNinstance                                        ! sanity checks
+   if (constitutive_j2_s0(i) < 0.0_pReal)          call IO_error(203)
+   if (constitutive_j2_gdot0(i) <= 0.0_pReal)      call IO_error(204)
+   if (constitutive_j2_n(i) <= 0.0_pReal)          call IO_error(205)
+   if (constitutive_j2_h0(i) <= 0.0_pReal)         call IO_error(206)
+   if (constitutive_j2_s_sat(i) <= 0.0_pReal)      call IO_error(207)
+   if (constitutive_j2_w0(i) <= 0.0_pReal)         call IO_error(208)
+   if (constitutive_j2_fTaylor(i) <= 0.0_pReal)    call IO_error(240)
+ enddo
+
+ do i = 1,maxNinstance
+   constitutive_j2_sizeDotState(i) = 1
+   constitutive_j2_sizeState(i)    = 1
+
+   do j = 1,maxval(phase_Noutput)
+     select case(constitutive_j2_output(j,i))
+       case('flowstress')
+         constitutive_j2_sizePostResults(i) = &
+         constitutive_j2_sizePostResults(i) + 1
+       case('strainrate')
+         constitutive_j2_sizePostResults(i) = &
+         constitutive_j2_sizePostResults(i) + 1
+     end select
+   enddo
+
+   forall(k=1:3)
+     forall(j=1:3) &
+       constitutive_j2_Cslip_66(k,j,i) =     constitutive_j2_C12(i)
+       constitutive_j2_Cslip_66(k,k,i) =     constitutive_j2_C11(i)
+       constitutive_j2_Cslip_66(k+3,k+3,i) = 0.5_pReal*(constitutive_j2_C11(i)-constitutive_j2_C12(i))
+   end forall
+   constitutive_j2_Cslip_66(:,:,i) = &
+     math_Mandel3333to66(math_Voigt66to3333(constitutive_j2_Cslip_66(:,:,i)))
+
+ enddo
+
+ return
+
+end subroutine
+
+
+function constitutive_j2_stateInit(ipc,ip,el)
+!*********************************************************************
+!* initial microstructural state                                     *
+!*********************************************************************
+ use prec, only: pReal,pInt
+ use material, only: material_phase, phase_constitutionInstance
+ implicit none
+
+!* Definition of variables
+ integer(pInt), intent(in) :: ipc,ip,el
+ integer(pInt) matID
+ real(pReal), dimension(1) :: &
+   constitutive_j2_stateInit
+
+ matID = phase_constitutionInstance(material_phase(ipc,ip,el))
+ constitutive_j2_stateInit = constitutive_j2_s0(matID)
+
+ return
+end function
+
+
+function constitutive_j2_homogenizedC(state,ipc,ip,el)
+!*********************************************************************
+!* homogenized elacticity matrix                                     *
+!* INPUT:                                                            *
+!*  - state           : state variables                              *
+!*  - ipc             : component-ID of current integration point    *
+!*  - ip              : current integration point                    *
+!*  - el              : current element                              *
+!*********************************************************************
+ use prec, only: pReal,pInt,p_vec
+ use mesh, only: mesh_NcpElems,mesh_maxNips
+ use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
+ implicit none
+
+!* Definition of variables
+ integer(pInt), intent(in) :: ipc,ip,el
+ integer(pInt) matID
+ real(pReal), dimension(6,6) :: constitutive_j2_homogenizedC
+ type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
+ 
+ matID = phase_constitutionInstance(material_phase(ipc,ip,el))
+ constitutive_j2_homogenizedC = constitutive_j2_Cslip_66(:,:,matID)
+
+ return
+
+end function
+
+
+subroutine constitutive_j2_microstructure(Temperature,state,ipc,ip,el)
+!*********************************************************************
+!* calculate derived quantities from state (not used here)           *
+!* INPUT:                                                            *
+!*  - Tp              : temperature                                  *
+!*  - ipc             : component-ID of current integration point    *
+!*  - ip              : current integration point                    *
+!*  - el              : current element                              *
+!*********************************************************************
+ use prec, only: pReal,pInt,p_vec
+ use mesh, only: mesh_NcpElems,mesh_maxNips
+ use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
+ implicit none
+
+!* Definition of variables
+ integer(pInt) ipc,ip,el, matID
+ real(pReal) Temperature
+ type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
+
+ matID = phase_constitutionInstance(material_phase(ipc,ip,el))
+
+end subroutine
+
+
+subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,state,ipc,ip,el)
+!*********************************************************************
+!* plastic velocity gradient and its tangent                         *
+!* INPUT:                                                            *
+!*  - Tstar_v         : 2nd Piola Kirchhoff stress tensor (Mandel)   *
+!*  - ipc             : component-ID at current integration point    *
+!*  - ip              : current integration point                    *
+!*  - el              : current element                              *
+!* OUTPUT:                                                           *
+!*  - Lp              : plastic velocity gradient                    *
+!*  - dLp_dTstar      : derivative of Lp (4th-rank tensor)           *
+!*********************************************************************
+ use prec, only: pReal,pInt,p_vec
+ use math, only: math_mul6x6,math_Mandel6to33,math_Plain3333to99
+ use lattice, only: lattice_Sslip,lattice_Sslip_v
+ use mesh, only: mesh_NcpElems,mesh_maxNips
+ use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
+
+ implicit none
+
+!* Definition of variables
+ integer(pInt) ipc,ip,el
+ integer(pInt) matID,i,k,l,m,n
+ real(pReal) Temperature
+ type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
+ real(pReal), dimension(6) :: Tstar_v
+ real(pReal), dimension(3,3) :: Lp
+ real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333
+ real(pReal), dimension(9,9) :: dLp_dTstar
+ real(pReal) norm_Tstar
+ matID = phase_constitutionInstance(material_phase(ipc,ip,el))
+
+ norm_Tstar = dsqrt(math_mul6x6(Tstar_v,Tstar_v))
+ 
+!* Calculation of Lp
+ Lp = math_Mandel6to33(Tstar_v)/norm_Tstar*constitutive_j2_gdot0(matID)/constitutive_j2_fTaylor(matID)* &
+      (dsqrt(1.5_pReal)/constitutive_j2_fTaylor(matID)*norm_Tstar/state(ipc,ip,el)%p(1))**constitutive_j2_n(matID)
+
+!* Calculation of the tangent of Lp
+ dLp_dTstar3333 = 0.0_pReal
+ dLp_dTstar = 0.0_pReal
+ dLp_dTstar = math_Plain3333to99(dLp_dTstar3333)
+
+ return
+end subroutine
+
+
+function constitutive_j2_dotState(Tstar_v,Temperature,state,ipc,ip,el)
+!*********************************************************************
+!* rate of change of microstructure                                  *
+!* INPUT:                                                            *
+!*  - Tstar_v         : 2nd Piola Kirchhoff stress tensor (Mandel)   *
+!*  - ipc             : component-ID at current integration point    *
+!*  - ip              : current integration point                    *
+!*  - el              : current element                              *
+!* OUTPUT:                                                           *
+!*  - constitutive_dotState : evolution of state variable            *
+!*********************************************************************
+ use prec, only: pReal,pInt,p_vec
+ use math, only: math_mul6x6
+ use lattice, only: lattice_Sslip_v
+ use mesh, only: mesh_NcpElems,mesh_maxNips
+ use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
+ implicit none
+
+!* Definition of variables
+ integer(pInt) ipc,ip,el
+ integer(pInt) matID
+ real(pReal) Temperature
+ type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
+ real(pReal), dimension(6) :: Tstar_v
+ real(pReal), dimension(1) :: constitutive_j2_dotState
+ real(pReal) norm_Tstar
+
+ matID = phase_constitutionInstance(material_phase(ipc,ip,el))
+
+ norm_Tstar = dsqrt(math_mul6x6(Tstar_v,Tstar_v))
+ constitutive_j2_dotState = constitutive_j2_gdot0(matID)/constitutive_j2_fTaylor(matID)* &
+                            (dsqrt(1.5_pReal)/constitutive_j2_fTaylor(matID)*norm_Tstar/state(ipc,ip,el)%p(1))** &
+							constitutive_j2_n(matID) * &
+							constitutive_j2_h0(matID)*(1.0_pReal-state(ipc,ip,el)%p(1)/constitutive_j2_s_sat(matID))** &
+							constitutive_j2_w0(matID)
+
+ return
+
+end function
+
+
+pure function constitutive_j2_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el)
+!*********************************************************************
+!* return array of constitutive results                              *
+!* INPUT:                                                            *
+!*  - Tstar_v         : 2nd Piola Kirchhoff stress tensor (Mandel)   *
+!*  - dt              : current time increment                       *
+!*  - ipc             : component-ID at current integration point    *
+!*  - ip              : current integration point                    *
+!*  - el              : current element                              *
+!*********************************************************************
+ use prec, only: pReal,pInt,p_vec
+ use math, only: math_mul6x6
+ use lattice, only: lattice_Sslip_v
+ use mesh, only: mesh_NcpElems,mesh_maxNips
+ use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput
+ implicit none
+
+!* Definition of variables
+ integer(pInt), intent(in) :: ipc,ip,el
+ real(pReal), intent(in) :: dt,Temperature
+ real(pReal), dimension(6), intent(in) :: Tstar_v
+ type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state
+ integer(pInt) matID,o,i,c,n
+ real(pReal) norm_Tstar
+ real(pReal), dimension(constitutive_j2_sizePostResults(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
+   constitutive_j2_postResults
+
+ matID = phase_constitutionInstance(material_phase(ipc,ip,el))
+ norm_Tstar = dsqrt(math_mul6x6(Tstar_v,Tstar_v))
+ c = 0_pInt
+ constitutive_j2_postResults = 0.0_pReal
+
+ do o = 1,phase_Noutput(material_phase(ipc,ip,el))
+   select case(constitutive_j2_output(o,matID))
+     case ('flowstress')
+       constitutive_j2_postResults(c+1) = state(ipc,ip,el)%p(1)
+       c = c + 1
+     case ('strainrate')
+	   constitutive_j2_postResults(c+1) = constitutive_j2_gdot0(matID)/constitutive_j2_fTaylor(matID)* &
+                                          (dsqrt(1.5_pReal)/constitutive_j2_fTaylor(matID)*norm_Tstar/state(ipc,ip,el)%p(1))** &
+										  constitutive_j2_n(matID)
+       c = c + 1
+   end select
+ enddo
+ 
+ return
+
+end function
+
+END MODULE

trunk/mpie_cpfem_marc2007r1.f90

@@ -36,6 +36,8 @@
  include "lattice.f90"          ! uses prec, math
  include "material.f90"         ! uses prec, math, IO, mesh
  include "constitutive_phenomenological.f90"     ! uses prec, math, IO, lattice, material, debug
+ include "constitutive_j2.f90"                   ! uses prec, math, IO, lattice, material, debug
+ include "constitutive_dislobased.f90"           ! uses prec, math, IO, lattice, material, debug
  include "constitutive.f90"     ! uses prec, IO, math, lattice, mesh, debug
  include "CPFEM.f90"            ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite
 

trunk/mpie_cpfem_marc2008r1.f90

@@ -36,6 +36,8 @@
  include "lattice.f90"          ! uses prec, math
  include "material.f90"         ! uses prec, math, IO, mesh
  include "constitutive_phenomenological.f90"     ! uses prec, math, IO, lattice, material, debug
+ include "constitutive_j2.f90"                   ! uses prec, math, IO, lattice, material, debug
+ include "constitutive_dislobased.f90"           ! uses prec, math, IO, lattice, material, debug
  include "constitutive.f90"     ! uses prec, IO, math, lattice, mesh, debug
  include "CPFEM.f90"            ! uses prec, math, mesh, constitutive, FEsolving, debug, lattice, IO, crystallite