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

changed all remaining routines to fortran-fast arrays (geometry reconstruction etc.) 

changed all remaining routines in f2py to more clever determination of array size (requires f2py >= 2.0)
enabled 3D visualize to work with odd resolution by switching to linear reconstruction
PLEASE NOTE: Redefinition of routines for f2py might cause trouble -> DELETE DAMASK_ROOT/lib/damask/core.so in this case
further changes: added pure statement where possible, polished, unified use of "Q" for "Quaternion" and reordered math to have similar routines together
This commit is contained in:
Martin Diehl 2013-01-31 16:28:08 +00:00
parent be06f34242
commit e74b5da19a
10 changed files with 1251 additions and 1439 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
code/constitutive_nonlocal.f90
View File

@ -2598,7 +2598,7 @@ subroutine constitutive_nonlocal_updateCompatibility(orientation,i,e)
use prec, only: pReal, & use prec, only: pReal, &
pInt pInt
use math, only: math_QuaternionDisorientation, & use math, only: math_qDisorientation, &
math_mul3x3, & math_mul3x3, &
math_qRot math_qRot
use material, only: material_phase, & use material, only: material_phase, &
@ -2727,7 +2727,7 @@ do n = 1_pInt,Nneighbors
!* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one.
!* All values below the threshold are set to zero. !* All values below the threshold are set to zero.
else else
absoluteMisorientation = math_QuaternionDisorientation(orientation(1:4,1,i,e), & absoluteMisorientation = math_qDisorientation(orientation(1:4,1,i,e), &
orientation(1:4,1,neighboring_i,neighboring_e), & orientation(1:4,1,neighboring_i,neighboring_e), &
0_pInt) ! no symmetry 0_pInt) ! no symmetry
do s1 = 1_pInt,ns ! my slip systems do s1 = 1_pInt,ns ! my slip systems

24
code/crystallite.f90
View File

