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

cleaned core module related stuff

This commit is contained in:
Martin Diehl 2016-03-24 11:49:23 +01:00
parent 4592db8dfb
commit 022b089fa7
11 changed files with 9 additions and 10433 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
code/Makefile
View File

@ -350,7 +350,7 @@ DAMASK_spectral.o: INTERFACENAME := spectral_interface.f90
SPECTRAL_SOLVER_FILES = spectral_mech_AL.o spectral_mech_Basic.o spectral_mech_Polarisation.o \
spectral_thermal.o spectral_damage.o
SPECTRAL_FILES = prec.o DAMASK_interface.o IO.o libs.o numerics.o debug.o math.o \
SPECTRAL_FILES = prec.o DAMASK_interface.o IO.o numerics.o debug.o math.o \
FEsolving.o mesh.o material.o lattice.o \
$(SOURCE_FILES) $(KINEMATICS_FILES) $(PLASTIC_FILES) constitutive.o \
crystallite.o \
@ -400,7 +400,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_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 \
$(SOURCE_FILES) $(KINEMATICS_FILES) $(PLASTIC_FILES) constitutive.o \
crystallite.o \
@ -614,9 +614,6 @@ debug.o: debug.f90 \
numerics.o
numerics.o: numerics.f90 \
libs.o
libs.o: libs.f90 \
IO.o
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
!--------------------------------------------------------------------------------------------------
#include "IO.f90"
#include "libs.f90"
#include "numerics.f90"
#include "debug.f90"
#include "math.f90"

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"

263
code/mesh.f90
View File

