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

Merge branch 'development' into 31_consistent_wd

This commit is contained in:
Martin Diehl 2018-08-03 00:08:16 +02:00
parent 7283ee6caa e0a39d202c
commit 3cceda5fe3
43 changed files with 183 additions and 1232 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
.gitlab-ci.yml
View File

@ -12,6 +12,7 @@ stages:
- compileAbaqus2017
- example
- performance
- createPackage
- createDocumentation
- saveDocumentation
- updateMaster
@ -411,6 +412,13 @@ SpectralRuntime:
- master
- release
###################################################################################################
createTar:
stage: createPackage
script:
- cd $(mktemp -d)
- $DAMASKROOT/PRIVATE/releasing/deployMe.sh $CI_COMMIT_SHA
###################################################################################################
AbaqusStd:
stage: createDocumentation

103
DAMASK_prerequisites.sh
View File

@ -12,21 +12,38 @@ echo + Send to damask@mpie.de for support
echo + view with \'cat $OUTFILE\'
echo ===========================================
function firstLevel {
echo -e '\n\n=============================================================================================='
echo $1
echo ==============================================================================================
}
function secondLevel {
echo ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
echo $1
echo ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
}
function thirdLevel {
echo -e '\n----------------------------------------------------------------------------------------------'
echo $1
echo ----------------------------------------------------------------------------------------------
}
function getDetails {
if which $1 &> /dev/null; then
echo ----------------------------------------------------------------------------------------------
echo $1:
echo ----------------------------------------------------------------------------------------------
secondLevel $1:
echo + location:
which $1
echo + $1 $2:
$1 $2
echo -e '\n'
else
echo $1 not found
fi
echo
}
# redirect STDOUT and STDERR to logfile
# https://stackoverflow.com/questions/11229385/redirect-all-output-in-a-bash-script-when-using-set-x^
exec > $OUTFILE 2>&1
@ -38,28 +55,18 @@ DAMASK_ROOT="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )"
echo XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
echo System report for \'$(hostname)\' created on $(date '+%Y-%m-%d %H:%M:%S') by \'$(whoami)\'
echo XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
echo
echo ==============================================================================================
echo DAMASK settings
echo ==============================================================================================
echo ----------------------------------------------------------------------------------------------
echo DAMASK_ROOT:
echo ----------------------------------------------------------------------------------------------
firstLevel "DAMASK settings"
secondLevel "DAMASK_ROOT"
echo $DAMASK_ROOT
echo
echo ----------------------------------------------------------------------------------------------
echo Version:
echo ----------------------------------------------------------------------------------------------
secondLevel "Version"
cat VERSION
echo
echo ----------------------------------------------------------------------------------------------
echo Settings in CONFIG:
echo ----------------------------------------------------------------------------------------------
secondLevel "Settings in CONFIG"
cat CONFIG
echo
echo ==============================================================================================
echo System
echo ==============================================================================================
firstLevel "System"
uname -a
echo
echo PATH: $PATH
@ -69,74 +76,52 @@ echo SHELL: $SHELL
echo PETSC_ARCH: $PETSC_ARCH
echo PETSC_DIR: $PETSC_DIR
ls $PETSC_DIR/lib
echo
echo ==============================================================================================
echo Python
echo ==============================================================================================
firstLevel "Python"
DEFAULT_PYTHON=python2.7
for executable in python python2 python3 python2.7; do
getDetails $executable '--version'
done
echo ----------------------------------------------------------------------------------------------
echo Details on $DEFAULT_PYTHON:
echo ----------------------------------------------------------------------------------------------
secondLevel "Details on $DEFAULT_PYTHON:"
echo $(ls -la $(which $DEFAULT_PYTHON))
for module in numpy scipy;do
echo -e '\n----------------------------------------------------------------------------------------------'
echo $module
echo ----------------------------------------------------------------------------------------------
thirdLevel $module
$DEFAULT_PYTHON -c "import $module; \
print('Version: {}'.format($module.__version__)); \
print('Location: {}'.format($module.__file__))"
done
echo ----------------------------------------------------------------------------------------------
echo vtk
echo ----------------------------------------------------------------------------------------------
thirdLevel vtk
$DEFAULT_PYTHON -c "import vtk; \
print('Version: {}'.format(vtk.vtkVersion.GetVTKVersion())); \
print('Location: {}'.format(vtk.__file__))"
echo ----------------------------------------------------------------------------------------------
echo h5py
echo ----------------------------------------------------------------------------------------------
thirdLevel h5py
$DEFAULT_PYTHON -c "import h5py; \
print('Version: {}'.format(h5py.version.version)); \
print('Location: {}'.format(h5py.__file__))"
echo
echo ==============================================================================================
echo GCC
echo ==============================================================================================
firstLevel "GNU Compiler Collection"
for executable in gcc g++ gfortran ;do
getDetails $executable '--version'
done
echo
echo ==============================================================================================
echo Intel Compiler Suite
echo ==============================================================================================
firstLevel "Intel Compiler Suite"
for executable in icc icpc ifort ;do
getDetails $executable '--version'
done
echo
echo ==============================================================================================
echo MPI Wrappers
echo ==============================================================================================
firstLevel "MPI Wrappers"
for executable in mpicc mpiCC mpic++ mpicpc mpicxx mpifort mpif90 mpif77; do
getDetails $executable '-show'
done
echo
echo ==============================================================================================
echo MPI Launchers
echo ==============================================================================================
firstLevel "MPI Launchers"
for executable in mpirun mpiexec; do
getDetails $executable '--version'
done
echo
echo ==============================================================================================
echo Abaqus
echo ==============================================================================================
firstLevel "Abaqus"
cd installation/mods_Abaqus # to have the right environment file
for executable in abaqus abq2016 abq2017; do
for executable in abaqus abq2017 abq2018; do
getDetails $executable 'information=all'
done
cd ../..

