store mapping MARC/FEM2DAMASK mapping

do not calculate the mapping for elements and nodes per call on
the fly, rather store it. Not memory efficient in the case that numbers
are not consequtive (order does not matter, but missing nodes/elements
would waste some 2 integers per missing number). However, this seem to cause
problems anyway when range indicators like '1 to 10' are used.
This commit is contained in:
Martin Diehl 2020-01-30 23:33:33 +01:00
parent d54b8714e1
commit 0b8ff64884
3 changed files with 41 additions and 91 deletions

View File

@ -152,7 +152,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
i, j, k, l, m, n, ph, homog, mySource i, j, k, l, m, n, ph, homog, mySource
logical updateJaco ! flag indicating if JAcobian has to be updated logical updateJaco ! flag indicating if JAcobian has to be updated
elCP = mesh_FEasCP('elem',elFE) elCP = mesh_FEM2DAMASK_elem(elFE)
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt & if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt &
.and. elCP == debug_e .and. ip == debug_i) then .and. elCP == debug_e .and. ip == debug_i) then

View File

@ -269,7 +269,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) & ! first after cutback if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) & ! first after cutback
computationMode = CPFEM_RESTOREJACOBIAN computationMode = CPFEM_RESTOREJACOBIAN
elseif (lovl == 6) then ! stress requested by marc elseif (lovl == 6) then ! stress requested by marc
cp_en = mesh_FEasCP('elem',m(1)) cp_en = mesh_FEM2DAMASK_elem(m(1))
if (cptim > theTime .or. inc /= theInc) then ! reached "convergence" if (cptim > theTime .or. inc /= theInc) then ! reached "convergence"
terminallyIll = .false. terminallyIll = .false.
cycleCounter = -1 ! first calc step increments this to cycle = 0 cycleCounter = -1 ! first calc step increments this to cycle = 0
@ -391,7 +391,7 @@ subroutine flux(f,ts,n,time)
real(pReal), dimension(2), intent(out) :: & real(pReal), dimension(2), intent(out) :: &
f f
call thermal_conduction_getSourceAndItsTangent(f(1), f(2), ts(3), n(3),mesh_FEasCP('elem',n(1))) call thermal_conduction_getSourceAndItsTangent(f(1), f(2), ts(3), n(3),mesh_FEM2DAMASK_elem(n(1)))
end subroutine flux end subroutine flux

View File

