Merge branch 'development' into even-more-HDF5-postprocessing

PRIVATE

@@ -1 +1 @@
-Subproject commit e8c935312c9ff22873f792282ff13e5a991752d3
+Subproject commit 5c5adbd8ccc0210fd6507431db8ec82ecec75352

VERSION

@@ -1 +1 @@
-v2.0.3-618-g4340c558
+v2.0.3-625-gbace77db

processing/post/addOrientations.py

@@ -27,6 +27,7 @@ Additional (globally fixed) rotations of the lab frame and/or crystal frame can
 representations = {
                   'quaternion': ['qu',4],
                   'rodrigues':  ['ro',4],
+                  'Rodrigues':  ['Ro',3],
                   'eulers':     ['eu',3],
                   'matrix':     ['om',9],
                   'angleaxis':  ['ax',4],
@@ -80,15 +81,20 @@ parser.add_option('-z',
                   dest = 'z',
                   metavar = 'string',
                   help = 'label of lab z vector (expressed in crystal coords)')
+parser.add_option('--lattice',
+                  dest = 'lattice',
+                  metavar = 'string',
+                  help = 'lattice structure to reduce rotation into fundamental zone')
 
 parser.set_defaults(output = [],
                     labrotation     = (1.,1.,1.,0.),                                                # no rotation about (1,1,1)
                     crystalrotation = (1.,1.,1.,0.),                                                # no rotation about (1,1,1)
+                    lattice = None,
                    )
 
 (options, filenames) = parser.parse_args()
 
-options.output = list(map(lambda x: x.lower(), options.output))
+#options.output = list(map(lambda x: x.lower(), options.output))
 if options.output == [] or (not set(options.output).issubset(set(representations))):
   parser.error('output must be chosen from {}.'.format(', '.join(representations)))
 
@@ -121,7 +127,7 @@ if filenames == []: filenames = [None]
 for name in filenames:
   try:    table = damask.ASCIItable(name = name,
                                     buffered = False)
-  except: continue
+  except Exception: continue
   damask.util.report(scriptName,name)
 
 # ------------------------------------------ read header ------------------------------------------
@@ -156,16 +162,16 @@ for name in filenames:
   outputAlive = True
   while outputAlive and table.data_read():                                                          # read next data line of ASCII table
     if   inputtype == 'eulers':
-      l = representations['eulers'][1]
+      d = representations['eulers'][1]
-      o = damask.Rotation.fromEulers(list(map(float,table.data[column:column+l])),options.degrees)
+      o = damask.Rotation.fromEulers(list(map(float,table.data[column:column+d])),options.degrees)
       
     elif inputtype == 'rodrigues':
-      l = representations['rodrigues'][1]
+      d = representations['rodrigues'][1]
-      o = damask.Rotation.fromRodrigues(list(map(float,table.data[column:column+l])))
+      o = damask.Rotation.fromRodrigues(list(map(float,table.data[column:column+d])))
       
     elif inputtype == 'matrix':
-      l = representations['matrix'][1]
+      d = representations['matrix'][1]
-      o = damask.Rotation.fromMatrix(list(map(float,table.data[column:column+l])))
+      o = damask.Rotation.fromMatrix(list(map(float,table.data[column:column+d])))
 
     elif inputtype == 'frame':
       M = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
@@ -174,14 +180,18 @@ for name in filenames:
       o = damask.Rotation.fromMatrix(M/np.linalg.norm(M,axis=0))
       
     elif inputtype == 'quaternion':
-      l = representations['quaternion'][1]
+      d = representations['quaternion'][1]
-      o = damask.Rotation.fromQuaternion(list(map(float,table.data[column:column+l])))
+      o = damask.Rotation.fromQuaternion(list(map(float,table.data[column:column+d])))
 
-    o= r*o*R                                                                                        # apply additional lab and crystal frame rotations
+    o = r*o*R                                                                                       # apply additional lab and crystal frame rotations
+
+    if options.lattice is not None:
+      o = damask.Orientation(rotation = o,lattice = options.lattice).reduced().rotation
 
     for output in options.output:
       if   output == 'quaternion': table.data_append(o.asQuaternion())
       elif output == 'rodrigues':  table.data_append(o.asRodrigues())
+      elif output == 'Rodrigues':  table.data_append(o.asRodrigues(vector=True))
       elif output == 'eulers':     table.data_append(o.asEulers(degrees=options.degrees))
       elif output == 'matrix':     table.data_append(o.asMatrix())
       elif output == 'angleaxis':  table.data_append(o.asAxisAngle(degrees=options.degrees))

python/damask/geom.py

@@ -73,14 +73,15 @@ class Geom():
     
     if size is not None and rescale:
       raise ValueError('Either set size explicitly or rescale automatically')
-    elif size is not None:
+
+    self.set_microstructure(microstructure)
+    self.set_origin(origin)
+
+    if size is not None:
       self.set_size(size)
     elif rescale:
       self.set_size(self.get_grid()/grid_old*self.size)
 
-    self.set_microstructure(microstructure)
-    self.set_origin(origin)
-   
     message = ['grid     a b c:      {}'.format(' x '.join(map(str,grid_old)))]
     if np.any(grid_old != self.get_grid()):
       message[-1] = util.delete(message[-1])

src/crystallite.f90

@@ -1274,7 +1274,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
       if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 &
           .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
                  .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then
-        write(6,'(a,i3,/)')                   '<< CRYST integrateStress >> iteration ', NiterationStressLp
+        write(6,'(a,i3,/)')                   '<< CRYST integrateStress >> Lp iteration ', NiterationStressLp
         write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Lpguess', transpose(Lpguess)
         write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Lp_constitutive', transpose(Lp_constitutive)
         write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Fi', transpose(Fi_new)
@@ -1377,6 +1377,7 @@ logical function integrateStress(ipc,ip,el,timeFraction)
       if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 &
           .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
                  .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then
+        write(6,'(a,i3,/)')                  '<< CRYST integrateStress >> Li iteration ', NiterationStressLi
         write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Li_constitutive', transpose(Li_constitutive)
         write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Liguess', transpose(Liguess)
       endif
@@ -1405,6 +1406,13 @@ logical function integrateStress(ipc,ip,el,timeFraction)
     else                                                                                            ! not converged and residuum not improved...
       steplengthLi   = num%subStepSizeLi * steplengthLi                                             ! ...try with smaller step length in same direction
       Liguess        = Liguess_old + steplengthLi * deltaLi
+#ifdef DEBUG
+        if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 &
+            .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
+                   .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then
+          write(6,'(a,1x,f7.4)') '<< CRYST integrateStress >> linear search for Liguess with step', steplengthLi
+        endif
+#endif
       cycle LiLoop
     endif
  
@@ -1450,6 +1458,13 @@ logical function integrateStress(ipc,ip,el,timeFraction)
     jacoCounterLi = jacoCounterLi + 1
  
     Liguess = Liguess + steplengthLi * deltaLi
+#ifdef DEBUG
+        if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 &
+            .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
+                   .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then
+          write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> corrected Liguess by', transpose(deltaLi)
+        endif
+#endif
   enddo LiLoop
  
   !* calculate new plastic and elastic deformation gradient

src/grid/DAMASK_grid.f90

@@ -118,6 +118,7 @@ program DAMASK_spectral
 ! assign mechanics solver depending on selected type
  select case (trim(config_numerics%getString('spectral_solver',defaultVal='basic')))
    case ('basic')
+     write(6,'(a)')   ' Basic solution algorithm'
      mech_init     => grid_mech_spectral_basic_init
      mech_forward  => grid_mech_spectral_basic_forward
      mech_solution => grid_mech_spectral_basic_solution
@@ -125,6 +126,7 @@ program DAMASK_spectral
    case ('polarisation')
      if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
        call IO_warning(42, ext_msg='debug Divergence')
+     write(6,'(a)')   ' Polarisation solution algorithm'
      mech_init     => grid_mech_spectral_polarisation_init
      mech_forward  => grid_mech_spectral_polarisation_forward
      mech_solution => grid_mech_spectral_polarisation_solution
@@ -132,6 +134,7 @@ program DAMASK_spectral
    case ('fem')
      if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
        call IO_warning(42, ext_msg='debug Divergence')
+     write(6,'(a)')   ' FEM solution algorithm'
      mech_init     => grid_mech_FEM_init
      mech_forward  => grid_mech_FEM_forward
      mech_solution => grid_mech_FEM_solution

src/plastic_isotropic.f90

@@ -283,49 +283,52 @@ end subroutine plastic_isotropic_LpAndItsTangent
 
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates plastic velocity gradient and its tangent
-! ToDo: Rename Tstar to Mi?
 !--------------------------------------------------------------------------------------------------
-subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar,Tstar,instance,of)
+subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
  
   real(pReal), dimension(3,3), intent(out) :: &
     Li                                                                                              !< inleastic velocity gradient
   real(pReal), dimension(3,3,3,3), intent(out)  :: &
