Merge branch 'development' into MiscImprovements

This commit is contained in:
Martin Diehl 2020-01-24 20:08:14 +01:00
commit bfe8585903
24 changed files with 247 additions and 920 deletions

1
.gitignore vendored
View File

@ -3,6 +3,7 @@
*.xdmf *.xdmf
*.bak *.bak
*~ *~
.DS_Store
bin bin
PRIVATE PRIVATE
build build

View File

@ -35,4 +35,3 @@ clean:
.PHONY: processing .PHONY: processing
processing: processing:
@./installation/symlink_Processing.py ${MAKEFLAGS} @./installation/symlink_Processing.py ${MAKEFLAGS}

@ -1 +1 @@
Subproject commit 4a8f3ccc2f9241b4afceece57ce71cae8025faa3 Subproject commit 83cb0b6ceea2e26c554ad4c6fbe6aabaee6fa27f

View File

@ -1 +1 @@
v2.0.3-1520-g701a9b18 v2.0.3-1531-g3d6ec695

View File

@ -103,9 +103,9 @@ set (COMPILE_FLAGS "${COMPILE_FLAGS} -ffpe-summary=all")
# print summary of floating point exeptions (invalid,zero,overflow,underflow,inexact,denormal) # print summary of floating point exeptions (invalid,zero,overflow,underflow,inexact,denormal)
# Additional options # Additional options
# -Warray-temporarieswarnings: because we have many temporary arrays (performance issue?): # -Warray-temporarieswarnings: because we have many temporary arrays (performance issue?)
# -Wimplicit-interface: no interfaces for lapack/MPI routines # -Wimplicit-interface: no interfaces for lapack/MPI routines
# -Wunsafe-loop-optimizations: warn if the loop cannot be optimized due to nontrivial assumptions. # -Wunsafe-loop-optimizations: warn if the loop cannot be optimized due to nontrivial assumptions
#------------------------------------------------------------------------------------------------ #------------------------------------------------------------------------------------------------
# Runtime debugging # Runtime debugging
@ -122,6 +122,7 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fcheck=all")
# Additional options # Additional options
# -ffpe-trap=precision,denormal,underflow # -ffpe-trap=precision,denormal,underflow
#------------------------------------------------------------------------------------------------ #------------------------------------------------------------------------------------------------
# precision settings # precision settings
set (PRECISION_FLAGS "${PRECISION_FLAGS} -fdefault-real-8") set (PRECISION_FLAGS "${PRECISION_FLAGS} -fdefault-real-8")

View File

@ -2,23 +2,51 @@
# PGI Compiler # PGI Compiler
################################################################################################### ###################################################################################################
if (OPENMP)
set (OPENMP_FLAGS "-mp")
else ()
set (OPENMP_FLAGS "-nomp")
endif ()
if (OPTIMIZATION STREQUAL "OFF") if (OPTIMIZATION STREQUAL "OFF")
set (OPTIMIZATION_FLAGS "-O0" ) set (OPTIMIZATION_FLAGS "-O0" )
elseif (OPTIMIZATION STREQUAL "DEFENSIVE") elseif (OPTIMIZATION STREQUAL "DEFENSIVE")
set (OPTIMIZATION_FLAGS "-O2") set (OPTIMIZATION_FLAGS "-O2 -fast")
elseif (OPTIMIZATION STREQUAL "AGGRESSIVE") elseif (OPTIMIZATION STREQUAL "AGGRESSIVE")
set (OPTIMIZATION_FLAGS "-O3") set (OPTIMIZATION_FLAGS "-O4 -fast -Mvect=sse")
endif () endif ()
set (STANDARD_CHECK "-Mallocatable=03 -Mstandard")
#------------------------------------------------------------------------------------------------ #------------------------------------------------------------------------------------------------
# Fine tuning compilation options # Fine tuning compilation options
set (COMPILE_FLAGS "${COMPILE_FLAGS} -Mpreprocess") set (COMPILE_FLAGS "${COMPILE_FLAGS} -Mpreprocess")
# preprocessor # preprocessor
set (STANDARD_CHECK "-Mallocatable=03") set (COMPILE_FLAGS "${COMPILE_FLAGS} -Minfo=all")
# instructs the compiler to produce information on standard error
#------------------------------------------------------------------------------------------------ set (COMPILE_FLAGS "${COMPILE_FLAGS} -Minform=warn")
# instructs the compiler to display error messages at the specified and higher levels
set (COMPILE_FLAGS "${COMPILE_FLAGS} -Mdclchk")
# instructs the compiler to require that all program variables be declared
#------------------------------------------------------------------------------------------------O
# Runtime debugging # Runtime debugging
set (DEBUG_FLAGS "${DEBUG_FLAGS} -g") set (DEBUG_FLAGS "${DEBUG_FLAGS} -g")
# Includes debugging information in the object module; sets the optimization level to zero unless a -O option is present on the command line # Includes debugging information in the object module; sets the optimization level to zero unless a -O option is present on the command line
set (DEBUG_FLAGS "${DEBUG_FLAGS} -C")
# Generates code to check array bounds
set (DEBUG_FLAGS "${DEBUG_FLAGS} -Mchkptr")
# Check for NULL pointers (pgf95, pgfortran only)
set (DEBUG_FLAGS "${DEBUG_FLAGS} -Mchkstk")
# Check the stack for available space upon entry to and before the start of a parallel region. Useful when many private variables are declared
set (DEBUG_FLAGS "${DEBUG_FLAGS} -Mbounds")
# Specifies whether array bounds checking is enabled or disabled
#------------------------------------------------------------------------------------------------
# precision settings
set (PRECISION_FLAGS "${PRECISION_FLAGS} -r8")
# Determines whether the compiler promotes REAL variables and constants to DOUBLE PRECISION

View File

11
env/DAMASK.csh vendored
View File

@ -2,15 +2,14 @@
# usage: source DAMASK_env.csh # usage: source DAMASK_env.csh
set CALLED=($_) set CALLED=($_)
set DIRNAME=`dirname $CALLED[2]` set ENV_ROOT=`dirname $CALLED[2]`
set DAMASK_ROOT=`python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" $DIRNAME"/../"` set DAMASK_ROOT=`python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" $ENV_ROOT"/../"`
source $DAMASK_ROOT/CONFIG source $ENV_ROOT/CONFIG
set path = ($DAMASK_ROOT/bin $path) set path = ($DAMASK_ROOT/bin $path)
set SOLVER=`which DAMASK_spectral` set SOLVER=`which DAMASK_spectral`
set PROCESSING=`which postResults`
if ( "x$DAMASK_NUM_THREADS" == "x" ) then if ( "x$DAMASK_NUM_THREADS" == "x" ) then
set DAMASK_NUM_THREADS=1 set DAMASK_NUM_THREADS=1
endif endif
@ -32,7 +31,6 @@ if ( $?prompt ) then
echo Using environment with ... echo Using environment with ...
echo "DAMASK $DAMASK_ROOT" echo "DAMASK $DAMASK_ROOT"
echo "Grid Solver $SOLVER" echo "Grid Solver $SOLVER"
echo "Post Processing $PROCESSING"
if ( $?PETSC_DIR) then if ( $?PETSC_DIR) then
echo "PETSc location $PETSC_DIR" echo "PETSc location $PETSC_DIR"
endif endif
@ -52,3 +50,6 @@ if ( ! $?PYTHONPATH ) then
else else
setenv PYTHONPATH $DAMASK_ROOT/python:$PYTHONPATH setenv PYTHONPATH $DAMASK_ROOT/python:$PYTHONPATH
endif endif
setenv MSC_ROOT
setenv ABAQUS_VERSION
setenv MARC_VERSION

22
env/DAMASK.sh vendored
View File

