Merge branch 'development' into HDF5-spectral-displacements

.gitlab-ci.yml

@@ -75,7 +75,7 @@ variables:
   MSC:                             "$MSC2019"
   IntelMarc:                       "$IntelCompiler17_8"
   IntelAbaqus:                     "$IntelCompiler16_4"
-  HDF5Marc:                        "HDF5/1.10.4/Intel-17.8"
+  HDF5Marc:                        "HDF5/1.10.5/Intel-17.8"
   # ++++++++++++ Documentation ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   Doxygen1_8_15:                   "Documentation/Doxygen/1.8.15"
   # ------------ Defaults ----------------------------------------------

PRIVATE

@@ -1 +1 @@
-Subproject commit cdf79281330b2cf6e11426685b1d0c999efc24cf
+Subproject commit 214c69be8b51adb39eb7ad25b139727c8b98afce

VERSION

@@ -1 +1 @@
-v2.0.3-860-g0dd14a06
+v2.0.3-922-gcebbfc90

installation/mods_MarcMentat/2019/Marc_tools/include_linux64

@@ -63,6 +63,7 @@ else
  INTEGER_PATH=/$MARC_INTEGER_SIZE
 fi
 
+FCOMP=ifort
 INTELPATH="/opt/intel/compilers_and_libraries_2017/linux"
 
 # find the root directory of the compiler installation:
@@ -103,9 +104,6 @@ if test "$DAMASK_HDF5" = "ON";then
   H5FC="$(h5fc -shlib -show)"
   HDF5_LIB=${H5FC//ifort/}
   FCOMP="$H5FC -DDAMASK_HDF5"
-  echo $FCOMP
-else
-  FCOMP=ifort
 fi
 
 # AEM
@@ -531,7 +529,7 @@ else
   FORT_OPT=" $FORT_OPT -save -zero"
 fi
 if test "$MARCHDF" = "HDF"; then
-  FORT_OPT="$FORT_OPT -DMARCHDF=$MARCHDF $HDF_INCLUDE"
+  FORT_OPT="$FORT_OPT -DMARCHDF=$MARCHDF"
 fi
 
 FORTLOW="$FCOMP $FORT_OPT $PROFILE -O0 $I8FFLAGS -I$MARC_SOURCE/common \
@@ -757,7 +755,7 @@ SECLIBS="-L$MARC_LIB -llapi"
 
 SOLVERLIBS="${BCSSOLVERLIBS} ${VKISOLVERLIBS} ${CASISOLVERLIBS} ${MF2SOLVERLIBS} \
   $MKLLIB -L$MARC_MKL -liomp5 \
-  $MARC_LIB/blas_src.a ${ACSI_LIB}/ACSI_MarcLib.a $KDTREE2_LIB/kdtree2.a $HDF_LIBS $HDF_LIB"
+  $MARC_LIB/blas_src.a ${ACSI_LIB}/ACSI_MarcLib.a $KDTREE2_LIB/kdtree2.a $HDF5_LIB"
 
 SOLVERLIBS_DLL=${SOLVERLIBS}
 if test "$AEM_DLL" -eq 1

processing/post/DADF5_vtk_cells.py

@@ -3,6 +3,7 @@
 import os
 import argparse
 
+import h5py
 import numpy as np
 import vtk
 from vtk.util import numpy_support
@@ -41,20 +42,29 @@ for filename in options.filenames:
   results = damask.DADF5(filename)
   
   if results.structured:                                                                            # for grid solvers use rectilinear grid
-    rGrid = vtk.vtkRectilinearGrid()
+    grid = vtk.vtkRectilineagrid()
     coordArray = [vtk.vtkDoubleArray(),
                   vtk.vtkDoubleArray(),
                   vtk.vtkDoubleArray(),
                  ]
 
-    rGrid.SetDimensions(*(results.grid+1))
+    grid.SetDimensions(*(results.grid+1))
     for dim in [0,1,2]:
       for c in np.linspace(0,results.size[dim],1+results.grid[dim]):
         coordArray[dim].InsertNextValue(c)
 
-    rGrid.SetXCoordinates(coordArray[0])
+    grid.SetXCoordinates(coordArray[0])
-    rGrid.SetYCoordinates(coordArray[1])
+    grid.SetYCoordinates(coordArray[1])
-    rGrid.SetZCoordinates(coordArray[2])
+    grid.SetZCoordinates(coordArray[2])
+  else:
+    nodes = vtk.vtkPoints()
+    with h5py.File(filename) as f:
+      nodes.SetData(numpy_support.numpy_to_vtk(f['/geometry/x_n'][()],deep=True))
+      grid = vtk.vtkUnstructuredGrid()
+      grid.SetPoints(nodes)
+      grid.Allocate(f['/geometry/T_c'].shape[0])
+      for i in f['/geometry/T_c']:
+        grid.InsertNextCell(vtk.VTK_HEXAHEDRON,8,i-1)
 
 
   for i,inc in enumerate(results.iter_visible('increments')):
@@ -75,7 +85,7 @@ for filename in options.filenames:
             shape = [array.shape[0],np.product(array.shape[1:])]
             vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE))
             vtk_data[-1].SetName('1_'+x[0].split('/',1)[1])
-            rGrid.GetCellData().AddArray(vtk_data[-1])
+            grid.GetCellData().AddArray(vtk_data[-1])
         else:
           x = results.get_dataset_location(label)
           if len(x) == 0:
@@ -84,7 +94,7 @@ for filename in options.filenames:
           shape = [array.shape[0],np.product(array.shape[1:])]
           vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE))
           vtk_data[-1].SetName('1_'+x[0].split('/',1)[1])
-          rGrid.GetCellData().AddArray(vtk_data[-1])
+          grid.GetCellData().AddArray(vtk_data[-1])
     
     results.set_visible('constituents',  False)
     results.set_visible('materialpoints',True)
@@ -99,7 +109,7 @@ for filename in options.filenames:
             shape = [array.shape[0],np.product(array.shape[1:])]
             vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE))
             vtk_data[-1].SetName('1_'+x[0].split('/',1)[1])
-            rGrid.GetCellData().AddArray(vtk_data[-1])
+            grid.GetCellData().AddArray(vtk_data[-1])
         else:
           x = results.get_dataset_location(label)
           if len(x) == 0:
@@ -108,10 +118,10 @@ for filename in options.filenames:
           shape = [array.shape[0],np.product(array.shape[1:])]
           vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE))
           vtk_data[-1].SetName('1_'+x[0].split('/',1)[1])
