Merge remote-tracking branch 'origin/development' into 2d-table-yaml

PRIVATE

@@ -1 +1 @@
-Subproject commit 48dd9972d9023caa8b04226112dcdd57fa0be6af
+Subproject commit 13dfa0ee9d702782f0b7999f3f7fb2384f58d768

VERSION

@@ -1 +1 @@
-v3.0.0-alpha2-572-g7b89cb41f
+v3.0.0-alpha2-602-ge2d4ab427

examples/ConfigFiles/Phase_Phenopowerlaw_cpTi.yaml

@@ -1,13 +1,14 @@
 # M. Levy, Handbook of Elastic Properties of Solids, Liquids, and Gases (2001)
 # C. Zambaldi, "Orientation informed nanoindentation of a-titanium: Indentation pileup in hexagonal metals deforming by prismatic slip", J. Mater. Res., Vol. 27, No. 1, Jan 14, 2012
-Ti-alpha:
+# Better use values from L. Wang, Z. Zheng, H. Phukan, P. Kenesei, J.-S. Park, J. Lind, R.M. Suter, T.R. Bieler, Direct measurement of critical resolved shear stress of prismatic and basal slip in polycrystalline Ti using high energy X-ray diffraction microscopy, Acta Mater 2017
+cpTi:
   lattice: hP
   c/a: 1.587
   mechanics:
     output: [F, P, F_e, F_p, L_p, O]
     elasticity: {C_11: 160.0e9, C_12: 90.0e9, C_13: 66.0e9, C_33: 181.7e9, C_44: 46.5e9, type: hooke}
     plasticity:
-      N_sl: [3, 3, 0, 0, 12]
+      N_sl: [3, 3, 0, 6, 12]
       a_sl: 2.0
       dot_gamma_0_sl: 0.001
       h_0_sl_sl: 200e6
@@ -15,5 +16,5 @@ Ti-alpha:
       n_sl: 20
       output: [gamma_sl]
       type: phenopowerlaw
-      xi_0_sl: [349e6, 150e6, 0, 0, 1107e6]
+      xi_0_sl: [0.15e9, 0.09e9, 0, 0.20e9, 0.25e9]
-      xi_inf_sl: [568e6, 1502e6, 0, 0, 3420e6]
+      xi_inf_sl: [0.24e9, 0.5e9, 0, 0.6e9, 0.8e9]

python/damask/_config.py

@@ -1,5 +1,6 @@
 import copy
 from io import StringIO
+from collections.abc import Iterable
 import abc
 
 import numpy as np
@@ -44,6 +45,42 @@ class Config(dict):
     copy = __copy__
 
 
+    def __or__(self,other):
+        """
+        Update configuration with contents of other.
+
+        Parameters
+        ----------
+        other : damask.Config or dict
+            Key-value pairs that update self.
+
+        """
+        duplicate = self.copy()
+        duplicate.update(other)
+        return duplicate
+
+
+    def __ior__(self,other):
+        """Update configuration with contents of other."""
+        return self.__or__(other)
+
+
+    def delete(self,keys):
+        """
+        Remove configuration keys.
+
+        Parameters
+        ----------
+        keys : iterable or scalar
+            Label of the key(s) to remove.
+
+        """
+        duplicate = self.copy()
+        for k in keys if isinstance(keys, Iterable) and not isinstance(keys, str) else [keys]:
+            del duplicate[k]
+        return duplicate
+
+
     @classmethod
     def load(cls,fname):
         """
@@ -99,30 +136,6 @@ class Config(dict):
             fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs))
 
 
-    def add(self,d):
-        """
-        Add dictionary.
-
-        d : dict
-            Dictionary to append.
-        """
-        duplicate = self.copy()
-        duplicate.update(d)
-        return duplicate
-
-
-    def delete(self,key):
-        """
-        Delete item.
-
-        key : str or scalar
-            Label of the key to remove.
-        """
-        duplicate = self.copy()
-        del duplicate[key]
-        return duplicate
-
-
     @property
     @abc.abstractmethod
     def is_complete(self):

python/damask/_result.py

@@ -371,7 +371,7 @@ class Result:
         with h5py.File(self.fname,'r') as f:
             for dataset in sets:
                 for group in self.groups_with_datasets(dataset):
-                    path = os.path.join(group,dataset)
+                    path = '/'.join([group,dataset])
                     inc,prop,name,cat,item = (path.split('/') + ['']*5)[:5]
                     key = '/'.join([prop,name+tag])
                     if key not in inGeom:
@@ -388,15 +388,15 @@ class Result:
                                    np.nan,
                                    dtype=np.dtype(f[path]))
                     data[inGeom[key]] = (f[path] if len(shape)>1 else np.expand_dims(f[path],1))[inData[key]]
-                    path = (os.path.join(*([prop,name]+([cat] if cat else [])+([item] if item else []))) if split else path)+tag
+                    path = ('/'.join([prop,name]+([cat] if cat else [])+([item] if item else [])) if split else path)+tag
                     if split:
                         try:
-                            tbl[inc].add(path,data)
+                            tbl[inc] = tbl[inc].add(path,data)
                         except KeyError:
                             tbl[inc] = Table(data.reshape(self.N_materialpoints,-1),{path:data.shape[1:]})
                     else:
                         try:
-                            tbl.add(path,data)
+                            tbl = tbl.add(path,data)
                         except AttributeError:
                             tbl = Table(data.reshape(self.N_materialpoints,-1),{path:data.shape[1:]})
 

python/damask/_rotation.py

@@ -475,15 +475,17 @@ class Rotation:
 
         Parameters
         ----------
-        vector : bool, optional
+        compact : bool, optional
-            Return as actual Rodrigues-Frank vector, i.e. axis
+            Return as actual Rodrigues-Frank vector,
-            and angle argument are not separated.
+            i.e. axis and angle argument are not separated.
 
         Returns
         -------
-        rho : numpy.ndarray of shape (...,4) unless vector == True:
+        rho : numpy.ndarray of shape (...,4) containing 
-            numpy.ndarray of shape (...,3)
+            [n_1, n_2, n_3, tan(ω/2)], ǀnǀ = 1 and ω ∈ [0,π]
-            Rodrigues-Frank vector: [n_1, n_2, n_3, tan(ω/2)], ǀnǀ = 1 and ω ∈ [0,π].
+            unless compact == True:
+            numpy.ndarray of shape (...,3) containing 
+            tan(ω/2) [n_1, n_2, n_3], ω ∈ [0,π].
 
         """
         ro = Rotation._qu2ro(self.quaternion)

python/damask/_table.py

@@ -189,6 +189,11 @@ class Table:
         label : str
             Column label.
 
+        Returns
+        -------
+        data : numpy.ndarray
+            Array of column data.
+
         """
         if re.match(r'[0-9]*?_',label):
             idx,key = label.split('_',1)
@@ -212,6 +217,11 @@ class Table:
         info : str, optional
             Human-readable information about the new data.
 
+        Returns
+        -------
+        table : Table
+            Updated table.
+
         """
         dup = self.copy()
         dup._add_comment(label,data.shape[1:],info)
