From 03f244995ed3f79240ae2c6e72d3da62af804e25 Mon Sep 17 00:00:00 2001 From: Franz Roters Date: Wed, 25 Apr 2007 07:33:24 +0000 Subject: [PATCH] corrected mesh_build_elements (unit) so that state variables are read correctly from input file --- trunk/mesh.f90 | 473 +++++++++++++++++++++++++------------------------ 1 file changed, 239 insertions(+), 234 deletions(-) diff --git a/trunk/mesh.f90 b/trunk/mesh.f90 index 94f93daa2..420705953 100644 --- a/trunk/mesh.f90 +++ b/trunk/mesh.f90 @@ -8,11 +8,11 @@ ! --------------------------- ! _Nelems : total number of elements in mesh -! _NcpElems : total number of CP elements in mesh +! _NcpElems : total number of CP elements in mesh ! _Nnodes : total number of nodes in mesh ! _maxNnodes : max number of nodes in any CP element -! _maxNips : max number of IPs in any CP element -! _maxNipNeighbors : max number of IP neighbors in any CP element +! _maxNips : max number of IPs in any CP element +! _maxNipNeighbors : max number of IP neighbors in any CP element ! _maxNsharedElems : max number of CP elements sharing a node ! ! _element : FEid, type, material, texture, node indices @@ -40,7 +40,7 @@ ! _ipNeighborhood : 6 or less neighboring IPs as [element_num, IP_index] ! order is +x,-x,+y,-y,+z,-z but meaning strongly depends on Elemtype ! --------------------------- - integer(pInt) mesh_Nelems,mesh_NcpElems + integer(pInt) mesh_Nelems,mesh_NcpElems integer(pInt) mesh_Nnodes,mesh_maxNnodes,mesh_maxNips,mesh_maxNipNeighbors,mesh_maxNsharedElems integer(pInt), dimension(:,:), allocatable, target :: mesh_mapFEtoCPelem,mesh_mapFEtoCPnode integer(pInt), dimension(:,:), allocatable :: mesh_element, mesh_sharedElem @@ -55,65 +55,65 @@ integer(pInt), parameter :: FE_maxNfaces = 6 integer(pInt), dimension(200):: FE_mapElemtype integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nnodes = & - (/8, & ! element 7 - 4 & ! element 134 + (/8, & ! element 7 + 4 & ! element 134 /) integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nips = & - (/8, & ! element 7 - 1 & ! element 134 - /) + (/8, & ! element 7 + 1 & ! element 134 + /) integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NipNeighbors = & - (/6, & ! element 7 - 4 & ! element 134 - /) + (/6, & ! element 7 + 4 & ! element 134 + /) integer(pInt), dimension(FE_maxNfaces,FE_Nelemtypes), parameter :: FE_NfaceNodes = & reshape((/& - 4,4,4,4,4,4, & ! element 7 - 3,3,3,3,0,0 & ! element 134 + 4,4,4,4,4,4, & ! element 7 + 3,3,3,3,0,0 & ! element 134 /),(/FE_maxNfaces,FE_Nelemtypes/)) integer(pInt), dimension(FE_maxNips,FE_Nelemtypes), parameter :: FE_nodeAtIP = & reshape((/& - 1,2,4,3,5,6,8,7, & ! element 7 - 1,0,0,0,0,0,0,0 & ! element 134 + 1,2,4,3,5,6,8,7, & ! element 7 + 1,0,0,0,0,0,0,0 & ! element 134 /),(/FE_maxNips,FE_Nelemtypes/)) integer(pInt), dimension(FE_maxNnodes,FE_Nelemtypes), parameter :: FE_ipAtNode = & reshape((/& - 1,2,4,3,5,6,8,7, & ! element 7 - 1,1,1,1,0,0,0,0 & ! element 134 + 1,2,4,3,5,6,8,7, & ! element 7 + 1,1,1,1,0,0,0,0 & ! element 134 /),(/FE_maxNnodes,FE_Nelemtypes/)) integer(pInt), dimension(FE_maxNfaceNodes,FE_maxNfaces,FE_Nelemtypes), parameter :: FE_nodeOnFace = & reshape((/& - 1,2,3,4 , & ! element 7 - 2,1,5,6 , & - 3,2,6,7 , & - 3,4,8,7 , & - 4,1,5,8 , & - 8,7,6,5 , & - 1,2,3,0 , & ! element 134 - 1,4,2,0 , & - 2,3,4,0 , & - 1,3,4,0 , & - 0,0,0,0 , & - 0,0,0,0 & + 1,2,3,4 , & ! element 7 + 2,1,5,6 , & + 3,2,6,7 , & + 3,4,8,7 , & + 4,1,5,8 , & + 8,7,6,5 , & + 1,2,3,0 , & ! element 134 + 1,4,2,0 , & + 2,3,4,0 , & + 1,3,4,0 , & + 0,0,0,0 , & + 0,0,0,0 & /),(/FE_maxNfaceNodes,FE_maxNfaces,FE_Nelemtypes/)) integer(pInt), dimension(FE_maxNneighbors,FE_maxNips,FE_Nelemtypes), parameter :: FE_ipNeighbor = & reshape((/& - 2,-5, 3,-2, 5,-1 , & ! element 7 - -3, 1, 4,-2, 6,-1 , & - 4,-5,-4, 1, 7,-1 , & - -3, 3,-4, 2, 8,-1 , & - 6,-5, 7,-2,-6, 1 , & - -3, 5, 8,-2,-6, 2 , & - 8,-5,-4, 5,-6, 3 , & - -3, 7,-4, 6,-6, 4 , & - -1,-2,-3,-4, 0, 0 , & ! element 134 - 0, 0, 0, 0, 0, 0 , & - 0, 0, 0, 0, 0, 0 , & - 0, 0, 0, 0, 0, 0 , & - 0, 0, 0, 0, 0, 0 , & - 0, 0, 0, 0, 0, 0 , & - 0, 0, 0, 0, 0, 0 , & - 0, 0, 0, 0, 0, 0 & + 2,-5, 3,-2, 5,-1 , & ! element 7 + -3, 1, 4,-2, 6,-1 , & + 4,-5,-4, 1, 7,-1 , & + -3, 3,-4, 2, 8,-1 , & + 6,-5, 7,-2,-6, 1 , & + -3, 5, 8,-2,-6, 2 , & + 8,-5,-4, 5,-6, 3 , & + -3, 7,-4, 6,-6, 4 , & + -1,-2,-3,-4, 0, 0 , & ! element 134 + 0, 0, 0, 0, 0, 0 , & + 0, 0, 0, 0, 0, 0 , & + 0, 0, 0, 0, 0, 0 , & + 0, 0, 0, 0, 0, 0 , & + 0, 0, 0, 0, 0, 0 , & + 0, 0, 0, 0, 0, 0 , & + 0, 0, 0, 0, 0, 0 & /),(/FE_maxNneighbors,FE_maxNips,FE_Nelemtypes/)) CONTAINS @@ -125,9 +125,9 @@ ! --------------------------- -!*********************************************************** -! initialization -!*********************************************************** +!*********************************************************** +! initialization +!*********************************************************** SUBROUTINE mesh_init () use prec, only: pInt @@ -135,23 +135,23 @@ implicit none integer(pInt), parameter :: fileUnit = 222 - - mesh_Nelems = 0_pInt - mesh_NcpElems = 0_pInt - mesh_Nnodes = 0_pInt - mesh_maxNips = 0_pInt - mesh_maxNnodes = 0_pInt - mesh_maxNsharedElems = 0_pInt - - FE_mapElemtype = 1 ! MISSING this should be zero... - FE_mapElemtype( 7) = 1 - FE_mapElemtype(134) = 2 - + + mesh_Nelems = 0_pInt + mesh_NcpElems = 0_pInt + mesh_Nnodes = 0_pInt + mesh_maxNips = 0_pInt + mesh_maxNnodes = 0_pInt + mesh_maxNsharedElems = 0_pInt + + FE_mapElemtype = 1 ! MISSING this should be zero... + FE_mapElemtype( 7) = 1 + FE_mapElemtype(134) = 2 + ! call to various subrountes to parse the stuff from the input file... if (IO_open_inputFile(fileUnit)) then call mesh_get_meshDimensions(fileUnit) - call mesh_build_nodeMapping(fileUnit) - call mesh_build_elemMapping(fileUnit) + call mesh_build_nodeMapping(fileUnit) + call mesh_build_elemMapping(fileUnit) call mesh_get_nodeElemDimensions(fileUnit) call mesh_build_nodes(fileUnit) call mesh_build_elements(fileUnit) @@ -162,7 +162,7 @@ call IO_error(100) endif - END SUBROUTINE + END SUBROUTINE !*********************************************************** @@ -232,38 +232,38 @@ integer(pInt) mesh_faceMatch integer(pInt), dimension(FE_NfaceNodes(face,FE_mapElemtype(mesh_element(2,elem)))) :: nodeMap integer(pInt) minN,NsharedElems,lonelyNode,faceNode,i,n,t - + minN = mesh_maxNsharedElems+1 ! init to worst case - mesh_faceMatch = 0_pInt ! intialize to "no match found" - t = FE_mapElemtype(mesh_element(2,elem)) ! figure elemType - + mesh_faceMatch = 0_pInt ! intialize to "no match found" + t = FE_mapElemtype(mesh_element(2,elem)) ! figure elemType + do faceNode=1,FE_NfaceNodes(face,t) ! loop over nodes on face nodeMap(faceNode) = mesh_FEasCP('node',mesh_element(4+FE_nodeOnFace(faceNode,face,t),elem)) ! CP id of face node - NsharedElems = mesh_sharedElem(1,nodeMap(faceNode)) ! figure # shared elements for this node + NsharedElems = mesh_sharedElem(1,nodeMap(faceNode)) ! figure # shared elements for this node if (NsharedElems < minN) then minN = NsharedElems ! remember min # shared elems lonelyNode = faceNode ! remember most lonely node endif - end do + end do candidate: do i=1,minN ! iterate over lonelyNode's shared elements - mesh_faceMatch = mesh_sharedElem(1+i,nodeMap(lonelyNode)) ! present candidate elem + mesh_faceMatch = mesh_sharedElem(1+i,nodeMap(lonelyNode)) ! present candidate elem if (mesh_faceMatch == elem) then ! my own element ? mesh_faceMatch = 0_pInt ! disregard cycle candidate endif - do faceNode=1,FE_NfaceNodes(face,t) ! check remaining face nodes to match - if (faceNode == lonelyNode) cycle ! disregard lonely node (matches anyway) + do faceNode=1,FE_NfaceNodes(face,t) ! check remaining face nodes to match + if (faceNode == lonelyNode) cycle ! disregard lonely node (matches anyway) n = nodeMap(faceNode) if (all(mesh_sharedElem(2:1+mesh_sharedElem(1,n),n) /= mesh_faceMatch)) then ! no ref to candidate elem? - mesh_faceMatch = 0_pInt ! set to "no match" (so far) + mesh_faceMatch = 0_pInt ! set to "no match" (so far) cycle candidate ! next candidate elem endif - end do + end do exit ! surviving candidate end do candidate - + return - + END FUNCTION @@ -372,51 +372,51 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements return END SUBROUTINE - -!******************************************************************** -! Build node mapping from FEM to CP -! + +!******************************************************************** +! Build node mapping from FEM to CP +! ! allocate globals: ! _mapFEtoCPnode -!******************************************************************** - SUBROUTINE mesh_build_nodeMapping (unit) - - use prec, only: pInt - use math, only: qsort - use IO - implicit none - - integer(pInt), dimension (mesh_Nnodes) :: node_count - integer(pInt) unit,i - integer(pInt), dimension (133) :: pos - character*300 line - -610 FORMAT(A300) - - allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt - node_count(:) = 0_pInt - - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,1) - if( IO_lc(IO_stringValue(line,pos,1)) == 'coordinates' ) then - read (unit,610,END=620) line ! skip crap line - do i=1,mesh_Nnodes - read (unit,610,END=620) line - mesh_mapFEtoCPnode(1,i) = IO_fixedIntValue (line,(/0,10/),1) - mesh_mapFEtoCPnode(2,i) = i - end do - exit - end if - end do - -620 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2)) - - return - END SUBROUTINE - - +!******************************************************************** + SUBROUTINE mesh_build_nodeMapping (unit) + + use prec, only: pInt + use math, only: qsort + use IO + implicit none + + integer(pInt), dimension (mesh_Nnodes) :: node_count + integer(pInt) unit,i + integer(pInt), dimension (133) :: pos + character*300 line + +610 FORMAT(A300) + + allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt + node_count(:) = 0_pInt + + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,1) + if( IO_lc(IO_stringValue(line,pos,1)) == 'coordinates' ) then + read (unit,610,END=620) line ! skip crap line + do i=1,mesh_Nnodes + read (unit,610,END=620) line + mesh_mapFEtoCPnode(1,i) = IO_fixedIntValue (line,(/0,10/),1) + mesh_mapFEtoCPnode(2,i) = i + end do + exit + end if + end do + +620 call qsort(mesh_mapFEtoCPnode,1,size(mesh_mapFEtoCPnode,2)) + + return + END SUBROUTINE + + !******************************************************************** ! Build element mapping from FEM to CP ! @@ -465,127 +465,132 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements END SUBROUTINE -!******************************************************************** +!******************************************************************** ! store x,y,z coordinates of all nodes in mesh ! ! allocate globals: ! _node -!******************************************************************** - SUBROUTINE mesh_build_nodes (unit) - - use prec, only: pInt - use IO - implicit none - - integer unit,i,j,m - integer(pInt), dimension(3) :: pos - integer(pInt), dimension(5), parameter :: node_ends = (/0,10,30,50,70/) - character*300 line - - allocate ( mesh_node (3,mesh_Nnodes) ) - mesh_node(:,:) = 0_pInt - -610 FORMAT(A300) - - rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,1) - if( IO_lc(IO_stringValue(line,pos,1)) == 'coordinates' ) then - read (unit,610,END=620) line ! skip crap line - do i=1,mesh_Nnodes - read (unit,610,END=620) line - m = mesh_FEasCP('node',IO_fixedIntValue (line,node_ends,1)) - do j=1,3 - mesh_node(j,m) = IO_fixedNoEFloatValue (line,node_ends,j+1) - end do - end do - exit - end if - end do - -620 return - - END SUBROUTINE - - -!******************************************************************** -! store FEid, type, mat, tex, and node list per element -! -! allocate globals: -! _element -!******************************************************************** - SUBROUTINE mesh_build_elements (unit) - - use prec, only: pInt - use IO - implicit none - - integer unit,i,j,sv,val,CP_elem - integer(pInt), dimension(133) :: pos - integer(pInt), dimension(1+mesh_NcpElems) :: contInts - character*300 line - - allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt - -610 FORMAT(A300) - +!******************************************************************** + SUBROUTINE mesh_build_nodes (unit) + + use prec, only: pInt + use IO + implicit none + + integer unit,i,j,m + integer(pInt), dimension(3) :: pos + integer(pInt), dimension(5), parameter :: node_ends = (/0,10,30,50,70/) + character*300 line + + allocate ( mesh_node (3,mesh_Nnodes) ) + mesh_node(:,:) = 0_pInt + +610 FORMAT(A300) + rewind(unit) - do - read (unit,610,END=620) line - pos = IO_stringPos(line,2) - - if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then - read (unit,610,END=620) line ! Garbage line - do i=1,mesh_Nelems - read (unit,610,END=620) line - pos = IO_stringPos(line,66) ! limit to 64 nodes max (plus ID, type) - CP_elem = mesh_FEasCP('elem',IO_intValue(line,pos,1)) - if (CP_elem /= 0) then ! disregard non CP elems - mesh_element (1,CP_elem) = IO_IntValue (line,pos,1) ! FE id - mesh_element (2,CP_elem) = IO_IntValue (line,pos,2) ! elem type - do j=1,FE_Nnodes(FE_mapElemtype(mesh_element(2,CP_elem))) - mesh_element(j+4,CP_elem) = IO_IntValue (line,pos,j+2) ! copy FE ids of nodes - end do - end if - end do - exit - endif - enddo + do + read (unit,610,END=620) line + pos = IO_stringPos(line,1) + if( IO_lc(IO_stringValue(line,pos,1)) == 'coordinates' ) then + read (unit,610,END=620) line ! skip crap line + do i=1,mesh_Nnodes + read (unit,610,END=620) line + m = mesh_FEasCP('node',IO_fixedIntValue (line,node_ends,1)) + do j=1,3 + mesh_node(j,m) = IO_fixedNoEFloatValue (line,node_ends,j+1) + end do + end do + exit + end if + end do - do ! fast forward to first "initial state" section - read (unit,610,END=620) line - if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial').and. & - (IO_lc(IO_stringValue(line,pos,2)) == 'state') ) exit - enddo - - do ! parse initial state section(s) - if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial').and. & - (IO_lc(IO_stringValue(line,pos,2)) == 'state') ) then - read (unit,610,END=620) line - pos = IO_stringPos(line,1) - sv = IO_IntValue (line,pos,1) ! figure state variable index - if( (sv == 2).or.(sv == 3) ) then ! only state vars 2 and 3 of interest - read (unit,610,END=620) line - do while (scan(IO_stringValue(line,pos,1),'+-',back=.true.)>1) - val = NINT(IO_fixedNoEFloatValue (line,(/0,20/),1)) ! state var's value - contInts = IO_continousIntValues(unit,mesh_Nelems) ! get affected elements - do i = 1,contInts(1) - CP_elem = mesh_FEasCP('elem',contInts(1+i)) - mesh_element(1+sv,CP_elem) = val - enddo - read (unit,610,END=620) line ! ignore IP range - read (unit,610,END=620) line ! read ahead (check in do loop) - enddo - endif - endif - enddo - -620 return - - END SUBROUTINE +620 return - + END SUBROUTINE + + +!******************************************************************** +! store FEid, type, mat, tex, and node list per element +! +! allocate globals: +! _element +!******************************************************************** + SUBROUTINE mesh_build_elements (unit) + + use prec, only: pInt + use IO + implicit none + + integer unit,i,j,sv,val,CP_elem + integer(pInt), dimension(133) :: pos + integer(pInt), dimension(1+mesh_NcpElems) :: contInts + character*300 line + + allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt + +610 FORMAT(A300) + + + rewind(unit) + do + read (unit,610,END=620) line + pos = IO_stringPos(line,2) + if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then + read (unit,610,END=620) line ! Garbage line + do i=1,mesh_Nelems + read (unit,610,END=620) line + pos = IO_stringPos(line,66) ! limit to 64 nodes max (plus ID, type) + CP_elem = mesh_FEasCP('elem',IO_intValue(line,pos,1)) + if (CP_elem /= 0) then ! disregard non CP elems + mesh_element (1,CP_elem) = IO_IntValue (line,pos,1) ! FE id + mesh_element (2,CP_elem) = IO_IntValue (line,pos,2) ! elem type + do j=1,FE_Nnodes(FE_mapElemtype(mesh_element(2,CP_elem))) + mesh_element(j+4,CP_elem) = IO_IntValue (line,pos,j+2) ! copy FE ids of nodes + end do + end if + end do + exit + endif + enddo + + rewind(unit) ! just in case "initial state" apears before "connectivity" + do ! fast forward to first "initial state" section + read (unit,610,END=620) line + pos = IO_stringPos(line,2) + if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial').and. & + (IO_lc(IO_stringValue(line,pos,2)) == 'state') ) exit + enddo + + do ! parse initial state section(s) + if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial').and. & + (IO_lc(IO_stringValue(line,pos,2)) == 'state') ) then + read (unit,610,END=620) line + pos = IO_stringPos(line,1) + sv = IO_IntValue (line,pos,1) ! figure state variable index + if( (sv == 2).or.(sv == 3) ) then ! only state vars 2 and 3 of interest + read (unit,610,END=620) line + pos = IO_stringPos(line,1) + do while (scan(IO_stringValue(line,pos,1),'+-',back=.true.)>1) + val = NINT(IO_fixedNoEFloatValue (line,(/0,20/),1)) ! state var's value + contInts = IO_continousIntValues(unit,mesh_Nelems) ! get affected elements + do i = 1,contInts(1) + CP_elem = mesh_FEasCP('elem',contInts(1+i)) + mesh_element(1+sv,CP_elem) = val + enddo + read (unit,610,END=620) line ! ignore IP range + pos = IO_stringPos(line,1) + enddo + endif + endif + read (unit,610,END=620) line ! read ahead (check in do loop) + pos = IO_stringPos(line,2) + enddo + +620 return + + END SUBROUTINE + + !******************************************************************** ! build list of elements shared by each node in mesh !