2
VERSION
View File

@ -1 +1 @@
v2.0.2-208-g798a9720
v2.0.2-261-gbc3f6ae9

1012
installation/patch/PETSc3.8
View File

File diff suppressed because it is too large Load Diff

3
installation/patch/README.md
View File

@ -14,9 +14,6 @@ patch -p1 < installation/patch/nameOfPatch
* **fwbw_derivative** switches the default spatial derivative from continuous to forward/backward difference.
This generally reduces spurious oscillations in the result as the spatial accuracy of the derivative is then compatible with the underlying solution grid.
* **PETSc-3.8** adjusts all includes and calls to PETSc to follow the 3.8.x API.
This allows to use the most recent version of PETSc.
## Create patch
commit your changes

8
installation/patch/python2to3.sh Executable file
View File

@ -0,0 +1,8 @@
#! /usr/bin/env bash
if [ $1x != 3to2x ]; then
echo 'python2.7 to python'
find . -name '*.py' | xargs sed -i 's/usr\/bin\/env python2.7/usr\/bin\/env python/g'
else
echo 'python to python2.7'
find . -name '*.py' | xargs sed -i 's/usr\/bin\/env python/usr\/bin\/env python2.7/g'
fi

24
lib/damask/util.py
View File

@ -1,5 +1,4 @@
# -*- coding: UTF-8 no BOM -*-
import sys,time,random,threading,os,subprocess,shlex
import numpy as np
from optparse import Option
@ -37,8 +36,8 @@ class bcolors:
def srepr(arg,glue = '\n'):
"""Joins arguments as individual lines"""
if (not hasattr(arg, "strip") and
hasattr(arg, "__getitem__") or
hasattr(arg, "__iter__")):
(hasattr(arg, "__getitem__") or
hasattr(arg, "__iter__"))):
return glue.join(str(x) for x in arg)
return arg if isinstance(arg,str) else repr(arg)
@ -60,9 +59,9 @@ def report_geom(info,
what = ['grid','size','origin','homogenization','microstructures']):
"""Reports (selected) geometry information"""
output = {
'grid' : 'grid a b c: {}'.format(' x '.join(map(str,info['grid' ]))),
'size' : 'size x y z: {}'.format(' x '.join(map(str,info['size' ]))),
'origin' : 'origin x y z: {}'.format(' : '.join(map(str,info['origin']))),
'grid' : 'grid a b c: {}'.format(' x '.join(list(map(str,info['grid' ])))),
'size' : 'size x y z: {}'.format(' x '.join(list(map(str,info['size' ])))),
'origin' : 'origin x y z: {}'.format(' : '.join(list(map(str,info['origin'])))),
'homogenization' : 'homogenization: {}'.format(info['homogenization']),
'microstructures' : 'microstructures: {}'.format(info['microstructures']),
}
@ -94,20 +93,21 @@ def execute(cmd,
stdout = subprocess.PIPE,
stderr = subprocess.PIPE,
stdin = subprocess.PIPE)
out,error = [i.replace("\x08","") for i in (process.communicate() if streamIn is None
else process.communicate(streamIn.read()))]
out,error = [i for i in (process.communicate() if streamIn is None
else process.communicate(streamIn.read().encode('utf-8')))]
out = out.decode('utf-8').replace('\x08','')
error = error.decode('utf-8').replace('\x08','')
os.chdir(initialPath)
if process.returncode != 0: raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode))
return out,error
def coordGridAndSize(coordinates):
"""Determines grid count and overall physical size along each dimension of an ordered array of coordinates"""
dim = coordinates.shape[1]
coords = [np.unique(coordinates[:,i]) for i in range(dim)]
mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords),'i')
mincorner = np.array(list(map(min,coords)))
maxcorner = np.array(list(map(max,coords)))
grid = np.array(list(map(len,coords)),'i')
size = grid/np.maximum(np.ones(dim,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other ones
return grid,size

4
processing/post/addCauchy.py
View File

@ -75,8 +75,8 @@ for name in filenames:
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
F = np.array(map(float,table.data[column[options.defgrad]:column[options.defgrad]+9]),'d').reshape(3,3)
P = np.array(map(float,table.data[column[options.stress ]:column[options.stress ]+9]),'d').reshape(3,3)
F = np.array(list(map(float,table.data[column[options.defgrad]:column[options.defgrad]+9])),'d').reshape(3,3)
P = np.array(list(map(float,table.data[column[options.stress ]:column[options.stress ]+9])),'d').reshape(3,3)
table.data_append(list(1.0/np.linalg.det(F)*np.dot(P,F.T).reshape(9))) # [Cauchy] = (1/det(F)) * [P].[F_transpose]
outputAlive = table.data_write() # output processed line

9
processing/post/addCompatibilityMismatch.py
View File

@ -282,19 +282,12 @@ for name in filenames:
table.data_readArray([options.defgrad,options.pos])
table.data_rewind()
if len(table.data.shape) < 2: table.data.shape += (1,) # expand to 2D shape
if table.data[:,9:].shape[1] < 3:
table.data = np.hstack((table.data,
np.zeros((table.data.shape[0],
3-table.data[:,9:].shape[1]),dtype='f'))) # fill coords up to 3D with zeros
coords = [np.unique(table.data[:,9+i]) for i in range(3)]
mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 set to smallest among other spacings
grid,size = damask.util.coordGridAndSize(table.data[:,9:12])
N = grid.prod()
if N != len(table.data): errors.append('data count {} does not match grid {}x{}x{}.'.format(N,*grid))

1
processing/post/addCurl.py
View File

@ -138,7 +138,6 @@ for name in filenames:
# --------------- figure out size and grid ---------------------------------------------------------
table.data_readArray()
grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)])
# ------------------------------------------ process value field -----------------------------------

