Merge branch 'development' of magit1.mpie.de:damask/DAMASK into development

.gitattributes

@@ -6,3 +6,4 @@
 # Denote all files that are truly binary and should not be modified.
 *.png binary
 *.jpg binary
+*.cae binary

code/IO.f90

@@ -6,10 +6,6 @@
 !> @brief  input/output functions, partly depending on chosen solver
 !--------------------------------------------------------------------------------------------------
 module IO
-#ifdef HDF
- use hdf5, only: &
-   HID_T
-#endif
  use prec, only: &
    pInt, &
    pReal
@@ -18,22 +14,8 @@ module IO
  private
  character(len=5), parameter, public :: &
    IO_EOF = '#EOF#'                                                                                 !< end of file string
-#ifdef HDF 
- integer(HID_T), public, protected :: tempCoordinates, tempResults
- integer(HID_T), private :: resultsFile, tempFile
- integer(pInt), private :: currentInc
-#endif
 
  public :: &
-#ifdef HDF
-   HDF5_mappingConstitutive, &
-   HDF5_mappingHomogenization, &
-   HDF5_mappingCells, &
-   HDF5_addGroup ,&
-   HDF5_forwardResults, &
-   HDF5_addScalarDataset, &
-   IO_formatIntToString ,&
-#endif
    IO_init, &
    IO_read, &
    IO_checkAndRewind, &
@@ -117,9 +99,6 @@ subroutine IO_init
 #include "compilation_info.f90"
  endif mainProcess
 
-#ifdef HDF 
- call HDF5_createJobFile
-#endif
 
 end subroutine IO_init
 
@@ -1944,526 +1923,4 @@ recursive function abaqus_assembleInputFile(unit1,unit2) result(createSuccess)
 end function abaqus_assembleInputFile
 #endif
 