-          rGrid.GetCellData().AddArray(vtk_data[-1])
+          grid.GetCellData().AddArray(vtk_data[-1])
           
-    if results.structured:
+    writer = vtk.vtkXMLRectilineagridWriter() if results.structured else \
-      writer = vtk.vtkXMLRectilinearGridWriter()
+             vtk.vtkXMLUnstructuredGridWriter()
 
     results.set_visible('constituents',  False)
     results.set_visible('materialpoints',False)
@@ -128,7 +138,6 @@ for filename in options.filenames:
     writer.SetCompressorTypeToZLib()
     writer.SetDataModeToBinary()
     writer.SetFileName(os.path.join(dirname,file_out))
-    if results.structured:
+    writer.SetInputData(grid)
-      writer.SetInputData(rGrid)
     
     writer.Write()

processing/post/DADF5_vtk_points.py

@@ -0,0 +1,121 @@
+#!/usr/bin/env python3
+
+import os
+import argparse
+
+import numpy as np
+import vtk
+from vtk.util import numpy_support
+
+import damask
+
+scriptName = os.path.splitext(os.path.basename(__file__))[0]
+scriptID   = ' '.join([scriptName,damask.version])
+
+# --------------------------------------------------------------------
+#                                MAIN
+# --------------------------------------------------------------------
+parser = argparse.ArgumentParser()
+
+#ToDo:  We need to decide on a way of handling arguments of variable lentght
+#https://stackoverflow.com/questions/15459997/passing-integer-lists-to-python
+
+#parser.add_argument('--version', action='version', version='%(prog)s {}'.format(scriptID))
+parser.add_argument('filenames', nargs='+',
+                    help='DADF5 files')
+parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string',
+                    help='name of subdirectory relative to the location of the DADF5 file to hold output')
+parser.add_argument('--mat', nargs='+',
+                    help='labels for materialpoint',dest='mat')
+parser.add_argument('--con', nargs='+',
+                    help='labels for constituent',dest='con')
+
+options = parser.parse_args()
+
+if options.mat is None: options.mat=[]
+if options.con is None: options.con=[]
+
+# --- loop over input files ------------------------------------------------------------------------
+
+for filename in options.filenames:
+  results = damask.DADF5(filename)
+  
+  Points   = vtk.vtkPoints()
+  Vertices = vtk.vtkCellArray()  
+  for c in results.cell_coordinates():
+    pointID = Points.InsertNextPoint(c)
+    Vertices.InsertNextCell(1)
+    Vertices.InsertCellPoint(pointID)
+
+  Polydata = vtk.vtkPolyData()
+  Polydata.SetPoints(Points)
+  Polydata.SetVerts(Vertices)
+  Polydata.Modified()
+  
+  for i,inc in enumerate(results.iter_visible('increments')):
+    print('Output step {}/{}'.format(i+1,len(results.increments)))
+    vtk_data = []
+    
+    results.set_visible('materialpoints',False)
+    results.set_visible('constituents',  True)
+    for label in options.con:
+      
+      for p in results.iter_visible('con_physics'):
+        if p != 'generic':
+          for c in results.iter_visible('constituents'):
+            x = results.get_dataset_location(label)
+            if len(x) == 0:
+              continue
+            array = results.read_dataset(x,0)
+            shape = [array.shape[0],np.product(array.shape[1:])]
+            vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE))
+            vtk_data[-1].SetName('1_'+x[0].split('/',1)[1])
+            Polydata.GetCellData().AddArray(vtk_data[-1])
+        else:
+          x = results.get_dataset_location(label)
+          if len(x) == 0:
+            continue
+          array = results.read_dataset(x,0)
+          shape = [array.shape[0],np.product(array.shape[1:])]
+          vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE))
+          vtk_data[-1].SetName('1_'+x[0].split('/',1)[1])
+          Polydata.GetCellData().AddArray(vtk_data[-1])
+    
+    results.set_visible('constituents',  False)
+    results.set_visible('materialpoints',True)
+    for label in options.mat:       
+      for p in results.iter_visible('mat_physics'):
+        if p != 'generic':
+          for m in results.iter_visible('materialpoints'):
+            x = results.get_dataset_location(label)
+            if len(x) == 0:
+              continue
+            array = results.read_dataset(x,0)
+            shape = [array.shape[0],np.product(array.shape[1:])]
+            vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE))
+            vtk_data[-1].SetName('1_'+x[0].split('/',1)[1])
+            Polydata.GetCellData().AddArray(vtk_data[-1])
+        else:
+          x = results.get_dataset_location(label)
+          if len(x) == 0:
+            continue
+          array = results.read_dataset(x,0)
+          shape = [array.shape[0],np.product(array.shape[1:])]
+          vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE))
+          vtk_data[-1].SetName('1_'+x[0].split('/',1)[1])
+          Polydata.GetCellData().AddArray(vtk_data[-1])
+          
+    writer = vtk.vtkXMLPolyDataWriter()
+
+
+    dirname  = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir))
+    if not os.path.isdir(dirname):
+      os.mkdir(dirname,0o755)
+    file_out = '{}_{}.{}'.format(os.path.splitext(os.path.split(filename)[-1])[0],inc,writer.GetDefaultFileExtension())
+    
+    writer.SetCompressorTypeToZLib()
+    writer.SetDataModeToBinary()
+    writer.SetFileName(os.path.join(dirname,file_out))
+    writer.SetInputData(Polydata)
+    
+    writer.Write()

processing/post/DADF5toDREAM3D.py