@ -10,14 +10,14 @@ function blink {
} }
if [ "$OSTYPE" == "linux-gnu" ] || [ "$OSTYPE" == 'linux' ]; then if [ "$OSTYPE" == "linux-gnu" ] || [ "$OSTYPE" == 'linux' ]; then
DAMASK_ROOT=$(dirname $BASH_SOURCE) ENV_ROOT=$(dirname $BASH_SOURCE)
else else
[[ "${BASH_SOURCE::1}" == "/" ]] && BASE="" || BASE="$(pwd)/" [[ "${BASH_SOURCE::1}" == "/" ]] && BASE="" || BASE="$(pwd)/"
STAT=$(stat "$(dirname $BASE$BASH_SOURCE)") STAT=$(stat "$(dirname $BASE$BASH_SOURCE)")
DAMASK_ROOT=${STAT##* } ENV_ROOT=${STAT##* }
fi fi
DAMASK_ROOT=$(canonicalPath "$DAMASK_ROOT/../") DAMASK_ROOT=$(canonicalPath "$ENV_ROOT/../")
# shorthand command to change to DAMASK_ROOT directory # shorthand command to change to DAMASK_ROOT directory
@ -27,7 +27,7 @@ eval "function DAMASK_root() { cd $DAMASK_ROOT; }"
set() { set() {
export $1$2$3 export $1$2$3
} }
source $DAMASK_ROOT/CONFIG source $ENV_ROOT/CONFIG
unset -f set unset -f set
# add BRANCH if DAMASK_ROOT is a git repository # add BRANCH if DAMASK_ROOT is a git repository
@ -38,9 +38,6 @@ PATH=${DAMASK_ROOT}/bin:$PATH
SOLVER=$(type -p DAMASK_spectral || true 2>/dev/null) SOLVER=$(type -p DAMASK_spectral || true 2>/dev/null)
[ "x$SOLVER" == "x" ] && SOLVER=$(blink 'Not found!') [ "x$SOLVER" == "x" ] && SOLVER=$(blink 'Not found!')
PROCESSING=$(type -p postResults || true 2>/dev/null)
[ "x$PROCESSING" == "x" ] && PROCESSING=$(blink 'Not found!')
[ "x$DAMASK_NUM_THREADS" == "x" ] && DAMASK_NUM_THREADS=1 [ "x$DAMASK_NUM_THREADS" == "x" ] && DAMASK_NUM_THREADS=1
# currently, there is no information that unlimited stack size causes problems # currently, there is no information that unlimited stack size causes problems
@ -60,7 +57,6 @@ if [ ! -z "$PS1" ]; then
echo Using environment with ... echo Using environment with ...
echo "DAMASK $DAMASK_ROOT $BRANCH" echo "DAMASK $DAMASK_ROOT $BRANCH"
echo "Grid Solver $SOLVER" echo "Grid Solver $SOLVER"
echo "Post Processing $PROCESSING"
if [ "x$PETSC_DIR" != "x" ]; then if [ "x$PETSC_DIR" != "x" ]; then
echo -n "PETSc location " echo -n "PETSc location "
[ -d $PETSC_DIR ] && echo $PETSC_DIR || blink $PETSC_DIR [ -d $PETSC_DIR ] && echo $PETSC_DIR || blink $PETSC_DIR
@ -93,12 +89,8 @@ fi
export DAMASK_NUM_THREADS export DAMASK_NUM_THREADS
export PYTHONPATH=$DAMASK_ROOT/python:$PYTHONPATH export PYTHONPATH=$DAMASK_ROOT/python:$PYTHONPATH
for var in BASE STAT SOLVER PROCESSING BRANCH; do for var in BASE STAT SOLVER BRANCH; do
unset "${var}" unset "${var}"
done done
for var in DAMASK MSC; do unset "ENV_ROOT"
unset "${var}_ROOT" unset "DAMASK_ROOT"
done
for var in ABAQUS MARC; do
unset "${var}_VERSION"
done

17
env/DAMASK.zsh vendored
View File

@ -9,6 +9,7 @@ function blink {
echo -e "\033[2;5m$1\033[0m" echo -e "\033[2;5m$1\033[0m"
} }
ENV_ROOT=$(canonicalPath "${0:a:h}")
DAMASK_ROOT=$(canonicalPath "${0:a:h}'/..") DAMASK_ROOT=$(canonicalPath "${0:a:h}'/..")
# shorthand command to change to DAMASK_ROOT directory # shorthand command to change to DAMASK_ROOT directory
@ -18,7 +19,7 @@ eval "function DAMASK_root() { cd $DAMASK_ROOT; }"
set() { set() {
export $1$2$3 export $1$2$3
} }
source $DAMASK_ROOT/CONFIG source $ENV_ROOT/CONFIG
unset -f set unset -f set
# add BRANCH if DAMASK_ROOT is a git repository # add BRANCH if DAMASK_ROOT is a git repository
@ -29,9 +30,6 @@ PATH=${DAMASK_ROOT}/bin:$PATH
SOLVER=$(which DAMASK_spectral || true 2>/dev/null) SOLVER=$(which DAMASK_spectral || true 2>/dev/null)
[[ "x$SOLVER" == "x" ]] && SOLVER=$(blink 'Not found!') [[ "x$SOLVER" == "x" ]] && SOLVER=$(blink 'Not found!')
PROCESSING=$(which postResults || true 2>/dev/null)
[[ "x$PROCESSING" == "x" ]] && PROCESSING=$(blink 'Not found!')
[[ "x$DAMASK_NUM_THREADS" == "x" ]] && DAMASK_NUM_THREADS=1 [[ "x$DAMASK_NUM_THREADS" == "x" ]] && DAMASK_NUM_THREADS=1
# currently, there is no information that unlimited stack size causes problems # currently, there is no information that unlimited stack size causes problems
@ -51,7 +49,6 @@ if [ ! -z "$PS1" ]; then
echo "Using environment with ..." echo "Using environment with ..."
echo "DAMASK $DAMASK_ROOT $BRANCH" echo "DAMASK $DAMASK_ROOT $BRANCH"
echo "Grid Solver $SOLVER" echo "Grid Solver $SOLVER"
echo "Post Processing $PROCESSING"
if [ "x$PETSC_DIR" != "x" ]; then if [ "x$PETSC_DIR" != "x" ]; then
echo -n "PETSc location " echo -n "PETSc location "
[ -d $PETSC_DIR ] && echo $PETSC_DIR || blink $PETSC_DIR [ -d $PETSC_DIR ] && echo $PETSC_DIR || blink $PETSC_DIR
@ -86,12 +83,8 @@ fi
export DAMASK_NUM_THREADS export DAMASK_NUM_THREADS
export PYTHONPATH=$DAMASK_ROOT/python:$PYTHONPATH export PYTHONPATH=$DAMASK_ROOT/python:$PYTHONPATH
for var in SOLVER PROCESSING BRANCH; do for var in SOLVER BRANCH; do
unset "${var}" unset "${var}"
done done
for var in DAMASK MSC; do unset "ENV_ROOT"
unset "${var}_ROOT" unset "DAMASK_ROOT"
done
for var in ABAQUS MARC; do
unset "${var}_VERSION"
done

View File

@ -2,11 +2,6 @@
SCRIPTLOCATION="$( cd "$( dirname "$0" )" && pwd )" SCRIPTLOCATION="$( cd "$( dirname "$0" )" && pwd )"
DAMASK_ROOT=$SCRIPTLOCATION/../../ DAMASK_ROOT=$SCRIPTLOCATION/../../
# defining set() allows to source the same file for tcsh and bash, with and without space around =
set() {
export $1$2$3
}
source $DAMASK_ROOT/CONFIG
if [ "x$MSC_ROOT" != "x" ]; then if [ "x$MSC_ROOT" != "x" ]; then
DEFAULT_DIR=$MSC_ROOT DEFAULT_DIR=$MSC_ROOT
@ -15,7 +10,7 @@ if [ "x$MARC_VERSION" != "x" ]; then
DEFAULT_VERSION=$MARC_VERSION DEFAULT_VERSION=$MARC_VERSION
fi fi
if [ "x$DAMASK_BIN" != "x" ]; then if [ "x$DAMASK_BIN" != "x" ]; then
BIN_DIR=$DAMASK_BIN BIN_DIR=$DAMASK_ROOT/bin
fi fi
while [ ! -d "$SCRIPTLOCATION/$VERSION" ] || [ -z "$VERSION" ] while [ ! -d "$SCRIPTLOCATION/$VERSION" ] || [ -z "$VERSION" ]

View File

@ -1 +1,2 @@
include damask/VERSION include damask/VERSION
include damask/LICENSE

13
python/README Normal file
View File

@ -0,0 +1,13 @@
DAMASK - The Düsseldorf Advanced Material Simulation Kit
Visit damask.mpie.de for installation and usage instructions
CONTACT INFORMATION
Max-Planck-Institut für Eisenforschung GmbH
Max-Planck-Str. 1
40237 Düsseldorf
Germany
Email: DAMASK@mpie.de
https://damask.mpie.de
https://magit1.mpie.de

View File

@ -1 +0,0 @@
../../README

View File

