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

Merge branch 'NoCoreModule' into development

This commit is contained in:
Martin Diehl 2016-07-01 15:07:10 +02:00
parent fa13aa30b0 7445aa0b32
commit af4307ea72
24 changed files with 380 additions and 12290 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
DAMASK_env.bat Normal file
View File

@ -0,0 +1,19 @@
:: sets up an environment for DAMASK on Windows
:: usage: call DAMASK_env.bat
@echo off
set LOCATION=%~dp0
set DAMASK_ROOT=%LOCATION%\DAMASK
set DAMASK_NUM_THREADS=2
chcp 1252
Title Düsseldorf Advanced Materials Simulation Kit - DAMASK, MPIE Düsseldorf
echo.
echo Düsseldorf Advanced Materials Simulation Kit - DAMASK
echo Max-Planck-Institut für Eisenforschung, Düsseldorf
echo http://damask.mpie.de
echo.
echo Preparing environment ...
echo DAMASK_ROOT=%DAMASK_ROOT%
echo DAMASK_NUM_THREADS=%DAMASK_NUM_THREADS%
set DAMASK_BIN=%DAMASK_ROOT%\bin
set PATH=%PATH%;%DAMASK_BIN%
set PYTHONPATH=%PYTHONPATH%;%DAMASK_ROOT%\lib

2
Makefile
View File

@ -21,11 +21,11 @@ marc:
processing: processing:
@if hash cython 2>/dev/null; then \ @if hash cython 2>/dev/null; then \
cd ./lib/damask; \ cd ./lib/damask; \
ln -s orientation.py corientation.pyx; \
CC=gcc python setup_corientation.py build_ext --inplace; \ CC=gcc python setup_corientation.py build_ext --inplace; \
rm -rv build; \ rm -rv build; \
rm *.c; \ rm *.c; \
fi fi
@./installation/compile_CoreModule.py ${MAKEFLAGS}
.PHONY: tidy .PHONY: tidy
tidy: tidy:

13
code/DAMASK_spectral.f90
View File

@ -617,12 +617,12 @@ program DAMASK_spectral
timeinc = timeinc/2.0_pReal timeinc = timeinc/2.0_pReal
elseif (solres(1)%termIll) then ! material point model cannot find a solution, exit in any casy elseif (solres(1)%termIll) then ! material point model cannot find a solution, exit in any casy
call IO_warning(850_pInt) call IO_warning(850_pInt)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written (e.g. for regridding) call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written
elseif (continueCalculation == 1_pInt) then elseif (continueCalculation == 1_pInt) then
guess = .true. ! accept non converged BVP solution guess = .true. ! accept non converged BVP solution
else ! default behavior, exit if spectral solver does not converge else ! default behavior, exit if spectral solver does not converge
call IO_warning(850_pInt) call IO_warning(850_pInt)
call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written (e.g. for regridding) call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! quit and provide information about last restart inc written
endif endif
else else
guess = .true. ! start guessing after first converged (sub)inc guess = .true. ! start guessing after first converged (sub)inc
@ -722,20 +722,17 @@ end program DAMASK_spectral
!> @brief quit subroutine to mimic behavior of FEM solvers !> @brief quit subroutine to mimic behavior of FEM solvers
!> @details exits the Spectral solver and reports time and duration. Exit code 0 signals !> @details exits the Spectral solver and reports time and duration. Exit code 0 signals
!> everything went fine. Exit code 1 signals an error, message according to IO_error. Exit code !> everything went fine. Exit code 1 signals an error, message according to IO_error. Exit code
!> 2 signals request for regridding, increment of last saved restart information is written to !> 2 signals no converged solution and increment of last saved restart information is written to
!> stderr. Exit code 3 signals no severe problems, but some increments did not converge !> stderr. Exit code 3 signals no severe problems, but some increments did not converge
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine quit(stop_id) subroutine quit(stop_id)
use prec, only: & use prec, only: &
pInt pInt
use numerics, only: &
worldrank
implicit none implicit none
integer(pInt), intent(in) :: stop_id integer(pInt), intent(in) :: stop_id
integer, dimension(8) :: dateAndTime ! type default integer integer, dimension(8) :: dateAndTime ! type default integer
if (worldrank == 0_pInt) then
call date_and_time(values = dateAndTime) call date_and_time(values = dateAndTime)
write(6,'(/,a)') 'DAMASK terminated on:' write(6,'(/,a)') 'DAMASK terminated on:'
write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',& write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',&
@ -744,11 +741,9 @@ subroutine quit(stop_id)
write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',& write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',&
dateAndTime(6),':',& dateAndTime(6),':',&
dateAndTime(7) dateAndTime(7)
endif
if (stop_id == 0_pInt) stop 0 ! normal termination if (stop_id == 0_pInt) stop 0 ! normal termination
if (stop_id < 0_pInt) then ! trigger regridding if (stop_id < 0_pInt) then ! terminally ill, restart might help
if (worldrank == 0_pInt) &
write(0,'(a,i6)') 'restart information available at ', stop_id*(-1_pInt) write(0,'(a,i6)') 'restart information available at ', stop_id*(-1_pInt)
stop 2 stop 2
endif endif

6
code/FEsolving.f90
View File

@ -72,11 +72,9 @@ subroutine FE_init
integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt), allocatable, dimension(:) :: chunkPos
#endif #endif
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- FEsolving init -+>>>' write(6,'(/,a)') ' <<<+- FEsolving init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess
modelName = getSolverJobName() modelName = getSolverJobName()
#ifdef Spectral #ifdef Spectral
@ -153,10 +151,6 @@ subroutine FE_init
200 close(FILEUNIT) 200 close(FILEUNIT)
endif endif
!--------------------------------------------------------------------------------------------------
! the following array are allocated by mesh.f90 and need to be deallocated in case of regridding
if (allocated(calcMode)) deallocate(calcMode)
if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP)
#endif #endif
if (iand(debug_level(debug_FEsolving),debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_FEsolving),debug_levelBasic) /= 0_pInt) then
write(6,'(a21,l1)') ' restart writing: ', restartWrite write(6,'(a21,l1)') ' restart writing: ', restartWrite

2
code/IO.f90
View File

@ -1625,8 +1625,6 @@ subroutine IO_error(error_ID,el,ip,g,ext_msg)
msg = 'update of gamma operator not possible when pre-calculated' msg = 'update of gamma operator not possible when pre-calculated'
case (880_pInt) case (880_pInt)
msg = 'mismatch of microstructure count and a*b*c in geom file' msg = 'mismatch of microstructure count and a*b*c in geom file'
case (890_pInt)
msg = 'invalid input for regridding'
case (891_pInt) case (891_pInt)
msg = 'unknown solver type selected' msg = 'unknown solver type selected'
case (892_pInt) case (892_pInt)

7
code/Makefile
View File

@ -351,7 +351,7 @@ DAMASK_spectral.o: INTERFACENAME := spectral_interface.f90
SPECTRAL_SOLVER_FILES = spectral_mech_AL.o spectral_mech_Basic.o spectral_mech_Polarisation.o \ SPECTRAL_SOLVER_FILES = spectral_mech_AL.o spectral_mech_Basic.o spectral_mech_Polarisation.o \
spectral_thermal.o spectral_damage.o spectral_thermal.o spectral_damage.o
SPECTRAL_FILES = C_routines.o system_routines.o prec.o DAMASK_interface.o IO.o libs.o numerics.o debug.o math.o \ SPECTRAL_FILES = C_routines.o system_routines.o prec.o DAMASK_interface.o IO.o numerics.o debug.o math.o \
FEsolving.o mesh.o material.o lattice.o \ FEsolving.o mesh.o material.o lattice.o \
$(SOURCE_FILES) $(KINEMATICS_FILES) $(PLASTIC_FILES) constitutive.o \ $(SOURCE_FILES) $(KINEMATICS_FILES) $(PLASTIC_FILES) constitutive.o \
crystallite.o \ crystallite.o \
@ -401,7 +401,7 @@ DAMASK_FEM.exe: INCLUDE_DIRS += -I./
FEM_SOLVER_FILES = FEM_mech.o FEM_thermal.o FEM_damage.o FEM_vacancyflux.o FEM_porosity.o FEM_hydrogenflux.o FEM_SOLVER_FILES = FEM_mech.o FEM_thermal.o FEM_damage.o FEM_vacancyflux.o FEM_porosity.o FEM_hydrogenflux.o
FEM_FILES = prec.o DAMASK_interface.o FEZoo.o IO.o libs.o numerics.o debug.o math.o \ FEM_FILES = prec.o DAMASK_interface.o FEZoo.o IO.o numerics.o debug.o math.o \
FEsolving.o mesh.o material.o lattice.o \ FEsolving.o mesh.o material.o lattice.o \
$(SOURCE_FILES) $(KINEMATICS_FILES) $(PLASTIC_FILES) constitutive.o \ $(SOURCE_FILES) $(KINEMATICS_FILES) $(PLASTIC_FILES) constitutive.o \
crystallite.o \ crystallite.o \
@ -612,9 +612,6 @@ debug.o: debug.f90 \
numerics.o numerics.o
numerics.o: numerics.f90 \ numerics.o: numerics.f90 \
libs.o
libs.o: libs.f90 \
IO.o IO.o
IO.o: IO.f90 \ IO.o: IO.f90 \

1
code/commercialFEM_fileList.f90
View File

@ -4,7 +4,6 @@
!> @details List of files needed by MSC.Marc, Abaqus/Explicit, and Abaqus/Standard !> @details List of files needed by MSC.Marc, Abaqus/Explicit, and Abaqus/Standard
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
#include "IO.f90" #include "IO.f90"
#include "libs.f90"
#include "numerics.f90" #include "numerics.f90"
#include "debug.f90" #include "debug.f90"
#include "math.f90" #include "math.f90"

12
code/core_quit.f90
View File

@ -1,12 +0,0 @@
!********************************************************************
! quit subroutine to satisfy IO_error for core module
!
!********************************************************************
subroutine quit(stop_id)
use prec, only: &
pInt
implicit none
integer(pInt), intent(in) :: stop_id
end subroutine

12
code/libs.f90
View File

@ -1,12 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief dummy source for inclusion of Library files
!--------------------------------------------------------------------------------------------------
module libs
!nothing in here
end module libs
#include "../lib/IR_Precision.f90"
#include "../lib/Lib_Base64.f90"
#include "../lib/Lib_VTK_IO.f90"

19
code/math.f90
View File

@ -72,10 +72,6 @@ module math
3_pInt,3_pInt & 3_pInt,3_pInt &
],[2,9]) !< arrangement in Plain notation ],[2,9]) !< arrangement in Plain notation
#ifdef Spectral
include 'fftw3.f03'
#endif
public :: & public :: &
math_init, & math_init, &
math_qsort, & math_qsort, &
@ -163,21 +159,6 @@ module math
math_rotate_forward33, & math_rotate_forward33, &
math_rotate_backward33, & math_rotate_backward33, &
math_rotate_forward3333 math_rotate_forward3333
#ifdef Spectral
public :: &
fftw_set_timelimit, &
fftw_plan_dft_3d, &
fftw_plan_many_dft_r2c, &
fftw_plan_many_dft_c2r, &
fftw_plan_with_nthreads, &
fftw_init_threads, &
fftw_alloc_complex, &
fftw_execute_dft, &
fftw_execute_dft_r2c, &
fftw_execute_dft_c2r, &
fftw_destroy_plan, &
math_tensorAvg
#endif
private :: & private :: &
math_partition, & math_partition, &
halton, & halton, &

591
code/mesh.f90
View File

