let lattice to the work

src/lattice.f90

@@ -499,6 +499,31 @@ module lattice
    lattice_structure, trans_lattice_structure
 
  
+ interface lattice_forestProjection ! DEPRECATED, use lattice_forestProjection_edge
+   module procedure slipProjection_transverse
+ end interface lattice_forestProjection
+ 
+ interface lattice_forestProjection_edge
+   module procedure slipProjection_transverse
+ end interface lattice_forestProjection_edge
+ 
+ interface lattice_forestProjection_screw
+   module procedure slipProjection_direction
+ end interface lattice_forestProjection_screw
+ 
+ interface lattice_slipProjection_modeI
+   module procedure slipProjection_normal
+ end interface lattice_slipProjection_modeI
+ 
+  interface lattice_slipProjection_modeII
+   module procedure slipProjection_direction
+ end interface lattice_slipProjection_modeII
+ 
+  interface lattice_slipProjection_modeIII
+   module procedure slipProjection_transverse
+ end interface lattice_slipProjection_modeIII
+
+
  public :: &
   lattice_init, &
   lattice_qDisorientation, &
@@ -517,10 +542,16 @@ module lattice
   lattice_interaction_SlipTwin, &
   lattice_interaction_SlipTrans, &
   lattice_interaction_TwinSlip, &
-  lattice_forestProjection, &
   lattice_characteristicShear_Twin, &
   lattice_C66_twin, &
-  lattice_C66_trans
+  lattice_C66_trans, &
+  lattice_forestProjection, &
+  lattice_forestProjection_edge, &
+  lattice_forestProjection_screw, &
+  lattice_slipProjection_modeI, &
+  lattice_slipProjection_modeII, &
+  lattice_slipProjection_modeIII
+
 
 contains
 
@@ -2181,9 +2212,11 @@ end function lattice_SchmidMatrix_cleavage
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Forest projection (for edge dislocations)
+!> @brief Projection of the transverse direction onto the slip plane
+!> @details: This projection is used to calculate forest hardening for edge dislocations and for 
+! mode III failure (ToDo: MD I am not 100% sure about mode III)
 !--------------------------------------------------------------------------------------------------
-function lattice_forestProjection(Nslip,structure,cOverA) result(projection)
+function slipProjection_transverse(Nslip,structure,cOverA) result(projection)
  use math, only: &
    math_mul3x3
  use IO, only: &
@@ -2231,7 +2264,118 @@ function lattice_forestProjection(Nslip,structure,cOverA) result(projection)
    projection(i,j) = abs(math_mul3x3(coordinateSystem(1:3,2,i),coordinateSystem(1:3,3,j)))
  enddo; enddo
 
