!############################################################## MODULE mesh !############################################################## use prec, only: pReal,pInt implicit none ! --------------------------- ! _Nelems : total number of elements in mesh ! _NcpElems : total number of CP elements in mesh ! _NelemTypes: total number of element types 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 ! _maxNsharedElems : max number of CP elements sharing a node ! ! _element : FEid, type, material, texture, node indices ! _node : x,y,z coordinates (initially!) ! _sharedElem : entryCount and list of elements containing node ! ! _mapFEtoCPelem : [sorted FEid, corresponding CPid] ! _mapFEtoCPnode : [sorted FEid, corresponding CPid] ! ! MISSING: these definitions should actually reside in the ! FE-solver specific part (different for MARC/ABAQUS)..! ! Hence, I suggest to prefix with "FE_" ! ! _mapElementtype : map MARC/ABAQUS elemtype to 1-maxN ! ! _Nnodes : # nodes in a specific type of element ! _Nips : # IPs in a specific type of element ! _NipNeighbors : # IP neighbors in a specific type of element ! _ipNeighbor : +x,-x,+y,-y,+z,-z list of intra-element IPs and ! (negative) neighbor faces per own IP in a specific type of element ! _NfaceNodes : # nodes per face in a specific type of element ! _nodeOnFace : list of node indices on each face of a specific type of element ! _ipAtNode : map node index to IP index in a specific type of element ! _nodeAtIP : map IP index to node index in a specific type of element ! _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,mesh_NelemTypes 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 integer(pInt), dimension(:,:,:,:), allocatable :: mesh_ipNeighborhood real(pReal), allocatable :: mesh_node (:,:) integer(pInt), parameter :: FE_Nelemtypes = 2 integer(pInt), parameter :: FE_maxNnodes = 8 integer(pInt), parameter :: FE_maxNips = 8 integer(pInt), parameter :: FE_maxNneighbors = 6 integer(pInt), parameter :: FE_maxNfaceNodes = 4 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 /) integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_Nips = & (/8, & ! element 7 1 & ! element 134 /) integer(pInt), dimension(FE_Nelemtypes), parameter :: FE_NipNeighbors = & (/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 /),(/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 /),(/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 /),(/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 & /),(/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 & /),(/FE_maxNneighbors,FE_maxNips,FE_Nelemtypes/)) CONTAINS ! --------------------------- ! subroutine mesh_init() ! function mesh_FEtoCPelement(FEid) ! function mesh_build_ipNeighorhood() ! subroutine mesh_parse_inputFile() ! --------------------------- !*********************************************************** ! initialization !*********************************************************** SUBROUTINE mesh_init () 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... END SUBROUTINE !*********************************************************** ! find face-matching element of same type ! ! !*********************************************************** FUNCTION mesh_faceMatch(face,elem) use prec, only: pInt implicit none integer(pInt) face,elem integer(pInt) mesh_faceMatch integer(pInt), dimension(FE_NfaceNodes(face,FE_mapElemtype(mesh_element(2,elem)))) :: nodeMap integer(pInt) minN,NsharedElems,lonelyNode,faceNode,i,j,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 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 if (NsharedElems < minN) then minN = NsharedElems ! remember min # shared elems lonelyNode = faceNode ! remember most lonely node endif end do candidate: do i=1,minN ! iterate over lonelyNode's shared elements 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) 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) cycle candidate ! next candidate elem endif end do exit ! surviving candidate end do candidate return END FUNCTION !*********************************************************** ! build up of IP neighborhood !*********************************************************** SUBROUTINE mesh_build_ipNeighborhood() use prec, only: pInt implicit none integer(pInt) e,t,i,j,k,n integer(pInt) neighbor,neighboringElem,neighboringIP,matchingElem,faceNode,linkingNode allocate(mesh_ipNeighborhood(2,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) do e = 1,mesh_NcpElems ! loop over cpElems t = FE_mapElemtype(mesh_element(2,e)) ! get elemType do i = 1,FE_Nips(t) ! loop over IPs of elem do n = 1,FE_NipNeighbors(t) ! loop over neighbors of IP neighbor = FE_ipNeighbor(n,i,t) if (neighbor > 0) then ! intra-element IP neighboringElem = e neighboringIP = neighbor else ! neighboring element's IP neighboringElem = 0_pInt neighboringIP = 0_pInt matchingElem = mesh_faceMatch(-neighbor,e) ! get CP elem id of face match if (matchingElem > 0 .and. & FE_mapElemtype(mesh_element(2,matchingElem)) == t) then ! found match of same type? matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matching face faceNode = FE_nodeOnFace(j,-neighbor,t) ! get face node id if (i == FE_ipAtNode(faceNode,t)) then ! ip linked to face node is me? linkingNode = mesh_element(4+faceNode,e) ! FE id of this facial node do k = 1,FE_Nnodes(t) ! loop over nodes in matching element if (linkingNode == mesh_element(4+k,matchingElem)) then neighboringElem = matchingElem neighboringIP = FE_ipAtNode(k,t) exit matchFace endif end do endif end do matchFace endif endif mesh_ipNeighborhood(1,n,i,e) = neighboringElem mesh_ipNeighborhood(2,n,i,e) = neighboringIP end do end do end do return END SUBROUTINE !*********************************************************** ! FE to CP id mapping by binary search thru lookup array ! ! valid questions are 'elem', 'node' !*********************************************************** FUNCTION mesh_FEasCP(what,id) use prec, only: pInt use IO, only: IO_lc implicit none character(len=*), intent(in) :: what integer(pInt), intent(in) :: id integer(pInt), dimension(:,:), pointer :: lookupMap integer(pInt) mesh_FEasCP, lower,upper,center mesh_FEasCP = 0_pInt select case(IO_lc(what(1:4))) case('elem') lookupMap => mesh_mapFEtoCPelem case('node') lookupMap => mesh_mapFEtoCPnode case default return end select lower = 1_pInt upper = size(lookupMap,2) ! check at bounds if (lookupMap(1,lower) == id) then mesh_FEasCP = lookupMap(2,lower) return elseif (lookupMap(1,upper) == id) then mesh_FEasCP = lookupMap(2,upper) return endif ! binary search in between bounds do while (upper-lower > 1) center = (lower+upper)/2 if (lookupMap(1,center) < id) then lower = center elseif (lookupMap(1,center) > id) then upper = center else mesh_FEasCP = lookupMap(2,center) exit end if end do return END FUNCTION !******************************************************************** ! Build node mapping from FEM to CP ! ! assign globals: ! _maxNsharedElems !******************************************************************** 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,j,Nnodes,cur_node 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 node mapping from FEM to CP ! ! assign globals: ! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsharedElems !******************************************************************** SUBROUTINE mesh_build_maxNofCPelems (unit) use prec, only: pInt use IO implicit none integer(pInt), dimension (mesh_Nnodes) :: node_count integer(pInt) unit,i,j,n,t,e integer(pInt), dimension (133) :: pos character*300 line 610 FORMAT(A300) node_count = 0_pInt rewind(unit) do read (unit,610,END=630) line pos = IO_stringPos(line,1) if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then read (unit,610,END=630) line ! Garbage line do i=1,mesh_Nelems ! read all elements read (unit,610,END=630) line pos = IO_stringPos(line,66) ! limit to 64 nodes max (plus ID, type) e = mesh_FEasCP('elem',IO_intValue(line,pos,1)) if (e /= 0) then t = FE_mapElemtype(IO_intValue(line,pos,2)) mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t)) mesh_maxNips = max(mesh_maxNips,FE_Nips(t)) mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t)) do j=1,FE_Nnodes(t) n = mesh_FEasCP('node',IO_IntValue (line,pos,j+2)) node_count(n) = node_count(n)+1 end do end if end do exit end if end do 630 mesh_maxNsharedElems = maxval(node_count) return END SUBROUTINE !******************************************************************** ! Build element mapping from FEM to CP !******************************************************************** SUBROUTINE mesh_build_elemMapping (unit) use prec, only: pInt use math, only: qsort use IO implicit none integer unit, i,CP_elem character*300 line integer(pInt), dimension (3) :: pos integer(pInt), dimension (1+mesh_NcpElems) :: contInts 610 FORMAT(A300) allocate ( mesh_mapFEtoCPelem(2,mesh_NcpElems) ) mesh_mapFEtoCPelem(:,:) = 0_pInt CP_elem = 0_pInt rewind(unit) do read (unit,610,END=620) line pos = IO_stringPos(line,1) if( IO_lc(IO_stringValue(line,pos,1)) == 'hypoelastic' ) then do i=1,3 ! skip three lines read (unit,610,END=620) line end do contInts = IO_continousIntValues(unit,mesh_NcpElems) do i = 1,contInts(1) CP_elem = CP_elem+1 mesh_mapFEtoCPelem(1,CP_elem) = contInts(1+i) mesh_mapFEtoCPelem(2,CP_elem) = CP_elem enddo end if end do 620 call qsort(mesh_mapFEtoCPelem,1,size(mesh_mapFEtoCPelem,2)) ! should be mesh_NcpElems return END SUBROUTINE !******************************************************************** !******************************************************************** SUBROUTINE mesh_build_sharedElems (unit) use prec, only: pInt use IO implicit none integer unit,i,j,CP_node,CP_elem integer(pInt), dimension (133) :: pos character*300 line 610 FORMAT(A300) allocate ( mesh_sharedElem( 1+mesh_maxNsharedElems,mesh_Nnodes) ) mesh_sharedElem(:,:) = 0_pInt rewind(unit) do read (unit,610,END=620) line pos = IO_stringPos(line,1) 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 do j = 1,FE_Nnodes(FE_mapElemtype(IO_intValue(line,pos,2))) CP_node = mesh_FEasCP('node',IO_IntValue (line,pos,j+2)) mesh_sharedElem(1,CP_node) = mesh_sharedElem(1,CP_node) + 1 mesh_sharedElem(1+mesh_sharedElem(1,CP_node),CP_node) = CP_elem enddo end if end do exit end if end do 620 return END SUBROUTINE !******************************************************************** !******************************************************************** SUBROUTINE mesh_build_nodeCoords (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 !******************************************************************** ! ! assign globals: ! _maxNnodes, _maxNips, _maxNipNeighbors !******************************************************************** SUBROUTINE mesh_build_elements (unit) use prec, only: pInt use IO implicit none integer unit,i,j,t,sv,val,CP_elem integer(pInt), dimension(133) :: pos integer(pInt), dimension(1+mesh_NcpElems) :: contInts character*300 line rewind(unit) allocate ( mesh_element (4+mesh_maxNnodes,mesh_NcpElems) ) write(*,*) 'allocated',4+mesh_maxNnodes,mesh_NcpElems mesh_element(:,:) = 0_pInt 610 FORMAT(A300) 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 write(*,*) 'done with connectivity.' do ! fast forward to "initial state" sections 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) ! read 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 to check in do loop enddo endif endif enddo 620 return END SUBROUTINE !******************************************************************** ! Get global variables ! ! assign globals: ! _Nelems, _Nnodes, NcpElem !******************************************************************** SUBROUTINE mesh_get_globals (unit) use prec, only: pInt use IO implicit none integer(pInt) unit,i,pos(41) character*300 line 610 FORMAT(A300) mesh_NelemTypes = 0_pInt rewind(unit) do read (unit,610,END=620) line pos = IO_stringPos(line,20) select case ( IO_lc(IO_Stringvalue(line,pos,1))) case('sizing') mesh_Nelems = IO_IntValue (line,pos,3) mesh_Nnodes = IO_IntValue (line,pos,4) case('hypoelastic') do i=1,4 read (unit,610,END=620) line end do pos = IO_stringPos(line,20) if( IO_lc(IO_Stringvalue(line,pos,2)).eq.'to' )then mesh_NcpElems = IO_IntValue(line,pos,3)-IO_IntValue(line,pos,1)+1 else mesh_NcpElems = mesh_NcpElems + pos(1) do while( IO_lc(IO_Stringvalue(line,pos,pos(1))).eq.'c' ) mesh_NcpElems = mesh_NcpElems - 1 ! Counted the c character from the line read (unit,610,END=620) line pos = IO_stringPos(line,20) mesh_NcpElems = mesh_NcpElems + pos(1) end do end if end select end do 620 return END SUBROUTINE END MODULE mesh