9
processing/post/addDeterminant.py
View File

@ -58,7 +58,7 @@ for name in filenames:
errors = []
remarks = []
for type, data in items.iteritems():
for type, data in items.items():
for what in data['labels']:
dim = table.label_dimension(what)
if dim != data['dim']: remarks.append('column {} is not a {}...'.format(what,type))
@ -81,12 +81,11 @@ for name in filenames:
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
for type, data in items.iteritems():
for type, data in items.items():
for column in data['column']:
table.data_append(determinant(map(float,table.data[column:
column+data['dim']])))
table.data_append(determinant(list(map(float,table.data[column: column+data['dim']]))))
outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization -----------------------------------
table.close() # close input ASCII table (works for stdin)
table.close() # close input ASCII table (works for stdin)

10
processing/post/addDeviator.py
View File

@ -66,7 +66,7 @@ for name in filenames:
remarks = []
column = {}
for type, data in items.iteritems():
for type, data in items.items():
for what in data['labels']:
dim = table.label_dimension(what)
if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type))
@ -83,7 +83,7 @@ for name in filenames:
# ------------------------------------------ assemble header --------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
for type, data in items.iteritems():
for type, data in items.items():
for label in data['active']:
table.labels_append(['{}_dev({})'.format(i+1,label) for i in range(data['dim'])] + \
(['sph({})'.format(label)] if options.spherical else [])) # extend ASCII header with new labels
@ -93,10 +93,10 @@ for name in filenames:
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
for type, data in items.iteritems():
for type, data in items.items():
for column in data['column']:
table.data_append(deviator(map(float,table.data[column:
column+data['dim']]),options.spherical))
table.data_append(deviator(list(map(float,table.data[column:
column+data['dim']])),options.spherical))
outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization -----------------------------------

8
processing/post/addDisplacement.py
View File

@ -168,13 +168,7 @@ for name in filenames:
np.zeros((table.data.shape[0],
3-table.data[:,9:].shape[1]),dtype='f'))) # fill coords up to 3D with zeros
coords = [np.unique(table.data[:,9+i]) for i in range(3)]
mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 set to smallest among other spacings
grid,size = damask.util.coordGridAndSize(table.data[:,9:12])
N = grid.prod()
if N != len(table.data): errors.append('data count {} does not match grid {}x{}x{}.'.format(N,*grid))

4
processing/post/addEhkl.py
View File

@ -88,9 +88,9 @@ for name in filenames:
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
for column in columns:
table.data_append(E_hkl(map(float,table.data[column:column+3]),options.hkl))
table.data_append(E_hkl(list(map(float,table.data[column:column+3])),options.hkl))
outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization -----------------------------------
table.close() # close ASCII tables
table.close() # close ASCII tables

12
processing/post/addEuclideanDistance.py
View File

@ -102,7 +102,7 @@ parser.add_option('-t',
help = 'feature type {{{}}} '.format(', '.join(map(lambda x:'/'.join(x['names']),features))) )
parser.add_option('-n',
'--neighborhood',
dest = 'neighborhood', choices = neighborhoods.keys(), metavar = 'string',
dest = 'neighborhood', choices = list(neighborhoods.keys()), metavar = 'string',
help = 'neighborhood type [neumann] {{{}}}'.format(', '.join(neighborhoods.keys())))
parser.add_option('-s',
'--scale',
@ -151,10 +151,8 @@ for name in filenames:
remarks = []
column = {}
coordDim = table.label_dimension(options.pos)
if not 3 >= coordDim >= 1:
if not 3 >= table.label_dimension(options.pos) >= 1:
errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos))
else: coordCol = table.label_index(options.pos)
if table.label_dimension(options.id) != 1: errors.append('grain identifier {} not found.'.format(options.id))
else: idCol = table.label_index(options.id)
@ -178,11 +176,7 @@ for name in filenames:
table.data_readArray()
coords = [np.unique(table.data[:,coordCol+i]) for i in range(coordDim)]
mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords)+[1]*(3-len(coords)),'i')
grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)])
N = grid.prod()
if N != len(table.data): errors.append('data count {} does not match grid {}.'.format(N,'x'.join(map(str,grid))))

6
processing/post/addGaussian.py
View File

