From 17b8205e3ff86289b4973e67355e988fb98438d8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 26 Jan 2011 12:56:52 +0000 Subject: [PATCH] reconstruction of geometry is now working. It is directly implemented in python (file spectral_post.py). reconstruction in fortran is not working (file reconstruct.f90) due to some problems with f2py --- processing/post/reconstruct.f90 | 448 +++++++++++++++++--------- processing/post/spectral_post.py | 524 ++++++++++++++++++++++++------- 2 files changed, 713 insertions(+), 259 deletions(-) diff --git a/processing/post/reconstruct.f90 b/processing/post/reconstruct.f90 index 753d4c57a..8b116858b 100644 --- a/processing/post/reconstruct.f90 +++ b/processing/post/reconstruct.f90 @@ -1,193 +1,306 @@ ! -*- f90 -*- -subroutine simple(defgrad,res_x,res_y,res_z,geomdimension,current_configuration) +function coordinates2(res_x,res_y,res_z,geomdimension,defgrad) implicit none integer, parameter :: pDouble = selected_real_kind(15,50) - integer i,j,k - integer res_x, res_y, res_z + integer i,j,k, l,m, s,o, loop, res_x, res_y, res_z + integer, dimension(3) :: res, init, oppo, me, rear + real*8, dimension(3,3) :: defgrad_av integer, dimension(3) :: resolution - real*8, dimension(3) :: geomdimension + real*8, dimension(3) :: geomdimension, myStep 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 + real*8, dimension(3,res_x, res_y, res_z) :: coordinates2 + real*8, dimension(3,3,res_x, res_y, res_z) :: defgrad + real*8, dimension(3,res_x,res_y,res_z,8) :: cornerCoords + real*8, dimension(3,2+res_x,2+res_y,2+res_z,6,8) :: coord + !f2py intent(in) res_x + !f2py intent(in) res_y + !f2py intent(in) 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 intent(out) coordinates2 + !f2py depend(res_x, res_y, res_z) coordinates2 !f2py depend(res_x, res_y, res_z) defgrad - resolution(1) = res_x; resolution(2) = res_y; resolution(3) = res_z + ! integer, dimension(3,8) :: corner = reshape((/ & + ! 0, 0, 0,& + ! 1, 0, 0,& + ! 1, 1, 0,& + ! 0, 1, 0,& + ! 1, 1, 1,& + ! 0, 1, 1,& + ! 0, 0, 1,& + ! 1, 0, 1 & + ! /), & + ! (/3,8/)) + + ! integer, dimension(3,8) :: step = reshape((/ & + ! 1, 1, 1,& + ! -1, 1, 1,& + ! -1,-1, 1,& + ! 1,-1, 1,& + ! -1,-1,-1,& + ! 1,-1,-1,& + ! 1, 1,-1,& + ! -1, 1,-1 & + ! /), & + ! (/3,8/)) + + ! integer, dimension(3,6) :: order = reshape((/ & + ! 1, 2, 3,& + ! 1, 3, 2,& + ! 2, 1, 3,& + ! 2, 3, 1,& + ! 3, 1, 2,& + ! 3, 2, 1 & + ! /), & + ! (/3,6/)) + + resolution = (/res_x,res_y,res_z/) + + write(6,*) 'defgrad', 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 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_pDouble + temp33_Real = real(0.0) else if((j==1).and.(i==1)) then - temp33_Real(1,:) = temp33_Real(1,:) + matmul(defgrad(i,j,k,:,:),& - (/0.0_pDouble,0.0_pDouble,geomdimension(3)/(real(resolution(3)))/)) + temp33_Real(1,:) = temp33_Real(1,:) + matmul(defgrad(:,:,i,j,k),& + (/real(0.0),real(0.0),real(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,:) + coordinates2(:,i,j,k) = temp33_Real(1,:) else if(i==1) then - temp33_Real(2,:) = temp33_Real(2,:) + matmul(defgrad(i,j,k,:,:),& - (/0.0_pDouble,geomdimension(2)/(real(resolution(2))),0.0_pDouble/)) + temp33_Real(2,:) = temp33_Real(2,:) + matmul(defgrad(:,:,i,j,k),& + (/real(0.0),real(geomdimension(2))/(real(resolution(2))),real(0.0)/)) temp33_Real(3,:) = temp33_Real(2,:) - current_configuration(i,j,k,:) = temp33_Real(2,:) + coordinates2(:,i,j,k) = temp33_Real(2,:) else - temp33_Real(3,:) = temp33_Real(3,:) + matmul(defgrad(i,j,k,:,:),& - (/geomdimension(1)/(real(resolution(1))),0.0_pDouble,0.0_pDouble/)) - current_configuration(i,j,k,:) = temp33_Real(3,:) + temp33_Real(3,:) = temp33_Real(3,:) + matmul(defgrad(:,:,i,j,k),& + (/real(geomdimension(1))/(real(resolution(1))),real(0.0),real(0.0)/)) + coordinates2(:,i,j,k) = temp33_Real(3,:) endif endif endif - !print*, current_configuration enddo; enddo; enddo -end subroutine simple + do i=1, res_x; do j = 1, res_y; do k = 1, res_z + coordinates2(:,i,j,k) = coordinates2(:,i,j,k)+ matmul(defgrad_av,(/geomdimension(1)/real(res_x),geomdimension(2)/real(res_y),geomdimension(3)/real(res_z)/)) + enddo; enddo; enddo + + res = (/res_x,res_y,res_z/) + do i=1,3; do j=1,3 + defgrad_av(i,j) = sum(defgrad(i,j,:,:,:)) /real(res(1)*res(2)*res(3)) + enddo; enddo -subroutine advanced(defgrad,res_x,res_y,res_z,geomdimension,current_configuration) + ! do s = 1,8 + ! init = corner(:,s)*(res-(/1,1,1/)) + ! oppo = corner(:,1+mod(s-1+4,8))*(res-(/1,1,1/)) + ! do o = 1,6 + ! do k = init(order(3,o)),oppo(order(3,o)),step(order(3,o),s) + ! rear(order(2,o)) = init(order(2,o)) + ! do j = init(order(2,o)),oppo(order(2,o)),step(order(2,o),s) + ! rear(order(1,o)) = init(order(1,o)) + ! do i = init(order(1,o)),oppo(order(1,o)),step(order(1,o),s) + ! ! print*, order(:,o) + ! me(order(1,o)) = i + ! me(order(2,o)) = j + ! me(order(3,o)) = k + ! ! print*, me +! ! if (all(me == init)) then +! ! coord(:,1+me(1),1+me(2),1+me(3),o,s) = 0.0 !& + ! ! geomdimension*(matmul(defgrad_av,real(corner(:,s),pDouble)) + & + ! ! matmul(defGrad(:,:,1+me(1),1+me(2),1+me(3)),step(:,s)/res/2.0_pDouble)) + ! ! else + ! ! myStep = (me-rear)*geomdimension/res + ! ! coord(:,1+me(1),1+me(2),1+me(3),o,s) = coord(:,1+rear(1),1+rear(2),1+rear(3),o,s) + & + ! ! 0.5_pDouble*matmul(defGrad(:,:,1+me(1),1+me(2),1+me(3)) + & + ! ! defGrad(:,:,1+rear(1),1+rear(2),1+rear(3)),& + ! ! myStep) + ! ! endif + ! ! rear = me + ! enddo + ! enddo + ! enddo + ! enddo ! orders + ! ! cornerCoords(:,:,:,:,s) = sum(coord(:,:,:,:,:,s),5)/6.0_pDouble + ! enddo ! corners + + ! coordinates = sum(cornerCoords(:,:,:,:,:),5)/8.0_pDouble ! plain average no shape functions... +end function + + + +! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine coordinates5(res_x,res_y,res_z,geomdimension,defgrad) +! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 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 + integer i,j,k, l,m, s,o, loop, res_x, res_y, res_z + integer, dimension(3) :: res, init, oppo, me, rear 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 + real*8, dimension(3) :: myStep + real*8, dimension(3,3) :: defGrad_av + integer, dimension(3,8) :: corner = reshape((/ & + 0, 0, 0,& + 1, 0, 0,& + 1, 1, 0,& + 0, 1, 0,& + 1, 1, 1,& + 0, 1, 1,& + 0, 0, 1,& + 1, 0, 1 & + /), & + (/3,8/)) + + integer, dimension(3,8) :: step = reshape((/ & + 1, 1, 1,& + -1, 1, 1,& + -1,-1, 1,& + 1,-1, 1,& + -1,-1,-1,& + 1,-1,-1,& + 1, 1,-1,& + -1, 1,-1 & + /), & + (/3,8/)) + + integer, dimension(3,6) :: order = reshape((/ & + 1, 2, 3,& + 1, 3, 2,& + 2, 1, 3,& + 2, 3, 1,& + 3, 1, 2,& + 3, 2, 1 & + /), & + (/3,6/)) + + real*8 defGrad(3,3,res_x,res_y,res_z) + real*8 coordinates(3,res_x,res_y,res_z) + real*8, dimension(3,res_x,res_y,res_z,8) :: cornerCoords + real*8, dimension(3,2+res_x,2+res_y,2+res_z,6,8) :: coord !f2py intent(in) res_x, res_y, res_z !f2py intent(in) geomdimension - !f2py intent(out) current_configuration + !f2py intent(out) coordinates !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 + !f2py depend(res_x, res_y, res_z) coordinates + !f2py depend(res_x, res_y, res_z) defgrad + !f2py depend(res_x, res_y, res_z) cornerCoords + !f2py depend(res_x, res_y, res_z) coord -subroutine mesh(inter,res_x,res_y,res_z,gdim,meshgeom) + res = (/res_x,res_y,res_z/) + do i=1,3; do j=1,3 + defgrad_av(i,j) = sum(defgrad(i,j,:,:,:)) /real(res(1)*res(2)*res(3)) + enddo; enddo + print*, 'defgra', defgrad + do s = 1,8 + init = corner(:,s)*(res-(/1,1,1/)) + oppo = corner(:,1+mod(s-1+4,8))*(res-(/1,1,1/)) + do o = 1,6 + do k = init(order(3,o)),oppo(order(3,o)),step(order(3,o),s) + rear(order(2,o)) = init(order(2,o)) + do j = init(order(2,o)),oppo(order(2,o)),step(order(2,o),s) + rear(order(1,o)) = init(order(1,o)) + do i = init(order(1,o)),oppo(order(1,o)),step(order(1,o),s) + me(order(1,o)) = i + me(order(2,o)) = j + me(order(3,o)) = k + if (all(me == init)) then + coord(:,1+me(1),1+me(2),1+me(3),o,s) = & + geomdimension*(matmul(defgrad_av,real(corner(:,s),pDouble)) + & + matmul(defGrad(:,:,1+me(1),1+me(2),1+me(3)),step(:,s)/res/2.0_pDouble)) + else + myStep = (me-rear)*geomdimension/res + coord(:,1+me(1),1+me(2),1+me(3),o,s) = coord(:,1+rear(1),1+rear(2),1+rear(3),o,s) + & + 0.5_pDouble*matmul(defGrad(:,:,1+me(1),1+me(2),1+me(3)) + & + defGrad(:,:,1+rear(1),1+rear(2),1+rear(3)),& + myStep) + endif + rear = me + enddo + enddo + enddo + enddo ! orders + cornerCoords(:,:,:,:,s) = sum(coord(:,:,:,:,:,s),5)/6.0_pDouble + enddo ! corners + + coordinates = sum(cornerCoords(:,:,:,:,:),5)/8.0_pDouble ! plain average no shape functions... + +end subroutine + +! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +subroutine mesh(res_x,res_y,res_z,geomdim,defgrad_av,centroids,nodes) +! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ implicit none integer, parameter :: pDouble = selected_real_kind(15,50) - integer i,j,k + integer i,j,k, n 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) + integer, dimension(3) :: res, shift, lookup, me + real*8, dimension(3) :: geomdim, diag = (/1,1,1/) + real*8, dimension(3,3) :: defgrad_av + real*8, dimension(3,res_x, res_y, res_z) :: centroids + real*8, dimension(3,res_x+2, res_y+2, res_z+2) :: wrappedCentroids + real*8, dimension(3,res_x+1, res_y+1, res_z+1) :: nodes + integer, dimension(3,8) :: neighbor = reshape((/ & + 0, 0, 0,& + 1, 0, 0,& + 1, 1, 0,& + 0, 1, 0,& + 0, 0, 1,& + 1, 0, 1,& + 1, 1, 1,& + 0, 1, 1 & + /), & + (/3,8/)) !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 + !f2py intent(in) centroids + !f2py intent(in) defgrad_av + !f2py intent(in) geomdim + !f2py intent(out) nodes + !f2py depend(res_x, res_y, res_z) centroids + !f2py depend(res_x, res_y, res_z) nodes + + res = (/res_x,res_y,res_z/) + wrappedCentroids(:,2:res(1)+1,2:res(2)+1,2:res(3)+1) = centroids + + do k = 0,res(3)+1 + do j = 0,res(2)+1 + do i = 0,res(1)+1 + if (& + k==0 .or. k==res(3)+1 .or. & + j==0 .or. j==res(2)+1 .or. & + i==0 .or. i==res(1)+1 & + ) then + me = (/i,j,k/) + shift = (res+diag-2*me)/(res+diag) + lookup = me+shift*res + wrappedCentroids(:,1+i,1+j,1+k) = centroids(:,lookup(1),lookup(2),lookup(3)) - & + matmul(defgrad_av,shift * geomdim) + endif + enddo + enddo + enddo + + nodes = 0.0_pDouble + + do k = 0,res(3) + do j = 0,res(2) + do i = 0,res(1) + + do n = 1,8 + nodes(:,1+i,1+j,1+k) = nodes(:,1+i,1+j,1+k) + wrappedCentroids(:,1+i+neighbor(1,n),& + 1+j+neighbor(2,n),& + 1+k+neighbor(3,n)) + enddo + nodes(:,1+i,1+j,1+k) = nodes(:,1+i,1+j,1+k) / 8.0_pDouble + + 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 @@ -277,3 +390,32 @@ end subroutine mesh ! enddo ! end looping over steps in current loadcase ! enddo ! end looping over loadcases ! close(539); close(538) + + + +! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +function coordinates3(res_x,res_y,res_z,geomdimension,defgrad) +! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + implicit none + integer i,j,k, l,m, s,o, loop, res_x, res_y, res_z + real*8 defgrad_av(3,3) + real*8 geomdimension(3) + integer res(3) + real*8 defgrad(3,3,res_x,res_y,res_z) + real*8 coordinates3(3,res_x,res_y,res_z) + + !f2py intent(in) res_x, res_y, res_z + !f2py depend(res_x, res_y, res_z) coordinates3 + !f2py depend(res_x, res_y, res_z) defgrad + + !f2py intent(in) geomdimension + !f2py intent(out) coordinates3 + !f2py intent(in) defgrad + + res = (/res_x,res_y,res_z/) + do i=1,3; do j=1,3 + defgrad_av(i,j) = sum(defgrad(i,j,:,:,:)) /real(res(1)*res(2)*res(3)) + enddo; enddo + print*, 'defgra', defgrad + coordinates3 = 0.0 +end function \ No newline at end of file diff --git a/processing/post/spectral_post.py b/processing/post/spectral_post.py index 353d86bb5..1894674a0 100644 --- a/processing/post/spectral_post.py +++ b/processing/post/spectral_post.py @@ -6,124 +6,436 @@ # 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 -import reconstruct +import os,sys,re,array,struct,numpy, time -#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 +class vector: + x,y,z = [None,None,None] + + def __init__(self,coords): + self.x = coords[0] + self.y = coords[1] + self.z = coords[2] -#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) - file.seek(begin,0) - value = struct.unpack('i',results.read(4))[0] - return value, results.tell() +class element: + items = [] + type = None -#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 + def __init__(self,nodes,type): + self.items = nodes + self.type = type -#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 '*********************************************************************************' +class element_scalar: + id = None + value = None + + def __init__(self,node,value): + self.id = node + self.value = value + + +class MPIEspectral_result: + + file = None + dataOffset = 0 + N_elemental_scalars = 0 + resolution = [0,0,0] + dimension = [0.0,0.0,0.0] + theTitle = '' + wd = '' + extrapolate = '' + N_increments = 0 + increment = 0 + N_nodes = 0 + N_node_scalars = 0 + N_elements = 0 + N_element_scalars = 0 + N_element_tensors = 0 + theNodes = [] + theElements = [] + + def __init__(self,filename): + + self.file = open(filename, 'rb') + + self.title = self._keyedString('load') + self.wd = self._keyedString('workingdir') + self.geometry = self._keyedString('geometry') + self.N_increments = self._keyedInt('increments') + self.N_element_scalars = self._keyedInt('materialpoint_sizeResults') + self.resolution = self._keyedPackedArray('resolution',3,'i') + self.N_nodes = (self.resolution[0]+1)*(self.resolution[1]+1)*(self.resolution[2]+1) + self.N_elements = self.resolution[0]*self.resolution[1]*self.resolution[2] + self.dimension = self._keyedPackedArray('dimension',3,'d') + a = self.resolution[0]+1 + b = self.resolution[1]+1 + c = self.resolution[2]+1 + self.file.seek(0) + self.dataOffset = self.file.read(2048).find('eoh')+7 + + + def __str__(self): + return '\n'.join([ + 'title: %s'%self.title, + 'workdir: %s'%self.wd, + 'extrapolation: %s'%self.extrapolate, + 'increments: %i'%self.N_increments, + 'increment: %i'%self.increment, + 'nodes: %i'%self.N_nodes, + 'resolution: %s'%(','.join(map(str,self.resolution))), + 'dimension: %s'%(','.join(map(str,self.dimension))), + 'elements: %i'%self.N_elements, + 'nodal_scalars: %i'%self.N_node_scalars, + 'elemental scalars: %i'%self.N_element_scalars, + 'end of header: %i'%self.dataOffset, + ] + ) + + def _keyedPackedArray(self,identifier,length = 3,type = 'd'): + match = {'d': 8,'i': 4} + self.file.seek(0) + m = re.search('%s%s'%(identifier,'(.{%i})'%(match[type])*length),self.file.read(2048)) + values = [] + if m: + for i in m.groups(): + values.append(struct.unpack(type,i)[0]) + return values + + def _keyedInt(self,identifier): + value = None + self.file.seek(0) + m = re.search('%s%s'%(identifier,'(.{4})'),self.file.read(2048)) + if m: + value = struct.unpack('i',m.group(1))[0] + return value + + def _keyedString(self,identifier): + value = None + self.file.seek(0) + m = re.search(r'(.{4})%s(.*?)\1'%identifier,self.file.read(2048)) + if m: + value = m.group(2) + return value + + def extrapolation(self,value): + self.extrapolate = value + + def element_scalar(self,elem,idx): + self.file.seek(self.dataOffset+(self.increment*(4+self.N_elements*self.N_element_scalars*8+4) + 4+(elem*self.N_element_scalars + idx)*8)) + value = struct.unpack('d',self.file.read(8))[0] + return [elemental_scalar(node,value) for node in self.theElements[elem].items] + +def readScalar(resolution,file,distance,startingPosition,offset): + currentPosition = startingPosition+offset*8+4 - distance*8 # we add distance later on + field = numpy.zeros([resolution[0],resolution[1],resolution[2]], 'd') + for x in range(0,resolution[0]): + for y in range(0,resolution[1]): + for z in range(0,resolution[2]): + currentPosition = currentPosition + distance*8 + p.file.seek(currentPosition) + field[x][y][z]=struct.unpack('d',p.file.read(8))[0] + return field + +def readTensor(resolution,file,distance,startingPosition,offset): + currentPosition = startingPosition+offset*8+4 - distance*8 # we add distance later on + field = numpy.zeros([resolution[0],resolution[1],resolution[2],3,3], 'd') + for x in range(0,resolution[0]): + for y in range(0,resolution[1]): + for z in range(0,resolution[2]): + currentPosition = currentPosition + distance*8 + p.file.seek(currentPosition) + for i in range(0,3): + for j in range(0,3): + field[x][y][z][i][j]=struct.unpack('d',p.file.read(8))[0] + return field + +#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +def mesh(res,geomdim,defgrad_av,centroids): +#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + neighbor = numpy.array([[0, 0, 0], + [1, 0, 0], + [1, 1, 0], + [0, 1, 0], + [0, 0, 1], + [1, 0, 1], + [1, 1, 1], + [0, 1, 1]]) + + wrappedCentroids = numpy.zeros([res[0]+2,res[1]+2,res[2]+2,3],'d') + nodes = numpy.zeros([res[0]+1,res[1]+1,res[2]+1,3],'d') + wrappedCentroids[1:-1,1:-1,1:-1] = centroids + diag = numpy.ones(3,'i') + shift = numpy.zeros(3,'i') + lookup = numpy.zeros(3,'i') + + for k in range(res[2]+2): + for j in range(res[1]+2): + for i in range(res[0]+2): + if (k==0 or k==res[2]+1 or \ + j==0 or j==res[1]+1 or \ + i==0 or i==res[0]+1 ): + me = numpy.array([i,j,k],'i') + shift = numpy.sign(res+diag-2*me)*(numpy.abs(res+diag-2*me)/(res+diag)) + lookup = me-diag+shift*res + wrappedCentroids[i,j,k] = centroids[lookup[0],lookup[1],lookup[2]]- \ + numpy.dot(defgrad_av, shift*geomdim) + for k in range(res[2]+1): + for j in range(res[1]+1): + for i in range(res[0]+1): + for n in range(8): + nodes[i,j,k] += wrappedCentroids[i+neighbor[n,0],j+neighbor[n,1],k+neighbor[n,2]] + nodes[:,:,:] /=8.0 + + return nodes + +# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +def centroids(res,geomdimension,defgrad): +# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + corner = numpy.array([[0, 0, 0], + [1, 0, 0], + [1, 1, 0], + [0, 1, 0], + [1, 1, 1], + [0, 1, 1], + [0, 0, 1], + [1, 0, 1]]) + step = numpy.array([[ 1, 1, 1], + [-1, 1, 1], + [-1,-1, 1], + [ 1,-1, 1], + [-1,-1,-1], + [ 1,-1,-1], + [ 1, 1,-1], + [-1, 1,-1]]) + + order = numpy.array([[0, 1, 2], + [0, 2, 1], + [1, 0, 2], + [1, 2, 0], + [2, 0, 1], + [2, 1, 0]]) + + cornerCoords = numpy.zeros([8,res[0],res[1],res[2],3], 'd') + coord = numpy.zeros([8,6,res[0],res[1],res[2],3], 'd') + centroids = numpy.zeros([res[0],res[1],res[2],3], 'd') + myStep = numpy.zeros(3,'d') + rear = numpy.zeros(3, 'i') + init = numpy.zeros(3, 'i') + oppo = numpy.zeros(3, 'i') + me = numpy.zeros(3, 'i') + ones = numpy.ones(3, 'i') + fones = numpy.ones(3, 'd') + defgrad_av=numpy.average(numpy.average(numpy.average(defgrad,0),0),0) + + for s in range(8):# corners + init = corner[s]*(res-ones) + oppo = corner[(s+4)%8]*(res-ones) + for o in range(6): # orders + for k in range(init[order[o,2]],oppo[order[o,2]]+step[s,order[o,2]],step[s,order[o,2]]): + rear[order[o,1]] = init[order[o,1]] + for j in range(init[order[o,1]],oppo[order[o,1]]+step[s,order[o,1]],step[s,order[o,1]]): + rear[order[o,0]] = init[order[o,0]] + for i in range(init[order[o,0]],oppo[order[o,0]]+step[s,order[o,0]],step[s,order[o,0]]): + me[order[o,0]] = i + me[order[o,1]] = j + me[order[o,2]] = k + if (numpy.all(me == init)): + coord[s,o,me[0],me[1],me[2]] = geomdimension * (numpy.dot(defgrad_av,corner[s]) + \ + numpy.dot(defgrad[me[0],me[1],me[2]],0.5*step[s]/res)) + else: + myStep = (me-rear)*geomdimension/res + coord[s,o,me[0],me[1],me[2]] = coord[s,o,rear[0],rear[1],rear[2]]+ \ + 0.5*numpy.dot(defgrad[me[0],me[1],me[2]] + \ + defgrad[rear[0],rear[1],rear[2]],myStep) + + rear[:] = me[:] + cornerCoords[s] = numpy.average(coord[s],0) + for k in range(res[2]): + for j in range(res[1]): + for i in range(res[0]): + parameter_coords=(2.0*numpy.array([i,j,k])-res+fones)/(res-fones) + pos = (fones + parameter_coords) + neg = (fones - parameter_coords) + if(k<3 and j<3 and i<3): + print i,j,k,':',pos, neg,'(',parameter_coords,')' + centroids[i,j,k] = ( cornerCoords[0,i,j,k] *neg[0]*neg[1]*neg[2]\ + + cornerCoords[1,i,j,k] *pos[0]*neg[1]*neg[2]\ + + cornerCoords[2,i,j,k] *pos[0]*pos[1]*neg[2]\ + + cornerCoords[3,i,j,k] *neg[0]*pos[1]*neg[2]\ + + cornerCoords[4,i,j,k] *pos[0]*pos[1]*pos[2]\ + + cornerCoords[5,i,j,k] *neg[0]*pos[1]*pos[2]\ + + cornerCoords[6,i,j,k] *neg[0]*neg[1]*pos[2]\ + + cornerCoords[7,i,j,k] *pos[0]*neg[1]*pos[2])*0.125 + return centroids, defgrad_av + +# function writes scalar values to a mesh (geometry) +def writeVtkAscii(filename,geometry,scalar,resolution): + prodnn=(p.resolution[0]+1)*(p.resolution[1]+1)*(p.resolution[2]+1) + vtk = open(filename, 'w') + vtk.write('# vtk DataFile Version 3.1\n') # header + vtk.write('just a test\n') # header + vtk.write('ASCII\n') # header + vtk.write('DATASET UNSTRUCTURED_GRID\n') # header + vtk.write('POINTS ') # header + vtk.write(str(prodnn)) # header + vtk.write(' FLOAT\n') # header + +# nodes + for k in range (resolution[2]+1): + for j in range (resolution[1]+1): + for i in range (resolution[0]+1): + vtk.write('\t'.join(map(str,geometry[i,j,k]))+'\n') + vtk.write('\n') + vtk.write('CELLS ') + vtk.write(str(resolution[0]*resolution[1]*resolution[2])) + vtk.write('\t') + vtk.write(str(resolution[0]*resolution[1]*resolution[2]*9)) + vtk.write('\n') + +# elements + for i in range (resolution[2]): + for j in range (resolution[1]): + for k in range (resolution[0]): + vtk.write('8') + vtk.write('\t') + base = i*(resolution[1]+1)*(resolution[2]+1)+j*(resolution[1]+1)+k + vtk.write(str(base)) + vtk.write('\t') + vtk.write(str(base+1)) + vtk.write('\t') + vtk.write(str(base+resolution[1]+2)) + vtk.write('\t') + vtk.write(str(base+resolution[1]+1)) + vtk.write('\t') + base = base + (resolution[1]+1)*(resolution[2]+1) + vtk.write(str(base)) + vtk.write('\t') + vtk.write(str(base+1)) + vtk.write('\t') + vtk.write(str(base+resolution[1]+2)) + vtk.write('\t') + vtk.write(str(base+resolution[1]+1)) + vtk.write('\n') + vtk.write('\n') + vtk.write('CELL_TYPES ') + vtk.write('\t') + vtk.write(str(resolution[0]*resolution[1]*resolution[2])) + vtk.write('\n') + for i in range (resolution[0]*resolution[1]*resolution[2]): + vtk.write('12\n') + vtk.write('\nCELL_DATA ') # header + vtk.write(str(resolution[0]*resolution[1]*resolution[2])) # header + vtk.write('\n') # header + vtk.write('SCALARS HorizontalSpeed float\n') # header + vtk.write('LOOKUP_TABLE default\n') # header + for k in range (resolution[2]): + for j in range (resolution[1]): + for i in range (resolution[0]): + vtk.write(str(scalar[i,j,k])) + vtk.write('\n') + return + +# function writes scalar values to a point field +def writeVtkAsciiDots(filename,coordinates,scalar,resolution): + prodnn=(p.resolution[0])*(p.resolution[1])*(p.resolution[2]) + vtk = open(filename, 'w') + vtk.write('# vtk DataFile Version 3.1\n') # header + vtk.write('just a test\n') # header + vtk.write('ASCII\n') # header + vtk.write('DATASET UNSTRUCTURED_GRID\n') # header + vtk.write('POINTS ') # header + vtk.write(str(prodnn)) # header + vtk.write(' FLOAT\n') # header +# points + for k in range (resolution[2]): + for j in range (resolution[1]): + for i in range (resolution[0]): + vtk.write('\t'.join(map(str,coordinates[i,j,k]))+'\n') + vtk.write('\n') + vtk.write('CELLS ') + vtk.write(str(prodnn)) + vtk.write('\t') + vtk.write(str(prodnn*2)) + vtk.write('\n') + for i in range(prodnn): + vtk.write('1\t' + str(i) + '\n') + vtk.write('CELL_TYPES ') + vtk.write('\t') + vtk.write(str(prodnn)) + vtk.write('\n') + for i in range (prodnn): + vtk.write('1\n') + vtk.write('\nPOINT_DATA ') # header + vtk.write(str(prodnn)) # header + vtk.write('\n') # header + vtk.write('SCALARS HorizontalSpeed float\n') # header + vtk.write('LOOKUP_TABLE default\n') # header + for k in range (resolution[2]): + for j in range (resolution[1]): + for i in range (resolution[0]): + vtk.write(str(scalar[i,j,k])) + vtk.write('\n') + return + +# functiongives the corner box for the average defgrad +def writeVtkAsciidefgrad_av(filename,diag,defgrad): + + points = numpy.array([\ + [0.0,0.0,0.0,],\ + [diag[0],0.0,0.0,],\ + [diag[0],diag[1],0.0,],\ + [0.0,diag[1],0.0,],\ + [0.0,0.0,diag[2],],\ + [diag[0],0.0,diag[2],],\ + [diag[0],diag[1],diag[2],],\ + [0.0,diag[1],diag[2],]]\ + ) + vtk = open(filename, 'w') + vtk.write('# vtk DataFile Version 3.1\n') # header + vtk.write('just a test\n') # header + vtk.write('ASCII\n') # header + vtk.write('DATASET UNSTRUCTURED_GRID\n') # header + vtk.write('POINTS 8') # header + vtk.write(' FLOAT\n') # header + + # points + for p in range (8): + vtk.write('\t'.join(map(str,numpy.dot(defgrad_av,points[p])))+'\n') + vtk.write('\n') + vtk.write('CELLS 8 16\n') + vtk.write('\n'.join(['1\t%i'%i for i in range(8)])+'\n') + vtk.write('CELL_TYPES 8\n') + vtk.write('\n'.join(['1']*8)+'\n') + + return + +print '*********************************************************************************' print 'Post Processing for Material subroutine for BVP solution using spectral method' print '*********************************************************************************\n' #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 +p = MPIEspectral_result('32x32x32x100.spectralOut') +p.extrapolation('') +print p # 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') +res_x=p.resolution[0] +res_y=p.resolution[1] +res_z=p.resolution[2] -# some funtions that might be useful -def InCharCount(location, character): - subj = file(location, "r") - body = subj.read() - subj.close() - return body.count(character) - - -def InCharCount(location, character): - subj = file(location, "r") - - nbr_of_char = 0 - for line in subj: - nbr_of_char = nbr_of_char + line.count(character) - - return nbr_of_char + +for i in range(120,121): + c_pos = p.dataOffset + i*(p.N_element_scalars*8*p.N_elements + 8) #8 accounts for header&footer + defgrad = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,16) + grain = readScalar(p.resolution,p.file,p.N_element_scalars,c_pos,7) + centroids, defgrad_av = centroids(p.resolution,p.dimension,defgrad) + print centroids.shape, defgrad_av.shape + ms = mesh(p.resolution,p.dimension,defgrad_av,centroids) + writeVtkAscii('mesh-%i.vtk'%i,ms,grain,p.resolution) + writeVtkAsciidefgrad_av('box-%i.vtk'%i,p.dimension,defgrad_av) + writeVtkAsciiDots('points-%i.vtk'%i,centroids,grain,p.resolution) + sys.stdout.flush()