-    dLi_dTstar                                                                                      !< derivative of Li with respect to the Mandel stress
+    dLi_dMi                                                                                         !< derivative of Li with respect to the Mandel stress
  
   real(pReal), dimension(3,3),   intent(in) :: &
-    Tstar                                                                                           !< Mandel stress ToDo: Mi?
+    Mi                                                                                              !< Mandel stress 
   integer,                       intent(in) :: &
     instance, &
     of
  
   real(pReal), dimension(3,3) :: &
-    Tstar_sph                                                                                       !< sphiatoric part of the Mandel stress
+    Mi_sph                                                                                          !< spherical part of the Mandel stress
   real(pReal) :: &
-    dot_gamma, &                                                                                    !< strainrate
+    dot_gamma, &                                                                                    !< shear rate
-    norm_Tstar_sph, &                                                                               !< euclidean norm of Tstar_sph
+    tr                                                                                               !< pressure
-    squarenorm_Tstar_sph                                                                            !< square of the euclidean norm of Tstar_sph
   integer :: &
     k, l, m, n
  
   associate(prm => param(instance), stt => state(instance))
  
-  Tstar_sph = math_spherical33(Tstar)
+  tr=math_trace33(math_spherical33(Mi))
-  squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph,Tstar_sph)
+
-  norm_Tstar_sph = sqrt(squarenorm_Tstar_sph)
+  if (prm%dilatation .and. abs(tr) > 0.0_pReal) then                                                 ! no stress or J2 plasticity --> Li and its derivative are zero
+    Li = math_I3 &
+       * prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(of))**(-prm%n) &
+       * tr * abs(tr)**(prm%n-1.0_pReal)
  