@ -83,7 +83,7 @@ for name in filenames:
if table.label_dimension(options.pos) != 3: errors.append('coordinates {} are not a vector.'.format(options.pos))
else: colCoord = table.label_index(options.pos)
for type, data in items.iteritems():
for type, data in items.items():
for what in (data['labels'] if data['labels'] is not None else []):
dim = table.label_dimension(what)
if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type))
@ -100,7 +100,7 @@ for name in filenames:
# ------------------------------------------ assemble header --------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
for type, data in items.iteritems():
for type, data in items.items():
for label in data['active']:
table.labels_append(['Gauss{}({})'.format(options.sigma,label)]) # extend ASCII header with new labels
table.head_write()
@ -114,7 +114,7 @@ for name in filenames:
# ------------------------------------------ process value field -----------------------------------
stack = [table.data]
for type, data in items.iteritems():
for type, data in items.items():
for i,label in enumerate(data['active']):
stack.append(ndimage.filters.gaussian_filter(table.data[:,data['column'][i]],
options.sigma,options.order,

10
processing/post/addIPFcolor.py
View File

@ -116,18 +116,18 @@ for name in filenames:
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
if inputtype == 'eulers':
o = damask.Orientation(Eulers = np.array(map(float,table.data[column:column+3]))*toRadians,
o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians,
symmetry = options.symmetry).reduced()
elif inputtype == 'matrix':
o = damask.Orientation(matrix = np.array(map(float,table.data[column:column+9])).reshape(3,3).transpose(),
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3).transpose(),
symmetry = options.symmetry).reduced()
elif inputtype == 'frame':
o = damask.Orientation(matrix = np.array(map(float,table.data[column[0]:column[0]+3] + \
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3])).reshape(3,3),
table.data[column[2]:column[2]+3]))).reshape(3,3),
symmetry = options.symmetry).reduced()
elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = np.array(map(float,table.data[column:column+4])),
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),
symmetry = options.symmetry).reduced()
table.data_append(o.IPFcolor(pole))

4
processing/post/addMises.py
View File

@ -70,7 +70,7 @@ for name in filenames:
errors = []
remarks = []
for type, data in items.iteritems():
for type, data in items.items():
for what in data['labels']:
dim = table.label_dimension(what)
if dim != data['dim']: remarks.append('column {} is not a {}...'.format(what,type))
@ -94,7 +94,7 @@ for name in filenames:
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
for type, data in items.iteritems():
for type, data in items.items():
for column in data['column']:
table.data_append(Mises(type,
np.array(table.data[column:column+data['dim']],'d').reshape(data['shape'])))

16
processing/post/addOrientations.py
View File

@ -80,7 +80,7 @@ parser.set_defaults(output = [],
(options, filenames) = parser.parse_args()
options.output = map(lambda x: x.lower(), options.output)
options.output = list(map(lambda x: x.lower(), options.output))
if options.output == [] or (not set(options.output).issubset(set(outputChoices))):
parser.error('output must be chosen from {}.'.format(', '.join(outputChoices)))
@ -147,21 +147,21 @@ for name in filenames:
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
if inputtype == 'eulers':
o = damask.Orientation(Eulers = np.array(map(float,table.data[column:column+3]))*toRadians,
o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians,
symmetry = options.symmetry).reduced()
elif inputtype == 'rodrigues':
o = damask.Orientation(Rodrigues= np.array(map(float,table.data[column:column+3])),
o = damask.Orientation(Rodrigues= np.array(list(map(float,table.data[column:column+3]))),
symmetry = options.symmetry).reduced()
elif inputtype == 'matrix':
o = damask.Orientation(matrix = np.array(map(float,table.data[column:column+9])).reshape(3,3).transpose(),
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3).transpose(),
symmetry = options.symmetry).reduced()
elif inputtype == 'frame':
o = damask.Orientation(matrix = np.array(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3])).reshape(3,3),
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3]))).reshape(3,3),
symmetry = options.symmetry).reduced()
elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = np.array(map(float,table.data[column:column+4])),
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),
symmetry = options.symmetry).reduced()
o.quaternion = r*o.quaternion*R # apply additional lab and crystal frame rotations

4
processing/post/addPK2.py
View File

@ -75,8 +75,8 @@ for name in filenames:
# ------------------------------------------ process data ------------------------------------------
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
F = np.array(map(float,table.data[column[options.defgrad]:column[options.defgrad]+9]),'d').reshape(3,3)
P = np.array(map(float,table.data[column[options.stress ]:column[options.stress ]+9]),'d').reshape(3,3)
F = np.array(list(map(float,table.data[column[options.defgrad]:column[options.defgrad]+9])),'d').reshape(3,3)
P = np.array(list(map(float,table.data[column[options.stress ]:column[options.stress ]+9])),'d').reshape(3,3)
table.data_append(list(np.dot(np.linalg.inv(F),P).reshape(9))) # [S] =[P].[F-1]
outputAlive = table.data_write() # output processed line

12
processing/post/addPole.py
View File

@ -120,15 +120,15 @@ for name in filenames:
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
if inputtype == 'eulers':
o = damask.Orientation(Eulers = np.array(map(float,table.data[column:column+3]))*toRadians)
o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians)
elif inputtype == 'matrix':
o = damask.Orientation(matrix = np.array(map(float,table.data[column:column+9])).reshape(3,3).transpose())
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3).transpose())
elif inputtype == 'frame':
o = damask.Orientation(matrix = np.array(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3])).reshape(3,3))
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3]))).reshape(3,3))
elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = np.array(map(float,table.data[column:column+4])))
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))))
rotatedPole = o.quaternion*pole # rotate pole according to crystal orientation
(x,y) = rotatedPole[0:2]/(1.+abs(pole[2])) # stereographic projection

12
processing/post/addSchmidfactors.py
View File

