diff --git a/VERSION b/VERSION index 99c96eda6..690c6cb1d 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-15-g3898fb2 +v2.0.0-43-ge39441f diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 53378a70d..71a5e4743 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -3569,11 +3569,7 @@ logical function crystallite_integrateStress(& maxticks external :: & -#if(FLOAT==8) dgesv -#elif(FLOAT==4) - sgesv -#endif !* be pessimistic crystallite_integrateStress = .false. @@ -3756,11 +3752,7 @@ logical function crystallite_integrateStress(& - math_Plain3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dT3333,dT_dFe3333),dFe_dLp3333)) dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine work = math_plain33to9(residuumLp) -#if(FLOAT==8) call dgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp -#elif(FLOAT==4) - call sgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp -#endif if (ierr /= 0_pInt) then #ifndef _OPENMP if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then @@ -3849,31 +3841,27 @@ logical function crystallite_integrateStress(& math_mul3333xx3333(dT_dFi3333, dFi_dLi3333))) & - math_Plain3333to99(math_mul3333xx3333(dLi_dFi3333, dFi_dLi3333)) work = math_plain33to9(residuumLi) -#if(FLOAT==8) call dgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li -#elif(FLOAT==4) - call sgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li -#endif - if (ierr /= 0_pInt) then + if (ierr /= 0_pInt) then #ifndef _OPENMP - if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then - write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLi inversion at el ip ipc ', & - el,mesh_element(1,el),ip,ipc - if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & - .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g)& - .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then - write(6,*) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLi',transpose(dRLi_dLi) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLi',transpose(math_Plain3333to99(dFe_dLi3333)) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFi_constitutive',transpose(math_Plain3333to99(dT_dFi3333)) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLi_dT_constitutive',transpose(math_Plain3333to99(dLi_dT3333)) - write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Li_constitutive',math_transpose33(Li_constitutive) - write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Liguess',math_transpose33(Liguess) - endif + if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then + write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLi inversion at el ip ipc ', & + el,mesh_element(1,el),ip,ipc + if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & + .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g)& + .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then + write(6,*) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLi',transpose(dRLi_dLi) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLi',transpose(math_Plain3333to99(dFe_dLi3333)) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFi_constitutive',transpose(math_Plain3333to99(dT_dFi3333)) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLi_dT_constitutive',transpose(math_Plain3333to99(dLi_dT3333)) + write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Li_constitutive',math_transpose33(Li_constitutive) + write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Liguess',math_transpose33(Liguess) endif -#endif - return endif +#endif + return + endif deltaLi = - math_plain9to33(work) endif diff --git a/code/lattice.f90 b/code/lattice.f90 index f5717e8e1..05a123125 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -2182,7 +2182,7 @@ pure function lattice_qDisorientation(Q1, Q2, struct) real(pReal), dimension(4) :: lattice_qDisorientation real(pReal), dimension(4), intent(in) :: & Q1, & ! 1st orientation - Q2 ! 2nd orientation + Q2 ! 2nd orientation integer(kind(LATTICE_undefined_ID)), optional, intent(in) :: & ! if given, symmetries between the two orientation will be considered struct diff --git a/code/math.f90 b/code/math.f90 index 2f76c1cb5..bf7460062 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -186,10 +186,6 @@ module math halton_seed_set, & i_to_halton, & prime - external :: & - dsyev, & - dgetrf, & - dgetri contains @@ -811,15 +807,13 @@ function math_invSym3333(A) integer(pInt), dimension(6) :: ipiv6 real(pReal), dimension(6,6) :: temp66_Real real(pReal), dimension(6) :: work6 + external :: & + dgetrf, & + dgetri temp66_real = math_Mandel3333to66(A) -#if(FLOAT==8) call dgetrf(6,6,temp66_real,6,ipiv6,ierr) call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr) -#elif(FLOAT==4) - call sgetrf(6,6,temp66_real,6,ipiv6,ierr) - call sgetri(6,temp66_real,6,ipiv6,work6,6,ierr) -#endif if (ierr == 0_pInt) then math_invSym3333 = math_Mandel66to3333(temp66_real) else @@ -847,13 +841,8 @@ subroutine math_invert(myDim,A, InvA, error) logical, intent(out) :: error invA = A -#if(FLOAT==8) call dgetrf(myDim,myDim,invA,myDim,ipiv,ierr) call dgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr) -#elif(FLOAT==4) - call sgetrf(myDim,myDim,invA,myDim,ipiv,ierr) - call sgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr) -#endif error = merge(.true.,.false., ierr /= 0_pInt) ! http://fortraninacworld.blogspot.de/2012/12/ternary-operator.html end subroutine math_invert @@ -1937,16 +1926,13 @@ subroutine math_eigenValuesVectorsSym(m,values,vectors,error) real(pReal), dimension(size(m,1)), intent(out) :: values real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: vectors logical, intent(out) :: error - integer(pInt) :: info real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f + external :: & + dsyev vectors = m ! copy matrix to input (doubles as output) array -#if(FLOAT==8) call dsyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info) -#elif(FLOAT==4) - call ssyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info) -#endif error = (info == 0_pInt) end subroutine math_eigenValuesVectorsSym @@ -2135,16 +2121,13 @@ function math_eigenvaluesSym(m) real(pReal), dimension(:,:), intent(in) :: m real(pReal), dimension(size(m,1)) :: math_eigenvaluesSym real(pReal), dimension(size(m,1),size(m,1)) :: vectors - integer(pInt) :: info real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f + external :: & + dsyev vectors = m ! copy matrix to input (doubles as output) array -#if(FLOAT==8) call dsyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info) -#elif(FLOAT==4) - call ssyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info) -#endif if (info /= 0_pInt) math_eigenvaluesSym = DAMASK_NaN end function math_eigenvaluesSym diff --git a/code/prec.f90 b/code/prec.f90 index 0f22cf649..efa0c76a3 100644 --- a/code/prec.f90 +++ b/code/prec.f90 @@ -4,9 +4,9 @@ !> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH -!> @brief setting precision for real and int type depending on makros "FLOAT" and "INT" +!> @brief setting precision for real and int type !> @details setting precision for real and int type and for DAMASK_NaN. Definition is made -!! depending on makros "FLOAT" and "INT" defined during compilation +!! depending on makro "INT" defined during compilation !! for details on NaN see https://software.intel.com/en-us/forums/topic/294680 !-------------------------------------------------------------------------------------------------- module prec @@ -18,18 +18,7 @@ module prec implicit none private -#if (FLOAT==4) -#if defined(Spectral) || defined(FEM) - SPECTRAL SOLVER AND OWN FEM DO NOT SUPPORT SINGLE PRECISION, STOPPING COMPILATION -#endif - integer, parameter, public :: pReal = 4 !< floating point single precition (was selected_real_kind(6,37), number with 6 significant digits, up to 1e+-37) -#ifdef __INTEL_COMPILER - real(pReal), parameter, public :: DAMASK_NaN = Z'7F800001' !< quiet NaN for single precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran) -#endif -#ifdef __GFORTRAN__ - real(pReal), parameter, public :: DAMASK_NaN = real(Z'7F800001', pReal) !< quiet NaN for single precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran) -#endif -#elif (FLOAT==8) +#if (FLOAT==8) integer, parameter, public :: pReal = 8 !< floating point double precision (was selected_real_kind(15,300), number with 15 significant digits, up to 1e+-300) #ifdef __INTEL_COMPILER real(pReal), parameter, public :: DAMASK_NaN = Z'7FF8000000000000' !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran) @@ -172,9 +161,9 @@ end subroutine prec_init !-------------------------------------------------------------------------------------------------- !> @brief figures out if a floating point number is NaN -! basically just a small wrapper, because gfortran < 4.9 does not have the IEEE module +! basically just a small wrapper, because gfortran < 5.0 does not have the IEEE module !-------------------------------------------------------------------------------------------------- -logical elemental function prec_isNaN(a) +logical elemental pure function prec_isNaN(a) implicit none real(pReal), intent(in) :: a @@ -187,4 +176,30 @@ logical elemental function prec_isNaN(a) #endif end function prec_isNaN + +!-------------------------------------------------------------------------------------------------- +!> @brief equality comparison for double precision +! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq +! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm +!-------------------------------------------------------------------------------------------------- +logical elemental pure function dEq(a,b,tol) + real(pReal), intent(in) :: a,b + real(pReal), intent(in), optional :: tol + real(pReal), parameter :: eps = 2.2204460492503131E-16 ! DBL_EPSILON in C + dEq = merge(.True., .False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) +end function dEq + + +!-------------------------------------------------------------------------------------------------- +!> @brief inequality comparison for double precision +! replaces "!=" but for certain (relative) tolerance. Counterpart to dEq +! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm +!-------------------------------------------------------------------------------------------------- +logical elemental pure function dNeq(a,b,tol) + real(pReal), intent(in) :: a,b + real(pReal), intent(in), optional :: tol + real(pReal), parameter :: eps = 2.2204460492503131E-16 ! DBL_EPSILON in C + dNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b]))) +end function dNeq + end module prec diff --git a/lib/damask/asciitable.py b/lib/damask/asciitable.py index 4abe4297d..f76cf6627 100644 --- a/lib/damask/asciitable.py +++ b/lib/damask/asciitable.py @@ -517,22 +517,27 @@ class ASCIItable(): # ------------------------------------------------------------------ def microstructure_read(self, - grid): + grid, + type = 'i', + strict = False): """read microstructure data (from .geom format)""" - N = grid.prod() # expected number of microstructure indices in data - microstructure = np.zeros(N,'i') # initialize as flat array + def datatype(item): + return int(item) if type.lower() == 'i' else float(item) + + N = grid.prod() # expected number of microstructure indices in data + microstructure = np.zeros(N,type) # initialize as flat array i = 0 while i < N and self.data_read(): items = self.data if len(items) > 2: - if items[1].lower() == 'of': items = [int(items[2])]*int(items[0]) - elif items[1].lower() == 'to': items = range(int(items[0]),1+int(items[2])) - else: items = map(int,items) - else: items = map(int,items) + if items[1].lower() == 'of': items = np.ones(datatype(items[0]))*datatype(items[2]) + elif items[1].lower() == 'to': items = np.arange(datatype(items[0]),1+datatype(items[2])) + else: items = map(datatype,items) + else: items = map(datatype,items) - s = min(len(items), N-i) # prevent overflow of microstructure array + s = min(len(items), N-i) # prevent overflow of microstructure array microstructure[i:i+s] = items[:s] - i += s + i += len(items) - return microstructure + return (microstructure, i == N and not self.data_read()) if strict else microstructure # check for proper point count and end of file diff --git a/lib/damask/test/test.py b/lib/damask/test/test.py index 523f25d60..826b190c1 100644 --- a/lib/damask/test/test.py +++ b/lib/damask/test/test.py @@ -22,11 +22,11 @@ class Test(): logger = logging.getLogger() logger.setLevel(0) - fh = logging.FileHandler('test.log') # create file handler which logs even debug messages + fh = logging.FileHandler('test.log') # create file handler which logs even debug messages fh.setLevel(logging.DEBUG) full = logging.Formatter('%(asctime)s - %(levelname)s: \n%(message)s') fh.setFormatter(full) - ch = logging.StreamHandler(stream=sys.stdout) # create console handler with a higher log level + ch = logging.StreamHandler(stream=sys.stdout) # create console handler with a higher log level ch.setLevel(logging.INFO) # create formatter and add it to the handlers plain = logging.Formatter('%(message)s') @@ -60,10 +60,10 @@ class Test(): try: self.postprocess(variant) if not self.compare(variant): - return variant+1 # return culprit + return variant+1 # return culprit except Exception as e : - logging.critical('\nWARNING:\n %s\n'%e) - return variant+1 # return culprit + logging.critical('\nWARNING:\n {}\n'.format(e)) + return variant+1 # return culprit return 0 else: if not self.testPossible(): return -1 @@ -74,13 +74,13 @@ class Test(): self.prepare(variant) self.run(variant) self.postprocess(variant) - if self.updateRequested: # update requested + if self.updateRequested: # update requested self.update(variant) - elif not (self.options.accept or self.compare(variant)): # no update, do comparison - return variant+1 # return culprit + elif not (self.options.accept or self.compare(variant)): # no update, do comparison + return variant+1 # return culprit except Exception as e : - logging.critical('\nWARNING:\n %s\n'%e) - return variant+1 # return culprit + logging.critical('\nWARNING:\n {}\n'.format(e)) + return variant+1 # return culprit return 0 def testPossible(self): @@ -94,13 +94,13 @@ class Test(): try: shutil.rmtree(self.dirCurrent()) except: - logging.warning('removal of directory "%s" not possible...'%(self.dirCurrent())) + logging.warning('removal of directory "{}" not possible...'.format(self.dirCurrent())) status = status and False try: os.mkdir(self.dirCurrent()) except: - logging.critical('creation of directory "%s" failed...'%(self.dirCurrent())) + logging.critical('creation of directory "{}" failed...'.format(self.dirCurrent())) status = status and False return status @@ -193,19 +193,19 @@ class Test(): try: shutil.copy2(self.fileInReference(file),self.fileInCurrent(targetfiles[i])) except: - logging.critical('Reference2Current: Unable to copy file %s'%file) + logging.critical('Reference2Current: Unable to copy file "{}"'.format(file)) def copy_Base2Current(self,sourceDir,sourcefiles=[],targetfiles=[]): - source=os.path.normpath(os.path.join(self.dirBase,'../../../'+sourceDir)) + source=os.path.normpath(os.path.join(self.dirBase,'../../..',sourceDir)) if len(targetfiles) == 0: targetfiles = sourcefiles for i,file in enumerate(sourcefiles): try: shutil.copy2(os.path.join(source,file),self.fileInCurrent(targetfiles[i])) except: logging.error(os.path.join(source,file)) - logging.critical('Base2Current: Unable to copy file %s'%file) + logging.critical('Base2Current: Unable to copy file "{}"'.format(file)) def copy_Current2Reference(self,sourcefiles=[],targetfiles=[]): @@ -215,7 +215,7 @@ class Test(): try: shutil.copy2(self.fileInCurrent(file),self.fileInReference(targetfiles[i])) except: - logging.critical('Current2Reference: Unable to copy file %s'%file) + logging.critical('Current2Reference: Unable to copy file "{}"'.format(file)) def copy_Proof2Current(self,sourcefiles=[],targetfiles=[]): @@ -225,7 +225,7 @@ class Test(): try: shutil.copy2(self.fileInProof(file),self.fileInCurrent(targetfiles[i])) except: - logging.critical('Proof2Current: Unable to copy file %s'%file) + logging.critical('Proof2Current: Unable to copy file "{}"'.format(file)) def copy_Current2Current(self,sourcefiles=[],targetfiles=[]): @@ -234,7 +234,7 @@ class Test(): try: shutil.copy2(self.fileInReference(file),self.fileInCurrent(targetfiles[i])) except: - logging.critical('Current2Current: Unable to copy file %s'%file) + logging.critical('Current2Current: Unable to copy file "{}"'.format(file)) def execute_inCurrentDir(self,cmd,streamIn=None): @@ -252,7 +252,7 @@ class Test(): def compare_Array(self,File1,File2): import numpy as np - logging.info('comparing\n '+File1+'\n '+File2) + logging.info('\n '.join(['comparing',File1,File2])) table1 = damask.ASCIItable(name=File1,readonly=True) table1.head_read() len1=len(table1.info)+2 @@ -270,8 +270,9 @@ class Test(): max_loc=np.argmax(abs(refArrayNonZero[curArray.nonzero()]/curArray[curArray.nonzero()]-1.)) refArrayNonZero = refArrayNonZero[curArray.nonzero()] curArray = curArray[curArray.nonzero()] - print(' ********\n * maximum relative error %e for %e and %e\n ********' - %(max_err, refArrayNonZero[max_loc],curArray[max_loc])) + print(' ********\n * maximum relative error {} between {} and {}\n ********'.format(max_err, + refArrayNonZero[max_loc], + curArray[max_loc])) return max_err else: raise Exception('mismatch in array size to compare') @@ -295,7 +296,7 @@ class Test(): absoluteTolerance=False,perLine=False,skipLines=[]): import numpy as np - logging.info('comparing ASCII Tables\n %s \n %s'%(file0,file1)) + logging.info('\n '.join(['comparing ASCII Tables',file0,file1])) if normHeadings == '': normHeadings = headings0 # check if comparison is possible and determine lenght of columns @@ -315,7 +316,7 @@ class Test(): for i in xrange(dataLength): if headings0[i]['shape'] != headings1[i]['shape']: - raise Exception('shape mismatch when comparing %s with %s '%(headings0[i]['label'],headings1[i]['label'])) + raise Exception('shape mismatch between {} and {} '.format(headings0[i]['label'],headings1[i]['label'])) shape[i] = headings0[i]['shape'] for j in xrange(np.shape(shape[i])[0]): length[i] *= shape[i][j] @@ -323,7 +324,9 @@ class Test(): for j in xrange(np.shape(normShape[i])[0]): normLength[i] *= normShape[i][j] else: - raise Exception('trying to compare %i with %i normed by %i data sets'%(len(headings0),len(headings1),len(normHeadings))) + raise Exception('trying to compare {} with {} normed by {} data sets'.format(len(headings0), + len(headings1), + len(normHeadings))) table0 = damask.ASCIItable(name=file0,readonly=True) table0.head_read() @@ -331,37 +334,34 @@ class Test(): table1.head_read() for i in xrange(dataLength): - key0 = {True :'1_%s', - False:'%s' }[length[i]>1]%headings0[i]['label'] - key1 = {True :'1_%s', - False:'%s' }[length[i]>1]%headings1[i]['label'] - normKey = {True :'1_%s', - False:'%s' }[normLength[i]>1]%normHeadings[i]['label'] + key0 = ('1_' if length[i]>1 else '') + headings0[i]['label'] + key1 = ('1_' if length[i]>1 else '') + headings1[i]['label'] + normKey = ('1_' if normLength[i]>1 else '') + normHeadings[i]['label'] if key0 not in table0.labels: - raise Exception('column %s not found in 1. table...\n'%key0) + raise Exception('column {} not found in 1. table...\n'.format(key0)) elif key1 not in table1.labels: - raise Exception('column %s not found in 2. table...\n'%key1) + raise Exception('column {} not found in 2. table...\n'.format(key1)) elif normKey not in table0.labels: - raise Exception('column %s not found in 1. table...\n'%normKey) + raise Exception('column {} not found in 1. table...\n'.format(normKey)) else: - column[0][i] = table0.labels.index(key0) # remember columns of requested data - column[1][i] = table1.labels.index(key1) # remember columns of requested data in second column - normColumn[i] = table0.labels.index(normKey) # remember columns of requested data in second column + column[0][i] = table0.labels.index(key0) + column[1][i] = table1.labels.index(key1) + normColumn[i] = table0.labels.index(normKey) line0 = 0 - while table0.data_read(): # read next data line of ASCII table + while table0.data_read(): # read next data line of ASCII table if line0 not in skipLines: for i in xrange(dataLength): myData = np.array(map(float,table0.data[column[0][i]:\ - column[0][i]+length[i]]),'d') + column[0][i]+length[i]]),'d') normData = np.array(map(float,table0.data[normColumn[i]:\ - normColumn[i]+normLength[i]]),'d') + normColumn[i]+normLength[i]]),'d') data[i] = np.append(data[i],np.reshape(myData,shape[i])) if normType == 'pInf': norm[i] = np.append(norm[i],np.max(np.abs(normData))) else: norm[i] = np.append(norm[i],np.linalg.norm(np.reshape(normData,normShape[i]),normType)) - line0 +=1 + line0 += 1 for i in xrange(dataLength): if not perLine: norm[i] = [np.max(norm[i]) for j in xrange(line0-len(skipLines))] @@ -370,12 +370,12 @@ class Test(): norm[i] = [1.0 for j in xrange(line0-len(skipLines))] absTol[i] = True if perLine: - logging.warning('At least one norm of %s in 1. table is 0.0, using absolute tolerance'%headings0[i]['label']) + logging.warning('At least one norm of {} in 1. table is 0.0, using absolute tolerance'.format(headings0[i]['label'])) else: - logging.warning('Maximum norm of %s in 1. table is 0.0, using absolute tolerance'%headings0[i]['label']) + logging.warning('Maximum norm of {} in 1. table is 0.0, using absolute tolerance'.format(headings0[i]['label'])) line1 = 0 - while table1.data_read(): # read next data line of ASCII table + while table1.data_read(): # read next data line of ASCII table if line1 not in skipLines: for i in xrange(dataLength): myData = np.array(map(float,table1.data[column[1][i]:\ @@ -384,21 +384,25 @@ class Test(): norm[i][line1-len(skipLines)]) line1 +=1 - if (line0 != line1): raise Exception('found %s lines in 1. table and %s in 2. table'%(line0,line1)) + if (line0 != line1): raise Exception('found {} lines in 1. table but {} in 2. table'.format(line0,line1)) logging.info(' ********') for i in xrange(dataLength): if absTol[i]: - logging.info(' * maximum absolute error %e for %s and %s'%(maxError[i],headings0[i]['label'],headings1[i]['label'])) + logging.info(' * maximum absolute error {} between {} and {}'.format(maxError[i], + headings0[i]['label'], + headings1[i]['label'])) else: - logging.info(' * maximum relative error %e for %s and %s'%(maxError[i],headings0[i]['label'],headings1[i]['label'])) + logging.info(' * maximum relative error {} between {} and {}'.format(maxError[i], + headings0[i]['label'], + headings1[i]['label'])) logging.info(' ********') return maxError def compare_TablesStatistically(self, - files = [None,None], # list of file names - columns = [None], # list of list of column labels (per file) + files = [None,None], # list of file names + columns = [None], # list of list of column labels (per file) meanTol = 1.0e-4, stdTol = 1.0e-6, preFilter = 1.0e-9): @@ -407,18 +411,18 @@ class Test(): threshold can be used to ignore small values (a negative number disables this feature) """ - if not (isinstance(files, Iterable) and not isinstance(files, str)): # check whether list of files is requested + if not (isinstance(files, Iterable) and not isinstance(files, str)): # check whether list of files is requested files = [str(files)] tables = [damask.ASCIItable(name = filename,readonly = True) for filename in files] for table in tables: table.head_read() - columns += [columns[0]]*(len(files)-len(columns)) # extend to same length as files - columns = columns[:len(files)] # truncate to same length as files + columns += [columns[0]]*(len(files)-len(columns)) # extend to same length as files + columns = columns[:len(files)] # truncate to same length as files for i,column in enumerate(columns): - if column is None: columns[i] = tables[i].labels # if no column is given, read all + if column is None: columns[i] = tables[i].labels # if no column is given, read all logging.info('comparing ASCIItables statistically') for i in xrange(len(columns)): @@ -428,7 +432,7 @@ class Test(): ) logging.info(files[i]+':'+','.join(columns[i])) - if len(files) < 2: return True # single table is always close to itself... + if len(files) < 2: return True # single table is always close to itself... data = [] for table,labels in zip(tables,columns): @@ -443,16 +447,16 @@ class Test(): normedDelta = np.where(normBy>preFilter,delta/normBy,0.0) mean = np.amax(np.abs(np.mean(normedDelta,0))) std = np.amax(np.std(normedDelta,0)) - logging.info('mean: %f'%mean) - logging.info('std: %f'%std) + logging.info('mean: {:f}'.format(mean)) + logging.info('std: {:f}'.format(std)) return (mean0.0, maximum, 1) # do not devide by zero for empty columns + maximum = np.where(maximum >0.0, maximum, 1) # avoid div by zero for empty columns for i in xrange(len(data)): data[i] /= maximum mask = np.zeros_like(table.data,dtype='bool') for table in data: - mask |= np.where(np.abs(table) 1])) + logging.critical(('The test' if len(self.variants) == 1 else 'All {} tests'.format(len(self.variants))) + ' passed') logging.critical('\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n') return 0 if culprit == -1: logging.warning('Warning: Could not start test') return 0 else: - logging.critical(' ********\n * Test %i failed...\n ********'%(culprit)) + logging.critical(' ********\n * Test {} failed...\n ********'.format(culprit)) logging.critical('\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n') return culprit diff --git a/lib/damask/util.py b/lib/damask/util.py index 77bc30403..244aaa25c 100644 --- a/lib/damask/util.py +++ b/lib/damask/util.py @@ -58,6 +58,24 @@ def emph(what): """emphasizes string on screen""" return bcolors.BOLD+srepr(what)+bcolors.ENDC +# ----------------------------- +def execute(cmd, + streamIn = None, + wd = './'): + """executes a command in given directory and returns stdout and stderr for optional stdin""" + initialPath = os.getcwd() + os.chdir(wd) + process = subprocess.Popen(shlex.split(cmd), + stdout = subprocess.PIPE, + stderr = subprocess.PIPE, + stdin = subprocess.PIPE) + out,error = [i.replace("\x08","") for i in (process.communicate() if streamIn is None + else process.communicate(streamIn.read()))] + os.chdir(initialPath) + if process.returncode != 0: raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode)) + return out,error + + # ----------------------------- # Matlab like trigonometric functions that take and return angles in degrees. # ----------------------------- @@ -75,7 +93,7 @@ def gridLocation(idx,res): # ----------------------------- def gridIndex(location,res): - return ( location[0] % res[0] + \ + return ( location[0] % res[0] + \ ( location[1] % res[1]) * res[0] + \ ( location[2] % res[2]) * res[1] * res[0] ) @@ -104,7 +122,9 @@ class extendableOption(Option): class backgroundMessage(threading.Thread): """reporting with animation to indicate progress""" - choices = {'bounce': ['_','o','O','°','¯','¯','°','O','o',], + choices = {'bounce': ['_', 'o', 'O', u'\u00B0', + u'\u203e',u'\u203e',u'\u00B0','O','o','_'], + 'spin': [u'\u25dc',u'\u25dd',u'\u25de',u'\u25df'], 'circle': [u'\u25f4',u'\u25f5',u'\u25f6',u'\u25f7'], 'hexagon': [u'\u2b22',u'\u2b23'], 'square': [u'\u2596',u'\u2598',u'\u259d',u'\u2597'], @@ -228,7 +248,7 @@ def leastsqBound(func, x0, args=(), bounds=None, Dfun=None, full_output=0, return shape(res), dt def _int2extGrad(p_int, bounds): - """Calculate the gradients of transforming the internal (unconstrained) to external (constained) parameter.""" + """Calculate the gradients of transforming the internal (unconstrained) to external (constrained) parameter.""" grad = np.empty_like(p_int) for i, (x, bound) in enumerate(zip(p_int, bounds)): lower, upper = bound @@ -430,22 +450,4 @@ def curve_fit_bound(f, xdata, ydata, p0=None, sigma=None, bounds=None, **kw): else: pcov = np.inf - if return_full: - return popt, pcov, infodict, errmsg, ier - else: - return popt, pcov - - -def execute(cmd,streamIn=None,wd='./'): - """executes a command in given directory and returns stdout and stderr for optional stdin""" - initialPath=os.getcwd() - os.chdir(wd) - process = subprocess.Popen(shlex.split(cmd),stdout=subprocess.PIPE,stderr = subprocess.PIPE,stdin=subprocess.PIPE) - if streamIn is not None: - out,error = [i.replace("\x08","") for i in process.communicate(streamIn.read())] - else: - out,error =[i.replace("\x08","") for i in process.communicate()] - os.chdir(initialPath) - if process.returncode !=0: raise RuntimeError(cmd+' failed with returncode '+str(process.returncode)) - return out,error - + return (popt, pcov, infodict, errmsg, ier) if return_full else (popt, pcov) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index ac9715c36..d8b1ee025 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -10,40 +10,35 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) def curlFFT(geomdim,field): - N = grid.prod() # field size - n = np.array(np.shape(field)[3:]).prod() # data size + grid = np.array(np.shape(field)[2::-1]) + N = grid.prod() # field size + n = np.array(np.shape(field)[3:]).prod() # data size - if n == 3: - dataType = 'vector' - elif n == 9: - dataType = 'tensor' + if n == 3: dataType = 'vector' + elif n == 9: dataType = 'tensor' field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2)) curl_fourier = np.zeros(field_fourier.shape,'c16') # differentiation in Fourier space k_s = np.zeros([3],'i') - TWOPIIMG = (0.0+2.0j*math.pi) + TWOPIIMG = 2.0j*math.pi for i in xrange(grid[2]): k_s[0] = i - if(grid[2]%2==0 and i == grid[2]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[0]=0 - elif (i > grid[2]//2): - k_s[0] = k_s[0] - grid[2] + if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif i > grid[2]//2: k_s[0] -= grid[2] for j in xrange(grid[1]): k_s[1] = j - if(grid[1]%2==0 and j == grid[1]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[1]=0 - elif (j > grid[1]//2): - k_s[1] = k_s[1] - grid[1] + if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif j > grid[1]//2: k_s[1] -= grid[1] for k in xrange(grid[0]//2+1): k_s[2] = k - if(grid[0]%2==0 and k == grid[0]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[2]=0 + if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + + xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field input order - xi = np.array([k_s[2]/geomdim[2]+0.0j,k_s[1]/geomdim[1]+0.j,k_s[0]/geomdim[0]+0.j],'c16') if dataType == 'tensor': for l in xrange(3): curl_fourier[i,j,k,0,l] = ( field_fourier[i,j,k,l,2]*xi[1]\ @@ -100,10 +95,8 @@ if options.vector is None and options.tensor is None: if filenames == []: filenames = [None] for name in filenames: - try: - table = damask.ASCIItable(name = name,buffered = False) - except: - continue + try: table = damask.ASCIItable(name = name,buffered = False) + except: continue damask.util.report(scriptName,name) # ------------------------------------------ read header ------------------------------------------ @@ -161,9 +154,10 @@ for name in filenames: stack = [table.data] for type, data in items.iteritems(): for i,label in enumerate(data['active']): - stack.append(curlFFT(size[::-1], # we need to reverse order here, because x - table.data[:,data['column'][i]:data['column'][i]+data['dim']]. # is fastest,ie rightmost, but leftmost in - reshape([grid[2],grid[1],grid[0]]+data['shape']))) # our x,y,z notation + # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation + stack.append(curlFFT(size[::-1], + table.data[:,data['column'][i]:data['column'][i]+data['dim']]. + reshape([grid[2],grid[1],grid[0]]+data['shape']))) # ------------------------------------------ output result ----------------------------------------- diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 75708a603..aadaceabf 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -10,39 +10,35 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) def divFFT(geomdim,field): - N = grid.prod() # field size - n = np.array(np.shape(field)[3:]).prod() # data size + grid = np.array(np.shape(field)[2::-1]) + N = grid.prod() # field size + n = np.array(np.shape(field)[3:]).prod() # data size field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2)) - div_fourier = np.zeros(field_fourier.shape[0:len(np.shape(field))-1],'c16') # size depents on whether tensor or vector + div_fourier = np.zeros(field_fourier.shape[0:len(np.shape(field))-1],'c16') # size depents on whether tensor or vector # differentiation in Fourier space k_s=np.zeros([3],'i') - TWOPIIMG = (0.0+2.0j*math.pi) + TWOPIIMG = 2.0j*math.pi for i in xrange(grid[2]): k_s[0] = i - if(grid[2]%2==0 and i == grid[2]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[0]=0 - elif (i > grid[2]//2): - k_s[0] = k_s[0] - grid[2] + if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif i > grid[2]//2: k_s[0] -= grid[2] for j in xrange(grid[1]): k_s[1] = j - if(grid[1]%2==0 and j == grid[1]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[1]=0 - elif (j > grid[1]//2): - k_s[1] = k_s[1] - grid[1] + if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif j > grid[1]//2: k_s[1] -= grid[1] for k in xrange(grid[0]//2+1): k_s[2] = k - if(grid[0]%2==0 and k == grid[0]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[2]=0 + if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - xi=np.array([k_s[2]/geomdim[2]+0.0j,k_s[1]/geomdim[1]+0.j,k_s[0]/geomdim[0]+0.j],'c16') - if n == 9: # tensor, 3x3 -> 3 + xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field input order + if n == 9: # tensor, 3x3 -> 3 for l in xrange(3): div_fourier[i,j,k,l] = sum(field_fourier[i,j,k,l,0:3]*xi) *TWOPIIMG - elif n == 3: # vector, 3 -> 1 + elif n == 3: # vector, 3 -> 1 div_fourier[i,j,k] = sum(field_fourier[i,j,k,0:3]*xi) *TWOPIIMG return np.fft.fftpack.irfftn(div_fourier,axes=(0,1,2)).reshape([N,n/3]) @@ -80,15 +76,13 @@ parser.set_defaults(coords = 'ipinitialcoord', if options.vector is None and options.tensor is None: parser.error('no data column specified.') -# --- loop over input files ------------------------------------------------------------------------- +# --- loop over input files ------------------------------------------------------------------------ if filenames == []: filenames = [None] for name in filenames: - try: - table = damask.ASCIItable(name = name,buffered = False) - except: - continue + try: table = damask.ASCIItable(name = name,buffered = False) + except: continue damask.util.report(scriptName,name) # ------------------------------------------ read header ------------------------------------------ @@ -140,16 +134,17 @@ for name in filenames: maxcorner = np.array(map(max,coords)) grid = np.array(map(len,coords),'i') size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1) - size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 set to smallest among other spacings + size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other ones # ------------------------------------------ process value field ----------------------------------- stack = [table.data] for type, data in items.iteritems(): for i,label in enumerate(data['active']): - stack.append(divFFT(size[::-1], # we need to reverse order here, because x - table.data[:,data['column'][i]:data['column'][i]+data['dim']]. # is fastest,ie rightmost, but leftmost in - reshape([grid[2],grid[1],grid[0]]+data['shape']))) # our x,y,z notation + # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation + stack.append(divFFT(size[::-1], + table.data[:,data['column'][i]:data['column'][i]+data['dim']]. + reshape([grid[2],grid[1],grid[0]]+data['shape']))) # ------------------------------------------ output result ----------------------------------------- @@ -158,4 +153,4 @@ for name in filenames: # ------------------------------------------ output finalization ----------------------------------- - table.close() # close input ASCII table (works for stdin) \ No newline at end of file + table.close() # close input ASCII table (works for stdin) diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py new file mode 100755 index 000000000..555e587de --- /dev/null +++ b/processing/post/addGradient.py @@ -0,0 +1,158 @@ +#!/usr/bin/env python +# -*- coding: UTF-8 no BOM -*- + +import os,sys,math +import numpy as np +from optparse import OptionParser +import damask + +scriptName = os.path.splitext(os.path.basename(__file__))[0] +scriptID = ' '.join([scriptName,damask.version]) + +def gradFFT(geomdim,field): + grid = np.array(np.shape(field)[2::-1]) + N = grid.prod() # field size + n = np.array(np.shape(field)[3:]).prod() # data size + if n == 3: dataType = 'vector' + elif n == 1: dataType = 'scalar' + + field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2)) + grad_fourier = np.zeros(field_fourier.shape+(3,),'c16') + +# differentiation in Fourier space + k_s = np.zeros([3],'i') + TWOPIIMG = 2.0j*math.pi + for i in xrange(grid[2]): + k_s[0] = i + if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif i > grid[2]//2: k_s[0] -= grid[2] + + for j in xrange(grid[1]): + k_s[1] = j + if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif j > grid[1]//2: k_s[1] -= grid[1] + + for k in xrange(grid[0]//2+1): + k_s[2] = k + if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + + xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field order + + grad_fourier[i,j,k,0,:] = field_fourier[i,j,k,0]*xi *TWOPIIMG # vector field from scalar data + + if dataType == 'vector': + grad_fourier[i,j,k,1,:] = field_fourier[i,j,k,1]*xi *TWOPIIMG # tensor field from vector data + grad_fourier[i,j,k,2,:] = field_fourier[i,j,k,2]*xi *TWOPIIMG + + return np.fft.fftpack.irfftn(grad_fourier,axes=(0,1,2)).reshape([N,3*n]) + + +# -------------------------------------------------------------------- +# MAIN +# -------------------------------------------------------------------- + +parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ +Add column(s) containing gradient of requested column(s). +Operates on periodic ordered three-dimensional data sets. +Deals with both vector- and scalar fields. + +""", version = scriptID) + +parser.add_option('-c','--coordinates', + dest = 'coords', + type = 'string', metavar='string', + help = 'column heading for coordinates [%default]') +parser.add_option('-v','--vector', + dest = 'vector', + action = 'extend', metavar = '', + help = 'heading of columns containing vector field values') +parser.add_option('-s','--scalar', + dest = 'scalar', + action = 'extend', metavar = '', + help = 'heading of columns containing scalar field values') + +parser.set_defaults(coords = 'ipinitialcoord', + ) + +(options,filenames) = parser.parse_args() + +if options.vector is None and options.scalar is None: + parser.error('no data column specified.') + +# --- loop over input files ------------------------------------------------------------------------ + +if filenames == []: filenames = [None] + +for name in filenames: + try: table = damask.ASCIItable(name = name,buffered = False) + except: continue + damask.util.report(scriptName,name) + +# ------------------------------------------ read header ------------------------------------------ + + table.head_read() + +# ------------------------------------------ sanity checks ---------------------------------------- + + items = { + 'scalar': {'dim': 1, 'shape': [1], 'labels':options.scalar, 'active':[], 'column': []}, + 'vector': {'dim': 3, 'shape': [3], 'labels':options.vector, 'active':[], 'column': []}, + } + errors = [] + remarks = [] + column = {} + + if table.label_dimension(options.coords) != 3: errors.append('coordinates {} are not a vector.'.format(options.coords)) + else: colCoord = table.label_index(options.coords) + + for type, data in items.iteritems(): + for what in (data['labels'] if data['labels'] is not None else []): + dim = table.label_dimension(what) + if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type)) + else: + items[type]['active'].append(what) + items[type]['column'].append(table.label_index(what)) + + if remarks != []: damask.util.croak(remarks) + if errors != []: + damask.util.croak(errors) + table.close(dismiss = True) + continue + +# ------------------------------------------ assemble header -------------------------------------- + + table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) + for type, data in items.iteritems(): + for label in data['active']: + table.labels_append(['{}_gradFFT({})'.format(i+1,label) for i in xrange(3 * data['dim'])]) # extend ASCII header with new labels + table.head_write() + +# --------------- figure out size and grid --------------------------------------------------------- + + table.data_readArray() + + coords = [np.unique(table.data[:,colCoord+i]) for i in xrange(3)] + mincorner = np.array(map(min,coords)) + maxcorner = np.array(map(max,coords)) + grid = np.array(map(len,coords),'i') + size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1) + size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) + +# ------------------------------------------ process value field ----------------------------------- + + stack = [table.data] + for type, data in items.iteritems(): + for i,label in enumerate(data['active']): + # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation + stack.append(gradFFT(size[::-1], + table.data[:,data['column'][i]:data['column'][i]+data['dim']]. + reshape([grid[2],grid[1],grid[0]]+data['shape']))) + +# ------------------------------------------ output result ----------------------------------------- + + if len(stack) > 1: table.data = np.hstack(tuple(stack)) + table.data_writeArray('%.12g') + +# ------------------------------------------ output finalization ----------------------------------- + + table.close() # close input ASCII table (works for stdin) diff --git a/processing/post/filterTable.py b/processing/post/filterTable.py index d375785dc..6723d9faf 100755 --- a/processing/post/filterTable.py +++ b/processing/post/filterTable.py @@ -35,7 +35,7 @@ Filter rows according to condition and columns by either white or black listing. Examples: Every odd row if x coordinate is positive -- " #ip.x# >= 0.0 and #_row_#%2 == 1 ). -All rows where label 'foo' equals 'bar' -- " #foo# == \"bar\" " +All rows where label 'foo' equals 'bar' -- " #s#foo# == 'bar' " """, version = scriptID) diff --git a/processing/post/vtk_addRectilinearGridData.py b/processing/post/vtk_addRectilinearGridData.py index 922b4b696..3ca54d84e 100755 --- a/processing/post/vtk_addRectilinearGridData.py +++ b/processing/post/vtk_addRectilinearGridData.py @@ -151,7 +151,7 @@ for name in filenames: writer = vtk.vtkXMLRectilinearGridWriter() writer.SetDataModeToBinary() writer.SetCompressorTypeToZLib() - writer.SetFileName(os.path.splitext(options.vtk)[0]+('' if options.inplace else '_added.vtr')) + writer.SetFileName(os.path.splitext(options.vtk)[0]+('.vtr' if options.inplace else '_added.vtr')) if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(rGrid) else: writer.SetInputData(rGrid) writer.Write() diff --git a/processing/pre/geom_canvas.py b/processing/pre/geom_canvas.py index 882f32c88..ec171f522 100755 --- a/processing/pre/geom_canvas.py +++ b/processing/pre/geom_canvas.py @@ -28,24 +28,32 @@ parser.add_option('-o', '--offset', help = 'a,b,c offset from old to new origin of grid [%default]') parser.add_option('-f', '--fill', dest = 'fill', - type = 'int', metavar = 'int', + type = 'float', metavar = 'float', help = '(background) canvas grain index. "0" selects maximum microstructure index + 1 [%default]') +parser.add_option('--float', + dest = 'real', + action = 'store_true', + help = 'input data is float [%default]') parser.set_defaults(grid = ['0','0','0'], offset = (0,0,0), fill = 0, + real = False, ) (options, filenames) = parser.parse_args() +datatype = 'f' if options.real else 'i' + + # --- loop over input files ------------------------------------------------------------------------- if filenames == []: filenames = [None] for name in filenames: - try: - table = damask.ASCIItable(name = name, - buffered = False, labeled = False) + try: table = damask.ASCIItable(name = name, + buffered = False, + labeled = False) except: continue damask.util.report(scriptName,name) @@ -71,7 +79,7 @@ for name in filenames: # --- read data ------------------------------------------------------------------------------------ - microstructure = table.microstructure_read(info['grid']).reshape(info['grid'],order='F') # read microstructure + microstructure = table.microstructure_read(info['grid'],datatype).reshape(info['grid'],order='F') # read microstructure # --- do work ------------------------------------------------------------------------------------ @@ -85,8 +93,8 @@ for name in filenames: else int(n) for o,n in zip(info['grid'],options.grid)],'i') newInfo['grid'] = np.where(newInfo['grid'] > 0, newInfo['grid'],info['grid']) - microstructure_cropped = np.zeros(newInfo['grid'],'i') - microstructure_cropped.fill(options.fill if options.fill > 0 else microstructure.max()+1) + microstructure_cropped = np.zeros(newInfo['grid'],datatype) + microstructure_cropped.fill(options.fill if options.real or options.fill > 0 else microstructure.max()+1) xindex = list(set(xrange(options.offset[0],options.offset[0]+newInfo['grid'][0])) & \ set(xrange(info['grid'][0]))) yindex = list(set(xrange(options.offset[1],options.offset[1]+newInfo['grid'][1])) & \ @@ -152,9 +160,9 @@ for name in filenames: # --- write microstructure information ------------------------------------------------------------ - formatwidth = int(math.floor(math.log10(microstructure_cropped.max())+1)) + format = '%g' if options.real else '%{}i'.format(int(math.floor(math.log10(microstructure_cropped.max())+1))) table.data = microstructure_cropped.reshape((newInfo['grid'][0],newInfo['grid'][1]*newInfo['grid'][2]),order='F').transpose() - table.data_writeArray('%%%ii'%(formatwidth),delimiter=' ') + table.data_writeArray(format,delimiter=' ') # --- output finalization -------------------------------------------------------------------------- diff --git a/processing/pre/geom_check.py b/processing/pre/geom_check.py index d607a1f73..518e8c3b0 100755 --- a/processing/pre/geom_check.py +++ b/processing/pre/geom_check.py @@ -54,12 +54,25 @@ for name in filenames: errors = [] if np.any(info['grid'] < 1): errors.append('invalid grid a b c.') if np.any(info['size'] <= 0.0): errors.append('invalid size x y z.') + +#--- read microstructure information -------------------------------------------------------------- + + if options.data: + microstructure,ok = table.microstructure_read(info['grid'],strict = True) # read microstructure + + if ok: + structure = vtk.vtkIntArray() + structure.SetName('Microstructures') + for idx in microstructure: structure.InsertNextValue(idx) + + else: errors.append('mismatch between data and grid dimension.') + if errors != []: damask.util.croak(errors) table.close(dismiss = True) continue -# --- generate VTK rectilinear grid -------------------------------------------------------------------------------- +# --- generate VTK rectilinear grid --------------------------------------------------------------- grid = vtk.vtkRectilinearGrid() grid.SetDimensions([x+1 for x in info['grid']]) @@ -72,18 +85,8 @@ for name in filenames: elif i == 1: grid.SetYCoordinates(temp) elif i == 2: grid.SetZCoordinates(temp) -#--- read microstructure information -------------------------------------------------------------- - if options.data: - microstructure = table.microstructure_read(info['grid']) # read microstructure - - structure = vtk.vtkIntArray() - structure.SetName('Microstructures') - - for idx in microstructure: - structure.InsertNextValue(idx) - - grid.GetCellData().AddArray(structure) + if options.data: grid.GetCellData().AddArray(structure) # --- write data ----------------------------------------------------------------------------------- if name: diff --git a/processing/pre/geom_unpack.py b/processing/pre/geom_unpack.py index a07166cdd..85c60d6e5 100755 --- a/processing/pre/geom_unpack.py +++ b/processing/pre/geom_unpack.py @@ -18,6 +18,14 @@ Unpack geometry files containing ranges "a to b" and/or "n of x" multiples (excl """, version = scriptID) +parser.add_option('-1', '--onedimensional', + dest = 'oneD', + action = 'store_true', + help = 'output geom file with one-dimensional data arrangement') + +parser.set_defaults(oneD = False, + ) + (options, filenames) = parser.parse_args() # --- loop over input files ------------------------------------------------------------------------- @@ -69,7 +77,8 @@ for name in filenames: microstructure = table.microstructure_read(info['grid']) # read microstructure formatwidth = int(math.floor(math.log10(microstructure.max())+1)) # efficient number printing format - table.data = microstructure.reshape((info['grid'][0],info['grid'][1]*info['grid'][2]),order='F').transpose() + table.data = microstructure if options.oneD else \ + microstructure.reshape((info['grid'][0],info['grid'][1]*info['grid'][2]),order='F').transpose() table.data_writeArray('%%%ii'%(formatwidth),delimiter = ' ') #--- output finalization --------------------------------------------------------------------------