@ -3357,8 +3357,8 @@ subroutine crystallite_orientations
!*** variables and functions from other modules ***! !*** variables and functions from other modules ***!
use math, only: math_pDecomposition, & use math, only: math_pDecomposition, &
math_RtoQuaternion, & math_RtoQ, &
math_QuaternionDisorientation, & math_qDisorientation, &
math_qConj math_qConj
use FEsolving, only: FEsolving_execElem, & use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP FEsolving_execIP
@ -3409,9 +3409,9 @@ logical error
call IO_warning(650_pInt, e, i, g) call IO_warning(650_pInt, e, i, g)
orientation = [1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal] ! fake orientation orientation = [1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal] ! fake orientation
else else
orientation = math_RtoQuaternion(transpose(R)) orientation = math_RtoQ(transpose(R))
endif endif
crystallite_rotation(1:4,g,i,e) = math_QuaternionDisorientation(crystallite_orientation0(1:4,g,i,e), & ! active rotation from ori0 crystallite_rotation(1:4,g,i,e) = math_qDisorientation(crystallite_orientation0(1:4,g,i,e), & ! active rotation from ori0
orientation, & ! to current orientation orientation, & ! to current orientation
0_pInt ) ! we don't want symmetry here 0_pInt ) ! we don't want symmetry here
crystallite_orientation(1:4,g,i,e) = orientation crystallite_orientation(1:4,g,i,e) = orientation
@ -3446,7 +3446,7 @@ logical error
neighboringStructure = constitutive_nonlocal_structure(neighboringInstance) ! get my neighbor's crystal structure neighboringStructure = constitutive_nonlocal_structure(neighboringInstance) ! get my neighbor's crystal structure
if (myStructure == neighboringStructure) then ! if my neighbor has same crystal structure like me if (myStructure == neighboringStructure) then ! if my neighbor has same crystal structure like me
crystallite_disorientation(:,n,1,i,e) = & crystallite_disorientation(:,n,1,i,e) = &
math_QuaternionDisorientation( crystallite_orientation(1:4,1,i,e), & math_qDisorientation( crystallite_orientation(1:4,1,i,e), &
crystallite_orientation(1:4,1,neighboring_i,neighboring_e), & crystallite_orientation(1:4,1,neighboring_i,neighboring_e), &
crystallite_symmetryID(1,i,e)) ! calculate disorientation crystallite_symmetryID(1,i,e)) ! calculate disorientation
else ! for neighbor with different phase else ! for neighbor with different phase
@ -3485,8 +3485,8 @@ function crystallite_postResults(&
) )
!*** variables and functions from other modules ***! !*** variables and functions from other modules ***!
use math, only: math_QuaternionToEuler, & use math, only: math_qToEuler, &
math_QuaternionToAxisAngle, & math_qToAxisAngle, &
math_mul33x33, & math_mul33x33, &
math_transpose33, & math_transpose33, &
math_det33, & math_det33, &
@ -3549,26 +3549,26 @@ function crystallite_postResults(&
crystallite_postResults(c+1:c+mySize) = crystallite_orientation(1:4,g,i,e) ! grain orientation as quaternion crystallite_postResults(c+1:c+mySize) = crystallite_orientation(1:4,g,i,e) ! grain orientation as quaternion
case ('eulerangles') case ('eulerangles')
mySize = 3_pInt mySize = 3_pInt
crystallite_postResults(c+1:c+mySize) = inDeg * math_QuaternionToEuler(crystallite_orientation(1:4,g,i,e)) ! grain orientation as Euler angles in degree crystallite_postResults(c+1:c+mySize) = inDeg * math_qToEuler(crystallite_orientation(1:4,g,i,e)) ! grain orientation as Euler angles in degree
case ('grainrotation') case ('grainrotation')
mySize = 4_pInt mySize = 4_pInt
crystallite_postResults(c+1:c+mySize) = math_QuaternionToAxisAngle(crystallite_rotation(1:4,g,i,e)) ! grain rotation away from initial orientation as axis-angle in crystal reference coordinates crystallite_postResults(c+1:c+mySize) = math_qToAxisAngle(crystallite_rotation(1:4,g,i,e)) ! grain rotation away from initial orientation as axis-angle in crystal reference coordinates
crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree
case ('grainrotationx') case ('grainrotationx')
mySize = 1_pInt mySize = 1_pInt
rotation = math_QuaternionToAxisAngle(math_qMul(math_qMul(crystallite_orientation(1:4,g,i,e), & rotation = math_qToAxisAngle(math_qMul(math_qMul(crystallite_orientation(1:4,g,i,e), &
crystallite_rotation(1:4,g,i,e)), & crystallite_rotation(1:4,g,i,e)), &
math_qConj(crystallite_orientation(1:4,g,i,e)))) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates math_qConj(crystallite_orientation(1:4,g,i,e)))) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates
crystallite_postResults(c+1) = inDeg * rotation(1) * rotation(4) ! angle in degree crystallite_postResults(c+1) = inDeg * rotation(1) * rotation(4) ! angle in degree
case ('grainrotationy') case ('grainrotationy')
mySize = 1_pInt mySize = 1_pInt
rotation = math_QuaternionToAxisAngle(math_qMul(math_qMul(crystallite_orientation(1:4,g,i,e), & rotation = math_qToAxisAngle(math_qMul(math_qMul(crystallite_orientation(1:4,g,i,e), &
crystallite_rotation(1:4,g,i,e)), & crystallite_rotation(1:4,g,i,e)), &
math_qConj(crystallite_orientation(1:4,g,i,e)))) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates math_qConj(crystallite_orientation(1:4,g,i,e)))) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates
crystallite_postResults(c+1) = inDeg * rotation(2) * rotation(4) ! angle in degree crystallite_postResults(c+1) = inDeg * rotation(2) * rotation(4) ! angle in degree
case ('grainrotationz') case ('grainrotationz')
mySize = 1_pInt mySize = 1_pInt
rotation = math_QuaternionToAxisAngle(math_qMul(math_qMul(crystallite_orientation(1:4,g,i,e), & rotation = math_qToAxisAngle(math_qMul(math_qMul(crystallite_orientation(1:4,g,i,e), &
crystallite_rotation(1:4,g,i,e)), & crystallite_rotation(1:4,g,i,e)), &
math_qConj(crystallite_orientation(1:4,g,i,e)))) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates math_qConj(crystallite_orientation(1:4,g,i,e)))) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates
crystallite_postResults(c+1) = inDeg * rotation(3) * rotation(4) ! angle in degree crystallite_postResults(c+1) = inDeg * rotation(3) * rotation(4) ! angle in degree

120
code/damask.core.pyf
View File

@ -9,7 +9,7 @@
! The auto-generated file is quite heavily corrected ! The auto-generated file is quite heavily corrected
! For modifying, notice the following hints: ! For modifying, notice the following hints:
! - if the dimension of an array depend on a array that is itself an input, use the C-Syntax: (1) becomes [0] etc. ! - if the dimension of an array depend on a array that is itself an input, use the C-Syntax: (1) becomes [0] etc.
! - be sure that the precision defined for math.f90 is integer, real*8, and complex*16 ! - be sure that the precision defined is integer, real*8, and complex*16
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
python module core ! in python module core ! in
interface ! in :core interface ! in :core
@ -74,20 +74,16 @@ python module core ! in
real*8, dimension(size(field,1),size(field,2),size(field,3),size(field,4)), depend(field) :: math_divergenceFDM real*8, dimension(size(field,1),size(field,2),size(field,3),size(field,4)), depend(field) :: math_divergenceFDM
end function math_divergenceFDM end function math_divergenceFDM
subroutine math_nearestNeighborSearch(spatialDim,Favg,geomdim,queryPoints,domainPoints,querySet,domainSet,indices) ! in :math:math.f90 function math_periodicNearestNeighbor(geomdim,Favg,querySet,domainSet) ! in :math:math.f90
! input variables ! input variables
integer, intent(in) :: spatialDim real*8, dimension(3), intent(in) :: geomdim
real*8, dimension(3,3), intent(in) :: Favg real*8, dimension(3,3), intent(in) :: Favg
real*8, dimension(3), intent(in) :: geomdim real*8, dimension(:,:), intent(in) :: querySet
integer, intent(in) :: queryPoints real*8, dimension(:,:), intent(in) :: domainSet
integer, intent(in) :: domainPoints integer, dimension(size(querySet,2)), depend(querySet) :: math_periodicNearestNeighbor
real*8, dimension(spatialDim,queryPoints), intent(in), depend(spatialDim,queryPoints) :: querySet
real*8, dimension(spatialDim,domainPoints), intent(in), depend(spatialDim,domainPoints) :: domainSet
! output variable
integer, dimension(queryPoints), intent(out), depend(queryPoints) :: indices
! depending on input ! depending on input
real*8, dimension(spatialDim,(3**spatialDim)*domainPoints), depend(spatialDim,domainPoints) :: domainSetLarge real*8, dimension(size(domainSet,1),(3_pInt**size(domainSet,1))*size(domainSet,2)), depend(domainSet) :: domainSetLarge
end subroutine math_nearestNeighborSearch end function math_periodicNearestNeighbor
function math_tensorAvg(field) ! in :math:math.f90 function math_tensorAvg(field) ! in :math:math.f90
! input variables ! input variables
@ -100,14 +96,14 @@ python module core ! in
! input variables ! input variables
real*8, dimension(:,:,:,:,:), intent(in) :: F real*8, dimension(:,:,:,:,:), intent(in) :: F
! output variables ! output variables
real*8, dimension(size(F,1),size(F,2),size(F,3),3,3), depend(F) :: math_logstrainSpat real*8, dimension(3,3,size(F,3),size(F,4),size(F,5)), depend(F) :: math_logstrainSpat
end function math_logstrainSpat end function math_logstrainSpat
function math_logstrainMat(F) ! in :math:math.f90 function math_logstrainMat(F) ! in :math:math.f90
! input variables ! input variables
real*8, dimension(:,:,:,:,:), intent(in) :: F real*8, dimension(:,:,:,:,:), intent(in) :: F
! function ! function
real*8, dimension(size(F,1),size(F,2),size(F,3),3,3), depend(F) :: math_logstrainMat real*8, dimension(3,3,size(F,3),size(F,4),size(F,5)), depend(F) :: math_logstrainMat
end function math_logstrainMat end function math_logstrainMat
function math_cauchy(F,P) ! in :math:math.f90 function math_cauchy(F,P) ! in :math:math.f90
@ -115,7 +111,7 @@ python module core ! in
real*8, dimension(:,:,:,:,:), intent(in) :: F real*8, dimension(:,:,:,:,:), intent(in) :: F
real*8, dimension(:,:,:,:,:), intent(in) :: P real*8, dimension(:,:,:,:,:), intent(in) :: P
! function ! function
real*8, dimension(size(F,1),size(F,2),size(F,3),3,3), depend(F) :: math_cauchy real*8, dimension(3,3,size(F,3),size(F,4),size(F,5)), depend(F) :: math_cauchy
end function math_cauchy end function math_cauchy
end module math end module math
@ -126,11 +122,11 @@ python module core ! in
end module fesolving end module fesolving
module mesh ! in :mesh:mesh.f90 module mesh ! in :mesh:mesh.f90
subroutine mesh_init(ip,element) subroutine mesh_init(ip,element)
integer, parameter :: ip = 1 integer, parameter :: ip = 1
integer, parameter :: element = 1 integer, parameter :: element = 1
end subroutine mesh_init end subroutine mesh_init
function mesh_regrid(adaptive,resNewInput,minRes) ! in :mesh:mesh.f90 function mesh_regrid(adaptive,resNewInput,minRes) ! in :mesh:mesh.f90
logical, intent(in) :: adaptive logical, intent(in) :: adaptive
integer, dimension(3) :: mesh_regrid integer, dimension(3) :: mesh_regrid
@ -138,62 +134,44 @@ python module core ! in
integer, dimension(3), intent(in), optional :: minRes = -1 integer, dimension(3), intent(in), optional :: minRes = -1
end function mesh_regrid end function mesh_regrid
subroutine mesh_regular_grid(res,geomdim,defgrad_av,centroids,nodes) ! in :math:math.f90 function mesh_nodesAroundCentres(gDim,Favg,centres) ! in :mesh:mesh.f90
! input variables real*8, dimension(:,:,:,:), intent(in) :: centres
integer, dimension(3), intent(in) :: res real*8, dimension(3), intent(in) :: gDim
real*8, dimension(3), intent(in) :: geomdim real*8, dimension(3,3), intent(in) :: Favg
real*8, dimension(3,3), intent(in) :: defgrad_av real*8, dimension(3,size(centres,2)+1,size(centres,3)+1,size(centres,4)+1), depend(centres) :: mesh_nodesAroundCentres
real*8, dimension(res[0], res[1], res[2], 3), intent(in), depend(res[0],res[1],res[2]) :: centroids real*8, dimension(3,size(centres,2)+1,size(centres,3)+1,size(centres,4)+1), depend(centres) :: wrappedCentres
! output variables end function mesh_nodesAroundCentres
real*8, dimension(res[0]+1,res[1]+1,res[2]+1,3), intent(out), depend(res[0],res[1],res[2]) :: nodes
! variables with dimension depending on input
real*8, dimension(res[0]+2,res[1]+2,res[2]+2,3), depend(res[0],res[1],res[2]) :: wrappedCentroids
end subroutine mesh_regular_grid
subroutine deformed_linear(res,geomdim,defgrad_av,defgrad,coord_avgCorner) ! in :math:math.f90 function mesh_deformedCoordsLinear(gDim,F,FavgIn) ! in :mesh:mesh.f90
! input variables real*8, dimension(:,:,:,:,:), intent(in) :: F
integer, dimension(3), intent(in) :: res real*8, dimension(3), intent(in) :: gDim
real*8, dimension(3), intent(in) :: geomdim real*8, dimension(3,3), intent(in), optional :: FavgIn = -1.0
real*8, dimension(3,3), intent(in) :: defgrad_av real*8, dimension(3,size(F,3),size(F,4),size(F,5)), depend(F) :: mesh_deformedCoordsLinear
real*8, dimension(res[0], res[1], res[2], 3,3),intent(in), depend(res[0],res[1],res[2]) :: defgrad real*8, dimension(3,0:size(F,3)-1,0:size(F,4)-1,0:size(F,5)-1,0:7), depend(F) :: coordsAvgOrder
! output variables end function mesh_deformedCoordsLinear
real*8, dimension(res[0], res[1], res[2], 3), intent(out), depend(res[0],res[1],res[2]) :: coord_avgCorner
! variables with dimension depending on input
real*8, dimension(8,6,res[0],res[1],res[2],3), depend(res[0],res[1],res[2]) :: coord
real*8, dimension(8,res[0],res[1],res[2],3), depend(res[0],res[1],res[2]) :: coord_avgOrder
end subroutine deformed_linear
subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords) ! in :math:math.f90 function mesh_deformedCoordsFFT(gDim,F,scalingIn,FavgIn) ! in :mesh:mesh.f90
! input variables real*8, dimension(:,:,:,:,:), intent(in) :: F
integer, dimension(3), intent(in) :: res real*8, dimension(3), intent(in) :: gDim
real*8, dimension(3), intent(in) :: geomdim real*8, intent(in), optional, :: scalingIn = -1.0
real*8, dimension(3,3), intent(in) :: defgrad_av real*8, dimension(3,3), intent(in), optional :: FavgIn = -1.0
real*8, intent(in) :: scaling real*8, dimension(3,size(F,3),size(F,4),size(F,5)), depend(F) :: mesh_deformedCoordsFFT
real*8, dimension(res[0], res[1], res[2], 3,3),intent(in), depend(res[0],res[1],res[2]) :: defgrad end function mesh_deformedCoordsFFT
! output variables
real*8, dimension(res[0], res[1], res[2], 3), intent(out), depend(res[0],res[1],res[2]) :: coords
end subroutine deformed_fft
subroutine volume_compare(res,geomdim,defgrad,nodes,volume_mismatch) ! in :math:math.f90 function mesh_volumeMismatch(gDim,F,nodes) ! in :mesh:mesh.f90
! input variables real*8, dimension(:,:,:,:,:), intent(in) :: F
integer, dimension(3), intent(in) :: res real*8, dimension(:,:,:,:), intent(in) :: nodes
real*8, dimension(3), intent(in) :: geomdim real*8, dimension(3), intent(in) :: gDim
real*8, dimension(res[0], res[1], res[2], 3,3),intent(in), depend(res[0],res[1],res[2]) :: defgrad real*8, dimension(size(F,3),size(F,4),size(F,5)), depend(F) :: mesh_volumeMismatch
real*8, dimension(res[0]+1,res[1]+1,res[2]+1,3), intent(in), depend(res[0],res[1],res[2]) :: nodes end function mesh_volumeMismatch
! output variables
real*8, dimension(res[0], res[1], res[2]), intent(out), depend(res[0],res[1],res[2])) :: volume_mismatch
end subroutine volume_compare
subroutine shape_compare(res,geomdim,defgrad,nodes,centroids,shape_mismatch) ! in :math:math.f90 function mesh_shapeMismatch(gDim,F,nodes,centres) ! in :mesh:mesh.f90
! input variables real*8, dimension(:,:,:,:,:), intent(in) :: F
integer, dimension(3), intent(in) :: res real*8, dimension(:,:,:,:), intent(in) :: nodes
real*8, dimension(3), intent(in) :: geomdim real*8, dimension(:,:,:,:), intent(in) :: centres
real*8, dimension(res[0], res[1], res[2], 3,3),intent(in), depend(res[0],res[1],res[2]) :: defgrad real*8, dimension(3), intent(in) :: gDim
real*8, dimension(res[0]+1,res[1]+1,res[2]+1,3), intent(in), depend(res[0],res[1],res[2]) :: nodes real*8, dimension(size(F,3),size(F,4),size(F,5)), depend(F) :: mesh_shapeMismatch
real*8, dimension(res[0], res[1], res[2], 3), intent(in), depend(res[0],res[1],res[2]) :: centroids end function mesh_shapeMismatch
! output variables
real*8, dimension(res[0], res[1], res[2]), intent(out), depend(res[0],res[1],res[2])) :: shape_mismatch
end subroutine shape_compare
end module mesh end module mesh
end interface end interface

