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

@@ -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 i,j,k, l,m, s,o, loop, res_x, res_y, res_z
- integer 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, dimension(3,res_x, res_y, res_z) :: coordinates2
- real*8 defgrad(res_x, res_y, res_z,3,3)
+ real*8, dimension(3,3,res_x, res_y, res_z) :: defgrad
- !f2py intent(in) 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
+ !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,:,:),&
+       temp33_Real(1,:) = temp33_Real(1,:) + matmul(defgrad(:,:,i,j,k),&
-                  (/0.0_pDouble,0.0_pDouble,geomdimension(3)/(real(resolution(3)))/))
+                  (/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,:,:),&
+         temp33_Real(2,:) = temp33_Real(2,:) + matmul(defgrad(:,:,i,j,k),&
-                   (/0.0_pDouble,geomdimension(2)/(real(resolution(2))),0.0_pDouble/))
+                   (/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,:,:),&
+         temp33_Real(3,:) = temp33_Real(3,:) + matmul(defgrad(:,:,i,j,k),&
-                  (/geomdimension(1)/(real(resolution(1))),0.0_pDouble,0.0_pDouble/))
+                  (/real(geomdimension(1))/(real(resolution(1))),real(0.0),real(0.0)/))
-         current_configuration(i,j,k,:) = temp33_Real(3,:)   
+         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 i,j,k, l,m, s,o, loop, res_x, res_y, res_z
- integer a,b,c
+ integer, dimension(3) :: res, init, oppo, me, rear
- 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, dimension(3) :: myStep
- real*8 current_configuration(res_x, res_y, res_z,3)
+ real*8, dimension(3,3) :: defGrad_av
- real*8 defgrad(res_x, res_y, res_z,3,3)
+ integer, dimension(3,8) :: corner = reshape((/ &
-! some declarations for the wrapping to python 
+                                     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) coordinates
- !f2py depend(res_x, res_y, res_z) defgrad
+ !f2py depend(res_x, res_y, res_z) defgrad 
- do i=1, 3; do j=1,3 !
+ !f2py depend(res_x, res_y, res_z) cornerCoords
-   defgrad_av(i,j) = sum(defgrad(:,:,:,i,j)) /real(res_x*res_y*res_z)
+ !f2py depend(res_x, res_y, res_z) coord
- 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)
+ 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
+ integer, dimension(3) :: res, shift, lookup, me
- real*8, dimension(3) :: gdim
+ real*8, dimension(3) :: geomdim, diag = (/1,1,1/)
- real*8 inter(res_x, res_y, res_z,3)
+ real*8, dimension(3,3) :: defgrad_av
- real*8 configuration(res_x+2, res_y+2, res_z+2,3)
+ real*8, dimension(3,res_x, res_y, res_z) :: centroids
- real*8 meshgeom(res_x+1, res_y+1, res_z+1,3)
+ 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(in) centroids
- !f2py intent(out) meshgeom
+ !f2py intent(in) defgrad_av
- !f2py intent(in) gdim
+ !f2py intent(in) geomdim
- !f2py depend(res_x, res_y, res_z) inter
+ !f2py intent(out) nodes
- !f2py depend(res_x, res_y, res_z) meshgeom
+ !f2py depend(res_x, res_y, res_z) centroids
- meshgeom = 0.0_pDouble
+ !f2py depend(res_x, res_y, res_z) nodes
- configuration = 0.0_pDouble
+
- 
+ res = (/res_x,res_y,res_z/)
- resolution(1) = res_x; resolution(2) = res_y; resolution(3) = res_z
+ wrappedCentroids(:,2:res(1)+1,2:res(2)+1,2:res(3)+1) = centroids
- 
+ 
- 
+ do k = 0,res(3)+1
- do k = 1, resolution(3)+1; do j = 1, resolution(2)+1; do i = 1, resolution(1)+1
+   do j = 0,res(2)+1
-   if((i==1).and.(j==1).and.(k==1)) then
+     do i = 0,res(1)+1
-  configuration(i,j,k,:)=inter(res_x,res_y,res_z,:)-gdim
+       if (&
-  configuration(res_x+2,j,k,:)=inter(i,res_y,res_z,:)+(/1.0,-1.0,-1.0/)*gdim
+         k==0 .or. k==res(3)+1 .or. &
-  configuration(i,res_y+2,k,:)=inter(res_x,j,res_z,:)+(/-1.0,1.0,-1.0/)*gdim
+         j==0 .or. j==res(2)+1 .or. &
-  configuration(i,j,res_z+2,:)=inter(res_x,res_y,k,:)+(/-1.0,-1.0,1.0/)*gdim
+         i==0 .or. i==res(1)+1      &
-  configuration(res_x+2,res_y+2,k,:)=inter(i,j,res_z,:)+(/1.0,1.0,-1.0/)*gdim
+         ) then
-  configuration(i,res_y+2,res_z+2,:)=inter(res_x,j,k,:)+(/-1.0,1.0,1.0/)*gdim
+         me = (/i,j,k/)
-  configuration(res_x+2,j,res_z+2,:)=inter(i,res_y,k,:)+(/1.0,-1.0,1.0/)*gdim
+         shift = (res+diag-2*me)/(res+diag)
-  configuration(res_x+2,res_y+2,res_z+2,:)=inter(i,j,k,:)+gdim
+         lookup = me+shift*res
-   endif       
+         wrappedCentroids(:,1+i,1+j,1+k) = centroids(:,lookup(1),lookup(2),lookup(3)) - &
-   if((i==1).and.(j==1).and.(k/=1)) then
+                                           matmul(defgrad_av,shift * geomdim)
-  configuration(1,1,k,:)=inter(res_x,res_y,k-1,:)+(/-1.0,-1.0,0.0/)*gdim
+       endif
-  configuration(res_x+2,1,k,:)=inter(1,res_y,k-1,:)+(/1.0,-1.0,0.0/)*gdim
+     enddo
-  configuration(1,res_y+2,k,:)=inter(res_x,1,k-1,:)+(/-1.0,1.0,0.0/)*gdim 
+   enddo
-  configuration(res_x+2,res_y+2,k,:)=inter(1,1,k-1,:)+(/1.0,1.0,0.0/)*gdim
+ enddo
-   endif
+ 
-   if((i==1).and.(j/=1).and.(k==1)) then
+ nodes = 0.0_pDouble
-  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
+ do k = 0,res(3)
-  configuration(1,j,res_z+2,:)=inter(res_x,j-1,1,:)+(/-1.0,0.0,1.0/)*gdim
+   do j = 0,res(2)
-  configuration(res_x+2,j,res_z+2,:)=inter(1,j-1,1,:)+(/1.0,0.0,1.0/)*gdim
+     do i = 0,res(1)
-   endif
+       
-   if((i/=1).and.(j==1).and.(k==1)) then
+       do n = 1,8
-     configuration(i,1,1,:)=inter(i-1,res_y,res_z,:)+(/0.0,-1.0,-1.0/)*gdim
+         nodes(:,1+i,1+j,1+k) = nodes(:,1+i,1+j,1+k) + wrappedCentroids(:,1+i+neighbor(1,n),&
-     configuration(i,1,res_z+2,:)=inter(i-1,res_y,1,:)+(/0.0,-1.0,1.0/)*gdim
+                                                                          1+j+neighbor(2,n),&
-     configuration(i,res_y+2,1,:)=inter(i-1,1,res_z,:)+(/0.0,1.0,-1.0/)*gdim
+                                                                          1+k+neighbor(3,n))
-     configuration(i,res_y+2,res_z+2,:)=inter(i-1,1,1,:)+(/0.0,1.0,1.0/)*gdim
+       enddo
-   endif 
+       nodes(:,1+i,1+j,1+k) = nodes(:,1+i,1+j,1+k) / 8.0_pDouble
-   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
+     enddo 
-     configuration(i,j,res_z+2,:)=inter(i-1,j-1,1,:)+(/0.0,0.0,1.0/)*gdim
+   enddo 
-   endif
+ enddo 
-   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
 
 
@@ -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

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 os,sys,re,array,struct,numpy, time
-import struct
-import numpy
-import reconstruct
 
