From 2dd0e06b249f69dbe570df22019b11982c3f5ca3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 21 Mar 2016 19:42:28 +0100 Subject: [PATCH 01/18] should reflect status of plasticity modules --- config/Phase_DisloKMC_Tungsten.config | 37 ------------------- ... Phase_Isotropic_AluminumIsotropic.config} | 2 +- 2 files changed, 1 insertion(+), 38 deletions(-) delete mode 100644 config/Phase_DisloKMC_Tungsten.config rename config/{Phase_J2_AluminumIsotropic.config => Phase_Isotropic_AluminumIsotropic.config} (94%) diff --git a/config/Phase_DisloKMC_Tungsten.config b/config/Phase_DisloKMC_Tungsten.config deleted file mode 100644 index 799c80e61..000000000 --- a/config/Phase_DisloKMC_Tungsten.config +++ /dev/null @@ -1,37 +0,0 @@ -### $Id$ ### - -[Tungsten] -elasticity hooke -plasticity dislokmc - -### Material parameters ### -lattice_structure bcc -C11 523.0e9 # From Marinica et al. Journal of Physics: Condensed Matter(2013) -C12 202.0e9 -C44 161.0e9 - -grainsize 2.0e-5 # Average grain size [m] 2.0e-5 -SolidSolutionStrength 0.0 # Strength due to elements in solid solution - -### Dislocation glide parameters ### -#per family -Nslip 12 0 -slipburgers 2.72e-10 # Burgers vector of slip system [m] -rhoedge0 1.0e12 # Initial edge dislocation density [m/m**3] -rhoedgedip0 1.0 # Initial edged dipole dislocation density [m/m**3] -Qedge 2.725e-19 # Activation energy for dislocation glide [J] -v0 3560.3 # Initial glide velocity [m/s](kmC) -p_slip 0.16 # p-exponent in glide velocity -q_slip 1.00 # q-exponent in glide velocity -u_slip 2.47 # u-exponent of stress pre-factor (kmC) -s_slip 0.97 # self hardening in glide velocity (kmC) -tau_peierls 2.03e9 # peierls stress [Pa] - -#hardening -dipoleformationfactor 0 # to have hardening due to dipole formation off -CLambdaSlip 10.0 # Adj. parameter controlling dislocation mean free path -D0 4.0e-5 # Vacancy diffusion prefactor [m**2/s] -Qsd 4.5e-19 # Activation energy for climb [J] -Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b] -Cedgedipmindistance 1.0 # Adj. parameter controlling the minimum dipole distance [in b] -interaction_slipslip 0.2 0.11 0.19 0.15 0.11 0.17 diff --git a/config/Phase_J2_AluminumIsotropic.config b/config/Phase_Isotropic_AluminumIsotropic.config similarity index 94% rename from config/Phase_J2_AluminumIsotropic.config rename to config/Phase_Isotropic_AluminumIsotropic.config index 7a41010d5..88a8ad579 100644 --- a/config/Phase_J2_AluminumIsotropic.config +++ b/config/Phase_Isotropic_AluminumIsotropic.config @@ -3,7 +3,7 @@ # Kuo, J. C., Mikrostrukturmechanik von Bikristallen mit Kippkorngrenzen. Shaker-Verlag 2004. http://edoc.mpg.de/204079 elasticity hooke -plasticity j2 +plasticity isotropic (output) flowstress (output) strainrate From a165b7b68cb36c2c796882937bdab2beff08728a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 23 Mar 2016 15:12:35 +0100 Subject: [PATCH 02/18] simplified statement --- code/math.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/code/math.f90 b/code/math.f90 index bf7460062..624990ac4 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -1645,14 +1645,14 @@ pure function math_qToAxisAngle(Q) real(pReal) :: halfAngle, sinHalfAngle real(pReal), dimension(4) :: math_qToAxisAngle - halfAngle = acos(max(-1.0_pReal, min(1.0_pReal, Q(1)))) ! limit to [-1,1] --> 0 to 180 deg + halfAngle = acos(math_limit(Q(1),-1.0_pReal,1.0_pReal)) sinHalfAngle = sin(halfAngle) - if (sinHalfAngle <= 1.0e-4_pReal) then ! very small rotation angle? + smallRotation: if (sinHalfAngle <= 1.0e-4_pReal) then math_qToAxisAngle = 0.0_pReal - else + else smallRotation math_qToAxisAngle= [ Q(2:4)/sinHalfAngle, halfAngle*2.0_pReal] - endif + endif smallRotation end function math_qToAxisAngle From ce99a071c2109d69118bf01a00f8dec7d381cac0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 23 Mar 2016 18:17:47 +0100 Subject: [PATCH 03/18] cleaning --- code/plastic_isotropic.f90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/code/plastic_isotropic.f90 b/code/plastic_isotropic.f90 index 81d3055cf..d8dc0f384 100644 --- a/code/plastic_isotropic.f90 +++ b/code/plastic_isotropic.f90 @@ -45,7 +45,7 @@ module plastic_isotropic gdot0, & n, & h0, & - h0_slopeLnRate, & + h0_slopeLnRate=0.0_pReal, & tausat, & a, & aTolFlowstress, & @@ -466,7 +466,7 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,e math_spherical33, & math_mul33xx33 use material, only: & - phaseAt, phasememberAt, & + phasememberAt, & plasticState, & material_phase, & phase_plasticityInstance @@ -495,7 +495,7 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,e k, l, m, n of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember - instance = phase_plasticityInstance(phaseAt(ipc,ip,el)) ! "phaseAt" equivalent to "material_phase" !! + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) Tstar_sph_33 = math_spherical33(math_Mandel6to33(Tstar_v)) ! spherical part of 2nd Piola-Kirchhoff stress squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph_33,Tstar_sph_33) @@ -535,7 +535,7 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el) use math, only: & math_mul6x6 use material, only: & - phaseAt, phasememberAt, & + phasememberAt, & plasticState, & material_phase, & phase_plasticityInstance @@ -559,7 +559,7 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el) of !< shortcut notation for offset position in state array of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember - instance = phase_plasticityInstance(phaseAt(ipc,ip,el)) ! "phaseAt" equivalent to "material_phase" !! + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) !-------------------------------------------------------------------------------------------------- ! norm of (deviatoric) 2nd Piola-Kirchhoff stress @@ -616,7 +616,7 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el) use material, only: & material_phase, & plasticState, & - phaseAt, phasememberAt, & + phasememberAt, & phase_plasticityInstance implicit none @@ -640,7 +640,7 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el) o of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember - instance = phase_plasticityInstance(phaseAt(ipc,ip,el)) ! "phaseAt" equivalent to "material_phase" !! + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) !-------------------------------------------------------------------------------------------------- ! norm of (deviatoric) 2nd Piola-Kirchhoff stress From 04d2148e61f624f96a3b7ef49ab2518c3e34c47a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 24 Mar 2016 17:55:15 +0100 Subject: [PATCH 04/18] will replace geom_fromAng (in combination with geom_fromTable) --- processing/pre/table_fromOIMgrainFile.py | 51 ++++++++++++++++++++++++ 1 file changed, 51 insertions(+) create mode 100755 processing/pre/table_fromOIMgrainFile.py diff --git a/processing/pre/table_fromOIMgrainFile.py b/processing/pre/table_fromOIMgrainFile.py new file mode 100755 index 000000000..c730c8892 --- /dev/null +++ b/processing/pre/table_fromOIMgrainFile.py @@ -0,0 +1,51 @@ +#!/usr/bin/env python +# -*- coding: UTF-8 no BOM -*- + +import os,sys +from optparse import OptionParser +import damask + +scriptName = os.path.splitext(os.path.basename(__file__))[0] +scriptID = ' '.join([scriptName,damask.version]) + +#-------------------------------------------------------------------------------------------------- +# MAIN +#-------------------------------------------------------------------------------------------------- + +parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ +Unpack geometry files containing ranges "a to b" and/or "n of x" multiples (exclusively in one line). + +""", version = scriptID) + +parser.add_option('-l', '--labels', + dest = 'labels', + help = 'output geom file with one-dimensional data arrangement') + +parser.set_defaults(labels = ['1_eulerangles','1_eulerangles','1_eulerangles', + '1_pos','2_pos', 'IQ', 'CI', 'Fit', 'GrainID',], + ) + +(options, filenames) = parser.parse_args() + +# --- loop over input files ------------------------------------------------------------------------- + +if filenames == []: filenames = [None] + +for name in filenames: + try: + table = damask.ASCIItable(name = name, + buffered = False, + labeled = False) + except: continue + damask.util.report(scriptName,name) + table.head_read() + data = [] + while table.data_read(): + data.append(table.data[0:len(options.labels)]) + + table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) + table.labels_append(options.labels) + table.head_write() + for i in data: + table.data = i + table.data_write() From 7423c1a06acefc7642784a38e60d2e6db0327281 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 24 Mar 2016 18:47:11 +0100 Subject: [PATCH 05/18] not needed anymore --- processing/post/addCauchyX.py | 61 ---------- processing/pre/geom_fromAng.py | 206 --------------------------------- 2 files changed, 267 deletions(-) delete mode 100755 processing/post/addCauchyX.py delete mode 100755 processing/pre/geom_fromAng.py diff --git a/processing/post/addCauchyX.py b/processing/post/addCauchyX.py deleted file mode 100755 index dc94e2a2e..000000000 --- a/processing/post/addCauchyX.py +++ /dev/null @@ -1,61 +0,0 @@ -#!/usr/bin/env python -# -*- coding: UTF-8 no BOM -*- - -import os,string,h5py -import numpy as np -from optparse import OptionParser -import damask - -# -------------------------------------------------------------------- -# MAIN -# -------------------------------------------------------------------- - -parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ -Add column(s) containing Cauchy stress based on given column(s) of -deformation gradient and first Piola--Kirchhoff stress. - -""" + string.replace('$Id$','\n','\\n') -) - - -parser.add_option('-f','--defgrad', dest='defgrad', \ - help='heading of columns containing deformation gradient [%default]') -parser.add_option('-p','--stress', dest='stress', \ - help='heading of columns containing first Piola--Kirchhoff stress [%default]') -parser.add_option('-o','--output', dest='output', \ - help='group containing requested data [%default]') -parser.set_defaults(defgrad = 'f') -parser.set_defaults(stress = 'p') -parser.set_defaults(output = 'crystallite') - -(options,filenames) = parser.parse_args() - -if options.defgrad is None or options.stress is None or options.output is None: - parser.error('missing data column...') - - -# ------------------------------------------ setup file handles --------------------------------------- - -files = [] -for name in filenames: - if os.path.exists(name): - files.append({'name':name, 'file':h5py.File(name,"a")}) - -# ------------------------------------------ loop over input files ------------------------------------ - -for myFile in files: - print(myFile['name']) - -# ------------------------------------------ loop over increments ------------------------------------- - for inc in myFile['file']['increments'].keys(): - print("Current Increment: "+inc) - for instance in myFile['file']['increments/'+inc+'/'+options.output].keys(): - dsets = myFile['file']['increments/'+inc+'/'+options.output+'/'+instance].keys() - if (options.defgrad in dsets and options.stress in dsets): - defgrad = myFile['file']['increments/'+inc+'/'+options.output+'/'+instance+'/'+options.defgrad] - stress = myFile['file']['increments/'+inc+'/'+options.output+'/'+instance+'/'+options.stress] - cauchy=np.zeros(np.shape(stress),'f') - for p in range(stress.shape[0]): - cauchy[p,...] = 1.0/np.linalg.det(defgrad[p,...])*np.dot(stress[p,...],defgrad[p,...].T) # [Cauchy] = (1/det(F)) * [P].[F_transpose] - cauchyFile = myFile['file']['increments/'+inc+'/'+options.output+'/'+instance].create_dataset('cauchy', data=cauchy) - cauchyFile.attrs['units'] = 'Pa' diff --git a/processing/pre/geom_fromAng.py b/processing/pre/geom_fromAng.py deleted file mode 100755 index 987ff8ac1..000000000 --- a/processing/pre/geom_fromAng.py +++ /dev/null @@ -1,206 +0,0 @@ -#!/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]) - -#-------------------------------------------------------------------------------------------------- -# MAIN -#-------------------------------------------------------------------------------------------------- -parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ - -Generate geometry description and material configuration from EBSD data in given square-gridded 'ang' file. -Two phases can be discriminated based on threshold value in a given data column. - -""", version = scriptID) - -parser.add_option('--column', dest='column', type='int', metavar = 'int', - help='data column to discriminate between both phases [%default]') -parser.add_option('-t','--threshold', dest='threshold', type='float', metavar = 'float', - help='threshold value for phase discrimination [%default]') -parser.add_option('--homogenization', dest='homogenization', type='int', metavar = 'int', - help='homogenization index for configuration [%default]') -parser.add_option('--phase', dest='phase', type='int', nargs = 2, metavar = 'int int', - help='phase indices for configuration %default') -parser.add_option('--crystallite', dest='crystallite', type='int', metavar = 'int', - help='crystallite index for configuration [%default]') -parser.add_option('-c','--compress', dest='compress', action='store_true', - help='lump identical microstructure and texture information [%default]') -parser.add_option('-a', '--axes', dest='axes', nargs = 3, metavar = 'string string string', - help='Euler angle coordinate system for configuration x,y,z = %default') -parser.add_option('-p', '--precision', dest='precision', choices=['0','1','2','3'], metavar = 'int', - help = 'euler angles decimal places for output format and compressing (0,1,2,3) [2]') - -parser.set_defaults(column = 11, - threshold = 0.5, - homogenization = 1, - phase = [1,2], - crystallite = 1, - compress = False, - axes = ['y','x','-z'], - precision = '2', - ) -(options,filenames) = parser.parse_args() - - -for i in options.axes: - if i.lower() not in ['x','+x','-x','y','+y','-y','z','+z','-z']: - parser.error('invalid axes %s %s %s' %(options.axes[0],options.axes[1],options.axes[2])) - -# --- loop over input files ------------------------------------------------------------------------- - -if filenames == []: filenames = [None] - -for name in filenames: - try: - table = damask.ASCIItable(name = name, - outname = os.path.splitext(name)[-2]+'.geom' if name else name, - buffered = False, labeled = False) - except: continue - damask.util.report(scriptName,name) - - info = { - 'grid': np.ones(3,'i'), - 'size': np.zeros(3,'d'), - 'origin': np.zeros(3,'d'), - 'microstructures': 0, - 'homogenization': options.homogenization - } - - coords = [{},{},{1:True}] - pos = {'min':[ float("inf"), float("inf")], - 'max':[-float("inf"),-float("inf")]} - - phase = [] - eulerangles = [] - -# --------------- read data ----------------------------------------------------------------------- - errors = [] - while table.data_read(): - words = table.data - if words[0] == '#': # process initial comments/header block - if len(words) > 2: - if words[2].lower() == 'hexgrid': - errors.append('The file has HexGrid format. Please first convert to SquareGrid...') - break - else: - currPos = words[3:5] - for i in xrange(2): - coords[i][currPos[i]] = True - currPos = map(float,currPos) - for i in xrange(2): - pos['min'][i] = min(pos['min'][i],currPos[i]) - pos['max'][i] = max(pos['max'][i],currPos[i]) - eulerangles.append(map(math.degrees,map(float,words[:3]))) - phase.append(options.phase[int(float(words[options.column-1]) > options.threshold)]) - - if errors != []: - damask.util.croak(errors) - continue - -# --------------- determine size and grid --------------------------------------------------------- - info['grid'] = np.array(map(len,coords),'i') - info['size'][0:2] = info['grid'][0:2]/(info['grid'][0:2]-1.0)* \ - np.array([pos['max'][0]-pos['min'][0], - pos['max'][1]-pos['min'][1]],'d') - info['size'][2]=info['size'][0]/info['grid'][0] - eulerangles = np.array(eulerangles,dtype='f').reshape(info['grid'].prod(),3) - phase = np.array(phase,dtype='i').reshape(info['grid'].prod()) - - limits = [360,180,360] - if any([np.any(eulerangles[:,i]>=limits[i]) for i in [0,1,2]]): - errors.append('Error: euler angles out of bound. Ang file might contain unidexed poins.') - for i,angle in enumerate(['phi1','PHI','phi2']): - for n in np.nditer(np.where(eulerangles[:,i]>=limits[i]),['zerosize_ok']): - errors.append('%s in line %i (%4.2f %4.2f %4.2f)\n' - %(angle,n,eulerangles[n,0],eulerangles[n,1],eulerangles[n,2])) - if errors != []: - damask.util.croak(errors) - continue - - eulerangles=np.around(eulerangles,int(options.precision)) # round to desired precision -# ensure, that rounded euler angles are not out of bounds (modulo by limits) - for i,angle in enumerate(['phi1','PHI','phi2']): - eulerangles[:,i]%=limits[i] - -# scale angles by desired precision and convert to int. create unique integer key from three euler angles by -# concatenating the string representation with leading zeros and store as integer and search unique euler angle keys. -# Texture IDs are the indices of the first occurrence, the inverse is used to construct the microstructure -# create a microstructure (texture/phase pair) for each point using unique texture IDs. -# Use longInt (64bit, i8) because the keys might be long - if options.compress: - formatString='{0:0>'+str(int(options.precision)+3)+'}' - euleranglesRadInt = (eulerangles*10**int(options.precision)).astype('int') - eulerKeys = np.array([int(''.join(map(formatString.format,euleranglesRadInt[i,:]))) \ - for i in xrange(info['grid'].prod())]) - devNull, texture, eulerKeys_idx = np.unique(eulerKeys, return_index = True, return_inverse=True) - msFull = np.array([[eulerKeys_idx[i],phase[i]] for i in xrange(info['grid'].prod())],'i8') - devNull,msUnique,matPoints = np.unique(msFull.view('c16'),True,True) - matPoints+=1 - microstructure = np.array([msFull[i] for i in msUnique]) # pick only unique microstructures - else: - texture = np.arange(info['grid'].prod()) - microstructure = np.hstack( zip(texture,phase) ).reshape(info['grid'].prod(),2) # create texture/phase pairs - formatOut = 1+int(math.log10(len(texture))) - textureOut =[''] - - eulerFormatOut='%%%i.%if'%(int(options.precision)+4,int(options.precision)) - outStringAngles='(gauss) phi1 '+eulerFormatOut+' Phi '+eulerFormatOut+' phi2 '+eulerFormatOut+' scatter 0.0 fraction 1.0' - for i in xrange(len(texture)): - textureOut += ['[Texture%s]'%str(i+1).zfill(formatOut), - 'axes %s %s %s'%(options.axes[0],options.axes[1],options.axes[2]), - outStringAngles%tuple(eulerangles[texture[i],...]) - ] - formatOut = 1+int(math.log10(len(microstructure))) - microstructureOut =[''] - for i in xrange(len(microstructure)): - microstructureOut += ['[Grain%s]'%str(i+1).zfill(formatOut), - 'crystallite\t%i'%options.crystallite, - '(constituent)\tphase %i\ttexture %i\tfraction 1.0'%(microstructure[i,1],microstructure[i,0]+1) - ] - - info['microstructures'] = len(microstructure) - -#--- report --------------------------------------------------------------------------------------- - damask.util.croak('grid a b c: %s\n'%(' x '.join(map(str,info['grid']))) + - 'size x y z: %s\n'%(' x '.join(map(str,info['size']))) + - 'origin x y z: %s\n'%(' : '.join(map(str,info['origin']))) + - 'homogenization: %i\n'%info['homogenization'] + - 'microstructures: %i\n\n'%info['microstructures']) - - if np.any(info['grid'] < 1): - damask.util.croak('invalid grid a b c.\n') - continue - if np.any(info['size'] <= 0.0): - damask.util.croak('invalid size x y z.\n') - continue - - -#--- write data/header -------------------------------------------------------------------------------- - table.info_clear() - table.info_append([' '.join([scriptID] + sys.argv[1:]), - "grid\ta %i\tb %i\tc %i"%(info['grid'][0],info['grid'][1],info['grid'][2],), - "size\tx %f\ty %f\tz %f"%(info['size'][0],info['size'][1],info['size'][2],), - "origin\tx %f\ty %f\tz %f"%(info['origin'][0],info['origin'][1],info['origin'][2],), - "microstructures\t%i"%info['microstructures'], - "homogenization\t%i"%info['homogenization'], - ] + - [line for line in microstructureOut + textureOut] - ) - table.head_write() - if options.compress: - matPoints = matPoints.reshape((info['grid'][1],info['grid'][0])) - table.data = matPoints - table.data_writeArray('%%%ii'%(1+int(math.log10(np.amax(matPoints)))),delimiter=' ') - else: - table.output_write("1 to %i\n"%(info['microstructures'])) - table.output_flush() - -# --- output finalization -------------------------------------------------------------------------- - - table.close() From 86e9615744fd8b18619afae5461b538ff50bd528 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 25 Mar 2016 11:03:56 +0100 Subject: [PATCH 06/18] avoid error prone definition of derived quantities --- code/lattice.f90 | 71 +++++++++++++++++++++++------------------------- 1 file changed, 34 insertions(+), 37 deletions(-) diff --git a/code/lattice.f90 b/code/lattice.f90 index 05a123125..c652f3941 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -94,11 +94,11 @@ module lattice LATTICE_fcc_NcleavageSystem = int([3, 4, 0],pInt) !< total # of cleavage systems per family for fcc integer(pInt), parameter, private :: & - LATTICE_fcc_Nslip = 12_pInt, & ! sum(lattice_fcc_NslipSystem), & !< total # of slip systems for fcc - LATTICE_fcc_Ntwin = 12_pInt, & ! sum(lattice_fcc_NtwinSystem) !< total # of twin systems for fcc + LATTICE_fcc_Nslip = sum(lattice_fcc_NslipSystem), & !< total # of slip systems for fcc + LATTICE_fcc_Ntwin = sum(lattice_fcc_NtwinSystem), & !< total # of twin systems for fcc LATTICE_fcc_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for fcc - LATTICE_fcc_Ntrans = 12_pInt, & !< total # of transformations for fcc - LATTICE_fcc_Ncleavage = 7_pInt !< total # of cleavage systems for fcc + LATTICE_fcc_Ntrans = sum(lattice_fcc_NtransSystem), & !< total # of transformation systems for fcc + LATTICE_fcc_Ncleavage = sum(lattice_fcc_NcleavageSystem) !< total # of cleavage systems for fcc real(pReal), dimension(3+3,LATTICE_fcc_Nslip), parameter, private :: & LATTICE_fcc_systemSlip = reshape(real([& @@ -365,7 +365,7 @@ module lattice !-------------------------------------------------------------------------------------------------- ! bcc integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: & - LATTICE_bcc_NslipSystem = int([ 12, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], pInt) !< total # of slip systems per family for bcc + LATTICE_bcc_NslipSystem = int([ 12, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], pInt) !< total # of slip systems per family for bcc integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: & LATTICE_bcc_NtwinSystem = int([ 12, 0, 0, 0], pInt) !< total # of twin systems per family for bcc @@ -374,16 +374,15 @@ module lattice LATTICE_bcc_NtransSystem = int([0,0],pInt) !< total # of transformation systems per family for bcc integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & - LATTICE_bcc_NcleavageSystem = int([3,6,0],pInt) !< total # of cleavage systems per family for bcc + LATTICE_bcc_NcleavageSystem = int([3,6,0],pInt) !< total # of cleavage systems per family for bcc integer(pInt), parameter, private :: & - LATTICE_bcc_Nslip = 24_pInt, & ! sum(lattice_bcc_NslipSystem), & !< total # of slip systems for bcc - LATTICE_bcc_Ntwin = 12_pInt, & ! sum(lattice_bcc_NtwinSystem) !< total # of twin systems for bcc - LATTICE_bcc_NnonSchmid = 6_pInt, & !< # of non-Schmid contributions for bcc. 6 known non schmid contributions for BCC (A. Koester, A. Ma, A. Hartmaier 2012) - LATTICE_bcc_Ntrans = 0_pInt, & !< total # of transformations for bcc - LATTICE_bcc_Ncleavage = 9_pInt !< total # of cleavage systems for bcc + LATTICE_bcc_Nslip = sum(lattice_bcc_NslipSystem), & !< total # of slip systems for bcc + LATTICE_bcc_Ntwin = sum(lattice_bcc_NtwinSystem), & !< total # of twin systems for bcc + LATTICE_bcc_NnonSchmid = 6_pInt, & !< total # of non-Schmid contributions for bcc (A. Koester, A. Ma, A. Hartmaier 2012) + LATTICE_bcc_Ntrans = sum(lattice_bcc_NtransSystem), & !< total # of transformation systems for bcc + LATTICE_bcc_Ncleavage = sum(lattice_bcc_NcleavageSystem) !< total # of cleavage systems for bcc - real(pReal), dimension(3+3,LATTICE_bcc_Nslip), parameter, private :: & LATTICE_bcc_systemSlip = reshape(real([& ! Slip direction Plane normal @@ -563,7 +562,7 @@ module lattice !-------------------------------------------------------------------------------------------------- ! hex integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: & - lattice_hex_NslipSystem = int([ 3, 3, 3, 6, 12, 6, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for hex + lattice_hex_NslipSystem = int([ 3, 3, 3, 6, 12, 6, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for hex integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: & lattice_hex_NtwinSystem = int([ 6, 6, 6, 6],pInt) !< # of slip systems per family for hex @@ -572,14 +571,14 @@ module lattice LATTICE_hex_NtransSystem = int([0,0],pInt) !< total # of transformation systems per family for hex integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & - LATTICE_hex_NcleavageSystem = int([3,0,0],pInt) !< total # of cleavage systems per family for hex + LATTICE_hex_NcleavageSystem = int([3,0,0],pInt) !< total # of cleavage systems per family for hex - integer(pInt), parameter , private :: & - LATTICE_hex_Nslip = 33_pInt, & ! sum(lattice_hex_NslipSystem), !< total # of slip systems for hex - LATTICE_hex_Ntwin = 24_pInt, & ! sum(lattice_hex_NtwinSystem) !< total # of twin systems for hex - LATTICE_hex_NnonSchmid = 0_pInt, & !< # of non-Schmid contributions for hex - LATTICE_hex_Ntrans = 0_pInt, & !< total # of transformations for hex - LATTICE_hex_Ncleavage = 3_pInt !< total # of transformations for hex + integer(pInt), parameter, private :: & + LATTICE_hex_Nslip = sum(lattice_hex_NslipSystem), & !< total # of slip systems for hex + LATTICE_hex_Ntwin = sum(lattice_hex_NtwinSystem), & !< total # of twin systems for hex + LATTICE_hex_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for hex + LATTICE_hex_Ntrans = sum(lattice_hex_NtransSystem), & !< total # of transformation systems for hex + LATTICE_hex_Ncleavage = sum(lattice_hex_NcleavageSystem) !< total # of cleavage systems for hex real(pReal), dimension(4+4,LATTICE_hex_Nslip), parameter, private :: & LATTICE_hex_systemSlip = reshape(real([& @@ -842,7 +841,6 @@ module lattice ],pReal),[ 4_pInt + 4_pInt,LATTICE_hex_Ncleavage]) - !-------------------------------------------------------------------------------------------------- ! bct integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: & @@ -856,14 +854,13 @@ module lattice integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & LATTICE_bct_NcleavageSystem = int([0,0,0],pInt) !< total # of cleavage systems per family for bct - - integer(pInt), parameter , private :: & - LATTICE_bct_Nslip = 52_pInt, & ! sum(lattice_bct_NslipSystem), !< total # of slip systems for bct - LATTICE_bct_Ntwin = 0_pInt, & ! sum(lattice_bcc_NtwinSystem) !< total # of twin systems for bct - LATTICE_bct_NnonSchmid = 0_pInt, & !< # of non-Schmid contributions for bct - LATTICE_bct_Ntrans = 0_pInt, & !< total # of transformations for bct - LATTICE_bct_Ncleavage = 0_pInt !< total # of transformations for bct + integer(pInt), parameter, private :: & + LATTICE_bct_Nslip = sum(lattice_bct_NslipSystem), & !< total # of slip systems for bct + LATTICE_bct_Ntwin = sum(lattice_bct_NtwinSystem), & !< total # of twin systems for bct + LATTICE_bct_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for bct + LATTICE_bct_Ntrans = sum(lattice_bct_NtransSystem), & !< total # of transformation systems for bct + LATTICE_bct_Ncleavage = sum(lattice_bct_NcleavageSystem) !< total # of cleavage systems for bct real(pReal), dimension(3+3,LATTICE_bct_Nslip), parameter, private :: & LATTICE_bct_systemSlip = reshape(real([& @@ -1007,10 +1004,10 @@ module lattice !-------------------------------------------------------------------------------------------------- ! isotropic integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & - LATTICE_iso_NcleavageSystem = int([3,0,0],pInt) !< total # of cleavage systems per family for isotropic + LATTICE_iso_NcleavageSystem = int([3,0,0],pInt) !< total # of cleavage systems per family for iso integer(pInt), parameter, private :: & - LATTICE_iso_Ncleavage = 3_pInt !< total # of cleavage systems for bcc + LATTICE_iso_Ncleavage = sum(LATTICE_iso_NcleavageSystem) !< total # of cleavage systems for iso real(pReal), dimension(3+3,LATTICE_iso_Ncleavage), parameter, private :: & LATTICE_iso_systemCleavage = reshape(real([& @@ -1023,10 +1020,10 @@ module lattice !-------------------------------------------------------------------------------------------------- ! orthorhombic integer(pInt), dimension(LATTICE_maxNcleavageFamily), parameter, public :: & - LATTICE_ortho_NcleavageSystem = int([1,1,1],pInt) !< total # of cleavage systems per family for orthotropic + LATTICE_ortho_NcleavageSystem = int([1,1,1],pInt) !< total # of cleavage systems per family for ortho integer(pInt), parameter, private :: & - LATTICE_ortho_Ncleavage = 3_pInt !< total # of cleavage systems for bcc + LATTICE_ortho_Ncleavage = sum(LATTICE_ortho_NcleavageSystem) !< total # of cleavage systems for ortho real(pReal), dimension(3+3,LATTICE_ortho_Ncleavage), parameter, private :: & LATTICE_ortho_systemCleavage = reshape(real([& @@ -1036,16 +1033,16 @@ module lattice 1, 0, 0, 0, 0, 1 & ],pReal),[ 3_pInt + 3_pInt,LATTICE_ortho_Ncleavage]) - real(pReal), dimension(:,:,:), allocatable, public, protected :: & + real(pReal), dimension(:,:,:), allocatable, public, protected :: & lattice_C66, lattice_trans_C66 - real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & + real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & lattice_C3333, lattice_trans_C3333 - real(pReal), dimension(:), allocatable, public, protected :: & + real(pReal), dimension(:), allocatable, public, protected :: & lattice_mu, & lattice_nu, & lattice_trans_mu, & lattice_trans_nu - real(pReal), dimension(:,:,:), allocatable, public, protected :: & + real(pReal), dimension(:,:,:), allocatable, public, protected :: & lattice_thermalConductivity33, & lattice_thermalExpansion33, & lattice_damageDiffusion33, & @@ -1054,7 +1051,7 @@ module lattice lattice_porosityDiffusion33, & lattice_hydrogenfluxDiffusion33, & lattice_hydrogenfluxMobility33 - real(pReal), dimension(:), allocatable, public, protected :: & + real(pReal), dimension(:), allocatable, public, protected :: & lattice_damageMobility, & lattice_porosityMobility, & lattice_massDensity, & From 97314619cc2cbba361461709562bdc2b384cb939 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 26 Mar 2016 19:59:05 +0100 Subject: [PATCH 07/18] corrections --- processing/pre/table_fromOIMgrainFile.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/processing/pre/table_fromOIMgrainFile.py b/processing/pre/table_fromOIMgrainFile.py index c730c8892..d6ac020bb 100755 --- a/processing/pre/table_fromOIMgrainFile.py +++ b/processing/pre/table_fromOIMgrainFile.py @@ -13,15 +13,15 @@ scriptID = ' '.join([scriptName,damask.version]) #-------------------------------------------------------------------------------------------------- parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ -Unpack geometry files containing ranges "a to b" and/or "n of x" multiples (exclusively in one line). +Adds header to OIM grain file to make it accesible as ASCII table """, version = scriptID) parser.add_option('-l', '--labels', dest = 'labels', - help = 'output geom file with one-dimensional data arrangement') + help = 'lables for requested columns') -parser.set_defaults(labels = ['1_eulerangles','1_eulerangles','1_eulerangles', +parser.set_defaults(labels = ['1_eulerangles','2_eulerangles','3_eulerangles', '1_pos','2_pos', 'IQ', 'CI', 'Fit', 'GrainID',], ) @@ -49,3 +49,7 @@ for name in filenames: for i in data: table.data = i table.data_write() + +# --- output finalization -------------------------------------------------------------------------- + + table.close() # close ASCII table From 752d2ddcc033a31f4bd867e9ed714af40c54bf47 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 10 Apr 2016 15:41:15 +0200 Subject: [PATCH 08/18] not really working and not used at all --- DAMASK_env.bat | 19 ------------------- 1 file changed, 19 deletions(-) delete mode 100644 DAMASK_env.bat diff --git a/DAMASK_env.bat b/DAMASK_env.bat deleted file mode 100644 index 202ddb00e..000000000 --- a/DAMASK_env.bat +++ /dev/null @@ -1,19 +0,0 @@ -:: sets up an environment for DAMASK on Windows -:: usage: call DAMASK_env.bat -@echo off -set LOCATION=%~dp0 -set DAMASK_ROOT=%LOCATION%\DAMASK -set DAMASK_NUM_THREADS=2 -chcp 1252 -Title Düsseldorf Advanced Materials Simulation Kit - DAMASK, MPIE Düsseldorf -echo. -echo Düsseldorf Advanced Materials Simulation Kit - DAMASK -echo Max-Planck-Institut für Eisenforschung, Düsseldorf -echo http://damask.mpie.de -echo. -echo Preparing environment ... -echo DAMASK_ROOT=%DAMASK_ROOT% -echo DAMASK_NUM_THREADS=%DAMASK_NUM_THREADS% -set DAMASK_BIN=%DAMASK_ROOT%\bin -set PATH=%PATH%;%DAMASK_BIN% -set PYTHONPATH=%PYTHONPATH%;%DAMASK_ROOT%\lib From 82063494fd8792c3d0431b347a85f8572507cb1c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 11 Apr 2016 19:47:07 +0200 Subject: [PATCH 09/18] adopted to new json format of paraview introduced "standard" colormaps for stress, strain, and orientation file is now implicit the first argument without key (in line with other scripts) --- lib/damask/colormaps.py | 40 ++++++++++++-------- processing/post/perceptualUniformColorMap.py | 16 ++++---- 2 files changed, 32 insertions(+), 24 deletions(-) diff --git a/lib/damask/colormaps.py b/lib/damask/colormaps.py index 314581471..45c88e8d5 100644 --- a/lib/damask/colormaps.py +++ b/lib/damask/colormaps.py @@ -303,36 +303,45 @@ class Colormap(): 'interpolate', ] __predefined__ = { - 'gray': {'left': Color('HSL',[0,1,1]), + 'gray': {'left': Color('HSL',[0,1,1]), 'right': Color('HSL',[0,0,0.15]), 'interpolate': 'perceptualuniform'}, - 'grey': {'left': Color('HSL',[0,1,1]), + 'grey': {'left': Color('HSL',[0,1,1]), 'right': Color('HSL',[0,0,0.15]), 'interpolate': 'perceptualuniform'}, - 'red': {'left': Color('HSL',[0,1,0.14]), + 'red': {'left': Color('HSL',[0,1,0.14]), 'right': Color('HSL',[0,0.35,0.91]), 'interpolate': 'perceptualuniform'}, - 'green': {'left': Color('HSL',[0.33333,1,0.14]), + 'green': {'left': Color('HSL',[0.33333,1,0.14]), 'right': Color('HSL',[0.33333,0.35,0.91]), 'interpolate': 'perceptualuniform'}, - 'blue': {'left': Color('HSL',[0.66,1,0.14]), + 'blue': {'left': Color('HSL',[0.66,1,0.14]), 'right': Color('HSL',[0.66,0.35,0.91]), 'interpolate': 'perceptualuniform'}, - 'seaweed': {'left': Color('HSL',[0.78,1.0,0.1]), + 'seaweed': {'left': Color('HSL',[0.78,1.0,0.1]), 'right': Color('HSL',[0.40000,0.1,0.9]), 'interpolate': 'perceptualuniform'}, - 'bluebrown': {'left': Color('HSL',[0.65,0.53,0.49]), + 'bluebrown': {'left': Color('HSL',[0.65,0.53,0.49]), 'right': Color('HSL',[0.11,0.75,0.38]), 'interpolate': 'perceptualuniform'}, - 'redgreen': {'left': Color('HSL',[0.97,0.96,0.36]), + 'redgreen': {'left': Color('HSL',[0.97,0.96,0.36]), 'right': Color('HSL',[0.33333,1.0,0.14]), 'interpolate': 'perceptualuniform'}, - 'bluered': {'left': Color('HSL',[0.65,0.53,0.49]), + 'bluered': {'left': Color('HSL',[0.65,0.53,0.49]), 'right': Color('HSL',[0.97,0.96,0.36]), 'interpolate': 'perceptualuniform'}, - 'blueredrainbow':{'left': Color('HSL',[2.0/3.0,1,0.5]), + 'blueredrainbow':{'left': Color('HSL',[2.0/3.0,1,0.5]), 'right': Color('HSL',[0,1,0.5]), 'interpolate': 'linear' }, + 'orientation': {'left': Color('RGB',[0.933334,0.878432,0.878431]), + 'right': Color('RGB',[0.250980,0.007843,0.000000]), + 'interpolate': 'perceptualuniform'}, + 'strain': {'left': Color('RGB',[0.941177,0.941177,0.870588]), + 'right': Color('RGB',[0.266667,0.266667,0.000000]), + 'interpolate': 'perceptualuniform'}, + 'stress': {'left': Color('RGB',[0.878432,0.874511,0.949019]), + 'right': Color('RGB',[0.000002,0.000000,0.286275]), + 'interpolate': 'perceptualuniform'}, } @@ -344,7 +353,7 @@ class Colormap(): predefined = None ): - if str(predefined).lower() in self.__predefined__: + if predefined is not None: left = self.__predefined__[predefined.lower()]['left'] right= self.__predefined__[predefined.lower()]['right'] interpolate = self.__predefined__[predefined.lower()]['interpolate'] @@ -442,11 +451,12 @@ class Colormap(): format = format.lower() # consistent comparison basis frac = 0.5*(np.array(crop) + 1.0) # rescale crop range to fractions colors = [self.color(float(i)/(steps-1)*(frac[1]-frac[0])+frac[0]).expressAs(model).color for i in xrange(steps)] - if format == 'paraview': - colormap = [''] \ - + [''%(color[0],color[1],color[2],) for i,color in colors] \ - + [''] + colormap = ['[\n {{\n "ColorSpace" : "RGB", "Name" : "{}",\n "RGBPoints" : ['.format(name)] \ + + [' {:4d},{:8.6f},{:8.6f},{:8.6f},'.format(i,color[0],color[1],color[2],) + for i,color in enumerate(colors[:-1])]\ + + [' {:4d},{:8.6f},{:8.6f},{:8.6f} '.format(i+1,colors[-1][0],colors[-1][1],colors[-1][2],)]\ + + [' ]\n }\n]'] elif format == 'gmsh': colormap = ['View.ColorTable = {'] \ diff --git a/processing/post/perceptualUniformColorMap.py b/processing/post/perceptualUniformColorMap.py index c2201f76b..bef6a2187 100755 --- a/processing/post/perceptualUniformColorMap.py +++ b/processing/post/perceptualUniformColorMap.py @@ -14,7 +14,7 @@ scriptID = ' '.join([scriptName,damask.version]) #Borland, D., & Taylor, R. M. (2007). Rainbow Color Map (Still) Considered Harmful. Computer Graphics and Applications, IEEE, 27(2), 14--17. #Moreland, K. (2009). Diverging Color Maps for Scientific Visualization. In Proc. 5th Int. Symp. Visual Computing (pp. 92--103). outtypes = ['paraview','gmsh','raw','GOM'] -extensions = ['.xml','.msh','.txt','.legend'] +extensions = ['.json','.msh','.txt','.legend'] colormodels = ['RGB','HSL','XYZ','CIELAB','MSH'] parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ @@ -34,11 +34,9 @@ parser.add_option('-r','--right', dest='right', type='float', nargs=3, metavar=' parser.add_option('-c','--colormodel', dest='colormodel', metavar='string', help='colormodel: '+', '.join(colormodels)+' [%default]') parser.add_option('-p','--predefined', dest='predefined', metavar='string', - help='predefined colormap [%default]') + help='predefined colormap') parser.add_option('-f','--format', dest='format', metavar='string', help='output format: '+', '.join(outtypes)+' [%default]') -parser.add_option('-b','--basename', dest='basename', metavar='string', - help='basename of output file [%default]') parser.set_defaults(colormodel = 'RGB') parser.set_defaults(predefined = None) parser.set_defaults(basename = None) @@ -48,7 +46,7 @@ parser.set_defaults(trim = (-1.0,1.0)) parser.set_defaults(left = (1.0,1.0,1.0)) parser.set_defaults(right = (0.0,0.0,0.0)) -(options,filenames) = parser.parse_args() +(options,filename) = parser.parse_args() if options.format not in outtypes: parser.error('invalid format: "%s" (can be %s).'%(options.format,', '.join(outtypes))) @@ -62,10 +60,10 @@ if options.trim[0] < -1.0 or \ parser.error('invalid trim range (-1 +1).') -name = options.format if options.basename is None\ - else options.basename -output = sys.stdout if options.basename is None\ - else open(os.path.basename(options.basename)+extensions[outtypes.index(options.format)],'w') +name = options.format if filename[0] is None\ + else filename[0] +output = sys.stdout if filename[0] is None\ + else open(os.path.basename(filename[0])+extensions[outtypes.index(options.format)],'w') colorLeft = damask.Color(options.colormodel.upper(), list(options.left)) colorRight = damask.Color(options.colormodel.upper(), list(options.right)) From 18f18aa4b93f31ca8d422277ec0f0ea149e1fb18 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 11 Apr 2016 19:55:55 +0200 Subject: [PATCH 10/18] detabbing --- code/lattice.f90 | 4 ++-- processing/pre/table_fromOIMgrainFile.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/code/lattice.f90 b/code/lattice.f90 index c652f3941..d08d42c06 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -576,7 +576,7 @@ module lattice integer(pInt), parameter, private :: & LATTICE_hex_Nslip = sum(lattice_hex_NslipSystem), & !< total # of slip systems for hex LATTICE_hex_Ntwin = sum(lattice_hex_NtwinSystem), & !< total # of twin systems for hex - LATTICE_hex_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for hex + LATTICE_hex_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for hex LATTICE_hex_Ntrans = sum(lattice_hex_NtransSystem), & !< total # of transformation systems for hex LATTICE_hex_Ncleavage = sum(lattice_hex_NcleavageSystem) !< total # of cleavage systems for hex @@ -858,7 +858,7 @@ module lattice integer(pInt), parameter, private :: & LATTICE_bct_Nslip = sum(lattice_bct_NslipSystem), & !< total # of slip systems for bct LATTICE_bct_Ntwin = sum(lattice_bct_NtwinSystem), & !< total # of twin systems for bct - LATTICE_bct_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for bct + LATTICE_bct_NnonSchmid = 0_pInt, & !< total # of non-Schmid contributions for bct LATTICE_bct_Ntrans = sum(lattice_bct_NtransSystem), & !< total # of transformation systems for bct LATTICE_bct_Ncleavage = sum(lattice_bct_NcleavageSystem) !< total # of cleavage systems for bct diff --git a/processing/pre/table_fromOIMgrainFile.py b/processing/pre/table_fromOIMgrainFile.py index d6ac020bb..9937af8c1 100755 --- a/processing/pre/table_fromOIMgrainFile.py +++ b/processing/pre/table_fromOIMgrainFile.py @@ -42,7 +42,7 @@ for name in filenames: data = [] while table.data_read(): data.append(table.data[0:len(options.labels)]) - + table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) table.labels_append(options.labels) table.head_write() From 70afa462b28fb23f6223a5061eef5c221b1b008b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 11 Apr 2016 20:25:24 +0200 Subject: [PATCH 11/18] should work now with odd resolution as well --- processing/post/addCurl.py | 7 ++++--- processing/post/addDivergence.py | 7 ++++--- processing/post/addGradient.py | 8 +++++--- 3 files changed, 13 insertions(+), 9 deletions(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index d8b1ee025..78377330c 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -10,6 +10,7 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) def curlFFT(geomdim,field): + shapeFFT = np.array(np.shape(field))[0:3] grid = np.array(np.shape(field)[2::-1]) N = grid.prod() # field size n = np.array(np.shape(field)[3:]).prod() # data size @@ -17,8 +18,8 @@ def curlFFT(geomdim,field): 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') + field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2),s=shapeFFT) + curl_fourier = np.empty(field_fourier.shape,'c16') # differentiation in Fourier space k_s = np.zeros([3],'i') @@ -55,7 +56,7 @@ def curlFFT(geomdim,field): curl_fourier[i,j,k,2] = ( field_fourier[i,j,k,1]*xi[0]\ -field_fourier[i,j,k,0]*xi[1]) *TWOPIIMG - return np.fft.fftpack.irfftn(curl_fourier,axes=(0,1,2)).reshape([N,n]) + return np.fft.fftpack.irfftn(curl_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n]) # -------------------------------------------------------------------- diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index aadaceabf..8037e6df5 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -10,12 +10,13 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) def divFFT(geomdim,field): + shapeFFT = np.array(np.shape(field))[0:3] 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 + field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2),s=shapeFFT) + div_fourier = np.empty(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') @@ -41,7 +42,7 @@ def divFFT(geomdim,field): 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]) + return np.fft.fftpack.irfftn(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n/3]) # -------------------------------------------------------------------- diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index 555e587de..44622d920 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -10,14 +10,16 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) def gradFFT(geomdim,field): + shapeFFT = np.array(np.shape(field))[0:3] 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') + field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2),s=shapeFFT) + grad_fourier = np.empty(field_fourier.shape+(3,),'c16') # differentiation in Fourier space k_s = np.zeros([3],'i') @@ -44,7 +46,7 @@ def gradFFT(geomdim,field): 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]) + return np.fft.fftpack.irfftn(grad_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,3*n]) # -------------------------------------------------------------------- From ce25acce77ae7f01b5b4865be7f207182dae1541 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 11 Apr 2016 21:00:43 +0200 Subject: [PATCH 12/18] no allocation for disorientation for local models only --- code/crystallite.f90 | 53 +++++++++++++++++++++++--------------------- 1 file changed, 28 insertions(+), 25 deletions(-) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 71a5e4743..0c15ea686 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -258,7 +258,8 @@ subroutine crystallite_init allocate(crystallite_orientation(4,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_orientation0(4,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_rotation(4,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_disorientation(4,nMax,cMax,iMax,eMax), source=0.0_pReal) + if (any(plasticState%nonLocal)) & + allocate(crystallite_disorientation(4,nMax,cMax,iMax,eMax),source=0.0_pReal) allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.) allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) allocate(crystallite_todo(cMax,iMax,eMax), source=.false.) @@ -3977,50 +3978,51 @@ subroutine crystallite_orientations ! --- CALCULATE ORIENTATION AND LATTICE ROTATION --- - !$OMP PARALLEL DO PRIVATE(orientation) + nonlocalPresent: if (any(plasticState%nonLocal)) then +!$OMP PARALLEL DO PRIVATE(orientation) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e)) ! somehow this subroutine is not threadsafe, so need critical statement here; not clear, what exactly the problem is - !$OMP CRITICAL (polarDecomp) +!$OMP CRITICAL (polarDecomp) orientation = math_RtoQ(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e)))) ! rotational part from polar decomposition as quaternion - !$OMP END CRITICAL (polarDecomp) +!$OMP END CRITICAL (polarDecomp) crystallite_rotation(1:4,c,i,e) = lattice_qDisorientation(crystallite_orientation0(1:4,c,i,e), & ! active rotation from ori0 - orientation) ! to current orientation (with no symmetry) + orientation) ! to current orientation (with no symmetry) crystallite_orientation(1:4,c,i,e) = orientation enddo; enddo; enddo - !$OMP END PARALLEL DO - - +!$OMP END PARALLEL DO + + ! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL --- ! --- we use crystallite_orientation from above, so need a separate loop - - !$OMP PARALLEL DO PRIVATE(myPhase,neighboring_e,neighboring_i,neighboringPhase) + +!$OMP PARALLEL DO PRIVATE(myPhase,neighboring_e,neighboring_i,neighboringPhase) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - myPhase = material_phase(1,i,e) ! get my phase (non-local models make no sense with more than one grain per material point) - if (plasticState(myPhase)%nonLocal) then ! if nonlocal model + myPhase = material_phase(1,i,e) ! get my phase (non-local models make no sense with more than one grain per material point) + if (plasticState(myPhase)%nonLocal) then ! if nonlocal model ! --- calculate disorientation between me and my neighbor --- - do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e)))) ! loop through my neighbors + do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e)))) ! loop through my neighbors neighboring_e = mesh_ipNeighborhood(1,n,i,e) neighboring_i = mesh_ipNeighborhood(2,n,i,e) - if (neighboring_e > 0 .and. neighboring_i > 0) then ! if neighbor exists - neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's phase - if (plasticState(neighboringPhase)%nonLocal) then ! neighbor got also nonlocal plasticity - if (lattice_structure(myPhase) == lattice_structure(neighboringPhase)) then ! if my neighbor has same crystal structure like me + if (neighboring_e > 0 .and. neighboring_i > 0) then ! if neighbor exists + neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's phase + if (plasticState(neighboringPhase)%nonLocal) then ! neighbor got also nonlocal plasticity + if (lattice_structure(myPhase) == lattice_structure(neighboringPhase)) then ! if my neighbor has same crystal structure like me crystallite_disorientation(:,n,1,i,e) = & lattice_qDisorientation( crystallite_orientation(1:4,1,i,e), & crystallite_orientation(1:4,1,neighboring_i,neighboring_e), & - lattice_structure(myPhase)) ! calculate disorientation for given symmetry - else ! for neighbor with different phase - crystallite_disorientation(:,n,1,i,e) = [0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal] ! 180 degree rotation about 100 axis + lattice_structure(myPhase)) ! calculate disorientation for given symmetry + else ! for neighbor with different phase + crystallite_disorientation(:,n,1,i,e) = [0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal]! 180 degree rotation about 100 axis endif - else ! for neighbor with local plasticity - crystallite_disorientation(:,n,1,i,e) = [-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal] ! homomorphic identity + else ! for neighbor with local plasticity + crystallite_disorientation(:,n,1,i,e) = [-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal]! homomorphic identity endif - else ! no existing neighbor - crystallite_disorientation(:,n,1,i,e) = [-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal] ! homomorphic identity + else ! no existing neighbor + crystallite_disorientation(:,n,1,i,e) = [-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal] ! homomorphic identity endif enddo @@ -4031,7 +4033,8 @@ subroutine crystallite_orientations endif enddo; enddo - !$OMP END PARALLEL DO +!$OMP END PARALLEL DO + endif nonlocalPresent end subroutine crystallite_orientations From 48e508b5add37e7dae91515664bc93198cb09ffe Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 12 Apr 2016 10:34:58 +0200 Subject: [PATCH 13/18] MPI_OFFSET_KIND is a long int --- code/DAMASK_spectral.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index 45451df08..371020fde 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -427,7 +427,7 @@ program DAMASK_spectral ! prepare MPI parallel out (including opening of file) allocate(outputSize(worldsize), source = 0_MPI_OFFSET_KIND) outputSize(worldrank+1) = size(materialpoint_results,kind=MPI_OFFSET_KIND)*int(pReal,MPI_OFFSET_KIND) - call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process + call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr)! get total output size over each process if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_allreduce') call MPI_file_open(PETSC_COMM_WORLD, & trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut', & From 43e74f00cf1f6637dabd65b779f6ed513f5aa80c Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 25 Apr 2016 12:14:05 +0200 Subject: [PATCH 14/18] updated version information after successful test of v2.0.0-160-g26f437b --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index c271e0939..833374b44 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-97-g8b27de7 +v2.0.0-160-g26f437b From 2eb7ad7432bbbb753ed8876f2d758129664e39d5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Apr 2016 12:57:38 +0200 Subject: [PATCH 15/18] polishing --- processing/post/addCurl.py | 16 ++++++++-------- processing/post/addDivergence.py | 12 ++++++------ processing/post/addGradient.py | 17 +++++++++-------- processing/pre/table_fromOIMgrainFile.py | 4 ++-- 4 files changed, 25 insertions(+), 24 deletions(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index 31a918316..f5db77234 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -63,25 +63,25 @@ def curlFFT(geomdim,field): # MAIN # -------------------------------------------------------------------- -parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ +parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """ Add column(s) containing curl of requested column(s). Operates on periodic ordered three-dimensional data sets. -Deals with both vector- and tensor-valued fields. +Deals with both vector- and tensor fields. """, version = scriptID) -parser.add_option('-c','--coordinates', +parser.add_option('-p','--pos','--periodiccellcenter' dest = 'coords', type = 'string', metavar = 'string', - help = 'column label of coordinates [%default]') + help = 'label of coordinates [%default]') parser.add_option('-v','--vector', dest = 'vector', action = 'extend', metavar = '', - help = 'column label(s) of vector field values') + help = 'label(s) of vector field values') parser.add_option('-t','--tensor', dest = 'tensor', action = 'extend', metavar = '', - help = 'column label(s) of tensor field values') + help = 'label(s) of tensor field values') parser.set_defaults(coords = 'pos', ) @@ -91,7 +91,7 @@ parser.set_defaults(coords = 'pos', 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] @@ -148,7 +148,7 @@ 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])) + 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 ----------------------------------- diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 935f7d83a..52ec2df75 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -19,7 +19,7 @@ def divFFT(geomdim,field): div_fourier = np.empty(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') + k_s = np.zeros([3],'i') TWOPIIMG = 2.0j*math.pi for i in xrange(grid[2]): k_s[0] = i @@ -49,25 +49,25 @@ def divFFT(geomdim,field): # MAIN # -------------------------------------------------------------------- -parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ +parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """ Add column(s) containing divergence of requested column(s). Operates on periodic ordered three-dimensional data sets. Deals with both vector- and tensor-valued fields. """, version = scriptID) -parser.add_option('-c','--coordinates', +parser.add_option('-p','--pos','--periodiccellcenter' dest = 'coords', type = 'string', metavar = 'string', - help = 'column label of coordinates [%default]') + help = 'label of coordinates [%default]') parser.add_option('-v','--vector', dest = 'vector', action = 'extend', metavar = '', - help = 'column label(s) of vector field values') + help = 'label(s) of vector field values') parser.add_option('-t','--tensor', dest = 'tensor', action = 'extend', metavar = '', - help = 'column label(s) of tensor field values') + help = 'label(s) of tensor field values') parser.set_defaults(coords = 'pos', ) diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index a741b3a70..7ef042044 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -14,6 +14,7 @@ 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' @@ -52,27 +53,27 @@ def gradFFT(geomdim,field): # MAIN # -------------------------------------------------------------------- -parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ +parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(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', +parser.add_option('-p','--pos','--periodiccellcenter' dest = 'coords', - type = 'string', metavar='string', - help = 'column heading for coordinates [%default]') + type = 'string', metavar = 'string', + help = 'label of coordinates [%default]') parser.add_option('-v','--vector', dest = 'vector', action = 'extend', metavar = '', - help = 'heading of columns containing vector field values') + help = 'label(s) of vector field values') parser.add_option('-s','--scalar', dest = 'scalar', action = 'extend', metavar = '', - help = 'heading of columns containing scalar field values') + help = 'label(s) of scalar field values') -parser.set_defaults(coords = 'ipinitialcoord', +parser.set_defaults(coords = 'pos', ) (options,filenames) = parser.parse_args() @@ -137,7 +138,7 @@ 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])) + 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 ----------------------------------- diff --git a/processing/pre/table_fromOIMgrainFile.py b/processing/pre/table_fromOIMgrainFile.py index 9937af8c1..d448ba896 100755 --- a/processing/pre/table_fromOIMgrainFile.py +++ b/processing/pre/table_fromOIMgrainFile.py @@ -19,9 +19,9 @@ Adds header to OIM grain file to make it accesible as ASCII table parser.add_option('-l', '--labels', dest = 'labels', - help = 'lables for requested columns') + help = 'lables of requested columns') -parser.set_defaults(labels = ['1_eulerangles','2_eulerangles','3_eulerangles', +parser.set_defaults(labels = ['1_euler','2_euler','3_euler', '1_pos','2_pos', 'IQ', 'CI', 'Fit', 'GrainID',], ) From 26e5f97ff3616e1e029fd83b05a010ea8c85150f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Apr 2016 13:22:34 +0200 Subject: [PATCH 16/18] missing comma --- processing/post/addCurl.py | 2 +- processing/post/addDivergence.py | 2 +- processing/post/addGradient.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index f5db77234..e368dc482 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -70,7 +70,7 @@ Deals with both vector- and tensor fields. """, version = scriptID) -parser.add_option('-p','--pos','--periodiccellcenter' +parser.add_option('-p','--pos','--periodiccellcenter', dest = 'coords', type = 'string', metavar = 'string', help = 'label of coordinates [%default]') diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 52ec2df75..d62c02fda 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -56,7 +56,7 @@ Deals with both vector- and tensor-valued fields. """, version = scriptID) -parser.add_option('-p','--pos','--periodiccellcenter' +parser.add_option('-p','--pos','--periodiccellcenter', dest = 'coords', type = 'string', metavar = 'string', help = 'label of coordinates [%default]') diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index 7ef042044..aaf1750a7 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -60,7 +60,7 @@ Deals with both vector- and scalar fields. """, version = scriptID) -parser.add_option('-p','--pos','--periodiccellcenter' +parser.add_option('-p','--pos','--periodiccellcenter', dest = 'coords', type = 'string', metavar = 'string', help = 'label of coordinates [%default]') From c30c7714ad7f134ba694f2c2b9ccae6a5641dc98 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Apr 2016 14:11:27 +0200 Subject: [PATCH 17/18] hickup, probably from merge --- processing/post/addGradient.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index aaf1750a7..b2d004bd2 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -18,7 +18,7 @@ def gradFFT(geomdim,field): if n == 3: dataType = 'vector' elif n == 1: dataType = 'scalar' - field_fourier = nps=shapeFFT).fft.fftpack.rfftn(field,axes=(0,1,2), + field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2),s=shapeFFT) grad_fourier = np.empty(field_fourier.shape+(3,),'c16') # differentiation in Fourier space From 7ed4ce859af63f7e01911ff96117cabd817ce9e6 Mon Sep 17 00:00:00 2001 From: zhangc43 Date: Mon, 25 Apr 2016 09:24:05 -0400 Subject: [PATCH 18/18] remove obsolete geom_frombarucentric script for microstructure reconstruction --- processing/pre/geom_fromBarycentric.py | 297 ------------------------- 1 file changed, 297 deletions(-) delete mode 100755 processing/pre/geom_fromBarycentric.py diff --git a/processing/pre/geom_fromBarycentric.py b/processing/pre/geom_fromBarycentric.py deleted file mode 100755 index fe687fd2a..000000000 --- a/processing/pre/geom_fromBarycentric.py +++ /dev/null @@ -1,297 +0,0 @@ -#!/usr/bin/env python - -## -# This script will read in all the seeds and partition the space -# using scipy.spatial.Delaunay triangulation. -# The unknown location will be then interpolated through Barycentric -# interpolation method, which relies on the triangulation. -# A rim will be automatically added to the patch, which will help -# improve the compatibility with the spectral solver as well as -# maintain meaningful microstructure(reduce artifacts). - - -import os -import numpy as np -import argparse -from scipy.spatial import Delaunay -import damask - -scriptName = os.path.splitext(os.path.basename(__file__))[0] -scriptID = ' '.join([scriptName,damask.version]) - - -OFFSET = 0.1 #resize the seeded volume to give space for rim/pan -PHANTOM_ID = -1 #grain ID for phantom seeds - -def d_print(info, data, separator=False): - """quickly print debug information""" - if(separator): print "*"*80 - print info - print data - - -def meshgrid2(*arrs): - """code inspired by http://stackoverflow.com/questions/1827489/numpy-meshgrid-in-3d""" - arrs = tuple(reversed(arrs)) - arrs = tuple(arrs) - lens = np.array(map(len, arrs)) - dim = len(arrs) - ans = [] - for i, arr in enumerate(arrs): - slc = np.ones(dim,'i') - slc[i] = lens[i] - arr2 = np.asarray(arr).reshape(slc) - for j, sz in enumerate(lens): - if j != i: - arr2 = arr2.repeat(sz, axis=j) - - ans.insert(0,arr2) - return tuple(ans) - - -#prepare command line interface -parser = argparse.ArgumentParser(prog="geoFromBarycentic", - description='''Generate geom file through \ - Barycentric interpolating seeds file.''', - epilog="requires numpy, and scipy.") -parser.add_argument("seeds", - help="seeds file in DAMASK format:\ - http://damask.mpie.de/Documentation/AsciiTableFormat", - default="test.seeds") -parser.add_argument("-v", "--version", - action="version", - version="%(prog)s 0.1") -parser.add_argument("-g", "--grid", - nargs=3, - help="grid size(mesh resolution, recommend using 2^x)", - default=[32,32,32], - type=int) -parser.add_argument("-s", "--size", - help="physical size of the target volume.", - nargs=3, - default=[1.0,1.0,1.0], - type=float) -parser.add_argument("-o", "--origin", - help="lower left corner of the patch.", - nargs=3, - default=[0.0,0.0,0.0], - type=float) -parser.add_argument('-m', '--homogenization', - help='homogenization index to be used', - default=1, - type=int) -parser.add_argument('-c', '--crystallite', - help='crystallite index to be used', - default=1, - type=int) -parser.add_argument('-p', '--phase', - help='phase index to be used', - default=1, - type=int) -parser.add_argument('-F', '--Favg', - help='reshape the periodicity, not useful for RIM method', - nargs=9, - default=[1.0,0.0,0.0, - 0.0,1.0,0.0, - 0.0,0.0,1.0], - type=float) -parser.add_argument("-G", "--geomFile", - help='the name of the output geom file', - default='seeds.geom', - type=str) -parser.add_argument("-C", "--configFile", - help='output dummy material.config file', - action='store_true', - default=False) -parser.add_argument("-d", "--debug", - help="start debugging script", - action='store_true', - default=False) -parser.add_argument("-S", "--seedsFile", - help="write out resized seeds file", - action='store_true', - default=False) -parser.add_argument("-r", '--addRim', - help="add rim and provide control of face lifting point", - action='store_true', - default=False) -args = parser.parse_args() # get all the arguments right after - -#quick help to user -print "*"*80 -parser.print_help() -print """Sample usage: -./geoFromBarycentic.py 20grains.seeds -g 128 128 128 -S -r; geom_check seeds.geom; seeds_check new_seed.seeds. -""" -print "*"*80 -if (args.debug): - d_print("args are:", parser.parse_args(),separator=True) - -#/\/\/\/\/# -# m a i n # -#\/\/\/\/\# -print "only work for 3D case now, 2D support coming soon..." -print "reading seeds file: {}".format(args.seeds) - -with open(args.seeds, 'r') as f: - rawtext = f.readlines() - n_header = int(rawtext.pop(0).split()[0]) - #record all the seeds position - if (args.addRim): - grid_shift = np.array(args.size) * np.array([OFFSET,OFFSET,OFFSET*2]) - s_coords = np.array([[np.array(float(item))*(1 - OFFSET*2) - for item in line.split()[:3]] + grid_shift - for line in rawtext[n_header:]]) - else: - #no need for shifting with periodicity - s_coords = np.array([[np.array(float(item)) - for item in line.split()[:3]] - for line in rawtext[n_header:]]) - - #record ID of the seeds: int/EulerAngles - if 'microstructure' in rawtext[n_header-1]: - s_id = [int(line.split()[-1]) for line in rawtext[n_header:]] - else: - print "WARNING:" - print "THIS SCRIPT DOES NOT UNDERSTAND HOW TO GROUP CRYSTALLITES." - print "ALL CRYSTAL ORIENTATIONS ARE CONSIDERED TO BE UNIQUE." - print "FOR MORE ACCURATE CONTROL OF SEEDS GROUPING, USE MICROSTRUCTURE ID." - s_id = range(len(s_coords)) - #s_eulers here is just a quick book keeping - s_eulers = np.array([[float(item) for item in line.split()[3:]] - for line in rawtext[n_header:]]) - -if(args.debug): - print d_print("resize point cloud to make space for rim/pan:", - s_coords) - -if(args.addRim): - #add binding box to create rim/pan for the volume where the ID of the seeds is - #unknown - print "Shrining the seeds to {}x in each direction".format(1 - OFFSET*2) - x,y,z = args.size[0],args.size[1],args.size[2] - print "Use circumscribed sphere to place phantom seeds." - r = np.sqrt(x**2+y**2+z**2)/2.0 - BINDBOX = [[0,0,0],[x,0,0],[0,y,0],[x,y,0], - [0,0,z],[x,0,z],[0,y,z],[x,y,z], - [x/2.0+r,y/2, z/2], [x/2.0-r, y/2, z/2], - [x/2, y/2.0+r, z/2], [x/2, y/2.0-r, z/2], - [x/2, y/2, z/2.0-r]] #8 corners + 5 face centers (no top) - print "Adding phantom seeds for RIM generation:" - for point in BINDBOX: - print point - s_coords = np.vstack([s_coords,point]) - s_id.append(PHANTOM_ID) -else: - #The idea here is that we read in each seed point, than duplicate in 3D (make a few copies), - #move on to the next seed point, repeat the same procedure. As for the ID list, we can just use the - #same one. The trick here is use the floor division to find the correct id since we pretty much duplicate - #the same point several times. - Favg = np.array(args.Favg).reshape((3,3)) - x,y,z = args.size[0],args.size[1],args.size[2] - tmp = [] - for seed in s_coords: - tmp += [np.dot(Favg, np.array(seed) + np.array([dx,dy,dz])) - for dz in [-z, 0, z] - for dy in [-y, 0, y] - for dx in [-x, 0, x]] - s_coords = tmp - -if (args.seedsFile): - with open("new_seed.seeds", "w") as f: - outstr = "4\theader\n" - outstr += "grid\ta {}\tb {}\tc {}\n".format(args.grid[0], - args.grid[1], - args.grid[2]) - outstr += "microstructures {}\n".format(len(set(s_id))) - outstr += "x\ty\tz\tmicrostructure" - if (args.addRim): - for i in range(len(s_id)): - outstr += "{}\t{}\t{}\t{}\n".format(s_coords[i][0], - s_coords[i][1], - s_coords[i][2], - s_id[i]) - else: - for i in range(len(s_coords)): - outstr += "{}\t{}\t{}\t{}\n".format(s_coords[i][0], - s_coords[i][1], - s_coords[i][2], - s_id[i//3**3]) - f.write(outstr) - -#triangulate space with given point-clouds -tri = Delaunay(s_coords) - -if(args.debug): - d_print("simplices:", tri.simplices, separator=True) - d_print("vertices:", s_coords[tri.simplices]) - -#populate grid points (only 3D for now) -''' -#populating grid using meshgrid2 -x = (np.arange(args.grid[0])+0.5)*args.size[0]/args.grid[0] -y = (np.arange(args.grid[1])+0.5)*args.size[1]/args.grid[1] -z = (np.arange(args.grid[2])+0.5)*args.size[2]/args.grid[2] -mesh_pts = np.transpose(np.vstack(map(np.ravel, meshgrid2(x, y, z)))) -print mesh_pts -''' -#this is actually faster than using meshgrid2 -mesh_pts = [[(i+0.5)*args.size[0]/args.grid[0], - (j+0.5)*args.size[1]/args.grid[1], - (k+0.5)*args.size[2]/args.grid[2]] - for k in range(args.grid[2]) - for j in range(args.grid[1]) - for i in range(args.grid[0])] - -mesh_ids = [PHANTOM_ID*2]*len(mesh_pts) #initialize grid - -#search ID for each grid point -s_id = np.array(s_id) #allow multi-indexing -mesh_idx = tri.find_simplex(mesh_pts) - -for i, pt in enumerate(mesh_pts): - if mesh_idx[i] < 0: - continue #didn't find any envelop tetrahedron --> something wrong! - #calculate Barycentric coordinates - bary_c = tri.transform[mesh_idx[i],:3,:3].dot(pt-tri.transform[mesh_idx[i],3,:]) - bary_c = np.append(bary_c, 1 - bary_c.sum()) - - if (args.addRim): - tmp_ids = s_id[tri.simplices[mesh_idx[i]]] #rim method - else: - tmp_ids = np.array(s_id[tri.simplices[mesh_idx[i]]//(3**3)]) #kill periodicity through floor division - #print tmp_ids - #print tri.simplices[mesh_idx[i]]//(3**3) - - max_weight = -1960 - for this_id in tmp_ids: - msk = [item==this_id for item in tmp_ids] #find vertex with the same id - tmp_weight = sum([bary_c[j] for j in range(len(bary_c)) if msk[j]]) - if tmp_weight > max_weight: - max_weight = tmp_weight - mesh_ids[i] = this_id - if (args.debug): - d_print("bary_c:",bary_c,separator=True) - d_print("vertex ID:", tmp_ids) - d_print("final ID:", mesh_ids[i]) - -mesh_ids = np.reshape(mesh_ids, (-1, args.grid[0])) - -#write to file -with open(args.geomFile, "w") as f: - outstr = "5\theader\n" - outstr += "grid\ta {}\tb {}\tc {}\n".format(args.grid[0], - args.grid[1], - args.grid[2]) - outstr += "size\tx {}\ty {}\tz {}\n".format(args.size[0], - args.size[1], - args.size[2]) - outstr += "origin\tx {}\ty {}\tz {}\n".format(args.origin[0], - args.origin[1], - args.origin[2]) - outstr += "homogenization\t{}\nmicrostructure\t{}\n".format(args.homogenization, - len(set(s_id))) - for row in mesh_ids: - row = [str(item) for item in list(row)] - outstr += "\t".join(row) + "\n" - f.write(outstr)