@@ -238,6 +248,11 @@ class Table:
         info : str, optional
             Human-readable information about the modified data.
 
+        Returns
+        -------
+        table : Table
+            Updated table.
+
         """
         dup = self.copy()
         dup._add_comment(label,data.shape[1:],info)
@@ -261,6 +276,11 @@ class Table:
         label : str
             Column label.
 
+        Returns
+        -------
+        table : Table
+            Updated table.
+
         """
         dup = self.copy()
         dup.data.drop(columns=label,inplace=True)
@@ -279,6 +299,11 @@ class Table:
         label_new : str or iterable of str
             New column label(s).
 
+        Returns
+        -------
+        table : Table
+            Updated table.
+
         """
         dup = self.copy()
         columns = dict(zip([old] if isinstance(old,str) else old,
@@ -300,6 +325,11 @@ class Table:
         ascending : bool or list, optional
             Set sort order.
 
+        Returns
+        -------
+        table : Table
+            Updated table.
+
         """
         dup = self.copy()
         dup._label_discrete()
@@ -320,6 +350,11 @@ class Table:
         other : Table
             Table to append.
 
+        Returns
+        -------
+        table : Table
+            Concatenated table.
+
         """
         if self.shapes != other.shapes or not self.data.columns.equals(other.data.columns):
             raise KeyError('Labels or shapes or order do not match')
@@ -340,6 +375,11 @@ class Table:
         other : Table
             Table to join.
 
+        Returns
+        -------
+        table : Table
+            Joined table.
+
         """
         if set(self.shapes) & set(other.shapes) or self.data.shape[0] != other.data.shape[0]:
             raise KeyError('Dublicated keys or row count mismatch')

python/tests/test_Config.py

@@ -23,8 +23,17 @@ class TestConfig:
             assert Config.load(f) == config
 
     def test_add_remove(self):
+        dummy = {'hello':'world','foo':'bar'}
         config = Config()
-        assert config.add({'hello':'world'}).delete('hello') == config
+        config |= dummy
+        assert config == Config() | dummy
+        config = config.delete(dummy)
+        assert config == Config()
+        assert (config |        dummy ).delete(        'hello'            ) == config | {'foo':'bar'}
+        assert (config |        dummy ).delete([       'hello',  'foo'   ]) == config
+        assert (config | Config(dummy)).delete({       'hello':1,'foo':2 }) == config
+        assert (config | Config(dummy)).delete(Config({'hello':1        })) == config | {'foo':'bar'}
+
 
     def test_repr(self,tmp_path):
         config = Config()

src/CPFEM.f90

@@ -178,11 +178,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
 
   if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward
 
-    chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP)))
+    !chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP)))
     !  case (THERMAL_conduction_ID) chosenThermal1
     !    temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = &
     !      temperature_inp
-    end select chosenThermal1
+    !end select chosenThermal1
     homogenization_F0(1:3,1:3,ma) = ffn
     homogenization_F(1:3,1:3,ma) = ffn1
 

src/grid/DAMASK_grid.f90

@@ -37,9 +37,10 @@ program DAMASK_grid
                                 f_out, &                                                            !< frequency of result writes
                                 f_restart                                                           !< frequency of restart writes
     logical ::                  estimate_rate                                                       !< follow trajectory of former loadcase
-    integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:)
   end type tLoadCase
 
+  integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:)
+
 !--------------------------------------------------------------------------------------------------
 ! variables related to information from load case and geom file
   real(pReal), dimension(9) :: temp_valueVector                                                     !< temporarily from loadcase file when reading in tensors (initialize to 0.0)
@@ -53,6 +54,7 @@ program DAMASK_grid
   integer, parameter :: &
     subStepFactor = 2                                                                               !< for each substep, divide the last time increment by 2.0
   real(pReal) :: &
+    T_0 = 300.0_pReal, &
     time = 0.0_pReal, &                                                                             !< elapsed time
     time0 = 0.0_pReal, &                                                                            !< begin of interval
     timeinc = 1.0_pReal, &                                                                          !< current time interval
@@ -78,8 +80,7 @@ program DAMASK_grid
     maxCutBack, &                                                                                   !< max number of cut backs
     stagItMax                                                                                       !< max number of field level staggered iterations
   character(len=pStringLen) :: &
-    incInfo, &
+    incInfo
-    loadcase_string
 
   type(tLoadCase), allocatable, dimension(:) :: loadCases                                           !< array of all load cases
   type(tSolutionState), allocatable, dimension(:) :: solres
@@ -98,10 +99,13 @@ program DAMASK_grid
     quit
   class (tNode), pointer :: &
     num_grid, &
-    debug_grid, &                                                                                   ! pointer to grid debug options
     config_load, &
     load_steps, &
     load_step, &
+    solver, &
+    initial_conditions, &
+    ic_thermal, &
+    thermal, &
     step_bc, &
     step_mech, &
     step_discretization, &
@@ -112,17 +116,11 @@ program DAMASK_grid
 ! init DAMASK (all modules)
 
   call CPFEM_initAll
-  print'(/,a)',   ' <<<+-  DAMASK_spectral init  -+>>>'; flush(IO_STDOUT)
+  print'(/,a)',   ' <<<+-  DAMASK_grid init  -+>>>'; flush(IO_STDOUT)
 
   print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019'
   print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
 
-!--------------------------------------------------------------------------------------------------
-! initialize field solver information
-  nActiveFields = 1
-  if (any(thermal_type  == THERMAL_conduction_ID  )) nActiveFields = nActiveFields + 1
-  if (any(damage_type   == DAMAGE_nonlocal_ID     )) nActiveFields = nActiveFields + 1
-  allocate(solres(nActiveFields))
 
 !-------------------------------------------------------------------------------------------------
 ! reading field paramters from numerics file and do sanity checks
@@ -133,19 +131,22 @@ program DAMASK_grid
   if (stagItMax < 0)    call IO_error(301,ext_msg='maxStaggeredIter')
   if (maxCutBack < 0)   call IO_error(301,ext_msg='maxCutBack')
 
+  config_load => YAML_parse_file(trim(interface_loadFile))
+  solver => config_load%get('solver')
+
 !--------------------------------------------------------------------------------------------------
 ! assign mechanics solver depending on selected type
 
-  debug_grid => config_debug%get('grid',defaultVal=emptyList)
+  nActiveFields = 1
-  select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic')))
+  select case (solver%get_asString('mechanical'))
-    case ('Basic')
+    case ('spectral_basic')
       mechanical_init         => grid_mechanical_spectral_basic_init
       mechanical_forward      => grid_mechanical_spectral_basic_forward
       mechanical_solution     => grid_mechanical_spectral_basic_solution
       mechanical_updateCoords => grid_mechanical_spectral_basic_updateCoords
       mechanical_restartWrite => grid_mechanical_spectral_basic_restartWrite
 
