From 322cbd2597b6a27be9e495971b93d6d3958ecc3e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 7 Jan 2011 14:37:05 +0000 Subject: [PATCH] changed tools for voronoi tessellation to match new specification (file extension and resolution). Also renamed the files, because the space characters cause trouble improved reconstruct.f90 and spectral_post.py, both files work now under linux added make_reconstruct.py, small shell script for using f2py --- processing/post/make_reconstruct.py | 8 + processing/post/reconstruct.f90 | 260 +++++++++++++++++- processing/post/spectral_post.py | 183 ++++++------ .../{voronoi fast.f90 => voronoi_fast.f90} | 22 +- ...ll memory.f90 => voronoi_small_memory.f90} | 22 +- 5 files changed, 375 insertions(+), 120 deletions(-) create mode 100644 processing/post/make_reconstruct.py rename processing/pre/{voronoi fast.f90 => voronoi_fast.f90} (92%) rename processing/pre/{voronoi small memory.f90 => voronoi_small_memory.f90} (90%) diff --git a/processing/post/make_reconstruct.py b/processing/post/make_reconstruct.py new file mode 100644 index 000000000..3cbafd3f7 --- /dev/null +++ b/processing/post/make_reconstruct.py @@ -0,0 +1,8 @@ +#!/bin/bash + +# This script is used to compile the python module used for geometry reconstruction. +# It uses the fortran wrapper f2py that is included in the numpy package to construct the +# module reconstruct.so out of the fortran code reconstruct.f90 +# written by M. Diehl, m.diehl@mpie.de + +f2py -c -m reconstruct reconstruct.f90 diff --git a/processing/post/reconstruct.f90 b/processing/post/reconstruct.f90 index 15541f61a..753d4c57a 100644 --- a/processing/post/reconstruct.f90 +++ b/processing/post/reconstruct.f90 @@ -1,14 +1,12 @@ ! -*- f90 -*- subroutine simple(defgrad,res_x,res_y,res_z,geomdimension,current_configuration) implicit none -! *** Precision of real and integer variables *** - integer, parameter :: pReal = selected_real_kind(15) ! 15 significant digits, up to 1e+-300 - integer, parameter :: pInt = selected_int_kind(9) ! up to +- 1e9 - integer(pInt) i,j,k - integer(pInt) res_x, res_y, res_z - integer(pInt), dimension(3) :: resolution - real(pReal) , dimension(3) :: geomdimension - real(pReal) , dimension(3,3) :: temp33_Real + integer, parameter :: pDouble = selected_real_kind(15,50) + integer i,j,k + integer res_x, res_y, res_z + integer, dimension(3) :: resolution + real*8, dimension(3) :: geomdimension + real*8, dimension(3,3) :: temp33_Real real*8 current_configuration(res_x, res_y, res_z,3) real*8 defgrad(res_x, res_y, res_z,3,3) !f2py intent(in) res_x, res_y, res_z @@ -19,29 +17,263 @@ subroutine simple(defgrad,res_x,res_y,res_z,geomdimension,current_configuration) !f2py depend(res_x, res_y, res_z) defgrad resolution(1) = res_x; resolution(2) = res_y; resolution(3) = res_z - print*, defgrad + do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) if((k==1).and.(j==1).and.(i==1)) then - temp33_Real =0.0_pReal + temp33_Real =0.0_pDouble else if((j==1).and.(i==1)) then temp33_Real(1,:) = temp33_Real(1,:) + matmul(defgrad(i,j,k,:,:),& - (/0.0_pReal,0.0_pReal,(real(resolution(3), pReal)/geomdimension(3))/)) + (/0.0_pDouble,0.0_pDouble,geomdimension(3)/(real(resolution(3)))/)) temp33_Real(2,:) = temp33_Real(1,:) temp33_Real(3,:) = temp33_Real(1,:) current_configuration(i,j,k,:) = temp33_Real(1,:) else if(i==1) then temp33_Real(2,:) = temp33_Real(2,:) + matmul(defgrad(i,j,k,:,:),& - (/0.0_pReal,(real(resolution(2),pReal)/geomdimension(2)),0.0_pReal/)) + (/0.0_pDouble,geomdimension(2)/(real(resolution(2))),0.0_pDouble/)) temp33_Real(3,:) = temp33_Real(2,:) current_configuration(i,j,k,:) = temp33_Real(2,:) else temp33_Real(3,:) = temp33_Real(3,:) + matmul(defgrad(i,j,k,:,:),& - (/(real(resolution(1), pReal)/geomdimension(1)),0.0_pReal,0.0_pReal/)) + (/geomdimension(1)/(real(resolution(1))),0.0_pDouble,0.0_pDouble/)) current_configuration(i,j,k,:) = temp33_Real(3,:) endif endif - endif + endif + !print*, current_configuration enddo; enddo; enddo end subroutine simple + +subroutine advanced(defgrad,res_x,res_y,res_z,geomdimension,current_configuration) + implicit none + integer, parameter :: pDouble = selected_real_kind(15,50) + integer i,j,k,l + integer a,b,c + integer res_x, res_y, res_z + integer, dimension(3) :: resolution + real*8, dimension(3) :: geomdimension + real*8, dimension(3,3) :: temp33_Real, defgrad_av !matrix, stores 3 times a 3dim position vector + real*8 current_configuration(res_x, res_y, res_z,3) + real*8 defgrad(res_x, res_y, res_z,3,3) +! some declarations for the wrapping to python + !f2py intent(in) res_x, res_y, res_z + !f2py intent(in) geomdimension + !f2py intent(out) current_configuration + !f2py intent(in) defgrad + !f2py depend(res_x, res_y, res_z) current_configuration + !f2py depend(res_x, res_y, res_z) defgrad + do i=1, 3; do j=1,3 ! + defgrad_av(i,j) = sum(defgrad(:,:,:,i,j)) /real(res_x*res_y*res_z) + enddo; enddo + + current_configuration(i,j,k,:) = 0.0_pDouble + do l=1,6 ! do 6 different paths, by using 3 'master edges' + select case(l) + case (1) + a = 1; b = 2; c = 3 + case (2) + a = 1; b = 3; c = 2 + case (3) + a = 2; b = 1; c = 3 + case (4) + a = 2; b = 3; c = 1 + case (5) + a = 3; b = 1; c = 2 + case (6) + a = 3; b = 2; c = 1 + end select + + resolution(a) = res_x; resolution(b) = res_y; resolution(c) = res_z + + do k = 1, resolution(c); do j = 1, resolution(b); do i = 1, resolution(a) + if((k==1).and.(j==1).and.(i==1)) then ! first FP + temp33_Real = 0.0_pDouble ! all positions set to zero + else + if((j==1).and.(i==1)) then ! all FPs on the 'master edge' + temp33_Real(1,:) = temp33_Real(1,:) + matmul((defgrad(i,j,k-1,:,:)+defgrad(i,j,k-1,:,:))/2.0_pDouble,& !using the average defgrad of the current step + (/0.0_pDouble,0.0_pDouble,geomdimension(c)/(real(resolution(c)))/)) + temp33_Real(2,:) = temp33_Real(1,:) + temp33_Real(3,:) = temp33_Real(1,:) + current_configuration(i,j,k,:) = current_configuration(i,j,k,:) + temp33_Real(1,:) + else + if(i==1) then + temp33_Real(2,:) = temp33_Real(2,:) + matmul((defgrad(i,j-1,k,:,:)+defgrad(i,j,k,:,:))/2.0_pDouble,& + (/0.0_pDouble,geomdimension(b)/(real(resolution(b))),0.0_pDouble/)) + temp33_Real(3,:) = temp33_Real(2,:) + current_configuration(i,j,k,:) = current_configuration(i,j,k,:) + temp33_Real(2,:) + else + temp33_Real(3,:) = temp33_Real(3,:) + matmul((defgrad(i-1,j,k,:,:)+defgrad(i,j,k,:,:))/2.0_pDouble,& + (/geomdimension(a)/(real(resolution(a))),0.0_pDouble,0.0_pDouble/)) + current_configuration(i,j,k,:) = current_configuration(i,j,k,:) + temp33_Real(3,:) + endif + endif + endif + enddo; enddo; enddo ! end of one reconstruction +enddo ! end of 6 reconstructions with different pathes +current_configuration = current_configuration/6.0_pDouble +end subroutine advanced + +subroutine mesh(inter,res_x,res_y,res_z,gdim,meshgeom) + implicit none + integer, parameter :: pDouble = selected_real_kind(15,50) + integer i,j,k + integer res_x, res_y, res_z + integer, dimension(3) :: resolution + real*8, dimension(3) :: gdim + real*8 inter(res_x, res_y, res_z,3) + real*8 configuration(res_x+2, res_y+2, res_z+2,3) + real*8 meshgeom(res_x+1, res_y+1, res_z+1,3) + !f2py intent(in) res_x, res_y, res_z + !f2py intent(in) inter + !f2py intent(out) meshgeom + !f2py intent(in) gdim + !f2py depend(res_x, res_y, res_z) inter + !f2py depend(res_x, res_y, res_z) meshgeom + meshgeom = 0.0_pDouble + configuration = 0.0_pDouble + + resolution(1) = res_x; resolution(2) = res_y; resolution(3) = res_z + + + do k = 1, resolution(3)+1; do j = 1, resolution(2)+1; do i = 1, resolution(1)+1 + if((i==1).and.(j==1).and.(k==1)) then + configuration(i,j,k,:)=inter(res_x,res_y,res_z,:)-gdim + configuration(res_x+2,j,k,:)=inter(i,res_y,res_z,:)+(/1.0,-1.0,-1.0/)*gdim + configuration(i,res_y+2,k,:)=inter(res_x,j,res_z,:)+(/-1.0,1.0,-1.0/)*gdim + configuration(i,j,res_z+2,:)=inter(res_x,res_y,k,:)+(/-1.0,-1.0,1.0/)*gdim + configuration(res_x+2,res_y+2,k,:)=inter(i,j,res_z,:)+(/1.0,1.0,-1.0/)*gdim + configuration(i,res_y+2,res_z+2,:)=inter(res_x,j,k,:)+(/-1.0,1.0,1.0/)*gdim + configuration(res_x+2,j,res_z+2,:)=inter(i,res_y,k,:)+(/1.0,-1.0,1.0/)*gdim + configuration(res_x+2,res_y+2,res_z+2,:)=inter(i,j,k,:)+gdim + endif + if((i==1).and.(j==1).and.(k/=1)) then + configuration(1,1,k,:)=inter(res_x,res_y,k-1,:)+(/-1.0,-1.0,0.0/)*gdim + configuration(res_x+2,1,k,:)=inter(1,res_y,k-1,:)+(/1.0,-1.0,0.0/)*gdim + configuration(1,res_y+2,k,:)=inter(res_x,1,k-1,:)+(/-1.0,1.0,0.0/)*gdim + configuration(res_x+2,res_y+2,k,:)=inter(1,1,k-1,:)+(/1.0,1.0,0.0/)*gdim + endif + if((i==1).and.(j/=1).and.(k==1)) then + configuration(1,j,1,:)=inter(res_x,j-1,res_z,:)+(/-1.0,0.0,-1.0/)*gdim + configuration(res_x+2,j,1,:)=inter(1,j-1,res_z,:)+(/1.0,0.0,-1.0/)*gdim + configuration(1,j,res_z+2,:)=inter(res_x,j-1,1,:)+(/-1.0,0.0,1.0/)*gdim + configuration(res_x+2,j,res_z+2,:)=inter(1,j-1,1,:)+(/1.0,0.0,1.0/)*gdim + endif + if((i/=1).and.(j==1).and.(k==1)) then + configuration(i,1,1,:)=inter(i-1,res_y,res_z,:)+(/0.0,-1.0,-1.0/)*gdim + configuration(i,1,res_z+2,:)=inter(i-1,res_y,1,:)+(/0.0,-1.0,1.0/)*gdim + configuration(i,res_y+2,1,:)=inter(i-1,1,res_z,:)+(/0.0,1.0,-1.0/)*gdim + configuration(i,res_y+2,res_z+2,:)=inter(i-1,1,1,:)+(/0.0,1.0,1.0/)*gdim + endif + if((i/=1).and.(j/=1).and.(k==1)) then + configuration(i,j,1,:)=inter(i-1,j-1,res_z,:)+(/0.0,0.0,-1.0/)*gdim + configuration(i,j,res_z+2,:)=inter(i-1,j-1,1,:)+(/0.0,0.0,1.0/)*gdim + endif + if((i==1).and.(j/=1).and.(k/=1)) then + configuration(1,j,k,:)=inter(res_x,j-1,k-1,:)+(/-1.0,0.0,0.0/)*gdim + configuration(res_x+2,j,k,:)=inter(i,j-1,k-1,:)+(/1.0,0.0,0.0/)*gdim + endif + if((i/=1).and.(j==1).and.(k/=1)) then + configuration(i,1,k,:)=inter(i-1,res_y,k-1,:)+(/0.0,-1.0,0.0/)*gdim + configuration(i,res_y+2,k,:)=inter(i-1,1,k-1,:)+(/0.0,1.0,0.0/)*gdim + endif + if((i/=1).and.(j/=1).and.(k/=1)) then + configuration(i,j,k,:)=inter(i-1,j-1,k-1,:) + endif + enddo; enddo; enddo + + do k = 1, resolution(3)+1; do j = 1, resolution(2)+1; do i = 1, resolution(1)+1 + meshgeom(i,j,k,:)=((configuration(i,j,k,:) +configuration(i+1,j,k,:))& + +(configuration(i,j+1,k,:) +configuration(i,j,k+1,:))& + +(configuration(i+1,j+1,k,:) +configuration(i,j+1,k+1,:))& + +(configuration(i+1,j,k+1,:) +configuration(i+1,j+1,k+1,:)))/8.0_pDouble + enddo; enddo; enddo +end subroutine mesh + + +!below some code I used for gmsh postprocessing. Might be helpful +!!gmsh output + ! character(len=1024) :: nriter + ! character(len=1024) :: nrstep + ! character(len=1024) :: nrloadcase + ! real(pReal), dimension(:,:,:,:), allocatable :: displacement + ! real(pReal), dimension(3,3) :: temp33_Real2 + ! real(pReal), dimension(3,3,3) :: Eigenvectorbasis + ! real(pReal), dimension(3) :: Eigenvalue + ! real(pReal) determinant +!!gmsh output +!!Postprocessing (gsmh output) + ! do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) + ! if((k==1).and.(j==1).and.(i==1)) then + ! temp33_Real =0.0_pReal + ! else + ! if((j==1).and.(i==1)) then + ! temp33_Real(1,:) = temp33_Real(1,:) + math_mul33x3(defgrad(i,j,k,:,:),& + ! (/0.0_pReal,0.0_pReal,(real(resolution(3))/meshdimension(3))/)) + ! temp33_Real(2,:) = temp33_Real(1,:) + ! temp33_Real(3,:) = temp33_Real(1,:) + ! displacement(i,j,k,:) = temp33_Real(1,:) + ! else + ! if(i==1) then + ! temp33_Real(2,:) = temp33_Real(2,:) + math_mul33x3(defgrad(i,j,k,:,:),& + ! (/0.0_pReal,(real(resolution(2))/meshdimension(2)),0.0_pReal/)) + ! temp33_Real(3,:) = temp33_Real(2,:) + ! displacement(i,j,k,:) = temp33_Real(2,:) + ! else + ! temp33_Real(3,:) = temp33_Real(3,:) + math_mul33x3(defgrad(i,j,k,:,:),& + ! (/(real(resolution(1))/meshdimension(1)),0.0_pReal,0.0_pReal/)) + ! displacement(i,j,k,:) = temp33_Real(3,:) + ! endif + ! endif + ! endif + ! enddo; enddo; enddo + + ! write(nrloadcase, *) loadcase; write(nriter, *) iter; write(nrstep, *) steps + ! open(589,file = 'stress' //trim(adjustl(nrloadcase))//'-'//trim(adjustl(nrstep))//'-'//trim(adjustl(nriter))//'_cpfem.msh') + ! open(588,file = 'logstrain'//trim(adjustl(nrloadcase))//'-'//trim(adjustl(nrstep))//'-'//trim(adjustl(nriter))//'_cpfem.msh') + ! write(589, '(4(A, /), I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', prodnn + ! write(588, '(4(A, /), I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', prodnn + + ! ielem = 0_pInt + ! do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) + ! ielem = ielem + 1 + ! write(589, '(I10, 3(tr2, E12.6))'), ielem, displacement(i,j,k,:) !for deformed configuration + ! write(588, '(I10, 3(tr2, E12.6))'), ielem, displacement(i,j,k,:) + !! write(589, '(4(I10,tr2))'), ielem, i-1,j-1,k-1 !for undeformed configuration + !!write(588, '(4(I10,tr2))'), ielem, i-1,j-1,k-1 + ! enddo; enddo; enddo + + ! write(589, '(2(A, /), I10)'), '$EndNodes', '$Elements', prodnn + ! write(588, '(2(A, /), I10)'), '$EndNodes', '$Elements', prodnn + + ! do i = 1, prodnn + ! write(589, '(I10, A, I10)'), i, ' 15 2 1 2', i + ! write(588, '(I10, A, I10)'), i, ' 15 2 1 2', i + ! enddo + + ! write(589, '(A)'), '$EndElements' + ! write(588, '(A)'), '$EndElements' + ! write(589, '(8(A, /), I10)'), '$NodeData', '1','"'//trim(adjustl('stress'//trim(adjustl(nrloadcase))//'-'//& + ! trim(adjustl(nrstep))//'-'//trim(adjustl(nriter))//'_cpfem.msh'))//'"','1','0.0', '3', '0', '9', prodnn + ! write(588, '(8(A, /), I10)'), '$NodeData', '1','"'//trim(adjustl('logstrain'//trim(adjustl(nrloadcase))//'-'//& + ! trim(adjustl(nrstep))//'-'//trim(adjustl(nriter))//'_cpfem.msh'))//'"','1','0.0', '3', '0', '9', prodnn + ! ielem = 0_pInt + ! do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) + ! ielem = ielem + 1 + ! write(589, '(i10, 9(tr2, E14.8))'), ielem, cstress_field(i,j,k,:,:) + ! call math_pDecomposition(defgrad(i,j,k,:,:),temp33_Real2,temp33_Real,errmatinv) !store R in temp33_Real + ! call math_invert3x3(temp33_Real, temp33_Real2, determinant, errmatinv) !inverse of R in temp33_Real2 + ! temp33_Real = math_mul33x33(defgrad(i,j,k,:,:), temp33_Real2) ! v = F o inv(R), store in temp33_Real + ! call math_spectral1(temp33_Real,Eigenvalue(1), Eigenvalue(2), Eigenvalue(3),& + ! Eigenvectorbasis(1,:,:),Eigenvectorbasis(2,:,:),Eigenvectorbasis(3,:,:)) + ! eigenvalue = log(sqrt(eigenvalue)) + ! temp33_Real = eigenvalue(1)*Eigenvectorbasis(1,:,:)+eigenvalue(2)*Eigenvectorbasis(2,:,:)+eigenvalue(3)*Eigenvectorbasis(3,:,:) + ! write(588, '(i10, 9(tr2, E14.8))'), ielem, temp33_Real + ! enddo; enddo; enddo + + ! write(589, *), '$EndNodeData' + ! write(588, *), '$EndNodeData' + ! close(589); close(588); close(540) + ! enddo ! end looping over steps in current loadcase + ! enddo ! end looping over loadcases +! close(539); close(538) diff --git a/processing/post/spectral_post.py b/processing/post/spectral_post.py index 04104a909..353d86bb5 100644 --- a/processing/post/spectral_post.py +++ b/processing/post/spectral_post.py @@ -1,100 +1,115 @@ #!/usr/bin/python # -*- coding: iso-8859-1 -*- + +# This script is used for the post processing of the results achieved by the spectral method. +# As it reads in the data coming from "materialpoint_results", it can be adopted to the data +# computed using the FEM solvers. Until now, its capable to handle elements with one IP in a regular order +# written by M. Diehl, m.diehl@mpie.de + import array import struct -#import numpy -print('post processing for mpie_spectral') +import numpy +import reconstruct -filename ='results.out' -results = open(filename, 'rb') -print('filename:', filename) -position=0 +#this funtion finds a string value in the given file for a given keyword, also returns the position +def searchNameForKeyword(searchstring, file): + file.seek(0,0) + begin = file.read(2048).find(searchstring) + len(searchstring) # function will return -1 if string is not found + file.seek(begin -len(searchstring) -4) # position of header (Fortran specific, 4 bit long on msuws4 + header = file.read(4) # store header. header will be repeated, thus it allows us to determine the length of the searchstring + file.seek(begin,0) # go back to starting postion + length = file.read(2048).find(header) # find header(footer) a second time + file.seek(begin,0) + return file.read(length), results.tell() # return searchstring and maximum position in file -#this funtion finds a string value in the given file for a given keyword -def searchNameForKeyword(searchstring, file, maxpos): +#this funtion finds an integer value in the given file for a given keyword, works similar as "searchNameForKeyword" +def searchValueForKeyword(searchstring, file): results.seek(0,0) - begin = file.read(2048).find(searchstring) + len(searchstring) # function will return -1 if string isnt found - file.seek(begin -len(searchstring) -4) #position of header - header = file.read(4) #store header - file.seek(begin,0) - length = file.read(2048).find(header) - file.seek(begin,0) - name = file.read(length) - return name, max(maxpos,results.tell()) - -#this funtion finds an integer value in the given file for a given keyword -def searchValueForKeyword(searchstring, file, maxpos): - results.seek(0,0) - begin = file.read(2048).find(searchstring) + len(searchstring) # function will return -1 if string isnt found - file.seek(begin,0) - value = array.array('i',[0]) + begin = file.read(2048).find(searchstring) + len(searchstring) + file.seek(begin,0) value = struct.unpack('i',results.read(4))[0] - return value, max(maxpos,results.tell()) + return value, results.tell() -#finds the value for the three integers followed the given keyword -def searchArrayForKeyword(searchstring, file, maxpos): - values = array.array('i',[0,0,0]) - results.seek(0,0) - begin = results.read(2048).find(searchstring) + len(searchstring) #end position of string "resolution" - pos = results.read(60).find(b'a') - results.seek(begin+pos+2,0) #why two, not 1?? - values[0]=struct.unpack('i',results.read(4))[0] - maxpos=max(maxpos,results.tell()) - results.seek(begin,0) - pos = results.read(60).find(b'b') - results.seek(begin+pos+1,0) - values[1]=struct.unpack('i',results.read(4))[0] - maxpos=max(maxpos,results.tell()) - results.seek(begin,0) - pos = results.read(60).find(b'c') - results.seek(begin+pos+1,0) - values[2]=struct.unpack('i',results.read(4))[0] - maxpos=max(maxpos,results.tell()) - return values, maxpos +#finds the value for the three integers (a,b,c) following the given keyword +def searchIntegerArrayForKeyword(searchstring, file): + values = array.array('i',[0,0,0]) + file.seek(0,0) + begin = file.read(2048).find(searchstring) + len(searchstring) #end position of string "searchstring" + pos = file.read(60).find('a') + file.seek(begin+pos+2,0) #why 2, not 1?? + values[0]=struct.unpack('i',file.read(4))[0] + maxpos=file.tell() + file.seek(begin,0) + pos = file.read(60).find('b') + file.seek(begin+pos+1,0) + values[1]=struct.unpack('i',file.read(4))[0] + maxpos = max(maxpos,file.tell()) + file.seek(begin,0) + pos = file.read(60).find('c') + file.seek(begin+pos+1,0) + values[2]=struct.unpack('i',file.read(4))[0] + maxpos=max(maxpos,file.tell()) + return values, maxpos -#finds the value for the three doubles followed the given keyword -def searchFarrayForKeyword(searchstring, file, maxpos): - values = array.array('d',[0,0,0]) - results.seek(0,0) - begin = results.read(2048).find(searchstring) + len(searchstring) #end position of string "resolution" - pos = results.read(60).find(b'x') - results.seek(begin+pos+2,0) #why two, not 1?? - values[0]=struct.unpack('d',results.read(8))[0] - maxpos=max(maxpos,results.tell()) - results.seek(begin,0) - pos = results.read(60).find(b'y') - results.seek(begin+pos+1,0) - values[1]=struct.unpack('d',results.read(8))[0] - maxpos=max(maxpos,results.tell()) - results.seek(begin,0) - pos = results.read(60).find(b'z') - results.seek(begin+pos+1,0) - values[2]=struct.unpack('d',results.read(8))[0] - maxpos=max(maxpos,results.tell()) - return values, maxpos +#finds the value for the three doubles (x,y,z) following the given keyword +def searchDoubleArrayForKeyword(searchstring, file): + values = array.array('d',[0,0,0]) + file.seek(0,0) + begin = file.read(2048).find(searchstring) + len(searchstring) + pos = file.read(60).find('x') + file.seek(begin+pos+2,0) + values[0]=struct.unpack('d',file.read(8))[0] + maxpos=file.tell() + file.seek(begin,0) + pos = file.read(60).find('y') + file.seek(begin+pos+1,0) + values[1] = struct.unpack('d',file.read(8))[0] + maxpos = max(maxpos,file.tell()) + file.seek(begin,0) + pos = file.read(60).find('z') + file.seek(begin+pos+1,0) + values[2] = struct.unpack('d',file.read(8))[0] + maxpos = max(maxpos,file.tell()) + return values, maxpos +print '*********************************************************************************' +print 'Post Processing for Material subroutine for BVP solution using spectral method' +print '*********************************************************************************\n' -loadcase, position = searchNameForKeyword(b'Loadcase ', results, position) -workingdir, position = searchNameForKeyword(b'Workingdir ', results, position) -jobname, position = searchNameForKeyword(b'JobName ', results, position) -totalincs, position = searchValueForKeyword(b'totalincs ', results, position) -materialpoint_sizeResults, position = searchValueForKeyword(b'materialpoint_sizeResults ', results, position) -resolution, position = searchArrayForKeyword(b'resolution ', results, position) -geomdimension, position = searchFarrayForKeyword(b'geomdimension ', results, position) -position=position+4 +13*8 -results.seek(position,0) -tnsr = array.array('d',[0,0,0,0,0,0,0,0,0]) -tnsr[0]=struct.unpack('d',results.read(8))[0] -tnsr[1]=struct.unpack('d',results.read(8))[0] -tnsr[2]=struct.unpack('d',results.read(8))[0] -tnsr[3]=struct.unpack('d',results.read(8))[0] -tnsr[4]=struct.unpack('d',results.read(8))[0] -print(tnsr) +#reading in the header of the results file +filename ='32x32x32x100.out' +print 'Results Filename:', filename +results = open(filename, 'rb') + +loadcase, position = searchNameForKeyword('Loadcase', results) +workingdir, temp = searchNameForKeyword('Workingdir', results) +position = max(temp, position) +jobname, temp = searchNameForKeyword('JobName', results) +position = max(temp, position) +totalincs, temp = searchValueForKeyword('totalincs', results) +position = max(temp, position) +materialpoint_sizeResults, temp = searchValueForKeyword('materialpoint_sizeResults', results) +position = max(temp, position) +resolution, temp = searchIntegerArrayForKeyword('resolution', results) +position = max(temp, position) +geomdimension, temp = searchDoubleArrayForKeyword('geomdimension', results) + +print 'Load case:', loadcase +print 'Workingdir:', workingdir +print 'Job Name:', jobname +print 'Total No. of Increments:', totalincs +print 'Materialpoint_sizeResults:', materialpoint_sizeResults +print 'Resolution:', resolution +print 'Geomdimension', geomdimension +print 'Position in File:', position + +# Ended reading of header +# Now starting to read information concerning output of materialpoint_results(:,1,:) + +filename = jobname[0:len(jobname)-5]+'.outputCrystallite' +outputCrystallite = open(filename, 'r') -# ended reading of header. now starting to read information concerning output -#filename = jobname[0:len(jobname)-5]+b'.outputCrystallite' -#outputCrystallite = open(filename, 'r') - +# some funtions that might be useful def InCharCount(location, character): subj = file(location, "r") body = subj.read() diff --git a/processing/pre/voronoi fast.f90 b/processing/pre/voronoi_fast.f90 similarity index 92% rename from processing/pre/voronoi fast.f90 rename to processing/pre/voronoi_fast.f90 index 0c676bd2e..e86dc9bd9 100644 --- a/processing/pre/voronoi fast.f90 +++ b/processing/pre/voronoi_fast.f90 @@ -18,7 +18,7 @@ program voronoi print*, '******************************************************************************' print*, '' print*, 'generates:' - print*, ' * mesh file "_GIVEN_NAME_.mesh": Geometrical information for solver' + 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:' @@ -37,7 +37,7 @@ program voronoi read(*, *), c write(*, '(A)', advance = 'NO') 'Please enter No. of Grains: ' read(*, *), N_Seeds - write(*, '(A)', advance = 'NO') 'Please enter name of mesh file: ' + 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 @@ -55,7 +55,7 @@ program voronoi 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, '(A)'), 'crystallite 1' write(20, trim(format2)), '(constituent) phase 1 texture ', i, ' fraction 1.0' end do @@ -77,8 +77,8 @@ program voronoi 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 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' @@ -163,16 +163,16 @@ program voronoi 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' + &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' + 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 mesh out +! 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 @@ -190,7 +190,7 @@ program voronoi 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 @@ -213,7 +213,7 @@ program voronoi write(20, *), '$EndNodeData' close(20) -! write 2d mesh out +! 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 diff --git a/processing/pre/voronoi small memory.f90 b/processing/pre/voronoi_small_memory.f90 similarity index 90% rename from processing/pre/voronoi small memory.f90 rename to processing/pre/voronoi_small_memory.f90 index ccb444670..24cd6eab7 100644 --- a/processing/pre/voronoi small memory.f90 +++ b/processing/pre/voronoi_small_memory.f90 @@ -18,7 +18,7 @@ program voronoi print*, '******************************************************************************' print*, '' print*, 'generates:' - print*, ' * mesh file "_GIVEN_NAME_.mesh": Geometrical information for solver' + 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' @@ -36,7 +36,7 @@ program voronoi read(*, *), c write(*, '(A)', advance = 'NO') 'Please enter No. of Grains: ' read(*, *), N_Seeds - write(*, '(A)', advance = 'NO') 'Please enter name of mesh file: ' + 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 @@ -53,7 +53,7 @@ program voronoi 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, '(A)'), 'crystallite 1' write(20, trim(format2)), '(constituent) phase 1 texture ', i, ' fraction 1.0' end do @@ -75,8 +75,8 @@ program voronoi 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 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' @@ -128,8 +128,8 @@ program voronoi 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. + &(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) @@ -168,13 +168,13 @@ program voronoi end do close(20) print*, 'voronoi tesselation finished' - print*, 'mesh file is written out' + 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 mesh out + +! 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 @@ -194,7 +194,7 @@ program voronoi write(20, *), '$EndNodeData' close(20) -! write 2d mesh out +! 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