@ -116,12 +116,8 @@ module mesh
#endif
#ifdef Spectral
#ifdef PETSc
#include <petsc/finclude/petscsys.h>
include 'fftw3-mpi.f03'
#else
include 'fftw3.f03'
#endif
#endif
! 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_ipCoordinates, &
mesh_cellCenterCoordinates, &
mesh_init_postprocessing, &
mesh_get_Ncellnodes, &
mesh_get_unitlength, &
mesh_get_nodeAtIP
#ifdef Spectral
public :: &
mesh_spectral_getGrid, &
mesh_spectral_getSize, &
mesh_nodesAroundCentres, &
mesh_deformedCoordsFFT, &
mesh_volumeMismatch, &
mesh_shapeMismatch
mesh_spectral_getSize
#endif
private :: &
@ -436,8 +427,7 @@ module mesh
mesh_spectral_build_nodes, &
mesh_spectral_build_elements, &
mesh_spectral_build_ipNeighborhood, &
#endif
#ifdef Marc4DAMASK
#elif defined Marc4DAMASK
mesh_marc_get_tableStyles, &
mesh_marc_count_nodesAndElements, &
mesh_marc_count_elementSets, &
@ -448,8 +438,7 @@ module mesh
mesh_marc_build_nodes, &
mesh_marc_count_cpSizes, &
mesh_marc_build_elements, &
#endif
#ifdef Abaqus
#elif defined Abaqus
mesh_abaqus_count_nodesAndElements, &
mesh_abaqus_count_elementSets, &
mesh_abaqus_count_materials, &
@ -473,11 +462,7 @@ module mesh
mesh_tell_statistics, &
FE_mapElemtype, &
mesh_faceMatch, &
mesh_build_FEdata, &
mesh_write_cellGeom, &
mesh_write_elemGeom, &
mesh_write_meshfile, &
mesh_read_meshfile
mesh_build_FEdata
contains
@ -487,9 +472,7 @@ contains
!! Order and routines strongly depend on type of solver
!--------------------------------------------------------------------------------------------------
subroutine mesh_init(ip,el)
#ifdef Spectral
use, intrinsic :: iso_c_binding
#endif
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 IO, only: &
@ -562,25 +545,18 @@ subroutine mesh_init(ip,el)
myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt)
#ifdef Spectral
#ifdef PETSc
call fftw_mpi_init()
#endif
call IO_open_file(FILEUNIT,geometryFile) ! parse info from geometry file...
if (myDebug) write(6,'(a)') ' Opened geometry file'; flush(6)
grid = mesh_spectral_getGrid(fileUnit)
geomSize = mesh_spectral_getSize(fileUnit)
#ifdef PETSc
gridMPI = int(grid,C_INTPTR_T)
alloc_local = fftw_mpi_local_size_3d(gridMPI(3), gridMPI(2), gridMPI(1)/2 +1, &
MPI_COMM_WORLD, local_K, local_K_offset)
grid3 = int(local_K,pInt)
grid3Offset = int(local_K_offset,pInt)
size3 = geomSize(3)*real(grid3,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)
call mesh_spectral_count()
if (myDebug) write(6,'(a)') ' Counted nodes/elements'; flush(6)
@ -592,8 +568,7 @@ subroutine mesh_init(ip,el)
if (myDebug) write(6,'(a)') ' Built nodes'; flush(6)
call mesh_spectral_build_elements(FILEUNIT)
if (myDebug) write(6,'(a)') ' Built elements'; flush(6)
#endif
#ifdef Marc4DAMASK
#elif defined Marc4DAMASK
call IO_open_inputFile(FILEUNIT,modelName) ! parse info from input file...
if (myDebug) write(6,'(a)') ' Opened input file'; flush(6)
call mesh_marc_get_tableStyles(FILEUNIT)
@ -616,8 +591,7 @@ subroutine mesh_init(ip,el)
if (myDebug) write(6,'(a)') ' Counted CP sizes'; flush(6)
call mesh_marc_build_elements(FILEUNIT)
if (myDebug) write(6,'(a)') ' Built elements'; flush(6)
#endif
#ifdef Abaqus
#elif defined Abaqus
call IO_open_inputFile(FILEUNIT,modelName) ! parse info from input file...
if (myDebug) write(6,'(a)') ' Opened input file'; flush(6)
noPart = IO_abaqus_hasNoPart(FILEUNIT)
@ -672,9 +646,6 @@ subroutine mesh_init(ip,el)
if (worldrank == 0_pInt) then
call mesh_tell_statistics
call mesh_write_meshfile
call mesh_write_cellGeom
call mesh_write_elemGeom
endif
if (usePingPong .and. (mesh_Nelems /= mesh_NcpElems)) &
@ -4511,228 +4482,6 @@ 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
#ifdef HDF
use IO, only: &
HDF5_mappingCells
#endif
implicit none
integer(I4P), dimension(1:mesh_Ncells) :: celltype
integer(I4P), dimension(mesh_Ncells*(1_pInt+FE_maxNcellnodesPerCell)) :: cellconnection
#ifdef HDF
integer(pInt), dimension(mesh_Ncells*FE_maxNcellnodesPerCell) :: cellconnectionHDF5
integer(pInt) :: j2=0_pInt
#endif
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
#ifdef HDF
cellconnectionHDF5(j2+1_pInt:j2+FE_NcellnodesPerCell(c)) &
= mesh_cell(1:FE_NcellnodesPerCell(c),i,e)-1_pInt
j2=j2 + FE_ncellnodesPerCell(c)
#endif
enddo
enddo
#ifdef HDF
call HDF5_mappingCells(cellconnectionHDF5(1:j2))
#endif
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
!--------------------------------------------------------------------------------------------------

159
configure vendored
View File

