Merge branch 'SimplifiedvtkRectilinear' into development

code/constitutive.f90

@@ -399,8 +399,6 @@ function constitutive_homogenizedC(ipc,ip,el)
    plastic_titanmod_homogenizedC
  use plastic_dislotwin, only: &
    plastic_dislotwin_homogenizedC
- use plastic_disloucla, only: &
-   plastic_disloucla_homogenizedC
  use lattice, only: &
    lattice_C66
 
@@ -414,8 +412,6 @@ function constitutive_homogenizedC(ipc,ip,el)
  plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
    case (PLASTICITY_DISLOTWIN_ID) plasticityType
      constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el)
-   case (PLASTICITY_DISLOUCLA_ID) plasticityType
-     constitutive_homogenizedC = plastic_disloucla_homogenizedC(ipc,ip,el)
    case (PLASTICITY_TITANMOD_ID) plasticityType
      constitutive_homogenizedC = plastic_titanmod_homogenizedC (ipc,ip,el)
    case default plasticityType

code/plastic_isotropic.f90

@@ -140,9 +140,10 @@ subroutine plastic_isotropic_init(fileUnit)
    sizeDeltaState
  character(len=65536) :: &
    tag       = '', &
-   outputtag = '', &
    line      = '', &
    extmsg    = ''
+ character(len=64) :: &
+   outputtag = ''
   integer(pInt) :: NipcMyPhase
 
  mainProcess: if (worldrank == 0) then 
@@ -382,8 +383,7 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
    math_mul33xx33, &
    math_transpose33
  use material, only: &
-   phaseAt, phasememberAt, &
+   phasememberAt, &
-   plasticState, &
    material_phase, &
    phase_plasticityInstance
 
@@ -413,7 +413,7 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
    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_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v))                                        ! deviatoric part of 2nd Piola-Kirchhoff stress
  squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33)
@@ -463,8 +463,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
 
@@ -491,7 +490,7 @@ real(pReal) :: &
    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)
@@ -529,8 +528,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
  
@@ -553,7 +551,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
@@ -609,8 +607,7 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
    math_mul6x6
  use material, only: &
    material_phase, &
-   plasticState, &
+   phasememberAt, &
-   phaseAt, phasememberAt, &
    phase_plasticityInstance
 
  implicit none
@@ -634,7 +631,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

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/addDisplacement.py

@@ -141,13 +141,13 @@ for name in filenames:
   if table.label_dimension(options.defgrad) != 9:
     errors.append('deformation gradient "{}" is not a 3x3 tensor.'.format(options.defgrad))
 
-  coordDim = table.label_dimension(options.pos)
+  coordDim = table.label_dimension(options.coords)
   if not 3 >= coordDim >= 1:
-    errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos))
+    errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.coords))
   elif coordDim < 3:
     remarks.append('appending {} dimension{} to coordinates "{}"...'.format(3-coordDim,
                                                                             's' if coordDim < 2 else '',
-                                                                            options.pos))
+                                                                            options.coords))
 
   if remarks != []: damask.util.croak(remarks)
   if errors  != []:
@@ -157,7 +157,7 @@ for name in filenames:
 
 # --------------- figure out size and grid ---------------------------------------------------------
 
-  table.data_readArray([options.defgrad,options.pos])
+  table.data_readArray([options.defgrad,options.coords])
   table.data_rewind()
 
   if len(table.data.shape) < 2: table.data.shape += (1,)                                            # expand to 2D shape
@@ -196,8 +196,8 @@ for name in filenames:
 
   table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
   table.labels_append((['{}_pos'         .format(i+1)   for i in xrange(3)] if options.nodal else []) +
-                       ['{}_avg({}).{}'  .format(i+1,options.defgrad,options.pos) for i in xrange(3)] +
+                       ['{}_avg({}).{}'  .format(i+1,options.defgrad,options.coords) for i in xrange(3)] +
-                       ['{}_fluct({}).{}'.format(i+1,options.defgrad,options.pos) for i in xrange(3)] )
+                       ['{}_fluct({}).{}'.format(i+1,options.defgrad,options.coords) for i in xrange(3)] )
   table.head_write()
 
 # ------------------------------------------ output data -------------------------------------------

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  != []:
@@ -82,16 +84,16 @@ for name in filenames:
  
 # ------------------------------------------ output result ---------------------------------------  
 
