diff --git a/processing/pre/voronoi_randomSeeding.f90 b/processing/pre/voronoi_randomSeeding.f90 index 11b97f9d4..4a0e94a29 100644 --- a/processing/pre/voronoi_randomSeeding.f90 +++ b/processing/pre/voronoi_randomSeeding.f90 @@ -18,18 +18,23 @@ program voronoi logical, dimension(:), allocatable :: seedmap character(len=1024) filename - integer(pInt) a, b, c, N_Seeds, seedpoint, i + integer(pInt), dimension(3) :: seedcoord + integer(pInt), dimension(:), allocatable :: rndInit + integer(pInt) a, b, c, N_Seeds, seedpoint, i, randomSeed, rndSize real(pReal), dimension(:,:), allocatable :: grainEuler, seeds real(pReal), parameter :: pi = 3.14159265358979323846264338327950288419716939937510_pReal - real(pReal) randomSeed print*, '******************************************************************************' print*, ' Voronoi description file' print*, '******************************************************************************' print*, '' print*, 'generates:' - print*, ' * description file "_GIVEN_NAME_.seeds":' + print*, ' * description file "_OUTPUT_.seeds":' print*, '' + write(*, '(A)', advance = 'NO') 'Please enter name of output seed file: ' + read(*, *), filename + write(*, '(A)', advance = 'NO') 'Please enter value random seed: ' + read(*, *), randomSeed; randomSeed = max(0_pInt,randomSeed) write(*, '(A)', advance = 'NO') 'Please enter value for first resolution: ' read(*, *), a write(*, '(A)', advance = 'NO') 'Please enter value for second resolution: ' @@ -38,43 +43,50 @@ program voronoi read(*, *), c write(*, '(A)', advance = 'NO') 'Please enter No. of Grains: ' read(*, *), N_Seeds - write(*, '(A)', advance = 'NO') 'Please enter name of geometry file: ' - read(*, *), filename allocate (seedmap(a*b*c)); seedmap = .false. ! logical to store information which position is occupied by a voronoi seed allocate (seeds(N_Seeds,3)) allocate (grainEuler(N_Seeds,3)) + + call random_seed(size=rndSize) + allocate(rndInit(rndSize)) + rndInit = randomSeed + call random_seed(put=rndInit) + call random_seed(get=rndInit) do i=1, N_Seeds call random_number(grainEuler(i,1)) call random_number(grainEuler(i,2)) call random_number(grainEuler(i,3)) - grainEuler(i,1) = (grainEuler(i,1))*360 - grainEuler(i,2) = acos(2.0_pReal*(grainEuler(i,2))-1.0_pReal)*180/pi - grainEuler(i,3) = grainEuler(i,3)*360 + grainEuler(i,1) = (grainEuler(i,1))*360.0 + grainEuler(i,2) = acos(2.0_pReal*(grainEuler(i,2))-1.0_pReal)*180.0/pi + grainEuler(i,3) = grainEuler(i,3)*360.0 enddo !generate random position of seeds for voronoi tessellation - i = 0 - do while (i /= N_Seeds) - call random_number(randomSeed) - seedpoint = int(randomSeed*(a*b*c)) - if (.not.seedmap(seedpoint+1)) then + i = 1 + do while (i <= N_Seeds) + call random_number(seeds(i,1)); seedcoord(1) = min(a,int(seeds(i,1)*a)+1_pInt)-1_pInt + call random_number(seeds(i,2)); seedcoord(2) = min(b,int(seeds(i,2)*b)+1_pInt)-1_pInt + call random_number(seeds(i,3)); seedcoord(3) = min(c,int(seeds(i,3)*c)+1_pInt)-1_pInt + seedpoint = seedcoord(1) + seedcoord(2)*a + seedcoord(3)*a*b + if (.not. seedmap(seedpoint+1)) then seedmap(seedpoint+1) = .true. i = i + 1 - seeds(i,1) = real(mod((seedpoint), a)+1)/real(a, pReal) - seeds(i,2) = real(mod(((seedpoint)/a), b)+1)/real(b,pReal) - seeds(i,3) = real(mod(((seedpoint)/(a*b)), c)+1)/real(c,pReal) end if end do ! write description file with orientation and position of each seed open(21, file = trim(filename)//('.seeds')) - write(21, '(A, I2, A, I2, A, I2)'), 'resolution a ', a, ' b ', b, ' c ', c - write(21,*), 'grains', N_Seeds + write(21, '(i1,a1,a6)') 4,achar(9),'header' + write(21, '(A, I2, A, I2, A, I2)') 'resolution a ', a, ' b ', b, ' c ', c + write(21,*) 'grains', N_Seeds + write(21,*) 'random seed ',rndInit(1) + write(21,'(6(a,a1))') 'x',achar(9),'y',achar(9),'z',achar(9),'phi1',achar(9),'Phi',achar(9),'phi2',achar(9) do i = 1, n_Seeds - write(21, '(6(F10.6,tr2))'),seeds(i,1), seeds(i,2), seeds(i,3),& - grainEuler(i,1), grainEuler(i,2), grainEuler(i,3) + write(21, '(6(F10.6,a1))'),seeds(i,1), achar(9), seeds(i,2), achar(9), seeds(i,3), achar(9), & + grainEuler(i,1),achar(9), grainEuler(i,2),achar(9), grainEuler(i,3),achar(9) end do close(21) + deallocate (rndInit) end program voronoi diff --git a/processing/pre/voronoi_tessellation.f90 b/processing/pre/voronoi_tessellation.f90 index 46370169a..68662b0d7 100644 --- a/processing/pre/voronoi_tessellation.f90 +++ b/processing/pre/voronoi_tessellation.f90 @@ -170,28 +170,28 @@ program voronoi implicit none logical gotN_Seeds, gotResolution - character(len=1024) input_name, output_name, format1, format2, N_Digits, line - integer(pInt) a, b, c, N_Seeds, seedPoint, minDistance, myDistance, i, j, k, l, m - integer(pInt), dimension(:), allocatable :: grainMap - integer(pInt) coordinates(3) - integer(pInt), dimension (15) :: posGeom + logical, dimension(:), allocatable :: grainCheck + character(len=1024) input_name, output_name, format1, format2, N_Digits, line, key + integer(pInt) a, b, c, N_Seeds, seedPoint, theGrain, minDistance, myDistance, i, j, k, l, m + integer(pInt), dimension(3) :: coordinates + integer(pInt), dimension (1+2*7) :: posGeom real(pReal), dimension(:,:), allocatable :: grainEuler, seeds + real(pReal), dimension(3) :: dim real(pReal), parameter :: pi = 3.14159265358979323846264338327950288419716939937510_pReal - real(pReal) scaling print*, '******************************************************************************' print*, ' Spectral Method Problem Set-up' print*, '******************************************************************************' print*, '' print*, 'generates:' - print*, ' * geom file "_GIVEN_NAME_.geom": Geometrical information for solver' + print*, ' * geom file "_OUTPUT_.geom": Geometrical information for solver' print*, ' * material file "material.config": Orientation information for solver' - print*, ' * "_GIVEN_NAME_.spectral": combined information for solver' + print*, ' * "_OUTPUT_.spectral": combined information for solver' print*, '' - write(*, '(A)', advance = 'NO') 'Enter filename of input file (extension .seeds): ' - read(*, *), input_name - write(*, '(A)', advance = 'NO') 'Enter filename of output file: ' + write(*, '(A)', advance = 'NO') 'Enter output filename: ' read(*, *), output_name + write(*, '(A)', advance = 'NO') 'Enter seed input file (extension .seeds): ' + read(*, *), input_name open(20, file = trim(input_name)//('.seeds'), status='old', action='read') rewind(20) @@ -221,19 +221,37 @@ program voronoi 100 allocate(grainEuler(N_Seeds,3)) allocate(Seeds(N_Seeds,3)) + allocate(grainCheck(N_Seeds)) + grainCheck = .false. print*, 'resolution: ' ,a,b,c - write(*, '(A)', advance = 'NO') 'Enter scaling factor: ' - read(*, *), scaling + write(*, '(A)', advance = 'NO') 'New first resolution: ' + read(*, *), a + write(*, '(A)', advance = 'NO') 'New second resolution: ' + read(*, *), b + write(*, '(A)', advance = 'NO') 'New third resolution: ' + read(*, *), c - a = int(a*scaling) - b = int(b*scaling) - c = int(c*scaling) - + write(*, '(A)', advance = 'NO') 'First dimension: ' + read(*, *), dim(1) + write(*, '(A)', advance = 'NO') 'Second dimension: ' + read(*, *), dim(2) + write(*, '(A)', advance = 'NO') 'Third dimension: ' + read(*, *), dim(3) + + rewind(20) + read(20,'(a1024)') line + posGeom = IO_stringPos(line,2) + key = IO_stringValue(line,posGeom,2) + if (IO_lc(key(1:4)) == 'head') then + do i=1,IO_intValue(line,posGeom,1); read(20,'(a1024)') line; enddo + else + rewind(20) + endif do i=1, N_seeds read(20,'(a1024)') line if (IO_isBlank(line)) cycle ! skip empty lines - posGeom = IO_stringPos(line,12) + posGeom = IO_stringPos(line,6) Seeds(i,1)=IO_floatValue(line,posGeom,1) Seeds(i,2)=IO_floatValue(line,posGeom,2) Seeds(i,3)=IO_floatValue(line,posGeom,3) @@ -247,7 +265,6 @@ program voronoi seeds(:,2) = seeds(:,2)*real(b, pReal) seeds(:,3) = seeds(:,3)*real(c, pReal) - allocate (grainMap(a*b*c)) ! calculate No. of digits needed for name of the grains i = 1 + int( log10(real( N_Seeds ))) write(N_Digits, *) i @@ -274,49 +291,51 @@ program voronoi end do close(20) print*, '' - print*, 'material config file is written out' + print*, 'material.config done.' !write header of geom file open(20, file = ((trim(output_name))//'.geom')) - write(20, '(A, I2, A, I2, A, I2)'), 'resolution a ', a, ' b ', b, ' c ', c - write(20, '(A, I4, A, I4, A, I4)'), 'dimension x ', a, ' y ', b, ' z ', c + open(21, file = ((trim(output_name))//'.spectral')) + write(20, '(A)'), '3 header' + write(20, '(A, I8, A, I8, A, I8)'), 'resolution a ', a, ' b ', b, ' c ', c + write(20, '(A, g15.10, A, g15.10, A, g15.10)'), 'dimension x ', dim(1), ' y ', dim(2), ' z ', dim(3) write(20, '(A)'), 'homogenization 1' -!initialize varibles, change values of some numbers for faster execution - format1 = '(I'//trim(N_Digits)//'.'//trim(N_Digits)//')' + format1 = '(I'//trim(N_Digits)//'.'//trim(N_Digits)//')' ! geom format + format2 = '(3(tr2, f6.2), 3(I10), I10, a)' ! spectral (Lebensohn) format -! perform voronoi tessellation and write result to file and to grainMap - do i = 1, a*b*c +! perform voronoi tessellation and write result to files + do i = 0, a*b*c-1 minDistance = a*a+b*b+c*c do j = 1, N_Seeds - do k = -1, 1 - do l = -1, 1 - do m = -1, 1 - myDistance = ((mod((i-1), a) +1-seeds(j,1)+m*a)**2+& - (mod(((i-1)/a), b) +1-seeds(j,2)+l*b)**2+& - (mod(((i-1)/(a*b)), c) +1-seeds(j,3)+k*c)**2) + do k = -1, 1 ! left, me, right image + do l = -1, 1 ! front, me, back image + do m = -1, 1 ! lower, me, upper image + myDistance = (( mod(i, a)+1 -seeds(j,1)-m*a)**2 + & + (mod((i/a), b)+1 -seeds(j,2)-l*b)**2 + & + (mod((i/(a*b)), c)+1 -seeds(j,3)-k*c)**2) if (myDistance < minDistance) then minDistance = myDistance - grainMap(i) = j + theGrain = j end if end do end do end do end do - write(20, format1), grainMap(i) + grainCheck(theGrain) = .true. + write(20, trim(format1)), theGrain + write(21, trim(format2)), grainEuler(theGrain,1), grainEuler(theGrain,2), grainEuler(theGrain,3), & + mod(i, a)+1, mod((i/a), b)+1, mod((i/(a*b)), c)+1, & + theGrain, ' 1' end do close(20) - print*, 'voronoi tesselation finished' - - open(20, file = ((trim(output_name))//'.spectral')) - format1 = '(3(tr2, f6.2), 3(I10), I10, a)' - do i = 1, a*b*c - j = grainMap(i) - write(20, trim(format1)), grainEuler(j,1), grainEuler(j,2), grainEuler(j,3), & - &mod((i-1), a)+1, mod(((i-1)/a), b)+1, mod(((i-1)/(a*b)), c)+1, & - &j, ' 1' - end do - print*, 'geometry files are written out' + close(21) + print*, 'voronoi tesselation done.' + if (all(grainCheck)) then + print*, 'all grains mapped!' + else + print*, 'only',count(grainCheck),'grains mapped!' + endif end program voronoi