From 9b4252da68628d6f3bf5644882cdee51bf35087c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 29 Nov 2019 16:53:40 +0100 Subject: [PATCH 01/19] smart handling of leading zeros --- processing/post/DADF5_postResults.py | 5 ++++- processing/post/DADF5_vtk_cells.py | 6 ++++-- processing/post/DADF5_vtk_points.py | 5 ++++- 3 files changed, 12 insertions(+), 4 deletions(-) diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py index a16ef147c..ab68fc76b 100755 --- a/processing/post/DADF5_postResults.py +++ b/processing/post/DADF5_postResults.py @@ -47,6 +47,8 @@ for filename in options.filenames: coords = np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) + N_digits = int(np.floor(np.log10(int(results.increments[-1][3:]))))+1 + N_digits = 5 # hack to keep test intact for i,inc in enumerate(results.iter_visible('increments')): print('Output step {}/{}'.format(i+1,len(results.increments))) @@ -92,5 +94,6 @@ for filename in options.filenames: dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) if not os.path.isdir(dirname): os.mkdir(dirname,0o755) - file_out = '{}_{}.txt'.format(os.path.splitext(os.path.split(filename)[-1])[0],inc) + file_out = '{}_{}.txt'.format(os.path.splitext(os.path.split(filename)[-1])[0], + inc[3:].zfill(N_digits)) np.savetxt(os.path.join(dirname,file_out),data,header=header,comments='') diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py index 1f5cc6686..58cb9771f 100755 --- a/processing/post/DADF5_vtk_cells.py +++ b/processing/post/DADF5_vtk_cells.py @@ -66,7 +66,7 @@ for filename in options.filenames: for i in f['/geometry/T_c']: grid.InsertNextCell(vtk.VTK_HEXAHEDRON,8,i-1) - + N_digits = int(np.floor(np.log10(int(results.increments[-1][3:]))))+1 for i,inc in enumerate(results.iter_visible('increments')): print('Output step {}/{}'.format(i+1,len(results.increments))) vtk_data = [] @@ -133,7 +133,9 @@ for filename in options.filenames: dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) if not os.path.isdir(dirname): os.mkdir(dirname,0o755) - file_out = '{}_{}.{}'.format(os.path.splitext(os.path.split(filename)[-1])[0],inc,writer.GetDefaultFileExtension()) + file_out = '{}_inc{}.{}'.format(os.path.splitext(os.path.split(filename)[-1])[0], + inc[3:].zfill(N_digits), + writer.GetDefaultFileExtension()) writer.SetCompressorTypeToZLib() writer.SetDataModeToBinary() diff --git a/processing/post/DADF5_vtk_points.py b/processing/post/DADF5_vtk_points.py index 87c1ad93e..9265cc3a0 100755 --- a/processing/post/DADF5_vtk_points.py +++ b/processing/post/DADF5_vtk_points.py @@ -52,6 +52,7 @@ for filename in options.filenames: Polydata.SetVerts(Vertices) Polydata.Modified() + N_digits = int(np.floor(np.log10(int(results.increments[-1][3:]))))+1 for i,inc in enumerate(results.iter_visible('increments')): print('Output step {}/{}'.format(i+1,len(results.increments))) vtk_data = [] @@ -111,7 +112,9 @@ for filename in options.filenames: dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) if not os.path.isdir(dirname): os.mkdir(dirname,0o755) - file_out = '{}_{}.{}'.format(os.path.splitext(os.path.split(filename)[-1])[0],inc,writer.GetDefaultFileExtension()) + file_out = '{}_inc{}.{}'.format(os.path.splitext(os.path.split(filename)[-1])[0], + inc[3:].zfill(N_digits), + writer.GetDefaultFileExtension()) writer.SetCompressorTypeToZLib() writer.SetDataModeToBinary() From dc3fc8f70cbd0fdb24af2d57d9bae294edba7366 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 29 Nov 2019 17:00:48 +0100 Subject: [PATCH 02/19] do not clutter with leading zeros --- python/damask/dadf5.py | 7 ++++--- src/results.f90 | 6 +++--- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index d879946eb..ee497a9ae 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -30,7 +30,7 @@ class DADF5(): """ with h5py.File(filename,'r') as f: - if f.attrs['DADF5-major'] != 0 or not 2 <= f.attrs['DADF5-minor'] <= 3: + if f.attrs['DADF5-major'] != 0 or not 2 <= f.attrs['DADF5-minor'] <= 4: raise TypeError('Unsupported DADF5 version {} '.format(f.attrs['DADF5-version'])) self.structured = 'grid' in f['geometry'].attrs.keys() @@ -40,8 +40,9 @@ class DADF5(): self.size = f['geometry'].attrs['size'] r=re.compile('inc[0-9]+') - self.increments = [i for i in f.keys() if r.match(i)] - self.times = [round(f[i].attrs['time/s'],12) for i in self.increments] + increments_unsorted = {int(i[3:]):i for i in f.keys() if r.match(i)} + self.increments = [increments_unsorted[i] for i in sorted(increments_unsorted)] + self.times = [round(f[i].attrs['time/s'],12) for i in self.increments] self.Nmaterialpoints, self.Nconstituents = np.shape(f['mapping/cellResults/constituent']) self.materialpoints = [m.decode() for m in np.unique(f['mapping/cellResults/materialpoint']['Name'])] diff --git a/src/results.f90 b/src/results.f90 index 471f994d6..e77600d82 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -68,9 +68,9 @@ subroutine results_init write(6,'(a)') ' https://doi.org/10.1007/s40192-017-0084-5' resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','w',.true.) - call HDF5_addAttribute(resultsFile,'DADF5-version',0.3_pReal) + call HDF5_addAttribute(resultsFile,'DADF5-version',0.4_pReal) call HDF5_addAttribute(resultsFile,'DADF5-major',0) - call HDF5_addAttribute(resultsFile,'DADF5-minor',3) + call HDF5_addAttribute(resultsFile,'DADF5-minor',4) call HDF5_addAttribute(resultsFile,'DAMASK',DAMASKVERSION) call get_command(commandLine) call HDF5_addAttribute(resultsFile,'call',trim(commandLine)) @@ -110,7 +110,7 @@ subroutine results_addIncrement(inc,time) real(pReal), intent(in) :: time character(len=pStringLen) :: incChar - write(incChar,'(i5.5)') inc ! allow up to 99999 increments + write(incChar,'(i10)') inc call HDF5_closeGroup(results_addGroup(trim('inc'//trim(adjustl(incChar))))) call results_setLink(trim('inc'//trim(adjustl(incChar))),'current') call HDF5_addAttribute(resultsFile,'time/s',time,trim('inc'//trim(adjustl(incChar)))) From e2b13a5ca3bb10042edc591edf994ecd3bae80d8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 29 Nov 2019 19:32:18 +0100 Subject: [PATCH 03/19] using de-facto standard naming --- processing/post/DADF5_postResults.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py index ab68fc76b..1d05289cb 100755 --- a/processing/post/DADF5_postResults.py +++ b/processing/post/DADF5_postResults.py @@ -94,6 +94,6 @@ for filename in options.filenames: dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) if not os.path.isdir(dirname): os.mkdir(dirname,0o755) - file_out = '{}_{}.txt'.format(os.path.splitext(os.path.split(filename)[-1])[0], - inc[3:].zfill(N_digits)) + file_out = '{}_inc{}.txt'.format(os.path.splitext(os.path.split(filename)[-1])[0], + inc[3:].zfill(N_digits)) np.savetxt(os.path.join(dirname,file_out),data,header=header,comments='') From 4185fcb4c3f3c8505dfae9ca2d0386668337a195 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 29 Nov 2019 21:40:54 +0100 Subject: [PATCH 04/19] adjust to new naming --- python/damask/dadf5.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index ee497a9ae..77ac5458b 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -166,7 +166,7 @@ class DADF5(): end increment (included) """ - self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','set') + self.__manage_visible(['inc{}'.format(i) for i in range(start,end+1)],'increments','set') def add_by_increment(self,start,end): @@ -181,7 +181,7 @@ class DADF5(): end increment (included) """ - self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','add') + self.__manage_visible(['inc{}'.format(i) for i in range(start,end+1)],'increments','add') def del_by_increment(self,start,end): @@ -196,7 +196,7 @@ class DADF5(): end increment (included) """ - self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','del') + self.__manage_visible(['inc{}'.format(i) for i in range(start,end+1)],'increments','del') def iter_visible(self,what): From 354c0123a1bb84f40443b03992c569e71201ead1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 30 Nov 2019 08:40:59 +0100 Subject: [PATCH 05/19] backward compatibility --- python/damask/dadf5.py | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 77ac5458b..7dddabfd4 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -30,7 +30,9 @@ class DADF5(): """ with h5py.File(filename,'r') as f: - if f.attrs['DADF5-major'] != 0 or not 2 <= f.attrs['DADF5-minor'] <= 4: + self.DADF5-major = f.attrs['DADF5-major'] + self.DADF5-minor = f.attrs['DADF5-minor'] + if self.DADF5-major != 0 or not 2 <= self.DADF5-minor <= 4: raise TypeError('Unsupported DADF5 version {} '.format(f.attrs['DADF5-version'])) self.structured = 'grid' in f['geometry'].attrs.keys() @@ -166,7 +168,10 @@ class DADF5(): end increment (included) """ - self.__manage_visible(['inc{}'.format(i) for i in range(start,end+1)],'increments','set') + if self.DADF5-minor >= 4: + self.__manage_visible([ 'inc{}'.format(i) for i in range(start,end+1)],'increments','set') + else: + self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','set') def add_by_increment(self,start,end): @@ -181,7 +186,10 @@ class DADF5(): end increment (included) """ - self.__manage_visible(['inc{}'.format(i) for i in range(start,end+1)],'increments','add') + if self.DADF5-minor >= 4: + self.__manage_visible([ 'inc{}'.format(i) for i in range(start,end+1)],'increments','add') + else: + self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','add') def del_by_increment(self,start,end): @@ -196,7 +204,10 @@ class DADF5(): end increment (included) """ - self.__manage_visible(['inc{}'.format(i) for i in range(start,end+1)],'increments','del') + if self.DADF5-minor >= 4: + self.__manage_visible([ 'inc{}'.format(i) for i in range(start,end+1)],'increments','del') + else: + self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','del') def iter_visible(self,what): From 8f77843b0f207ecf731d0f79e730cb4573c4ec99 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 30 Nov 2019 09:01:37 +0100 Subject: [PATCH 06/19] improved naming - clearer variable names - no '-' in attributes - single source of truth --- python/damask/dadf5.py | 17 +++++++++++------ src/results.f90 | 7 +++---- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 7dddabfd4..d03617752 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -30,9 +30,14 @@ class DADF5(): """ with h5py.File(filename,'r') as f: - self.DADF5-major = f.attrs['DADF5-major'] - self.DADF5-minor = f.attrs['DADF5-minor'] - if self.DADF5-major != 0 or not 2 <= self.DADF5-minor <= 4: + try: + self.version_major = f.attrs['DADF5_version_major'] + self.version_minor = f.attrs['DADF5_version_minor'] + except KeyError: + self.version_major = f.attrs['DADF5-major'] + self.version_minor = f.attrs['DADF5-minor'] + + if self.version_major != 0 or not 2 <= self.version_minor <= 4: raise TypeError('Unsupported DADF5 version {} '.format(f.attrs['DADF5-version'])) self.structured = 'grid' in f['geometry'].attrs.keys() @@ -168,7 +173,7 @@ class DADF5(): end increment (included) """ - if self.DADF5-minor >= 4: + if self.version_minor >= 4: self.__manage_visible([ 'inc{}'.format(i) for i in range(start,end+1)],'increments','set') else: self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','set') @@ -186,7 +191,7 @@ class DADF5(): end increment (included) """ - if self.DADF5-minor >= 4: + if self.version_minor >= 4: self.__manage_visible([ 'inc{}'.format(i) for i in range(start,end+1)],'increments','add') else: self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','add') @@ -204,7 +209,7 @@ class DADF5(): end increment (included) """ - if self.DADF5-minor >= 4: + if self.version_minor >= 4: self.__manage_visible([ 'inc{}'.format(i) for i in range(start,end+1)],'increments','del') else: self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','del') diff --git a/src/results.f90 b/src/results.f90 index e77600d82..a6b01790d 100644 --- a/src/results.f90 +++ b/src/results.f90 @@ -68,10 +68,9 @@ subroutine results_init write(6,'(a)') ' https://doi.org/10.1007/s40192-017-0084-5' resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','w',.true.) - call HDF5_addAttribute(resultsFile,'DADF5-version',0.4_pReal) - call HDF5_addAttribute(resultsFile,'DADF5-major',0) - call HDF5_addAttribute(resultsFile,'DADF5-minor',4) - call HDF5_addAttribute(resultsFile,'DAMASK',DAMASKVERSION) + call HDF5_addAttribute(resultsFile,'DADF5_version_major',0) + call HDF5_addAttribute(resultsFile,'DADF5_version_minor',4) + call HDF5_addAttribute(resultsFile,'DAMASK_version',DAMASKVERSION) call get_command(commandLine) call HDF5_addAttribute(resultsFile,'call',trim(commandLine)) call HDF5_closeGroup(results_addGroup('mapping')) From 8189b50509962ce4a68bd1fbc85054a121e06c76 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 2 Dec 2019 12:58:23 +0100 Subject: [PATCH 07/19] same functionality but tested --- src/math.f90 | 27 --------------------------- src/plastic_nonlocal.f90 | 15 ++++++--------- 2 files changed, 6 insertions(+), 36 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 7eba379d2..6b83e17a3 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -838,33 +838,6 @@ pure function math_Voigt66to3333(m66) end function math_Voigt66to3333 -!-------------------------------------------------------------------------------------------------- -!> @brief action of a quaternion on a vector (rotate vector v with Q) -!> @details deprecated -!-------------------------------------------------------------------------------------------------- -pure function math_qRot(Q,v) - - real(pReal), dimension(4), intent(in) :: Q - real(pReal), dimension(3), intent(in) :: v - real(pReal), dimension(3) :: math_qRot - real(pReal), dimension(4,4) :: T - integer :: i, j - - do i = 1,4 - do j = 1,i - T(i,j) = Q(i) * Q(j) - enddo - enddo - - math_qRot = [-v(1)*(T(3,3)+T(4,4)) + v(2)*(T(3,2)-T(4,1)) + v(3)*(T(4,2)+T(3,1)), & - v(1)*(T(3,2)+T(4,1)) - v(2)*(T(2,2)+T(4,4)) + v(3)*(T(4,3)-T(2,1)), & - v(1)*(T(4,2)-T(3,1)) + v(2)*(T(4,3)+T(2,1)) - v(3)*(T(2,2)+T(3,3))] - - math_qRot = 2.0_pReal * math_qRot + v - -end function math_qRot - - !-------------------------------------------------------------------------------------------------- !> @brief rotation matrix from Bunge-Euler (3-1-3) angles (in radians) !> @details deprecated diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90 index 23bfb50aa..5469515dc 100644 --- a/src/plastic_nonlocal.f90 +++ b/src/plastic_nonlocal.f90 @@ -1836,8 +1836,6 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) ns, & ! number of active slip systems s1, & ! slip system index (me) s2 ! slip system index (my neighbor) - real(pReal), dimension(4) :: & - absoluteMisorientation ! absolute misorientation (without symmetry) between me and my neighbor real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phaseAt(1,e))),& totalNslip(phase_plasticityInstance(material_phaseAt(1,e))),& nIPneighbors) :: & @@ -1848,7 +1846,7 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) nThresholdValues logical, dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,e)))) :: & belowThreshold - type(rotation) :: rot + type(rotation) :: mis Nneighbors = nIPneighbors ph = material_phaseAt(1,e) @@ -1914,18 +1912,17 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. !* All values below the threshold are set to zero. else - rot = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e)) - absoluteMisorientation = rot%asQuaternion() + mis = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e)) mySlipSystems: do s1 = 1,ns neighborSlipSystems: do s2 = 1,ns my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), & - math_qRot(absoluteMisorientation, prm%slip_normal(1:3,s2))) & + mis%rotate(prm%slip_normal(1:3,s2))) & * abs(math_inner(prm%slip_direction(1:3,s1), & - math_qRot(absoluteMisorientation, prm%slip_direction(1:3,s2)))) + mis%rotate(prm%slip_direction(1:3,s2)))) my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), & - math_qRot(absoluteMisorientation, prm%slip_normal(1:3,s2)))) & + mis%rotate(prm%slip_normal(1:3,s2)))) & * abs(math_inner(prm%slip_direction(1:3,s1), & - math_qRot(absoluteMisorientation, prm%slip_direction(1:3,s2)))) + mis%rotate(prm%slip_direction(1:3,s2)))) enddo neighborSlipSystems my_compatibilitySum = 0.0_pReal From be099e38c2fcc451b6774ebe3da1b053d899e0f6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 2 Dec 2019 16:22:27 +0100 Subject: [PATCH 08/19] might be of use --- src/rotations.f90 | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/src/rotations.f90 b/src/rotations.f90 index e042cff21..cd4a64547 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -64,6 +64,7 @@ module rotations procedure, public :: asRodrigues procedure, public :: asMatrix !------------------------------------------ + procedure, public :: fromQuaternion procedure, public :: fromEulers procedure, public :: fromAxisAngle procedure, public :: fromMatrix @@ -157,6 +158,18 @@ end function asHomochoric !--------------------------------------------------------------------------------------------------- ! Initialize rotation from different representations !--------------------------------------------------------------------------------------------------- +subroutine fromQuaternion(self,qu) + + class(rotation), intent(out) :: self + real(pReal), dimension(4), intent(in) :: qu + + if (dNeq(norm2(qu),1.0)) & + call IO_error(402,ext_msg='fromQuaternion') + + self%q = qu + +end subroutine fromQuaternion +!--------------------------------------------------------------------------------------------------- subroutine fromEulers(self,eu,degrees) class(rotation), intent(out) :: self From 83453d10ef25e079c562faed6082e3e13e7fe9bf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 2 Dec 2019 16:37:22 +0100 Subject: [PATCH 09/19] use rotation class for consistent handling of rotations --- src/grid/DAMASK_grid.f90 | 24 ++++++++++-------------- src/grid/spectral_utilities.f90 | 3 ++- 2 files changed, 12 insertions(+), 15 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index e83cf3283..d71dcac18 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -28,7 +28,6 @@ program DAMASK_spectral use grid_damage_spectral use grid_thermal_spectral use results - use rotations implicit none @@ -78,7 +77,6 @@ program DAMASK_spectral character(len=6) :: loadcase_string character(len=1024) :: & incInfo - type(rotation) :: R type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tLoadCase) :: newLoadCase type(tSolutionState), allocatable, dimension(:) :: solres @@ -189,6 +187,7 @@ program DAMASK_spectral newLoadCase%ID(field) = FIELD_DAMAGE_ID endif damageActive + call newLoadCase%rot%fromEulers(real([0.0,0.0,0.0],pReal)) readIn: do i = 1, chunkPos(1) select case (IO_lc(IO_stringValue(line,chunkPos,i))) case('fdot','dotf','l','f') ! assign values for the deformation BC matrix @@ -244,14 +243,13 @@ program DAMASK_spectral do j = 1, 3 temp_valueVector(j) = IO_floatValue(line,chunkPos,i+k+j) enddo - call R%fromEulers(temp_valueVector(1:3),degrees=(l==1)) - newLoadCase%rotation = R%asMatrix() + call newLoadCase%rot%fromEulers(temp_valueVector(1:3),degrees=(l==1)) case('rotation','rot') ! assign values for the rotation matrix temp_valueVector = 0.0_pReal do j = 1, 9 temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) enddo - newLoadCase%rotation = math_9to33(temp_valueVector) + call newLoadCase%rot%fromMatrix(math_9to33(temp_valueVector)) end select enddo readIn @@ -295,14 +293,12 @@ program DAMASK_spectral endif enddo; write(6,'(/)',advance='no') enddo - if (any(abs(matmul(newLoadCase%rotation, & - transpose(newLoadCase%rotation))-math_I3) > & - reshape(spread(tol_math_check,1,9),[ 3,3]))& - .or. abs(math_det33(newLoadCase%rotation)) > & - 1.0_pReal + tol_math_check) errorID = 846 ! given rotation matrix contains strain - if (any(dNeq(newLoadCase%rotation, math_I3))) & + if (any(abs(matmul(newLoadCase%rot%asMatrix(), & + transpose(newLoadCase%rot%asMatrix()))-math_I3) > & + reshape(spread(tol_math_check,1,9),[ 3,3]))) errorID = 846 ! given rotation matrix contains strain + if (any(dNeq(newLoadCase%rot%asMatrix(), math_I3))) & write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',& - transpose(newLoadCase%rotation) + transpose(newLoadCase%rot%asMatrix()) if (newLoadCase%time < 0.0_pReal) errorID = 834 ! negative time increment write(6,'(2x,a,f12.6)') 'time: ', newLoadCase%time if (newLoadCase%incs < 1) errorID = 835 ! non-positive incs count @@ -469,7 +465,7 @@ program DAMASK_spectral cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, & deformation_BC = loadCases(currentLoadCase)%deformation, & stress_BC = loadCases(currentLoadCase)%stress, & - rotation_BC = loadCases(currentLoadCase)%rotation) + rotation_BC = loadCases(currentLoadCase)%rot%asMatrix()) case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack) case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack) @@ -488,7 +484,7 @@ program DAMASK_spectral solres(field) = mech_solution (& incInfo,timeinc,timeIncOld, & stress_BC = loadCases(currentLoadCase)%stress, & - rotation_BC = loadCases(currentLoadCase)%rotation) + rotation_BC = loadCases(currentLoadCase)%rot%asMatrix()) case(FIELD_THERMAL_ID) solres(field) = grid_thermal_spectral_solution(timeinc,timeIncOld) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 8b0e430f3..5107bd900 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -10,6 +10,7 @@ module spectral_utilities use prec use math + use rotations use IO use mesh_grid use numerics @@ -90,7 +91,7 @@ module spectral_utilities end type tBoundaryCondition type, public :: tLoadCase - real(pReal), dimension (3,3) :: rotation = math_I3 !< rotation of BC + type(rotation) :: rot !< rotation of BC type(tBoundaryCondition) :: stress, & !< stress BC deformation !< deformation BC (Fdot or L) real(pReal) :: time = 0.0_pReal !< length of increment From 8a9d3f8d6d48d579714b39d5901a44e70c732797 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 2 Dec 2019 20:06:58 +0100 Subject: [PATCH 10/19] avoid code duplication --- src/grid/DAMASK_grid.f90 | 4 +-- src/grid/grid_mech_FEM.f90 | 13 ++++----- src/grid/grid_mech_spectral_basic.f90 | 21 ++++++++------- src/grid/grid_mech_spectral_polarisation.f90 | 28 +++++++++++--------- src/grid/spectral_utilities.f90 | 3 ++- src/math.f90 | 14 ---------- src/rotations.f90 | 2 +- 7 files changed, 38 insertions(+), 47 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index d71dcac18..51c97456b 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -465,7 +465,7 @@ program DAMASK_spectral cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, & deformation_BC = loadCases(currentLoadCase)%deformation, & stress_BC = loadCases(currentLoadCase)%stress, & - rotation_BC = loadCases(currentLoadCase)%rot%asMatrix()) + rotation_BC = loadCases(currentLoadCase)%rot) case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack) case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack) @@ -484,7 +484,7 @@ program DAMASK_spectral solres(field) = mech_solution (& incInfo,timeinc,timeIncOld, & stress_BC = loadCases(currentLoadCase)%stress, & - rotation_BC = loadCases(currentLoadCase)%rot%asMatrix()) + rotation_BC = loadCases(currentLoadCase)%rot) case(FIELD_THERMAL_ID) solres(field) = grid_thermal_spectral_solution(timeinc,timeIncOld) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 0c3844fcf..be79af835 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -242,7 +242,8 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation timeinc_old !< time increment of last successful increment type(tBoundaryCondition), intent(in) :: & stress_BC - real(pReal), dimension(3,3), intent(in) :: rotation_BC + type(rotation), intent(in) :: & + rotation_BC type(tSolutionState) :: & solution !-------------------------------------------------------------------------------------------------- @@ -254,7 +255,7 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + S = utilities_maskedCompliance(rotation_BC%asMatrix(),stress_BC%maskLogical,C_volAvg) !-------------------------------------------------------------------------------------------------- ! set module wide available data params%stress_mask = stress_BC%maskFloat @@ -297,7 +298,7 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - real(pReal), dimension(3,3), intent(in) :: & + type(rotation), intent(in) :: & rotation_BC PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: & @@ -482,7 +483,7 @@ subroutine formResidual(da_local,x_local, & trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) + ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' deformation gradient aim =', transpose(F_aim) flush(6) @@ -498,7 +499,7 @@ subroutine formResidual(da_local,x_local, & x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) enddo; enddo; enddo ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 - F(1:3,1:3,ii,jj,kk) = math_rotate_backward33(F_aim,params%rotation_BC) + transpose(matmul(BMat,x_elem)) + F(1:3,1:3,ii,jj,kk) = params%rotation_BC%rotTensor2(F_aim,active=.true.) + transpose(matmul(BMat,x_elem)) enddo; enddo; enddo call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr) @@ -506,7 +507,7 @@ subroutine formResidual(da_local,x_local, & ! evaluate constitutive response call Utilities_constitutiveResponse(P_current,& P_av,C_volAvg,devNull, & - F,params%timeinc,params%rotation_BC) + F,params%timeinc,params%rotation_BC%asMatrix()) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 3dc978dc6..74b56db61 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -212,7 +212,8 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ timeinc_old !< time increment of last successful increment type(tBoundaryCondition), intent(in) :: & stress_BC - real(pReal), dimension(3,3), intent(in) :: rotation_BC + type(rotation), intent(in) :: & + rotation_BC type(tSolutionState) :: & solution !-------------------------------------------------------------------------------------------------- @@ -224,7 +225,7 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + S = utilities_maskedCompliance(rotation_BC%asMatrix(),stress_BC%maskLogical,C_volAvg) if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg) !-------------------------------------------------------------------------------------------------- @@ -269,7 +270,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - real(pReal), dimension(3,3), intent(in) :: & + type(rotation), intent(in) :: & rotation_BC PetscErrorCode :: ierr PetscScalar, dimension(:,:,:,:), pointer :: F @@ -299,9 +300,9 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime endif - Fdot = utilities_calculateRate(guess, & - F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, & - math_rotate_backward33(F_aimDot,rotation_BC)) + Fdot = utilities_calculateRate(guess, & + F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, & + rotation_BC%rotTensor2(F_aimDot,active=.true.)) F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3]) materialpoint_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3]) @@ -311,7 +312,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average - math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3]) + rotation_BC%rotTensor2(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) end subroutine grid_mech_spectral_basic_forward @@ -446,7 +447,7 @@ subroutine formResidual(in, F, & trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) + ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' deformation gradient aim =', transpose(F_aim) flush(6) @@ -456,7 +457,7 @@ subroutine formResidual(in, F, & ! evaluate constitutive response call utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory) P_av,C_volAvg,C_minMaxAvg, & - F,params%timeinc,params%rotation_BC) + F,params%timeinc,params%rotation_BC%asMatrix()) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- @@ -471,7 +472,7 @@ subroutine formResidual(in, F, & tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform call utilities_FFTtensorForward ! FFT forward of global "tensorField_real" err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use - call utilities_fourierGammaConvolution(math_rotate_backward33(deltaF_aim,params%rotation_BC)) ! convolution of Gamma and tensorField_fourier, with arg + call utilities_fourierGammaConvolution(params%rotation_BC%rotTensor2(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier call utilities_FFTtensorBackward ! FFT backward of global tensorField_fourier !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index c1b5d79c9..5c5a9f10c 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -14,6 +14,7 @@ module grid_mech_spectral_polarisation use DAMASK_interface use HDF5_utilities use math + use rotations use spectral_utilities use IO use FEsolving @@ -222,12 +223,13 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old, ! input data for solution character(len=*), intent(in) :: & incInfoIn - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & timeinc, & !< time increment of current solution timeinc_old !< time increment of last successful increment type(tBoundaryCondition), intent(in) :: & stress_BC - real(pReal), dimension(3,3), intent(in) :: rotation_BC + type(rotation), intent(in) :: & + rotation_BC type(tSolutionState) :: & solution !-------------------------------------------------------------------------------------------------- @@ -239,7 +241,7 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old, !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + S = utilities_maskedCompliance(rotation_BC%asMatrix(),stress_BC%maskLogical,C_volAvg) if (num%update_gamma) then call utilities_updateGamma(C_minMaxAvg) C_scale = C_minMaxAvg @@ -288,7 +290,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - real(pReal), dimension(3,3), intent(in) ::& + type(rotation), intent(in) :: & rotation_BC PetscErrorCode :: ierr PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau @@ -324,10 +326,10 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc Fdot = utilities_calculateRate(guess, & F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, & - math_rotate_backward33(F_aimDot,rotation_BC)) + rotation_BC%rotTensor2(F_aimDot,active=.true.)) F_tauDot = utilities_calculateRate(guess, & F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, & - math_rotate_backward33(F_aimDot,rotation_BC)) + rotation_BC%rotTensor2(F_aimDot,active=.true.)) F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) @@ -338,7 +340,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average - math_rotate_backward33(F_aim,rotation_BC)),& + rotation_BC%rotTensor2(F_aim,active=.true.)),& [9,grid(1),grid(2),grid3]) if (guess) then F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & @@ -349,8 +351,8 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc F_lambda33 = math_mul3333xx33(S_scale,matmul(F_lambda33, & math_mul3333xx33(C_scale,& matmul(transpose(F_lambda33),& - F_lambda33)-math_I3))*0.5_pReal)& - + math_I3 + F_lambda33)-math_I3))*0.5_pReal) & + + math_I3 F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) enddo; enddo; enddo endif @@ -514,7 +516,7 @@ subroutine formResidual(in, FandF_tau, & trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) + ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' deformation gradient aim =', transpose(F_aim) flush(6) @@ -533,7 +535,7 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! doing convolution in Fourier space call utilities_FFTtensorForward - call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) + call utilities_fourierGammaConvolution(params%rotation_BC%rotTensor2(polarBeta*F_aim,active=.true.)) call utilities_FFTtensorBackward !-------------------------------------------------------------------------------------------------- @@ -544,14 +546,14 @@ subroutine formResidual(in, FandF_tau, & ! evaluate constitutive response call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) P_av,C_volAvg,C_minMaxAvg, & - F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC) + F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC%asMatrix()) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- ! stress BC handling F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & - -math_rotate_forward33(F_av,params%rotation_BC)) + & + -params%rotation_BC%rotTensor2(F_av)) + & params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc ! calculate divergence tensorField_real = 0.0_pReal diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 5107bd900..7b826033b 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -104,7 +104,8 @@ module spectral_utilities end type tLoadCase type, public :: tSolutionParams !< @todo use here the type definition for a full loadcase - real(pReal), dimension(3,3) :: stress_mask, stress_BC, rotation_BC + real(pReal), dimension(3,3) :: stress_mask, stress_BC + type(rotation) :: rotation_BC real(pReal) :: timeinc real(pReal) :: timeincOld end type tSolutionParams diff --git a/src/math.f90 b/src/math.f90 index 6b83e17a3..bac2f8192 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1315,20 +1315,6 @@ pure function math_rotate_forward33(tensor,R) end function math_rotate_forward33 -!-------------------------------------------------------------------------------------------------- -!> @brief rotate 33 tensor backward -!> @details deprecated -!-------------------------------------------------------------------------------------------------- -pure function math_rotate_backward33(tensor,R) - - real(pReal), dimension(3,3) :: math_rotate_backward33 - real(pReal), dimension(3,3), intent(in) :: tensor, R - - math_rotate_backward33 = matmul(transpose(R),matmul(tensor,R)) - -end function math_rotate_backward33 - - !-------------------------------------------------------------------------------------------------- !> @brief rotate 3333 tensor C'_ijkl=g_im*g_jn*g_ko*g_lp*C_mnop !> @details deprecated diff --git a/src/rotations.f90 b/src/rotations.f90 index cd4a64547..a4a0bac88 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -163,7 +163,7 @@ subroutine fromQuaternion(self,qu) class(rotation), intent(out) :: self real(pReal), dimension(4), intent(in) :: qu - if (dNeq(norm2(qu),1.0)) & + if (dNeq(norm2(qu),1.0_pReal)) & call IO_error(402,ext_msg='fromQuaternion') self%q = qu From f5292019e5820c0f7f9a3f43c342560e12e66368 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 2 Dec 2019 20:23:50 +0100 Subject: [PATCH 11/19] use rotation class --- src/grid/grid_mech_FEM.f90 | 7 ++-- src/grid/grid_mech_spectral_basic.f90 | 7 ++-- src/grid/grid_mech_spectral_polarisation.f90 | 7 ++-- src/grid/spectral_utilities.f90 | 26 ++++++++------- src/math.f90 | 35 -------------------- 5 files changed, 23 insertions(+), 59 deletions(-) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index be79af835..e87e02d09 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -207,8 +207,7 @@ subroutine grid_mech_FEM_init call utilities_updateCoords(F) call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 F, & ! target F - 0.0_pReal, & ! time increment - math_I3) ! no rotation of boundary condition + 0.0_pReal) ! time increment call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr) CHKERRQ(ierr) call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr) @@ -255,7 +254,7 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(rotation_BC%asMatrix(),stress_BC%maskLogical,C_volAvg) + S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) !-------------------------------------------------------------------------------------------------- ! set module wide available data params%stress_mask = stress_BC%maskFloat @@ -507,7 +506,7 @@ subroutine formResidual(da_local,x_local, & ! evaluate constitutive response call Utilities_constitutiveResponse(P_current,& P_av,C_volAvg,devNull, & - F,params%timeinc,params%rotation_BC%asMatrix()) + F,params%timeinc,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 74b56db61..36caabf90 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -173,8 +173,7 @@ subroutine grid_mech_spectral_basic_init call Utilities_updateCoords(reshape(F,shape(F_lastInc))) call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F - 0.0_pReal, & ! time increment - math_I3) ! no rotation of boundary condition + 0.0_pReal) ! time increment call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! deassociate pointer restartRead2: if (interface_restartInc > 0) then @@ -225,7 +224,7 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(rotation_BC%asMatrix(),stress_BC%maskLogical,C_volAvg) + S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg) !-------------------------------------------------------------------------------------------------- @@ -457,7 +456,7 @@ subroutine formResidual(in, F, & ! evaluate constitutive response call utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory) P_av,C_volAvg,C_minMaxAvg, & - F,params%timeinc,params%rotation_BC%asMatrix()) + F,params%timeinc,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 5c5a9f10c..07d9a5afc 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -187,8 +187,7 @@ subroutine grid_mech_spectral_polarisation_init call Utilities_updateCoords(reshape(F,shape(F_lastInc))) call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F - 0.0_pReal, & ! time increment - math_I3) ! no rotation of boundary condition + 0.0_pReal) ! time increment call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer restartRead2: if (interface_restartInc > 0) then @@ -241,7 +240,7 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old, !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(rotation_BC%asMatrix(),stress_BC%maskLogical,C_volAvg) + S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) if (num%update_gamma) then call utilities_updateGamma(C_minMaxAvg) C_scale = C_minMaxAvg @@ -546,7 +545,7 @@ subroutine formResidual(in, FandF_tau, & ! evaluate constitutive response call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) P_av,C_volAvg,C_minMaxAvg, & - F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC%asMatrix()) + F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 7b826033b..a1265556a 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -686,10 +686,11 @@ end function utilities_curlRMS !-------------------------------------------------------------------------------------------------- function utilities_maskedCompliance(rot_BC,mask_stress,C) - real(pReal), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance - real(pReal), intent(in) , dimension(3,3,3,3) :: C !< current average stiffness - real(pReal), intent(in) , dimension(3,3) :: rot_BC !< rotation of load frame - logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC + real(pReal), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance + real(pReal), intent(in), dimension(3,3,3,3) :: C !< current average stiffness + type(rotation), intent(in) :: rot_BC !< rotation of load frame + logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC + integer :: j, k, m, n logical, dimension(9) :: mask_stressVector real(pReal), dimension(9,9) :: temp99_Real @@ -707,7 +708,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) allocate (c_reduced(size_reduced,size_reduced), source =0.0_pReal) allocate (s_reduced(size_reduced,size_reduced), source =0.0_pReal) allocate (sTimesC(size_reduced,size_reduced), source =0.0_pReal) - temp99_Real = math_3333to99(math_rotate_forward3333(C,rot_BC)) + temp99_Real = math_3333to99(rot_BC%rotTensor4(C)) if(debugGeneral) then write(6,'(/,a)') ' ... updating masked compliance ............................................' @@ -836,12 +837,12 @@ end subroutine utilities_fourierTensorDivergence subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& F,timeinc,rotation_BC) - real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness - real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress - real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress - real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target - real(pReal), intent(in) :: timeinc !< loading time - real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame + real(pReal), intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness + real(pReal), intent(out), dimension(3,3) :: P_av !< average PK stress + real(pReal), intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress + real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target + real(pReal), intent(in) :: timeinc !< loading time + type(rotation), intent(in), optional :: rotation_BC !< rotation of load frame integer :: & @@ -863,7 +864,8 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& if (debugRotation) & write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',& transpose(P_av)*1.e-6_pReal - P_av = math_rotate_forward33(P_av,rotation_BC) + if(present(rotation_BC)) & + P_av = rotation_BC%rotTensor2(P_av) write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& transpose(P_av)*1.e-6_pReal flush(6) diff --git a/src/math.f90 b/src/math.f90 index bac2f8192..0b06c9186 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1301,41 +1301,6 @@ real(pReal) pure function math_areaTriangle(v1,v2,v3) end function math_areaTriangle -!-------------------------------------------------------------------------------------------------- -!> @brief rotate 33 tensor forward -!> @details deprecated -!-------------------------------------------------------------------------------------------------- -pure function math_rotate_forward33(tensor,R) - - real(pReal), dimension(3,3) :: math_rotate_forward33 - real(pReal), dimension(3,3), intent(in) :: tensor, R - - math_rotate_forward33 = matmul(R,matmul(tensor,transpose(R))) - -end function math_rotate_forward33 - - -!-------------------------------------------------------------------------------------------------- -!> @brief rotate 3333 tensor C'_ijkl=g_im*g_jn*g_ko*g_lp*C_mnop -!> @details deprecated -!-------------------------------------------------------------------------------------------------- -pure function math_rotate_forward3333(tensor,R) - - real(pReal), dimension(3,3,3,3) :: math_rotate_forward3333 - real(pReal), dimension(3,3), intent(in) :: R - real(pReal), dimension(3,3,3,3), intent(in) :: tensor - integer :: i,j,k,l,m,n,o,p - - math_rotate_forward3333 = 0.0_pReal - do i = 1,3;do j = 1,3;do k = 1,3;do l = 1,3 - do m = 1,3;do n = 1,3;do o = 1,3;do p = 1,3 - math_rotate_forward3333(i,j,k,l) = math_rotate_forward3333(i,j,k,l) & - + R(i,m) * R(j,n) * R(k,o) * R(l,p) * tensor(m,n,o,p) - enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo - -end function math_rotate_forward3333 - - !-------------------------------------------------------------------------------------------------- !> @brief limits a scalar value to a certain range (either one or two sided) ! Will return NaN if left > right From fecd4632b45f9949542306cfd311d9c001269550 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 4 Dec 2019 05:49:17 +0100 Subject: [PATCH 12/19] correct reporting of time --- python/damask/dadf5.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index d03617752..4790a1b38 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -333,8 +333,8 @@ class DADF5(): """Return information on all active datasets in the file.""" message = '' with h5py.File(self.filename,'r') as f: - for s,i in enumerate(self.iter_visible('increments')): - message+='\n{} ({}s)\n'.format(i,self.times[s]) + for i in self.iter_visible('increments'): + message+='\n{} ({}s)\n'.format(i,self.times[self.increments.index(i)]) for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']): for oo in self.iter_visible(o): message+=' {}\n'.format(oo) From 32ebcea207d23c10c76b8fa97a5fc1b0a281b500 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 4 Dec 2019 05:49:43 +0100 Subject: [PATCH 13/19] more tests --- python/tests/test_DADF5.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/python/tests/test_DADF5.py b/python/tests/test_DADF5.py index 5a6478f03..37c5a9e87 100644 --- a/python/tests/test_DADF5.py +++ b/python/tests/test_DADF5.py @@ -42,6 +42,22 @@ class TestDADF5: in_file = default.read_dataset(loc['sigma'],0) assert np.allclose(in_memory,in_file) + def test_add_determinant(self,default): + default.add_determinant('P') + loc = {'P': default.get_dataset_location('P'), + 'det(P)':default.get_dataset_location('det(P)')} + in_memory = np.linalg.det(default.read_dataset(loc['P'],0)) + in_file = default.read_dataset(loc['det(P)'],0) + assert np.allclose(in_memory,in_file) + + def test_add_norm(self,default): + default.add_norm('F',1) + loc = {'F': default.get_dataset_location('F'), + '|F|_1':default.get_dataset_location('|F|_1')} + in_memory = np.linalg.norm(default.read_dataset(loc['F'],0),ord=1,axis=(1,2),keepdims=True) + in_file = default.read_dataset(loc['|F|_1'],0) + assert np.allclose(in_memory,in_file) + def test_add_absolute(self,default): default.add_absolute('Fe') loc = {'Fe': default.get_dataset_location('Fe'), From 6c31e228553a3c453fba8905ed70f12a939307ea Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 4 Dec 2019 05:59:52 +0100 Subject: [PATCH 14/19] testing new style DADF5 file --- .../DADF5/12grains6x7x8_tensionY.hdf5 | Bin 1935952 -> 1935528 bytes python/tests/test_DADF5.py | 8 +++++++- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/python/tests/reference/DADF5/12grains6x7x8_tensionY.hdf5 b/python/tests/reference/DADF5/12grains6x7x8_tensionY.hdf5 index caafce4784e2e9e3f29f8bd6dcad4e338a238153..39c17fadb5cc8fe724606eda02d7a62ac7ee8c6b 100644 GIT binary patch delta 19434 zcmc&+dz@3%ozJ}|6K;~p%$-MOk{O1{$T|~ifsCv*tyMB08bmZuYom5?&{_r*W~joX zVkKa^j@mj&(|S_v+TivPL2Cxo8+`5H`k10@Gd?gX#f%n=57w}*>?qjodFI|b+_3() z@PXXl?|#qkobx;9{J!^h&P{T^*pQBIh|ZIq&>4ui9|!4{li1zxJxc-QcV30Rk0FMu0V|)+T-h+LxGF zMb##bu0_VMxg|Z{TFmX4Obw#RtNl9|jVcAL9IzdiBCL!CPVFXx3sVV_IPTN2{5wVR^V)C@P5X7X4K)qCiNX=slg zmtk`5C93x`T!{IF+f)5S{S{!iz&Un&)mZLFA61x!=Go)pze}X@E*mpt8nr>2+0rG{ z!KZHptP3TZODxzDHv3w4gjz0JDi_-sw|5AMGhS#C1GPIZeBpdiEY>9UOe{~7{v_Gd zjwZxC*~E_OyePR6KXrJ(tKWPt@vx&ql;rsMj6}&HiM{nB<26DeQ(d2Esm+SPbr?4d z;r@O^99Ux_>RF1L;C zrG6;>Ik9Jdm+-p2XM&XI&h@WOF%b8qT8I|*(hpOtn$C2N%Q87%L-l-Nq{nFij>|E* z^h$jEKd{H`s>TH0+9q}G?0C6Z6y(I$Umdi?UZtfo`x1MkPQkZx;?-tR%8CoZ$f0hz6mKNj~z48qklndEKTj0%Z=wO~54@`{)yL$agH+EZ zfN6jg=B-GUjg7ojE7;Y_Oz^{FmQG=!yd>c0W`|oikg&bI!QT5ewHlo;DOaus=Zc9x zP4v%~pP^Q{Z9kL4&r%IpeDq7+xHZAyVvFxD*0RxX|1kM#kF&5`s!q)RZiNop-){jtG3n*ozZu7UZn- z1}{9xWEYEed;yLz101I@7$5%>=Gq(F5-+NB*uXqd@WU0=V9>6Ko6sJSJH>hWZh*-e z0H4@i^}5gtkJSP?Z0#~wYApt6+aj>rEVRRn3}W~{mooaOBWeEsm1E{a?^`lz^29te zG@)ALYKH(bf2mrjA71#=9&n40g7*2WXLII*DKh|5UEUlo+^?(+g3o7&g;j9%N-!Ht zUJu~E-3r+1@>wh{VF-?_0xPY-mDDkPuwym2%HF%0iqY%Qm`a_>$ZCN2B?bR0=`S$@ydfHcX!_$NW-ZqVSNsrlDR#@0E}_V{ z5%qj>Bem>OI&NmG>12Cn-R>XG+8qVaqo z2|lnqIY=|UMVa{rNk2q1*sFcY z#5HE?LLmk}Oq*|>(g93$Ad@k}H4d?&b1?$}e?bE>a4rFj)>f&O44g|q!Fd{xG2so1 zYc&`n!kexXGUn04_UL^yw4_p1C~lpK{5O7I6`HK0@cmR(zdOus^Qm5jE96JJ_w@aqs(G;(lQU_m=PRj~73VnB+}7V!p_Qa5hbZ%05rwYL z(^VVnn!pH8kAhBR>jSvaF}Q1c)f-^)!7BJgEuiCpSxks3)`zP4=YYw-RNB2!2j1O=$OxpR3?lnWfcOp3lS(pLiZ zPGvOZZI<_}s<_#1=>nN^r1JVocwTkG-44s3ebhnA6lyJrCtmFo8!alRu+I%Q%a;BH zNe1E-X$p~s^9H4z`Ud5gTxzz^fbO={QADx9LLZqVR3YY!f$939NhW7!P(7buDQ18- zWlS~5IK!@~HSyy9PO(+frP2X6jJr_3V=ihp(}Ejvt|+$PcC4L$Oy7(?CJ%j~Z$=+0 z&YPhzhns=2Lz~(G`h<3-+6?kk5Ark6G}!{y3n68^-BQ~GieIxNt~~A$Jr+VB`wR#q ze-`p2CJ-xJ@ijE4qu(I1=L@U&8NPz}`5yScNwfoBSQX(ci!p(O&au!?j;g6hAi?W3 zAcGi(!u_y#uLfff;xK;z-)9zhdJ4<|{yi#M=|9JE0)u?@kMgi6OjFx-Vz7hXd`4)q zsUrhlywTCA)E&o*E&YWBp1s1n9%N0^l-B2j*#aE@yJcpxkuwVmr(u_+o)@gY?4E`m zn1VzuAz55ZE{Qh-7n1A13xf60@oCnD!0Ve zpax{%)o<1zLvLv?#=LPC=(9^|i%#@qe0;i+TQUt1hwE&*#{B0<@WsdLJH;MNyg-=| z2%{SxrmEJsOW)APr*cVU#kmsV5YIe%KcH-#Cm{|&&+|PJdhT4HFXm-(kcT%KPEC{N z>T#;M_0ULiw_Z2$5si=&+PpYB9CjW zuUg}>x%5XgV6nVMd)AC;C z%FH`9a>a3ZCQtsA>UrQE`-9%)B__+;;yJqArfMV-&Xk@s?v8JplsNm<_4dp*Y9RNZ z1h?N;)u{;E@$J%sQsUUYHS>qqZQc$ITt^_u9r|{RGP&?bJa-hv?yqv&{zwt|j6?=s zoXo%v6g0bb=}QhVIsF9H)BjDT7%R_Pk}caqHIS@ctrPXsZ0S_$cH*v>9+eXIx7XXF zk5W@9MN&2{wj<_HB+Ov*eMzg%*I_0{?xA|#f-p0{Tc8@meeAVc+!Aj)9NC3tLF{F9 z8q@tC^icMMG$;g=RZmKYH}bDYh&S3n?^6;c4(ICz#4!fPUq&m!zYk4tn7Oge3Oimw zwzGFzQ3$+yJ(AZUS@b*^&sE%d1A#Xs8kbQG6*UxP7-ZY1aMz>znk zZs1+1!Q%ZlW#CO|rD5XXxlaQGoz$gh$0bdjjYrkfA;{0mZ*?m;TH z=@z=riZ{;qAR_vFYGE%+0OzUUv$N{^f=XEB5pqlO9G0$yr|Msup!n8!_Dv52JM=V ziT8m-r(NF*7wkD`djmud*(jsnfg&l2q{I)|`sWByC^Xs8=J`R?or;ul3m9#)Q%=DK zPR1!m5U?F>UmrR{vgrbJ zbhekoTa7C5)}Bsr(4jgQZEgY21$Im^*3<>m)F@ncfgSCVU#t!UITC>*i|qY95r*Ow zmP;+7rFp}KF3~qE!{nT@;zB!m$Ve}y_s5pl=|e_wiJc@Wq096IBTP<)sh-ccC^Ntt zqA|!fx5{PDi$`!H!ZfiIMYLuf)-pIlW3KDp8E2$j0(Scm*A9Rq^q&3F?e70)}jPPR{(+{Z;&qjxG zY-Yo`BCOl!=wDE{nW}kF=|@!dve3pCR6e4|Y0RO84GZts+%~l+4IaTw4w9kRHE}X= zY`jz1y$Ls8+Tx&j$-u3SDJ(LnoOY`NZFCJh;=pSQ4|1=vdxvAg)LmhGBmy#zJH8x= z+>U_lz!X4>d?lDAq?M&_;)W#mBR@Xk+>BlpFqu}Zz2+a;r=NXxpF#EAU!$LA_}C`~ zBK#Q`h(6~s++2qHM~vf+^}~rp_A@t`O+lo@e|(IpIG{{j>_wIkg_-&yGP&Zo7?U$M zQGJ-s>1Q52Bf{mproAR25zh zV9Y7&xifswLAW!2kd~cskOp; z&qmK?pG@hNx8%S%n>Agcw#Sr?5^9Q2ge(5b(ccq(g}yn&tB5i=^%~XlrbJnBj#JIm zed0hAs_6pRIZ>AWH(pGoqYgr$pwCgWbMc{51x1>46dQr5wuYo@Lz50drFnyNtyEjX zs1)?p;Y*TMfP6hs(?c^ZsC z*eBE%zCi;r2>XOigX1-nPDA&rXk}AljWL}D-%wGNkBl?}1n=F1_WU|tsX>k&eXWKP zYEBt>t!9n@dfRGg#%AQ&$+lYT&(h11Rct4h_+=DO@ZjqNc5M#EUq`MDpF%X)wXS{0 z%WITSzUH3uKtQUkNUpu7)>5{s;ABFH#Y;6H11}TMz>hT`11}TMADpI%zGaH7LtZ-vq?e(W=qSS`m8~J!m|=@ynKCRkaI!RQ6#}+fy&IW3pm4M7 zRk|a%HMyH)vIkh^FSyHMxyqhq*s&j>-k?)v@h9>Wt|_8BBu$GlD0h;kC}9~d&X?zC7GOklInRpAKIiJoD7q5n~_gwnx;{agQ6Yvda5aHmI>7iY@w!NaO)Pi3&^*t zlmufz$qgrdEcf>W8HzWoTL-Oin5BpU60S`Dc1RImD)0c1VUB57^eu zkLwEtncSD4dOqDk%m8nQ#vq$+)^eS2Hv|HQcML`+%%fQ1#jRRl=)@?^UEoO)Y2f#uCT~_86(+>Hh z+(LL`vY-JOggpW()z+z&48k5Em%^GlN-o87RkZTw5%~mi+5c-WEUs5UMyc9I0%)_T z!$JawM_D&ehx{A)*K1H~p4pGbGuDZ(95_`N3hY+|vxF<8j(XLSfoG3j(KsxIG?+?! z;)-fO#(d&?rVjDR*uyl+6i4de#s|=$Ku4w?`Bd6bPx&Op?(;dX!I$|Yx0Bi^Jz7ur z1Z{lvJXYVoAfxCbTb{|2yQrS8t28Um>(JQ4kC9|gOf5<&%k_9Yr4voKq*QX zKUq)NB=B_o6gCk*p-KCY*1VpWbDt~chA=%@8xj!SkSU^{w)W9ts*9FT8m zP-%f|`GLOX20|$1LlsE$yi@8~lo zACp6uQ$3H{p6U7&EXd@(84Wmj@f1*!tg6T>&f>*RJ);cGXdq`xQaiQfMa{q0(by@T z+Moh)_7d!9AlXY{hQ3lwLDC7{xKgJc=jx)e|1{J%lD*jFnN;OFwSjFoqe^l5lAWb5 zTVS%uNA-N(6;ILQGzRg!!yC>vwOVq5Eq--Mzrv?AkRzx3Y1G!_=?!dxA|%PBsvC9W zl3{YLhU$4^GRy!UU5!EBVNE|sDp@>VZ47-~!Q*S7=}hwc2Aa;K&uBpDj4#+g(wS^` zLq*E*EJB{8rzEip@g=1k=^~;=(Gke@Myw_W?$ISgjd>^HLzk6#CwXy$$`%HDS%gwz zKhl5*0t@eWPQ#G!ED zawnyx{1whA)YN;Yvh7MIqNew1^@hKAIrkiD+ov!(PjcJ-jY7tDv#D6zO{ZdLkG^b# z$;p?fo=?RnGr-%DRt@s0sHv9bq<_-%3!hPm&p8Q44Lnb6$HrJ=GGk7fh&|Ax8;l^6 z`?ga(Z%mLG;EmB3#2pS@V0YWpevrp8@m)=0@E2u3_9Zj{(JOVp_pXz2)a0L?h@-~- z>LeT$`VS}Js5cC@+_DGBvcPGP_dR?`Px4c6#X)B~-WLMtKQv)A;V9P!1dW5Iycb8N_q7;fXcvC<5u_or=81MhEnH_d9BwsR`k|TUS z&UBbs^L*S)54%KU4)ax^qz{tP5h{CGu;lA9|B)W2F^JbCIqp)6Qoag)(#(##)X0mE zyTndcdE($#K!9|#nH_xnPcx+}JYr;DZn8`He6ktQ)qvN9=L((xy+@l*N>`cLE=pI} zO6AQtvl8~XlI@@rLf>ga$6r9Y)HRi^k|7tNtMJ8?Eg~VfZJBGiePAh-W8ZR-JyP)s zO;JKu8FYuxRggjX@Mwza=*k~!(M8B!CP#lx^_;Fk-=5r4dq9kyoFG z$T!95H>ojQWxwMhbmb3IQ(@G4^4*K&QwAn-g^SRY%z(U3O+g)9d3yD23o^OyYScEG zQU$(CA-&68Y}Xc5if3Z+`}(r-N-FzS(N8`TUCaP)i^d?Hi6o~~tECz5wNv^vbhV4n zm2?fY?Y-8;rYu6@#$*?|b%Rl0vMEUQyfKBV^^MjT-ZgFdW9)yB}U8Cr!VAfhXA zY?HG3I#*2-=-uohL{(hvsz_9!Eo4F9U8YeuzS-4IlPUjCh#C`B;5;`cD)hQY#}-#R zSw87i2vM1+#%^<|9AxlpPdI4wF%8HlYJP4R77uDL20@eXRCHhLtyW9l@4MY) z{C;10u?AW(q0$T)O=TCk%Tb4=vOzbcE8k)_qpQ(jcV)T?$FM&MU5&+wYe`BKNDwtn zsnE-P7rEQXqk(|l{afs&9*+bOg1_vy#Eo9|OWowINL!=Otum6q01&bXzMuga4FG<_ z{HO+FFaX3*xU)sA&49JY8|L5EV2nwt?^kX_TCtzdSe3c0h4|d}Q%&T-kv?}@3ni^S zhUBYQQ<9Jt{_4c-)L8b;7P>g|cOt83aYlZokCZ7Uj}B8kUz`KX0Iwsd8sv*pQCH3TyC?mgYxnvHIo@?RMHTm=lI_;qf6^X@VRu`DS3J@(mA`z4Ti|v39oKBCGWC&whF-H^RK~rc zQX7kf4Ci5Agyp#M7h_k^FIf&S+}H>D(I_!lK1}sKj9ZBQW!X2SxWNMd<`nz>Q+e3f zuYhn|`By>Ew^?o>zl5Py0q-caEVqX~rd16bX(7)SlNvt|GL9k75Hj{MB)6?8NxI`D z|D)f?nB7LMI4;2CULVy{S0QJ9qVIf|$&pXdkne@Qb!dN2Rngq3v#lJ5Nb;Ecw1u#- z>nOG5N6jbSyf{8(#)T(Q_t5JZrBC&hY6_A@^TwI}rpI}i95l7^`S1ChLITHIlvWG+ z