-    case ('Polarisation')
+    case ('spectral_polarization')
       mechanical_init         => grid_mechanical_spectral_polarisation_init
       mechanical_forward      => grid_mechanical_spectral_polarisation_forward
       mechanical_solution     => grid_mechanical_spectral_polarisation_solution
@@ -160,32 +161,36 @@ program DAMASK_grid
       mechanical_restartWrite => grid_mechanical_FEM_restartWrite
 
     case default
-      call IO_error(error_ID = 891, ext_msg = trim(num_grid%get_asString('solver')))
+      call IO_error(error_ID = 891, ext_msg = trim(solver%get_asString('mechanical')))
 
   end select
 
+!--------------------------------------------------------------------------------------------------
+! initialize field solver information
+  if (solver%get_asString('thermal',defaultVal = 'n/a') == 'spectral') nActiveFields = nActiveFields + 1
+  if (solver%get_asString('damage', defaultVal = 'n/a') == 'spectral') nActiveFields = nActiveFields + 1
+
+  allocate(solres(nActiveFields))
+  allocate(    ID(nActiveFields))
+
+  field = 1
+  ID(field) = FIELD_MECH_ID                                                                         ! mechanical active by default
+  thermalActive: if (solver%get_asString('thermal',defaultVal = 'n/a') == 'spectral') then
+    field = field + 1
+    ID(field) = FIELD_THERMAL_ID
+  endif thermalActive
+  damageActive: if (solver%get_asString('damage',defaultVal = 'n/a') == 'spectral') then
+    field = field + 1
+    ID(field) = FIELD_DAMAGE_ID
+  endif damageActive
+
 
 !--------------------------------------------------------------------------------------------------
-! reading information from load case file and to sanity checks
-  config_load => YAML_parse_file(trim(interface_loadFile))
-
   load_steps => config_load%get('loadstep')
   allocate(loadCases(load_steps%length))                                                            ! array of load cases
 
   do l = 1, load_steps%length
 
-    allocate(loadCases(l)%ID(nActiveFields))
-    field = 1
-    loadCases(l)%ID(field) = FIELD_MECH_ID                                                          ! mechanical active by default
-    thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then
-      field = field + 1
-      loadCases(l)%ID(field) = FIELD_THERMAL_ID
-    endif thermalActive
-    damageActive: if (any(damage_type == DAMAGE_nonlocal_ID)) then
-      field = field + 1
-      loadCases(l)%ID(field) = FIELD_DAMAGE_ID
-    endif damageActive
-
     load_step => load_steps%get(l)
     step_bc   => load_step%get('boundary_conditions')
     step_mech => step_bc%get('mechanical')
@@ -228,11 +233,9 @@ program DAMASK_grid
 
     loadCases(l)%f_restart = load_step%get_asInt('f_restart', defaultVal=huge(0))
     loadCases(l)%f_out     = load_step%get_asInt('f_out',     defaultVal=1)
-    loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. & 
+    loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. l>1)
-                                                       merge(.true.,.false.,l > 1))
 
     reportAndCheck: if (worldrank == 0) then
-      write (loadcase_string, '(i0)' ) l
       print'(/,a,i0)', ' load case: ', l
       print*, ' estimate_rate:', loadCases(l)%estimate_rate
       if (loadCases(l)%deformation%myType == 'L') then
@@ -292,7 +295,7 @@ program DAMASK_grid
       if (loadCases(l)%f_restart < huge(0)) &
         print'(a,i0)',   '  f_restart: ', loadCases(l)%f_restart
 
-      if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string)                 ! exit with error message
+      if (errorID > 0) call IO_error(error_ID = errorID, el = l)
 
     endif reportAndCheck
   enddo
@@ -301,12 +304,14 @@ program DAMASK_grid
 ! doing initialization depending on active solvers
   call spectral_Utilities_init
   do field = 1, nActiveFields
-    select case (loadCases(1)%ID(field))
+    select case (ID(field))
       case(FIELD_MECH_ID)
         call mechanical_init
 
       case(FIELD_THERMAL_ID)
-        call grid_thermal_spectral_init
+        initial_conditions => config_load%get('initial_conditions',defaultVal=emptyDict)
+        thermal => initial_conditions%get('thermal',defaultVal=emptyDict)
+        call grid_thermal_spectral_init(thermal%get_asFloat('T',defaultVal = T_0))
 
       case(FIELD_DAMAGE_ID)
         call grid_damage_spectral_init
@@ -377,7 +382,7 @@ program DAMASK_grid
 !--------------------------------------------------------------------------------------------------
 ! forward fields
           do field = 1, nActiveFields
-            select case(loadCases(l)%ID(field))
+            select case(ID(field))
               case(FIELD_MECH_ID)
                 call mechanical_forward (&
                         cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, &
@@ -397,7 +402,7 @@ program DAMASK_grid
           stagIterate = .true.
           do while (stagIterate)
             do field = 1, nActiveFields
-              select case(loadCases(l)%ID(field))
+              select case(ID(field))
                 case(FIELD_MECH_ID)
                   solres(field) = mechanical_solution(incInfo)
                 case(FIELD_THERMAL_ID)

src/grid/grid_mech_spectral_polarisation.f90

@@ -116,7 +116,7 @@ subroutine grid_mechanical_spectral_polarisation_init
     num_grid, &
     debug_grid
 
-  print'(/,a)', ' <<<+-  grid_mechanical_spectral_polarisation init  -+>>>'; flush(IO_STDOUT)
+  print'(/,a)', ' <<<+-  grid_mechanical_spectral_polarization init  -+>>>'; flush(IO_STDOUT)
 
   print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
   print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'

src/grid/grid_thermal_spectral.f90

@@ -61,7 +61,9 @@ contains
 !> @brief allocates all neccessary fields and fills them with data
 ! ToDo: Restart not implemented
 !--------------------------------------------------------------------------------------------------
-subroutine grid_thermal_spectral_init
+subroutine grid_thermal_spectral_init(T_0)
+
+  real(pReal), intent(in) :: T_0
 
   PetscInt, dimension(0:worldsize-1) :: localK
   integer :: i, j, k, ce
@@ -131,9 +133,10 @@ subroutine grid_thermal_spectral_init
   ce = 0
   do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
     ce = ce + 1
-    T_current(i,j,k) = homogenization_thermal_T(ce)
+    T_current(i,j,k) = T_0
     T_lastInc(i,j,k) = T_current(i,j,k)
     T_stagInc(i,j,k) = T_current(i,j,k)
+    call homogenization_thermal_setField(T_0,0.0_pReal,ce)
   enddo; enddo; enddo
   call DMDAVecGetArrayF90(thermal_grid,solution_vec,x_scal,ierr); CHKERRQ(ierr)                     !< get the data out of PETSc to work with
   x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current

src/homogenization.f90

@@ -35,15 +35,11 @@ module homogenization
     homogState, &
     damageState_h
 
-
+  integer(kind(THERMAL_isothermal_ID)),       dimension(:),   allocatable :: &
-  real(pReal), dimension(:), allocatable, public, protected :: &
-    thermal_initialT
-
-  integer(kind(THERMAL_isothermal_ID)),       dimension(:),   allocatable, public, protected :: &
     thermal_type                                                                                    !< thermal transport model
-  integer(kind(DAMAGE_none_ID)),              dimension(:),   allocatable, public, protected :: &
+  integer(kind(DAMAGE_none_ID)),              dimension(:),   allocatable :: &
     damage_type                                                                                     !< nonlocal damage model
-  integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:),   allocatable, public, protected :: &
+  integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:),   allocatable :: &
     homogenization_type                                                                             !< type of each homogenization
 
   type, private :: tNumerics_damage