-  if (prm%dilatation .and. norm_Tstar_sph > 0.0_pReal) then                                         ! no stress or J2 plastitiy --> Li and its derivative are zero
+#ifdef DEBUG
-    dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Tstar_sph /(prm%M*stt%xi(of))) **prm%n
+    if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 &
- 
+        .and. (of == prm%of_debug .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then
-    Li = math_I3/sqrt(3.0_pReal) * dot_gamma/prm%M
+      write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> pressure / MPa', tr/3.0_pReal*1.0e-6_pReal
+      write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> gdot', prm%dot_gamma_0 * (3.0_pReal*prm%M*stt%xi(of))**(-prm%n) &
+                                                           * tr * abs(tr)**(prm%n-1.0_pReal)
+    end if
+#endif
+
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
-      dLi_dTstar(k,l,m,n) = (prm%n-1.0_pReal) * Tstar_sph(k,l)*Tstar_sph(m,n) / squarenorm_Tstar_sph
+      dLi_dMi(k,l,m,n) = n / tr * Li(k,l) * math_I3(m,n)
-    forall (k=1:3,l=1:3) &
+
-      dLi_dTstar(k,l,k,l) = dLi_dTstar(k,l,k,l) + 1.0_pReal
- 
-    dLi_dTstar = dot_gamma / prm%M * dLi_dTstar / norm_Tstar_sph
   else
-    Li = 0.0_pReal
+    Li      = 0.0_pReal
-    dLi_dTstar = 0.0_pReal
+    dLi_dMi = 0.0_pReal
   endif
  
   end associate

src/plastic_phenopowerlaw.f90

@@ -250,6 +250,7 @@ subroutine plastic_phenopowerlaw_init
 !--------------------------------------------------------------------------------------------------
 ! slip-twin related parameters
    slipAndTwinActive: if (prm%totalNslip > 0 .and. prm%totalNtwin > 0) then
+     prm%h0_TwinSlip          = config%getFloat('h0_twinslip')
      prm%interaction_SlipTwin = lattice_interaction_SlipByTwin(prm%Nslip,prm%Ntwin,&
                                                                config%getFloats('interaction_sliptwin'), &
                                                                config%getString('lattice_structure'))