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

tried to implement an inverse reconstruction, meaning to calculate an average deformation gradient from the 8 cornes of a node in reference and current configuration. Not working yet. 

polishing, mainly in fft reconstruction.
added fftw library file and fortran file for type specification
added folder references with literature for visualization/geometry reconstruction
This commit is contained in:
Martin Diehl 2011-02-21 16:30:18 +00:00
parent 8dd1a694a3
commit c399a06c97
24 changed files with 59072 additions and 106 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
documentation/References/ACS_1986_Comparison of an Integration Procedure to Fourier Transform and Data Averaging Procedures in Chromatographic Data Analysis.pdf Normal file
View File

Binary file not shown.

BIN
documentation/References/Brachman_1999_Numerical integration required during Fourier integral analysis.pdf Normal file
View File

Binary file not shown.

BIN
documentation/References/Bronstein_1992_Surface reconstruction from derivatives.pdf Normal file
View File

Binary file not shown.

BIN
documentation/References/Campos_2002_Integration in the Fourier domain for restoration of a function from its slope.pdf Normal file
View File

Binary file not shown.

BIN
documentation/References/Elster_2002_High-accuracy reconstruction of a function f(x) when only derivative is known at discrete measurement points.pdf Normal file
View File

Binary file not shown.

BIN
documentation/References/Fileformats vtk.pdf Normal file
View File

Binary file not shown.

BIN
documentation/References/Huang_Interactive Shape Interpolation through Controllable Dynamic Deformation.pdf Normal file
View File

Binary file not shown.

BIN
documentation/References/Li_2002_A Fast Spectral Method for Active 3D Shape Reconstruction.pdf Normal file
View File

Binary file not shown.

BIN
documentation/References/Mai_Wavefront Reconstruction from Slope Data Within Pupils of Arbitrary Shapes Using Iterative Fourier Transform.pdf Normal file
View File

Binary file not shown.

BIN
documentation/References/ParaView UsersGuide.pdf Normal file
View File

Binary file not shown.

BIN
documentation/References/ParaView Visualization Course.pdf Normal file
View File

Binary file not shown.

BIN
documentation/References/Rainbow Color Map (Still) Considered Harmful.pdf Normal file
View File

Binary file not shown.

BIN
documentation/References/Scientific Visualization with VTK.pdf Normal file
View File

Binary file not shown.

BIN
documentation/References/Shen_2006_Fast-Fourier-transform based numerical integration method for the Rayleigh–Sommerfeld diffraction formula.pdf Normal file
View File

Binary file not shown.

58645
documentation/References/Visualization in the geosciences.pdf Normal file
View File

File diff suppressed because one or more lines are too long

BIN
documentation/References/Visualizing with VTK - A Tutorial.pdf Normal file
View File

Binary file not shown.

BIN
documentation/References/gmsh Manual.pdf Normal file
View File

Binary file not shown.

126
documentation/Visualization/Geometry_Reconstruction.tex Normal file
View File