@ -32,17 +32,12 @@ module mesh
real(pReal), public, protected :: & real(pReal), public, protected :: &
mesh_unitlength !< physical length of one unit in mesh mesh_unitlength !< physical length of one unit in mesh
integer, dimension(:,:), allocatable, target :: & integer, dimension(:), allocatable, public :: &
mesh_mapFEtoCPelem, & !< [sorted FEid, corresponding CPid] mesh_FEM2DAMASK_elem, & !< DAMASK element ID for Marc element ID
mesh_mapFEtoCPnode !< [sorted FEnode, corresponding CPnode] mesh_FEM2DAMASK_node !< DAMASK node ID for Marc node ID
integer, dimension(:), allocatable :: &
mapMarc2DAMASK_elem, & !< DAMASK element ID for Marc element ID
mapMarc2DAMASK_node !< DAMASK node ID for Marc node ID
public :: & public :: &
mesh_init, & mesh_init
mesh_FEasCP
contains contains
@ -89,7 +84,7 @@ subroutine mesh_init(ip,el)
FEsolving_execIP = [1,elem%nIPs] FEsolving_execIP = [1,elem%nIPs]
allocate(calcMode(elem%nIPs,nElems),source=.false.) ! pretend to have collected what first call is asking (F = I) allocate(calcMode(elem%nIPs,nElems),source=.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_FEM2DAMASK_elem(el)) = .true. ! first ip,el needs to be already pingponged to "calc"
allocate(cellNodeDefinition(elem%nNodes-1)) allocate(cellNodeDefinition(elem%nNodes-1))
@ -215,13 +210,11 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni
call inputRead_elemType(elem, & call inputRead_elemType(elem, &
nElems,inputFile) nElems,inputFile)
allocate (mesh_mapFEtoCPelem(2,nElems), source=0) call inputRead_mapElems(mesh_FEM2DAMASK_elem,&
call inputRead_mapElems(mapMarc2DAMASK_elem,& nElems,elem%nNodes,inputFile)
elem%nNodes,inputFile)
allocate (mesh_mapFEtoCPnode(2,Nnodes), source=0) call inputRead_mapNodes(mesh_FEM2DAMASK_node,&
call inputRead_mapNodes(mapMarc2DAMASK_node,& nNodes,inputFile)
inputFile)
call inputRead_elemNodes(node0_elem, & call inputRead_elemNodes(node0_elem, &
Nnodes,inputFile) Nnodes,inputFile)
@ -432,14 +425,16 @@ end subroutine inputRead_mapElemSets
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Maps elements from FE ID to internal (consecutive) representation. !> @brief Maps elements from FE ID to internal (consecutive) representation.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine inputRead_mapElems(map, & subroutine inputRead_mapElems(FEM2DAMASK, &
nNodes,fileContent) nElems,nNodesPerElem,fileContent)
integer, allocatable, dimension(:), intent(out) :: map integer, allocatable, dimension(:), intent(out) :: FEM2DAMASK
integer, intent(in) :: nNodes !< number of nodes per element integer, intent(in) :: nElems, & !< number of elements
nNodesPerElem !< number of nodes per element
character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, dimension(2,nElems) :: map_unsorted
integer, allocatable, dimension(:) :: chunkPos integer, allocatable, dimension(:) :: chunkPos
integer :: i,j,l,nNodesAlreadyRead integer :: i,j,l,nNodesAlreadyRead
@ -448,11 +443,11 @@ subroutine inputRead_mapElems(map, &
if(chunkPos(1) < 1) cycle if(chunkPos(1) < 1) cycle
if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'connectivity') then if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'connectivity') then
j = 0 j = 0
do i = 1,size(mesh_mapFEtoCPelem,2) do i = 1,nElems
chunkPos = IO_stringPos(fileContent(l+1+i+j)) chunkPos = IO_stringPos(fileContent(l+1+i+j))
mesh_mapFEtoCPelem(:,i) = [IO_intValue(fileContent(l+1+i+j),chunkPos,1),i] map_unsorted(:,i) = [IO_intValue(fileContent(l+1+i+j),chunkPos,1),i]
nNodesAlreadyRead = chunkPos(1) - 2 nNodesAlreadyRead = chunkPos(1) - 2
do while(nNodesAlreadyRead < nNodes) ! read on if not all nodes in one line do while(nNodesAlreadyRead < nNodesPerElem) ! read on if not all nodes in one line
j = j + 1 j = j + 1
chunkPos = IO_stringPos(fileContent(l+1+i+j)) chunkPos = IO_stringPos(fileContent(l+1+i+j))
nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1) nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1)
@ -462,10 +457,10 @@ subroutine inputRead_mapElems(map, &
endif endif
enddo enddo
call math_sort(mesh_mapFEtoCPelem) call math_sort(map_unsorted)
allocate(map(minval(mesh_mapFEtoCPelem(1,:)):maxval(mesh_mapFEtoCPelem(1,:))),source=-1) allocate(FEM2DAMASK(minval(map_unsorted(1,:)):maxval(map_unsorted(1,:))),source=-1)
do i = 1,size(mesh_mapFEtoCPelem,2) do i = 1, nElems
map(mesh_mapFEtoCPelem(1,i)) = mesh_mapFEtoCPelem(2,i) FEM2DAMASK(map_unsorted(1,i)) = map_unsorted(2,i)
enddo enddo
end subroutine inputRead_mapElems end subroutine inputRead_mapElems
@ -474,13 +469,15 @@ end subroutine inputRead_mapElems
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Maps node from FE ID to internal (consecutive) representation. !> @brief Maps node from FE ID to internal (consecutive) representation.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine inputRead_mapNodes(map, & subroutine inputRead_mapNodes(FEM2DAMASK, &
fileContent) nNodes,fileContent)
integer, allocatable, dimension(:), intent(out) :: map integer, allocatable, dimension(:), intent(out) :: FEM2DAMASK
integer, intent(in) :: nNodes !< number of nodes
character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines character(len=*), dimension(:), intent(in) :: fileContent !< file content, separated per lines
integer, dimension(2,nNodes) :: map_unsorted
integer, allocatable, dimension(:) :: chunkPos integer, allocatable, dimension(:) :: chunkPos
integer :: i, l integer :: i, l
@ -488,18 +485,18 @@ subroutine inputRead_mapNodes(map, &
chunkPos = IO_stringPos(fileContent(l)) chunkPos = IO_stringPos(fileContent(l))
if(chunkPos(1) < 1) cycle if(chunkPos(1) < 1) cycle
if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates') then if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates') then
do i = 1,size(mesh_mapFEtoCPnode,2) do i = 1,nNodes
chunkPos = IO_stringPos(fileContent(l+1+i)) chunkPos = IO_stringPos(fileContent(l+1+i))
mesh_mapFEtoCPnode(:,i) = [IO_intValue(fileContent(l+1+i),chunkPos,1),i] map_unsorted(:,i) = [IO_intValue(fileContent(l+1+i),chunkPos,1),i]
enddo enddo
exit exit
endif endif
enddo enddo
call math_sort(mesh_mapFEtoCPnode) call math_sort(map_unsorted)
allocate(map(minval(mesh_mapFEtoCPnode(1,:)):maxval(mesh_mapFEtoCPnode(1,:))),source=-1) allocate(FEM2DAMASK(minval(map_unsorted(1,:)):maxval(map_unsorted(1,:))),source=-1)
do i = 1,size(mesh_mapFEtoCPnode,2) do i = 1, nNodes
map(mesh_mapFEtoCPnode(1,i)) = mesh_mapFEtoCPnode(2,i) FEM2DAMASK(map_unsorted(1,i)) = map_unsorted(2,i)
enddo enddo
end subroutine inputRead_mapNodes end subroutine inputRead_mapNodes
@ -526,7 +523,7 @@ subroutine inputRead_elemNodes(nodes, &
if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates') then if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates') then
do i=1,nNode do i=1,nNode
chunkPos = IO_stringPos(fileContent(l+1+i)) chunkPos = IO_stringPos(fileContent(l+1+i))
m = mesh_FEasCP('node',IO_intValue(fileContent(l+1+i),chunkPos,1)) m = mesh_FEM2DAMASK_node(IO_intValue(fileContent(l+1+i),chunkPos,1))
do j = 1,3 do j = 1,3
nodes(j,m) = mesh_unitlength * IO_floatValue(fileContent(l+1+i),chunkPos,j+1) nodes(j,m) = mesh_unitlength * IO_floatValue(fileContent(l+1+i),chunkPos,j+1)
enddo enddo
@ -650,11 +647,11 @@ function inputRead_connectivityElem(nElem,nNodes,fileContent)
j = 0 j = 0
do i = 1,nElem do i = 1,nElem
chunkPos = IO_stringPos(fileContent(l+1+i+j)) chunkPos = IO_stringPos(fileContent(l+1+i+j))
e = mesh_FEasCP('elem',IO_intValue(fileContent(l+1+i+j),chunkPos,1)) e = mesh_FEM2DAMASK_elem(IO_intValue(fileContent(l+1+i+j),chunkPos,1))
if (e /= 0) then ! disregard non CP elems if (e /= 0) then ! disregard non CP elems
do k = 1,chunkPos(1)-2 do k = 1,chunkPos(1)-2
inputRead_connectivityElem(k,e) = & inputRead_connectivityElem(k,e) = &
mesh_FEasCP('node',IO_IntValue(fileContent(l+1+i+j),chunkPos,k+2)) mesh_FEM2DAMASK_node(IO_IntValue(fileContent(l+1+i+j),chunkPos,k+2))
enddo enddo
nNodesAlreadyRead = chunkPos(1) - 2 nNodesAlreadyRead = chunkPos(1) - 2
do while(nNodesAlreadyRead < nNodes) ! read on if not all nodes in one line do while(nNodesAlreadyRead < nNodes) ! read on if not all nodes in one line
@ -662,7 +659,7 @@ function inputRead_connectivityElem(nElem,nNodes,fileContent)
chunkPos = IO_stringPos(fileContent(l+1+i+j)) chunkPos = IO_stringPos(fileContent(l+1+i+j))
do k = 1,chunkPos(1) do k = 1,chunkPos(1)
inputRead_connectivityElem(nNodesAlreadyRead+k,e) = & inputRead_connectivityElem(nNodesAlreadyRead+k,e) = &
mesh_FEasCP('node',IO_IntValue(fileContent(l+1+i+j),chunkPos,k)) mesh_FEM2DAMASK_node(IO_IntValue(fileContent(l+1+i+j),chunkPos,k))
enddo enddo
nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1) nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1)
enddo enddo
@ -717,7 +714,7 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza
if (initialcondTableStyle == 2) m = m + 2 if (initialcondTableStyle == 2) m = m + 2
contInts = continuousIntValues(fileContent(l+k+m+1:),nElem,nameElemSet,mapElemSet,size(nameElemSet)) ! get affected elements contInts = continuousIntValues(fileContent(l+k+m+1:),nElem,nameElemSet,mapElemSet,size(nameElemSet)) ! get affected elements
do i = 1,contInts(1) do i = 1,contInts(1)
e = mesh_FEasCP('elem',contInts(1+i)) e = mesh_FEM2DAMASK_elem(contInts(1+i))
if (sv == 2) microstructureAt(e) = myVal if (sv == 2) microstructureAt(e) = myVal
if (sv == 3) homogenizationAt(e) = myVal if (sv == 3) homogenizationAt(e) = myVal
enddo enddo
@ -1050,53 +1047,6 @@ function IPareaNormal(elem,nElem,connectivity,node)
end function IPareaNormal end function IPareaNormal
!--------------------------------------------------------------------------------------------------
!> @brief Gives the FE to CP ID mapping by binary search through lookup array
!! valid questions (what) are 'elem', 'node'
!--------------------------------------------------------------------------------------------------
integer function mesh_FEasCP(what,myID)
character(len=*), intent(in) :: what
integer, intent(in) :: myID
integer, dimension(:,:), pointer :: lookupMap
integer :: lower,upper,center
mesh_FEasCP = 0
select case(IO_lc(what(1:4)))
case('elem')
lookupMap => mesh_mapFEtoCPelem
case('node')
lookupMap => mesh_mapFEtoCPnode
case default
return
endselect
lower = 1
upper = int(size(lookupMap,2),pInt)
if (lookupMap(1,lower) == myID) then ! check at bounds QUESTION is it valid to extend bounds by 1 and just do binary search w/o init check at bounds?
mesh_FEasCP = lookupMap(2,lower)
return
elseif (lookupMap(1,upper) == myID) then
mesh_FEasCP = lookupMap(2,upper)
return
endif
binarySearch: do while (upper-lower > 1)
center = (lower+upper)/2
if (lookupMap(1,center) < myID) then
lower = center
elseif (lookupMap(1,center) > myID) then
upper = center
else
mesh_FEasCP = lookupMap(2,center)
exit
endif
enddo binarySearch
end function mesh_FEasCP
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return integer list corresponding to items in consecutive lines. !> @brief return integer list corresponding to items in consecutive lines.
!! First integer in array is counter !! First integer in array is counter