223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

cleaning

This commit is contained in:
Martin Diehl 2019-06-13 00:39:35 +02:00
parent 2d52560b6d
commit e117ffbc0c
1 changed files with 256 additions and 756 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

530
src/mesh_marc.f90
View File

@ -52,9 +52,9 @@ module mesh
real(pReal), dimension(:,:,:), allocatable:: & real(pReal), dimension(:,:,:), allocatable:: &
mesh_ipArea !< area of interface to neighboring IP (initially!) mesh_ipArea !< area of interface to neighboring IP (initially!)
real(pReal),dimension(:,:,:,:), allocatable :: & real(pReal),dimension(:,:,:,:), allocatable :: &
mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!)
! -------------------------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------------------------
type(tMesh) :: theMesh type(tMesh) :: theMesh
@ -67,12 +67,6 @@ module mesh
mesh_Ncells, & !< total number of cells in mesh mesh_Ncells, & !< total number of cells in mesh
mesh_maxNsharedElems !< max number of CP elements sharing a node mesh_maxNsharedElems !< max number of CP elements sharing a node
integer, dimension(:,:), allocatable :: &
mesh_sharedElem, & !< entryCount and list of elements containing node
mesh_nodeTwins !< node twins are surface nodes that lie exactly on opposite sides of the mesh (surfaces nodes with equal coordinate values in two dimensions)
logical, dimension(3) :: mesh_periodicSurface !< flag indicating periodic outer surfaces (used for fluxes)
integer, dimension(:,:), allocatable :: & integer, dimension(:,:), allocatable :: &
mesh_cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID mesh_cellnodeParent !< cellnode's parent element ID, cellnode's intra-element ID
@ -89,14 +83,8 @@ integer, dimension(:,:), allocatable :: &
! Hence, I suggest to prefix with "FE_" ! Hence, I suggest to prefix with "FE_"
integer, parameter :: & integer, parameter :: &
FE_Nelemtypes = 13, &
FE_Ngeomtypes = 10, & FE_Ngeomtypes = 10, &
FE_Ncelltypes = 4, & FE_Ncelltypes = 4, &
FE_maxNipNeighbors = 6, &
FE_maxmaxNnodesAtIP = 8, & !< max number of (equivalent) nodes attached to an IP
FE_maxNmatchingNodesPerFace = 4, &
FE_maxNfaces = 6, &
FE_maxNcellnodes = 64, &
FE_maxNcellnodesPerCell = 8, & FE_maxNcellnodesPerCell = 8, &
FE_maxNcellfaces = 6, & FE_maxNcellfaces = 6, &
FE_maxNcellnodesPerCellface = 4 FE_maxNcellnodesPerCellface = 4
@ -115,84 +103,6 @@ integer, dimension(:,:), allocatable :: &
8 & ! element 21 (3D 20node 27ip) 8 & ! element 21 (3D 20node 27ip)
],pInt) ],pInt)
integer, dimension(FE_maxNfaces,FE_Ngeomtypes), parameter :: FE_NmatchingNodesPerFace = & !< number of matching nodes per face in a specific type of element geometry
reshape(int([ &
2,2,2,0,0,0, & ! element 6 (2D 3node 1ip)
2,2,2,0,0,0, & ! element 125 (2D 6node 3ip)
2,2,2,2,0,0, & ! element 11 (2D 4node 4ip)
2,2,2,2,0,0, & ! element 27 (2D 8node 9ip)
3,3,3,3,0,0, & ! element 134 (3D 4node 1ip)
3,3,3,3,0,0, & ! element 127 (3D 10node 4ip)
3,4,4,4,3,0, & ! element 136 (3D 6node 6ip)
4,4,4,4,4,4, & ! element 117 (3D 8node 1ip)
4,4,4,4,4,4, & ! element 7 (3D 8node 8ip)
4,4,4,4,4,4 & ! element 21 (3D 20node 27ip)
],pInt),[FE_maxNipNeighbors,FE_Ngeomtypes])
integer, dimension(FE_maxNmatchingNodesPerFace,FE_maxNfaces,FE_Ngeomtypes), parameter :: FE_face = & !< List of node indices on each face of a specific type of element geometry
reshape(int([&
1,2,0,0 , & ! element 6 (2D 3node 1ip)
2,3,0,0 , &
3,1,0,0 , &
0,0,0,0 , &
0,0,0,0 , &
0,0,0,0 , &
1,2,0,0 , & ! element 125 (2D 6node 3ip)
2,3,0,0 , &
3,1,0,0 , &
0,0,0,0 , &
0,0,0,0 , &
0,0,0,0 , &
1,2,0,0 , & ! element 11 (2D 4node 4ip)
2,3,0,0 , &
3,4,0,0 , &
4,1,0,0 , &
0,0,0,0 , &
0,0,0,0 , &
1,2,0,0 , & ! element 27 (2D 8node 9ip)
2,3,0,0 , &
3,4,0,0 , &
4,1,0,0 , &
0,0,0,0 , &
0,0,0,0 , &
1,2,3,0 , & ! element 134 (3D 4node 1ip)
1,4,2,0 , &
2,3,4,0 , &
1,3,4,0 , &
0,0,0,0 , &
0,0,0,0 , &
1,2,3,0 , & ! element 127 (3D 10node 4ip)
1,4,2,0 , &
2,4,3,0 , &
1,3,4,0 , &
0,0,0,0 , &
0,0,0,0 , &
1,2,3,0 , & ! element 136 (3D 6node 6ip)
1,4,5,2 , &
2,5,6,3 , &
1,3,6,4 , &
4,6,5,0 , &
0,0,0,0 , &
1,2,3,4 , & ! element 117 (3D 8node 1ip)
2,1,5,6 , &
3,2,6,7 , &
4,3,7,8 , &
4,1,5,8 , &
8,7,6,5 , &
1,2,3,4 , & ! element 7 (3D 8node 8ip)
2,1,5,6 , &
3,2,6,7 , &
4,3,7,8 , &
4,1,5,8 , &
8,7,6,5 , &
1,2,3,4 , & ! element 21 (3D 20node 27ip)
2,1,5,6 , &
3,2,6,7 , &
4,3,7,8 , &
4,1,5,8 , &
8,7,6,5 &
],pInt),[FE_maxNmatchingNodesPerFace,FE_maxNfaces,FE_Ngeomtypes])
integer, dimension(FE_Ncelltypes), parameter :: FE_NcellnodesPerCellface = & !< number of cell nodes per cell face in a specific cell type integer, dimension(FE_Ncelltypes), parameter :: FE_NcellnodesPerCellface = & !< number of cell nodes per cell face in a specific cell type
int([& int([&
@ -260,9 +170,7 @@ subroutine mesh_init(ip,el)
myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0) myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0)
call IO_open_inputFile(FILEUNIT,modelName) ! parse info from input file... call IO_open_inputFile(FILEUNIT,modelName)
if (myDebug) write(6,'(a)') ' Opened input file'; flush(6)
fileFormatVersion = mesh_marc_get_fileFormat(FILEUNIT) fileFormatVersion = mesh_marc_get_fileFormat(FILEUNIT)
if (myDebug) write(6,'(a)') ' Got input file format'; flush(6) if (myDebug) write(6,'(a)') ' Got input file format'; flush(6)
@ -329,14 +237,7 @@ subroutine mesh_init(ip,el)
call mesh_build_ipAreas call mesh_build_ipAreas
if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6)
call mesh_build_nodeTwins
if (myDebug) write(6,'(a)') ' Built node twins'; flush(6)
call mesh_build_sharedElems
if (myDebug) write(6,'(a)') ' Built shared elements'; flush(6)
call mesh_build_ipNeighborhood
call IP_neighborhood2 call IP_neighborhood2
if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6) if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6)
if (usePingPong .and. (mesh_Nelems /= theMesh%nElems)) & if (usePingPong .and. (mesh_Nelems /= theMesh%nElems)) &
@ -354,14 +255,11 @@ subroutine mesh_init(ip,el)
calcMode = .false. ! pretend to have collected what first call is asking (F = I) calcMode = .false. ! pretend to have collected what first call is asking (F = I)
calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc" calcMode(ip,mesh_FEasCP('elem',el)) = .true. ! first ip,el needs to be already pingponged to "calc"
theMesh%homogenizationAt = mesh_element(3,:)
theMesh%microstructureAt = mesh_element(4,:)
call discretization_init(mesh_element(3,:),mesh_element(4,:),& call discretization_init(mesh_element(3,:),mesh_element(4,:),&
reshape(mesh_ipCoordinates,[3,theMesh%elem%nIPs*theMesh%nElems]),& reshape(mesh_ipCoordinates,[3,theMesh%elem%nIPs*theMesh%nElems]),&
mesh_node0) mesh_node0)
call geometry_plastic_nonlocal_setIPvolume(mesh_ipVolume) call geometry_plastic_nonlocal_setIPvolume(mesh_ipVolume)
call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood) call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood2)
call geometry_plastic_nonlocal_setIParea(mesh_IParea) call geometry_plastic_nonlocal_setIParea(mesh_IParea)
call geometry_plastic_nonlocal_setIPareaNormal(mesh_IPareaNormal) call geometry_plastic_nonlocal_setIPareaNormal(mesh_IPareaNormal)
@ -396,7 +294,7 @@ integer function mesh_marc_get_fileFormat(fileUnit)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Figures out table styles for initial cond and hypoelastic !> @brief Figures out table styles for initial cond and hypoelastic
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_marc_get_tableStyles(initialcond, hypoelastic,fileUnit) subroutine mesh_marc_get_tableStyles(initialcond,hypoelastic,fileUnit)
integer, intent(out) :: initialcond, hypoelastic integer, intent(out) :: initialcond, hypoelastic
integer, intent(in) :: fileUnit integer, intent(in) :: fileUnit
@ -499,8 +397,8 @@ subroutine mesh_marc_count_nodesAndElements(nNodes, nElems, fileUnit)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_marc_count_elementSets(nElemSets,maxNelemInSet,fileUnit) subroutine mesh_marc_count_elementSets(nElemSets,maxNelemInSet,fileUnit)
integer, intent(in) :: fileUnit
integer, intent(out) :: nElemSets, maxNelemInSet integer, intent(out) :: nElemSets, maxNelemInSet
integer, intent(in) :: fileUnit
integer, allocatable, dimension(:) :: chunkPos integer, allocatable, dimension(:) :: chunkPos
character(len=300) :: line character(len=300) :: line
@ -528,9 +426,9 @@ subroutine mesh_marc_count_elementSets(nElemSets,maxNelemInSet,fileUnit)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_marc_map_elementSets(nameElemSet,mapElemSet,fileUnit) subroutine mesh_marc_map_elementSets(nameElemSet,mapElemSet,fileUnit)
integer, intent(in) :: fileUnit
character(len=64), dimension(:), intent(out) :: nameElemSet character(len=64), dimension(:), intent(out) :: nameElemSet
integer, dimension(:,:), intent(out) :: mapElemSet integer, dimension(:,:), intent(out) :: mapElemSet
integer, intent(in) :: fileUnit
integer, allocatable, dimension(:) :: chunkPos integer, allocatable, dimension(:) :: chunkPos
character(len=300) :: line character(len=300) :: line
@ -877,8 +775,8 @@ subroutine buildCells(thisMesh,elem,connectivity_elem)
integer,dimension(:), allocatable :: candidates_local integer,dimension(:), allocatable :: candidates_local
integer,dimension(:,:,:), allocatable :: connectivity_cell integer,dimension(:,:,:), allocatable :: connectivity_cell
integer,dimension(:,:), allocatable :: connectivity_cell_reshape integer,dimension(:,:), allocatable :: connectivity_cell_reshape
real(pReal), dimension(:,:), allocatable :: coordinates,nodes5 real(pReal), dimension(:,:), allocatable :: nodes_new,nodes
integer :: e, n, c, p, s,u,i,m,j,nParentNodes,nCellNode,ierr,Nelem,candidateID integer :: e, n, c, p, s,i,m,j,nParentNodes,nCellNode,Nelem,candidateID
Nelem = thisMesh%Nelems Nelem = thisMesh%Nelems
@ -976,7 +874,7 @@ subroutine buildCells(thisMesh,elem,connectivity_elem)
! calculate coordinates of cell nodes and insert their ID into the cell conectivity ! calculate coordinates of cell nodes and insert their ID into the cell conectivity
coordinates = reshape([(0.0_pReal,j = 1, 3*i)], [3,i]) nodes_new = reshape([(0.0_pReal,j = 1, 3*i)], [3,i])
i = 1 i = 1
n = 1 n = 1
@ -987,10 +885,10 @@ subroutine buildCells(thisMesh,elem,connectivity_elem)
e = candidates_global(nParentNodes*2+1,n+j) e = candidates_global(nParentNodes*2+1,n+j)
c = candidates_global(nParentNodes*2+2,n+j) c = candidates_global(nParentNodes*2+2,n+j)
do m = 1, nParentNodes do m = 1, nParentNodes
coordinates(:,i) = coordinates(:,i) & nodes_new(:,i) = nodes_new(:,i) &
+ thisMesh%node_0(:,parentsAndWeights(m,1)) * real(parentsAndWeights(m,2),pReal) + thisMesh%node_0(:,parentsAndWeights(m,1)) * real(parentsAndWeights(m,2),pReal)
enddo enddo
coordinates(:,i) = coordinates(:,i)/real(sum(parentsAndWeights(:,2)),pReal) nodes_new(:,i) = nodes_new(:,i)/real(sum(parentsAndWeights(:,2)),pReal)
do while (n+j<= size(candidates_local)*Nelem) do while (n+j<= size(candidates_local)*Nelem)
if (any(candidates_global(1:2*nParentNodes,n+j)/=candidates_global(1:2*nParentNodes,n))) exit if (any(candidates_global(1:2*nParentNodes,n+j)/=candidates_global(1:2*nParentNodes,n))) exit
@ -1005,11 +903,11 @@ subroutine buildCells(thisMesh,elem,connectivity_elem)
enddo enddo
nCellNode = nCellNode + i nCellNode = nCellNode + i
if (i/=0) nodes5 = reshape([nodes5,coordinates],[3,nCellNode]) if (i/=0) nodes = reshape([nodes,nodes_new],[3,nCellNode])
enddo enddo
thisMesh%node_0 = nodes5 thisMesh%node_0 = nodes
mesh_cell2 = connectivity_cell mesh_cell2 = connectivity_cell
connectivity_cell_reshape = reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*thisMesh%Nelems]) connectivity_cell_reshape = reshape(connectivity_cell,[elem%NcellNodesPerCell,elem%nIPs*Nelem])
#if defined(DAMASK_HDF5) #if defined(DAMASK_HDF5)
call results_openJobFile call results_openJobFile
@ -1017,6 +915,7 @@ subroutine buildCells(thisMesh,elem,connectivity_elem)
'connectivity of the cells','-') 'connectivity of the cells','-')
call results_closeJobFile call results_closeJobFile
#endif #endif
end subroutine buildCells end subroutine buildCells
@ -1252,17 +1151,6 @@ subroutine IP_neighborhood2
endif endif
f = f +1 f = f +1
enddo enddo
call geometry_plastic_nonlocal_setIPneighborhood(mesh_ipNeighborhood2)
do e = 1,theMesh%nElems
do i = 1,theMesh%elem%nIPs
do n = 1, theMesh%elem%nIPneighbors
write(6,'(a,i1.1,x,i1.1,x,i1.1)') 'e,i,n ',e,i,n
write(6,'(i1.1,x,i1.1,x,i3.2)') mesh_ipNeighborhood(1:3,n,i,e)
write(6,'(i1.1,x,i1.1,x,i3.2)') mesh_ipNeighborhood2(1:3,n,i,e)
enddo
enddo
enddo
end subroutine IP_neighborhood2 end subroutine IP_neighborhood2
@ -1369,394 +1257,6 @@ subroutine mesh_build_ipAreas
end subroutine mesh_build_ipAreas end subroutine mesh_build_ipAreas
!--------------------------------------------------------------------------------------------------
!> @brief assignment of twin nodes for each cp node, allocate globals '_nodeTwins'
!--------------------------------------------------------------------------------------------------
subroutine mesh_build_nodeTwins
integer dir, & ! direction of periodicity
node, &
minimumNode, &
maximumNode, &
n1, &
n2
integer, dimension(mesh_Nnodes+1) :: minimumNodes, maximumNodes ! list of surface nodes (minimum and maximum coordinate value) with first entry giving the number of nodes
real(pReal) minCoord, maxCoord, & ! extreme positions in one dimension
tolerance ! tolerance below which positions are assumed identical
real(pReal), dimension(3) :: distance ! distance between two nodes in all three coordinates
logical, dimension(mesh_Nnodes) :: unpaired
allocate(mesh_nodeTwins(3,mesh_Nnodes))
mesh_nodeTwins = 0
tolerance = 0.001_pReal * minval(mesh_ipVolume) ** 0.333_pReal
do dir = 1,3 ! check periodicity in directions of x,y,z
if (mesh_periodicSurface(dir)) then ! only if periodicity is requested
!*** find out which nodes sit on the surface
!*** and have a minimum or maximum position in this dimension
minimumNodes = 0
maximumNodes = 0
minCoord = minval(mesh_node0(dir,:))
maxCoord = maxval(mesh_node0(dir,:))
do node = 1,mesh_Nnodes ! loop through all nodes and find surface nodes
if (abs(mesh_node0(dir,node) - minCoord) <= tolerance) then
minimumNodes(1) = minimumNodes(1) + 1
minimumNodes(minimumNodes(1)+1) = node
elseif (abs(mesh_node0(dir,node) - maxCoord) <= tolerance) then
maximumNodes(1) = maximumNodes(1) + 1
maximumNodes(maximumNodes(1)+1) = node
endif
enddo
!*** find the corresponding node on the other side with the same position in this dimension
unpaired = .true.
do n1 = 1,minimumNodes(1)
minimumNode = minimumNodes(n1+1)
if (unpaired(minimumNode)) then
do n2 = 1,maximumNodes(1)
maximumNode = maximumNodes(n2+1)
distance = abs(mesh_node0(:,minimumNode) - mesh_node0(:,maximumNode))
if (sum(distance) - distance(dir) <= tolerance) then ! minimum possible distance (within tolerance)
mesh_nodeTwins(dir,minimumNode) = maximumNode
mesh_nodeTwins(dir,maximumNode) = minimumNode
unpaired(maximumNode) = .false. ! remember this node, we don't have to look for his partner again
exit
endif
enddo
endif
enddo
endif
enddo
end subroutine mesh_build_nodeTwins
!--------------------------------------------------------------------------------------------------
!> @brief get maximum count of shared elements among cpElements and build list of elements shared
!! by each node in mesh. Allocate globals '_maxNsharedElems' and '_sharedElem'
!--------------------------------------------------------------------------------------------------
subroutine mesh_build_sharedElems
integer(pint) e, & ! element index
g, & ! element type
node, & ! CP node index
n, & ! node index per element
myDim, & ! dimension index
nodeTwin ! node twin in the specified dimension
integer, dimension (mesh_Nnodes) :: node_count
integer, dimension(:), allocatable :: node_seen
allocate(node_seen(maxval(FE_NmatchingNodes)))
node_count = 0
do e = 1,theMesh%nElems
g = theMesh%elem%geomType
node_seen = 0 ! reset node duplicates
do n = 1,FE_NmatchingNodes(g) ! check each node of element
node = mesh_element(4+n,e)
if (all(node_seen /= node)) then
node_count(node) = node_count(node) + 1 ! if FE node not yet encountered -> count it
do myDim = 1,3 ! check in each dimension...
nodeTwin = mesh_nodeTwins(myDim,node)
if (nodeTwin > 0) & ! if I am a twin of some node...
node_count(nodeTwin) = node_count(nodeTwin) + 1 ! -> count me again for the twin node
enddo
endif
node_seen(n) = node ! remember this node to be counted already
enddo
enddo
mesh_maxNsharedElems = int(maxval(node_count),pInt) ! most shared node
allocate(mesh_sharedElem(1+mesh_maxNsharedElems,mesh_Nnodes),source=0)
do e = 1,theMesh%nElems
g = theMesh%elem%geomType
node_seen = 0
do n = 1,FE_NmatchingNodes(g)
node = mesh_element(4+n,e)
if (all(node_seen /= node)) then
mesh_sharedElem(1,node) = mesh_sharedElem(1,node) + 1 ! count for each node the connected elements
mesh_sharedElem(mesh_sharedElem(1,node)+1,node) = e ! store the respective element id
do myDim = 1,3 ! check in each dimension...
nodeTwin = mesh_nodeTwins(myDim,node)
if (nodeTwin > 0) then ! if i am a twin of some node...
mesh_sharedElem(1,nodeTwin) = mesh_sharedElem(1,nodeTwin) + 1 ! ...count me again for the twin
mesh_sharedElem(mesh_sharedElem(1,nodeTwin)+1,nodeTwin) = e ! store the respective element id
endif
enddo
endif
node_seen(n) = node
enddo
enddo
end subroutine mesh_build_sharedElems
!--------------------------------------------------------------------------------------------------
!> @brief build up of IP neighborhood, allocate globals '_ipNeighborhood'
!--------------------------------------------------------------------------------------------------
subroutine mesh_build_ipNeighborhood
integer :: myElem, & ! my CP element index
myIP, &
myType, & ! my element type
myFace, &
neighbor, & ! neighor index
neighboringIPkey, & ! positive integer indicating the neighboring IP (for intra-element) and negative integer indicating the face towards neighbor (for neighboring element)
candidateIP, &
neighboringType, & ! element type of neighbor
NlinkedNodes, & ! number of linked nodes
twin_of_linkedNode, & ! node twin of a specific linkedNode
NmatchingNodes, & ! number of matching nodes
dir, & ! direction of periodicity
matchingElem, & ! CP elem number of matching element
matchingFace, & ! face ID of matching element
a, anchor, &
neighboringIP, &
neighboringElem, &
pointingToMe
integer, dimension(FE_maxmaxNnodesAtIP) :: &
linkedNodes = 0, &
matchingNodes
logical checkTwins
allocate(mesh_ipNeighborhood(3,theMesh%elem%nIPneighbors,theMesh%elem%nIPs,theMesh%nElems))
mesh_ipNeighborhood = 0
do myElem = 1,theMesh%nElems ! loop over cpElems
myType = theMesh%elem%geomType
do myIP = 1,theMesh%elem%nIPs
do neighbor = 1,FE_NipNeighbors(theMesh%elem%cellType) ! loop over neighbors of IP
neighboringIPkey = theMesh%elem%IPneighbor(neighbor,myIP)
!*** if the key is positive, the neighbor is inside the element
!*** that means, we have already found our neighboring IP
if (neighboringIPkey > 0) then
mesh_ipNeighborhood(1,neighbor,myIP,myElem) = myElem
mesh_ipNeighborhood(2,neighbor,myIP,myElem) = neighboringIPkey
!*** if the key is negative, the neighbor resides in a neighboring element
!*** that means, we have to look through the face indicated by the key and see which element is behind that face
elseif (neighboringIPkey < 0) then ! neighboring element's IP
myFace = -neighboringIPkey
call mesh_faceMatch(myElem, myFace, matchingElem, matchingFace) ! get face and CP elem id of face match
if (matchingElem > 0) then ! found match?
neighboringType = theMesh%elem%geomType
!*** trivial solution if neighbor has only one IP
if (theMesh%elem%nIPs == 1) then
mesh_ipNeighborhood(1,neighbor,myIP,myElem) = matchingElem
mesh_ipNeighborhood(2,neighbor,myIP,myElem) = 1
cycle
endif
!*** find those nodes which build the link to the neighbor
NlinkedNodes = 0
linkedNodes = 0
do a = 1,theMesh%elem%maxNnodeAtIP
anchor = theMesh%elem%NnodeAtIP(a,myIP)
if (anchor /= 0) then ! valid anchor node
if (any(FE_face(:,myFace,myType) == anchor)) then ! ip anchor sits on face?
NlinkedNodes = NlinkedNodes + 1
linkedNodes(NlinkedNodes) = mesh_element(4+anchor,myElem) ! CP id of anchor node
else ! something went wrong with the linkage, since not all anchors sit on my face
NlinkedNodes = 0
linkedNodes = 0
exit
endif
endif
enddo
!*** loop through the ips of my neighbor
!*** and try to find an ip with matching nodes
!*** also try to match with node twins
checkCandidateIP: do candidateIP = 1,theMesh%elem%nIPs
NmatchingNodes = 0
matchingNodes = 0
do a = 1,theMesh%elem%maxNnodeAtIP
anchor = theMesh%elem%NnodeAtIP(a,candidateIP)
if (anchor /= 0) then ! valid anchor node
if (any(FE_face(:,matchingFace,neighboringType) == anchor)) then ! sits on matching face?
NmatchingNodes = NmatchingNodes + 1
matchingNodes(NmatchingNodes) = mesh_element(4+anchor,matchingElem) ! CP id of neighbor's anchor node
else ! no matching, because not all nodes sit on the matching face
NmatchingNodes = 0
matchingNodes = 0
exit
endif
endif
enddo
if (NmatchingNodes /= NlinkedNodes) & ! this ip has wrong count of anchors on face
cycle checkCandidateIP
!*** check "normal" nodes whether they match or not
checkTwins = .false.
do a = 1,NlinkedNodes
if (all(matchingNodes /= linkedNodes(a))) then ! this linkedNode does not match any matchingNode
checkTwins = .true.
exit ! no need to search further
endif
enddo
!*** if no match found, then also check node twins
if(checkTwins) then
dir = int(maxloc(abs(mesh_ipAreaNormal(1:3,neighbor,myIP,myElem)),1),pInt) ! check for twins only in direction of the surface normal
do a = 1,NlinkedNodes
twin_of_linkedNode = mesh_nodeTwins(dir,linkedNodes(a))
if (twin_of_linkedNode == 0 .or. & ! twin of linkedNode does not exist...
all(matchingNodes /= twin_of_linkedNode)) then ! ... or it does not match any matchingNode
cycle checkCandidateIP ! ... then check next candidateIP
endif
enddo
endif
!*** we found a match !!!
mesh_ipNeighborhood(1,neighbor,myIP,myElem) = matchingElem
mesh_ipNeighborhood(2,neighbor,myIP,myElem) = candidateIP
exit checkCandidateIP
enddo checkCandidateIP
endif ! end of valid external matching
endif ! end of internal/external matching
enddo
enddo
enddo
do myElem = 1,theMesh%nElems ! loop over cpElems
myType = theMesh%elem%geomType
do myIP = 1,theMesh%elem%nIPs
do neighbor = 1,FE_NipNeighbors(theMesh%elem%cellType) ! loop over neighbors of IP
neighboringElem = mesh_ipNeighborhood(1,neighbor,myIP,myElem)
neighboringIP = mesh_ipNeighborhood(2,neighbor,myIP,myElem)
if (neighboringElem > 0 .and. neighboringIP > 0) then ! if neighbor exists ...
neighboringType = theMesh%elem%geomType
do pointingToMe = 1,FE_NipNeighbors(theMesh%elem%cellType) ! find neighboring index that points from my neighbor to myself
if ( myElem == mesh_ipNeighborhood(1,pointingToMe,neighboringIP,neighboringElem) &
.and. myIP == mesh_ipNeighborhood(2,pointingToMe,neighboringIP,neighboringElem)) then ! possible candidate
if (math_mul3x3(mesh_ipAreaNormal(1:3,neighbor,myIP,myElem),&
mesh_ipAreaNormal(1:3,pointingToMe,neighboringIP,neighboringElem)) < 0.0_pReal) then ! area normals have opposite orientation (we have to check that because of special case for single element with two ips and periodicity. In this case the neighbor is identical in two different directions.)
mesh_ipNeighborhood(3,neighbor,myIP,myElem) = pointingToMe ! found match
exit ! so no need to search further
endif
endif
enddo
endif
enddo
enddo
enddo
contains
!--------------------------------------------------------------------------------------------------
!> @brief find face-matching element of same type
!--------------------------------------------------------------------------------------------------
subroutine mesh_faceMatch(elem, face ,matchingElem, matchingFace)
integer, intent(out) :: matchingElem, & ! matching CP element ID
matchingFace ! matching face ID
integer, intent(in) :: face, & ! face ID
elem ! CP elem ID
integer, dimension(FE_NmatchingNodesPerFace(face,theMesh%elem%geomType)) :: &
myFaceNodes ! global node ids on my face
integer :: myType, &
candidateType, &
candidateElem, &
candidateFace, &
candidateFaceNode, &
minNsharedElems, &
NsharedElems, &
lonelyNode = 0, &
i, &
n, &
dir ! periodicity direction
integer, dimension(:), allocatable :: element_seen
logical checkTwins
matchingElem = 0
matchingFace = 0
minNsharedElems = mesh_maxNsharedElems + 1 ! init to worst case
myType =theMesh%elem%geomType
do n = 1,FE_NmatchingNodesPerFace(face,myType) ! loop over nodes on face
myFaceNodes(n) = mesh_element(4+FE_face(n,face,myType),elem) ! CP id of face node
NsharedElems = mesh_sharedElem(1,myFaceNodes(n)) ! figure # shared elements for this node
if (NsharedElems < minNsharedElems) then
minNsharedElems = NsharedElems ! remember min # shared elems
lonelyNode = n ! remember most lonely node
endif
enddo
allocate(element_seen(minNsharedElems))
element_seen = 0
checkCandidate: do i = 1,minNsharedElems ! iterate over lonelyNode's shared elements
candidateElem = mesh_sharedElem(1+i,myFaceNodes(lonelyNode)) ! present candidate elem
if (all(element_seen /= candidateElem)) then ! element seen for the first time?
element_seen(i) = candidateElem
candidateType = theMesh%elem%geomType
checkCandidateFace: do candidateFace = 1,FE_maxNipNeighbors ! check each face of candidate
if (FE_NmatchingNodesPerFace(candidateFace,candidateType) &
/= FE_NmatchingNodesPerFace(face,myType) & ! incompatible face
.or. (candidateElem == elem .and. candidateFace == face)) then ! this is my face
cycle checkCandidateFace
endif
checkTwins = .false.
do n = 1,FE_NmatchingNodesPerFace(candidateFace,candidateType) ! loop through nodes on face
candidateFaceNode = mesh_element(4+FE_face(n,candidateFace,candidateType),candidateElem)
if (all(myFaceNodes /= candidateFaceNode)) then ! candidate node does not match any of my face nodes
checkTwins = .true. ! perhaps the twin nodes do match
exit
endif
enddo
if(checkTwins) then
checkCandidateFaceTwins: do dir = 1,3
do n = 1,FE_NmatchingNodesPerFace(candidateFace,candidateType) ! loop through nodes on face
candidateFaceNode = mesh_element(4+FE_face(n,candidateFace,candidateType),candidateElem)
if (all(myFaceNodes /= mesh_nodeTwins(dir,candidateFaceNode))) then ! node twin does not match either
if (dir == 3) then
cycle checkCandidateFace
else
cycle checkCandidateFaceTwins ! try twins in next dimension
endif
endif
enddo
exit checkCandidateFaceTwins
enddo checkCandidateFaceTwins
endif
matchingFace = candidateFace
matchingElem = candidateElem
exit checkCandidate ! found my matching candidate
enddo checkCandidateFace
endif
enddo checkCandidate
end subroutine mesh_faceMatch
end subroutine mesh_build_ipNeighborhood
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief get properties of different types of finite elements !> @brief get properties of different types of finite elements
!> @details assign globals FE_cellface !> @details assign globals FE_cellface