@@ -133,10 +129,8 @@ module homogenization
 
 
     module function thermal_conduction_getConductivity(ce) result(K)
-
       integer, intent(in) :: ce
       real(pReal), dimension(3,3) :: K
-
     end function thermal_conduction_getConductivity
 
     module function thermal_conduction_getSpecificHeat(ce) result(c_P)
@@ -177,7 +171,6 @@ module homogenization
     end function damage_nonlocal_getMobility
 
     module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce)
-
       integer, intent(in) :: ce
       real(pReal),   intent(in) :: &
         phi
@@ -185,21 +178,17 @@ module homogenization
         phiDot, dPhiDot_dPhi
     end subroutine damage_nonlocal_getSourceAndItsTangent
 
-
     module subroutine damage_nonlocal_putNonLocalDamage(phi,ce)
-
       integer, intent(in) :: ce
       real(pReal),   intent(in) :: &
         phi
-
     end subroutine damage_nonlocal_putNonLocalDamage
 
     module subroutine damage_nonlocal_results(ho,group)
-
       integer,          intent(in) :: ho
       character(len=*), intent(in) :: group
-
     end subroutine damage_nonlocal_results
+
   end interface
 
   public ::  &
@@ -242,8 +231,7 @@ subroutine homogenization_init()
 
   allocate(homogState      (size(material_name_homogenization)))
   allocate(damageState_h   (size(material_name_homogenization)))
-  call material_parseHomogenization
+  call material_parseHomogenization()
-
 
   num_homog        => config_numerics%get('homogenization',defaultVal=emptyDict)
   num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict)
@@ -251,12 +239,10 @@ subroutine homogenization_init()
   num%nMPstate  = num_homogGeneric%get_asInt('nMPstate',defaultVal=10)
   if (num%nMPstate < 1) call IO_error(301,ext_msg='nMPstate')
 
-
   call mechanical_init(num_homog)
   call thermal_init()
   call damage_init()
-
+  call damage_nonlocal_init()
-  call damage_nonlocal_init
 
 
 end subroutine homogenization_init
@@ -341,8 +327,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
               print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
             terminallyIll = .true.                                                                  ! ...and kills all others
           endif
-         call thermal_homogenize(ip,el)
         enddo
+        call thermal_homogenize(ip,el)
       enddo
     enddo
     !$OMP END DO
@@ -556,7 +542,6 @@ subroutine material_parseHomogenization
   allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID)
   allocate(thermal_type(size(material_name_homogenization)),        source=THERMAL_isothermal_ID)
   allocate(damage_type (size(material_name_homogenization)),        source=DAMAGE_none_ID)
-  allocate(thermal_initialT(size(material_name_homogenization)),              source=300.0_pReal)
 
   do h=1, size(material_name_homogenization)
     homog => material_homogenization%get(h)
@@ -572,15 +557,10 @@ subroutine material_parseHomogenization
         call IO_error(500,ext_msg=homogMech%get_asString('type'))
     end select
 
-
     if (homog%contains('thermal')) then
       homogThermal => homog%get('thermal')
-        thermal_initialT(h) =  homogThermal%get_asFloat('T_0',defaultVal=300.0_pReal)
-
         select case (homogThermal%get_asString('type'))
-          case('isothermal')
+          case('pass')
-            thermal_type(h) = THERMAL_isothermal_ID
-          case('conduction')
             thermal_type(h) = THERMAL_conduction_ID
           case default
             call IO_error(500,ext_msg=homogThermal%get_asString('type'))
@@ -590,9 +570,7 @@ subroutine material_parseHomogenization
     if (homog%contains('damage')) then
       homogDamage => homog%get('damage')
         select case (homogDamage%get_asString('type'))
-          case('none')
+          case('pass')
-            damage_type(h) = DAMAGE_none_ID
-          case('nonlocal')
             damage_type(h) = DAMAGE_nonlocal_ID
           case default
             call IO_error(500,ext_msg=homogDamage%get_asString('type'))
@@ -600,7 +578,6 @@ subroutine material_parseHomogenization
     endif
   enddo
 
-
 end subroutine material_parseHomogenization
 
 

src/homogenization_thermal.f90

@@ -44,7 +44,7 @@ module subroutine thermal_init()
   allocate(current(configHomogenizations%length))
 
   do ho = 1, configHomogenizations%length
-    allocate(current(ho)%T(count(material_homogenizationAt2==ho)), source=thermal_initialT(ho))
+    allocate(current(ho)%T(count(material_homogenizationAt2==ho)), source=300.0_pReal)
     allocate(current(ho)%dot_T(count(material_homogenizationAt2==ho)), source=0.0_pReal)
     configHomogenization => configHomogenizations%get(ho)
     associate(prm => param(ho))

src/material.f90

@@ -30,7 +30,7 @@ module material
     material_homogenizationAt, &                                                                    !< homogenization ID of each element
     material_homogenizationAt2, &                                                                   !< per cell
     material_homogenizationMemberAt2                                                                !< cell
-  integer, dimension(:,:),   allocatable, public, protected :: &                                    ! (ip,elem)
+  integer, dimension(:,:),   allocatable :: &                                                       ! (ip,elem)
     material_homogenizationMemberAt                                                                 !< position of the element within its homogenization instance
   integer, dimension(:,:),   allocatable, public, protected :: &                                    ! (constituent,elem)
     material_phaseAt, &                                                                             !< phase ID of each element

src/phase_damage.f90

@@ -133,7 +133,7 @@ module subroutine damage_init
 
   integer :: &
     ph, &                                                                                           !< counter in phase loop
-    Nconstituents
+    Nmembers
   class(tNode), pointer :: &
    phases, &
    phase, &
@@ -151,10 +151,10 @@ module subroutine damage_init
 
   do ph = 1,phases%length
 
-    Nconstituents = count(material_phaseAt2 == ph)
+    Nmembers = count(material_phaseAt2 == ph)
 