@ -1,20 +1,20 @@
"""Main aggregator.""" """Main aggregator."""
import os import os
import re
with open(os.path.join(os.path.dirname(__file__),'VERSION')) as f:
version = f.readline()[1:-1]
name = 'damask' name = 'damask'
with open(os.path.join(os.path.dirname(__file__),'VERSION')) as f:
version = re.sub(r'^v','',f.readline().strip())
# classes # classes
from .environment import Environment # noqa from .environment import Environment # noqa
from .asciitable import ASCIItable # noqa
from .table import Table # noqa from .table import Table # noqa
from .asciitable import ASCIItable # noqa
from .config import Material # noqa from .config import Material # noqa
from .colormaps import Colormap, Color # noqa from .colormaps import Colormap, Color # noqa
from .orientation import Symmetry, Lattice, Rotation, Orientation # noqa from .orientation import Symmetry, Lattice, Rotation, Orientation # noqa
from .dadf5 import DADF5 # noqa from .dadf5 import DADF5 # noqa
from .geom import Geom # noqa from .geom import Geom # noqa
from .solver import Solver # noqa from .solver import Solver # noqa
@ -22,6 +22,10 @@ from .test import Test # noqa
from .util import extendableOption # noqa from .util import extendableOption # noqa
# functions in modules # functions in modules
from . import mechanics # noqa from . import mechanics # noqa
from . import grid_filters # noqa from . import grid_filters # noqa
# clean temporary variables
del os
del re
del f

View File

@ -1,27 +1,25 @@
import os import os
import re
class Environment(): class Environment():
__slots__ = [ \ __slots__ = [ \
'options', 'options',
] ]
def __init__(self): def __init__(self):
self.options = {} """Read and provide values of DAMASK configuration."""
self.get_options() self.options = {}
self.get_options()
def relPath(self,relative = '.'): def relPath(self,relative = '.'):
return os.path.join(self.rootDir(),relative) return os.path.join(self.rootDir(),relative)
def rootDir(self): def rootDir(self):
return os.path.normpath(os.path.join(os.path.realpath(__file__),'../../../')) return os.path.normpath(os.path.join(os.path.realpath(__file__),'../../../'))
def get_options(self): def get_options(self):
with open(self.relPath(self.rootDir()+'/CONFIG')) as configFile: for item in ['DAMASK_NUM_THREADS',
for line in configFile: 'MSC_ROOT',
l = re.sub('^set ', '', line).strip() # remove "set" (tcsh) when setting variables 'MARC_VERSION',
if l and not l.startswith('#'): 'ABAQUS_VERSION',
items = re.split(r'\s*=\s*',l) ]:
if len(items) == 2: self.options[item] = os.environ[item] if item in os.environ else None
self.options[items[0].upper()] = \
re.sub('\$\{*DAMASK_ROOT\}*',self.rootDir(),os.path.expandvars(items[1])) # expand all shell variables and DAMASK_ROOT

View File