@ -0,0 +1,126 @@
\documentclass[a4paper,12pt]{article}
\usepackage{bm}
\usepackage{mathbbol} %poor substitute for mathbb. Changed to mathbbg to have an option for greek bb
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{graphicx}
\newcommand{\term}[1]{\textsc{#1}}
\newcommand{\eref}[1]{eq.~\eqref{#1}}
\newcommand{\Eref}[1]{Eq.~\eqref{#1}}
\newcommand{\fref}[1]{fig.~\ref{#1}}
\newcommand{\Fref}[1]{Fig.~\ref{#1}}
\newcommand{\tref}[1]{tab.~\ref{#1}}
\newcommand{\Tref}[1]{Tab.~\ref{#1}}
\newcommand{\cref}[1]{chapter~\ref{#1}}
\newcommand{\Cref}[1]{Chapter~\ref{#1}}
\newcommand{\sref}[1]{section~\ref{#1}}
\newcommand{\Sref}[1]{Section~\ref{#1}}
\newcommand{\lref}[1]{listing~\ref{#1}}
\newcommand{\Lref}[1]{Listing~\ref{#1}}
\newcommand{\ie}{i.e.}
\newcommand{\eg}{e.g.}
\newcommand{\cf}{cf.}
\newcommand{\field}[1]{\ensuremath{\mathcal{#1}}}
\newcommand{\tnsr}[1]{\ensuremath{\bm{{#1}}}}
\newcommand{\tnsrfour}[1]{\ensuremath{\mathbb{#1}}}
\newcommand{\gammaop}{\ensuremath{\mathbbg{\Gamma}}}
\newcommand{\abs}[1]{\ensuremath{\left|{#1}\right|}}
\newcommand{\norm}[1]{\ensuremath{\left|\left|{#1}\right|\right|}}
\newcommand{\vctr}[1]{\ensuremath{\bm{#1}}}
\newcommand{\inc}[1]{\ensuremath{{\rm d}{#1}}}
\newcommand{\sign}[1]{\ensuremath{\operatorname{sign}\left({#1}\right)}}
\newcommand{\grad}[1][]{\ensuremath{\operatorname{grad}{#1}}}
\newcommand{\divergence}[1][]{\ensuremath{\operatorname{div}{#1}}}
\newcommand{\totalder}[2]{\ensuremath{\frac{\inc{#1}}{\inc{#2}}}}
\newcommand{\partialder}[2]{\ensuremath{\frac{\partial{#1}}{\partial{#2}}}}
\newcommand{\identity}{\ensuremath{\tnsr I}}
\newcommand{\F}[1][]{\ensuremath{\tnsr F^{{\rm #1}}}}
\newcommand{\Favg}{{\ensuremath{\overline{\F}^{(n)}}}}
\newcommand{\Fp}[1][]{\ensuremath{\tnsr F_\text{p}^{#1}}}
\newcommand{\Fe}[1][]{\ensuremath{\tnsr F_\text{e}^{#1}}}
\newcommand{\Lp}{\ensuremath{\tnsr L_\text{p}}}
%\newcommand{\Q}[1]{\ensuremath{\tnsr Q^{(#1)}}}
%\newcommand{\x}[2][]{\ensuremath{\vctr x^{(#2)}_\text{#1}}}
%\newcommand{\y}[2][]{\ensuremath{\vctr x^{(#2)}_\text{#1}}}
\newcommand{\dg}[2][]{\ensuremath{\Delta\vctr g^{(#2)}_\text{#1}}}
\newcommand{\g}[1][]{\ensuremath{\vctr g_\text{#1}}}
%\newcommand{\A}[2][]{\ensuremath{A^{(#2)}_\text{#1}}}
\newcommand{\N}[2]{\ensuremath{\varrho^{(#1)}_\text{#2}}}
\newcommand{\Burgers}[1]{\ensuremath{\vctr s^{(#1)}}}
\newcommand{\n}[1]{\ensuremath{\vctr n^{(#1)}}}
\newcommand{\m}[2]{\ensuremath{\vctr m^{(#1)}_{#2}}}
\newcommand{\ld}[1]{\ensuremath{\vctr p^{(#1)}}}
\newcommand{\velocity}[2]{\ensuremath{v^{(#1)}_\text{#2}}}
\newcommand{\avgvelocity}[2]{\ensuremath{{\overline v}^{(#1)}_ \text{#2}}}
\newcommand{\flux}[2]{\ensuremath{\vctr f^{(#1)}_ \text{#2}}}
\newcommand{\averageflux}[2]{\ensuremath{\overline{\vctr f}^{(#1)}_ \text{#2}}}
\newcommand{\interfaceflux}[2]{\ensuremath{\tilde{\vctr f}^{(#1)}_ \text{#2}}}
\newcommand{\transmissivity}[1]{\ensuremath{\chi^{(#1)}}}
\newcommand{\galpha}{\ensuremath{\gamma^{(\alpha)}}}
\newcommand{\dotgalpha}{\ensuremath{\dot{\gamma}^{(\alpha)}}}
\newcommand{\taualpha}{\ensuremath{\tau^{(\alpha)}}}
\newcommand{\taualphamax}{\ensuremath{\hat\tau^{(\alpha)}}}
\newcommand{\density}[2]{\ensuremath{\varrho^{(#1)}_ \text{#2}}}
\newcommand{\densityfunc}[2]{\ensuremath{{\tilde\varrho}^{(#1)}_ \text{#2}}}
\newcommand{\avgdensity}[2]{\ensuremath{{\overline\varrho}^{(#1)}_ \text{#2}}}
\newcommand{\dotdensity}[2]{\ensuremath{\dot{\varrho}^{(#1)}_ \text{#2}}}
\newcommand{\densityexcess}[2]{\ensuremath{\Delta\varrho^{(#1)}_ \text{#2}}}
\newcommand{\cs}[2][]{\ensuremath{\sigma^{(#1)}_ \text{#2}}}
\title{Geometry reconstruction using Fast Fourier transform}
\date{\today}
\author{Martin Diehl}
\begin{document}
\maketitle
The presented method allows the shape reconstruction of a volume element with periodic boundary conditions.
The deformation gradient on each point of a regular, three-dimensional lattice in undeformed configuration must be known.
The deformation gradient maps each point (or a line in the infinitesimal neighborhood of the point) into the current configuration.
It is defined as:
\begin{equation}
\F(\vctr x) = \left(
\begin{array}{ccc}
\frac{\partial y_1}{\partial x_1} & \frac{\partial y_1}{\partial x_2} & \frac{\partial y_1}{\partial x_3} \\
\frac{\partial y_2}{\partial x_1} & \frac{\partial y_2}{\partial x_2} & \frac{\partial y_2}{\partial x_3} \\
\frac{\partial y_3}{\partial x_1} & \frac{\partial y_3}{\partial x_2} & \frac{\partial y_3}{\partial x_3} \\
\end{array}
\right)
\end{equation}
with \vctr y\ are the coordinates in current configuration and \vctr x\ are the coordinates in reference configuration.
The three-dimensional field of second order tensors is transformed to the Fourier space, giving three-dimensional field of second order tensors that depend on the three dimensional wave vector and not on the vector of spatial coordinates:
\begin{equation}
\mathcal F \left( \F(\vctr x) \right)= \hat{\F}(\vctr k)
\end{equation}
Integration in Fourier space works is defined for the one dimensional case as:
\begin{equation}
\mathcal{F} \left( \int_{-\infty}^{x} g (\tau) d \tau \right) = \frac{\mathcal{F}{g(x)}}{i2 \pi k} + c \delta(k)
\end{equation}
where the last term is the integration constant.
Constant (or linear after integration) terms cannot properly handled when using the integration property of the Fourier domain, as a division by the ``constant wave'' ($k=0$) is not possible.
Thus, to carry out the integration the function is separated into an average and a locally fluctuating part.
The locally fluctuating part is integrated in Fourier space, while the integration of the constant part is easier in the real domain.
The fluctuation field of the position vector in deformed configuration in Fourier space reads as:
\begin{equation}
\hat{\tilde{y}}_{\rm j}(\vctr k) = \hat{F_{\rm ji}}(\vctr k) \left(k_{\rm i} i 2 \pi \right)^{-1} \forall k_i \neq 0
\end{equation}
The average part is set to zero in Fourier space.
The inverse Fourier transform gives the locally fluctuating part of each position in current configuration:
\begin{equation}
\mathcal{F}^{-1}\left(\hat{\tilde{\vctr y}}(\vctr k) \right) = \tilde{\vctr y}(\vctr x)
\end{equation}
and the position vector in undeformed configuration is given as:
\begin{equation}
\vctr y(\vctr x) = \overline{\F}: \vctr x + \tilde{\vctr y}(\vctr x)
\end{equation}
\end{document}

72
processing/post/fftw3.f Normal file
View File

@ -0,0 +1,72 @@
INTEGER FFTW_R2HC
PARAMETER (FFTW_R2HC=0)
INTEGER FFTW_HC2R
PARAMETER (FFTW_HC2R=1)
INTEGER FFTW_DHT
PARAMETER (FFTW_DHT=2)
INTEGER FFTW_REDFT00
PARAMETER (FFTW_REDFT00=3)
INTEGER FFTW_REDFT01
PARAMETER (FFTW_REDFT01=4)
INTEGER FFTW_REDFT10
PARAMETER (FFTW_REDFT10=5)
INTEGER FFTW_REDFT11
PARAMETER (FFTW_REDFT11=6)
INTEGER FFTW_RODFT00
PARAMETER (FFTW_RODFT00=7)
INTEGER FFTW_RODFT01
PARAMETER (FFTW_RODFT01=8)
INTEGER FFTW_RODFT10
PARAMETER (FFTW_RODFT10=9)
INTEGER FFTW_RODFT11
PARAMETER (FFTW_RODFT11=10)
INTEGER FFTW_FORWARD
PARAMETER (FFTW_FORWARD=-1)
INTEGER FFTW_BACKWARD
PARAMETER (FFTW_BACKWARD=+1)
INTEGER FFTW_MEASURE
PARAMETER (FFTW_MEASURE=0)
INTEGER FFTW_DESTROY_INPUT
PARAMETER (FFTW_DESTROY_INPUT=1)
INTEGER FFTW_UNALIGNED
PARAMETER (FFTW_UNALIGNED=2)
INTEGER FFTW_CONSERVE_MEMORY
PARAMETER (FFTW_CONSERVE_MEMORY=4)
INTEGER FFTW_EXHAUSTIVE
PARAMETER (FFTW_EXHAUSTIVE=8)
INTEGER FFTW_PRESERVE_INPUT
PARAMETER (FFTW_PRESERVE_INPUT=16)
INTEGER FFTW_PATIENT
PARAMETER (FFTW_PATIENT=32)
INTEGER FFTW_ESTIMATE
PARAMETER (FFTW_ESTIMATE=64)
INTEGER FFTW_ESTIMATE_PATIENT
PARAMETER (FFTW_ESTIMATE_PATIENT=128)
INTEGER FFTW_BELIEVE_PCOST
PARAMETER (FFTW_BELIEVE_PCOST=256)
INTEGER FFTW_NO_DFT_R2HC
PARAMETER (FFTW_NO_DFT_R2HC=512)
INTEGER FFTW_NO_NONTHREADED
PARAMETER (FFTW_NO_NONTHREADED=1024)
INTEGER FFTW_NO_BUFFERING
PARAMETER (FFTW_NO_BUFFERING=2048)
INTEGER FFTW_NO_INDIRECT_OP
PARAMETER (FFTW_NO_INDIRECT_OP=4096)
INTEGER FFTW_ALLOW_LARGE_GENERIC
PARAMETER (FFTW_ALLOW_LARGE_GENERIC=8192)
INTEGER FFTW_NO_RANK_SPLITS
PARAMETER (FFTW_NO_RANK_SPLITS=16384)
INTEGER FFTW_NO_VRANK_SPLITS
PARAMETER (FFTW_NO_VRANK_SPLITS=32768)
INTEGER FFTW_NO_VRECURSE
PARAMETER (FFTW_NO_VRECURSE=65536)
INTEGER FFTW_NO_SIMD
PARAMETER (FFTW_NO_SIMD=131072)
INTEGER FFTW_NO_SLOW
PARAMETER (FFTW_NO_SLOW=262144)
INTEGER FFTW_NO_FIXED_RADIX_LARGE_N
PARAMETER (FFTW_NO_FIXED_RADIX_LARGE_N=524288)
INTEGER FFTW_ALLOW_PRUNING
PARAMETER (FFTW_ALLOW_PRUNING=1048576)
INTEGER FFTW_WISDOM_ONLY
PARAMETER (FFTW_WISDOM_ONLY=2097152)

BIN
processing/post/libfftw3.a Normal file
View File

Binary file not shown.

2
processing/post/postResults
View File

@ -75,6 +75,7 @@ class MPIEspectral_result: # mimic py_post result object
dimension = [0.0,0.0,0.0] dimension = [0.0,0.0,0.0]
theTitle = '' theTitle = ''
wd = '' wd = ''
geometry = ''
extrapolate = '' extrapolate = ''
N_increments = 0 N_increments = 0
increment = 0 increment = 0
@ -123,6 +124,7 @@ class MPIEspectral_result: # mimic py_post result object
return '\n'.join([ return '\n'.join([
'title: %s'%self.theTitle, 'title: %s'%self.theTitle,
'workdir: %s'%self.wd, 'workdir: %s'%self.wd,
'geometry: %s'%self.geometry,
'extrapolation: %s'%self.extrapolate, 'extrapolation: %s'%self.extrapolate,
'increments: %i'%self.N_increments, 'increments: %i'%self.N_increments,
'increment: %i'%self.increment, 'increment: %i'%self.increment,

212
processing/post/postprocessingMath.f90
View File

@ -435,37 +435,19 @@ subroutine deformed_fft(res_x,res_y,res_z,geomdim,defgrad,defgrad_av,scaling,coo
implicit none implicit none
integer res_x, res_y, res_z integer res_x, res_y, res_z
real*8 geomdim(3) real*8 geomdim(3)
real*8 defgrad(res_x,res_y,res_z,3,3) real*8 defgrad(res_x,res_y,res_z,3,3)
complex*16 defgrad_fft(res_x,res_y,res_z,3,3)
real*8 defgrad_av(3,3) real*8 defgrad_av(3,3)
real*8 scaling real*8 scaling
real*8 coords(res_x, res_y,res_z,3) real*8 coords(res_x,res_y,res_z,3)
complex*16 coords_fft(res_x/2+1,res_y,res_z,3) complex*16 coords_fft(res_x/2+1,res_y,res_z,3)
real*8 waves(res_x/2+1,res_y,res_z,3) complex*16 defgrad_fft(res_x,res_y,res_z,3,3)
include 'fftw3.f' !header file for fftw3 (declaring variables). Library files are also needed integer i, j, k
integer*8 :: plan_fft(2) integer k_s(3)
real*8, parameter :: pi = 3.14159265358979323846264338327950288419716939937510
real*8 zero
real*8 temp33_Real(3,3)
integer i, j, k, ierr
integer k_s(3)
real*8 step(3) real*8 step(3)
complex*16 img real*8 offset_coords(3)
real*8, parameter :: pi = 3.14159265358979323846264338327950288419716939937510
img = cmplx(0.0,1.0) integer*8 :: plan_fft(2)
include 'fftw3.f' !header file for fftw3 (declaring variables). Library files are also needed
do k = 1, res_z
k_s(3) = k-1
if(k > res_z/2+1) k_s(3) = k_s(3)-res_z
do j = 1, res_y
k_s(2) = j-1
if(j > res_y/2+1) k_s(2) = k_s(2)-res_y
do i = 1, res_x/2+1
k_s(1) = i-1
waves(i,j,k,:) = real(k_s)/geomdim
enddo; enddo; enddo
call dfftw_plan_many_dft(plan_fft(1),3,(/res_x,res_y,res_z/),9,& call dfftw_plan_many_dft(plan_fft(1),3,(/res_x,res_y,res_z/),9,&
defgrad_fft,(/res_x,res_y,res_z/),1,res_x*res_y*res_z,& defgrad_fft,(/res_x,res_y,res_z/),1,res_x*res_y*res_z,&
@ -475,68 +457,142 @@ subroutine deformed_fft(res_x,res_y,res_z,geomdim,defgrad,defgrad_av,scaling,coo
coords_fft,(/res_x/2+1,res_y,res_z/),1,(res_x/2+1)*res_y*res_z,& coords_fft,(/res_x/2+1,res_y,res_z/),1,(res_x/2+1)*res_y*res_z,&
coords, (/res_x, res_y,res_z/),1, res_x* res_y*res_z,FFTW_PATIENT) coords, (/res_x, res_y,res_z/),1, res_x* res_y*res_z,FFTW_PATIENT)
coords_fft=0.0 coords_fft = 0.0
defgrad_fft = defgrad defgrad_fft = defgrad
call dfftw_execute_dft(plan_fft(1), defgrad_fft, defgrad_fft)
do k = 1, res_z; do j = 1, res_y; do i = 1, res_x/2+1
if(i/=1) then
coords_fft(i,j,k,1) = defgrad_fft(i,j,k,1,1)/(waves(i,j,k,1)*img)
coords_fft(i,j,k,2) = defgrad_fft(i,j,k,2,1)/(waves(i,j,k,1)*img)
coords_fft(i,j,k,3) = defgrad_fft(i,j,k,3,1)/(waves(i,j,k,1)*img)
endif
if(j/=1) then
coords_fft(i,j,k,1) = coords_fft(i,j,k,1) + defgrad_fft(i,j,k,1,2)/(waves(i,j,k,2)*img)
coords_fft(i,j,k,2) = coords_fft(i,j,k,2) + defgrad_fft(i,j,k,2,2)/(waves(i,j,k,2)*img)
coords_fft(i,j,k,3) = coords_fft(i,j,k,3) + defgrad_fft(i,j,k,3,2)/(waves(i,j,k,2)*img)
endif
if(k/=1) then
coords_fft(i,j,k,1) = coords_fft(i,j,k,1) + defgrad_fft(i,j,k,1,3)/(waves(i,j,k,3)*img)
coords_fft(i,j,k,2) = coords_fft(i,j,k,2) + defgrad_fft(i,j,k,2,3)/(waves(i,j,k,3)*img)
coords_fft(i,j,k,3) = coords_fft(i,j,k,3) + defgrad_fft(i,j,k,3,3)/(waves(i,j,k,3)*img)
endif
enddo; enddo; enddo
call dfftw_execute_dft_c2r(plan_fft(2), coords_fft, coords)
coords = coords/real(res_x*res_y*res_z)
do k = 1, res_z; do j = 1, res_y; do i = 1, res_x
coords(i,j,k,:) = coords(i,j,k,:)/(geomdim*2.0*pi)
enddo; enddo; enddo
step(1) = geomdim(1)/real(res_x) step(1) = geomdim(1)/real(res_x)
step(2) = geomdim(2)/real(res_y) step(2) = geomdim(2)/real(res_y)
step(3) = geomdim(3)/real(res_z) step(3) = geomdim(3)/real(res_z)
call dfftw_execute_dft(plan_fft(1), defgrad_fft, defgrad_fft)
temp33_Real(1,:) = matmul(defgrad_av,step/2.0) & do k = 1, res_z
- matmul(defgrad_av,(/zero,zero,step(3)/)) ! start below origin k_s(3) = k-1
if(k > res_z/2+1) k_s(3) = k_s(3)-res_z
do k = 1, res_z; do j = 1, res_y; do i = 1, res_x do j = 1, res_y
if((j==1).and.(i==1)) then k_s(2) = j-1
temp33_Real(1,:) = temp33_Real(1,:) + matmul(defgrad_av,& if(j > res_y/2+1) k_s(2) = k_s(2)-res_y
(/zero,zero,step(3)/)) do i = 1, res_x/2+1
temp33_Real(2,:) = temp33_Real(1,:) k_s(1) = i-1
temp33_Real(3,:) = temp33_Real(1,:) if(i/=1) coords_fft(i,j,k,:) = coords_fft(i,j,k,:)&
coords(i,j,k,:) = coords(i,j,k,:) + temp33_Real(1,:) + defgrad_fft(i,j,k,:,1)*geomdim(1)/(real(k_s(1))*cmplx(0.0,1.0)*pi*2.0)
else if(j/=1) coords_fft(i,j,k,:) = coords_fft(i,j,k,:)&
if(i==1) then + defgrad_fft(i,j,k,:,2)*geomdim(2)/(real(k_s(2))*cmplx(0.0,1.0)*pi*2.0)
temp33_Real(2,:) = temp33_Real(2,:) + matmul(defgrad_av,& if(k/=1) coords_fft(i,j,k,:) = coords_fft(i,j,k,:)&
(/zero,step(2),zero/)) + defgrad_fft(i,j,k,:,3)*geomdim(3)/(real(k_s(3))*cmplx(0.0,1.0)*pi*2.0)
temp33_Real(3,:) = temp33_Real(2,:)
coords(i,j,k,:) = coords(i,j,k,:) + temp33_Real(2,:)
else
temp33_Real(3,:) = temp33_Real(3,:) + matmul(defgrad_av,&
(/step(1),zero,zero/))
coords(i,j,k,:) = coords(i,j,k,:) + temp33_Real(3,:)
endif
endif
enddo; enddo; enddo enddo; enddo; enddo
call dfftw_execute_dft_c2r(plan_fft(2), coords_fft, coords)
coords = coords/real(res_x*res_y*res_z)
offset_coords = matmul(defgrad(1,1,1,:,:),step/2.0) - scaling*coords(1,1,1,:)
do k = 1, res_z; do j = 1, res_y; do i = 1, res_x
coords(i,j,k,:) = scaling*coords(i,j,k,:) + offset_coords + matmul(defgrad_av,&
(/step(1)*real(i-1),&
step(2)*real(j-1),&
step(3)*real(k-1)/))
enddo; enddo; enddo
end subroutine deformed_fft end subroutine deformed_fft
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine inverse_reconstruction(res_x,res_y,res_z,reference_configuration,current_configuration,defgrad)
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
use math
implicit none
integer res_x, res_y, res_z
real*8 reference_configuration(res_x+1,res_y+1,res_z+1,3)
real*8 current_configuration(res_x+1,res_y+1,res_z+1,3)
real*8 defgrad(res_x,res_y,res_z,3,3)
real*8 delta, tolerance, res, res_center
real*8 reference(8,3)
real*8 current(8,3)
real*8 defgrad_temp(3,3)
real*8 dres_dF(3,3)
real*8 identity(3,3)
real*8 ref_bar(3)
real*8 current_bar(3)
real*8 r(8)
real*8 differentiate(9,3,3)
integer i, j, k, m, l, x, y, o
identity = 0.0
identity(1,1) = 1.0
identity(2,2) = 1.0
identity(3,3) = 1.0
differentiate = 0.0
tolerance = 1e-10
delta = 1e-9
k = 0
do j = 1, 3; do i = 1, 3
k = k+1
differentiate(k,i,j) = 1.0
enddo; enddo
do k = 1, res_z
do j = 1, res_y
do i = 1, res_x
reference(1,:) = reference_configuration(i ,j ,k ,:)
reference(2,:) = reference_configuration(i+1,j ,k ,:)
reference(3,:) = reference_configuration(i+1,j+1,k ,:)
reference(4,:) = reference_configuration(i ,j+1,k ,:)
reference(5,:) = reference_configuration(i ,j ,k+1,:)
reference(6,:) = reference_configuration(i+1,j ,k+1,:)
reference(7,:) = reference_configuration(i+1,j+1,k+1,:)
reference(8,:) = reference_configuration(i ,j+1,k+1,:)
current(1,:) = current_configuration(i ,j ,k ,:)
current(2,:) = current_configuration(i+1,j ,k ,:)
current(3,:) = current_configuration(i+1,j+1,k ,:)
current(4,:) = current_configuration(i ,j+1,k ,:)
current(5,:) = current_configuration(i ,j ,k+1,:)
current(6,:) = current_configuration(i+1,j ,k+1,:)
current(7,:) = current_configuration(i+1,j+1,k+1,:)
current(8,:) = current_configuration(i ,j+1,k+1,:)
do o=1,3
ref_bar(o) = sum(reference(:,o))/8.0
current_bar(o) = sum(current(:,o))/8.0
enddo
do o=1,8
reference(o,:) = reference(o,:) -ref_bar
current(o,:) = current(o,:) -current_bar
enddo
defgrad_temp = identity
res_center = 2.0*tolerance
o=0
do while(res_center >= tolerance)
o = o + 1
do l = 1,8 ! loop over corners
r(l) = sqrt(sum((current(l,:)-matmul(defgrad_temp,reference(l,:)))**2)) ! corner distance
enddo
res_center = sum(r*r) ! current residuum
print*, 'res_center', res_center
m=0
do y=1,3; do x=1,3 ! numerical differentiation
m = m+1
do l = 1,8
r(l) = sqrt(sum((current(l,:)-matmul((defgrad_temp+differentiate(m,:,:)*delta),reference(l,:)))**2)) ! corner distance
enddo
res = sum(r*r)
print*,'res step', m, res
dres_dF(x,y) = (res-res_center)/delta
enddo; enddo
print*, 'dres_dF', dres_dF
print*, 'deltadef', math_inv3x3(dres_dF)*res_center
defgrad_temp = defgrad_temp - math_inv3x3(dres_dF)*res_center ! Newton--Raphson
print*, o, res_center
pause
enddo
defgrad(i,j,k,:,:) = defgrad_temp
enddo; enddo; enddo
end subroutine inverse_reconstruction
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine tensor_avg(res_x,res_y,res_z,tensor,avg) subroutine tensor_avg(res_x,res_y,res_z,tensor,avg)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

54
processing/post/postprocessingMath.pyf
View File

@ -3,6 +3,50 @@
python module postprocessingMath ! in python module postprocessingMath ! in
interface ! in :postprocessingMath interface ! in :postprocessingMath
module math ! in :postprocessingMath:postprocessingMath.f90
real*8 parameter,optional,dimension(3,3) :: math_i3=reshape((/1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0/),(/3,3/))
real*8 parameter,optional :: pi=3.14159265359
function math_mul33x33(a,b) ! in :postprocessingMath:postprocessingMath.f90:math
real*8 dimension(3,3),intent(in) :: a
real*8 dimension(3,3),intent(in) :: b
real*8 dimension(3,3) :: math_mul33x33
end function math_mul33x33
subroutine math_invert3x3(a,inva,deta,error) ! in :postprocessingMath:postprocessingMath.f90:math
real*8 dimension(3,3),intent(in) :: a
real*8 dimension(3,3),intent(out) :: inva
real*8 intent(out) :: deta
logical intent(out) :: error
end subroutine math_invert3x3
function math_det3x3(m) ! in :postprocessingMath:postprocessingMath.f90:math
real*8 dimension(3,3),intent(in) :: m
real*8 :: math_det3x3
end function math_det3x3
subroutine math_pdecomposition(fe,u,r,error) ! in :postprocessingMath:postprocessingMath.f90:math
real*8 dimension(3,3),intent(in) :: fe
real*8 dimension(3,3),intent(out) :: u
real*8 dimension(3,3),intent(out) :: r
logical intent(out) :: error
end subroutine math_pdecomposition
function math_inv3x3(a) ! in :postprocessingMath:postprocessingMath.f90:math
real*8 dimension(3,3),intent(in) :: a
real*8 dimension(3,3) :: math_inv3x3
end function math_inv3x3
subroutine math_hi(m,hi1m,hi2m,hi3m) ! in :postprocessingMath:postprocessingMath.f90:math
real*8 dimension(3,3),intent(in) :: m
real*8 intent(out) :: hi1m
real*8 intent(out) :: hi2m
real*8 intent(out) :: hi3m
end subroutine math_hi
subroutine math_spectral1(m,ew1,ew2,ew3,eb1,eb2,eb3) ! in :postprocessingMath:postprocessingMath.f90:math
real*8 dimension(3,3),intent(in) :: m
real*8 intent(out) :: ew1
real*8 intent(out) :: ew2
real*8 intent(out) :: ew3
real*8 dimension(3,3),intent(out) :: eb1
real*8 dimension(3,3),intent(out) :: eb2
real*8 dimension(3,3),intent(out) :: eb3
end subroutine math_spectral1
end module math
subroutine mesh(res_x,res_y,res_z,geomdim,defgrad_av,centroids,nodes) ! in :postprocessingMath:postprocessingMath.f90 subroutine mesh(res_x,res_y,res_z,geomdim,defgrad_av,centroids,nodes) ! in :postprocessingMath:postprocessingMath.f90
integer intent(in) :: res_x integer intent(in) :: res_x
integer intent(in) :: res_y integer intent(in) :: res_y
@ -35,8 +79,16 @@ python module postprocessingMath ! in
complex*16 dimension(res_x/2+1,res_y,res_z,3),depend(res_x,res_y,res_z) :: coords_fft complex*16 dimension(res_x/2+1,res_y,res_z,3),depend(res_x,res_y,res_z) :: coords_fft
real*8 dimension(res_x,res_y,res_z,3,3),intent(in),depend(res_x,res_y,res_z) :: defgrad real*8 dimension(res_x,res_y,res_z,3,3),intent(in),depend(res_x,res_y,res_z) :: defgrad
complex*16 dimension(res_x,res_y,res_z,3,3),depend(res_x,res_y,res_z) :: defgrad_fft complex*16 dimension(res_x,res_y,res_z,3,3),depend(res_x,res_y,res_z) :: defgrad_fft
real*8 dimension(res_x/2+1,res_y,res_z,3),depend(res_x,res_y,res_z) :: waves complex*16 dimension(res_x/2+1,res_y,res_z,3),depend(res_x,res_y,res_z) :: integrator
end subroutine deformed_fft end subroutine deformed_fft
subroutine inverse_reconstruction(res_x,res_y,res_z,reference_configuration,current_configuration,defgrad) ! in :postprocessingMath:postprocessingMath.f90
integer intent(in) :: res_x
integer intent(in) :: res_y
integer intent(in) :: res_z
real*8 dimension(res_x+1,res_y+1,res_z+1,3),intent(in),depend(res_x,res_y,res_z) :: reference_configuration
real*8 dimension(res_x+1,res_y+1,res_z+1,3),intent(in),depend(res_x,res_y,res_z) :: current_configuration
real*8 dimension(res_x,res_y,res_z,3,3),intent(out),depend(res_x,res_y,res_z) :: defgrad
end subroutine inverse_reconstruction
subroutine tensor_avg(res_x,res_y,res_z,tensor,avg) ! in :postprocessingMath:postprocessingMath.f90 subroutine tensor_avg(res_x,res_y,res_z,tensor,avg) ! in :postprocessingMath:postprocessingMath.f90
integer intent(in) :: res_x integer intent(in) :: res_x
integer intent(in) :: res_y integer intent(in) :: res_y

63
processing/post/spectral_post.py
View File

@ -6,7 +6,7 @@
# computed using the FEM solvers. Until now, its capable to handle elements with one IP in a regular order # 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 # written by M. Diehl, m.diehl@mpie.de
import os,sys,re,array,struct,numpy, time, postprocessingMath import os,sys,re,array,struct,numpy, time, postprocessingMath, math
class vector: class vector:
x,y,z = [None,None,None] x,y,z = [None,None,None]
@ -306,7 +306,7 @@ print 'Post Processing for Material subroutine for BVP solution using spectral m
print '*********************************************************************************\n' print '*********************************************************************************\n'
#reading in the header of the results file #reading in the header of the results file
name = 'dipl10' name = 'dipl32_shear'
p = MPIEspectral_result(name+'.spectralOut') p = MPIEspectral_result(name+'.spectralOut')
p.extrapolation('') p.extrapolation('')
print p print p
@ -317,37 +317,50 @@ res_x=p.resolution[0]
res_y=p.resolution[1] res_y=p.resolution[1]
res_z=p.resolution[2] res_z=p.resolution[2]
# for i in range(1,3):
# print('new step')
# c_pos = p.dataOffset + i*(p.N_element_scalars*8*p.N_elements + 8) #8 accounts for header&footer
# for j in range(p.N_element_scalars):
# #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 z in range(0,resolution[2]):
# #for y in range(0,resolution[1]):
# #for x in range(0,resolution[0]):
# currentPosition = c_pos + j*8 +4
# p.file.seek(currentPosition)
# print(struct.unpack('d',p.file.read(8)))
ms=numpy.zeros([res_x,res_y,res_z,3], 'd') ms=numpy.zeros([res_x,res_y,res_z,3], 'd')
for i in range(249,250): print 'data structure'
for i in range(p.N_element_scalars):
c_pos = p.dataOffset + i*8.0 + 4.0
p.file.seek(c_pos)
print(i, struct.unpack('d',p.file.read(8)))
for i in range(1,2):
c_pos = p.dataOffset + i*(p.N_element_scalars*8*p.N_elements + 8) #8 accounts for header&footer c_pos = p.dataOffset + i*(p.N_element_scalars*8*p.N_elements + 8) #8 accounts for header&footer
defgrad = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,16) defgrad = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,7)
logstrain = postprocessingMath.logstrain_mat(res_x,res_y,res_z,defgrad) defgrad_av = postprocessingMath.tensor_avg(res_x,res_y,res_z,defgrad)
logstrain2 = postprocessingMath.logstrain_spat(res_x,res_y,res_z,defgrad) centroids_coord = postprocessingMath.deformed(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av)
p_stress = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,58) undeformed = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord)
c_stress = postprocessingMath.calculate_cauchy(res_x,res_y,res_z,defgrad,p_stress) #writeVtkAscii(name+'-mesh-undeformed.vtk',undeformed,defgrad[:,:,:,1,2],p.resolution)
vm = postprocessingMath.calculate_mises(res_x,res_y,res_z,c_stress)
for i in range(240,241):
c_pos = p.dataOffset + i*(p.N_element_scalars*8*p.N_elements + 8) #8 accounts for header&footer
defgrad = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,7)
defgrad_av = postprocessingMath.tensor_avg(res_x,res_y,res_z,defgrad)
#defgrad = numpy.zeros([p.resolution[0],p.resolution[1],p.resolution[2],3,3], 'd')
#for z in range(p.resolution[2]):
# defgrad[:,:,z,1,2] = (2.0*z)/(p.resolution[2]-1.0)+ 3.8*math.sin(z*20.0/(p.resolution[2]-1.0)*2*math.pi)
# defgrad[:,:,z,0,2] = (2.0*z)/(p.resolution[2]-1.0)+ 5.0*math.cos(z/(p.resolution[2]-1.0)*2*math.pi)
#defgrad[:,:,:,0,0] = 1.0
#defgrad[:,:,:,1,1] = 1.0
#defgrad[:,:,:,2,2] = 1.0
#logstrain = postprocessingMath.logstrain_mat(res_x,res_y,res_z,defgrad)
#logstrain2 = postprocessingMath.logstrain_spat(res_x,res_y,res_z,defgrad)
#p_stress = readTensor(p.resolution,p.file,p.N_element_scalars,c_pos,52)
#c_stress = postprocessingMath.calculate_cauchy(res_x,res_y,res_z,defgrad,p_stress)
#vm = postprocessingMath.calculate_mises(res_x,res_y,res_z,c_stress)
#defgrad_av = postprocessingMath.tensor_avg(res_x,res_y,res_z,defgrad)
#subroutine inverse_reconstruction(res_x,res_y,res_z,reference_configuration,current_configuration,defgrad)
centroids_coord = postprocessingMath.deformed(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av)
deformed = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord)
#writeVtkAscii(name+'-mesh-deformed.vtk',deformed,defgrad[:,:,:,1,2],p.resolution)
defgrad = postprocessingMath.inverse_reconstruction(res_x,res_y,res_z,undeformed,deformed)
defgrad_av = postprocessingMath.tensor_avg(res_x,res_y,res_z,defgrad) defgrad_av = postprocessingMath.tensor_avg(res_x,res_y,res_z,defgrad)
centroids_coord = postprocessingMath.deformed(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av) centroids_coord = postprocessingMath.deformed(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av)
centroids_coord2 = postprocessingMath.deformed_fft(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av,1.0) centroids_coord2 = postprocessingMath.deformed_fft(res_x,res_y,res_z,p.dimension,defgrad,defgrad_av,1.0)
ms = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord) ms = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord)
ms2 = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord2) ms2 = postprocessingMath.mesh(p.resolution[0],p.resolution[1],p.resolution[2],p.dimension,defgrad_av,centroids_coord2)
writeVtkAscii(name+'-mesh-usual-%s.vtk'%i,ms,logstrain[:,:,:,1,2],p.resolution) writeVtkAscii(name+'-mesh-usual-%s.vtk'%i,ms,defgrad[:,:,:,1,2],p.resolution)
writeVtkAscii(name+'-mesh-fft-%s.vtk'%i,ms2,logstrain[:,:,:,1,2],p.resolution) writeVtkAscii(name+'-mesh-fft-%s.vtk'%i,ms2,defgrad[:,:,:,1,2],p.resolution)
#writeVtkAsciidefgrad_av(name+'-box-%i.vtk'%i,p.dimension,defgrad_av) #writeVtkAsciidefgrad_av(name+'-box-%i.vtk'%i,p.dimension,defgrad_av)
#writeVtkAsciiDots(name+'-points-%i.vtk'%i,centroids_coord,grain,p.resolution) #writeVtkAsciiDots(name+'-points-%i.vtk'%i,centroids_coord,grain,p.resolution)
sys.stdout.flush() sys.stdout.flush()