From e80e055c756c1c3fa3f2d3614141fa8aff6be88f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 31 Oct 2010 16:06:10 +0000 Subject: [PATCH] added three small programs useful for spectral method (fortran sourececode) voronoi fast.f90 and voronoi small memory.f90 are two variants to do a voronoi tessellation and write the result to a mesh file that can be interpreted by mpie_spectral.f90. Difference is the memory management resulting in one faster and one memory saving version. voronoi fast.f90 has also the ability to write out a file for Ricardo Lebensohns spectral code. colormap.f90 is a simple code that can be used to generate colormaps for gmsh. --- processing/post/colormap.f90 | 69 ++++++++ processing/pre/voronoi fast.f90 | 218 ++++++++++++++++++++++++ processing/pre/voronoi small memory.f90 | 218 ++++++++++++++++++++++++ 3 files changed, 505 insertions(+) create mode 100644 processing/post/colormap.f90 create mode 100644 processing/pre/voronoi fast.f90 create mode 100644 processing/pre/voronoi small memory.f90 diff --git a/processing/post/colormap.f90 b/processing/post/colormap.f90 new file mode 100644 index 000000000..17a9ebc38 --- /dev/null +++ b/processing/post/colormap.f90 @@ -0,0 +1,69 @@ +program colormap + +implicit none + real startH, endH, startS, endS, startL, endL, h_strich,x,c,m + integer steps,i,j + character(len=100) name + print*, '******************************************************************************' + print*, ' write colormap for gmsh' + print*, '******************************************************************************' + print*, '' + write(*, '(A)', advance = 'NO') 'Please enter startvalue for L: ' + read(*, *), startL + write(*, '(A)', advance = 'NO') 'Please enter endvalue for L: ' + read(*, *), endL + write(*, '(A)', advance = 'NO') 'Please enter startvalue for S: ' + read(*, *), startS + write(*, '(A)', advance = 'NO') 'Please enter endvalue for S: ' + read(*, *), endS + write(*, '(A)', advance = 'NO') 'Please enter steps: ' + read(*, *), steps + do j=0,360 + write(name, *) j + name=adjustl(name) + startH = real(j) + endH =real(j) + open(20, file = ('colormap_')//trim(name)//('.map')) + write(20,*),'View.ColorTable = {' + if(endH6.0) h_strich = h_strich-6.0 + c = (1- abs(2*(startL + real(i)*(endL-startL)/real(steps))-1))*sqrt(startS + real(i)*(endS-startS)/real(steps)) + x = c*(1- abs(mod(h_strich, real(2))-1)) + m = (startL + real(i)*(endL-startL)/real(steps)) -.5*c + if ((0.0 <= h_strich).and.(h_strich<1.0)) then + write(20,*),'{',(c+m)*255,',',(x+m)*255,',',(0.0+m)*255,'},' + else if ((1.0 <= h_strich).and.(h_strich<2.0)) then + write(20,*),'{',(x+m)*255,',',(c+m)*255,',',(0.0+m)*255,'},' + else if ((2.0 <= h_strich).and.(h_strich<3.0)) then + write(20,*),'{',(0.0+m)*255,',',(c+m)*255,',',(x+m)*255,'},' + else if ((3.0 <= h_strich).and.(h_strich<4.0)) then + write(20,*),'{',(0.0+m)*255,',',(x+m)*255,',',(c+m)*255,'},' + else if ((4.0 <= h_strich).and.(h_strich<5.0)) then + write(20,*),'{',(x+m)*255,',',(0.0+m)*255,',',(c+m)*255,'},' + else if ((5.0 <= h_strich).and.(h_strich<=6.0)) then + write(20,*),'{',(c+m)*255,',',(0.0+m)*255,',',(x+m)*255,'},' + endif + enddo + H_strich = (startH + real(i)*(endH-startH)/real(steps))/60 + if(h_strich>6.0) h_strich = h_strich-6.0 + c = (1- abs(2*(startL + real(i)*(endL-startL)/real(steps))-1))*(startS + real(i)*(endS-startS)/real(steps)) + x = c*(1- abs(mod(h_strich, real(2))-1)) + m = (startL + real(i)*(endL-startL)/real(steps)) -.5*c + if ((0.0 <= h_strich).and.(h_strich<1.0)) then + write(20,*),'{',(c+m)*255,',',(x+m)*255,',',(0.0+m)*255,'}' + else if ((1.0 <= h_strich).and.(h_strich<2.0)) then + write(20,*),'{',(x+m)*255,',',(c+m)*255,',',(0.0+m)*255,'}' + else if ((2.0 <= h_strich).and.(h_strich<3.0)) then + write(20,*),'{',(0.0+m)*255,',',(c+m)*255,',',(x+m)*255,'}' + else if ((3.0 <= h_strich).and.(h_strich<4.0)) then + write(20,*),'{',(0.0+m)*255,',',(x+m)*255,',',c+m,'}' + else if ((4.0 <= h_strich).and.(h_strich<5.0)) then + write(20,*),'{',(x+m)*255,',',(0.0+m)*255,',',(c+m)*255,'}' + else if ((5.0 <= h_strich).and.(h_strich<=6.0)) then + write(20,*),'{',(c+m)*255,',',(0.0+m)*255,',',(x+m)*255,'}' + endif + write(20,*),'};' + enddo +end program \ No newline at end of file diff --git a/processing/pre/voronoi fast.f90 b/processing/pre/voronoi fast.f90 new file mode 100644 index 000000000..699dccbe7 --- /dev/null +++ b/processing/pre/voronoi fast.f90 @@ -0,0 +1,218 @@ +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*, ' * mesh file "_GIVEN_NAME_.mesh": 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 mesh 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 mesh file, should be done before the following change of variables + open(20, file = ((trim(name))//'.mesh')) + 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*, 'mesh 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 mesh 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 mesh 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 small memory.f90 b/processing/pre/voronoi small memory.f90 new file mode 100644 index 000000000..ccb444670 --- /dev/null +++ b/processing/pre/voronoi small memory.f90 @@ -0,0 +1,218 @@ +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*, ' * mesh file "_GIVEN_NAME_.mesh": 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 mesh 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 mesh file, should be done before the following change of variables + open(20, file = ((trim(name))//'.mesh')) + 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*, 'mesh 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 mesh 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 mesh 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