@@ -0,0 +1,149 @@
+#!/usr/bin/env python3
+
+import argparse
+import os
+
+import h5py
+import numpy as np
+
+import damask
+
+class AttributeManagerNullterm(h5py.AttributeManager): 
+  """
+  Attribute management for DREAM.3D hdf5 files.
+  
+  String attribute values are stored as fixed-length string with NULLTERM
+  
+  References
+  ----------
+    https://stackoverflow.com/questions/38267076
+    https://stackoverflow.com/questions/52750232
+
+  """ 
+
+  def create(self, name, data, shape=None, dtype=None):
+    if isinstance(data,str):
+      tid = h5py.h5t.C_S1.copy()
+      tid.set_size(len(data + ' '))
+      super().create(name=name,data=data+' ',dtype = h5py.Datatype(tid))
+    else:
+      super().create(name=name,data=data,shape=shape,dtype=dtype)
+     
+
+h5py._hl.attrs.AttributeManager = AttributeManagerNullterm # 'Monkey patch'
+
+
+# --------------------------------------------------------------------
+# Crystal structure specifications
+# --------------------------------------------------------------------
+Crystal_structures = {'fcc': 1,
+                      'bcc': 1,
+                      'hcp': 0,
+                      'bct': 7,
+                      'ort': 6} #TODO: is bct Tetragonal low/Tetragonal high?
+Phase_types = {'Primary': 0} #further additions to these can be done by looking at 'Create Ensemble Info' filter
+
+
+# --------------------------------------------------------------------
+#                                MAIN
+# --------------------------------------------------------------------
+parser = argparse.ArgumentParser(description='Creating a file for DREAM3D from DAMASK data')
+parser.add_argument('filenames',nargs='+',help='HDF5 based output file')
+parser.add_argument('--inc',nargs='+',help='Increment for which DREAM3D to be used, eg. 00025',type=int)
+parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string',
+                    help='name of subdirectory to hold output')
+
+options = parser.parse_args()
+
+# --------------------------------------------------------------------
+# loop over input files
+for filename in options.filenames:
+  f = damask.DADF5(filename)  #DAMASK output file
+  count = 0
+  for increment in f.increments:
+    if int(increment[3:]) not in options.inc: 
+      count = count + 1
+      continue
+
+    #-------output file creation-------------------------------------
+    dirname  = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir))
+    print(dirname)
+    try:
+      os.mkdir(dirname)
+    except FileExistsError:
+      pass
+
+    o = h5py.File(dirname + '/' + os.path.splitext(filename)[0] + '_{}.dream3D'.format(increment),'w')
+    #-----------------------------------------------------------------
+    o.attrs['DADF5toDREAM3D'] = '1.0'
+    o.attrs['FileVersion']    = '7.0' 
+    #-----------------------------------------------------------------
+    
+
+    for g in ['DataContainerBundles','Pipeline']: # empty groups (needed)
+      o.create_group(g)
+
+    data_container_label = 'DataContainers/ImageDataContainer'        
+    cell_data_label      = data_container_label + '/CellData'
+
+
+    # Phase information of DREAM.3D is constituent ID in DAMASK
+    o[cell_data_label + '/Phases'] = f.get_constituent_ID().reshape(tuple(f.grid)+(1,))  
+    # Data quaternions
+    DAMASK_quaternion = f.read_dataset(f.get_dataset_location('orientation'),0)
+    DREAM_3D_quaternion = np.empty((np.prod(f.grid),4),dtype=np.float32)
+    # Convert: DAMASK uses P = -1, DREAM.3D uses P = +1. Also change position of imagninary part
+    DREAM_3D_quaternion = np.hstack((-DAMASK_quaternion['x'],-DAMASK_quaternion['y'],-DAMASK_quaternion['z'],
+                                      DAMASK_quaternion['w']))
+    o[cell_data_label + '/Quats'] = DREAM_3D_quaternion.reshape(tuple(f.grid)+(4,))
+    
+    # Attributes to CellData group
+    o[cell_data_label].attrs['AttributeMatrixType'] = np.array([3],np.uint32)
+    o[cell_data_label].attrs['TupleDimensions']     = f.grid.astype(np.uint64)
+  
+    # Common Attributes for groups in CellData
+    for group in ['/Phases','/Quats']:
+      o[cell_data_label + group].attrs['DataArrayVersion']      = np.array([2],np.int32)
+      o[cell_data_label + group].attrs['Tuple Axis Dimensions'] = 'x={},y={},z={}'.format(*f.grid)
+      
+    # phase attributes
+    o[cell_data_label + '/Phases'].attrs['ComponentDimensions'] = np.array([1],np.uint64)
+    o[cell_data_label + '/Phases'].attrs['ObjectType']          = 'DataArray<int32_t>'
+    
+    # Quats attributes
+    o[cell_data_label + '/Quats'].attrs['ComponentDimensions'] = np.array([4],np.uint64)
+    o[cell_data_label + '/Quats'].attrs['ObjectType']          = 'DataArray<float>'        
+
+    # Create EnsembleAttributeMatrix
+    ensemble_label = data_container_label + '/EnsembleAttributeMatrix' 
+    
+    # Data CrystalStructures
+    o[ensemble_label + '/CrystalStructures'] = np.uint32(np.array([999,\
+                                                   Crystal_structures[f.get_crystal_structure()]])).reshape((2,1))
+    o[ensemble_label + '/PhaseTypes']        = np.uint32(np.array([999,Phase_types['Primary']])).reshape((2,1))    # ToDo
+   
+    # Attributes Ensemble Matrix
+    o[ensemble_label].attrs['AttributeMatrixType'] = np.array([11],np.uint32)
+    o[ensemble_label].attrs['TupleDimensions']     = np.array([2], np.uint64)
+    
+    # Attributes for data in Ensemble matrix
+    for group in ['CrystalStructures','PhaseTypes']: # 'PhaseName' not required MD: But would be nice to take the phase name mapping
+      o[ensemble_label+'/'+group].attrs['ComponentDimensions']   = np.array([1],np.uint64)
+      o[ensemble_label+'/'+group].attrs['Tuple Axis Dimensions'] = 'x=2'
+      o[ensemble_label+'/'+group].attrs['DataArrayVersion']      = np.array([2],np.int32)
+      o[ensemble_label+'/'+group].attrs['ObjectType']            = 'DataArray<uint32_t>'
+      o[ensemble_label+'/'+group].attrs['TupleDimensions']       = np.array([2],np.uint64)
+      
+  
+    # Create geometry info
+    geom_label = data_container_label + '/_SIMPL_GEOMETRY'
+    
+    o[geom_label + '/DIMENSIONS'] = np.int64(f.grid)
+    o[geom_label + '/ORIGIN']     = np.float32(np.zeros(3))
+    o[geom_label + '/SPACING']    = np.float32(f.size)
+  
+    o[geom_label].attrs['GeometryName']     = 'ImageGeometry'
+    o[geom_label].attrs['GeometryTypeName'] = 'ImageGeometry'
+    o[geom_label].attrs['GeometryType']          = np.array([0],np.uint32) 
+    o[geom_label].attrs['SpatialDimensionality'] = np.array([3],np.uint32) 
+    o[geom_label].attrs['UnitDimensionality']    = np.array([3],np.uint32) 

