223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

cleaning

This commit is contained in:
Martin Diehl 2019-06-11 09:48:07 +02:00
parent 9dfe71aa06
commit b2409d6998
13 changed files with 1587 additions and 1700 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/DAMASK_interface.f90
View File

@ -40,12 +40,6 @@ module DAMASK_interface
setSIGTERM, &
setSIGUSR1, &
setSIGUSR2
private :: &
setWorkingDirectory, &
getGeometryFile, &
getLoadCaseFile, &
rectifyPath, &
makeRelativePath
contains

306
src/HDF5_utilities.f90
View File

@ -9,9 +9,7 @@ module HDF5_utilities
use IO
use rotations
use numerics
#if defined(PETSc) || defined(DAMASK_HDF5)
use HDF5
#endif
#ifdef PETSc
use PETSC
#endif
@ -21,7 +19,7 @@ module HDF5_utilities
#if defined(PETSc) || defined(DAMASK_HDF5)
!--------------------------------------------------------------------------------------------------
!> @brief reads pInt or pReal data of defined shape from file ! ToDo: order of arguments wrong
!> @brief reads integer or float data of defined shape from file ! ToDo: order of arguments wrong
!> @details for parallel IO, all dimension except for the last need to match
!--------------------------------------------------------------------------------------------------
interface HDF5_read
@ -44,7 +42,7 @@ module HDF5_utilities
end interface HDF5_read
!--------------------------------------------------------------------------------------------------
!> @brief writes pInt or pReal data of defined shape to file ! ToDo: order of arguments wrong
!> @brief writes integer or real data of defined shape to file ! ToDo: order of arguments wrong
!> @details for parallel IO, all dimension except for the last need to match
!--------------------------------------------------------------------------------------------------
interface HDF5_write
@ -69,7 +67,7 @@ module HDF5_utilities
end interface HDF5_write
!--------------------------------------------------------------------------------------------------
!> @brief attached attributes of type char,pInt or pReal to a file/dataset/group
!> @brief attached attributes of type char, integer or real to a file/dataset/group
!--------------------------------------------------------------------------------------------------
interface HDF5_addAttribute
module procedure HDF5_addAttribute_str
@ -114,7 +112,7 @@ subroutine HDF5_utilities_init
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_pInt,ext_msg='Default integer size does not match H5T_NATIVE_INTEGER')
call IO_error(0,ext_msg='Default integer size does not match H5T_NATIVE_INTEGER')
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)')
@ -144,30 +142,30 @@ integer(HID_T) function HDF5_openFile(fileName,mode,parallel) ! ToDo: simply "op
endif
call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5pcreate_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_openFile: h5pcreate_f')
#ifdef PETSc
if (present(parallel)) then; if (parallel) then
call h5pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5pset_fapl_mpio_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_openFile: h5pset_fapl_mpio_f')
endif; endif
#endif
if (m == 'w') then
call h5fcreate_f(fileName,H5F_ACC_TRUNC_F,HDF5_openFile,hdferr,access_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5fcreate_f (w)')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_openFile: h5fcreate_f (w)')
elseif(m == 'a') then
call h5fopen_f(fileName,H5F_ACC_RDWR_F,HDF5_openFile,hdferr,access_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5fopen_f (a)')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_openFile: h5fopen_f (a)')
elseif(m == 'r') then
call h5fopen_f(fileName,H5F_ACC_RDONLY_F,HDF5_openFile,hdferr,access_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5fopen_f (r)')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_openFile: h5fopen_f (r)')
else
call IO_error(1_pInt,ext_msg='HDF5_openFile: h5fopen_f unknown access mode: '//trim(m))
call IO_error(1,ext_msg='HDF5_openFile: h5fopen_f unknown access mode: '//trim(m))
endif
call h5pclose_f(plist_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_openFile: h5pclose_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_openFile: h5pclose_f')
end function HDF5_openFile
@ -182,7 +180,7 @@ subroutine HDF5_closeFile(fileHandle)
integer :: hdferr
call h5fclose_f(fileHandle,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_closeFile: h5fclose_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_closeFile: h5fclose_f')
end subroutine HDF5_closeFile
@ -201,19 +199,19 @@ integer(HID_T) function HDF5_addGroup(fileHandle,groupName)
!-------------------------------------------------------------------------------------------------
! creating a property list for data access properties
call h5pcreate_f(H5P_GROUP_ACCESS_F, aplist_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_addGroup: h5pcreate_f ('//trim(groupName)//')')
if (hdferr < 0) call IO_error(1,ext_msg = 'HDF5_addGroup: h5pcreate_f ('//trim(groupName)//')')
!-------------------------------------------------------------------------------------------------
! setting I/O mode to collective
#ifdef PETSc
call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_addGroup: h5pset_all_coll_metadata_ops_f ('//trim(groupName)//')')
if (hdferr < 0) call IO_error(1,ext_msg = 'HDF5_addGroup: h5pset_all_coll_metadata_ops_f ('//trim(groupName)//')')
#endif
!-------------------------------------------------------------------------------------------------
! Create group
call h5gcreate_f(fileHandle, trim(groupName), HDF5_addGroup, hdferr, OBJECT_NAMELEN_DEFAULT_F,gapl_id = aplist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_addGroup: h5gcreate_f ('//trim(groupName)//')')
if (hdferr < 0) call IO_error(1,ext_msg = 'HDF5_addGroup: h5gcreate_f ('//trim(groupName)//')')
call h5pclose_f(aplist_id,hdferr)
@ -237,19 +235,19 @@ integer(HID_T) function HDF5_openGroup(fileHandle,groupName)
!-------------------------------------------------------------------------------------------------
! creating a property list for data access properties
call h5pcreate_f(H5P_GROUP_ACCESS_F, aplist_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_openGroup: h5pcreate_f ('//trim(groupName)//')')
if (hdferr < 0) call IO_error(1,ext_msg = 'HDF5_openGroup: h5pcreate_f ('//trim(groupName)//')')
!-------------------------------------------------------------------------------------------------
! setting I/O mode to collective
#ifdef PETSc
call h5pget_all_coll_metadata_ops_f(aplist_id, is_collective, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_openGroup: h5pset_all_coll_metadata_ops_f ('//trim(groupName)//')')
if (hdferr < 0) call IO_error(1,ext_msg = 'HDF5_openGroup: h5pset_all_coll_metadata_ops_f ('//trim(groupName)//')')
#endif
!-------------------------------------------------------------------------------------------------
! opening the group
call h5gopen_f(fileHandle, trim(groupName), HDF5_openGroup, hdferr, gapl_id = aplist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_openGroup: h5gopen_f ('//trim(groupName)//')')
if (hdferr < 0) call IO_error(1,ext_msg = 'HDF5_openGroup: h5gopen_f ('//trim(groupName)//')')
call h5pclose_f(aplist_id,hdferr)
@ -265,7 +263,7 @@ subroutine HDF5_closeGroup(group_id)
integer :: hdferr
call h5gclose_f(group_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_closeGroup: h5gclose_f (el is ID)', el = int(group_id,pInt))
if (hdferr < 0) call IO_error(1,ext_msg = 'HDF5_closeGroup: h5gclose_f (el is ID)', el = int(group_id))
end subroutine HDF5_closeGroup
@ -288,11 +286,11 @@ logical function HDF5_objectExists(loc_id,path)
endif
call h5lexists_f(loc_id, p, HDF5_objectExists, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_objectExists: h5oexists_by_name_f')
if (hdferr < 0) call IO_error(1,ext_msg = 'HDF5_objectExists: h5oexists_by_name_f')
if(HDF5_objectExists) then
call h5oexists_by_name_f(loc_id, p, HDF5_objectExists, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_objectExists: h5oexists_by_name_f')
if (hdferr < 0) call IO_error(1,ext_msg = 'HDF5_objectExists: h5oexists_by_name_f')
endif
end function HDF5_objectExists
@ -319,27 +317,27 @@ subroutine HDF5_addAttribute_str(loc_id,attrLabel,attrValue,path)
endif
call h5screate_f(H5S_SCALAR_F,space_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5screate_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_str: h5screate_f')
call h5tcopy_f(H5T_NATIVE_CHARACTER, type_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5tcopy_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_str: 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='HDF5_addAttribute_str: h5tset_size_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_str: h5tset_size_f')
call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5aexists_by_name_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_str: h5aexists_by_name_f')
if (attrExists) then
call h5adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5adelete_by_name_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_str: h5adelete_by_name_f')
endif
call h5acreate_by_name_f(loc_id,trim(p),trim(attrLabel),type_id,space_id,attr_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5acreate_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_str: 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='HDF5_addAttribute_str: h5awrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_str: h5awrite_f')
call h5aclose_f(attr_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5aclose_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_str: h5aclose_f')
call h5tclose_f(type_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5tclose_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_str: h5tclose_f')
call h5sclose_f(space_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_str: h5sclose_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_str: h5sclose_f')
end subroutine HDF5_addAttribute_str
@ -351,7 +349,7 @@ subroutine HDF5_addAttribute_int(loc_id,attrLabel,attrValue,path)
integer(HID_T), intent(in) :: loc_id
character(len=*), intent(in) :: attrLabel
integer(pInt), intent(in) :: attrValue
integer, intent(in) :: attrValue
character(len=*), intent(in), optional :: path
integer :: hdferr
@ -366,21 +364,21 @@ subroutine HDF5_addAttribute_int(loc_id,attrLabel,attrValue,path)
endif
call h5screate_f(H5S_SCALAR_F,space_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int: h5screate_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int: h5screate_f')
call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int: h5aexists_by_name_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int: h5aexists_by_name_f')
if (attrExists) then
call h5adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int: h5adelete_by_name_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int: h5adelete_by_name_f')
endif
call h5acreate_by_name_f(loc_id,trim(p),trim(attrLabel),H5T_NATIVE_INTEGER,space_id,attr_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int: h5acreate_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int: h5acreate_f')
call h5awrite_f(attr_id, H5T_NATIVE_INTEGER, attrValue, int([1],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int: h5awrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int: h5awrite_f')
call h5aclose_f(attr_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int: h5tclose_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int: h5tclose_f')
call h5sclose_f(space_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int: h5sclose_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int: h5sclose_f')
end subroutine HDF5_addAttribute_int
@ -407,21 +405,21 @@ subroutine HDF5_addAttribute_real(loc_id,attrLabel,attrValue,path)
endif
call h5screate_f(H5S_SCALAR_F,space_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_real: h5screate_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_real: h5screate_f')
call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_real: h5aexists_by_name_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_real: h5aexists_by_name_f')
if (attrExists) then
call h5adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_real: h5adelete_by_name_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_real: h5adelete_by_name_f')
endif
call h5acreate_by_name_f(loc_id,trim(p),trim(attrLabel),H5T_NATIVE_DOUBLE,space_id,attr_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_real: h5acreate_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_real: h5acreate_f')
call h5awrite_f(attr_id, H5T_NATIVE_DOUBLE, attrValue, int([1],HSIZE_T), hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_real: h5awrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_real: h5awrite_f')
call h5aclose_f(attr_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_real: h5tclose_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_real: h5tclose_f')
call h5sclose_f(space_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_real: h5sclose_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_real: h5sclose_f')
end subroutine HDF5_addAttribute_real
@ -433,7 +431,7 @@ subroutine HDF5_addAttribute_int_array(loc_id,attrLabel,attrValue,path)
integer(HID_T), intent(in) :: loc_id
character(len=*), intent(in) :: attrLabel
integer(pInt), intent(in), dimension(:) :: attrValue
integer, intent(in), dimension(:) :: attrValue
character(len=*), intent(in), optional :: path
integer :: hdferr
@ -451,21 +449,21 @@ subroutine HDF5_addAttribute_int_array(loc_id,attrLabel,attrValue,path)
array_size = size(attrValue,kind=HSIZE_T)
call h5screate_simple_f(1, array_size, space_id, hdferr, array_size)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int_array: h5screate_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int_array: h5screate_f')
call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int_array: h5aexists_by_name_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int_array: h5aexists_by_name_f')
if (attrExists) then
call h5adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int_array: h5adelete_by_name_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int_array: h5adelete_by_name_f')
endif
call h5acreate_by_name_f(loc_id,trim(p),trim(attrLabel),H5T_NATIVE_INTEGER,space_id,attr_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int_array: h5acreate_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int_array: h5acreate_f')
call h5awrite_f(attr_id, H5T_NATIVE_INTEGER, attrValue, array_size, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int_array: h5awrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int_array: h5awrite_f')
call h5aclose_f(attr_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int_array: h5tclose_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int_array: h5tclose_f')
call h5sclose_f(space_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int_array: h5sclose_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int_array: h5sclose_f')
end subroutine HDF5_addAttribute_int_array
@ -495,21 +493,21 @@ subroutine HDF5_addAttribute_real_array(loc_id,attrLabel,attrValue,path)
array_size = size(attrValue,kind=HSIZE_T)
call h5screate_simple_f(1, array_size, space_id, hdferr, array_size)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int_array: h5screate_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int_array: h5screate_f')
call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int_array: h5aexists_by_name_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int_array: h5aexists_by_name_f')
if (attrExists) then
call h5adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int_array: h5adelete_by_name_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int_array: h5adelete_by_name_f')
endif
call h5acreate_by_name_f(loc_id,trim(p),trim(attrLabel),H5T_NATIVE_DOUBLE,space_id,attr_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int_array: h5acreate_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int_array: h5acreate_f')
call h5awrite_f(attr_id, H5T_NATIVE_DOUBLE, attrValue, array_size, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int_array: h5awrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int_array: h5awrite_f')
call h5aclose_f(attr_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int_array: h5tclose_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int_array: h5tclose_f')
call h5sclose_f(space_id,hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_addAttribute_int_array: h5sclose_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_addAttribute_int_array: h5sclose_f')
end subroutine HDF5_addAttribute_real_array
@ -525,19 +523,19 @@ subroutine HDF5_setLink(loc_id,target_name,link_name)
logical :: linkExists
call h5lexists_f(loc_id, link_name,linkExists, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_setLink: h5lexists_soft_f ('//trim(link_name)//')')
if (hdferr < 0) call IO_error(1,ext_msg = 'HDF5_setLink: h5lexists_soft_f ('//trim(link_name)//')')
if (linkExists) then
call h5ldelete_f(loc_id,link_name, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_setLink: h5ldelete_soft_f ('//trim(link_name)//')')
if (hdferr < 0) call IO_error(1,ext_msg = 'HDF5_setLink: h5ldelete_soft_f ('//trim(link_name)//')')
endif
call h5lcreate_soft_f(target_name, loc_id, link_name, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg = 'HDF5_setLink: h5lcreate_soft_f ('//trim(target_name)//' '//trim(link_name)//')')
if (hdferr < 0) call IO_error(1,ext_msg = 'HDF5_setLink: h5lcreate_soft_f ('//trim(target_name)//' '//trim(link_name)//')')
end subroutine HDF5_setLink
!--------------------------------------------------------------------------------------------------
!> @brief read dataset of type pReal with 1 dimension
!> @brief read dataset of type real with 1 dimension
!--------------------------------------------------------------------------------------------------
subroutine HDF5_read_real1(loc_id,dataset,datasetName,parallel)
@ -570,14 +568,14 @@ subroutine HDF5_read_real1(loc_id,dataset,datasetName,parallel)
call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_real1: h5dread_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_read_real1: h5dread_f')
call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id)
end subroutine HDF5_read_real1
!--------------------------------------------------------------------------------------------------
!> @brief read dataset of type pReal with 2 dimensions
!> @brief read dataset of type real with 2 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_read_real2(loc_id,dataset,datasetName,parallel)
@ -610,14 +608,14 @@ subroutine HDF5_read_real2(loc_id,dataset,datasetName,parallel)
call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_real2: h5dread_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_read_real2: h5dread_f')
call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id)
end subroutine HDF5_read_real2
!--------------------------------------------------------------------------------------------------
!> @brief read dataset of type pReal with 2 dimensions
!> @brief read dataset of type real with 2 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_read_real3(loc_id,dataset,datasetName,parallel)
@ -650,14 +648,14 @@ subroutine HDF5_read_real3(loc_id,dataset,datasetName,parallel)
call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_real3: h5dread_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_read_real3: h5dread_f')
call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id)
end subroutine HDF5_read_real3
!--------------------------------------------------------------------------------------------------
!> @brief read dataset of type pReal with 4 dimensions
!> @brief read dataset of type real with 4 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_read_real4(loc_id,dataset,datasetName,parallel)
@ -690,14 +688,14 @@ subroutine HDF5_read_real4(loc_id,dataset,datasetName,parallel)
call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_real4: h5dread_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_read_real4: h5dread_f')
call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id)
end subroutine HDF5_read_real4
!--------------------------------------------------------------------------------------------------
!> @brief read dataset of type pReal with 5 dimensions
!> @brief read dataset of type real with 5 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_read_real5(loc_id,dataset,datasetName,parallel)
@ -730,14 +728,14 @@ subroutine HDF5_read_real5(loc_id,dataset,datasetName,parallel)
call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_real5: h5dread_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_read_real5: h5dread_f')
call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id)
end subroutine HDF5_read_real5
!--------------------------------------------------------------------------------------------------
!> @brief read dataset of type pReal with 6 dimensions
!> @brief read dataset of type real with 6 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_read_real6(loc_id,dataset,datasetName,parallel)
@ -770,14 +768,14 @@ subroutine HDF5_read_real6(loc_id,dataset,datasetName,parallel)
call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_real6: h5dread_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_read_real6: h5dread_f')
call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id)
end subroutine HDF5_read_real6
!--------------------------------------------------------------------------------------------------
!> @brief read dataset of type pReal with 7 dimensions
!> @brief read dataset of type real with 7 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_read_real7(loc_id,dataset,datasetName,parallel)
@ -810,7 +808,7 @@ subroutine HDF5_read_real7(loc_id,dataset,datasetName,parallel)
call h5dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_real7: h5dread_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_read_real7: h5dread_f')
call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id)
@ -818,7 +816,7 @@ end subroutine HDF5_read_real7
!--------------------------------------------------------------------------------------------------
!> @brief read dataset of type pInt with 1 dimension
!> @brief read dataset of type integer with 1 dimension
!--------------------------------------------------------------------------------------------------
subroutine HDF5_read_int1(loc_id,dataset,datasetName,parallel)
@ -851,14 +849,14 @@ subroutine HDF5_read_int1(loc_id,dataset,datasetName,parallel)
call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_int1: h5dread_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_read_int1: h5dread_f')
call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id)
end subroutine HDF5_read_int1
!--------------------------------------------------------------------------------------------------
!> @brief read dataset of type pInt with 2 dimensions
!> @brief read dataset of type integer with 2 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_read_int2(loc_id,dataset,datasetName,parallel)
@ -891,14 +889,14 @@ subroutine HDF5_read_int2(loc_id,dataset,datasetName,parallel)
call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_int2: h5dread_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_read_int2: h5dread_f')
call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id)
end subroutine HDF5_read_int2
!--------------------------------------------------------------------------------------------------
!> @brief read dataset of type pInt with 3 dimensions
!> @brief read dataset of type integer with 3 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_read_int3(loc_id,dataset,datasetName,parallel)
@ -931,14 +929,14 @@ subroutine HDF5_read_int3(loc_id,dataset,datasetName,parallel)
call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_int3: h5dread_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_read_int3: h5dread_f')
call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id)
end subroutine HDF5_read_int3
!--------------------------------------------------------------------------------------------------
!> @brief read dataset of type pInt withh 4 dimensions
!> @brief read dataset of type integer withh 4 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_read_int4(loc_id,dataset,datasetName,parallel)
@ -971,14 +969,14 @@ subroutine HDF5_read_int4(loc_id,dataset,datasetName,parallel)
call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_int4: h5dread_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_read_int4: h5dread_f')
call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id)
end subroutine HDF5_read_int4
!--------------------------------------------------------------------------------------------------
!> @brief read dataset of type pInt with 5 dimensions
!> @brief read dataset of type integer with 5 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_read_int5(loc_id,dataset,datasetName,parallel)
@ -1011,14 +1009,14 @@ subroutine HDF5_read_int5(loc_id,dataset,datasetName,parallel)
call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_int5: h5dread_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_read_int5: h5dread_f')
call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id)
end subroutine HDF5_read_int5
!--------------------------------------------------------------------------------------------------
!> @brief read dataset of type pInt with 6 dimensions
!> @brief read dataset of type integer with 6 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_read_int6(loc_id,dataset,datasetName,parallel)
@ -1051,14 +1049,14 @@ subroutine HDF5_read_int6(loc_id,dataset,datasetName,parallel)
call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_int6: h5dread_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_read_int6: h5dread_f')
call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id)
end subroutine HDF5_read_int6
!--------------------------------------------------------------------------------------------------
!> @brief read dataset of type pInt with 7 dimensions
!> @brief read dataset of type integer with 7 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_read_int7(loc_id,dataset,datasetName,parallel)
@ -1091,7 +1089,7 @@ subroutine HDF5_read_int7(loc_id,dataset,datasetName,parallel)
call h5dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_read_int7: h5dread_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_read_int7: h5dread_f')
call finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id)
@ -1099,7 +1097,7 @@ end subroutine HDF5_read_int7
!--------------------------------------------------------------------------------------------------
!> @brief write dataset of type pReal with 1 dimension
!> @brief write dataset of type real with 1 dimension
!--------------------------------------------------------------------------------------------------
subroutine HDF5_write_real1(loc_id,dataset,datasetName,parallel)
@ -1132,7 +1130,7 @@ subroutine HDF5_write_real1(loc_id,dataset,datasetName,parallel)
if (product(totalShape) /= 0) then
call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_real1: h5dwrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_write_real1: h5dwrite_f')
endif
call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
@ -1140,7 +1138,7 @@ subroutine HDF5_write_real1(loc_id,dataset,datasetName,parallel)
end subroutine HDF5_write_real1
!--------------------------------------------------------------------------------------------------
!> @brief write dataset of type pReal with 2 dimensions
!> @brief write dataset of type real with 2 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_write_real2(loc_id,dataset,datasetName,parallel)
@ -1173,7 +1171,7 @@ subroutine HDF5_write_real2(loc_id,dataset,datasetName,parallel)
if (product(totalShape) /= 0) then
call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_real2: h5dwrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_write_real2: h5dwrite_f')
endif
call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
@ -1181,7 +1179,7 @@ subroutine HDF5_write_real2(loc_id,dataset,datasetName,parallel)
end subroutine HDF5_write_real2
!--------------------------------------------------------------------------------------------------
!> @brief write dataset of type pReal with 3 dimensions
!> @brief write dataset of type real with 3 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_write_real3(loc_id,dataset,datasetName,parallel)
@ -1214,7 +1212,7 @@ subroutine HDF5_write_real3(loc_id,dataset,datasetName,parallel)
if (product(totalShape) /= 0) then
call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_real3: h5dwrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_write_real3: h5dwrite_f')
endif
call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
@ -1222,7 +1220,7 @@ subroutine HDF5_write_real3(loc_id,dataset,datasetName,parallel)
end subroutine HDF5_write_real3
!--------------------------------------------------------------------------------------------------
!> @brief write dataset of type pReal with 4 dimensions
!> @brief write dataset of type real with 4 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_write_real4(loc_id,dataset,datasetName,parallel)
@ -1255,7 +1253,7 @@ subroutine HDF5_write_real4(loc_id,dataset,datasetName,parallel)
if (product(totalShape) /= 0) then
call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_real4: h5dwrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_write_real4: h5dwrite_f')
endif
call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
@ -1264,7 +1262,7 @@ end subroutine HDF5_write_real4
!--------------------------------------------------------------------------------------------------
!> @brief write dataset of type pReal with 5 dimensions
!> @brief write dataset of type real with 5 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_write_real5(loc_id,dataset,datasetName,parallel)
@ -1297,7 +1295,7 @@ subroutine HDF5_write_real5(loc_id,dataset,datasetName,parallel)
if (product(totalShape) /= 0) then
call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_real5: h5dwrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_write_real5: h5dwrite_f')
endif
call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
@ -1305,7 +1303,7 @@ subroutine HDF5_write_real5(loc_id,dataset,datasetName,parallel)
end subroutine HDF5_write_real5
!--------------------------------------------------------------------------------------------------
!> @brief write dataset of type pReal with 6 dimensions
!> @brief write dataset of type real with 6 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_write_real6(loc_id,dataset,datasetName,parallel)
@ -1338,7 +1336,7 @@ subroutine HDF5_write_real6(loc_id,dataset,datasetName,parallel)
if (product(totalShape) /= 0) then
call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_real6: h5dwrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_write_real6: h5dwrite_f')
endif
call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
@ -1346,7 +1344,7 @@ subroutine HDF5_write_real6(loc_id,dataset,datasetName,parallel)
end subroutine HDF5_write_real6
!--------------------------------------------------------------------------------------------------
!> @brief write dataset of type pReal with 7 dimensions
!> @brief write dataset of type real with 7 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_write_real7(loc_id,dataset,datasetName,parallel)
@ -1379,7 +1377,7 @@ subroutine HDF5_write_real7(loc_id,dataset,datasetName,parallel)
if (product(totalShape) /= 0) then
call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(totalShape,HSIZE_T), hdferr,&
file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_real7: h5dwrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_write_real7: h5dwrite_f')
endif
call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
@ -1388,7 +1386,7 @@ end subroutine HDF5_write_real7
!--------------------------------------------------------------------------------------------------
!> @brief write dataset of type pInt with 1 dimension
!> @brief write dataset of type integer with 1 dimension
!--------------------------------------------------------------------------------------------------
subroutine HDF5_write_int1(loc_id,dataset,datasetName,parallel)
@ -1421,7 +1419,7 @@ subroutine HDF5_write_int1(loc_id,dataset,datasetName,parallel)
if (product(totalShape) /= 0) then
call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_int1: h5dwrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_write_int1: h5dwrite_f')
endif
call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
@ -1429,7 +1427,7 @@ subroutine HDF5_write_int1(loc_id,dataset,datasetName,parallel)
end subroutine HDF5_write_int1
!--------------------------------------------------------------------------------------------------
!> @brief write dataset of type pInt with 2 dimensions
!> @brief write dataset of type integer with 2 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_write_int2(loc_id,dataset,datasetName,parallel)
@ -1462,7 +1460,7 @@ subroutine HDF5_write_int2(loc_id,dataset,datasetName,parallel)
if (product(totalShape) /= 0) then
call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_int2: h5dwrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_write_int2: h5dwrite_f')
endif
call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
@ -1470,7 +1468,7 @@ subroutine HDF5_write_int2(loc_id,dataset,datasetName,parallel)
end subroutine HDF5_write_int2
!--------------------------------------------------------------------------------------------------
!> @brief write dataset of type pInt with 3 dimensions
!> @brief write dataset of type integer with 3 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_write_int3(loc_id,dataset,datasetName,parallel)
@ -1503,7 +1501,7 @@ subroutine HDF5_write_int3(loc_id,dataset,datasetName,parallel)
if (product(totalShape) /= 0) then
call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_int3: h5dwrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_write_int3: h5dwrite_f')
endif
call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
@ -1511,7 +1509,7 @@ subroutine HDF5_write_int3(loc_id,dataset,datasetName,parallel)
end subroutine HDF5_write_int3
!--------------------------------------------------------------------------------------------------
!> @brief write dataset of type pInt with 4 dimensions
!> @brief write dataset of type integer with 4 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_write_int4(loc_id,dataset,datasetName,parallel)
@ -1544,7 +1542,7 @@ subroutine HDF5_write_int4(loc_id,dataset,datasetName,parallel)
if (product(totalShape) /= 0) then
call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_int4: h5dwrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_write_int4: h5dwrite_f')
endif
call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
@ -1552,7 +1550,7 @@ subroutine HDF5_write_int4(loc_id,dataset,datasetName,parallel)
end subroutine HDF5_write_int4
!--------------------------------------------------------------------------------------------------
!> @brief write dataset of type pInt with 5 dimensions
!> @brief write dataset of type integer with 5 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_write_int5(loc_id,dataset,datasetName,parallel)
@ -1585,7 +1583,7 @@ subroutine HDF5_write_int5(loc_id,dataset,datasetName,parallel)
if (product(totalShape) /= 0) then
call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_int5: h5dwrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_write_int5: h5dwrite_f')
endif
call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
@ -1593,7 +1591,7 @@ subroutine HDF5_write_int5(loc_id,dataset,datasetName,parallel)
end subroutine HDF5_write_int5
!--------------------------------------------------------------------------------------------------
!> @brief write dataset of type pInt with 6 dimensions
!> @brief write dataset of type integer with 6 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_write_int6(loc_id,dataset,datasetName,parallel)
@ -1626,7 +1624,7 @@ subroutine HDF5_write_int6(loc_id,dataset,datasetName,parallel)
if (product(totalShape) /= 0) then
call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_int6: h5dwrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_write_int6: h5dwrite_f')
endif
call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
@ -1634,7 +1632,7 @@ subroutine HDF5_write_int6(loc_id,dataset,datasetName,parallel)
end subroutine HDF5_write_int6
!--------------------------------------------------------------------------------------------------
!> @brief write dataset of type pInt with 7 dimensions
!> @brief write dataset of type integer with 7 dimensions
!--------------------------------------------------------------------------------------------------
subroutine HDF5_write_int7(loc_id,dataset,datasetName,parallel)
@ -1667,7 +1665,7 @@ subroutine HDF5_write_int7(loc_id,dataset,datasetName,parallel)
if (product(totalShape) /= 0) then
call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER,dataset,int(totalShape,HSIZE_T), hdferr,&
file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_int7: h5dwrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_write_int7: h5dwrite_f')
endif
call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
@ -1743,7 +1741,7 @@ subroutine HDF5_write_rotation(loc_id,dataset,datasetName,parallel)
file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
call h5dwrite_f(dset_id, z_id,dataset_asArray(4,:),int(totalShape,HSIZE_T), hdferr,&
file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_rotation: h5dwrite_f')
if (hdferr < 0) call IO_error(1,ext_msg='HDF5_write_rotation: h5dwrite_f')
endif
call finalize_write(plist_id, dset_id, filespace_id, memspace_id)
@ -1768,7 +1766,7 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_
globalShape !< shape of the dataset (all processes)
integer(HID_T), intent(out) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
integer(pInt), dimension(worldsize) :: &
integer, dimension(worldsize) :: &
readSize !< contribution of all processes
integer :: ierr
integer :: hdferr
@ -1776,17 +1774,17 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_
!-------------------------------------------------------------------------------------------------
! creating a property list for transfer properties (is collective for MPI)
call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_read: h5pcreate_f')
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5pcreate_f')
!--------------------------------------------------------------------------------------------------
readSize = 0_pInt
readSize(worldrank+1) = int(localShape(ubound(localShape,1)),pInt)
readSize = 0
readSize(worldrank+1) = int(localShape(ubound(localShape,1)))
#ifdef PETSc
if (parallel) then
call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_read: h5pset_dxpl_mpio_f')
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5pset_dxpl_mpio_f')
call MPI_allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process
if (ierr /= 0) call IO_error(894_pInt,ext_msg='initialize_read: MPI_allreduce')
if (ierr /= 0) call IO_error(894,ext_msg='initialize_read: MPI_allreduce')
endif
#endif
myStart = int(0,HSIZE_T)
@ -1796,28 +1794,28 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_
!--------------------------------------------------------------------------------------------------
! create dataspace in memory (local shape)
call h5screate_simple_f(size(localShape), localShape, memspace_id, hdferr, localShape)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_read: h5screate_simple_f/memspace_id')
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5screate_simple_f/memspace_id')
!--------------------------------------------------------------------------------------------------
! creating a property list for IO and set it to collective
call h5pcreate_f(H5P_DATASET_ACCESS_F, aplist_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_read: h5pcreate_f')
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5pcreate_f')
#ifdef PETSc
call h5pset_all_coll_metadata_ops_f(aplist_id, .true., hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_read: h5pset_all_coll_metadata_ops_f')
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5pset_all_coll_metadata_ops_f')
#endif
!--------------------------------------------------------------------------------------------------
! open the dataset in the file and get the space ID
call h5dopen_f(loc_id,datasetName,dset_id,hdferr, dapl_id = aplist_id)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_read: h5dopen_f')
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5dopen_f')
call h5dget_space_f(dset_id, filespace_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_read: h5dget_space_f')
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5dget_space_f')
!--------------------------------------------------------------------------------------------------
! select a hyperslab (the portion of the current process) in the file
call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, localShape, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_read: h5sselect_hyperslab_f')
if (hdferr < 0) call IO_error(1,ext_msg='initialize_read: h5sselect_hyperslab_f')
end subroutine initialize_read
@ -1831,15 +1829,15 @@ subroutine finalize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id
integer :: hdferr
call h5pclose_f(plist_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='finalize_read: plist_id')
if (hdferr < 0) call IO_error(1,ext_msg='finalize_read: plist_id')
call h5pclose_f(aplist_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='finalize_read: aplist_id')
if (hdferr < 0) call IO_error(1,ext_msg='finalize_read: aplist_id')
call h5dclose_f(dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='finalize_read: h5dclose_f')
if (hdferr < 0) call IO_error(1,ext_msg='finalize_read: h5dclose_f')
call h5sclose_f(filespace_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='finalize_read: h5sclose_f/filespace_id')
if (hdferr < 0) call IO_error(1,ext_msg='finalize_read: h5sclose_f/filespace_id')
call h5sclose_f(memspace_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='finalize_read: h5sclose_f/memspace_id')
if (hdferr < 0) call IO_error(1,ext_msg='finalize_read: h5sclose_f/memspace_id')
end subroutine finalize_read
@ -1870,22 +1868,22 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
!-------------------------------------------------------------------------------------------------
! creating a property list for transfer properties (is collective when reading in parallel)
call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_write: h5pcreate_f')
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5pcreate_f')
#ifdef PETSc
if (parallel) then
call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_write: h5pset_dxpl_mpio_f')
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5pset_dxpl_mpio_f')
endif
#endif
!--------------------------------------------------------------------------------------------------
! determine the global data layout among all processes
writeSize = 0_pInt
writeSize(worldrank+1) = int(myShape(ubound(myShape,1)),pInt)
writeSize = 0
writeSize(worldrank+1) = int(myShape(ubound(myShape,1)))
#ifdef PETSc
if (parallel) then
call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get total output size over each process
if (ierr /= 0) call IO_error(894_pInt,ext_msg='initialize_write: MPI_allreduce')
if (ierr /= 0) call IO_error(894,ext_msg='initialize_write: MPI_allreduce')
endif
#endif
myStart = int(0,HSIZE_T)
@ -1895,17 +1893,16 @@ if (parallel) then
!--------------------------------------------------------------------------------------------------
! create dataspace in memory (local shape) and in file (global shape)
call h5screate_simple_f(size(myShape), myShape, memspace_id, hdferr, myShape)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_write: h5dopen_f')
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5dopen_f')
call h5screate_simple_f(size(totalShape), totalShape, filespace_id, hdferr, totalShape)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_write: h5dget_space_f')
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5dget_space_f')
!--------------------------------------------------------------------------------------------------
! create dataset in the file and select a hyperslab from it (the portion of the current process)
call h5dcreate_f(loc_id, trim(datasetName), datatype, filespace_id, dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_write: h5dcreate_f')
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5dcreate_f')
call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, myShape, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='initialize_write: h5sselect_hyperslab_f')
if (hdferr < 0) call IO_error(1,ext_msg='initialize_write: h5sselect_hyperslab_f')
end subroutine initialize_write
@ -1919,14 +1916,15 @@ subroutine finalize_write(plist_id, dset_id, filespace_id, memspace_id)
integer :: hdferr
call h5pclose_f(plist_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='finalize_write: plist_id')
if (hdferr < 0) call IO_error(1,ext_msg='finalize_write: plist_id')
call h5dclose_f(dset_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='finalize_write: h5dclose_f')
if (hdferr < 0) call IO_error(1,ext_msg='finalize_write: h5dclose_f')
call h5sclose_f(filespace_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='finalize_write: h5sclose_f/filespace_id')
if (hdferr < 0) call IO_error(1,ext_msg='finalize_write: h5sclose_f/filespace_id')
call h5sclose_f(memspace_id, hdferr)
if (hdferr < 0) call IO_error(1_pInt,ext_msg='finalize_write: h5sclose_f/memspace_id')
if (hdferr < 0) call IO_error(1,ext_msg='finalize_write: h5sclose_f/memspace_id')
end subroutine finalize_write
#endif
end module HDF5_Utilities

366
src/damage_local.f90
View File

@ -3,43 +3,46 @@
!> @brief material subroutine for locally evolving damage field
!--------------------------------------------------------------------------------------------------
module damage_local
use prec
use material
use numerics
use config
use prec
use material
use numerics
use config
use source_damage_isoBrittle
use source_damage_isoDuctile
use source_damage_anisoBrittle
use source_damage_anisoDuctile
implicit none
private
integer, dimension(:,:), allocatable, target, public :: &
damage_local_sizePostResult !< size of each post result output
implicit none
private
integer, dimension(:,:), allocatable, target, public :: &
damage_local_sizePostResult
character(len=64), dimension(:,:), allocatable, target, public :: &
damage_local_output
integer, dimension(:), allocatable, target, public :: &
damage_local_Noutput
character(len=64), dimension(:,:), allocatable, target, public :: &
damage_local_output !< name of each post result output
integer, dimension(:), allocatable, target, public :: &
damage_local_Noutput !< number of outputs per instance of this damage
enum, bind(c)
enumerator :: undefined_ID, &
damage_ID
end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable :: &
damage_local_outputID !< ID of each post result output
type :: tParameters
integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID
end type tParameters
type(tparameters), dimension(:), allocatable :: &
param
public :: &
damage_local_init, &
damage_local_updateState, &
damage_local_postResults
enum, bind(c)
enumerator :: &
undefined_ID, &
damage_ID
end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable :: &
damage_local_outputID !< ID of each post result output
type :: tParameters
integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID
end type tParameters
type(tparameters), dimension(:), allocatable :: &
param
public :: &
damage_local_init, &
damage_local_updateState, &
damage_local_postResults
contains
@ -49,167 +52,160 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine damage_local_init
integer :: maxNinstance,homog,instance,i
integer :: sizeState
integer :: NofMyHomog, h
integer :: maxNinstance,homog,instance,i
integer :: sizeState
integer :: NofMyHomog, h
integer(kind(undefined_ID)) :: &
outputID
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
outputID
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
character(len=65536), dimension(:), allocatable :: &
outputs
write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_local_label//' init -+>>>'
maxNinstance = count(damage_type == DAMAGE_local_ID)
if (maxNinstance == 0) return
allocate(damage_local_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
allocate(damage_local_output (maxval(homogenization_Noutput),maxNinstance))
damage_local_output = ''
allocate(damage_local_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
allocate(damage_local_Noutput (maxNinstance), source=0)
allocate(param(maxNinstance))
outputs
do h = 1, size(damage_type)
if (damage_type(h) /= DAMAGE_LOCAL_ID) cycle
associate(prm => param(damage_typeInstance(h)), &
config => config_homogenization(h))
write(6,'(/,a)') ' <<<+- damage_'//DAMAGE_local_label//' init -+>>>'
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
maxNinstance = count(damage_type == DAMAGE_local_ID)
if (maxNinstance == 0) return
allocate(damage_local_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
allocate(damage_local_output (maxval(homogenization_Noutput),maxNinstance))
damage_local_output = ''
allocate(damage_local_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
allocate(damage_local_Noutput (maxNinstance), source=0)
allocate(param(maxNinstance))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('damage')
damage_local_output(i,damage_typeInstance(h)) = outputs(i)
damage_local_Noutput(instance) = damage_local_Noutput(instance) + 1
damage_local_sizePostResult(i,damage_typeInstance(h)) = 1
prm%outputID = [prm%outputID , damage_ID]
end select
enddo
do h = 1, size(damage_type)
if (damage_type(h) /= DAMAGE_LOCAL_ID) cycle
associate(prm => param(damage_typeInstance(h)), &
config => config_homogenization(h))
outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
allocate(prm%outputID(0))
do i=1, size(outputs)
outputID = undefined_ID
select case(outputs(i))
case ('damage')
damage_local_output(i,damage_typeInstance(h)) = outputs(i)
damage_local_Noutput(instance) = damage_local_Noutput(instance) + 1
damage_local_sizePostResult(i,damage_typeInstance(h)) = 1
prm%outputID = [prm%outputID , damage_ID]
end select
enddo
homog = h
homog = h
NofMyHomog = count(material_homogenizationAt == homog)
instance = damage_typeInstance(homog)
NofMyHomog = count(material_homogenizationAt == homog)
instance = damage_typeInstance(homog)
! allocate state arrays
sizeState = 1
damageState(homog)%sizeState = sizeState
damageState(homog)%sizePostResults = sum(damage_local_sizePostResult(:,instance))
allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%subState0(sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%state (sizeState,NofMyHomog), source=damage_initialPhi(homog))
nullify(damageMapping(homog)%p)
damageMapping(homog)%p => mappingHomogenization(1,:,:)
deallocate(damage(homog)%p)
damage(homog)%p => damageState(homog)%state(1,:)
end associate
enddo
sizeState = 1
damageState(homog)%sizeState = sizeState
damageState(homog)%sizePostResults = sum(damage_local_sizePostResult(:,instance))
allocate(damageState(homog)%state0 (sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%subState0(sizeState,NofMyHomog), source=damage_initialPhi(homog))
allocate(damageState(homog)%state (sizeState,NofMyHomog), source=damage_initialPhi(homog))
nullify(damageMapping(homog)%p)
damageMapping(homog)%p => mappingHomogenization(1,:,:)
deallocate(damage(homog)%p)
damage(homog)%p => damageState(homog)%state(1,:)
end associate
enddo
end subroutine damage_local_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates local change in damage field
!--------------------------------------------------------------------------------------------------
function damage_local_updateState(subdt, ip, el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
subdt
logical, dimension(2) :: &
damage_local_updateState
integer :: &
homog, &
offset
real(pReal) :: &
phi, phiDot, dPhiDot_dPhi
homog = material_homogenizationAt(el)
offset = mappingHomogenization(1,ip,el)
phi = damageState(homog)%subState0(1,offset)
call damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el)
phi = max(residualStiffness,min(1.0_pReal,phi + subdt*phiDot))
damage_local_updateState = [ abs(phi - damageState(homog)%state(1,offset)) &
<= err_damage_tolAbs &
.or. abs(phi - damageState(homog)%state(1,offset)) &
<= err_damage_tolRel*abs(damageState(homog)%state(1,offset)), &
.true.]
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
subdt
logical, dimension(2) :: &
damage_local_updateState
integer :: &
homog, &
offset
real(pReal) :: &
phi, phiDot, dPhiDot_dPhi
homog = material_homogenizationAt(el)
offset = mappingHomogenization(1,ip,el)
phi = damageState(homog)%subState0(1,offset)
call damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el)
phi = max(residualStiffness,min(1.0_pReal,phi + subdt*phiDot))
damage_local_updateState = [ abs(phi - damageState(homog)%state(1,offset)) &
<= err_damage_tolAbs &
.or. abs(phi - damageState(homog)%state(1,offset)) &
<= err_damage_tolRel*abs(damageState(homog)%state(1,offset)), &
.true.]
damageState(homog)%state(1,offset) = phi
damageState(homog)%state(1,offset) = phi
end function damage_local_updateState
!--------------------------------------------------------------------------------------------------
!> @brief calculates homogenized local damage driving forces
!--------------------------------------------------------------------------------------------------
subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el)
use source_damage_isoBrittle, only: &
source_damage_isobrittle_getRateAndItsTangent
use source_damage_isoDuctile, only: &
source_damage_isoductile_getRateAndItsTangent
use source_damage_anisoBrittle, only: &
source_damage_anisobrittle_getRateAndItsTangent
use source_damage_anisoDuctile, only: &
source_damage_anisoductile_getRateAndItsTangent
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi
integer :: &
phase, &
grain, &
source, &
constituent
real(pReal) :: &
phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi
integer :: &
phase, &
grain, &
source, &
constituent
real(pReal) :: &
phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi
phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
phase = phaseAt(grain,ip,el)
constituent = phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_damage_isoBrittle_ID)
call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
phase = phaseAt(grain,ip,el)
constituent = phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_damage_isoBrittle_ID)
call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_isoDuctile_ID)
call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_isoDuctile_ID)
call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoBrittle_ID)
call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoBrittle_ID)
call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoDuctile_ID)
call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoDuctile_ID)
call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case default
localphiDot = 0.0_pReal
dLocalphiDot_dPhi = 0.0_pReal
case default
localphiDot = 0.0_pReal
dLocalphiDot_dPhi = 0.0_pReal
end select
phiDot = phiDot + localphiDot
dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi
enddo
enddo
phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal)
end select
phiDot = phiDot + localphiDot
dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi
enddo
enddo
phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal)
end subroutine damage_local_getSourceAndItsTangent
@ -219,31 +215,31 @@ end subroutine damage_local_getSourceAndItsTangent
!--------------------------------------------------------------------------------------------------
function damage_local_postResults(ip,el)
integer, intent(in) :: &
ip, & !< integration point
el !< element
real(pReal), dimension(sum(damage_local_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: &
damage_local_postResults
integer, intent(in) :: &
ip, & !< integration point
el !< element
real(pReal), dimension(sum(damage_local_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: &
damage_local_postResults
integer :: &
instance, homog, offset, o, c
homog = material_homogenizationAt(el)
offset = damageMapping(homog)%p(ip,el)
instance = damage_typeInstance(homog)
associate(prm => param(instance))
c = 0
integer :: instance, homog, offset, o, c
homog = material_homogenizationAt(el)
offset = damageMapping(homog)%p(ip,el)
instance = damage_typeInstance(homog)
associate(prm => param(instance))
c = 0
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (damage_ID)
damage_local_postResults(c+1) = damage(homog)%p(offset)
c = c + 1
end select
enddo outputsLoop
end associate
outputsLoop: do o = 1,size(prm%outputID)
select case(prm%outputID(o))
case (damage_ID)
damage_local_postResults(c+1) = damage(homog)%p(offset)
c = c + 1
end select
enddo outputsLoop
end associate
end function damage_local_postResults
end module damage_local

23
src/damage_nonlocal.f90
View File

@ -19,26 +19,25 @@ module damage_nonlocal
implicit none
private
integer, dimension(:,:), allocatable, target, public :: &
damage_nonlocal_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
damage_nonlocal_output !< name of each post result output
integer, dimension(:), allocatable, target, public :: &
damage_nonlocal_Noutput !< number of outputs per instance of this damage
integer, dimension(:,:), allocatable, target, public :: &
damage_nonlocal_sizePostResult
character(len=64), dimension(:,:), allocatable, target, public :: &
damage_nonlocal_output
integer, dimension(:), allocatable, target, public :: &
damage_nonlocal_Noutput
enum, bind(c)
enumerator :: undefined_ID, &
damage_ID
enumerator :: &
undefined_ID, &
damage_ID
end enum
type :: tParameters
integer(kind(undefined_ID)), dimension(:), allocatable :: &
integer(kind(undefined_ID)), dimension(:), allocatable :: &
outputID
end type tParameters
type(tparameters), dimension(:), allocatable :: &
type(tparameters), dimension(:), allocatable :: &
param
public :: &

10
src/future.f90
View File

@ -3,7 +3,11 @@
!> @brief New fortran functions for compiler versions that do not support them
!--------------------------------------------------------------------------------------------------
module future
use prec
implicit none
public
contains
#if defined(__GFORTRAN__) || __INTEL_COMPILER < 1800
@ -11,6 +15,7 @@ contains
!> @brief substitute for the findloc intrinsic (only for integer, dimension(:) at the moment)
!--------------------------------------------------------------------------------------------------
function findloc(a,v)
integer, intent(in), dimension(:) :: a
integer, intent(in) :: v
integer :: i,j
@ -29,13 +34,10 @@ end function findloc
#if defined(__PGI)
!--------------------------------------------------------------------------------------------------
!> @brief substitute for the norm2 intrinsic (only for real,dimension(3) at the moment)
!> @brief substitute for the norm2 intrinsic (only for real, dimension(3) at the moment)
!--------------------------------------------------------------------------------------------------
real(pReal) pure function norm2(v)
use prec, only: &
pReal
implicit none
real(pReal), intent(in), dimension(3) :: v
norm2 = sqrt(sum(v**2))

17
src/homogenization.f90
View File

@ -23,16 +23,15 @@ module homogenization
use damage_none
use damage_local
use damage_nonlocal
#if defined(PETSc) || defined(DAMASK_HDF5)
use results
use HDF5_utilities
#endif
implicit none
private
!--------------------------------------------------------------------------------------------------
! General variables for the homogenization at a material point
implicit none
private
real(pReal), dimension(:,:,:,:), allocatable, public :: &
real(pReal), dimension(:,:,:,:), allocatable, public :: &
materialpoint_F0, & !< def grad of IP at start of FE increment
materialpoint_F, & !< def grad of IP to be reached at end of FE increment
materialpoint_P !< first P--K stress of IP
@ -45,17 +44,17 @@ module homogenization
thermal_maxSizePostResults, &
damage_maxSizePostResults
real(pReal), dimension(:,:,:,:), allocatable, private :: &
real(pReal), dimension(:,:,:,:), allocatable :: &
materialpoint_subF0, & !< def grad of IP at beginning of homogenization increment
materialpoint_subF !< def grad of IP to be reached at end of homog inc
real(pReal), dimension(:,:), allocatable, private :: &
real(pReal), dimension(:,:), allocatable :: &
materialpoint_subFrac, &
materialpoint_subStep, &
materialpoint_subdt
logical, dimension(:,:), allocatable, private :: &
logical, dimension(:,:), allocatable :: &
materialpoint_requested, &
materialpoint_converged
logical, dimension(:,:,:), allocatable, private :: &
logical, dimension(:,:,:), allocatable :: &
materialpoint_doneAndHappy
interface

724
src/mesh/DAMASK_FEM.f90
View File

@ -8,452 +8,418 @@
!--------------------------------------------------------------------------------------------------
program DAMASK_FEM
#include <petsc/finclude/petscsys.h>
use PetscDM
use prec, only: &
pInt, &
pReal, &
tol_math_check
use DAMASK_interface, only: &
DAMASK_interface_init, &
loadCaseFile, &
getSolverJobName
use IO, only: &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_error, &
IO_lc, &
IO_intOut, &
IO_warning
use math ! need to include the whole module for FFTW
use CPFEM2
use FEsolving, only: &
restartWrite, &
restartInc
use numerics, only: &
worldrank, &
maxCutBack, &
stagItMax
use mesh, only: &
mesh_Nboundaries, &
mesh_boundaries, &
geomMesh
use FEM_Utilities, only: &
utilities_init, &
tSolutionState, &
tLoadCase, &
cutBack, &
maxFields, &
nActiveFields, &
FIELD_MECH_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID, &
FIELD_MECH_label
use FEM_mech
implicit none
use PetscDM
use prec
use DAMASK_interface
use IO
use math
use CPFEM2
use FEsolving
use numerics
use mesh
use FEM_Utilities
use FEM_mech
implicit none
!--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file
integer(pInt), allocatable, dimension(:) :: chunkPos ! this is longer than needed for geometry parsing
integer(pInt) :: &
N_def = 0_pInt !< # of rate of deformation specifiers found in load case file
character(len=65536) :: &
line
integer, allocatable, dimension(:) :: chunkPos ! this is longer than needed for geometry parsing
integer :: &
N_def = 0 !< # of rate of deformation specifiers found in load case file
character(len=65536) :: &
line
!--------------------------------------------------------------------------------------------------
! loop variables, convergence etc.
integer(pInt), parameter :: &
subStepFactor = 2_pInt !< for each substep, divide the last time increment by 2.0
real(pReal) :: &
time = 0.0_pReal, & !< elapsed time
time0 = 0.0_pReal, & !< begin of interval
timeinc = 0.0_pReal, & !< current time interval
timeIncOld = 0.0_pReal, & !< previous time interval
remainingLoadCaseTime = 0.0_pReal !< remaining time of current load case
logical :: &
guess, & !< guess along former trajectory
stagIterate
integer(pInt) :: &
i, &
errorID, &
cutBackLevel = 0_pInt, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$
stepFraction = 0_pInt !< fraction of current time interval
integer(pInt) :: &
currentLoadcase = 0_pInt, & !< current load case
currentFace = 0_pInt, &
inc, & !< current increment in current load case
totalIncsCounter = 0_pInt, & !< total # of increments
convergedCounter = 0_pInt, & !< # of converged increments
notConvergedCounter = 0_pInt, & !< # of non-converged increments
fileUnit = 0_pInt, & !< file unit for reading load case and writing results
myStat, &
statUnit = 0_pInt, & !< file unit for statistics output
lastRestartWritten = 0_pInt, & !< total increment No. at which last restart information was written
stagIter, &
component
character(len=6) :: loadcase_string
character(len=1024) :: &
incInfo
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
type(tSolutionState), allocatable, dimension(:) :: solres
PetscInt :: faceSet, currentFaceSet
PetscInt :: field, dimPlex
PetscErrorCode :: ierr
external :: &
quit
integer, parameter :: &
subStepFactor = 2 !< for each substep, divide the last time increment by 2.0
real(pReal) :: &
time = 0.0_pReal, & !< elapsed time
time0 = 0.0_pReal, & !< begin of interval
timeinc = 0.0_pReal, & !< current time interval
timeIncOld = 0.0_pReal, & !< previous time interval
remainingLoadCaseTime = 0.0_pReal !< remaining time of current load case
logical :: &
guess, & !< guess along former trajectory
stagIterate
integer :: &
i, &
errorID, &
cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$
stepFraction = 0 !< fraction of current time interval
integer :: &
currentLoadcase = 0, & !< current load case
currentFace = 0, &
inc, & !< current increment in current load case
totalIncsCounter = 0, & !< total # of increments
convergedCounter = 0, & !< # of converged increments
notConvergedCounter = 0, & !< # of non-converged increments
fileUnit = 0, & !< file unit for reading load case and writing results
myStat, &
statUnit = 0, & !< file unit for statistics output
lastRestartWritten = 0, & !< total increment No. at which last restart information was written
stagIter, &
component
character(len=6) :: loadcase_string
character(len=1024) :: &
incInfo
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
type(tSolutionState), allocatable, dimension(:) :: solres
PetscInt :: faceSet, currentFaceSet
PetscInt :: field, dimPlex
PetscErrorCode :: ierr
external :: &
quit
!--------------------------------------------------------------------------------------------------
! init DAMASK (all modules)
call CPFEM_initAll
write(6,'(/,a)') ' <<<+- DAMASK_FEM init -+>>>'
call CPFEM_initAll
write(6,'(/,a)') ' <<<+- DAMASK_FEM init -+>>>'
! reading basic information from load case file and allocate data structure containing load cases
call DMGetDimension(geomMesh,dimPlex,ierr); CHKERRA(ierr) !< dimension of mesh (2D or 3D)
nActiveFields = 1
allocate(solres(nActiveFields))
call DMGetDimension(geomMesh,dimPlex,ierr); CHKERRA(ierr) !< dimension of mesh (2D or 3D)
nActiveFields = 1
allocate(solres(nActiveFields))
!--------------------------------------------------------------------------------------------------
! reading basic information from load case file and allocate data structure containing load cases
open(newunit=fileunit,iostat=myStat,file=trim(loadCaseFile),action='read')
if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=trim(loadCaseFile))
do
read(fileUnit, '(A)', iostat=myStat) line
if ( myStat /= 0_pInt) exit
if (IO_isBlank(line)) cycle ! skip empty lines
chunkPos = IO_stringPos(line)
do i = 1_pInt, chunkPos(1) ! reading compulsory parameters for loadcase
select case (IO_lc(IO_stringValue(line,chunkPos,i)))
case('$loadcase')
N_def = N_def + 1_pInt
end select
enddo ! count all identifiers to allocate memory and do sanity check
enddo
allocate (loadCases(N_def))
do i = 1, size(loadCases)
allocate(loadCases(i)%fieldBC(nActiveFields))
field = 1
loadCases(i)%fieldBC(field)%ID = FIELD_MECH_ID
enddo
do i = 1, size(loadCases)
do field = 1, nActiveFields
select case (loadCases(i)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
loadCases(i)%fieldBC(field)%nComponents = dimPlex !< X, Y (, Z) displacements
allocate(loadCases(i)%fieldBC(field)%componentBC(loadCases(i)%fieldBC(field)%nComponents))
do component = 1, loadCases(i)%fieldBC(field)%nComponents
select case (component)
case (1)
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_X_ID
case (2)
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Y_ID
case (3)
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Z_ID
end select
enddo
end select
do component = 1, loadCases(i)%fieldBC(field)%nComponents
allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal)
allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.)
enddo
enddo
enddo
open(newunit=fileunit,iostat=myStat,file=trim(loadCaseFile),action='read')
if (myStat /= 0) call IO_error(100,el=myStat,ext_msg=trim(loadCaseFile))
do
read(fileUnit, '(A)', iostat=myStat) line
if ( myStat /= 0) exit
if (IO_isBlank(line)) cycle ! skip empty lines
chunkPos = IO_stringPos(line)
do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase
select case (IO_lc(IO_stringValue(line,chunkPos,i)))
case('$loadcase')
N_def = N_def + 1
end select
enddo ! count all identifiers to allocate memory and do sanity check
enddo
allocate (loadCases(N_def))
do i = 1, size(loadCases)
allocate(loadCases(i)%fieldBC(nActiveFields))
field = 1
loadCases(i)%fieldBC(field)%ID = FIELD_MECH_ID
enddo
do i = 1, size(loadCases)
do field = 1, nActiveFields
select case (loadCases(i)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
loadCases(i)%fieldBC(field)%nComponents = dimPlex !< X, Y (, Z) displacements
allocate(loadCases(i)%fieldBC(field)%componentBC(loadCases(i)%fieldBC(field)%nComponents))
do component = 1, loadCases(i)%fieldBC(field)%nComponents
select case (component)
case (1)
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_X_ID
case (2)
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Y_ID
case (3)
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Z_ID
end select
enddo
end select
do component = 1, loadCases(i)%fieldBC(field)%nComponents
allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal)
allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.)
enddo
enddo
enddo
!--------------------------------------------------------------------------------------------------
! reading the load case and assign values to the allocated data structure
rewind(fileUnit)
do
read(fileUnit, '(A)', iostat=myStat) line
if ( myStat /= 0_pInt) exit
if (IO_isBlank(line)) cycle ! skip empty lines
chunkPos = IO_stringPos(line)
do i = 1_pInt, chunkPos(1)
select case (IO_lc(IO_stringValue(line,chunkPos,i)))
rewind(fileUnit)
do
read(fileUnit, '(A)', iostat=myStat) line
if ( myStat /= 0) exit
if (IO_isBlank(line)) cycle ! skip empty lines
chunkPos = IO_stringPos(line)
do i = 1, chunkPos(1)
select case (IO_lc(IO_stringValue(line,chunkPos,i)))
!--------------------------------------------------------------------------------------------------
! loadcase information
case('$loadcase')
currentLoadCase = IO_intValue(line,chunkPos,i+1_pInt)
case('face')
currentFace = IO_intValue(line,chunkPos,i+1_pInt)
currentFaceSet = -1_pInt
do faceSet = 1, mesh_Nboundaries
if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet
enddo
if (currentFaceSet < 0_pInt) call IO_error(error_ID = errorID, ext_msg = 'invalid BC')
case('t','time','delta') ! increment time
loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1_pInt)
case('n','incs','increments','steps') ! number of increments
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1_pInt)
case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling)
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1_pInt)
loadCases(currentLoadCase)%logscale = 1_pInt
case('freq','frequency','outputfreq') ! frequency of result writings
loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1_pInt)
case('r','restart','restartwrite') ! frequency of writing restart information
loadCases(currentLoadCase)%restartfrequency = &
max(0_pInt,IO_intValue(line,chunkPos,i+1_pInt))
case('guessreset','dropguessing')
loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
case('$loadcase')
currentLoadCase = IO_intValue(line,chunkPos,i+1)
case('face')
currentFace = IO_intValue(line,chunkPos,i+1)
currentFaceSet = -1
do faceSet = 1, mesh_Nboundaries
if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet
enddo
if (currentFaceSet < 0) call IO_error(error_ID = errorID, ext_msg = 'invalid BC')
case('t','time','delta') ! increment time
loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1)
case('n','incs','increments','steps') ! number of increments
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1)
case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling)
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1)
loadCases(currentLoadCase)%logscale = 1
case('freq','frequency','outputfreq') ! frequency of result writings
loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1)
case('r','restart','restartwrite') ! frequency of writing restart information
loadCases(currentLoadCase)%restartfrequency = &
max(0,IO_intValue(line,chunkPos,i+1))
case('guessreset','dropguessing')
loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
!--------------------------------------------------------------------------------------------------
! boundary condition information
case('x') ! X displacement field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_X_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('y') ! Y displacement field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_Y_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
case('z') ! Z displacement field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_Z_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1_pInt)
endif
enddo
endif
enddo
end select
enddo; enddo
close(fileUnit)
case('x') ! X displacement field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_X_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1)
endif
enddo
endif
enddo
case('y') ! Y displacement field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_Y_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1)
endif
enddo
endif
enddo
case('z') ! Z displacement field
do field = 1, nActiveFields
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_Z_ID) then
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1)
endif
enddo
endif
enddo
end select
enddo; enddo
close(fileUnit)
!--------------------------------------------------------------------------------------------------
! consistency checks and output of load case
loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase
errorID = 0_pInt
checkLoadcases: do currentLoadCase = 1_pInt, size(loadCases)
write (loadcase_string, '(i6)' ) currentLoadCase
write(6,'(1x,a,i6)') 'load case: ', currentLoadCase
if (.not. loadCases(currentLoadCase)%followFormerTrajectory) &
write(6,'(2x,a)') 'drop guessing along trajectory'
do field = 1_pInt, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_MECH_label)
end select
do faceSet = 1_pInt, mesh_Nboundaries
do component = 1_pInt, loadCases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask(faceSet)) &
write(6,'(4x,a,i2,a,i2,a,f12.7)') 'Face ', mesh_boundaries(faceSet), &
' Component ', component, &
' Value ', loadCases(currentLoadCase)%fieldBC(field)% &
componentBC(component)%Value(faceSet)
enddo
enddo
enddo
write(6,'(2x,a,f12.6)') 'time: ', loadCases(currentLoadCase)%time
if (loadCases(currentLoadCase)%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count
write(6,'(2x,a,i5)') 'increments: ', loadCases(currentLoadCase)%incs
if (loadCases(currentLoadCase)%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency
write(6,'(2x,a,i5)') 'output frequency: ', &
loadCases(currentLoadCase)%outputfrequency
write(6,'(2x,a,i5,/)') 'restart frequency: ', &
loadCases(currentLoadCase)%restartfrequency
if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
enddo checkLoadcases
loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase
errorID = 0
checkLoadcases: do currentLoadCase = 1, size(loadCases)
write (loadcase_string, '(i6)' ) currentLoadCase
write(6,'(1x,a,i6)') 'load case: ', currentLoadCase
if (.not. loadCases(currentLoadCase)%followFormerTrajectory) &
write(6,'(2x,a)') 'drop guessing along trajectory'
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
write(6,'(2x,a)') 'Field '//trim(FIELD_MECH_label)
end select
do faceSet = 1, mesh_Nboundaries
do component = 1, loadCases(currentLoadCase)%fieldBC(field)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask(faceSet)) &
write(6,'(4x,a,i2,a,i2,a,f12.7)') 'Face ', mesh_boundaries(faceSet), &
' Component ', component, &
' Value ', loadCases(currentLoadCase)%fieldBC(field)% &
componentBC(component)%Value(faceSet)
enddo
enddo
enddo
write(6,'(2x,a,f12.6)') 'time: ', loadCases(currentLoadCase)%time
if (loadCases(currentLoadCase)%incs < 1) errorID = 835 ! non-positive incs count
write(6,'(2x,a,i5)') 'increments: ', loadCases(currentLoadCase)%incs
if (loadCases(currentLoadCase)%outputfrequency < 1) errorID = 836 ! non-positive result frequency
write(6,'(2x,a,i5)') 'output frequency: ', &
loadCases(currentLoadCase)%outputfrequency
write(6,'(2x,a,i5,/)') 'restart frequency: ', &
loadCases(currentLoadCase)%restartfrequency
if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
enddo checkLoadcases
!--------------------------------------------------------------------------------------------------
! doing initialization depending on active solvers
call Utilities_init()
do field = 1, nActiveFields
select case (loadCases(1)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
call FEM_mech_init(loadCases(1)%fieldBC(field))
end select
enddo
call Utilities_init
do field = 1, nActiveFields
select case (loadCases(1)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
call FEM_mech_init(loadCases(1)%fieldBC(field))
end select
enddo
loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases)
time0 = time ! load case start time
guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc
incLooping: do inc = 1_pInt, loadCases(currentLoadCase)%incs
totalIncsCounter = totalIncsCounter + 1_pInt
loadCaseLooping: do currentLoadCase = 1, size(loadCases)
time0 = time ! load case start time
guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc
incLooping: do inc = 1, loadCases(currentLoadCase)%incs
totalIncsCounter = totalIncsCounter + 1
!--------------------------------------------------------------------------------------------------
! forwarding time
timeIncOld = timeinc ! last timeinc that brought former inc to an end
if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale
timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
else
if (currentLoadCase == 1_pInt) then ! 1st load case of logarithmic scale
if (inc == 1_pInt) then ! 1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real( 1_pInt-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1_pInt-loadCases(1)%incs ,pReal))
endif
else ! not-1st load case of logarithmic scale
timeinc = time0 * &
( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal))&
-(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc-1_pInt ,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal)))
endif
endif
timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step
timeIncOld = timeinc ! last timeinc that brought former inc to an end
if (loadCases(currentLoadCase)%logscale == 0) then ! linear scale
timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
else
if (currentLoadCase == 1) then ! 1st load case of logarithmic scale
if (inc == 1) then ! 1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real( 1-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1-loadCases(1)%incs ,pReal))
endif
else ! not-1st load case of logarithmic scale
timeinc = time0 * &
( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal))&
-(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc-1 ,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal)))
endif
endif
timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step
skipping: if (totalIncsCounter <= restartInc) then ! not yet at restart inc?
time = time + timeinc ! just advance time, skip already performed calculation
guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference
else skipping
stepFraction = 0_pInt ! fraction scaled by stepFactor**cutLevel
skipping: if (totalIncsCounter <= restartInc) then ! not yet at restart inc?
time = time + timeinc ! just advance time, skip already performed calculation
guess = .true.
else skipping
stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time
time = time + timeinc ! forward target time
stepFraction = stepFraction + 1_pInt ! count step
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time
time = time + timeinc ! forward target time
stepFraction = stepFraction + 1 ! count step
!--------------------------------------------------------------------------------------------------
! report begin of new step
write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,a,es12.5'//&
',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//&
',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') &
'Time', time, &
's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
'-', stepFraction, '/', subStepFactor**cutBackLevel,&
' of load case ', currentLoadCase,'/',size(loadCases)
write(incInfo,&
write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,a,es12.5'//&
',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//&
',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') &
'Time', time, &
's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
'-', stepFraction, '/', subStepFactor**cutBackLevel,&
' of load case ', currentLoadCase,'/',size(loadCases)
write(incInfo,&
'(a,'//IO_intOut(totalIncsCounter)//&
',a,'//IO_intOut(sum(loadCases%incs))//&
',a,'//IO_intOut(stepFraction)//&
',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') &
'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
'-',stepFraction, '/', subStepFactor**cutBackLevel
flush(6)
flush(6)
!--------------------------------------------------------------------------------------------------
! forward fields
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
call FEM_mech_forward (&
guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
call FEM_mech_forward (&
guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select
enddo
end select
enddo
!--------------------------------------------------------------------------------------------------
! solve fields
stagIter = 0_pInt
stagIterate = .true.
do while (stagIterate)
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
solres(field) = FEM_mech_solution (&
incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
stagIter = 0
stagIterate = .true.
do while (stagIterate)
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
solres(field) = FEM_mech_solution (&
incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select
end select
if(.not. solres(field)%converged) exit ! no solution found
if(.not. solres(field)%converged) exit ! no solution found
enddo
stagIter = stagIter + 1_pInt
stagIterate = stagIter < stagItMax &
.and. all(solres(:)%converged) &
.and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration
enddo
enddo
stagIter = stagIter + 1
stagIterate = stagIter < stagItMax &
.and. all(solres(:)%converged) &
.and. .not. all(solres(:)%stagConverged) ! stationary with respect to staggered iteration
enddo
! check solution
cutBack = .False.
if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
if (cutBackLevel < maxCutBack) then ! do cut back
write(6,'(/,a)') ' cut back detected'
cutBack = .True.
stepFraction = (stepFraction - 1_pInt) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1_pInt
time = time - timeinc ! rewind time
timeinc = timeinc/2.0_pReal
else ! default behavior, exit if spectral solver does not converge
call IO_warning(850_pInt)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written
endif
else
guess = .true. ! start guessing after first converged (sub)inc
timeIncOld = timeinc
endif
if (.not. cutBack) then
if (worldrank == 0) write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
solres%converged, solres%iterationsNeeded ! write statistics about accepted solution
endif
enddo subStepLooping
cutBack = .False.
if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
if (cutBackLevel < maxCutBack) then ! do cut back
write(6,'(/,a)') ' cut back detected'
cutBack = .True.
stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1
time = time - timeinc ! rewind time
timeinc = timeinc/2.0_pReal
else ! default behavior, exit if spectral solver does not converge
call IO_warning(850)
call quit(-1*(lastRestartWritten+1)) ! quit and provide information about last restart inc written
endif
else
guess = .true. ! start guessing after first converged (sub)inc
timeIncOld = timeinc
endif
if (.not. cutBack) then
if (worldrank == 0) write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
solres%converged, solres%iterationsNeeded ! write statistics about accepted solution
endif
enddo subStepLooping
cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc
cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc
if (all(solres(:)%converged)) then
convergedCounter = convergedCounter + 1_pInt
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report converged inc
' increment ', totalIncsCounter, ' converged'
else
notConvergedCounter = notConvergedCounter + 1_pInt
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged'
endif; flush(6)
if (all(solres(:)%converged)) then
convergedCounter = convergedCounter + 1
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report converged inc
' increment ', totalIncsCounter, ' converged'
else
notConvergedCounter = notConvergedCounter + 1
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged'
endif; flush(6)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
write(6,'(1/,a)') ' ... writing results to file ......................................'
call CPFEM_results(totalIncsCounter,time)
endif
if ( loadCases(currentLoadCase)%restartFrequency > 0_pInt & ! writing of restart info requested ...
.and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! ... and at frequency of writing restart information
restartWrite = .true. ! set restart parameter for FEsolving
lastRestartWritten = inc ! QUESTION: first call to CPFEM_general will write?
endif
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency
write(6,'(1/,a)') ' ... writing results to file ......................................'
call CPFEM_results(totalIncsCounter,time)
endif
if ( loadCases(currentLoadCase)%restartFrequency > 0 & ! writing of restart info requested ...
.and. mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0) then ! ... and at frequency of writing restart information
restartWrite = .true. ! set restart parameter for FEsolving
lastRestartWritten = inc ! first call to CPFEM_general will write
endif
endif skipping
endif skipping
enddo incLooping
enddo loadCaseLooping
enddo loadCaseLooping
!--------------------------------------------------------------------------------------------------
! report summary of whole calculation
write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,'//IO_intOut(convergedCounter)//',a,'//IO_intOut(notConvergedCounter + convergedCounter)//',a,f5.1,a)') &
convergedCounter, ' out of ', &
notConvergedCounter + convergedCounter, ' (', &
real(convergedCounter, pReal)/&
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, ' %) increments converged!'
flush(6)
close(statUnit)
write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,'//IO_intOut(convergedCounter)//',a,'//IO_intOut(notConvergedCounter + convergedCounter)//',a,f5.1,a)') &
convergedCounter, ' out of ', &
notConvergedCounter + convergedCounter, ' (', &
real(convergedCounter, pReal)/&
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, ' %) increments converged!'
flush(6)
close(statUnit)
if (notConvergedCounter > 0_pInt) call quit(2_pInt) ! error if some are not converged
call quit(0_pInt) ! no complains ;)
if (notConvergedCounter > 0) call quit(2) ! error if some are not converged
call quit(0) ! no complains ;)
end program DAMASK_FEM

1054
src/mesh/FEM_mech.f90
View File

File diff suppressed because it is too large Load Diff

265
src/mesh/FEM_utilities.f90
View File

@ -6,85 +6,92 @@ module FEM_utilities
#include <petsc/finclude/petscdmplex.h>
#include <petsc/finclude/petscdmda.h>
#include <petsc/finclude/petscis.h>
use prec, only: pReal, pInt
use PETScdmplex
use PETScdmda
use PETScis
use PETScdmplex
use PETScdmda
use PETScis
use prec
use FEsolving
use homogenization
use numerics
use debug
use math
use mesh
implicit none
private
implicit none
private
!--------------------------------------------------------------------------------------------------
!
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
integer(pInt), public, parameter :: maxFields = 6_pInt
integer(pInt), public :: nActiveFields = 0_pInt
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
integer, public, parameter :: maxFields = 6
integer, public :: nActiveFields = 0
!--------------------------------------------------------------------------------------------------
! grid related information information
real(pReal), public :: wgt !< weighting factor 1/Nelems
real(pReal), public :: wgt !< weighting factor 1/Nelems
!--------------------------------------------------------------------------------------------------
! field labels information
character(len=*), parameter, public :: &
FIELD_MECH_label = 'mechanical'
enum, bind(c)
enumerator :: FIELD_UNDEFINED_ID, &
FIELD_MECH_ID
end enum
enum, bind(c)
enumerator :: COMPONENT_UNDEFINED_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID
end enum
character(len=*), parameter, public :: &
FIELD_MECH_label = 'mechanical'
enum, bind(c)
enumerator :: FIELD_UNDEFINED_ID, &
FIELD_MECH_ID
end enum
enum, bind(c)
enumerator :: COMPONENT_UNDEFINED_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID
end enum
!--------------------------------------------------------------------------------------------------
! variables controlling debugging
logical, private :: &
logical :: &
debugPETSc !< use some in debug defined options for more verbose PETSc solution
!--------------------------------------------------------------------------------------------------
! derived types
type, public :: tSolutionState !< return type of solution from FEM solver variants
logical :: converged = .true.
logical :: stagConverged = .true.
integer(pInt) :: iterationsNeeded = 0_pInt
end type tSolutionState
type, public :: tComponentBC
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
real(pReal), allocatable :: Value(:)
logical, allocatable :: Mask(:)
end type tComponentBC
type, public :: tFieldBC
integer(kind(FIELD_UNDEFINED_ID)) :: ID
integer(pInt) :: nComponents = 0_pInt
type(tComponentBC), allocatable :: componentBC(:)
end type tFieldBC
type, public :: tLoadCase
real(pReal) :: time = 0.0_pReal !< length of increment
integer(pInt) :: incs = 0_pInt, & !< number of increments
outputfrequency = 1_pInt, & !< frequency of result writes
restartfrequency = 0_pInt, & !< frequency of restart writes
logscale = 0_pInt !< linear/logarithmic time inc flag
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer(pInt), allocatable :: faceID(:)
type(tFieldBC), allocatable :: fieldBC(:)
end type tLoadCase
type, public :: tSolutionState !< return type of solution from FEM solver variants
logical :: converged = .true.
logical :: stagConverged = .true.
integer :: iterationsNeeded = 0
end type tSolutionState
public :: &
utilities_init, &
utilities_constitutiveResponse, &
utilities_projectBCValues, &
FIELD_MECH_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID
type, public :: tComponentBC
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
real(pReal), allocatable :: Value(:)
logical, allocatable :: Mask(:)
end type tComponentBC
type, public :: tFieldBC
integer(kind(FIELD_UNDEFINED_ID)) :: ID
integer :: nComponents = 0
type(tComponentBC), allocatable :: componentBC(:)
end type tFieldBC
type, public :: tLoadCase
real(pReal) :: time = 0.0_pReal !< length of increment
integer :: incs = 0, & !< number of increments
outputfrequency = 1, & !< frequency of result writes
restartfrequency = 0, & !< frequency of restart writes
logscale = 0 !< linear/logarithmic time inc flag
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer, allocatable :: faceID(:)
type(tFieldBC), allocatable :: fieldBC(:)
end type tLoadCase
public :: &
utilities_init, &
utilities_constitutiveResponse, &
utilities_projectBCValues, &
FIELD_MECH_ID, &
COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, &
COMPONENT_MECH_Z_ID
contains
@ -92,45 +99,32 @@ contains
!> @brief allocates all neccessary fields, sets debug flags
!--------------------------------------------------------------------------------------------------
subroutine utilities_init
use numerics, only: &
structOrder, &
petsc_defaultOptions, &
petsc_options
use debug, only: &
debug_level, &
debug_SPECTRAL, &
debug_SPECTRALPETSC,&
PETSCDEBUG
use math ! must use the whole module for use of FFTW
use mesh, only: &
mesh_NcpElemsGlobal, &
mesh_maxNips
character(len=1024) :: petsc_optionsPhysics
PetscErrorCode :: ierr
character(len=1024) :: petsc_optionsPhysics
PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>'
write(6,'(/,a)') ' <<<+- DAMASK_FEM_utilities init -+>>>'
!--------------------------------------------------------------------------------------------------
! set debugging parameters
debugPETSc = iand(debug_level(debug_SPECTRAL),debug_SPECTRALPETSC) /= 0
if(debugPETSc) write(6,'(3(/,a),/)') &
' Initializing PETSc with debug options: ', &
trim(PETScDebug), &
' add more using the PETSc_Options keyword in numerics.config '
flush(6)
call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr)
CHKERRQ(ierr)
if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr)
CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_defaultOptions),ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
CHKERRQ(ierr)
write(petsc_optionsPhysics,'(a,i0)') '-mechFE_petscspace_degree ' , structOrder
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsPhysics),ierr)
CHKERRQ(ierr)
wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal)
debugPETSc = iand(debug_level(debug_SPECTRAL),debug_SPECTRALPETSC) /= 0
if(debugPETSc) write(6,'(3(/,a),/)') &
' Initializing PETSc with debug options: ', &
trim(PETScDebug), &
' add more using the PETSc_Options keyword in numerics.config '
flush(6)
call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr)
CHKERRQ(ierr)
if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr)
CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_defaultOptions),ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
CHKERRQ(ierr)
write(petsc_optionsPhysics,'(a,i0)') '-mechFE_petscspace_degree ' , structOrder
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsPhysics),ierr)
CHKERRQ(ierr)
wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal)
end subroutine utilities_init
@ -139,28 +133,23 @@ end subroutine utilities_init
!> @brief calculates constitutive response
!--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
use FEsolving, only: &
restartWrite
use homogenization, only: &
materialpoint_P, &
materialpoint_stressAndItsTangent
real(pReal), intent(in) :: timeinc !< loading time
logical, intent(in) :: forwardData !< age results
real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress
PetscErrorCode :: ierr
real(pReal), intent(in) :: timeinc !< loading time
logical, intent(in) :: forwardData !< age results
real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress
PetscErrorCode :: ierr
write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field
call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field
restartWrite = .false. ! reset restartWrite status
cutBack = .false. ! reset cutBack status
P_av = sum(sum(materialpoint_P,dim=4),dim=3) * wgt ! average of P
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
restartWrite = .false. ! reset restartWrite status
cutBack = .false. ! reset cutBack status
P_av = sum(sum(materialpoint_P,dim=4),dim=3) * wgt ! average of P
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
end subroutine utilities_constitutiveResponse
@ -170,32 +159,32 @@ end subroutine utilities_constitutiveResponse
!--------------------------------------------------------------------------------------------------
subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCValue,BCDotValue,timeinc)
Vec :: localVec
PetscInt :: field, comp, nBcPoints, point, dof, numDof, numComp, offset
PetscSection :: section
IS :: bcPointsIS
PetscInt, pointer :: bcPoints(:)
PetscScalar, pointer :: localArray(:)
PetscScalar :: BCValue,BCDotValue,timeinc
PetscErrorCode :: ierr
Vec :: localVec
PetscInt :: field, comp, nBcPoints, point, dof, numDof, numComp, offset
PetscSection :: section
IS :: bcPointsIS
PetscInt, pointer :: bcPoints(:)
PetscScalar, pointer :: localArray(:)
PetscScalar :: BCValue,BCDotValue,timeinc
PetscErrorCode :: ierr
call PetscSectionGetFieldComponents(section,field,numComp,ierr); CHKERRQ(ierr)
call ISGetSize(bcPointsIS,nBcPoints,ierr); CHKERRQ(ierr)
if (nBcPoints > 0) call ISGetIndicesF90(bcPointsIS,bcPoints,ierr)
call VecGetArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)
do point = 1, nBcPoints
call PetscSectionGetFieldDof(section,bcPoints(point),field,numDof,ierr)
CHKERRQ(ierr)
call PetscSectionGetFieldOffset(section,bcPoints(point),field,offset,ierr)
CHKERRQ(ierr)
do dof = offset+comp+1, offset+numDof, numComp
localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc
enddo
enddo
call VecRestoreArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)
call VecAssemblyBegin(localVec, ierr); CHKERRQ(ierr)
call VecAssemblyEnd (localVec, ierr); CHKERRQ(ierr)
if (nBcPoints > 0) call ISRestoreIndicesF90(bcPointsIS,bcPoints,ierr)
call PetscSectionGetFieldComponents(section,field,numComp,ierr); CHKERRQ(ierr)
call ISGetSize(bcPointsIS,nBcPoints,ierr); CHKERRQ(ierr)
if (nBcPoints > 0) call ISGetIndicesF90(bcPointsIS,bcPoints,ierr)
call VecGetArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)
do point = 1, nBcPoints
call PetscSectionGetFieldDof(section,bcPoints(point),field,numDof,ierr)
CHKERRQ(ierr)
call PetscSectionGetFieldOffset(section,bcPoints(point),field,offset,ierr)
CHKERRQ(ierr)
do dof = offset+comp+1, offset+numDof, numComp
localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc
enddo
enddo
call VecRestoreArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)
call VecAssemblyBegin(localVec, ierr); CHKERRQ(ierr)
call VecAssemblyEnd (localVec, ierr); CHKERRQ(ierr)
if (nBcPoints > 0) call ISRestoreIndicesF90(bcPointsIS,bcPoints,ierr)
end subroutine utilities_projectBCValues

