Merge branch 'python-polishing' into lattice-into-orientation

PRIVATE

@@ -1 +1 @@
-Subproject commit ab64793bb04c506d815ebc850672ed0f2d013e67
+Subproject commit 95f7faea920dd6956884e4a55f72e5d5b1ffcdc8

VERSION

@@ -1 +1 @@
-v3.0.0-alpha3-40-ga1c46b445
+v3.0.0-alpha3-70-g7774ca721

env/DAMASK.sh

@@ -30,7 +30,6 @@ PATH=${DAMASK_ROOT}/bin:$PATH
 SOLVER=$(type -p DAMASK_grid || true 2>/dev/null)
 [ "x$SOLVER" == "x" ] && SOLVER=$(blink 'Not found!')
 
-[ "x$OMP_NUM_THREADS" == "x" ] && OMP_NUM_THREADS=1
 
 # currently, there is no information that unlimited stack size causes problems
 # still, http://software.intel.com/en-us/forums/topic/501500 suggest to fix it
@@ -55,9 +54,8 @@ if [ ! -z "$PS1" ]; then
     [[ $(canonicalPath "$PETSC_DIR") == $PETSC_DIR ]] \
     || echo "               ~~> "$(canonicalPath "$PETSC_DIR")
   fi
-  echo -n "MSC.Marc/Mentat    "
+  [ "x$PETSC_ARCH" != "x" ] && echo "PETSc architecture $PETSC_ARCH"
-  [ -d $MSC_ROOT ] && echo $MSC_ROOT || blink $MSC_ROOT
+  [ "x$OMP_NUM_THREADS" == "x" ] && export OMP_NUM_THREADS=4
-  echo
   echo "Multithreading     OMP_NUM_THREADS=$OMP_NUM_THREADS"
   echo -n "heap  size         "
    [[ "$(ulimit -d)" == "unlimited" ]] \

env/DAMASK.zsh

@@ -20,7 +20,6 @@ PATH=${DAMASK_ROOT}/bin:$PATH
 SOLVER=$(which DAMASK_grid || true 2>/dev/null)
 [[ "x$SOLVER" == "x" ]] && SOLVER=$(blink 'Not found!')
 
-[[ "x$OMP_NUM_THREADS" == "x" ]] && OMP_NUM_THREADS=1
 
 # currently, there is no information that unlimited stack size causes problems
 # still, http://software.intel.com/en-us/forums/topic/501500 suggest to fix it
@@ -45,11 +44,8 @@ if [ ! -z "$PS1" ]; then
     [[ $(canonicalPath "$PETSC_DIR") == $PETSC_DIR ]] \
     || echo "               ~~> "$(canonicalPath "$PETSC_DIR")
   fi
-  [[ "x$PETSC_ARCH" == "x" ]] \
+  [[ "x$PETSC_ARCH" != "x" ]] && echo "PETSc architecture $PETSC_ARCH"
-  || echo "PETSc architecture $PETSC_ARCH"
+  [[ "x$OMP_NUM_THREADS" == "x" ]] && export OMP_NUM_THREADS=4
-  echo -n "MSC.Marc/Mentat    "
-  [ -d $MSC_ROOT ] && echo $MSC_ROOT || blink $MSC_ROOT
-  echo
   echo "Multithreading     OMP_NUM_THREADS=$OMP_NUM_THREADS"
   echo -n "heap  size         "
    [[ "$(ulimit -d)" == "unlimited" ]] \

examples/config/phase/thermal/Al.yaml

@@ -0,0 +1,7 @@
+references:
+  - www.engineeringtoolbox.com/thermal-conductivity-metals-d_858.html
+  - www.engineeringtoolbox.com/specific-heat-metals-d_152.html
+c_p: 0.91e3
+K_11: 236.0
+K_22: 236.0
+K_33: 236.0

examples/config/phase/thermal/Steel-0.5C.yaml

@@ -0,0 +1,7 @@
+references:
+  - www.engineeringtoolbox.com/thermal-conductivity-metals-d_858.html
+  - www.engineeringtoolbox.com/specific-heat-metals-d_152.html
+c_p: 0.49e3
+K_11: 54.0
+K_22: 54.0
+K_33: 54.0

installation/mods_MarcMentat/apply_DAMASK_modifications.py

@@ -7,6 +7,11 @@ from pathlib import Path
 
 import damask
 