@ -1,743 +0,0 @@
# -*- coding: UTF-8 no BOM -*-
# This tool converts a msc.marc result file into the vtk format that
# can be viewed by Paraview software (Kitware), or MayaVi (needs xml-vtk, or ...
#
# About the vtk format: http://www.vtk.org/VTK/project/about.html
# Some example vtk files: http://people.sc.fsu.edu/~jburkardt/data/vtk/vtk.html
# www.paraview.org
import os,sys,re
import numpy as np
import py_post # MSC closed source module to access marc result files
class MARC_POST():
def __init__(self):
self.projdir='./'
def opent16(self,incr=None):
self.fpath=os.path.join(self.projdir,self.postname)
print('Trying to open ',self.fpath,' ...')
self.p=py_post.post_open(self.fpath)
if self.p is None:
print('Could not open %s.'%self.postname); #return 'err'#; sys.exit(1)
raise Exception('Could not open t16')
print('Postfile %s%s is open ...'%(self.projdir,self.postname))
self.maxincr=self.p.increments()
print('and has %i increments'%self.maxincr)
if incr is None:
self.p.moveto(self.maxincr-1)
else:
self.p.moveto(incr+1)
print('moved to increment ', self.p.increment)
self.p.extrapolation('translate') # linear, translate, average. query with p.extrapolate
print('extrapolation method is ', self.p.extrapolate)
print('collecting model information')
self.t16info(printFlag=0)
print('t16 is open')
self.p
def t16info(self, printFlag=1
):
if not self.p:
self.p=self.opent16()#post_open(self.postname)
print(self.p)
oldincr=self.p.position
if oldincr==0:
self.p.moveto(1)
self.nnodes=self.p.nodes() #self.p.node(self.nnodes) crashes; only to nnodes-1 possible
self.nodes=range(0,self.nnodes)
self.nscals=self.p.node_scalars(); #print 'nscals', nscals
self.nscal_list=['']*self.nscals
self.nel=self.p.elements()
self.elscals=self.p.element_scalars()
self.elscal_list=['']*self.elscals
self.eltens=self.p.element_tensors()
self.elten_list=['']*self.eltens
for i in range(0,self.nscals):
self.nscal_list[i]=self.p.node_scalar_label(i)
for i in range (0,self.elscals):
self.elscal_list[i]=self.p.element_scalar_label(i)
if printFlag==1: print(i, self.elscal_list[i])
for i in range (0,self.eltens):
self.elten_list[i]=self.p.element_tensor_label(i)
if printFlag==1: print(i, self.elten_list[i])
for i in range(0,self.p.element_tensors()):
if printFlag==1: print('Element Tensor: ', i, self.p.element_tensor_label(i))
if printFlag==1: print('')
for i in range(0,self.p.element_scalars()):
if printFlag==1: print('Element Scalar: ', i, self.p.element_scalar_label(i))
if oldincr==0:
self.p.moveto(0)
def closet16(self):
print('Trying to close FEM result file ...')
try:
if self.p:
self.p.close()
print('FEM result file closed.')
self.p=None
else:
print('post object not open?')
except:
print('ERROR. Could not close FEM result file.')
def getLabelNr(self, label=None, type='Scalar'):
if type[0]=='S' or type[0]=='s': # element scalar
labelNr=self.elscal_list.index(label)
elif type[0]=='N': # node scalar
labelNr=self.nscal_list.index(label)
elif type[0]=='T' or type[0]=='t': # tensor
labelNr=self.elten_list.index(label)
print('Found label %s at index %i'%(label,labelNr))
return labelNr
def writeNodes2VTK(self, fobj):
self.points=[]
self.VTKcnt=200 # number of values per line in vtk file
fobj.write('POINTS %i'%self.p.nodes()+' float\n')
self.nodes_dict={} # store the node IDs in case of holes in the numbering
for iNd in self.nodes:
nd=self.p.node(iNd)
disp=self.p.node_displacement(iNd)
nd_xyz=[nd.x+disp[0], nd.y+disp[1], nd.z+disp[2]]
self.points.append(nd_xyz) # for pyvtk
fobj.write('%f %f %f \n'%
(nd.x+disp[0], nd.y+disp[1], nd.z+disp[2]))
self.nodes_dict[nd.id-1]=iNd
fobj.write('\n')
print('Nodes written to VTK: %i'%self.p.nodes())
def writeElements2VTK(self, fobj):
fobj.write('\nCELLS %i %i'%(self.p.elements(),self.p.elements()*9)+'\n')
self.cells=[] #for pyvtk
for iEl in range(0,self.nel):
el=self.p.element(iEl)
cell_nodes=[] # for pyvtk
ndlist=el.items
for k in [0, 1, 2, 3, 4, 5, 6, 7]: # FOR CELL TYPE VTK_HEXAHEDRON
node=ndlist[k]-1
cell_nodes.append(self.nodes_dict[node])
self.cells.append(cell_nodes) # for pyvtk
for e in self.cells:
fobj.write('8 ')
for n in e:
fobj.write('%6i '%n)
fobj.write('\n')
fobj.write('\nCELL_TYPES %i'%self.p.elements()+'\n')
cnt=0
for iEl in range(0,self.nel):
cnt=cnt+1
#fobj.write('11\n') #VTK_VOXEL
fobj.write('12 ') #VTK_HEXAHEDRON
if cnt>self.VTKcnt:
fobj.write('\n');cnt=0
fobj.write('\n')
print('Elements written to VTK: %i'%self.p.elements())
def writeElScalars2NodesVTK(self,fobj):
fobj.write('\nPOINT_DATA %i\n'%self.p.nodes())
nScal=12
nComponents=1+nScal
fobj.write('SCALARS scalars float %i\n'%nComponents)
fobj.write('LOOKUP_TABLE default\n')
idxScal=self.nscal_list.index('Displacement Z')
for iNd in self.nodes:
fobj.write('%f '%(self.p.node_scalar(iNd,idxScal)))
for iEl in range(0,self.nel):
el=self.p.element(iEl)
ndlist=el.items
if (iNd+1) in ndlist:
idx=ndlist.index(iNd+1)
for iV in range(0,nScal):
elData=self.p.element_scalar(iEl,35+iV)
fobj.write('%f '%(elData[idx].value))
break
fobj.write('\n')
fobj.write('\n')
def writeNodeScalars2VTK(self,fobj):
fobj.write('\nPOINT_DATA %i\n'%self.p.nodes())
self.pointDataScalars=[]
for idxNdScal in range(-3,self.nscals): #now include node x,y,z
if idxNdScal>=0:
datalabel=self.nscal_list[idxNdScal]
datalabel=re.sub("\s",'_',datalabel)
else:
if idxNdScal==-3: datalabel='node.x'
if idxNdScal==-2: datalabel='node.y'
if idxNdScal==-1: datalabel='node.z'
fobj.write('SCALARS %s float %i\n'%(datalabel,1))#nComponents))
fobj.write('LOOKUP_TABLE default\n')
cnt=0
for iNd in range(0,self.nnodes):
cnt=cnt+1
if idxNdScal>=0:
ndData=self.p.node_scalar(iNd,idxNdScal)
else:
nd=self.p.node(iNd)
if idxNdScal==-3: ndData=nd.x
if idxNdScal==-2: ndData=nd.y
if idxNdScal==-1: ndData=nd.z
fobj.write('%E '%(ndData))
if cnt>self.VTKcnt:
fobj.write('\n')
cnt=0
fobj.write('\n')
fobj.write('\n')
def writeElementData2VTK(self,fobj):
self.sig_vMises=[]
self.sig33=[]
idx_sig_vMises=self.getLabelNr('Equivalent Von Mises Stress')
idx_sig33=self.getLabelNr('Comp 33 of Cauchy Stress')
fobj.write('\nCELL_DATA %i\n'%self.p.elements())
for idxElScal in range(0,self.elscals):
datalabel=self.elscal_list[idxElScal]
datalabel=re.sub("\s",'_',datalabel)
fobj.write('\n\nSCALARS %s float %i\n'%(datalabel,1))
fobj.write('LOOKUP_TABLE default\n')
cnt=0
for iEl in range(0,self.nel):
cnt=cnt+1
elData=self.p.element_scalar(iEl,idxElScal)
avgScal=0.0
if datalabel in ['phi1', 'PHI','phi2']: # Euler angles should not be averaged
avgScal=avgScal+elData[0].value
else:
for IP in range(0,8):
avgScal=avgScal+elData[IP].value
avgScal=avgScal/8.
fobj.write('%E '%(avgScal))
if idxElScal==idx_sig_vMises:
self.sig_vMises.append(avgScal)
elif idxElScal==idx_sig33:
self.sig33.append(avgScal)
if cnt>self.VTKcnt:
fobj.write('\n')
cnt=0
fobj.write('\n')
def get_avg_el_scal(self,idxElScal):
result=[]
datalabel=self.elscal_list[idxElScal]
print('Collecting %s from all elements'%datalabel)
for iEl in range(0,self.nel):
elData=self.p.element_scalar(iEl,idxElScal)
avgScal=0.0
for IP in range(0,8):
avgScal=avgScal+elData[IP].value
avgScal=avgScal/8.
result.append(avgScal)
return result
def writeUniaxiality2VTK(self,fobj):
datalabel='uniaxiality_sig_vMises_durch_sig33'
fobj.write('SCALARS %s float %i\n'%(datalabel,1))
fobj.write('LOOKUP_TABLE default\n')
cnt=0
for iEl in range(0,self.nel):
cnt=cnt+1
if abs(self.sig_vMises[iEl])<1e-5:
datum=0.
else:
datum=self.sig33[iEl]/self.sig_vMises[iEl]
fobj.write('%E '%(datum))
if cnt>self.VTKcnt:
fobj.write('\n')
cnt=0
fobj.write('\n')
def stress_per_element(self):
self.stress=[]
for iEl in range(0,self.nel):
sig=self.avg_elten(2,elID=iEl)
self.stress.append(sig[0])
def mean_stress_per_element(self):
self.mean_stress=[]
for iEl in range(0,self.nel):
sig=self.stress[iEl]
self.mean_stress.append(self.meanStress(sig))
def triaxiality_per_element(self):
# classical triaxiality
# 1/3 : uniax tension
self.triaxiality=[]
for iEl in range(0,self.nel):
t=self.mean_stress[iEl]/self.sig_vMises[iEl]
self.triaxiality.append(t)
def moreElData2VTK(self,fobj,data=[],label='datalabel'):
fobj.write('SCALARS %s float %i\n'%(label,1))
fobj.write('LOOKUP_TABLE default\n')
cnt=0
for iEl in range(0,self.nel):
cnt=cnt+1
fobj.write('%E '%(data[iEl]))
if cnt>self.VTKcnt:
fobj.write('\n')
cnt=0
fobj.write('\n')
def calc_lode_parameter(self):
self.lode=[]
try:
self.stress
except:
self.stress_per_element()
for iEl in range(0,self.nel):
sig=self.stress[iEl]
lode=self.stress2lode(sig)
self.lode.append(lode)
def stress2lode(self,stress):
[pStress,pAxes]=self.princStress(stress)
s1=pStress[0]
s2=pStress[1]
s3=pStress[2]
lode=(2*s2-s1-s3) / ( s1 - s3 )
return lode
def princStress(self, stress):
"""
Function to compute 3D principal stresses and sort them.
from: http://geodynamics.org/svn/cig/short/3D/PyLith/trunk/playpen/postproc/vtkcff.py
"""
stressMat=np.array(stress)
(princStress, princAxes) = np.linalg.eigh(stressMat)
idx = princStress.argsort()
princStressOrdered = princStress[idx]
princAxesOrdered = princAxes[:,idx]
return princStressOrdered, princAxesOrdered
def avg_elten(self,
idxElTen, mat=0, elID=None):
tensum=np.zeros((3,3));
T=np.zeros((3,3));
pts=0;
avg=np.zeros((3,3));
if elID is None:
averaged_elements=range(0,self.nel)
else:
averaged_elements=[elID]
for i in averaged_elements:
if mat==0 or int(self.p.element_scalar(i,4)[0].value)==mat:
T=self.p.element_tensor(i,idxElTen)
for k in range (0,8):
tensum[0][0] = tensum[0][0] + T[k].t11
tensum[0][1] = tensum[0][1] + T[k].t12
tensum[0][2] = tensum[0][2] + T[k].t13
tensum[1][1] = tensum[1][1] + T[k].t22
tensum[1][2] = tensum[1][2] + T[k].t23
tensum[2][2] = tensum[2][2] + T[k].t33
pts=pts+1
avg=tensum/pts
avg=self.fillComponents(avg)
del [T]
return (avg,tensum,pts)
def fillComponents(self,
halftensor
):
halftensor[1][0]=halftensor[0][1]
halftensor[2][0]=halftensor[0][2]
halftensor[2][1]=halftensor[1][2]
return halftensor
def vMises(self,tensor33):
t=tensor33
s=(t[0,0]-t[1,1])**2+(t[1,1]-t[2,2])**2+(t[0,0]-t[2,2])**2+\
6*(t[0,1]**2+t[1,2]**2+t[2,0]**2)
vM=np.sqrt(s/2.)
return vM
def meanStress(self,tensor33):
t=tensor33
s=t[0,0]+t[1,1]+t[2,2]
ms=s/3.
return ms
def invariants(self,tensor33):
t=tensor33
I1=t[0,0]+t[1,1]+t[2,2]
I2=t[0,0]*t[1,1]+t[1,1]*t[2,2]+t[0,0]*t[2,2]-\
t[0,1]**2-t[1,2]**2-t[0,2]**2
I3=t[0,0]*t[1,1]*t[2,2]+\
2*t[0,1]*t[1,2]*t[2,0]-\
t[2,2]*t[0,1]**2-t[0,0]*t[1,2]**2-t[1,1]*t[0,2]**2
return [ I1, I2, I3 ]
class VTK_WRITER():
"""
The resulting vtk-file can be imported in Paraview 3.12
Then use Filters: Cell Data to Point Data + Contour
to plot semi-transparent iso-surfaces.
"""
def __init__(self):
self.p=MARC_POST() # self.p
def openFile(self, filename='test.vtp'):
self.f=open(filename,'w+')
self.fname=filename
def writeFirstLines(self,
vtkFile=None,
version='2.0',
comment='Test',
dformat='ASCII', # BINARY | [ASCII]
dtype='UNSTRUCTURED_GRID' # UNSTRUCTURED GRID
):
if vtkFile is None:
vtkFile=self.f
# First Line contains Data format version
self.versionVTK=version
vtkFile.write('# vtk DataFile Version %s\n'%self.versionVTK)
# Comment goes to 2nd line and has maximum 256 chars
vtkFile.write(comment+'\n')
vtkFile.write(dformat+'\n')
vtkFile.write('DATASET '+dtype+'\n')
def marc2vtkBatch(self):
for iori in range(1,63):
self.p.postname='indent_fric0.3_R2.70_cA146.0_h0.320_ori%03i_OST_h19d.t16'%(iori)
if os.path.exists(self.p.postname):
self.marc2vtk(mode='fast', batchMode=1)
def marc2vtk(self, label=None, mode='fast',
batchMode=0,
incRange=None,
incStepMult=1.):
if batchMode==0:
try:
self.p
except:
self.p=MARC_POST()
### ---- CHANGE dir/job/model to process here
os.chdir('M:/nicu');
jobname='ori001'
self.p.postname='indent_fric0.3_R0.25_cA90.0_h0.010_4320els_'+jobname+'.t16'
### ----
self.p.opent16()
self.p.t16info()
incMax=self.p.p.increments();
if incRange is None:
incStep=5
incRange=range(0,incMax+1,incStep)
self.vtkPath=os.getcwd()+'/vtk_%s/'%self.p.postname
if not os.path.exists(self.vtkPath):
os.mkdir(self.vtkPath)
for inc in incRange:
print('Increment: %i'%inc)
self.p.p.moveto(inc)
t=self.p.p.time
sys.stdout.write('inc:%i, time:%.3f\n'%(self.p.p.increment,t))
self.incStr='inc%03i'%(inc*incStepMult)
self.openFile(filename=self.vtkPath+self.p.postname[0:-4]+'_'+
self.incStr+'.vtk')
self.writeFirstLines(comment=self.p.postname,
dtype='UNSTRUCTURED_GRID')
self.p.writeNodes2VTK(fobj=self.f)
self.p.writeElements2VTK(fobj=self.f)
self.p.writeNodeScalars2VTK(fobj=self.f)
self.p.writeElementData2VTK(fobj=self.f)
# insert generation and writing of derived post values
# *here*
self.f.close()
print('Increment (self.p.p.increment): %i'%self.p.p.increment)
print('Data written.')
print(self.p.postname+' ready.')
def scaleBar(self, length=1.0, posXYZ=[0., 0., 0.]):
self.fsb=open('micronbar_l%.1f.vtp'%length,'w+')
self.writeFirstLines(self.fsb, comment='micronbar')
pts=np.array([])
width=length*1.
height=length*1.
wVec=np.array([0., width, 0.])
lVec=np.array([length,0.,0.])
hVec=np.array([0.,0.,height])
posXYZ=posXYZ-0.5*wVec-0.5*lVec#-0.5*hVec # CENTERING Y/N
posXYZ=np.array(posXYZ)
pts=[posXYZ, posXYZ+lVec,
posXYZ+wVec,
posXYZ+wVec+lVec]
pts.extend([pts[0]+hVec,pts[1]+hVec,pts[2]+hVec,pts[3]+hVec])
print(len(pts), pts)
self.fsb.write('POINTS %i float\n'%len(pts))
for npts in range(0,len(pts)):
self.fsb.write('%f %f %f\n'%(pts[npts][0], pts[npts][1], pts[npts][2]))
if 1: #Triad
nCells=3
ptsPerCell=2 # Lines (Type=3)
cellSize=(ptsPerCell+1)*nCells
self.fsb.write('CELLS %i %i\n'%(nCells,cellSize))
self.fsb.write('2 0 1\n') #X-Line
self.fsb.write('2 0 2\n') #Y-Line
self.fsb.write('2 0 4\n') #Z-Line
self.fsb.write('CELL_TYPES %i\n'%(nCells))
self.fsb.write('3\n3\n3\n')#Line
else: # Cube, change posXYZ
nCells=1
ptsPerCell=2 # Lines (Type=3)
cellSize=(ptsPerCell+1)*nCells
self.fsb.write('CELLS %i %i\n'%(nCells,cellSize))
self.fsb.write('2 0 1\n') #Line
self.fsb.write('CELL_TYPES %i\n'%(nCells))
self.fsb.write('3\n')#Line
self.fsb.write('\n')
self.fsb.close()
print(self.fsb)
def example_unstructured(self):
self.openFile(filename='example_unstructured_grid.vtk')
self.f.write("""
# vtk DataFile Version 2.0
example_unstruct_grid
ASCII
POINTS 12 float
0 0 0
1 0 0
1 1 0
0 1 0
0 0 1
1 0 1
1 1 1
0 1 1
0 0 1.9
1 0 1.9
1 1 1.9
0 1 1.9
CELLS 2 18
8 0 1 2 3 4 5 6 7
8 4 5 6
7 8 9 10 11
CELL_TYPES 2
12
12
POINT_DATA 12
SCALARS nodex float 1
LOOKUP_TABLE default
2.34E+12
2.00
0.00
1.62
5.03
1.02
1.50
0.00
3 5 6 23423423423423423423.23423423""")
self.f.close()
def writeNodes2VTK(self, fobj):
self.VTKcnt=200 # how many numbers per line in vtk file
fobj.write('POINTS %i'%self.p.nodes()+' float\n')
for iNd in self.nodes:
nd=self.p.node(iNd)
disp=self.p.node_displacement(iNd)
fobj.write('%f %f %f \n'%
(nd.x+disp[0], nd.y+disp[1], nd.z+disp[2]))
fobj.write('\n')
print('Nodes written to VTK: %i'%self.p.nodes())
def writeElements2VTK(self, fobj):
fobj.write('\nCELLS %i %i'%(self.p.elements(),self.p.elements()*9)+'\n')
for iEl in range(0,self.nel):
el=self.p.element(iEl)
fobj.write('8 ')
ndlist=el.items
for k in [0, 1, 2, 3, 4, 5, 6, 7]: # FOR CELL TYPE VTK_HEXAHEDRON
fobj.write('%6i '%(ndlist[k]-1))
fobj.write('\n')
fobj.write('\nCELL_TYPES %i'%self.p.elements()+'\n')
cnt=0
for iEl in range(0,self.nel):
cnt=cnt+1
fobj.write('12 ') #VTK_HEXAHEDRON
if cnt>self.VTKcnt:
fobj.write('\n');cnt=0
fobj.write('\n')
print('Elements written to VTK: %i'%self.p.elements())
def writeElScalars2NodesVTK(self,fobj):
fobj.write('\nPOINT_DATA %i\n'%self.p.nodes())
nScal=12
nComponents=1+nScal
fobj.write('SCALARS scalars float %i\n'%nComponents)
fobj.write('LOOKUP_TABLE default\n')
idxScal=self.nscal_list.index('Displacement Z')
for iNd in self.nodes:
fobj.write('%f '%(self.p.node_scalar(iNd,idxScal)))
for iEl in range(0,self.nel):
el=self.p.element(iEl)
ndlist=el.items
if (iNd+1) in ndlist:
idx=ndlist.index(iNd+1)
for iV in range(0,nScal):
elData=self.p.element_scalar(iEl,35+iV)
fobj.write('%f '%(elData[idx].value))
break
fobj.write('\n')
fobj.write('\n')
def writeNodeScalars2VTK(self,fobj):
fobj.write('\nPOINT_DATA %i\n'%self.p.nodes())
for idxNdScal in range(-3,self.nscals): # include node x,y,z
if idxNdScal>=0:
datalabel=self.nscal_list[idxNdScal]
datalabel=re.sub("\s",'_',datalabel)
else:
if idxNdScal==-3: datalabel='node.x'
if idxNdScal==-2: datalabel='node.y'
if idxNdScal==-1: datalabel='node.z'
fobj.write('SCALARS %s float %i\n'%(datalabel,1))#nComponents))
fobj.write('LOOKUP_TABLE default\n')
cnt=0
for iNd in range(0,self.nnodes):
cnt=cnt+1
if idxNdScal>=0:
ndData=self.p.node_scalar(iNd,idxNdScal)
else:
nd=self.p.node(iNd)
if idxNdScal==-3: ndData=nd.x
if idxNdScal==-2: ndData=nd.y
if idxNdScal==-1: ndData=nd.z
fobj.write('%E '%(ndData))
if cnt>self.VTKcnt:
fobj.write('\n')
cnt=0
fobj.write('\n')
fobj.write('\n')
def writeElementData2VTK(self,fobj):
fobj.write('\nCELL_DATA %i\n'%self.p.elements())
for idxElScal in range(0,self.elscals):
datalabel=self.elscal_list[idxElScal]
datalabel=re.sub("\s",'_',datalabel)
fobj.write('\n\nSCALARS %s float %i\n'%(datalabel,1))#nComponents))
fobj.write('LOOKUP_TABLE default\n')
cnt=0
for iEl in range(0,self.nel):
cnt=cnt+1
elData=self.p.element_scalar(iEl,idxElScal)
avgScal=0.0
if datalabel in ['phi1', 'PHI','phi2']:
avgScal=avgScal+elData[0].value
else:
for IP in range(0,8):
avgScal=avgScal+elData[IP].value
avgScal=avgScal/8.
fobj.write('%E '%(avgScal))
if cnt>self.VTKcnt:
fobj.write('\n')
cnt=0
fobj.write('\n')
def example1(self):
self.openFile()
self.writeFirstLines()
self.f.write("""DATASET POLYDATA
POINTS 8 float
0.0 0.0 0.0
1.0 0.0 0.0
1.0 1.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0
1.0 0.0 1.0
1.0 1.0 1.0
0.0 1.0 1.0
POLYGONS 6 30
4 0 1 2 3
4 4 5 6 7
4 0 1 5 4
4 2 3 7 6
4 0 4 7 3
4 1 2 6 5
CELL_DATA 6
SCALARS cell_scalars int 1
LOOKUP_TABLE default
0
1
2
3
4
5
NORMALS cell_normals float
0 0 -1
0 0 1
0 -1 0
0 1 0
-1 0 0
1 0 0
FIELD FieldData 2
cellIds 1 6 int
0 1 2 3 4 5
faceAttributes 2 6 float
0.0 1.0 1.0 2.0 2.0 3.0 3.0 4.0 4.0 5.0 5.0 6.0
POINT_DATA 8
SCALARS sample_scalars float 1
LOOKUP_TABLE my_table
0.0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
LOOKUP_TABLE my_table 8
0.0 0.0 0.0 1.0
1.0 0.0 0.0 1.0
0.0 1.0 0.0 1.0
1.0 1.0 0.0 1.0
0.0 0.0 1.0 1.0
1.0 0.0 1.0 1.0
0.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0""")
self.f.close()
import pyvtk
class marc_to_vtk():
"""
Anybody wants to implement it with pyvtk?
The advantage would be that pyvtk can also wirte the
<xml>-VTK format and binary.
These can be plotted with mayavi.
"""
def __init__(self):
self.p=[]#MARC_POST() # self.p
def run(self):
vtk = pyvtk.VtkData(\
pyvtk.UnstructuredGrid(self.p.points,
hexahedron=self.p.cells),
'm2v output')
vtk.tofile('m2v_file')