469
src/mesh/FEM_zoo.f90
View File

@ -3,29 +3,30 @@
!> @brief Interpolation data used by the FEM solver
!--------------------------------------------------------------------------------------------------
module FEM_Zoo
use prec, only: pReal, pInt, group_float
implicit none
private
integer(pInt), parameter, public:: &
maxOrder = 5 !< current max interpolation set at cubic (intended to be arbitrary)
real(pReal), dimension(2,3), private, parameter :: &
triangle = reshape([-1.0_pReal, -1.0_pReal, &
1.0_pReal, -1.0_pReal, &
-1.0_pReal, 1.0_pReal], shape=[2,3])
real(pReal), dimension(3,4), private, parameter :: &
tetrahedron = reshape([-1.0_pReal, -1.0_pReal, -1.0_pReal, &
1.0_pReal, -1.0_pReal, -1.0_pReal, &
-1.0_pReal, 1.0_pReal, -1.0_pReal, &
-1.0_pReal, -1.0_pReal, 1.0_pReal], shape=[3,4])
integer(pInt), dimension(3,maxOrder), public, protected :: &
FEM_Zoo_nQuadrature !< number of quadrature points for a given spatial dimension(1-3) and interpolation order(1-maxOrder)
type(group_float), dimension(3,maxOrder), public, protected :: &
FEM_Zoo_QuadratureWeights, & !< quadrature weights for each quadrature rule
FEM_Zoo_QuadraturePoints !< quadrature point coordinates (in simplical system) for each quadrature rule
public :: &
FEM_Zoo_init
use prec
implicit none
private
integer, parameter, public:: &
maxOrder = 5 !< current max interpolation set at cubic (intended to be arbitrary)
real(pReal), dimension(2,3), private, parameter :: &
triangle = reshape([-1.0_pReal, -1.0_pReal, &
1.0_pReal, -1.0_pReal, &
-1.0_pReal, 1.0_pReal], shape=[2,3])
real(pReal), dimension(3,4), private, parameter :: &
tetrahedron = reshape([-1.0_pReal, -1.0_pReal, -1.0_pReal, &
1.0_pReal, -1.0_pReal, -1.0_pReal, &
-1.0_pReal, 1.0_pReal, -1.0_pReal, &
-1.0_pReal, -1.0_pReal, 1.0_pReal], shape=[3,4])
integer, dimension(3,maxOrder), public, protected :: &
FEM_Zoo_nQuadrature !< number of quadrature points for a given spatial dimension(1-3) and interpolation order(1-maxOrder)
type(group_float), dimension(3,maxOrder), public, protected :: &
FEM_Zoo_QuadratureWeights, & !< quadrature weights for each quadrature rule
FEM_Zoo_QuadraturePoints !< quadrature point coordinates (in simplical system) for each quadrature rule
public :: &
FEM_Zoo_init
contains
@ -34,306 +35,336 @@ contains
!> @brief initializes FEM interpolation data
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_init
implicit none
write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>'
write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>'
!--------------------------------------------------------------------------------------------------
! 2D linear
FEM_Zoo_nQuadrature(2,1) = 1
allocate(FEM_Zoo_QuadratureWeights(2,1)%p(1))
allocate(FEM_Zoo_QuadraturePoints (2,1)%p(2))
FEM_Zoo_QuadratureWeights(2,1)%p(1) = 1.0_pReal
call FEM_Zoo_permutationStar3([1.0_pReal/3.0_pReal], &
FEM_Zoo_QuadraturePoints(2,1)%p(1:2))
FEM_Zoo_nQuadrature(2,1) = 1
allocate(FEM_Zoo_QuadratureWeights(2,1)%p(1))
allocate(FEM_Zoo_QuadraturePoints (2,1)%p(2))
FEM_Zoo_QuadratureWeights(2,1)%p(1) = 1.0_pReal
call FEM_Zoo_permutationStar3([1.0_pReal/3.0_pReal], &
FEM_Zoo_QuadraturePoints(2,1)%p(1:2))
!--------------------------------------------------------------------------------------------------
! 2D quadratic
FEM_Zoo_nQuadrature(2,2) = 3
allocate(FEM_Zoo_QuadratureWeights(2,2)%p(3))
allocate(FEM_Zoo_QuadraturePoints (2,2)%p(6))
FEM_Zoo_QuadratureWeights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal
call FEM_Zoo_permutationStar21([1.0_pReal/6.0_pReal], &
FEM_Zoo_QuadraturePoints(2,2)%p(1:6))
FEM_Zoo_nQuadrature(2,2) = 3
allocate(FEM_Zoo_QuadratureWeights(2,2)%p(3))
allocate(FEM_Zoo_QuadraturePoints (2,2)%p(6))
FEM_Zoo_QuadratureWeights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal
call FEM_Zoo_permutationStar21([1.0_pReal/6.0_pReal], &
FEM_Zoo_QuadraturePoints(2,2)%p(1:6))
!--------------------------------------------------------------------------------------------------
! 2D cubic
FEM_Zoo_nQuadrature(2,3) = 6
allocate(FEM_Zoo_QuadratureWeights(2,3)%p(6 ))
allocate(FEM_Zoo_QuadraturePoints (2,3)%p(12))
FEM_Zoo_QuadratureWeights(2,3)%p(1:3) = 0.22338158967801146570_pReal
call FEM_Zoo_permutationStar21([0.44594849091596488632_pReal], &
FEM_Zoo_QuadraturePoints(2,3)%p(1:6))
FEM_Zoo_QuadratureWeights(2,3)%p(4:6) = 0.10995174365532186764_pReal
call FEM_Zoo_permutationStar21([0.091576213509770743460_pReal], &
FEM_Zoo_QuadraturePoints(2,3)%p(7:12))
FEM_Zoo_nQuadrature(2,3) = 6
allocate(FEM_Zoo_QuadratureWeights(2,3)%p(6 ))
allocate(FEM_Zoo_QuadraturePoints (2,3)%p(12))
FEM_Zoo_QuadratureWeights(2,3)%p(1:3) = 0.22338158967801146570_pReal
call FEM_Zoo_permutationStar21([0.44594849091596488632_pReal], &
FEM_Zoo_QuadraturePoints(2,3)%p(1:6))
FEM_Zoo_QuadratureWeights(2,3)%p(4:6) = 0.10995174365532186764_pReal
call FEM_Zoo_permutationStar21([0.091576213509770743460_pReal], &
FEM_Zoo_QuadraturePoints(2,3)%p(7:12))
!--------------------------------------------------------------------------------------------------
! 2D quartic
FEM_Zoo_nQuadrature(2,4) = 12
allocate(FEM_Zoo_QuadratureWeights(2,4)%p(12))
allocate(FEM_Zoo_QuadraturePoints (2,4)%p(24))
FEM_Zoo_QuadratureWeights(2,4)%p(1:3) = 0.11678627572638_pReal
call FEM_Zoo_permutationStar21([0.24928674517091_pReal], &
FEM_Zoo_QuadraturePoints(2,4)%p(1:6))
FEM_Zoo_QuadratureWeights(2,4)%p(4:6) = 0.05084490637021_pReal
call FEM_Zoo_permutationStar21([0.06308901449150_pReal], &
FEM_Zoo_QuadraturePoints(2,4)%p(7:12))
FEM_Zoo_QuadratureWeights(2,4)%p(7:12) = 0.08285107561837_pReal
call FEM_Zoo_permutationStar111([0.31035245103378_pReal, 0.63650249912140_pReal], &
FEM_Zoo_QuadraturePoints(2,4)%p(13:24))
FEM_Zoo_nQuadrature(2,4) = 12
allocate(FEM_Zoo_QuadratureWeights(2,4)%p(12))
allocate(FEM_Zoo_QuadraturePoints (2,4)%p(24))
FEM_Zoo_QuadratureWeights(2,4)%p(1:3) = 0.11678627572638_pReal
call FEM_Zoo_permutationStar21([0.24928674517091_pReal], &
FEM_Zoo_QuadraturePoints(2,4)%p(1:6))
FEM_Zoo_QuadratureWeights(2,4)%p(4:6) = 0.05084490637021_pReal
call FEM_Zoo_permutationStar21([0.06308901449150_pReal], &
FEM_Zoo_QuadraturePoints(2,4)%p(7:12))
FEM_Zoo_QuadratureWeights(2,4)%p(7:12) = 0.08285107561837_pReal
call FEM_Zoo_permutationStar111([0.31035245103378_pReal, 0.63650249912140_pReal], &
FEM_Zoo_QuadraturePoints(2,4)%p(13:24))
!--------------------------------------------------------------------------------------------------
! 2D order 5
FEM_Zoo_nQuadrature(2,5) = 16
allocate(FEM_Zoo_QuadratureWeights(2,5)%p(16))
allocate(FEM_Zoo_QuadraturePoints (2,5)%p(32))
FEM_Zoo_QuadratureWeights(2,5)%p(1 ) = 0.14431560767779_pReal
call FEM_Zoo_permutationStar3([0.33333333333333_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(1:2))
FEM_Zoo_QuadratureWeights(2,5)%p(2:4) = 0.09509163426728_pReal
call FEM_Zoo_permutationStar21([0.45929258829272_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(3:8))
FEM_Zoo_QuadratureWeights(2,5)%p(5:7) = 0.10321737053472_pReal
call FEM_Zoo_permutationStar21([0.17056930775176_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(9:14))
FEM_Zoo_QuadratureWeights(2,5)%p(8:10) = 0.03245849762320_pReal
call FEM_Zoo_permutationStar21([0.05054722831703_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(15:20))
FEM_Zoo_QuadratureWeights(2,5)%p(11:16) = 0.02723031417443_pReal
call FEM_Zoo_permutationStar111([0.26311282963464_pReal, 0.72849239295540_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(21:32))
FEM_Zoo_nQuadrature(2,5) = 16
allocate(FEM_Zoo_QuadratureWeights(2,5)%p(16))
allocate(FEM_Zoo_QuadraturePoints (2,5)%p(32))
FEM_Zoo_QuadratureWeights(2,5)%p(1 ) = 0.14431560767779_pReal
call FEM_Zoo_permutationStar3([0.33333333333333_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(1:2))
FEM_Zoo_QuadratureWeights(2,5)%p(2:4) = 0.09509163426728_pReal
call FEM_Zoo_permutationStar21([0.45929258829272_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(3:8))
FEM_Zoo_QuadratureWeights(2,5)%p(5:7) = 0.10321737053472_pReal
call FEM_Zoo_permutationStar21([0.17056930775176_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(9:14))
FEM_Zoo_QuadratureWeights(2,5)%p(8:10) = 0.03245849762320_pReal
call FEM_Zoo_permutationStar21([0.05054722831703_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(15:20))
FEM_Zoo_QuadratureWeights(2,5)%p(11:16) = 0.02723031417443_pReal
call FEM_Zoo_permutationStar111([0.26311282963464_pReal, 0.72849239295540_pReal], &
FEM_Zoo_QuadraturePoints(2,5)%p(21:32))
!--------------------------------------------------------------------------------------------------
! 3D linear
FEM_Zoo_nQuadrature(3,1) = 1
allocate(FEM_Zoo_QuadratureWeights(3,1)%p(1))
allocate(FEM_Zoo_QuadraturePoints (3,1)%p(3))
FEM_Zoo_QuadratureWeights(3,1)%p(1) = 1.0_pReal
call FEM_Zoo_permutationStar4([0.25_pReal], &
FEM_Zoo_QuadraturePoints(3,1)%p(1:3))
FEM_Zoo_nQuadrature(3,1) = 1
allocate(FEM_Zoo_QuadratureWeights(3,1)%p(1))
allocate(FEM_Zoo_QuadraturePoints (3,1)%p(3))
FEM_Zoo_QuadratureWeights(3,1)%p(1) = 1.0_pReal
call FEM_Zoo_permutationStar4([0.25_pReal], &
FEM_Zoo_QuadraturePoints(3,1)%p(1:3))
!--------------------------------------------------------------------------------------------------
! 3D quadratic
FEM_Zoo_nQuadrature(3,2) = 4
allocate(FEM_Zoo_QuadratureWeights(3,2)%p(4 ))
allocate(FEM_Zoo_QuadraturePoints (3,2)%p(12))
FEM_Zoo_QuadratureWeights(3,2)%p(1:4) = 0.25_pReal
call FEM_Zoo_permutationStar31([0.13819660112501051518_pReal], &
FEM_Zoo_QuadraturePoints(3,2)%p(1:12))
FEM_Zoo_nQuadrature(3,2) = 4
allocate(FEM_Zoo_QuadratureWeights(3,2)%p(4 ))
allocate(FEM_Zoo_QuadraturePoints (3,2)%p(12))
FEM_Zoo_QuadratureWeights(3,2)%p(1:4) = 0.25_pReal
call FEM_Zoo_permutationStar31([0.13819660112501051518_pReal], &
FEM_Zoo_QuadraturePoints(3,2)%p(1:12))
!--------------------------------------------------------------------------------------------------
! 3D cubic
FEM_Zoo_nQuadrature(3,3) = 14
allocate(FEM_Zoo_QuadratureWeights(3,3)%p(14))
allocate(FEM_Zoo_QuadraturePoints (3,3)%p(42))
FEM_Zoo_QuadratureWeights(3,3)%p(1:4) = 0.073493043116361949544_pReal
call FEM_Zoo_permutationStar31([0.092735250310891226402_pReal], &
FEM_Zoo_QuadraturePoints(3,3)%p(1:12))
FEM_Zoo_QuadratureWeights(3,3)%p(5:8) = 0.11268792571801585080_pReal
call FEM_Zoo_permutationStar31([0.31088591926330060980_pReal], &
FEM_Zoo_QuadraturePoints(3,3)%p(13:24))
FEM_Zoo_QuadratureWeights(3,3)%p(9:14) = 0.042546020777081466438_pReal
call FEM_Zoo_permutationStar22([0.045503704125649649492_pReal], &
FEM_Zoo_QuadraturePoints(3,3)%p(25:42))
FEM_Zoo_nQuadrature(3,3) = 14
allocate(FEM_Zoo_QuadratureWeights(3,3)%p(14))
allocate(FEM_Zoo_QuadraturePoints (3,3)%p(42))
FEM_Zoo_QuadratureWeights(3,3)%p(1:4) = 0.073493043116361949544_pReal
call FEM_Zoo_permutationStar31([0.092735250310891226402_pReal], &
FEM_Zoo_QuadraturePoints(3,3)%p(1:12))
FEM_Zoo_QuadratureWeights(3,3)%p(5:8) = 0.11268792571801585080_pReal
call FEM_Zoo_permutationStar31([0.31088591926330060980_pReal], &
FEM_Zoo_QuadraturePoints(3,3)%p(13:24))
FEM_Zoo_QuadratureWeights(3,3)%p(9:14) = 0.042546020777081466438_pReal
call FEM_Zoo_permutationStar22([0.045503704125649649492_pReal], &
FEM_Zoo_QuadraturePoints(3,3)%p(25:42))
!--------------------------------------------------------------------------------------------------
! 3D quartic
FEM_Zoo_nQuadrature(3,4) = 35
allocate(FEM_Zoo_QuadratureWeights(3,4)%p(35))
allocate(FEM_Zoo_QuadraturePoints (3,4)%p(105))
FEM_Zoo_QuadratureWeights(3,4)%p(1:4) = 0.0021900463965388_pReal
call FEM_Zoo_permutationStar31([0.0267367755543735_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(1:12))
FEM_Zoo_QuadratureWeights(3,4)%p(5:16) = 0.0143395670177665_pReal
call FEM_Zoo_permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(13:48))
FEM_Zoo_QuadratureWeights(3,4)%p(17:22) = 0.0250305395686746_pReal
call FEM_Zoo_permutationStar22([0.4547545999844830_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(49:66))
FEM_Zoo_QuadratureWeights(3,4)%p(23:34) = 0.0479839333057554_pReal
call FEM_Zoo_permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(67:102))
FEM_Zoo_QuadratureWeights(3,4)%p(35) = 0.0931745731195340_pReal
call FEM_Zoo_permutationStar4([0.25_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(103:105))
FEM_Zoo_nQuadrature(3,4) = 35
allocate(FEM_Zoo_QuadratureWeights(3,4)%p(35))
allocate(FEM_Zoo_QuadraturePoints (3,4)%p(105))
FEM_Zoo_QuadratureWeights(3,4)%p(1:4) = 0.0021900463965388_pReal
call FEM_Zoo_permutationStar31([0.0267367755543735_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(1:12))
FEM_Zoo_QuadratureWeights(3,4)%p(5:16) = 0.0143395670177665_pReal
call FEM_Zoo_permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(13:48))
FEM_Zoo_QuadratureWeights(3,4)%p(17:22) = 0.0250305395686746_pReal
call FEM_Zoo_permutationStar22([0.4547545999844830_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(49:66))
FEM_Zoo_QuadratureWeights(3,4)%p(23:34) = 0.0479839333057554_pReal
call FEM_Zoo_permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(67:102))
FEM_Zoo_QuadratureWeights(3,4)%p(35) = 0.0931745731195340_pReal
call FEM_Zoo_permutationStar4([0.25_pReal], &
FEM_Zoo_QuadraturePoints(3,4)%p(103:105))
!--------------------------------------------------------------------------------------------------
! 3D quintic
FEM_Zoo_nQuadrature(3,5) = 56
allocate(FEM_Zoo_QuadratureWeights(3,5)%p(56))
allocate(FEM_Zoo_QuadraturePoints (3,5)%p(168))
FEM_Zoo_QuadratureWeights(3,5)%p(1:4) = 0.0010373112336140_pReal
call FEM_Zoo_permutationStar31([0.0149520651530592_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(1:12))
FEM_Zoo_QuadratureWeights(3,5)%p(5:16) = 0.0096016645399480_pReal
call FEM_Zoo_permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(13:48))
FEM_Zoo_QuadratureWeights(3,5)%p(17:28) = 0.0164493976798232_pReal
call FEM_Zoo_permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(49:84))
FEM_Zoo_QuadratureWeights(3,5)%p(29:40) = 0.0153747766513310_pReal
call FEM_Zoo_permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(85:120))
FEM_Zoo_QuadratureWeights(3,5)%p(41:52) = 0.0293520118375230_pReal
call FEM_Zoo_permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(121:156))
FEM_Zoo_QuadratureWeights(3,5)%p(53:56) = 0.0366291366405108_pReal
call FEM_Zoo_permutationStar31([0.1344783347929940_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(157:168))
FEM_Zoo_nQuadrature(3,5) = 56
allocate(FEM_Zoo_QuadratureWeights(3,5)%p(56))
allocate(FEM_Zoo_QuadraturePoints (3,5)%p(168))
FEM_Zoo_QuadratureWeights(3,5)%p(1:4) = 0.0010373112336140_pReal
call FEM_Zoo_permutationStar31([0.0149520651530592_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(1:12))
FEM_Zoo_QuadratureWeights(3,5)%p(5:16) = 0.0096016645399480_pReal
call FEM_Zoo_permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(13:48))
FEM_Zoo_QuadratureWeights(3,5)%p(17:28) = 0.0164493976798232_pReal
call FEM_Zoo_permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(49:84))
FEM_Zoo_QuadratureWeights(3,5)%p(29:40) = 0.0153747766513310_pReal
call FEM_Zoo_permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(85:120))
FEM_Zoo_QuadratureWeights(3,5)%p(41:52) = 0.0293520118375230_pReal
call FEM_Zoo_permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(121:156))
FEM_Zoo_QuadratureWeights(3,5)%p(53:56) = 0.0366291366405108_pReal
call FEM_Zoo_permutationStar31([0.1344783347929940_pReal], &
FEM_Zoo_QuadraturePoints(3,5)%p(157:168))
end subroutine FEM_Zoo_init
!--------------------------------------------------------------------------------------------------
!> @brief star 3 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar3(point,qPt)
implicit none
real(pReal) :: point(1), qPt(2,1), temp(3,1)
real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(2,1), intent(out) :: qPt
real(pReal), dimension(3,1) :: temp
temp(:,1) = [point(1), point(1), point(1)]
qPt = matmul(triangle, temp)
temp(:,1) = [point(1), point(1), point(1)]
qPt = matmul(triangle, temp)
end subroutine FEM_Zoo_permutationStar3
!--------------------------------------------------------------------------------------------------
!> @brief star 21 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar21(point,qPt)
implicit none
real(pReal) :: point(1), qPt(2,3), temp(3,3)
real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(2,3), intent(out) :: qPt
real(pReal), dimension(3,3) :: temp
temp(:,1) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1)]
temp(:,2) = [point(1), 1.0_pReal - 2.0_pReal*point(1), point(1)]
temp(:,3) = [1.0_pReal - 2.0_pReal*point(1), point(1), point(1)]
temp(:,1) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1)]
temp(:,2) = [point(1), 1.0_pReal - 2.0_pReal*point(1), point(1)]
temp(:,3) = [1.0_pReal - 2.0_pReal*point(1), point(1), point(1)]
qPt = matmul(triangle, temp)
end subroutine FEM_Zoo_permutationStar21
!--------------------------------------------------------------------------------------------------
!> @brief star 111 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar111(point,qPt)
implicit none
real(pReal) :: point(2), qPt(2,6), temp(3,6)
real(pReal), dimension(2), intent(in) :: point
real(pReal), dimension(2,6), intent(out) :: qPt
real(pReal), dimension(3,6) :: temp
temp(:,1) = [point(1), point(2), 1.0_pReal - point(1) - point(2)]
temp(:,2) = [point(1), 1.0_pReal - point(1) - point(2), point(2)]
temp(:,4) = [point(2), 1.0_pReal - point(1) - point(2), point(1)]
temp(:,5) = [1.0_pReal - point(1) - point(2), point(2), point(1)]
temp(:,6) = [1.0_pReal - point(1) - point(2), point(1), point(2)]
qPt = matmul(triangle, temp)
temp(:,1) = [point(1), point(2), 1.0_pReal - point(1) - point(2)]
temp(:,2) = [point(1), 1.0_pReal - point(1) - point(2), point(2)]
temp(:,4) = [point(2), 1.0_pReal - point(1) - point(2), point(1)]
temp(:,5) = [1.0_pReal - point(1) - point(2), point(2), point(1)]
temp(:,6) = [1.0_pReal - point(1) - point(2), point(1), point(2)]
qPt = matmul(triangle, temp)
end subroutine FEM_Zoo_permutationStar111
!--------------------------------------------------------------------------------------------------
!> @brief star 4 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar4(point,qPt)
implicit none
real(pReal) :: point(1), qPt(3,1), temp(4,1)
real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(3,1), intent(out) :: qPt
real(pReal), dimension(4,1) :: temp
temp(:,1) = [point(1), point(1), point(1), point(1)]
qPt = matmul(tetrahedron, temp)
temp(:,1) = [point(1), point(1), point(1), point(1)]
qPt = matmul(tetrahedron, temp)
end subroutine FEM_Zoo_permutationStar4
!--------------------------------------------------------------------------------------------------
!> @brief star 31 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar31(point,qPt)
implicit none
real(pReal) :: point(1), qPt(3,4), temp(4,4)
real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(3,4), intent(out) :: qPt
real(pReal), dimension(4,4) :: temp
temp(:,1) = [point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1)]
temp(:,2) = [point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1)]
temp(:,3) = [point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1)]
temp(:,4) = [1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)]
qPt = matmul(tetrahedron, temp)
temp(:,1) = [point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1)]
temp(:,2) = [point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1)]
temp(:,3) = [point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1)]
temp(:,4) = [1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)]
qPt = matmul(tetrahedron, temp)
end subroutine FEM_Zoo_permutationStar31
!--------------------------------------------------------------------------------------------------
!> @brief star 22 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar22(point,qPt)
implicit none
real(pReal) :: point(1), qPt(3,6), temp(4,6)
real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(3,6), intent(out) :: qPt
real(pReal), dimension(4,6) :: temp
temp(:,1) = [point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1)]
temp(:,2) = [point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1)]
temp(:,3) = [0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1)]
temp(:,4) = [0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1)]
temp(:,5) = [0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1)]
temp(:,6) = [point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)]
qPt = matmul(tetrahedron, temp)
temp(:,1) = [point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1)]
temp(:,2) = [point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1)]
temp(:,3) = [0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1)]
temp(:,4) = [0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1)]
temp(:,5) = [0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1)]
temp(:,6) = [point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)]
qPt = matmul(tetrahedron, temp)
end subroutine FEM_Zoo_permutationStar22
!--------------------------------------------------------------------------------------------------
!> @brief star 211 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar211(point,qPt)
implicit none
real(pReal) :: point(2), qPt(3,12), temp(4,12)
real(pReal), dimension(2), intent(in) :: point
real(pReal), dimension(3,12), intent(out) :: qPt
real(pReal), dimension(4,12) :: temp
temp(:,1 ) = [point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
temp(:,2 ) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2)]
temp(:,3 ) = [point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
temp(:,4 ) = [point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)]
temp(:,5 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2)]
temp(:,6 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1)]
temp(:,7 ) = [point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
temp(:,8 ) = [point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)]
temp(:,9 ) = [point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1)]
temp(:,10) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2)]
temp(:,11) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1)]
temp(:,12) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)]
qPt = matmul(tetrahedron, temp)
temp(:,1 ) = [point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
temp(:,2 ) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2)]
temp(:,3 ) = [point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
temp(:,4 ) = [point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)]
temp(:,5 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2)]
temp(:,6 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1)]
temp(:,7 ) = [point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
temp(:,8 ) = [point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)]
temp(:,9 ) = [point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1)]
temp(:,10) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2)]
temp(:,11) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1)]
temp(:,12) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)]
qPt = matmul(tetrahedron, temp)
end subroutine FEM_Zoo_permutationStar211
!--------------------------------------------------------------------------------------------------
!> @brief star 1111 permutation of input
!--------------------------------------------------------------------------------------------------
subroutine FEM_Zoo_permutationStar1111(point,qPt)
implicit none
real(pReal) :: point(3), qPt(3,24), temp(4,24)
real(pReal), dimension(3), intent(in) :: point
real(pReal), dimension(3,24), intent(out) :: qPt
real(pReal), dimension(4,24) :: temp
temp(:,1 ) = [point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,2 ) = [point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3)]
temp(:,3 ) = [point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,4 ) = [point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2)]
temp(:,5 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3)]
temp(:,6 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2)]
temp(:,7 ) = [point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,8 ) = [point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3)]
temp(:,9 ) = [point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,10) = [point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1)]
temp(:,11) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3)]
temp(:,12) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1)]
temp(:,13) = [point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,14) = [point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2)]
temp(:,15) = [point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,16) = [point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1)]
temp(:,17) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2)]
temp(:,18) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1)]
temp(:,19) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3)]
temp(:,20) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2)]
temp(:,21) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3)]
temp(:,22) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1)]
temp(:,23) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2)]
temp(:,24) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)]
temp(:,1 ) = [point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,2 ) = [point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3)]
temp(:,3 ) = [point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,4 ) = [point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2)]
temp(:,5 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3)]
temp(:,6 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2)]
temp(:,7 ) = [point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,8 ) = [point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3)]
temp(:,9 ) = [point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,10) = [point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1)]
temp(:,11) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3)]
temp(:,12) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1)]
temp(:,13) = [point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,14) = [point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2)]
temp(:,15) = [point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3)]
temp(:,16) = [point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1)]
temp(:,17) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2)]
temp(:,18) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1)]
temp(:,19) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3)]
temp(:,20) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2)]
temp(:,21) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3)]
temp(:,22) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1)]
temp(:,23) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2)]
temp(:,24) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)]
qPt = matmul(tetrahedron, temp)
qPt = matmul(tetrahedron, temp)
end subroutine FEM_Zoo_permutationStar1111
end module FEM_Zoo