@ -252,15 +252,15 @@ for name in filenames:
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
if inputtype == 'eulers':
o = damask.Orientation(Eulers = np.array(map(float,table.data[column:column+3]))*toRadians,)
o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians,)
elif inputtype == 'matrix':
o = damask.Orientation(matrix = np.array(map(float,table.data[column:column+9])).reshape(3,3).transpose(),)
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3).transpose(),)
elif inputtype == 'frame':
o = damask.Orientation(matrix = np.array(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3])).reshape(3,3),)
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3]))).reshape(3,3),)
elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = np.array(map(float,table.data[column:column+4])),)
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),)
rotForce = o.quaternion.conjugated() * force
rotNormal = o.quaternion.conjugated() * normal

6
processing/post/addSpectralDecomposition.py
View File

@ -58,7 +58,7 @@ for name in filenames:
errors = []
remarks = []
for type, data in items.iteritems():
for type, data in items.items():
for what in data['labels']:
dim = table.label_dimension(what)
if dim != data['dim']: remarks.append('column {} is not a {}...'.format(what,type))
@ -84,9 +84,9 @@ for name in filenames:
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
for type, data in items.iteritems():
for type, data in items.items():
for column in data['column']:
(u,v) = np.linalg.eigh(np.array(map(float,table.data[column:column+data['dim']])).reshape(data['shape']))
(u,v) = np.linalg.eigh(np.array(list(map(float,table.data[column:column+data['dim']]))).reshape(data['shape']))
if options.rh and np.dot(np.cross(v[:,0], v[:,1]), v[:,2]) < 0.0 : v[:, 2] *= -1.0 # ensure right-handed eigenvector basis
table.data_append(list(u)) # vector of max,mid,min eigval
table.data_append(list(v.transpose().reshape(data['dim']))) # 3x3=9 combo vector of max,mid,min eigvec coordinates

4
processing/post/addStrainTensors.py
View File

@ -101,7 +101,7 @@ for name in filenames:
errors = []
remarks = []
for type, data in items.iteritems():
for type, data in items.items():
for what in data['labels']:
dim = table.label_dimension(what)
if dim != data['dim']: remarks.append('column {} is not a {}...'.format(what,type))
@ -132,7 +132,7 @@ for name in filenames:
while outputAlive and table.data_read(): # read next data line of ASCII table
for column in items['tensor']['column']: # loop over all requested defgrads
F = np.array(map(float,table.data[column:column+items['tensor']['dim']]),'d').reshape(items['tensor']['shape'])
F = np.array(list(map(float,table.data[column:column+items['tensor']['dim']])),'d').reshape(items['tensor']['shape'])
(U,S,Vh) = np.linalg.svd(F) # singular value decomposition
R = np.dot(U,Vh) # rotation of polar decomposition
stretch['U'] = np.dot(np.linalg.inv(R),F) # F = RU

25
processing/post/averageDown.py
View File

@ -76,7 +76,6 @@ for name in filenames:
remarks = []
if table.label_dimension(options.pos) != 3: errors.append('coordinates {} are not a vector.'.format(options.pos))
else: colCoord = table.label_index(options.pos)
if remarks != []: damask.util.croak(remarks)
if errors != []:
@ -94,14 +93,7 @@ for name in filenames:
table.data_readArray()
if (any(options.grid) == 0 or any(options.size) == 0.0):
coords = [np.unique(table.data[:,colCoord+i]) for i in range(3)]
mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 set to smallest among other spacings
delta = size/np.maximum(np.ones(3,'d'), grid)
origin = mincorner - 0.5*delta # shift from cell center to corner
grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)])
else:
grid = np.array(options.grid,'i')
@ -129,16 +121,15 @@ for name in filenames:
#--- generate grid --------------------------------------------------------------------------------
if colCoord:
x = (0.5 + shift[0] + np.arange(packedGrid[0],dtype=float))/packedGrid[0]*size[0] + origin[0]
y = (0.5 + shift[1] + np.arange(packedGrid[1],dtype=float))/packedGrid[1]*size[1] + origin[1]
z = (0.5 + shift[2] + np.arange(packedGrid[2],dtype=float))/packedGrid[2]*size[2] + origin[2]
x = (0.5 + shift[0] + np.arange(packedGrid[0],dtype=float))/packedGrid[0]*size[0] + origin[0]
y = (0.5 + shift[1] + np.arange(packedGrid[1],dtype=float))/packedGrid[1]*size[1] + origin[1]
z = (0.5 + shift[2] + np.arange(packedGrid[2],dtype=float))/packedGrid[2]*size[2] + origin[2]
xx = np.tile( x, packedGrid[1]* packedGrid[2])
yy = np.tile(np.repeat(y,packedGrid[0] ),packedGrid[2])
zz = np.repeat(z,packedGrid[0]*packedGrid[1])
xx = np.tile( x, packedGrid[1]* packedGrid[2])
yy = np.tile(np.repeat(y,packedGrid[0] ),packedGrid[2])
zz = np.repeat(z,packedGrid[0]*packedGrid[1])
table.data[:,colCoord:colCoord+3] = np.squeeze(np.dstack((xx,yy,zz)))
table.data[:,table.label_indexragen(options.pos)] = np.squeeze(np.dstack((xx,yy,zz)))
# ------------------------------------------ output result -----------------------------------------

10
processing/post/blowUp.py
View File

