Merge remote-tracking branch 'origin/20-NewStyleDislotwin' into 20-NewStyleDislotwin

VERSION

@@ -1 +1 @@
-v2.0.2-514-gbfa56e9b
+v2.0.2-540-gce6e6679

src/math.f90

@@ -12,7 +12,7 @@ module math
 
  implicit none
  private
- real(pReal),    parameter, public :: PI = 3.141592653589793_pReal                                  !< ratio of a circle's circumference to its diameter
+ real(pReal),    parameter, public :: PI = acos(-1.0_pReal)                                         !< ratio of a circle's circumference to its diameter
  real(pReal),    parameter, public :: INDEG = 180.0_pReal/PI                                        !< conversion from radian into degree
  real(pReal),    parameter, public :: INRAD = PI/180.0_pReal                                        !< conversion from degree into radian
  complex(pReal), parameter, public :: TWOPIIMG = (0.0_pReal,2.0_pReal)*(PI,0.0_pReal)               !< Re(0.0), Im(2xPi)
@@ -36,13 +36,13 @@ module math
 
  real(pReal), dimension(6), parameter, private :: &
    nrmMandel = [&
-     1.0_pReal,                1.0_pReal,                1.0_pReal,&
+     1.0_pReal,       1.0_pReal,       1.0_pReal, &
-     1.414213562373095_pReal,  1.414213562373095_pReal,  1.414213562373095_pReal ]                  !< weighting for Mandel notation (forward)
+     sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal) ]                                            !< weighting for Mandel notation (forward)
 
  real(pReal), dimension(6), parameter , public :: &
    invnrmMandel = [&
-     1.0_pReal,                1.0_pReal,                1.0_pReal,&
+     1.0_pReal,                 1.0_pReal,                 1.0_pReal, &
-     0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal ]                 !< weighting for Mandel notation (backward)
+     1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal) ]              !< weighting for Mandel notation (backward)
 
  integer(pInt), dimension (2,6), parameter, private :: &
    mapVoigt = reshape([&
@@ -160,7 +160,7 @@ module math
    math_rotate_forward33, &
    math_rotate_backward33, &
    math_rotate_forward3333, &
-   math_limit
+   math_clip
  private :: &
    math_check, &
    halton
@@ -1366,16 +1366,16 @@ pure function math_RtoEuler(R)
  sqhk =sqrt(R(1,3)*R(1,3)+R(2,3)*R(2,3))
 
 ! calculate PHI
- math_RtoEuler(2) = acos(math_limit(R(3,3)/sqhkl,-1.0_pReal, 1.0_pReal))
+ math_RtoEuler(2) = acos(math_clip(R(3,3)/sqhkl,-1.0_pReal, 1.0_pReal))
 
  if((math_RtoEuler(2) < 1.0e-8_pReal) .or. (pi-math_RtoEuler(2) < 1.0e-8_pReal)) then
    math_RtoEuler(3) = 0.0_pReal
-   math_RtoEuler(1) = acos(math_limit(R(1,1)/squvw, -1.0_pReal, 1.0_pReal))
+   math_RtoEuler(1) = acos(math_clip(R(1,1)/squvw, -1.0_pReal, 1.0_pReal))
    if(R(2,1) > 0.0_pReal) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1)
  else
-   math_RtoEuler(3) = acos(math_limit(R(2,3)/sqhk, -1.0_pReal, 1.0_pReal))
+   math_RtoEuler(3) = acos(math_clip(R(2,3)/sqhk, -1.0_pReal, 1.0_pReal))
    if(R(1,3) < 0.0) math_RtoEuler(3) = 2.0_pReal*pi-math_RtoEuler(3)
-   math_RtoEuler(1) = acos(math_limit(-R(3,2)/sin(math_RtoEuler(2)), -1.0_pReal, 1.0_pReal))
+   math_RtoEuler(1) = acos(math_clip(-R(3,2)/sin(math_RtoEuler(2)), -1.0_pReal, 1.0_pReal))
    if(R(3,1) < 0.0) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1)
  end if
 
@@ -1657,7 +1657,7 @@ pure function math_qToEuler(qPassive)
  math_qToEuler(2) = acos(1.0_pReal-2.0_pReal*(q(2)**2+q(3)**2))
 
  if (abs(math_qToEuler(2)) < 1.0e-6_pReal) then