8
src/mesh_FEM.f90
View File

@ -75,10 +75,9 @@ contains
subroutine tMesh_FEM_init(self,dimen,order,nodes)
implicit none
integer, intent(in) :: dimen
integer, intent(in) :: dimen
integer, intent(in) :: order
real(pReal), intent(in), dimension(:,:) :: nodes
real(pReal), intent(in), dimension(:,:) :: nodes
class(tMesh_FEM) :: self
if (dimen == 2) then
@ -238,9 +237,8 @@ end subroutine mesh_init
!--------------------------------------------------------------------------------------------------
pure function mesh_cellCenterCoordinates(ip,el)
implicit none
integer, intent(in) :: el, & !< element number
ip !< integration point number
ip !< integration point number
real(pReal), dimension(3) :: mesh_cellCenterCoordinates !< x,y,z coordinates of the cell center of the requested IP cell
end function mesh_cellCenterCoordinates

2
src/results.f90
View File

@ -107,7 +107,7 @@ end subroutine results_closeJobFile
!--------------------------------------------------------------------------------------------------
subroutine results_addIncrement(inc,time)
integer(pInt), intent(in) :: inc
integer, intent(in) :: inc
real(pReal), intent(in) :: time
character(len=pStringLen) :: incChar

37
src/source_damage_anisoDuctile.f90
View File