@ -64,7 +64,6 @@ for name in filenames:
remarks = []
if table.label_dimension(options.pos) != 3: errors.append('coordinates "{}" are not a vector.'.format(options.pos))
else: colCoord = table.label_index(options.pos)
colElem = table.label_index('elem')
@ -79,12 +78,7 @@ for name in filenames:
table.data_readArray(options.pos)
table.data_rewind()
coords = [np.unique(table.data[:,i]) for i in range(3)]
mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords),'i')
size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 set to smallest among other spacings
grid,size = damask.util.coordGridAndSize(table.data)
packing = np.array(options.packing,'i')
outSize = grid*packing
@ -113,7 +107,7 @@ for name in filenames:
for c in range(outSize[2]):
for b in range(outSize[1]):
for a in range(outSize[0]):
data[a,b,c,colCoord:colCoord+3] = [a+0.5,b+0.5,c+0.5]*elementSize
data[a,b,c,table.label_indexrange(options.pos)] = [a+0.5,b+0.5,c+0.5]*elementSize
if colElem != -1: data[a,b,c,colElem] = elem
table.data = data[a,b,c,:].tolist()
outputAlive = table.data_write() # output processed line

6
processing/post/rotateData.py
View File

@ -73,7 +73,7 @@ for name in filenames:
remarks = []
column = {}
for type, data in items.iteritems():
for type, data in items.items():
for what in data['labels']:
dim = table.label_dimension(what)
if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type))
@ -100,13 +100,13 @@ for name in filenames:
for column in items[datatype]['column']: # loop over all requested labels
table.data[column:column+items[datatype]['dim']] = \
q * np.array(map(float,table.data[column:column+items[datatype]['dim']]))
q * np.array(list(map(float,table.data[column:column+items[datatype]['dim']])))
datatype = 'tensor'
for column in items[datatype]['column']: # loop over all requested labels
table.data[column:column+items[datatype]['dim']] = \
np.dot(R,np.dot(np.array(map(float,table.data[column:column+items[datatype]['dim']])).\
np.dot(R,np.dot(np.array(list(map(float,table.data[column:column+items[datatype]['dim']]))).\
reshape(items[datatype]['shape']),R.transpose())).reshape(items[datatype]['dim'])
outputAlive = table.data_write() # output processed line

2
processing/post/vtk2ang.py
View File

@ -421,8 +421,6 @@ for filename in filenames:
meshActor.GetProperty().SetOpacity(0.2)
meshActor.GetProperty().SetColor(1.0,1.0,0)
meshActor.GetProperty().BackfaceCullingOn()
# meshActor.GetProperty().SetEdgeColor(1,1,0.5)
# meshActor.GetProperty().EdgeVisibilityOn()
boxpoints = vtk.vtkPoints()
for n in range(8):

2
processing/post/vtk_rectilinearGrid.py
View File

@ -82,7 +82,7 @@ for name in filenames:
[coords[i][j-1] + coords[i][j] for j in range(1,len(coords[i]))] + \
[3.0 * coords[i][-1] - coords[i][-1 - int(len(coords[i]) > 1)]]) for i in range(3)]
grid = np.array(map(len,coords),'i')
grid = np.array(list(map(len,coords)),'i')
N = grid.prod() if options.mode == 'point' else (grid-1).prod()
if N != len(table.data):

12
processing/pre/geom_fromOsteonGeometry.py
View File

@ -55,9 +55,9 @@ parser.set_defaults(canal = 25e-6,
(options,filename) = parser.parse_args()
if np.any(options.grid < 2):
if np.any(np.array(options.grid) < 2):
parser('invalid grid a b c.')
if np.any(options.size <= 0.0):
if np.any(np.array(options.size) <= 0.0):
parser('invalid size x y z.')
# --- open input files ----------------------------------------------------------------------------
@ -114,12 +114,8 @@ for y in range(info['grid'][1]):
info['microstructures'] += 1
#--- report ---------------------------------------------------------------------------------------
damask.util.croak(['grid a b c: %s'%(' x '.join(map(str,info['grid']))),
'size x y z: %s'%(' x '.join(map(str,info['size']))),
'origin x y z: %s'%(' : '.join(map(str,info['origin']))),
'homogenization: %i'%info['homogenization'],
'microstructures: %i'%info['microstructures']])
# -------------------------------------- switch according to task ----------------------------------
damask.util.report_geom(info,['grid','size','origin','homogenization','microstructures'])
formatwidth = 1+int(math.floor(math.log10(info['microstructures']-1)))
header = [scriptID + ' ' + ' '.join(sys.argv[1:])]
header.append('<microstructure>')

10
processing/pre/geom_fromTable.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*-
import os,sys,math,types,time
import os,sys,math,time
import scipy.spatial, numpy as np
from optparse import OptionParser
import damask
@ -152,7 +152,7 @@ for name in filenames:
continue
table.data_readArray([options.pos] \
+ ([label] if isinstance(label, types.StringTypes) else label) \
+ (label if isinstance(label, list) else [label]) \
+ ([options.phase] if options.phase else []))
if coordDim == 2:
@ -165,9 +165,9 @@ for name in filenames:
# --------------- figure out size and grid ---------------------------------------------------------
coords = [np.unique(table.data[:,i]) for i in range(3)]
mincorner = np.array(map(min,coords))
maxcorner = np.array(map(max,coords))
grid = np.array(map(len,coords),'i')
mincorner = np.array(list(map(min,coords)))
maxcorner = np.array(list(map(max,coords)))
grid = np.array(list(map(len,coords)),'i')
size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1)
size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 set to smallest among other spacings
delta = size/np.maximum(np.ones(3,'d'), grid)