-   math_qToEuler(1) = sign(2.0_pReal*acos(math_limit(q(1),-1.0_pReal, 1.0_pReal)),q(4))
+   math_qToEuler(1) = sign(2.0_pReal*acos(math_clip(q(1),-1.0_pReal, 1.0_pReal)),q(4))
    math_qToEuler(3) = 0.0_pReal
  else
    math_qToEuler(1) = atan2(+q(1)*q(3)+q(2)*q(4), q(1)*q(2)-q(3)*q(4))
@@ -1684,7 +1684,7 @@ pure function math_qToAxisAngle(Q)
  real(pReal) :: halfAngle, sinHalfAngle
  real(pReal), dimension(4) :: math_qToAxisAngle
 
- halfAngle = acos(math_limit(Q(1),-1.0_pReal,1.0_pReal))
+ halfAngle = acos(math_clip(Q(1),-1.0_pReal,1.0_pReal))
  sinHalfAngle = sin(halfAngle)
 
  smallRotation: if (sinHalfAngle <= 1.0e-4_pReal) then
@@ -1744,7 +1744,7 @@ real(pReal) pure function math_EulerMisorientation(EulerA,EulerB)
  cosTheta = (math_trace33(math_mul33x33(math_EulerToR(EulerB), &
                               transpose(math_EulerToR(EulerA)))) - 1.0_pReal) * 0.5_pReal
 
- math_EulerMisorientation = acos(math_limit(cosTheta,-1.0_pReal,1.0_pReal))
+ math_EulerMisorientation = acos(math_clip(cosTheta,-1.0_pReal,1.0_pReal))
 
 end function math_EulerMisorientation
 
@@ -2055,7 +2055,7 @@ function math_eigenvectorBasisSym33(m)
    EB(3,3,3)=1.0_pReal
  else threeSimilarEigenvalues
    rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