+def copy_and_replace(in_file,dst,msc_root,editor):
+    with open(in_file) as f_in, open(dst/Path(in_file).name,'w') as f_out:
+        f_out.write(f_in.read().replace('%INSTALLDIR%',msc_root).replace('%EDITOR%',editor))
+
+
 parser = argparse.ArgumentParser(
                   description='Apply DAMASK modification to MSC.Marc/Mentat',
                   prog = Path(__file__).name,
@@ -17,34 +22,25 @@ parser.add_argument('--editor', dest='editor', metavar='string', default='vi',
 parser.add_argument('--msc-root', dest='msc_root', metavar='string',
                     default=damask.solver._marc._msc_root,
                     help='MSC.Marc/Mentat root directory')
-parser.add_argument('--msc-version', dest='msc_version', metavar='string',
+parser.add_argument('--msc-version', dest='msc_version', type=float, metavar='float',
                     default=damask.solver._marc._msc_version,
                     help='MSC.Marc/Mentat version')
-parser.add_argument('--damask-root', dest='damask_root', type=float, metavar = 'float',
+parser.add_argument('--damask-root', dest='damask_root', metavar = 'string',
                     default=damask.solver._marc._damask_root,
                     help='DAMASK root directory')
 
-
 args = parser.parse_args()
 msc_root    = Path(args.msc_root)
 damask_root = Path(args.damask_root)
 msc_version = args.msc_version
 
-def copy_and_replace(in_file,dst):
-    with open(in_file) as f:
-        content = f.read()
-    content = content.replace('%INSTALLDIR%',str(msc_root))
-    content = content.replace('%EDITOR%', args.editor)
-    with open(dst/Path(in_file).name,'w') as f:
-        f.write(content)
-
 
 print('adapting Marc tools...\n')
 
 src = damask_root/f'installation/mods_MarcMentat/{msc_version}/Marc_tools'
 dst = msc_root/f'marc{msc_version}/tools'
 for in_file in glob.glob(str(src/'*damask*')) + [str(src/'include_linux64')]:
-    copy_and_replace(in_file,dst)
+    copy_and_replace(in_file,dst,args.msc_root,args.editor)
 
 
 print('adapting Mentat scripts and menus...\n')
@@ -52,12 +48,12 @@ print('adapting Mentat scripts and menus...\n')
 src = damask_root/f'installation/mods_MarcMentat/{msc_version}/Mentat_bin'
 dst = msc_root/f'mentat{msc_version}/bin'
 for in_file in glob.glob(str(src/'*[!.original]')):
-    copy_and_replace(in_file,dst)
+    copy_and_replace(in_file,dst,args.msc_root,args.editor)
 
 src = damask_root/f'installation/mods_MarcMentat/{msc_version}/Mentat_menus'
 dst = msc_root/f'mentat{msc_version}/menus'
 for in_file in glob.glob(str(src/'job_run.ms')):
-    copy_and_replace(in_file,dst)
+    copy_and_replace(in_file,dst,args.msc_root,args.editor)
 
 
 print('compiling Mentat menu binaries...')

python/damask/_colormap.py

@@ -110,13 +110,13 @@ class Colormap(mpl.colors.ListedColormap):
             low_,high_ = map(Colormap._rgb2msh,low_high)
 
         elif model.lower() == 'hsv':
-            if np.any(low_high<0) or np.any(low_high[:,1:3]>1) or np.any(low_high[:,0]>360):
+            if np.any(low_high<0) or np.any(low_high>[360,1,1]):
                 raise ValueError(f'HSV color {low} | {high} are out of range.')
 
             low_,high_ = map(Colormap._hsv2msh,low_high)
 
         elif model.lower() == 'hsl':
-            if np.any(low_high<0) or np.any(low_high[:,1:3]>1) or np.any(low_high[:,0]>360):
+            if np.any(low_high<0) or np.any(low_high>[360,1,1]):
                 raise ValueError(f'HSL color {low} | {high} are out of range.')
 
             low_,high_ = map(Colormap._hsl2msh,low_high)

python/damask/_grid.py

@@ -392,7 +392,7 @@ class Grid:
             seeds_p   = seeds
             coords    = grid_filters.coordinates0_point(cells,size).reshape(-1,3)
 
-        pool = mp.Pool(int(os.environ.get('OMP_NUM_THREADS',1)))
+        pool = mp.Pool(int(os.environ.get('OMP_NUM_THREADS',4)))
         result = pool.map_async(partial(Grid._find_closest_seed,seeds_p,weights_p), [coord for coord in coords])
         pool.close()
         pool.join()
@@ -436,8 +436,12 @@ class Grid:
 
         """
         coords = grid_filters.coordinates0_point(cells,size).reshape(-1,3)
-        KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
+        tree = spatial.cKDTree(seeds,boxsize=size) if periodic else \
-        devNull,material_ = KDTree.query(coords)
+               spatial.cKDTree(seeds)
+        try:
+            material_ = tree.query(coords, workers = int(os.environ.get('OMP_NUM_THREADS',4)))[1]
+        except TypeError:
+            material_ = tree.query(coords, n_jobs = int(os.environ.get('OMP_NUM_THREADS',4)))[1]    # scipy <1.6
 
         return Grid(material = (material_ if material is None else material[material_]).reshape(cells),
                     size     = size,

python/damask/_orientation.py

@@ -134,30 +134,15 @@ class Orientation(Rotation):
         """
         Rotation.__init__(self) if rotation is None else Rotation.__init__(self,rotation=rotation)
 
-        if (    lattice not in _lattice_symmetries
-            and lattice not in _crystal_families):
-            raise KeyError(f'Lattice "{lattice}" is unknown')
-
-        self.family     = None
-        self.lattice    = None
-        self.a          = None
-        self.b          = None
-        self.c          = None
-        self.alpha      = None
-        self.beta       = None
-        self.gamma      = None
         self.kinematics = None
 
         if lattice in _lattice_symmetries:
             self.family  = _lattice_symmetries[lattice]
             self.lattice = lattice
+
             self.a = 1 if a is None else a
             self.b = b
             self.c = c
-            self.alpha = (np.radians(alpha) if degrees else alpha) if alpha is not None else None
-            self.beta  = (np.radians(beta)  if degrees else beta)  if beta  is not None else None
-            self.gamma = (np.radians(gamma) if degrees else gamma) if gamma is not None else None
-
             self.a = float(self.a) if self.a is not None else \
                      (self.b / self.ratio['b'] if self.b is not None and self.ratio['b'] is not None else
                       self.c / self.ratio['c'] if self.c is not None and self.ratio['c'] is not None else None)
@@ -169,9 +154,13 @@ class Orientation(Rotation):
                      (self.a * self.ratio['c'] if self.a is not None and self.ratio['c'] is not None else
                       self.b / self.ratio['b'] * self.ratio['c']
                       if self.c is not None and self.ratio['b'] is not None and self.ratio['c'] is not None else None)
-            self.alpha = self.alpha if self.alpha is not None else self.immutable['alpha'] if 'alpha' in self.immutable else None
+
-            self.beta  = self.beta  if self.beta  is not None else self.immutable['beta']  if 'beta'  in self.immutable else None
+            self.alpha = np.radians(alpha) if degrees and alpha is not None else alpha
-            self.gamma = self.gamma if self.gamma is not None else self.immutable['gamma'] if 'gamma' in self.immutable else None
+            self.beta  = np.radians(beta)  if degrees and beta  is not None else beta
+            self.gamma = np.radians(gamma) if degrees and gamma is not None else gamma
+            if self.alpha is None and 'alpha' in self.immutable: self.alpha = self.immutable['alpha']
+            if self.beta  is None and 'beta'  in self.immutable: self.beta  = self.immutable['beta']
+            if self.gamma is None and 'gamma' in self.immutable: self.gamma = self.immutable['gamma']
 
             if \
                 (self.a     is None) \
@@ -197,7 +186,13 @@ class Orientation(Rotation):
                                          {'direction':self.Bravais_to_Miller(uvtw=master[m][:,0:4]),
                                           'plane':    self.Bravais_to_Miller(hkil=master[m][:,4:8])}
         elif lattice in _crystal_families:
-            self.family = lattice
+            self.family  = lattice
+            self.lattice = None
+
+            self.a = self.b = self.c = None
+            self.alpha = self.beta = self.gamma = None
+        else:
+            raise KeyError(f'Lattice "{lattice}" is unknown')
 
 
     def __repr__(self):

python/damask/_result.py

@@ -1179,7 +1179,7 @@ class Result:
 
         """
         chunk_size = 1024**2//8
-        pool = mp.Pool(int(os.environ.get('OMP_NUM_THREADS',1)))
+        pool = mp.Pool(int(os.environ.get('OMP_NUM_THREADS',4)))
         lock = mp.Manager().Lock()
 
         groups = []

python/damask/_vtk.py

@@ -150,9 +150,8 @@ class VTK:
         ----------
         fname : str or pathlib.Path
             Filename for reading. Valid extensions are .vtr, .vtu, .vtp, and .vtk.
-        dataset_type : str, optional
+        dataset_type : {'vtkRectilinearGrid', 'vtkUnstructuredGrid', 'vtkPolyData'}, optional
             Name of the vtk.vtkDataSet subclass when opening a .vtk file.
-            Valid types are vtkRectilinearGrid, vtkUnstructuredGrid, and vtkPolyData.
 
         Returns
         -------

python/damask/mechanics.py

@@ -124,9 +124,6 @@ def strain(F,t,m):
     """
     Calculate strain tensor (Seth–Hill family).
 
-    For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and
-    https://de.wikipedia.org/wiki/Verzerrungstensor
-
     Parameters
     ----------
     F : numpy.ndarray of shape (...,3,3)
@@ -142,6 +139,11 @@ def strain(F,t,m):
     epsilon : numpy.ndarray of shape (...,3,3)
         Strain of F.
 
+    References
+    ----------
+    https://en.wikipedia.org/wiki/Finite_strain_theory
+    https://de.wikipedia.org/wiki/Verzerrungstensor
+
     """
     if   t == 'V':
         w,n = _np.linalg.eigh(deformation_Cauchy_Green_left(F))
@@ -150,7 +152,6 @@ def strain(F,t,m):
 
     if   m > 0.0:
         eps = 1.0/(2.0*abs(m)) * (+ _np.einsum('...j,...kj,...lj',w**m,n,n) - _np.eye(3))
-
     elif m < 0.0:
         eps = 1.0/(2.0*abs(m)) * (- _np.einsum('...j,...kj,...lj',w**m,n,n) + _np.eye(3))
     else:

python/damask/seeds.py

@@ -3,8 +3,8 @@
 from scipy import spatial as _spatial
 import numpy as _np
 
-from . import util
+from . import util as _util
-from . import grid_filters
+from . import grid_filters as _grid_filters
 
 
 def from_random(size,N_seeds,cells=None,rng_seed=None):
@@ -34,7 +34,7 @@ def from_random(size,N_seeds,cells=None,rng_seed=None):
     if cells is None:
         coords = rng.random((N_seeds,3)) * size
     else:
-        grid_coords = grid_filters.coordinates0_point(cells,size).reshape(-1,3,order='F')
+        grid_coords = _grid_filters.coordinates0_point(cells,size).reshape(-1,3,order='F')
         coords = grid_coords[rng.choice(_np.prod(cells),N_seeds, replace=False)] \
                + _np.broadcast_to(size/cells,(N_seeds,3))*(rng.random((N_seeds,3))*.5-.25)          # wobble without leaving cells
 
@@ -73,12 +73,12 @@ def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,rng_seed=
 
     s = 1
     i = 0
-    progress = util._ProgressBar(N_seeds+1,'',50)
+    progress = _util._ProgressBar(N_seeds+1,'',50)
     while s < N_seeds:
         candidates = rng.random((N_candidates,3))*_np.broadcast_to(size,(N_candidates,3))
         tree = _spatial.cKDTree(coords[:s],boxsize=size) if periodic else \
                _spatial.cKDTree(coords[:s])
-        distances, dev_null = tree.query(candidates)
+        distances = tree.query(candidates)[0]
         best = distances.argmax()
         if distances[best] > distance:                                                              # require minimum separation
             coords[s] = candidates[best]                                                            # maximum separation to existing point cloud
@@ -120,7 +120,7 @@ def from_grid(grid,selection=None,invert=False,average=False,periodic=True):
     material = grid.material.reshape((-1,1),order='F')
     mask = _np.full(grid.cells.prod(),True,dtype=bool) if selection is None else \
            _np.isin(material,selection,invert=invert).flatten()
-    coords = grid_filters.coordinates0_point(grid.cells,grid.size).reshape(-1,3,order='F')
+    coords = _grid_filters.coordinates0_point(grid.cells,grid.size).reshape(-1,3,order='F')
 
     if not average:
         return (coords[mask],material[mask])

python/damask/util.py

@@ -1,3 +1,5 @@
+"""Miscellaneous helper functionality."""
+
 import sys
 import datetime
 import os
@@ -177,6 +179,16 @@ def execute(cmd,wd='./',env=None):
 
 
 def natural_sort(key):
+    """
+    Natural sort.
+
+    For use in python's 'sorted'.
+
+    References
+    ----------
+    https://en.wikipedia.org/wiki/Natural_sort_order
+
+    """
     convert = lambda text: int(text) if text.isdigit() else text
     return [ convert(c) for c in re.split('([0-9]+)', key) ]
 
@@ -191,9 +203,9 @@ def show_progress(iterable,N_iter=None,prefix='',bar_length=50):
     ----------
     iterable : iterable/function with yield statement
         Iterable (or function with yield statement) to be decorated.
-    N_iter : int
+    N_iter : int, optional
         Total # of iterations. Needed if number of iterations can not be obtained as len(iterable).
-    prefix : str, optional.
+    prefix : str, optional
         Prefix string.
     bar_length : int, optional
         Character length of bar. Defaults to 50.
@@ -509,6 +521,7 @@ def dict_prune(d):
             v = dict_prune(v)
         if not isinstance(v,dict) or v != {}:
             new[k] = v
+
     return new
 
 

src/HDF5_utilities.f90

@@ -1775,7 +1775,7 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
 
   integer,          dimension(worldsize)      :: writeSize                                          !< contribution of all processes
   integer(HID_T) ::  dcpl
-  integer :: ierr, hdferr
+  integer :: ierr, hdferr, HDF5_major, HDF5_minor, HDF5_release
   integer(HSIZE_T), parameter :: chunkSize = 1024_HSIZE_T**2/8_HSIZE_T
 
 !-------------------------------------------------------------------------------------------------
@@ -1808,14 +1808,18 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
   call h5pcreate_f(H5P_DATASET_CREATE_F, dcpl, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
   if(product(totalShape) >= chunkSize*2_HSIZE_T) then
-    call h5pset_chunk_f(dcpl, size(totalShape), getChunks(totalShape,chunkSize), hdferr)
+    call H5get_libversion_f(HDF5_major,HDF5_minor,HDF5_release,hdferr)
-    if(hdferr < 0) error stop 'HDF5 error'
+    if (hdferr < 0) error stop 'HDF5 error'
-    call h5pset_shuffle_f(dcpl, hdferr)
+    if (HDF5_major == 1 .and. HDF5_minor >= 12) then                                                ! https://forum.hdfgroup.org/t/6186
-    if(hdferr < 0) error stop 'HDF5 error'
+      call h5pset_chunk_f(dcpl, size(totalShape), getChunks(totalShape,chunkSize), hdferr)
-    call h5pset_deflate_f(dcpl, 6, hdferr)
+      if (hdferr < 0) error stop 'HDF5 error'
-    if(hdferr < 0) error stop 'HDF5 error'
+      call h5pset_shuffle_f(dcpl, hdferr)
-    call h5pset_Fletcher32_f(dcpl,hdferr)
+      if (hdferr < 0) error stop 'HDF5 error'
-    if(hdferr < 0) error stop 'HDF5 error'
+      call h5pset_deflate_f(dcpl, 6, hdferr)
+      if (hdferr < 0) error stop 'HDF5 error'
+      call h5pset_Fletcher32_f(dcpl,hdferr)
+      if (hdferr < 0) error stop 'HDF5 error'
+    endif
   endif
 
 !--------------------------------------------------------------------------------------------------

src/commercialFEM_fileList.f90

@@ -21,6 +21,7 @@
 #include "lattice.f90"
 #include "phase.f90"
 #include "phase_mechanical.f90"
+#include "phase_mechanical_elastic.f90"
 #include "phase_mechanical_plastic.f90"
 #include "phase_mechanical_plastic_none.f90"
 #include "phase_mechanical_plastic_isotropic.f90"

src/grid/grid_thermal_spectral.f90

@@ -90,7 +90,8 @@ subroutine grid_thermal_spectral_init(T_0)
 
 !--------------------------------------------------------------------------------------------------
 ! set default and user defined options for PETSc
- call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-thermal_snes_type ngmres',ierr)
+ call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-thermal_snes_type newtonls -thermal_snes_mf &
+                               &-thermal_snes_ksp_ew -thermal_ksp_type fgmres',ierr)
  CHKERRQ(ierr)
  call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr)
  CHKERRQ(ierr)

src/phase.f90

@@ -58,10 +58,6 @@ module phase
       grain
   end type tDebugOptions
 
-  integer, dimension(:), allocatable, public :: &                                                   !< ToDo: should be protected (bug in Intel compiler)
-    phase_elasticityInstance, &
-    phase_NstiffnessDegradations
-
   logical, dimension(:), allocatable, public :: &                                                   ! ToDo: should be protected (bug in Intel Compiler)
     phase_localPlasticity                                                                           !< flags phases with local constitutive law
 
@@ -298,7 +294,6 @@ module phase
   end interface
 
 
-
   type(tDebugOptions) :: debugConstitutive
 #if __INTEL_COMPILER >= 1900
   public :: &

src/phase_mechanical.f90

@@ -5,10 +5,6 @@ submodule(phase) mechanical
 
 
   enum, bind(c); enumerator :: &
-    ELASTICITY_UNDEFINED_ID, &
-    ELASTICITY_HOOKE_ID, &
-    STIFFNESS_DEGRADATION_UNDEFINED_ID, &
-    STIFFNESS_DEGRADATION_DAMAGE_ID, &
     PLASTICITY_UNDEFINED_ID, &
     PLASTICITY_NONE_ID, &
     PLASTICITY_ISOTROPIC_ID, &
@@ -23,11 +19,6 @@ submodule(phase) mechanical
     KINEMATICS_THERMAL_EXPANSION_ID
   end enum
 
-  integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:),   allocatable :: &
-    phase_elasticity                                                                                !< elasticity of each phase
-  integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)),     dimension(:,:), allocatable :: &
-    phase_stiffnessDegradation                                                                      !< active stiffness degradation mechanisms of each phase
-
   type(tTensorContainer), dimension(:), allocatable :: &
     ! current value
     phase_mechanical_Fe, &
@@ -57,9 +48,27 @@ submodule(phase) mechanical
       class(tNode), pointer :: phases
     end subroutine eigendeformation_init
 
+    module subroutine elastic_init(phases)
+      class(tNode), pointer :: phases
+    end subroutine elastic_init
+
     module subroutine plastic_init
     end subroutine plastic_init
 
+    module subroutine phase_hooke_SandItsTangents(S,dS_dFe,dS_dFi,Fe,Fi,ph,en)
+      integer, intent(in) :: &
+        ph, &
+        en
+      real(pReal),   intent(in),  dimension(3,3) :: &
+        Fe, &                                                                                       !< elastic deformation gradient
+        Fi                                                                                          !< intermediate deformation gradient
+      real(pReal),   intent(out), dimension(3,3) :: &
+        S                                                                                           !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
+      real(pReal),   intent(out), dimension(3,3,3,3) :: &
+        dS_dFe, &                                                                                   !< derivative of 2nd P-K stress with respect to elastic deformation gradient
+        dS_dFi                                                                                      !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
+    end subroutine phase_hooke_SandItsTangents
+    
     module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en)
       real(pReal), dimension(3,3),     intent(out) :: &
         Li                                                                                          !< inleastic velocity gradient
@@ -73,7 +82,6 @@ submodule(phase) mechanical
     end subroutine plastic_isotropic_LiAndItsTangent
 
     module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken)
-
       integer, intent(in) :: &
         co, &                                                                                           !< component-ID of integration point
         ip, &                                                                                           !< integration point
@@ -198,17 +206,11 @@ module subroutine mechanical_init(materials,phases)
     constituents, &
     constituent, &
     phase, &
-    mech, &
+    mech
-    elastic, &
-    stiffDegradation
 
   print'(/,a)', ' <<<+-  phase:mechanical init  -+>>>'
 
 !-------------------------------------------------------------------------------------------------
-! initialize elasticity (hooke)                         !ToDO: Maybe move to elastic submodule along with function homogenizedC?
-  allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID)
-  allocate(phase_elasticityInstance(phases%length), source = 0)
-  allocate(phase_NstiffnessDegradations(phases%length),source=0)
   allocate(output_constituent(phases%length))
 
   allocate(phase_mechanical_Fe(phases%length))
@@ -253,32 +255,8 @@ module subroutine mechanical_init(materials,phases)
 #else
     output_constituent(ph)%label  = mech%get_as1dString('output',defaultVal=emptyStringArray)
 #endif
-    elastic => mech%get('elastic')
-    if (IO_lc(elastic%get_asString('type')) == 'hooke') then ! accept small letter h for the moment
-      phase_elasticity(ph) = ELASTICITY_HOOKE_ID
-    else
-      call IO_error(200,ext_msg=elastic%get_asString('type'))
-    endif
-    stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList)                      ! check for stiffness degradation mechanisms
-    phase_NstiffnessDegradations(ph) = stiffDegradation%length
   enddo
 
-  allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),phases%length), &
-                        source=STIFFNESS_DEGRADATION_undefined_ID)
-
-  if(maxVal(phase_NstiffnessDegradations)/=0) then
-    do ph = 1, phases%length
-      phase => phases%get(ph)
-      mech    => phase%get('mechanical')
-      stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList)
-      do stiffDegradationCtr = 1, stiffDegradation%length
-        if(stiffDegradation%get_asString(stiffDegradationCtr) == 'damage') &
-            phase_stiffnessDegradation(stiffDegradationCtr,ph) = STIFFNESS_DEGRADATION_damage_ID
-      enddo
-    enddo
-  endif
-
-
   do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2)
     do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
       material     => materials%get(discretization_materialAt(el))
@@ -306,6 +284,9 @@ module subroutine mechanical_init(materials,phases)
   enddo; enddo
 
 
+! initialize elasticity
+  call elastic_init(phases)
+
 ! initialize plasticity
   allocate(plasticState(phases%length))
   allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID)
@@ -313,9 +294,6 @@ module subroutine mechanical_init(materials,phases)
 
   call plastic_init()
 
-  do ph = 1, phases%length
-    phase_elasticityInstance(ph) = count(phase_elasticity(1:ph) == phase_elasticity(ph))
-  enddo
 
   num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
 
@@ -348,51 +326,6 @@ module subroutine mechanical_init(materials,phases)
 end subroutine mechanical_init
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
-!> the elastic and intermediate deformation gradients using Hooke's law
-!--------------------------------------------------------------------------------------------------
-subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
-                                              Fe, Fi, ph, en)
-
-  integer, intent(in) :: &
-    ph, &
-    en
-  real(pReal),   intent(in),  dimension(3,3) :: &
-    Fe, &                                                                                           !< elastic deformation gradient
-    Fi                                                                                              !< intermediate deformation gradient
-  real(pReal),   intent(out), dimension(3,3) :: &
-    S                                                                                               !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
-  real(pReal),   intent(out), dimension(3,3,3,3) :: &
-    dS_dFe, &                                                                                       !< derivative of 2nd P-K stress with respect to elastic deformation gradient
-    dS_dFi                                                                                          !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
-
-  real(pReal), dimension(3,3) :: E
-  real(pReal), dimension(3,3,3,3) :: C
-  integer :: &
-    d, &                                                                                            !< counter in degradation loop
-    i, j
-
-  C = math_66toSym3333(phase_homogenizedC(ph,en))
-
-  DegradationLoop: do d = 1, phase_NstiffnessDegradations(ph)
-    degradationType: select case(phase_stiffnessDegradation(d,ph))
-      case (STIFFNESS_DEGRADATION_damage_ID) degradationType
-        C = C * damage_phi(ph,en)**2
-    end select degradationType
-  enddo DegradationLoop
-
-  E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3)                                                  !< Green-Lagrange strain in unloaded configuration
-  S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi))                                        !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
-
-  do i =1, 3;do j=1,3
-    dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi)))                !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
-    dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3))                             !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
-  enddo; enddo
-
-end subroutine phase_hooke_SandItsTangents
-
-
 module subroutine mechanical_results(group,ph)
 
   character(len=*), intent(in) :: group
@@ -1082,26 +1015,6 @@ module subroutine mechanical_forward()
 end subroutine mechanical_forward
 
 
-
-!--------------------------------------------------------------------------------------------------
-!> @brief returns the homogenize elasticity matrix
-!> ToDo: homogenizedC66 would be more consistent
-!--------------------------------------------------------------------------------------------------
-module function phase_homogenizedC(ph,en) result(C)
-
-  real(pReal), dimension(6,6) :: C
-  integer,      intent(in)    :: ph, en
-
-  plasticType: select case (phase_plasticity(ph))
-    case (PLASTICITY_DISLOTWIN_ID) plasticType
-     C = plastic_dislotwin_homogenizedC(ph,en)
-    case default plasticType
-     C = lattice_C66(1:6,1:6,ph)
-  end select plasticType
-
-end function phase_homogenizedC
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief calculate stress (P)
 !--------------------------------------------------------------------------------------------------

src/phase_mechanical_elastic.f90

@@ -0,0 +1,135 @@
+submodule(phase:mechanical) elastic
+
+  enum, bind(c); enumerator :: &
+    ELASTICITY_UNDEFINED_ID, &
+    ELASTICITY_HOOKE_ID, &
+    STIFFNESS_DEGRADATION_UNDEFINED_ID, &
+    STIFFNESS_DEGRADATION_DAMAGE_ID
+  end enum
+
+  integer, dimension(:), allocatable :: &
+    phase_NstiffnessDegradations
+ 
+  integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:),   allocatable :: &
+    phase_elasticity                                                                                !< elasticity of each phase
+  integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)),     dimension(:,:), allocatable :: &
+    phase_stiffnessDegradation                                                                      !< active stiffness degradation mechanisms of each phase
+
+contains
+
+
+module subroutine elastic_init(phases)
+
+  class(tNode), pointer :: &
+    phases
+
+  integer :: &
+    ph, &
+    stiffDegradationCtr
+  class(tNode), pointer :: &
+    phase, &
+    mech, &
+    elastic, &
+    stiffDegradation
+ 
+  print'(/,a)', ' <<<+-  phase:mechanical:elastic init  -+>>>'
+
+  allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID)
+  allocate(phase_NstiffnessDegradations(phases%length),source=0)
+  
+  do ph = 1, phases%length
+    phase   => phases%get(ph)
+    mech    => phase%get('mechanical')
+    elastic => mech%get('elastic')
+    if(IO_lc(elastic%get_asString('type')) == 'hooke') then ! accept small letter h for the moment
+      phase_elasticity(ph) = ELASTICITY_HOOKE_ID
+    else
+      call IO_error(200,ext_msg=elastic%get_asString('type'))
+    endif
+    stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList)                      ! check for stiffness degradation mechanisms
+    phase_NstiffnessDegradations(ph) = stiffDegradation%length
+  enddo
+
+  allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),phases%length), &
+                        source=STIFFNESS_DEGRADATION_undefined_ID)
+
+  if(maxVal(phase_NstiffnessDegradations)/=0) then
+    do ph = 1, phases%length
+      phase => phases%get(ph)
+      mech    => phase%get('mechanical')
+      stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList)
+      do stiffDegradationCtr = 1, stiffDegradation%length
+        if(stiffDegradation%get_asString(stiffDegradationCtr) == 'damage') &
+            phase_stiffnessDegradation(stiffDegradationCtr,ph) = STIFFNESS_DEGRADATION_damage_ID
+      enddo
+    enddo
+  endif
+
+end subroutine elastic_init
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
+!> the elastic and intermediate deformation gradients using Hooke's law
+!--------------------------------------------------------------------------------------------------
+module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
+                                              Fe, Fi, ph, en)
+
+  integer, intent(in) :: &
+    ph, &
+    en
+  real(pReal),   intent(in),  dimension(3,3) :: &
+    Fe, &                                                                                           !< elastic deformation gradient
+    Fi                                                                                              !< intermediate deformation gradient
+  real(pReal),   intent(out), dimension(3,3) :: &
+    S                                                                                               !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
+  real(pReal),   intent(out), dimension(3,3,3,3) :: &
+    dS_dFe, &                                                                                       !< derivative of 2nd P-K stress with respect to elastic deformation gradient
+    dS_dFi                                                                                          !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
+
+  real(pReal), dimension(3,3) :: E
+  real(pReal), dimension(3,3,3,3) :: C
+  integer :: &
+    d, &                                                                                            !< counter in degradation loop
+    i, j
+
+  C = math_66toSym3333(phase_homogenizedC(ph,en))
+
+  DegradationLoop: do d = 1, phase_NstiffnessDegradations(ph)
+    degradationType: select case(phase_stiffnessDegradation(d,ph))
+      case (STIFFNESS_DEGRADATION_damage_ID) degradationType
+        C = C * damage_phi(ph,en)**2
+    end select degradationType
+  enddo DegradationLoop
+
+  E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3)                                                  !< Green-Lagrange strain in unloaded configuration
+  S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi))                                        !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
+
+  do i =1, 3;do j=1,3
+    dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi)))                !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
+    dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3))                             !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
+  enddo; enddo
+
+end subroutine phase_hooke_SandItsTangents
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns the homogenized elasticity matrix
+!> ToDo: homogenizedC66 would be more consistent
+!--------------------------------------------------------------------------------------------------
+module function phase_homogenizedC(ph,en) result(C)
+
+  real(pReal), dimension(6,6) :: C
+  integer,      intent(in)    :: ph, en
+
+  plasticType: select case (phase_plasticity(ph))
+    case (PLASTICITY_DISLOTWIN_ID) plasticType
+     C = plastic_dislotwin_homogenizedC(ph,en)
+    case default plasticType
+     C = lattice_C66(1:6,1:6,ph)
+  end select plasticType
+
+end function phase_homogenizedC
+
+
+end submodule elastic