-end function lattice_forestProjection
+end function slipProjection_transverse
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Projection of the slip direction onto the slip plane
+!> @details: This projection is used to calculate forest hardening for screw dislocations and for 
+! mode II failure (ToDo: MD I am not 100% sure about mode II)
+!--------------------------------------------------------------------------------------------------
+function slipProjection_direction(Nslip,structure,cOverA) result(projection)
+ use math, only: &
+   math_mul3x3
+ use IO, only: &
+   IO_error
+
+ implicit none
+ integer(pInt),    dimension(:),            intent(in) :: Nslip                                     !< number of active slip systems per family
+ character(len=*),                          intent(in) :: structure                                 !< lattice structure
+ real(pReal),                               intent(in) :: cOverA                                    !< c/a ratio
+ real(pReal),     dimension(sum(Nslip),sum(Nslip))     :: projection
+
+ real(pReal),     dimension(3,3,sum(Nslip))            :: coordinateSystem
+ real(pReal),     dimension(:,:), allocatable          :: slipSystems
+ integer(pInt),   dimension(:), allocatable            :: NslipMax
+ integer(pInt) :: i, j
+
+ if (len_trim(structure) /= 3_pInt) &
+   call IO_error(137_pInt,ext_msg='lattice_forestProjection: '//trim(structure))
+
+ select case(structure(1:3))
+   case('fcc')
+     NslipMax    = LATTICE_FCC_NSLIPSYSTEM
+     slipSystems = LATTICE_FCC_SYSTEMSLIP
+   case('bcc')
+     NslipMax    = LATTICE_BCC_NSLIPSYSTEM
+     slipSystems = LATTICE_BCC_SYSTEMSLIP
+   case('hex')
+     NslipMax    = LATTICE_HEX_NSLIPSYSTEM
+     slipSystems = LATTICE_HEX_SYSTEMSLIP
+   case('bct')
+     NslipMax    = LATTICE_BCT_NSLIPSYSTEM
+     slipSystems = LATTICE_BCT_SYSTEMSLIP
+   case default
+     call IO_error(137_pInt,ext_msg='lattice_forestProjection: '//trim(structure))
+ end select
+
+ if (any(NslipMax(1:size(Nslip)) - Nslip < 0_pInt)) &
+   call IO_error(145_pInt,ext_msg='Nslip '//trim(structure))
+ if (any(Nslip < 0_pInt)) &
+   call IO_error(144_pInt,ext_msg='Nslip '//trim(structure))
+
+ coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure,cOverA)
+
+ do i=1_pInt, sum(Nslip); do j=1_pInt, sum(Nslip)
+   projection(i,j) = abs(math_mul3x3(coordinateSystem(1:3,2,i),coordinateSystem(1:3,1,j)))
+ enddo; enddo
+
+end function slipProjection_direction
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Projection of the slip plane onto itself
+!> @details: This projection is used for mode I failure
+!--------------------------------------------------------------------------------------------------
+function slipProjection_normal(Nslip,structure,cOverA) result(projection)
+ use math, only: &
+   math_mul3x3
+ use IO, only: &
+   IO_error
+
+ implicit none
+ integer(pInt),    dimension(:),            intent(in) :: Nslip                                     !< number of active slip systems per family
+ character(len=*),                          intent(in) :: structure                                 !< lattice structure
+ real(pReal),                               intent(in) :: cOverA                                    !< c/a ratio
+ real(pReal),     dimension(sum(Nslip),sum(Nslip))     :: projection
+
+ real(pReal),     dimension(3,3,sum(Nslip))            :: coordinateSystem
+ real(pReal),     dimension(:,:), allocatable          :: slipSystems
+ integer(pInt),   dimension(:), allocatable            :: NslipMax
+ integer(pInt) :: i, j
+
+ if (len_trim(structure) /= 3_pInt) &
+   call IO_error(137_pInt,ext_msg='lattice_forestProjection: '//trim(structure))
+
+ select case(structure(1:3))
+   case('fcc')
+     NslipMax    = LATTICE_FCC_NSLIPSYSTEM
+     slipSystems = LATTICE_FCC_SYSTEMSLIP
+   case('bcc')
+     NslipMax    = LATTICE_BCC_NSLIPSYSTEM
+     slipSystems = LATTICE_BCC_SYSTEMSLIP
+   case('hex')
+     NslipMax    = LATTICE_HEX_NSLIPSYSTEM
+     slipSystems = LATTICE_HEX_SYSTEMSLIP
+   case('bct')
+     NslipMax    = LATTICE_BCT_NSLIPSYSTEM
+     slipSystems = LATTICE_BCT_SYSTEMSLIP
+   case default
+     call IO_error(137_pInt,ext_msg='lattice_forestProjection: '//trim(structure))
+ end select
+
+ if (any(NslipMax(1:size(Nslip)) - Nslip < 0_pInt)) &
+   call IO_error(145_pInt,ext_msg='Nslip '//trim(structure))
+ if (any(Nslip < 0_pInt)) &
+   call IO_error(144_pInt,ext_msg='Nslip '//trim(structure))
+
+ coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure,cOverA)
+
+ do i=1_pInt, sum(Nslip); do j=1_pInt, sum(Nslip)
+   projection(i,j) = abs(math_mul3x3(coordinateSystem(1:3,2,i),coordinateSystem(1:3,2,j)))
+ enddo; enddo
+
+end function slipProjection_normal
 
 
 !--------------------------------------------------------------------------------------------------

src/plastic_nonlocal.f90

@@ -64,9 +64,7 @@ module plastic_nonlocal
    minDipoleHeightPerSlipFamily, &                                                                  !< minimum stable edge/screw dipole height for each family and instance
    minDipoleHeight, &                                                                               !< minimum stable edge/screw dipole height for each slip system and instance
    peierlsStressPerSlipFamily, &                                                                    !< Peierls stress (edge and screw) 
-   peierlsStress, &                                                                                 !< Peierls stress (edge and screw) 
+   peierlsStress                                                                                 !< Peierls stress (edge and screw) 
-   forestProjectionEdge, &                                                                          !< matrix of forest projections of edge dislocations for each instance
-   forestProjectionScrew                                                                          !< matrix of forest projections of screw dislocations for each instance
 
  real(pReal), dimension(:,:,:,:), allocatable, private :: &
    rhoDotEdgeJogsOutput, &
@@ -172,8 +170,8 @@ module plastic_nonlocal
    
    real(pReal), dimension(:,:), allocatable     :: &
    interactionSlipSlip ,&                                                                              !< coefficients for slip-slip interaction for each interaction type and instance
-   forestProjectionEdge, &                                                                          !< matrix of forest projections of edge dislocations for each instance
+   forestProjection_Edge, &                                                                          !< matrix of forest projections of edge dislocations for each instance
-   forestProjectionScrew                                                                          !< matrix of forest projections of screw dislocations for each instance
+   forestProjection_Screw                                                                          !< matrix of forest projections of screw dislocations for each instance
      real(pReal), dimension(:), allocatable, private :: &
    nonSchmidCoeff
    integer(pInt) :: totalNslip
@@ -531,8 +529,6 @@ allocate(iTauB(maxTotalNslip,maxNinstances),   source=0_pInt)
 
 allocate(lambda0(maxTotalNslip,maxNinstances),                                        source=0.0_pReal)
 allocate(minDipoleHeight(maxTotalNslip,2,maxNinstances),                              source=-1.0_pReal)
-allocate(forestProjectionEdge(maxTotalNslip,maxTotalNslip,maxNinstances),             source=0.0_pReal)
-allocate(forestProjectionScrew(maxTotalNslip,maxTotalNslip,maxNinstances),            source=0.0_pReal)
 allocate(sourceProbability(maxTotalNslip,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%nElems), &
                                                                                    source=2.0_pReal)
 
@@ -681,17 +677,6 @@ allocate(colinearSystem(maxTotalNslip,maxNinstances),
    
        do s2 = 1_pInt,ns
 
-         !*** calculation of forest projections for edge and screw dislocations. s2 acts as forest for s1
-   
-         forestProjectionEdge(s1,s2,instance) &
-             = abs(math_mul3x3(lattice_sn(1:3,slipSystemLattice(s1,instance),phase), &
-                               lattice_st(1:3,slipSystemLattice(s2,instance),phase)))                   ! forest projection of edge dislocations is the projection of (t = b x n) onto the slip normal of the respective slip plane
-         
-         forestProjectionScrew(s1,s2,instance) &
-             = abs(math_mul3x3(lattice_sn(1:3,slipSystemLattice(s1,instance),phase), &
-                               lattice_sd(1:3,slipSystemLattice(s2,instance),phase)))                   ! forest projection of screw dislocations is the projection of b onto the slip normal of the respective splip plane
-     
-         
          !*** colinear slip system (only makes sense for fcc like it is defined here)
          
          if ((all(dEq(lattice_sd(1:3,slipSystemLattice(s1,instance),phase), &
@@ -769,8 +754,10 @@ param(instance)%probabilisticMultiplication = .false.
    
 
    prm%burgers = math_expand(prm%burgers,prm%Nslip)
-   
+        prm%forestProjection_edge     = lattice_forestProjection_edge (prm%Nslip,config%getString('lattice_structure'),&
-  
+                                                          config%getFloat('c/a',defaultVal=0.0_pReal))
+          prm%forestProjection_screw     = lattice_forestProjection_screw (prm%Nslip,config%getString('lattice_structure'),&
+                                                          config%getFloat('c/a',defaultVal=0.0_pReal))
    minDipoleHeightPerSlipFamily(:,1_pInt,instance) = config_phase(p)%getFloats('minimumdipoleheightedge')!,'ddipminedge')
    minDipoleHeightPerSlipFamily(:,2_pInt,instance) = config_phase(p)%getFloats('minimumdipoleheightscrew')!,'ddipminscrew')
    peierlsStressPerSlipFamily(:,1_pInt,instance) = config_phase(p)%getFloat('peierlsstressedge')!,'peierlsstress_edge')
@@ -1293,9 +1280,9 @@ where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < prm%significantN &
 
 forall (s = 1_pInt:ns) &
   rhoForest(s) = dot_product((sum(abs(rhoSgl(1:ns,[1,2,5,6])),2) + rhoDip(1:ns,1)), &
-                              forestProjectionEdge(s,1:ns,instance)) & 
+                              prm%forestProjection_Edge(s,1:ns)) & 
                + dot_product((sum(abs(rhoSgl(1:ns,[3,4,7,8])),2) + rhoDip(1:ns,2)), &
-                              forestProjectionScrew(s,1:ns,instance))
+                              prm%forestProjection_Screw(s,1:ns))
 
 
 !*** calculate the threshold shear stress for dislocation slip