View File

@ -5,30 +5,32 @@ import damask
class Abaqus(Solver): class Abaqus(Solver):
"""Wrapper to run DAMASK with Abaqus.""" """Wrapper to run DAMASK with Abaqus."""
def __init__(self,version=int(damask.Environment().options['ABAQUS_VERSION'])): def __init__(self,version=damask.Environment().options['ABAQUS_VERSION']):
""" """
Create a Abaqus solver object. Create a Abaqus solver object.
Parameters Parameters
---------- ----------
version : integer version : integer
Abaqus version Abaqus version
""" """
self.solver ='Abaqus' self.solver = 'Abaqus'
self.version = int(version) try:
self.version = int(version)
except TypeError:
self.version = -1
def return_run_command(self,model): def return_run_command(self,model):
env=damask.Environment() try:
try: cmd = 'abq{}'.format(self.version)
cmd='abq'+self.version subprocess.check_output([cmd,'information=release'])
subprocess.check_output([cmd,'information=release']) except OSError: # link to abqXXX not existing
except OSError: # link to abqXXX not existing cmd = 'abaqus'
cmd='abaqus' process = subprocess.Popen([cmd,'information=release'],stdout = subprocess.PIPE,stderr = subprocess.PIPE)
process = subprocess.Popen(['abaqus','information=release'],stdout = subprocess.PIPE,stderr = subprocess.PIPE) detectedVersion = int(process.stdout.readlines()[1].split()[1].decode('utf-8'))
detectedVersion = int(process.stdout.readlines()[1].split()[1].decode('utf-8')) if self.version not in [detectedVersion,-1]:
if self.version != detectedVersion: raise Exception('found Abaqus version {}, but requested {}'.format(detectedVersion,self.version))
raise Exception('found Abaqus version {}, but requested {}'.format(detectedVersion,self.version)) return '{} -job {} -user {}/src/DAMASK_abaqus interactive'.format(cmd,model,damask.Environment().rootDir())
return '{} -job {} -user {}/src/DAMASK_abaqus interactive'.format(cmd,model,env.rootDir())

