added function calls in init()

This commit is contained in:
Philip Eisenlohr 2007-04-10 11:22:53 +00:00
parent 35266bd513
commit df3f327ece
1 changed files with 341 additions and 312 deletions

View File

@ -9,7 +9,6 @@
! --------------------------- ! ---------------------------
! _Nelems : total number of elements in mesh ! _Nelems : total number of elements in mesh
! _NcpElems : total number of CP 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 ! _Nnodes : total number of nodes in mesh
! _maxNnodes : max number of nodes in any CP element ! _maxNnodes : max number of nodes in any CP element
! _maxNips : max number of IPs 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] ! _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 ! 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) mesh_Nnodes,mesh_maxNnodes,mesh_maxNips,mesh_maxNipNeighbors,mesh_maxNsharedElems
integer(pInt), dimension(:,:), allocatable, target :: mesh_mapFEtoCPelem,mesh_mapFEtoCPnode integer(pInt), dimension(:,:), allocatable, target :: mesh_mapFEtoCPelem,mesh_mapFEtoCPnode
integer(pInt), dimension(:,:), allocatable :: mesh_element, mesh_sharedElem integer(pInt), dimension(:,:), allocatable :: mesh_element, mesh_sharedElem
@ -129,7 +128,13 @@
!*********************************************************** !***********************************************************
! initialization ! 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_Nelems = 0_pInt
mesh_NcpElems = 0_pInt mesh_NcpElems = 0_pInt
@ -142,14 +147,82 @@
FE_mapElemtype( 7) = 1 FE_mapElemtype( 7) = 1
FE_mapElemtype(134) = 2 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 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 ! find face-matching element of same type
! !!***********************************************************
!
!***********************************************************
FUNCTION mesh_faceMatch(face,elem) FUNCTION mesh_faceMatch(face,elem)
use prec, only: pInt use prec, only: pInt
@ -194,120 +267,117 @@ candidate: do i=1,minN ! iterate over lonelyNode's shared elements
END FUNCTION END FUNCTION
!*********************************************************** !********************************************************************
! build up of IP neighborhood ! get count of elements, nodes, and cp elements in mesh
!*********************************************************** ! for subsequent array allocations
SUBROUTINE mesh_build_ipNeighborhood() !
! assign globals:
! _Nelems, _Nnodes, _NcpElems
!********************************************************************
SUBROUTINE mesh_get_meshDimensions (unit)
use prec, only: pInt use prec, only: pInt
use IO
implicit none implicit none
integer(pInt) e,t,i,j,k,n integer(pInt) unit,i,pos(41)
integer(pInt) neighbor,neighboringElem,neighboringIP,matchingElem,faceNode,linkingNode character*300 line
allocate(mesh_ipNeighborhood(2,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) 610 FORMAT(A300)
do e = 1,mesh_NcpElems ! loop over cpElems rewind(unit)
t = FE_mapElemtype(mesh_element(2,e)) ! get elemType do
do i = 1,FE_Nips(t) ! loop over IPs of elem read (unit,610,END=620) line
do n = 1,FE_NipNeighbors(t) ! loop over neighbors of IP pos = IO_stringPos(line,20)
neighbor = FE_ipNeighbor(n,i,t)
if (neighbor > 0) then ! intra-element IP select case ( IO_lc(IO_Stringvalue(line,pos,1)))
neighboringElem = e case('sizing')
neighboringIP = neighbor mesh_Nelems = IO_IntValue (line,pos,3)
else ! neighboring element's IP mesh_Nnodes = IO_IntValue (line,pos,4)
neighboringElem = 0_pInt case('hypoelastic')
neighboringIP = 0_pInt do i=1,4
matchingElem = mesh_faceMatch(-neighbor,e) ! get CP elem id of face match read (unit,610,END=620) line
if (matchingElem > 0 .and. & end do
FE_mapElemtype(mesh_element(2,matchingElem)) == t) then ! found match of same type? pos = IO_stringPos(line,20)
matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matching face if( IO_lc(IO_Stringvalue(line,pos,2)).eq.'to' )then
faceNode = FE_nodeOnFace(j,-neighbor,t) ! get face node id mesh_NcpElems = IO_IntValue(line,pos,3)-IO_IntValue(line,pos,1)+1
if (i == FE_ipAtNode(faceNode,t)) then ! ip linked to face node is me? else
linkingNode = mesh_element(4+faceNode,e) ! FE id of this facial node mesh_NcpElems = mesh_NcpElems + pos(1)
do k = 1,FE_Nnodes(t) ! loop over nodes in matching element do while( IO_lc(IO_Stringvalue(line,pos,pos(1))).eq.'c' )
if (linkingNode == mesh_element(4+k,matchingElem)) then mesh_NcpElems = mesh_NcpElems - 1 ! Counted the c character from the line
neighboringElem = matchingElem read (unit,610,END=620) line
neighboringIP = FE_ipAtNode(k,t) pos = IO_stringPos(line,20)
exit matchFace mesh_NcpElems = mesh_NcpElems + pos(1)
endif end do
end do end if
endif
end do matchFace end select
endif
endif
mesh_ipNeighborhood(1,n,i,e) = neighboringElem
mesh_ipNeighborhood(2,n,i,e) = neighboringIP
end do
end do
end do end do
return 620 return
END SUBROUTINE END SUBROUTINE
!********************************************************************
!*********************************************************** ! get maximum count of nodes, IPs, IP neighbors, and shared elements
! FE to CP id mapping by binary search thru lookup array ! for subsequent array allocations
! !
! valid questions are 'elem', 'node' ! assign globals:
!*********************************************************** ! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsharedElems
FUNCTION mesh_FEasCP(what,id) !********************************************************************
SUBROUTINE mesh_get_nodeElemDimensions (unit)
use prec, only: pInt use prec, only: pInt
use IO, only: IO_lc use IO
implicit none implicit none
character(len=*), intent(in) :: what integer(pInt), dimension (mesh_Nnodes) :: node_count
integer(pInt), intent(in) :: id integer(pInt) unit,i,j,n,t,e
integer(pInt), dimension(:,:), pointer :: lookupMap integer(pInt), dimension (133) :: pos
integer(pInt) mesh_FEasCP, lower,upper,center character*300 line
mesh_FEasCP = 0_pInt 610 FORMAT(A300)
select case(IO_lc(what(1:4)))
case('elem') node_count = 0_pInt
lookupMap => mesh_mapFEtoCPelem
case('node') rewind(unit)
lookupMap => mesh_mapFEtoCPnode do
case default read (unit,610,END=630) line
return pos = IO_stringPos(line,1)
end select if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then
read (unit,610,END=630) line ! Garbage line
lower = 1_pInt do i=1,mesh_Nelems ! read all elements
upper = size(lookupMap,2) read (unit,610,END=630) line
pos = IO_stringPos(line,66) ! limit to 64 nodes max (plus ID, type)
! check at bounds e = mesh_FEasCP('elem',IO_intValue(line,pos,1))
if (lookupMap(1,lower) == id) then if (e /= 0) then
mesh_FEasCP = lookupMap(2,lower) t = FE_mapElemtype(IO_intValue(line,pos,2))
return mesh_maxNnodes = max(mesh_maxNnodes,FE_Nnodes(t))
elseif (lookupMap(1,upper) == id) then mesh_maxNips = max(mesh_maxNips,FE_Nips(t))
mesh_FEasCP = lookupMap(2,upper) mesh_maxNipNeighbors = max(mesh_maxNipNeighbors,FE_NipNeighbors(t))
return do j=1,FE_Nnodes(t)
endif n = mesh_FEasCP('node',IO_IntValue (line,pos,j+2))
node_count(n) = node_count(n)+1
! binary search in between bounds end do
do while (upper-lower > 1) end if
center = (lower+upper)/2 end do
if (lookupMap(1,center) < id) then exit
lower = center
elseif (lookupMap(1,center) > id) then
upper = center
else
mesh_FEasCP = lookupMap(2,center)
exit
end if end if
end do end do
return
END FUNCTION
630 mesh_maxNsharedElems = maxval(node_count)
return
END SUBROUTINE
!******************************************************************** !********************************************************************
! Build node mapping from FEM to CP ! Build node mapping from FEM to CP
! !
! assign globals: ! allocate globals:
! _maxNsharedElems ! _mapFEtoCPnode
!******************************************************************** !********************************************************************
SUBROUTINE mesh_build_nodeMapping (unit) 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) 610 FORMAT(A300)
allocate ( mesh_mapFEtoCPnode(2,mesh_Nnodes) ) allocate (mesh_mapFEtoCPnode(2,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt
mesh_mapFEtoCPnode(:,:) = 0_pInt
node_count(:) = 0_pInt node_count(:) = 0_pInt
rewind(unit) rewind(unit)
@ -348,150 +417,61 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
END SUBROUTINE 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 ! store x,y,z coordinates of all nodes in mesh
! !
! assign globals: ! allocate globals:
! _maxNnodes, _maxNips, _maxNipNeighbors, _maxNsharedElems ! _node
!******************************************************************** !********************************************************************
SUBROUTINE mesh_build_maxNofCPelems (unit) SUBROUTINE mesh_build_nodes (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 prec, only: pInt
use IO use IO
@ -530,9 +510,10 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
!******************************************************************** !********************************************************************
! ! store FEid, type, mat, tex, and node list per element
! assign globals: !
! _maxNnodes, _maxNips, _maxNipNeighbors ! allocate globals:
! _element
!******************************************************************** !********************************************************************
SUBROUTINE mesh_build_elements (unit) 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 integer(pInt), dimension(1+mesh_NcpElems) :: contInts
character*300 line character*300 line
rewind(unit) allocate (mesh_element (4+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt
allocate ( mesh_element (4+mesh_maxNnodes,mesh_NcpElems) )
write(*,*) 'allocated',4+mesh_maxNnodes,mesh_NcpElems
mesh_element(:,:) = 0_pInt
610 FORMAT(A300) 610 FORMAT(A300)
rewind(unit)
do do
read (unit,610,END=620) line read (unit,610,END=620) line
pos = IO_stringPos(line,2) pos = IO_stringPos(line,2)
@ -573,8 +552,8 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
exit exit
endif endif
enddo 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 read (unit,610,END=620) line
if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial').and. & if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial').and. &
(IO_lc(IO_stringValue(line,pos,2)) == 'state') ) exit (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 read (unit,610,END=620) line
do while (scan(IO_stringValue(line,pos,1),'+-',back=.true.)>1) do while (scan(IO_stringValue(line,pos,1),'+-',back=.true.)>1)
val = NINT(IO_fixedNoEFloatValue (line,(/0,20/),1)) ! state var's value 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) do i = 1,contInts(1)
CP_elem = mesh_FEasCP('elem',contInts(1+i)) CP_elem = mesh_FEasCP('elem',contInts(1+i))
mesh_element(1+sv,CP_elem) = val mesh_element(1+sv,CP_elem) = val
enddo enddo
read (unit,610,END=620) line ! ignore IP range 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 enddo
endif endif
endif endif
@ -607,59 +586,109 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
END SUBROUTINE END SUBROUTINE
!******************************************************************** !********************************************************************
! Get global variables ! build list of elements shared by each node in mesh
! !
! assign globals: ! allocate globals:
! _Nelems, _Nnodes, NcpElem ! _sharedElem
!******************************************************************** !********************************************************************
SUBROUTINE mesh_get_globals (unit) SUBROUTINE mesh_build_sharedElems (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
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 END MODULE mesh