-   phi=acos(math_limit(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
+   phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
    values = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* &
                              [cos(phi/3.0_pReal), &
                               cos((phi+2.0_pReal*PI)/3.0_pReal), &
@@ -2120,7 +2120,7 @@ function math_eigenvectorBasisSym33_log(m)
    EB(3,3,3)=1.0_pReal
  else threeSimilarEigenvalues
    rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
-   phi=acos(math_limit(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
+   phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
    values = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* &
                              [cos(phi/3.0_pReal), &
                               cos((phi+2.0_pReal*PI)/3.0_pReal), &
@@ -2232,7 +2232,7 @@ function math_eigenvaluesSym33(m)
    math_eigenvaluesSym33 = math_eigenvaluesSym(m)
  else
    rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
-   phi=acos(math_limit(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
+   phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
    math_eigenvaluesSym33 = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* &
                          [cos(phi/3.0_pReal), &
                           cos((phi+2.0_pReal*PI)/3.0_pReal), &
@@ -2617,7 +2617,7 @@ end function math_rotate_forward3333
 !> @brief limits a scalar value to a certain range (either one or two sided)
 ! Will return NaN if left > right
 !--------------------------------------------------------------------------------------------------
-real(pReal) pure function math_limit(a, left, right)
+real(pReal) pure function math_clip(a, left, right)
  use, intrinsic :: &
    IEEE_arithmetic
 
@@ -2626,14 +2626,14 @@ real(pReal) pure function math_limit(a, left, right)
  real(pReal), intent(in), optional :: left, right
 
    
- math_limit = min ( &
+ math_clip = min ( &
                    max (merge(left, -huge(a), present(left)), a), &
                         merge(right, huge(a), present(right)) &
                   )
 
  if (present(left) .and. present(right)) &
-   math_limit = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_limit, left>right)
+   math_clip = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_clip, left>right)
 
-end function math_limit
+end function math_clip
 
 end module math

src/plastic_dislotwin.f90

@@ -448,50 +448,50 @@ subroutine plastic_dislotwin_init(fileUnit)
      outputID = undefined_ID
      select case(outputs(i))
        case ('edge_density')
-         outputID = edge_density_ID
+         outputID = merge(edge_density_ID,undefined_ID,prm%totalNslip > 0_pInt)
          outputSize = prm%totalNslip
        case ('dipole_density')
-         outputID = dipole_density_ID
+         outputID = merge(dipole_density_ID,undefined_ID,prm%totalNslip > 0_pInt)
          outputSize = prm%totalNslip
        case ('shear_rate_slip','shearrate_slip')
-         outputID = shear_rate_slip_ID
+         outputID = merge(shear_rate_slip_ID,undefined_ID,prm%totalNslip > 0_pInt)
          outputSize = prm%totalNslip
        case ('accumulated_shear_slip')
-         outputID = accumulated_shear_slip_ID
+         outputID = merge(accumulated_shear_slip_ID,undefined_ID,prm%totalNslip > 0_pInt)
          outputSize = prm%totalNslip
        case ('mfp_slip')
-         outputID = mfp_slip_ID
+         outputID = merge(mfp_slip_ID,undefined_ID,prm%totalNslip > 0_pInt)
          outputSize = prm%totalNslip
        case ('resolved_stress_slip')
-         outputID = resolved_stress_slip_ID
+         outputID = merge(resolved_stress_slip_ID,undefined_ID,prm%totalNslip > 0_pInt)
          outputSize = prm%totalNslip
        case ('threshold_stress_slip')
-         outputID= threshold_stress_slip_ID
+         outputID= merge(threshold_stress_slip_ID,undefined_ID,prm%totalNslip > 0_pInt)
          outputSize = prm%totalNslip
        case ('edge_dipole_distance')
-         outputID = edge_dipole_distance_ID
+         outputID = merge(edge_dipole_distance_ID,undefined_ID,prm%totalNslip > 0_pInt)
          outputSize = prm%totalNslip
        case ('stress_exponent')
-         outputID = stress_exponent_ID
+         outputID = merge(stress_exponent_ID,undefined_ID,prm%totalNslip > 0_pInt)
          outputSize = prm%totalNslip
 
        case ('twin_fraction')
-         outputID = twin_fraction_ID
+         outputID = merge(twin_fraction_ID,undefined_ID,prm%totalNtwin >0_pInt)
          outputSize = prm%totalNtwin
        case ('shear_rate_twin','shearrate_twin')
-         outputID = shear_rate_twin_ID
+         outputID = merge(shear_rate_twin_ID,undefined_ID,prm%totalNtwin >0_pInt)
          outputSize = prm%totalNtwin
        case ('accumulated_shear_twin')
-         outputID = accumulated_shear_twin_ID
+         outputID = merge(accumulated_shear_twin_ID,undefined_ID,prm%totalNtwin >0_pInt)
          outputSize = prm%totalNtwin
        case ('mfp_twin')
-         outputID = mfp_twin_ID
+         outputID = merge(mfp_twin_ID,undefined_ID,prm%totalNtwin >0_pInt)
          outputSize = prm%totalNtwin
        case ('resolved_stress_twin')
-         outputID = resolved_stress_twin_ID
+         outputID = merge(resolved_stress_twin_ID,undefined_ID,prm%totalNtwin >0_pInt)
          outputSize = prm%totalNtwin
        case ('threshold_stress_twin')
-         outputID = threshold_stress_twin_ID
+         outputID = merge(threshold_stress_twin_ID,undefined_ID,prm%totalNtwin >0_pInt)
          outputSize = prm%totalNtwin
          
        case ('resolved_stress_shearband')
@@ -840,54 +840,28 @@ subroutine plastic_dislotwin_init(fileUnit)
    dst%strainTransFraction=>plasticState(p)%dotState(startIndex:endIndex,:)
    plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolTransFrac
 
-   invLambdaSlip0 = spread(0.0_pReal,1,prm%totalNslip)                                   ! calculation required? Seems to be the same as in microstructure
+   dst%whole => plasticState(p)%dotState
-   forall (i = 1_pInt:prm%totalNslip) &
-     invLambdaSlip0(i) = sqrt(dot_product(math_expand(prm%rho0,prm%Nslip)+ &
-     math_expand(prm%rhoDip0,prm%Nslip),prm%forestProjectionEdge(1:prm%totalNslip,i)))/ &
-                         prm%CLambdaSlip(i)
-   mse%invLambdaSlip = spread(math_expand(invLambdaSlip0,prm%Nslip),2, NofMyPhase)
-
 
+   allocate(mse%invLambdaSlip(prm%totalNslip,NofMyPhase),source=0.0_pReal)
    allocate(mse%invLambdaSlipTwin(prm%totalNslip,NofMyPhase),source=0.0_pReal)
    allocate(mse%invLambdaTwin(prm%totalNtwin,NofMyPhase),source=0.0_pReal)
    allocate(mse%invLambdaSlipTrans(prm%totalNtrans,NofMyPhase),source=0.0_pReal)
    allocate(mse%invLambdaTrans(prm%totalNtrans,NofMyPhase),source=0.0_pReal)
 
-   MeanFreePathSlip0  = prm%GrainSize/(1.0_pReal+invLambdaSlip0*prm%GrainSize)
+   allocate(mse%mfp_slip(prm%totalNslip,NofMyPhase),  source=0.0_pReal)
-   mse%mfp_slip       = spread(math_expand(MeanFreePathSlip0,prm%Nslip),2, NofMyPhase)
+   allocate(mse%mfp_twin(prm%totalNtwin,NofMyPhase),  source=0.0_pReal)
-   
+   allocate(mse%mfp_trans(prm%totalNtrans,NofMyPhase),source=0.0_pReal)
-   MeanFreePathTwin0  = spread(prm%GrainSize,1,prm%totalNtwin)
-   mse%mfp_twin       = spread(math_expand(MeanFreePathTwin0,prm%Ntwin),2, NofMyPhase)
-
-   MeanFreePathTrans0 = spread(prm%GrainSize,1,prm%totalNtrans)
-   mse%mfp_trans      = spread(math_expand(MeanFreePathTrans0,prm%Ntrans),2, NofMyPhase)
-
-   tauSlipThreshold0  = spread(0.0_pReal,1,prm%totalNslip)
-   forall (i = 1_pInt:prm%totalNslip) tauSlipThreshold0(i) = &
-     prm%mu*prm%burgers_slip(i) * sqrt(dot_product(math_expand(prm%rho0 + prm%rhoDip0,prm%Nslip),&
-                                                   prm%interaction_SlipSlip(i,1:prm%totalNslip)))
-   mse%threshold_stress_slip = spread(math_expand(tauSlipThreshold0,prm%Nslip),2, NofMyPhase)
 
+   allocate(mse%threshold_stress_slip(prm%totalNslip,NofMyPhase),  source=0.0_pReal)
    allocate(mse%threshold_stress_twin(prm%totalNtwin,NofMyPhase),  source=0.0_pReal)
    allocate(mse%threshold_stress_trans(prm%totalNtrans,NofMyPhase),source=0.0_pReal)
 
-   TwinVolume0= spread(0.0_pReal,1,prm%totalNtwin)
-   forall (i = 1_pInt:prm%totalNtwin) TwinVolume0(i) = &
-     (PI/4.0_pReal)*prm%twinsize(i)*MeanFreePathTwin0(i)**2.0_pReal
-   mse%twinVolume = &
-      spread(math_expand(TwinVolume0,prm%Ntwin),2, NofMyPhase)
-
-   MartensiteVolume0= spread(0.0_pReal,1,prm%totalNtrans)
-   forall (i = 1_pInt:prm%totalNtrans) MartensiteVolume0(i) = &
-     (PI/4.0_pReal)*prm%lamellarsizePerTransSystem(i)*MeanFreePathTrans0(i)**2.0_pReal
-   mse%martensiteVolume = &
-     spread(math_expand(MartensiteVolume0,prm%Ntrans),2, NofMyPhase)
-
-   dst%whole => plasticState(p)%dotState
-
    allocate(mse%tau_r_twin(prm%totalNtwin,NofMyPhase),   source=0.0_pReal)
    allocate(mse%tau_r_trans(prm%totalNtrans,NofMyPhase), source=0.0_pReal)
 
+   allocate(mse%twinVolume(prm%totalNtwin,NofMyPhase),   source=0.0_pReal)
+   allocate(mse%martensiteVolume(prm%totalNtrans,NofMyPhase), source=0.0_pReal)
+
    end associate
  enddo
  
@@ -916,7 +890,6 @@ function plastic_dislotwin_homogenizedC(ipc,ip,el)
                   of
  real(pReal) :: sumf_twin, sumf_trans
 
- !* Shortened notation
  of = phasememberAt(ipc,ip,el)
  associate(prm => param(phase_plasticityInstance(material_phase(ipc,ip,el))),&
            stt => state(phase_plasticityInstance(material_phase(ipc,ip,el))))