python/damask/dadf5.py

@@ -309,6 +309,19 @@ class DADF5():
     return path
     
     
+  def get_constituent_ID(self,c=0):
+    """Pointwise constituent ID."""
+    with h5py.File(self.filename,'r') as f:
+      names = f['/mapping/cellResults/constituent']['Name'][:,c].astype('str')
+    return np.array([int(n.split('_')[0]) for n in names.tolist()],dtype=np.int32)
+
+
+  def get_crystal_structure(self):                                                                  # ToDo: extension to multi constituents/phase
+    """Info about the crystal structure."""
+    with h5py.File(self.filename,'r') as f:
+      return f[self.get_dataset_location('orientation')[0]].attrs['Lattice'].astype('str')          # np.bytes_ to string
+
+
   def read_dataset(self,path,c):
     """
     Dataset for all points/cells.
@@ -318,7 +331,7 @@ class DADF5():
     with h5py.File(self.filename,'r') as f:
       shape = (self.Nmaterialpoints,) + np.shape(f[path[0]])[1:]
       if len(shape) == 1: shape = shape +(1,)
-      dataset = np.full(shape,np.nan)
+      dataset = np.full(shape,np.nan,dtype=np.dtype(f[path[0]]))
       for pa in path:
         label = pa.split('/')[2]
         
@@ -345,6 +358,20 @@ class DADF5():
     return dataset
 
 
+  def cell_coordinates(self):
+    """Initial coordinates of the cell centers."""
+    if self.structured:
+      delta = self.size/self.grid*0.5
+      z, y, x = np.meshgrid(np.linspace(delta[2],self.size[2]-delta[2],self.grid[2]),
+                            np.linspace(delta[1],self.size[1]-delta[1],self.grid[1]),
+                            np.linspace(delta[0],self.size[0]-delta[0],self.grid[0]),
+                           )
+      return np.concatenate((x[:,:,:,None],y[:,:,:,None],y[:,:,:,None]),axis = 3).reshape([np.product(self.grid),3])
+    else:
+      with h5py.File(self.filename,'r') as f:
+        return f['geometry/x_c'][()]
+
+
   def add_Cauchy(self,P='P',F='F'):
     """
     Adds Cauchy stress calculated from 1st Piola-Kirchhoff stress and deformation gradient.

src/CPFEM.f90

@@ -381,6 +381,7 @@ subroutine CPFEM_results(inc,time)
   call results_addIncrement(inc,time)
   call constitutive_results
   call crystallite_results
+  call homogenization_results
   call results_removeLink('current') ! ToDo: put this into closeJobFile
   call results_closeJobFile
 #endif

src/DAMASK_abaqus.f

@@ -313,8 +313,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
    
  call CPFEM_general(computationMode,usePingPong,dfgrd0,dfgrd1,temperature,dtime,noel,npt,stress_h,ddsdde_h)
 
-!     Mandel:              11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
+!     DAMASK:            11, 22, 33, 12, 23, 13
-!     straight:            11, 22, 33, 12, 23, 13
 !     ABAQUS explicit:   11, 22, 33, 12, 23, 13
 !     ABAQUS implicit:   11, 22, 33, 12, 13, 23
 !     ABAQUS implicit:   11, 22, 33, 12

src/DAMASK_interface.f90

@@ -13,7 +13,7 @@
 #define INTEL_MIN 1600
 #define PETSC_MAJOR 3
 #define PETSC_MINOR_MIN 10
-#define PETSC_MINOR_MAX 11
+#define PETSC_MINOR_MAX 12
 
 module DAMASK_interface
   use, intrinsic :: iso_fortran_env

src/DAMASK_marc.f90

@@ -265,8 +265,8 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
          call debug_reset()                                                                         ! resets debugging
          outdatedFFN1  = .false.
          cycleCounter  = cycleCounter + 1
-         mesh_cellnode = mesh_build_cellnodes()                                                     ! update cell node coordinates
+         !mesh_cellnode = mesh_build_cellnodes()                                                     ! update cell node coordinates
-         call mesh_build_ipCoordinates()                                                            ! update ip coordinates
+         !call mesh_build_ipCoordinates()                                                            ! update ip coordinates
        endif
        if (outdatedByNewInc) then
          computationMode = ior(computationMode,CPFEM_AGERESULTS)                                    ! calc and age results
@@ -315,9 +315,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
  lastLovl = lovl                                                                                    ! record lovl
 
  call CPFEM_general(computationMode,usePingPong,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde)
-!     Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
-!     Marc:   11, 22, 33, 12, 23, 13
-!     Marc:   11, 22, 33, 12
 
  d = ddsdde(1:ngens,1:ngens)
  s = stress(1:ndi+nshear)

src/HDF5_utilities.f90

@@ -112,14 +112,10 @@ subroutine HDF5_utilities_init
   call h5open_f(hdferr)
   if (hdferr < 0) call IO_error(1,ext_msg='HDF5_Utilities_init: h5open_f')
 
-#ifndef Marc4DAMASK
-! This test should ensure that integer size matches. For some reasons, the HDF5 libraries
-! that come with MSC.Marc>=2019 seem to be of 4byte even though it is a 8byte Marc version
   call h5tget_size_f(H5T_NATIVE_INTEGER,typeSize, hdferr)
   if (hdferr < 0) call IO_error(1,ext_msg='HDF5_Utilities_init: h5tget_size_f (int)')
   if (int(bit_size(0),SIZE_T)/=typeSize*8) &
     call IO_error(0,ext_msg='Default integer size does not match H5T_NATIVE_INTEGER')
-#endif
 
   call h5tget_size_f(H5T_NATIVE_DOUBLE,typeSize, hdferr)
   if (hdferr < 0) call IO_error(1,ext_msg='HDF5_Utilities_init: h5tget_size_f (double)')

src/IO.f90

@@ -44,7 +44,6 @@ module IO
     IO_extractValue, &
     IO_countDataLines
 #elif defined(Marc4DAMASK)
-    IO_skipChunks, &
     IO_fixedNoEFloatValue, &
     IO_fixedIntValue, &
     IO_countNumericalDataLines
