centralized function for ones/zeros

src/misc.f90

@@ -21,7 +21,9 @@ module misc
     misc_init, &
     misc_selfTest, &
     misc_optional, &
-    misc_prefixOptions
+    misc_prefixOptions, &
+    misc_ones, &
+    misc_zeros
 
 contains
 
@@ -136,6 +138,34 @@ pure function misc_prefixOptions(string,prefix) result(prefixed)
 end function misc_prefixOptions
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief 1D array of zeros.
+!--------------------------------------------------------------------------------------------------
+pure function misc_zeros(N)
+
+  integer, intent(in) :: N                                                                          !< number of ones
+  real(pREAL), dimension(N) :: misc_zeros
+
+
+  misc_zeros = 0._pREAL
+
+end function misc_zeros
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief 1D array of ones.
+!--------------------------------------------------------------------------------------------------
+pure function misc_ones(N)
+
+  integer, intent(in) :: N                                                                          !< number of ones
+  real(pREAL), dimension(N) :: misc_ones
+
+
+  misc_ones = 1._pREAL
+
+end function misc_ones
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief Check correctness of some misc functions.
 !--------------------------------------------------------------------------------------------------
@@ -143,7 +173,7 @@ subroutine misc_selfTest()
 
   real(pREAL) :: r
   character(len=:),      allocatable :: str,out
-
+  integer :: N
 
   call random_number(r)
   if (test_str('DAMASK') /= 'DAMASK')                        error stop 'optional_str, present'
@@ -163,6 +193,12 @@ subroutine misc_selfTest()
   out=misc_prefixOptions(str,'p_')
   if (out /= '-p_a -1 -p_more 123 -p_flag -')                error stop 'misc_prefixOptions'
 
+  N = int(r*99._pReal)
+  if (size(misc_zeros(N)) /= N)           error stop 'shape zeros'
+  if (size(misc_ones(N)) /= N)            error stop 'shape ones'
+  if (any(dNeq(misc_zeros(N),0.0_pReal))) error stop 'value zeros'
+  if (any(dNeq(misc_ones(N),1.0_pReal)))  error stop 'value ones'
+
 contains
 
   function test_str(str_in) result(str_out)

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -85,8 +85,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
     N_tw                                                                                            !< number of twin-systems for a given twin family
   real(pREAL), dimension(:), allocatable :: &
     xi_0_sl, &                                                                                      !< initial critical shear stress for slip
-    xi_0_tw, &                                                                                      !< initial critical shear stress for twin
+    xi_0_tw                                                                                          !< initial critical shear stress for twin
-    ones
   real(pREAL), dimension(:,:), allocatable :: &
     a_nS                                                                                            !< non-Schmid coefficients
   character(len=:), allocatable :: &
@@ -141,15 +140,18 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
 !--------------------------------------------------------------------------------------------------
 ! slip related parameters
     slipActive: if (prm%sum_N_sl > 0) then
-      prm%dot_gamma_0_sl = math_expand(pl%get_as1dReal('dot_gamma_0_sl',requiredSize=size(N_sl)                          ), N_sl)
+      prm%dot_gamma_0_sl = math_expand(pl%get_as1dReal('dot_gamma_0_sl',requiredSize=size(N_sl)), N_sl)
-      prm%n_sl           = math_expand(pl%get_as1dReal('n_sl',          requiredSize=size(N_sl)                          ), N_sl)
+      prm%n_sl           = math_expand(pl%get_as1dReal('n_sl',          requiredSize=size(N_sl)), N_sl)
-      prm%a_sl           = math_expand(pl%get_as1dReal('a_sl',          requiredSize=size(N_sl)                          ), N_sl)
+      prm%a_sl           = math_expand(pl%get_as1dReal('a_sl',          requiredSize=size(N_sl)), N_sl)
-      prm%h_0_sl_sl      = math_expand(pl%get_as1dReal('h_0_sl-sl',     requiredSize=size(N_sl)                          ), N_sl)
+      prm%h_0_sl_sl      = math_expand(pl%get_as1dReal('h_0_sl-sl',     requiredSize=size(N_sl)), N_sl)
-      prm%c_1            = math_expand(pl%get_as1dReal('c_1',           requiredSize=size(N_sl),defaultVal=0.0_pREAL*ones), N_sl)
-      prm%c_2            = math_expand(pl%get_as1dReal('c_2',           requiredSize=size(N_sl),defaultVal=          ones), N_sl)
       xi_0_sl            = math_expand(pl%get_as1dReal('xi_0_sl',       requiredSize=size(N_sl)), N_sl)
       prm%xi_inf_sl      = math_expand(pl%get_as1dReal('xi_inf_sl',     requiredSize=size(N_sl)), N_sl)
