!-------------------------------------------------------------------------------------------------- !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief Parse geometry file to set up discretization and geometry for nonlocal model !-------------------------------------------------------------------------------------------------- module mesh #include use PETScsys use prec use system_routines use DAMASK_interface use IO use debug use numerics use discretization use geometry_plastic_nonlocal use FEsolving implicit none private real(pReal), dimension(:,:,:), allocatable, public :: & mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) integer, dimension(3), public, protected :: & grid !< (global) grid integer, public, protected :: & grid3, & !< (local) grid in 3rd direction grid3Offset !< (local) grid offset in 3rd direction real(pReal), dimension(3), public, protected :: & geomSize real(pReal), public, protected :: & size3, & !< (local) size in 3rd direction size3offset !< (local) size offset in 3rd direction public :: & mesh_init contains !-------------------------------------------------------------------------------------------------- !> @brief reads the geometry file to obtain information on discretization !-------------------------------------------------------------------------------------------------- subroutine mesh_init(ip,el) integer, intent(in), optional :: el, ip ! for compatibility reasons include 'fftw3-mpi.f03' real(pReal), dimension(3) :: & mySize !< domain size of this process integer, dimension(3) :: & myGrid !< domain grid of this process integer, dimension(:), allocatable :: & microstructureAt, & homogenizationAt integer :: j integer(C_INTPTR_T) :: & devNull, z, z_offset write(6,'(/,a)') ' <<<+- mesh init -+>>>' 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 devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T), & int(grid(2),C_INTPTR_T), & int(grid(1),C_INTPTR_T)/2+1, & PETSC_COMM_WORLD, & z, & ! domain grid size along z z_offset) ! domain grid offset along z grid3 = int(z) grid3Offset = int(z_offset) size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal) size3Offset = geomSize(3)*real(grid3Offset,pReal)/real(grid(3),pReal) 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 = 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 !-------------------------------------------------------------------------------------------------- ! 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 !-------------------------------------------------------------------------------------------------- !> @brief Parses geometry file !> @details important variables have an implicit "save" attribute. Therefore, this function is ! supposed to be called only once! !-------------------------------------------------------------------------------------------------- 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 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, & myStat, & l, & !< line counter c, & !< counter for # microstructures in line o, & !< order of "to" packing e, & !< "element", i.e. spectral collocation point i, j grid = -1 geomSize = -1.0_pReal !-------------------------------------------------------------------------------------------------- ! 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) call IO_error(100,ext_msg=trim(geometryFile)) allocate(character(len=fileLength)::rawData) read(fileUnit) rawData close(fileUnit) !-------------------------------------------------------------------------------------------------- ! get header length endPos = index(rawData,new_line('')) if(endPos <= index(rawData,'head')) then startPos = len(rawData) call IO_error(error_ID=841, ext_msg='readGeom') else chunkPos = IO_stringPos(rawData(1:endPos)) 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 if (endPos < startPos) endPos = len(rawData) ! end of file without new line line = rawData(startPos:endPos) 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,.true.)) ) case ('grid') if (chunkPos(1) > 6) then do j = 2,6,2 select case (IO_lc(IO_stringValue(line,chunkPos,j))) case('a') grid(1) = IO_intValue(line,chunkPos,j+1) case('b') grid(2) = IO_intValue(line,chunkPos,j+1) case('c') grid(3) = IO_intValue(line,chunkPos,j+1) end select enddo endif case ('size') if (chunkPos(1) > 6) then do j = 2,6,2 select case (IO_lc(IO_stringValue(line,chunkPos,j))) case('x') geomSize(1) = IO_floatValue(line,chunkPos,j+1) case('y') geomSize(2) = IO_floatValue(line,chunkPos,j+1) case('z') geomSize(3) = IO_floatValue(line,chunkPos,j+1) end select enddo endif case ('homogenization') if (chunkPos(1) > 1) h = IO_intValue(line,chunkPos,2) end select enddo !-------------------------------------------------------------------------------------------------- ! sanity checks 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, 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 do while (startPos < len(rawData)) 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 l = l + 1 chunkPos = IO_stringPos(trim(line)) noCompression: if (chunkPos(1) /= 3) then c = chunkPos(1) 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) = [(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 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) = [(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, el=e) end subroutine readGeom !--------------------------------------------------------------------------------------------------- !> @brief Calculate position of IPs/cell centres (pretend to be an element) !--------------------------------------------------------------------------------------------------- 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,1,product(grid)) :: ipCoordinates 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 = 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 nodes !-------------------------------------------------------------------------------------------------- !> @brief Calculate IP interface areas !-------------------------------------------------------------------------------------------------- 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(6,1,product(grid)) :: cellEdgeArea 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 cellEdgeArea !-------------------------------------------------------------------------------------------------- !> @brief Calculate IP interface areas normals !-------------------------------------------------------------------------------------------------- pure function cellEdgeNormal(nElems) integer, intent(in) :: nElems real, dimension(3,6,1,nElems) :: cellEdgeNormal 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 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 !-------------------------------------------------------------------------------------------------- !> @brief builds mesh of (distorted) cubes for given coordinates (= center of the cubes) !-------------------------------------------------------------------------------------------------- function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes) real(pReal), intent(in), dimension(:,:,:,:) :: & centres real(pReal), dimension(3,size(centres,2)+1,size(centres,3)+1,size(centres,4)+1) :: & nodes real(pReal), intent(in), dimension(3) :: & gDim real(pReal), intent(in), dimension(3,3) :: & Favg real(pReal), dimension(3,size(centres,2)+2,size(centres,3)+2,size(centres,4)+2) :: & wrappedCentres integer(pInt) :: & i,j,k,n integer(pInt), dimension(3), parameter :: & diag = 1_pInt integer(pInt), dimension(3) :: & shift = 0_pInt, & lookup = 0_pInt, & me = 0_pInt, & iRes = 0_pInt integer(pInt), dimension(3,8) :: & neighbor = reshape([ & 0_pInt, 0_pInt, 0_pInt, & 1_pInt, 0_pInt, 0_pInt, & 1_pInt, 1_pInt, 0_pInt, & 0_pInt, 1_pInt, 0_pInt, & 0_pInt, 0_pInt, 1_pInt, & 1_pInt, 0_pInt, 1_pInt, & 1_pInt, 1_pInt, 1_pInt, & 0_pInt, 1_pInt, 1_pInt ], [3,8]) !-------------------------------------------------------------------------------------------------- ! initializing variables iRes = [size(centres,2),size(centres,3),size(centres,4)] nodes = 0.0_pReal wrappedCentres = 0.0_pReal !-------------------------------------------------------------------------------------------------- ! building wrappedCentres = centroids + ghosts wrappedCentres(1:3,2_pInt:iRes(1)+1_pInt,2_pInt:iRes(2)+1_pInt,2_pInt:iRes(3)+1_pInt) = centres do k = 0_pInt,iRes(3)+1_pInt do j = 0_pInt,iRes(2)+1_pInt do i = 0_pInt,iRes(1)+1_pInt if (k==0_pInt .or. k==iRes(3)+1_pInt .or. & ! z skin j==0_pInt .or. j==iRes(2)+1_pInt .or. & ! y skin i==0_pInt .or. i==iRes(1)+1_pInt ) then ! x skin me = [i,j,k] ! me on skin shift = sign(abs(iRes+diag-2_pInt*me)/(iRes+diag),iRes+diag-2_pInt*me) lookup = me-diag+shift*iRes wrappedCentres(1:3,i+1_pInt, j+1_pInt, k+1_pInt) = & centres(1:3,lookup(1)+1_pInt,lookup(2)+1_pInt,lookup(3)+1_pInt) & - matmul(Favg, real(shift,pReal)*gDim) endif enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! averaging do k = 0_pInt,iRes(3); do j = 0_pInt,iRes(2); do i = 0_pInt,iRes(1) do n = 1_pInt,8_pInt nodes(1:3,i+1_pInt,j+1_pInt,k+1_pInt) = & nodes(1:3,i+1_pInt,j+1_pInt,k+1_pInt) + wrappedCentres(1:3,i+1_pInt+neighbor(1,n), & j+1_pInt+neighbor(2,n), & k+1_pInt+neighbor(3,n) ) enddo enddo; enddo; enddo nodes = nodes/8.0_pReal end function mesh_nodesAroundCentres end module mesh