adopted improvements done by Arun Prakash.

mesh:
elemType identification based on lower case
Abaqus now reports more errors

IO:
new function to inquire whether inputfile contains "parts"
new function to assemble multiply included inputfile into a flat one
awareness of range generation in element numbers
error reporting
This commit is contained in:
Philip Eisenlohr 2010-07-13 10:26:07 +00:00
parent eb7830dc8f
commit 4d110126da
2 changed files with 229 additions and 56 deletions

View File

@ -5,6 +5,7 @@
CONTAINS
!---------------------------
! function IO_abaqus_assembleInputFile
! function IO_open_file(unit,relPath)
! function IO_open_inputFile(unit)
! function IO_hybridIA(Nast,ODFfileName)
@ -40,6 +41,90 @@ subroutine IO_init ()
return
endsubroutine
!********************************************************************
! AP: 12.07.10
! create a new input file for abaqus simulations
! by removing all comment lines and including "include"s
!********************************************************************
recursive function IO_abaqus_assembleInputFile(unit1,unit2) result(createSuccess)
use prec
use mpie_interface
implicit none
character(len=300) line,fname
integer(pInt), intent(in) :: unit1, unit2
logical createSuccess,fexist
integer(pInt) i
do
read(unit2,'(A300)',END=220) line
line = IO_lc(trim(line))
! call IO_lcInPlace(line)
if (line(1:8)=='*include') then
fname = trim(getSolverWorkingDirectoryName())//trim(line(9+scan(line(9:),'='):))
inquire(file=fname, exist=fexist)
if (.not.(fexist)) then
write(6,*)'ERROR: file does not exist error in IO_abaqus_assembleInputFile'
write(6,*)'filename: ', trim(fname)
createSuccess = .false.
return
endif
open(unit2+1,err=200,status='old',file=fname)
if (IO_abaqus_assembleInputFile(unit1,unit2+1)) then
createSuccess=.true.
close(unit2+1)
else
createSuccess=.false.
return
endif
else if (line(1:2) /= '**') then
write(unit1,'(A)') trim(line)
endif
enddo
220 createSuccess = .true.
return
200 createSuccess =.false.
return
end function
!***********************************************************
! check if the input file for Abaqus contains part info
!***********************************************************
function IO_abaqus_hasNoPart(unit)
use prec, only: pInt
implicit none
integer(pInt) unit
integer(pInt), parameter :: maxNchunks = 1
integer(pInt), dimension(1+2*maxNchunks) :: pos
logical IO_abaqus_hasNoPart
character(len=300) line
IO_abaqus_hasNoPart = .true.
610 FORMAT(A300)
rewind(unit)
do
read(unit,610,END=620) line
pos = IO_stringPos(line,maxNchunks)
if (IO_lc(IO_stringValue(line,pos,1)) == '*part' ) then
IO_abaqus_hasNoPart = .false.
exit
endif
enddo
620 return
endfunction
!********************************************************************
! open existing file to given unit
! path to file is relative to working directory
@ -54,17 +139,21 @@ endsubroutine
character(len=*) relPath
integer(pInt) unit
IO_open_file = .false.
open(unit,status='old',err=100,file=trim(getSolverWorkingDirectoryName())//relPath)
IO_open_file = .true.
return
100 IO_open_file = .false.
return
100 return
endfunction
!********************************************************************
! open FEM inputfile to given unit
! AP: 12.07.10
! : changed the function to open *.inp_assembly, which is basically
! the input file without comment lines and possibly assembled includes
!********************************************************************
logical function IO_open_inputFile(unit)
@ -74,12 +163,23 @@ endsubroutine
integer(pInt), intent(in) :: unit
open(unit,status='old',err=100,file=trim(getSolverWorkingDirectoryName())//&
trim(getSolverJobName())//InputFileExtension)
IO_open_inputFile = .true.
return
100 IO_open_inputFile = .false.
return
IO_open_inputFile = .false.
if (FEsolver == 'Abaqus') then
open(unit+1,status='old',err=100,&
file=trim(getSolverWorkingDirectoryName())//&
trim(getSolverJobName())//InputFileExtension)
open(unit,err=100,file=trim(getSolverWorkingDirectoryName())//&
trim(getSolverJobName())//InputFileExtension//'_assembly')
IO_open_inputFile = IO_abaqus_assembleInputFile(unit,unit+1) ! strip comments and concatenate any "include"s
close(unit+1)
else
open(unit,status='old',err=100,file=trim(getSolverWorkingDirectoryName())//&
trim(getSolverJobName())//InputFileExtension)
IO_open_inputFile = .true.
endif
100 return
endfunction
@ -97,12 +197,13 @@ endsubroutine
integer(pInt), intent(in) :: unit
character(*), intent(in) :: newExt
IO_open_jobFile = .false.
open(unit,status='replace',err=100,file=trim(getSolverWorkingDirectoryName())//&
trim(getSolverJobName())//'.'//newExt)
IO_open_jobFile = .true.
return
100 IO_open_jobFile = .false.
return
100 return
endfunction
@ -139,9 +240,11 @@ endsubroutine
function IO_hybridIA(Nast,ODFfileName)
use prec, only: pReal, pInt
implicit none
real(pReal), parameter :: pi = 3.14159265358979323846264338327950288419716939937510_pReal
real(pReal), parameter :: inRad = pi/180.0_pReal
character(len=*) ODFfileName
character(len=80) line
character(len=*), parameter :: fileFormat = '(A80)'
@ -154,9 +257,6 @@ endsubroutine
real(pReal), dimension(:,:,:), allocatable :: dV_V
real(pReal), dimension(3,Nast) :: IO_hybridIA
real(pReal), parameter :: pi = 3.14159265358979323846264338327950288419716939937510_pReal
real(pReal), parameter :: inRad = pi/180.0_pReal
if (.not. IO_open_file(999,ODFfileName)) goto 100
!--- parse header of ODF file ---
@ -732,6 +832,9 @@ endfunction
!********************************************************************
! count lines containig data up to next *keyword
! AP: changed the function to neglect comment lines between keyword definitions.
! : is not changed back to the original version since *.inp_assembly does not
! : contain any comment lines (12.07.2010)
!********************************************************************
function IO_countDataLines (unit)
@ -739,7 +842,7 @@ endfunction
implicit none
integer(pInt) IO_countDataLines,unit
integer(pInt), parameter :: maxNchunks = 64
integer(pInt), parameter :: maxNchunks = 1
integer(pInt), dimension(1+2*maxNchunks) :: pos
character(len=300) line,tmp
@ -748,10 +851,10 @@ endfunction
read(unit,'(A300)',end=100) line
pos = IO_stringPos(line,maxNchunks)
tmp = IO_lc(IO_stringValue(line,pos,1))
if (tmp(1:1) == '*' ) then ! found keyword
if (tmp(1:1) == '*' .and. tmp(2:2) /= '*') then ! found keyword
exit
else
IO_countDataLines = IO_countDataLines + 1_pInt
if (tmp(2:2) /= '*') IO_countDataLines = IO_countDataLines + 1_pInt
endif
enddo
100 backspace(unit)
@ -837,8 +940,10 @@ endfunction
integer(pInt) :: lookupMaxN
integer(pInt), dimension(:,:) :: lookupMap
character(len=300) line
logical rangeGeneration
IO_continousIntValues = 0
rangeGeneration = .false.
select case (FEsolver)
case ('Marc')
@ -880,6 +985,14 @@ endfunction
backspace(unit)
enddo
! check if the element values in the elset are auto generated
backspace(unit)
read(unit,'(A300)',end=100) line
pos = IO_stringPos(line,maxNchunks)
do i = 1,pos(1)
if (IO_lc(IO_stringValue(line,pos,i)) == 'generate') rangeGeneration = .true.
enddo
do l = 1,count
read(unit,'(A300)',end=100) line
pos = IO_stringPos(line,maxNchunks)
@ -894,11 +1007,17 @@ endfunction
endif
enddo
enddo
else ! assuming range generation
else if (rangeGeneration) then ! range generation
do i = IO_intValue(line,pos,1),IO_intValue(line,pos,2),max(1,IO_intValue(line,pos,3))
IO_continousIntValues(1) = IO_continousIntValues(1) + 1
IO_continousIntValues(1+IO_continousIntValues(1)) = i
enddo
else ! read individual elem nums
do i = 1,pos(1)
! write(*,*)'IO_CIV-int',IO_intValue(line,pos,i)
IO_continousIntValues(1) = IO_continousIntValues(1) + 1
IO_continousIntValues(1+IO_continousIntValues(1)) = IO_intValue(line,pos,i)
enddo
endif
enddo
@ -923,7 +1042,7 @@ endfunction
integer(pInt), intent(in) :: ID
integer(pInt), optional, intent(in) :: e,i,g
character(len=*), optional, intent(in) :: ext_msg
character(len=80) msg
character(len=120) msg
select case (ID)
case (40)
@ -1098,6 +1217,32 @@ endfunction
case (700)
msg = 'Singular matrix in stress iteration'
! Error messages related to parsing of Abaqus input file
case (900)
msg = 'PARSE ERROR: Improper definition of nodes in input file (Nnodes < 2)'
case (901)
msg = 'PARSE ERROR: No Elements defined in input file (Nelems = 0)'
case (902)
msg = 'PARSE ERROR: No Element sets defined in input file (Atleast one *Elset must exist)'
case (903)
msg = 'PARSE ERROR: No Materials defined in input file (Look into section assigments)'
case (904)
msg = 'PARSE ERROR: No elements could be assigned for Elset: '//ext_msg
case (905)
msg = 'PARSE ERROR: Error in mesh_abaqus_map_materials'
case (906)
msg = 'PARSE ERROR: Error in mesh_abaqus_count_cpElements'
case (907)
msg = 'PARSE ERROR: Incorrect size of mesh_mapFEtoCPelem in mesh_abaqus_map_elements; Size cannot be zero'
case (908)
msg = 'PARSE ERROR: Incorrect size of mesh_mapFEtoCPnode in mesh_abaqus_map_nodes; Size cannot be zero'
case (909)
msg = 'PARSE ERROR: Incorrect size of mesh_node in mesh_abaqus_build_nodes; must be equal to mesh_Nnodes'
case(910)
msg = 'PARSE ERROR: Incorrect element type mapping in '//ext_msg
case default
msg = 'Unknown error number...'
end select

View File

@ -65,6 +65,8 @@
integer(pInt), dimension(:,:,:), allocatable :: FE_ipNeighbor
integer(pInt), dimension(:,:,:), allocatable :: FE_subNodeParent
integer(pInt), dimension(:,:,:,:), allocatable :: FE_subNodeOnIPFace
logical :: noPart ! for cases where the ABAQUS input file does not use part/assembly information
integer(pInt) :: hypoelasticTableStyle
integer(pInt) :: initialcondTableStyle
@ -213,7 +215,7 @@
use mpie_interface
use prec, only: pInt
use IO, only: IO_error,IO_open_InputFile
use IO, only: IO_error,IO_open_InputFile,IO_abaqus_hasNoPart
use FEsolving, only: parallelExecution, FEsolving_execElem, FEsolving_execIP, calcMode, lastMode
implicit none
@ -228,7 +230,6 @@
call mesh_build_FEdata() ! --- get properties of the different types of elements
if (IO_open_inputFile(fileUnit)) then ! --- parse info from input file...
select case (FEsolver)
@ -253,6 +254,7 @@
call mesh_marc_count_cpSizes(fileunit)
call mesh_marc_build_elements(fileUnit)
case ('Abaqus')
noPart = IO_abaqus_hasNoPart(fileUnit)
call mesh_abaqus_count_nodesAndElements(fileUnit)
call mesh_abaqus_count_elementSets(fileUnit)
call mesh_abaqus_count_materials(fileUnit)
@ -273,7 +275,7 @@
call mesh_build_ipVolumes()
call mesh_build_ipAreas()
call mesh_tell_statistics()
parallelExecution = (parallelExecution .and. (mesh_Nelems == mesh_NcpElems)) ! plus potential killer from non-local constitutive
else
call IO_error(101) ! cannot open input file
@ -296,36 +298,38 @@
! mapping of FE element types to internal representation
!***********************************************************
function FE_mapElemtype(what)
use IO, only: IO_lc
implicit none
character(len=*), intent(in) :: what
integer(pInt) FE_mapElemtype
select case (what)
select case (IO_lc(what))
case ( '7', &
'C3D8')
'c3d8')
FE_mapElemtype = 1 ! Three-dimensional Arbitrarily Distorted Brick
case ('134', &
'C3D4')
'c3d4')
FE_mapElemtype = 2 ! Three-dimensional Four-node Tetrahedron
case ( '11', &
'CPE4')
'cpe4')
FE_mapElemtype = 3 ! Arbitrary Quadrilateral Plane-strain
case ( '27', &
'CPE8')
'cpe8')
FE_mapElemtype = 4 ! Plane Strain, Eight-node Distorted Quadrilateral
case ('157')
FE_mapElemtype = 5 ! Three-dimensional, Low-order, Tetrahedron, Herrmann Formulations
case ('136', &
'C3D6')
'c3d6')
FE_mapElemtype = 6 ! Three-dimensional Arbitrarily Distorted Pentahedral
case ( '21', &
'C3D20')
'c3d20')
FE_mapElemtype = 7 ! Three-dimensional Arbitrarily Distorted quadratic hexahedral
case ( '117', &
'123', &
'C3D8R')
'c3d8r')
FE_mapElemtype = 8 ! Three-dimensional Arbitrarily Distorted linear hexahedral with reduced integration
case default
FE_mapElemtype = 0 ! unknown element --> should raise an error upstream..!
@ -1256,7 +1260,7 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. &
IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false.
if (inPart) then
if (inPart .or. noPart) then
select case ( IO_lc(IO_stringValue(line,pos,1)))
case('*node')
if( &
@ -1277,8 +1281,11 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
endselect
endif
enddo
620 if (mesh_Nnodes < 2) call IO_error(900)
if (mesh_Nelems == 0) call IO_error(901)
620 return
return
endsubroutine
@ -1355,12 +1362,14 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. &
IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false.
if ( inPart .and. IO_lc(IO_stringValue(line,pos,1)) == '*elset' ) &
if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,pos,1)) == '*elset' ) &
mesh_NelemSets = mesh_NelemSets + 1_pInt
enddo
620 return
620 continue
if (mesh_NelemSets == 0) call IO_error(902)
return
endsubroutine
@ -1395,14 +1404,15 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. &
IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false.
if ( inPart .and. &
if ( (inPart .or. noPart) .and. &
IO_lc(IO_StringValue(line,pos,1)) == '*solid' .and. &
IO_lc(IO_StringValue(line,pos,2)) == 'section' ) &
mesh_Nmaterials = mesh_Nmaterials + 1_pInt
enddo
620 return
620 if (mesh_Nmaterials == 0) call IO_error(903)
return
endsubroutine
@ -1510,8 +1520,9 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
endselect
enddo
620 return
620 if (mesh_NcpElems == 0) call IO_error(906)
return
endsubroutine
@ -1572,7 +1583,7 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
integer(pInt), dimension (1+2*maxNchunks) :: pos
character(len=300) line
integer(pInt) unit,elemSet
integer(pInt) unit,elemSet,i
logical inPart
allocate (mesh_nameElemSet(mesh_NelemSets)) ; mesh_nameElemSet = ''
@ -1590,15 +1601,20 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. &
IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false.
if ( inPart .and. IO_lc(IO_stringValue(line,pos,1)) == '*elset' ) then
if ( (inPart .or. noPart) .and. IO_lc(IO_stringValue(line,pos,1)) == '*elset' ) then
elemSet = elemSet + 1_pInt
mesh_nameElemSet(elemSet) = IO_extractValue(IO_lc(IO_stringValue(line,pos,2)),'elset')
mesh_mapElemSet(:,elemSet) = IO_continousIntValues(unit,mesh_Nelems,mesh_nameElemSet,mesh_mapElemSet,elemSet-1)
endif
enddo
640 return
640 do i = 1,elemSet
! write(6,*)'elemSetName: ',mesh_nameElemSet(i)
! write(6,*)'elems in Elset',mesh_mapElemSet(:,i)
if (mesh_mapElemSet(1,i) == 0) call IO_error(ID=904,ext_msg=mesh_nameElemSet(i))
enddo
return
endsubroutine
@ -1636,7 +1652,7 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. &
IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false.
if ( inPart .and. &
if ( (inPart .or. noPart) .and. &
IO_lc(IO_StringValue(line,pos,1)) == '*solid' .and. &
IO_lc(IO_StringValue(line,pos,2)) == 'section' ) then
@ -1658,8 +1674,14 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
endif
enddo
620 return
620 if (count==0) call IO_error(905)
do i=1,count
! write(6,*)'name of materials: ',i,mesh_nameMaterial(i)
! write(6,*)'name of elemSets: ',i,mesh_mapMaterial(i)
if (mesh_nameMaterial(i)=='' .or. mesh_mapMaterial(i)=='') call IO_error(905)
enddo
return
endsubroutine
@ -1682,7 +1704,7 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
mesh_mapFEtoCPnode(:,i) = i
return
endsubroutine
@ -1770,7 +1792,7 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. &
IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false.
if( inPart .and. &
if( (inPart .or. noPart) .and. &
IO_lc(IO_stringValue(line,pos,1)) == '*node' .and. &
( IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. &
IO_lc(IO_stringValue(line,pos,2)) /= 'print' .and. &
@ -1793,6 +1815,7 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
650 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2))
if (size(mesh_mapFEtoCPnode) == 0) call IO_error(908)
return
endsubroutine
@ -1927,6 +1950,8 @@ FE_ipNeighbor(:FE_NipNeighbors(8),:FE_Nips(8),8) = & ! element 117
660 call qsort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2)) ! should be mesh_NcpElems
if (size(mesh_mapFEtoCPelem) < 2) call IO_error(907)
return
endsubroutine
@ -2040,13 +2065,14 @@ subroutine mesh_marc_count_cpSizes (unit)
if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. &
IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false.
if( inPart .and. &
if( (inPart .or. noPart) .and. &
IO_lc(IO_stringValue(line,pos,1)) == '*element' .and. &
( IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. &
IO_lc(IO_stringValue(line,pos,2)) /= 'matrix' .and. &
IO_lc(IO_stringValue(line,pos,2)) /= 'response' ) &
) then
t = FE_mapElemtype(IO_extractValue(IO_stringValue(line,pos,2),'type')) ! remember elem type
t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,pos,2)),'type')) ! remember elem type
if (t==0) call IO_error(ID=910,ext_msg='mesh_abaqus_count_cpSizes')
count = IO_countDataLines(unit)
do i = 1,count
backspace(unit)
@ -2232,7 +2258,7 @@ subroutine mesh_marc_count_cpSizes (unit)
if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. &
IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false.
if( inPart .and. &
if( (inPart .or. noPart) .and. &
IO_lc(IO_stringValue(line,pos,1)) == '*node' .and. &
( IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. &
IO_lc(IO_stringValue(line,pos,2)) /= 'print' .and. &
@ -2252,7 +2278,8 @@ subroutine mesh_marc_count_cpSizes (unit)
endif
enddo
670 return
670 if (size(mesh_node,2) /= mesh_Nnodes) call IO_error(909)
return
endsubroutine
@ -2328,7 +2355,7 @@ subroutine mesh_marc_count_cpSizes (unit)
mesh_element ( 2,e) = FE_mapElemtype('C3D8R') ! elem type
mesh_element ( 3,e) = homog ! homogenization
mesh_element ( 4,e) = IO_IntValue(line,pos,1) ! microstructure
mesh_element ( 5,e) = e + (e-1)/a + (e-1)/a/b*(a+1) ! base node
mesh_element ( 5,e) = e + (e-1)/a + (e-1)/a/b*(a+1) ! base node
mesh_element ( 6,e) = mesh_element ( 5,e) + 1
mesh_element ( 7,e) = mesh_element ( 5,e) + (a+1) + 1
mesh_element ( 8,e) = mesh_element ( 5,e) + (a+1)
@ -2465,13 +2492,14 @@ subroutine mesh_marc_count_cpSizes (unit)
if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. &
IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false.
if( inPart .and. &
if( (inPart .or. noPart) .and. &
IO_lc(IO_stringValue(line,pos,1)) == '*element' .and. &
( IO_lc(IO_stringValue(line,pos,2)) /= 'output' .and. &
IO_lc(IO_stringValue(line,pos,2)) /= 'matrix' .and. &
IO_lc(IO_stringValue(line,pos,2)) /= 'response' ) &
) then
t = FE_mapElemtype(IO_extractValue(IO_stringValue(line,pos,2),'type')) ! remember elem type
t = FE_mapElemtype(IO_extractValue(IO_lc(IO_stringValue(line,pos,2)),'type')) ! remember elem type
if (t==0) call IO_error(ID=910,ext_msg='mesh_abaqus_build_elements')
count = IO_countDataLines(unit)
do i = 1,count
backspace(unit)