@@ -1023,27 +1022,6 @@ integer function IO_countNumericalDataLines(fileUnit)
   backspace(fileUnit)
 
 end function IO_countNumericalDataLines
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief reads file to skip (at least) N chunks (may be over multiple lines)
-!--------------------------------------------------------------------------------------------------
-subroutine IO_skipChunks(fileUnit,N)
-
-  integer, intent(in)  :: fileUnit, &                                                               !< file handle
-                          N                                                                         !< minimum number of chunks to skip
- 
-  integer              :: remainingChunks
-  character(len=65536) :: line
- 
-  line = ''
-  remainingChunks = N
- 
-  do while (trim(line) /= IO_EOF .and. remainingChunks > 0)
-    line = IO_read(fileUnit)
-    remainingChunks = remainingChunks - (size(IO_stringPos(line))-1)/2
-  enddo
-end subroutine IO_skipChunks
 #endif
 
 

src/commercialFEM_fileList.f90

@@ -14,11 +14,11 @@
 #include "Lambert.f90"
 #include "rotations.f90"
 #include "FEsolving.f90"
-#include "geometry_plastic_nonlocal.f90"
 #include "element.f90"
 #include "mesh_base.f90"
 #include "HDF5_utilities.f90"
 #include "results.f90"
+#include "geometry_plastic_nonlocal.f90"
 #include "discretization.f90"
 #ifdef Abaqus
 #include "mesh_abaqus.f90"

src/crystallite.f90

@@ -2090,7 +2090,7 @@ end subroutine setConvergenceFlag
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Standard forwarding of state as state = state0 + dotState * (delta t)
+!> @brief Standard forwarding of state as state = state0 + dotState * (delta t)                        comment seems wrong!
 !--------------------------------------------------------------------------------------------------
 subroutine update_stress(timeFraction)
 

src/element.f90

@@ -146,11 +146,11 @@ module element
        8  & ! 3D 8node
     ]                                                                                               !< number of cell nodes in a specific cell type
 
-  ! *** FE_ipNeighbor ***
+  ! *** IPneighbor ***
-  ! is a list of the neighborhood of each IP.
+  ! list of the neighborhood of each IP.
   ! It is sorted in (local) +x,-x, +y,-y, +z,-z direction.
-  ! Positive integers denote an intra-FE IP identifier.
+  ! Positive integers denote an intra-element IP identifier.
-  ! Negative integers denote the interface behind which the neighboring (extra-FE) IP will be located.
+  ! Negative integers denote the interface behind which the neighboring (extra-element) IP will be located.
  
   integer, dimension(nIPneighbor(cellType(1)),nIP(1)), parameter, private :: IPneighbor1 = &
     reshape([&
@@ -256,10 +256,6 @@ module element
       -3,26,-4,24,-6,18  &
      ],[nIPneighbor(cellType(10)),nIP(10)])
 
