Merge remote branch 'origin/development' into pheno+

examples/SpectralMethod/example_postprocess.py

@@ -1,33 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: UTF-8 no BOM -*-
-
-import os
-import glob
-from subprocess import call
-
-geom_name = '20grains16x16x16_tensionX'
-postResults = 'postResults --cr f,p --split --separation x,y,z '+geom_name+'.spectralOut'
-
-sts = call(postResults, shell=True)
-
-os.chdir('./postProc/')
-ascii_files = glob.glob(geom_name+'_inc*.txt')
-print ascii_files
-
-showTable = "showTable -a "
-addCauchy = 'addCauchy '
-addMises = 'addMises -s Cauchy '
-addStrainTensors = "addStrainTensors -0 -v "
-visualize3D = "3Dvisualize -s 'Mises(Cauchy)',1_p  Cauchy "
-
-
-postProc = [addCauchy, addMises, addStrainTensors, visualize3D]
-
-
-for f in ascii_files:
-    print f
-    for p in postProc:
-        p = p+f
-        print p
-        sts = call(p,shell=True)
-

examples/SpectralMethod/example_run.py

@@ -1,17 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: UTF-8 no BOM -*-
-import sys
-
-resolutions = [16,32,64]
-resolution = resolutions[0]
-
-try:
-  resolution = int(sys.argv[1])
-except:
-  pass
-
-if resolution not in resolutions:
-  resolution = resolutions[0]
-  
-from subprocess import call
-call('make run%s'%('x'.join([str(resolution)]*3)), shell=True)

processing/post/vtk_pointcloud.py

@@ -49,7 +49,9 @@ for name in filenames:
   remarks = []
   coordDim = table.label_dimension(options.pos)
   if not 3 >= coordDim >= 1: errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos))
-  elif coordDim < 3:         remarks.append('appending {} dimensions to coordinates "{}"...'.format(3-coordDim,options.pos))
+  elif coordDim < 3:        remarks.append('appending {} dimension{} to coordinates "{}"...'.format(3-coordDim,
+                                                                                                    's' if coordDim < 2 else '',
+                                                                                                    options.pos))
 
   if remarks != []: damask.util.croak(remarks)
   if errors  != []:

processing/post/vtk_rectilinearGrid.py

@@ -14,7 +14,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 # --------------------------------------------------------------------
 
 parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
