get F from crystallite module

code/constitutive.f90

@@ -488,7 +488,7 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
    case (PLASTICITY_NONLOCAL_ID) plasticityType
      call plastic_nonlocal_microstructure (Fe,Fp,ip,el)
    case (PLASTICITY_PHENOPLUS_ID) plasticityType
-     call plastic_phenoplus_microstructure(orientations,ipc,ip,el)
+     call plastic_phenoplus_microstructure(orientations,ipc,ip,el,Fe,Fp)
  end select plasticityType
 
 end subroutine constitutive_microstructure

code/plastic_phenoplus.f90

@@ -557,7 +557,7 @@ subroutine plastic_phenoplus_init(fileUnit)
 ! allocate state arrays
      sizeState = plastic_phenoplus_totalNslip(instance) &                         ! s_slip
                + plastic_phenoplus_totalNtwin(instance) &                         ! s_twin
-               + 2_pInt &                                                         ! sum(gamma) + sum(f)
+               + 2_pInt &                                                         ! sum(gamma) + sum(twinVolFrac)
                + plastic_phenoplus_totalNslip(instance) &                         ! accshear_slip
                + plastic_phenoplus_totalNtwin(instance) &                         ! accshear_twin
                + plastic_phenoplus_totalNslip(instance)                           ! kappa
@@ -568,7 +568,7 @@ subroutine plastic_phenoplus_init(fileUnit)
      !           memory leak issue.
      sizeDotState = plastic_phenoplus_totalNslip(instance) &                      ! s_slip
                   + plastic_phenoplus_totalNtwin(instance) &                      ! s_twin
-                  + 2_pInt &                                                      ! sum(gamma) + sum(f)
+                  + 2_pInt &                                                      ! sum(gamma) + sum(twinVolFrac)
                   + plastic_phenoplus_totalNslip(instance) &                      ! accshear_slip
                   + plastic_phenoplus_totalNtwin(instance)                        ! accshear_twin
 
@@ -739,7 +739,7 @@ end subroutine plastic_phenoplus_aTolState
 !--------------------------------------------------------------------------------------------------
 !> @brief calculate push-up factors (kappa) for each voxel based on its neighbors
 !--------------------------------------------------------------------------------------------------
-subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el)
+subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el,Fe,Fp)
  use math, only:  pi, &
                   math_mul33x33, &
                   math_mul3x3, &
@@ -772,7 +772,10 @@ subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el)
  integer(pInt), intent(in) :: &
    ipc, &                                                     !< component-ID of integration point
    ip, &                                                      !< integration point
-   el                                                                                     !< element
+   el
+ real(pReal), dimension(3,3), intent(in) :: &
+                                Fe, &                         ! elastic deformation gradient
+                                Fp                            ! elastic deformation gradient                                                                                     !< element
  real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
                                                 orientation   ! crystal orientation in quaternions
 
@@ -802,10 +805,18 @@ subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el)
  real(pReal)    kappa_max, &                  !
                 tmp_myshear_slip, &           !temp storage for accumulative shear for me
                 mprime_cut, &                 !m' cutoff to consider neighboring effect
+                dtaylor_cut, &                !threshold for determine high contrast interface using Taylor factor
                 avg_acshear_ne, &             !the average accumulative shear from my neighbor
+                taylor_me, &                  !Taylor factor for me
+                taylor_ne, &                  !Taylor factor for my current neighbor
                 tmp_mprime, &                 !temp holder for m' value
                 tmp_acshear                   !temp holder for accumulative shear for m'
 
+ real(pReal), dimension(3,3) :: &
+                Fe_me, &                      !my elastic deformation gradient
+                Fp_me, &                      !my plastic deformation gradient
+                Fe_ne, &                      !elastic deformation gradient of my current neighbor
+                Fp_ne                         !plastic deformation gradient of my current neighbor
 
  real(pReal), dimension(plastic_phenoplus_totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: &
                       m_primes, &                   !m' between me_alpha(one) and neighbor beta(all)
@@ -824,6 +835,10 @@ subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el)
                                                 ne_mprimes                                !m' between each neighbor
 
  !***Get my properties
+ !@TODO
+ ! still need to know how to access the total strain for current material point
+ ! also, need to figure out an efficient way to calculate gamma_dot for the material
+ ! point and its neighbors
  Nneighbors          = FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el))))
  ph                  = phaseAt(ipc,ip,el)                                   !get my phase
  of                  = phasememberAt(ipc,ip,el)                             !get my spatial location offset in memory
@@ -833,88 +848,114 @@ subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el)
  nt                  = plastic_phenoplus_totalNtwin(instance)
  offset_acshear_slip = ns + nt + 2_pInt
  index_kappa         = ns + nt + 2_pInt + ns + nt                           !location of kappa in plasticState
+
+ !***init calculation for given voxel
  mprime_cut          = 0.7_pReal                                            !set by Dr.Bieler
+ dtaylor_cut         = 1.0_pReal                                            !set by Chen, quick test only
 