3
processing/pre/geom_fromVoronoiTessellation.py
View File

@ -15,8 +15,7 @@ scriptID = ' '.join([scriptName,damask.version])
def meshgrid2(*arrs):
"""Code inspired by http://stackoverflow.com/questions/1827489/numpy-meshgrid-in-3d"""
arrs = tuple(reversed(arrs))
arrs = tuple(arrs)
lens = np.array(map(len, arrs))
lens = np.array(list(map(len, arrs)))
dim = len(arrs)
ans = []
for i, arr in enumerate(arrs):

2
processing/pre/geom_translate.py
View File

@ -92,7 +92,7 @@ for name in filenames:
}
substituted = np.copy(microstructure)
for k, v in sub.iteritems(): substituted[microstructure==k] = v # substitute microstructure indices
for k, v in sub.items(): substituted[microstructure==k] = v # substitute microstructure indices
substituted += options.microstructure # shift microstructure indices

2
processing/pre/hybridIA_linODFsampling.py
View File

@ -270,7 +270,7 @@ for name in filenames:
ODF['limit'] = np.radians(limits[1,:]) # right hand limits in radians
ODF['center'] = 0.0 if all(limits[0,:]<1e-8) else 0.5 # vertex or cell centered
ODF['interval'] = np.array(map(len,[np.unique(table.data[:,i]) for i in range(3)]),'i') # steps are number of distict values
ODF['interval'] = np.array(list(map(len,[np.unique(table.data[:,i]) for i in range(3)])),'i') # steps are number of distict values
ODF['nBins'] = ODF['interval'].prod()
ODF['delta'] = np.radians(np.array(limits[1,0:3]-limits[0,0:3])/(ODF['interval']-1)) # step size

13
processing/pre/mentat_spectralBox.py
View File

@ -77,7 +77,14 @@ def mesh(r,d):
"%f %f %f"%(-d[0],d[1],d[2]),
"%f %f %f"%(-d[0],d[1],0.0),
"*add_elements",
range(1,9),
"1",
"2",
"3",
"4",
"5",
"6",
"7",
"8",
"*sub_divisions",
"%i %i %i"%(r[2],r[1],r[0]),
"*subdivide_elements",
@ -201,7 +208,7 @@ if options.port:
except:
parser.error('no valid Mentat release found.')
# --- loop over input files -------------------------------------------------------------------------
# --- loop over input files ------------------------------------------------------------------------
if filenames == []: filenames = [None]
@ -236,7 +243,7 @@ for name in filenames:
# --- read data ------------------------------------------------------------------------------------
microstructure = table.microstructure_read(info['grid']).reshape(info['grid'].prod(),order='F') # read microstructure
microstructure = table.microstructure_read(info['grid']).reshape(info['grid'].prod(),order='F') # read microstructure
cmds = [\
init(),

2
processing/pre/patchFromReconstructedBoundaries.py
View File

@ -344,7 +344,7 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn):
else:
myNeighbors[grainNeighbors[leg][side]] = 1
if myNeighbors: # do I have any neighbors (i.e., non-bounding box segment)
candidateGrains = sorted(myNeighbors.iteritems(), key=lambda p: (p[1],p[0]), reverse=True) # sort grain counting
candidateGrains = sorted(myNeighbors.items(), key=lambda p: (p[1],p[0]), reverse=True) # sort grain counting
# most frequent one not yet seen?
rcData['grainMapping'].append(candidateGrains[0 if candidateGrains[0][0] not in rcData['grainMapping'] else 1][0]) # must be me then
# special case of bi-crystal situation...

4
src/DAMASK_marc.f90
View File

@ -216,8 +216,8 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
! Marc common blocks are in fixed format so they have to be reformated to free format (f90)
! Beware of changes in newer Marc versions
#include QUOTE(PASTE(../lib/MarcInclude/concom,Marc4DAMASK)) ! concom is needed for inc, lovl
#include QUOTE(PASTE(../lib/MarcInclude/creeps,Marc4DAMASK)) ! creeps is needed for timinc (time increment)
#include QUOTE(PASTE(./MarcInclude/concom,Marc4DAMASK)) ! concom is needed for inc, lovl
#include QUOTE(PASTE(./MarcInclude/creeps,Marc4DAMASK)) ! creeps is needed for timinc (time increment)
logical :: cutBack
real(pReal), dimension(6) :: stress

0
lib/MarcInclude/concom2016 → src/MarcInclude/concom2016
View File

0
lib/MarcInclude/concom2017 → src/MarcInclude/concom2017
View File

0
lib/MarcInclude/creeps2016 → src/MarcInclude/creeps2016
View File

0
lib/MarcInclude/creeps2017 → src/MarcInclude/creeps2017
View File

1
src/config.f90
View File

@ -670,6 +670,7 @@ function getStrings(this,key,defaultVal,requiredShape,raw)
endif
else notAllocated
if (whole) then
str = item%string%val(item%string%pos(4):)
getStrings = [getStrings,str]
else
do i=2_pInt,item%string%pos(1)

30
src/math.f90
View File