-! MD: probably not needed END
-! --------------------------------------------------------------------------------------------------
-
-
  
   integer, dimension(nNode(1),NcellNode(geomType(1))), parameter :: cellNodeParentNodeWeights1 = &
     reshape([&
@@ -757,7 +753,7 @@ subroutine tElement_init(self,elemType)
        self%cell        = CELL10
    end select
 
-   self%NcellNodesPerCell = NCELLNODEPERCELL(self%cellType)
+   self%NcellnodesPerCell = NCELLNODEPERCELL(self%cellType)
   
    select case(self%cellType)
      case(1)

src/geometry_plastic_nonlocal.f90

@@ -7,6 +7,7 @@
 !--------------------------------------------------------------------------------------------------
 module geometry_plastic_nonlocal
   use prec
+  use results
 
   implicit none
   private
@@ -32,6 +33,7 @@ module geometry_plastic_nonlocal
     geometry_plastic_nonlocal_setIPvolume, &
     geometry_plastic_nonlocal_setIParea, &
     geometry_plastic_nonlocal_setIPareaNormal, &
+    geometry_plastic_nonlocal_results, &
     geometry_plastic_nonlocal_disable
     
 contains
@@ -112,4 +114,45 @@ subroutine geometry_plastic_nonlocal_disable
   
 end subroutine geometry_plastic_nonlocal_disable
 
+
+!---------------------------------------------------------------------------------------------------
+!> @brief Writes geometry data to results file
+!---------------------------------------------------------------------------------------------------
+subroutine geometry_plastic_nonlocal_results
+  
+  integer,     dimension(:),   allocatable :: shp
+
+#if defined(DAMASK_HDF5)
+  call results_openJobFile
+
+  writeVolume: block
+    real(pReal), dimension(:), allocatable :: temp
+    shp = shape(geometry_plastic_nonlocal_IPvolume0)
+    temp = reshape(geometry_plastic_nonlocal_IPvolume0,[shp(1)*shp(2)])
+    call results_writeDataset('geometry',temp,'v_0',&
+                              'initial cell volume','m³')
+  end block writeVolume
+
+  writeAreas: block
+    real(pReal), dimension(:,:), allocatable :: temp
+    shp = shape(geometry_plastic_nonlocal_IParea0)
+    temp = reshape(geometry_plastic_nonlocal_IParea0,[shp(1),shp(2)*shp(3)])
+    call results_writeDataset('geometry',temp,'a_0',&
+                              'initial cell face area','m²')
+  end block writeAreas
+
+  writeNormals: block
+    real(pReal), dimension(:,:,:), allocatable :: temp
+    shp = shape(geometry_plastic_nonlocal_IPareaNormal0)
+    temp = reshape(geometry_plastic_nonlocal_IPareaNormal0,[shp(1),shp(2),shp(3)*shp(4)])
+    call results_writeDataset('geometry',temp,'n_0',&
+                              'initial cell face normals','-',transposed=.false.)
+  end block writeNormals
+
+
+  call results_closeJobFile
+#endif
+  
+end subroutine geometry_plastic_nonlocal_results
+
 end module geometry_plastic_nonlocal

src/lattice.f90

@@ -228,86 +228,86 @@ module lattice
   real(pReal), dimension(4+4,LATTICE_HEX_NSLIP), parameter :: &
     LATTICE_HEX_SYSTEMSLIP = reshape(real([&
      ! Slip direction     Plane normal
-     ! Basal systems <11.0>{00.1} (independent of c/a-ratio, Bravais notation (4 coordinate base))
+     ! Basal systems <-1-1.0>{00.1} (independent of c/a-ratio, Bravais notation (4 coordinate base))
        2, -1, -1,  0,     0,  0,  0,  1, &
       -1,  2, -1,  0,     0,  0,  0,  1, &
       -1, -1,  2,  0,     0,  0,  0,  1, &
-     ! 1st type prismatic systems <11.0>{10.0}  (independent of c/a-ratio)
+     ! 1st type prismatic systems <-1-1.0>{1-1.0}  (independent of c/a-ratio)
        2, -1, -1,  0,     0,  1, -1,  0, &
       -1,  2, -1,  0,    -1,  0,  1,  0, &
       -1, -1,  2,  0,     1, -1,  0,  0, &
-     ! 2nd type prismatic systems <10.0>{11.0} -- a slip; plane normals independent of c/a-ratio
+     ! 2nd type prismatic systems <-11.0>{11.0} -- a slip; plane normals independent of c/a-ratio
-       0,  1,  -1, 0,     2, -1, -1,  0, &
+      -1,  1,  0,  0,     1,  1, -2,  0, &
-      -1,  0,  1,  0,    -1,  2, -1,  0, &
+       0, -1,  1,  0,    -2,  1,  1,  0, &
-       1, -1,  0,  0,    -1, -1,  2,  0,  &
+       1,  0, -1,  0,     1, -2,  1,  0, &
-     ! 1st type 1st order pyramidal systems <11.0>{-11.1} -- plane normals depend on the c/a-ratio
+     ! 1st type 1st order pyramidal systems <-1-1.0>{-11.1} -- plane normals depend on the c/a-ratio
-       2, -1, -1,  0,     0,  1, -1,  1, &
+      -1,  2, -1,  0,     1,  0, -1,  1, &
-      -1,  2, -1,  0,    -1,  0,  1,  1, &
+      -2,  1,  1,  0,     0,  1, -1,  1, &
-      -1, -1,  2,  0,     1, -1,  0,  1, &
+      -1, -1,  2,  0,    -1,  1,  0,  1, &
-       1,  1, -2,  0,    -1,  1,  0,  1, &
+       1, -2,  1,  0,    -1,  0,  1,  1, &
-      -2,  1,  1,  0,     0, -1,  1,  1, &
+       2, -1, -1,  0,     0, -1,  1,  1, &
-       1, -2,  1,  0,     1,  0, -1,  1, &
+       1,  1, -2,  0,     1, -1,  0,  1, &
      ! pyramidal system: c+a slip <11.3>{-10.1} -- plane normals depend on the c/a-ratio
-       2, -1, -1,  3,    -1,  1,  0,  1, &
-       1, -2,  1,  3,    -1,  1,  0,  1, &
-      -1, -1,  2,  3,     1,  0, -1,  1, &
       -2,  1,  1,  3,     1,  0, -1,  1, &
-      -1,  2, -1,  3,     0, -1,  1,  1, &
+      -1, -1,  2,  3,     1,  0, -1,  1, &
-       1,  1, -2,  3,     0, -1,  1,  1, &
-      -2,  1,  1,  3,     1, -1,  0,  1, &
-      -1,  2, -1,  3,     1, -1,  0,  1, &
-       1,  1, -2,  3,    -1,  0,  1,  1, &
-       2, -1, -1,  3,    -1,  0,  1,  1, &
-       1, -2,  1,  3,     0,  1, -1,  1, &
       -1, -1,  2,  3,     0,  1, -1,  1, &
+       1, -2,  1,  3,     0,  1, -1,  1, &
+       1, -2,  1,  3,    -1,  1,  0,  1, &
+       2, -1, -1,  3,    -1,  1,  0,  1, &
+       2, -1, -1,  3,    -1,  0,  1,  1, &
+       1,  1, -2,  3,    -1,  0,  1,  1, &
+       1,  1, -2,  3,     0, -1,  1,  1, &
+      -1,  2, -1,  3,     0, -1,  1,  1, &
+      -1,  2, -1,  3,     1, -1,  0,  1, &
+      -2,  1,  1,  3,     1, -1,  0,  1, &
      ! pyramidal system: c+a slip <11.3>{-1-1.2} -- as for hexagonal ice (Castelnau et al. 1996, similar to twin system found below)
-       2, -1, -1,  3,    -2,  1,  1,  2, & ! sorted according to similar twin system
+      -1, -1,  2,  3,     1,  1, -2,  2, & ! <11.3>{-1-1.2} shear = 2((c/a)^2-2)/(3 c/a) 
-      -1,  2, -1,  3,     1, -2,  1,  2, & ! <11.3>{-1-1.2} shear = 2((c/a)^2-2)/(3 c/a)
-      -1, -1,  2,  3,     1,  1, -2,  2, &
-      -2,  1,  1,  3,     2, -1, -1,  2, &
        1, -2,  1,  3,    -1,  2, -1,  2, &
-       1,  1, -2,  3,    -1, -1,  2,  2  &
+       2, -1, -1,  3,    -2,  1,  1,  2, &
-      ],pReal),shape(LATTICE_HEX_SYSTEMSLIP))                                                       !< slip systems for hex sorted by A. Alankar & P. Eisenlohr
+       1,  1, -2,  3,    -1, -1,  2,  2, &
+      -1,  2, -1,  3,     1, -2,  1,  2, &
+      -2,  1,  1,  3,     2, -1, -1,  2  &
+      ],pReal),shape(LATTICE_HEX_SYSTEMSLIP))                                                       !< slip systems for hex, sorted by P. Eisenlohr CCW around <c> starting next to a_1 axis
  
   character(len=*), dimension(6), parameter :: LATTICE_HEX_SLIPFAMILY_NAME = &
