diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index 137d48562..7ddcfa06a 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -59,17 +59,16 @@ subroutine mesh_init(ip,el) microstructureAt, & homogenizationAt - logical :: myDebug integer :: j integer(C_INTPTR_T) :: & devNull, z, z_offset write(6,'(/,a)') ' <<<+- mesh init -+>>>' - myDebug = iand(debug_level(debug_mesh),debug_levelBasic) /= 0 - - call mesh_spectral_read_grid(grid,geomSize,microstructureAt,homogenizationAt) + call readGeom(grid,geomSize,microstructureAt,homogenizationAt) +!-------------------------------------------------------------------------------------------------- +! grid solver specific quantities if(worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)') call fftw_mpi_init @@ -86,33 +85,33 @@ subroutine mesh_init(ip,el) myGrid = [grid(1:2),grid3] mySize = [geomSize(1:2),size3] +!-------------------------------------------------------------------------------------------------- +! general discretization microstructureAt = microstructureAt(product(grid(1:2))*grid3Offset+1: & product(grid(1:2))*(grid3Offset+grid3)) ! reallocate/shrink in case of MPI homogenizationAt = homogenizationAt(product(grid(1:2))*grid3Offset+1: & product(grid(1:2))*(grid3Offset+grid3)) ! reallocate/shrink in case of MPI - - mesh_ipCoordinates = mesh_build_ipCoordinates(myGrid,mySize,grid3Offset) - if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) - - call geometry_plastic_nonlocal_setIPvolume( & - reshape([(product(mySize/real(myGrid,pReal)),j=1,product(myGrid))],[1,product(myGrid)])) - call geometry_plastic_nonlocal_setIParea(mesh_build_ipAreas(mySize,myGrid)) - call geometry_plastic_nonlocal_setIPareaNormal(mesh_build_ipNormals(product(myGrid))) - call geometry_plastic_nonlocal_setIPneighborhood(mesh_spectral_build_ipNeighborhood(myGrid)) - if (myDebug) write(6,'(a)') ' Built nonlocal geometry'; flush(6) - - if (debug_e < 1 .or. debug_e > product(myGrid)) & - call IO_error(602,ext_msg='element') ! selected element does not exist - if (debug_i /= 1) & - call IO_error(602,ext_msg='IP') ! selected element does not have requested IP + mesh_ipCoordinates = IPcoordinates(myGrid,mySize,grid3Offset) + call discretization_init(homogenizationAt,microstructureAt, & + reshape(mesh_ipCoordinates,[3,product(myGrid)]), & + Nodes(myGrid,mySize,grid3Offset)) FEsolving_execElem = [1,product(myGrid)] ! parallel loop bounds set to comprise all elements allocate(FEsolving_execIP(2,product(myGrid)),source=1) ! parallel loop bounds set to comprise the only IP - call discretization_init(homogenizationAt,microstructureAt, & - reshape(mesh_ipCoordinates,[3,product(myGrid)]), & - mesh_spectral_build_nodes(myGrid,mySize,grid3Offset)) +!-------------------------------------------------------------------------------------------------- +! geometry information required by the nonlocal CP model + call geometry_plastic_nonlocal_setIPvolume(reshape([(product(mySize/real(myGrid,pReal)),j=1,product(myGrid))], & + [1,product(myGrid)])) + call geometry_plastic_nonlocal_setIParea (cellEdgeArea(mySize,myGrid)) + call geometry_plastic_nonlocal_setIPareaNormal (cellEdgeNormal(product(myGrid))) + call geometry_plastic_nonlocal_setIPneighborhood(IPneighborhood(myGrid)) + +!-------------------------------------------------------------------------------------------------- +! sanity checks for debugging + if (debug_e < 1 .or. debug_e > product(myGrid)) call IO_error(602,ext_msg='element') ! selected element does not exist + if (debug_i /= 1) call IO_error(602,ext_msg='IP') ! selected IP does not exist end subroutine mesh_init @@ -122,20 +121,20 @@ end subroutine mesh_init !> @details important variables have an implicit "save" attribute. Therefore, this function is ! supposed to be called only once! !-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_read_grid(grid,geomSize,microstructure,homogenization) +subroutine readGeom(grid,geomSize,microstructure,homogenization) integer, dimension(3), intent(out) :: grid ! grid (for all processes!) real(pReal), dimension(3), intent(out) :: geomSize ! size (for all processes!) integer, dimension(:), intent(out), allocatable :: & - microstructure, & - homogenization + microstructure, & + homogenization - character(len=:), allocatable :: rawData - character(len=65536) :: line - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: h =- 1_pInt - integer(pInt) :: & - headerLength = -1_pInt, & !< length of header (in lines) + character(len=:), allocatable :: rawData + character(len=65536) :: line + integer, allocatable, dimension(:) :: chunkPos + integer :: & + h =- 1, & + headerLength = -1, & !< length of header (in lines) fileLength, & !< length of the geom file (in characters) fileUnit, & startPos, endPos, & @@ -146,15 +145,15 @@ subroutine mesh_spectral_read_grid(grid,geomSize,microstructure,homogenization) e, & !< "element", i.e. spectral collocation point i, j - grid = -1_pInt + grid = -1 geomSize = -1.0_pReal !-------------------------------------------------------------------------------------------------- -! read data as stream +! read raw data as stream inquire(file = trim(geometryFile), size=fileLength) open(newunit=fileUnit, file=trim(geometryFile), access='stream',& status='old', position='rewind', action='read',iostat=myStat) - if(myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=trim(geometryFile)) + if(myStat /= 0) call IO_error(100,ext_msg=trim(geometryFile)) allocate(character(len=fileLength)::rawData) read(fileUnit) rawData close(fileUnit) @@ -164,251 +163,244 @@ subroutine mesh_spectral_read_grid(grid,geomSize,microstructure,homogenization) endPos = index(rawData,new_line('')) if(endPos <= index(rawData,'head')) then startPos = len(rawData) - call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_read_grid') + call IO_error(error_ID=841, ext_msg='readGeom') else chunkPos = IO_stringPos(rawData(1:endPos)) - if (chunkPos(1) < 2_pInt) call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_read_grid') - headerLength = IO_intValue(rawData(1:endPos),chunkPos,1_pInt) - startPos = endPos + 1_pInt + if (chunkPos(1) < 2) call IO_error(error_ID=841, ext_msg='readGeom') + headerLength = IO_intValue(rawData(1:endPos),chunkPos,1) + startPos = endPos + 1 endif !-------------------------------------------------------------------------------------------------- ! read and interprete header l = 0 do while (l < headerLength .and. startPos < len(rawData)) - endPos = startPos + index(rawData(startPos:),new_line('')) - 1_pInt + endPos = startPos + index(rawData(startPos:),new_line('')) - 1 if (endPos < startPos) endPos = len(rawData) ! end of file without new line line = rawData(startPos:endPos) - startPos = endPos + 1_pInt - l = l + 1_pInt + startPos = endPos + 1 + l = l + 1 chunkPos = IO_stringPos(trim(line)) if (chunkPos(1) < 2) cycle ! need at least one keyword value pair - select case ( IO_lc(IO_StringValue(trim(line),chunkPos,1_pInt,.true.)) ) + select case ( IO_lc(IO_StringValue(trim(line),chunkPos,1,.true.)) ) case ('grid') if (chunkPos(1) > 6) then - do j = 2_pInt,6_pInt,2_pInt + do j = 2,6,2 select case (IO_lc(IO_stringValue(line,chunkPos,j))) case('a') - grid(1) = IO_intValue(line,chunkPos,j+1_pInt) + grid(1) = IO_intValue(line,chunkPos,j+1) case('b') - grid(2) = IO_intValue(line,chunkPos,j+1_pInt) + grid(2) = IO_intValue(line,chunkPos,j+1) case('c') - grid(3) = IO_intValue(line,chunkPos,j+1_pInt) + grid(3) = IO_intValue(line,chunkPos,j+1) end select enddo endif case ('size') if (chunkPos(1) > 6) then - do j = 2_pInt,6_pInt,2_pInt + do j = 2,6,2 select case (IO_lc(IO_stringValue(line,chunkPos,j))) case('x') - geomSize(1) = IO_floatValue(line,chunkPos,j+1_pInt) + geomSize(1) = IO_floatValue(line,chunkPos,j+1) case('y') - geomSize(2) = IO_floatValue(line,chunkPos,j+1_pInt) + geomSize(2) = IO_floatValue(line,chunkPos,j+1) case('z') - geomSize(3) = IO_floatValue(line,chunkPos,j+1_pInt) + geomSize(3) = IO_floatValue(line,chunkPos,j+1) end select enddo endif case ('homogenization') - if (chunkPos(1) > 1) h = IO_intValue(line,chunkPos,2_pInt) + if (chunkPos(1) > 1) h = IO_intValue(line,chunkPos,2) end select enddo !-------------------------------------------------------------------------------------------------- ! sanity checks - if(h < 1_pInt) & - call IO_error(error_ID = 842_pInt, ext_msg='homogenization (mesh_spectral_read_grid)') - if(any(grid < 1_pInt)) & - call IO_error(error_ID = 842_pInt, ext_msg='grid (mesh_spectral_read_grid)') + if(h < 1) & + call IO_error(error_ID = 842, ext_msg='homogenization (readGeom)') + if(any(grid < 1)) & + call IO_error(error_ID = 842, ext_msg='grid (readGeom)') if(any(geomSize < 0.0_pReal)) & - call IO_error(error_ID = 842_pInt, ext_msg='size (mesh_spectral_read_grid)') + call IO_error(error_ID = 842, ext_msg='size (readGeom)') allocate(microstructure(product(grid)), source = -1) ! too large in case of MPI (shrink later, not very elegant) allocate(homogenization(product(grid)), source = h) ! too large in case of MPI (shrink later, not very elegant) !-------------------------------------------------------------------------------------------------- ! read and interpret content - e = 1_pInt + e = 1 do while (startPos < len(rawData)) - endPos = startPos + index(rawData(startPos:),new_line('')) - 1_pInt + endPos = startPos + index(rawData(startPos:),new_line('')) - 1 if (endPos < startPos) endPos = len(rawData) ! end of file without new line line = rawData(startPos:endPos) - startPos = endPos + 1_pInt - l = l + 1_pInt + startPos = endPos + 1 + l = l + 1 chunkPos = IO_stringPos(trim(line)) noCompression: if (chunkPos(1) /= 3) then c = chunkPos(1) - microstructure(e:e+c-1_pInt) = [(IO_intValue(line,chunkPos,i+1_pInt), i=0_pInt, c-1_pInt)] + microstructure(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)] else noCompression compression: if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'of') then c = IO_intValue(line,chunkPos,1) - microstructure(e:e+c-1_pInt) = [(IO_intValue(line,chunkPos,3),i = 1_pInt,IO_intValue(line,chunkPos,1))] + microstructure(e:e+c-1) = [(IO_intValue(line,chunkPos,3),i = 1,IO_intValue(line,chunkPos,1))] else if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'to') then compression - c = abs(IO_intValue(line,chunkPos,3) - IO_intValue(line,chunkPos,1)) + 1_pInt - o = merge(+1_pInt, -1_pInt, IO_intValue(line,chunkPos,3) > IO_intValue(line,chunkPos,1)) - microstructure(e:e+c-1_pInt) = [(i, i = IO_intValue(line,chunkPos,1),IO_intValue(line,chunkPos,3),o)] + c = abs(IO_intValue(line,chunkPos,3) - IO_intValue(line,chunkPos,1)) + 1 + o = merge(+1, -1, IO_intValue(line,chunkPos,3) > IO_intValue(line,chunkPos,1)) + microstructure(e:e+c-1) = [(i, i = IO_intValue(line,chunkPos,1),IO_intValue(line,chunkPos,3),o)] else compression c = chunkPos(1) - microstructure(e:e+c-1_pInt) = [(IO_intValue(line,chunkPos,i+1_pInt), i=0_pInt, c-1_pInt)] + microstructure(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)] endif compression endif noCompression e = e+c end do - if (e-1 /= product(grid)) call IO_error(error_ID = 843_pInt, el=e) + if (e-1 /= product(grid)) call IO_error(error_ID = 843, el=e) -end subroutine mesh_spectral_read_grid +end subroutine readGeom !--------------------------------------------------------------------------------------------------- -!> @brief Calculates position of nodes (pretend to be an element) +!> @brief Calculate position of IPs/cell centres (pretend to be an element) !--------------------------------------------------------------------------------------------------- -pure function mesh_spectral_build_nodes(grid,geomSize,grid3Offset) result(nodes) +function IPcoordinates(grid,geomSize,grid3Offset) - integer, dimension(3), intent(in) :: grid ! grid (for this process!) - real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!) - integer, intent(in) :: grid3Offset ! grid(3) offset - real(pReal), dimension(3,product(grid+1)) :: nodes - integer :: n,a,b,c + integer, dimension(3), intent(in) :: grid ! grid (for this process!) + real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!) + integer, intent(in) :: grid3Offset ! grid(3) offset - n = 0 - do c = 0, grid3 - do b = 0, grid(2) - do a = 0, grid(1) - n = n + 1 - nodes(1:3,n) = geomSize/real(grid,pReal) * real([a,b,grid3Offset+c],pReal) - enddo - enddo - enddo - -end function mesh_spectral_build_nodes - - -!--------------------------------------------------------------------------------------------------- -!> @brief Calculates position of IPs/cell centres (pretend to be an element) -!--------------------------------------------------------------------------------------------------- -function mesh_build_ipCoordinates(grid,geomSize,grid3Offset) result(ipCoordinates) - - integer, dimension(3), intent(in) :: grid ! grid (for this process!) - real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!) - integer, intent(in) :: grid3Offset ! grid(3) offset real(pReal), dimension(3,1,product(grid)) :: ipCoordinates - integer :: n,a,b,c + + integer :: & + a,b,c, & + i + i = 0 + do c = 1, grid(3); do b = 1, grid(2); do a = 1, grid(1) + i = i + 1 + IPcoordinates(1:3,1,i) = geomSize/real(grid,pReal) * (real([a,b,grid3Offset+c],pReal) -0.5_pReal) + enddo; enddo; enddo + +end function IPcoordinates + + +!--------------------------------------------------------------------------------------------------- +!> @brief Calculate position of nodes (pretend to be an element) +!--------------------------------------------------------------------------------------------------- +pure function nodes(grid,geomSize,grid3Offset) + + integer, dimension(3), intent(in) :: grid ! grid (for this process!) + real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!) + integer, intent(in) :: grid3Offset ! grid(3) offset + + real(pReal), dimension(3,product(grid+1)) :: nodes + + integer :: & + a,b,c, & + n + n = 0 - do c = 1, grid(3) - do b = 1, grid(2) - do a = 1, grid(1) - n = n + 1 - ipCoordinates(1:3,1,n) = geomSize/real(grid,pReal) * (real([a,b,grid3Offset+c],pReal) -0.5_pReal) - enddo - enddo - enddo + do c = 0, grid3; do b = 0, grid(2); do a = 0, grid(1) + n = n + 1 + nodes(1:3,n) = geomSize/real(grid,pReal) * real([a,b,grid3Offset+c],pReal) + enddo; enddo; enddo -end function mesh_build_ipCoordinates +end function nodes !-------------------------------------------------------------------------------------------------- -!> @brief build neighborhood relations for spectral -!> @details assign globals: mesh_ipNeighborhood +!> @brief Calculate IP interface areas !-------------------------------------------------------------------------------------------------- -pure function mesh_spectral_build_ipNeighborhood(grid) result(IPneighborhood) - - integer, dimension(3), intent(in) :: grid ! grid (for this process!) - - integer, dimension(3,6,1,product(grid)) :: IPneighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] - - integer :: & - x,y,z, & - e - - e = 0 - do z = 0,grid(3)-1 - do y = 0,grid(2)-1 - do x = 0,grid(1)-1 - e = e + 1 - IPneighborhood(1,1,1,e) = z * grid(1) * grid(2) & - + y * grid(1) & - + modulo(x+1,grid(1)) & - + 1 - IPneighborhood(1,2,1,e) = z * grid(1) * grid(2) & - + y * grid(1) & - + modulo(x-1,grid(1)) & - + 1 - IPneighborhood(1,3,1,e) = z * grid(1) * grid(2) & - + modulo(y+1,grid(2)) * grid(1) & - + x & - + 1 - IPneighborhood(1,4,1,e) = z * grid(1) * grid(2) & - + modulo(y-1,grid(2)) * grid(1) & - + x & - + 1 - IPneighborhood(1,5,1,e) = modulo(z+1,grid(3)) * grid(1) * grid(2) & - + y * grid(1) & - + x & - + 1 - IPneighborhood(1,6,1,e) = modulo(z-1,grid(3)) * grid(1) * grid(2) & - + y * grid(1) & - + x & - + 1 - IPneighborhood(2,1:6,1,e) = 1 - IPneighborhood(3,1,1,e) = 2 - IPneighborhood(3,2,1,e) = 1 - IPneighborhood(3,3,1,e) = 4 - IPneighborhood(3,4,1,e) = 3 - IPneighborhood(3,5,1,e) = 6 - IPneighborhood(3,6,1,e) = 5 - enddo - enddo - enddo - -end function mesh_spectral_build_ipNeighborhood - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculation of IP interface areas -!-------------------------------------------------------------------------------------------------- -pure function mesh_build_ipAreas(geomSize,grid) result(IPareas) +pure function cellEdgeArea(geomSize,grid) - real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!) - integer, dimension(3), intent(in) :: grid ! grid (for this process!) + real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!) + integer, dimension(3), intent(in) :: grid ! grid (for this process!) - real(pReal), dimension(6,1,product(grid)) :: IPareas + real(pReal), dimension(6,1,product(grid)) :: cellEdgeArea - IPareas(1:2,1,:) = geomSize(2)/real(grid(2)) * geomSize(3)/real(grid(3)) - IPareas(3:4,1,:) = geomSize(3)/real(grid(3)) * geomSize(1)/real(grid(1)) - IPareas(5:6,1,:) = geomSize(1)/real(grid(1)) * geomSize(2)/real(grid(2)) + cellEdgeArea(1:2,1,:) = geomSize(2)/real(grid(2)) * geomSize(3)/real(grid(3)) + cellEdgeArea(3:4,1,:) = geomSize(3)/real(grid(3)) * geomSize(1)/real(grid(1)) + cellEdgeArea(5:6,1,:) = geomSize(1)/real(grid(1)) * geomSize(2)/real(grid(2)) -end function mesh_build_ipAreas +end function cellEdgeArea !-------------------------------------------------------------------------------------------------- -!> @brief calculation of IP interface areas normals +!> @brief Calculate IP interface areas normals !-------------------------------------------------------------------------------------------------- -pure function mesh_build_ipNormals(nElems) result(IPnormals) +pure function cellEdgeNormal(nElems) integer, intent(in) :: nElems - real, dimension(3,6,1,nElems) :: IPnormals + real, dimension(3,6,1,nElems) :: cellEdgeNormal - IPnormals(1:3,1,1,:) = spread([+1.0_pReal, 0.0_pReal, 0.0_pReal],2,nElems) - IPnormals(1:3,2,1,:) = spread([-1.0_pReal, 0.0_pReal, 0.0_pReal],2,nElems) - IPnormals(1:3,3,1,:) = spread([ 0.0_pReal,+1.0_pReal, 0.0_pReal],2,nElems) - IPnormals(1:3,4,1,:) = spread([ 0.0_pReal,-1.0_pReal, 0.0_pReal],2,nElems) - IPnormals(1:3,5,1,:) = spread([ 0.0_pReal, 0.0_pReal,+1.0_pReal],2,nElems) - IPnormals(1:3,6,1,:) = spread([ 0.0_pReal, 0.0_pReal,-1.0_pReal],2,nElems) + cellEdgeNormal(1:3,1,1,:) = spread([+1.0_pReal, 0.0_pReal, 0.0_pReal],2,nElems) + cellEdgeNormal(1:3,2,1,:) = spread([-1.0_pReal, 0.0_pReal, 0.0_pReal],2,nElems) + cellEdgeNormal(1:3,3,1,:) = spread([ 0.0_pReal,+1.0_pReal, 0.0_pReal],2,nElems) + cellEdgeNormal(1:3,4,1,:) = spread([ 0.0_pReal,-1.0_pReal, 0.0_pReal],2,nElems) + cellEdgeNormal(1:3,5,1,:) = spread([ 0.0_pReal, 0.0_pReal,+1.0_pReal],2,nElems) + cellEdgeNormal(1:3,6,1,:) = spread([ 0.0_pReal, 0.0_pReal,-1.0_pReal],2,nElems) -end function mesh_build_ipNormals +end function cellEdgeNormal +!-------------------------------------------------------------------------------------------------- +!> @brief Build IP neighborhood relations +!-------------------------------------------------------------------------------------------------- +pure function IPneighborhood(grid) + + integer, dimension(3), intent(in) :: grid ! grid (for this process!) + + integer, dimension(3,6,1,product(grid)) :: IPneighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] + + integer :: & + x,y,z, & + e + + e = 0 + do z = 0,grid(3)-1; do y = 0,grid(2)-1; do x = 0,grid(1)-1 + e = e + 1 + IPneighborhood(1,1,1,e) = z * grid(1) * grid(2) & + + y * grid(1) & + + modulo(x+1,grid(1)) & + + 1 + IPneighborhood(1,2,1,e) = z * grid(1) * grid(2) & + + y * grid(1) & + + modulo(x-1,grid(1)) & + + 1 + IPneighborhood(1,3,1,e) = z * grid(1) * grid(2) & + + modulo(y+1,grid(2)) * grid(1) & + + x & + + 1 + IPneighborhood(1,4,1,e) = z * grid(1) * grid(2) & + + modulo(y-1,grid(2)) * grid(1) & + + x & + + 1 + IPneighborhood(1,5,1,e) = modulo(z+1,grid(3)) * grid(1) * grid(2) & + + y * grid(1) & + + x & + + 1 + IPneighborhood(1,6,1,e) = modulo(z-1,grid(3)) * grid(1) * grid(2) & + + y * grid(1) & + + x & + + 1 + IPneighborhood(2,1:6,1,e) = 1 + IPneighborhood(3,1,1,e) = 2 + IPneighborhood(3,2,1,e) = 1 + IPneighborhood(3,3,1,e) = 4 + IPneighborhood(3,4,1,e) = 3 + IPneighborhood(3,5,1,e) = 6 + IPneighborhood(3,6,1,e) = 5 + enddo; enddo; enddo + +end function IPneighborhood !--------------------------------------------------------------------------------------------------