@ -160,7 +160,7 @@ module math
math_rotate_forward33, &
math_rotate_backward33, &
math_rotate_forward3333, &
math_limit
math_clip
private :: &
math_check, &
halton
@ -1363,16 +1363,16 @@ pure function math_RtoEuler(R)
sqhk =sqrt(R(1,3)*R(1,3)+R(2,3)*R(2,3))
! calculate PHI
math_RtoEuler(2) = acos(math_limit(R(3,3)/sqhkl,-1.0_pReal, 1.0_pReal))
math_RtoEuler(2) = acos(math_clip(R(3,3)/sqhkl,-1.0_pReal, 1.0_pReal))
if((math_RtoEuler(2) < 1.0e-8_pReal) .or. (pi-math_RtoEuler(2) < 1.0e-8_pReal)) then
math_RtoEuler(3) = 0.0_pReal
math_RtoEuler(1) = acos(math_limit(R(1,1)/squvw, -1.0_pReal, 1.0_pReal))
math_RtoEuler(1) = acos(math_clip(R(1,1)/squvw, -1.0_pReal, 1.0_pReal))
if(R(2,1) > 0.0_pReal) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1)
else
math_RtoEuler(3) = acos(math_limit(R(2,3)/sqhk, -1.0_pReal, 1.0_pReal))
math_RtoEuler(3) = acos(math_clip(R(2,3)/sqhk, -1.0_pReal, 1.0_pReal))
if(R(1,3) < 0.0) math_RtoEuler(3) = 2.0_pReal*pi-math_RtoEuler(3)
math_RtoEuler(1) = acos(math_limit(-R(3,2)/sin(math_RtoEuler(2)), -1.0_pReal, 1.0_pReal))
math_RtoEuler(1) = acos(math_clip(-R(3,2)/sin(math_RtoEuler(2)), -1.0_pReal, 1.0_pReal))
if(R(3,1) < 0.0) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1)
end if
@ -1654,7 +1654,7 @@ pure function math_qToEuler(qPassive)
math_qToEuler(2) = acos(1.0_pReal-2.0_pReal*(q(2)**2+q(3)**2))
if (abs(math_qToEuler(2)) < 1.0e-6_pReal) then
math_qToEuler(1) = sign(2.0_pReal*acos(math_limit(q(1),-1.0_pReal, 1.0_pReal)),q(4))
math_qToEuler(1) = sign(2.0_pReal*acos(math_clip(q(1),-1.0_pReal, 1.0_pReal)),q(4))
math_qToEuler(3) = 0.0_pReal
else
math_qToEuler(1) = atan2(+q(1)*q(3)+q(2)*q(4), q(1)*q(2)-q(3)*q(4))
@ -1681,7 +1681,7 @@ pure function math_qToAxisAngle(Q)
real(pReal) :: halfAngle, sinHalfAngle
real(pReal), dimension(4) :: math_qToAxisAngle
halfAngle = acos(math_limit(Q(1),-1.0_pReal,1.0_pReal))
halfAngle = acos(math_clip(Q(1),-1.0_pReal,1.0_pReal))
sinHalfAngle = sin(halfAngle)
smallRotation: if (sinHalfAngle <= 1.0e-4_pReal) then
@ -1741,7 +1741,7 @@ real(pReal) pure function math_EulerMisorientation(EulerA,EulerB)
cosTheta = (math_trace33(math_mul33x33(math_EulerToR(EulerB), &
transpose(math_EulerToR(EulerA)))) - 1.0_pReal) * 0.5_pReal
math_EulerMisorientation = acos(math_limit(cosTheta,-1.0_pReal,1.0_pReal))
math_EulerMisorientation = acos(math_clip(cosTheta,-1.0_pReal,1.0_pReal))
end function math_EulerMisorientation
@ -2052,7 +2052,7 @@ function math_eigenvectorBasisSym33(m)
EB(3,3,3)=1.0_pReal
else threeSimilarEigenvalues
rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
phi=acos(math_limit(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
values = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* &
[cos(phi/3.0_pReal), &
cos((phi+2.0_pReal*PI)/3.0_pReal), &
@ -2117,7 +2117,7 @@ function math_eigenvectorBasisSym33_log(m)
EB(3,3,3)=1.0_pReal
else threeSimilarEigenvalues
rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
phi=acos(math_limit(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
values = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* &
[cos(phi/3.0_pReal), &
cos((phi+2.0_pReal*PI)/3.0_pReal), &
@ -2229,7 +2229,7 @@ function math_eigenvaluesSym33(m)
math_eigenvaluesSym33 = math_eigenvaluesSym(m)
else
rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
phi=acos(math_limit(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
math_eigenvaluesSym33 = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* &
[cos(phi/3.0_pReal), &
cos((phi+2.0_pReal*PI)/3.0_pReal), &
@ -2614,7 +2614,7 @@ end function math_rotate_forward3333
!> @brief limits a scalar value to a certain range (either one or two sided)
! Will return NaN if left > right
!--------------------------------------------------------------------------------------------------
real(pReal) pure function math_limit(a, left, right)
real(pReal) pure function math_clip(a, left, right)
use, intrinsic :: &
IEEE_arithmetic
@ -2623,14 +2623,14 @@ real(pReal) pure function math_limit(a, left, right)
real(pReal), intent(in), optional :: left, right
math_limit = min ( &
math_clip = min ( &
max (merge(left, -huge(a), present(left)), a), &
merge(right, huge(a), present(right)) &
)
if (present(left) .and. present(right)) &
math_limit = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_limit, left>right)
math_clip = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_clip, left>right)
end function math_limit
end function math_clip
end module math