Merge branch 'Marc-No-nonlocal' into 'development'

Marc no nonlocal

See merge request damask/DAMASK!98

.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 ----------------------------------------------

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
@@ -40,21 +41,30 @@ if options.con is None: options.con=[]
 for filename in options.filenames:
   results = damask.DADF5(filename)
   
-  if results.structured:                                                                               # for grid solvers use rectilinear grid
+  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()
 
 
     dirname  = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir))
@@ -122,7 +132,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()

python/damask/dadf5.py

@@ -348,6 +348,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,11 +313,10 @@ 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 explicit:     11, 22, 33, 12, 23, 13
+!     ABAQUS implicit:   11, 22, 33, 12, 13, 23
-!     ABAQUS implicit:     11, 22, 33, 12, 13, 23
+!     ABAQUS implicit:   11, 22, 33, 12
-!     ABAQUS implicit:     11, 22, 33, 12
 
  ddsdde = ddsdde_h(1:ntens,1:ntens)
  stress = stress_h(1:ntens)

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

@@ -1072,8 +1072,8 @@ integer function IO_countContinuousIntValues(fileUnit)
    if (chunkPos(1) < 1) then                                                                        ! empty line
      exit
    elseif (IO_lc(IO_stringValue(line,chunkPos,2)) == 'to' ) then                                    ! found range indicator
-     IO_countContinuousIntValues = 1 + abs(  IO_intValue(line,chunkPos,3) &
+     IO_countContinuousIntValues = 1 + abs( IO_intValue(line,chunkPos,3) &
-                                                - IO_intValue(line,chunkPos,1))
+                                           -IO_intValue(line,chunkPos,1))
      exit                                                                                           ! only one single range indicator allowed                              
    else
      IO_countContinuousIntValues = IO_countContinuousIntValues+chunkPos(1)-1                        ! add line's count when assuming 'c'

src/mesh_FEM.f90

@@ -221,9 +221,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)

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