diff --git a/code/IO.f90 b/code/IO.f90 index 0038d767d..85b4014ef 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -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 diff --git a/code/mesh.f90 b/code/mesh.f90 index e9c623e87..88572043c 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -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)