1309
code/math.f90
View File

File diff suppressed because it is too large Load Diff

1082
code/mesh.f90
View File

File diff suppressed because it is too large Load Diff

14
lib/damask/__init__.py
View File

@ -10,7 +10,6 @@ from .result import Result # one class with subclasses
from .geometry import Geometry # one class with subclasses from .geometry import Geometry # one class with subclasses
from .solver import Solver # one class with subclasses from .solver import Solver # one class with subclasses
from .test import Test from .test import Test
from util import *
try: try:
from . import core from . import core
@ -30,7 +29,7 @@ try:
core.math.curlFFT = core.math.math_curlFFT core.math.curlFFT = core.math.math_curlFFT
core.math.divergenceFFT = core.math.math_divergenceFFT core.math.divergenceFFT = core.math.math_divergenceFFT
core.math.divergenceFDM = core.math.math_divergenceFDM core.math.divergenceFDM = core.math.math_divergenceFDM
#core.math.periodicNearestNeighbor = core.math.math_periodicNearestNeighbor core.math.periodicNearestNeighbor = core.math.math_periodicNearestNeighbor
core.math.tensorAvg = core.math.math_tensorAvg core.math.tensorAvg = core.math.math_tensorAvg
core.math.logstrainSpat = core.math.math_logstrainSpat core.math.logstrainSpat = core.math.math_logstrainSpat
core.math.logstrainMat = core.math.math_logstrainMat core.math.logstrainMat = core.math.math_logstrainMat
@ -38,11 +37,12 @@ try:
core.FEsolving.init = core.FEsolving.FE_init core.FEsolving.init = core.FEsolving.FE_init
core.mesh.init = core.mesh.mesh_init core.mesh.init = core.mesh.mesh_init
core.mesh.regrid = core.mesh.mesh_regrid core.mesh.regrid = core.mesh.mesh_regrid
#core.mesh.volumeMismatch = core.mesh.mesh_volumeMismatch core.mesh.nodesAroundCentres = core.mesh.mesh_nodesAroundCentres
#core.mesh.shapeMismatch = core.mesh.mesh_shapeMismatch core.mesh.deformedCoordsLinear = core.mesh.mesh_deformedCoordsLinear
#core.mesh.nodesAroundCentroids = core.mesh.mesh_spectral_nodesAroundCentroids core.mesh.deformedCoordsFFT = core.mesh.mesh_deformedCoordsFFT
#core.mesh.deformedCoordsLin = core.mesh.mesh_deformedCoordsLin core.mesh.volumeMismatch = core.mesh.mesh_volumeMismatch
#core.mesh.deformedCoordsFFT = core.mesh.mesh_deformedCoordsFFT core.mesh.shapeMismatch = core.mesh.mesh_shapeMismatch
except ImportError, e: except ImportError, e:
core = None # from http://www.python.org/dev/peps/pep-0008/ core = None # from http://www.python.org/dev/peps/pep-0008/
if('setup_processing' not in sys.argv[0]): if('setup_processing' not in sys.argv[0]):