c~ir3hfxmD+Y5|ia=RL=)%@-uxM8iS~gyKSidYOOR~1FdTRaZKnlwGuKenwV+O zO4p(nA#q!h@9Nt!%H+a9eOpGE0X}FNgS@|*ZqVK4HnlNO{+q=uPHUxHoNiMtw6-EH z4z#xtE*_IwD{`^VgFH%@0ev-mq`j4Kv1=w#V=nf*gD)vF#Gv45rOYri4XG=0ar)F& zl@txAfN*i@QVqy}3J559s|IAi1%!-;CR!;O$98DYgp5O7t;S?L^u7jKk&HD%Mu=hT zJ&oT<&Tcx?O37G$zm<`3YS#41WbB_kox3~@H^0wa?&&6KOvdQvEj~bM^42~2^$P?c zk3Zo<#O3toElYnAOkig zpwUk>AY(cXo;4lO@xV+PtzxIP!R-Nir}z~$^pI0{5#1o~bOafaFC9(E5xnO~_FcND zwd84SbbV$|L#2o0wewHc2g@jvQ?sa^$BGa$z$=NU2KoBb6lGtlnbk)5SQCLsPP0Po z3dtd(zBWq7#qKse9a{&krtKiVbdsgGVg0UDJgqbwS5^ps!HxeXe@^@#`SZ&3Ll}w< zw;kP620{QngsHs_dNssk|GN#b4Tp}I%p}BE(xx2SLp0*Q6nqF>kGSYn6;tjku~6BM zZnjt!$Gq}UH5Edwpa^Gmo^|#~C#?IAtFKwv-xu3*s6JLMm$%HLf6p!`7dmJ)Q7Fg% zarElmRk1B{^K8}@zGH@p-6J2hi1kSaDZ^@rqH{YhSSkoLk{}REnfb!HxhMKhoOeRI z|CD(rtT_3~d6$3Ne@c7IQL^4TRRYarW6_c}DQu73aHw_X4kazd;?KJ!r7ZT0rEi1( zDFAVGAtn$09?HG&AtZJiaTz8Hd!XFUaKRPO0j_V5$zv;_Jiu^SnYyY^88WmPRwygWKR+s4O~QaL$YA;x;DcT2g|9Jgax&3rYsvAj;o{eH*Y zk1~1a!n^Y29&n8h0Y_oPg##EbFY#0Y3A5Z0qi}9h1-Dd%_DCS!wC( zWcRr4WF8B3xD69~jW(p1oZA1HJDQFsSpy$9N)@lP~6yc=b=n2v2ZI3Nah zcn($8$=QG4VSkdv8(aEGiwlq3Z(VTSva>|Nl8=Sy*zhN9a+e@Rr<+jEU9uu3W7D`~ zDoIP?q@*L(@{cyje<(WLhXT!wuQ&{`0-e5DzFi8IpjnbgX9v2>A;k+Yd#5WliSsnb zpDLt80p#)1ESpfHB({(lh89R-A>b!Jl`j$fs8<%-rH-hkYuUKG5(TZ|iR$b!@#Evw zp)xTgsqT4V|7`_LJoJLHNu1JZzVfY-f<|JeA81negF)@T5!>X*lOhTHH1Z+p;R&?n zx8i`DcuG_oKee1FBF~fJ(Wc(bpzbHJieQD5ckdM1LW&aAw21mQfDV%);dj{OtC9DL zlEIx*;!$VGz7Iw8#6-z;1-;n2sB9CSR1!GK#M!fcQL;%q6P@mm;CYxhdsc;XyA--c z!1f;QDRyUsgU^ z;vo*^sW4M&#OV}9r{BMGjqFi)P4%!<#*PZ3_Huc^u8WoE_!?)cw7}R4*;gr}TOKTZ zK^#R-SIY3XnwKc&E1l5dAQDT&W;EoKN69~KHI;MBi12STGx7iA#(2#CvJ<;tashtQ zLzVc+IYcz{m(mpl!}`|ilnvx1| ztK|NbVRoAjIUSsi7fzg^2TG~TaYIZVe-O$!X+`>sqic-G3AARl{JK>Wd;73^y0(;R z6Cw934(4b}L>{n)*68 zSH@h|vli|T4ahr|zk;U}Zg2Dc)n_V0#(I7x2RA}Fr?3uIpWC6+ho=+@1YI75H;LM> z#XU;iP|zbBW7p%3`>qGgUFzZWZ0g;BJDR=}-c5M}!etWb21L>%) z18C_3rLAhq5E%;*bm~TVK+LL>8|BJ6aR}`i!^2Yfo*X;(lji#Gn9yisvy3BrZ2S{N z%%IM@aF0Uw(2^W3=|NFalKLqv3H(D5r_}v-$>)g$Gk_4OoE?`@lO#6FxEV7hy;T%% zmotCTflR_X1*JdLflT5%<M`$p}x}aCLjp0 zzLEyJ?qQ%K8{}5$@{;_ZCIsq1*4`5QhN2Glmdp{czb3z3ahfPEnb4oa7aa+OZZF+ET%ZFLp-)-p*MZFVbCw#9mn7e|$M1pAk{T~X4{lXj)h+jw;=!6j zzZ@^^_l&T+e5mKqYUh^HQ$%axUZ|Z`hwm-L2G5{f|6SVOId;D>jwP9#nSgRWlSWx} zu3l5O=4YjNCh6+oka3122@jy%6Q%3yy1+=|ZCG9N05Lm>Zho-zMcKEl6z!>$;fcIk zDR+x;^~Hxu`{&5MeWm2#;|D9dA1W=z8kcm+H^|}FO3`0BnUXJNj%a3|dMbuDqwu)m zm)$mtku~I?#S(k-5C36x4#FaVADz11d6`|g9z@>_PV|^4&$evGy6=FE<$01(N3iUa^be#>J=eP)yN4ub$hw>gf)Z@4j zCQmq68oUhk+*|6g<%=ye&`Doa@^dnViAJ7sSYr0)*4n!&KwlVl;7F6RRo(Bf;9PCQ zX^B1k{OW~XcAGb+gHt->X8G%kW`&qM=w@aO(3iJhHNYaZ7CLJ1Yl2v7>0dZA%UIpp z1ahDmKKZEhvHIK`O&@Ot`TySFv1+YU>^N7)k9Fv`&8m-8=jW|sqw}N!mrgZ}Cc3Ik_lP&p2K=TT@r`dw$Zp&*4k`hq@MNhZ1%gNn3EXnND;2D+` zrfifD%QORvvXSng#WwG4l|IvX-C5S?kv_E`dvD2FIUts)Z~j(1L6lRA%24BzVoSZ5 zzr$yjF;0@xi^(PRwj5kSi?Oriz}ILo_4U+`Tguv`kx#+ZB(+ohVj0fu;{GzT+>WwT z>fAEmmLi^~+!Fdv9mvG0UOm!+WGRWCZ1vS3q=<7Ta zZt%KsCZ`rbIs7XKLrek3P2?3uQ)IQsZOzwITh3IlcXFtAkuqS-M_ROh5$-|gTCE3Y z^m|SZYUxs5l-=hjz+2=QQQjI?*DO)`=Sc4L(bQ~0TavP5GF4%k>rmDT4J*J}M8!X6 zt?3~JHOpi7DQNx*`2;zBzf!?P@JlsirFDE|#bIN!bx_i1pO zt}bOKoN(W-^sma@2RQCDH?Ma;l>1oDoO#t|$bOUqAo~u&C#Qj#|1s+GvSrr5)@TZr z4Yd-5Cedougo0yp&sK%fQD!SB5Vk3?r>9%h>8<3hmq?h666`KFSBEyR(bRUM89639 zAHib*B}Lg|`#B(#NjBVklL|?B7LNYXST(`q^kb|7^GK2@;HG5r3fXY8YmF*s`Dg=X z*6pM({znwbLh(mHVM=8dB1b!1H*}qs$8x zR&*H4a`sLu361_VF6ev~1A>o7Q^mO^@{$6f95+#vYryL1U%1G-85hCwq3XTI ziaDlhkF_Yrgx|12KriBPNOC+j4P@eRNOG3wKqfAyAmLgas2ItNOz1$L_i}W4mfG=# z6_ZRy8p2EniD?^l#y(r#)pYrA5Ouy}>u+!}7@raOYA@kUo!PK_7}yJ{diL4W`HfDz z*0Ilq9{jq!uQLNTc;zl8hgs<3)6UHlaGa(v)@%O&d81B`!mFxP85{P{2-HelBU zO_CYOIhWnojV@Gm(DtGnKV$={A&-huO^4ch$ksncjH79n9d7{~#QpgKqRBPcfoO<= ze@-;JFuYmz2kbyJU1!>jL_-#(7T8x<$Ih_Rz?48O_bH2GkrGU)0VqnL>43dWcAlfF zN_a-jEwrP}SIMitGs2MEZTXrDL4%LQ)xiz04a;%kOim6%IcJ*Dv!EWwC77H(+a8_X zL&KjkG^K$I zwMfTR=Ih3BB~NwRDdlvX4~oNIx5GQu2u6~1zL~F-8z{$n63Bs5;FDY8WeT_@I)%Ki z`A(oU-^3##)UVkBn<;o`(J8thTFX*G!U zlcJ${cvH3sk+3M|){@(F{U95{m#d5IUCGLTT_YBg{S=zS2>0tiCb5u$W}ef5OkyFS zA2~a>93AmijKl zpWHfRd&SWvt#@dyCfmWNeZbM*IKfbSf@rlIlXwFQ#b32~Y=fX_<~j$=Lic(HT6#=P z-@@ySp)bQ^frUUm3*$@yZzQcL7W0=W8-k1|(&OL{?gcFig7O}Qi!j#WxIBw3fM6J5IOQdy_x(%`8qWpI z60!rZ637+i7Y?s;MyG!k3sfrBd|gx(uv4bAkM0BY@r(n=|8-p-DOIL24xq{~gYpJ+ zRjI8$W^t&=Xvh$g$D5&?GgRbNW6#EzoOsQFsgj&9uxTYURgqlI{OG2uh3VHEgfgw6 zgP>{=Sw+T4E@Aa5R85Ov zv#PVA0&iB;&=RuxKJ@{1esT89)>ddNW)j9JF=o!xflR_U1*NanflR_UWzVTp1+Zso zi-uM|tF0)qWtG-YU$AARJt4YfrEYn)0(&~}=L)(ui1(`6pRbrB%9&Y}#WS|Srb_PH z>}+xgzU}eSV%V~(N@L%qTWAUK?bx5On^&u;bj96g$;q+CszMt= zlM_8XQM%Tyi=`}gXqipAkZg7WMfzHtMv9c;-RcipoPBsB>wHWE^aV+Dt1aEmb%p!B zu!j|L{yP-;R}60!C(yKPnIn#?lY_*X;B`*Ach%#np#f+JP1o04Rvz+zIM`Stpz_~{ zO|~W|2sPJ4r=LJwD`c-6_`XvkSxP`%c&m%%YC^%Dm2NN0KtL0ZRvY48ipjahp`1r~ zUxB`?wwFe7HDV6l%bngpy5!sFq{-RDMtGW2-*dw1S#%@0o&BDZ63h_0%^$7~PEF!X zkN(g&SQ1Q5--@3RWM}ao03@ZZ`GFJ9j~tpDarSqHqQ;s-OdgLxIUi;ttU9+uQ`dX5 zvr@EaEuvK$oo&*nj#>$2&fMyx$=Tp-pm%H&8}ErtFy4KUDub7UOpaU)<=l`UQ@{<; zDdb(v_W^ZAo%sE4t7Z!t|8s>04ab}VVu$)s%!yO8)PqjSFQFehX=>JIvK*L1y<4#+ z|If*gr|}Pu(MeJi@2et3G&S>1@);q!@4%PE$)#hv6UYTOX0|X$PV`Cpbr_SVN4dr4 ztkM*jL_Nwa?!GGE7ViQbnq+1>w>!-fGxs_jvS?zKe@3V(QiZeN?^fmCr+60JeRb7< zJifMyC1(Cds~|Bm^YRpG|1MV3@i=uUP0ZqGm=ujCH81;qO|BsUz8YV$#4NM63MXdH zpH&s*8vk#qfNP34oiaq?unuJ6bqWf!RBM_{yiP$=K^>?V*93=jAakyn*;R$PMui&lGFO0!TtGA=Os#1(*DiO|_Jhq1c=G$2fl+d-a z8m+(7-Pfra>pK%5`&rQA(=fBch*R>{`S6LNk5=PD zl5cArffL71R0H=6{;C?dhdeK`Kqw3E;#)1w*=pyns+Y|X2T}CxYP?AQ1n$}wtP>{Q zuLcStdj9!3;V_0b%M;f*Axq0O)@TGm58#cX-ujvqw)l4VW<=^?wCZ8krS^%PAWtP+ zX!*URvz?EEJn=Yu;((w6CygGZQ%N~#3^xLtwB8tw^&CgSqmRZ&L72G^l=3L(4S^iE z7(O{M_zLuO=qVC7NzN{FX+wz|niy=Zq2x5v41DT7s)i8M8Qa6%b=r{_KWvpoXN>bDCf3} zGX)%{Q^@

jrxI(vntbv1Sa)^^|y|_8LfwMi$pVT9jz7!D*3qNe!i~iMcgJY0HUP z7S}Y>v}kB4b}@T(pZOZOj*(4$p}aOvUy<`@F?oLvx8uv=^fj`qMq?e5w{}WRLpSI^ zCXuwU*5WK_=)&6K+~r1xoweZVk&CEL34h7La&k>vjbFGrt7`Ek^(C}~c)7!g zFIVrZuEj}G&%)ZG#5HtTtwuK{?xu93Y|?>D-2K+Csy0c#(qT;8O%=J{(Sgi)Ym?OE z(kwEktM9Ug?JRfF@_yW{+bKF;EpbHL5}2O06s-@6(&iJM9w+cxz?$^A=sSBtfv} zVi)k$)N+@RucY9u>Z2FCaIQ3UtBVW+^3W3(yS{KuQHiT?Z+5W-;GV?hLw(kzYs1)L_n?UGLXZnJbSl zdGtdl=WI1tpl?JYbe_-29dFm?!v-8UeB4F3v+D^Mnqh1&IdCZbZ^5t``L&C3r{b0>7D}J9lMz^q22hs2p>OOxV6sU{v4B?+n}5qlY7c& zv`!)Ku&y8U%WKj>%@_#cJ%4mj_LS1EsmJYcRp7&ihg_65Xa3}(%ZX-*NYCqZx{_n> zDYXB9s~J+hF4X>-tC{j<@PWEKZzkW=CifrlWpUn2zva@X$s~|d-b^*tX&{q8PC--W z=s+fsoDyeZq7H~Nd94l&J92Nj%oD#vTn8=6p1Nm7PapogEV&`i|6LCLMnR zqEqs%dbD}5y;VA_i#%jk3zySX$3X_2w8ur39%L^m$7x`K;dmc(^$b)t)oZjZ$7JUf*tYlK;NrIcqV%|7@c3|}gvv(4 z#<~e6r&qBGv?W1EG6me0tfr6;P+haI9k`;PU;Zobv&1Lq<5z(0)RoMb+?8NV&ti|^ zVT>?&VhNP<*)qZuaAR}|Nr&k^fYuM#oEz4R0V%GtiDsO>XlvFnF$pqB~>- zUz{Rm8*SS2LNCjfR%xk63q<7Ts64gA14Neonui~6jCz2`%p5KJi1-nXH^!(-X&RWl zgBBA{rx7h?M;zCjn+T^^vPB2K5&lbh0kOq@DA^DWi z)tuztkq+Df+IntlV22!UYoLc57jA+ZymA+l!z_66dFf^fI8LXK%u8~}(WX^}Kr|A? zrjxxhUF4;~J@?y!zI5Hu04x^&aRaazd1PcEQ|fq9eIwR@6T|qE4P?BKSDEoU8orRZ za(`$5<|11D*)rx|)UKx+Fmny&8i2V*|K4C^E-C#EI_~X;D>sa?h}pq**SItCyKVK5 zR*f^<2wV5$IBgvi-q+&{mtJO!)H;Twh$P?^3X!*rPiu_HiFcr!Kdo{25^VmbQSv|J zLcWD^lU-MhGH~m=4Gs309O!WWy#Zb|KGwN`GH?fW3}xUHgYtgqs!~gRV6Tvz(?qX0 zF3V(PE|f#B1tI;Ov0K9XAbUT=j~p0SR|LsaTT?#%&}=gBA%1-LcSdXh>QS4yZ4;jvZA&pZCpH2*dTd7kk^M&ljC%#QNeDP% zX*6_yjL8WF%JWZ)Dd46|X$m<3>sp2Mt+X-Ueb@?+8sW!BB+wl%X=Fb>Qqo9&d}QH4 zqcLG7kGyU)Cd?FYV@5TF#29imNb3gduQcmyp;VAo8i9X<&5gjnQ%WP|->kQh{uYU3 zYb?sYY1BTuu^FXKvRI{~v5VP4Z1PL=+IuoCK0q(&K4M2Jy^fG#E50nwz`?f0JQbQT zF!_BF>6E--aXL*jf2YF~qu=oRI*>X2j-NCO)35s&2vHqvP3Xb*u)~kn#T!DuA$)_0b0s=X#fBK diff --git a/python/tests/test_DADF5.py b/python/tests/test_DADF5.py index 37c5a9e87..1dd93d003 100644 --- a/python/tests/test_DADF5.py +++ b/python/tests/test_DADF5.py @@ -23,7 +23,13 @@ def reference_dir(reference_dir_base): class TestDADF5: - + + def test_time_increments(self,default): + shape = default.read_dataset(default.get_dataset_location('F'),0).shape + default.set_by_time(0.0,20.0) + for i in default.iter_visible('increments'): + assert shape == default.read_dataset(default.get_dataset_location('F'),0).shape + def test_add_deviator(self,default): default.add_deviator('P') loc = {'P' :default.get_dataset_location('P'), From 6902d3dd0add631fda7dede02030bcfd59c94a4b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 4 Dec 2019 06:04:06 +0100 Subject: [PATCH 15/19] repeated test better sort alphabetically --- python/tests/test_DADF5.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/python/tests/test_DADF5.py b/python/tests/test_DADF5.py index 1dd93d003..2235321b5 100644 --- a/python/tests/test_DADF5.py +++ b/python/tests/test_DADF5.py @@ -72,14 +72,6 @@ class TestDADF5: in_file = default.read_dataset(loc['|Fe|'],0) assert np.allclose(in_memory,in_file) - def test_add_determinant(self,default): - default.add_determinant('P') - loc = {'P': default.get_dataset_location('P'), - 'det(P)': default.get_dataset_location('det(P)')} - in_memory = np.linalg.det(default.read_dataset(loc['P'],0)).reshape(-1,1) - in_file = default.read_dataset(loc['det(P)'],0) - assert np.allclose(in_memory,in_file) - def test_add_spherical(self,default): default.add_spherical('P') loc = {'P': default.get_dataset_location('P'), From 285075bb6b35a0e6083223109fa0b056f19f8196 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 4 Dec 2019 06:15:32 +0100 Subject: [PATCH 16/19] ordered alphabetically --- python/damask/dadf5.py | 332 ++++++++++++++++++------------------- python/tests/test_DADF5.py | 33 ++-- 2 files changed, 183 insertions(+), 182 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 4790a1b38..5c23636db 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -437,6 +437,76 @@ class DADF5(): else: with h5py.File(self.filename,'r') as f: return f['geometry/x_c'][()] + + + def add_absolute(self,x): + """ + Add absolute value. + + Parameters + ---------- + x : str + Label of the dataset containing a scalar, vector, or tensor. + + """ + def __add_absolute(x): + + return { + 'data': np.abs(x['data']), + 'label': '|{}|'.format(x['label']), + 'meta': { + 'Unit': x['meta']['Unit'], + 'Description': 'Absolute value of {} ({})'.format(x['label'],x['meta']['Description']), + 'Creator': 'dadf5.py:add_abs v{}'.format(version) + } + } + + requested = [{'label':x,'arg':'x'}] + + self.__add_generic_pointwise(__add_absolute,requested) + + + def add_calculation(self,formula,label,unit='n/a',description=None,vectorized=True): + """ + Add result of a general formula. + + Parameters + ---------- + formula : str + Formula, refer to datasets by ‘#Label#‘. + label : str + Label of the dataset containing the result of the calculation. + unit : str, optional + Physical unit of the result. + description : str, optional + Human readable description of the result. + vectorized : bool, optional + Indicate whether the formula is written in vectorized form. Default is ‘True’. + + """ + if vectorized is not True: + raise NotImplementedError + + def __add_calculation(**kwargs): + + formula = kwargs['formula'] + for d in re.findall(r'#(.*?)#',formula): + formula = formula.replace('#{}#'.format(d),"kwargs['{}']['data']".format(d)) + + return { + 'data': eval(formula), + 'label': kwargs['label'], + 'meta': { + 'Unit': kwargs['unit'], + 'Description': '{} (formula: {})'.format(kwargs['description'],kwargs['formula']), + 'Creator': 'dadf5.py:add_calculation v{}'.format(version) + } + } + + requested = [{'label':d,'arg':d} for d in set(re.findall(r'#(.*?)#',formula))] # datasets used in the formula + pass_through = {'formula':formula,'label':label,'unit':unit,'description':description} + + self.__add_generic_pointwise(__add_calculation,requested,pass_through) def add_Cauchy(self,P='P',F='F'): @@ -470,6 +540,90 @@ class DADF5(): self.__add_generic_pointwise(__add_Cauchy,requested) + def add_determinant(self,x): + """ + Add the determinant of a tensor. + + Parameters + ---------- + x : str + Label of the dataset containing a tensor. + + """ + def __add_determinant(x): + + return { + 'data': np.linalg.det(x['data']), + 'label': 'det({})'.format(x['label']), + 'meta': { + 'Unit': x['meta']['Unit'], + 'Description': 'Determinant of tensor {} ({})'.format(x['label'],x['meta']['Description']), + 'Creator': 'dadf5.py:add_determinant v{}'.format(version) + } + } + + requested = [{'label':x,'arg':'x'}] + + self.__add_generic_pointwise(__add_determinant,requested) + + + def add_deviator(self,x): + """ + Add the deviatoric part of a tensor. + + Parameters + ---------- + x : str + Label of the dataset containing a tensor. + + """ + def __add_deviator(x): + + if not np.all(np.array(x['data'].shape[1:]) == np.array([3,3])): + raise ValueError + + return { + 'data': mechanics.deviatoric_part(x['data']), + 'label': 's_{}'.format(x['label']), + 'meta': { + 'Unit': x['meta']['Unit'], + 'Description': 'Deviator of tensor {} ({})'.format(x['label'],x['meta']['Description']), + 'Creator': 'dadf5.py:add_deviator v{}'.format(version) + } + } + + requested = [{'label':x,'arg':'x'}] + + self.__add_generic_pointwise(__add_deviator,requested) + + + def add_maximum_shear(self,x): + """ + Add maximum shear components of symmetric tensor. + + Parameters + ---------- + x : str + Label of the dataset containing a symmetric tensor. + + """ + def __add_maximum_shear(x): + + return { + 'data': mechanics.maximum_shear(x['data']), + 'label': 'max_shear({})'.format(x['label']), + 'meta': { + 'Unit': x['meta']['Unit'], + 'Description': 'Maximum shear component of of {} ({})'.format(x['label'],x['meta']['Description']), + 'Creator': 'dadf5.py:add_maximum_shear v{}'.format(version) + } + } + + requested = [{'label':x,'arg':'x'}] + + self.__add_generic_pointwise(__add_maximum_shear,requested) + + def add_Mises(self,x): """ Add the equivalent Mises stress or strain of a symmetric tensor. @@ -540,58 +694,33 @@ class DADF5(): self.__add_generic_pointwise(__add_norm,requested,{'ord':ord}) - def add_absolute(self,x): + def add_principal_components(self,x): """ - Add absolute value. - + Add principal components of symmetric tensor. + + The principal components are sorted in descending order, each repeated according to its multiplicity. + Parameters ---------- x : str - Label of the dataset containing a scalar, vector, or tensor. + Label of the dataset containing a symmetric tensor. """ - def __add_absolute(x): + def __add_principal_components(x): return { - 'data': np.abs(x['data']), - 'label': '|{}|'.format(x['label']), + 'data': mechanics.principal_components(x['data']), + 'label': 'lambda_{}'.format(x['label']), 'meta': { 'Unit': x['meta']['Unit'], - 'Description': 'Absolute value of {} ({})'.format(x['label'],x['meta']['Description']), - 'Creator': 'dadf5.py:add_abs v{}'.format(version) + 'Description': 'Pricipal components of {} ({})'.format(x['label'],x['meta']['Description']), + 'Creator': 'dadf5.py:add_principal_components v{}'.format(version) } } requested = [{'label':x,'arg':'x'}] - - self.__add_generic_pointwise(__add_absolute,requested) - - def add_determinant(self,x): - """ - Add the determinant of a tensor. - - Parameters - ---------- - x : str - Label of the dataset containing a tensor. - - """ - def __add_determinant(x): - - return { - 'data': np.linalg.det(x['data']), - 'label': 'det({})'.format(x['label']), - 'meta': { - 'Unit': x['meta']['Unit'], - 'Description': 'Determinant of tensor {} ({})'.format(x['label'],x['meta']['Description']), - 'Creator': 'dadf5.py:add_determinant v{}'.format(version) - } - } - - requested = [{'label':x,'arg':'x'}] - - self.__add_generic_pointwise(__add_determinant,requested) + self.__add_generic_pointwise(__add_principal_components,requested) def add_spherical(self,x): @@ -624,79 +753,6 @@ class DADF5(): self.__add_generic_pointwise(__add_spherical,requested) - def add_deviator(self,x): - """ - Add the deviatoric part of a tensor. - - Parameters - ---------- - x : str - Label of the dataset containing a tensor. - - """ - def __add_deviator(x): - - if not np.all(np.array(x['data'].shape[1:]) == np.array([3,3])): - raise ValueError - - return { - 'data': mechanics.deviatoric_part(x['data']), - 'label': 's_{}'.format(x['label']), - 'meta': { - 'Unit': x['meta']['Unit'], - 'Description': 'Deviator of tensor {} ({})'.format(x['label'],x['meta']['Description']), - 'Creator': 'dadf5.py:add_deviator v{}'.format(version) - } - } - - requested = [{'label':x,'arg':'x'}] - - self.__add_generic_pointwise(__add_deviator,requested) - - - def add_calculation(self,formula,label,unit='n/a',description=None,vectorized=True): - """ - Add result of a general formula. - - Parameters - ---------- - formula : str - Formula, refer to datasets by ‘#Label#‘. - label : str - Label of the dataset containing the result of the calculation. - unit : str, optional - Physical unit of the result. - description : str, optional - Human readable description of the result. - vectorized : bool, optional - Indicate whether the formula is written in vectorized form. Default is ‘True’. - - """ - if vectorized is not True: - raise NotImplementedError - - def __add_calculation(**kwargs): - - formula = kwargs['formula'] - for d in re.findall(r'#(.*?)#',formula): - formula = formula.replace('#{}#'.format(d),"kwargs['{}']['data']".format(d)) - - return { - 'data': eval(formula), - 'label': kwargs['label'], - 'meta': { - 'Unit': kwargs['unit'], - 'Description': '{} (formula: {})'.format(kwargs['description'],kwargs['formula']), - 'Creator': 'dadf5.py:add_calculation v{}'.format(version) - } - } - - requested = [{'label':d,'arg':d} for d in set(re.findall(r'#(.*?)#',formula))] # datasets used in the formula - pass_through = {'formula':formula,'label':label,'unit':unit,'description':description} - - self.__add_generic_pointwise(__add_calculation,requested,pass_through) - - def add_strain_tensor(self,F='F',t='U',m=0): """ Add strain tensor calculated from a deformation gradient. @@ -729,62 +785,6 @@ class DADF5(): requested = [{'label':F,'arg':'F'}] self.__add_generic_pointwise(__add_strain_tensor,requested,{'t':t,'m':m}) - - - def add_principal_components(self,x): - """ - Add principal components of symmetric tensor. - - The principal components are sorted in descending order, each repeated according to its multiplicity. - - Parameters - ---------- - x : str - Label of the dataset containing a symmetric tensor. - - """ - def __add_principal_components(x): - - return { - 'data': mechanics.principal_components(x['data']), - 'label': 'lambda_{}'.format(x['label']), - 'meta': { - 'Unit': x['meta']['Unit'], - 'Description': 'Pricipal components of {} ({})'.format(x['label'],x['meta']['Description']), - 'Creator': 'dadf5.py:add_principal_components v{}'.format(version) - } - } - - requested = [{'label':x,'arg':'x'}] - - self.__add_generic_pointwise(__add_principal_components,requested) - - - def add_maximum_shear(self,x): - """ - Add maximum shear components of symmetric tensor. - - Parameters - ---------- - x : str - Label of the dataset containing a symmetric tensor. - - """ - def __add_maximum_shear(x): - - return { - 'data': mechanics.maximum_shear(x['data']), - 'label': 'max_shear({})'.format(x['label']), - 'meta': { - 'Unit': x['meta']['Unit'], - 'Description': 'Maximum shear component of of {} ({})'.format(x['label'],x['meta']['Description']), - 'Creator': 'dadf5.py:add_maximum_shear v{}'.format(version) - } - } - - requested = [{'label':x,'arg':'x'}] - - self.__add_generic_pointwise(__add_maximum_shear,requested) def __add_generic_pointwise(self,func,datasets_requested,extra_args={}): diff --git a/python/tests/test_DADF5.py b/python/tests/test_DADF5.py index 2235321b5..8aa8ec174 100644 --- a/python/tests/test_DADF5.py +++ b/python/tests/test_DADF5.py @@ -28,14 +28,15 @@ class TestDADF5: shape = default.read_dataset(default.get_dataset_location('F'),0).shape default.set_by_time(0.0,20.0) for i in default.iter_visible('increments'): - assert shape == default.read_dataset(default.get_dataset_location('F'),0).shape + assert shape == default.read_dataset(default.get_dataset_location('F'),0).shape - def test_add_deviator(self,default): - default.add_deviator('P') - loc = {'P' :default.get_dataset_location('P'), - 's_P':default.get_dataset_location('s_P')} - in_memory = mechanics.deviatoric_part(default.read_dataset(loc['P'],0)) - in_file = default.read_dataset(loc['s_P'],0) + + def test_add_absolute(self,default): + default.add_absolute('Fe') + loc = {'Fe': default.get_dataset_location('Fe'), + '|Fe|': default.get_dataset_location('|Fe|')} + in_memory = np.abs(default.read_dataset(loc['Fe'],0)) + in_file = default.read_dataset(loc['|Fe|'],0) assert np.allclose(in_memory,in_file) def test_add_Cauchy(self,default): @@ -52,10 +53,18 @@ class TestDADF5: default.add_determinant('P') loc = {'P': default.get_dataset_location('P'), 'det(P)':default.get_dataset_location('det(P)')} - in_memory = np.linalg.det(default.read_dataset(loc['P'],0)) + in_memory = np.linalg.det(default.read_dataset(loc['P'],0)).reshape((-1,1)) in_file = default.read_dataset(loc['det(P)'],0) assert np.allclose(in_memory,in_file) + def test_add_deviator(self,default): + default.add_deviator('P') + loc = {'P' :default.get_dataset_location('P'), + 's_P':default.get_dataset_location('s_P')} + in_memory = mechanics.deviatoric_part(default.read_dataset(loc['P'],0)) + in_file = default.read_dataset(loc['s_P'],0) + assert np.allclose(in_memory,in_file) + def test_add_norm(self,default): default.add_norm('F',1) loc = {'F': default.get_dataset_location('F'), @@ -64,14 +73,6 @@ class TestDADF5: in_file = default.read_dataset(loc['|F|_1'],0) assert np.allclose(in_memory,in_file) - def test_add_absolute(self,default): - default.add_absolute('Fe') - loc = {'Fe': default.get_dataset_location('Fe'), - '|Fe|': default.get_dataset_location('|Fe|')} - in_memory = np.abs(default.read_dataset(loc['Fe'],0)) - in_file = default.read_dataset(loc['|Fe|'],0) - assert np.allclose(in_memory,in_file) - def test_add_spherical(self,default): default.add_spherical('P') loc = {'P': default.get_dataset_location('P'), From cb0d39eee6cb806e14eb2b3382a24611a4162ed5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 4 Dec 2019 19:00:56 +0100 Subject: [PATCH 17/19] not needed anymore --- src/plastic_nonlocal.f90 | 40 +++++++++++++++++----------------------- 1 file changed, 17 insertions(+), 23 deletions(-) diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90 index 5469515dc..5a3f8c173 100644 --- a/src/plastic_nonlocal.f90 +++ b/src/plastic_nonlocal.f90 @@ -24,13 +24,10 @@ module plastic_nonlocal implicit none private - real(pReal), parameter, private :: & + real(pReal), parameter :: & KB = 1.38e-23_pReal !< Physical parameter, Boltzmann constant in J/Kelvin - - integer, dimension(:,:), allocatable, target, public :: & - plastic_nonlocal_sizePostResult !< size of each post result output - character(len=64), dimension(:,:), allocatable, target, public :: & + character(len=64), dimension(:,:), allocatable :: & plastic_nonlocal_output !< name of each post result output @@ -54,18 +51,18 @@ module plastic_nonlocal mob_scr_neg = 4 !< mobile screw positive ! BEGIN DEPRECATES - integer, dimension(:,:,:), allocatable, private :: & + integer, dimension(:,:,:), allocatable :: & iRhoU, & !< state indices for unblocked density iRhoB, & !< state indices for blocked density iRhoD, & !< state indices for dipole density iV, & !< state indices for dislcation velocities iD !< state indices for stable dipole height - integer, dimension(:), allocatable, private, protected :: & + integer, dimension(:), allocatable :: & totalNslip !< total number of active slip systems for each instance !END DEPRECATED - real(pReal), dimension(:,:,:,:,:,:), allocatable, private :: & - compatibility !< slip system compatibility between me and my neighbors + real(pReal), dimension(:,:,:,:,:,:), allocatable :: & + compatibility !< slip system compatibility between me and my neighbors enum, bind(c) enumerator :: & @@ -93,7 +90,7 @@ module plastic_nonlocal gamma_ID end enum - type, private :: tParameters !< container type for internal constitutive parameters + type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & atomicVolume, & !< atomic volume Dsd0, & !< prefactor for self-diffusion coefficient @@ -143,19 +140,19 @@ module plastic_nonlocal interactionSlipSlip ,& !< coefficients for slip-slip interaction forestProjection_Edge, & !< matrix of forest projections of edge dislocations forestProjection_Screw !< matrix of forest projections of screw dislocations - real(pReal), dimension(:), allocatable, private :: & + real(pReal), dimension(:), allocatable :: & nonSchmidCoeff - real(pReal), dimension(:,:,:), allocatable, private :: & + real(pReal), dimension(:,:,:), allocatable :: & Schmid, & !< Schmid contribution nonSchmid_pos, & nonSchmid_neg !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws) integer :: & totalNslip - integer, dimension(:) ,allocatable , public:: & + integer, dimension(:) ,allocatable :: & Nslip,& colinearSystem !< colinear system to the active slip system (only valid for fcc!) - logical, private :: & + logical :: & shortRangeStressCorrection, & !< flag indicating the use of the short range stress correction by a excess density gradient term probabilisticMultiplication @@ -164,13 +161,13 @@ module plastic_nonlocal end type tParameters - type, private :: tNonlocalMicrostructure + type :: tNonlocalMicrostructure real(pReal), allocatable, dimension(:,:) :: & tau_pass, & tau_Back end type tNonlocalMicrostructure - type, private :: tNonlocalState + type :: tNonlocalState real(pReal), pointer, dimension(:,:) :: & rho, & ! < all dislocations rhoSgl, & @@ -196,16 +193,16 @@ module plastic_nonlocal v_scr_neg end type tNonlocalState - type(tNonlocalState), allocatable, dimension(:), private :: & + type(tNonlocalState), allocatable, dimension(:) :: & deltaState, & dotState, & state - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - type(tNonlocalMicrostructure), dimension(:), allocatable, private :: microstructure + type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure - integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & + integer(kind(undefined_ID)), dimension(:,:), allocatable :: & plastic_nonlocal_outputID !< ID of each post result output public :: & @@ -268,7 +265,6 @@ subroutine plastic_nonlocal_init allocate(deltaState(maxNinstances)) allocate(microstructure(maxNinstances)) - allocate(plastic_nonlocal_sizePostResult(maxval(phase_Noutput), maxNinstances), source=0) allocate(plastic_nonlocal_output(maxval(phase_Noutput), maxNinstances)) plastic_nonlocal_output = '' allocate(plastic_nonlocal_outputID(maxval(phase_Noutput), maxNinstances), source=undefined_ID) @@ -498,7 +494,6 @@ subroutine plastic_nonlocal_init if (outputID /= undefined_ID) then plastic_nonlocal_output(i,phase_plasticityInstance(p)) = outputs(i) - plastic_nonlocal_sizePostResult(i,phase_plasticityInstance(p)) = prm%totalNslip prm%outputID = [prm%outputID , outputID] endif @@ -524,7 +519,6 @@ subroutine plastic_nonlocal_init prm%totalNslip,0,0) plasticState(p)%nonlocal = .true. plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention - plasticState(p)%sizePostResults = sum(plastic_nonlocal_sizePostResult(:,phase_plasticityInstance(p))) totalNslip(phase_plasticityInstance(p)) = prm%totalNslip From 0f4f415c5be0ce182c4608eab98dae8f1471f847 Mon Sep 17 00:00:00 2001 From: Test User Date: Sat, 7 Dec 2019 13:47:55 +0100 Subject: [PATCH 18/19] [skip ci] updated version information after successful test of v2.0.3-1218-g5a6111ec --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index ae832e588..88e96b85c 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-1177-ga41871e2 +v2.0.3-1218-g5a6111ec From a6e636a1c3e12615c2b2f2a5e958c3c781a48733 Mon Sep 17 00:00:00 2001 From: Test User Date: Sat, 7 Dec 2019 19:50:17 +0100 Subject: [PATCH 19/19] [skip ci] updated version information after successful test of v2.0.3-1228-g3e269f04 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 88e96b85c..35ae2dc0f 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-1218-g5a6111ec +v2.0.3-1228-g3e269f04