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
This commit is contained in:
Martin Diehl
parent
ccc6aac10b
commit
322cbd2597
|
|
@ -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
|
||||||
|
|
@ -1,14 +1,12 @@
|
||||||
! -*- f90 -*-
|
! -*- f90 -*-
|
||||||
subroutine simple(defgrad,res_x,res_y,res_z,geomdimension,current_configuration)
|
subroutine simple(defgrad,res_x,res_y,res_z,geomdimension,current_configuration)
|
||||||
implicit none
|
implicit none
|
||||||
! *** Precision of real and integer variables ***
|
integer, parameter :: pDouble = selected_real_kind(15,50)
|
||||||
integer, parameter :: pReal = selected_real_kind(15) ! 15 significant digits, up to 1e+-300
|
integer i,j,k
|
||||||
integer, parameter :: pInt = selected_int_kind(9) ! up to +- 1e9
|
integer res_x, res_y, res_z
|
||||||
integer(pInt) i,j,k
|
integer, dimension(3) :: resolution
|
||||||
integer(pInt) res_x, res_y, res_z
|
real*8, dimension(3) :: geomdimension
|
||||||
integer(pInt), dimension(3) :: resolution
|
real*8, dimension(3,3) :: temp33_Real
|
||||||
real(pReal) , dimension(3) :: geomdimension
|
|
||||||
real(pReal) , dimension(3,3) :: temp33_Real
|
|
||||||
real*8 current_configuration(res_x, res_y, res_z,3)
|
real*8 current_configuration(res_x, res_y, res_z,3)
|
||||||
real*8 defgrad(res_x, res_y, res_z,3,3)
|
real*8 defgrad(res_x, res_y, res_z,3,3)
|
||||||
!f2py intent(in) res_x, res_y, res_z
|
!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
|
!f2py depend(res_x, res_y, res_z) defgrad
|
||||||
|
|
||||||
resolution(1) = res_x; resolution(2) = res_y; resolution(3) = res_z
|
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)
|
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
|
if((k==1).and.(j==1).and.(i==1)) then
|
||||||
temp33_Real =0.0_pReal
|
temp33_Real =0.0_pDouble
|
||||||
else
|
else
|
||||||
if((j==1).and.(i==1)) then
|
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_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(2,:) = temp33_Real(1,:)
|
||||||
temp33_Real(3,:) = temp33_Real(1,:)
|
temp33_Real(3,:) = temp33_Real(1,:)
|
||||||
current_configuration(i,j,k,:) = temp33_Real(1,:)
|
current_configuration(i,j,k,:) = temp33_Real(1,:)
|
||||||
else
|
else
|
||||||
if(i==1) then
|
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_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,:)
|
temp33_Real(3,:) = temp33_Real(2,:)
|
||||||
current_configuration(i,j,k,:) = temp33_Real(2,:)
|
current_configuration(i,j,k,:) = temp33_Real(2,:)
|
||||||
else
|
else
|
||||||
temp33_Real(3,:) = temp33_Real(3,:) + matmul(defgrad(i,j,k,:,:),&
|
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,:)
|
current_configuration(i,j,k,:) = temp33_Real(3,:)
|
||||||
endif
|
endif
|
||||||
endif
|
endif
|
||||||
endif
|
endif
|
||||||
|
!print*, current_configuration
|
||||||
enddo; enddo; enddo
|
enddo; enddo; enddo
|
||||||
end subroutine simple
|
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)
|
||||||
|
|
|
||||||
|
|
@ -1,100 +1,115 @@
|
||||||
#!/usr/bin/python
|
#!/usr/bin/python
|
||||||
# -*- coding: iso-8859-1 -*-
|
# -*- 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 array
|
||||||
import struct
|
import struct
|
||||||
#import numpy
|
import numpy
|
||||||
print('post processing for mpie_spectral')
|
import reconstruct
|
||||||
|
|
||||||
filename ='results.out'
|
#this funtion finds a string value in the given file for a given keyword, also returns the position
|
||||||
results = open(filename, 'rb')
|
def searchNameForKeyword(searchstring, file):
|
||||||
print('filename:', filename)
|
file.seek(0,0)
|
||||||
position=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
|
#this funtion finds an integer value in the given file for a given keyword, works similar as "searchNameForKeyword"
|
||||||
def searchNameForKeyword(searchstring, file, maxpos):
|
def searchValueForKeyword(searchstring, file):
|
||||||
results.seek(0,0)
|
results.seek(0,0)
|
||||||
begin = file.read(2048).find(searchstring) + len(searchstring) # function will return -1 if string isnt found
|
begin = file.read(2048).find(searchstring) + len(searchstring)
|
||||||
file.seek(begin -len(searchstring) -4) #position of header
|
|
||||||
header = file.read(4) #store header
|
|
||||||
file.seek(begin,0)
|
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])
|
|
||||||
value = struct.unpack('i',results.read(4))[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
|
#finds the value for the three integers (a,b,c) following the given keyword
|
||||||
def searchArrayForKeyword(searchstring, file, maxpos):
|
def searchIntegerArrayForKeyword(searchstring, file):
|
||||||
values = array.array('i',[0,0,0])
|
values = array.array('i',[0,0,0])
|
||||||
results.seek(0,0)
|
file.seek(0,0)
|
||||||
begin = results.read(2048).find(searchstring) + len(searchstring) #end position of string "resolution"
|
begin = file.read(2048).find(searchstring) + len(searchstring) #end position of string "searchstring"
|
||||||
pos = results.read(60).find(b'a')
|
pos = file.read(60).find('a')
|
||||||
results.seek(begin+pos+2,0) #why two, not 1??
|
file.seek(begin+pos+2,0) #why 2, not 1??
|
||||||
values[0]=struct.unpack('i',results.read(4))[0]
|
values[0]=struct.unpack('i',file.read(4))[0]
|
||||||
maxpos=max(maxpos,results.tell())
|
maxpos=file.tell()
|
||||||
results.seek(begin,0)
|
file.seek(begin,0)
|
||||||
pos = results.read(60).find(b'b')
|
pos = file.read(60).find('b')
|
||||||
results.seek(begin+pos+1,0)
|
file.seek(begin+pos+1,0)
|
||||||
values[1]=struct.unpack('i',results.read(4))[0]
|
values[1]=struct.unpack('i',file.read(4))[0]
|
||||||
maxpos=max(maxpos,results.tell())
|
maxpos = max(maxpos,file.tell())
|
||||||
results.seek(begin,0)
|
file.seek(begin,0)
|
||||||
pos = results.read(60).find(b'c')
|
pos = file.read(60).find('c')
|
||||||
results.seek(begin+pos+1,0)
|
file.seek(begin+pos+1,0)
|
||||||
values[2]=struct.unpack('i',results.read(4))[0]
|
values[2]=struct.unpack('i',file.read(4))[0]
|
||||||
maxpos=max(maxpos,results.tell())
|
maxpos=max(maxpos,file.tell())
|
||||||
return values, maxpos
|
return values, maxpos
|
||||||
|
|
||||||
#finds the value for the three doubles followed the given keyword
|
#finds the value for the three doubles (x,y,z) following the given keyword
|
||||||
def searchFarrayForKeyword(searchstring, file, maxpos):
|
def searchDoubleArrayForKeyword(searchstring, file):
|
||||||
values = array.array('d',[0,0,0])
|
values = array.array('d',[0,0,0])
|
||||||
results.seek(0,0)
|
file.seek(0,0)
|
||||||
begin = results.read(2048).find(searchstring) + len(searchstring) #end position of string "resolution"
|
begin = file.read(2048).find(searchstring) + len(searchstring)
|
||||||
pos = results.read(60).find(b'x')
|
pos = file.read(60).find('x')
|
||||||
results.seek(begin+pos+2,0) #why two, not 1??
|
file.seek(begin+pos+2,0)
|
||||||
values[0]=struct.unpack('d',results.read(8))[0]
|
values[0]=struct.unpack('d',file.read(8))[0]
|
||||||
maxpos=max(maxpos,results.tell())
|
maxpos=file.tell()
|
||||||
results.seek(begin,0)
|
file.seek(begin,0)
|
||||||
pos = results.read(60).find(b'y')
|
pos = file.read(60).find('y')
|
||||||
results.seek(begin+pos+1,0)
|
file.seek(begin+pos+1,0)
|
||||||
values[1]=struct.unpack('d',results.read(8))[0]
|
values[1] = struct.unpack('d',file.read(8))[0]
|
||||||
maxpos=max(maxpos,results.tell())
|
maxpos = max(maxpos,file.tell())
|
||||||
results.seek(begin,0)
|
file.seek(begin,0)
|
||||||
pos = results.read(60).find(b'z')
|
pos = file.read(60).find('z')
|
||||||
results.seek(begin+pos+1,0)
|
file.seek(begin+pos+1,0)
|
||||||
values[2]=struct.unpack('d',results.read(8))[0]
|
values[2] = struct.unpack('d',file.read(8))[0]
|
||||||
maxpos=max(maxpos,results.tell())
|
maxpos = max(maxpos,file.tell())
|
||||||
return values, maxpos
|
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)
|
#reading in the header of the results file
|
||||||
workingdir, position = searchNameForKeyword(b'Workingdir ', results, position)
|
filename ='32x32x32x100.out'
|
||||||
jobname, position = searchNameForKeyword(b'JobName ', results, position)
|
print 'Results Filename:', filename
|
||||||
totalincs, position = searchValueForKeyword(b'totalincs ', results, position)
|
results = open(filename, 'rb')
|
||||||
materialpoint_sizeResults, position = searchValueForKeyword(b'materialpoint_sizeResults ', results, position)
|
|
||||||
resolution, position = searchArrayForKeyword(b'resolution ', results, position)
|
loadcase, position = searchNameForKeyword('Loadcase', results)
|
||||||
geomdimension, position = searchFarrayForKeyword(b'geomdimension ', results, position)
|
workingdir, temp = searchNameForKeyword('Workingdir', results)
|
||||||
position=position+4 +13*8
|
position = max(temp, position)
|
||||||
results.seek(position,0)
|
jobname, temp = searchNameForKeyword('JobName', results)
|
||||||
tnsr = array.array('d',[0,0,0,0,0,0,0,0,0])
|
position = max(temp, position)
|
||||||
tnsr[0]=struct.unpack('d',results.read(8))[0]
|
totalincs, temp = searchValueForKeyword('totalincs', results)
|
||||||
tnsr[1]=struct.unpack('d',results.read(8))[0]
|
position = max(temp, position)
|
||||||
tnsr[2]=struct.unpack('d',results.read(8))[0]
|
materialpoint_sizeResults, temp = searchValueForKeyword('materialpoint_sizeResults', results)
|
||||||
tnsr[3]=struct.unpack('d',results.read(8))[0]
|
position = max(temp, position)
|
||||||
tnsr[4]=struct.unpack('d',results.read(8))[0]
|
resolution, temp = searchIntegerArrayForKeyword('resolution', results)
|
||||||
print(tnsr)
|
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
|
# some funtions that might be useful
|
||||||
#filename = jobname[0:len(jobname)-5]+b'.outputCrystallite'
|
|
||||||
#outputCrystallite = open(filename, 'r')
|
|
||||||
|
|
||||||
def InCharCount(location, character):
|
def InCharCount(location, character):
|
||||||
subj = file(location, "r")
|
subj = file(location, "r")
|
||||||
body = subj.read()
|
body = subj.read()
|
||||||
|
|
|
||||||
|
|
@ -18,7 +18,7 @@ program voronoi
|
||||||
print*, '******************************************************************************'
|
print*, '******************************************************************************'
|
||||||
print*, ''
|
print*, ''
|
||||||
print*, 'generates:'
|
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*, ' * material file "material.config": Orientation information for solver'
|
||||||
print*, ' * "_GIVEN_NAME_.spectral": combined information for solver'
|
print*, ' * "_GIVEN_NAME_.spectral": combined information for solver'
|
||||||
print*, 'optional output:'
|
print*, 'optional output:'
|
||||||
|
|
@ -37,7 +37,7 @@ program voronoi
|
||||||
read(*, *), c
|
read(*, *), c
|
||||||
write(*, '(A)', advance = 'NO') 'Please enter No. of Grains: '
|
write(*, '(A)', advance = 'NO') 'Please enter No. of Grains: '
|
||||||
read(*, *), N_Seeds
|
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
|
read(*, *), name
|
||||||
write(*, '(A)', advance = 'NO') 'Should the visualization files be generated (y/n)? '
|
write(*, '(A)', advance = 'NO') 'Should the visualization files be generated (y/n)? '
|
||||||
read(*, *), choice
|
read(*, *), choice
|
||||||
|
|
@ -77,8 +77,8 @@ program voronoi
|
||||||
print*, ''
|
print*, ''
|
||||||
print*, 'material config file is written out'
|
print*, 'material config file is written out'
|
||||||
|
|
||||||
!write header of mesh file, should be done before the following change of variables
|
!write header of geom file, should be done before the following change of variables
|
||||||
open(20, file = ((trim(name))//'.mesh'))
|
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, 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, I4, A, I4, A, I4)'), 'dimension x ', 2**a, ' y ', 2**b, ' z ', 2**c
|
||||||
write(20, '(A)'), 'homogenization 1'
|
write(20, '(A)'), 'homogenization 1'
|
||||||
|
|
@ -166,13 +166,13 @@ program voronoi
|
||||||
&mod((i-1), a)+1, mod(((i-1)/a), b)+1, mod(((i-1)/(ab)), c)+1, &
|
&mod((i-1), a)+1, mod(((i-1)/a), b)+1, mod(((i-1)/(ab)), c)+1, &
|
||||||
&j, ' 1'
|
&j, ' 1'
|
||||||
end do
|
end do
|
||||||
print*, 'mesh files are written out'
|
print*, 'geometry files are written out'
|
||||||
|
|
||||||
!write visualization files (in case wanted)
|
!write visualization files (in case wanted)
|
||||||
if (choice == 'y' .or. choice == 'Y') then
|
if (choice == 'y' .or. choice == 'Y') then
|
||||||
print*, 'for more information on gmsh: http://geuz.org/gmsh/'
|
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'))
|
open(20, file = ((trim(name))//'_3Dfull.msh'))
|
||||||
write(20, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', abc
|
write(20, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', abc
|
||||||
do i = 1, abc
|
do i = 1, abc
|
||||||
|
|
@ -213,7 +213,7 @@ program voronoi
|
||||||
write(20, *), '$EndNodeData'
|
write(20, *), '$EndNodeData'
|
||||||
close(20)
|
close(20)
|
||||||
|
|
||||||
! write 2d mesh out
|
! write 2d geometry out
|
||||||
open(20, file = ((trim(name))//'_2D.msh'))
|
open(20, file = ((trim(name))//'_2D.msh'))
|
||||||
write(20, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', ab
|
write(20, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', ab
|
||||||
do i = 1, ab
|
do i = 1, ab
|
||||||
|
|
@ -18,7 +18,7 @@ program voronoi
|
||||||
print*, '******************************************************************************'
|
print*, '******************************************************************************'
|
||||||
print*, ''
|
print*, ''
|
||||||
print*, 'generates:'
|
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*, ' * material file "material.config": Orientation information for solver'
|
||||||
print*, 'optional output:'
|
print*, 'optional output:'
|
||||||
print*, ' * view file "_GIVEN_NAME_2D.msh": Information for visualization in gmsh'
|
print*, ' * view file "_GIVEN_NAME_2D.msh": Information for visualization in gmsh'
|
||||||
|
|
@ -36,7 +36,7 @@ program voronoi
|
||||||
read(*, *), c
|
read(*, *), c
|
||||||
write(*, '(A)', advance = 'NO') 'Please enter No. of Grains: '
|
write(*, '(A)', advance = 'NO') 'Please enter No. of Grains: '
|
||||||
read(*, *), N_Seeds
|
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
|
read(*, *), name
|
||||||
write(*, '(A)', advance = 'NO') 'Should the visualization files be generated (y/n)? '
|
write(*, '(A)', advance = 'NO') 'Should the visualization files be generated (y/n)? '
|
||||||
read(*, *), choice
|
read(*, *), choice
|
||||||
|
|
@ -75,8 +75,8 @@ program voronoi
|
||||||
print*, ''
|
print*, ''
|
||||||
print*, 'material config file is written out'
|
print*, 'material config file is written out'
|
||||||
|
|
||||||
!write header of mesh file, should be done before the following change of variables
|
!write header of geometry file, should be done before the following change of variables
|
||||||
open(20, file = ((trim(name))//'.mesh'))
|
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, 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, I4, A, I4, A, I4)'), 'dimension x ', 2**a, ' y ', 2**b, ' z ', 2**c
|
||||||
write(20, '(A)'), 'homogenization 1'
|
write(20, '(A)'), 'homogenization 1'
|
||||||
|
|
@ -168,13 +168,13 @@ program voronoi
|
||||||
end do
|
end do
|
||||||
close(20)
|
close(20)
|
||||||
print*, 'voronoi tesselation finished'
|
print*, 'voronoi tesselation finished'
|
||||||
print*, 'mesh file is written out'
|
print*, 'geometry file is written out'
|
||||||
|
|
||||||
!write visualization files (in case wanted)
|
!write visualization files (in case wanted)
|
||||||
if (choice == 'y' .or. choice == 'Y') then
|
if (choice == 'y' .or. choice == 'Y') then
|
||||||
print*, 'for more information on gmsh: http://geuz.org/gmsh/'
|
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'))
|
open(20, file = ((trim(name))//'_3D.msh'))
|
||||||
write(20, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', abc_Red
|
write(20, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', abc_Red
|
||||||
do j = 1, abc_Red
|
do j = 1, abc_Red
|
||||||
|
|
@ -194,7 +194,7 @@ program voronoi
|
||||||
write(20, *), '$EndNodeData'
|
write(20, *), '$EndNodeData'
|
||||||
close(20)
|
close(20)
|
||||||
|
|
||||||
! write 2d mesh out
|
! write 2d geometry out
|
||||||
open(20, file = ((trim(name))//'_2D.msh'))
|
open(20, file = ((trim(name))//'_2D.msh'))
|
||||||
write(20, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', ab
|
write(20, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', ab
|
||||||
do i = 1, ab
|
do i = 1, ab
|
||||||
Loading…
Reference in New Issue