67
processing/post/3Dvisualize.py
View File

@ -132,7 +132,7 @@ def vtk_writeASCII_mesh(mesh,data,res,sep):
'ASCII', 'ASCII',
'DATASET UNSTRUCTURED_GRID', 'DATASET UNSTRUCTURED_GRID',
'POINTS %i float'%N1, 'POINTS %i float'%N1,
[[['\t'.join(map(str,mesh[i,j,k])) for i in range(res[0]+1)] for j in range(res[1]+1)] for k in range(res[2]+1)], [[['\t'.join(map(str,mesh[:,i,j,k])) for i in range(res[0]+1)] for j in range(res[1]+1)] for k in range(res[2]+1)],
'CELLS %i %i'%(N,N*9), 'CELLS %i %i'%(N,N*9),
] ]
@ -189,7 +189,7 @@ def gmsh_writeASCII_mesh(mesh,data,res,sep):
'$EndMeshFormat', '$EndMeshFormat',
'$Nodes', '$Nodes',
'%i float'%N1, '%i float'%N1,
[[['\t'.join(map(str,l,mesh[i,j,k])) for l in range(1,N1+1) for i in range(res[0]+1)] for j in range(res[1]+1)] for k in range(res[2]+1)], [[['\t'.join(map(str,l,mesh[:,i,j,k])) for l in range(1,N1+1) for i in range(res[0]+1)] for j in range(res[1]+1)] for k in range(res[2]+1)],
'$EndNodes', '$EndNodes',
'$Elements', '$Elements',
'%i'%N, '%i'%N,
@ -346,7 +346,7 @@ 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('--filenodalcoords', dest='filenodalcoords', type='string', \ parser.add_option('--filenodalcoords', dest='filenodalcoords', type='string', \
help='ASCII table containing nodal coords') help='ASCII table containing nodal coords')
parser.add_option('--labelnodalcoords', dest='nodalcoords', type='string', \ parser.add_option('--labelnodalcoords', dest='nodalcoords', type='string', nargs=3, \
help='labels of nodal coords in ASCII table') help='labels of nodal coords in ASCII table')
parser.add_option('-l', '--linear', dest='linearreconstruction', action='store_true',\ parser.add_option('-l', '--linear', dest='linearreconstruction', action='store_true',\
help='use linear reconstruction of geometry [%default]') help='use linear reconstruction of geometry [%default]')
@ -367,14 +367,14 @@ 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(filenodalcoords = '') parser.set_defaults(filenodalcoords = '')
parser.set_defaults(labelnodalcoords = 'coord') parser.set_defaults(labelnodalcoords = ['coord.x','coord.y','coord.z'])
parser.set_defaults(linearreconstruction = False) parser.set_defaults(linearreconstruction = False)
sep = {'n': '\n', 't': '\t', 's': ' '} sep = {'n': '\n', 't': '\t', 's': ' '}
(options, args) = parser.parse_args() (options, args) = parser.parse_args()
if options.scaling !=1.0 and options.linearreconstruction: print 'cannot scale for linear reconstruction' if options.scaling !=1.0 and options.linearreconstruction: print 'cannot scale for linear reconstruction'
if options.scaling !=1.0 and options.filenodalcoords!='': print 'cannot scale when reading coordinate from file' if options.scaling !=1.0 and options.filenodalcoords!='': print 'cannot scale when reading coordinate from file'
options.separator = options.separator.lower() options.separator = options.separator.lower()
for filename in args: for filename in args:
if not os.path.exists(filename): if not os.path.exists(filename):
@ -440,6 +440,7 @@ for filename in args:
key=lambda x:(x[locol+0],x[locol+1],x[locol+2])),'d') # sort with z as fastest and x as slowest index key=lambda x:(x[locol+0],x[locol+1],x[locol+2])),'d') # sort with z as fastest and x as slowest index
N = len(values) N = len(values)
grid = [{},{},{}] grid = [{},{},{}]
for j in xrange(3): for j in xrange(3):
for i in xrange(N): for i in xrange(N):
@ -448,7 +449,9 @@ for filename in args:
res = numpy.array([len(grid[0]),\ res = numpy.array([len(grid[0]),\
len(grid[1]),\ len(grid[1]),\
len(grid[2]),],'i') len(grid[2]),],'i')
if (res[0]%2 != 0 or res[1]%2 != 0 or (res[2] != 1 and res[2]%2 !=0)):
print 'using linear reconstruction for uneven resolution'
options.linearreconstruction = True
dim = numpy.ones(3) dim = numpy.ones(3)
for i,r in enumerate(res): for i,r in enumerate(res):
@ -460,34 +463,42 @@ for filename in args:
else: else:
dim[2] = options.unitlength dim[2] = options.unitlength
if options.filenodalcoords: if options.undeformed:
Favg = numpy.eye(3)
else:
Favg = damask.core.math.tensorAvg(
numpy.reshape(numpy.transpose(values[:,column['tensor'][options.defgrad]:
column['tensor'][options.defgrad]+9]),
(3,3,res[0],res[1],res[2])))
if not options.filenodalcoords:
F = numpy.reshape(numpy.transpose(values[:,column['tensor'][options.defgrad]:
column['tensor'][options.defgrad]+9]),
(3,3,res[0],res[1],res[2]))
if options.linearreconstruction:
centroids = damask.core.mesh.deformedCoordsLinear(dim,F,Favg)
else:
centroids = damask.core.mesh.deformedCoordsFFT(dim,F,options.scaling,Favg)
mesh = damask.core.mesh.nodesAroundCentres(dim,Favg,centroids)
else:
mesh = numpy.zeros(((res[0]+1)*(res[1]+1)*(res[2]+1),3),'d') mesh = numpy.zeros(((res[0]+1)*(res[1]+1)*(res[2]+1),3),'d')
mesh=mesh.reshape(res[0]+1,res[1]+1,res[2]+1,3)
filenodalcoords = open(options.filenodalcoords) filenodalcoords = open(options.filenodalcoords)
tablenodalcoords = damask.ASCIItable(filenodalcoords) tablenodalcoords = damask.ASCIItable(filenodalcoords)
tablenodalcoords.head_read() tablenodalcoords.head_read()
coord = tablenodalcoords.labels.index(options.labelnodalcoords+'.x') columns = [tablenodalcoords.labels.index(options.labelnodalcoords[0]),
tablenodalcoords.labels.index(options.labelnodalcoords[1]),
tablenodalcoords.labels.index(options.labelnodalcoords[2])]
i = 0 i = 0
while tablenodalcoords.data_read(): while tablenodalcoords.data_read(): # read next data line of ASCII table
mesh[i%(res[0]+1),(i//(res[0]+1))%(res[1]+1),(i//(res[0]+1)//(res[1]+1)) % (res[2]+1),:]=\
[float(tablenodalcoords.data[coord]),
float(tablenodalcoords.data[coord+1]), mesh[i,:]=float(tablenodalcoords.data[column[:]])
float(tablenodalcoords.data[coord+2])]
i += 1
i += 1
mesh=mesh.reshape(res[0]+1,res[1]+1,res[2]+1,3)
else:
F = numpy.reshape(values[:,column['tensor'][options.defgrad]: column['tensor'][options.defgrad]+9],(res[0],res[1],res[2],3,3))
if options.undeformed:
Favg = numpy.eye(3)
else:
Favg = damask.core.math.tensorAvg(F)
if options.linearreconstruction:
centroids = damask.core.mesh.deformed_linear(res,dim,Favg,F)
else:
centroids = damask.core.mesh.deformed_fft(res,dim,Favg,options.scaling,F)
mesh = damask.core.mesh.mesh_regular_grid(res,dim,Favg,centroids)
fields = {\ fields = {\
'tensor': {},\ 'tensor': {},\
'vector': {},\ 'vector': {},\

43
processing/post/addCompatibilityMismatch.py
View File

@ -50,28 +50,23 @@ parser.add_option('--no-volume','-v', dest='noVolume', action='store_false', \
help='do not calculate volume mismatch [%default]') help='do not calculate volume mismatch [%default]')
parser.add_option('-c','--coordinates', dest='coords', type='string',\ parser.add_option('-c','--coordinates', dest='coords', type='string',\
help='column heading for coordinates [%default]') help='column heading for coordinates [%default]')
parser.add_option('-f','--deformation', dest='defgrad', action='extend', type='string', \ parser.add_option('-f','--deformation', dest='F', action='extend', type='string', \
help='heading(s) of columns containing deformation tensor values %default') help='heading(s) of columns containing deformation tensor values %default')
parser.set_defaults(noVolume = False) parser.set_defaults(noVolume = False)
parser.set_defaults(noShape = False) parser.set_defaults(noShape = False)
parser.set_defaults(coords = 'ip') parser.set_defaults(coords = 'ip')
parser.set_defaults(defgrad = 'f') parser.set_defaults(F = 'f')
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()
defgrad_av = {}
centroids = {}
nodes = {}
shape_mismatch = {}
volume_mismatch= {}
datainfo = { # list of requested labels per datatype datainfo = { # list of requested labels per datatype
'defgrad': {'len':9, 'F': {'len':9,
'label':[]}, 'label':[]},
} }
if options.defgrad != None: datainfo['defgrad']['label'] += options.defgrad if options.F != None: datainfo['F']['label'] += options.F
# ------------------------------------------ setup file handles --------------------------------------- # ------------------------------------------ setup file handles ---------------------------------------
@ -121,13 +116,13 @@ for file in files:
# --------------- figure out columns to process # --------------- figure out columns to process
key = '1_%s' %options.defgrad key = '1_%s' %options.F
if key not in table.labels: if key not in table.labels:
sys.stderr.write('column %s not found...\n'%key) sys.stderr.write('column %s not found...\n'%key)
else: else:
defgrad = numpy.array([0.0 for i in xrange(N*9)]).reshape(list(res)+[3,3]) F = numpy.array([0.0 for i in xrange(N*9)]).reshape([3,3]+list(res))
if not options.noShape: table.labels_append(['mismatch_shape(%s)' %options.defgrad]) if not options.noShape: table.labels_append(['shapeMismatch(%s)' %options.F])
if not options.noVolume: table.labels_append(['mismatch_volume(%s)'%options.defgrad]) if not options.noVolume: table.labels_append(['volMismatch(%s)'%options.F])
column = table.labels.index(key) column = table.labels.index(key)
# ------------------------------------------ assemble header --------------------------------------- # ------------------------------------------ assemble header ---------------------------------------
@ -142,13 +137,19 @@ for file in files:
while table.data_read(): # read next data line of ASCII table while table.data_read(): # read next data line of ASCII table
(x,y,z) = location(idx,res) # figure out (x,y,z) position from line count (x,y,z) = location(idx,res) # figure out (x,y,z) position from line count
idx += 1 idx += 1
defgrad[x,y,z] = numpy.array(map(float,table.data[column:column+9]),'d').reshape(3,3) F[0:3,0:3,x,y,z] = numpy.array(map(float,table.data[column:column+9]),'d').reshape(3,3)
defgrad_av = damask.core.math.tensorAvg(defgrad) Favg = damask.core.math.tensorAvg(F)
centroids = damask.core.mesh.deformed_fft(res,geomdim,defgrad_av,1.0,defgrad)
nodes = damask.core.mesh.mesh_regular_grid(res,geomdim,defgrad_av,centroids) if (res[0]%2 != 0 or res[1]%2 != 0 or (res[2] != 1 and res[2]%2 !=0)):
if not options.noShape: shape_mismatch = damask.core.mesh.shape_compare( res,geomdim,defgrad,nodes,centroids) print 'using linear reconstruction for uneven resolution'
if not options.noVolume: volume_mismatch = damask.core.mesh.volume_compare(res,geomdim,defgrad,nodes) centres = damask.core.mesh.deformedCoordsLin(geomdim,F,Favg)
else:
centres = damask.core.mesh.deformedCoordsFFT(geomdim,F,1.0,Favg)
nodes = damask.core.mesh.nodesAroundCentres(geomdim,Favg,centres)
if not options.noShape: shapeMismatch = damask.core.mesh.shapeMismatch( geomdim,F,nodes,centres)
if not options.noVolume: volumeMismatch = damask.core.mesh.volumeMismatch(geomdim,F,nodes)
# ------------------------------------------ process data --------------------------------------- # ------------------------------------------ process data ---------------------------------------
@ -157,8 +158,8 @@ for file in files:
while table.data_read(): # read next data line of ASCII table while table.data_read(): # read next data line of ASCII table
(x,y,z) = location(idx,res) # figure out (x,y,z) position from line count (x,y,z) = location(idx,res) # figure out (x,y,z) position from line count
idx += 1 idx += 1
if not options.noShape: table.data_append(shape_mismatch[x,y,z]) if not options.noShape: table.data_append( shapeMismatch[x,y,z])
if not options.noVolume: table.data_append(volume_mismatch[x,y,z]) if not options.noVolume: table.data_append(volumeMismatch[x,y,z])
table.data_write() # output processed line table.data_write() # output processed line

18
processing/post/addDeformedConfiguration.py
View File

@ -117,8 +117,8 @@ for file in files:
if key not in table.labels: if key not in table.labels:
sys.stderr.write('column %s not found...\n'%key) sys.stderr.write('column %s not found...\n'%key)
else: else:
defgrad = numpy.array([0.0 for i in xrange(N*9)]).reshape(list(res)+[3,3]) F = numpy.array([0.0 for i in xrange(N*9)]).reshape([3,3]+list(res))
table.labels_append(['%s_coords'%(coord+1) for coord in xrange(3)]) # extend ASCII header with new labels table.labels_append(['%s_coordsMod'%(coord+1) for coord in xrange(3)]) # extend ASCII header with new labels
column = table.labels.index(key) column = table.labels.index(key)
@ -134,14 +134,16 @@ for file in files:
while table.data_read(): # read next data line of ASCII table while table.data_read(): # read next data line of ASCII table
(x,y,z) = location(idx,res) # figure out (x,y,z) position from line count (x,y,z) = location(idx,res) # figure out (x,y,z) position from line count
idx += 1 idx += 1
defgrad[x,y,z] = numpy.array(map(float,table.data[column:column+9]),'d').reshape(3,3) F[:,:,x,y,z] = numpy.array(map(float,table.data[column:column+9]),'d').reshape(3,3)
# ------------------------------------------ process value field ---------------------------- # ------------------------------------------ process value field ----------------------------
defgrad_av = damask.core.math.tensorAvg(defgrad) Favg = damask.core.math.tensorAvg(F)
if options.linearreconstruction: if options.linearreconstruction:
centroids = damask.core.mesh.deformed_fft(res,geomdim,defgrad_av,1.0,defgrad) centroids = damask.core.mesh.deformedCoordsLin(geomdim,F,Favg)
else: else:
centroids = damask.core.mesh.deformed_fft(res,geomdim,defgrad_av,1.0,defgrad) print 'ddd'
centroids = damask.core.mesh.deformedCoordsFFT(geomdim,F,1.0,Favg)
# ------------------------------------------ process data --------------------------------------- # ------------------------------------------ process data ---------------------------------------
table.data_rewind() table.data_rewind()
@ -149,7 +151,7 @@ for file in files:
while table.data_read(): # read next data line of ASCII table while table.data_read(): # read next data line of ASCII table
(x,y,z) = location(idx,res) # figure out (x,y,z) position from line count (x,y,z) = location(idx,res) # figure out (x,y,z) position from line count
idx += 1 idx += 1
table.data_append(list(centroids[x,y,z])) table.data_append(list(centroids[:,x,y,z]))
table.data_write() # output processed line table.data_write() # output processed line

9
processing/pre/geom_fromVoronoiTessellation.py
View File

@ -206,13 +206,10 @@ for file in files:
/float(info['resolution'][2]) /float(info['resolution'][2])
undeformed[:,i] += shift undeformed[:,i] += shift
indices = damask.core.math.math_nearestNeighborSearch(3,\ indices = damask.core.math.periodicNearestNeighbor(\
numpy.array(([1.0,0.0,0.0],\
[0.0,1.0,0.0],\
[0.0,0.0,1.0]),'d'),\
info['dimension'],\ info['dimension'],\
N,info['grains'],undeformed,coords)//3**3 + 1 # floor division to kill periodic images numpy.eye(3),\
undeformed,coords)//3**3 + 1 # floor division to kill periodic images
for n in xrange(info['resolution'][1:3].prod()): # loop over 2nd and 3rd dimension for n in xrange(info['resolution'][1:3].prod()): # loop over 2nd and 3rd dimension
file['output'].write({ True: ' ', file['output'].write({ True: ' ',
False:'\n'}[options.twoD].\ False:'\n'}[options.twoD].\