View File

@ -1,89 +1,91 @@
import os import os
import subprocess import subprocess
import shlex import shlex
import string
from .solver import Solver from .solver import Solver
import damask import damask
class Marc(Solver): class Marc(Solver):
"""Wrapper to run DAMASK with MSCMarc.""" """Wrapper to run DAMASK with MSCMarc."""
def __init__(self,version=damask.Environment().options['MARC_VERSION']): def __init__(self,version=damask.Environment().options['MARC_VERSION']):
""" """
Create a Marc solver object. Create a Marc solver object.
Parameters Parameters
---------- ----------
version : float version : float
Marc version Marc version
""" """
self.solver ='Marc' self.solver = 'Marc'
self.version = damask.environment.Environment().options['MARC_VERSION'] try:
self.version = int(version)
except TypeError:
self.version = -1
#-------------------------- #--------------------------
def libraryPath(self): def libraryPath(self):
path_MSC = damask.environment.Environment().options['MSC_ROOT'] path_MSC = damask.Environment().options['MSC_ROOT']
path_lib = '{}/mentat{}/shlib/linux64'.format(path_MSC,self.version) path_lib = '{}/mentat{}/shlib/linux64'.format(path_MSC,self.version)
return path_lib if os.path.exists(path_lib) else '' return path_lib if os.path.exists(path_lib) else ''
#-------------------------- #--------------------------
def toolsPath(self): def toolsPath(self):
path_MSC = damask.environment.Environment().options['MSC_ROOT'] path_MSC = damask.Environment().options['MSC_ROOT']
path_tools = '{}/marc{}/tools'.format(path_MSC,self.version) path_tools = '{}/marc{}/tools'.format(path_MSC,self.version)
return path_tools if os.path.exists(path_tools) else '' return path_tools if os.path.exists(path_tools) else ''
#-------------------------- #--------------------------
def submit_job(self, def submit_job(self,
model, model,
job = 'job1', job = 'job1',
logfile = False, logfile = False,
compile = False, compile = False,
optimization ='', optimization = '',
): ):
damaskEnv = damask.environment.Environment() damaskEnv = damask.Environment()
user = os.path.join(damaskEnv.relPath('src'),'DAMASK_marc{}.{}'.format(self.version,'f90' if compile else 'marc')) user = os.path.join(damaskEnv.relPath('src'),'DAMASK_marc{}.{}'.format(self.version,'f90' if compile else 'marc'))
if not os.path.isfile(user): if not os.path.isfile(user):
raise FileNotFoundError("DAMASK4Marc ({}) '{}' not found".format(('source' if compile else 'binary'),user)) raise FileNotFoundError("DAMASK4Marc ({}) '{}' not found".format(('source' if compile else 'binary'),user))
# Define options [see Marc Installation and Operation Guide, pp 23] # Define options [see Marc Installation and Operation Guide, pp 23]
script = 'run_damask_{}mp'.format(optimization) script = 'run_damask_{}mp'.format(optimization)
cmd = os.path.join(self.toolsPath(),script) + \ cmd = os.path.join(self.toolsPath(),script) + \
' -jid ' + model + '_' + job + \ ' -jid ' + model + '_' + job + \
' -nprocd 1 -autorst 0 -ci n -cr n -dcoup 0 -b no -v no' ' -nprocd 1 -autorst 0 -ci n -cr n -dcoup 0 -b no -v no'
if compile: cmd += ' -u ' + user + ' -save y' if compile: cmd += ' -u ' + user + ' -save y'
else: cmd += ' -prog ' + os.path.splitext(user)[0] else: cmd += ' -prog ' + os.path.splitext(user)[0]
print('job submission with{} compilation: {}'.format('' if compile else 'out',user)) print('job submission {} compilation: {}'.format('with' if compile else 'without',user))
if logfile: log = open(logfile, 'w') if logfile: log = open(logfile, 'w')
print(cmd) print(cmd)
process = subprocess.Popen(shlex.split(cmd),stdout = log,stderr = subprocess.STDOUT) process = subprocess.Popen(shlex.split(cmd),stdout = log,stderr = subprocess.STDOUT)
log.close() log.close()
process.wait() process.wait()
#-------------------------- #--------------------------
def exit_number_from_outFile(self,outFile=None): def exit_number_from_outFile(self,outFile=None):
import string exitnumber = -1
exitnumber = -1 with open(outFile,'r') as fid_out:
fid_out = open(outFile,'r') for line in fid_out:
for line in fid_out: if (string.find(line,'tress iteration') != -1):
if (string.find(line,'tress iteration') != -1): print(line)
print(line) elif (string.find(line,'Exit number') != -1):
elif (string.find(line,'Exit number') != -1): substr = line[string.find(line,'Exit number'):len(line)]
substr = line[string.find(line,'Exit number'):len(line)] exitnumber = int(substr[12:16])
exitnumber = int(substr[12:16])
fid_out.close() return exitnumber
return exitnumber