- !***gather my accumulative shear from palsticState
+ !***calculate my Taylor factor
- FINDMYSHEAR: do j = 1_pInt,ns
-  me_acshear(j) = plasticState(ph)%state(offset_acshear_slip+j, of)
- enddo FINDMYSHEAR
-
  !***gather my orientation and slip systems
  my_orientation      = orientation(1:4, ipc, ip, el)
+ Fe_me               = Fe(ipc,ip,el)
+ Fp_me               = Fp(ipc,ip,el)
  slipNormal(1:3, 1:ns) = lattice_sn(1:3, 1:ns, ph)
  slipDirect(1:3, 1:ns) = lattice_sd(1:3, 1:ns, ph)
- kappa_max             = plastic_phenoplus_kappa_max(instance)                            !maximum pushups allowed (READIN)
 
- !***calculate kappa between me and all my neighbors
+ !***loop into the geometry to figure out who is my closest neighbor
- LOOPMYSLIP: DO me_slip=1_pInt,ns
-  vld_Nneighbors   = Nneighbors
-  tmp_myshear_slip = me_acshear(me_slip)
-  tmp_mprime       = 0.0_pReal                                                            !highest m' from all neighbors
-  tmp_acshear      = 0.0_pReal                                                            !accumulative shear from highest m'
-
-  !***go through my neighbors to find highest m'
  LOOPNEIGHBORS: DO n=1_pInt, Nneighbors
-   neighbor_el          = mesh_ipNeighborhood(1,n,ip,el)
+ !******for each of my neighbor, calculate the Taylor factor
-   neighbor_ip          = mesh_ipNeighborhood(2,n,ip,el)
+  ne_taylor = 1.0
-   neighbor_n           = 1                                                               !It is ipc
+ !*********for the high contrast interface
-   neighbor_of          = phasememberAt( neighbor_n, neighbor_ip, neighbor_el)
+  IF (abs(taylor_ne - taylor_me) > dtaylor_cut) THEN
-   neighbor_ph          = phaseAt( neighbor_n, neighbor_ip, neighbor_el)
+ !*********  gather neighbor orientation and slip systems
-   neighbor_tex         = material_texture(1,neighbor_ip,neighbor_el)
+ !*********  calculate m' (need to loop through all my slip systems as well)
-   neighbor_orientation = orientation(1:4, neighbor_n, neighbor_ip, neighbor_el)          !ipc is always 1.
+ !*********  if m'>mprime_cut kappa=1.5 else 1.0
-   absMisorientation    = lattice_qDisorientation(my_orientation, &
+ !******
-                                                  neighbor_orientation, &
-                                                  0_pInt)                                 !no need for explicit calculation of symmetry
-
-   !***find the accumulative shear for this neighbor
-   LOOPFINDNEISHEAR: DO ne_slip_ac=1_pInt, ns
-    ne_acshear(ne_slip_ac) = plasticState(ph)%state(offset_acshear_slip+ne_slip_ac, &
-                                                    neighbor_of)
-   ENDDO LOOPFINDNEISHEAR
-
-   !***calculate the average accumulative shear and use it as cutoff
-   avg_acshear_ne = SUM(ne_acshear)/ns
-
-   !***
-   IF (ph==neighbor_ph) THEN
-    !***walk through all the
-    LOOPNEIGHBORSLIP: DO ne_slip=1_pInt,ns
-      !***only consider slip system that is active (above average accumulative shear)
-      IF (ne_acshear(ne_slip) > avg_acshear_ne) THEN
-        m_primes(ne_slip) = abs(math_mul3x3(slipNormal(1:3,me_slip), &
-                                math_qRot(absMisorientation, slipNormal(1:3,ne_slip)))) &
-                           *abs(math_mul3x3(slipDirect(1:3,me_slip), &
-                                math_qRot(absMisorientation, slipDirect(1:3,ne_slip))))
-        !***find the highest m' and corresponding accumulative shear
-        IF (m_primes(ne_slip) > tmp_mprime) THEN
-          tmp_mprime  = m_primes(ne_slip)
-          tmp_acshear = ne_acshear(ne_slip)
-        ENDIF
-      ENDIF
-    ENDDO LOOPNEIGHBORSLIP
-
   ELSE
-    ne_mprimes(n)  = 0.0_pReal
-    vld_Nneighbors = vld_Nneighbors - 1_pInt
   ENDIF
-
+ !***end of search
  ENDDO LOOPNEIGHBORS
 
-  !***check if this element close to rim
+ ! !***gather my accumulative shear from palsticState
-  IF (vld_Nneighbors < Nneighbors) THEN
+ ! FINDMYSHEAR: do j = 1_pInt,ns
-   !***rim voxel, no modification allowed
+ !  me_acshear(j) = plasticState(ph)%state(offset_acshear_slip+j, of)
-   plasticState(ph)%state(index_kappa+me_slip, of) = 1.0_pReal
+ ! enddo FINDMYSHEAR
-  ELSE
-   !***patch voxel, started to calculate push up factor for gamma_dot
-   IF ((tmp_mprime > mprime_cut) .AND. (tmp_acshear > tmp_myshear_slip)) THEN
-    plasticState(ph)%state(index_kappa+me_slip, of) = 1.0_pReal / tmp_mprime
-   ELSE
-    !***minimum damping factor is 0.5
-    plasticState(ph)%state(index_kappa+me_slip, of) = 0.5_pReal + tmp_mprime * 0.5_pReal
-   ENDIF
-  ENDIF
 