-    ['<1 1 . 1>{0 0 . 1}  ', &
+    ['< 1 1 . 0>{ 0  0 . 1}', &
-     '<1 1 . 1>{1 0 . 0}  ', &
+     '< 1 1 . 0>{ 1  0 . 0}', &
-     '<1 0 . 0>{1 1 . 0}  ', &
+     '<-1 1 . 0>{ 1  1 . 0}', &
-     '<1 1 . 0>{-1 1 . 1} ', &
+     '< 1 1 . 0>{ 1 -1 . 1}', &
      '< 1 1 . 3>{-1  0 . 1}', &
      '< 1 1 . 3>{-1 -1 . 2}']
  
   real(pReal), dimension(4+4,LATTICE_HEX_NTWIN), parameter :: &
     LATTICE_HEX_SYSTEMTWIN =  reshape(real([&
      ! Compression or Tension = f(twinning shear=f(c/a)) for each metal ! (according to Yoo 1981)
-       1, -1,  0,  1,    -1,  1,  0,  2, & ! <-10.1>{10.2} shear = (3-(c/a)^2)/(sqrt(3) c/a)
+      -1,  0,  1,  1,     1,  0, -1,  2, & ! <-10.1>{10.2} shear = (3-(c/a)^2)/(sqrt(3) c/a)
-      -1,  0,  1,  1,     1,  0, -1,  2, &
+       0, -1,  1,  1,     0,  1, -1,  2, &
+       1, -1,  0,  1,    -1,  1,  0,  2, &
+       1,  0, -1,  1,    -1,  0,  1,  2, &
        0,  1, -1,  1,     0, -1,  1,  2, &
       -1,  1,  0,  1,     1, -1,  0,  2, &
-       1,  0, -1,  1,    -1,  0,  1,  2, &
-       0, -1,  1,  1,     0,  1, -1,  2, &
 !
-       2, -1, -1,  6,    -2,  1,  1,  1, & ! <11.6>{-1-1.1} shear = 1/(c/a)
+      -1, -1,  2,  6,     1,  1, -2,  1, & ! <11.6>{-1-1.1} shear = 1/(c/a)
-      -1,  2, -1,  6,     1, -2,  1,  1, &
-      -1, -1,  2,  6,     1,  1, -2,  1, &
-      -2,  1,  1,  6,     2, -1, -1,  1, &
        1, -2,  1,  6,    -1,  2, -1,  1, &
+       2, -1, -1,  6,    -2,  1,  1,  1, &
        1,  1, -2,  6,    -1, -1,  2,  1, &
+      -1,  2, -1,  6,     1, -2,  1,  1, &
+      -2,  1,  1,  6,     2, -1, -1,  1, &
 !
-      -1,  1,  0, -2,    -1,  1,  0,  1, & !! <10.-2>{10.1} shear = (4(c/a)^2-9)/(4 sqrt(3) c/a)
+       1,  0, -1, -2,     1,  0, -1,  1, & ! <10.-2>{10.1} shear = (4(c/a)^2-9)/(4 sqrt(3) c/a)
-       1,  0, -1, -2,     1,  0, -1,  1, &
+       0,  1, -1, -2,     0,  1, -1,  1, &
+      -1,  1,  0, -2,    -1,  1,  0,  1, &
+      -1,  0,  1, -2,    -1,  0,  1,  1, &
        0, -1,  1, -2,     0, -1,  1,  1, &
        1, -1,  0, -2,     1, -1,  0,  1, &
-      -1,  0,  1, -2,    -1,  0,  1,  1, &
-       0,  1, -1, -2,     0,  1, -1,  1, &
 !
-       2, -1, -1, -3,     2, -1, -1,  2, & ! <11.-3>{11.2} shear = 2((c/a)^2-2)/(3 c/a)
+       1,  1, -2, -3,     1,  1, -2,  2, & ! <11.-3>{11.2} shear = 2((c/a)^2-2)/(3 c/a)
       -1,  2, -1, -3,    -1,  2, -1,  2, &
-      -1, -1,  2, -3,    -1, -1,  2,  2, &
       -2,  1,  1, -3,    -2,  1,  1,  2, &
+      -1, -1,  2, -3,    -1, -1,  2,  2, &
        1, -2,  1, -3,     1, -2,  1,  2, &
-       1,  1, -2, -3,     1,  1, -2,  2  &
+       2, -1, -1, -3,     2, -1, -1,  2  &
-      ],pReal),shape(LATTICE_HEX_SYSTEMTWIN))                                                       !< twin systems for hex, order follows Prof. Tom Bieler's scheme
+      ],pReal),shape(LATTICE_HEX_SYSTEMTWIN))                                                       !< twin systems for hex, sorted by P. Eisenlohr CCW around <c> starting next to a_1 axis
  
   character(len=*), dimension(4), parameter :: LATTICE_HEX_TWINFAMILY_NAME = &
     ['<-1 0 .  1>{ 1  0 . 2}', &

src/math.f90

@@ -1307,10 +1307,10 @@ real(pReal) pure function math_volTetrahedron(v1,v2,v3,v4)
   real(pReal), dimension (3,3) :: m
 
   m(1:3,1) = v1-v2
-  m(1:3,2) = v2-v3
+  m(1:3,2) = v1-v3
-  m(1:3,3) = v3-v4
+  m(1:3,3) = v1-v4
 
-  math_volTetrahedron = math_det33(m)/6.0_pReal
+  math_volTetrahedron = abs(math_det33(m))/6.0_pReal
 
 end function math_volTetrahedron
 
@@ -1404,6 +1404,7 @@ subroutine unitTest
   integer,     dimension(5) :: range_out_ = [1,2,3,4,5]
 
   real(pReal)                 :: det
+  real(pReal), dimension(3)   :: v3_1,v3_2,v3_3,v3_4
   real(pReal), dimension(6)   :: v6
   real(pReal), dimension(9)   :: v9
   real(pReal), dimension(3,3) :: t33,t33_2
@@ -1452,6 +1453,15 @@ subroutine unitTest
   if(any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) &
     call IO_error(401,ext_msg='math_symmetric33')
  
+  call random_number(v3_1)
+  call random_number(v3_2)
+  call random_number(v3_3)
+  call random_number(v3_4)
+  
+  if(dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, &
+          math_volTetrahedron(v3_1,v3_2,v3_3,v3_4),tol=1.0e-12_pReal)) &
+  call IO_error(401,ext_msg='math_volTetrahedron')
+    
   call random_number(t33)
   if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) &
     call IO_error(401,ext_msg='math_det33/math_detSym33')

src/mesh_FEM.f90

@@ -28,8 +28,7 @@ module mesh
  integer, public, protected :: &
    mesh_Nboundaries, &
    mesh_NcpElems, &                                                                                 !< total number of CP elements in mesh
-   mesh_NcpElemsGlobal, &
+   mesh_NcpElemsGlobal
-   mesh_Nnodes                                                                                      !< total number of nodes in mesh
 
 !!!! BEGIN DEPRECATED !!!!!
  integer, public, protected :: &
@@ -69,8 +68,7 @@ module mesh
  public :: &
    mesh_init, &
    mesh_FEM_build_ipVolumes, &
-   mesh_FEM_build_ipCoordinates, &
+   mesh_FEM_build_ipCoordinates
-   mesh_cellCenterCoordinates
 
 contains
 
