From 0b8ff64884f2c0509edefc520192deb2f5c64ecf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 30 Jan 2020 23:33:33 +0100 Subject: [PATCH] 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. --- src/CPFEM.f90 | 2 +- src/DAMASK_marc.f90 | 4 +- src/mesh_marc.f90 | 126 +++++++++++++------------------------------- 3 files changed, 41 insertions(+), 91 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 1ef5f1389..fd406219f 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -152,7 +152,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt i, j, k, l, m, n, ph, homog, mySource 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 & .and. elCP == debug_e .and. ip == debug_i) then diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index ba7751372..f7efb5279 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -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 computationMode = CPFEM_RESTOREJACOBIAN 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" terminallyIll = .false. 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) :: & 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 diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90 index 05bb77b8a..14e001eaf 100644 --- a/src/mesh_marc.f90 +++ b/src/mesh_marc.f90 @@ -32,17 +32,12 @@ module mesh real(pReal), public, protected :: & mesh_unitlength !< physical length of one unit in mesh - integer, dimension(:,:), allocatable, target :: & - mesh_mapFEtoCPelem, & !< [sorted FEid, corresponding CPid] - mesh_mapFEtoCPnode !< [sorted FEnode, corresponding CPnode] - - integer, dimension(:), allocatable :: & - mapMarc2DAMASK_elem, & !< DAMASK element ID for Marc element ID - mapMarc2DAMASK_node !< DAMASK node ID for Marc node ID + integer, dimension(:), allocatable, public :: & + mesh_FEM2DAMASK_elem, & !< DAMASK element ID for Marc element ID + mesh_FEM2DAMASK_node !< DAMASK node ID for Marc node ID public :: & - mesh_init, & - mesh_FEasCP + mesh_init contains @@ -89,7 +84,7 @@ subroutine mesh_init(ip,el) FEsolving_execIP = [1,elem%nIPs] 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)) @@ -215,13 +210,11 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,microstructureAt,homogeni call inputRead_elemType(elem, & nElems,inputFile) - allocate (mesh_mapFEtoCPelem(2,nElems), source=0) - call inputRead_mapElems(mapMarc2DAMASK_elem,& - elem%nNodes,inputFile) + call inputRead_mapElems(mesh_FEM2DAMASK_elem,& + nElems,elem%nNodes,inputFile) - allocate (mesh_mapFEtoCPnode(2,Nnodes), source=0) - call inputRead_mapNodes(mapMarc2DAMASK_node,& - inputFile) + call inputRead_mapNodes(mesh_FEM2DAMASK_node,& + nNodes,inputFile) call inputRead_elemNodes(node0_elem, & Nnodes,inputFile) @@ -432,14 +425,16 @@ end subroutine inputRead_mapElemSets !-------------------------------------------------------------------------------------------------- !> @brief Maps elements from FE ID to internal (consecutive) representation. !-------------------------------------------------------------------------------------------------- -subroutine inputRead_mapElems(map, & - nNodes,fileContent) +subroutine inputRead_mapElems(FEM2DAMASK, & + 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 + integer, dimension(2,nElems) :: map_unsorted integer, allocatable, dimension(:) :: chunkPos integer :: i,j,l,nNodesAlreadyRead @@ -448,11 +443,11 @@ subroutine inputRead_mapElems(map, & if(chunkPos(1) < 1) cycle if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'connectivity') then j = 0 - do i = 1,size(mesh_mapFEtoCPelem,2) + do i = 1,nElems 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 - 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 chunkPos = IO_stringPos(fileContent(l+1+i+j)) nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1) @@ -462,10 +457,10 @@ subroutine inputRead_mapElems(map, & endif enddo - call math_sort(mesh_mapFEtoCPelem) - allocate(map(minval(mesh_mapFEtoCPelem(1,:)):maxval(mesh_mapFEtoCPelem(1,:))),source=-1) - do i = 1,size(mesh_mapFEtoCPelem,2) - map(mesh_mapFEtoCPelem(1,i)) = mesh_mapFEtoCPelem(2,i) + call math_sort(map_unsorted) + allocate(FEM2DAMASK(minval(map_unsorted(1,:)):maxval(map_unsorted(1,:))),source=-1) + do i = 1, nElems + FEM2DAMASK(map_unsorted(1,i)) = map_unsorted(2,i) enddo end subroutine inputRead_mapElems @@ -474,13 +469,15 @@ end subroutine inputRead_mapElems !-------------------------------------------------------------------------------------------------- !> @brief Maps node from FE ID to internal (consecutive) representation. !-------------------------------------------------------------------------------------------------- -subroutine inputRead_mapNodes(map, & - fileContent) +subroutine inputRead_mapNodes(FEM2DAMASK, & + 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 + integer, dimension(2,nNodes) :: map_unsorted integer, allocatable, dimension(:) :: chunkPos integer :: i, l @@ -488,18 +485,18 @@ subroutine inputRead_mapNodes(map, & chunkPos = IO_stringPos(fileContent(l)) if(chunkPos(1) < 1) cycle 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)) - 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 exit endif enddo - call math_sort(mesh_mapFEtoCPnode) - allocate(map(minval(mesh_mapFEtoCPnode(1,:)):maxval(mesh_mapFEtoCPnode(1,:))),source=-1) - do i = 1,size(mesh_mapFEtoCPnode,2) - map(mesh_mapFEtoCPnode(1,i)) = mesh_mapFEtoCPnode(2,i) + call math_sort(map_unsorted) + allocate(FEM2DAMASK(minval(map_unsorted(1,:)):maxval(map_unsorted(1,:))),source=-1) + do i = 1, nNodes + FEM2DAMASK(map_unsorted(1,i)) = map_unsorted(2,i) enddo end subroutine inputRead_mapNodes @@ -526,7 +523,7 @@ subroutine inputRead_elemNodes(nodes, & if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates') then do i=1,nNode 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 nodes(j,m) = mesh_unitlength * IO_floatValue(fileContent(l+1+i),chunkPos,j+1) enddo @@ -650,11 +647,11 @@ function inputRead_connectivityElem(nElem,nNodes,fileContent) j = 0 do i = 1,nElem 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 do k = 1,chunkPos(1)-2 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 nNodesAlreadyRead = chunkPos(1) - 2 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)) do k = 1,chunkPos(1) 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 nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1) enddo @@ -717,7 +714,7 @@ subroutine inputRead_microstructureAndHomogenization(microstructureAt,homogeniza if (initialcondTableStyle == 2) m = m + 2 contInts = continuousIntValues(fileContent(l+k+m+1:),nElem,nameElemSet,mapElemSet,size(nameElemSet)) ! get affected elements 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 == 3) homogenizationAt(e) = myVal enddo @@ -1050,53 +1047,6 @@ function IPareaNormal(elem,nElem,connectivity,node) 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. !! First integer in array is counter