-#this funtion finds a string value in the given file for a given keyword, also returns the position
+class vector:
-def searchNameForKeyword(searchstring, file):
+  x,y,z = [None,None,None]
-    file.seek(0,0)
+  
-    begin = file.read(2048).find(searchstring) + len(searchstring) # function will return -1 if string is not found
+  def __init__(self,coords):
-    file.seek(begin -len(searchstring) -4)                         # position of header (Fortran specific, 4 bit long on msuws4
+    self.x = coords[0]
-    header = file.read(4)                                          # store header. header will be repeated, thus it allows us to determine the length of the searchstring
+    self.y = coords[1]
-    file.seek(begin,0)                                             # go back to starting postion
+    self.z = coords[2]
-    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 an integer value in the given file for a given keyword, works similar as "searchNameForKeyword"
+class element:
-def searchValueForKeyword(searchstring, file):
+  items = []
-    results.seek(0,0)
+  type = None
-    begin = file.read(2048).find(searchstring) + len(searchstring) 
-    file.seek(begin,0)                                                                     
-    value = struct.unpack('i',results.read(4))[0]
-    return value, results.tell()
 
-#finds the value for the three integers (a,b,c) following the given keyword
+  def __init__(self,nodes,type):
-def searchIntegerArrayForKeyword(searchstring, file):
+    self.items = nodes
-	values = array.array('i',[0,0,0])
+    self.type = type
-	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 (x,y,z) following the given keyword
+class element_scalar:
-def searchDoubleArrayForKeyword(searchstring, file):
+  id = None
-	values = array.array('d',[0,0,0])
+  value = None
-	file.seek(0,0)
+
-	begin = file.read(2048).find(searchstring) + len(searchstring) 
+  def __init__(self,node,value):
-	pos = file.read(60).find('x')
+    self.id = node
-	file.seek(begin+pos+2,0) 
+    self.value = value
-	values[0]=struct.unpack('d',file.read(8))[0]
+
-	maxpos=file.tell()
+
-	file.seek(begin,0)
+class MPIEspectral_result:
-	pos = file.read(60).find('y')
+
-	file.seek(begin+pos+1,0)
+  file = None
-	values[1] = struct.unpack('d',file.read(8))[0]
+  dataOffset = 0
-	maxpos = max(maxpos,file.tell())
+  N_elemental_scalars = 0
-	file.seek(begin,0)
+  resolution = [0,0,0]
-	pos = file.read(60).find('z')
+  dimension = [0.0,0.0,0.0]
-	file.seek(begin+pos+1,0)
+  theTitle = ''
-	values[2] = struct.unpack('d',file.read(8))[0]
+  wd = ''
-	maxpos = max(maxpos,file.tell())
+  extrapolate = ''
-	return values, maxpos
+  N_increments = 0
-print '*********************************************************************************'	
+  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'
+p = MPIEspectral_result('32x32x32x100.spectralOut')
-print 'Results Filename:', filename 
+p.extrapolation('')
-results = open(filename, 'rb')
+print p
-
-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'
+res_x=p.resolution[0]
-outputCrystallite = open(filename, 'r')
+res_y=p.resolution[1]
+res_z=p.resolution[2]
 
 
-# some funtions that might be useful
+
-def InCharCount(location, character):
+for i in range(120,121):
-    subj = file(location, "r")
+  c_pos = p.dataOffset + i*(p.N_element_scalars*8*p.N_elements + 8) #8 accounts for header&footer
-    body = subj.read()
+  defgrad = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,16)         
-    subj.close()
+  grain = readScalar(p.resolution,p.file,p.N_element_scalars,c_pos,7)
-    return body.count(character)
+  centroids, defgrad_av = centroids(p.resolution,p.dimension,defgrad)
-
+  print centroids.shape, defgrad_av.shape
-
+  ms = mesh(p.resolution,p.dimension,defgrad_av,centroids)
-def InCharCount(location, character):
+  writeVtkAscii('mesh-%i.vtk'%i,ms,grain,p.resolution)
-    subj = file(location, "r")
+  writeVtkAsciidefgrad_av('box-%i.vtk'%i,p.dimension,defgrad_av)
-
+  writeVtkAsciiDots('points-%i.vtk'%i,centroids,grain,p.resolution)
-    nbr_of_char = 0
+  sys.stdout.flush()
-    for line in subj:
-        nbr_of_char = nbr_of_char + line.count(character)
-
-    return nbr_of_char