@ -116,12 +116,8 @@ module mesh
#endif #endif
#ifdef Spectral #ifdef Spectral
#ifdef PETSc
#include <petsc/finclude/petscsys.h> #include <petsc/finclude/petscsys.h>
include 'fftw3-mpi.f03' include 'fftw3-mpi.f03'
#else
include 'fftw3.f03'
#endif
#endif #endif
! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS) ! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS)
@ -413,18 +409,13 @@ module mesh
mesh_build_ipVolumes, & mesh_build_ipVolumes, &
mesh_build_ipCoordinates, & mesh_build_ipCoordinates, &
mesh_cellCenterCoordinates, & mesh_cellCenterCoordinates, &
mesh_init_postprocessing, &
mesh_get_Ncellnodes, & mesh_get_Ncellnodes, &
mesh_get_unitlength, & mesh_get_unitlength, &
mesh_get_nodeAtIP mesh_get_nodeAtIP
#ifdef Spectral #ifdef Spectral
public :: & public :: &
mesh_spectral_getGrid, & mesh_spectral_getGrid, &
mesh_spectral_getSize, & mesh_spectral_getSize
mesh_nodesAroundCentres, &
mesh_deformedCoordsFFT, &
mesh_volumeMismatch, &
mesh_shapeMismatch
#endif #endif
private :: & private :: &
@ -436,8 +427,7 @@ module mesh
mesh_spectral_build_nodes, & mesh_spectral_build_nodes, &
mesh_spectral_build_elements, & mesh_spectral_build_elements, &
mesh_spectral_build_ipNeighborhood, & mesh_spectral_build_ipNeighborhood, &
#endif #elif defined Marc4DAMASK
#ifdef Marc4DAMASK
mesh_marc_get_tableStyles, & mesh_marc_get_tableStyles, &
mesh_marc_count_nodesAndElements, & mesh_marc_count_nodesAndElements, &
mesh_marc_count_elementSets, & mesh_marc_count_elementSets, &
@ -448,8 +438,7 @@ module mesh
mesh_marc_build_nodes, & mesh_marc_build_nodes, &
mesh_marc_count_cpSizes, & mesh_marc_count_cpSizes, &
mesh_marc_build_elements, & mesh_marc_build_elements, &
#endif #elif defined Abaqus
#ifdef Abaqus
mesh_abaqus_count_nodesAndElements, & mesh_abaqus_count_nodesAndElements, &
mesh_abaqus_count_elementSets, & mesh_abaqus_count_elementSets, &
mesh_abaqus_count_materials, & mesh_abaqus_count_materials, &
@ -473,11 +462,7 @@ module mesh
mesh_tell_statistics, & mesh_tell_statistics, &
FE_mapElemtype, & FE_mapElemtype, &
mesh_faceMatch, & mesh_faceMatch, &
mesh_build_FEdata, & mesh_build_FEdata
mesh_write_cellGeom, &
mesh_write_elemGeom, &
mesh_write_meshfile, &
mesh_read_meshfile
contains contains
@ -487,9 +472,7 @@ contains
!! Order and routines strongly depend on type of solver !! Order and routines strongly depend on type of solver
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_init(ip,el) subroutine mesh_init(ip,el)
#ifdef Spectral
use, intrinsic :: iso_c_binding use, intrinsic :: iso_c_binding
#endif
use DAMASK_interface use DAMASK_interface
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)
use IO, only: & use IO, only: &
@ -531,11 +514,9 @@ subroutine mesh_init(ip,el)
integer(pInt) :: j integer(pInt) :: j
logical :: myDebug logical :: myDebug
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- mesh init -+>>>' write(6,'(/,a)') ' <<<+- mesh init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess
if (allocated(mesh_mapFEtoCPelem)) deallocate(mesh_mapFEtoCPelem) if (allocated(mesh_mapFEtoCPelem)) deallocate(mesh_mapFEtoCPelem)
if (allocated(mesh_mapFEtoCPnode)) deallocate(mesh_mapFEtoCPnode) if (allocated(mesh_mapFEtoCPnode)) deallocate(mesh_mapFEtoCPnode)
@ -562,25 +543,18 @@ subroutine mesh_init(ip,el)
myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt)
#ifdef Spectral #ifdef Spectral
#ifdef PETSc
call fftw_mpi_init() call fftw_mpi_init()
#endif
call IO_open_file(FILEUNIT,geometryFile) ! parse info from geometry file... call IO_open_file(FILEUNIT,geometryFile) ! parse info from geometry file...
if (myDebug) write(6,'(a)') ' Opened geometry file'; flush(6) if (myDebug) write(6,'(a)') ' Opened geometry file'; flush(6)
grid = mesh_spectral_getGrid(fileUnit) grid = mesh_spectral_getGrid(fileUnit)
geomSize = mesh_spectral_getSize(fileUnit) geomSize = mesh_spectral_getSize(fileUnit)
#ifdef PETSc
gridMPI = int(grid,C_INTPTR_T) gridMPI = int(grid,C_INTPTR_T)
alloc_local = fftw_mpi_local_size_3d(gridMPI(3), gridMPI(2), gridMPI(1)/2 +1, & alloc_local = fftw_mpi_local_size_3d(gridMPI(3), gridMPI(2), gridMPI(1)/2 +1, &
MPI_COMM_WORLD, local_K, local_K_offset) MPI_COMM_WORLD, local_K, local_K_offset)
grid3 = int(local_K,pInt) grid3 = int(local_K,pInt)
grid3Offset = int(local_K_offset,pInt) grid3Offset = int(local_K_offset,pInt)
size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal) size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal)
size3Offset = geomSize(3)*real(grid3Offset,pReal)/real(grid(3),pReal) size3Offset = geomSize(3)*real(grid3Offset,pReal)/real(grid(3),pReal)
#endif
if (myDebug) write(6,'(a)') ' Grid partitioned'; flush(6) if (myDebug) write(6,'(a)') ' Grid partitioned'; flush(6)
call mesh_spectral_count() call mesh_spectral_count()
if (myDebug) write(6,'(a)') ' Counted nodes/elements'; flush(6) if (myDebug) write(6,'(a)') ' Counted nodes/elements'; flush(6)
@ -592,8 +566,7 @@ subroutine mesh_init(ip,el)
if (myDebug) write(6,'(a)') ' Built nodes'; flush(6) if (myDebug) write(6,'(a)') ' Built nodes'; flush(6)
call mesh_spectral_build_elements(FILEUNIT) call mesh_spectral_build_elements(FILEUNIT)
if (myDebug) write(6,'(a)') ' Built elements'; flush(6) if (myDebug) write(6,'(a)') ' Built elements'; flush(6)
#endif #elif defined Marc4DAMASK
#ifdef Marc4DAMASK
call IO_open_inputFile(FILEUNIT,modelName) ! parse info from input file... call IO_open_inputFile(FILEUNIT,modelName) ! parse info from input file...
if (myDebug) write(6,'(a)') ' Opened input file'; flush(6) if (myDebug) write(6,'(a)') ' Opened input file'; flush(6)
call mesh_marc_get_tableStyles(FILEUNIT) call mesh_marc_get_tableStyles(FILEUNIT)
@ -616,8 +589,7 @@ subroutine mesh_init(ip,el)
if (myDebug) write(6,'(a)') ' Counted CP sizes'; flush(6) if (myDebug) write(6,'(a)') ' Counted CP sizes'; flush(6)
call mesh_marc_build_elements(FILEUNIT) call mesh_marc_build_elements(FILEUNIT)
if (myDebug) write(6,'(a)') ' Built elements'; flush(6) if (myDebug) write(6,'(a)') ' Built elements'; flush(6)
#endif #elif defined Abaqus
#ifdef Abaqus
call IO_open_inputFile(FILEUNIT,modelName) ! parse info from input file... call IO_open_inputFile(FILEUNIT,modelName) ! parse info from input file...
if (myDebug) write(6,'(a)') ' Opened input file'; flush(6) if (myDebug) write(6,'(a)') ' Opened input file'; flush(6)
noPart = IO_abaqus_hasNoPart(FILEUNIT) noPart = IO_abaqus_hasNoPart(FILEUNIT)
@ -666,15 +638,12 @@ subroutine mesh_init(ip,el)
if (myDebug) write(6,'(a)') ' Built shared elements'; flush(6) if (myDebug) write(6,'(a)') ' Built shared elements'; flush(6)
call mesh_build_ipNeighborhood call mesh_build_ipNeighborhood
#else #else
call mesh_spectral_build_ipNeighborhood(FILEUNIT) call mesh_spectral_build_ipNeighborhood
#endif #endif
if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6) if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6)
if (worldrank == 0_pInt) then if (worldrank == 0_pInt) then
call mesh_tell_statistics call mesh_tell_statistics
call mesh_write_meshfile
call mesh_write_cellGeom
call mesh_write_elemGeom
endif endif
if (usePingPong .and. (mesh_Nelems /= mesh_NcpElems)) & if (usePingPong .and. (mesh_Nelems /= mesh_NcpElems)) &
@ -963,7 +932,7 @@ subroutine mesh_build_ipCoordinates
do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell
myCoords = myCoords + mesh_cellnode(1:3,mesh_cell(n,i,e)) myCoords = myCoords + mesh_cellnode(1:3,mesh_cell(n,i,e))
enddo enddo
mesh_ipCoordinates(1:3,i,e) = myCoords / real(FE_NcellnodesPerCell(c),pReal) mesh_ipCoordinates(1:3,i,e) = myCoords / FE_NcellnodesPerCell(c)
enddo enddo
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -990,7 +959,7 @@ pure function mesh_cellCenterCoordinates(ip,el)
do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates + mesh_cellnode(1:3,mesh_cell(n,ip,el)) mesh_cellCenterCoordinates = mesh_cellCenterCoordinates + mesh_cellnode(1:3,mesh_cell(n,ip,el))
enddo enddo
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates / real(FE_NcellnodesPerCell(c),pReal) mesh_cellCenterCoordinates = mesh_cellCenterCoordinates / FE_NcellnodesPerCell(c)
end function mesh_cellCenterCoordinates end function mesh_cellCenterCoordinates
@ -1417,11 +1386,9 @@ end subroutine mesh_spectral_build_elements
!> @brief build neighborhood relations for spectral !> @brief build neighborhood relations for spectral
!> @details assign globals: mesh_ipNeighborhood !> @details assign globals: mesh_ipNeighborhood
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_build_ipNeighborhood(fileUnit) subroutine mesh_spectral_build_ipNeighborhood
implicit none implicit none
integer(pInt), intent(in) :: &
fileUnit
integer(pInt) :: & integer(pInt) :: &
x,y,z, & x,y,z, &
e e
@ -1558,332 +1525,8 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes)
nodes = nodes/8.0_pReal nodes = nodes/8.0_pReal
end function mesh_nodesAroundCentres end function mesh_nodesAroundCentres
!--------------------------------------------------------------------------------------------------
!> @brief calculate coordinates in current configuration for given defgrad
! using integration in Fourier space
!--------------------------------------------------------------------------------------------------
function mesh_deformedCoordsFFT(gDim,F,FavgIn,scalingIn) result(coords)
use IO, only: &
IO_error
use numerics, only: &
fftw_timelimit, &
fftw_planner_flag
use debug, only: &
debug_mesh, &
debug_level, &
debug_levelBasic
use math, only: &
PI, &
math_mul33x3
implicit none
real(pReal), intent(in), dimension(:,:,:,:,:) :: F
real(pReal), dimension(3,size(F,3),size(F,4),size(F,5)) :: coords
real(pReal), intent(in), dimension(3) :: gDim
real(pReal), intent(in), dimension(3,3), optional :: FavgIn
real(pReal), intent(in), dimension(3), optional :: scalingIn
! allocatable arrays for fftw c routines
type(C_PTR) :: planForth, planBack
type(C_PTR) :: coords_fftw, defgrad_fftw
real(pReal), dimension(:,:,:,:,:), pointer :: F_real
complex(pReal), dimension(:,:,:,:,:), pointer :: F_fourier
real(pReal), dimension(:,:,:,:), pointer :: coords_real
complex(pReal), dimension(:,:,:,:), pointer :: coords_fourier
! other variables
integer(pInt) :: i, j, k, m, res1Red
integer(pInt), dimension(3) :: k_s, iRes
real(pReal), dimension(3) :: scaling, step, offset_coords, integrator
real(pReal), dimension(3,3) :: Favg
integer(pInt), dimension(2:3,2) :: Nyquist ! highest frequencies to be removed (1 if even, 2 if odd)
if (present(scalingIn)) then
where (scalingIn < 0.0_pReal) ! invalid values. in case of f2py -1 if not present
scaling = [1.0_pReal,1.0_pReal,1.0_pReal]
elsewhere
scaling = scalingIn
end where
else
scaling = 1.0_pReal
endif
iRes = [size(F,3),size(F,4),size(F,5)]
integrator = gDim / 2.0_pReal / PI ! see notes where it is used
res1Red = iRes(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
step = gDim/real(iRes, pReal)
Nyquist(2,1:2) = [iRes(2)/2_pInt + 1_pInt, iRes(2)/2_pInt + 1_pInt + mod(iRes(2),2_pInt)]
Nyquist(3,1:2) = [iRes(3)/2_pInt + 1_pInt, iRes(3)/2_pInt + 1_pInt + mod(iRes(3),2_pInt)]
!--------------------------------------------------------------------------------------------------
! report
if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then
write(6,'(a)') ' Restore geometry using FFT-based integration'
write(6,'(a,3(i12 ))') ' grid a b c: ', iRes
write(6,'(a,3(es12.5))') ' size x y z: ', gDim
endif
!--------------------------------------------------------------------------------------------------
! sanity check
if (pReal /= C_DOUBLE .or. pInt /= C_INT) &
call IO_error(0_pInt,ext_msg='Fortran to C in mesh_deformedCoordsFFT')
!--------------------------------------------------------------------------------------------------
! allocation and FFTW initialization
defgrad_fftw = fftw_alloc_complex(int(res1Red *iRes(2)*iRes(3)*9_pInt,C_SIZE_T)) ! C_SIZE_T is of type integer(8)
coords_fftw = fftw_alloc_complex(int(res1Red *iRes(2)*iRes(3)*3_pInt,C_SIZE_T)) ! C_SIZE_T is of type integer(8)
call c_f_pointer(defgrad_fftw, F_real, &
[iRes(1)+2_pInt-mod(iRes(1),2_pInt),iRes(2),iRes(3),3_pInt,3_pInt])
call c_f_pointer(defgrad_fftw, F_fourier, &
[res1Red, iRes(2),iRes(3),3_pInt,3_pInt])
call c_f_pointer(coords_fftw, coords_real, &
[iRes(1)+2_pInt-mod(iRes(1),2_pInt),iRes(2),iRes(3),3_pInt])
call c_f_pointer(coords_fftw, coords_fourier, &
[res1Red, iRes(2),iRes(3),3_pInt])
call fftw_set_timelimit(fftw_timelimit)
planForth = fftw_plan_many_dft_r2c(3_pInt,[iRes(3),iRes(2) ,iRes(1)],9_pInt,& ! dimensions , length in each dimension in reversed order
F_real,[iRes(3),iRes(2) ,iRes(1)+2_pInt-mod(iRes(1),2_pInt)],& ! input data , physical length in each dimension in reversed order
1_pInt, iRes(3)*iRes(2)*(iRes(1)+2_pInt-mod(iRes(1),2_pInt)),& ! striding , product of physical lenght in the 3 dimensions
F_fourier,[iRes(3),iRes(2) ,res1Red],&
1_pInt, iRes(3)*iRes(2)* res1Red,fftw_planner_flag)
planBack = fftw_plan_many_dft_c2r(3_pInt,[iRes(3),iRes(2) ,iRes(1)],3_pInt,&
coords_fourier,[iRes(3),iRes(2) ,res1Red],&
1_pInt, iRes(3)*iRes(2)* res1Red,&
coords_real,[iRes(3),iRes(2) ,iRes(1)+2_pInt-mod(iRes(1),2_pInt)],&
1_pInt, iRes(3)*iRes(2)*(iRes(1)+2_pInt-mod(iRes(1),2_pInt)),&
fftw_planner_flag)
F_real(1:iRes(1),1:iRes(2),1:iRes(3),1:3,1:3) = &
reshape(F,[iRes(1),iRes(2),iRes(3),3,3], order = [4,5,1,2,3]) ! F_real is overwritten during plan creatio, is larger (padding) and has different order
!--------------------------------------------------------------------------------------------------
! FFT
call fftw_execute_dft_r2c(planForth, F_real, F_fourier)
!--------------------------------------------------------------------------------------------------
! if no average F is given, compute it in Fourier space
if (present(FavgIn)) then
if (all(FavgIn < 0.0_pReal)) then
Favg = real(F_fourier(1,1,1,1:3,1:3),pReal)/real(product(iRes),pReal) !the f2py way to tell it is not present
else
Favg = FavgIn
endif
else
Favg = real(F_fourier(1,1,1,1:3,1:3),pReal)/real(product(iRes),pReal)
endif
!--------------------------------------------------------------------------------------------------
! remove highest frequency in each direction, in third direction only if not 2D
if(iRes(1)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
F_fourier (res1Red, 1:iRes(2), 1:iRes(3), 1:3,1:3) &
= cmplx(0.0_pReal,0.0_pReal,pReal)
if(iRes(2)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
F_fourier (1:res1Red,Nyquist(2,1):Nyquist(2,2),1:iRes(3), 1:3,1:3) &
= cmplx(0.0_pReal,0.0_pReal,pReal)
if(iRes(3)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
F_fourier (1:res1Red,1:iRes(2), Nyquist(3,1):Nyquist(3,2),1:3,1:3) &
= cmplx(0.0_pReal,0.0_pReal,pReal)
!--------------------------------------------------------------------------------------------------
! integration in Fourier space
coords_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
do k = 1_pInt, iRes(3)
k_s(3) = k-1_pInt
if(k > iRes(3)/2_pInt+1_pInt) k_s(3) = k_s(3)-iRes(3)
do j = 1_pInt, iRes(2)
k_s(2) = j-1_pInt
if(j > iRes(2)/2_pInt+1_pInt) k_s(2) = k_s(2)-iRes(2)
do i = 1_pInt, res1Red
k_s(1) = i-1_pInt
do m = 1_pInt,3_pInt
coords_fourier(i,j,k,m) = sum(F_fourier(i,j,k,m,1:3)*&
cmplx(0.0_pReal,real(k_s,pReal)*integrator,pReal))
enddo
if (any(k_s /= 0_pInt)) coords_fourier(i,j,k,1:3) = &
coords_fourier(i,j,k,1:3) / cmplx(-sum(k_s*k_s),0.0_pReal,pReal)
enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! iFFT and freeing memory
call fftw_execute_dft_c2r(planBack,coords_fourier,coords_real)
coords = reshape(coords_real(1:iRes(1),1:iRes(2),1:iRes(3),1:3), [3,iRes(1),iRes(2),iRes(3)], &
order = [2,3,4,1])/real(product(iRes),pReal) ! weight and change order
call fftw_destroy_plan(planForth)
call fftw_destroy_plan(planBack)
call fftw_free(defgrad_fftw)
call fftw_free(coords_fftw)
!--------------------------------------------------------------------------------------------------
! add average to scaled fluctuation and put (0,0,0) on (0,0,0)
offset_coords = math_mul33x3(F(1:3,1:3,1,1,1),step/2.0_pReal) - scaling*coords(1:3,1,1,1)
forall(k = 1_pInt:iRes(3), j = 1_pInt:iRes(2), i = 1_pInt:iRes(1)) &
coords(1:3,i,j,k) = scaling(1:3)*coords(1:3,i,j,k) &
+ offset_coords &
+ math_mul33x3(Favg,step*real([i,j,k]-1_pInt,pReal))
end function mesh_deformedCoordsFFT
!--------------------------------------------------------------------------------------------------
!> @brief calculates the mismatch between volume of reconstructed (compatible) cube and
! determinant of defgrad at the FP
!--------------------------------------------------------------------------------------------------
function mesh_volumeMismatch(gDim,F,nodes) result(vMismatch)
use IO, only: &
IO_error
use debug, only: &
debug_mesh, &
debug_level, &
debug_levelBasic
use math, only: &
math_det33, &
math_volTetrahedron
implicit none
real(pReal), intent(in), dimension(:,:,:,:,:) :: &
F
real(pReal), dimension(size(F,3),size(F,4),size(F,5)) :: &
vMismatch
real(pReal), intent(in), dimension(:,:,:,:) :: &
nodes
real(pReal), dimension(3) :: &
gDim
integer(pInt), dimension(3) :: &
iRes
real(pReal), dimension(3,8) :: coords
integer(pInt) :: i,j,k
real(pReal) :: volInitial
iRes = [size(F,3),size(F,4),size(F,5)]
volInitial = product(gDim)/real(product(iRes), pReal)
!--------------------------------------------------------------------------------------------------
! report and check
if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then
write(6,'(a)') ' Calculating volume mismatch'
write(6,'(a,3(i12 ))') ' grid a b c: ', iRes
write(6,'(a,3(es12.5))') ' size x y z: ', gDim
endif
if (any([iRes/=size(nodes,2)-1_pInt,iRes/=size(nodes,3)-1_pInt,iRes/=size(nodes,4)-1_pInt]))&
call IO_error(0_pInt,ext_msg='Arrays F and nodes in mesh_volumeMismatch')
!--------------------------------------------------------------------------------------------------
! calculate actual volume and volume resulting from deformation gradient
do k = 1_pInt,iRes(3)
do j = 1_pInt,iRes(2)
do i = 1_pInt,iRes(1)
coords(1:3,1) = nodes(1:3,i, j, k )
coords(1:3,2) = nodes(1:3,i+1_pInt,j, k )
coords(1:3,3) = nodes(1:3,i+1_pInt,j+1_pInt,k )
coords(1:3,4) = nodes(1:3,i, j+1_pInt,k )
coords(1:3,5) = nodes(1:3,i, j, k+1_pInt)
coords(1:3,6) = nodes(1:3,i+1_pInt,j, k+1_pInt)
coords(1:3,7) = nodes(1:3,i+1_pInt,j+1_pInt,k+1_pInt)
coords(1:3,8) = nodes(1:3,i, j+1_pInt,k+1_pInt)
vMismatch(i,j,k) = &
abs(math_volTetrahedron(coords(1:3,7),coords(1:3,1),coords(1:3,8),coords(1:3,4))) &
+ abs(math_volTetrahedron(coords(1:3,7),coords(1:3,1),coords(1:3,8),coords(1:3,5))) &
+ abs(math_volTetrahedron(coords(1:3,7),coords(1:3,1),coords(1:3,3),coords(1:3,4))) &
+ abs(math_volTetrahedron(coords(1:3,7),coords(1:3,1),coords(1:3,3),coords(1:3,2))) &
+ abs(math_volTetrahedron(coords(1:3,7),coords(1:3,5),coords(1:3,2),coords(1:3,6))) &
+ abs(math_volTetrahedron(coords(1:3,7),coords(1:3,5),coords(1:3,2),coords(1:3,1)))
vMismatch(i,j,k) = vMismatch(i,j,k)/math_det33(F(1:3,1:3,i,j,k))
enddo; enddo; enddo
vMismatch = vMismatch/volInitial
end function mesh_volumeMismatch
!--------------------------------------------------------------------------------------------------
!> @brief Routine to calculate the mismatch between the vectors from the central point to
! the corners of reconstructed (combatible) volume element and the vectors calculated by deforming
! the initial volume element with the current deformation gradient
!--------------------------------------------------------------------------------------------------
function mesh_shapeMismatch(gDim,F,nodes,centres) result(sMismatch)
use IO, only: &
IO_error
use debug, only: &
debug_mesh, &
debug_level, &
debug_levelBasic
use math, only: &
math_mul33x3
implicit none
real(pReal), intent(in), dimension(:,:,:,:,:) :: &
F
real(pReal), dimension(size(F,3),size(F,4),size(F,5)) :: &
sMismatch
real(pReal), intent(in), dimension(:,:,:,:) :: &
nodes, &
centres
real(pReal), dimension(3) :: &
gDim, &
fRes
integer(pInt), dimension(3) :: &
iRes
real(pReal), dimension(3,8) :: coordsInitial
integer(pInt) i,j,k
iRes = [size(F,3),size(F,4),size(F,5)]
fRes = real(iRes,pReal)
!--------------------------------------------------------------------------------------------------
! report and check
if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then
write(6,'(a)') ' Calculating shape mismatch'
write(6,'(a,3(i12 ))') ' grid a b c: ', iRes
write(6,'(a,3(es12.5))') ' size x y z: ', gDim
endif
if(any([iRes/=size(nodes,2)-1_pInt,iRes/=size(nodes,3)-1_pInt,iRes/=size(nodes,4)-1_pInt]) .or.&
any([iRes/=size(centres,2), iRes/=size(centres,3), iRes/=size(centres,4)]))&
call IO_error(0_pInt,ext_msg='Arrays F and nodes/centres in mesh_shapeMismatch')
!--------------------------------------------------------------------------------------------------
! initial positions
coordsInitial(1:3,1) = [-gDim(1)/fRes(1),-gDim(2)/fRes(2),-gDim(3)/fRes(3)]
coordsInitial(1:3,2) = [+gDim(1)/fRes(1),-gDim(2)/fRes(2),-gDim(3)/fRes(3)]
coordsInitial(1:3,3) = [+gDim(1)/fRes(1),+gDim(2)/fRes(2),-gDim(3)/fRes(3)]
coordsInitial(1:3,4) = [-gDim(1)/fRes(1),+gDim(2)/fRes(2),-gDim(3)/fRes(3)]
coordsInitial(1:3,5) = [-gDim(1)/fRes(1),-gDim(2)/fRes(2),+gDim(3)/fRes(3)]
coordsInitial(1:3,6) = [+gDim(1)/fRes(1),-gDim(2)/fRes(2),+gDim(3)/fRes(3)]
coordsInitial(1:3,7) = [+gDim(1)/fRes(1),+gDim(2)/fRes(2),+gDim(3)/fRes(3)]
coordsInitial(1:3,8) = [-gDim(1)/fRes(1),+gDim(2)/fRes(2),+gDim(3)/fRes(3)]
coordsInitial = coordsInitial/2.0_pReal
!--------------------------------------------------------------------------------------------------
! compare deformed original and deformed positions to actual positions
do k = 1_pInt,iRes(3)
do j = 1_pInt,iRes(2)
do i = 1_pInt,iRes(1)
sMismatch(i,j,k) = &
sqrt(sum((nodes(1:3,i, j, k ) - centres(1:3,i,j,k)&
- math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,1)))**2.0_pReal))&
+ sqrt(sum((nodes(1:3,i+1_pInt,j, k ) - centres(1:3,i,j,k)&
- math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,2)))**2.0_pReal))&
+ sqrt(sum((nodes(1:3,i+1_pInt,j+1_pInt,k ) - centres(1:3,i,j,k)&
- math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,3)))**2.0_pReal))&
+ sqrt(sum((nodes(1:3,i, j+1_pInt,k ) - centres(1:3,i,j,k)&
- math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,4)))**2.0_pReal))&
+ sqrt(sum((nodes(1:3,i, j, k+1_pInt) - centres(1:3,i,j,k)&
- math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,5)))**2.0_pReal))&
+ sqrt(sum((nodes(1:3,i+1_pInt,j, k+1_pInt) - centres(1:3,i,j,k)&
- math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,6)))**2.0_pReal))&
+ sqrt(sum((nodes(1:3,i+1_pInt,j+1_pInt,k+1_pInt) - centres(1:3,i,j,k)&
- math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,7)))**2.0_pReal))&
+ sqrt(sum((nodes(1:3,i, j+1_pInt,k+1_pInt) - centres(1:3,i,j,k)&
- math_mul33x3(F(1:3,1:3,i,j,k), coordsInitial(1:3,8)))**2.0_pReal))
enddo; enddo; enddo
end function mesh_shapeMismatch
#endif #endif
#ifdef Marc4DAMASK #ifdef Marc4DAMASK
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Figures out table styles (Marc only) and stores to 'initialcondTableStyle' and !> @brief Figures out table styles (Marc only) and stores to 'initialcondTableStyle' and
@ -3070,6 +2713,7 @@ use IO, only: &
implicit none implicit none
integer(pInt), intent(in) :: fileUnit integer(pInt), intent(in) :: fileUnit
#ifndef Spectral #ifndef Spectral
integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) chunk, Nchunks integer(pInt) chunk, Nchunks
@ -3081,9 +2725,10 @@ use IO, only: &
mesh_periodicSurface = .true. mesh_periodicSurface = .true.
#else #else
mesh_periodicSurface = .false. mesh_periodicSurface = .false.
#if defined(Marc4DAMASK) #ifdef Marc4DAMASK
keyword = '$damask' keyword = '$damask'
#elif defined(Abaqus) #endif
#ifdef Abaqus
keyword = '**damask' keyword = '**damask'
#endif #endif
@ -3691,7 +3336,6 @@ integer(pInt) function FE_mapElemtype(what)
'c3d20t') 'c3d20t')
FE_mapElemtype = 13_pInt ! Three-dimensional Arbitrarily Distorted quadratic hexahedral FE_mapElemtype = 13_pInt ! Three-dimensional Arbitrarily Distorted quadratic hexahedral
case default case default
FE_mapElemtype = -1_pInt ! error return
call IO_error(error_ID=190_pInt,ext_msg=IO_lc(what)) call IO_error(error_ID=190_pInt,ext_msg=IO_lc(what))
end select end select
@ -3700,7 +3344,6 @@ end function FE_mapElemtype
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief find face-matching element of same type !> @brief find face-matching element of same type
!> @details currently not used, check if needed for HDF5 output, otherwise delete
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_faceMatch(elem, face ,matchingElem, matchingFace) subroutine mesh_faceMatch(elem, face ,matchingElem, matchingFace)
@ -4511,212 +4154,6 @@ subroutine mesh_build_FEdata
end subroutine mesh_build_FEdata end subroutine mesh_build_FEdata
!--------------------------------------------------------------------------------------------------
!> @brief writes out initial cell geometry
!--------------------------------------------------------------------------------------------------
subroutine mesh_write_cellGeom
use DAMASK_interface, only: &
getSolverJobName, &
getSolverWorkingDirectoryName
use IR_Precision, only: &
I4P
use Lib_VTK_IO, only: &
VTK_ini, &
VTK_geo, &
VTK_con, &
VTK_end
implicit none
integer(I4P), dimension(1:mesh_Ncells) :: celltype
integer(I4P), dimension(mesh_Ncells*(1_pInt+FE_maxNcellnodesPerCell)) :: cellconnection
integer(I4P):: error
integer(I4P):: g, c, e, CellID, i, j
cellID = 0_pInt
j = 0_pInt
do e = 1_pInt, mesh_NcpElems ! loop over cpElems
g = FE_geomtype(mesh_element(2_pInt,e)) ! get geometry type
c = FE_celltype(g) ! get cell type
do i = 1_pInt,FE_Nips(g) ! loop over ips=cells in this element
cellID = cellID + 1_pInt
celltype(cellID) = MESH_VTKCELLTYPE(c)
cellconnection(j+1_pInt:j+FE_NcellnodesPerCell(c)+1_pInt) &
= [FE_NcellnodesPerCell(c),mesh_cell(1:FE_NcellnodesPerCell(c),i,e)-1_pInt] ! number of cellnodes per cell & list of global cellnode IDs belnging to this cell (cellnode counting starts at 0)
j = j + FE_NcellnodesPerCell(c) + 1_pInt
enddo
enddo
error=VTK_ini(output_format = 'ASCII', &
title=trim(getSolverJobName())//' cell mesh', &
filename = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'_ipbased.vtk', &
mesh_topology = 'UNSTRUCTURED_GRID')
!ToDo: check error here
error=VTK_geo(NN = int(mesh_Ncellnodes,I4P), &
X = mesh_cellnode(1,1:mesh_Ncellnodes), &
Y = mesh_cellnode(2,1:mesh_Ncellnodes), &
Z = mesh_cellnode(3,1:mesh_Ncellnodes))
!ToDo: check error here
error=VTK_con(NC = int(mesh_Ncells,I4P), &
connect = cellconnection(1:j), &
!ToDo: check error here
cell_type = celltype)
error=VTK_end()
!ToDo: check error here
end subroutine mesh_write_cellGeom
!--------------------------------------------------------------------------------------------------
!> @brief writes out initial element geometry
!--------------------------------------------------------------------------------------------------
subroutine mesh_write_elemGeom
use DAMASK_interface, only: &
getSolverJobName, &
getSolverWorkingDirectoryName
use IR_Precision, only: &
I4P
use Lib_VTK_IO, only: &
VTK_ini, &
VTK_geo, &
VTK_con, &
VTK_end
implicit none
integer(I4P), dimension(1:mesh_NcpElems) :: elemtype
integer(I4P), dimension(mesh_NcpElems*(1_pInt+FE_maxNnodes)) :: elementconnection
integer(I4P):: error
integer(pInt):: e, t, n, i
i = 0_pInt
do e = 1_pInt, mesh_NcpElems ! loop over cpElems
t = mesh_element(2,e) ! get element type
elemtype(e) = MESH_VTKELEMTYPE(t)
elementconnection(i+1_pInt) = FE_Nnodes(t) ! number of nodes per element
do n = 1_pInt,FE_Nnodes(t)
elementconnection(i+1_pInt+n) = mesh_element(4_pInt+n,e) - 1_pInt ! global node ID of node that belongs to this element (node counting starts at 0)
enddo
i = i + 1_pInt + FE_Nnodes(t)
enddo
error=VTK_ini(output_format = 'ASCII', &
title=trim(getSolverJobName())//' element mesh', &
filename = trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'_nodebased.vtk', &
mesh_topology = 'UNSTRUCTURED_GRID')
!ToDo: check error here
error=VTK_geo(NN = int(mesh_Nnodes,I4P), &
X = mesh_node0(1,1:mesh_Nnodes), &
Y = mesh_node0(2,1:mesh_Nnodes), &
Z = mesh_node0(3,1:mesh_Nnodes))
!ToDo: check error here
error=VTK_con(NC = int(mesh_Nelems,I4P), &
connect = elementconnection(1:i), &
cell_type = elemtype)
!ToDo: check error here
error =VTK_end()
!ToDo: check error here
end subroutine mesh_write_elemGeom
!--------------------------------------------------------------------------------------------------
!> @brief writes description file for mesh
!--------------------------------------------------------------------------------------------------
subroutine mesh_write_meshfile
use IO, only: &
IO_write_jobFile
implicit none
integer(pInt), parameter :: fileUnit = 223_pInt
integer(pInt) :: e,i,t,g,c,n
call IO_write_jobFile(fileUnit,'mesh')
write(fileUnit,'(A16,E10.3)') 'unitlength', mesh_unitlength
write(fileUnit,'(A16,I10)') 'maxNcellnodes', mesh_maxNcellnodes
write(fileUnit,'(A16,I10)') 'maxNips', mesh_maxNips
write(fileUnit,'(A16,I10)') 'maxNnodes', mesh_maxNnodes
write(fileUnit,'(A16,I10)') 'Nnodes', mesh_Nnodes
write(fileUnit,'(A16,I10)') 'NcpElems', mesh_NcpElems
do e = 1_pInt,mesh_NcpElems
t = mesh_element(2,e)
write(fileUnit,'(20(I10))') mesh_element(1_pInt:4_pInt+FE_Nnodes(t),e)
enddo
write(fileUnit,'(A16,I10)') 'Ncellnodes', mesh_Ncellnodes
do n = 1_pInt,mesh_Ncellnodes
write(fileUnit,'(2(I10))') mesh_cellnodeParent(1:2,n)
enddo
write(fileUnit,'(A16,I10)') 'Ncells', mesh_Ncells
do e = 1_pInt,mesh_NcpElems
t = mesh_element(2,e)
g = FE_geomtype(t)
c = FE_celltype(g)
do i = 1_pInt,FE_Nips(g)
write(fileUnit,'(8(I10))') &
mesh_cell(1_pInt:FE_NcellnodesPerCell(c),i,e)
enddo
enddo
close(fileUnit)
end subroutine mesh_write_meshfile
!--------------------------------------------------------------------------------------------------
!> @brief reads mesh description file
!--------------------------------------------------------------------------------------------------
integer function mesh_read_meshfile(filepath)
implicit none
character(len=*), intent(in) :: filepath
integer(pInt), parameter :: fileUnit = 223_pInt
integer(pInt) :: e,i,t,g,n
open(fileUnit,status='old',err=100,iostat=mesh_read_meshfile,action='read',file=filepath)
read(fileUnit,'(TR16,E10.3)',err=100,iostat=mesh_read_meshfile) mesh_unitlength
read(fileUnit,'(TR16,I10)',err=100,iostat=mesh_read_meshfile) mesh_maxNcellnodes
read(fileUnit,'(TR16,I10)',err=100,iostat=mesh_read_meshfile) mesh_maxNips
read(fileUnit,'(TR16,I10)',err=100,iostat=mesh_read_meshfile) mesh_maxNnodes
read(fileUnit,'(TR16,I10)',err=100,iostat=mesh_read_meshfile) mesh_Nnodes
read(fileUnit,'(TR16,I10)',err=100,iostat=mesh_read_meshfile) mesh_NcpElems
if (.not. allocated(mesh_element)) allocate(mesh_element(4_pInt+mesh_maxNnodes,mesh_NcpElems))
mesh_element = 0_pInt
do e = 1_pInt,mesh_NcpElems
read(fileUnit,'(20(I10))',err=100,iostat=mesh_read_meshfile) &
mesh_element(:,e)
enddo
read(fileUnit,'(TR16,I10)',err=100,iostat=mesh_read_meshfile) mesh_Ncellnodes
if (.not. allocated(mesh_cellnodeParent)) allocate(mesh_cellnodeParent(2_pInt,mesh_Ncellnodes))
do n = 1_pInt,mesh_Ncellnodes
read(fileUnit,'(2(I10))',err=100,iostat=mesh_read_meshfile) mesh_cellnodeParent(1:2,n)
enddo
read(fileUnit,'(TR16,I10)',err=100,iostat=mesh_read_meshfile) mesh_Ncells
if (.not. allocated(mesh_cell)) allocate(mesh_cell(FE_maxNcellnodesPerCell,mesh_maxNips,mesh_NcpElems))
do e = 1_pInt,mesh_NcpElems
t = mesh_element(2,e)
g = FE_geomtype(t)
do i = 1_pInt,FE_Nips(g)
read(fileUnit,'(8(I10))',err=100,iostat=mesh_read_meshfile) mesh_cell(:,i,e)
enddo
enddo
close(fileUnit)
mesh_read_meshfile = 0 ! successfully read data
100 continue
end function mesh_read_meshfile
!--------------------------------------------------------------------------------------------------
!> @brief initializes mesh data for use in post processing
!--------------------------------------------------------------------------------------------------
integer function mesh_init_postprocessing(filepath)
implicit none
character(len=*), intent(in) :: filepath
call mesh_build_FEdata
mesh_init_postprocessing = mesh_read_meshfile(filepath)
end function mesh_init_postprocessing
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns global variable mesh_Ncellnodes !> @brief returns global variable mesh_Ncellnodes
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