-Create regular voxel grid from points in an ASCIItable (or geom file).
+Create regular voxel grid from points in an ASCIItable.
 
 """, version = scriptID)
 
@@ -28,15 +28,9 @@ parser.add_option('-p',
                   dest = 'pos',
                   type = 'string', metavar = 'string',
                   help = 'label of coordinates [%default]')
-parser.add_option('-g',
-                  '--geom',
-                  dest = 'geom',
-                  action = 'store_true',
-                  help = 'geom input format')
 
 parser.set_defaults(mode   = 'cell',
                     pos    = 'pos',
-                    geom   = False,
                    )
 
 (options, filenames) = parser.parse_args()
@@ -46,10 +40,9 @@ parser.set_defaults(mode   = 'cell',
 if filenames == []: filenames = [None]
 
 for name in filenames:
-  isGeom = options.geom or (name is not None and name.endswith('.geom'))
   try:    table = damask.ASCIItable(name = name,
                                     buffered = False,
-                                    labeled  = not isGeom,
+                                    labeled  = True,
                                     readonly = True,
                                    )
   except: continue
@@ -61,9 +54,11 @@ for name in filenames:
 
   remarks = []
   errors  = []
-  coordDim = 3 if isGeom else table.label_dimension(options.pos)
+  coordDim = table.label_dimension(options.pos)
   if not 3 >= coordDim >= 1: errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos))
-  elif coordDim < 3:         remarks.append('appending {} dimensions to coordinates "{}"...'.format(3-coordDim,options.pos))
+  elif coordDim < 3:        remarks.append('appending {} dimension{} to coordinates "{}"...'.format(3-coordDim,
+                                                                                                    's' if coordDim < 2 else '',
+                                                                                                    options.pos))
 
   if remarks != []: damask.util.croak(remarks)
   if errors  != []:
@@ -73,14 +68,6 @@ for name in filenames:
 
 # --------------- figure out size and grid ---------------------------------------------------------
 
-  if isGeom:
-    info,extra_header = table.head_getGeom()
-    coords = [np.linspace(info['origin'][i],
-                          info['origin'][i]+info['size'][i],
-                          num = info['grid'][i]+1,
-                          endpoint = True,
-                         ) for i in xrange(3)]
-  else:
   table.data_readArray(options.pos)
   if len(table.data.shape) < 2: table.data.shape += (1,)                                            # expand to 2D shape
   if table.data.shape[1] < 3:
@@ -96,9 +83,9 @@ for name in filenames:
                              [3.0 * coords[i][-1] - coords[i][-1 - (len(coords[i]) > 1)]]) for i in xrange(3)]
 
   grid = np.array(map(len,coords),'i')
-  N = grid.prod() if options.mode == 'point' or isGeom else (grid-1).prod()
+  N = grid.prod() if options.mode == 'point' else (grid-1).prod()
 
-  if not isGeom and N != len(table.data):
+  if N != len(table.data):
     errors.append('data count {} does not match grid {}x{}x{}.'.format(N,*(grid - (options.mode == 'cell')) ))
   if errors  != []:
     damask.util.croak(errors)
@@ -131,7 +118,7 @@ for name in filenames:
     writer.SetDataModeToBinary()
     writer.SetFileName(os.path.join(os.path.split(name)[0],
                                     os.path.splitext(os.path.split(name)[1])[0] +
-                                    ('' if isGeom else '_{}({})'.format(options.pos, options.mode)) +
+                                    '_{}({})'.format(options.pos, options.mode) +
                                     '.' + writer.GetDefaultFileExtension()))
   else:
     writer = vtk.vtkDataSetWriter()

processing/pre/geom_check.sh

@@ -2,8 +2,10 @@
 
 for geom in "$@"
 do
-  vtk_rectilinearGrid \
+  geom_toTable \
-    --geom $geom
+  < $geom \
+  | \
+  vtk_rectilinearGrid > ${geom%.*}.vtk
 
   geom_toTable \
   < $geom \
@@ -11,5 +13,6 @@ do
   vtk_addRectilinearGridData \
     --scalar microstructure \
     --inplace \
-    --vtk ${geom%.*}.vtr
+    --vtk ${geom%.*}.vtk
+  rm ${geom%.*}.vtk
 done

processing/pre/geom_fromEuclideanDistance.py

@@ -1,219 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: UTF-8 no BOM -*-
-
-import os,sys,math,itertools
-import numpy as np
-from scipy import ndimage
-from optparse import OptionParser
-import damask
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-def periodic_3Dpad(array, rimdim=(1,1,1)):
-
-  rimdim = np.array(rimdim,'i')
-  size = np.array(array.shape,'i')
-  padded = np.empty(size+2*rimdim,array.dtype)
-  padded[rimdim[0]:rimdim[0]+size[0],
-         rimdim[1]:rimdim[1]+size[1],
-         rimdim[2]:rimdim[2]+size[2]] = array
-
-  p = np.zeros(3,'i')
-  for side in xrange(3):
-    for p[(side+2)%3] in xrange(padded.shape[(side+2)%3]):
-      for p[(side+1)%3] in xrange(padded.shape[(side+1)%3]):
-        for p[side%3] in xrange(rimdim[side%3]):
-          spot = (p-rimdim)%size
-          padded[p[0],p[1],p[2]] = array[spot[0],spot[1],spot[2]]
-        for p[side%3] in xrange(rimdim[side%3]+size[side%3],size[side%3]+2*rimdim[side%3]):
-          spot = (p-rimdim)%size
-          padded[p[0],p[1],p[2]] = array[spot[0],spot[1],spot[2]]
-  return padded
-
-#--------------------------------------------------------------------------------------------------
-#                                MAIN
-#--------------------------------------------------------------------------------------------------
-
-features = [
-            {'aliens': 1, 'alias': ['boundary','biplane'],},
-            {'aliens': 2, 'alias': ['tripleline',],},
-            {'aliens': 3, 'alias': ['quadruplepoint',],}
-           ]
-
-neighborhoods = {
-                  'neumann':np.array([
-                                        [-1, 0, 0],
-                                        [ 1, 0, 0],
-                                        [ 0,-1, 0],
-                                        [ 0, 1, 0],
-                                        [ 0, 0,-1],
-                                        [ 0, 0, 1],
-                                      ]),
-                  'moore':np.array([
-                                        [-1,-1,-1],
-                                        [ 0,-1,-1],
-                                        [ 1,-1,-1],
-                                        [-1, 0,-1],
-                                        [ 0, 0,-1],
-                                        [ 1, 0,-1],
-                                        [-1, 1,-1],
-                                        [ 0, 1,-1],
-                                        [ 1, 1,-1],
-#
-                                        [-1,-1, 0],
-                                        [ 0,-1, 0],
-                                        [ 1,-1, 0],
-                                        [-1, 0, 0],
-#
-                                        [ 1, 0, 0],
-                                        [-1, 1, 0],
-                                        [ 0, 1, 0],
-                                        [ 1, 1, 0],
-#
-                                        [-1,-1, 1],
-                                        [ 0,-1, 1],
-                                        [ 1,-1, 1],
-                                        [-1, 0, 1],
-                                        [ 0, 0, 1],
-                                        [ 1, 0, 1],
-                                        [-1, 1, 1],
-                                        [ 0, 1, 1],
-                                        [ 1, 1, 1],
-                                      ])
-                }
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
-Produce geom files containing Euclidean distance to grain structural features:
-boundaries, triple lines, and quadruple points.
-
-""", version = scriptID)
-
-parser.add_option('-t','--type',
-                  dest = 'type',
-                  action = 'extend', metavar = '<string LIST>',
-                  help = 'feature type {%s} '%(', '.join(map(lambda x:'|'.join(x['alias']),features))) )
-parser.add_option('-n','--neighborhood',
-                  dest = 'neighborhood',
-                  choices = neighborhoods.keys(), metavar = 'string',
-                  help = 'type of neighborhood {%s} [neumann]'%(', '.join(neighborhoods.keys())))
-parser.add_option('-s', '--scale',
-                  dest = 'scale',
-                  type = 'float', metavar = 'float',
-                  help = 'voxel size [%default]')
-
-parser.set_defaults(type = [],
-                    neighborhood = 'neumann',
-                    scale = 1.0,
-                   )
-
-(options,filenames) = parser.parse_args()
-
-if len(options.type) == 0 or \
-   not set(options.type).issubset(set(list(itertools.chain(*map(lambda x: x['alias'],features))))):
-  parser.error('sleect feature type from (%s).'%(', '.join(map(lambda x:'|'.join(x['alias']),features))) )
-if 'biplane' in options.type and 'boundary' in options.type:
-  parser.error("only one alias out 'biplane' and 'boundary' required")
-  
-feature_list = []
-for i,feature in enumerate(features):
-  for name in feature['alias']:
-    for myType in options.type:
-      if name.startswith(myType):
-        feature_list.append(i)                                                                      # remember selected features
-        break
-
-# --- loop over input files -------------------------------------------------------------------------
-
-if filenames == []: filenames = [None]
-
-for name in filenames:
-  try:
-    table = damask.ASCIItable(name = name,
-                              buffered = False, labeled = False, readonly = True)
-  except: continue
-  damask.util.report(scriptName,name)
-
-# --- interpret header ----------------------------------------------------------------------------
-
-  table.head_read()
-  info,extra_header = table.head_getGeom()
-  
-  damask.util.croak(['grid     a b c:  %s'%(' x '.join(map(str,info['grid']))),
-               'size     x y z:  %s'%(' x '.join(map(str,info['size']))),
-               'origin   x y z:  %s'%(' : '.join(map(str,info['origin']))),
-               'homogenization:  %i'%info['homogenization'],
-               'microstructures: %i'%info['microstructures'],
-              ])
-
-  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.')
-  if errors != []:
-    damask.util.croak(errors)
-    table.close(dismiss = True)
-    continue
-
-# --- read data ------------------------------------------------------------------------------------
-
-  microstructure = table.microstructure_read(info['grid']).reshape(info['grid'],order='F')          # read microstructure
-
-  table.close()
-  
-  neighborhood = neighborhoods[options.neighborhood]
-  convoluted = np.empty([len(neighborhood)]+list(info['grid']+2),'i')
-  structure = periodic_3Dpad(microstructure)
-  
-  for i,p in enumerate(neighborhood):
-    stencil = np.zeros((3,3,3),'i')
-    stencil[1,1,1] = -1
-    stencil[p[0]+1,
-            p[1]+1,
-            p[2]+1] = 1
-    convoluted[i,:,:,:] = ndimage.convolve(structure,stencil)
-  
-  convoluted = np.sort(convoluted,axis = 0)
-  uniques = np.where(convoluted[0,1:-1,1:-1,1:-1] != 0, 1,0)                                        # initialize unique value counter (exclude myself [= 0])
-
-  for i in xrange(1,len(neighborhood)):                                                             # check remaining points in neighborhood
-    uniques += np.where(np.logical_and(
-                           convoluted[i,1:-1,1:-1,1:-1] != convoluted[i-1,1:-1,1:-1,1:-1],          # flip of ID difference detected?
-                           convoluted[i,1:-1,1:-1,1:-1] != 0),                                      # not myself?
-                           1,0)                                                                     # count flip
-
-  for feature in feature_list:
-    try:
-      table = damask.ASCIItable(outname = features[feature]['alias'][0]+'_'+name if name else name,
-                                buffered = False, labeled = False)
-    except: continue
-
-    damask.util.croak(features[feature]['alias'][0])
-      
-    distance = np.where(uniques >= features[feature]['aliens'],0.0,1.0)                             # seed with 0.0 when enough unique neighbor IDs are present
-    distance = ndimage.morphology.distance_transform_edt(distance)*[options.scale]*3
-
-    info['microstructures'] = int(math.ceil(distance.max()))
-
-#--- write header ---------------------------------------------------------------------------------
-
-    table.info_clear()
-    table.info_append(extra_header+[
-      scriptID + ' ' + ' '.join(sys.argv[1:]),
-      "grid\ta {grid[0]}\tb {grid[1]}\tc {grid[2]}".format(grid=info['grid']),
-      "size\tx {size[0]}\ty {size[1]}\tz {size[2]}".format(size=info['size']),
-      "origin\tx {origin[0]}\ty {origin[1]}\tz {origin[2]}".format(origin=info['origin']),
-      "homogenization\t{homog}".format(homog=info['homogenization']),
-      "microstructures\t{microstructures}".format(microstructures=info['microstructures']),
-      ])
-    table.labels_clear()
-    table.head_write()
-    
-# --- write microstructure information ------------------------------------------------------------
-
-    formatwidth = int(math.floor(math.log10(distance.max())+1))
-    table.data = distance.reshape((info['grid'][0],info['grid'][1]*info['grid'][2]),order='F').transpose()
-    table.data_writeArray('%%%ii'%(formatwidth),delimiter=' ')
-    
-#--- output finalization --------------------------------------------------------------------------
-
-    table.close()

