223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

made linear shape reconstruction working again, pretty similar results like corrected FFT reconstruction

This commit is contained in:
Martin Diehl 2012-07-18 18:39:59 +00:00
parent b2a7f85101
commit 13caf2d389
2 changed files with 77 additions and 79 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

108
code/math.f90
View File

@ -16,16 +16,19 @@
! You should have received a copy of the GNU General Public License ! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>. ! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
! !
!############################################################## !--------------------------------------------------------------------------------------------------
!* $Id$ !* $Id$
!############################################################## !--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Mathematical library, including random number generation and tensor represenations
!--------------------------------------------------------------------------------------------------
#ifdef Spectral #ifdef Spectral
#include "kdtree2.f90" #include "kdtree2.f90"
#endif #endif
module math module math
!##############################################################
use, intrinsic :: iso_c_binding use, intrinsic :: iso_c_binding
use prec, only: pReal,pInt use prec, only: pReal,pInt
@ -35,15 +38,13 @@ module math
real(pReal), parameter, public :: INRAD = pi/180.0_pReal real(pReal), parameter, public :: INRAD = pi/180.0_pReal
complex(pReal), parameter, public :: TWOPIIMG = (0.0_pReal,2.0_pReal)* pi complex(pReal), parameter, public :: TWOPIIMG = (0.0_pReal,2.0_pReal)* pi
! *** 3x3 Identity ***
real(pReal), dimension(3,3), parameter, public :: & real(pReal), dimension(3,3), parameter, public :: &
math_I3 = reshape([& math_I3 = reshape([&
1.0_pReal,0.0_pReal,0.0_pReal, & 1.0_pReal,0.0_pReal,0.0_pReal, &
0.0_pReal,1.0_pReal,0.0_pReal, & 0.0_pReal,1.0_pReal,0.0_pReal, &
0.0_pReal,0.0_pReal,1.0_pReal & 0.0_pReal,0.0_pReal,1.0_pReal &
],[3,3]) ],[3,3]) !< 3x3 Identity
! *** Mandel notation ***
integer(pInt), dimension (2,6), parameter, private :: & integer(pInt), dimension (2,6), parameter, private :: &
mapMandel = reshape([& mapMandel = reshape([&
1_pInt,1_pInt, & 1_pInt,1_pInt, &
@ -52,7 +53,7 @@ module math
1_pInt,2_pInt, & 1_pInt,2_pInt, &
2_pInt,3_pInt, & 2_pInt,3_pInt, &
1_pInt,3_pInt & 1_pInt,3_pInt &
],[2,6]) ],[2,6]) !< Mandel notation
real(pReal), dimension(6), parameter, private :: & real(pReal), dimension(6), parameter, private :: &
nrmMandel = [& nrmMandel = [&
@ -64,7 +65,6 @@ module math
1.0_pReal, 1.0_pReal, 1.0_pReal,& 1.0_pReal, 1.0_pReal, 1.0_pReal,&
0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal] 0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal]
! *** Voigt notation ***
integer(pInt), dimension (2,6), parameter, private :: & integer(pInt), dimension (2,6), parameter, private :: &
mapVoigt = reshape([& mapVoigt = reshape([&
1_pInt,1_pInt, & 1_pInt,1_pInt, &
@ -73,7 +73,7 @@ module math
2_pInt,3_pInt, & 2_pInt,3_pInt, &
1_pInt,3_pInt, & 1_pInt,3_pInt, &
1_pInt,2_pInt & 1_pInt,2_pInt &
],[2,6]) ],[2,6]) !< Voigt notation
real(pReal), dimension(6), parameter, private :: & real(pReal), dimension(6), parameter, private :: &
nrmVoigt = 1.0_pReal, & nrmVoigt = 1.0_pReal, &
@ -154,9 +154,9 @@ real(pReal), dimension(4,36), parameter, private :: &
contains contains
!************************************************************************** !--------------------------------------------------------------------------------------------------
! initialization of module !> @brief initialization of random seed generator
!************************************************************************** !--------------------------------------------------------------------------------------------------
subroutine math_init subroutine math_init
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
@ -3172,15 +3172,11 @@ end subroutine mesh_regular_grid
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine deformed_linear(res,geomdim,defgrad_av,defgrad,coord_avgCorner) subroutine deformed_linear(res,geomdim,defgrad_av,defgrad,coord)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Routine to calculate coordinates in current configuration for given defgrad ! Routine to calculate coordinates in current configuration for given defgrad
! using linear interpolation (blurres out high frequency defomation) ! using linear interpolation (blurres out high frequency defomation)
! !
use debug, only: debug_math, &
debug_level, &
debug_levelBasic
implicit none implicit none
! input variables ! input variables
integer(pInt), intent(in), dimension(3) :: res integer(pInt), intent(in), dimension(3) :: res
@ -3188,15 +3184,14 @@ subroutine deformed_linear(res,geomdim,defgrad_av,defgrad,coord_avgCorner)
real(pReal), intent(in), dimension(3,3) :: defgrad_av real(pReal), intent(in), dimension(3,3) :: defgrad_av
real(pReal), intent(in), dimension( res(1),res(2),res(3),3,3) :: defgrad real(pReal), intent(in), dimension( res(1),res(2),res(3),3,3) :: defgrad
! output variables ! output variables
real(pReal), intent(out), dimension( res(1),res(2),res(3),3) :: coord_avgCorner real(pReal), intent(out), dimension( res(1),res(2),res(3),3) :: coord
! variables with dimension depending on input ! variables with dimension depending on input
real(pReal), dimension(8,6,res(1),res(2),res(3),3) :: coord
real(pReal), dimension( 8,res(1),res(2),res(3),3) :: coord_avgOrder real(pReal), dimension( 8,res(1),res(2),res(3),3) :: coord_avgOrder
! other variables ! other variables
real(pReal), dimension(3) :: myStep, fones = 1.0_pReal, parameter_coords, negative, positive real(pReal), dimension(3) :: myStep, fones = 1.0_pReal, parameter_coords, negative, positive, offset_coords
integer(pInt), dimension(3) :: rear, init, ones = 1_pInt, oppo, me integer(pInt), dimension(3) :: rear, init, ones = 1_pInt, oppo, me
integer(pInt) i, j, k, s, o integer(pInt) i, j, k, s, o
integer(pInt), dimension(3,8) :: corner = reshape((/ & integer(pInt), dimension(3,8) :: corner = reshape([ &
0_pInt, 0_pInt, 0_pInt,& 0_pInt, 0_pInt, 0_pInt,&
1_pInt, 0_pInt, 0_pInt,& 1_pInt, 0_pInt, 0_pInt,&
1_pInt, 1_pInt, 0_pInt,& 1_pInt, 1_pInt, 0_pInt,&
@ -3205,9 +3200,8 @@ subroutine deformed_linear(res,geomdim,defgrad_av,defgrad,coord_avgCorner)
0_pInt, 1_pInt, 1_pInt,& 0_pInt, 1_pInt, 1_pInt,&
0_pInt, 0_pInt, 1_pInt,& 0_pInt, 0_pInt, 1_pInt,&
1_pInt, 0_pInt, 1_pInt & 1_pInt, 0_pInt, 1_pInt &
/), & ],[3,8])
(/3,8/)) integer(pInt), dimension(3,8) :: step = reshape([&
integer(pInt), dimension(3,8) :: step = reshape((/ &
1_pInt, 1_pInt, 1_pInt,& 1_pInt, 1_pInt, 1_pInt,&
-1_pInt, 1_pInt, 1_pInt,& -1_pInt, 1_pInt, 1_pInt,&
-1_pInt,-1_pInt, 1_pInt,& -1_pInt,-1_pInt, 1_pInt,&
@ -3216,64 +3210,57 @@ subroutine deformed_linear(res,geomdim,defgrad_av,defgrad,coord_avgCorner)
1_pInt,-1_pInt,-1_pInt,& 1_pInt,-1_pInt,-1_pInt,&
1_pInt, 1_pInt,-1_pInt,& 1_pInt, 1_pInt,-1_pInt,&
-1_pInt, 1_pInt,-1_pInt & -1_pInt, 1_pInt,-1_pInt &
/), & ], [3,8])
(/3,8/)) integer(pInt), dimension(3,6) :: order = reshape([ &
integer(pInt), dimension(3,6) :: order = reshape((/ &
1_pInt, 2_pInt, 3_pInt,& 1_pInt, 2_pInt, 3_pInt,&
1_pInt, 3_pInt, 2_pInt,& 1_pInt, 3_pInt, 2_pInt,&
2_pInt, 1_pInt, 3_pInt,& 2_pInt, 1_pInt, 3_pInt,&
2_pInt, 3_pInt, 1_pInt,& 2_pInt, 3_pInt, 1_pInt,&
3_pInt, 1_pInt, 2_pInt,& 3_pInt, 1_pInt, 2_pInt,&
3_pInt, 2_pInt, 1_pInt & 3_pInt, 2_pInt, 1_pInt &
/), & ], [3,6])
(/3,6/))
if (iand(debug_level(debug_math),debug_levelBasic) /= 0_pInt) then
print*, 'Restore geometry using linear integration'
print '(a,3(e12.5))', ' Dimension: ', geomdim
print '(a,3(i5))', ' Resolution:', res
endif
coord_avgOrder = 0.0_pReal coord_avgOrder = 0.0_pReal
do s = 0_pInt, 7_pInt ! corners (from 0 to 7) do s = 0_pInt, 7_pInt ! corners (from 0 to 7)
init = corner(:,s+1_pInt)*(res-ones) +ones init = corner(:,s+1_pInt)*(res-ones) +ones
oppo = corner(:,mod((s+4_pInt),8_pInt)+1_pInt)*(res-ones) +ones oppo = corner(:,mod((s+4_pInt),8_pInt)+1_pInt)*(res-ones) +ones
do o=1_pInt,6_pInt ! orders (from 1 to 6) do o=1_pInt,6_pInt ! orders (from 1 to 6)
coord = 0_pReal
do k = init(order(3,o)), oppo(order(3,o)), step(order(3,o),s+1_pInt) do k = init(order(3,o)), oppo(order(3,o)), step(order(3,o),s+1_pInt)
rear(order(2,o)) = init(order(2,o)) rear(order(2,o)) = init(order(2,o))
do j = init(order(2,o)), oppo(order(2,o)), step(order(2,o),s+1_pInt) do j = init(order(2,o)), oppo(order(2,o)), step(order(2,o),s+1_pInt)
rear(order(1,o)) = init(order(1,o)) rear(order(1,o)) = init(order(1,o))
do i = init(order(1,o)), oppo(order(1,o)), step(order(1,o),s+1_pInt) do i = init(order(1,o)), oppo(order(1,o)), step(order(1,o),s+1_pInt)
me(order(1,o)) = i me(order(1:3,o)) = [i,j,k]
me(order(2,o)) = j if ( all(me==init)) then
me(order(3,o)) = k coord(me(1),me(2),me(3),1:3) = geomdim * (matmul(defgrad_av,real(corner(1:3,s+1),pReal)) + &
if ( (me(1)==init(1)).and.(me(2)==init(2)).and. (me(3)==init(3)) ) then
coord(s+1_pInt,o,me(1),me(2),me(3),1:3) = geomdim * (matmul(defgrad_av,real(corner(1:3,s+1),pReal)) + &
matmul(defgrad(me(1),me(2),me(3),1:3,1:3),0.5_pReal*real(step(1:3,s+1_pInt)/res,pReal))) matmul(defgrad(me(1),me(2),me(3),1:3,1:3),0.5_pReal*real(step(1:3,s+1_pInt)/res,pReal)))
else else
myStep = (me-rear)*geomdim/res myStep = (me-rear)*geomdim/res
coord(s+1_pInt,o,me(1),me(2),me(3),1:3) = coord(s+1_pInt,o,rear(1),rear(2),rear(3),1:3) + & coord(me(1),me(2),me(3),1:3) = coord(rear(1),rear(2),rear(3),1:3) + &
0.5_pReal*matmul(defgrad(me(1),me(2),me(3),1:3,1:3) + & 0.5_pReal*matmul(defgrad(me(1),me(2),me(3),1:3,1:3) + &
defgrad(rear(1),rear(2),rear(3),1:3,1:3),myStep) defgrad(rear(1),rear(2),rear(3),1:3,1:3),myStep)
endif endif
rear = me rear = me
enddo; enddo; enddo; enddo enddo; enddo; enddo
do i = 1_pInt,6_pInt coord_avgOrder(s+1_pInt,1:res(1),1:res(2),1:res(3),1:3) = &
coord_avgOrder(s+1_pInt,1:res(1),1:res(2),1:res(3),1:3) = coord_avgOrder(s+1_pInt, 1:res(1),1:res(2),1:res(3),1:3)& coord_avgOrder(s+1_pInt,1:res(1),1:res(2),1:res(3),1:3) + coord/6.0_pReal
+ coord(s+1_pInt,i,1:res(1),1:res(2),1:res(3),1:3)/6.0_pReal
enddo enddo
offset_coords = coord_avgOrder(s+1,1,1,1,1:3)
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
coord_avgOrder(s+1,i,j,k,1:3) = coord_avgOrder(s+1,i,j,k,1:3) - offset_coords
enddo; enddo; enddo
enddo enddo
do k = 0_pInt, res(3)-1_pInt do k = 0_pInt, res(3)-1_pInt
do j = 0_pInt, res(2)-1_pInt do j = 0_pInt, res(2)-1_pInt
do i = 0_pInt, res(1)-1_pInt do i = 0_pInt, res(1)-1_pInt
parameter_coords = (2.0_pReal*(/real(i,pReal)+0.0_pReal,real(j,pReal)+0.0_pReal,real(k,pReal)+0.0_pReal/)& parameter_coords = (2.0_pReal*real([i,j,k]+1,pReal)-real(res,pReal))/(real(res,pReal))
-real(res,pReal)+fones)/(real(res,pReal)-fones)
positive = fones + parameter_coords positive = fones + parameter_coords
negative = fones - parameter_coords negative = fones - parameter_coords
coord_avgCorner(i+1_pInt,j+1_pInt,k+1_pInt,1:3)& coord(i+1_pInt,j+1_pInt,k+1_pInt,1:3)&
=(coord_avgOrder(1,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *negative(1)*negative(2)*negative(3)& =(coord_avgOrder(1,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *negative(1)*negative(2)*negative(3)&
+ coord_avgOrder(2,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *positive(1)*negative(2)*negative(3)& + coord_avgOrder(2,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *positive(1)*negative(2)*negative(3)&
+ coord_avgOrder(3,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *positive(1)*positive(2)*negative(3)& + coord_avgOrder(3,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *positive(1)*positive(2)*negative(3)&
@ -3284,6 +3271,13 @@ subroutine deformed_linear(res,geomdim,defgrad_av,defgrad,coord_avgCorner)
+ coord_avgOrder(8,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *positive(1)*negative(2)*positive(3))*0.125_pReal + coord_avgOrder(8,i+1_pInt,j+1_pInt,k+1_pInt,1:3) *positive(1)*negative(2)*positive(3))*0.125_pReal
enddo; enddo; enddo enddo; enddo; enddo
offset_coords = matmul(defgrad(1,1,1,1:3,1:3),geomdim/real(res, pReal)/2.0_pReal) - coord(1,1,1,1:3)
do k = 1_pInt, res(3)
do j = 1_pInt, res(2)
do i = 1_pInt, res(1)
coord(i,j,k,1:3) = coord(i,j,k,1:3)+ offset_coords
enddo; enddo; enddo
end subroutine deformed_linear end subroutine deformed_linear
#ifdef Spectral #ifdef Spectral
@ -3407,15 +3401,11 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords)
enddo; enddo; enddo enddo; enddo; enddo
call fftw_destroy_plan(fftw_forth); call fftw_destroy_plan(fftw_back) call fftw_destroy_plan(fftw_forth)
call c_f_pointer(C_NULL_PTR, defgrad_real, [res(1)+2_pInt,res(2),res(3),3_pInt,3_pInt]) ! let all pointers point on NULL-Type call fftw_destroy_plan(fftw_back)
call c_f_pointer(C_NULL_PTR, defgrad_fourier, [res1_red ,res(2),res(3),3_pInt,3_pInt]) call fftw_free(defgrad_fftw)
call c_f_pointer(C_NULL_PTR, coords_real, [res(1)+2_pInt,res(2),res(3),3_pInt]) call fftw_free(coords_fftw)
call c_f_pointer(C_NULL_PTR, coords_fourier,[res1_red ,res(2),res(3),3_pInt])
if(.not. (c_associated(C_LOC(defgrad_real(1,1,1,1,1))) .and. c_associated(C_LOC(defgrad_fourier(1,1,1,1,1)))))& ! Check if pointers are deassociated and free memory
call fftw_free(defgrad_fftw) ! This procedure ensures that optimization do not mix-up lines, because a
if(.not.(c_associated(C_LOC(coords_real(1,1,1,1))) .and. c_associated(C_LOC(coords_fourier(1,1,1,1)))))& ! simple fftw_free(field_fftw) could be done immediately after the last line where field_fftw appears, e.g:
call fftw_free(coords_fftw) ! call c_f_pointer(field_fftw, field_fourier, [res1_red ,res(2),res(3),vec_tens,3])
end subroutine deformed_fft end subroutine deformed_fft

12
processing/post/3Dvisualize.py
View File

@ -344,6 +344,9 @@ parser.add_option('--scaling', dest='scaling', type='float', \
help='scaling of fluctuation [%default]') help='scaling of fluctuation [%default]')
parser.add_option('-u', '--unitlength', dest='unitlength', type='float', \ parser.add_option('-u', '--unitlength', dest='unitlength', type='float', \
help='set unit length for 2D model [%default]') help='set unit length for 2D model [%default]')
parser.add_option('-l', '--linear', dest='linearreconstruction', action='store_true',\
help='use linear reconstruction of geometry [%default]')
parser.set_defaults(defgrad = 'f') parser.set_defaults(defgrad = 'f')
parser.set_defaults(separator = 't') parser.set_defaults(separator = 't')
parser.set_defaults(scalar = []) parser.set_defaults(scalar = [])
@ -359,6 +362,7 @@ parser.set_defaults(scaling = 1.0)
parser.set_defaults(undeformed = False) parser.set_defaults(undeformed = False)
parser.set_defaults(unitlength = 0.0) parser.set_defaults(unitlength = 0.0)
parser.set_defaults(cell = True) parser.set_defaults(cell = True)
parser.set_defaults(linearreconstruction = False)
sep = {'n': '\n', 't': '\t', 's': ' '} sep = {'n': '\n', 't': '\t', 's': ' '}
@ -456,8 +460,12 @@ for filename in args:
defgrad_av = damask.core.math.tensor_avg(res,numpy.reshape(values[:,column['tensor'][options.defgrad]: defgrad_av = damask.core.math.tensor_avg(res,numpy.reshape(values[:,column['tensor'][options.defgrad]:
column['tensor'][options.defgrad]+9], column['tensor'][options.defgrad]+9],
(res[0],res[1],res[2],3,3))) (res[0],res[1],res[2],3,3)))
if options.linearreconstruction:
#centroids = damask.core.math.deformed_linear(res,dim,defgrad_av, centroids = damask.core.math.deformed_linear(res,dim,defgrad_av,
numpy.reshape(values[:,column['tensor'][options.defgrad]:
column['tensor'][options.defgrad]+9],
(res[0],res[1],res[2],3,3)))
else:
centroids = damask.core.math.deformed_fft(res,dim,defgrad_av,options.scaling, centroids = damask.core.math.deformed_fft(res,dim,defgrad_av,options.scaling,
numpy.reshape(values[:,column['tensor'][options.defgrad]: numpy.reshape(values[:,column['tensor'][options.defgrad]:
column['tensor'][options.defgrad]+9], column['tensor'][options.defgrad]+9],