28
code/spectral_interface.f90
View File

@ -13,9 +13,7 @@ module DAMASK_interface
pInt pInt
implicit none implicit none
private private
#ifdef PETSc
#include <petsc/finclude/petscsys.h> #include <petsc/finclude/petscsys.h>
#endif
logical, public, protected :: appendToOutFile = .false. !< Append to existing spectralOut file (in case of restart, not in case of regridding) logical, public, protected :: appendToOutFile = .false. !< Append to existing spectralOut file (in case of restart, not in case of regridding)
integer(pInt), public, protected :: spectralRestartInc = 1_pInt !< Increment at which calculation starts integer(pInt), public, protected :: spectralRestartInc = 1_pInt !< Increment at which calculation starts
character(len=1024), public, protected :: & character(len=1024), public, protected :: &
@ -44,13 +42,10 @@ contains
!> @brief initializes the solver by interpreting the command line arguments. Also writes !> @brief initializes the solver by interpreting the command line arguments. Also writes
!! information on computation to screen !! information on computation to screen
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn) subroutine DAMASK_interface_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)
implicit none implicit none
character(len=1024), optional, intent(in) :: &
loadCaseParameterIn, & !< if using the f2py variant, the -l argument of DAMASK_spectral.exe
geometryParameterIn !< if using the f2py variant, the -g argument of DAMASK_spectral.exe
character(len=1024) :: & character(len=1024) :: &
commandLine, & !< command line call as string commandLine, & !< command line call as string
loadCaseArg ='', & !< -l argument given to DAMASK_spectral.exe loadCaseArg ='', & !< -l argument given to DAMASK_spectral.exe
@ -67,9 +62,7 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn)
chunkPos chunkPos
integer, dimension(8) :: & integer, dimension(8) :: &
dateAndTime ! type default integer dateAndTime ! type default integer
#ifdef PETSc
PetscErrorCode :: ierr PetscErrorCode :: ierr
#endif
external :: & external :: &
quit,& quit,&
MPI_Comm_rank,& MPI_Comm_rank,&
@ -81,7 +74,6 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc Init ! PETSc Init
#ifdef PETSc
#ifdef _OPENMP #ifdef _OPENMP
call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,ierr);CHKERRQ(ierr) ! in case of OpenMP, don't rely on PETScInitialize doing MPI init call MPI_Init_Thread(MPI_THREAD_FUNNELED,threadLevel,ierr);CHKERRQ(ierr) ! in case of OpenMP, don't rely on PETScInitialize doing MPI init
if (threadLevel<MPI_THREAD_FUNNELED) then if (threadLevel<MPI_THREAD_FUNNELED) then
@ -92,7 +84,6 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn)
call PetscInitialize(PETSC_NULL_CHARACTER,ierr) ! according to PETSc manual, that should be the first line in the code call PetscInitialize(PETSC_NULL_CHARACTER,ierr) ! according to PETSc manual, that should be the first line in the code
CHKERRQ(ierr) ! this is a macro definition, it is case sensitive CHKERRQ(ierr) ! this is a macro definition, it is case sensitive
call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr) call MPI_Comm_rank(PETSC_COMM_WORLD,worldrank,ierr);CHKERRQ(ierr)
#endif
mainProcess: if (worldrank == 0) then mainProcess: if (worldrank == 0) then
if (output_unit /= 6) then if (output_unit /= 6) then
@ -116,11 +107,6 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn)
write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>'
#include "compilation_info.f90" #include "compilation_info.f90"
if ( present(loadcaseParameterIn) .and. present(geometryParameterIn)) then ! both mandatory parameters given in function call
geometryArg = geometryParameterIn
loadcaseArg = loadcaseParameterIn
commandLine = 'n/a'
else if ( .not.( present(loadcaseParameterIn) .and. present(geometryParameterIn))) then ! none parameters given in function call, trying to get them from command line
call get_command(commandLine) call get_command(commandLine)
chunkPos = IIO_stringPos(commandLine) chunkPos = IIO_stringPos(commandLine)
do i = 1, chunkPos(1) do i = 1, chunkPos(1)
@ -138,7 +124,6 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn)
write(6,'(a)') ' --load (-l, --loadcase)' write(6,'(a)') ' --load (-l, --loadcase)'
write(6,'(a)') ' --workingdir (-w, --wd, --workingdirectory, -d, --directory)' write(6,'(a)') ' --workingdir (-w, --wd, --workingdirectory, -d, --directory)'
write(6,'(a)') ' --restart (-r, --rs)' write(6,'(a)') ' --restart (-r, --rs)'
write(6,'(a)') ' --regrid (--rg)'
write(6,'(a)') ' --help (-h)' write(6,'(a)') ' --help (-h)'
write(6,'(/,a)')' -----------------------------------------------------------------------' write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(a)') ' Mandatory arguments:' write(6,'(a)') ' Mandatory arguments:'
@ -170,13 +155,6 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn)
write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".' write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".'
write(6,'(a)') ' Works only if the restart information for total increment' write(6,'(a)') ' Works only if the restart information for total increment'
write(6,'(a)') ' No. XX-1 is available in the working directory.' write(6,'(a)') ' No. XX-1 is available in the working directory.'
write(6,'(/,a)')' --regrid XX'
write(6,'(a)') ' Reads in total increment No. XX-1 and continues to'
write(6,'(a)') ' calculate total increment No. XX.'
write(6,'(a)') ' Attention: Overwrites existing results file '
write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".'
write(6,'(a)') ' Works only if the restart information for total increment'
write(6,'(a)') ' No. XX-1 is available in the working directory.'
write(6,'(/,a)')' -----------------------------------------------------------------------' write(6,'(/,a)')' -----------------------------------------------------------------------'
write(6,'(a)') ' Help:' write(6,'(a)') ' Help:'
write(6,'(/,a)')' --help' write(6,'(/,a)')' --help'
@ -192,12 +170,8 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn)
case ('-r', '--rs', '--restart') case ('-r', '--rs', '--restart')
spectralRestartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt) spectralRestartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt)
appendToOutFile = .true. appendToOutFile = .true.
case ('--rg', '--regrid')
spectralRestartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt)
appendToOutFile = .false.
end select end select
enddo enddo
endif
if (len(trim(loadcaseArg)) == 0 .or. len(trim(geometryArg)) == 0) then if (len(trim(loadcaseArg)) == 0 .or. len(trim(geometryArg)) == 0) then
write(6,'(a)') ' Please specify geometry AND load case (-h for help)' write(6,'(a)') ' Please specify geometry AND load case (-h for help)'