-
-#ifdef HDF
-!--------------------------------------------------------------------------------------------------
-!> @brief creates and initializes HDF5 output files
-!--------------------------------------------------------------------------------------------------
-subroutine HDF5_createJobFile
- use hdf5
- use DAMASK_interface, only: &
-   getSolverWorkingDirectoryName, &
-   getSolverJobName
-
- implicit none
- integer              :: hdferr
- integer(SIZE_T)      :: typeSize 
- character(len=1024)  :: path
- integer(HID_T) :: prp_id
- integer(SIZE_T), parameter :: increment = 104857600 ! increase temp file in memory in 100MB steps
-
-
-!--------------------------------------------------------------------------------------------------
-! initialize HDF5 library and check if integer and float type size match
- call h5open_f(hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5open_f')
- call h5tget_size_f(H5T_NATIVE_INTEGER,typeSize, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5tget_size_f (int)')
- if (int(pInt,SIZE_T)/=typeSize) call IO_error(0_pInt,ext_msg='pInt does not match H5T_NATIVE_INTEGER')
- call h5tget_size_f(H5T_NATIVE_DOUBLE,typeSize, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5tget_size_f (double)')
- if (int(pReal,SIZE_T)/=typeSize) call IO_error(0_pInt,ext_msg='pReal does not match H5T_NATIVE_DOUBLE')
-
-!--------------------------------------------------------------------------------------------------
-! open file
- path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//'DAMASKout'
- call h5fcreate_f(path,H5F_ACC_TRUNC_F,resultsFile,hdferr)
- if (hdferr < 0) call IO_error(100_pInt,ext_msg=path)
- call HDF5_addStringAttribute(resultsFile,'createdBy','$Id$')
-
-!--------------------------------------------------------------------------------------------------
-! open temp file
- path = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.'//'DAMASKoutTemp'
- call h5pcreate_f(H5P_FILE_ACCESS_F, prp_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5pcreate_f')
- call h5pset_fapl_core_f(prp_id, increment, .false., hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_createJobFile: h5pset_fapl_core_f')
- call h5fcreate_f(path,H5F_ACC_TRUNC_F,tempFile,hdferr)
- if (hdferr < 0) call IO_error(100_pInt,ext_msg=path)
-
-!--------------------------------------------------------------------------------------------------
-! create mapping groups in out file
- call HDF5_closeGroup(HDF5_addGroup("mapping"))
- call HDF5_closeGroup(HDF5_addGroup("results"))
- call HDF5_closeGroup(HDF5_addGroup("coordinates"))
-
-!--------------------------------------------------------------------------------------------------
-! create results group in temp file
- tempResults     = HDF5_addGroup("results",tempFile)
- tempCoordinates = HDF5_addGroup("coordinates",tempFile)
-
-end subroutine HDF5_createJobFile
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief creates and initializes HDF5 output file
-!--------------------------------------------------------------------------------------------------
-subroutine HDF5_closeJobFile()
- use hdf5
- 
- implicit none
- integer :: hdferr
- call h5fclose_f(resultsFile,hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_closeJobFile: h5fclose_f')
-
-end subroutine HDF5_closeJobFile
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief adds a new group to the results file, or if loc is present at the given location
-!--------------------------------------------------------------------------------------------------
-integer(HID_T) function HDF5_addGroup(path,loc)
- use hdf5
-
- implicit none
- character(len=*), intent(in) :: path
- integer(HID_T), intent(in),optional :: loc
- integer                      :: hdferr
-
- if (present(loc)) then
-   call h5gcreate_f(loc, trim(path), HDF5_addGroup, hdferr)
- else
-   call h5gcreate_f(resultsFile, trim(path), HDF5_addGroup, hdferr)
- endif
- if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_addGroup: h5gcreate_f ('//trim(path)//' )')
-
-end function HDF5_addGroup
-
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief adds a new group to the results file
-!--------------------------------------------------------------------------------------------------
-integer(HID_T) function HDF5_openGroup(path)
- use hdf5
-
- implicit none
- character(len=*), intent(in) :: path
- integer                      :: hdferr
-
- call h5gopen_f(resultsFile, trim(path), HDF5_openGroup, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_openGroup: h5gopen_f ('//trim(path)//' )')
-
-end function HDF5_openGroup
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief closes a group
-!--------------------------------------------------------------------------------------------------
-subroutine HDF5_closeGroup(ID)
- use hdf5
-
- implicit none
- integer(HID_T), intent(in) :: ID
- integer                    :: hdferr
-
- call h5gclose_f(ID, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_closeGroup: h5gclose_f')
-
-end subroutine HDF5_closeGroup
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief adds a new group to the results file
-!--------------------------------------------------------------------------------------------------
-subroutine HDF5_addStringAttribute(entity,attrLabel,attrValue)
- use hdf5
-
- implicit none
- integer(HID_T),   intent(in)  :: entity
- character(len=*), intent(in)  :: attrLabel, attrValue
- integer                       :: hdferr
- integer(HID_T)                :: attr_id, space_id, type_id
-
- call h5screate_f(H5S_SCALAR_F,space_id,hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5screate_f')
- call h5tcopy_f(H5T_NATIVE_CHARACTER, type_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5tcopy_f')
- call h5tset_size_f(type_id, int(len(trim(attrValue)),HSIZE_T), hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5tset_size_f')
- call h5acreate_f(entity, trim(attrLabel),type_id,space_id,attr_id,hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5acreate_f')
- call h5awrite_f(attr_id, type_id, trim(attrValue), int([1],HSIZE_T), hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5awrite_f')
- call h5aclose_f(attr_id,hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5aclose_f')
- call h5sclose_f(space_id,hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_addStringAttribute: h5sclose_f')
-
-end subroutine HDF5_addStringAttribute
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief adds the unique mapping from spatial position and constituent ID to results
-!--------------------------------------------------------------------------------------------------
-subroutine HDF5_mappingConstitutive(mapping)
- use hdf5
-
- implicit none
- integer(pInt), intent(in), dimension(:,:,:) :: mapping
-
- integer        :: hdferr, NmatPoints,Nconstituents
- integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id,instance_id,position_id
-
- Nconstituents=size(mapping,1)
- NmatPoints=size(mapping,2)
- mapping_ID = HDF5_openGroup("mapping")
-
-!--------------------------------------------------------------------------------------------------
-! create dataspace
- call h5screate_simple_f(2, int([Nconstituents,NmatPoints],HSIZE_T), space_id, hdferr, &
-                            int([Nconstituents,NmatPoints],HSIZE_T))
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive')
-
-!--------------------------------------------------------------------------------------------------
-! compound type
- call h5tcreate_f(H5T_COMPOUND_F, 6_SIZE_T, dtype_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tcreate_f dtype_id')
-
- call h5tinsert_f(dtype_id, "Constitutive Instance",         0_SIZE_T, H5T_STD_U16LE, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tinsert_f 0')
- call h5tinsert_f(dtype_id, "Position in Instance Results",  2_SIZE_T, H5T_STD_U32LE, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tinsert_f 2')
-
-!--------------------------------------------------------------------------------------------------
-! create Dataset
- call h5dcreate_f(mapping_id, "Constitutive", dtype_id, space_id, dset_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive')
-
-!--------------------------------------------------------------------------------------------------
-! Create memory types (one compound datatype for each member)
- call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), instance_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tcreate_f instance_id')
- call h5tinsert_f(instance_id, "Constitutive Instance",        0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tinsert_f instance_id')
-
- call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tcreate_f position_id')
- call h5tinsert_f(position_id, "Position in Instance Results", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tinsert_f position_id')
-
-!--------------------------------------------------------------------------------------------------
-! write data by fields in the datatype. Fields order is not important.
- call h5dwrite_f(dset_id, position_id, mapping(1:Nconstituents,1:NmatPoints,1), &
-                                            int([Nconstituents,  NmatPoints],HSIZE_T), hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5dwrite_f position_id')
- 
- call h5dwrite_f(dset_id, instance_id, mapping(1:Nconstituents,1:NmatPoints,2), &
-                                            int([Nconstituents,  NmatPoints],HSIZE_T), hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5dwrite_f instance_id')
-
-!--------------------------------------------------------------------------------------------------
-!close types, dataspaces
- call h5tclose_f(dtype_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tclose_f dtype_id')
- call h5tclose_f(position_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tclose_f position_id')
- call h5tclose_f(instance_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5tclose_f instance_id')
- call h5dclose_f(dset_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5dclose_f')
- call h5sclose_f(space_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5sclose_f')
- call HDF5_closeGroup(mapping_ID)
-
-end subroutine HDF5_mappingConstitutive
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief adds the unique mapping from spatial position and constituent ID to results
-!--------------------------------------------------------------------------------------------------
-subroutine HDF5_mappingCrystallite(mapping)
- use hdf5
-
- implicit none
- integer(pInt), intent(in), dimension(:,:,:) :: mapping
-
- integer        :: hdferr, NmatPoints,Nconstituents
- integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id,instance_id,position_id
-
- Nconstituents=size(mapping,1)
- NmatPoints=size(mapping,2)
- mapping_ID = HDF5_openGroup("mapping")
-
-!--------------------------------------------------------------------------------------------------
-! create dataspace
- call h5screate_simple_f(2, int([Nconstituents,NmatPoints],HSIZE_T), space_id, hdferr, &
-                            int([Nconstituents,NmatPoints],HSIZE_T))
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite')
-
-!--------------------------------------------------------------------------------------------------
-! compound type
- call h5tcreate_f(H5T_COMPOUND_F, 6_SIZE_T, dtype_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tcreate_f dtype_id')
-
- call h5tinsert_f(dtype_id, "Crystallite Instance",          0_SIZE_T, H5T_STD_U16LE, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f 0')
- call h5tinsert_f(dtype_id, "Position in Instance Results",  2_SIZE_T, H5T_STD_U32LE, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f 2')
-
-!--------------------------------------------------------------------------------------------------
-! create Dataset
- call h5dcreate_f(mapping_id, "Crystallite", dtype_id, space_id, dset_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite')
-
-!--------------------------------------------------------------------------------------------------
-! Create memory types (one compound datatype for each member)
- call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), instance_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tcreate_f instance_id')
- call h5tinsert_f(instance_id, "Crystallite Instance",        0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f instance_id')
-
- call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tcreate_f position_id')
- call h5tinsert_f(position_id, "Position in Instance Results", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tinsert_f position_id')
-
-!--------------------------------------------------------------------------------------------------
-! write data by fields in the datatype. Fields order is not important.
- call h5dwrite_f(dset_id, position_id, mapping(1:Nconstituents,1:NmatPoints,1), &
-                                            int([Nconstituents,  NmatPoints],HSIZE_T), hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dwrite_f position_id')
- 
- call h5dwrite_f(dset_id, instance_id, mapping(1:Nconstituents,1:NmatPoints,2), &
-                                            int([Nconstituents,  NmatPoints],HSIZE_T), hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dwrite_f instance_id')
-
-!--------------------------------------------------------------------------------------------------
-!close types, dataspaces
- call h5tclose_f(dtype_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f dtype_id')
- call h5tclose_f(position_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f position_id')
- call h5tclose_f(instance_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5tclose_f instance_id')
- call h5dclose_f(dset_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5dclose_f')
- call h5sclose_f(space_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCrystallite: h5sclose_f')
- call HDF5_closeGroup(mapping_ID)
-
-end subroutine HDF5_mappingCrystallite
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief adds the unique mapping from spatial position to results
-!--------------------------------------------------------------------------------------------------
-subroutine HDF5_mappingHomogenization(mapping)
- use hdf5
-
- implicit none
- integer(pInt), intent(in), dimension(:,:) :: mapping
-
- integer        :: hdferr, NmatPoints
- integer(HID_T) :: mapping_id, dtype_id, dset_id, space_id,instance_id,position_id,elem_id,ip_id
-
- NmatPoints=size(mapping,1)
- mapping_ID = HDF5_openGroup("mapping")
-
-!--------------------------------------------------------------------------------------------------
-! create dataspace
- call h5screate_simple_f(1, int([NmatPoints],HSIZE_T), space_id, hdferr, &
-                            int([NmatPoints],HSIZE_T))
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization')
-
-!--------------------------------------------------------------------------------------------------
-! compound type
- call h5tcreate_f(H5T_COMPOUND_F, 11_SIZE_T, dtype_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tcreate_f dtype_id')
- 
- call h5tinsert_f(dtype_id, "Homogenization Instance",      0_SIZE_T, H5T_STD_U16LE, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f 0')
- call h5tinsert_f(dtype_id, "Position in Instance Results", 2_SIZE_T, H5T_STD_U32LE, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f 2')
- call h5tinsert_f(dtype_id, "Element Number",               6_SIZE_T, H5T_STD_U32LE, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f 6')
- call h5tinsert_f(dtype_id, "Material Point Number",       10_SIZE_T, H5T_STD_U8LE, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f 10')
-
-!--------------------------------------------------------------------------------------------------
-! create Dataset
- call h5dcreate_f(mapping_id, "Homogenization", dtype_id, space_id, dset_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization')
-
-!--------------------------------------------------------------------------------------------------
-! Create memory types (one compound datatype for each member)
- call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), instance_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tcreate_f instance_id')
- call h5tinsert_f(instance_id, "Homogenization Instance",      0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f instance_id')
-
- call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), position_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tcreate_f position_id')
- call h5tinsert_f(position_id, "Position in Instance Results", 0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f position_id')
-
- call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), elem_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tcreate_f elem_id')
- call h5tinsert_f(elem_id,     "Element Number",               0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f elem_id')
-
- call h5tcreate_f(H5T_COMPOUND_F, int(pInt,SIZE_T), ip_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tcreate_f ip_id')
- call h5tinsert_f(ip_id,       "Material Point Number",        0_SIZE_T, H5T_NATIVE_INTEGER, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tinsert_f ip_id')
-
-!--------------------------------------------------------------------------------------------------
-! write data by fields in the datatype. Fields order is not important.
- call h5dwrite_f(dset_id, position_id, mapping(1:NmatPoints,1), &
-                                            int([NmatPoints],HSIZE_T), hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5dwrite_f position_id')
- 
- call h5dwrite_f(dset_id, instance_id, mapping(1:NmatPoints,2), &
-                                            int([NmatPoints],HSIZE_T), hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5dwrite_f position_id')
-
- call h5dwrite_f(dset_id, elem_id,     mapping(1:NmatPoints,3), &
-                                            int([NmatPoints],HSIZE_T), hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5dwrite_f elem_id')
-
- call h5dwrite_f(dset_id, ip_id,       mapping(1:NmatPoints,4), &
-                                            int([NmatPoints],HSIZE_T), hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5dwrite_f ip_id')
-
-!--------------------------------------------------------------------------------------------------
-!close types, dataspaces
- call h5tclose_f(dtype_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tclose_f dtype_id')
- call h5tclose_f(position_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tclose_f position_id')
- call h5tclose_f(instance_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tclose_f instance_id')
- call h5tclose_f(ip_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tclose_f ip_id')
- call h5tclose_f(elem_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5tclose_f elem_id')
- call h5dclose_f(dset_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5dclose_f')
- call h5sclose_f(space_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingHomogenization: h5sclose_f')
- call HDF5_closeGroup(mapping_ID)
-
-end subroutine HDF5_mappingHomogenization
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief adds the unique cell to node mapping
-!--------------------------------------------------------------------------------------------------
-subroutine HDF5_mappingCells(mapping)
- use hdf5
-
- implicit none
- integer(pInt), intent(in), dimension(:) :: mapping
-
- integer        :: hdferr, Nnodes
- integer(HID_T) :: mapping_id, dset_id, space_id
-
- Nnodes=size(mapping)
- mapping_ID = HDF5_openGroup("mapping")
-
-!--------------------------------------------------------------------------------------------------
-! create dataspace
- call h5screate_simple_f(1, int([Nnodes],HSIZE_T), space_id, hdferr, &
-                            int([Nnodes],HSIZE_T))
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCells: h5screate_simple_f')
-
-!--------------------------------------------------------------------------------------------------
-! create Dataset
- call h5dcreate_f(mapping_id, "Cell",H5T_NATIVE_INTEGER, space_id, dset_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCells')
-
-!--------------------------------------------------------------------------------------------------
-! write data by fields in the datatype. Fields order is not important.
- call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER, mapping, int([Nnodes],HSIZE_T), hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingCells: h5dwrite_f instance_id')
-
-!--------------------------------------------------------------------------------------------------
-!close types, dataspaces
- call h5dclose_f(dset_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5dclose_f')
- call h5sclose_f(space_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='IO_mappingConstitutive: h5sclose_f')
- call HDF5_closeGroup(mapping_ID)
-
-end subroutine HDF5_mappingCells
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief creates a new scalar dataset in the given group location
-!--------------------------------------------------------------------------------------------------
-subroutine HDF5_addScalarDataset(group,nnodes,label,SIunit)
- use hdf5
-
- implicit none
- integer(HID_T),   intent(in) :: group
- integer(pInt),    intent(in) :: nnodes
- character(len=*), intent(in)  :: SIunit,label
-
- integer        :: hdferr
- integer(HID_T) :: dset_id, space_id
-
-!--------------------------------------------------------------------------------------------------
-! create dataspace
- call h5screate_simple_f(1, int([Nnodes],HSIZE_T), space_id, hdferr, &
-                            int([Nnodes],HSIZE_T))
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addScalarDataset: h5screate_simple_f')
-
-!--------------------------------------------------------------------------------------------------
-! create Dataset
- call h5dcreate_f(group, trim(label),H5T_NATIVE_DOUBLE, space_id, dset_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addScalarDataset: h5dcreate_f')
- call HDF5_addStringAttribute(dset_id,'unit',trim(SIunit))
-
-!--------------------------------------------------------------------------------------------------
-!close types, dataspaces
- call h5dclose_f(dset_id, hdferr)
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addScalarDataset: h5dclose_f')
- call h5sclose_f(space_id, hdferr)
-
-end subroutine HDF5_addScalarDataset
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief returns nicely formatted string of integer value
-!--------------------------------------------------------------------------------------------------
-function IO_formatIntToString(myInt)
-
- implicit none
- integer(pInt), intent(in) :: myInt
- character(len=1_pInt + int(log10(real(myInt)),pInt)) :: IO_formatIntToString
- write(IO_formatIntToString,'('//IO_intOut(myInt)//')') myInt
- 
- end function
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief copies the current temp results to the actual results file
-!--------------------------------------------------------------------------------------------------
-subroutine HDF5_forwardResults
- use hdf5
-
- implicit none
- integer        :: hdferr
- integer(HID_T)   :: new_loc_id
-    
- new_loc_id = HDF5_openGroup("results")
- currentInc = currentInc + 1_pInt
- call h5ocopy_f(tempFile, 'results', new_loc_id,dst_name=IO_formatIntToString(currentInc), hdferr=hdferr)
- if (hdferr < 0_pInt) call IO_error(1_pInt,ext_msg='HDF5_forwardResults: h5ocopy_f')
- call HDF5_closeGroup(new_loc_id)
-
-end subroutine HDF5_forwardResults
-
-
-#endif
 end module IO

code/Makefile

@@ -309,7 +309,7 @@ KINEMATICS_FILES = \
    kinematics_vacancy_strain.o kinematics_hydrogen_strain.o
 
 PLASTIC_FILES = \
-   plastic_dislotwin.o plastic_disloUCLA.o plastic_isotropic.o plastic_j2.o \
+   plastic_dislotwin.o plastic_disloUCLA.o plastic_isotropic.o \
    plastic_phenopowerlaw.o plastic_titanmod.o plastic_nonlocal.o plastic_none.o \
    plastic_phenoplus.o
 
@@ -579,9 +579,6 @@ plastic_phenoplus.o:           plastic_phenoplus.f90 \
 plastic_isotropic.o:           plastic_isotropic.f90 \
                                lattice.o
 
-plastic_j2.o:                  plastic_j2.f90 \
-                               lattice.o
-
 plastic_none.o:                plastic_none.f90 \
                                lattice.o
 ifeq "$(F90)" "gfortran"

code/commercialFEM_fileList.f90

@@ -28,7 +28,6 @@
 #include "kinematics_hydrogen_strain.f90"
 #include "plastic_none.f90"
 #include "plastic_isotropic.f90"
-#include "plastic_j2.f90"
 #include "plastic_phenopowerlaw.f90"
 #include "plastic_phenoplus.f90"
 #include "plastic_titanmod.f90"

code/constitutive.f90

@@ -69,7 +69,6 @@ subroutine constitutive_init()
    ELASTICITY_hooke_ID, &
    PLASTICITY_none_ID, &
    PLASTICITY_isotropic_ID, &
-   PLASTICITY_j2_ID, &
    PLASTICITY_phenopowerlaw_ID, &
    PLASTICITY_phenoplus_ID, &
    PLASTICITY_dislotwin_ID, &
@@ -93,7 +92,6 @@ subroutine constitutive_init()
    ELASTICITY_HOOKE_label, &
    PLASTICITY_NONE_label, &
    PLASTICITY_ISOTROPIC_label, &
-   PLASTICITY_J2_label, &
    PLASTICITY_PHENOPOWERLAW_label, &
    PLASTICITY_PHENOPLUS_label, &
    PLASTICITY_DISLOTWIN_label, &
@@ -114,7 +112,6 @@ subroutine constitutive_init()
 
  use plastic_none
  use plastic_isotropic
- use plastic_j2
  use plastic_phenopowerlaw
  use plastic_phenoplus
  use plastic_dislotwin
@@ -160,7 +157,6 @@ subroutine constitutive_init()
 ! parse plasticities from config file
  if (any(phase_plasticity == PLASTICITY_NONE_ID))          call plastic_none_init
  if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID))     call plastic_isotropic_init(FILEUNIT)
- if (any(phase_plasticity == PLASTICITY_J2_ID))            call plastic_j2_init(FILEUNIT)
  if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init(FILEUNIT)
  if (any(phase_plasticity == PLASTICITY_PHENOPLUS_ID))     call plastic_phenoplus_init(FILEUNIT)
  if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID))     call plastic_dislotwin_init(FILEUNIT)
@@ -217,11 +213,6 @@ subroutine constitutive_init()
            thisNoutput => plastic_isotropic_Noutput
            thisOutput => plastic_isotropic_output
            thisSize   => plastic_isotropic_sizePostResult
-         case (PLASTICITY_J2_ID) plasticityType
-           outputName = PLASTICITY_J2_label
-           thisNoutput => plastic_j2_Noutput
-           thisOutput => plastic_j2_output
-           thisSize   => plastic_j2_sizePostResult
          case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
            outputName = PLASTICITY_PHENOPOWERLAW_label
            thisNoutput => plastic_phenopowerlaw_Noutput
@@ -513,7 +504,6 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
    thermalMapping, &
    PLASTICITY_NONE_ID, &
    PLASTICITY_ISOTROPIC_ID, &
-   PLASTICITY_J2_ID, &
    PLASTICITY_PHENOPOWERLAW_ID, &
    PLASTICITY_PHENOPLUS_ID, &
    PLASTICITY_DISLOTWIN_ID, &
@@ -522,8 +512,6 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
    PLASTICITY_NONLOCAL_ID
  use plastic_isotropic, only: &
    plastic_isotropic_LpAndItsTangent
- use plastic_j2, only: &
-   plastic_j2_LpAndItsTangent
  use plastic_phenopowerlaw, only: &
    plastic_phenopowerlaw_LpAndItsTangent
  use plastic_phenoplus, only: &
@@ -574,8 +562,6 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
      dLp_dMstar = 0.0_pReal
    case (PLASTICITY_ISOTROPIC_ID) plasticityType
      call plastic_isotropic_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
-   case (PLASTICITY_J2_ID) plasticityType
-     call plastic_j2_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
    case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
      call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
    case (PLASTICITY_PHENOPLUS_ID) plasticityType
@@ -903,7 +889,6 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
    homogenization_maxNgrains, &
    PLASTICITY_none_ID, &
    PLASTICITY_isotropic_ID, &
-   PLASTICITY_j2_ID, &
    PLASTICITY_phenopowerlaw_ID, &
    PLASTICITY_phenoplus_ID, &
    PLASTICITY_dislotwin_ID, &
@@ -916,8 +901,6 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
    SOURCE_thermal_externalheat_ID
  use plastic_isotropic, only:  &
    plastic_isotropic_dotState
- use plastic_j2, only:  &
-   plastic_j2_dotState
  use plastic_phenopowerlaw, only: &
    plastic_phenopowerlaw_dotState
  use plastic_phenoplus, only: &
@@ -971,8 +954,6 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
  plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
    case (PLASTICITY_ISOTROPIC_ID) plasticityType
      call plastic_isotropic_dotState    (Tstar_v,ipc,ip,el)
-   case (PLASTICITY_J2_ID) plasticityType
-     call plastic_j2_dotState           (Tstar_v,ipc,ip,el)
    case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
      call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
    case (PLASTICITY_PHENOPLUS_ID) plasticityType
@@ -1117,7 +1098,6 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
    homogenization_maxNgrains, &
    PLASTICITY_NONE_ID, &
    PLASTICITY_ISOTROPIC_ID, &
-   PLASTICITY_J2_ID, &
    PLASTICITY_PHENOPOWERLAW_ID, &
    PLASTICITY_PHENOPLUS_ID, &
    PLASTICITY_DISLOTWIN_ID, &
@@ -1130,8 +1110,6 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
    SOURCE_damage_anisoDuctile_ID
  use plastic_isotropic, only: &
    plastic_isotropic_postResults
- use plastic_j2, only: &
-   plastic_j2_postResults
  use plastic_phenopowerlaw, only: &
    plastic_phenopowerlaw_postResults
  use plastic_phenoplus, only: &
@@ -1185,8 +1163,6 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
      constitutive_postResults(startPos:endPos) = plastic_titanmod_postResults(ipc,ip,el)
    case (PLASTICITY_ISOTROPIC_ID) plasticityType
      constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
-   case (PLASTICITY_J2_ID) plasticityType
-     constitutive_postResults(startPos:endPos) = plastic_j2_postResults(Tstar_v,ipc,ip,el)
    case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
      constitutive_postResults(startPos:endPos) = &
        plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)

code/homogenization.f90

@@ -71,12 +71,6 @@ contains
 !> @brief module initialization
 !--------------------------------------------------------------------------------------------------
 subroutine homogenization_init
-#ifdef HDF
- use hdf5, only: &
-   HID_T
- use IO, only : &
-   HDF5_mappingHomogenization
-#endif
  use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
  use math, only: &
    math_I3
@@ -131,12 +125,6 @@ subroutine homogenization_init
  character(len=64), dimension(:,:), pointer :: thisOutput
  character(len=32) :: outputName                                                                    !< name of output, intermediate fix until HDF5 output is ready
  logical :: knownHomogenization, knownThermal, knownDamage, knownVacancyflux, knownPorosity, knownHydrogenflux
-#ifdef HDF
- integer(pInt), dimension(:,:), allocatable :: mapping
- integer(pInt), dimension(:), allocatable :: InstancePosition
- allocate(mapping(mesh_ncpelems,4),source=0_pInt)
- allocate(InstancePosition(material_Nhomogenization),source=0_pInt)
-#endif
 
 
 !--------------------------------------------------------------------------------------------------
@@ -396,17 +384,6 @@ subroutine homogenization_init
 
 !--------------------------------------------------------------------------------------------------
 ! allocate and initialize global state and postresutls variables
-#ifdef HDF
- elementLooping: do e = 1,mesh_NcpElems
-   myInstance = homogenization_typeInstance(mesh_element(3,e))
-   IpLooping: do i = 1,FE_Nips(FE_geomtype(mesh_element(2,e)))
-       InstancePosition(myInstance) = InstancePosition(myInstance)+1_pInt
-       mapping(e,1:4) = [instancePosition(myinstance),myinstance,e,i]
-   enddo IpLooping
- enddo elementLooping
- call  HDF5_mappingHomogenization(mapping)
-#endif
-
  homogenization_maxSizePostResults = 0_pInt
  thermal_maxSizePostResults        = 0_pInt
  damage_maxSizePostResults         = 0_pInt

code/material.f90

@@ -24,7 +24,6 @@ module material
    ELASTICITY_hooke_label               = 'hooke', &
    PLASTICITY_none_label                = 'none', &
    PLASTICITY_isotropic_label           = 'isotropic', &
-   PLASTICITY_j2_label                  = 'j2', &
    PLASTICITY_phenopowerlaw_label       = 'phenopowerlaw', &
    PLASTICITY_phenoplus_label           = 'phenoplus', &
    PLASTICITY_dislotwin_label           = 'dislotwin', &
@@ -74,7 +73,6 @@ module material
    enumerator :: PLASTICITY_undefined_ID, &
                  PLASTICITY_none_ID, &
                  PLASTICITY_isotropic_ID, &
-                 PLASTICITY_j2_ID, &
                  PLASTICITY_phenopowerlaw_ID, &
                  PLASTICITY_phenoplus_ID, &
                  PLASTICITY_dislotwin_ID, &
@@ -313,7 +311,6 @@ module material
    ELASTICITY_hooke_ID ,&
    PLASTICITY_none_ID, &
    PLASTICITY_isotropic_ID, &
-   PLASTICITY_J2_ID, &
    PLASTICITY_phenopowerlaw_ID, &
    PLASTICITY_phenoplus_ID, &
    PLASTICITY_dislotwin_ID, &
@@ -351,9 +348,6 @@ module material
    HYDROGENFLUX_cahnhilliard_ID, &
    HOMOGENIZATION_none_ID, &
    HOMOGENIZATION_isostrain_ID, &
-#ifdef HDF
-   material_NconstituentsPhase, &
-#endif
    HOMOGENIZATION_RGC_ID
 
  private :: &
@@ -982,8 +976,6 @@ subroutine material_parsePhase(fileUnit,myPart)
              phase_plasticity(section) = PLASTICITY_NONE_ID
            case (PLASTICITY_ISOTROPIC_label)
              phase_plasticity(section) = PLASTICITY_ISOTROPIC_ID
-           case (PLASTICITY_J2_label)
-             phase_plasticity(section) = PLASTICITY_J2_ID
            case (PLASTICITY_PHENOPOWERLAW_label)
              phase_plasticity(section) = PLASTICITY_PHENOPOWERLAW_ID
            case (PLASTICITY_PHENOPLUS_label)
@@ -1603,14 +1595,4 @@ subroutine material_populateGrains
 
 end subroutine material_populateGrains
 
-#ifdef HDF
-integer(pInt) pure function material_NconstituentsPhase(matID)
-
- implicit none
- integer(pInt), intent(in) :: matID
-
- material_NconstituentsPhase = count(microstructure_phase == matID)
-end function
-#endif
-
 end module material

code/mesh.f90

@@ -4525,17 +4525,9 @@ subroutine mesh_write_cellGeom
    VTK_geo, &
    VTK_con, &
    VTK_end 
-#ifdef HDF
- use IO, only: &
-   HDF5_mappingCells 
-#endif
  implicit none
  integer(I4P), dimension(1:mesh_Ncells)                                :: celltype
  integer(I4P), dimension(mesh_Ncells*(1_pInt+FE_maxNcellnodesPerCell)) :: cellconnection
-#ifdef HDF
- integer(pInt), dimension(mesh_Ncells*FE_maxNcellnodesPerCell) :: cellconnectionHDF5
- integer(pInt) :: j2=0_pInt
-#endif
  integer(I4P):: error
  integer(I4P):: g, c, e, CellID, i, j
 
@@ -4550,16 +4542,8 @@ subroutine mesh_write_cellGeom
      cellconnection(j+1_pInt:j+FE_NcellnodesPerCell(c)+1_pInt) &
        = [FE_NcellnodesPerCell(c),mesh_cell(1:FE_NcellnodesPerCell(c),i,e)-1_pInt]                 ! number of cellnodes per cell & list of global cellnode IDs belnging to this cell (cellnode counting starts at 0)
      j = j + FE_NcellnodesPerCell(c) + 1_pInt
-#ifdef HDF
-     cellconnectionHDF5(j2+1_pInt:j2+FE_NcellnodesPerCell(c)) &
-       = mesh_cell(1:FE_NcellnodesPerCell(c),i,e)-1_pInt
-  j2=j2 + FE_ncellnodesPerCell(c)
-#endif
    enddo
  enddo
-#ifdef HDF
- call HDF5_mappingCells(cellconnectionHDF5(1:j2))
-#endif
 
  error=VTK_ini(output_format = 'ASCII', &
        title=trim(getSolverJobName())//' cell mesh', &

code/plastic_isotropic.f90

@@ -497,28 +497,23 @@ real(pReal) :: &
  squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph_33,Tstar_sph_33)
  norm_Tstar_sph = sqrt(squarenorm_Tstar_sph) 
 
- if (param(instance)%dilatation) then
+ if (param(instance)%dilatation .and. norm_Tstar_sph > 0.0_pReal) then                              ! Tstar == 0 or J2 plascitiy --> both Li and dLi_dTstar are zero
-     if (norm_Tstar_sph <= 0.0_pReal) then                                                          ! Tstar == 0 --> both Li and dLi_dTstar are zero
+   gamma_dot = param(instance)%gdot0 &
-       Li = 0.0_pReal
+               * (sqrt(1.5_pReal) * norm_Tstar_sph / param(instance)%fTaylor / state(instance)%flowstress(of) ) &
-       dLi_dTstar_3333 = 0.0_pReal
+               **param(instance)%n
-     else
-       gamma_dot = param(instance)%gdot0 &
-                   * (sqrt(1.5_pReal) * norm_Tstar_sph / param(instance)%fTaylor / state(instance)%flowstress(of) ) &
-                   **param(instance)%n
 
-       Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/param(instance)%fTaylor
+   Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/param(instance)%fTaylor
 
-       !--------------------------------------------------------------------------------------------------
+   !--------------------------------------------------------------------------------------------------
-       ! Calculation of the tangent of Li
+   ! Calculation of the tangent of Li
-       forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
+   forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
-         dLi_dTstar_3333(k,l,m,n) = (param(instance)%n-1.0_pReal) * &
+     dLi_dTstar_3333(k,l,m,n) = (param(instance)%n-1.0_pReal) * &
-                                          Tstar_sph_33(k,l)*Tstar_sph_33(m,n) / squarenorm_Tstar_sph
+                                      Tstar_sph_33(k,l)*Tstar_sph_33(m,n) / squarenorm_Tstar_sph
-       forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
+   forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
-         dLi_dTstar_3333(k,l,k,l) = dLi_dTstar_3333(k,l,k,l) + 1.0_pReal
+     dLi_dTstar_3333(k,l,k,l) = dLi_dTstar_3333(k,l,k,l) + 1.0_pReal
 
-       dLi_dTstar_3333 = gamma_dot / param(instance)%fTaylor * &
+   dLi_dTstar_3333 = gamma_dot / param(instance)%fTaylor * &
-                                          dLi_dTstar_3333 / norm_Tstar_sph
+                                      dLi_dTstar_3333 / norm_Tstar_sph
-     endif
  else
   Li = 0.0_pReal
   dLi_dTstar_3333 = 0.0_pReal

code/plastic_j2.f90

@@ -1,564 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
-!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine for isotropic (J2) plasticity
-!> @details Isotropic (J2) Plasticity which resembles the phenopowerlaw plasticity without
-!! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an
-!! untextured polycrystal
-!--------------------------------------------------------------------------------------------------
-module plastic_j2
-#ifdef HDF
- use hdf5, only: &
-   HID_T
-#endif
-
- use prec, only: &
-   pReal,&
-   pInt
- 
- implicit none
- private
- integer(pInt),                       dimension(:),     allocatable,         public, protected :: &
-   plastic_j2_sizePostResults                                                                  !< cumulative size of post results
-   
- integer(pInt),                       dimension(:,:),   allocatable, target, public :: &
-   plastic_j2_sizePostResult                                                                   !< size of each post result output
-   
- character(len=64),                   dimension(:,:),   allocatable, target, public :: &
-   plastic_j2_output                                                                           !< name of each post result output
- 
- integer(pInt),                       dimension(:),     allocatable, target, public :: &
-   plastic_j2_Noutput                                                                          !< number of outputs per instance
- real(pReal),                         dimension(:),     allocatable,         private :: &
-   plastic_j2_fTaylor, &                                                                       !< Taylor factor
-   plastic_j2_tau0, &                                                                          !< initial plastic stress
-   plastic_j2_gdot0, &                                                                         !< reference velocity
-   plastic_j2_n, &                                                                             !< Visco-plastic parameter
-!--------------------------------------------------------------------------------------------------
-! h0 as function of h0 = A + B log (gammadot) 
-   plastic_j2_h0, &
-   plastic_j2_h0_slopeLnRate, &
-   plastic_j2_tausat, &                                                                        !< final plastic stress
-   plastic_j2_a, &
-   plastic_j2_aTolResistance, &
-   plastic_j2_aTolShear, &
-!--------------------------------------------------------------------------------------------------
-! tausat += (asinh((gammadot / SinhFitA)**(1 / SinhFitD)))**(1 / SinhFitC) / (SinhFitB * (gammadot / gammadot0)**(1/n))
-   plastic_j2_tausat_SinhFitA, &                                                               !< fitting parameter for normalized strain rate vs. stress function
-   plastic_j2_tausat_SinhFitB, &                                                               !< fitting parameter for normalized strain rate vs. stress function
-   plastic_j2_tausat_SinhFitC, &                                                               !< fitting parameter for normalized strain rate vs. stress function
-   plastic_j2_tausat_SinhFitD                                                                  !< fitting parameter for normalized strain rate vs. stress function
-
- enum, bind(c) 
-   enumerator :: undefined_ID, &
-                 flowstress_ID, &
-                 strainrate_ID
- end enum
- integer(kind(undefined_ID)),           dimension(:,:),   allocatable,        private :: & 
-   plastic_j2_outputID                                                                         !< ID of each post result output
- 
-
-#ifdef HDF 
- type plastic_j2_tOutput
-   real(pReal),                         dimension(:),     allocatable,         private :: &
-     flowstress, &
-     strainrate
- logical :: flowstressActive = .false., strainrateActive = .false.                ! if we can write the output block wise, this is not needed anymore because we can do an if(allocated(xxx))                                 
- end type plastic_j2_tOutput
- type(plastic_j2_tOutput), allocatable, dimension(:) :: plastic_j2_Output2
-integer(HID_T), allocatable, dimension(:) :: outID
-#endif
-
-
- public  :: &
-   plastic_j2_init, &
-   plastic_j2_LpAndItsTangent, &
-   plastic_j2_dotState, &
-   plastic_j2_postResults
-
-contains
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief module initialization
-!> @details reads in material parameters, allocates arrays, and does sanity checks
-!--------------------------------------------------------------------------------------------------
-subroutine plastic_j2_init(fileUnit)
- use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
-#ifdef HDF 
- use hdf5
-#endif 
- use debug, only: &
-   debug_level, &
-   debug_constitutive, &
-   debug_levelBasic
- use numerics, only: &
-   analyticJaco, &
-   worldrank, &
-   numerics_integrator
- use math, only: &
-   math_Mandel3333to66, &
-   math_Voigt66to3333
- use IO, only: &
-   IO_read, &
-   IO_lc, &
-   IO_getTag, &
-   IO_isBlank, &
-   IO_stringPos, &
-   IO_stringValue, &
-   IO_floatValue, &
-   IO_error, &
-   IO_timeStamp, &
-#ifdef HDF 
-   tempResults, &
-   HDF5_addGroup, &
-   HDF5_addScalarDataset,&
-#endif
-   IO_EOF
- use material, only: &
-   phase_plasticity, &
-   phase_plasticityInstance, &
-   phase_Noutput, &
-   PLASTICITY_J2_label, &
-   PLASTICITY_J2_ID, &
-   material_phase, &
-   plasticState, &
-   MATERIAL_partPhase
-   
- use lattice  
-
- implicit none
- integer(pInt), intent(in) :: fileUnit
- 
- 
- integer(pInt), allocatable, dimension(:) :: chunkPos
- integer(pInt) :: &
-   o, &
-   phase, & 
-   maxNinstance, &
-   instance, &
-   mySize, &
-   sizeDotState, &
-   sizeState, &
-   sizeDeltaState
- character(len=65536) :: &
-   tag  = '', &
-   line = ''
-  integer(pInt) :: NofMyPhase
-
-#ifdef HDF 
- character(len=5) :: &
-   str1
- integer(HID_T) :: ID,ID2,ID4
-#endif
-
- mainProcess: if (worldrank == 0) then 
-   write(6,'(/,a)')   ' <<<+-  constitutive_'//PLASTICITY_J2_label//' init  -+>>>'
-   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
-#include "compilation_info.f90"
- endif mainProcess
- 
- maxNinstance = int(count(phase_plasticity == PLASTICITY_J2_ID),pInt)
- if (maxNinstance == 0_pInt) return
-
- if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
-   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
-
-#ifdef HDF 
-  allocate(plastic_j2_Output2(maxNinstance))
-  allocate(outID(maxNinstance))
-#endif
-
- allocate(plastic_j2_sizePostResults(maxNinstance),                      source=0_pInt)
- allocate(plastic_j2_sizePostResult(maxval(phase_Noutput), maxNinstance),source=0_pInt)
- allocate(plastic_j2_output(maxval(phase_Noutput), maxNinstance))
-          plastic_j2_output = ''
- allocate(plastic_j2_outputID(maxval(phase_Noutput),maxNinstance),       source=undefined_ID)
- allocate(plastic_j2_Noutput(maxNinstance),                              source=0_pInt)
- allocate(plastic_j2_fTaylor(maxNinstance),                              source=0.0_pReal)
- allocate(plastic_j2_tau0(maxNinstance),                                 source=0.0_pReal)
- allocate(plastic_j2_gdot0(maxNinstance),                                source=0.0_pReal)
- allocate(plastic_j2_n(maxNinstance),                                    source=0.0_pReal)
- allocate(plastic_j2_h0(maxNinstance),                                   source=0.0_pReal)
- allocate(plastic_j2_h0_slopeLnRate(maxNinstance),                       source=0.0_pReal)
- allocate(plastic_j2_tausat(maxNinstance),                               source=0.0_pReal)
- allocate(plastic_j2_a(maxNinstance),                                    source=0.0_pReal)
- allocate(plastic_j2_aTolResistance(maxNinstance),                       source=0.0_pReal)
- allocate(plastic_j2_aTolShear     (maxNinstance),                       source=0.0_pReal)
- allocate(plastic_j2_tausat_SinhFitA(maxNinstance),                      source=0.0_pReal)
- allocate(plastic_j2_tausat_SinhFitB(maxNinstance),                      source=0.0_pReal)
- allocate(plastic_j2_tausat_SinhFitC(maxNinstance),                      source=0.0_pReal)
- allocate(plastic_j2_tausat_SinhFitD(maxNinstance),                      source=0.0_pReal)
- 
- rewind(fileUnit)
- phase = 0_pInt
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase)         ! wind forward to <phase>
-   line = IO_read(fileUnit)
- enddo
- 
- parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit                                                                                           
-   endif
-   if (IO_getTag(line,'[',']') /= '') then                                                          ! next section
-     phase = phase + 1_pInt                                                                         ! advance section counter
-     if (phase_plasticity(phase) == PLASTICITY_J2_ID) then
-       instance = phase_plasticityInstance(phase)
-     endif
-     cycle                                                                                          ! skip to next line
-   endif
-   if (phase > 0_pInt ) then; if (phase_plasticity(phase) == PLASTICITY_J2_ID) then                 ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran
-     instance = phase_plasticityInstance(phase)                                                     ! which instance of my plasticity is present phase
-     chunkPos = IO_stringPos(line) 
-     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
-
-     select case(tag)
-       case ('(output)')
-         select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
-           case ('flowstress')
-             plastic_j2_Noutput(instance) = plastic_j2_Noutput(instance) + 1_pInt
-             plastic_j2_outputID(plastic_j2_Noutput(instance),instance) = flowstress_ID
-             plastic_j2_output(plastic_j2_Noutput(instance),instance) = &
-                                                IO_lc(IO_stringValue(line,chunkPos,2_pInt))
-           case ('strainrate')
-             plastic_j2_Noutput(instance) = plastic_j2_Noutput(instance) + 1_pInt
-             plastic_j2_outputID(plastic_j2_Noutput(instance),instance) = strainrate_ID
-             plastic_j2_output(plastic_j2_Noutput(instance),instance) = &
-                                                IO_lc(IO_stringValue(line,chunkPos,2_pInt))
-           case default
-
-         end select
-       case ('tau0')
-         plastic_j2_tau0(instance)         = IO_floatValue(line,chunkPos,2_pInt)
-         if (plastic_j2_tau0(instance) < 0.0_pReal) &
-           call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
-       case ('gdot0')
-         plastic_j2_gdot0(instance)        = IO_floatValue(line,chunkPos,2_pInt)
-         if (plastic_j2_gdot0(instance) <= 0.0_pReal) &
-           call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
-       case ('n')
-         plastic_j2_n(instance)            = IO_floatValue(line,chunkPos,2_pInt)
-         if (plastic_j2_n(instance) <= 0.0_pReal) &
-           call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
-       case ('h0')
-         plastic_j2_h0(instance)           = IO_floatValue(line,chunkPos,2_pInt)
-       case ('h0_slope','slopelnrate')
-         plastic_j2_h0_slopeLnRate(instance)  = IO_floatValue(line,chunkPos,2_pInt)
-       case ('tausat')
-         plastic_j2_tausat(instance)          = IO_floatValue(line,chunkPos,2_pInt)
-         if (plastic_j2_tausat(instance) <= 0.0_pReal) &
-           call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
-       case ('tausat_sinhfita')
-         plastic_j2_tausat_SinhFitA(instance) = IO_floatValue(line,chunkPos,2_pInt)
-       case ('tausat_sinhfitb')
-         plastic_j2_tausat_SinhFitB(instance) = IO_floatValue(line,chunkPos,2_pInt)
-       case ('tausat_sinhfitc')
-         plastic_j2_tausat_SinhFitC(instance) = IO_floatValue(line,chunkPos,2_pInt)
-       case ('tausat_sinhfitd')
-         plastic_j2_tausat_SinhFitD(instance) = IO_floatValue(line,chunkPos,2_pInt)
-       case ('a', 'w0')
-         plastic_j2_a(instance)               = IO_floatValue(line,chunkPos,2_pInt)
-         if (plastic_j2_a(instance) <= 0.0_pReal) &
-           call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
-       case ('taylorfactor')
-         plastic_j2_fTaylor(instance)         = IO_floatValue(line,chunkPos,2_pInt)
-         if (plastic_j2_fTaylor(instance) <= 0.0_pReal) &
-           call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
-       case ('atol_resistance')
-         plastic_j2_aTolResistance(instance)  = IO_floatValue(line,chunkPos,2_pInt)
-         if (plastic_j2_aTolResistance(instance) <= 0.0_pReal) &
-           call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
-       case ('atol_shear')
-         plastic_j2_aTolShear(instance)  = IO_floatValue(line,chunkPos,2_pInt)
-
-       case default
-
-     end select
-   endif; endif
- enddo parsingFile
-
- initializeInstances: do phase = 1_pInt, size(phase_plasticity)
-   myPhase: if (phase_plasticity(phase) == PLASTICITY_j2_ID) then
-     NofMyPhase=count(material_phase==phase)
-     instance = phase_plasticityInstance(phase)
-!--------------------------------------------------------------------------------------------------
-!  sanity checks
-     if (plastic_j2_aTolShear(instance) <= 0.0_pReal) &
-       plastic_j2_aTolShear(instance) = 1.0e-6_pReal                                         ! default absolute tolerance 1e-6
-
-!--------------------------------------------------------------------------------------------------
-!  Determine size of postResults array
-     outputsLoop: do o = 1_pInt,plastic_j2_Noutput(instance)
-       select case(plastic_j2_outputID(o,instance))
-         case(flowstress_ID,strainrate_ID)
-           mySize = 1_pInt
-         case default
-       end select
-  
-       outputFound: if (mySize > 0_pInt) then
-         plastic_j2_sizePostResult(o,instance) = mySize
-         plastic_j2_sizePostResults(instance) = &
-         plastic_j2_sizePostResults(instance) + mySize
-       endif outputFound
-     enddo outputsLoop
-
-!--------------------------------------------------------------------------------------------------
-! allocate state arrays
-     sizeState    = 2_pInt
-     sizeDotState = sizeState
-     sizeDeltaState = 0_pInt
-     plasticState(phase)%sizeState = sizeState
-     plasticState(phase)%sizeDotState = sizeDotState
-     plasticState(phase)%sizeDeltaState = sizeDeltaState
-     plasticState(phase)%sizePostResults = plastic_j2_sizePostResults(instance)
-     plasticState(phase)%nSlip = 1
-     plasticState(phase)%nTwin = 0
-     plasticState(phase)%nTrans= 0
-     allocate(plasticState(phase)%aTolState          (   sizeState))
-     plasticState(phase)%aTolState(1) = plastic_j2_aTolResistance(instance)
-     plasticState(phase)%aTolState(2) = plastic_j2_aTolShear(instance)
-     allocate(plasticState(phase)%state0             (   sizeState,NofMyPhase))
-     plasticState(phase)%state0(1,1:NofMyPhase) = plastic_j2_tau0(instance)
-     plasticState(phase)%state0(2,1:NofMyPhase) = 0.0_pReal
-     allocate(plasticState(phase)%partionedState0    (   sizeState,NofMyPhase),source=0.0_pReal)
-     allocate(plasticState(phase)%subState0          (   sizeState,NofMyPhase),source=0.0_pReal)
-     allocate(plasticState(phase)%state              (   sizeState,NofMyPhase),source=0.0_pReal)
-     allocate(plasticState(phase)%dotState           (sizeDotState,NofMyPhase),source=0.0_pReal)
-     allocate(plasticState(phase)%deltaState       (sizeDeltaState,NofMyPhase),source=0.0_pReal)
-     if (.not. analyticJaco) then
-       allocate(plasticState(phase)%state_backup     (   sizeState,NofMyPhase),source=0.0_pReal)
-       allocate(plasticState(phase)%dotState_backup  (sizeDotState,NofMyPhase),source=0.0_pReal)
-     endif
-     if (any(numerics_integrator == 1_pInt)) then
-       allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase),source=0.0_pReal)
-       allocate(plasticState(phase)%previousDotState2(sizeDotState,NofMyPhase),source=0.0_pReal)
-     endif
-     if (any(numerics_integrator == 4_pInt)) &
-       allocate(plasticState(phase)%RK4dotState      (sizeDotState,NofMyPhase),source=0.0_pReal)
-     if (any(numerics_integrator == 5_pInt)) &
-       allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
-     plasticState(phase)%slipRate        => plasticState(phase)%dotState(2:2,1:NofMyPhase)
-     plasticState(phase)%accumulatedSlip => plasticState(phase)%state   (2:2,1:NofMyPhase)
-   endif myPhase
- enddo initializeInstances
-
-end subroutine plastic_j2_init
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief calculates plastic velocity gradient and its tangent
-!--------------------------------------------------------------------------------------------------
-subroutine plastic_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
- use math, only: &
-   math_mul6x6, &
-   math_Mandel6to33, &
-   math_Plain3333to99, &
-   math_deviatoric33, &
-   math_mul33xx33
- use material, only: &
-   phaseAt, phasememberAt, &
-   plasticState, &
-   material_phase, &
-   phase_plasticityInstance
-
- implicit none
- real(pReal), dimension(3,3), intent(out) :: &
-   Lp                                                                                               !< plastic velocity gradient
- real(pReal), dimension(9,9), intent(out) :: &
-   dLp_dTstar99                                                                                     !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
-
- real(pReal), dimension(6),   intent(in) :: &
-   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor in Mandel notation
- integer(pInt),               intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
-   ip, &                                                                                            !< integration point
-   el                                                                                               !< element
-
- real(pReal), dimension(3,3) :: &
-   Tstar_dev_33                                                                                     !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
- real(pReal), dimension(3,3,3,3) :: &
-   dLp_dTstar_3333                                                                                  !< derivative of Lp with respect to Tstar as 4th order tensor
- real(pReal) :: &
-   gamma_dot, &                                                                                     !< strainrate
-   norm_Tstar_dev, &                                                                                !< euclidean norm of Tstar_dev
-   squarenorm_Tstar_dev                                                                             !< square of the euclidean norm of Tstar_dev
- integer(pInt) :: &
-   instance, &
-   k, l, m, n
-
- 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)
- norm_Tstar_dev = sqrt(squarenorm_Tstar_dev) 
-
- if (norm_Tstar_dev <= 0.0_pReal) then                                                              ! Tstar == 0 --> both Lp and dLp_dTstar are zero
-   Lp = 0.0_pReal
-   dLp_dTstar99 = 0.0_pReal
- else
-   gamma_dot = plastic_j2_gdot0(instance) &
-             * (sqrt(1.5_pReal) * norm_Tstar_dev / (plastic_j2_fTaylor(instance) * &
-               plasticState(phaseAt(ipc,ip,el))%state(1,phasememberAt(ipc,ip,el)))) &
-                                                  **plastic_j2_n(instance)
-
-   Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/plastic_j2_fTaylor(instance)
-
-!--------------------------------------------------------------------------------------------------
-! Calculation of the tangent of Lp
-   forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
-     dLp_dTstar_3333(k,l,m,n) = (plastic_j2_n(instance)-1.0_pReal) * &
-                                      Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev
-   forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
-     dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal
-   forall (k=1_pInt:3_pInt,m=1_pInt:3_pInt) &
-     dLp_dTstar_3333(k,k,m,m) = dLp_dTstar_3333(k,k,m,m) - 1.0_pReal/3.0_pReal
-   dLp_dTstar99 = math_Plain3333to99(gamma_dot / plastic_j2_fTaylor(instance) * &
-                                      dLp_dTstar_3333 / norm_Tstar_dev)
- end if
-end subroutine plastic_j2_LpAndItsTangent
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief calculates the rate of change of microstructure
-!--------------------------------------------------------------------------------------------------
-subroutine plastic_j2_dotState(Tstar_v,ipc,ip,el)
- use math, only: &
-   math_mul6x6
- use material, only: &
-   phaseAt, phasememberAt, &
-   plasticState, &
-   material_phase, &
-   phase_plasticityInstance
- 
- implicit none
- real(pReal), dimension(6), intent(in):: &
-   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor in Mandel notation
- integer(pInt),             intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
-   ip, &                                                                                            !< integration point
-   el                                                                                               !< element
- real(pReal), dimension(6) :: &
-   Tstar_dev_v                                                                                      !< deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation
- real(pReal) :: &
-   gamma_dot, &                                                                                     !< strainrate
-   hardening, &                                                                                     !< hardening coefficient
-   saturation, &                                                                                    !< saturation resistance
-   norm_Tstar_dev                                                                                   !< euclidean norm of Tstar_dev
- integer(pInt) :: &
-   instance, &                                                                                      !< instance of my instance (unique number of my constitutive model)
-   of, &                                                                                            !< shortcut notation for offset position in state array
-   ph                                                                                               !< shortcut notation for phase ID (unique number of all phases, regardless of constitutive model)
-
- of = phasememberAt(ipc,ip,el)
- ph = phaseAt(ipc,ip,el)
- instance = phase_plasticityInstance(material_phase(ipc,ip,el))
-
-!--------------------------------------------------------------------------------------------------
-! norm of deviatoric part of 2nd Piola-Kirchhoff stress
- Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
- Tstar_dev_v(4:6) = Tstar_v(4:6)
- norm_Tstar_dev = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
-
-!--------------------------------------------------------------------------------------------------
-! strain rate 
- gamma_dot = plastic_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev & 
-            / &!-----------------------------------------------------------------------------------
-           (plastic_j2_fTaylor(instance)*plasticState(ph)%state(1,of)) )**plastic_j2_n(instance)
- 
-!--------------------------------------------------------------------------------------------------
-! hardening coefficient
- if (abs(gamma_dot) > 1e-12_pReal) then
-   if (abs(plastic_j2_tausat_SinhFitA(instance)) <= tiny(0.0_pReal)) then
-     saturation = plastic_j2_tausat(instance)
-   else
-     saturation = (  plastic_j2_tausat(instance) &
-                   + ( log(  ( gamma_dot / plastic_j2_tausat_SinhFitA(instance)&
-                               )**(1.0_pReal / plastic_j2_tausat_SinhFitD(instance))&
-                            + sqrt(  ( gamma_dot / plastic_j2_tausat_SinhFitA(instance) &
-                                      )**(2.0_pReal / plastic_j2_tausat_SinhFitD(instance)) &
-                                   + 1.0_pReal ) &
-                            ) & ! asinh(K) = ln(K + sqrt(K^2 +1))
-                       )**(1.0_pReal / plastic_j2_tausat_SinhFitC(instance)) &
-                   / (  plastic_j2_tausat_SinhFitB(instance) &
-                      * (gamma_dot / plastic_j2_gdot0(instance))**(1.0_pReal / plastic_j2_n(instance)) &
-                      ) &
-                   )
-   endif
-   hardening = ( plastic_j2_h0(instance) + plastic_j2_h0_slopeLnRate(instance) * log(gamma_dot) ) &
-               * abs( 1.0_pReal - plasticState(ph)%state(1,of)/saturation )**plastic_j2_a(instance) &
-               * sign(1.0_pReal, 1.0_pReal - plasticState(ph)%state(1,of)/saturation)
- else
-   hardening = 0.0_pReal
- endif
-
-  plasticState(ph)%dotState(1,of) = hardening * gamma_dot
-  plasticState(ph)%dotState(2,of) =             gamma_dot
-
-end subroutine plastic_j2_dotState
-
-!--------------------------------------------------------------------------------------------------
-!> @brief return array of constitutive results
-!--------------------------------------------------------------------------------------------------
-function plastic_j2_postResults(Tstar_v,ipc,ip,el)
- use math, only: &
-   math_mul6x6
- use material, only: &
-   material_phase, &
-   plasticState, &
-   phaseAt, phasememberAt, &
-   phase_plasticityInstance
-
- implicit none
- real(pReal), dimension(6),  intent(in) :: &
-   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor in Mandel notation
- integer(pInt),              intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
-   ip, &                                                                                            !< integration point
-   el                                                                                               !< element
- real(pReal), dimension(plastic_j2_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
-                                           plastic_j2_postResults
-
- real(pReal), dimension(6) :: &
-   Tstar_dev_v                                                                                      ! deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation
- real(pReal) :: &
-   norm_Tstar_dev                                                                                   ! euclidean norm of Tstar_dev
- integer(pInt) :: &
-   instance, &                                                                                      !< instance of my instance (unique number of my constitutive model)
-   of, &                                                                                            !< shortcut notation for offset position in state array
-   ph, &                                                                                            !< shortcut notation for phase ID (unique number of all phases, regardless of constitutive model)
-   c, &
-   o
-
- of = phasememberAt(ipc,ip,el)
- ph = phaseAt(ipc,ip,el)
- instance = phase_plasticityInstance(material_phase(ipc,ip,el))
- 
-!--------------------------------------------------------------------------------------------------
-! calculate deviatoric part of 2nd Piola-Kirchhoff stress and its norm
- Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
- Tstar_dev_v(4:6) = Tstar_v(4:6)
- norm_Tstar_dev = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
- 
- c = 0_pInt
- plastic_j2_postResults = 0.0_pReal
-
- outputsLoop: do o = 1_pInt,plastic_j2_Noutput(instance)
-   select case(plastic_j2_outputID(o,instance))
-     case (flowstress_ID)
-       plastic_j2_postResults(c+1_pInt) = plasticState(ph)%state(1,of)
-       c = c + 1_pInt
-     case (strainrate_ID)
-       plastic_j2_postResults(c+1_pInt) = &
-                plastic_j2_gdot0(instance) * (            sqrt(1.5_pReal) * norm_Tstar_dev & 
-             / &!----------------------------------------------------------------------------------
-              (plastic_j2_fTaylor(instance) * plasticState(ph)%state(1,of)) ) ** plastic_j2_n(instance)
-       c = c + 1_pInt
-   end select
- enddo outputsLoop
-
-end function plastic_j2_postResults
-
-
-end module plastic_j2

config/Phase_Isotropic_AluminumIsotropic.config

@@ -3,7 +3,7 @@
 # Kuo, J. C., Mikrostrukturmechanik von Bikristallen mit Kippkorngrenzen. Shaker-Verlag 2004. http://edoc.mpg.de/204079
 
 elasticity	hooke
-plasticity	j2
+plasticity	isotropic
 
 (output) flowstress
 (output) strainrate

processing/post/addCompatibilityMismatch.py

@@ -36,7 +36,7 @@ parser.add_option('--no-volume','-v',
                   dest = 'volume',
                   action = 'store_false',
                   help = 'omit volume mismatch')
-parser.set_defaults(coords   = 'ipinitialcoord',
+parser.set_defaults(coords   = 'pos',
                     defgrad  = 'f',
                     shape = True,
                     volume = True,

processing/post/blowUp.py

@@ -76,9 +76,9 @@ for name in filenames:
 
 # --------------- figure out size and grid ---------------------------------------------------------
 
-  table.data_readArray(options.coords)
+  table.data_readArray()
 
-  coords = [np.unique(table.data[:,i]) for i in xrange(3)]
+  coords = [np.unique(table.data[:,colCoord+i]) for i in xrange(3)]
   mincorner = np.array(map(min,coords))
   maxcorner = np.array(map(max,coords))
   grid   = np.array(map(len,coords),'i')

processing/post/postResults.py

@@ -1003,7 +1003,8 @@ fileOpen = False
 assembleHeader = True
 header = []
 standard = ['inc'] + \
-           ['time'] if options.time else [] + \
+           {True: ['time'],
+            False:[]}[options.time] + \
            ['elem','node','ip','grain','1_pos','2_pos','3_pos']
 
 # ---------------------------   loop over positions   --------------------------------

processing/post/vtk_rectilinearGrid.py

@@ -38,7 +38,7 @@ parser.set_defaults(coords = 'pos',
 if filenames == []: filenames = [None]
 
 for name in filenames:
-  isGeom = name.endswith('.geom')
+  isGeom = name is not None and name.endswith('.geom')
   try:    table = damask.ASCIItable(name = name,
                                     buffered = False,
                                     labeled  = not isGeom,