@@ -107,7 +105,8 @@ subroutine mesh_init
  integer, parameter :: FILEUNIT = 222
  integer :: j
  integer, allocatable, dimension(:) :: chunkPos
- integer :: dimPlex
+ integer :: dimPlex, &
+   mesh_Nnodes                                                                                      !< total number of nodes in mesh
  integer, parameter :: &
    mesh_ElemType=1                                                                             !< Element type of the mesh (only support homogeneous meshes)
  character(len=512) :: &
@@ -221,9 +220,6 @@ subroutine mesh_init
  call theMesh%init(dimplex,integrationOrder,mesh_node0)
  call theMesh%setNelems(mesh_NcpElems)
 
- theMesh%homogenizationAt  = mesh_element(3,:)
- theMesh%microstructureAt  = mesh_element(4,:)
- 
    call discretization_init(mesh_element(3,:),mesh_element(4,:),&
                            reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), &
                            mesh_node0)
@@ -231,26 +227,8 @@ subroutine mesh_init
 end subroutine mesh_init
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief Calculates cell center coordinates.
-!--------------------------------------------------------------------------------------------------
-pure function mesh_cellCenterCoordinates(ip,el)
- 
- integer, intent(in) :: el, &                                                                  !< element number
-                        ip                                                                     !< integration point number
- real(pReal), dimension(3) :: mesh_cellCenterCoordinates                                             !< x,y,z coordinates of the cell center of the requested IP cell
-
-end function mesh_cellCenterCoordinates
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume'
-!> @details The IP volume is calculated differently depending on the cell type.
-!> 2D cells assume an element depth of one in order to calculate the volume.
-!> For the hexahedral cell we subdivide the cell into subvolumes of pyramidal
-!> shape with a cell face as basis and the central ip at the tip. This subvolume is
-!> calculated as an average of four tetrahedals with three corners on the cell face 
-!> and one corner at the central ip.
 !--------------------------------------------------------------------------------------------------
 subroutine mesh_FEM_build_ipVolumes(dimPlex)
  

src/mesh_base.f90

@@ -1,4 +1,3 @@
-
 !--------------------------------------------------------------------------------------------------
 !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
 !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
@@ -8,14 +7,13 @@
 !--------------------------------------------------------------------------------------------------
 module mesh_base
 
-  use, intrinsic :: iso_c_binding
   use prec
   use element
 
   implicit none
 
 !---------------------------------------------------------------------------------------------------
-!> Properties of a the whole mesh (consisting of one type of elements)
+!> Properties of a whole mesh (consisting of one type of elements)
 !---------------------------------------------------------------------------------------------------
  type, public :: tMesh
    type(tElement) :: &
@@ -33,11 +31,7 @@ module mesh_base
      elemType, &
      Ncells, &
      nIPneighbors, &
-     NcellNodes, &
+     NcellNodes
-     maxElemsPerNode
-   integer(pInt), dimension(:), allocatable, public :: &
-     homogenizationAt, &
-     microstructureAt
    integer(pInt), dimension(:,:), allocatable, public :: &
      connectivity
    contains
@@ -47,6 +41,7 @@ module mesh_base
  end type tMesh
 
 contains
+
 subroutine tMesh_base_init(self,meshType,elemType,nodes)
  
  class(tMesh) :: self

src/plastic_isotropic.f90

@@ -434,7 +434,7 @@ end function plastic_isotropic_postResults
 !> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_isotropic_results(instance,group)
-#if defined(PETSc) || defined(DAMASKHDF5)
+#if defined(PETSc) || defined(DAMASK_HDF5)
 
   integer, intent(in) :: instance
   character(len=*), intent(in) :: group

src/plastic_nonlocal.f90

@@ -1452,9 +1452,9 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
     opposite_n, &                                                                                   !< neighbor index pointing to me when looking from my opposite neighbor
     t, &                                                                                            !< type of dislocation
     o,&                                                                                             !< offset shortcut
-    no,&                                                                                            !< neighbour offset shortcut
+    no,&                                                                                            !< neighbor offset shortcut
     p,&                                                                                             !< phase shortcut
-    np,&                                                                                            !< neighbour phase shortcut
+    np,&                                                                                            !< neighbor phase shortcut
     topp, &                                                                                         !< type of dislocation with opposite sign to t
     s                                                                                               !< index of my current slip system
   real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,el))),10) :: &

src/results.f90

@@ -65,7 +65,7 @@ subroutine results_init
   write(6,'(/,a)') ' <<<+-  results init  -+>>>'
 
   write(6,'(/,a)') ' Diehl et al., Integrating Materials and Manufacturing Innovation 6(1):83–91, 2017'
-  write(6,'(a)')   ' https://doi.org/10.1007/s40192-018-0118-7'
+  write(6,'(a)')   ' https://doi.org/10.1007/s40192-017-0084-5'
 
   resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','w',.true.)
   call HDF5_addAttribute(resultsFile,'DADF5-version',0.3_pReal)
@@ -296,21 +296,34 @@ end subroutine results_writeVectorDataset_real
 !--------------------------------------------------------------------------------------------------
 !> @brief stores a tensor dataset in a group
 !--------------------------------------------------------------------------------------------------
-subroutine results_writeTensorDataset_real(group,dataset,label,description,SIunit)
+subroutine results_writeTensorDataset_real(group,dataset,label,description,SIunit,transposed)
 
   character(len=*), intent(in)                   :: label,group,description
   character(len=*), intent(in), optional         :: SIunit
+  logical,          intent(in), optional         :: transposed
   real(pReal),      intent(in), dimension(:,:,:) :: dataset
   
   integer :: i 
+  logical :: T
   integer(HID_T) :: groupHandle
   real(pReal), dimension(:,:,:), allocatable :: dataset_transposed
 
   
+  if(present(transposed)) then
+    T = transposed
+  else
+    T = .true.
+  endif
+
+  if(T) then
+    if(size(dataset,1) /= size(dataset,2)) call IO_error(0,ext_msg='transpose non-symmetric tensor')
     allocate(dataset_transposed,mold=dataset)
-  do i=1,size(dataset,3)
+    do i=1,size(dataset_transposed,3)
-    dataset_transposed(1:3,1:3,i) = transpose(dataset(1:3,1:3,i))
+      dataset_transposed(:,:,i) = transpose(dataset(:,:,i))
     enddo
+  else
+    allocate(dataset_transposed,source=dataset)
+  endif
 
   groupHandle = results_openGroup(group)