-  writer = vtk.vtkXMLPolyDataWriter()
-  writer.SetCompressorTypeToZLib()
-
   if name:
+    writer = vtk.vtkXMLPolyDataWriter()
+    writer.SetCompressorTypeToZLib()
     writer.SetDataModeToBinary()
     writer.SetFileName(os.path.join(os.path.split(name)[0],
                                     os.path.splitext(os.path.split(name)[1])[0] +
                                     '.' + writer.GetDefaultFileExtension()))
   else:
-    writer.SetDataModeToAscii()
+    writer = vtk.vtkDataSetWriter()
+    writer.SetHeader('# powered by '+scriptID)
     writer.WriteToOutputStringOn()
   
   if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(Polydata)
@@ -99,6 +101,6 @@ for name in filenames:
 
   writer.Write()
 
-  if name is None:  sys.stdout.write(writer.GetOutputString())
+  if name is None: sys.stdout.write(writer.GetOutputString()[:writer.GetOutputStringLength()])      # limiting of outputString is fix for vtk <7.0
 
   table.close()

processing/post/vtk_rectilinearGrid.py

@@ -6,7 +6,6 @@ import numpy as np
 import damask
 from optparse import OptionParser
 
-
 scriptName = os.path.splitext(os.path.basename(__file__))[0]
 scriptID   = ' '.join([scriptName,damask.version])
 
@@ -15,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)
 
@@ -29,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()
@@ -47,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
@@ -62,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  != []:
@@ -74,32 +68,24 @@ for name in filenames:
 
 # --------------- figure out size and grid ---------------------------------------------------------
 
-  if isGeom:
+  table.data_readArray(options.pos)
-    info,extra_header = table.head_getGeom()
+  if len(table.data.shape) < 2: table.data.shape += (1,)                                            # expand to 2D shape
-    coords = [np.linspace(info['origin'][i],
+  if table.data.shape[1] < 3:
-                          info['origin'][i]+info['size'][i],
+    table.data = np.hstack((table.data,
-                          num = info['grid'][i]+1,
+                            np.zeros((table.data.shape[0],
-                          endpoint = True,
+                                      3-table.data.shape[1]),dtype='f')))                           # fill coords up to 3D with zeros
-                         ) 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:
-      table.data = np.hstack((table.data,
-                              np.zeros((table.data.shape[0],
-                                        3-table.data.shape[1]),dtype='f')))                           # fill coords up to 3D with zeros
 
-    coords = [np.unique(table.data[:,i]) for i in xrange(3)]
+  coords = [np.unique(table.data[:,i]) for i in xrange(3)]
 
-    if options.mode == 'cell':
+  if options.mode == 'cell':
-      coords = [0.5 * np.array([3.0 * coords[i][0] - coords[i][0 + len(coords[i]) > 1]] + \
+    coords = [0.5 * np.array([3.0 * coords[i][0] - coords[i][0 + len(coords[i]) > 1]] + \
-                               [coords[i][j-1] + coords[i][j] for j in xrange(1,len(coords[i]))] + \
+                             [coords[i][j-1] + coords[i][j] for j in xrange(1,len(coords[i]))] + \
-                               [3.0 * coords[i][-1] - coords[i][-1 - (len(coords[i]) > 1)]]) for i in xrange(3)]
+                             [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)
@@ -132,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()
@@ -144,7 +130,6 @@ for name in filenames:
 
   writer.Write()
 
-  if name is None:  sys.stdout.write(writer.GetOutputString() if name else
+  if name is None: sys.stdout.write(writer.GetOutputString()[:writer.GetOutputStringLength()])      # limiting of outputString is fix for vtk <7.0
-                                     writer.GetOutputString()[0:writer.GetOutputStringLength()])
 
   table.close()

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