diff --git a/trunk/mesh.f90 b/trunk/mesh.f90 index 0133338c0..5f13294de 100644 --- a/trunk/mesh.f90 +++ b/trunk/mesh.f90 @@ -9,7 +9,6 @@ ! --------------------------- ! _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 @@ -41,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,mesh_NelemTypes + 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 @@ -129,7 +128,13 @@ !*********************************************************** ! initialization !*********************************************************** - SUBROUTINE mesh_init () + SUBROUTINE mesh_init () + + use prec, only: pInt + use IO, only: IO_open_InputFile + implicit none + + integer(pInt), parameter :: fileUnit = 222 mesh_Nelems = 0_pInt mesh_NcpElems = 0_pInt @@ -142,14 +147,82 @@ FE_mapElemtype( 7) = 1 FE_mapElemtype(134) = 2 -! call to various subrountes to parse the stuff from the input file... +! 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_get_nodeElemDimensions(fileUnit) + call mesh_build_nodeMapping(fileUnit) + call mesh_build_elemMapping(fileUnit) + call mesh_build_nodes(fileUnit) + call mesh_build_elements(fileUnit) + call mesh_build_sharedElems(fileUnit) + call mesh_build_ipNeighborhood() + close (fileUnit) + else + call IO_error(100) + endif + 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 QUESTION is it valid to extend bounds by 1 and just do binary search w/o init 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 + + !*********************************************************** ! find face-matching element of same type -! -! -!*********************************************************** +!!*********************************************************** FUNCTION mesh_faceMatch(face,elem) use prec, only: pInt @@ -194,120 +267,117 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements END FUNCTION -!*********************************************************** -! build up of IP neighborhood -!*********************************************************** - SUBROUTINE mesh_build_ipNeighborhood() - +!******************************************************************** +! get count of elements, nodes, and cp elements in mesh +! for subsequent array allocations +! +! assign globals: +! _Nelems, _Nnodes, _NcpElems +!******************************************************************** + SUBROUTINE mesh_get_meshDimensions (unit) + use prec, only: pInt + use IO 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 + + integer(pInt) unit,i,pos(41) + character*300 line + +610 FORMAT(A300) + + 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 - - return + +620 return END SUBROUTINE + - -!*********************************************************** -! FE to CP id mapping by binary search thru lookup array +!******************************************************************** +! get maximum count of nodes, IPs, IP neighbors, and shared elements +! for subsequent array allocations ! -! valid questions are 'elem', 'node' -!*********************************************************** - FUNCTION mesh_FEasCP(what,id) - +! assign globals: +! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsharedElems +!******************************************************************** + SUBROUTINE mesh_get_nodeElemDimensions (unit) + use prec, only: pInt - use IO, only: IO_lc + use IO 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 + + 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 - return - - END FUNCTION - +630 mesh_maxNsharedElems = maxval(node_count) + + return + END SUBROUTINE + !******************************************************************** ! Build node mapping from FEM to CP ! -! assign globals: -! _maxNsharedElems +! allocate globals: +! _mapFEtoCPnode !******************************************************************** SUBROUTINE mesh_build_nodeMapping (unit) @@ -323,8 +393,7 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc 610 FORMAT(A300) - allocate ( mesh_mapFEtoCPnode(2,mesh_Nnodes) ) - mesh_mapFEtoCPnode(:,:) = 0_pInt + allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt node_count(:) = 0_pInt rewind(unit) @@ -348,150 +417,61 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc END SUBROUTINE +!******************************************************************** +! Build element mapping from FEM to CP +! +! allocate globals: +! _mapFEtoCPelem +!******************************************************************** + 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 + + !******************************************************************** -! Build node mapping from FEM to CP -! -! assign globals: -! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsharedElems +! store x,y,z coordinates of all nodes in mesh +! +! allocate globals: +! _node !******************************************************************** - 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) + SUBROUTINE mesh_build_nodes (unit) use prec, only: pInt use IO @@ -530,9 +510,10 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc !******************************************************************** -! -! assign globals: -! _maxNnodes, _maxNips, _maxNipNeighbors +! store FEid, type, mat, tex, and node list per element +! +! allocate globals: +! _element !******************************************************************** SUBROUTINE mesh_build_elements (unit) @@ -545,13 +526,11 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc 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 + 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) @@ -573,8 +552,8 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc exit endif enddo - write(*,*) 'done with connectivity.' - do ! fast forward to "initial state" sections + + 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 @@ -590,13 +569,13 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc 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 + 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 to check in do loop + read (unit,610,END=620) line ! read ahead (check in do loop) enddo endif endif @@ -607,59 +586,109 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc 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 +!******************************************************************** +! build list of elements shared by each node in mesh +! +! allocate globals: +! _sharedElem +!******************************************************************** + 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 + + +!*********************************************************** +! build up of IP neighborhood +! +! allocate globals +! _ipNeighborhood +!*********************************************************** + 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)) ; mesh_ipNeighborhood = 0_pInt + + 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 + END MODULE mesh \ No newline at end of file