159
configure vendored
View File

@ -22,19 +22,6 @@ class extendableOption(Option):
Option.take_action(self, action, dest, opt, value, values, parser) Option.take_action(self, action, dest, opt, value, values, parser)
# -----------------------------
def filePresent(paths,files,warning=False):
for path in paths:
for file in files:
if os.path.isfile(os.path.join(path,file)): return True
if warning: print "Warning: %s not found in %s"%(','.join(files),','.join(paths))
return False
######################################################## ########################################################
# MAIN # MAIN
######################################################## ########################################################
@ -42,35 +29,15 @@ def filePresent(paths,files,warning=False):
parser = OptionParser(option_class=extendableOption, usage='%prog options', description = """ parser = OptionParser(option_class=extendableOption, usage='%prog options', description = """
Configures the compilation and installation of DAMASK Configures the compilation and installation of DAMASK
""" + string.replace('$Id$','\n','\\n') """)
)
#--- determine default compiler ----------------------------------------------------------------------
compiler = os.getenv('F90')
if compiler == None:
compiler = 'ifort' if subprocess.call(['which', 'ifort'], stdout=subprocess.PIPE, stderr=subprocess.PIPE) == 0 \
else 'gfortran'
#--- default option values --------------------------------------------------------------------------
BLAS_order = ['IMKL','ACML','LAPACK','OPENBLAS']
defaults={'DAMASK_BIN':'depending on access rights', defaults={'DAMASK_BIN':'depending on access rights',
'F90':compiler,
'FFTW_ROOT':'/usr', 'FFTW_ROOT':'/usr',
'MSC_ROOT' :'/msc', 'MSC_ROOT' :'/msc',
'DAMASK_NUM_THREADS':4, 'DAMASK_NUM_THREADS':4,
'MARC_VERSION':'2015', 'MARC_VERSION':'2015',
'blasType':'LAPACK',
'blasRoot':{'LAPACK' :'/usr',
'ACML' :'/opt/acml6.1.0',
'IMKL' : os.getenv('MKLROOT') if os.getenv('MKLROOT') else '/opt/intel/composerxe/mkl',
'OPENBLAS' :'/usr',
},
'spectralOptions':{}, 'spectralOptions':{},
} }
#--- if local config file exists, read, otherwise assume global config file ------------------------ #--- if local config file exists, read, otherwise assume global config file ------------------------
configFile = os.path.join(os.getenv('HOME'),'.damask/damask.conf') \ configFile = os.path.join(os.getenv('HOME'),'.damask/damask.conf') \
if os.path.isfile(os.path.join(os.getenv('HOME'),'.damask/damask.conf')) \ if os.path.isfile(os.path.join(os.getenv('HOME'),'.damask/damask.conf')) \
@ -91,129 +58,25 @@ try:
defaults['DAMASK_NUM_THREADS'] = int(value) defaults['DAMASK_NUM_THREADS'] = int(value)
if key == 'DAMASK_BIN': if key == 'DAMASK_BIN':
defaults['DAMASK_BIN'] = value defaults['DAMASK_BIN'] = value
if key in ['F90','FFTW_ROOT','MSC_ROOT','spectralOptions','MARC_VERSION']:
defaults[key] = value
for theKey in reversed(BLAS_order):
if key == theKey+'_ROOT' and value != None and value != '':
defaults['blasType'] = theKey
defaults['blasRoot'][theKey] = value
except IOError: except IOError:
pass pass
parser.add_option('--prefix', dest='prefix', metavar='string', parser.add_option('--prefix', dest='prefix', metavar='string',
help='location of (links to) DAMASK executables [%default]') help='location of (links to) DAMASK executables [%default]')
parser.add_option('--with-FC','--with-fc',
dest='compiler', metavar='string',
help='F90 compiler [%default]')
parser.add_option('--with-FFTW-dir','--with-fftw-dir',
dest='fftwRoot', metavar='string',
help='root directory of FFTW [%default]')
parser.add_option('--with-MSC-dir','--with-msc-dir',
dest='mscRoot', metavar='string',
help='root directory of MSC.Marc/Mentat [%default]')
parser.add_option('--with-MARC-version','--with-marc-version',
dest='marcVersion', metavar='string',
help='version of MSC.Marc/Mentat [%default]')
parser.add_option('--with-OMP-threads','--with-omp-threads', parser.add_option('--with-OMP-threads','--with-omp-threads',
dest='threads', type='int', metavar='int', dest='threads', type='int', metavar='int',
help='number of openMP threads [%default]') help='number of openMP threads [%default]')
parser.add_option('--with-BLAS-type','--with-blas-type',
dest='blasType', metavar='string',
help='type of BLAS/LAPACK library [%default] {{{}}}'.format(','.join(BLAS_order)))
parser.add_option('--with-BLAS-dir','--with-blas-dir',
dest='blasRoot',metavar='string',
help='root directory of BLAS/LAPACK library [%default]')
parser.add_option('--with-spectral-options', dest='spectraloptions', action='extend', metavar='<string LIST>', parser.add_option('--with-spectral-options', dest='spectraloptions', action='extend', metavar='<string LIST>',
help='options for compilation of spectral solver') help='options for compilation of spectral solver')
parser.set_defaults(prefix = defaults['DAMASK_BIN']) parser.set_defaults(prefix = defaults['DAMASK_BIN'])
parser.set_defaults(compiler = defaults['F90'])
parser.set_defaults(fftwRoot = defaults['FFTW_ROOT'])
parser.set_defaults(mscRoot = defaults['MSC_ROOT']) parser.set_defaults(mscRoot = defaults['MSC_ROOT'])
parser.set_defaults(marcVersion = defaults['MARC_VERSION']) parser.set_defaults(marcVersion = defaults['MARC_VERSION'])
parser.set_defaults(threads = defaults['DAMASK_NUM_THREADS']) parser.set_defaults(threads = defaults['DAMASK_NUM_THREADS'])
parser.set_defaults(blasType = defaults['blasType'])
#--- set default for blasRoot depending on current option (or default) for blasType --------------------
blasType = defaults['blasType'].upper()
for i, arg in enumerate(sys.argv):
if arg.lower().startswith('--with-blas-type'):
if arg.lower().endswith('--with-blas-type'):
blasType = sys.argv[i+1].upper()
else:
blasType = sys.argv[i][17:].upper()
if blasType not in BLAS_order:
blasType = defaults['blasType'].upper()
parser.set_defaults(blasRoot = defaults['blasRoot'][blasType])
parser.set_defaults(spectraloptions = []) parser.set_defaults(spectraloptions = [])
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()
#--- consistency checks --------------------------------------------------------------------------------
options.compiler = options.compiler.lower()
options.blasType = options.blasType.upper()
options.fftwRoot = os.path.normpath(options.fftwRoot)
options.mscRoot = os.path.normpath(options.mscRoot)
options.blasRoot = os.path.normpath(options.blasRoot)
locations = {
'FFTW' : [os.path.join(options.fftwRoot,'lib64'),os.path.join(options.fftwRoot,'lib')],
'LAPACK' : [os.path.join(options.blasRoot,'lib64'),os.path.join(options.blasRoot,'lib')],
'OPENBLAS': [os.path.join(options.blasRoot,'lib64'),os.path.join(options.blasRoot,'lib')],
'ACML' : [os.path.join(options.blasRoot,'%s64/lib'%options.compiler)],
'ACML_mp' : [os.path.join(options.blasRoot,'%s64_mp/lib'%options.compiler)],
'IMKL' : [os.path.join(options.blasRoot,'lib/intel64')],
}
libraries = {
'FFTW' : [
'libfftw3.so','libfftw3.a',
'libfftw3_threads.so','libfftw3_threads.a',
],
'LAPACK' : [
'liblapack.so','liblapack.a','liblapack.dylib',
],
'OPENBLAS' : [
'libopenblas.so','libopenblas.a','libopenblas.dylib',
],
'ACML' : [
'libacml.so','libacml.a',
],
'ACML_mp' : [
'libacml_mp.so','libacml_mp.a',
],
'IMKL' : [
'libmkl_core.so','libmkl_core.a',
'libmkl_sequential.so','libmkl_sequential.a',
'libmkl_intel_thread.so','libmkl_intel_thread.a',
'libmkl_intel_lp64.so','libmkl_intel_lp64.a',
'libmkl_gnu_thread.so','libmkl_gnu_thread.a',
'libmkl_gf_lp64.so','libmkl_gf_lp64.a',
],
}
if options.compiler not in ['ifort','gfortran']:
print('Error: Unknown compiler option: %s'%options.compiler)
sys.exit(1)
if not subprocess.call(['which', options.compiler], stdout=subprocess.PIPE, stderr=subprocess.PIPE) == 0:
print('Compiler Warning: executable %s not found!'%options.compiler)
if not os.path.isdir(options.mscRoot):
print('Warning: MSC root directory %s not found!'%options.mscRoot)
filePresent(locations['FFTW'],libraries['FFTW'],warning=True)
if options.blasType in ['LAPACK','OPENBLAS','IMKL']:
filePresent(locations[options.blasType],libraries[options.blasType],warning=True)
elif options.blasType == 'ACML':
filePresent(locations[options.blasType],libraries[options.blasType],warning=True)
filePresent(locations[options.blasType+'_mp'],libraries[options.blasType+'_mp'],warning=True)
else:
print('Error: Unknown BLAS/LAPACK library: %s'%options.blasType)
sys.exit(1)
#--- read config file if present to keep comments and order --------------------------------------- #--- read config file if present to keep comments and order ---------------------------------------
output = [] output = []
@ -228,12 +91,6 @@ try:
if items[0] == 'DAMASK_BIN': if items[0] == 'DAMASK_BIN':
line = '%s=%s'%(items[0],options.prefix) line = '%s=%s'%(items[0],options.prefix)
options.prefix ='depending on access rights' options.prefix ='depending on access rights'
if items[0] == 'F90':
line = '%s=%s'%(items[0],options.compiler)
options.compiler =''
if items[0] == 'FFTW_ROOT':
line = '%s=%s'%(items[0],options.fftwRoot)
options.fftwRoot =''
if items[0] == 'MSC_ROOT': if items[0] == 'MSC_ROOT':
line = '%s=%s'%(items[0],options.mscRoot) line = '%s=%s'%(items[0],options.mscRoot)
options.mscRoot ='' options.mscRoot =''
@ -243,14 +100,6 @@ try:
if items[0] == 'DAMASK_NUM_THREADS': if items[0] == 'DAMASK_NUM_THREADS':
line = '%s=%s'%(items[0],options.threads) line = '%s=%s'%(items[0],options.threads)
options.threads ='' options.threads =''
for blasType in defaults['blasRoot'].keys():
if items[0] == '%s_ROOT'%blasType and items[0] == '%s_ROOT'%options.blasType:
line = '%s=%s'%(items[0],options.blasRoot)
options.blasType=''
elif items[0] == '#%s_ROOT'%blasType and items[0] == '#%s_ROOT'%options.blasType:
line = '%s=%s'%(items[0][1:],options.blasRoot)
options.blasType=''
elif items[0] == '%s_ROOT'%blasType: line = '#'+line
for spectralOption in options.spectraloptions: for spectralOption in options.spectraloptions:
[key,value] = re.split('[= ]',spectralOption)[0:2] [key,value] = re.split('[= ]',spectralOption)[0:2]
if key == items[0]: if key == items[0]:
@ -264,18 +113,12 @@ except IOError:
for opt, value in options.__dict__.items(): for opt, value in options.__dict__.items():
if opt == 'prefix' and value != 'depending on access rights': if opt == 'prefix' and value != 'depending on access rights':
output.append('DAMASK_BIN=%s'%value) output.append('DAMASK_BIN=%s'%value)
if opt == 'compiler' and value != '':
output.append('F90=%s'%value)
if opt == 'fftwRoot' and value != '':
output.append('FFTW_ROOT=%s'%value)
if opt == 'mscRoot' and value != '': if opt == 'mscRoot' and value != '':
output.append('MSC_ROOT=%s'%value) output.append('MSC_ROOT=%s'%value)
if opt == 'marcVersion' and value != '': if opt == 'marcVersion' and value != '':
output.append('MARC_VERSION=%s'%value) output.append('MARC_VERSION=%s'%value)
if opt == 'threads' and value != '': if opt == 'threads' and value != '':
output.append('DAMASK_NUM_THREADS=%s'%value) output.append('DAMASK_NUM_THREADS=%s'%value)
if opt == 'blasType' and value != '':
output.append('%s_ROOT=%s'%(options.blasType,options.blasRoot))
for spectralOption in options.spectraloptions: for spectralOption in options.spectraloptions:
output.append(spectralOption) output.append(spectralOption)

