diff --git a/processing/pre/voronoi_fast.f90 b/processing/pre/voronoi_fast.f90 deleted file mode 100644 index e86dc9bd9..000000000 --- a/processing/pre/voronoi_fast.f90 +++ /dev/null @@ -1,237 +0,0 @@ -program voronoi - use prec, only: pReal, pInt - - implicit none - - logical, dimension(:), allocatable :: seedmap - character(len=1024) name, format, format2, N_Digits - character choice - integer(pInt) a, b, c, ab, abc, abc_Red, N_Seeds, seedPoint, minDistance, myDistance, i, j, k, l, m - integer(pInt), dimension(:), allocatable :: seeds, grainMap, visual_Case - integer(pInt) coordinates(3) - real(pReal), dimension(:), allocatable :: grainEuler - real(pReal), parameter :: pi = 3.14159265358979323846264338327950288419716939937510_pReal - real(pReal) randomSeed - - print*, '******************************************************************************' - print*, ' Spectral Method Problem Set-up' - print*, '******************************************************************************' - print*, '' - print*, 'generates:' - print*, ' * geom file "_GIVEN_NAME_.geom": Geometrical information for solver' - print*, ' * material file "material.config": Orientation information for solver' - print*, ' * "_GIVEN_NAME_.spectral": combined information for solver' - print*, 'optional output:' - print*, ' * view file "_GIVEN_NAME_2D.msh": Information for visualization in gmsh' - print*, ' * view file "_GIVEN_NAME_3D.msh": Information for visualization in gmsh' - print*, '' - print*, 'hints:' - print*, ' * a+b+c should not exeed 30' - print*, ' * file extension is added to given name' - print*, '' - write(*, '(A)', advance = 'NO') 'Please enter value for a: ' - read(*, *), a - write(*, '(A)', advance = 'NO') 'Please enter value for b: ' - read(*, *), b - write(*, '(A)', advance = 'NO') 'Please enter value for c: ' - 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(*, *), name - write(*, '(A)', advance = 'NO') 'Should the visualization files be generated (y/n)? ' - read(*, *), choice - -! calculate No. of digits needed for name of the grains - i = 1 + int( log10(real( N_Seeds ))) - write(N_Digits, *) i - N_Digits = adjustl( N_Digits ) - allocate(grainEuler(3*N_Seeds)) - -!write material.config header and add a microstructure entry for every grain - open(20, file = trim(name)//('_material.config')) - write(20, '(A)'), '' - format = '(A, I'//trim(N_Digits)//'.'//trim(N_Digits)//', A)' - format2 = '(A, I'//trim(N_Digits)//', A)' - do i = 1, N_Seeds - write(20, trim(format)), '[Grain', i, ']' - write(20, '(A)'), 'crystallite 1' - write(20, trim(format2)), '(constituent) phase 1 texture ', i, ' fraction 1.0' - end do - -! get random euler angles for every grain, store them in grainEuler and write them to the material.config file - format2 = '(A, F10.6, A, F10.6, A, F10.6, A)' - write(20, '(/, A)'), '' - do i = 1, N_Seeds - call random_number(grainEuler(i*3 -2)) - call random_number(grainEuler(i*3 -1)) - call random_number(grainEuler(i*3 -0)) - grainEuler(i*3 -2) = (grainEuler(i*3 -2))*360 - grainEuler(i*3 -1) = acos(2.0_pReal*(grainEuler(i*3 -1))-1.0_pReal)*180/pi - grainEuler(i*3 -0) = grainEuler(i*3)*360 - write(20, trim(format)), '[Grain', i, ']' - write(20, trim(format2)), '(gauss) phi1 ', grainEuler(i*3-2), ' Phi ', grainEuler(i*3-1), & - &' Phi2 ', grainEuler(i*3), ' scatter 0 fraction 1' - end do - close(20) - print*, '' - print*, 'material config file is written out' - -!write header of geom file, should be done before the following change of variables - open(20, file = ((trim(name))//'.geom')) - write(20, '(A, I2, A, I2, A, I2)'), 'resolution a ', 2**a, ' b ', 2**b, ' c ', 2**c - write(20, '(A, I4, A, I4, A, I4)'), 'dimension x ', 2**a, ' y ', 2**b, ' z ', 2**c - write(20, '(A)'), 'homogenization 1' - -!initialize varibles, change values of some numbers for faster execution - a = 2**a - b = 2**b - c = 2**c - ab = a * b - abc = a * b * c - abc_Red = abc -(a-1)*(b-1)*(c-1) - format = '(I'//trim(N_Digits)//'.'//trim(N_Digits)//')' - - allocate (seedmap(abc)); seedmap = .false. - allocate (seeds(3*N_Seeds)) - allocate (grainMap(abc)) - allocate (visual_Case(abc_Red)) - k = 1 - -!build array with x-y-z-coordinates of each point - do i = 1, abc - coordinates(1) = mod((i-1), a) +1 - coordinates(2) = mod(((i-1)/a), b) +1 - coordinates(3) = mod(((i-1)/(ab)), c) +1 - if((coordinates(3) == 1)) then - visual_Case(k) = i - k = k +1 - else - if((coordinates(2) == 1)) then - visual_Case(k) = i - k = k +1 - else - if((coordinates(1) == 1)) then - visual_Case(k) = i - k = k +1 - else - end if - end if - end if - end do - -!generate random position of seeds for voronoi tessellation - i = 0 - do while (i /= N_Seeds) - call random_number(randomSeed) - seedpoint = int(randomSeed*(abc)) - if (.not.seedmap(seedpoint+1)) then - seedmap(seedpoint+1) = .true. - seeds(i*3+1) = (mod((seedpoint), a)+1) - seeds(i*3+2) = (mod(((seedpoint)/a), b)+1) - seeds(i*3+3) = (mod(((seedpoint)/(ab)), c)+1) - i = i +1 - else - end if - end do - -! perform voronoi tessellation and write result to file and to grainMap - do i = 1, abc - 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*3-2)+m*a)**2+& - (mod(((i-1)/a), b) +1-seeds(j*3-1)+l*b)**2+& - (mod(((i-1)/(ab)), c) +1-seeds(j*3-0)+k*c)**2) - if (myDistance < minDistance) then - minDistance = myDistance - grainMap(i) = j - end if - end do - end do - end do - end do - write(20, format), grainMap(i) - end do - close(20) - print*, 'voronoi tesselation finished' - - open(20, file = ((trim(name))//'.spectral')) - format = '(tr2, f6.2, tr2, f6.2, tr2, f6.2, I10, I10, I10, I10, a)' - do i = 1, abc - j = grainMap(i) - write(20, trim(format)), grainEuler(j*3-2), grainEuler(j*3-1), grainEuler(j*3), & - &mod((i-1), a)+1, mod(((i-1)/a), b)+1, mod(((i-1)/(ab)), c)+1, & - &j, ' 1' - end do - print*, 'geometry files are written out' - -!write visualization files (in case wanted) - if (choice == 'y' .or. choice == 'Y') then - print*, 'for more information on gmsh: http://geuz.org/gmsh/' - -! write full geometry out - open(20, file = ((trim(name))//'_3Dfull.msh')) - write(20, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', abc - do i = 1, abc - write(20, '(I10, I10, I10, I10)'), i, mod((i-1), a) +1, mod(((i-1)/a), b) +1, mod(((i-1)/(ab)), c) +1 - end do - write(20, '(A, /, A, /, I10)'), '$EndNodes', '$Elements', abc - do i = 1, abc - write(20, '(I10, A, I10, A, I10)'), i, ' 15 2', grainMap(i), ' 2', i - end do - write(20, '(A)'), '$EndElements' - write(20, '(A, /, A, /, A, /, A, /, A, /, A, /, A, /, A, /, I10)'), '$NodeData', '1', '"Grain No."', '1', & - &'0.0', '3', '0', '1', abc - do i = 1, abc - write(20, '(I10, tr2, I10)'), i, grainMap(i) - end do - write(20, *), '$EndNodeData' - close(20) - -! write 3d skin out - open(20, file = ((trim(name))//'_3D.msh')) - write(20, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', abc_Red - do j = 1, abc_Red - i = visual_Case(j) - write(20, '(I10, I10, I10, I10)'), i, mod((i-1), a) +1, mod(((i-1)/a), b) +1, mod(((i-1)/(ab)), c) +1 - end do - write(20, '(A, /, A, /, I10)'), '$EndNodes', '$Elements', abc_Red - do j = 1, abc_Red - i = visual_case(j) - write(20, '(I10, A, I10, A, I10)'), i, ' 15 2', grainMap(i), ' 2', i - end do - write(20, '(A)'), '$EndElements' - write(20, '(A, /, A, /, A, /, A, /, A, /, A, /, A, /, A, /, I10)'), '$NodeData', '1', '"Grain No."', '1', & - &'0.0', '3', '0', '1', abc_Red - do j = 1, abc_Red - i = visual_case(j) - write(20, '(I10, tr2, I10)'), i, grainMap(i) - end do - write(20, *), '$EndNodeData' - close(20) - -! write 2d geometry out - open(20, file = ((trim(name))//'_2D.msh')) - write(20, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', ab - do i = 1, ab - write(20, '(I10, I10, I10, I10)'), i, mod((i-1), a) +1, mod(((i-1)/a), b) +1, mod(((i-1)/(ab)), c) +1 - end do - write(20, '(A, /, A, /, I10)'), '$EndNodes', '$Elements', ab - do j = 1, ab - write(20, '(I10, A, I10, A, I10)'), j, ' 15 2', grainMap(j), ' 2', j - end do - write(20, '(A)'), '$EndElements' - write(20, '(A, /, A, /, A, /, A, /, A, /, A, /, A, /, A, /, I10)'), '$NodeData', '1', '"Grain No."', & - &'1', '0.0', '3', '0', '1', ab - do j = 1, ab - write(20, '(I10, tr2, I10)'), j, grainMap(j) - end do - write(20, *), '$EndNodeData' - close(20) - print*, 'visualization files are written out' - else - end if -end program voronoi diff --git a/processing/pre/voronoi_randomSeeding.f90 b/processing/pre/voronoi_randomSeeding.f90 new file mode 100644 index 000000000..11b97f9d4 --- /dev/null +++ b/processing/pre/voronoi_randomSeeding.f90 @@ -0,0 +1,80 @@ +!prec.f90 407 2009-08-31 15:09:15Z MPIE\f.roters +!############################################################## + MODULE prec +!############################################################## + + implicit none + +! *** Precision of real and integer variables *** + integer, parameter :: pReal = selected_real_kind(15,300) ! 15 significant digits, up to 1e+-300 + integer, parameter :: pInt = selected_int_kind(9) ! up to +- 1e9 + integer, parameter :: pLongInt = 8 ! should be 64bit + + END MODULE prec + +program voronoi + use prec, only: pReal, pInt + implicit none + + logical, dimension(:), allocatable :: seedmap + character(len=1024) filename + integer(pInt) a, b, c, N_Seeds, seedpoint, i + 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*, '' + write(*, '(A)', advance = 'NO') 'Please enter value for first resolution: ' + read(*, *), a + write(*, '(A)', advance = 'NO') 'Please enter value for second resolution: ' + read(*, *), b + write(*, '(A)', advance = 'NO') 'Please enter value for third resolution: ' + 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)) + + 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 + 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 + 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 + 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) + end do + close(21) +end program voronoi diff --git a/processing/pre/voronoi_small_memory.f90 b/processing/pre/voronoi_small_memory.f90 deleted file mode 100644 index 24cd6eab7..000000000 --- a/processing/pre/voronoi_small_memory.f90 +++ /dev/null @@ -1,218 +0,0 @@ -program voronoi - use prec, only: pReal, pInt - - implicit none - - logical seed - character(len=1024) name, format, format2, N_Digits - character choice - integer(pInt) a, b, c, ab, abc, abc_Red, N_Seeds, seedPoint, minDistance, myDistance, i, j, k, l, m, n, z - integer(pInt), dimension(:), allocatable :: seedmap, grainMap, visual_Case - integer(pInt) coordinates(3) - real(pReal) grainEuler(3) - real(pReal), parameter :: pi = 3.14159265358979323846264338327950288419716939937510_pReal - real(pReal) randomSeed - - print*, '******************************************************************************' - print*, ' Spectral Method Problem Set-up' - print*, '******************************************************************************' - print*, '' - print*, 'generates:' - print*, ' * geometry file "_GIVEN_NAME_.geom": Geometrical information for solver' - print*, ' * material file "material.config": Orientation information for solver' - print*, 'optional output:' - print*, ' * view file "_GIVEN_NAME_2D.msh": Information for visualization in gmsh' - print*, ' * view file "_GIVEN_NAME_3D.msh": Information for visualization in gmsh' - print*, '' - print*, 'hints:' - print*, ' * a+b+c should not exeed 30' - print*, ' * file extension is added to given name' - print*, '' - write(*, '(A)', advance = 'NO') 'Please enter value for a: ' - read(*, *), a - write(*, '(A)', advance = 'NO') 'Please enter value for b: ' - read(*, *), b - write(*, '(A)', advance = 'NO') 'Please enter value for c: ' - 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(*, *), name - write(*, '(A)', advance = 'NO') 'Should the visualization files be generated (y/n)? ' - read(*, *), choice - -! calculate No. of digits needed for name of the grains - i = 1 + int( log10(real( N_Seeds ))) - write(N_Digits, *) i - N_Digits = adjustl( N_Digits ) - -!write material.config header and add a microstructure entry for every grain - open(20, file = trim(name)//('_material.config')) - write(20, '(A)'), '' - format = '(A, I'//trim(N_Digits)//'.'//trim(N_Digits)//', A)' - format2 = '(A, I'//trim(N_Digits)//', A)' - do i = 1, N_Seeds - write(20, trim(format)), '[Grain', i, ']' - write(20, '(A)'), 'crystallite 1' - write(20, trim(format2)), '(constituent) phase 1 texture ', i, ' fraction 1.0' - end do - -! get random euler angles for every grain, store them in grainEuler and write them to the material.config file - format2 = '(A, F10.6, A, F10.6, A, F10.6, A)' - write(20, '(/, A)'), '' - do i = 1, N_Seeds - call random_number(grainEuler(1)) - call random_number(grainEuler(2)) - call random_number(grainEuler(3)) - grainEuler(1) = (grainEuler(1)) * 360 - grainEuler(2) = acos(2.0_pReal * (grainEuler(2))-1.0_pReal) * 180/pi - grainEuler(3) = grainEuler(3) * 360 - write(20, trim(format)), '[Grain', i, ']' - write(20, trim(format2)), '(gauss) phi1 ', grainEuler(1), ' Phi ', grainEuler(2), & - &' Phi2 ', grainEuler(3), ' scatter 0 fraction 1' - end do - close(20) - print*, '' - print*, 'material config file is written out' - -!write header of geometry file, should be done before the following change of variables - open(20, file = ((trim(name))//'.geom')) - write(20, '(A, I2, A, I2, A, I2)'), 'resolution a ', 2**a, ' b ', 2**b, ' c ', 2**c - write(20, '(A, I4, A, I4, A, I4)'), 'dimension x ', 2**a, ' y ', 2**b, ' z ', 2**c - write(20, '(A)'), 'homogenization 1' - -!initialize varibles, change values of some numbers for faster execution - a = 2**a - b = 2**b - c = 2**c - ab = a * b - abc = a * b * c - abc_Red = abc -(a-1)*(b-1)*(c-1) - format = '(I'//trim(N_Digits)//'.'//trim(N_Digits)//')' - - allocate (seedmap(3 * N_Seeds)); seedmap = 0 - allocate (grainMap(abc_Red)) - allocate (visual_Case(abc_Red)) - k = 1 - -!build array with x-y-z-coordinates of each point - do i = 1, abc - coordinates(1) = mod((i-1), a) +1 - coordinates(2) = mod(((i-1)/a), b) +1 - coordinates(3) = mod(((i-1)/(ab)), c) +1 - if((coordinates(3) == 1)) then - visual_Case(k) = i - k = k +1 - else - if((coordinates(2) == 1)) then - visual_Case(k) = i - k = k +1 - else - if((coordinates(1) == 1)) then - visual_Case(k) = i - k = k +1 - else - end if - end if - end if - end do - -!generate random position of seeds for voronoi tessellation - i = 0 - do while (i /= N_Seeds) - call random_number(randomSeed) - seedpoint = int(randomSeed*(abc)) - seed = .false. - coordinates(1) = mod((seedpoint), a) +1 - coordinates(2) = mod(((seedpoint)/a), b) +1 - coordinates(3) = mod(((seedpoint)/(ab)), c) +1 - do j = 1, i - if((coordinates(1) == seedmap(j*3 -2)) .and. & - &(coordinates(2) == seedmap(j*3 -1)) .and. & - &(coordinates(3) == seedmap(j*3 -0))) seed = .true. - end do - if (.not.seed) then - seedmap(i*3 +1) = coordinates(1) - seedmap(i*3 +2) = coordinates(2) - seedmap(i*3 +3) = coordinates(3) - i = i +1 - end if - end do - -! perform voronoi tessellation and write result to file and to grainMap - n = 1 - do i = 1, abc - 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-seedmap(j*3-2)+m*a)**2+& - (mod(((i-1)/a), b) +1-seedmap(j*3-1)+l*b)**2+& - (mod(((i-1)/(ab)), c) +1-seedmap(j*3-0)+k*c)**2) - if (myDistance < minDistance) then - minDistance = myDistance - z = j - end if - end do - end do - end do - end do - write(20, format), z - do k = 1, abc_Red - if(visual_Case(k) == i) then - grainMap(n) = z - n = n +1 - end if - end do - end do - close(20) - print*, 'voronoi tesselation finished' - print*, 'geometry file is written out' - -!write visualization files (in case wanted) - if (choice == 'y' .or. choice == 'Y') then - print*, 'for more information on gmsh: http://geuz.org/gmsh/' - -! write full geometry out - open(20, file = ((trim(name))//'_3D.msh')) - write(20, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', abc_Red - do j = 1, abc_Red - i = visual_Case(j) - write(20, '(I10, I10, I10, I10)'), i, mod((i-1), a) +1, mod(((i-1)/a), b) +1, mod(((i-1)/(ab)), c) +1 - end do - write(20, '(A, /, A, /, I10)'), '$EndNodes', '$Elements', abc_Red - do j = 1, abc_Red - write(20, '(I10, A, I10, A, I10)'), visual_Case(j), ' 15 2', grainMap(j), ' 2', visual_Case(j) - end do - write(20, '(A)'), '$EndElements' - write(20, '(A, /, A, /, A, /, A, /, A, /, A, /, A, /, A, /, I10)'), '$NodeData', '1', '"Grain No."', '1', & - &'0.0', '3', '0', '1', abc_Red - do j = 1, abc_Red - write(20, '(I10, tr2, I10)'), visual_Case(j), grainMap(j) - end do - write(20, *), '$EndNodeData' - close(20) - -! write 2d geometry out - open(20, file = ((trim(name))//'_2D.msh')) - write(20, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', ab - do i = 1, ab - write(20, '(I10, I10, I10, I10)'), i, mod((i-1), a) +1, mod(((i-1)/a), b) +1, mod(((i-1)/(ab)), c) +1 - end do - write(20, '(A, /, A, /, I10)'), '$EndNodes', '$Elements', ab - do j = 1, ab - write(20, '(I10, A, I10, A, I10)'), j, ' 15 2', grainMap(j), ' 2', j - end do - write(20, '(A)'), '$EndElements' - write(20, '(A, /, A, /, A, /, A, /, A, /, A, /, A, /, A, /, I10)'), '$NodeData', '1', '"Grain No."', & - &'1', '0.0', '3', '0', '1', ab - do j = 1, ab - write(20, '(I10, tr2, I10)'), j, grainMap(j) - end do - write(20, *), '$EndNodeData' - close(20) - print*, 'visualization files are written out' - else - end if -end program voronoi diff --git a/processing/pre/voronoi_tessellation.f90 b/processing/pre/voronoi_tessellation.f90 new file mode 100644 index 000000000..46370169a --- /dev/null +++ b/processing/pre/voronoi_tessellation.f90 @@ -0,0 +1,322 @@ +!prec.f90 407 2009-08-31 15:09:15Z MPIE\f.roters +!############################################################## + MODULE prec +!############################################################## + + implicit none + +! *** Precision of real and integer variables *** + integer, parameter :: pReal = selected_real_kind(15,300) ! 15 significant digits, up to 1e+-300 + integer, parameter :: pInt = selected_int_kind(9) ! up to +- 1e9 + integer, parameter :: pLongInt = 8 ! should be 64bit + + END MODULE prec + +!IO.f90 693 2010-11-04 18:18:01Z MPIE\c.kords +!############################################################## + MODULE IO +!############################################################## + + CONTAINS +!******************************************************************** +! identifies lines without content +!******************************************************************** + pure function IO_isBlank (line) + + use prec, only: pInt + implicit none + + character(len=*), intent(in) :: line + character(len=*), parameter :: blank = achar(32)//achar(9)//achar(10)//achar(13) ! whitespaces + character(len=*), parameter :: comment = achar(35) ! comment id '#' + integer(pInt) posNonBlank, posComment + logical IO_isBlank + + posNonBlank = verify(line,blank) + posComment = scan(line,comment) + IO_isBlank = posNonBlank == 0 .or. posNonBlank == posComment + + return + + endfunction + +!******************************************************************** +! read string value at pos from line +!******************************************************************** + pure function IO_stringValue (line,positions,pos) + + use prec, only: pReal,pInt + implicit none + + character(len=*), intent(in) :: line + integer(pInt), intent(in) :: positions(*),pos + character(len=1+positions(pos*2+1)-positions(pos*2)) IO_stringValue + + if (positions(1) < pos) then + IO_stringValue = '' + else + IO_stringValue = line(positions(pos*2):positions(pos*2+1)) + endif + return + + endfunction + +!******************************************************************** +! read float value at pos from line +!******************************************************************** + pure function IO_floatValue (line,positions,pos) + + use prec, only: pReal,pInt + implicit none + + character(len=*), intent(in) :: line + integer(pInt), intent(in) :: positions(*),pos + real(pReal) IO_floatValue + + if (positions(1) < pos) then + IO_floatValue = 0.0_pReal + else + read(UNIT=line(positions(pos*2):positions(pos*2+1)),ERR=100,FMT=*) IO_floatValue + endif + return +100 IO_floatValue = huge(1.0_pReal) + return + + endfunction + +!******************************************************************** +! read int value at pos from line +!******************************************************************** + pure function IO_intValue (line,positions,pos) + + use prec, only: pReal,pInt + implicit none + + character(len=*), intent(in) :: line + integer(pInt), intent(in) :: positions(*),pos + integer(pInt) IO_intValue + + if (positions(1) < pos) then + IO_intValue = 0_pInt + else + read(UNIT=line(positions(pos*2):positions(pos*2+1)),ERR=100,FMT=*) IO_intValue + endif + return +100 IO_intValue = huge(1_pInt) + return + + endfunction + +!******************************************************************** +! change character in line to lower case +!******************************************************************** + pure function IO_lc (line) + + use prec, only: pInt + implicit none + + character (len=*), intent(in) :: line + character (len=len(line)) IO_lc + integer(pInt) i + + IO_lc = line + do i=1,len(line) + if(640) + left = right + verify(line(right+1:),sep) + right = left + scan(line(left:),sep) - 2 + if ( IO_stringPos(1)' + format1 = '(A, I'//trim(N_Digits)//'.'//trim(N_Digits)//', A)' + format2 = '(A, I'//trim(N_Digits)//', A)' + do i = 1, N_Seeds + write(20, trim(format1)), '[Grain', i, ']' + write(20, '(A)'), 'crystallite 1' + write(20, trim(format2)), '(constituent) phase 1 texture ', i, ' fraction 1.0' + end do + +! get random euler angles for every grain, store them in grainEuler and write them to the material.config file + format2 = '(6(A, F10.6))' + write(20, '(/, A)'), '' + do i = 1, N_Seeds + write(20, trim(format1)), '[Grain', i, ']' + write(20, trim(format2)), '(gauss) phi1 ', grainEuler(i,1), ' Phi ', grainEuler(i,2), & + &' Phi2 ', grainEuler(i,3), ' scatter 0 fraction 1' + end do + close(20) + print*, '' + print*, 'material config file is written out' + +!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 + write(20, '(A)'), 'homogenization 1' + +!initialize varibles, change values of some numbers for faster execution + format1 = '(I'//trim(N_Digits)//'.'//trim(N_Digits)//')' + + +! perform voronoi tessellation and write result to file and to grainMap + do i = 1, a*b*c + 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) + if (myDistance < minDistance) then + minDistance = myDistance + grainMap(i) = j + end if + end do + end do + end do + end do + write(20, format1), grainMap(i) + 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' + +end program voronoi