View File

@ -1,8 +1,9 @@
import setuptools import setuptools
import os import os
import re
with open(os.path.join(os.path.dirname(__file__),'damask/VERSION')) as f: with open(os.path.join(os.path.dirname(__file__),'damask/VERSION')) as f:
version = f.readline()[1:-1] version = re.sub(r'(-([^-]*)).*$',r'.\2',re.sub(r'^v(\d+\.\d+(\.\d+)?)',r'\1',f.readline().strip()))
setuptools.setup( setuptools.setup(
name="damask", name="damask",
@ -18,12 +19,13 @@ setuptools.setup(
"pandas", "pandas",
"scipy", "scipy",
"h5py", "h5py",
"vtk" "vtk",
], ],
license = 'GPL3',
classifiers = [ classifiers = [
"Intended Audience :: Science/Research",
"Topic :: Scientific/Engineering",
"Programming Language :: Python :: 3", "Programming Language :: Python :: 3",
"License :: OSI Approved :: GPL3", "License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)",
"Operating System :: OS Independent", "Operating System :: OS Independent",
], ],
) )

View File

@ -2,8 +2,6 @@ import os
import pytest import pytest
import damask
def pytest_addoption(parser): def pytest_addoption(parser):
parser.addoption("--update", parser.addoption("--update",
action="store_true", action="store_true",
@ -17,5 +15,4 @@ def update(request):
@pytest.fixture @pytest.fixture
def reference_dir_base(): def reference_dir_base():
"""Directory containing reference results.""" """Directory containing reference results."""
env = damask.Environment() return os.path.join(os.path.dirname(__file__),'reference')
return os.path.join(env.rootDir(),'python','tests','reference')

View File

@ -433,10 +433,10 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
forall (i=1:3, j=1:3) do i =1, 3;do j=1,3
dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
end forall enddo; enddo
end subroutine constitutive_hooke_SandItsTangents end subroutine constitutive_hooke_SandItsTangents

View File

@ -43,10 +43,17 @@ module lattice
LATTICE_FCC_NCLEAVAGESYSTEM = [3, 4] !< # of cleavage systems per family for fcc LATTICE_FCC_NCLEAVAGESYSTEM = [3, 4] !< # of cleavage systems per family for fcc
integer, parameter :: & integer, parameter :: &
#ifndef __PGI
LATTICE_FCC_NSLIP = sum(LATTICE_FCC_NSLIPSYSTEM), & !< total # of slip systems for fcc LATTICE_FCC_NSLIP = sum(LATTICE_FCC_NSLIPSYSTEM), & !< total # of slip systems for fcc
LATTICE_FCC_NTWIN = sum(LATTICE_FCC_NTWINSYSTEM), & !< total # of twin systems for fcc LATTICE_FCC_NTWIN = sum(LATTICE_FCC_NTWINSYSTEM), & !< total # of twin systems for fcc
LATTICE_FCC_NTRANS = sum(LATTICE_FCC_NTRANSSYSTEM), & !< total # of transformation systems for fcc LATTICE_FCC_NTRANS = sum(LATTICE_FCC_NTRANSSYSTEM), & !< total # of transformation systems for fcc
LATTICE_FCC_NCLEAVAGE = sum(LATTICE_FCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for fcc LATTICE_FCC_NCLEAVAGE = sum(LATTICE_FCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for fcc
#else
LATTICE_FCC_NSLIP = 18, &
LATTICE_FCC_NTWIN = 12, &
LATTICE_FCC_NTRANS = 12, &
LATTICE_FCC_NCLEAVAGE = 7
#endif
real(pReal), dimension(3+3,LATTICE_FCC_NSLIP), parameter :: & real(pReal), dimension(3+3,LATTICE_FCC_NSLIP), parameter :: &
LATTICE_FCC_SYSTEMSLIP = reshape(real([& LATTICE_FCC_SYSTEMSLIP = reshape(real([&
@ -128,9 +135,15 @@ module lattice
LATTICE_BCC_NCLEAVAGESYSTEM = [3, 6] !< # of cleavage systems per family for bcc LATTICE_BCC_NCLEAVAGESYSTEM = [3, 6] !< # of cleavage systems per family for bcc
integer, parameter :: & integer, parameter :: &
#ifndef __PGI
LATTICE_BCC_NSLIP = sum(LATTICE_BCC_NSLIPSYSTEM), & !< total # of slip systems for bcc LATTICE_BCC_NSLIP = sum(LATTICE_BCC_NSLIPSYSTEM), & !< total # of slip systems for bcc
LATTICE_BCC_NTWIN = sum(LATTICE_BCC_NTWINSYSTEM), & !< total # of twin systems for bcc LATTICE_BCC_NTWIN = sum(LATTICE_BCC_NTWINSYSTEM), & !< total # of twin systems for bcc
LATTICE_BCC_NCLEAVAGE = sum(LATTICE_BCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for bcc LATTICE_BCC_NCLEAVAGE = sum(LATTICE_BCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for bcc
#else
LATTICE_BCC_NSLIP = 24, &
LATTICE_BCC_NTWIN = 12, &
LATTICE_BCC_NCLEAVAGE = 9
#endif
real(pReal), dimension(3+3,LATTICE_BCC_NSLIP), parameter :: & real(pReal), dimension(3+3,LATTICE_BCC_NSLIP), parameter :: &
LATTICE_BCC_SYSTEMSLIP = reshape(real([& LATTICE_BCC_SYSTEMSLIP = reshape(real([&
@ -206,9 +219,15 @@ module lattice
LATTICE_HEX_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for hex LATTICE_HEX_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for hex
integer, parameter :: & integer, parameter :: &
#ifndef __PGI
LATTICE_HEX_NSLIP = sum(LATTICE_HEX_NSLIPSYSTEM), & !< total # of slip systems for hex LATTICE_HEX_NSLIP = sum(LATTICE_HEX_NSLIPSYSTEM), & !< total # of slip systems for hex
LATTICE_HEX_NTWIN = sum(LATTICE_HEX_NTWINSYSTEM), & !< total # of twin systems for hex LATTICE_HEX_NTWIN = sum(LATTICE_HEX_NTWINSYSTEM), & !< total # of twin systems for hex
LATTICE_HEX_NCLEAVAGE = sum(LATTICE_HEX_NCLEAVAGESYSTEM) !< total # of cleavage systems for hex LATTICE_HEX_NCLEAVAGE = sum(LATTICE_HEX_NCLEAVAGESYSTEM) !< total # of cleavage systems for hex
#else
LATTICE_HEX_NSLIP = 33, &
LATTICE_HEX_NTWIN = 24, &
LATTICE_HEX_NCLEAVAGE = 3
#endif
real(pReal), dimension(4+4,LATTICE_HEX_NSLIP), parameter :: & real(pReal), dimension(4+4,LATTICE_HEX_NSLIP), parameter :: &
LATTICE_HEX_SYSTEMSLIP = reshape(real([& LATTICE_HEX_SYSTEMSLIP = reshape(real([&
@ -301,7 +320,11 @@ module lattice
LATTICE_BCT_NSLIPSYSTEM = [2, 2, 2, 4, 2, 4, 2, 2, 4, 8, 4, 8, 8 ] !< # of slip systems per family for bct (Sn) Bieler J. Electr Mater 2009 LATTICE_BCT_NSLIPSYSTEM = [2, 2, 2, 4, 2, 4, 2, 2, 4, 8, 4, 8, 8 ] !< # of slip systems per family for bct (Sn) Bieler J. Electr Mater 2009
integer, parameter :: & integer, parameter :: &
#ifndef __PGI
LATTICE_BCT_NSLIP = sum(LATTICE_BCT_NSLIPSYSTEM) !< total # of slip systems for bct LATTICE_BCT_NSLIP = sum(LATTICE_BCT_NSLIPSYSTEM) !< total # of slip systems for bct
#else
LATTICE_BCT_NSLIP = 52
#endif
real(pReal), dimension(3+3,LATTICE_BCT_NSLIP), parameter :: & real(pReal), dimension(3+3,LATTICE_BCT_NSLIP), parameter :: &
LATTICE_BCT_SYSTEMSLIP = reshape(real([& LATTICE_BCT_SYSTEMSLIP = reshape(real([&
@ -379,7 +402,11 @@ module lattice
LATTICE_ISO_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for iso LATTICE_ISO_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for iso
integer, parameter :: & integer, parameter :: &
#ifndef __PGI
LATTICE_ISO_NCLEAVAGE = sum(LATTICE_ISO_NCLEAVAGESYSTEM) !< total # of cleavage systems for iso LATTICE_ISO_NCLEAVAGE = sum(LATTICE_ISO_NCLEAVAGESYSTEM) !< total # of cleavage systems for iso
#else
LATTICE_ISO_NCLEAVAGE = 3
#endif
real(pReal), dimension(3+3,LATTICE_ISO_NCLEAVAGE), parameter :: & real(pReal), dimension(3+3,LATTICE_ISO_NCLEAVAGE), parameter :: &
LATTICE_ISO_SYSTEMCLEAVAGE= reshape(real([& LATTICE_ISO_SYSTEMCLEAVAGE= reshape(real([&
@ -396,7 +423,11 @@ module lattice
LATTICE_ORT_NCLEAVAGESYSTEM = [1, 1, 1] !< # of cleavage systems per family for ortho LATTICE_ORT_NCLEAVAGESYSTEM = [1, 1, 1] !< # of cleavage systems per family for ortho
integer, parameter :: & integer, parameter :: &
#ifndef __PGI
LATTICE_ORT_NCLEAVAGE = sum(LATTICE_ORT_NCLEAVAGESYSTEM) !< total # of cleavage systems for ortho LATTICE_ORT_NCLEAVAGE = sum(LATTICE_ORT_NCLEAVAGESYSTEM) !< total # of cleavage systems for ortho
#else
LATTICE_ORT_NCLEAVAGE = 3
#endif
real(pReal), dimension(3+3,LATTICE_ORT_NCLEAVAGE), parameter :: & real(pReal), dimension(3+3,LATTICE_ORT_NCLEAVAGE), parameter :: &
LATTICE_ORT_SYSTEMCLEAVAGE = reshape(real([& LATTICE_ORT_SYSTEMCLEAVAGE = reshape(real([&

View File

@ -610,7 +610,7 @@ function om2ax(om) result(ax)
else else
call dgeev('N','V',3,om_,3,Wr,Wi,devNull,3,VR,3,work,size(work,1),ierr) call dgeev('N','V',3,om_,3,Wr,Wi,devNull,3,VR,3,work,size(work,1),ierr)
if (ierr /= 0) call IO_error(0,ext_msg='Error in om2ax: DGEEV return not zero') if (ierr /= 0) call IO_error(0,ext_msg='Error in om2ax: DGEEV return not zero')
#if defined(__GFORTRAN__) && __GNUC__<9 || defined(__INTEL_COMPILER) && INTEL_COMPILER<1800 #if defined(__GFORTRAN__) && __GNUC__<9 || defined(__INTEL_COMPILER) && INTEL_COMPILER<1800 || defined(__PGI)
i = maxloc(merge(1,0,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1) i = maxloc(merge(1,0,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1)
#else #else
i = findloc(cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal),.true.,dim=1) !find eigenvalue (1,0) i = findloc(cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal),.true.,dim=1) !find eigenvalue (1,0)
@ -1266,38 +1266,49 @@ subroutine unitTest
sin(2.0_pReal*PI*x(1))*A] sin(2.0_pReal*PI*x(1))*A]
if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal) if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
endif endif
#ifndef __PGI
if(dNeq0(norm2(om2qu(qu2om(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'om2qu/qu2om,' if(dNeq0(norm2(om2qu(qu2om(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'om2qu/qu2om,'
if(dNeq0(norm2(eu2qu(qu2eu(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'eu2qu/qu2eu,' if(dNeq0(norm2(eu2qu(qu2eu(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'eu2qu/qu2eu,'
if(dNeq0(norm2(ax2qu(qu2ax(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'ax2qu/qu2ax,' if(dNeq0(norm2(ax2qu(qu2ax(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'ax2qu/qu2ax,'
if(dNeq0(norm2(ro2qu(qu2ro(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'ro2qu/qu2ro,' if(dNeq0(norm2(ro2qu(qu2ro(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'ro2qu/qu2ro,'
if(dNeq0(norm2(ho2qu(qu2ho(qu))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2qu/qu2ho,' if(dNeq0(norm2(ho2qu(qu2ho(qu))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2qu/qu2ho,'
if(dNeq0(norm2(cu2qu(qu2cu(qu))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2qu/qu2cu,' if(dNeq0(norm2(cu2qu(qu2cu(qu))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2qu/qu2cu,'
#endif
om = qu2om(qu) om = qu2om(qu)
#ifndef __PGI
if(dNeq0(norm2(om2qu(eu2om(om2eu(om)))-qu),1.0e-7_pReal)) msg = trim(msg)//'eu2om/om2eu,' if(dNeq0(norm2(om2qu(eu2om(om2eu(om)))-qu),1.0e-7_pReal)) msg = trim(msg)//'eu2om/om2eu,'
if(dNeq0(norm2(om2qu(ax2om(om2ax(om)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ax2om/om2ax,' if(dNeq0(norm2(om2qu(ax2om(om2ax(om)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ax2om/om2ax,'
if(dNeq0(norm2(om2qu(ro2om(om2ro(om)))-qu),1.0e-12_pReal)) msg = trim(msg)//'ro2om/om2ro,' if(dNeq0(norm2(om2qu(ro2om(om2ro(om)))-qu),1.0e-12_pReal)) msg = trim(msg)//'ro2om/om2ro,'
if(dNeq0(norm2(om2qu(ho2om(om2ho(om)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2om/om2ho,' if(dNeq0(norm2(om2qu(ho2om(om2ho(om)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2om/om2ho,'
if(dNeq0(norm2(om2qu(cu2om(om2cu(om)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2om/om2cu,' if(dNeq0(norm2(om2qu(cu2om(om2cu(om)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2om/om2cu,'
#endif
eu = qu2eu(qu) eu = qu2eu(qu)
#ifndef __PGI
if(dNeq0(norm2(eu2qu(ax2eu(eu2ax(eu)))-qu),1.0e-12_pReal)) msg = trim(msg)//'ax2eu/eu2ax,' if(dNeq0(norm2(eu2qu(ax2eu(eu2ax(eu)))-qu),1.0e-12_pReal)) msg = trim(msg)//'ax2eu/eu2ax,'
if(dNeq0(norm2(eu2qu(ro2eu(eu2ro(eu)))-qu),1.0e-12_pReal)) msg = trim(msg)//'ro2eu/eu2ro,' if(dNeq0(norm2(eu2qu(ro2eu(eu2ro(eu)))-qu),1.0e-12_pReal)) msg = trim(msg)//'ro2eu/eu2ro,'
if(dNeq0(norm2(eu2qu(ho2eu(eu2ho(eu)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2eu/eu2ho,' if(dNeq0(norm2(eu2qu(ho2eu(eu2ho(eu)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2eu/eu2ho,'
if(dNeq0(norm2(eu2qu(cu2eu(eu2cu(eu)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2eu/eu2cu,' if(dNeq0(norm2(eu2qu(cu2eu(eu2cu(eu)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2eu/eu2cu,'
#endif
ax = qu2ax(qu) ax = qu2ax(qu)
#ifndef __PGI
if(dNeq0(norm2(ax2qu(ro2ax(ax2ro(ax)))-qu),1.0e-12_pReal)) msg = trim(msg)//'ro2ax/ax2ro,' if(dNeq0(norm2(ax2qu(ro2ax(ax2ro(ax)))-qu),1.0e-12_pReal)) msg = trim(msg)//'ro2ax/ax2ro,'
if(dNeq0(norm2(ax2qu(ho2ax(ax2ho(ax)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2ax/ax2ho,' if(dNeq0(norm2(ax2qu(ho2ax(ax2ho(ax)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2ax/ax2ho,'
if(dNeq0(norm2(ax2qu(cu2ax(ax2cu(ax)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2ax/ax2cu,' if(dNeq0(norm2(ax2qu(cu2ax(ax2cu(ax)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2ax/ax2cu,'
#endif
ro = qu2ro(qu) ro = qu2ro(qu)
#ifndef __PGI
if(dNeq0(norm2(ro2qu(ho2ro(ro2ho(ro)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2ro/ro2ho,' if(dNeq0(norm2(ro2qu(ho2ro(ro2ho(ro)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2ro/ro2ho,'
if(dNeq0(norm2(ro2qu(cu2ro(ro2cu(ro)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2ro/ro2cu,' if(dNeq0(norm2(ro2qu(cu2ro(ro2cu(ro)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2ro/ro2cu,'
#endif
ho = qu2ho(qu) ho = qu2ho(qu)
#ifndef __PGI
if(dNeq0(norm2(ho2qu(cu2ho(ho2cu(ho)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2ho/ho2cu,' if(dNeq0(norm2(ho2qu(cu2ho(ho2cu(ho)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2ho/ho2cu,'
#endif
call R%fromMatrix(om) call R%fromMatrix(om)