141
installation/compile_CoreModule.py
View File

@ -1,141 +0,0 @@
#!/usr/bin/env python2
# -*- coding: UTF-8 no BOM -*-
import os,sys,glob,subprocess,shlex
from damask import Environment
from damask import version as DAMASKVERSION
# compiles fortran code for Python
scriptID = '$Id$'
damaskEnv = Environment()
baseDir = damaskEnv.relPath('installation/')
codeDir = damaskEnv.relPath('code/')
keywords=['IMKL_ROOT','ACML_ROOT','LAPACK_ROOT','FFTW_ROOT','F90']
options={}
#--- getting options from damask.conf or, if not present, from envinronment -----------------------
for option in keywords:
try:
value = damaskEnv.options[option]
except:
value = os.getenv(option)
if value is None: value = '' # env not set
options[option]=value
#--- overwrite default options with keyword=value pair from argument list to mimic make behavior --
for i, arg in enumerate(sys.argv):
for option in keywords:
if arg.startswith(option):
options[option] = sys.argv[i][len(option)+1:]
#--- check for valid compiler and set options -----------------------------------------------------
compilers = ['ifort','gfortran']
if options['F90'] not in compilers:
sys.exit('compiler "F90" (in installation/options or as Shell variable) has to be one out of: %s'%(', '.join(compilers)))
compiler = {
'gfortran': '--fcompiler=gnu95 --f90flags="-fPIC -fno-range-check -xf95-cpp-input -std=f2008 -fall-intrinsics'+\
' -fdefault-real-8 -fdefault-double-8"',
'ifort': '--fcompiler=intelem --f90flags="-fPIC -fpp -stand f08 -diag-disable 5268 -assume byterecl'+\
' -real-size 64 -integer-size 32 -shared-intel"',
}[options['F90']]
#--- option not depending on compiler -------------------------------------------------------------
compileOptions = ' -DSpectral -DFLOAT=8 -DINT=4 -I%s/lib -DDAMASKVERSION=\\\\\"\\\"%s\\\\\"\\\"'%(damaskEnv.rootDir(),DAMASKVERSION)
#--- this saves the path of libraries to core.so, hence it is known during runtime ----------------
if options['F90'] == 'gfortran':
# solved error: Undefined symbols for architecture x86_64: "_PyArg_ParseTupleAndKeywords"
# as found on https://lists.macosforge.org/pipermail/macports-dev/2013-May/022735.html
LDFLAGS = '-shared -Wl,-undefined,dynamic_lookup'
else:
# some f2py versions/configurations compile with openMP, so linking against openMP is needed
# to prevent errors during loading of core module
LDFLAGS = ' -openmp -Wl'
#--- run path of for fftw during runtime ----------------------------------------------------------
LDFLAGS += ',-rpath,%s/lib,-rpath,%s/lib64'%(options['FFTW_ROOT'],options['FFTW_ROOT'])
# see http://cens.ioc.ee/pipermail/f2py-users/2003-December/000621.html
if options['IMKL_ROOT']:
if options['F90'] == 'gfortran':
arch = 'gf'
elif options['F90'] == 'ifort':
arch = 'intel'
lib_lapack = '-L%s/lib/intel64 -lmkl_%s_lp64 -lmkl_core -lmkl_sequential -lm'\
%(options['IMKL_ROOT'],arch)
LDFLAGS +=',-rpath,%s/lib/intel64'%(options['IMKL_ROOT'])
elif options['ACML_ROOT'] != '':
lib_lapack = '-L%s/%s64/lib -lacml'%(options['ACML_ROOT'],options['F90'])
LDFLAGS +=',-rpath,%s/%s64/lib'%(options['ACML_ROOT'],options['F90'])
elif options['LAPACK_ROOT'] != '':
lib_lapack = '-L%s/lib -L%s/lib64 -llapack'%(options['LAPACK_ROOT'],options['LAPACK_ROOT'])
LDFLAGS +=',-rpath,%s/lib,-rpath,%s/lib64'%(options['LAPACK_ROOT'],options['LAPACK_ROOT'])
#--------------------------------------------------------------------------------------------------
# f2py does not (yet) support setting of special flags for the linker, hence they must be set via
# environment variable ----------------------------------------------------------------------------
my_env = os.environ
my_env["LDFLAGS"] = LDFLAGS
#--- delete old file ------------------------------------------------------------------------------
try:
os.remove(os.path.join(damaskEnv.relPath('lib/damask'),'core.so'))
except OSError, e:
print ("Error when deleting: %s - %s." % (e.filename,e.strerror))
# The following command is used to compile the fortran files and make the functions defined
# in damask.core.pyf available for python in the module core.so
# It uses the fortran wrapper f2py that is included in the numpy package to construct the
# module core.so out of the fortran code in the f90 files
# For the generation of the pyf file use the following lines:
###########################################################################
#'f2py -h damask.core.pyf' +\
#' --overwrite-signature --no-lower prec.f90 DAMASK_spectral_interface.f90 math.f90 mesh.f90,...'
###########################################################################
os.chdir(codeDir)
cmd = 'f2py damask.core.pyf' +\
' -c --no-lower %s'%(compiler) +\
compileOptions+\
' C_routines.c'+\
' system_routines.f90'+\
' prec.f90'+\
' spectral_interface.f90'+\
' IO.f90'+\
' libs.f90'+\
' numerics.f90'+\
' debug.f90'+\
' math.f90'+\
' FEsolving.f90'+\
' mesh.f90'+\
' core_quit.f90'+\
' -L%s/lib -lfftw3'%(options['FFTW_ROOT'])+\
' %s'%lib_lapack
print('Executing: '+cmd)
try:
subprocess.call(shlex.split(cmd),env=my_env)
except subprocess.CalledProcessError:
print('build failed')
except OSError:
print ('f2py not found')
try:
os.rename(os.path.join(codeDir,'core.so'),\
os.path.join(damaskEnv.relPath('lib/damask'),'core.so'))
except:
pass
modules = glob.glob('*.mod')
for module in modules:
print 'removing', module
os.remove(module)
#--- check if compilation of core module was successful -------------------------------------------
try:
with open(damaskEnv.relPath('lib/damask/core.so')) as f: pass
except IOError as e:
print '*********\n* core.so not found, compilation of core modules was not successful\n*********'

1230
lib/IR_Precision.f90
View File

File diff suppressed because it is too large Load Diff

909
lib/Lib_Base64.f90
View File

