diff --git a/PRIVATE b/PRIVATE index d96bfb329..3a2f89547 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit d96bfb32920c96a8a43958f76a209d34c6bd841a +Subproject commit 3a2f89547c264044a7bfab9d33aee78eec495a76 diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py index caf09d536..fa47805bb 100755 --- a/processing/post/DADF5_postResults.py +++ b/processing/post/DADF5_postResults.py @@ -22,10 +22,15 @@ parser.add_argument('filenames', nargs='+', help='DADF5 files') parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string', help='name of subdirectory to hold output') +parser.add_argument('--mat', nargs='+', + help='labels for materialpoint/homogenization',dest='mat') +parser.add_argument('--con', nargs='+', + help='labels for constituent/crystallite/constitutive',dest='con') options = parser.parse_args() -options.labels = ['Fe','Fp','xi_sl'] +if options.mat is None: options.mat=[] +if options.con is None: options.con=[] # --- loop over input files ------------------------------------------------------------------------ @@ -48,16 +53,13 @@ for filename in options.filenames: data = np.array([inc['inc'] for j in range(np.product(results.grid))]).reshape([np.product(results.grid),1]) header+= 'inc' - - data = np.concatenate((data,np.array([j+1 for j in range(np.product(results.grid))]).reshape([np.product(results.grid),1])),1) - header+=' node' coords = coords.reshape([np.product(results.grid),3]) data = np.concatenate((data,coords),1) header+=' 1_pos 2_pos 3_pos' results.active['increments'] = [inc] - for label in options.labels: + for label in options.con: for o in results.c_output_types: results.active['c_output_types'] = [o] for c in results.constituents: @@ -67,12 +69,33 @@ for filename in options.filenames: continue label = x[0].split('/')[-1] array = results.read_dataset(x,0) - d = np.product(np.shape(array)[1:]) + d = int(np.product(np.shape(array)[1:])) array = np.reshape(array,[np.product(results.grid),d]) data = np.concatenate((data,array),1) - header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)]) - + if d>1: + header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)]) + else: + header+=' '+label + + for label in options.mat: + for o in results.m_output_types: + results.active['m_output_types'] = [o] + for m in results.materialpoints: + results.active['materialpoints'] = [m] + x = results.get_dataset_location(label) + if len(x) == 0: + continue + label = x[0].split('/')[-1] + array = results.read_dataset(x,0) + d = int(np.product(np.shape(array)[1:])) + array = np.reshape(array,[np.product(results.grid),d]) + data = np.concatenate((data,array),1) + + if d>1: + header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)]) + else: + header+=' '+label dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) try: diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py index dc1177488..ef27a973c 100755 --- a/processing/post/DADF5_vtk_cells.py +++ b/processing/post/DADF5_vtk_cells.py @@ -23,10 +23,15 @@ parser.add_argument('filenames', nargs='+', help='DADF5 files') parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string', help='name of subdirectory to hold output') +parser.add_argument('--mat', nargs='+', + help='labels for materialpoint/homogenization',dest='mat') +parser.add_argument('--con', nargs='+', + help='labels for constituent/crystallite/constitutive',dest='con') options = parser.parse_args() -options.labels = ['Fe','Fp','xi_sl'] +if options.mat is None: options.mat=[] +if options.con is None: options.con=[] # --- loop over input files ------------------------------------------------------------------------ @@ -54,7 +59,9 @@ for filename in options.filenames: print('Output step {}/{}'.format(i+1,len(results.increments))) vtk_data = [] results.active['increments'] = [inc] - for label in options.labels: + + for label in options.con: + for o in results.c_output_types: results.active['c_output_types'] = [o] if o != 'generic': diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 841029af8..f5fc81b16 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -48,16 +48,24 @@ class DADF5(): for o in f['inc{:05}/constituent/{}'.format(self.increments[0]['inc'],c)].keys(): self.c_output_types.append(o) self.c_output_types = list(set(self.c_output_types)) # make unique + + self.m_output_types = [] + for m in self.materialpoints: + for o in f['inc{:05}/materialpoint/{}'.format(self.increments[0]['inc'],m)].keys(): + self.m_output_types.append(o) + self.m_output_types = list(set(self.m_output_types)) # make unique self.active= {'increments': self.increments, 'constituents': self.constituents, 'materialpoints': self.materialpoints, 'constituent': self.Nconstituents, - 'c_output_types': self.c_output_types} + 'c_output_types': self.c_output_types, + 'm_output_types': self.m_output_types} self.filename = filename self.mode = mode - + + def list_data(self): """Shows information on all datasets in the file""" with h5py.File(self.filename,'r') as f: @@ -73,6 +81,16 @@ class DADF5(): print(' {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode())) except: pass + for m in self.active['materialpoints']: + group_materialpoint = group_inc+'/materialpoint/'+m + for t in self.active['m_output_types']: + print(' {}'.format(t)) + group_output_types = group_materialpoint+'/'+t + try: + for x in f[group_output_types].keys(): + print(' {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode())) + except: + pass def get_dataset_location(self,label): @@ -81,14 +99,25 @@ class DADF5(): with h5py.File(self.filename,'r') as f: for i in self.active['increments']: group_inc = 'inc{:05}'.format(i['inc']) + for c in self.active['constituents']: group_constituent = group_inc+'/constituent/'+c for t in self.active['c_output_types']: try: f[group_constituent+'/'+t+'/'+label] path.append(group_constituent+'/'+t+'/'+label) - except: - pass + except Exception as e: + print('unable to locate constituents dataset: '+ str(e)) + + for m in self.active['materialpoints']: + group_materialpoint = group_inc+'/materialpoint/'+m + for t in self.active['m_output_types']: + try: + f[group_materialpoint+'/'+t+'/'+label] + path.append(group_materialpoint+'/'+t+'/'+label) + except Exception as e: + print('unable to locate materialpoints dataset: '+ str(e)) + return path @@ -100,13 +129,29 @@ class DADF5(): """ with h5py.File(self.filename,'r') as f: shape = (self.Nmaterialpoints,) + np.shape(f[path[0]])[1:] + if len(shape) == 1: shape = shape +(1,) dataset = np.full(shape,np.nan) for pa in path: label = pa.split('/')[2] - p = np.where(f['mapping/cellResults/constituent'][:,c]['Name'] == str.encode(label))[0] - u = (f['mapping/cellResults/constituent'][p,c]['Position']) - dataset[p,:] = f[pa][u,:] - + try: + p = np.where(f['mapping/cellResults/constituent'][:,c]['Name'] == str.encode(label))[0] + u = (f['mapping/cellResults/constituent'][p,c]['Position']) + a = np.array(f[pa]) + if len(a.shape) == 1: + a=a.reshape([a.shape[0],1]) + dataset[p,:] = a[u,:] + except Exception as e: + print('unable to read constituent: '+ str(e)) + try: + p = np.where(f['mapping/cellResults/materialpoint']['Name'] == str.encode(label))[0] + u = (f['mapping/cellResults/materialpoint'][p.tolist()]['Position']) + a = np.array(f[pa]) + if len(a.shape) == 1: + a=a.reshape([a.shape[0],1]) + dataset[p,:] = a[u,:] + except Exception as e: + print('unable to read materialpoint: '+ str(e)) + return dataset diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 0ae1323f0..5e1cd71eb 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -46,7 +46,6 @@ #include "plastic_nonlocal.f90" #include "constitutive.f90" #include "crystallite.f90" -#include "homogenization_mech_RGC.f90" #include "thermal_isothermal.f90" #include "thermal_adiabatic.f90" #include "thermal_conduction.f90" @@ -56,4 +55,5 @@ #include "homogenization.f90" #include "homogenization_mech_none.f90" #include "homogenization_mech_isostrain.f90" +#include "homogenization_mech_RGC.f90" #include "CPFEM.f90" diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 605cf5094..3210f02d4 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -6,7 +6,20 @@ !-------------------------------------------------------------------------------------------------- module homogenization use prec + use IO + use config + use debug + use math use material + use numerics + use constitutive + use crystallite + use mesh + use FEsolving +#if defined(PETSc) || defined(DAMASK_HDF5) + use results + use HDF5_utilities +#endif !-------------------------------------------------------------------------------------------------- ! General variables for the homogenization at a material point @@ -22,7 +35,6 @@ module homogenization materialpoint_results !< results array of material point integer, public, protected :: & materialpoint_sizeResults, & - homogenization_maxSizePostResults, & thermal_maxSizePostResults, & damage_maxSizePostResults @@ -47,11 +59,24 @@ module homogenization module subroutine mech_isostrain_init end subroutine mech_isostrain_init + module subroutine mech_RGC_init + end subroutine mech_RGC_init + + module subroutine mech_isostrain_partitionDeformation(F,avgF) real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point end subroutine mech_isostrain_partitionDeformation + module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of) + real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient + real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point + integer, intent(in) :: & + instance, & + of + end subroutine mech_RGC_partitionDeformation + + module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point @@ -60,7 +85,37 @@ module homogenization real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses integer, intent(in) :: instance end subroutine mech_isostrain_averageStressAndItsTangent + + module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) + real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point + real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point + + real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses + real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses + integer, intent(in) :: instance + end subroutine mech_RGC_averageStressAndItsTangent + + + module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) + logical, dimension(2) :: mech_RGC_updateState + real(pReal), dimension(:,:,:), intent(in) :: & + P,& !< partitioned stresses + F,& !< partitioned deformation gradients + F0 !< partitioned initial deformation gradients + real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses + real(pReal), dimension(3,3), intent(in) :: avgF !< average F + real(pReal), intent(in) :: dt !< time increment + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + end function mech_RGC_updateState + + module subroutine mech_RGC_results(instance,group) + integer, intent(in) :: instance !< homogenization instance + character(len=*), intent(in) :: group !< group name in HDF5 file + end subroutine mech_RGC_results + end interface public :: & @@ -68,11 +123,6 @@ module homogenization materialpoint_stressAndItsTangent, & materialpoint_postResults, & homogenization_results - private :: & - partitionDeformation, & - updateState, & - averageStressAndItsTangent, & - postResults contains @@ -81,36 +131,12 @@ contains !> @brief module initialization !-------------------------------------------------------------------------------------------------- subroutine homogenization_init - use math, only: & - math_I3 - use debug, only: & - debug_level, & - debug_homogenization, & - debug_levelBasic, & - debug_e, & - debug_g - use mesh, only: & - theMesh, & - mesh_element - use constitutive, only: & - constitutive_plasticity_maxSizePostResults, & - constitutive_source_maxSizePostResults - use crystallite, only: & - crystallite_maxSizePostResults - use config, only: & - config_deallocate, & - config_homogenization, & - homogenization_name - use homogenization_mech_RGC use thermal_isothermal use thermal_adiabatic use thermal_conduction use damage_none use damage_local use damage_nonlocal - use IO - use numerics, only: & - worldrank integer, parameter :: FILEUNIT = 200 integer :: e,i,p @@ -120,10 +146,9 @@ subroutine homogenization_init character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready logical :: valid - if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init - if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call homogenization_RGC_init + if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mech_RGC_init if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init if (any(thermal_type == THERMAL_adiabatic_ID)) call thermal_adiabatic_init @@ -135,39 +160,14 @@ subroutine homogenization_init !-------------------------------------------------------------------------------------------------- ! write description file for homogenization output - mainProcess2: if (worldrank == 0) then + mainProcess: if (worldrank == 0) then call IO_write_jobFile(FILEUNIT,'outputHomogenization') do p = 1,size(config_homogenization) if (any(material_homogenizationAt == p)) then - i = homogenization_typeInstance(p) ! which instance of this homogenization type - valid = .true. ! assume valid - select case(homogenization_type(p)) ! split per homogenization type - case (HOMOGENIZATION_NONE_ID) - outputName = HOMOGENIZATION_NONE_label - thisOutput => null() - thisSize => null() - case (HOMOGENIZATION_ISOSTRAIN_ID) - outputName = HOMOGENIZATION_ISOSTRAIN_label - thisOutput => null() - thisSize => null() - case (HOMOGENIZATION_RGC_ID) - outputName = HOMOGENIZATION_RGC_label - thisOutput => homogenization_RGC_output - thisSize => homogenization_RGC_sizePostResult - case default - valid = .false. - end select write(FILEUNIT,'(/,a,/)') '['//trim(homogenization_name(p))//']' - if (valid) then - write(FILEUNIT,'(a)') '(type)'//char(9)//trim(outputName) - write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p) - if (homogenization_type(p) /= HOMOGENIZATION_NONE_ID .and. & - homogenization_type(p) /= HOMOGENIZATION_ISOSTRAIN_ID) then - do e = 1,size(thisOutput(:,i)) - write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) - enddo - endif - endif + write(FILEUNIT,'(a)') '(type) n/a' + write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p) + i = thermal_typeInstance(p) ! which instance of this thermal type valid = .true. ! assume valid select case(thermal_type(p)) ! split per thermal type @@ -197,6 +197,7 @@ subroutine homogenization_init enddo endif endif + i = damage_typeInstance(p) ! which instance of this damage type valid = .true. ! assume valid select case(damage_type(p)) ! split per damage type @@ -229,7 +230,7 @@ subroutine homogenization_init endif enddo close(FILEUNIT) - endif mainProcess2 + endif mainProcess call config_deallocate('material.config/homogenization') @@ -237,7 +238,7 @@ subroutine homogenization_init ! allocate and initialize global variables allocate(materialpoint_dPdF(3,3,3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) allocate(materialpoint_F0(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) - materialpoint_F0 = spread(spread(math_I3,3,theMesh%elem%nIPs),4,theMesh%nElems) ! initialize to identity + materialpoint_F0 = spread(spread(math_I3,3,theMesh%elem%nIPs),4,theMesh%nElems) ! initialize to identity allocate(materialpoint_F(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) materialpoint_F = materialpoint_F0 ! initialize to identity allocate(materialpoint_subF0(3,3,theMesh%elem%nIPs,theMesh%nElems), source=0.0_pReal) @@ -252,18 +253,15 @@ subroutine homogenization_init !-------------------------------------------------------------------------------------------------- ! allocate and initialize global state and postresutls variables - homogenization_maxSizePostResults = 0 thermal_maxSizePostResults = 0 damage_maxSizePostResults = 0 do p = 1,size(config_homogenization) - homogenization_maxSizePostResults = max(homogenization_maxSizePostResults,homogState (p)%sizePostResults) thermal_maxSizePostResults = max(thermal_maxSizePostResults, thermalState (p)%sizePostResults) damage_maxSizePostResults = max(damage_maxSizePostResults ,damageState (p)%sizePostResults) enddo materialpoint_sizeResults = 1 & ! grain count - + 1 + homogenization_maxSizePostResults & ! homogSize & homogResult - + thermal_maxSizePostResults & + + 1 + thermal_maxSizePostResults & + damage_maxSizePostResults & + homogenization_maxNgrains * (1 + crystallite_maxSizePostResults & ! crystallite size & crystallite results + 1 + constitutive_plasticity_maxSizePostResults & ! constitutive size & constitutive results @@ -273,11 +271,6 @@ subroutine homogenization_init write(6,'(/,a)') ' <<<+- homogenization init -+>>>' if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then -#ifdef TODO - write(6,'(a32,1x,7(i8,1x))') 'homogenization_state0: ', shape(homogenization_state0) - write(6,'(a32,1x,7(i8,1x))') 'homogenization_subState0: ', shape(homogenization_subState0) - write(6,'(a32,1x,7(i8,1x))') 'homogenization_state: ', shape(homogenization_state) -#endif write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F0: ', shape(materialpoint_F0) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F: ', shape(materialpoint_F) @@ -303,51 +296,6 @@ end subroutine homogenization_init !> @brief parallelized calculation of stress and corresponding tangent at material points !-------------------------------------------------------------------------------------------------- subroutine materialpoint_stressAndItsTangent(updateJaco,dt) - use numerics, only: & - subStepMinHomog, & - subStepSizeHomog, & - stepIncreaseHomog, & - nMPstate - use FEsolving, only: & - FEsolving_execElem, & - FEsolving_execIP, & - terminallyIll - use mesh, only: & - mesh_element - use crystallite, only: & - crystallite_F0, & - crystallite_Fp0, & - crystallite_Fp, & - crystallite_Fi0, & - crystallite_Fi, & - crystallite_Lp0, & - crystallite_Lp, & - crystallite_Li0, & - crystallite_Li, & - crystallite_S0, & - crystallite_S, & - crystallite_partionedF0, & - crystallite_partionedF, & - crystallite_partionedFp0, & - crystallite_partionedLp0, & - crystallite_partionedFi0, & - crystallite_partionedLi0, & - crystallite_partionedS0, & - crystallite_dt, & - crystallite_requested, & - crystallite_stress, & - crystallite_stressTangent, & - crystallite_orientations -#ifdef DEBUG - use debug, only: & - debug_level, & - debug_homogenization, & - debug_levelBasic, & - debug_levelExtensive, & - debug_levelSelective, & - debug_e, & - debug_i -#endif real(pReal), intent(in) :: dt !< time increment logical, intent(in) :: updateJaco !< initiating Jacobian update @@ -640,14 +588,6 @@ end subroutine materialpoint_stressAndItsTangent !> @brief parallelized calculation of result array at material points !-------------------------------------------------------------------------------------------------- subroutine materialpoint_postResults - use FEsolving, only: & - FEsolving_execElem, & - FEsolving_execIP - use mesh, only: & - mesh_element - use crystallite, only: & - crystallite_sizePostResults, & - crystallite_postResults integer :: & thePos, & @@ -697,12 +637,6 @@ end subroutine materialpoint_postResults !> @brief partition material point def grad onto constituents !-------------------------------------------------------------------------------------------------- subroutine partitionDeformation(ip,el) - use mesh, only: & - mesh_element - use crystallite, only: & - crystallite_partionedF - use homogenization_mech_RGC, only: & - homogenization_RGC_partitionDeformation integer, intent(in) :: & ip, & !< integration point @@ -719,7 +653,7 @@ subroutine partitionDeformation(ip,el) materialpoint_subF(1:3,1:3,ip,el)) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - call homogenization_RGC_partitionDeformation(& + call mech_RGC_partitionDeformation(& crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & materialpoint_subF(1:3,1:3,ip,el),& ip, & @@ -734,15 +668,6 @@ end subroutine partitionDeformation !> "happy" with result !-------------------------------------------------------------------------------------------------- function updateState(ip,el) - use mesh, only: & - mesh_element - use crystallite, only: & - crystallite_P, & - crystallite_dPdF, & - crystallite_partionedF,& - crystallite_partionedF0 - use homogenization_mech_RGC, only: & - homogenization_RGC_updateState use thermal_adiabatic, only: & thermal_adiabatic_updateState use damage_local, only: & @@ -758,14 +683,14 @@ function updateState(ip,el) case (HOMOGENIZATION_RGC_ID) chosenHomogenization updateState = & updateState .and. & - homogenization_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & - crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & - crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el),& - materialpoint_subF(1:3,1:3,ip,el),& - materialpoint_subdt(ip,el), & - crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & - ip, & - el) + mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & + crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & + crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el),& + materialpoint_subF(1:3,1:3,ip,el),& + materialpoint_subdt(ip,el), & + crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & + ip, & + el) end select chosenHomogenization chosenThermal: select case (thermal_type(mesh_element(3,el))) @@ -793,12 +718,6 @@ end function updateState !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- subroutine averageStressAndItsTangent(ip,el) - use mesh, only: & - mesh_element - use crystallite, only: & - crystallite_P,crystallite_dPdF - use homogenization_mech_RGC, only: & - homogenization_RGC_averageStressAndItsTangent integer, intent(in) :: & ip, & !< integration point @@ -818,7 +737,7 @@ subroutine averageStressAndItsTangent(ip,el) homogenization_typeInstance(mesh_element(3,el))) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - call homogenization_RGC_averageStressAndItsTangent(& + call mech_RGC_averageStressAndItsTangent(& materialpoint_P(1:3,1:3,ip,el), & materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& crystallite_P(1:3,1:3,1:homogenization_Ngrains(mesh_element(3,el)),ip,el), & @@ -834,10 +753,6 @@ end subroutine averageStressAndItsTangent !> if homogenization_sizePostResults(i,e) > 0 !! !-------------------------------------------------------------------------------------------------- function postResults(ip,el) - use mesh, only: & - mesh_element - use homogenization_mech_RGC, only: & - homogenization_RGC_postResults use thermal_adiabatic, only: & thermal_adiabatic_postResults use thermal_conduction, only: & @@ -856,23 +771,12 @@ function postResults(ip,el) postResults integer :: & startPos, endPos ,& - of, instance, homog + homog postResults = 0.0_pReal startPos = 1 - endPos = homogState(material_homogenizationAt(el))%sizePostResults - chosenHomogenization: select case (homogenization_type(mesh_element(3,el))) - - case (HOMOGENIZATION_RGC_ID) chosenHomogenization - instance = homogenization_typeInstance(material_homogenizationAt(el)) - of = mappingHomogenization(1,ip,el) - postResults(startPos:endPos) = homogenization_RGC_postResults(instance,of) - - end select chosenHomogenization - - startPos = endPos + 1 - endPos = endPos + thermalState(material_homogenizationAt(el))%sizePostResults + endPos = thermalState(material_homogenizationAt(el))%sizePostResults chosenThermal: select case (thermal_type(mesh_element(3,el))) case (THERMAL_adiabatic_ID) chosenThermal @@ -905,14 +809,9 @@ end function postResults !-------------------------------------------------------------------------------------------------- subroutine homogenization_results #if defined(PETSc) || defined(DAMASK_HDF5) - use results - use homogenization_mech_RGC - use HDF5_utilities use config, only: & config_name_homogenization => homogenization_name ! anticipate logical name - use material, only: & - homogenization_typeInstance, & material_homogenization_type => homogenization_type integer :: p diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index d3feac1eb..1cbe837d2 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -6,16 +6,7 @@ !> @brief Relaxed grain cluster (RGC) homogenization scheme !> Nconstituents is defined as p x q x r (cluster) !-------------------------------------------------------------------------------------------------- -module homogenization_mech_RGC - use prec - use material - - implicit none - private - integer, dimension(:,:), allocatable,target, public :: & - homogenization_RGC_sizePostResult - character(len=64), dimension(:,:), allocatable,target, public :: & - homogenization_RGC_output ! name of each post result output +submodule(homogenization) homogenization_mech_RGC enum, bind(c) enumerator :: & @@ -28,7 +19,7 @@ module homogenization_mech_RGC magnitudemismatch_ID end enum - type, private :: tParameters + type :: tParameters integer, dimension(:), allocatable :: & Nconstituents real(pReal) :: & @@ -39,11 +30,11 @@ module homogenization_mech_RGC angles integer :: & of_debug = 0 - integer(kind(undefined_ID)), dimension(:), allocatable :: & + integer(kind(undefined_ID)), dimension(:), allocatable :: & outputID end type tParameters - type, private :: tRGCstate + type :: tRGCstate real(pReal), pointer, dimension(:) :: & work, & penaltyEnergy @@ -51,7 +42,7 @@ module homogenization_mech_RGC relaxationVector end type tRGCstate - type, private :: tRGCdependentState + type :: tRGCdependentState real(pReal), allocatable, dimension(:) :: & volumeDiscrepancy, & relaxationRate_avg, & @@ -62,56 +53,25 @@ module homogenization_mech_RGC orientation end type tRGCdependentState - type(tparameters), dimension(:), allocatable, private :: & + type(tparameters), dimension(:), allocatable :: & param - type(tRGCstate), dimension(:), allocatable, private :: & + type(tRGCstate), dimension(:), allocatable :: & state, & state0 - type(tRGCdependentState), dimension(:), allocatable, private :: & + type(tRGCdependentState), dimension(:), allocatable :: & dependentState - public :: & - homogenization_RGC_init, & - homogenization_RGC_partitionDeformation, & - homogenization_RGC_averageStressAndItsTangent, & - homogenization_RGC_updateState, & - homogenization_RGC_postResults, & - mech_RGC_results ! name suited for planned submodule situation - private :: & - relaxationVector, & - interfaceNormal, & - getInterface, & - grain1to3, & - grain3to1, & - interface4to1, & - interface1to4 - contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all necessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_init() - use debug, only: & -#ifdef DEBUG - debug_i, & - debug_e, & -#endif - debug_level, & - debug_homogenization, & - debug_levelBasic - use math, only: & - math_EulerToR, & - INRAD - use IO, only: & - IO_error - use config, only: & - config_homogenization +module subroutine mech_RGC_init integer :: & Ninstance, & h, i, & - NofMyHomog, outputSize, & + NofMyHomog, & sizeState, nIntFaceTot character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] @@ -139,9 +99,6 @@ subroutine homogenization_RGC_init() allocate(state0(Ninstance)) allocate(dependentState(Ninstance)) - allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput),Ninstance),source=0) - allocate(homogenization_RGC_output(maxval(homogenization_Noutput),Ninstance)) - homogenization_RGC_output='' do h = 1, size(homogenization_type) if (homogenization_type(h) /= HOMOGENIZATION_RGC_ID) cycle @@ -176,28 +133,20 @@ subroutine homogenization_RGC_init() case('constitutivework') outputID = constitutivework_ID - outputSize = 1 case('penaltyenergy') outputID = penaltyenergy_ID - outputSize = 1 case('volumediscrepancy') outputID = volumediscrepancy_ID - outputSize = 1 case('averagerelaxrate') outputID = averagerelaxrate_ID - outputSize = 1 case('maximumrelaxrate') outputID = maximumrelaxrate_ID - outputSize = 1 case('magnitudemismatch') outputID = magnitudemismatch_ID - outputSize = 3 end select if (outputID /= undefined_ID) then - homogenization_RGC_output(i,homogenization_typeInstance(h)) = outputs(i) - homogenization_RGC_sizePostResult(i,homogenization_typeInstance(h)) = outputSize prm%outputID = [prm%outputID , outputID] endif @@ -211,7 +160,7 @@ subroutine homogenization_RGC_init() + size(['avg constitutive work ','average penalty energy']) homogState(h)%sizeState = sizeState - homogState(h)%sizePostResults = sum(homogenization_RGC_sizePostResult(:,homogenization_typeInstance(h))) + homogState(h)%sizePostResults = 0 allocate(homogState(h)%state0 (sizeState,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%subState0(sizeState,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%state (sizeState,NofMyHomog), source=0.0_pReal) @@ -235,23 +184,17 @@ subroutine homogenization_RGC_init() enddo -end subroutine homogenization_RGC_init +end subroutine mech_RGC_init !-------------------------------------------------------------------------------------------------- !> @brief partitions the deformation gradient onto the constituents !-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_partitionDeformation(F,avgF,instance,of) -#ifdef DEBUG - use debug, only: & - debug_level, & - debug_homogenization, & - debug_levelExtensive -#endif +module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of) real(pReal), dimension (:,:,:), intent(out) :: F !< partioned F per grain - real(pReal), dimension (:,:), intent(in) :: avgF !< averaged F + real(pReal), dimension (3,3), intent(in) :: avgF !< averaged F integer, intent(in) :: & instance, & of @@ -291,45 +234,14 @@ subroutine homogenization_RGC_partitionDeformation(F,avgF,instance,of) end associate -end subroutine homogenization_RGC_partitionDeformation +end subroutine mech_RGC_partitionDeformation !-------------------------------------------------------------------------------------------------- !> @brief update the internal state of the homogenization scheme and tell whether "done" and ! "happy" with result !-------------------------------------------------------------------------------------------------- -function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) -#ifdef DEBUG - use debug, only: & - debug_level, & - debug_homogenization,& - debug_levelExtensive -#endif - use math, only: & - math_invert2 - use numerics, only: & - absTol_RGC, & - relTol_RGC, & - absMax_RGC, & - relMax_RGC, & - pPert_RGC, & - maxdRelax_RGC, & - viscPower_RGC, & - viscModus_RGC, & - refRelaxRate_RGC - - real(pReal), dimension(:,:,:), intent(in) :: & - P,& !< array of P - F,& !< array of F - F0 !< array of initial F - real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< array of current grain stiffness - real(pReal), dimension(3,3), intent(in) :: avgF !< average F - real(pReal), intent(in) :: dt !< time increment - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - - logical, dimension(2) :: homogenization_RGC_updateState +module procedure mech_RGC_updateState integer, dimension(4) :: intFaceN,intFaceP,faceID integer, dimension(3) :: nGDim,iGr3N,iGr3P @@ -347,7 +259,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) #endif zeroTimeStep: if(dEq0(dt)) then - homogenization_RGC_updateState = .true. ! pretend everything is fine and return + mech_RGC_updateState = .true. ! pretend everything is fine and return return endif zeroTimeStep @@ -466,12 +378,12 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) endif #endif - homogenization_RGC_updateState = .false. + mech_RGC_updateState = .false. !-------------------------------------------------------------------------------------------------- ! If convergence reached => done and happy if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then - homogenization_RGC_updateState = .true. + mech_RGC_updateState = .true. #ifdef DEBUG if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 & .and. prm%of_debug == of) write(6,'(1x,a55,/)')'... done and happy' @@ -513,7 +425,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !-------------------------------------------------------------------------------------------------- ! if residual blows-up => done but unhappy elseif (residMax > relMax_RGC*stresMax .or. residMax > absMax_RGC) then ! try to restart when residual blows up exceeding maximum bound - homogenization_RGC_updateState = [.true.,.false.] ! with direct cut-back + mech_RGC_updateState = [.true.,.false.] ! with direct cut-back #ifdef DEBUG if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0 & @@ -656,9 +568,10 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !-------------------------------------------------------------------------------------------------- ! ... of the numerical viscosity traction "rmatrix" allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal) - forall (i=1:3*nIntFaceTot) & + do i=1,3*nIntFaceTot rmatrix(i,i) = viscModus_RGC*viscPower_RGC/(refRelaxRate_RGC*dt)* & ! tangent due to numerical viscosity traction appears (abs(drelax(i))/(refRelaxRate_RGC*dt))**(viscPower_RGC - 1.0_pReal) ! only in the main diagonal term + enddo #ifdef DEBUG if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0) then @@ -710,7 +623,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) enddo; enddo stt%relaxationVector(:,of) = relax + drelax ! Updateing the state variable for the next iteration if (any(abs(drelax) > maxdRelax_RGC)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large - homogenization_RGC_updateState = [.true.,.false.] + mech_RGC_updateState = [.true.,.false.] !$OMP CRITICAL (write2out) write(6,'(1x,a,1x,i3,1x,a,1x,i3,1x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback' write(6,'(1x,a,1x,e15.8)')'due to large relaxation change =',maxval(abs(drelax)) @@ -736,10 +649,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !> @brief calculate stress-like penalty due to deformation mismatch !-------------------------------------------------------------------------------------------------- subroutine stressPenalty(rPen,nMis,avgF,fDef,ip,el,instance,of) - use math, only: & - math_civita - use numerics, only: & - xSmoo_RGC real(pReal), dimension (:,:,:), intent(out) :: rPen !< stress-like penalty real(pReal), dimension (:,:), intent(out) :: nMis !< total amount of mismatch @@ -852,13 +761,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !> @brief calculate stress-like penalty due to volume discrepancy !-------------------------------------------------------------------------------------------------- subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain,instance,of) - use math, only: & - math_det33, & - math_inv33 - use numerics, only: & - maxVolDiscr_RGC,& - volDiscrMod_RGC,& - volDiscrPow_RGC real(pReal), dimension (:,:,:), intent(out) :: vPen ! stress-like penalty due to volume real(pReal), intent(out) :: vDiscrep ! total volume discrepancy @@ -907,8 +809,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) ! deformation !-------------------------------------------------------------------------------------------------- function surfaceCorrection(avgF,instance,of) - use math, only: & - math_invert33 real(pReal), dimension(3) :: surfaceCorrection @@ -940,8 +840,6 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor !-------------------------------------------------------------------------------------------------- function equivalentModuli(grainID,ip,el) - use constitutive, only: & - constitutive_homogenizedC real(pReal), dimension(2) :: equivalentModuli @@ -1012,13 +910,13 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) end subroutine grainDeformation -end function homogenization_RGC_updateState +end procedure mech_RGC_updateState !-------------------------------------------------------------------------------------------------- !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- -subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) +module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point @@ -1030,68 +928,19 @@ subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF, avgP = sum(P,3) /real(product(param(instance)%Nconstituents),pReal) dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%Nconstituents),pReal) -end subroutine homogenization_RGC_averageStressAndItsTangent - - -!-------------------------------------------------------------------------------------------------- -!> @brief return array of homogenization results for post file inclusion -!-------------------------------------------------------------------------------------------------- -pure function homogenization_RGC_postResults(instance,of) result(postResults) - - integer, intent(in) :: & - instance, & - of - - integer :: & - o,c - real(pReal), dimension(sum(homogenization_RGC_sizePostResult(:,instance))) :: & - postResults - - associate(stt => state(instance), dst => dependentState(instance), prm => param(instance)) - - c = 0 - - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - - case (constitutivework_ID) - postResults(c+1) = stt%work(of) - c = c + 1 - case (magnitudemismatch_ID) - postResults(c+1:c+3) = dst%mismatch(1:3,of) - c = c + 3 - case (penaltyenergy_ID) - postResults(c+1) = stt%penaltyEnergy(of) - c = c + 1 - case (volumediscrepancy_ID) - postResults(c+1) = dst%volumeDiscrepancy(of) - c = c + 1 - case (averagerelaxrate_ID) - postResults(c+1) = dst%relaxationrate_avg(of) - c = c + 1 - case (maximumrelaxrate_ID) - postResults(c+1) = dst%relaxationrate_max(of) - c = c + 1 - end select - - enddo outputsLoop - - end associate - -end function homogenization_RGC_postResults +end subroutine mech_RGC_averageStressAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file ! ToDo: check wheter units are correct !-------------------------------------------------------------------------------------------------- -subroutine mech_RGC_results(instance,group) +module subroutine mech_RGC_results(instance,group) #if defined(PETSc) || defined(DAMASK_HDF5) - use results, only: & - results_writeDataset - integer, intent(in) :: instance - character(len=*) :: group + integer, intent(in) :: instance + character(len=*), intent(in) :: group + integer :: o associate(stt => state(instance), dst => dependentState(instance), prm => param(instance)) @@ -1122,8 +971,8 @@ subroutine mech_RGC_results(instance,group) end associate #else - integer, intent(in) :: instance - character(len=*) :: group + integer, intent(in) :: instance + character(len=*), intent(in) :: group #endif end subroutine mech_RGC_results @@ -1262,7 +1111,7 @@ integer pure function interface4to1(iFace4D, nGDim) else interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1) & + nGDim(3)*nGDim(1)*(iFace4D(3)-1) & - + (nGDim(1)-1)*nGDim(2)*nGDim(3) ! total number of interfaces normal //e1 + + (nGDim(1)-1)*nGDim(2)*nGDim(3) ! total # of interfaces normal || e1 endif case(3) @@ -1271,8 +1120,8 @@ integer pure function interface4to1(iFace4D, nGDim) else interface4to1 = iFace4D(2) + nGDim(1)*(iFace4D(3)-1) & + nGDim(1)*nGDim(2)*(iFace4D(4)-1) & - + (nGDim(1)-1)*nGDim(2)*nGDim(3) & ! total number of interfaces normal //e1 - + nGDim(1)*(nGDim(2)-1)*nGDim(3) ! total number of interfaces normal //e2 + + (nGDim(1)-1)*nGDim(2)*nGDim(3) & ! total # of interfaces normal || e1 + + nGDim(1)*(nGDim(2)-1)*nGDim(3) ! total # of interfaces normal || e2 endif case default @@ -1296,23 +1145,23 @@ pure function interface1to4(iFace1D, nGDim) !-------------------------------------------------------------------------------------------------- ! compute the total number of interfaces, which ... - nIntFace = [(nGDim(1)-1)*nGDim(2)*nGDim(3), & ! ... normal //e1 - nGDim(1)*(nGDim(2)-1)*nGDim(3), & ! ... normal //e2 - nGDim(1)*nGDim(2)*(nGDim(3)-1)] ! ... normal //e3 + nIntFace = [(nGDim(1)-1)*nGDim(2)*nGDim(3), & ! ... normal || e1 + nGDim(1)*(nGDim(2)-1)*nGDim(3), & ! ... normal || e2 + nGDim(1)*nGDim(2)*(nGDim(3)-1)] ! ... normal || e3 !-------------------------------------------------------------------------------------------------- ! get the corresponding interface ID in 4D (normal and local position) - if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then ! interface with normal //e1 + if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then ! interface with normal || e1 interface1to4(1) = 1 interface1to4(3) = mod((iFace1D-1),nGDim(2))+1 interface1to4(4) = mod(int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)),nGDim(3))+1 interface1to4(2) = int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)/real(nGDim(3),pReal))+1 - elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then ! interface with normal //e2 + elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then ! interface with normal || e2 interface1to4(1) = 2 interface1to4(4) = mod((iFace1D-nIntFace(1)-1),nGDim(3))+1 interface1to4(2) = mod(int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)),nGDim(1))+1 interface1to4(3) = int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)/real(nGDim(1),pReal))+1 - elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then ! interface with normal //e3 + elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then ! interface with normal || e3 interface1to4(1) = 3 interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1),nGDim(1))+1 interface1to4(3) = mod(int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)),nGDim(2))+1 @@ -1322,4 +1171,4 @@ pure function interface1to4(iFace1D, nGDim) end function interface1to4 -end module homogenization_mech_RGC +end submodule homogenization_mech_RGC diff --git a/src/homogenization_mech_isostrain.f90 b/src/homogenization_mech_isostrain.f90 index 7dd7bad7d..c1eba821c 100644 --- a/src/homogenization_mech_isostrain.f90 +++ b/src/homogenization_mech_isostrain.f90 @@ -6,8 +6,6 @@ !-------------------------------------------------------------------------------------------------- submodule(homogenization) homogenization_mech_isostrain - implicit none - enum, bind(c) enumerator :: & parallel_ID, & @@ -30,14 +28,6 @@ contains !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- module subroutine mech_isostrain_init - use debug, only: & - debug_HOMOGENIZATION, & - debug_level, & - debug_levelBasic - use IO, only: & - IO_error - use config, only: & - config_homogenization integer :: & Ninstance, & diff --git a/src/homogenization_mech_none.f90 b/src/homogenization_mech_none.f90 index e7a5a12e6..d5b24e5d1 100644 --- a/src/homogenization_mech_none.f90 +++ b/src/homogenization_mech_none.f90 @@ -6,20 +6,12 @@ !-------------------------------------------------------------------------------------------------- submodule(homogenization) homogenization_mech_none - implicit none - contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- module subroutine mech_none_init - use debug, only: & - debug_HOMOGENIZATION, & - debug_level, & - debug_levelBasic - use config, only: & - config_homogenization integer :: & Ninstance, & @@ -27,7 +19,7 @@ module subroutine mech_none_init NofMyHomog write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_NONE_label//' init -+>>>' - + Ninstance = count(homogenization_type == HOMOGENIZATION_NONE_ID) if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90 index 7db1f5f7f..cb4b3cbae 100644 --- a/src/plastic_dislotwin.f90 +++ b/src/plastic_dislotwin.f90 @@ -1067,8 +1067,8 @@ subroutine plastic_dislotwin_results(instance,group) case (rho_dip_ID) call results_writeDataset(group,stt%rho_dip,'rho_dip',& 'dislocation dipole density''1/m²') - case (dot_gamma_sl_ID) - call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',& + case (gamma_sl_ID) + call results_writeDataset(group,stt%gamma_sl,'gamma_sl',& 'plastic shear','1') case (Lambda_sl_ID) call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',& diff --git a/src/results.f90 b/src/results.f90 index c5582b927..05db831f7 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -299,18 +299,26 @@ end subroutine results_writeVectorDataset_real !-------------------------------------------------------------------------------------------------- subroutine results_writeTensorDataset_real(group,dataset,label,description,SIunit) - character(len=*), intent(in) :: label,group,description - character(len=*), intent(in), optional :: SIunit - real(pReal), intent(inout), dimension(:,:,:) :: dataset + character(len=*), intent(in) :: label,group,description + character(len=*), intent(in), optional :: SIunit + real(pReal), intent(in), dimension(:,:,:) :: dataset + integer :: i integer(HID_T) :: groupHandle - + real(pReal), dimension(:,:,:), allocatable :: dataset_transposed + + + allocate(dataset_transposed,mold=dataset) + do i=1,size(dataset,3) + dataset_transposed(1:3,1:3,i) = transpose(dataset(1:3,1:3,i)) + enddo + groupHandle = results_openGroup(group) #ifdef PETSc - call HDF5_write(groupHandle,dataset,label,.true.) + call HDF5_write(groupHandle,dataset_transposed,label,.true.) #else - call HDF5_write(groupHandle,dataset,label,.false.) + call HDF5_write(groupHandle,dataset_transposed,label,.false.) #endif if (HDF5_objectExists(groupHandle,label)) &