-    allocate(current(ph)%phi(Nconstituents),source=1.0_pReal)
+    allocate(current(ph)%phi(Nmembers),source=1.0_pReal)
-    allocate(current(ph)%d_phi_d_dot_phi(Nconstituents),source=0.0_pReal)
+    allocate(current(ph)%d_phi_d_dot_phi(Nmembers),source=0.0_pReal)
 
     phase => phases%get(ph)
     sources => phase%get('damage',defaultVal=emptyList)

src/phase_damage_anisobrittle.f90

@@ -40,7 +40,7 @@ module function anisobrittle_init() result(mySources)
     phase, &
     sources, &
     src
-  integer :: Nconstituents,p
+  integer :: Nmembers,p
   integer, dimension(:), allocatable :: N_cl
   character(len=pStringLen) :: extmsg = ''
 
@@ -92,8 +92,8 @@ module function anisobrittle_init() result(mySources)
         if (any(prm%g_crit <  0.0_pReal)) extmsg = trim(extmsg)//' g_crit'
         if (any(prm%s_crit <  0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
 
-        Nconstituents = count(material_phaseAt==p) * discretization_nIPs
+        Nmembers = count(material_phaseAt==p) * discretization_nIPs
-        call phase_allocateState(damageState(p),Nconstituents,1,1,0)
+        call phase_allocateState(damageState(p),Nmembers,1,1,0)
         damageState(p)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
         if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'
 

src/phase_damage_anisoductile.f90

@@ -35,7 +35,7 @@ module function anisoductile_init() result(mySources)
     pl, &
     sources, &
     src
-  integer :: Ninstances,Nconstituents,p
+  integer :: Ninstances,Nmembers,p
   integer, dimension(:), allocatable :: N_sl
   character(len=pStringLen) :: extmsg = ''
 
@@ -78,8 +78,8 @@ module function anisoductile_init() result(mySources)
         if (prm%q              <= 0.0_pReal)  extmsg = trim(extmsg)//' q'
         if (any(prm%gamma_crit <  0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
 
-        Nconstituents=count(material_phaseAt2==p)
+        Nmembers=count(material_phaseAt2==p)
-        call phase_allocateState(damageState(p),Nconstituents,1,1,0)
+        call phase_allocateState(damageState(p),Nmembers,1,1,0)
         damageState(p)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
         if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'
 

src/phase_damage_isobrittle.f90

@@ -31,7 +31,7 @@ module function isobrittle_init() result(mySources)
     phase, &
     sources, &
     src
-  integer :: Nconstituents,p
+  integer :: Nmembers,p
   character(len=pStringLen) :: extmsg = ''
 
 
@@ -64,8 +64,8 @@ module function isobrittle_init() result(mySources)
         ! sanity checks
         if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit'
 
-        Nconstituents = count(material_phaseAt2==p)
+        Nmembers = count(material_phaseAt2==p)
-        call phase_allocateState(damageState(p),Nconstituents,1,1,1)
+        call phase_allocateState(damageState(p),Nmembers,1,1,1)
         damageState(p)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
         if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'
 

src/phase_damage_isoductile.f90

@@ -33,7 +33,7 @@ module function isoductile_init() result(mySources)
     phase, &
     sources, &
     src
-  integer :: Ninstances,Nconstituents,p
+  integer :: Ninstances,Nmembers,p
   character(len=pStringLen) :: extmsg = ''
 
 
@@ -68,8 +68,8 @@ module function isoductile_init() result(mySources)
         if (prm%q          <= 0.0_pReal) extmsg = trim(extmsg)//' q'
         if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
 
-        Nconstituents=count(material_phaseAt2==p)
+        Nmembers=count(material_phaseAt2==p)
-        call phase_allocateState(damageState(p),Nconstituents,1,1,0)
+        call phase_allocateState(damageState(p),Nmembers,1,1,0)
         damageState(p)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
         if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'
 

src/phase_mechanical.f90

@@ -191,7 +191,7 @@ module subroutine mechanical_init(materials,phases)
     ph, &
     me, &
     stiffDegradationCtr, &
-    Nconstituents
+    Nmembers
   class(tNode), pointer :: &
     num_crystallite, &
     material, &
@@ -229,22 +229,22 @@ module subroutine mechanical_init(materials,phases)
   allocate(material_orientation0(homogenization_maxNconstituents,phases%length,maxVal(material_phaseMemberAt)))
 
   do ph = 1, phases%length
-    Nconstituents = count(material_phaseAt == ph) * discretization_nIPs
+    Nmembers = count(material_phaseAt == ph) * discretization_nIPs
 
-    allocate(phase_mechanical_Fi(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Fe(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_Fe(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Fi0(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_Fi0(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Fp(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_Fp(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Fp0(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_Fp0(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Li(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Li0(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Lp0(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_Lp0(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Lp(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_S(ph)%data(3,3,Nconstituents),source=0.0_pReal)
+    allocate(phase_mechanical_S(ph)%data(3,3,Nmembers),source=0.0_pReal)
-    allocate(phase_mechanical_P(ph)%data(3,3,Nconstituents),source=0.0_pReal)
+    allocate(phase_mechanical_P(ph)%data(3,3,Nmembers),source=0.0_pReal)
-    allocate(phase_mechanical_S0(ph)%data(3,3,Nconstituents),source=0.0_pReal)
+    allocate(phase_mechanical_S0(ph)%data(3,3,Nmembers),source=0.0_pReal)
-    allocate(phase_mechanical_F(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_F(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_F0(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers))
 
     phase   => phases%get(ph)
     mech    => phase%get('mechanics')
@@ -279,7 +279,6 @@ module subroutine mechanical_init(materials,phases)
   endif
 
 
-  !$OMP PARALLEL DO PRIVATE(ph,me,material,constituents,constituent) 
   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))
@@ -305,7 +304,6 @@ module subroutine mechanical_init(materials,phases)
 
     enddo
   enddo; enddo
-  !$OMP END PARALLEL DO
 
 
 ! initialize plasticity

src/phase_mechanical_plastic_dislotungsten.f90

@@ -79,7 +79,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
   logical, dimension(:), allocatable :: myPlasticity
   integer :: &
     ph, i, &
-    Nconstituents, &
+    Nmembers, &
     sizeState, sizeDotState, &
     startIndex, endIndex
   integer,    dimension(:), allocatable :: &
@@ -220,18 +220,18 @@ module function plastic_dislotungsten_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-    Nconstituents = count(material_phaseAt2 == ph)
+    Nmembers = count(material_phaseAt2 == ph)
     sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
     sizeState = sizeDotState
 
-    call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0)
+    call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
 
 !--------------------------------------------------------------------------------------------------
 ! state aliases and initialization
     startIndex = 1
     endIndex   = prm%sum_N_sl
     stt%rho_mob => plasticState(ph)%state(startIndex:endIndex,:)
-    stt%rho_mob =  spread(rho_mob_0,2,Nconstituents)
+    stt%rho_mob =  spread(rho_mob_0,2,Nmembers)
     dot%rho_mob => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
     if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
@@ -239,7 +239,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_sl
     stt%rho_dip => plasticState(ph)%state(startIndex:endIndex,:)
-    stt%rho_dip =  spread(rho_dip_0,2,Nconstituents)
+    stt%rho_dip =  spread(rho_dip_0,2,Nmembers)
     dot%rho_dip => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
 
@@ -251,8 +251,8 @@ module function plastic_dislotungsten_init() result(myPlasticity)
     ! global alias
     plasticState(ph)%slipRate        => plasticState(ph)%dotState(startIndex:endIndex,:)
 
-    allocate(dst%Lambda_sl(prm%sum_N_sl,Nconstituents),         source=0.0_pReal)
+    allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers),         source=0.0_pReal)
-    allocate(dst%threshold_stress(prm%sum_N_sl,Nconstituents),  source=0.0_pReal)
+    allocate(dst%threshold_stress(prm%sum_N_sl,Nmembers),  source=0.0_pReal)
 
     plasticState(ph)%state0 = plasticState(ph)%state                                                ! ToDo: this could be done centrally
 

src/phase_mechanical_plastic_dislotwin.f90

@@ -127,7 +127,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
   logical, dimension(:), allocatable :: myPlasticity
   integer :: &
     ph, i, &
-    Nconstituents, &
+    Nmembers, &
     sizeState, sizeDotState, &
     startIndex, endIndex
   integer,     dimension(:), allocatable :: &
@@ -406,21 +406,21 @@ module function plastic_dislotwin_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-    Nconstituents  = count(material_phaseAt2 == ph)
+    Nmembers  = count(material_phaseAt2 == ph)
     sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl &
                  + size(['f_tw'])                           * prm%sum_N_tw &
                  + size(['f_tr'])                           * prm%sum_N_tr
     sizeState = sizeDotState
 
 
-    call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0)
+    call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
 
 !--------------------------------------------------------------------------------------------------
 ! locally defined state aliases and initialization of state0 and atol
     startIndex = 1
     endIndex   = prm%sum_N_sl
     stt%rho_mob=>plasticState(ph)%state(startIndex:endIndex,:)
-    stt%rho_mob= spread(rho_mob_0,2,Nconstituents)
+    stt%rho_mob= spread(rho_mob_0,2,Nmembers)
     dot%rho_mob=>plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
     if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
@@ -428,7 +428,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_sl
     stt%rho_dip=>plasticState(ph)%state(startIndex:endIndex,:)
-    stt%rho_dip= spread(rho_dip_0,2,Nconstituents)
+    stt%rho_dip= spread(rho_dip_0,2,Nmembers)
     dot%rho_dip=>plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
 
@@ -454,18 +454,18 @@ module function plastic_dislotwin_init() result(myPlasticity)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tr',defaultVal=1.0e-6_pReal)
     if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tr'
 
-    allocate(dst%Lambda_sl             (prm%sum_N_sl,Nconstituents),source=0.0_pReal)
+    allocate(dst%Lambda_sl             (prm%sum_N_sl,Nmembers),source=0.0_pReal)
-    allocate(dst%tau_pass              (prm%sum_N_sl,Nconstituents),source=0.0_pReal)
+    allocate(dst%tau_pass              (prm%sum_N_sl,Nmembers),source=0.0_pReal)
 
-    allocate(dst%Lambda_tw             (prm%sum_N_tw,Nconstituents),source=0.0_pReal)
+    allocate(dst%Lambda_tw             (prm%sum_N_tw,Nmembers),source=0.0_pReal)
-    allocate(dst%tau_hat_tw            (prm%sum_N_tw,Nconstituents),source=0.0_pReal)
+    allocate(dst%tau_hat_tw            (prm%sum_N_tw,Nmembers),source=0.0_pReal)
-    allocate(dst%tau_r_tw              (prm%sum_N_tw,Nconstituents),source=0.0_pReal)
+    allocate(dst%tau_r_tw              (prm%sum_N_tw,Nmembers),source=0.0_pReal)
-    allocate(dst%V_tw                  (prm%sum_N_tw,Nconstituents),source=0.0_pReal)
+    allocate(dst%V_tw                  (prm%sum_N_tw,Nmembers),source=0.0_pReal)
 
-    allocate(dst%Lambda_tr             (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
+    allocate(dst%Lambda_tr             (prm%sum_N_tr,Nmembers),source=0.0_pReal)
-    allocate(dst%tau_hat_tr            (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
+    allocate(dst%tau_hat_tr            (prm%sum_N_tr,Nmembers),source=0.0_pReal)
-    allocate(dst%tau_r_tr              (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
+    allocate(dst%tau_r_tr              (prm%sum_N_tr,Nmembers),source=0.0_pReal)
-    allocate(dst%V_tr                  (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
+    allocate(dst%V_tr                  (prm%sum_N_tr,Nmembers),source=0.0_pReal)
 
     plasticState(ph)%state0 = plasticState(ph)%state                                                ! ToDo: this could be done centrally
 

src/phase_mechanical_plastic_isotropic.f90

@@ -52,7 +52,7 @@ module function plastic_isotropic_init() result(myPlasticity)
   logical, dimension(:), allocatable :: myPlasticity
   integer :: &
     ph, &
-    Nconstituents, &
+    Nmembers, &
     sizeState, sizeDotState
   real(pReal) :: &
     xi_0                                                                                            !< initial critical stress
@@ -119,11 +119,11 @@ module function plastic_isotropic_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-    Nconstituents = count(material_phaseAt2 == ph)
+    Nmembers = count(material_phaseAt2 == ph)
     sizeDotState = size(['xi   ','gamma'])
     sizeState = sizeDotState
 
-    call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0)
+    call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
 
 !--------------------------------------------------------------------------------------------------
 ! state aliases and initialization

src/phase_mechanical_plastic_kinehardening.f90

@@ -62,7 +62,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
   logical, dimension(:), allocatable :: myPlasticity
   integer :: &
     ph, o,  &
-    Nconstituents, &
+    Nmembers, &
     sizeState, sizeDeltaState, sizeDotState, &
     startIndex, endIndex
   integer,     dimension(:), allocatable :: &
@@ -165,19 +165,19 @@ module function plastic_kinehardening_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-    Nconstituents = count(material_phaseAt2 == ph)
+    Nmembers = count(material_phaseAt2 == ph)
     sizeDotState   = size(['crss     ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml
     sizeDeltaState = size(['sense ',   'chi0  ',    'gamma0'   ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml
     sizeState = sizeDotState + sizeDeltaState
 
-    call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
+    call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
 
 !--------------------------------------------------------------------------------------------------
 ! state aliases and initialization
     startIndex = 1
     endIndex   = prm%sum_N_sl
     stt%crss => plasticState(ph)%state   (startIndex:endIndex,:)
-    stt%crss = spread(xi_0, 2, Nconstituents)
+    stt%crss = spread(xi_0, 2, Nmembers)
     dot%crss => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
     if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'

src/phase_mechanical_plastic_nonlocal.f90

@@ -177,7 +177,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
   integer :: &
     Ninstances, &
     ph, &
-    Nconstituents, &
+    Nmembers, &
     sizeState, sizeDotState, sizeDependentState, sizeDeltaState, &
     s1, s2, &
     s, t, l
@@ -398,7 +398,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-    Nconstituents  = count(material_phaseAt2 == ph)
+    Nmembers  = count(material_phaseAt2 == ph)
     sizeDotState = size([   'rhoSglEdgePosMobile   ','rhoSglEdgeNegMobile   ', &
                             'rhoSglScrewPosMobile  ','rhoSglScrewNegMobile  ', &
                             'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', &
@@ -412,9 +412,9 @@ module function plastic_nonlocal_init() result(myPlasticity)
                                 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl     !< other dependent state variables that are not updated by microstructure
     sizeDeltaState            = sizeDotState
 
-    call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
+    call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
 
-    allocate(geom(ph)%V_0(Nconstituents))
+    allocate(geom(ph)%V_0(Nmembers))
     call storeGeometry(ph)
 
     plasticState(ph)%nonlocal         = pl%get_asBool('nonlocal')
@@ -486,26 +486,26 @@ module function plastic_nonlocal_init() result(myPlasticity)
         dot%rho_dip_scr => plasticState(ph)%dotState               (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
         del%rho_dip_scr => plasticState(ph)%deltaState             (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
 
-    stt%gamma => plasticState(ph)%state                      (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
+    stt%gamma => plasticState(ph)%state                      (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers)
-    dot%gamma => plasticState(ph)%dotState                   (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
+    dot%gamma => plasticState(ph)%dotState                   (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers)
-    del%gamma => plasticState(ph)%deltaState                 (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
+    del%gamma => plasticState(ph)%deltaState                 (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers)
     plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl )  = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-2_pReal)
     if(any(plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) &
       extmsg = trim(extmsg)//' atol_gamma'
-    plasticState(ph)%slipRate => plasticState(ph)%dotState    (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
+    plasticState(ph)%slipRate => plasticState(ph)%dotState    (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers)
 
-    stt%rho_forest => plasticState(ph)%state                 (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nconstituents)
+    stt%rho_forest => plasticState(ph)%state                 (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nmembers)
-    stt%v          => plasticState(ph)%state                 (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents)
+    stt%v          => plasticState(ph)%state                 (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers)
-        stt%v_edg_pos  => plasticState(ph)%state             (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nconstituents)
+        stt%v_edg_pos  => plasticState(ph)%state             (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nmembers)
-        stt%v_edg_neg  => plasticState(ph)%state             (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nconstituents)
+        stt%v_edg_neg  => plasticState(ph)%state             (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nmembers)
-        stt%v_scr_pos  => plasticState(ph)%state             (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nconstituents)
+        stt%v_scr_pos  => plasticState(ph)%state             (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nmembers)
-        stt%v_scr_neg  => plasticState(ph)%state             (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents)
+        stt%v_scr_neg  => plasticState(ph)%state             (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers)
 
-    allocate(dst%tau_pass(prm%sum_N_sl,Nconstituents),source=0.0_pReal)
+    allocate(dst%tau_pass(prm%sum_N_sl,Nmembers),source=0.0_pReal)
-    allocate(dst%tau_back(prm%sum_N_sl,Nconstituents),source=0.0_pReal)
+    allocate(dst%tau_back(prm%sum_N_sl,Nmembers),source=0.0_pReal)
     end associate
 
-    if (Nconstituents > 0) call stateInit(ini,ph,Nconstituents)
+    if (Nmembers > 0) call stateInit(ini,ph,Nmembers)
     plasticState(ph)%state0 = plasticState(ph)%state
 
 !--------------------------------------------------------------------------------------------------
@@ -527,7 +527,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
     if(.not. myPlasticity(ph)) cycle
 
     phase => phases%get(ph)
-    Nconstituents = count(material_phaseAt2 == ph)
+    Nmembers = count(material_phaseAt2 == ph)
     l = 0
     do t = 1,4
       do s = 1,param(ph)%sum_N_sl
@@ -1579,13 +1579,13 @@ end subroutine plastic_nonlocal_results
 !--------------------------------------------------------------------------------------------------
 !> @brief populates the initial dislocation density
 !--------------------------------------------------------------------------------------------------
-subroutine stateInit(ini,phase,Nconstituents)
+subroutine stateInit(ini,phase,Nmembers)
 
   type(tInitialParameters) :: &
     ini
   integer,intent(in) :: &
     phase, &
-    Nconstituents
+    Nmembers
   integer :: &
     i, &
     e, &
@@ -1602,7 +1602,7 @@ subroutine stateInit(ini,phase,Nconstituents)
     totalVolume, &
     densityBinning, &
     minimumIpVolume
-  real(pReal), dimension(Nconstituents) :: &
+  real(pReal), dimension(Nmembers) :: &
     volume
 
 
@@ -1622,13 +1622,13 @@ subroutine stateInit(ini,phase,Nconstituents)
     meanDensity = 0.0_pReal
     do while(meanDensity < ini%random_rho_u)
       call random_number(rnd)
-      phasemember = nint(rnd(1)*real(Nconstituents,pReal) + 0.5_pReal)
+      phasemember = nint(rnd(1)*real(Nmembers,pReal) + 0.5_pReal)
       s           = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal)
       meanDensity = meanDensity + densityBinning * volume(phasemember) / totalVolume
       stt%rhoSglMobile(s,phasemember) = densityBinning
     enddo
   else                                ! homogeneous distribution with noise
-    do e = 1, Nconstituents
+    do e = 1, Nmembers
       do f = 1,size(ini%N_sl,1)
         from = 1 + sum(ini%N_sl(1:f-1))
         upto = sum(ini%N_sl(1:f))
@@ -1822,16 +1822,13 @@ subroutine storeGeometry(ph)
   integer, intent(in) :: ph
 
   integer :: ip, el, ce, co
+  real(pReal), dimension(:), allocatable :: V
 
-  ce = 0
+
-  do el = 1, size(material_homogenizationMemberAt,2)
+  V = reshape(IPvolume,[product(shape(IPvolume))])
-    do ip = 1, size(material_homogenizationMemberAt,1)
+  do ce = 1, size(material_homogenizationMemberAt2,1)
-      ce = ce + 1
     do co = 1, homogenization_maxNconstituents
-        if(material_phaseAt2(co,ce) == ph) then
+      if (material_phaseAt2(co,ce) == ph) geom(ph)%V_0(material_phaseMemberAt2(co,ce)) = V(ce)
-           geom(ph)%V_0(material_phaseMemberAt2(co,ce)) = IPvolume(ip,el)
-        endif
-      enddo
     enddo
   enddo
 

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -71,7 +71,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
   logical, dimension(:), allocatable :: myPlasticity
   integer :: &
     ph, i, &
-    Nconstituents, &
+    Nmembers, &
     sizeState, sizeDotState, &
     startIndex, endIndex
   integer,     dimension(:), allocatable :: &
@@ -223,20 +223,20 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-    Nconstituents = count(material_phaseAt2 == ph)
+    Nmembers = count(material_phaseAt2 == ph)
     sizeDotState = size(['xi_sl   ','gamma_sl']) * prm%sum_N_sl &
                  + size(['xi_tw   ','gamma_tw']) * prm%sum_N_tw
     sizeState = sizeDotState
 
 
-    call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0)
+    call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
 
 !--------------------------------------------------------------------------------------------------
 ! state aliases and initialization
     startIndex = 1
     endIndex   = prm%sum_N_sl
     stt%xi_slip => plasticState(ph)%state   (startIndex:endIndex,:)
-    stt%xi_slip =  spread(xi_0_sl, 2, Nconstituents)
+    stt%xi_slip =  spread(xi_0_sl, 2, Nmembers)
     dot%xi_slip => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
     if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
@@ -244,7 +244,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_tw
     stt%xi_twin => plasticState(ph)%state   (startIndex:endIndex,:)
-    stt%xi_twin =  spread(xi_0_tw, 2, Nconstituents)
+    stt%xi_twin =  spread(xi_0_tw, 2, Nmembers)
     dot%xi_twin => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
     if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'

src/phase_thermal.f90

@@ -15,13 +15,13 @@ submodule(phase) thermal
     THERMAL_EXTERNALHEAT_ID
   end enum
 
-  type :: tDataContainer
+  type :: tDataContainer             ! ?? not very telling name. Better: "fieldQuantities" ??
     real(pReal), dimension(:), allocatable :: T, dot_T
   end type tDataContainer
   integer(kind(THERMAL_UNDEFINED_ID)),  dimension(:,:), allocatable :: &
     thermal_source
 
-  type(tDataContainer), dimension(:), allocatable :: current
+  type(tDataContainer), dimension(:), allocatable :: current          ! ?? not very telling name. Better: "field" ??
 
   integer :: thermal_source_maxSizeDotState
 
@@ -78,7 +78,7 @@ module subroutine thermal_init(phases)
 
   integer :: &
     ph, so, &
-    Nconstituents
+    Nmembers
 
 
   print'(/,a)', ' <<<+-  phase:thermal init  -+>>>'
@@ -89,18 +89,13 @@ module subroutine thermal_init(phases)
   allocate(thermal_Nsources(phases%length),source = 0)
 
   do ph = 1, phases%length
-
+    Nmembers = count(material_phaseAt2 == ph)
-    Nconstituents = count(material_phaseAt2 == ph)
+    allocate(current(ph)%T(Nmembers),source=300.0_pReal)
-
+    allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal)
-    allocate(current(ph)%T(Nconstituents),source=300.0_pReal)
-    allocate(current(ph)%dot_T(Nconstituents),source=0.0_pReal)
     phase => phases%get(ph)
-    if(phase%contains('thermal')) then
+    thermal => phase%get('thermal',defaultVal=emptyDict)
-      thermal => phase%get('thermal')
     sources => thermal%get('source',defaultVal=emptyList)
-
     thermal_Nsources(ph) = sources%length
-    endif
     allocate(thermalstate(ph)%p(thermal_Nsources(ph)))
   enddo
 
@@ -112,7 +107,7 @@ module subroutine thermal_init(phases)
   endif
 
   thermal_source_maxSizeDotState = 0
-  PhaseLoop2:do ph = 1,phases%length
+  do ph = 1,phases%length
 
     do so = 1,thermal_Nsources(ph)
       thermalState(ph)%p(so)%state  = thermalState(ph)%p(so)%state0
@@ -120,7 +115,7 @@ module subroutine thermal_init(phases)
 
     thermal_source_maxSizeDotState  = max(thermal_source_maxSizeDotState, &
                                           maxval(thermalState(ph)%p%sizeDotState))
-  enddo PhaseLoop2
+  enddo
 
 end subroutine thermal_init
 
@@ -156,7 +151,6 @@ module subroutine phase_thermal_getRate(TDot, ph,me)
    Tdot = Tdot + my_Tdot
   enddo
 
-
 end subroutine phase_thermal_getRate
 
 
@@ -185,7 +179,7 @@ function phase_thermal_collectDotState(ph,me) result(broken)
 end function phase_thermal_collectDotState
 
 
-module function thermal_stress(Delta_t,ph,me) result(converged_)
+module function thermal_stress(Delta_t,ph,me) result(converged_)           ! ?? why is this called "stress" when it seems closer to "updateState" ??
 
   real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: ph, me
@@ -301,14 +295,12 @@ function thermal_active(source_label,src_length)  result(active_source)
   allocate(active_source(src_length,phases%length), source = .false. )
   do p = 1, phases%length
     phase => phases%get(p)
-    if (phase%contains('thermal')) then
+    thermal =>  phase%get('thermal',defaultVal=emptyDict)
-      thermal =>  phase%get('thermal',defaultVal=emptyList)
     sources =>  thermal%get('source',defaultVal=emptyList)
     do s = 1, sources%length
       src => sources%get(s)
-        if(src%get_asString('type') == source_label) active_source(s,p) = .true.
+      active_source(s,p) = src%get_asString('type') == source_label
     enddo
-    endif
   enddo
 
 

src/phase_thermal_dissipation.f90

@@ -31,7 +31,7 @@ module function dissipation_init(source_length) result(mySources)
     phase, &
     sources, thermal, &
     src
-  integer :: so,Nconstituents,ph
+  integer :: so,Nmembers,ph
 
 
   mySources = thermal_active('dissipation',source_length)
@@ -54,8 +54,8 @@ module function dissipation_init(source_length) result(mySources)
           src => sources%get(so)
 
           prm%kappa = src%get_asFloat('kappa')
-          Nconstituents = count(material_phaseAt2 == ph)
+          Nmembers = count(material_phaseAt2 == ph)
-          call phase_allocateState(thermalState(ph)%p(so),Nconstituents,0,0,0)
+          call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0)
 
         end associate
       endif

src/phase_thermal_externalheat.f90

@@ -38,7 +38,7 @@ module function externalheat_init(source_length) result(mySources)
     phase, &
     sources, thermal, &
     src
-  integer :: so,Nconstituents,ph
+  integer :: so,Nmembers,ph
 
 
   mySources = thermal_active('externalheat',source_length)
@@ -67,8 +67,8 @@ module function externalheat_init(source_length) result(mySources)
 
           prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n))
 
-          Nconstituents = count(material_phaseAt2 == ph)
+          Nmembers = count(material_phaseAt2 == ph)
-          call phase_allocateState(thermalState(ph)%p(so),Nconstituents,1,1,0)
+          call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0)
         end associate
       endif
     enddo