- ENDDO LOOPMYSLIP
+ ! !***gather my orientation and slip systems
+ ! my_orientation        = orientation(1:4, ipc, ip, el)
+ ! slipNormal(1:3, 1:ns) = lattice_sn(1:3, 1:ns, ph)
+ ! slipDirect(1:3, 1:ns) = lattice_sd(1:3, 1:ns, ph)
+ ! kappa_max             = plastic_phenoplus_kappa_max(instance)                            !maximum pushups allowed (READIN)
+
+ ! !***calculate kappa between me and all my neighbors
+ ! LOOPMYSLIP: DO me_slip=1_pInt,ns
+ !  vld_Nneighbors   = Nneighbors
+ !  tmp_myshear_slip = me_acshear(me_slip)
+ !  tmp_mprime       = 0.0_pReal                                                            !highest m' from all neighbors
+ !  tmp_acshear      = 0.0_pReal                                                            !accumulative shear from highest m'
+
+ !  !***go through my neighbors to find highest m'
+ !  LOOPNEIGHBORS: DO n=1_pInt,Nneighbors
+ !   neighbor_el          = mesh_ipNeighborhood(1,n,ip,el)
+ !   neighbor_ip          = mesh_ipNeighborhood(2,n,ip,el)
+ !   neighbor_n           = 1                                                               !It is ipc
+ !   neighbor_of          = phasememberAt( neighbor_n, neighbor_ip, neighbor_el)
+ !   neighbor_ph          = phaseAt( neighbor_n, neighbor_ip, neighbor_el)
+ !   neighbor_tex         = material_texture(1,neighbor_ip,neighbor_el)
+ !   neighbor_orientation = orientation(1:4, neighbor_n, neighbor_ip, neighbor_el)          !ipc is always 1.
+ !   absMisorientation    = lattice_qDisorientation(my_orientation, &
+ !                                                  neighbor_orientation, &
+ !                                                  0_pInt)                                 !no need for explicit calculation of symmetry
+
+ !   !***find the accumulative shear for this neighbor
+ !   LOOPFINDNEISHEAR: DO ne_slip_ac=1_pInt, ns
+ !    ne_acshear(ne_slip_ac) = plasticState(ph)%state(offset_acshear_slip+ne_slip_ac, &
+ !                                                    neighbor_of)
+ !   ENDDO LOOPFINDNEISHEAR
+
+ !   !***calculate the average accumulative shear and use it as cutoff
+ !   avg_acshear_ne = SUM(ne_acshear)/ns
+
+ !   !***
+ !   IF (ph==neighbor_ph) THEN
+ !    !***walk through all the
+ !    LOOPNEIGHBORSLIP: DO ne_slip=1_pInt,ns
+ !      !***only consider slip system that is active (above average accumulative shear)
+ !      IF (ne_acshear(ne_slip) > avg_acshear_ne) THEN
+ !        m_primes(ne_slip) = abs(math_mul3x3(slipNormal(1:3,me_slip), &
+ !                                math_qRot(absMisorientation, slipNormal(1:3,ne_slip)))) &
+ !                           *abs(math_mul3x3(slipDirect(1:3,me_slip), &
+ !                                math_qRot(absMisorientation, slipDirect(1:3,ne_slip))))
+ !        !***find the highest m' and corresponding accumulative shear
+ !        IF (m_primes(ne_slip) > tmp_mprime) THEN
+ !          tmp_mprime  = m_primes(ne_slip)
+ !          tmp_acshear = ne_acshear(ne_slip)
+ !        ENDIF
+ !      ENDIF
+ !    ENDDO LOOPNEIGHBORSLIP
+
+ !   ELSE
+ !    ne_mprimes(n)  = 0.0_pReal
+ !    vld_Nneighbors = vld_Nneighbors - 1_pInt
+ !   ENDIF
+
+ !  ENDDO LOOPNEIGHBORS
+
+ !  !***check if this element close to rim
+ !  IF (vld_Nneighbors < Nneighbors) THEN
+ !   !***rim voxel, no modification allowed
+ !   plasticState(ph)%state(index_kappa+me_slip, of) = 1.0_pReal
+ !  ELSE
+ !   !***patch voxel, started to calculate push up factor for gamma_dot
+ !   IF ((tmp_mprime > mprime_cut) .AND. (tmp_acshear > tmp_myshear_slip)) THEN
+ !    plasticState(ph)%state(index_kappa+me_slip, of) = 1.5_pReal
+ !   ELSE
+ !    !***minimum damping factor is 0.5
+ !    plasticState(ph)%state(index_kappa+me_slip, of) = 1.0_pReal
+ !   ENDIF
+ !  ENDIF
+
+ ! ENDDO LOOPMYSLIP
 
 end subroutine plastic_phenoplus_microstructure