processing/pre/seeds_fromTable.py

@@ -1,138 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: UTF-8 no BOM -*-
-
-import os,itertools
-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 = """
-Create seed file by taking microstructure indices from given ASCIItable column.
-White and black-listing of microstructure indices is possible.
-
-Examples:
---white 1,2,5 --index grainID isolates grainID entries of value 1, 2, and 5;
---black 1 --index grainID takes all grainID entries except for value 1.
-
-""", version = scriptID)
-
-parser.add_option('-p',
-                  '--pos', '--seedposition',
-                   dest = 'pos',
-                   type = 'string', metavar = 'string',
-                   help = 'label of coordinates [%default]')
-parser.add_option('--boundingbox',
-                  dest = 'box',
-                  type = 'float', nargs = 6, metavar = ' '.join(['float']*6),
-                  help = 'min (x,y,z) and max (x,y,z) coordinates of bounding box [tight]')
-parser.add_option('-m',
-                  '--microstructure',
-                  dest = 'microstructure',
-                  type = 'string', metavar = 'string',
-                  help = 'label of microstructures [%default]')
-parser.add_option('--weight',
-                  dest = 'weight',
-                  type = 'string', metavar = 'string',
-                  help = 'label of weights [%default]')
-parser.add_option('-w',
-                  '--white',
-                  dest = 'whitelist',
-                  action = 'extend', metavar = '<int LIST>',
-                  help = 'whitelist of microstructure indices')
-parser.add_option('-b',
-                  '--black',
-                  dest = 'blacklist',
-                  action = 'extend', metavar = '<int LIST>',
-                  help = 'blacklist of microstructure indices')
-
-parser.set_defaults(pos            = 'pos',
-                    microstructure = 'microstructure',
-                    weight         = None,
-                   )
-
-(options,filenames) = parser.parse_args()
-
-if options.whitelist is not None: options.whitelist = map(int,options.whitelist)
-if options.blacklist is not None: options.blacklist = map(int,options.blacklist)
-
-# --- loop over input files -------------------------------------------------------------------------
-
-if filenames == []: filenames = [None]
-
-for name in filenames:
-  try:    table = damask.ASCIItable(name = name,
-                                    outname = os.path.splitext(name)[0]+'.seeds' if name else name,
-                                    buffered = False)
-  except: continue
-  damask.util.report(scriptName,name)
-
-  table.head_read()                                                                                 # read ASCII header info
-
-# ------------------------------------------ sanity checks ---------------------------------------  
-
-  missing_labels = table.data_readArray([options.pos,options.microstructure] +
-                                        ([options.weight] if options.weight else []))
-
-  errors = []
-  if len(missing_labels) > 0:
-    errors.append('column{} {} not found'.format('s' if len(missing_labels) > 1 else '',
-                                                ', '.join(missing_labels)))
-  input = {options.pos: 3,
-           options.microstructure: 1,}
-  if options.weight: input.update({options.weight: 1})
-  for label, dim in input.iteritems():
-    if table.label_dimension(label) != dim:
-      errors.append('column {} has wrong dimension'.format(label))
-  
-  if errors != []:
-    damask.util.croak(errors)
-    table.close(dismiss = True)                                                                     # close ASCII table file handles and delete output file
-    continue
-
-# ------------------------------------------ process data ------------------------------------------
-  
-# --- finding bounding box -------------------------------------------------------------------------
-
-  boundingBox = np.array((np.amin(table.data[:,0:3],axis = 0),np.amax(table.data[:,0:3],axis = 0)))
-  if options.box:
-    boundingBox[0,:] = np.minimum(options.box[0:3],boundingBox[0,:])
-    boundingBox[1,:] = np.maximum(options.box[3:6],boundingBox[1,:])
-
-# --- rescaling coordinates ------------------------------------------------------------------------
-
-  table.data[:,0:3] -= boundingBox[0,:]
-  table.data[:,0:3] /= boundingBox[1,:]-boundingBox[0,:]
-
-# --- filtering of grain voxels --------------------------------------------------------------------
-
-  mask = np.logical_and(
-           np.ones_like(table.data[:,3],bool) if options.whitelist is None \
-           else np.in1d(table.data[:,3].ravel(), options.whitelist).reshape(table.data[:,3].shape),
-           np.ones_like(table.data[:,3],bool) if options.blacklist is None \
-           else    np.invert(np.in1d(table.data[:,3].ravel(), options.blacklist).reshape(table.data[:,3].shape))
-          )
-  table.data = table.data[mask]
-
-# ------------------------------------------ assemble header ---------------------------------------  
-
-  table.info = [
-                scriptID,
-                'size {}'.format(' '.join(list(itertools.chain.from_iterable(zip(['x','y','z'],
-                                                                       map(str,boundingBox[1,:]-boundingBox[0,:])))))),
-               ]
-  table.labels_clear()
-  table.labels_append(['1_pos','2_pos','3_pos','microstructure'] +
-                      ['weight'] if options.weight else [])                                        # implicitly switching label processing/writing on
-  table.head_write()
-  
-# ------------------------------------------ output result ---------------------------------------  
-
-  table.data_writeArray()
-  table.close()                                                                                     # close ASCII tables