-      prm%f_sat_sl_tw    = math_expand(pl%get_as1dReal('f_sat_sl-tw',   requiredSize=size(N_sl),defaultVal=0.0_pREAL*ones), N_sl)
+      prm%c_1            = math_expand(pl%get_as1dReal('c_1',           requiredSize=size(N_sl), &
+                                                                        defaultVal=misc_ones(size(N_sl))), N_sl)
+      prm%c_2            = math_expand(pl%get_as1dReal('c_2',           requiredSize=size(N_sl), &
+                                                                        defaultVal=misc_ones(size(N_sl))), N_sl)
+      prm%f_sat_sl_tw    = math_expand(pl%get_as1dReal('f_sat_sl-tw',   requiredSize=size(N_sl), &
+                                                                        defaultVal=misc_zeros(size(N_sl))), N_sl)
 
       prm%h_sl_sl = crystal_interaction_SlipBySlip(N_sl,pl%get_as1dReal('h_sl-sl'),phase_lattice(ph))
 
@@ -191,13 +193,14 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
 !--------------------------------------------------------------------------------------------------
 ! twin related parameters
     twinActive: if (prm%sum_N_tw > 0) then
-      ones = [(1.0_pREAL,i=1,size(N_tw))]
+      prm%dot_gamma_0_tw = math_expand(pl%get_as1dReal('dot_gamma_0_tw', requiredSize=size(N_tw)), N_tw)
-      prm%dot_gamma_0_tw = math_expand(pl%get_as1dReal('dot_gamma_0_tw', requiredSize=size(N_tw)                          ), N_tw)
+      prm%n_tw           = math_expand(pl%get_as1dReal('n_tw',           requiredSize=size(N_tw)), N_tw)
-      prm%n_tw           = math_expand(pl%get_as1dReal('n_tw',           requiredSize=size(N_tw)                          ), N_tw)
+      prm%h_0_tw_tw      = math_expand(pl%get_as1dReal('h_0_tw-tw',      requiredSize=size(N_tw)), N_tw)
-      prm%c_3            = math_expand(pl%get_as1dReal('c_3',            requiredSize=size(N_tw),defaultVal=0.0_pREAL*ones), N_tw)
-      prm%c_4            = math_expand(pl%get_as1dReal('c_4',            requiredSize=size(N_tw),defaultVal=0.0_pREAL*ones), N_tw)
-      prm%h_0_tw_tw      = math_expand(pl%get_as1dReal('h_0_tw-tw',      requiredSize=size(N_tw)                          ), N_tw)
       xi_0_tw            = math_expand(pl%get_as1dReal('xi_0_tw',        requiredSize=size(N_tw)), N_tw)
+      prm%c_3            = math_expand(pl%get_as1dReal('c_3',            requiredSize=size(N_tw), &
+                                                                         defaultVal=misc_ones(size(N_tw))), N_tw)
+      prm%c_4            = math_expand(pl%get_as1dReal('c_4',            requiredSize=size(N_tw), &
+                                                                         defaultVal=misc_zeros(size(N_tw))), N_tw)
 
       prm%gamma_char = crystal_characteristicShear_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
       prm%h_tw_tw    = crystal_interaction_TwinByTwin(N_tw,pl%get_as1dReal('h_tw-tw'),phase_lattice(ph))