@ -1,909 +0,0 @@
!> @ingroup Library
!> @{
!> @defgroup Lib_Base64Library Lib_Base64
!> base64 encoding/decoding library
!> @}
!> @ingroup Interface
!> @{
!> @defgroup Lib_Base64Interface Lib_Base64
!> base64 encoding/decoding library
!> @}
!> @ingroup PublicProcedure
!> @{
!> @defgroup Lib_Base64PublicProcedure Lib_Base64
!> base64 encoding/decoding library
!> @}
!> @ingroup PrivateProcedure
!> @{
!> @defgroup Lib_Base64PrivateProcedure Lib_Base64
!> base64 encoding/decoding library
!> @}
!> @ingroup GlobalVarPar
!> @{
!> @defgroup Lib_Base64GlobalVarPar Lib_Base64
!> base64 encoding/decoding library
!> @}
!> @ingroup PrivateVarPar
!> @{
!> @defgroup Lib_Base64PrivateVarPar Lib_Base64
!> base64 encoding/decoding library
!> @}
!> This module contains base64 encoding/decoding procedures.
!> @todo \b Decoding: Implement decoding functions.
!> @todo \b DocComplete: Complete the documentation.
!> @ingroup Lib_Base64Library
module Lib_Base64
!-----------------------------------------------------------------------------------------------------------------------------------
USE IR_Precision ! Integers and reals precision definition.
!-----------------------------------------------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------------------------------------------
implicit none
private
public:: b64_encode
!public:: b64_decode
public:: pack_data
!-----------------------------------------------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------------------------------------------
!> @ingroup Lib_Base64GlobalVarPar
!> @{
!> @}
!> @ingroup Lib_Base64PrivateVarPar
!> @{
character(64):: base64="ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/" !< Base64 alphabet.
!> @}
!-----------------------------------------------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------------------------------------------
!> @brief Subroutine for encoding numbers (integer and real) to base64.
!> @ingroup Lib_Base64Interface
interface b64_encode
module procedure b64_encode_R8_a, &
b64_encode_R4_a, &
b64_encode_I8_a, &
b64_encode_I4_a, &
b64_encode_I2_a, &
b64_encode_I1_a
endinterface
!!> @brief Subroutine for decoding numbers (integer and real) from base64.
!!> @ingroup Lib_Base64Interface
!interface b64_decode
! module procedure b64_decode_R8_a, &
! b64_decode_R4_a, &
! b64_decode_I8_a, &
! b64_decode_I4_a, &
! b64_decode_I2_a, &
! b64_decode_I1_a
!endinterface
!> @brief Subroutine for packing different kinds of data into single I1P array. This is useful for encoding different kinds
!> variables into a single stream of bits.
!> @ingroup Lib_Base64Interface
interface pack_data
module procedure pack_data_R8_R4,pack_data_R8_I8,pack_data_R8_I4,pack_data_R8_I2,pack_data_R8_I1, &
pack_data_R4_R8,pack_data_R4_I8,pack_data_R4_I4,pack_data_R4_I2,pack_data_R4_I1, &
pack_data_I8_R8,pack_data_I8_R4,pack_data_I8_I4,pack_data_I8_I2,pack_data_I8_I1, &
pack_data_I4_R8,pack_data_I4_R4,pack_data_I4_I8,pack_data_I4_I2,pack_data_I4_I1, &
pack_data_I2_R8,pack_data_I2_R4,pack_data_I2_I8,pack_data_I2_I4,pack_data_I2_I1, &
pack_data_I1_R8,pack_data_I1_R4,pack_data_I1_I8,pack_data_I1_I4,pack_data_I1_I2
endinterface
!-----------------------------------------------------------------------------------------------------------------------------------
contains
!> @ingroup Lib_Base64PrivateProcedure
!> @{
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R8_R4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R8P), intent(IN):: a1(1:) !< Firs data stream.
real(R4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R8_R4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R8_I8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R8P), intent(IN):: a1(1:) !< First data stream.
integer(I8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R8_I8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R8_I4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R8P), intent(IN):: a1(1:) !< First data stream.
integer(I4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R8_I4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R8_I2(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R8P), intent(IN):: a1(1:) !< First data stream.
integer(I2P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R8_I2
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R8_I1(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R8P), intent(IN):: a1(1:) !< First data stream.
integer(I1P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R8_I1
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R4_R8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R4P), intent(IN):: a1(1:) !< Firs data stream.
real(R8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R4_R8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R4_I8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R4P), intent(IN):: a1(1:) !< First data stream.
integer(I8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R4_I8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R4_I4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R4P), intent(IN):: a1(1:) !< First data stream.
integer(I4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R4_I4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R4_I2(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R4P), intent(IN):: a1(1:) !< First data stream.
integer(I2P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R4_I2
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_R4_I1(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
real(R4P), intent(IN):: a1(1:) !< First data stream.
integer(I1P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_R4_I1
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I8_R8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I8P), intent(IN):: a1(1:) !< First data stream.
real(R8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I8_R8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I8_R4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I8P), intent(IN):: a1(1:) !< First data stream.
real(R4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I8_R4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I8_I4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I8P), intent(IN):: a1(1:) !< First data stream.
integer(I4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I8_I4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I8_I2(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I8P), intent(IN):: a1(1:) !< First data stream.
integer(I2P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I8_I2
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I8_I1(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I8P), intent(IN):: a1(1:) !< First data stream.
integer(I1P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I8_I1
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I4_R8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: a1(1:) !< First data stream.
real(R8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I4_R8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I4_R4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: a1(1:) !< First data stream.
real(R4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I4_R4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I4_I8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: a1(1:) !< First data stream.
integer(I8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I4_I8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I4_I2(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: a1(1:) !< First data stream.
integer(I2P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I4_I2
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I4_I1(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: a1(1:) !< First data stream.
integer(I1P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I4_I1
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I2_R8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I2P), intent(IN):: a1(1:) !< First data stream.
real(R8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I2_R8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I2_R4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I2P), intent(IN):: a1(1:) !< First data stream.
real(R4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I2_R4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I2_I8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I2P), intent(IN):: a1(1:) !< First data stream.
integer(I8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I2_I8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I2_I4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I2P), intent(IN):: a1(1:) !< First data stream.
integer(I4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I2_I4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I2_I1(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I2P), intent(IN):: a1(1:) !< First data stream.
integer(I1P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I2_I1
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I1_R8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I1P), intent(IN):: a1(1:) !< First data stream.
real(R8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I1_R8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I1_R4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I1P), intent(IN):: a1(1:) !< First data stream.
real(R4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I1_R4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I1_I8(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I1P), intent(IN):: a1(1:) !< First data stream.
integer(I8P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I1_I8
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I1_I4(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I1P), intent(IN):: a1(1:) !< First data stream.
integer(I4P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I1_I4
!> @brief Subroutine for packing different kinds of data into single I1P array.
pure subroutine pack_data_I1_I2(a1,a2,packed)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I1P), intent(IN):: a1(1:) !< First data stream.
integer(I2P), intent(IN):: a2(1:) !< Second data stream.
integer(I1P), allocatable, intent(INOUT):: packed(:) !< Packed data into I1P array.
integer(I1P), allocatable:: p1(:) !< Temporary packed data of first stream.
integer(I1P), allocatable:: p2(:) !< Temporary packed data of second stream.
integer(I4P):: np !< Size of temporary packed data.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
np = size(transfer(a1,p1)) ; allocate(p1(1:np)) ; p1 = transfer(a1,p1)
np = size(transfer(a2,p2)) ; allocate(p2(1:np)) ; p2 = transfer(a2,p2)
if (allocated(packed)) deallocate(packed) ; allocate(packed(1:size(p1,dim=1)+size(p2,dim=1))) ; packed = [p1,p2]
deallocate(p1,p2)
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine pack_data_I1_I2
!> @brief Subroutine for encoding bits (must be multiple of 24 bits) into base64 charcaters code (of length multiple of 4).
!> @note The bits stream are encoded in chunks of 24 bits as the following example (in little endian order):
!> @code
!> +--first octet--+-second octet--+--third octet--+
!> |7 6 5 4 3 2 1 0|7 6 5 4 3 2 1 0|7 6 5 4 3 2 1 0|
!> +-----------+---+-------+-------+---+-----------+
!> |5 4 3 2 1 0|5 4 3 2 1 0|5 4 3 2 1 0|5 4 3 2 1 0|
!> +--1.index--+--2.index--+--3.index--+--4.index--+
!> @endcode
!> The 4 indexes are stored into 4 elements 8 bits array, thus 2 bits of each array element are not used.
pure subroutine encode_bits(bits,padd,code)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I1P), intent(IN):: bits(1:) !< Bits to be encoded.
integer(I4P), intent(IN):: padd !< Number of padding characters ('=').
character(1), intent(OUT):: code(1:) !< Characters code.
integer(I1P):: sixb(1:4) !< 6 bits slices (stored into 8 bits integer) of 24 bits input.
integer(I8P):: c,e !< Counters.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
c = 1_I8P
do e=1_I8P,size(bits,dim=1,kind=I8P),3_I8P ! loop over array elements: 3 bytes (24 bits) scanning
sixb = 0_I1P
call mvbits(bits(e ),2,6,sixb(1),0)
call mvbits(bits(e ),0,2,sixb(2),4) ; call mvbits(bits(e+1),4,4,sixb(2),0)
call mvbits(bits(e+1),0,4,sixb(3),2) ; call mvbits(bits(e+2),6,2,sixb(3),0)
call mvbits(bits(e+2),0,6,sixb(4),0)
sixb = sixb + 1_I1P
code(c :c )(1:1) = base64(sixb(1):sixb(1))
code(c+1:c+1)(1:1) = base64(sixb(2):sixb(2))
code(c+2:c+2)(1:1) = base64(sixb(3):sixb(3))
code(c+3:c+3)(1:1) = base64(sixb(4):sixb(4))
c = c + 4_I8P
enddo
if (padd>0) code(size(code,dim=1)-padd+1:)(1:1)='='
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine encode_bits
!> @brief Subroutine for encoding array numbers to base64 (R8P).
pure subroutine b64_encode_R8_a(nB,n,code)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: nB !< Number of bytes of single element of n.
real(R8P), intent(IN):: n(1:) !< Array of numbers to be encoded.
character(1), allocatable, intent(OUT):: code(:) !< Encoded array.
integer(I1P):: nI1P(1:((size(n,dim=1)*nB+2)/3)*3) !< One byte integer array containing n.
integer(I4P):: padd !< Number of padding characters ('=').
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
if (allocated(code)) deallocate(code) ; allocate(code(1:((size(n,dim=1)*nB+2)/3)*4)) ! allocating code chars
nI1P = transfer(n,nI1P) ! casting n to integer array of 1 byte elem
padd = mod((size(n,dim=1)*nB),3_I4P) ; if (padd>0_I4P) padd = 3_I4P - padd ! computing the number of padding characters
call encode_bits(bits=nI1P,padd=padd,code=code) ! encoding bits
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine b64_encode_R8_a
!> @brief Subroutine for encoding array numbers to base64 (R4P).
pure subroutine b64_encode_R4_a(nB,n,code)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: nB !< Number of bytes of single element of n.
real(R4P), intent(IN):: n(1:) !< Array of numbers to be encoded.
character(1), allocatable, intent(OUT):: code(:) !< Encoded array.
integer(I1P):: nI1P(1:((size(n,dim=1)*nB+2)/3)*3) !< One byte integer array containing n.
integer(I4P):: padd !< Number of padding characters ('=').
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
if (allocated(code)) deallocate(code) ; allocate(code(1:((size(n,dim=1)*nB+2)/3)*4)) ! allocating code chars
nI1P = transfer(n,nI1P) ! casting n to integer array of 1 byte elem
padd = mod((size(n,dim=1)*nB),3_I4P) ; if (padd>0_I4P) padd = 3_I4P - padd ! computing the number of padding characters
call encode_bits(bits=nI1P,padd=padd,code=code) ! encoding bits
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine b64_encode_R4_a
!> @brief Subroutine for encoding array numbers to base64 (I8P).
pure subroutine b64_encode_I8_a(nB,n,code)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: nB !< Number of bytes of single element of n.
integer(I8P), intent(IN):: n(1:) !< Array of numbers to be encoded.
character(1), allocatable, intent(OUT):: code(:) !< Encoded array.
integer(I1P):: nI1P(1:((size(n,dim=1)*nB+2)/3)*3) !< One byte integer array containing n.
integer(I4P):: padd !< Number of padding characters ('=').
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
if (allocated(code)) deallocate(code) ; allocate(code(1:((size(n,dim=1)*nB+2)/3)*4)) ! allocating code chars
nI1P = transfer(n,nI1P) ! casting n to integer array of 1 byte elem
padd = mod((size(n,dim=1)*nB),3_I4P) ; if (padd>0_I4P) padd = 3_I4P - padd ! computing the number of padding characters
call encode_bits(bits=nI1P,padd=padd,code=code) ! encoding bits
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine b64_encode_I8_a
!> @brief Subroutine for encoding array numbers to base64 (I4P).
pure subroutine b64_encode_I4_a(nB,n,code)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: nB !< Number of bytes of single element of n.
integer(I4P), intent(IN):: n(1:) !< Array of numbers to be encoded.
character(1), allocatable, intent(OUT):: code(:) !< Encoded array.
integer(I1P):: nI1P(1:((size(n,dim=1)*nB+2)/3)*3) !< One byte integer array containing n.
integer(I4P):: padd !< Number of padding characters ('=').
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
if (allocated(code)) deallocate(code) ; allocate(code(1:((size(n,dim=1)*nB+2)/3)*4)) ! allocating code chars
nI1P = transfer(n,nI1P) ! casting n to integer array of 1 byte elem
padd = mod((size(n,dim=1)*nB),3_I4P) ; if (padd>0_I4P) padd = 3_I4P - padd ! computing the number of padding characters
call encode_bits(bits=nI1P,padd=padd,code=code) ! encoding bits
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine b64_encode_I4_a
!> @brief Subroutine for encoding array numbers to base64 (I2P).
pure subroutine b64_encode_I2_a(nB,n,code)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: nB !< Number of bytes of single element of n.
integer(I2P), intent(IN):: n(1:) !< Array of numbers to be encoded.
character(1), allocatable, intent(OUT):: code(:) !< Encoded array.
integer(I1P):: nI1P(1:((size(n,dim=1)*nB+2)/3)*3) !< One byte integer array containing n.
integer(I4P):: padd !< Number of padding characters ('=').
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
if (allocated(code)) deallocate(code) ; allocate(code(1:((size(n,dim=1)*nB+2)/3)*4)) ! allocating code chars
nI1P = transfer(n,nI1P) ! casting n to integer array of 1 byte elem
padd = mod((size(n,dim=1)*nB),3_I4P) ; if (padd>0_I4P) padd = 3_I4P - padd ! computing the number of padding characters
call encode_bits(bits=nI1P,padd=padd,code=code) ! encoding bits
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine b64_encode_I2_a
!> @brief Subroutine for encoding array numbers to base64 (I1P).
pure subroutine b64_encode_I1_a(nB,n,code)
!---------------------------------------------------------------------------------------------------------------------------------
implicit none
integer(I4P), intent(IN):: nB !< Number of bytes of single element of n.
integer(I1P), intent(IN):: n(1:) !< Array of numbers to be encoded.
character(1), allocatable, intent(OUT):: code(:) !< Encoded array.
integer(I1P):: nI1P(1:((size(n,dim=1)*nB+2)/3)*3) !< One byte integer array containing n.
integer(I4P):: padd !< Number of padding characters ('=').
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
if (allocated(code)) deallocate(code) ; allocate(code(1:((size(n,dim=1)*nB+2)/3)*4)) ! allocating code chars
nI1P = transfer(n,nI1P) ! casting n to integer array of 1 byte elem
padd = mod((size(n,dim=1)*nB),3_I4P) ; if (padd>0_I4P) padd = 3_I4P - padd ! computing the number of padding characters
call encode_bits(bits=nI1P,padd=padd,code=code) ! encoding bits
return
!---------------------------------------------------------------------------------------------------------------------------------
endsubroutine b64_encode_I1_a
!!> @brief Subroutine for decoding array numbers from base64 (R8P).
!pure subroutine b64_decode_R8_a(code,n)
!!--------------------------------------------------------------------------------------------------------------------------------
!implicit none
!real(R8P), intent(OUT):: n(1:) !< Number to be decoded.
!character(ncR8P*size(n,dim=1)), intent(IN):: code !< Encoded number.
!integer(I4P):: c,d !< Counters.
!!--------------------------------------------------------------------------------------------------------------------------------
!!--------------------------------------------------------------------------------------------------------------------------------
!d = 1_I4P
!do c=1,len(code),ncR8P
! call b64_decode_R8_s(code=code(c:c+ncR8P-1),n=n(d))
! d = d + 1_I4P
!enddo
!return
!!--------------------------------------------------------------------------------------------------------------------------------
!endsubroutine b64_decode_R8_a
!> @}
endmodule Lib_Base64

6070
lib/Lib_VTK_IO.f90
View File

File diff suppressed because it is too large Load Diff

1
lib/damask/.gitignore vendored
View File

@ -1,2 +1,3 @@
core.so core.so
corientation.so corientation.so
*.pyx

43
lib/damask/__init__.py
View File

@ -1,7 +1,7 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
"""Main aggregator""" """Main aggregator"""
import sys, os import os
with open(os.path.join(os.path.dirname(__file__),'../../VERSION')) as f: with open(os.path.join(os.path.dirname(__file__),'../../VERSION')) as f:
version = f.readline()[:-1] version = f.readline()[:-1]
@ -10,47 +10,14 @@ from .environment import Environment # noqa
from .asciitable import ASCIItable # noqa from .asciitable import ASCIItable # noqa
from .config import Material # noqa from .config import Material # noqa
from .colormaps import Colormap, Color # noqa from .colormaps import Colormap, Color # noqa
try:
from .corientation import Quaternion, Rodrigues, Symmetry, Orientation # noqa
print "Import Cython version of Orientation module"
except:
from .orientation import Quaternion, Rodrigues, Symmetry, Orientation # noqa from .orientation import Quaternion, Rodrigues, Symmetry, Orientation # noqa
# try:
# from .corientation import Quaternion, Rodrigues, Symmetry, Orientation
# print "Import Cython version of Orientation module"
# except:
# from .orientation import Quaternion, Rodrigues, Symmetry, Orientation
#from .block import Block # only one class #from .block import Block # only one class
from .result import Result # noqa from .result import Result # noqa
from .geometry import Geometry # noqa from .geometry import Geometry # noqa
from .solver import Solver # noqa from .solver import Solver # noqa
from .test import Test # noqa from .test import Test # noqa
from .util import extendableOption # noqa from .util import extendableOption # noqa
try:
from . import core
# cleaning up namespace
###################################################################################################
# capitalize according to convention
core.IO = core.io
core.FEsolving = core.fesolving
core.DAMASK_interface = core.damask_interface
# remove modulePrefix_
core.prec.init = core.prec.prec_init
core.DAMASK_interface.init = core.DAMASK_interface.DAMASK_interface_init
core.IO.init = core.IO.IO_init
core.numerics.init = core.numerics.numerics_init
core.debug.init = core.debug.debug_init
core.math.init = core.math.math_init
core.math.tensorAvg = core.math.math_tensorAvg
core.FEsolving.init = core.FEsolving.FE_init
core.mesh.init = core.mesh.mesh_init
core.mesh.nodesAroundCentres = core.mesh.mesh_nodesAroundCentres
core.mesh.deformedCoordsFFT = core.mesh.mesh_deformedCoordsFFT
core.mesh.volumeMismatch = core.mesh.mesh_volumeMismatch
core.mesh.shapeMismatch = core.mesh.mesh_shapeMismatch
except (ImportError,AttributeError) as e:
core = None # from http://www.python.org/dev/peps/pep-0008/
if os.path.split(sys.argv[0])[1] not in ('symLink_Processing.py',
'compile_CoreModule.py',
):
sys.stderr.write('\nWARNING: Core module (Fortran code) not available, \n'\
"try to run 'make processing'\n"\
'Error message when importing core.so: %s\n\n'%e)

1277
lib/damask/corientation.pyx
View File

File diff suppressed because it is too large Load Diff

1226
lib/fftw3.f03
View File

File diff suppressed because it is too large Load Diff

433
processing/post/3Dvisualize.py
View File

@ -1,433 +0,0 @@
#!/usr/bin/env python2
# -*- coding: UTF-8 no BOM -*-
import os,sys,re,fnmatch,vtk
import numpy as np
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
def outStdout(cmd,locals):
if cmd[0:3] == '(!)':
exec(cmd[3:])
elif cmd[0:3] == '(?)':
cmd = eval(cmd[3:])
print cmd
else:
print cmd
return
def outFile(cmd,locals):
if cmd[0:3] == '(!)':
exec(cmd[3:])
elif cmd[0:3] == '(?)':
cmd = eval(cmd[3:])
locals['filepointer'].write(cmd+'\n')
else:
locals['filepointer'].write(cmd+'\n')
return
def output(cmds,locals,dest):
for cmd in cmds:
if isinstance(cmd,list):
output(cmd,locals,dest)
else:
{\
'File': outFile,\
'Stdout': outStdout,\
}[dest](str(cmd),locals)
return
def transliterateToFloat(x):
try:
return float(x)
except:
return 0.0
def unravel(item):
if hasattr(item,'__contains__'): return ' '.join(map(unravel,item))
else: return str(item)
# ++++++++++++++++++++++++++++++++++++++++++++++++++++
def vtk_writeASCII_mesh(mesh,data,res,sep):
"""function writes data array defined on a hexahedral mesh (geometry)"""
info = {\
'tensor': {'name':'tensor','len':9},\
'vector': {'name':'vector','len':3},\
'scalar': {'name':'scalar','len':1},\
'double': {'name':'scalar','len':2},\
'triple': {'name':'scalar','len':3},\
'quadruple': {'name':'scalar','len':4},\
}
N1 = (res[0]+1)*(res[1]+1)*(res[2]+1)
N = res[0]*res[1]*res[2]
cmds = [\
'# vtk DataFile Version 3.1',
'powered by %s'%scriptID,
'ASCII',
'DATASET UNSTRUCTURED_GRID',
'POINTS %i double'%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)],
'CELLS %i %i'%(N,N*9),
]
# cells
for z in range (res[2]):
for y in range (res[1]):
for x in range (res[0]):
base = z*(res[1]+1)*(res[0]+1)+y*(res[0]+1)+x
cmds.append('8 '+'\t'.join(map(str,[ \
base,
base+1,
base+res[0]+2,
base+res[0]+1,
base+(res[1]+1)*(res[0]+1),
base+(res[1]+1)*(res[0]+1)+1,
base+(res[1]+1)*(res[0]+1)+res[0]+2,
base+(res[1]+1)*(res[0]+1)+res[0]+1,
])))
cmds += [\
'CELL_TYPES %i'%N,
['12']*N,
'CELL_DATA %i'%N,
]
for type in data:
plural = {True:'',False:'S'}[type.lower().endswith('s')]
for item in data[type]['_order_']:
cmds += [\
'%s %s double'%(info[type]['name'].upper()+plural,item),
{True:'LOOKUP_TABLE default',False:''}[info[type]['name'][:3]=='sca'],
[[[sep.join(map(unravel,data[type][item][:,j,k]))] for j in range(res[1])] for k in range(res[2])],
]
return cmds
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def vtk_writeASCII_points(coordinates,data,res,sep):
"""function writes data array defined on a point field"""
N = res[0]*res[1]*res[2]
cmds = [\
'# vtk DataFile Version 3.1',
'powered by %s'%scriptID,
'ASCII',
'DATASET UNSTRUCTURED_GRID',
'POINTS %i double'%N,
[[['\t'.join(map(str,coordinates[i,j,k])) for i in range(res[0])] for j in range(res[1])] for k in range(res[2])],
'CELLS %i %i'%(N,N*2),
['1\t%i'%i for i in range(N)],
'CELL_TYPES %i'%N,
['1']*N,
'POINT_DATA %i'%N,
]
for type in data:
plural = {True:'',False:'S'}[type.lower().endswith('s')]
for item in data[type]:
cmds += [\
'%s %s double'%(type.upper()+plural,item),
{True:'LOOKUP_TABLE default',False:''}[type.lower()[:3]=='sca'],
[[[sep.join(map(unravel,data[type][item][:,j,k]))] for j in range(res[1])] for k in range(res[2])],
]
return cmds
# ----------------------- MAIN -------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog [options] datafile[s]', description = """
Produce VTK file from data field.
Coordinates are taken from (consecutive) x, y, and z columns.
""", version = scriptID)
sepChoices = ['n','t','s']
parser.add_option('-s', '--scalar', dest='scalar', action='extend', metavar = '<string LIST>',
help='list of single scalars to visualize')
parser.add_option( '--double', dest='double', action='extend', metavar = '<string LIST>',
help='list of two scalars to visualize')
parser.add_option( '--triple', dest='triple', action='extend', metavar = '<string LIST>',
help='list of three scalars to visualize')
parser.add_option( '--quadruple', dest='quadruple', action='extend', metavar = '<string LIST>',
help='list of four scalars to visualize')
parser.add_option('-v', '--vector', dest='vector', action='extend', metavar = '<string LIST>',
help='list of vectors to visualize')
parser.add_option('-t', '--tensor', dest='tensor', action='extend', metavar = '<string LIST>',
help='list of tensors to visualize')
parser.add_option('-d', '--deformation', dest='defgrad', metavar = 'string',
help='heading of deformation gradient columns [%default]')
parser.add_option('--reference', dest='undeformed', action='store_true',
help='map results to reference (undeformed) configuration [%default]')
parser.add_option('-c','--cell', dest='cell', action='store_true',
help='data is cell-centered [%default]')
parser.add_option('-p','--vertex', dest='cell', action='store_false',
help='data is vertex-centered')
parser.add_option('--mesh', dest='output_mesh', action='store_true',
help='produce VTK mesh file [%default]')
parser.add_option('--nomesh', dest='output_mesh', action='store_false',
help='omit VTK mesh file')
parser.add_option('--points', dest='output_points', action='store_true',
help='produce VTK points file [%default]')
parser.add_option('--nopoints', dest='output_points', action='store_false',
help='omit VTK points file')
parser.add_option('--separator', dest='separator', type='choice', choices=sepChoices, metavar='string',
help='data separator {%s} [t]'%(' '.join(map(str,sepChoices))))
parser.add_option('--scaling', dest='scaling', action='extend', metavar = '<float LIST>',
help='scaling of fluctuation')
parser.add_option('-u', '--unitlength', dest='unitlength', type='float', metavar = 'float',
help='set unit length for 2D model [%default]')
parser.set_defaults(defgrad = 'f')
parser.set_defaults(separator = 't')
parser.set_defaults(scalar = [])
parser.set_defaults(double = [])
parser.set_defaults(triple = [])
parser.set_defaults(quadruple = [])
parser.set_defaults(vector = [])
parser.set_defaults(tensor = [])
parser.set_defaults(output_mesh = True)
parser.set_defaults(output_points = False)
parser.set_defaults(scaling = [])
parser.set_defaults(undeformed = False)
parser.set_defaults(unitlength = 0.0)
parser.set_defaults(cell = True)
sep = {'n': '\n', 't': '\t', 's': ' '}
(options, args) = parser.parse_args()
options.scaling += [1.0 for i in xrange(max(0,3-len(options.scaling)))]
options.scaling = map(float, options.scaling)
for filename in args:
if not os.path.exists(filename):
continue
file = open(filename)
content = file.readlines()
file.close()
m = re.search('(\d+)\s*head', content[0].lower())
if m is None:
continue
print filename,'\n'
sys.stdout.flush()
headrow = int(m.group(1))
headings = content[headrow].split()
column = {}
matches = {}
maxcol = 0
locol = -1
for col,head in enumerate(headings):
if head == {True:'1_pos',False:'1_nodeinitialcoord'}[options.cell]:
locol = col
maxcol = max(maxcol,col+3)
break
if locol < 0:
print 'missing coordinates..!'
continue
column['tensor'] = {}
matches['tensor'] = {}
for label in [options.defgrad] + options.tensor:
column['tensor'][label] = -1
for col,head in enumerate(headings):
if head == label or head == '1_'+label:
column['tensor'][label] = col
maxcol = max(maxcol,col+9)
matches['tensor'][label] = [label]
break
if not options.undeformed and column['tensor'][options.defgrad] < 0:
print 'missing deformation gradient "%s"..!'%options.defgrad
continue
column['vector'] = {}
matches['vector'] = {}
for label in options.vector:
column['vector'][label] = -1
for col,head in enumerate(headings):
if head == label or head == '1_'+label:
column['vector'][label] = col
maxcol = max(maxcol,col+3)
matches['vector'][label] = [label]
break
for length,what in enumerate(['scalar','double','triple','quadruple']):
column[what] = {}
labels = eval("options.%s"%what)
matches[what] = {}
for col,head in enumerate(headings):
for needle in labels:
if fnmatch.fnmatch(head,needle):
column[what][head] = col
maxcol = max(maxcol,col+1+length)
if needle not in matches[what]:
matches[what][needle] = [head]
else:
matches[what][needle] += [head]
values = np.array(sorted([map(transliterateToFloat,line.split()[:maxcol]) for line in content[headrow+1:]],
key=lambda x:(x[locol+0],x[locol+1],x[locol+2])),'d') # sort with z as fastest and x as slowest index
values2 = np.array([map(transliterateToFloat,line.split()[:maxcol]) for line in content[headrow+1:]],'d') # sort with x as fastest and z as slowest index
N = len(values)
tempGrid = [{},{},{}]
for j in xrange(3):
for i in xrange(N):
tempGrid[j][str(values[i,locol+j])] = True
grid = np.array([len(tempGrid[0]),\
len(tempGrid[1]),\
len(tempGrid[2]),],'i')
dim = np.ones(3)
for i,r in enumerate(grid):
if r > 1:
dim[i] = (max(map(float,tempGrid[i].keys()))-min(map(float,tempGrid[i].keys())))*r/(r-1.0)
if grid[2]==1: # for 2D case set undefined dimension to given unitlength or alternatively give it the length of the smallest element
if options.unitlength == 0.0:
dim[2] = min(dim/grid)
else:
dim[2] = options.unitlength
print dim
if options.undeformed:
Favg = np.eye(3)
else:
Favg = damask.core.math.tensorAvg(
np.reshape(np.transpose(values[:,column['tensor'][options.defgrad]:
column['tensor'][options.defgrad]+9]),
(3,3,grid[0],grid[1],grid[2])))
F = np.reshape(np.transpose(values[:,column['tensor'][options.defgrad]:
column['tensor'][options.defgrad]+9]),
(3,3,grid[0],grid[1],grid[2]))
centroids = damask.core.mesh.deformedCoordsFFT(dim,F,Favg,options.scaling)
nodes = damask.core.mesh.nodesAroundCentres(dim,Favg,centroids)
fields = {\
'tensor': {},\
'vector': {},\
'scalar': {},\
'double': {},\
'triple': {},\
'quadruple': {},\
}
reshape = {\
'tensor': [3,3],\
'vector': [3],\
'scalar': [],\
'double': [2],\
'triple': [3],\
'quadruple': [4],\
}
length = {\
'tensor': 9,\
'vector': 3,\
'scalar': 1,\
'double': 2,\
'triple': 3,\
'quadruple': 4,\
}
# vtk lib out
if False:
points = vtk.vtkPoints()
for z in range (grid[2]+1):
for y in range (grid[1]+1):
for x in range (grid[0]+1):
points.InsertNextPoint(nodes[:,x,y,z])
data=[]
j=0
for datatype in fields.keys():
for what in eval('options.'+datatype):
for label in matches[datatype][what]:
col = column[datatype][label]
if col != -1:
data.append(vtk.vtkFloatArray())
data[j].SetNumberOfComponents(length[datatype])
for i in xrange(grid[2]*grid[1]*grid[0]):
for k in xrange(length[datatype]):
data[j].InsertNextValue(values2[i,col+k])
data[j].SetName(label)
j+=1
if options.output_mesh:
hexs = vtk.vtkCellArray()
i = 0
elems=[]
for z in range (grid[2]):
for y in range (grid[1]):
for x in range (grid[0]):
elems.append(vtk.vtkHexahedron())
base = z*(grid[1]+1)*(grid[0]+1)+y*(grid[0]+1)+x
elems[i].GetPointIds().SetId(0, base)
elems[i].GetPointIds().SetId(1, base+1)
elems[i].GetPointIds().SetId(2, base+grid[0]+2)
elems[i].GetPointIds().SetId(3, base+grid[0]+1)
elems[i].GetPointIds().SetId(4, base+(grid[1]+1)*(grid[0]+1))
elems[i].GetPointIds().SetId(5, base+(grid[1]+1)*(grid[0]+1)+1)
elems[i].GetPointIds().SetId(6, base+(grid[1]+1)*(grid[0]+1)+grid[0]+2)
elems[i].GetPointIds().SetId(7, base+(grid[1]+1)*(grid[0]+1)+grid[0]+1)
hexs.InsertNextCell(elems[i])
i+=1
uGrid = vtk.vtkUnstructuredGrid()
uGrid.SetPoints(points)
i = 0
for z in range (grid[2]):
for y in range (grid[1]):
for x in range (grid[0]):
uGrid.InsertNextCell(elems[i].GetCellType(), elems[i].GetPointIds())
i+=1
for i in xrange(len(data)):
uGrid.GetCellData().AddArray(data[i])
outWriter = vtk.vtkXMLUnstructuredGridWriter()
outWriter.SetDataModeToBinary()
outWriter.SetCompressorTypeToZLib()
(head,tail) = os.path.split(filename)
outWriter.SetFileName(os.path.join(head,'mesh_'+os.path.splitext(tail)[0]+'.vtu'))
outWriter.SetInput(uGrid)
outWriter.Write()
for datatype in fields.keys():
print '\n%s:'%datatype,
fields[datatype]['_order_'] = []
for what in eval('options.'+datatype):
for label in matches[datatype][what]:
col = column[datatype][label]
if col != -1:
print label,
fields[datatype][label] = np.reshape(values[:,col:col+length[datatype]],[grid[0],grid[1],grid[2]]+reshape[datatype])
fields[datatype]['_order_'] += [label]
print '\n'
out = {}
if options.output_mesh: out['mesh'] = vtk_writeASCII_mesh(nodes,fields,grid,sep[options.separator])
if options.output_points: out['points'] = vtk_writeASCII_points(centroids,fields,grid,sep[options.separator])
for what in out.keys():
print what
(head,tail) = os.path.split(filename)
vtk = open(os.path.join(head,what+'_'+os.path.splitext(tail)[0]+'.vtk'), 'w')
output(out[what],{'filepointer':vtk},'File')
vtk.close()
print

297
processing/post/addCompatibilityMismatch.py
View File

@ -1,14 +1,210 @@
#!/usr/bin/env python2 #!/usr/bin/env python2
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
import os,sys import os
import math
import numpy as np import numpy as np
import scipy.ndimage
from optparse import OptionParser from optparse import OptionParser
import damask import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version]) scriptID = ' '.join([scriptName,damask.version])
#--------------------------------------------------------------------------------------------------
def cell2node(cellData,grid):
nodeData = 0.0
datalen = np.array(cellData.shape[3:]).prod()
for i in xrange(datalen):
node = scipy.ndimage.convolve(cellData.reshape(tuple(grid[::-1])+(datalen,))[...,i],
np.ones((2,2,2))/8., # 2x2x2 neighborhood of cells
mode = 'wrap',
origin = -1, # offset to have cell origin as center
) # now averaged at cell origins
node = np.append(node,node[np.newaxis,0,:,:,...],axis=0) # wrap along z
node = np.append(node,node[:,0,np.newaxis,:,...],axis=1) # wrap along y
node = np.append(node,node[:,:,0,np.newaxis,...],axis=2) # wrap along x
nodeData = node[...,np.newaxis] if i==0 else np.concatenate((nodeData,node[...,np.newaxis]),axis=-1)
return nodeData
#--------------------------------------------------------------------------------------------------
def deformationAvgFFT(F,grid,size,nodal=False,transformed=False):
"""calculate average cell center (or nodal) deformation for deformation gradient field specified in each grid cell"""
if nodal:
x, y, z = np.meshgrid(np.linspace(0,size[2],1+grid[2]),
np.linspace(0,size[1],1+grid[1]),
np.linspace(0,size[0],1+grid[0]),
indexing = 'ij')
else:
x, y, z = np.meshgrid(np.linspace(size[2]/grid[2]/2.,size[2]-size[2]/grid[2]/2.,grid[2]),
np.linspace(size[1]/grid[1]/2.,size[1]-size[1]/grid[1]/2.,grid[1]),
np.linspace(size[0]/grid[0]/2.,size[0]-size[0]/grid[0]/2.,grid[0]),
indexing = 'ij')
origCoords = np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3)
F_fourier = F if transformed else np.fft.rfftn(F,axes=(0,1,2)) # transform or use provided data
Favg = np.real(F_fourier[0,0,0,:,:])/grid.prod() # take zero freq for average
avgDeformation = np.einsum('ml,ijkl->ijkm',Favg,origCoords) # dX = Favg.X
return avgDeformation
#--------------------------------------------------------------------------------------------------
def displacementFluctFFT(F,grid,size,nodal=False,transformed=False):
"""calculate cell center (or nodal) displacement for deformation gradient field specified in each grid cell"""
integrator = 0.5j * size / math.pi
kk, kj, ki = np.meshgrid(np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2])),
np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1])),
np.arange(grid[0]//2+1),
indexing = 'ij')
k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3)
k_sSquared = np.einsum('...l,...l',k_s,k_s)
k_sSquared[0,0,0] = 1.0 # ignore global average frequency
#--------------------------------------------------------------------------------------------------
# integration in Fourier space
displacement_fourier = -np.einsum('ijkml,ijkl,l->ijkm',
F if transformed else np.fft.rfftn(F,axes=(0,1,2)),
k_s,
integrator,
) / k_sSquared[...,np.newaxis]
#--------------------------------------------------------------------------------------------------
# backtransformation to real space
displacement = np.fft.irfftn(displacement_fourier,grid[::-1],axes=(0,1,2))
return cell2node(displacement,grid) if nodal else displacement
def volTetrahedron(coords):
"""
Return the volume of the tetrahedron with given vertices or sides.
Ifvertices are given they must be in a NumPy array with shape (4,3): the
position vectors of the 4 vertices in 3 dimensions; if the six sides are
given, they must be an array of length 6. If both are given, the sides
will be used in the calculation.
This method implements
Tartaglia's formula using the Cayley-Menger determinant:
|0 1 1 1 1 |
|1 0 s1^2 s2^2 s3^2|
288 V^2 = |1 s1^2 0 s4^2 s5^2|
|1 s2^2 s4^2 0 s6^2|
|1 s3^2 s5^2 s6^2 0 |
where s1, s2, ..., s6 are the tetrahedron side lengths.
from http://codereview.stackexchange.com/questions/77593/calculating-the-volume-of-a-tetrahedron
"""
# The indexes of rows in the vertices array corresponding to all
# possible pairs of vertices
vertex_pair_indexes = np.array(((0, 1), (0, 2), (0, 3),
(1, 2), (1, 3), (2, 3)))
# Get all the squares of all side lengths from the differences between
# the 6 different pairs of vertex positions
vertices = np.concatenate((coords[0],coords[1],coords[2],coords[3])).reshape([4,3])
vertex1, vertex2 = vertex_pair_indexes[:,0], vertex_pair_indexes[:,1]
sides_squared = np.sum((vertices[vertex1] - vertices[vertex2])**2,axis=-1)
# Set up the Cayley-Menger determinant
M = np.zeros((5,5))
# Fill in the upper triangle of the matrix
M[0,1:] = 1
# The squared-side length elements can be indexed using the vertex
# pair indices (compare with the determinant illustrated above)
M[tuple(zip(*(vertex_pair_indexes + 1)))] = sides_squared
# The matrix is symmetric, so we can fill in the lower triangle by
# adding the transpose
M = M + M.T
return np.sqrt(np.linalg.det(M) / 288)
def volumeMismatch(size,F,nodes):
"""
calculates the volume mismatch
volume mismatch is defined as the difference between volume of reconstructed
(compatible) cube and determinant of defgrad at the FP
"""
coords = np.empty([8,3])
vMismatch = np.empty(grid[::-1])
volInitial = size.prod()/grid.prod()
#--------------------------------------------------------------------------------------------------
# calculate actual volume and volume resulting from deformation gradient
for k in xrange(grid[2]):
for j in xrange(grid[1]):
for i in xrange(grid[0]):
coords[0,0:3] = nodes[k, j, i ,0:3]
coords[1,0:3] = nodes[k ,j, i+1,0:3]
coords[2,0:3] = nodes[k ,j+1,i+1,0:3]
coords[3,0:3] = nodes[k, j+1,i ,0:3]
coords[4,0:3] = nodes[k+1,j, i ,0:3]
coords[5,0:3] = nodes[k+1,j, i+1,0:3]
coords[6,0:3] = nodes[k+1,j+1,i+1,0:3]
coords[7,0:3] = nodes[k+1,j+1,i ,0:3]
vMismatch[k,j,i] = \
( abs(volTetrahedron([coords[6,0:3],coords[0,0:3],coords[7,0:3],coords[3,0:3]])) \
+ abs(volTetrahedron([coords[6,0:3],coords[0,0:3],coords[7,0:3],coords[4,0:3]])) \
+ abs(volTetrahedron([coords[6,0:3],coords[0,0:3],coords[2,0:3],coords[3,0:3]])) \
+ abs(volTetrahedron([coords[6,0:3],coords[0,0:3],coords[2,0:3],coords[1,0:3]])) \
+ abs(volTetrahedron([coords[6,0:3],coords[4,0:3],coords[1,0:3],coords[5,0:3]])) \
+ abs(volTetrahedron([coords[6,0:3],coords[4,0:3],coords[1,0:3],coords[0,0:3]]))) \
/np.linalg.det(F[k,j,i,0:3,0:3])
return vMismatch/volInitial
def shapeMismatch(size,F,nodes,centres):
"""
Routine to calculate the shape mismatch
shape mismatch is defined as difference between the vectors from the central point to
the corners of reconstructed (combatible) volume element and the vectors calculated by deforming
the initial volume element with the current deformation gradient
"""
coordsInitial = np.empty([8,3])
sMismatch = np.empty(grid[::-1])
#--------------------------------------------------------------------------------------------------
# initial positions
coordsInitial[0,0:3] = [-size[0]/grid[0],-size[1]/grid[1],-size[2]/grid[2]]
coordsInitial[1,0:3] = [+size[0]/grid[0],-size[1]/grid[1],-size[2]/grid[2]]
coordsInitial[2,0:3] = [+size[0]/grid[0],+size[1]/grid[1],-size[2]/grid[2]]
coordsInitial[3,0:3] = [-size[0]/grid[0],+size[1]/grid[1],-size[2]/grid[2]]
coordsInitial[4,0:3] = [-size[0]/grid[0],-size[1]/grid[1],+size[2]/grid[2]]
coordsInitial[5,0:3] = [+size[0]/grid[0],-size[1]/grid[1],+size[2]/grid[2]]
coordsInitial[6,0:3] = [+size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]]
coordsInitial[7,0:3] = [-size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]]
coordsInitial = coordsInitial/2.0
#--------------------------------------------------------------------------------------------------
# compare deformed original and deformed positions to actual positions
for k in xrange(grid[2]):
for j in xrange(grid[1]):
for i in xrange(grid[0]):
sMismatch[k,j,i] = \
+ np.linalg.norm(nodes[k, j, i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[0,0:3]))\
+ np.linalg.norm(nodes[k, j, i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[1,0:3]))\
+ np.linalg.norm(nodes[k, j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[2,0:3]))\
+ np.linalg.norm(nodes[k, j+1,i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[3,0:3]))\
+ np.linalg.norm(nodes[k+1,j, i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[4,0:3]))\
+ np.linalg.norm(nodes[k+1,j, i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[5,0:3]))\
+ np.linalg.norm(nodes[k+1,j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[6,0:3]))\
+ np.linalg.norm(nodes[k+1,j+1,i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[7,0:3]))
return sMismatch
# -------------------------------------------------------------------- # --------------------------------------------------------------------
# MAIN # MAIN
# -------------------------------------------------------------------- # --------------------------------------------------------------------
@ -64,13 +260,16 @@ for name in filenames:
errors = [] errors = []
remarks = [] remarks = []
if table.label_dimension(options.pos) != 3:
errors.append('coordinates "{}" are not a vector.'.format(options.pos))
else: colCoord = table.label_index(options.pos)
if table.label_dimension(options.defgrad) != 9: if table.label_dimension(options.defgrad) != 9:
errors.append('deformation gradient "{}" is not a tensor.'.format(options.defgrad)) errors.append('deformation gradient "{}" is not a 3x3 tensor.'.format(options.defgrad))
else: colF = table.label_index(options.defgrad)
coordDim = table.label_dimension(options.pos)
if not 3 >= coordDim >= 1:
errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos))
elif coordDim < 3:
remarks.append('appending {} dimension{} to coordinates "{}"...'.format(3-coordDim,
's' if coordDim < 2 else '',
options.pos))
if remarks != []: damask.util.croak(remarks) if remarks != []: damask.util.croak(remarks)
if errors != []: if errors != []:
@ -78,64 +277,60 @@ for name in filenames:
table.close(dismiss=True) table.close(dismiss=True)
continue continue
# ------------------------------------------ assemble header --------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
if options.shape: table.labels_append('shapeMismatch({})'.format(options.defgrad))
if options.volume: table.labels_append('volMismatch({})'.format(options.defgrad))
# --------------- figure out size and grid --------------------------------------------------------- # --------------- figure out size and grid ---------------------------------------------------------
table.data_readArray() table.data_readArray([options.defgrad,options.pos])
table.data_rewind()
coords = [np.unique(table.data[:,colCoord+i]) for i in xrange(3)] if len(table.data.shape) < 2: table.data.shape += (1,) # expand to 2D shape
if table.data[:,9:].shape[1] < 3:
table.data = np.hstack((table.data,
np.zeros((table.data.shape[0],
3-table.data[:,9:].shape[1]),dtype='f'))) # fill coords up to 3D with zeros
coords = [np.unique(table.data[:,9+i]) for i in xrange(3)]
mincorner = np.array(map(min,coords)) mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords)) maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords),'i') grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1) size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # grid==1 spacing set to smallest among other ones size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 set to smallest among other spacings
N = grid.prod() N = grid.prod()
# --------------- figure out columns to process --------------------------------------------------- if N != len(table.data): errors.append('data count {} does not match grid {}x{}x{}.'.format(N,*grid))
key = '1_'+options.defgrad if errors != []:
if table.label_index(key) == -1: damask.util.croak(errors)
damask.util.croak('column "{}" not found...'.format(key)) table.close(dismiss = True)
continue continue
else:
column = table.label_index(key) # remember columns of requested data
# ------------------------------------------ assemble header --------------------------------------- # -----------------------------process data and assemble header -------------------------------------
if options.shape: table.labels_append(['shapeMismatch({})'.format(options.defgrad)])
if options.volume: table.labels_append(['volMismatch({})'.format(options.defgrad)]) F_fourier = np.fft.rfftn(table.data[:,:9].reshape(grid[2],grid[1],grid[0],3,3),axes=(0,1,2)) # perform transform only once...
nodes = displacementFluctFFT(F_fourier,grid,size,True,transformed=True)\
+ deformationAvgFFT (F_fourier,grid,size,True,transformed=True)
if options.shape:
table.labels_append(['shapeMismatch({})'.format(options.defgrad)])
centres = displacementFluctFFT(F_fourier,grid,size,False,transformed=True)\
+ deformationAvgFFT (F_fourier,grid,size,False,transformed=True)
if options.volume:
table.labels_append(['volMismatch({})'.format(options.defgrad)])
table.head_write() table.head_write()
if options.shape:
shapeMismatch = shapeMismatch( size,table.data[:,:9].reshape(grid[2],grid[1],grid[0],3,3),nodes,centres)
if options.volume:
volumeMismatch = volumeMismatch(size,table.data[:,:9].reshape(grid[2],grid[1],grid[0],3,3),nodes)
# ------------------------------------------ read deformation gradient field ----------------------- # ------------------------------------------ output data -------------------------------------------
table.data_rewind() for i in xrange(grid[2]):
F = np.zeros(N*9,'d').reshape([3,3]+list(grid)) for j in xrange(grid[1]):
idx = 0 for k in xrange(grid[0]):
while table.data_read(): table.data_read()
(x,y,z) = damask.util.gridLocation(idx,grid) # figure out (x,y,z) position from line count if options.shape: table.data_append(shapeMismatch[i,j,k])
idx += 1 if options.volume: table.data_append(volumeMismatch[i,j,k])
F[0:3,0:3,x,y,z] = np.array(map(float,table.data[column:column+9]),'d').reshape(3,3) table.data_write()
Favg = damask.core.math.tensorAvg(F)
centres = damask.core.mesh.deformedCoordsFFT(size,F,Favg,[1.0,1.0,1.0])
nodes = damask.core.mesh.nodesAroundCentres(size,Favg,centres)
if options.shape: shapeMismatch = damask.core.mesh.shapeMismatch( size,F,nodes,centres)
if options.volume: volumeMismatch = damask.core.mesh.volumeMismatch(size,F,nodes)
# ------------------------------------------ process data ------------------------------------------
table.data_rewind()
idx = 0
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
(x,y,z) = damask.util.gridLocation(idx,grid) # figure out (x,y,z) position from line count
idx += 1
if options.shape: table.data_append( shapeMismatch[x,y,z])
if options.volume: table.data_append(volumeMismatch[x,y,z])
outputAlive = table.data_write()
# ------------------------------------------ output finalization ----------------------------------- # ------------------------------------------ output finalization -----------------------------------