src/crystallite.f90

@@ -3990,7 +3990,6 @@ subroutine crystallite_orientations
  use plastic_nonlocal, only: &
    plastic_nonlocal_updateCompatibility
 
-
  implicit none
  integer(pInt) &
    c, &                                                                                             !< counter in integration point component loop
@@ -4006,16 +4005,15 @@ subroutine crystallite_orientations
 
  ! --- CALCULATE ORIENTATION AND LATTICE ROTATION ---
 
- 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)
-         orientation = math_RtoQ(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e))))          ! rotational part from polar decomposition as quaternion
+       orientation = math_RtoQ(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e))))
 !$OMP END CRITICAL (polarDecomp)
-         crystallite_rotation(1:4,c,i,e) = lattice_qDisorientation(crystallite_orientation0(1:4,c,i,e), &  ! active rotation from ori0
+       crystallite_rotation(1:4,c,i,e) = lattice_qDisorientation(crystallite_orientation0(1:4,c,i,e), &! active rotation from initial
                                                                   orientation)                         ! to current orientation (with no symmetry)
        crystallite_orientation(1:4,c,i,e) = orientation
  enddo; enddo; enddo
@@ -4025,6 +4023,7 @@ subroutine crystallite_orientations
  ! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL ---
  ! --- we use crystallite_orientation from above, so need a separate loop
   
+ nonlocalPresent: if (any(plasticState%nonLocal)) then
 !$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)