@ -6,9 +6,15 @@
!--------------------------------------------------------------------------------------------------
module source_damage_anisoDuctile
use prec
use debug
use IO
use math
use material
use config
implicit none
private
integer, dimension(:), allocatable, public, protected :: &
source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism?
source_damage_anisoDuctile_instance !< instance of damage source mechanism
@ -57,26 +63,6 @@ contains
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_init
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_error
use math, only: &
math_expand
use material, only: &
material_allocateSourceState, &
phase_source, &
phase_Nsources, &
phase_Noutput, &
SOURCE_damage_anisoDuctile_label, &
SOURCE_damage_anisoDuctile_ID, &
material_phase, &
sourceState
use config, only: &
config_phase
integer :: Ninstance,phase,instance,source,sourceOffset
integer :: NofMyPhase,p ,i
@ -181,13 +167,6 @@ end subroutine source_damage_anisoDuctile_init
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
use material, only: &
phaseAt, phasememberAt, &
plasticState, &
sourceState, &
material_homogenizationAt, &
damage, &
damageMapping
integer, intent(in) :: &
ipc, & !< component-ID of integration point
@ -222,8 +201,6 @@ end subroutine source_damage_anisoDuctile_dotState
!> @brief returns local part of nonlocal damage driving force
!--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
use material, only: &
sourceState
integer, intent(in) :: &
phase, &
@ -249,8 +226,6 @@ end subroutine source_damage_anisoDuctile_getRateAndItsTangent
!> @brief return array of local damage results
!--------------------------------------------------------------------------------------------------
function source_damage_anisoDuctile_postResults(phase, constituent)
use material, only: &
sourceState
integer, intent(in) :: &
phase, &