@ -22,19 +22,6 @@ class extendableOption(Option):
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
########################################################
@ -42,35 +29,15 @@ def filePresent(paths,files,warning=False):
parser = OptionParser(option_class=extendableOption, usage='%prog options', description = """
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',
'F90':compiler,
'FFTW_ROOT':'/usr',
'MSC_ROOT' :'/msc',
'DAMASK_NUM_THREADS':4,
'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':{},
}
#--- if local config file exists, read, otherwise assume global config file ------------------------
configFile = 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)
if key == 'DAMASK_BIN':
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:
pass
parser.add_option('--prefix', dest='prefix', metavar='string',
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',
dest='threads', type='int', metavar='int',
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>',
help='options for compilation of spectral solver')
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(marcVersion = defaults['MARC_VERSION'])
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 = [])
(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 ---------------------------------------
output = []
@ -228,12 +91,6 @@ try:
if items[0] == 'DAMASK_BIN':
line = '%s=%s'%(items[0],options.prefix)
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':
line = '%s=%s'%(items[0],options.mscRoot)
options.mscRoot =''
@ -243,14 +100,6 @@ try:
if items[0] == 'DAMASK_NUM_THREADS':
line = '%s=%s'%(items[0],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:
[key,value] = re.split('[= ]',spectralOption)[0:2]
if key == items[0]:
@ -264,18 +113,12 @@ except IOError:
for opt, value in options.__dict__.items():
if opt == 'prefix' and value != 'depending on access rights':
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 != '':
output.append('MSC_ROOT=%s'%value)
if opt == 'marcVersion' and value != '':
output.append('MARC_VERSION=%s'%value)
if opt == 'threads' and 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:
output.append(spectralOption)

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

1226
lib/fftw3.f03
View File

File diff suppressed because it is too large Load Diff

144
processing/post/marc_deformedGeometry.py
View File

@ -1,144 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os,sys,shutil
import numpy as np
import damask
from optparse import OptionParser
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# -----------------------------
# MAIN FUNCTION STARTS HERE
# -----------------------------
# --- input parsing
parser = OptionParser(usage='%prog [options] resultfile', description = """
Create vtk files for the (deformed) geometry that belongs to a .t16 (MSC.Marc) results file.
""", version = scriptID)
parser.add_option('-d','--dir', dest='dir', \
help='name of subdirectory to hold output [%default]')
parser.add_option('-r','--range', dest='range', type='int', nargs=3, \
help='range of positions (or increments) to output (start, end, step) [all]')
parser.add_option('--increments', action='store_true', dest='getIncrements', \
help='switch to increment range [%default]')
parser.add_option('-t','--type', dest='type', type='choice', choices=['ipbased','nodebased'], \
help='processed geometry type [ipbased and nodebased]')
parser.set_defaults(dir = 'vtk')
parser.set_defaults(getIncrements= False)
(options, files) = parser.parse_args()
# --- basic sanity checks
if files == []:
parser.print_help()
parser.error('no file specified...')
filename = os.path.splitext(files[0])[0]
if not os.path.exists(filename+'.t16'):
parser.print_help()
parser.error('invalid file "%s" specified...'%filename+'.t16')
if not options.type :
options.type = ['nodebased', 'ipbased']
else:
options.type = [options.type]
# --- more sanity checks
sys.path.append(damask.solver.Marc().libraryPath('../../'))
try:
import py_post
except:
print('error: no valid Mentat release found')
sys.exit(-1)
# --------------------------- open results file and initialize mesh ----------
p = py_post.post_open(filename+'.t16')
p.moveto(0)
Nnodes = p.nodes()
Nincrements = p.increments() - 1 # t16 contains one "virtual" increment (at 0)
if damask.core.mesh.mesh_init_postprocessing(filename+'.mesh') > 0:
print('error: init not successful')
sys.exit(-1)
Ncellnodes = damask.core.mesh.mesh_get_Ncellnodes()
unitlength = damask.core.mesh.mesh_get_unitlength()
# --------------------------- create output dir --------------------------------
dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir))
if not os.path.isdir(dirname):
os.mkdir(dirname,0755)
# --------------------------- get positions --------------------------------
incAtPosition = {}
positionOfInc = {}
for position in range(Nincrements):
p.moveto(position+1)
incAtPosition[position] = p.increment # remember "real" increment at this position
positionOfInc[p.increment] = position # remember position of "real" increment
if not options.range:
options.getIncrements = False
locations = range(Nincrements) # process all positions
else:
options.range = list(options.range) # convert to list
if options.getIncrements:
locations = [positionOfInc[x] for x in range(options.range[0],options.range[1]+1,options.range[2])
if x in positionOfInc]
else:
locations = range( max(0,options.range[0]),
min(Nincrements,options.range[1]+1),
options.range[2] )
increments = [incAtPosition[x] for x in locations] # build list of increments to process
# --------------------------- loop over positions --------------------------------
for incCount,position in enumerate(locations): # walk through locations
p.moveto(position+1) # wind to correct position
# --- get displacements
node_displacement = [[0,0,0] for i in range(Nnodes)]
for n in range(Nnodes):
if p.node_displacements():
node_displacement[n] = map(lambda x:x*unitlength,list(p.node_displacement(n)))
c = damask.core.mesh.mesh_build_cellnodes(np.array(node_displacement).T,Ncellnodes)
cellnode_displacement = [[c[i][n] for i in range(3)] for n in range(Ncellnodes)]
# --- append displacements to corresponding files
for geomtype in options.type:
outFilename = eval('"'+eval("'%%s_%%s_inc%%0%ii.vtk'%(math.log10(max(increments+[1]))+1)")\
+'"%(dirname + os.sep + os.path.split(filename)[1],geomtype,increments[incCount])')
print outFilename
shutil.copyfile('%s_%s.vtk'%(filename,geomtype),outFilename)
with open(outFilename,'a') as myfile:
myfile.write("POINT_DATA %i\n"%{'nodebased':Nnodes,'ipbased':Ncellnodes}[geomtype])
myfile.write("VECTORS displacement double\n")
coordinates = {'nodebased':node_displacement,'ipbased':cellnode_displacement}[geomtype]
for n in range({'nodebased':Nnodes,'ipbased':Ncellnodes}[geomtype]):
myfile.write("%.8e %.8e %.8e\n"%(coordinates[n][0],coordinates[n][1],coordinates[n][2]))
# --------------------------- DONE --------------------------------

421
processing/post/marc_extractData.py
View File

@ -1,421 +0,0 @@
#!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*-
import os,sys,string,re,time
from optparse import OptionParser, OptionGroup
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# -----------------------------
def ParseOutputFormat(filename,homogID,crystID,phaseID):
"""parse .output* files in order to get a list of outputs"""
myID = {'Homogenization': homogID,
'Crystallite': crystID,
'Constitutive': phaseID,
}
format = {}
for what in ['Homogenization','Crystallite','Constitutive']:
content = []
format[what] = {'outputs':{},'specials':{'brothers':[]}}
for prefix in ['']+map(str,range(1,17)):
if os.path.exists(prefix+filename+'.output'+what):
try:
file = open(prefix+filename+'.output'+what)
content = file.readlines()
file.close()
break
except:
pass
if content == []: continue # nothing found...
tag = ''
tagID = 0
for line in content:
if re.match("\s*$",line) or re.match("#",line): # skip blank lines and comments
continue
m = re.match("\[(.+)\]",line) # look for block indicator
if m: # next section
tag = m.group(1)
tagID += 1
format[what]['specials']['brothers'].append(tag)
if tag == myID[what] or (myID[what].isdigit() and tagID == int(myID[what])):
format[what]['specials']['_id'] = tagID
format[what]['outputs'] = []
tag = myID[what]
else: # data from section
if tag == myID[what]:
(output,length) = line.split()
output.lower()
if length.isdigit():
length = int(length)
if re.match("\((.+)\)",output): # special data, e.g. (Ngrains)
format[what]['specials'][output] = length
elif length > 0:
format[what]['outputs'].append([output,length])
if '_id' not in format[what]['specials']:
print "\nsection '%s' not found in <%s>"%(myID[what], what)
print '\n'.join(map(lambda x:' [%s]'%x, format[what]['specials']['brothers']))
return format
# -----------------------------
def ParsePostfile(p,filename, outputFormat, legacyFormat):
"""
parse postfile in order to get position and labels of outputs
needs "outputFormat" for mapping of output names to postfile output indices
"""
startVar = {True: 'GrainCount',
False:'HomogenizationCount'}
# --- build statistics
stat = { \
'IndexOfLabel': {}, \
'Title': p.title(), \
'Extrapolation': p.extrapolate, \
'NumberOfIncrements': p.increments() - 1, \
'NumberOfNodes': p.nodes(), \
'NumberOfNodalScalars': p.node_scalars(), \
'LabelOfNodalScalar': [None]*p.node_scalars() , \
'NumberOfElements': p.elements(), \
'NumberOfElementalScalars': p.element_scalars(), \
'LabelOfElementalScalar': [None]*p.element_scalars() , \
'NumberOfElementalTensors': p.element_tensors(), \
'LabelOfElementalTensor': [None]*p.element_tensors(), \
}
# --- find labels
for labelIndex in range(stat['NumberOfNodalScalars']):
label = p.node_scalar_label(labelIndex)
stat['IndexOfLabel'][label] = labelIndex
stat['LabelOfNodalScalar'][labelIndex] = label
for labelIndex in range(stat['NumberOfElementalScalars']):
label = p.element_scalar_label(labelIndex)
stat['IndexOfLabel'][label] = labelIndex
stat['LabelOfElementalScalar'][labelIndex] = label
for labelIndex in range(stat['NumberOfElementalTensors']):
label = p.element_tensor_label(labelIndex)
stat['IndexOfLabel'][label] = labelIndex
stat['LabelOfElementalTensor'][labelIndex] = label
if 'User Defined Variable 1' in stat['IndexOfLabel']: # output format without dedicated names?
stat['IndexOfLabel'][startVar[legacyFormat]] = stat['IndexOfLabel']['User Defined Variable 1'] # adjust first named entry
if startVar[legacyFormat] in stat['IndexOfLabel']: # does the result file contain relevant user defined output at all?
startIndex = stat['IndexOfLabel'][startVar[legacyFormat]]
stat['LabelOfElementalScalar'][startIndex] = startVar[legacyFormat]
# We now have to find a mapping for each output label as defined in the .output* files to the output position in the post file
# Since we know where the user defined outputs start ("startIndex"), we can simply assign increasing indices to the labels
# given in the .output* file
offset = 1
if legacyFormat:
stat['LabelOfElementalScalar'][startIndex + offset] = startVar[not legacyFormat] # add HomogenizationCount as second
offset += 1
for (name,N) in outputFormat['Homogenization']['outputs']:
for i in range(N):
label = {False: '%s'%( name),
True:'%i_%s'%(i+1,name)}[N > 1]
stat['IndexOfLabel'][label] = startIndex + offset
stat['LabelOfElementalScalar'][startIndex + offset] = label
offset += 1
if not legacyFormat:
stat['IndexOfLabel'][startVar[not legacyFormat]] = startIndex + offset
stat['LabelOfElementalScalar'][startIndex + offset] = startVar[not legacyFormat] # add GrainCount
offset += 1
if '(ngrains)' in outputFormat['Homogenization']['specials']:
for grain in range(outputFormat['Homogenization']['specials']['(ngrains)']):
stat['IndexOfLabel']['%i_CrystalliteCount'%(grain+1)] = startIndex + offset # report crystallite count
stat['LabelOfElementalScalar'][startIndex + offset] = '%i_CrystalliteCount'%(grain+1) # add GrainCount
offset += 1
for (name,N) in outputFormat['Crystallite']['outputs']: # add crystallite outputs
for i in range(N):
label = {False: '%i_%s'%(grain+1, name),
True:'%i_%i_%s'%(grain+1,i+1,name)}[N > 1]
stat['IndexOfLabel'][label] = startIndex + offset
stat['LabelOfElementalScalar'][startIndex + offset] = label
offset += 1
stat['IndexOfLabel']['%i_ConstitutiveCount'%(grain+1)] = startIndex + offset # report constitutive count
stat['LabelOfElementalScalar'][startIndex + offset] = '%i_ConstitutiveCount'%(grain+1) # add GrainCount
offset += 1
for (name,N) in outputFormat['Constitutive']['outputs']: # add constitutive outputs
for i in range(N):
label = {False: '%i_%s'%(grain+1, name),
True:'%i_%i_%s'%(grain+1,i+1,name)}[N > 1]
stat['IndexOfLabel'][label] = startIndex + offset
try:
stat['LabelOfElementalScalar'][startIndex + offset] = label
except IndexError:
print 'trying to assign %s at position %i+%i'%(label,startIndex,offset)
sys.exit(1)
offset += 1
return stat
# -----------------------------
def GetIncrementLocations(p,Nincrements,options):
"""get mapping between positions in postfile and increment number"""
incAtPosition = {}
positionOfInc = {}
for position in range(Nincrements):
p.moveto(position+1)
incAtPosition[position] = p.increment # remember "real" increment at this position
positionOfInc[p.increment] = position # remember position of "real" increment
if not options.range:
options.getIncrements = False
locations = range(Nincrements) # process all positions
else:
options.range = list(options.range) # convert to list
if options.getIncrements:
locations = [positionOfInc[x] for x in range(options.range[0],options.range[1]+1,options.range[2])
if x in positionOfInc]
else:
locations = range( max(0,options.range[0]),
min(Nincrements,options.range[1]+1),
options.range[2] )
increments = [incAtPosition[x] for x in locations] # build list of increments to process
return [increments,locations]
# -----------------------------
def SummarizePostfile(stat,where=sys.stdout):
where.write('\n\n')
where.write('title:\t%s'%stat['Title'] + '\n\n')
where.write('extraplation:\t%s'%stat['Extrapolation'] + '\n\n')
where.write('increments:\t%i'%(stat['NumberOfIncrements']) + '\n\n')
where.write('nodes:\t%i'%stat['NumberOfNodes'] + '\n\n')
where.write('elements:\t%i'%stat['NumberOfElements'] + '\n\n')
where.write('nodal scalars:\t%i'%stat['NumberOfNodalScalars'] + '\n\n '\
+'\n '.join(stat['LabelOfNodalScalar']) + '\n\n')
where.write('elemental scalars:\t%i'%stat['NumberOfElementalScalars'] + '\n\n '\
+ '\n '.join(stat['LabelOfElementalScalar']) + '\n\n')
where.write('elemental tensors:\t%i'%stat['NumberOfElementalTensors'] + '\n\n '\
+ '\n '.join(stat['LabelOfElementalTensor']) + '\n\n')
return True
# -----------------------------
def SummarizeOutputfile(format,where=sys.stdout):
where.write('\nUser Defined Outputs')
for what in format.keys():
where.write('\n\n %s:'%what)
for output in format[what]['outputs']:
where.write('\n %s'%output)
return True
# -----------------------------
def writeHeader(myfile,stat,geomtype):
myfile.write('2\theader\n')
myfile.write(string.replace('$Id$','\n','\\n')+
'\t' + ' '.join(sys.argv[1:]) + '\n')
if geomtype == 'nodebased':
myfile.write('node')
for i in range(stat['NumberOfNodalScalars']):
myfile.write('\t%s'%''.join(stat['LabelOfNodalScalar'][i].split()))
elif geomtype == 'ipbased':
myfile.write('elem\tip')
for i in range(stat['NumberOfElementalScalars']):
myfile.write('\t%s'%''.join(stat['LabelOfElementalScalar'][i].split()))
myfile.write('\n')
return True
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Extract data from a .t16 (MSC.Marc) results file.
""", version = scriptID)
parser.add_option('-i','--info', action='store_true', dest='info', \
help='list contents of resultfile [%default]')
parser.add_option('-l','--legacy', action='store_true', dest='legacy', \
help='legacy user result block (starts with GrainCount) [%default]')
parser.add_option('-d','--dir', dest='dir', \
help='name of subdirectory to hold output [%default]')
parser.add_option('-r','--range', dest='range', type='int', nargs=3, \
help='range of positions (or increments) to output (start, end, step) [all]')
parser.add_option('--increments', action='store_true', dest='getIncrements', \
help='switch to increment range [%default]')
parser.add_option('-t','--type', dest='type', type='choice', choices=['ipbased','nodebased'], \
help='processed geometry type [ipbased and nodebased]')
group_material = OptionGroup(parser,'Material identifier')
group_material.add_option('--homogenization', dest='homog', \
help='homogenization identifier (as string or integer [%default])', metavar='<ID>')
group_material.add_option('--crystallite', dest='cryst', \
help='crystallite identifier (as string or integer [%default])', metavar='<ID>')
group_material.add_option('--phase', dest='phase', \
help='phase identifier (as string or integer [%default])', metavar='<ID>')
parser.add_option_group(group_material)
parser.set_defaults(info = False)
parser.set_defaults(legacy = False)
parser.set_defaults(dir = 'vtk')
parser.set_defaults(getIncrements= False)
parser.set_defaults(homog = '1')
parser.set_defaults(cryst = '1')
parser.set_defaults(phase = '1')
(options, files) = parser.parse_args()
# --- sanity checks
if files == []:
parser.print_help()
parser.error('no file specified...')
filename = os.path.splitext(files[0])[0]
if not os.path.exists(filename+'.t16'):
parser.print_help()
parser.error('invalid file "%s" specified...'%filename+'.t16')
sys.path.append(damask.solver.Marc().libraryPath('../../'))
try:
import py_post
except:
print('error: no valid Mentat release found')
sys.exit(-1)
if not options.type :
options.type = ['nodebased', 'ipbased']
else:
options.type = [options.type]
# --- initialize mesh data
if damask.core.mesh.mesh_init_postprocessing(filename+'.mesh'):
print('error: init not successful')
sys.exit(-1)
# --- check if ip data available for all elements; if not, then .t19 file is required
p = py_post.post_open(filename+'.t16')
asciiFile = False
p.moveto(1)
for e in range(p.elements()):
if not damask.core.mesh.mesh_get_nodeAtIP(str(p.element(e).type),1):
if os.path.exists(filename+'.t19'):
p.close()
p = py_post.post_open(filename+'.t19')
asciiFile = True
break
# --- parse *.output and *.t16 file
outputFormat = ParseOutputFormat(filename,options.homog,options.cryst,options.phase)
p.moveto(1)
p.extrapolation('translate')
stat = ParsePostfile(p,filename,outputFormat,options.legacy)
# --- output info
if options.info:
print '\n\nMentat release %s'%damask.solver.Marc().version('../../')
SummarizePostfile(stat)
SummarizeOutputfile(outputFormat)
sys.exit(0)
# --- create output dir
dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir))
if not os.path.isdir(dirname):
os.mkdir(dirname,0755)
# --- get positions
[increments,locations] = GetIncrementLocations(p,stat['NumberOfIncrements'],options)
# --- loop over positions
time_start = time.time()
for incCount,position in enumerate(locations): # walk through locations
p.moveto(position+1) # wind to correct position
time_delta = (float(len(locations)) / float(incCount+1) - 1.0) * (time.time() - time_start)
sys.stdout.write("\r(%02i:%02i:%02i) processing increment %i of %i..."\
%(time_delta//3600,time_delta%3600//60,time_delta%60,incCount+1,len(locations)))
sys.stdout.flush()
# --- write header
outFilename = {}
for geomtype in options.type:
outFilename[geomtype] = eval('"'+eval("'%%s_%%s_inc%%0%ii.txt'%(math.log10(max(increments+[1]))+1)")\
+'"%(dirname + os.sep + os.path.split(filename)[1],geomtype,increments[incCount])')
with open(outFilename[geomtype],'w') as myfile:
writeHeader(myfile,stat,geomtype)
# --- write node based data
if geomtype == 'nodebased':
for n in range(stat['NumberOfNodes']):
myfile.write(str(n))
for l in range(stat['NumberOfNodalScalars']):
myfile.write('\t'+str(p.node_scalar(n,l)))
myfile.write('\n')
# --- write ip based data
elif geomtype == 'ipbased':
for e in range(stat['NumberOfElements']):
if asciiFile:
print 'ascii postfile not yet supported'
sys.exit(-1)
else:
ipData = [[]]
for l in range(stat['NumberOfElementalScalars']):
data = p.element_scalar(e,l)
for i in range(len(data)): # at least as many nodes as ips
node = damask.core.mesh.mesh_get_nodeAtIP(str(p.element(e).type),i+1) # fortran indexing starts at 1
if not node: break # no more ips
while i >= len(ipData): ipData.append([])
ipData[i].extend([data[node-1].value]) # python indexing starts at 0
for i in range(len(ipData)):
myfile.write('\t'.join(map(str,[e,i]+ipData[i]))+'\n')
p.close()
sys.stdout.write("\n")