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

Merge remote-tracking branch 'origin/development' into clean-and-polish-damage

This commit is contained in:
Martin Diehl 2020-03-10 13:49:11 +01:00
parent aa75591ea9 0dea7a3c7e
commit 705ee908a2
34 changed files with 3751 additions and 4959 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
.gitattributes vendored
View File

@ -8,3 +8,10 @@
*.jpg binary *.jpg binary
*.hdf5 binary *.hdf5 binary
*.pdf binary *.pdf binary
# ignore files from MSC.Marc in language statistics
installation/mods_MarcMentat/* linguist-vendored
src/MarcInclude/* linguist-vendored
# ignore reference files for tests in language statistics
python/tests/reference/* linguist-vendored

7
.gitlab-ci.yml
View File

@ -142,13 +142,6 @@ Pre_General:
- master - master
- release - release
grid_geometryPacking:
stage: preprocessing
script: grid_geometryPacking/test.py
except:
- master
- release
################################################################################################### ###################################################################################################
Post_AverageDown: Post_AverageDown:
stage: postprocessing stage: postprocessing

2
VERSION
View File

@ -1 +1 @@
v2.0.3-1747-ga2e8e5b1 v2.0.3-1862-g0b340a6d

8
processing/post/DADF5_postResults.py
View File

@ -52,15 +52,15 @@ for filename in options.filenames:
table = damask.Table(np.ones(np.product(results.grid),dtype=int)*int(inc[3:]),{'inc':(1,)}) table = damask.Table(np.ones(np.product(results.grid),dtype=int)*int(inc[3:]),{'inc':(1,)})
table.add('pos',coords.reshape((-1,3))) table.add('pos',coords.reshape((-1,3)))
results.set_visible('materialpoints',False) results.pick('materialpoints',False)
results.set_visible('constituents', True) results.pick('constituents', True)
for label in options.con: for label in options.con:
x = results.get_dataset_location(label) x = results.get_dataset_location(label)
if len(x) != 0: if len(x) != 0:
table.add(label,results.read_dataset(x,0,plain=True).reshape((results.grid.prod(),-1))) table.add(label,results.read_dataset(x,0,plain=True).reshape((results.grid.prod(),-1)))
results.set_visible('constituents', False) results.pick('constituents', False)
results.set_visible('materialpoints',True) results.pick('materialpoints',True)
for label in options.mat: for label in options.mat:
x = results.get_dataset_location(label) x = results.get_dataset_location(label)
if len(x) != 0: if len(x) != 0:

1139
processing/post/postResults.py
View File

File diff suppressed because it is too large Load Diff

21
processing/pre/seeds_fromDistribution.py
View File

@ -221,13 +221,9 @@ baseFile=os.path.splitext(os.path.basename(options.seedFile))[0]
points = np.array(options.grid).prod().astype('float') points = np.array(options.grid).prod().astype('float')
# ----------- calculate target distribution and bin edges # ----------- calculate target distribution and bin edges
targetGeomFile = os.path.splitext(os.path.basename(options.target))[0]+'.geom' targetGeom = damask.Geom.from_file(os.path.splitext(os.path.basename(options.target))[0]+'.geom')
targetGeomTable = damask.ASCIItable(targetGeomFile,None,labeled=False,readonly=True) nMicrostructures = len(np.unique(targetGeom.microstructure))
targetGeomTable.head_read() targetVolFrac = np.bincount(targetGeom.microstructure.flatten())/targetGeom.grid.prod().astype(np.float)
info,devNull = targetGeomTable.head_getGeom()
nMicrostructures = info['microstructures']
targetVolFrac = np.bincount(targetGeomTable.microstructure_read(info['grid']))[1:nMicrostructures+1]/\
float(info['grid'].prod())
target=[] target=[]
for i in range(1,nMicrostructures+1): for i in range(1,nMicrostructures+1):
targetHist,targetBins = np.histogram(targetVolFrac,bins=i) #bin boundaries targetHist,targetBins = np.histogram(targetVolFrac,bins=i) #bin boundaries
@ -251,13 +247,12 @@ initialGeomVFile = StringIO()
initialGeomVFile.write(damask.util.execute('geom_fromVoronoiTessellation '+ initialGeomVFile.write(damask.util.execute('geom_fromVoronoiTessellation '+
' -g '+' '.join(list(map(str, options.grid))),bestSeedsVFile)[0]) ' -g '+' '.join(list(map(str, options.grid))),bestSeedsVFile)[0])
initialGeomVFile.seek(0) initialGeomVFile.seek(0)
initialGeomTable = damask.ASCIItable(initialGeomVFile,None,labeled=False,readonly=True) initialGeom = damask.Geom.from_file(initialGeomVFile)
initialGeomTable.head_read()
info,devNull = initialGeomTable.head_getGeom()
if info['microstructures'] != nMicrostructures: damask.util.croak('error. Microstructure count mismatch') if len(np.unique(targetGeom.microstructure)) != nMicrostructures:
damask.util.croak('error. Microstructure count mismatch')
initialData = np.bincount(initialGeomTable.microstructure_read(info['grid']))/points initialData = np.bincount(initialGeom.microstructure.flatten())/points
for i in range(nMicrostructures): for i in range(nMicrostructures):
initialHist = np.histogram(initialData,bins=target[i]['bins'])[0] initialHist = np.histogram(initialData,bins=target[i]['bins'])[0]
target[i]['error']=np.sqrt(np.square(np.array(target[i]['histogram']-initialHist)).sum()) target[i]['error']=np.sqrt(np.square(np.array(target[i]['histogram']-initialHist)).sum())
@ -273,7 +268,7 @@ for i in range(nMicrostructures):
if options.maxseeds < 1: if options.maxseeds < 1:
maxSeeds = info['microstructures'] maxSeeds = len(np.unique(initialGeom.microstructure))
else: else:
maxSeeds = options.maxseeds maxSeeds = options.maxseeds

13
python/README
View File

@ -1,13 +0,0 @@
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

1
python/README Symbolic link
View File

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

148
python/damask/Lambert.py
View File

@ -6,27 +6,27 @@
# Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University # Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University
# All rights reserved. # All rights reserved.
# #
# Redistribution and use in source and binary forms, with or without modification, are # Redistribution and use in source and binary forms, with or without modification, are
# permitted provided that the following conditions are met: # permitted provided that the following conditions are met:
# #
# - Redistributions of source code must retain the above copyright notice, this list # - Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer. # of conditions and the following disclaimer.
# - Redistributions in binary form must reproduce the above copyright notice, this # - Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or # list of conditions and the following disclaimer in the documentation and/or
# other materials provided with the distribution. # other materials provided with the distribution.
# - Neither the names of Marc De Graef, Carnegie Mellon University nor the names # - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
# of its contributors may be used to endorse or promote products derived from # of its contributors may be used to endorse or promote products derived from
# this software without specific prior written permission. # this software without specific prior written permission.
# #
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE # ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL # LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR # DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER # SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, # CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE # OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
# USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. # USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#################################################################################################### ####################################################################################################
@ -44,7 +44,7 @@ def CubeToBall(cube):
---------- ----------
cube : numpy.ndarray cube : numpy.ndarray
coordinates of a point in a uniform refinable cubical grid. coordinates of a point in a uniform refinable cubical grid.
References References
---------- ----------
D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014 D. Roşca et al., Modelling and Simulation in Materials Science and Engineering 22:075013, 2014
@ -52,37 +52,37 @@ def CubeToBall(cube):
""" """
if np.abs(np.max(cube))>np.pi**(2./3.) * 0.5: if np.abs(np.max(cube))>np.pi**(2./3.) * 0.5:
raise ValueError raise ValueError
# transform to the sphere grid via the curved square, and intercept the zero point # transform to the sphere grid via the curved square, and intercept the zero point
if np.allclose(cube,0.0,rtol=0.0,atol=1.0e-300): if np.allclose(cube,0.0,rtol=0.0,atol=1.0e-300):
ball = np.zeros(3) ball = np.zeros(3)
else: else:
# get pyramide and scale by grid parameter ratio # get pyramide and scale by grid parameter ratio
p = get_order(cube) p = get_order(cube)
XYZ = cube[p] * sc XYZ = cube[p] * sc
# intercept all the points along the z-axis # intercept all the points along the z-axis
if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-300): if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-300):
ball = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]]) ball = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]])
else: else:
order = [1,0] if np.abs(XYZ[1]) <= np.abs(XYZ[0]) else [0,1] order = [1,0] if np.abs(XYZ[1]) <= np.abs(XYZ[0]) else [0,1]
q = np.pi/12.0 * XYZ[order[0]]/XYZ[order[1]] q = np.pi/12.0 * XYZ[order[0]]/XYZ[order[1]]
c = np.cos(q) c = np.cos(q)
s = np.sin(q) s = np.sin(q)
q = R1*2.0**0.25/beta * XYZ[order[1]] / np.sqrt(np.sqrt(2.0)-c) q = R1*2.0**0.25/beta * XYZ[order[1]] / np.sqrt(np.sqrt(2.0)-c)
T = np.array([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q T = np.array([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q
# transform to sphere grid (inverse Lambert) # transform to sphere grid (inverse Lambert)
# note that there is no need to worry about dividing by zero, since XYZ[2] can not become zero # note that there is no need to worry about dividing by zero, since XYZ[2] can not become zero
c = np.sum(T**2) c = np.sum(T**2)
s = c * np.pi/24.0 /XYZ[2]**2 s = c * np.pi/24.0 /XYZ[2]**2
c = c * np.sqrt(np.pi/24.0)/XYZ[2] c = c * np.sqrt(np.pi/24.0)/XYZ[2]
q = np.sqrt( 1.0 - s ) q = np.sqrt( 1.0 - s )
ball = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ]) ball = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ])
# reverse the coordinates back to the regular order according to the original pyramid number # reverse the coordinates back to the regular order according to the original pyramid number
ball = ball[p] ball = ball[p]
return ball return ball
@ -103,46 +103,46 @@ def BallToCube(ball):
""" """
rs = np.linalg.norm(ball) rs = np.linalg.norm(ball)
if rs > R1: if rs > R1:
raise ValueError raise ValueError
if np.allclose(ball,0.0,rtol=0.0,atol=1.0e-300): if np.allclose(ball,0.0,rtol=0.0,atol=1.0e-300):
cube = np.zeros(3) cube = np.zeros(3)
else: else:
p = get_order(ball) p = get_order(ball)
xyz3 = ball[p] xyz3 = ball[p]
# inverse M_3 # inverse M_3
xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) ) xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) )
# inverse M_2 # inverse M_2
qxy = np.sum(xyz2**2) qxy = np.sum(xyz2**2)
if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-300): if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-300):
Tinv = np.zeros(2) Tinv = np.zeros(2)
else: else:
q2 = qxy + np.max(np.abs(xyz2))**2 q2 = qxy + np.max(np.abs(xyz2))**2
sq2 = np.sqrt(q2) sq2 = np.sqrt(q2)
q = (beta/np.sqrt(2.0)/R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2))*sq2)) q = (beta/np.sqrt(2.0)/R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2))*sq2))
tt = np.clip((np.min(np.abs(xyz2))**2+np.max(np.abs(xyz2))*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0) tt = np.clip((np.min(np.abs(xyz2))**2+np.max(np.abs(xyz2))*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0)
Tinv = np.array([1.0,np.arccos(tt)/np.pi*12.0]) if np.abs(xyz2[1]) <= np.abs(xyz2[0]) else \ Tinv = np.array([1.0,np.arccos(tt)/np.pi*12.0]) if np.abs(xyz2[1]) <= np.abs(xyz2[0]) else \
np.array([np.arccos(tt)/np.pi*12.0,1.0]) np.array([np.arccos(tt)/np.pi*12.0,1.0])
Tinv = q * np.where(xyz2<0.0,-Tinv,Tinv) Tinv = q * np.where(xyz2<0.0,-Tinv,Tinv)
# inverse M_1 # inverse M_1
cube = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /sc cube = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /sc
# reverse the coordinates back to the regular order according to the original pyramid number
cube = cube[p]
# reverse the coordinates back to the regular order according to the original pyramid number
cube = cube[p]
return cube return cube
def get_order(xyz): def get_order(xyz):
""" """
Get order of the coordinates. Get order of the coordinates.
Depending on the pyramid in which the point is located, the order need to be adjusted. Depending on the pyramid in which the point is located, the order need to be adjusted.
Parameters Parameters
---------- ----------
xyz : numpy.ndarray xyz : numpy.ndarray
@ -157,10 +157,10 @@ def get_order(xyz):
""" """
if (abs(xyz[0])<= xyz[2]) and (abs(xyz[1])<= xyz[2]) or \ if (abs(xyz[0])<= xyz[2]) and (abs(xyz[1])<= xyz[2]) or \
(abs(xyz[0])<=-xyz[2]) and (abs(xyz[1])<=-xyz[2]): (abs(xyz[0])<=-xyz[2]) and (abs(xyz[1])<=-xyz[2]):
return [0,1,2] return [0,1,2]
elif (abs(xyz[2])<= xyz[0]) and (abs(xyz[1])<= xyz[0]) or \ elif (abs(xyz[2])<= xyz[0]) and (abs(xyz[1])<= xyz[0]) or \
(abs(xyz[2])<=-xyz[0]) and (abs(xyz[1])<=-xyz[0]): (abs(xyz[2])<=-xyz[0]) and (abs(xyz[1])<=-xyz[0]):
return [1,2,0] return [1,2,0]
elif (abs(xyz[0])<= xyz[1]) and (abs(xyz[2])<= xyz[1]) or \ elif (abs(xyz[0])<= xyz[1]) and (abs(xyz[2])<= xyz[1]) or \
(abs(xyz[0])<=-xyz[1]) and (abs(xyz[2])<=-xyz[1]): (abs(xyz[0])<=-xyz[1]) and (abs(xyz[2])<=-xyz[1]):
return [2,0,1] return [2,0,1]

8
python/damask/__init__.py
View File

@ -13,9 +13,11 @@ 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 .rotation import Rotation # noqa
from .dadf5 import DADF5 # noqa from .lattice import Symmetry, Lattice# noqa
from .dadf5 import DADF5 as Result # noqa from .orientation import Orientation # noqa
from .result import Result # noqa
from .result import Result as DADF5 # noqa
from .geom import Geom # noqa from .geom import Geom # noqa
from .solver import Solver # noqa from .solver import Solver # noqa

69
python/damask/asciitable.py
View File

@ -11,12 +11,12 @@ class ASCIItable():
"""Read and write to ASCII tables.""" """Read and write to ASCII tables."""
tmpext = '_tmp' # filename extension for in-place access tmpext = '_tmp' # filename extension for in-place access
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def __init__(self, def __init__(self,
name = None, name = None,
outname = None, outname = None,
buffered = False, # flush writes buffered = False, # is ignored, only exists for compatibility reasons
labeled = True, # assume table has labels labeled = True, # assume table has labels
readonly = False, # no reading from file readonly = False, # no reading from file
): ):
@ -52,7 +52,7 @@ class ASCIItable():
if self.__IO__['in'] is None \ if self.__IO__['in'] is None \
or self.__IO__['out'] is None: raise IOError # complain if any required file access not possible or self.__IO__['out'] is None: raise IOError # complain if any required file access not possible
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def _removeCRLF(self, def _removeCRLF(self,
@ -63,7 +63,6 @@ class ASCIItable():
except AttributeError: except AttributeError:
return str(string) return str(string)
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def _quote(self, def _quote(self,
what): what):
@ -71,6 +70,7 @@ class ASCIItable():
return '{quote}{content}{quote}'.format( return '{quote}{content}{quote}'.format(
quote = ('"' if str(what)=='' or re.search(r"\s",str(what)) else ''), quote = ('"' if str(what)=='' or re.search(r"\s",str(what)) else ''),
content = what) content = what)
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def close(self, def close(self,
dismiss = False): dismiss = False):
@ -110,7 +110,7 @@ class ASCIItable():
def head_read(self): def head_read(self):
""" """
Get column labels. Get column labels.
by either reading the first row or, by either reading the first row or,
if keyword "head[*]" is present, the last line of the header if keyword "head[*]" is present, the last line of the header
""" """
@ -121,7 +121,7 @@ class ASCIItable():
firstline = self.__IO__['in'].readline().strip() firstline = self.__IO__['in'].readline().strip()
m = re.search(r'(\d+)\s+head', firstline.lower()) # search for "head" keyword m = re.search(r'(\d+)\s+head', firstline.lower()) # search for "head" keyword
if m: # proper ASCIItable format if m: # proper ASCIItable format
if self.__IO__['labeled']: # table features labels if self.__IO__['labeled']: # table features labels
@ -145,7 +145,7 @@ class ASCIItable():
if self.__IO__['labeled']: # table features labels if self.__IO__['labeled']: # table features labels
self.tags = self.data # get tags from last line in "header"... self.tags = self.data # get tags from last line in "header"...
self.data_read() # ...and remove from buffer self.data_read() # ...and remove from buffer
if self.__IO__['labeled']: # table features tags if self.__IO__['labeled']: # table features tags
self.__IO__['tags'] = list(self.tags) # backup tags (make COPY, not link) self.__IO__['tags'] = list(self.tags) # backup tags (make COPY, not link)
@ -163,7 +163,7 @@ class ASCIItable():
if self.__IO__['labeled']: if self.__IO__['labeled']:
head.append('\t'.join(map(self._quote,self.tags))) head.append('\t'.join(map(self._quote,self.tags)))
if len(self.tags) == 0: raise ValueError('no labels present.') if len(self.tags) == 0: raise ValueError('no labels present.')
return self.output_write(head) return self.output_write(head)
# ------------------------------------------------------------------ # ------------------------------------------------------------------
@ -178,15 +178,11 @@ class ASCIItable():
'grid': lambda x: int(x), 'grid': lambda x: int(x),
'size': lambda x: float(x), 'size': lambda x: float(x),
'origin': lambda x: float(x), 'origin': lambda x: float(x),
'homogenization': lambda x: int(x),
'microstructures': lambda x: int(x),
} }
info = { info = {
'grid': np.zeros(3,'i'), 'grid': np.zeros(3,'i'),
'size': np.zeros(3,'d'), 'size': np.zeros(3,'d'),
'origin': np.zeros(3,'d'), 'origin': np.zeros(3,'d'),
'homogenization': 0,
'microstructures': 0,
} }
extra_header = [] extra_header = []
@ -234,7 +230,7 @@ class ASCIItable():
raw = False): raw = False):
""" """
Tell abstract labels. Tell abstract labels.
"x" for "1_x","2_x",... unless raw output is requested. "x" for "1_x","2_x",... unless raw output is requested.
operates on object tags or given list. operates on object tags or given list.
""" """
@ -347,7 +343,7 @@ class ASCIItable():
""" """
start = self.label_index(labels) start = self.label_index(labels)
dim = self.label_dimension(labels) dim = self.label_dimension(labels)
return np.hstack([range(s,s+d) for s,d in zip(start,dim)]).astype(int) \ return np.hstack([range(s,s+d) for s,d in zip(start,dim)]).astype(int) \
if isinstance(labels, Iterable) and not isinstance(labels, str) \ if isinstance(labels, Iterable) and not isinstance(labels, str) \
else range(start,start+dim) else range(start,start+dim)
@ -375,15 +371,6 @@ class ASCIItable():
self.tags = list(self.__IO__['tags']) # restore label info found in header (as COPY, not link) self.tags = list(self.__IO__['tags']) # restore label info found in header (as COPY, not link)
self.__IO__['labeled'] = len(self.tags) > 0 self.__IO__['labeled'] = len(self.tags) > 0
# ------------------------------------------------------------------
def data_skipLines(self,
count):
"""Wind forward by count number of lines."""
for i in range(count):
alive = self.data_read()
return alive
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def data_read(self, def data_read(self,
advance = True, advance = True,
@ -425,8 +412,8 @@ class ASCIItable():
columns = [] columns = []
for i,(c,d) in enumerate(zip(indices[present],dimensions[present])): # for all valid labels ... for i,(c,d) in enumerate(zip(indices[present],dimensions[present])): # for all valid labels ...
# ... transparently add all components unless column referenced by number or with explicit dimension # ... transparently add all components unless column referenced by number or with explicit dimension
columns += list(range(c,c + columns += list(range(c,c +
(d if str(c) != str(labels[present[i]]) else (d if str(c) != str(labels[present[i]]) else
1))) 1)))
use = np.array(columns) if len(columns) > 0 else None use = np.array(columns) if len(columns) > 0 else None
@ -457,7 +444,7 @@ class ASCIItable():
output = [fmt % value for value in row] if fmt else list(map(repr,row)) output = [fmt % value for value in row] if fmt else list(map(repr,row))
except Exception: except Exception:
output = [fmt % row] if fmt else [repr(row)] output = [fmt % row] if fmt else [repr(row)]
try: try:
self.__IO__['out'].write(delimiter.join(output) + '\n') self.__IO__['out'].write(delimiter.join(output) + '\n')
except Exception: except Exception:
@ -473,33 +460,3 @@ class ASCIItable():
for item in what: self.data_append(item) for item in what: self.data_append(item)
except TypeError: except TypeError:
self.data += [str(what)] self.data += [str(what)]
# ------------------------------------------------------------------
def microstructure_read(self,
grid,
type = 'i',
strict = False):
"""Read microstructure data (from .geom format)."""
def datatype(item):
return int(item) if type.lower() == 'i' else float(item)
N = grid.prod() # expected number of microstructure indices in data
microstructure = np.zeros(N,type) # initialize as flat array
i = 0
while i < N and self.data_read():
items = self.data
if len(items) > 2:
if items[1].lower() == 'of':
items = np.ones(datatype(items[0]))*datatype(items[2])
elif items[1].lower() == 'to':
items = np.linspace(datatype(items[0]),datatype(items[2]),
abs(datatype(items[2])-datatype(items[0]))+1,dtype=int)
else: items = list(map(datatype,items))
else: items = list(map(datatype,items))
s = min(len(items), N-i) # prevent overflow of microstructure array
microstructure[i:i+s] = items[:s]
i += len(items)
return (microstructure, i == N and not self.data_read()) if strict else microstructure # check for proper point count and end of file

558
python/damask/colormaps.py
View File

@ -1,359 +1,355 @@
import numpy as np import numpy as np
class Color(): class Color():
"""Color representation in and conversion between different color-spaces.""" """Color representation in and conversion between different color-spaces."""
__slots__ = [ __slots__ = [
'model', 'model',
'color', 'color',
'__dict__', '__dict__',
] ]
# ------------------------------------------------------------------ def __init__(self,
def __init__(self, model = 'RGB',
model = 'RGB', color = np.zeros(3,'d')):
color = np.zeros(3,'d')): """
""" Create a Color object.
Create a Color object.
Parameters
----------
model : string
color model
color : numpy.ndarray
vector representing the color according to the selected model
""" Parameters
self.__transforms__ = \ ----------
{'HSV': {'index': 0, 'next': self._HSV2HSL}, model : string
'HSL': {'index': 1, 'next': self._HSL2RGB, 'prev': self._HSL2HSV}, color model
'RGB': {'index': 2, 'next': self._RGB2XYZ, 'prev': self._RGB2HSL}, color : numpy.ndarray
'XYZ': {'index': 3, 'next': self._XYZ2CIELAB, 'prev': self._XYZ2RGB}, vector representing the color according to the selected model
'CIELAB': {'index': 4, 'next': self._CIELAB2MSH, 'prev': self._CIELAB2XYZ},
'MSH': {'index': 5, 'prev': self._MSH2CIELAB},
}
model = model.upper() """
if model not in list(self.__transforms__.keys()): model = 'RGB' self.__transforms__ = \
if model == 'RGB' and max(color) > 1.0: # are we RGB255 ? {'HSV': {'index': 0, 'next': self._HSV2HSL},
for i in range(3): 'HSL': {'index': 1, 'next': self._HSL2RGB, 'prev': self._HSL2HSV},
color[i] /= 255.0 # rescale to RGB 'RGB': {'index': 2, 'next': self._RGB2XYZ, 'prev': self._RGB2HSL},
'XYZ': {'index': 3, 'next': self._XYZ2CIELAB, 'prev': self._XYZ2RGB},
'CIELAB': {'index': 4, 'next': self._CIELAB2MSH, 'prev': self._CIELAB2XYZ},
'MSH': {'index': 5, 'prev': self._MSH2CIELAB},
}
if model == 'HSL': # are we HSL ? model = model.upper()
if abs(color[0]) > 1.0: color[0] /= 360.0 # with angular hue? if model not in list(self.__transforms__.keys()): model = 'RGB'
while color[0] >= 1.0: color[0] -= 1.0 # rewind to proper range if model == 'RGB' and max(color) > 1.0: # are we RGB255 ?
while color[0] < 0.0: color[0] += 1.0 # rewind to proper range for i in range(3):
color[i] /= 255.0 # rescale to RGB
self.model = model if model == 'HSL': # are we HSL ?
self.color = np.array(color,'d') if abs(color[0]) > 1.0: color[0] /= 360.0 # with angular hue?
while color[0] >= 1.0: color[0] -= 1.0 # rewind to proper range
while color[0] < 0.0: color[0] += 1.0 # rewind to proper range
self.model = model
self.color = np.array(color,'d')
# ------------------------------------------------------------------ def __repr__(self):
def __repr__(self): """Color model and values."""
"""Color model and values.""" return 'Model: %s Color: %s'%(self.model,str(self.color))
return 'Model: %s Color: %s'%(self.model,str(self.color))
# ------------------------------------------------------------------ def __str__(self):
def __str__(self): """Color model and values."""
"""Color model and values.""" return self.__repr__()
return self.__repr__()
# ------------------------------------------------------------------ def convert_to(self,toModel = 'RGB'):
def convertTo(self,toModel = 'RGB'): """
""" Change the color model permanently.
Change the color model permanently.
Parameters
----------
toModel : string
color model
""" Parameters
toModel = toModel.upper() ----------
if toModel not in list(self.__transforms__.keys()): return toModel : string
color model
sourcePos = self.__transforms__[self.model]['index'] """
targetPos = self.__transforms__[toModel]['index'] toModel = toModel.upper()
if toModel not in list(self.__transforms__.keys()): return
while sourcePos < targetPos: sourcePos = self.__transforms__[self.model]['index']
self.__transforms__[self.model]['next']() targetPos = self.__transforms__[toModel]['index']
sourcePos += 1
while sourcePos > targetPos: while sourcePos < targetPos:
self.__transforms__[self.model]['prev']() self.__transforms__[self.model]['next']()
sourcePos -= 1 sourcePos += 1
return self
while sourcePos > targetPos:
self.__transforms__[self.model]['prev']()
sourcePos -= 1
return self
# ------------------------------------------------------------------ def express_as(self,asModel = 'RGB'):
def expressAs(self,asModel = 'RGB'): """
""" Return the color in a different model.
Return the color in a different model.
Parameters
----------
asModel : string
color model
""" Parameters
return self.__class__(self.model,self.color).convertTo(asModel) ----------
asModel : string
color model
"""
return self.__class__(self.model,self.color).convert_to(asModel)
def _HSV2HSL(self): def _HSV2HSL(self):
""" """
Convert H(ue) S(aturation) V(alue or brightness) to H(ue) S(aturation) L(uminance). Convert H(ue) S(aturation) V(alue or brightness) to H(ue) S(aturation) L(uminance).
All values are in the range [0,1] All values are in the range [0,1]
http://codeitdown.com/hsl-hsb-hsv-color http://codeitdown.com/hsl-hsb-hsv-color
""" """
if self.model != 'HSV': return if self.model != 'HSV': return
converted = Color('HSL',np.array([ converted = Color('HSL',np.array([
self.color[0], self.color[0],
1. if self.color[2] == 0.0 or (self.color[1] == 0.0 and self.color[2] == 1.0) \ 1. if self.color[2] == 0.0 or (self.color[1] == 0.0 and self.color[2] == 1.0) \
else self.color[1]*self.color[2]/(1.-abs(self.color[2]*(2.-self.color[1])-1.)), else self.color[1]*self.color[2]/(1.-abs(self.color[2]*(2.-self.color[1])-1.)),
0.5*self.color[2]*(2.-self.color[1]), 0.5*self.color[2]*(2.-self.color[1]),
])) ]))
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _HSL2HSV(self): def _HSL2HSV(self):
""" """
Convert H(ue) S(aturation) L(uminance) to H(ue) S(aturation) V(alue or brightness). Convert H(ue) S(aturation) L(uminance) to H(ue) S(aturation) V(alue or brightness).
All values are in the range [0,1] All values are in the range [0,1]
http://codeitdown.com/hsl-hsb-hsv-color http://codeitdown.com/hsl-hsb-hsv-color
""" """
if self.model != 'HSL': return if self.model != 'HSL': return
h = self.color[0] h = self.color[0]
b = self.color[2]+0.5*(self.color[1]*(1.-abs(2*self.color[2]-1))) b = self.color[2]+0.5*(self.color[1]*(1.-abs(2*self.color[2]-1)))
s = 1.0 if b == 0.0 else 2.*(b-self.color[2])/b s = 1.0 if b == 0.0 else 2.*(b-self.color[2])/b
converted = Color('HSV',np.array([h,s,b])) converted = Color('HSV',np.array([h,s,b]))
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _HSL2RGB(self): def _HSL2RGB(self):
""" """
Convert H(ue) S(aturation) L(uminance) to R(red) G(reen) B(lue). Convert H(ue) S(aturation) L(uminance) to R(red) G(reen) B(lue).
All values are in the range [0,1] All values are in the range [0,1]
from http://en.wikipedia.org/wiki/HSL_and_HSV from http://en.wikipedia.org/wiki/HSL_and_HSV
""" """
if self.model != 'HSL': return if self.model != 'HSL': return
sextant = self.color[0]*6.0 sextant = self.color[0]*6.0
c = (1.0 - abs(2.0 * self.color[2] - 1.0))*self.color[1] c = (1.0 - abs(2.0 * self.color[2] - 1.0))*self.color[1]
x = c*(1.0 - abs(sextant%2 - 1.0)) x = c*(1.0 - abs(sextant%2 - 1.0))
m = self.color[2] - 0.5*c m = self.color[2] - 0.5*c
converted = Color('RGB',np.array([ converted = Color('RGB',np.array([
[c+m, x+m, m], [c+m, x+m, m],
[x+m, c+m, m], [x+m, c+m, m],
[m, c+m, x+m], [m, c+m, x+m],
[m, x+m, c+m], [m, x+m, c+m],
[x+m, m, c+m], [x+m, m, c+m],
[c+m, m, x+m], [c+m, m, x+m],
][int(sextant)],'d')) ][int(sextant)],'d'))
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _RGB2HSL(self): def _RGB2HSL(self):
""" """
Convert R(ed) G(reen) B(lue) to H(ue) S(aturation) L(uminance). Convert R(ed) G(reen) B(lue) to H(ue) S(aturation) L(uminance).
All values are in the range [0,1] All values are in the range [0,1]
from http://130.113.54.154/~monger/hsl-rgb.html from http://130.113.54.154/~monger/hsl-rgb.html
""" """
if self.model != 'RGB': return if self.model != 'RGB': return
HSL = np.zeros(3,'d') HSL = np.zeros(3,'d')
maxcolor = self.color.max() maxcolor = self.color.max()
mincolor = self.color.min() mincolor = self.color.min()
HSL[2] = (maxcolor + mincolor)/2.0 HSL[2] = (maxcolor + mincolor)/2.0
if(mincolor == maxcolor): if(mincolor == maxcolor):
HSL[0] = 0.0 HSL[0] = 0.0
HSL[1] = 0.0 HSL[1] = 0.0
else:
if (HSL[2]<0.5):
HSL[1] = (maxcolor - mincolor)/(maxcolor + mincolor)
else: else:
HSL[1] = (maxcolor - mincolor)/(2.0 - maxcolor - mincolor) if (HSL[2]<0.5):
if (maxcolor == self.color[0]): HSL[1] = (maxcolor - mincolor)/(maxcolor + mincolor)
HSL[0] = 0.0 + (self.color[1] - self.color[2])/(maxcolor - mincolor) else:
elif (maxcolor == self.color[1]): HSL[1] = (maxcolor - mincolor)/(2.0 - maxcolor - mincolor)
HSL[0] = 2.0 + (self.color[2] - self.color[0])/(maxcolor - mincolor) if (maxcolor == self.color[0]):
elif (maxcolor == self.color[2]): HSL[0] = 0.0 + (self.color[1] - self.color[2])/(maxcolor - mincolor)
HSL[0] = 4.0 + (self.color[0] - self.color[1])/(maxcolor - mincolor) elif (maxcolor == self.color[1]):
HSL[0] = HSL[0]*60.0 # scaling to 360 might be dangerous for small values HSL[0] = 2.0 + (self.color[2] - self.color[0])/(maxcolor - mincolor)
if (HSL[0] < 0.0): elif (maxcolor == self.color[2]):
HSL[0] = HSL[0] + 360.0 HSL[0] = 4.0 + (self.color[0] - self.color[1])/(maxcolor - mincolor)
for i in range(2): HSL[0] = HSL[0]*60.0 # scaling to 360 might be dangerous for small values
HSL[i+1] = min(HSL[i+1],1.0) if (HSL[0] < 0.0):
HSL[i+1] = max(HSL[i+1],0.0) HSL[0] = HSL[0] + 360.0
for i in range(2):
HSL[i+1] = min(HSL[i+1],1.0)
HSL[i+1] = max(HSL[i+1],0.0)
converted = Color('HSL', HSL) converted = Color('HSL', HSL)
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _RGB2XYZ(self): def _RGB2XYZ(self):
""" """
Convert R(ed) G(reen) B(lue) to CIE XYZ. Convert R(ed) G(reen) B(lue) to CIE XYZ.
All values are in the range [0,1] All values are in the range [0,1]
from http://www.cs.rit.edu/~ncs/color/t_convert.html from http://www.cs.rit.edu/~ncs/color/t_convert.html
""" """
if self.model != 'RGB': return if self.model != 'RGB': return
XYZ = np.zeros(3,'d') XYZ = np.zeros(3,'d')
RGB_lin = np.zeros(3,'d') RGB_lin = np.zeros(3,'d')
convert = np.array([[0.412453,0.357580,0.180423], convert = np.array([[0.412453,0.357580,0.180423],
[0.212671,0.715160,0.072169], [0.212671,0.715160,0.072169],
[0.019334,0.119193,0.950227]]) [0.019334,0.119193,0.950227]])
for i in range(3): for i in range(3):
if (self.color[i] > 0.04045): RGB_lin[i] = ((self.color[i]+0.0555)/1.0555)**2.4 if (self.color[i] > 0.04045): RGB_lin[i] = ((self.color[i]+0.0555)/1.0555)**2.4
else: RGB_lin[i] = self.color[i] /12.92 else: RGB_lin[i] = self.color[i] /12.92
XYZ = np.dot(convert,RGB_lin) XYZ = np.dot(convert,RGB_lin)
for i in range(3): for i in range(3):
XYZ[i] = max(XYZ[i],0.0) XYZ[i] = max(XYZ[i],0.0)
converted = Color('XYZ', XYZ) converted = Color('XYZ', XYZ)
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _XYZ2RGB(self): def _XYZ2RGB(self):
""" """
Convert CIE XYZ to R(ed) G(reen) B(lue). Convert CIE XYZ to R(ed) G(reen) B(lue).
All values are in the range [0,1] All values are in the range [0,1]
from http://www.cs.rit.edu/~ncs/color/t_convert.html from http://www.cs.rit.edu/~ncs/color/t_convert.html
""" """
if self.model != 'XYZ': if self.model != 'XYZ':
return return
convert = np.array([[ 3.240479,-1.537150,-0.498535], convert = np.array([[ 3.240479,-1.537150,-0.498535],
[-0.969256, 1.875992, 0.041556], [-0.969256, 1.875992, 0.041556],
[ 0.055648,-0.204043, 1.057311]]) [ 0.055648,-0.204043, 1.057311]])
RGB_lin = np.dot(convert,self.color) RGB_lin = np.dot(convert,self.color)
RGB = np.zeros(3,'d') RGB = np.zeros(3,'d')
for i in range(3): for i in range(3):
if (RGB_lin[i] > 0.0031308): RGB[i] = ((RGB_lin[i])**(1.0/2.4))*1.0555-0.0555 if (RGB_lin[i] > 0.0031308): RGB[i] = ((RGB_lin[i])**(1.0/2.4))*1.0555-0.0555
else: RGB[i] = RGB_lin[i] *12.92 else: RGB[i] = RGB_lin[i] *12.92
for i in range(3): for i in range(3):
RGB[i] = min(RGB[i],1.0) RGB[i] = min(RGB[i],1.0)
RGB[i] = max(RGB[i],0.0) RGB[i] = max(RGB[i],0.0)
maxVal = max(RGB) # clipping colors according to the display gamut maxVal = max(RGB) # clipping colors according to the display gamut
if (maxVal > 1.0): RGB /= maxVal if (maxVal > 1.0): RGB /= maxVal
converted = Color('RGB', RGB) converted = Color('RGB', RGB)
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _CIELAB2XYZ(self): def _CIELAB2XYZ(self):
""" """
Convert CIE Lab to CIE XYZ. Convert CIE Lab to CIE XYZ.
All values are in the range [0,1] All values are in the range [0,1]
from http://www.easyrgb.com/index.php?X=MATH&H=07#text7 from http://www.easyrgb.com/index.php?X=MATH&H=07#text7
""" """
if self.model != 'CIELAB': return if self.model != 'CIELAB': return
ref_white = np.array([.95047, 1.00000, 1.08883],'d') # Observer = 2, Illuminant = D65 ref_white = np.array([.95047, 1.00000, 1.08883],'d') # Observer = 2, Illuminant = D65
XYZ = np.zeros(3,'d') XYZ = np.zeros(3,'d')
XYZ[1] = (self.color[0] + 16.0 ) / 116.0 XYZ[1] = (self.color[0] + 16.0 ) / 116.0
XYZ[0] = XYZ[1] + self.color[1]/ 500.0 XYZ[0] = XYZ[1] + self.color[1]/ 500.0
XYZ[2] = XYZ[1] - self.color[2]/ 200.0 XYZ[2] = XYZ[1] - self.color[2]/ 200.0
for i in range(len(XYZ)): for i in range(len(XYZ)):
if (XYZ[i] > 6./29. ): XYZ[i] = XYZ[i]**3. if (XYZ[i] > 6./29. ): XYZ[i] = XYZ[i]**3.
else: XYZ[i] = 108./841. * (XYZ[i] - 4./29.) else: XYZ[i] = 108./841. * (XYZ[i] - 4./29.)
converted = Color('XYZ', XYZ*ref_white) converted = Color('XYZ', XYZ*ref_white)
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _XYZ2CIELAB(self):
"""
Convert CIE XYZ to CIE Lab.
All values are in the range [0,1]
from http://en.wikipedia.org/wiki/Lab_color_space,
http://www.cs.rit.edu/~ncs/color/t_convert.html
"""
if self.model != 'XYZ': return
ref_white = np.array([.95047, 1.00000, 1.08883],'d') # Observer = 2, Illuminant = D65
XYZ = self.color/ref_white
for i in range(len(XYZ)):
if (XYZ[i] > 216./24389 ): XYZ[i] = XYZ[i]**(1.0/3.0)
else: XYZ[i] = (841./108. * XYZ[i]) + 16.0/116.0
converted = Color('CIELAB', np.array([ 116.0 * XYZ[1] - 16.0,
500.0 * (XYZ[0] - XYZ[1]),
200.0 * (XYZ[1] - XYZ[2]) ]))
self.model = converted.model
self.color = converted.color
def _CIELAB2MSH(self): def _XYZ2CIELAB(self):
""" """
Convert CIE Lab to Msh colorspace. Convert CIE XYZ to CIE Lab.
from http://www.cs.unm.edu/~kmorel/documents/ColorMaps/DivergingColorMapWorkshop.xls All values are in the range [0,1]
""" from http://en.wikipedia.org/wiki/Lab_color_space,
if self.model != 'CIELAB': return http://www.cs.rit.edu/~ncs/color/t_convert.html
"""
if self.model != 'XYZ': return
Msh = np.zeros(3,'d') ref_white = np.array([.95047, 1.00000, 1.08883],'d') # Observer = 2, Illuminant = D65
Msh[0] = np.sqrt(np.dot(self.color,self.color)) XYZ = self.color/ref_white
if (Msh[0] > 0.001):
Msh[1] = np.arccos(self.color[0]/Msh[0])
if (self.color[1] != 0.0):
Msh[2] = np.arctan2(self.color[2],self.color[1])
converted = Color('MSH', Msh) for i in range(len(XYZ)):
self.model = converted.model if (XYZ[i] > 216./24389 ): XYZ[i] = XYZ[i]**(1.0/3.0)
self.color = converted.color else: XYZ[i] = (841./108. * XYZ[i]) + 16.0/116.0
converted = Color('CIELAB', np.array([ 116.0 * XYZ[1] - 16.0,
500.0 * (XYZ[0] - XYZ[1]),
200.0 * (XYZ[1] - XYZ[2]) ]))
self.model = converted.model
self.color = converted.color
def _MSH2CIELAB(self): def _CIELAB2MSH(self):
""" """
Convert Msh colorspace to CIE Lab. Convert CIE Lab to Msh colorspace.
with s,h in radians from http://www.cs.unm.edu/~kmorel/documents/ColorMaps/DivergingColorMapWorkshop.xls
from http://www.cs.unm.edu/~kmorel/documents/ColorMaps/DivergingColorMapWorkshop.xls """
""" if self.model != 'CIELAB': return
if self.model != 'MSH': return
Lab = np.zeros(3,'d') Msh = np.zeros(3,'d')
Lab[0] = self.color[0] * np.cos(self.color[1]) Msh[0] = np.sqrt(np.dot(self.color,self.color))
Lab[1] = self.color[0] * np.sin(self.color[1]) * np.cos(self.color[2]) if (Msh[0] > 0.001):
Lab[2] = self.color[0] * np.sin(self.color[1]) * np.sin(self.color[2]) Msh[1] = np.arccos(self.color[0]/Msh[0])
if (self.color[1] != 0.0):
Msh[2] = np.arctan2(self.color[2],self.color[1])
converted = Color('CIELAB', Lab) converted = Color('MSH', Msh)
self.model = converted.model self.model = converted.model
self.color = converted.color self.color = converted.color
def _MSH2CIELAB(self):
"""
Convert Msh colorspace to CIE Lab.
with s,h in radians
from http://www.cs.unm.edu/~kmorel/documents/ColorMaps/DivergingColorMapWorkshop.xls
"""
if self.model != 'MSH': return
Lab = np.zeros(3,'d')
Lab[0] = self.color[0] * np.cos(self.color[1])
Lab[1] = self.color[0] * np.sin(self.color[1]) * np.cos(self.color[2])
Lab[2] = self.color[0] * np.sin(self.color[1]) * np.sin(self.color[2])
converted = Color('CIELAB', Lab)
self.model = converted.model
self.color = converted.color
class Colormap(): class Colormap():
@ -416,7 +412,7 @@ class Colormap():
): ):
""" """
Create a Colormap object. Create a Colormap object.
Parameters Parameters
---------- ----------
left : Color left : Color
@ -498,13 +494,13 @@ class Colormap():
def interpolate_linear(lo, hi, frac): def interpolate_linear(lo, hi, frac):
"""Linear interpolation between lo and hi color at given fraction; output in model of lo color.""" """Linear interpolation between lo and hi color at given fraction; output in model of lo color."""
interpolation = (1.0 - frac) * np.array(lo.color[:]) \ interpolation = (1.0 - frac) * np.array(lo.color[:]) \
+ frac * np.array(hi.expressAs(lo.model).color[:]) + frac * np.array(hi.express_as(lo.model).color[:])
return Color(lo.model,interpolation) return Color(lo.model,interpolation)
if self.interpolate == 'perceptualuniform': if self.interpolate == 'perceptualuniform':
return interpolate_Msh(self.left.expressAs('MSH').color, return interpolate_Msh(self.left.express_as('MSH').color,
self.right.expressAs('MSH').color,fraction) self.right.express_as('MSH').color,fraction)
elif self.interpolate == 'linear': elif self.interpolate == 'linear':
return interpolate_linear(self.left, return interpolate_linear(self.left,
self.right,fraction) self.right,fraction)
@ -528,7 +524,7 @@ class Colormap():
""" """
format = format.lower() # consistent comparison basis format = format.lower() # consistent comparison basis
frac = 0.5*(np.array(crop) + 1.0) # rescale crop range to fractions frac = 0.5*(np.array(crop) + 1.0) # rescale crop range to fractions
colors = [self.color(float(i)/(steps-1)*(frac[1]-frac[0])+frac[0]).expressAs(model).color for i in range(steps)] colors = [self.color(float(i)/(steps-1)*(frac[1]-frac[0])+frac[0]).express_as(model).color for i in range(steps)]
if format == 'paraview': if format == 'paraview':
colormap = ['[\n {{\n "ColorSpace": "RGB", "Name": "{}", "DefaultMap": true,\n "RGBPoints" : ['.format(name)] \ colormap = ['[\n {{\n "ColorSpace": "RGB", "Name": "{}", "DefaultMap": true,\n "RGBPoints" : ['.format(name)] \
+ [' {:4d},{:8.6f},{:8.6f},{:8.6f},'.format(i,color[0],color[1],color[2],) \ + [' {:4d},{:8.6f},{:8.6f},{:8.6f},'.format(i,color[0],color[1],color[2],) \

1175
python/damask/dadf5.py
View File

File diff suppressed because it is too large Load Diff

6
python/damask/geom.py
View File

@ -80,7 +80,7 @@ class Geom():
if size is not None: if size is not None:
self.set_size(size) self.set_size(size)
elif rescale: elif rescale:
self.set_size(self.get_grid()/grid_old*self.size) self.set_size(self.get_grid()/grid_old*self.size)
message = ['grid a b c: {}'.format(' x '.join(map(str,grid_old)))] message = ['grid a b c: {}'.format(' x '.join(map(str,grid_old)))]
if np.any(grid_old != self.get_grid()): if np.any(grid_old != self.get_grid()):
@ -269,7 +269,7 @@ class Geom():
comments = [] comments = []
for i,line in enumerate(content[:header_length]): for i,line in enumerate(content[:header_length]):
items = line.lower().strip().split() items = line.lower().strip().split()
key = items[0] if len(items) > 0 else '' key = items[0] if items else ''
if key == 'grid': if key == 'grid':
grid = np.array([ int(dict(zip(items[1::2],items[2::2]))[i]) for i in ['a','b','c']]) grid = np.array([ int(dict(zip(items[1::2],items[2::2]))[i]) for i in ['a','b','c']])
elif key == 'size': elif key == 'size':
@ -524,7 +524,7 @@ class Geom():
"""Renumber sorted microstructure indices to 1,...,N.""" """Renumber sorted microstructure indices to 1,...,N."""
renumbered = np.empty(self.get_grid(),dtype=self.microstructure.dtype) renumbered = np.empty(self.get_grid(),dtype=self.microstructure.dtype)
for i, oldID in enumerate(np.unique(self.microstructure)): for i, oldID in enumerate(np.unique(self.microstructure)):
renumbered = np.where(self.microstructure == oldID, i+1, renumbered) renumbered = np.where(self.microstructure == oldID, i+1, renumbered)
return self.update(renumbered) return self.update(renumbered)
#self.add_comments('tbd') #self.add_comments('tbd')

70
python/damask/grid_filters.py
View File

@ -1,14 +1,14 @@
from scipy import spatial from scipy import spatial
import numpy as np import numpy as np
def __ks(size,grid,first_order=False): def _ks(size,grid,first_order=False):
""" """
Get wave numbers operator. Get wave numbers operator.
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
""" """
k_sk = np.where(np.arange(grid[0])>grid[0]//2,np.arange(grid[0])-grid[0],np.arange(grid[0]))/size[0] k_sk = np.where(np.arange(grid[0])>grid[0]//2,np.arange(grid[0])-grid[0],np.arange(grid[0]))/size[0]
@ -30,14 +30,14 @@ def curl(size,field):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
""" """
n = np.prod(field.shape[3:]) n = np.prod(field.shape[3:])
k_s = __ks(size,field.shape[:3],True) k_s = _ks(size,field.shape[:3],True)
e = np.zeros((3, 3, 3)) e = np.zeros((3, 3, 3))
e[0, 1, 2] = e[1, 2, 0] = e[2, 0, 1] = +1.0 # Levi-Civita symbol e[0, 1, 2] = e[1, 2, 0] = e[2, 0, 1] = +1.0 # Levi-Civita symbol
e[0, 2, 1] = e[2, 1, 0] = e[1, 0, 2] = -1.0 e[0, 2, 1] = e[2, 1, 0] = e[1, 0, 2] = -1.0
field_fourier = np.fft.rfftn(field,axes=(0,1,2)) field_fourier = np.fft.rfftn(field,axes=(0,1,2))
@ -54,11 +54,11 @@ def divergence(size,field):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
""" """
n = np.prod(field.shape[3:]) n = np.prod(field.shape[3:])
k_s = __ks(size,field.shape[:3],True) k_s = _ks(size,field.shape[:3],True)
field_fourier = np.fft.rfftn(field,axes=(0,1,2)) field_fourier = np.fft.rfftn(field,axes=(0,1,2))
divergence = (np.einsum('ijkl,ijkl ->ijk', k_s,field_fourier)*2.0j*np.pi if n == 3 else # vector, 3 -> 1 divergence = (np.einsum('ijkl,ijkl ->ijk', k_s,field_fourier)*2.0j*np.pi if n == 3 else # vector, 3 -> 1
@ -74,11 +74,11 @@ def gradient(size,field):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
""" """
n = np.prod(field.shape[3:]) n = np.prod(field.shape[3:])
k_s = __ks(size,field.shape[:3],True) k_s = _ks(size,field.shape[:3],True)
field_fourier = np.fft.rfftn(field,axes=(0,1,2)) field_fourier = np.fft.rfftn(field,axes=(0,1,2))
gradient = (np.einsum('ijkl,ijkm->ijkm', field_fourier,k_s)*2.0j*np.pi if n == 1 else # scalar, 1 -> 3 gradient = (np.einsum('ijkl,ijkm->ijkm', field_fourier,k_s)*2.0j*np.pi if n == 1 else # scalar, 1 -> 3
@ -96,7 +96,7 @@ def cell_coord0(grid,size,origin=np.zeros(3)):
grid : numpy.ndarray grid : numpy.ndarray
number of grid points. number of grid points.
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
origin : numpy.ndarray, optional origin : numpy.ndarray, optional
physical origin of the periodic field. Default is [0.0,0.0,0.0]. physical origin of the periodic field. Default is [0.0,0.0,0.0].
@ -108,7 +108,8 @@ def cell_coord0(grid,size,origin=np.zeros(3)):
np.linspace(start[0],end[0],grid[0]), np.linspace(start[0],end[0],grid[0]),
indexing = 'ij') indexing = 'ij')
return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3)
def cell_displacement_fluct(size,F): def cell_displacement_fluct(size,F):
""" """
@ -117,14 +118,14 @@ def cell_displacement_fluct(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
""" """
integrator = 0.5j*size/np.pi integrator = 0.5j*size/np.pi
k_s = __ks(size,F.shape[:3],False) k_s = _ks(size,F.shape[:3],False)
k_s_squared = np.einsum('...l,...l',k_s,k_s) k_s_squared = np.einsum('...l,...l',k_s,k_s)
k_s_squared[0,0,0] = 1.0 k_s_squared[0,0,0] = 1.0
@ -136,6 +137,7 @@ def cell_displacement_fluct(size,F):
return np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3]) return np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3])
def cell_displacement_avg(size,F): def cell_displacement_avg(size,F):
""" """
Cell center displacement field from average part of the deformation gradient field. Cell center displacement field from average part of the deformation gradient field.
@ -143,7 +145,7 @@ def cell_displacement_avg(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
@ -151,6 +153,7 @@ def cell_displacement_avg(size,F):
F_avg = np.average(F,axis=(0,1,2)) F_avg = np.average(F,axis=(0,1,2))
return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),cell_coord0(F.shape[:3][::-1],size)) return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),cell_coord0(F.shape[:3][::-1],size))
def cell_displacement(size,F): def cell_displacement(size,F):
""" """
Cell center displacement field from deformation gradient field. Cell center displacement field from deformation gradient field.
@ -158,13 +161,14 @@ def cell_displacement(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
""" """
return cell_displacement_avg(size,F) + cell_displacement_fluct(size,F) return cell_displacement_avg(size,F) + cell_displacement_fluct(size,F)
def cell_coord(size,F,origin=np.zeros(3)): def cell_coord(size,F,origin=np.zeros(3)):
""" """
Cell center positions. Cell center positions.
@ -172,7 +176,7 @@ def cell_coord(size,F,origin=np.zeros(3)):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
origin : numpy.ndarray, optional origin : numpy.ndarray, optional
@ -181,6 +185,7 @@ def cell_coord(size,F,origin=np.zeros(3)):
""" """
return cell_coord0(F.shape[:3][::-1],size,origin) + cell_displacement(size,F) return cell_coord0(F.shape[:3][::-1],size,origin) + cell_displacement(size,F)
def cell_coord0_gridSizeOrigin(coord0,ordered=True): def cell_coord0_gridSizeOrigin(coord0,ordered=True):
""" """
Return grid 'DNA', i.e. grid, size, and origin from array of cell positions. Return grid 'DNA', i.e. grid, size, and origin from array of cell positions.
@ -200,11 +205,11 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True):
size = grid/np.maximum(grid-1,1) * (maxcorner-mincorner) size = grid/np.maximum(grid-1,1) * (maxcorner-mincorner)
delta = size/grid delta = size/grid
origin = mincorner - delta*.5 origin = mincorner - delta*.5
# 1D/2D: size/origin combination undefined, set origin to 0.0 # 1D/2D: size/origin combination undefined, set origin to 0.0
size [np.where(grid==1)] = origin[np.where(grid==1)]*2. size [np.where(grid==1)] = origin[np.where(grid==1)]*2.
origin[np.where(grid==1)] = 0.0 origin[np.where(grid==1)] = 0.0
if grid.prod() != len(coord0): if grid.prod() != len(coord0):
raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid)) raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid))
@ -221,6 +226,7 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True):
return (grid,size,origin) return (grid,size,origin)
def coord0_check(coord0): def coord0_check(coord0):
""" """
Check whether coordinates lie on a regular grid. Check whether coordinates lie on a regular grid.
@ -234,7 +240,6 @@ def coord0_check(coord0):
cell_coord0_gridSizeOrigin(coord0,ordered=True) cell_coord0_gridSizeOrigin(coord0,ordered=True)
def node_coord0(grid,size,origin=np.zeros(3)): def node_coord0(grid,size,origin=np.zeros(3)):
""" """
Nodal positions (undeformed). Nodal positions (undeformed).
@ -244,7 +249,7 @@ def node_coord0(grid,size,origin=np.zeros(3)):
grid : numpy.ndarray grid : numpy.ndarray
number of grid points. number of grid points.
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
origin : numpy.ndarray, optional origin : numpy.ndarray, optional
physical origin of the periodic field. Default is [0.0,0.0,0.0]. physical origin of the periodic field. Default is [0.0,0.0,0.0].
@ -253,8 +258,9 @@ def node_coord0(grid,size,origin=np.zeros(3)):
np.linspace(origin[1],size[1]+origin[1],1+grid[1]), np.linspace(origin[1],size[1]+origin[1],1+grid[1]),
np.linspace(origin[0],size[0]+origin[0],1+grid[0]), np.linspace(origin[0],size[0]+origin[0],1+grid[0]),
indexing = 'ij') indexing = 'ij')
return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3)
def node_displacement_fluct(size,F): def node_displacement_fluct(size,F):
""" """
@ -263,13 +269,14 @@ def node_displacement_fluct(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
""" """
return cell_2_node(cell_displacement_fluct(size,F)) return cell_2_node(cell_displacement_fluct(size,F))
def node_displacement_avg(size,F): def node_displacement_avg(size,F):
""" """
Nodal displacement field from average part of the deformation gradient field. Nodal displacement field from average part of the deformation gradient field.
@ -277,7 +284,7 @@ def node_displacement_avg(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
@ -285,6 +292,7 @@ def node_displacement_avg(size,F):
F_avg = np.average(F,axis=(0,1,2)) F_avg = np.average(F,axis=(0,1,2))
return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),node_coord0(F.shape[:3][::-1],size)) return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),node_coord0(F.shape[:3][::-1],size))
def node_displacement(size,F): def node_displacement(size,F):
""" """
Nodal displacement field from deformation gradient field. Nodal displacement field from deformation gradient field.
@ -292,13 +300,14 @@ def node_displacement(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
""" """
return node_displacement_avg(size,F) + node_displacement_fluct(size,F) return node_displacement_avg(size,F) + node_displacement_fluct(size,F)
def node_coord(size,F,origin=np.zeros(3)): def node_coord(size,F,origin=np.zeros(3)):
""" """
Nodal positions. Nodal positions.
@ -306,7 +315,7 @@ def node_coord(size,F,origin=np.zeros(3)):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
origin : numpy.ndarray, optional origin : numpy.ndarray, optional
@ -315,22 +324,25 @@ def node_coord(size,F,origin=np.zeros(3)):
""" """
return node_coord0(F.shape[:3][::-1],size,origin) + node_displacement(size,F) return node_coord0(F.shape[:3][::-1],size,origin) + node_displacement(size,F)
def cell_2_node(cell_data): def cell_2_node(cell_data):
"""Interpolate periodic cell data to nodal data.""" """Interpolate periodic cell data to nodal data."""
n = ( cell_data + np.roll(cell_data,1,(0,1,2)) n = ( cell_data + np.roll(cell_data,1,(0,1,2))
+ np.roll(cell_data,1,(0,)) + np.roll(cell_data,1,(1,)) + np.roll(cell_data,1,(2,)) + np.roll(cell_data,1,(0,)) + np.roll(cell_data,1,(1,)) + np.roll(cell_data,1,(2,))
+ np.roll(cell_data,1,(0,1)) + np.roll(cell_data,1,(1,2)) + np.roll(cell_data,1,(2,0)))*0.125 + np.roll(cell_data,1,(0,1)) + np.roll(cell_data,1,(1,2)) + np.roll(cell_data,1,(2,0)))*0.125
return np.pad(n,((0,1),(0,1),(0,1))+((0,0),)*len(cell_data.shape[3:]),mode='wrap') return np.pad(n,((0,1),(0,1),(0,1))+((0,0),)*len(cell_data.shape[3:]),mode='wrap')
def node_2_cell(node_data): def node_2_cell(node_data):
"""Interpolate periodic nodal data to cell data.""" """Interpolate periodic nodal data to cell data."""
c = ( node_data + np.roll(node_data,1,(0,1,2)) c = ( node_data + np.roll(node_data,1,(0,1,2))
+ np.roll(node_data,1,(0,)) + np.roll(node_data,1,(1,)) + np.roll(node_data,1,(2,)) + np.roll(node_data,1,(0,)) + np.roll(node_data,1,(1,)) + np.roll(node_data,1,(2,))
+ np.roll(node_data,1,(0,1)) + np.roll(node_data,1,(1,2)) + np.roll(node_data,1,(2,0)))*0.125 + np.roll(node_data,1,(0,1)) + np.roll(node_data,1,(1,2)) + np.roll(node_data,1,(2,0)))*0.125
return c[:-1,:-1,:-1] return c[:-1,:-1,:-1]
def node_coord0_gridSizeOrigin(coord0,ordered=False): def node_coord0_gridSizeOrigin(coord0,ordered=False):
""" """
Return grid 'DNA', i.e. grid, size, and origin from array of nodal positions. Return grid 'DNA', i.e. grid, size, and origin from array of nodal positions.
@ -349,7 +361,7 @@ def node_coord0_gridSizeOrigin(coord0,ordered=False):
grid = np.array(list(map(len,coords)),'i') - 1 grid = np.array(list(map(len,coords)),'i') - 1
size = maxcorner-mincorner size = maxcorner-mincorner
origin = mincorner origin = mincorner
if (grid+1).prod() != len(coord0): if (grid+1).prod() != len(coord0):
raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid)) raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid))

641
python/damask/lattice.py Normal file
View File

@ -0,0 +1,641 @@
import numpy as np
from .rotation import Rotation
P = -1
# ******************************************************************************************
class Symmetry:
"""
Symmetry operations for lattice systems.
References
----------
https://en.wikipedia.org/wiki/Crystal_system
"""
lattices = [None,'orthorhombic','tetragonal','hexagonal','cubic',]
def __init__(self, symmetry = None):
"""
Symmetry Definition.
Parameters
----------
symmetry : str, optional
label of the crystal system
"""
if symmetry is not None and symmetry.lower() not in Symmetry.lattices:
raise KeyError('Symmetry/crystal system "{}" is unknown'.format(symmetry))
self.lattice = symmetry.lower() if isinstance(symmetry,str) else symmetry
def __copy__(self):
"""Copy."""
return self.__class__(self.lattice)
copy = __copy__
def __repr__(self):
"""Readable string."""
return '{}'.format(self.lattice)
def __eq__(self, other):
"""
Equal to other.
Parameters
----------
other : Symmetry
Symmetry to check for equality.
"""
return self.lattice == other.lattice
def __neq__(self, other):
"""
Not Equal to other.
Parameters
----------
other : Symmetry
Symmetry to check for inequality.
"""
return not self.__eq__(other)
def __cmp__(self,other):
"""
Linear ordering.
Parameters
----------
other : Symmetry
Symmetry to check for for order.
"""
myOrder = Symmetry.lattices.index(self.lattice)
otherOrder = Symmetry.lattices.index(other.lattice)
return (myOrder > otherOrder) - (myOrder < otherOrder)
def symmetryOperations(self,members=[]):
"""List (or single element) of symmetry operations as rotations."""
if self.lattice == 'cubic':
symQuats = [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, 0.0, 0.0, 1.0 ],
[ 0.0, 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2) ],
[ 0.0, 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ],
[ 0.0, -0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ],
[ 0.5, 0.5, 0.5, 0.5 ],
[-0.5, 0.5, 0.5, 0.5 ],
[-0.5, 0.5, 0.5, -0.5 ],
[-0.5, 0.5, -0.5, 0.5 ],
[-0.5, -0.5, 0.5, 0.5 ],
[-0.5, -0.5, 0.5, -0.5 ],
[-0.5, -0.5, -0.5, 0.5 ],
[-0.5, 0.5, -0.5, -0.5 ],
[-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[-0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2), 0.0 ],
[-0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2), 0.0 ],
[-0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0, 0.0 ],
[-0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0, 0.0 ],
]
elif self.lattice == 'hexagonal':
symQuats = [
[ 1.0, 0.0, 0.0, 0.0 ],
[-0.5*np.sqrt(3), 0.0, 0.0, -0.5 ],
[ 0.5, 0.0, 0.0, 0.5*np.sqrt(3) ],
[ 0.0, 0.0, 0.0, 1.0 ],
[-0.5, 0.0, 0.0, 0.5*np.sqrt(3) ],
[-0.5*np.sqrt(3), 0.0, 0.0, 0.5 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, -0.5*np.sqrt(3), 0.5, 0.0 ],
[ 0.0, 0.5, -0.5*np.sqrt(3), 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ],
[ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ],
]
elif self.lattice == 'tetragonal':
symQuats = [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, 0.0, 0.0, 1.0 ],
[ 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ],
[ 0.0, -0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ],
[ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
]
elif self.lattice == 'orthorhombic':
symQuats = [
[ 1.0,0.0,0.0,0.0 ],
[ 0.0,1.0,0.0,0.0 ],
[ 0.0,0.0,1.0,0.0 ],
[ 0.0,0.0,0.0,1.0 ],
]
else:
symQuats = [
[ 1.0,0.0,0.0,0.0 ],
]
symOps = list(map(Rotation,
np.array(symQuats)[np.atleast_1d(members) if members != [] else range(len(symQuats))]))
try:
iter(members) # asking for (even empty) list of members?
except TypeError:
return symOps[0] # no, return rotation object
else:
return symOps # yes, return list of rotations
def inFZ(self,rodrigues):
"""
Check whether given Rodrigues-Frank vector falls into fundamental zone of own symmetry.
Fundamental zone in Rodrigues space is point symmetric around origin.
"""
if (len(rodrigues) != 3):
raise ValueError('Input is not a Rodrigues-Frank vector.\n')
if np.any(rodrigues == np.inf): return False
Rabs = abs(rodrigues)
if self.lattice == 'cubic':
return np.sqrt(2.0)-1.0 >= Rabs[0] \
and np.sqrt(2.0)-1.0 >= Rabs[1] \
and np.sqrt(2.0)-1.0 >= Rabs[2] \
and 1.0 >= Rabs[0] + Rabs[1] + Rabs[2]
elif self.lattice == 'hexagonal':
return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] and 1.0 >= Rabs[2] \
and 2.0 >= np.sqrt(3)*Rabs[0] + Rabs[1] \
and 2.0 >= np.sqrt(3)*Rabs[1] + Rabs[0] \
and 2.0 >= np.sqrt(3) + Rabs[2]
elif self.lattice == 'tetragonal':
return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] \
and np.sqrt(2.0) >= Rabs[0] + Rabs[1] \
and np.sqrt(2.0) >= Rabs[2] + 1.0
elif self.lattice == 'orthorhombic':
return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] and 1.0 >= Rabs[2]
else:
return True
def inDisorientationSST(self,rodrigues):
"""
Check whether given Rodrigues-Frank vector (of misorientation) falls into standard stereographic triangle of own symmetry.
References
----------
A. Heinz and P. Neumann, Acta Crystallographica Section A 47:780-789, 1991
https://doi.org/10.1107/S0108767391006864
"""
if (len(rodrigues) != 3):
raise ValueError('Input is not a Rodrigues-Frank vector.\n')
R = rodrigues
epsilon = 0.0
if self.lattice == 'cubic':
return R[0] >= R[1]+epsilon and R[1] >= R[2]+epsilon and R[2] >= epsilon
elif self.lattice == 'hexagonal':
return R[0] >= np.sqrt(3)*(R[1]-epsilon) and R[1] >= epsilon and R[2] >= epsilon
elif self.lattice == 'tetragonal':
return R[0] >= R[1]-epsilon and R[1] >= epsilon and R[2] >= epsilon
elif self.lattice == 'orthorhombic':
return R[0] >= epsilon and R[1] >= epsilon and R[2] >= epsilon
else:
return True
def inSST(self,
vector,
proper = False,
color = False):
"""
Check whether given vector falls into standard stereographic triangle of own symmetry.
proper considers only vectors with z >= 0, hence uses two neighboring SSTs.
Return inverse pole figure color if requested.
Bases are computed from
basis = {'cubic' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
[1.,0.,1.]/np.sqrt(2.), # direction of green
[1.,1.,1.]/np.sqrt(3.)]).T), # direction of blue
'hexagonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
[1.,0.,0.], # direction of green
[np.sqrt(3.),1.,0.]/np.sqrt(4.)]).T), # direction of blue
'tetragonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
[1.,0.,0.], # direction of green
[1.,1.,0.]/np.sqrt(2.)]).T), # direction of blue
'orthorhombic' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
[1.,0.,0.], # direction of green
[0.,1.,0.]]).T), # direction of blue
}
"""
if self.lattice == 'cubic':
basis = {'improper':np.array([ [-1. , 0. , 1. ],
[ np.sqrt(2.) , -np.sqrt(2.) , 0. ],
[ 0. , np.sqrt(3.) , 0. ] ]),
'proper':np.array([ [ 0. , -1. , 1. ],
[-np.sqrt(2.) , np.sqrt(2.) , 0. ],
[ np.sqrt(3.) , 0. , 0. ] ]),
}
elif self.lattice == 'hexagonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -np.sqrt(3.) , 0. ],
[ 0. , 2. , 0. ] ]),
'proper':np.array([ [ 0. , 0. , 1. ],
[-1. , np.sqrt(3.) , 0. ],
[ np.sqrt(3.) , -1. , 0. ] ]),
}
elif self.lattice == 'tetragonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -1. , 0. ],
[ 0. , np.sqrt(2.) , 0. ] ]),
'proper':np.array([ [ 0. , 0. , 1. ],
[-1. , 1. , 0. ],
[ np.sqrt(2.) , 0. , 0. ] ]),
}
elif self.lattice == 'orthorhombic':
basis = {'improper':np.array([ [ 0., 0., 1.],
[ 1., 0., 0.],
[ 0., 1., 0.] ]),
'proper':np.array([ [ 0., 0., 1.],
[-1., 0., 0.],
[ 0., 1., 0.] ]),
}
else: # direct exit for unspecified symmetry
if color:
return (True,np.zeros(3,'d'))
else:
return True
v = np.array(vector,dtype=float)
if proper: # check both improper ...
theComponents = np.around(np.dot(basis['improper'],v),12)
inSST = np.all(theComponents >= 0.0)
if not inSST: # ... and proper SST
theComponents = np.around(np.dot(basis['proper'],v),12)
inSST = np.all(theComponents >= 0.0)
else:
v[2] = abs(v[2]) # z component projects identical
theComponents = np.around(np.dot(basis['improper'],v),12) # for positive and negative values
inSST = np.all(theComponents >= 0.0)
if color: # have to return color array
if inSST:
rgb = np.power(theComponents/np.linalg.norm(theComponents),0.5) # smoothen color ramps
rgb = np.minimum(np.ones(3,dtype=float),rgb) # limit to maximum intensity
rgb /= max(rgb) # normalize to (HS)V = 1
else:
rgb = np.zeros(3,dtype=float)
return (inSST,rgb)
else:
return inSST
# code derived from https://github.com/ezag/pyeuclid
# suggested reading: http://web.mit.edu/2.998/www/QuaternionReport1.pdf
# ******************************************************************************************
class Lattice:
"""
Lattice system.
Currently, this contains only a mapping from Bravais lattice to symmetry
and orientation relationships. It could include twin and slip systems.
References
----------
https://en.wikipedia.org/wiki/Bravais_lattice
"""
lattices = {
'triclinic':{'symmetry':None},
'bct':{'symmetry':'tetragonal'},
'hex':{'symmetry':'hexagonal'},
'fcc':{'symmetry':'cubic','c/a':1.0},
'bcc':{'symmetry':'cubic','c/a':1.0},
}
def __init__(self, lattice):
"""
New lattice of given type.
Parameters
----------
lattice : str
Bravais lattice.
"""
self.lattice = lattice
self.symmetry = Symmetry(self.lattices[lattice]['symmetry'])
def __repr__(self):
"""Report basic lattice information."""
return 'Bravais lattice {} ({} symmetry)'.format(self.lattice,self.symmetry)
# Kurdjomov--Sachs orientation relationship for fcc <-> bcc transformation
# from S. Morito et al., Journal of Alloys and Compounds 577:s587-s592, 2013
# also see K. Kitahara et al., Acta Materialia 54:1279-1288, 2006
KS = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]]],dtype='float'),
'directions': np.array([
[[ -1, 0, 1],[ -1, -1, 1]],
[[ -1, 0, 1],[ -1, 1, -1]],
[[ 0, 1, -1],[ -1, -1, 1]],
[[ 0, 1, -1],[ -1, 1, -1]],
[[ 1, -1, 0],[ -1, -1, 1]],
[[ 1, -1, 0],[ -1, 1, -1]],
[[ 1, 0, -1],[ -1, -1, 1]],
[[ 1, 0, -1],[ -1, 1, -1]],
[[ -1, -1, 0],[ -1, -1, 1]],
[[ -1, -1, 0],[ -1, 1, -1]],
[[ 0, 1, 1],[ -1, -1, 1]],
[[ 0, 1, 1],[ -1, 1, -1]],
[[ 0, -1, 1],[ -1, -1, 1]],
[[ 0, -1, 1],[ -1, 1, -1]],
[[ -1, 0, -1],[ -1, -1, 1]],
[[ -1, 0, -1],[ -1, 1, -1]],
[[ 1, 1, 0],[ -1, -1, 1]],
[[ 1, 1, 0],[ -1, 1, -1]],
[[ -1, 1, 0],[ -1, -1, 1]],
[[ -1, 1, 0],[ -1, 1, -1]],
[[ 0, -1, -1],[ -1, -1, 1]],
[[ 0, -1, -1],[ -1, 1, -1]],
[[ 1, 0, 1],[ -1, -1, 1]],
[[ 1, 0, 1],[ -1, 1, -1]]],dtype='float')}
# Greninger--Troiano orientation relationship for fcc <-> bcc transformation
# from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
GT = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 1, 1, 1],[ 1, 0, 1]],
[[ 1, 1, 1],[ 1, 1, 0]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ -1, -1, 1],[ -1, 0, 1]],
[[ -1, -1, 1],[ -1, -1, 0]],
[[ -1, -1, 1],[ 0, -1, 1]],
[[ -1, 1, 1],[ -1, 0, 1]],
[[ -1, 1, 1],[ -1, 1, 0]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 1, 0, 1]],
[[ 1, -1, 1],[ 1, -1, 0]],
[[ 1, -1, 1],[ 0, -1, 1]],
[[ 1, 1, 1],[ 1, 1, 0]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 1, 0, 1]],
[[ -1, -1, 1],[ -1, -1, 0]],
[[ -1, -1, 1],[ 0, -1, 1]],
[[ -1, -1, 1],[ -1, 0, 1]],
[[ -1, 1, 1],[ -1, 1, 0]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ -1, 0, 1]],
[[ 1, -1, 1],[ 1, -1, 0]],
[[ 1, -1, 1],[ 0, -1, 1]],
[[ 1, -1, 1],[ 1, 0, 1]]],dtype='float'),
'directions': np.array([
[[ -5,-12, 17],[-17, -7, 17]],
[[ 17, -5,-12],[ 17,-17, -7]],
[[-12, 17, -5],[ -7, 17,-17]],
[[ 5, 12, 17],[ 17, 7, 17]],
[[-17, 5,-12],[-17, 17, -7]],
[[ 12,-17, -5],[ 7,-17,-17]],
[[ -5, 12,-17],[-17, 7,-17]],
[[ 17, 5, 12],[ 17, 17, 7]],
[[-12,-17, 5],[ -7,-17, 17]],
[[ 5,-12,-17],[ 17, -7,-17]],
[[-17, -5, 12],[-17,-17, 7]],
[[ 12, 17, 5],[ 7, 17, 17]],
[[ -5, 17,-12],[-17, 17, -7]],
[[-12, -5, 17],[ -7,-17, 17]],
[[ 17,-12, -5],[ 17, -7,-17]],
[[ 5,-17,-12],[ 17,-17, -7]],
[[ 12, 5, 17],[ 7, 17, 17]],
[[-17, 12, -5],[-17, 7,-17]],
[[ -5,-17, 12],[-17,-17, 7]],
[[-12, 5,-17],[ -7, 17,-17]],
[[ 17, 12, 5],[ 17, 7, 17]],
[[ 5, 17, 12],[ 17, 17, 7]],
[[ 12, -5,-17],[ 7,-17,-17]],
[[-17,-12, 5],[-17,-7, 17]]],dtype='float')}
# Greninger--Troiano' orientation relationship for fcc <-> bcc transformation
# from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
GTprime = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 7, 17, 17],[ 12, 5, 17]],
[[ 17, 7, 17],[ 17, 12, 5]],
[[ 17, 17, 7],[ 5, 17, 12]],
[[ -7,-17, 17],[-12, -5, 17]],
[[-17, -7, 17],[-17,-12, 5]],
[[-17,-17, 7],[ -5,-17, 12]],
[[ 7,-17,-17],[ 12, -5,-17]],
[[ 17, -7,-17],[ 17,-12, -5]],
[[ 17,-17, -7],[ 5,-17,-12]],
[[ -7, 17,-17],[-12, 5,-17]],
[[-17, 7,-17],[-17, 12, -5]],
[[-17, 17, -7],[ -5, 17,-12]],
[[ 7, 17, 17],[ 12, 17, 5]],
[[ 17, 7, 17],[ 5, 12, 17]],
[[ 17, 17, 7],[ 17, 5, 12]],
[[ -7,-17, 17],[-12,-17, 5]],
[[-17, -7, 17],[ -5,-12, 17]],
[[-17,-17, 7],[-17, -5, 12]],
[[ 7,-17,-17],[ 12,-17, -5]],
[[ 17, -7,-17],[ 5, -12,-17]],
[[ 17,-17, -7],[ 17, -5,-12]],
[[ -7, 17,-17],[-12, 17, -5]],
[[-17, 7,-17],[ -5, 12,-17]],
[[-17, 17, -7],[-17, 5,-12]]],dtype='float'),
'directions': np.array([
[[ 0, 1, -1],[ 1, 1, -1]],
[[ -1, 0, 1],[ -1, 1, 1]],
[[ 1, -1, 0],[ 1, -1, 1]],
[[ 0, -1, -1],[ -1, -1, -1]],
[[ 1, 0, 1],[ 1, -1, 1]],
[[ 1, -1, 0],[ 1, -1, -1]],
[[ 0, 1, -1],[ -1, 1, -1]],
[[ 1, 0, 1],[ 1, 1, 1]],
[[ -1, -1, 0],[ -1, -1, 1]],
[[ 0, -1, -1],[ 1, -1, -1]],
[[ -1, 0, 1],[ -1, -1, 1]],
[[ -1, -1, 0],[ -1, -1, -1]],
[[ 0, -1, 1],[ 1, -1, 1]],
[[ 1, 0, -1],[ 1, 1, -1]],
[[ -1, 1, 0],[ -1, 1, 1]],
[[ 0, 1, 1],[ -1, 1, 1]],
[[ -1, 0, -1],[ -1, -1, -1]],
[[ -1, 1, 0],[ -1, 1, -1]],
[[ 0, -1, 1],[ -1, -1, 1]],
[[ -1, 0, -1],[ -1, 1, -1]],
[[ 1, 1, 0],[ 1, 1, 1]],
[[ 0, 1, 1],[ 1, 1, 1]],
[[ 1, 0, -1],[ 1, -1, -1]],
[[ 1, 1, 0],[ 1, 1, -1]]],dtype='float')}
# Nishiyama--Wassermann orientation relationship for fcc <-> bcc transformation
# from H. Kitahara et al., Materials Characterization 54:378-386, 2005
NW = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ 1, -1, 1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]]],dtype='float'),
'directions': np.array([
[[ 2, -1, -1],[ 0, -1, 1]],
[[ -1, 2, -1],[ 0, -1, 1]],
[[ -1, -1, 2],[ 0, -1, 1]],
[[ -2, -1, -1],[ 0, -1, 1]],
[[ 1, 2, -1],[ 0, -1, 1]],
[[ 1, -1, 2],[ 0, -1, 1]],
[[ 2, 1, -1],[ 0, -1, 1]],
[[ -1, -2, -1],[ 0, -1, 1]],
[[ -1, 1, 2],[ 0, -1, 1]],
[[ 2, -1, 1],[ 0, -1, 1]], #It is wrong in the paper, but matrix is correct
[[ -1, 2, 1],[ 0, -1, 1]],
[[ -1, -1, -2],[ 0, -1, 1]]],dtype='float')}
# Pitsch orientation relationship for fcc <-> bcc transformation
# from Y. He et al., Acta Materialia 53:1179-1190, 2005
Pitsch = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 0, 1, 0],[ -1, 0, 1]],
[[ 0, 0, 1],[ 1, -1, 0]],
[[ 1, 0, 0],[ 0, 1, -1]],
[[ 1, 0, 0],[ 0, -1, -1]],
[[ 0, 1, 0],[ -1, 0, -1]],
[[ 0, 0, 1],[ -1, -1, 0]],
[[ 0, 1, 0],[ -1, 0, -1]],
[[ 0, 0, 1],[ -1, -1, 0]],
[[ 1, 0, 0],[ 0, -1, -1]],
[[ 1, 0, 0],[ 0, -1, 1]],
[[ 0, 1, 0],[ 1, 0, -1]],
[[ 0, 0, 1],[ -1, 1, 0]]],dtype='float'),
'directions': np.array([
[[ 1, 0, 1],[ 1, -1, 1]],
[[ 1, 1, 0],[ 1, 1, -1]],
[[ 0, 1, 1],[ -1, 1, 1]],
[[ 0, 1, -1],[ -1, 1, -1]],
[[ -1, 0, 1],[ -1, -1, 1]],
[[ 1, -1, 0],[ 1, -1, -1]],
[[ 1, 0, -1],[ 1, -1, -1]],
[[ -1, 1, 0],[ -1, 1, -1]],
[[ 0, -1, 1],[ -1, -1, 1]],
[[ 0, 1, 1],[ -1, 1, 1]],
[[ 1, 0, 1],[ 1, -1, 1]],
[[ 1, 1, 0],[ 1, 1, -1]]],dtype='float')}
# Bain orientation relationship for fcc <-> bcc transformation
# from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
Bain = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 1, 0, 0],[ 1, 0, 0]],
[[ 0, 1, 0],[ 0, 1, 0]],
[[ 0, 0, 1],[ 0, 0, 1]]],dtype='float'),
'directions': np.array([
[[ 0, 1, 0],[ 0, 1, 1]],
[[ 0, 0, 1],[ 1, 0, 1]],
[[ 1, 0, 0],[ 1, 1, 0]]],dtype='float')}
def relationOperations(self,model):
"""
Crystallographic orientation relationships for phase transformations.
References
----------
S. Morito et al., Journal of Alloys and Compounds 577:s587-s592, 2013
https://doi.org/10.1016/j.jallcom.2012.02.004
K. Kitahara et al., Acta Materialia 54(5):1279-1288, 2006
https://doi.org/10.1016/j.actamat.2005.11.001
Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
https://doi.org/10.1107/S0021889805038276
H. Kitahara et al., Materials Characterization 54(4-5):378-386, 2005
https://doi.org/10.1016/j.matchar.2004.12.015
Y. He et al., Acta Materialia 53(4):1179-1190, 2005
https://doi.org/10.1016/j.actamat.2004.11.021
"""
models={'KS':self.KS, 'GT':self.GT, 'GT_prime':self.GTprime,
'NW':self.NW, 'Pitsch': self.Pitsch, 'Bain':self.Bain}
try:
relationship = models[model]
except KeyError :
raise KeyError('Orientation relationship "{}" is unknown'.format(model))
if self.lattice not in relationship['mapping']:
raise ValueError('Relationship "{}" not supported for lattice "{}"'.format(model,self.lattice))
r = {'lattice':Lattice((set(relationship['mapping'])-{self.lattice}).pop()), # target lattice
'rotations':[] }
myPlane_id = relationship['mapping'][self.lattice]
otherPlane_id = (myPlane_id+1)%2
myDir_id = myPlane_id +2
otherDir_id = otherPlane_id +2
for miller in np.hstack((relationship['planes'],relationship['directions'])):
myPlane = miller[myPlane_id]/ np.linalg.norm(miller[myPlane_id])
myDir = miller[myDir_id]/ np.linalg.norm(miller[myDir_id])
myMatrix = np.array([myDir,np.cross(myPlane,myDir),myPlane])
otherPlane = miller[otherPlane_id]/ np.linalg.norm(miller[otherPlane_id])
otherDir = miller[otherDir_id]/ np.linalg.norm(miller[otherDir_id])
otherMatrix = np.array([otherDir,np.cross(otherPlane,otherDir),otherPlane])
r['rotations'].append(Rotation.fromMatrix(np.dot(otherMatrix.T,myMatrix)))
return r

97
python/damask/mechanics.py
View File

@ -21,21 +21,21 @@ def Cauchy(P,F):
return symmetric(sigma) return symmetric(sigma)
def deviatoric_part(x): def deviatoric_part(T):
""" """
Return deviatoric part of a tensor. Return deviatoric part of a tensor.
Parameters Parameters
---------- ----------
x : numpy.array of shape (:,3,3) or (3,3) T : numpy.array of shape (:,3,3) or (3,3)
Tensor of which the deviatoric part is computed. Tensor of which the deviatoric part is computed.
""" """
return x - np.eye(3)*spherical_part(x) if np.shape(x) == (3,3) else \ return T - np.eye(3)*spherical_part(T) if np.shape(T) == (3,3) else \
x - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[x.shape[0],3,3]),spherical_part(x)) T - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[T.shape[0],3,3]),spherical_part(T))
def eigenvalues(x): def eigenvalues(T_sym):
""" """
Return the eigenvalues, i.e. principal components, of a symmetric tensor. Return the eigenvalues, i.e. principal components, of a symmetric tensor.
@ -44,14 +44,14 @@ def eigenvalues(x):
Parameters Parameters
---------- ----------
x : numpy.array of shape (:,3,3) or (3,3) T_sym : numpy.array of shape (:,3,3) or (3,3)
Symmetric tensor of which the eigenvalues are computed. Symmetric tensor of which the eigenvalues are computed.
""" """
return np.linalg.eigvalsh(symmetric(x)) return np.linalg.eigvalsh(symmetric(T_sym))
def eigenvectors(x,RHS=False): def eigenvectors(T_sym,RHS=False):
""" """
Return eigenvectors of a symmetric tensor. Return eigenvectors of a symmetric tensor.
@ -59,47 +59,47 @@ def eigenvectors(x,RHS=False):
Parameters Parameters
---------- ----------
x : numpy.array of shape (:,3,3) or (3,3) T_sym : numpy.array of shape (:,3,3) or (3,3)
Symmetric tensor of which the eigenvectors are computed. Symmetric tensor of which the eigenvectors are computed.
RHS: bool, optional RHS: bool, optional
Enforce right-handed coordinate system. Default is False. Enforce right-handed coordinate system. Default is False.
""" """
(u,v) = np.linalg.eigh(symmetric(x)) (u,v) = np.linalg.eigh(symmetric(T_sym))
if RHS: if RHS:
if np.shape(x) == (3,3): if np.shape(T_sym) == (3,3):
if np.linalg.det(v) < 0.0: v[:,2] *= -1.0 if np.linalg.det(v) < 0.0: v[:,2] *= -1.0
else: else:
v[np.linalg.det(v) < 0.0,:,2] *= -1.0 v[np.linalg.det(v) < 0.0,:,2] *= -1.0
return v return v
def left_stretch(x): def left_stretch(T):
""" """
Return the left stretch of a tensor. Return the left stretch of a tensor.
Parameters Parameters
---------- ----------
x : numpy.array of shape (:,3,3) or (3,3) T : numpy.array of shape (:,3,3) or (3,3)
Tensor of which the left stretch is computed. Tensor of which the left stretch is computed.
""" """
return __polar_decomposition(x,'V')[0] return __polar_decomposition(T,'V')[0]
def maximum_shear(x): def maximum_shear(T_sym):
""" """
Return the maximum shear component of a symmetric tensor. Return the maximum shear component of a symmetric tensor.
Parameters Parameters
---------- ----------
x : numpy.array of shape (:,3,3) or (3,3) T_sym : numpy.array of shape (:,3,3) or (3,3)
Symmetric tensor of which the maximum shear is computed. Symmetric tensor of which the maximum shear is computed.
""" """
w = eigenvalues(x) w = eigenvalues(T_sym)
return (w[0] - w[2])*0.5 if np.shape(x) == (3,3) else \ return (w[0] - w[2])*0.5 if np.shape(T_sym) == (3,3) else \
(w[:,0] - w[:,2])*0.5 (w[:,0] - w[:,2])*0.5
@ -147,53 +147,54 @@ def PK2(P,F):
S = np.einsum('ijk,ikl->ijl',np.linalg.inv(F),P) S = np.einsum('ijk,ikl->ijl',np.linalg.inv(F),P)
return symmetric(S) return symmetric(S)
def right_stretch(x):
def right_stretch(T):
""" """
Return the right stretch of a tensor. Return the right stretch of a tensor.
Parameters Parameters
---------- ----------
x : numpy.array of shape (:,3,3) or (3,3) T : numpy.array of shape (:,3,3) or (3,3)
Tensor of which the right stretch is computed. Tensor of which the right stretch is computed.
""" """
return __polar_decomposition(x,'U')[0] return __polar_decomposition(T,'U')[0]
def rotational_part(x): def rotational_part(T):
""" """
Return the rotational part of a tensor. Return the rotational part of a tensor.
Parameters Parameters
---------- ----------
x : numpy.array of shape (:,3,3) or (3,3) T : numpy.array of shape (:,3,3) or (3,3)
Tensor of which the rotational part is computed. Tensor of which the rotational part is computed.
""" """
return __polar_decomposition(x,'R')[0] return __polar_decomposition(T,'R')[0]
def spherical_part(x,tensor=False): def spherical_part(T,tensor=False):
""" """
Return spherical (hydrostatic) part of a tensor. Return spherical (hydrostatic) part of a tensor.
Parameters Parameters
---------- ----------
x : numpy.array of shape (:,3,3) or (3,3) T : numpy.array of shape (:,3,3) or (3,3)
Tensor of which the hydrostatic part is computed. Tensor of which the hydrostatic part is computed.
tensor : bool, optional tensor : bool, optional
Map spherical part onto identity tensor. Default is false Map spherical part onto identity tensor. Default is false
""" """
if x.shape == (3,3): if T.shape == (3,3):
sph = np.trace(x)/3.0 sph = np.trace(T)/3.0
return sph if not tensor else np.eye(3)*sph return sph if not tensor else np.eye(3)*sph
else: else:
sph = np.trace(x,axis1=1,axis2=2)/3.0 sph = np.trace(T,axis1=1,axis2=2)/3.0
if not tensor: if not tensor:
return sph return sph
else: else:
return np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),(x.shape[0],3,3)),sph) return np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),(T.shape[0],3,3)),sph)
def strain_tensor(F,t,m): def strain_tensor(F,t,m):
@ -234,73 +235,73 @@ def strain_tensor(F,t,m):
eps eps
def symmetric(x): def symmetric(T):
""" """
Return the symmetrized tensor. Return the symmetrized tensor.
Parameters Parameters
---------- ----------
x : numpy.array of shape (:,3,3) or (3,3) T : numpy.array of shape (:,3,3) or (3,3)
Tensor of which the symmetrized values are computed. Tensor of which the symmetrized values are computed.
""" """
return (x+transpose(x))*0.5 return (T+transpose(T))*0.5
def transpose(x): def transpose(T):
""" """
Return the transpose of a tensor. Return the transpose of a tensor.
Parameters Parameters
---------- ----------
x : numpy.array of shape (:,3,3) or (3,3) T : numpy.array of shape (:,3,3) or (3,3)
Tensor of which the transpose is computed. Tensor of which the transpose is computed.
""" """
return x.T if np.shape(x) == (3,3) else \ return T.T if np.shape(T) == (3,3) else \
np.transpose(x,(0,2,1)) np.transpose(T,(0,2,1))
def __polar_decomposition(x,requested): def __polar_decomposition(T,requested):
""" """
Singular value decomposition. Singular value decomposition.
Parameters Parameters
---------- ----------
x : numpy.array of shape (:,3,3) or (3,3) T : numpy.array of shape (:,3,3) or (3,3)
Tensor of which the singular values are computed. Tensor of which the singular values are computed.
requested : iterable of str requested : iterable of str
Requested outputs: ‘R’ for the rotation tensor, Requested outputs: ‘R’ for the rotation tensor,
‘V’ for left stretch tensor and ‘U’ for right stretch tensor. ‘V’ for left stretch tensor and ‘U’ for right stretch tensor.
""" """
u, s, vh = np.linalg.svd(x) u, s, vh = np.linalg.svd(T)
R = np.dot(u,vh) if np.shape(x) == (3,3) else \ R = np.dot(u,vh) if np.shape(T) == (3,3) else \
np.einsum('ijk,ikl->ijl',u,vh) np.einsum('ijk,ikl->ijl',u,vh)
output = [] output = []
if 'R' in requested: if 'R' in requested:
output.append(R) output.append(R)
if 'V' in requested: if 'V' in requested:
output.append(np.dot(x,R.T) if np.shape(x) == (3,3) else np.einsum('ijk,ilk->ijl',x,R)) output.append(np.dot(T,R.T) if np.shape(T) == (3,3) else np.einsum('ijk,ilk->ijl',T,R))
if 'U' in requested: if 'U' in requested:
output.append(np.dot(R.T,x) if np.shape(x) == (3,3) else np.einsum('ikj,ikl->ijl',R,x)) output.append(np.dot(R.T,T) if np.shape(T) == (3,3) else np.einsum('ikj,ikl->ijl',R,T))
return tuple(output) return tuple(output)
def __Mises(x,s): def __Mises(T_sym,s):
""" """
Base equation for Mises equivalent of a stres or strain tensor. Base equation for Mises equivalent of a stres or strain tensor.
Parameters Parameters
---------- ----------
x : numpy.array of shape (:,3,3) or (3,3) T_sym : numpy.array of shape (:,3,3) or (3,3)
Symmetric tensor of which the von Mises equivalent is computed. Symmetric tensor of which the von Mises equivalent is computed.
s : float s : float
Scaling factor (2/3 for strain, 3/2 for stress). Scaling factor (2/3 for strain, 3/2 for stress).
""" """
d = deviatoric_part(x) d = deviatoric_part(T_sym)
return np.sqrt(s*(np.sum(d**2.0))) if np.shape(x) == (3,3) else \ return np.sqrt(s*(np.sum(d**2.0))) if np.shape(T_sym) == (3,3) else \
np.sqrt(s*np.einsum('ijk->i',d**2.0)) np.sqrt(s*np.einsum('ijk->i',d**2.0))

1689
python/damask/orientation.py
View File

File diff suppressed because it is too large Load Diff

1139
python/damask/result.py Normal file
View File

File diff suppressed because it is too large Load Diff

837
python/damask/rotation.py Normal file
View File

@ -0,0 +1,837 @@
import numpy as np
from . import Lambert
P = -1
def iszero(a):
return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0)
class Rotation:
u"""
Orientation stored with functionality for conversion to different representations.
References
----------
D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015
https://doi.org/10.1088/0965-0393/23/8/083501
Conventions
-----------
Convention 1: Coordinate frames are right-handed.
Convention 2: A rotation angle ω is taken to be positive for a counterclockwise rotation
when viewing from the end point of the rotation axis towards the origin.
Convention 3: Rotations will be interpreted in the passive sense.
Convention 4: Euler angle triplets are implemented using the Bunge convention,
with the angular ranges as [0, 2π],[0, π],[0, 2π].
Convention 5: The rotation angle ω is limited to the interval [0, π].
Convention 6: the real part of a quaternion is positive, Re(q) > 0
Convention 7: P = -1 (as default).
Usage
-----
Vector "a" (defined in coordinate system "A") is passively rotated
resulting in new coordinates "b" when expressed in system "B".
b = Q * a
b = np.dot(Q.asMatrix(),a)
"""
__slots__ = ['quaternion']
def __init__(self,quaternion = np.array([1.0,0.0,0.0,0.0])):
"""
Initializes to identity unless specified.
Parameters
----------
quaternion : numpy.ndarray, optional
Unit quaternion that follows the conventions. Use .fromQuaternion to perform a sanity check.
"""
self.quaternion = quaternion.copy()
def __copy__(self):
"""Copy."""
return self.__class__(self.quaternion)
copy = __copy__
def __repr__(self):
"""Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles."""
return '\n'.join([
'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)),
'Matrix:\n{}'.format(self.asMatrix()),
'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.asEulers(degrees=True)),
])
def __mul__(self, other):
"""
Multiplication.
Parameters
----------
other : numpy.ndarray or Rotation
Vector, second or fourth order tensor, or rotation object that is rotated.
Todo
----
Document details active/passive)
considere rotation of (3,3,3,3)-matrix
"""
if isinstance(other, Rotation): # rotate a rotation
self_q = self.quaternion[0]
self_p = self.quaternion[1:]
other_q = other.quaternion[0]
other_p = other.quaternion[1:]
R = self.__class__(np.append(self_q*other_q - np.dot(self_p,other_p),
self_q*other_p + other_q*self_p + P * np.cross(self_p,other_p)))
return R.standardize()
elif isinstance(other, (tuple,np.ndarray)):
if isinstance(other,tuple) or other.shape == (3,): # rotate a single (3)-vector or meshgrid
A = self.quaternion[0]**2.0 - np.dot(self.quaternion[1:],self.quaternion[1:])
B = 2.0 * ( self.quaternion[1]*other[0]
+ self.quaternion[2]*other[1]
+ self.quaternion[3]*other[2])
C = 2.0 * P*self.quaternion[0]
return np.array([
A*other[0] + B*self.quaternion[1] + C*(self.quaternion[2]*other[2] - self.quaternion[3]*other[1]),
A*other[1] + B*self.quaternion[2] + C*(self.quaternion[3]*other[0] - self.quaternion[1]*other[2]),
A*other[2] + B*self.quaternion[3] + C*(self.quaternion[1]*other[1] - self.quaternion[2]*other[0]),
])
elif other.shape == (3,3,): # rotate a single (3x3)-matrix
return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T))
elif other.shape == (3,3,3,3,):
raise NotImplementedError
else:
return NotImplemented
else:
return NotImplemented
def inverse(self):
"""In-place inverse rotation/backward rotation."""
self.quaternion[1:] *= -1
return self
def inversed(self):
"""Inverse rotation/backward rotation."""
return self.copy().inverse()
def standardize(self):
"""In-place quaternion representation with positive q."""
if self.quaternion[0] < 0.0: self.quaternion*=-1
return self
def standardized(self):
"""Quaternion representation with positive q."""
return self.copy().standardize()
def misorientation(self,other):
"""
Get Misorientation.
Parameters
----------
other : Rotation
Rotation to which the misorientation is computed.
"""
return other*self.inversed()
def average(self,other):
"""
Calculate the average rotation.
Parameters
----------
other : Rotation
Rotation from which the average is rotated.
"""
return Rotation.fromAverage([self,other])
################################################################################################
# convert to different orientation representations (numpy arrays)
def asQuaternion(self):
"""
Unit quaternion [q, p_1, p_2, p_3] unless quaternion == True: damask.quaternion object.
Parameters
----------
quaternion : bool, optional
return quaternion as DAMASK object.
"""
return self.quaternion
def asEulers(self,
degrees = False):
"""
Bunge-Euler angles: (φ_1, ϕ, φ_2).
Parameters
----------
degrees : bool, optional
return angles in degrees.
"""
eu = Rotation.qu2eu(self.quaternion)
if degrees: eu = np.degrees(eu)
return eu
def asAxisAngle(self,
degrees = False,
pair = False):
"""
Axis angle representation [n_1, n_2, n_3, ω] unless pair == True: ([n_1, n_2, n_3], ω).
Parameters
----------
degrees : bool, optional
return rotation angle in degrees.
pair : bool, optional
return tuple of axis and angle.
"""
ax = Rotation.qu2ax(self.quaternion)
if degrees: ax[3] = np.degrees(ax[3])
return (ax[:3],np.degrees(ax[3])) if pair else ax
def asMatrix(self):
"""Rotation matrix."""
return Rotation.qu2om(self.quaternion)
def asRodrigues(self,
vector = False):
"""
Rodrigues-Frank vector representation [n_1, n_2, n_3, tan(ω/2)] unless vector == True: [n_1, n_2, n_3] * tan(ω/2).
Parameters
----------
vector : bool, optional
return as actual Rodrigues--Frank vector, i.e. rotation axis scaled by tan(ω/2).
"""
ro = Rotation.qu2ro(self.quaternion)
return ro[:3]*ro[3] if vector else ro
def asHomochoric(self):
"""Homochoric vector: (h_1, h_2, h_3)."""
return Rotation.qu2ho(self.quaternion)
def asCubochoric(self):
"""Cubochoric vector: (c_1, c_2, c_3)."""
return Rotation.qu2cu(self.quaternion)
def asM(self):
"""
Intermediate representation supporting quaternion averaging.
References
----------
F. Landis Markley et al., Journal of Guidance, Control, and Dynamics 30(4):1193-1197, 2007
https://doi.org/10.2514/1.28949
"""
return np.outer(self.quaternion,self.quaternion)
################################################################################################
# static constructors. The input data needs to follow the convention, options allow to
# relax these convections
@staticmethod
def fromQuaternion(quaternion,
acceptHomomorph = False,
P = -1):
qu = quaternion if isinstance(quaternion,np.ndarray) and quaternion.dtype == np.dtype(float) \
else np.array(quaternion,dtype=float)
if P > 0: qu[1:4] *= -1 # convert from P=1 to P=-1
if qu[0] < 0.0:
if acceptHomomorph:
qu *= -1.
else:
raise ValueError('Quaternion has negative first component.\n{}'.format(qu[0]))
if not np.isclose(np.linalg.norm(qu), 1.0):
raise ValueError('Quaternion is not of unit length.\n{} {} {} {}'.format(*qu))
return Rotation(qu)
@staticmethod
def fromEulers(eulers,
degrees = False):
eu = eulers if isinstance(eulers, np.ndarray) and eulers.dtype == np.dtype(float) \
else np.array(eulers,dtype=float)
eu = np.radians(eu) if degrees else eu
if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or eu[1] > np.pi:
raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].\n{} {} {}.'.format(*eu))
return Rotation(Rotation.eu2qu(eu))
@staticmethod
def fromAxisAngle(angleAxis,
degrees = False,
normalise = False,
P = -1):
ax = angleAxis if isinstance(angleAxis, np.ndarray) and angleAxis.dtype == np.dtype(float) \
else np.array(angleAxis,dtype=float)
if P > 0: ax[0:3] *= -1 # convert from P=1 to P=-1
if degrees: ax[ 3] = np.radians(ax[3])
if normalise: ax[0:3] /= np.linalg.norm(ax[0:3])
if ax[3] < 0.0 or ax[3] > np.pi:
raise ValueError('Axis angle rotation angle outside of [0..π].\n'.format(ax[3]))
if not np.isclose(np.linalg.norm(ax[0:3]), 1.0):
raise ValueError('Axis angle rotation axis is not of unit length.\n{} {} {}'.format(*ax[0:3]))
return Rotation(Rotation.ax2qu(ax))
@staticmethod
def fromBasis(basis,
orthonormal = True,
reciprocal = False,
):
om = basis if isinstance(basis, np.ndarray) else np.array(basis).reshape((3,3))
if reciprocal:
om = np.linalg.inv(om.T/np.pi) # transform reciprocal basis set
orthonormal = False # contains stretch
if not orthonormal:
(U,S,Vh) = np.linalg.svd(om) # singular value decomposition
om = np.dot(U,Vh)
if not np.isclose(np.linalg.det(om),1.0):
raise ValueError('matrix is not a proper rotation.\n{}'.format(om))
if not np.isclose(np.dot(om[0],om[1]), 0.0) \
or not np.isclose(np.dot(om[1],om[2]), 0.0) \
or not np.isclose(np.dot(om[2],om[0]), 0.0):
raise ValueError('matrix is not orthogonal.\n{}'.format(om))
return Rotation(Rotation.om2qu(om))
@staticmethod
def fromMatrix(om,
):
return Rotation.fromBasis(om)
@staticmethod
def fromRodrigues(rodrigues,
normalise = False,
P = -1):
ro = rodrigues if isinstance(rodrigues, np.ndarray) and rodrigues.dtype == np.dtype(float) \
else np.array(rodrigues,dtype=float)
if P > 0: ro[0:3] *= -1 # convert from P=1 to P=-1
if normalise: ro[0:3] /= np.linalg.norm(ro[0:3])
if not np.isclose(np.linalg.norm(ro[0:3]), 1.0):
raise ValueError('Rodrigues rotation axis is not of unit length.\n{} {} {}'.format(*ro[0:3]))
if ro[3] < 0.0:
raise ValueError('Rodriques rotation angle not positive.\n'.format(ro[3]))
return Rotation(Rotation.ro2qu(ro))
@staticmethod
def fromHomochoric(homochoric,
P = -1):
ho = homochoric if isinstance(homochoric, np.ndarray) and homochoric.dtype == np.dtype(float) \
else np.array(homochoric,dtype=float)
if P > 0: ho *= -1 # convert from P=1 to P=-1
return Rotation(Rotation.ho2qu(ho))
@staticmethod
def fromCubochoric(cubochoric,
P = -1):
cu = cubochoric if isinstance(cubochoric, np.ndarray) and cubochoric.dtype == np.dtype(float) \
else np.array(cubochoric,dtype=float)
ho = Rotation.cu2ho(cu)
if P > 0: ho *= -1 # convert from P=1 to P=-1
return Rotation(Rotation.ho2qu(ho))
@staticmethod
def fromAverage(rotations,
weights = []):
"""
Average rotation.
References
----------
F. Landis Markley et al., Journal of Guidance, Control, and Dynamics 30(4):1193-1197, 2007
https://doi.org/10.2514/1.28949
Parameters
----------
rotations : list of Rotations
Rotations to average from
weights : list of floats, optional
Weights for each rotation used for averaging
"""
if not all(isinstance(item, Rotation) for item in rotations):
raise TypeError("Only instances of Rotation can be averaged.")
N = len(rotations)
if weights == [] or not weights:
weights = np.ones(N,dtype='i')
for i,(r,n) in enumerate(zip(rotations,weights)):
M = r.asM() * n if i == 0 \
else M + r.asM() * n # noqa add (multiples) of this rotation to average noqa
eig, vec = np.linalg.eig(M/N)
return Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True)
@staticmethod
def fromRandom():
r = np.random.random(3)
A = np.sqrt(r[2])
B = np.sqrt(1.0-r[2])
return Rotation(np.array([np.cos(2.0*np.pi*r[0])*A,
np.sin(2.0*np.pi*r[1])*B,
np.cos(2.0*np.pi*r[1])*B,
np.sin(2.0*np.pi*r[0])*A])).standardize()
####################################################################################################
# Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations
####################################################################################################
# Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
# Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification, are
# permitted provided that the following conditions are met:
#
# - Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
# - Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or
# other materials provided with the distribution.
# - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
# of its contributors may be used to endorse or promote products derived from
# this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
# USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
####################################################################################################
#---------- Quaternion ----------
@staticmethod
def qu2om(qu):
"""Quaternion to rotation matrix."""
qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2)
om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2)
om[1,0] = 2.0*(qu[2]*qu[1]+qu[0]*qu[3])
om[0,1] = 2.0*(qu[1]*qu[2]-qu[0]*qu[3])
om[2,1] = 2.0*(qu[3]*qu[2]+qu[0]*qu[1])
om[1,2] = 2.0*(qu[2]*qu[3]-qu[0]*qu[1])
om[0,2] = 2.0*(qu[1]*qu[3]+qu[0]*qu[2])
om[2,0] = 2.0*(qu[3]*qu[1]-qu[0]*qu[2])
return om if P > 0.0 else om.T
@staticmethod
def qu2eu(qu):
"""Quaternion to Bunge-Euler angles."""
q03 = qu[0]**2+qu[3]**2
q12 = qu[1]**2+qu[2]**2
chi = np.sqrt(q03*q12)
if iszero(chi):
eu = np.array([np.arctan2(-P*2.0*qu[0]*qu[3],qu[0]**2-qu[3]**2), 0.0, 0.0]) if iszero(q12) else \
np.array([np.arctan2(2.0*qu[1]*qu[2],qu[1]**2-qu[2]**2), np.pi, 0.0])
else:
eu = np.array([np.arctan2((-P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]-qu[2]*qu[3])*chi ),
np.arctan2( 2.0*chi, q03-q12 ),
np.arctan2(( P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]+qu[2]*qu[3])*chi )])
# reduce Euler angles to definition range, i.e a lower limit of 0.0
eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu)
return eu
@staticmethod
def qu2ax(qu):
"""
Quaternion to axis angle pair.
Modified version of the original formulation, should be numerically more stable
"""
if iszero(qu[1]**2+qu[2]**2+qu[3]**2): # set axis to [001] if the angle is 0/360
ax = [ 0.0, 0.0, 1.0, 0.0 ]
elif not iszero(qu[0]):
s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2)
omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0))
ax = [ qu[1]*s, qu[2]*s, qu[3]*s, omega ]
else:
ax = [ qu[1], qu[2], qu[3], np.pi]
return np.array(ax)
@staticmethod
def qu2ro(qu):
"""Quaternion to Rodriques-Frank vector."""
if iszero(qu[0]):
ro = [qu[1], qu[2], qu[3], np.inf]
else:
s = np.linalg.norm([qu[1],qu[2],qu[3]])
ro = [0.0,0.0,P,0.0] if iszero(s) else \
[ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(np.clip(qu[0],-1.0,1.0)))]
return np.array(ro)
@staticmethod
def qu2ho(qu):
"""Quaternion to homochoric vector."""
omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0))
if iszero(omega):
ho = np.array([ 0.0, 0.0, 0.0 ])
else:
ho = np.array([qu[1], qu[2], qu[3]])
f = 0.75 * ( omega - np.sin(omega) )
ho = ho/np.linalg.norm(ho) * f**(1./3.)
return ho
@staticmethod
def qu2cu(qu):
"""Quaternion to cubochoric vector."""
return Rotation.ho2cu(Rotation.qu2ho(qu))
#---------- Rotation matrix ----------
@staticmethod
def om2qu(om):
"""
Rotation matrix to quaternion.
The original formulation (direct conversion) had (numerical?) issues
"""
return Rotation.eu2qu(Rotation.om2eu(om))
@staticmethod
def om2eu(om):
"""Rotation matrix to Bunge-Euler angles."""
if abs(om[2,2]) < 1.0:
zeta = 1.0/np.sqrt(1.0-om[2,2]**2)
eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta),
np.arccos(om[2,2]),
np.arctan2(om[0,2]*zeta, om[1,2]*zeta)])
else:
eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation
# reduce Euler angles to definition range, i.e a lower limit of 0.0
eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu)
return eu
@staticmethod
def om2ax(om):
"""Rotation matrix to axis angle pair."""
ax=np.empty(4)
# first get the rotation angle
t = 0.5*(om.trace() -1.0)
ax[3] = np.arccos(np.clip(t,-1.0,1.0))
if iszero(ax[3]):
ax = [ 0.0, 0.0, 1.0, 0.0]
else:
w,vr = np.linalg.eig(om)
# next, find the eigenvalue (1,0j)
i = np.where(np.isclose(w,1.0+0.0j))[0][0]
ax[0:3] = np.real(vr[0:3,i])
diagDelta = np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]])
ax[0:3] = np.where(iszero(diagDelta), ax[0:3],np.abs(ax[0:3])*np.sign(-P*diagDelta))
return np.array(ax)
@staticmethod
def om2ro(om):
"""Rotation matrix to Rodriques-Frank vector."""
return Rotation.eu2ro(Rotation.om2eu(om))
@staticmethod
def om2ho(om):
"""Rotation matrix to homochoric vector."""
return Rotation.ax2ho(Rotation.om2ax(om))
@staticmethod
def om2cu(om):
"""Rotation matrix to cubochoric vector."""
return Rotation.ho2cu(Rotation.om2ho(om))
#---------- Bunge-Euler angles ----------
@staticmethod
def eu2qu(eu):
"""Bunge-Euler angles to quaternion."""
ee = 0.5*eu
cPhi = np.cos(ee[1])
sPhi = np.sin(ee[1])
qu = np.array([ cPhi*np.cos(ee[0]+ee[2]),
-P*sPhi*np.cos(ee[0]-ee[2]),
-P*sPhi*np.sin(ee[0]-ee[2]),
-P*cPhi*np.sin(ee[0]+ee[2]) ])
if qu[0] < 0.0: qu*=-1
return qu
@staticmethod
def eu2om(eu):
"""Bunge-Euler angles to rotation matrix."""
c = np.cos(eu)
s = np.sin(eu)
om = np.array([[+c[0]*c[2]-s[0]*s[2]*c[1], +s[0]*c[2]+c[0]*s[2]*c[1], +s[2]*s[1]],
[-c[0]*s[2]-s[0]*c[2]*c[1], -s[0]*s[2]+c[0]*c[2]*c[1], +c[2]*s[1]],
[+s[0]*s[1], -c[0]*s[1], +c[1] ]])
om[np.where(iszero(om))] = 0.0
return om
@staticmethod
def eu2ax(eu):
"""Bunge-Euler angles to axis angle pair."""
t = np.tan(eu[1]*0.5)
sigma = 0.5*(eu[0]+eu[2])
delta = 0.5*(eu[0]-eu[2])
tau = np.linalg.norm([t,np.sin(sigma)])
alpha = np.pi if iszero(np.cos(sigma)) else \
2.0*np.arctan(tau/np.cos(sigma))
if iszero(alpha):
ax = np.array([ 0.0, 0.0, 1.0, 0.0 ])
else:
ax = -P/tau * np.array([ t*np.cos(delta), t*np.sin(delta), np.sin(sigma) ]) # passive axis angle pair so a minus sign in front
ax = np.append(ax,alpha)
if alpha < 0.0: ax *= -1.0 # ensure alpha is positive
return ax
@staticmethod
def eu2ro(eu):
"""Bunge-Euler angles to Rodriques-Frank vector."""
ro = Rotation.eu2ax(eu) # convert to axis angle pair representation
if ro[3] >= np.pi: # Differs from original implementation. check convention 5
ro[3] = np.inf
elif iszero(ro[3]):
ro = np.array([ 0.0, 0.0, P, 0.0 ])
else:
ro[3] = np.tan(ro[3]*0.5)
return ro
@staticmethod
def eu2ho(eu):
"""Bunge-Euler angles to homochoric vector."""
return Rotation.ax2ho(Rotation.eu2ax(eu))
@staticmethod
def eu2cu(eu):
"""Bunge-Euler angles to cubochoric vector."""
return Rotation.ho2cu(Rotation.eu2ho(eu))
#---------- Axis angle pair ----------
@staticmethod
def ax2qu(ax):
"""Axis angle pair to quaternion."""
if iszero(ax[3]):
qu = np.array([ 1.0, 0.0, 0.0, 0.0 ])
else:
c = np.cos(ax[3]*0.5)
s = np.sin(ax[3]*0.5)
qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ])
return qu
@staticmethod
def ax2om(ax):
"""Axis angle pair to rotation matrix."""
c = np.cos(ax[3])
s = np.sin(ax[3])
omc = 1.0-c
om=np.diag(ax[0:3]**2*omc + c)
for idx in [[0,1,2],[1,2,0],[2,0,1]]:
q = omc*ax[idx[0]] * ax[idx[1]]
om[idx[0],idx[1]] = q + s*ax[idx[2]]
om[idx[1],idx[0]] = q - s*ax[idx[2]]
return om if P < 0.0 else om.T
@staticmethod
def ax2eu(ax):
"""Rotation matrix to Bunge Euler angles."""
return Rotation.om2eu(Rotation.ax2om(ax))
@staticmethod
def ax2ro(ax):
"""Axis angle pair to Rodriques-Frank vector."""
if iszero(ax[3]):
ro = [ 0.0, 0.0, P, 0.0 ]
else:
ro = [ax[0], ax[1], ax[2]]
# 180 degree case
ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \
[np.tan(ax[3]*0.5)]
return np.array(ro)
@staticmethod
def ax2ho(ax):
"""Axis angle pair to homochoric vector."""
f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0)
ho = ax[0:3] * f
return ho
@staticmethod
def ax2cu(ax):
"""Axis angle pair to cubochoric vector."""
return Rotation.ho2cu(Rotation.ax2ho(ax))
#---------- Rodrigues-Frank vector ----------
@staticmethod
def ro2qu(ro):
"""Rodriques-Frank vector to quaternion."""
return Rotation.ax2qu(Rotation.ro2ax(ro))
@staticmethod
def ro2om(ro):
"""Rodgrigues-Frank vector to rotation matrix."""
return Rotation.ax2om(Rotation.ro2ax(ro))
@staticmethod
def ro2eu(ro):
"""Rodriques-Frank vector to Bunge-Euler angles."""
return Rotation.om2eu(Rotation.ro2om(ro))
@staticmethod
def ro2ax(ro):
"""Rodriques-Frank vector to axis angle pair."""
ta = ro[3]
if iszero(ta):
ax = [ 0.0, 0.0, 1.0, 0.0 ]
elif not np.isfinite(ta):
ax = [ ro[0], ro[1], ro[2], np.pi ]
else:
angle = 2.0*np.arctan(ta)
ta = 1.0/np.linalg.norm(ro[0:3])
ax = [ ro[0]/ta, ro[1]/ta, ro[2]/ta, angle ]
return np.array(ax)
@staticmethod
def ro2ho(ro):
"""Rodriques-Frank vector to homochoric vector."""
if iszero(np.sum(ro[0:3]**2.0)):
ho = [ 0.0, 0.0, 0.0 ]
else:
f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi
ho = ro[0:3] * (0.75*f)**(1.0/3.0)
return np.array(ho)
@staticmethod
def ro2cu(ro):
"""Rodriques-Frank vector to cubochoric vector."""
return Rotation.ho2cu(Rotation.ro2ho(ro))
#---------- Homochoric vector----------
@staticmethod
def ho2qu(ho):
"""Homochoric vector to quaternion."""
return Rotation.ax2qu(Rotation.ho2ax(ho))
@staticmethod
def ho2om(ho):
"""Homochoric vector to rotation matrix."""
return Rotation.ax2om(Rotation.ho2ax(ho))
@staticmethod
def ho2eu(ho):
"""Homochoric vector to Bunge-Euler angles."""
return Rotation.ax2eu(Rotation.ho2ax(ho))
@staticmethod
def ho2ax(ho):
"""Homochoric vector to axis angle pair."""
tfit = np.array([+1.0000000000018852, -0.5000000002194847,
-0.024999992127593126, -0.003928701544781374,
-0.0008152701535450438, -0.0002009500426119712,
-0.00002397986776071756, -0.00008202868926605841,
+0.00012448715042090092, -0.0001749114214822577,
+0.0001703481934140054, -0.00012062065004116828,
+0.000059719705868660826, -0.00001980756723965647,
+0.000003953714684212874, -0.00000036555001439719544])
# normalize h and store the magnitude
hmag_squared = np.sum(ho**2.)
if iszero(hmag_squared):
ax = np.array([ 0.0, 0.0, 1.0, 0.0 ])
else:
hm = hmag_squared
# convert the magnitude to the rotation angle
s = tfit[0] + tfit[1] * hmag_squared
for i in range(2,16):
hm *= hmag_squared
s += tfit[i] * hm
ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0)))
return ax
@staticmethod
def ho2ro(ho):
"""Axis angle pair to Rodriques-Frank vector."""
return Rotation.ax2ro(Rotation.ho2ax(ho))
@staticmethod
def ho2cu(ho):
"""Homochoric vector to cubochoric vector."""
return Lambert.BallToCube(ho)
#---------- Cubochoric ----------
@staticmethod
def cu2qu(cu):
"""Cubochoric vector to quaternion."""
return Rotation.ho2qu(Rotation.cu2ho(cu))
@staticmethod
def cu2om(cu):
"""Cubochoric vector to rotation matrix."""
return Rotation.ho2om(Rotation.cu2ho(cu))
@staticmethod
def cu2eu(cu):
"""Cubochoric vector to Bunge-Euler angles."""
return Rotation.ho2eu(Rotation.cu2ho(cu))
@staticmethod
def cu2ax(cu):
"""Cubochoric vector to axis angle pair."""
return Rotation.ho2ax(Rotation.cu2ho(cu))
@staticmethod
def cu2ro(cu):
"""Cubochoric vector to Rodriques-Frank vector."""
return Rotation.ho2ro(Rotation.cu2ho(cu))
@staticmethod
def cu2ho(cu):
"""Cubochoric vector to homochoric vector."""
return Lambert.CubeToBall(cu)

6
python/damask/table.py
View File

@ -32,7 +32,7 @@ class Table():
"""Label data individually, e.g. v v v ==> 1_v 2_v 3_v.""" """Label data individually, e.g. v v v ==> 1_v 2_v 3_v."""
labels = [] labels = []
for label,shape in self.shapes.items(): for label,shape in self.shapes.items():
size = np.prod(shape) size = int(np.prod(shape))
labels += ['{}{}'.format('' if size == 1 else '{}_'.format(i+1),label) for i in range(size)] labels += ['{}{}'.format('' if size == 1 else '{}_'.format(i+1),label) for i in range(size)]
self.data.columns = labels self.data.columns = labels
@ -41,14 +41,14 @@ class Table():
"""Label data condensed, e.g. 1_v 2_v 3_v ==> v v v.""" """Label data condensed, e.g. 1_v 2_v 3_v ==> v v v."""
labels = [] labels = []
for label,shape in self.shapes.items(): for label,shape in self.shapes.items():
labels += [label] * np.prod(shape) labels += [label] * int(np.prod(shape))
self.data.columns = labels self.data.columns = labels
def __add_comment(self,label,shape,info): def __add_comment(self,label,shape,info):
if info is not None: if info is not None:
self.comments.append('{}{}: {}'.format(label, self.comments.append('{}{}: {}'.format(label,
' '+str(shape) if np.prod(shape) > 1 else '', ' '+str(shape) if np.prod(shape,dtype=int) > 1 else '',
info)) info))

344
python/damask/util.py
View File

@ -6,8 +6,6 @@ import shlex
from fractions import Fraction from fractions import Fraction
from functools import reduce from functools import reduce
from optparse import Option from optparse import Option
from queue import Queue
from threading import Thread
import numpy as np import numpy as np
@ -40,218 +38,212 @@ class bcolors:
self.BOLD = '' self.BOLD = ''
self.UNDERLINE = '' self.UNDERLINE = ''
self.CROSSOUT = '' self.CROSSOUT = ''
# -----------------------------
def srepr(arg,glue = '\n'): def srepr(arg,glue = '\n'):
"""Joins arguments as individual lines.""" r"""
if (not hasattr(arg, "strip") and Join arguments as individual lines.
(hasattr(arg, "__getitem__") or
hasattr(arg, "__iter__"))): Parameters
return glue.join(str(x) for x in arg) ----------
return arg if isinstance(arg,str) else repr(arg) arg : iterable
Items to join.
glue : str, optional
Defaults to \n.
"""
if (not hasattr(arg, "strip") and
(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)
# -----------------------------
def croak(what, newline = True): def croak(what, newline = True):
"""Writes formated to stderr.""" """
if what is not None: Write formated to stderr.
sys.stderr.write(srepr(what,glue = '\n') + ('\n' if newline else ''))
sys.stderr.flush() Parameters
----------
what : str or iterable
Content to be displayed
newline : bool, optional
Separate items of what by newline. Defaults to True.
"""
if not what:
sys.stderr.write(srepr(what,glue = '\n') + ('\n' if newline else ''))
sys.stderr.flush()
# -----------------------------
def report(who = None, def report(who = None,
what = None): what = None):
"""Reports script and file name.""" """
croak( (emph(who)+': ' if who is not None else '') + (what if what is not None else '') + '\n' ) Reports script and file name.
DEPRECATED
"""
croak( (emph(who)+': ' if who is not None else '') + (what if what is not None else '') + '\n' )
# -----------------------------
def emph(what): def emph(what):
"""Formats string with emphasis.""" """Formats string with emphasis."""
return bcolors.BOLD+srepr(what)+bcolors.ENDC return bcolors.BOLD+srepr(what)+bcolors.ENDC
# -----------------------------
def deemph(what): def deemph(what):
"""Formats string with deemphasis.""" """Formats string with deemphasis."""
return bcolors.DIM+srepr(what)+bcolors.ENDC return bcolors.DIM+srepr(what)+bcolors.ENDC
# -----------------------------
def delete(what): def delete(what):
"""Formats string as deleted.""" """Formats string as deleted."""
return bcolors.DIM+srepr(what)+bcolors.ENDC return bcolors.DIM+srepr(what)+bcolors.ENDC
# -----------------------------
def strikeout(what): def strikeout(what):
"""Formats string as strikeout.""" """Formats string as strikeout."""
return bcolors.CROSSOUT+srepr(what)+bcolors.ENDC return bcolors.CROSSOUT+srepr(what)+bcolors.ENDC
# -----------------------------
def execute(cmd, def execute(cmd,
streamIn = None, streamIn = None,
wd = './', wd = './',
env = None): env = None):
"""Executes a command in given directory and returns stdout and stderr for optional stdin.""" """
initialPath = os.getcwd() Execute command.
os.chdir(wd)
myEnv = os.environ if env is None else env Parameters
process = subprocess.Popen(shlex.split(cmd), ----------
stdout = subprocess.PIPE, cmd : str
stderr = subprocess.PIPE, Command to be executed.
stdin = subprocess.PIPE, streanIn :, optional
env = myEnv) Input (via pipe) for executed process.
out,error = [i for i in (process.communicate() if streamIn is None wd : str, optional
else process.communicate(streamIn.read().encode('utf-8')))] Working directory of process. Defaults to ./ .
out = out.decode('utf-8').replace('\x08','') env :
error = error.decode('utf-8').replace('\x08','') Environment
os.chdir(initialPath)
if process.returncode != 0: raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode)) """
return out,error initialPath = os.getcwd()
os.chdir(wd)
myEnv = os.environ if env is None else env
process = subprocess.Popen(shlex.split(cmd),
stdout = subprocess.PIPE,
stderr = subprocess.PIPE,
stdin = subprocess.PIPE,
env = myEnv)
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
# -----------------------------
class extendableOption(Option): class extendableOption(Option):
""" """
Used for definition of new option parser action 'extend', which enables to take multiple option arguments. Used for definition of new option parser action 'extend', which enables to take multiple option arguments.
Adopted from online tutorial http://docs.python.org/library/optparse.html Adopted from online tutorial http://docs.python.org/library/optparse.html
""" DEPRECATED
"""
ACTIONS = Option.ACTIONS + ("extend",) ACTIONS = Option.ACTIONS + ("extend",)
STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",) STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",)
TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",) TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",)
ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",) ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",)
def take_action(self, action, dest, opt, value, values, parser):
if action == "extend":
lvalue = value.split(",")
values.ensure_value(dest, []).extend(lvalue)
else:
Option.take_action(self, action, dest, opt, value, values, parser)
def take_action(self, action, dest, opt, value, values, parser):
if action == "extend":
lvalue = value.split(",")
values.ensure_value(dest, []).extend(lvalue)
else:
Option.take_action(self, action, dest, opt, value, values, parser)
# Print iterations progress
# from https://gist.github.com/aubricus/f91fb55dc6ba5557fbab06119420dd6a
def progressBar(iteration, total, prefix='', bar_length=50): def progressBar(iteration, total, prefix='', bar_length=50):
""" """
Call in a loop to create terminal progress bar. Call in a loop to create terminal progress bar.
From https://gist.github.com/aubricus/f91fb55dc6ba5557fbab06119420dd6a
@params: Parameters
iteration - Required : current iteration (Int) ----------
total - Required : total iterations (Int) iteration : int
prefix - Optional : prefix string (Str) Current iteration.
bar_length - Optional : character length of bar (Int) total : int
""" Total iterations.
fraction = iteration / float(total) prefix : str, optional
if not hasattr(progressBar, "last_fraction"): # first call to function Prefix string.
progressBar.start_time = time.time() bar_length : int, optional
progressBar.last_fraction = -1.0 Character length of bar. Defaults to 50.
remaining_time = ' n/a'
else: """
if fraction <= progressBar.last_fraction or iteration == 0: # reset: called within a new loop fraction = iteration / float(total)
progressBar.start_time = time.time() if not hasattr(progressBar, "last_fraction"): # first call to function
progressBar.last_fraction = -1.0 progressBar.start_time = time.time()
remaining_time = ' n/a' progressBar.last_fraction = -1.0
remaining_time = ' n/a'
else: else:
progressBar.last_fraction = fraction if fraction <= progressBar.last_fraction or iteration == 0: # reset: called within a new loop
remainder = (total - iteration) * (time.time()-progressBar.start_time)/iteration progressBar.start_time = time.time()
remaining_time = '{: 3d}:'.format(int( remainder//3600)) + \ progressBar.last_fraction = -1.0
'{:02d}:'.format(int((remainder//60)%60)) + \ remaining_time = ' n/a'
'{:02d}' .format(int( remainder %60)) else:
progressBar.last_fraction = fraction
remainder = (total - iteration) * (time.time()-progressBar.start_time)/iteration
remaining_time = '{: 3d}:'.format(int( remainder//3600)) + \
'{:02d}:'.format(int((remainder//60)%60)) + \
'{:02d}' .format(int( remainder %60))
filled_length = int(round(bar_length * fraction)) filled_length = int(round(bar_length * fraction))
bar = '█' * filled_length + '░' * (bar_length - filled_length) bar = '█' * filled_length + '░' * (bar_length - filled_length)
sys.stderr.write('\r{} {} {}'.format(prefix, bar, remaining_time)), sys.stderr.write('\r{} {} {}'.format(prefix, bar, remaining_time)),
if iteration == total: sys.stderr.write('\n')
sys.stderr.flush()
if iteration == total:
sys.stderr.write('\n')
sys.stderr.flush()
def scale_to_coprime(v): def scale_to_coprime(v):
"""Scale vector to co-prime (relatively prime) integers.""" """Scale vector to co-prime (relatively prime) integers."""
MAX_DENOMINATOR = 1000 MAX_DENOMINATOR = 1000
def get_square_denominator(x): def get_square_denominator(x):
"""Denominator of the square of a number.""" """Denominator of the square of a number."""
return Fraction(x ** 2).limit_denominator(MAX_DENOMINATOR).denominator return Fraction(x ** 2).limit_denominator(MAX_DENOMINATOR).denominator
def lcm(a, b): def lcm(a, b):
"""Least common multiple.""" """Least common multiple."""
return a * b // np.gcd(a, b) return a * b // np.gcd(a, b)
denominators = [int(get_square_denominator(i)) for i in v] denominators = [int(get_square_denominator(i)) for i in v]
s = reduce(lcm, denominators) ** 0.5 s = reduce(lcm, denominators) ** 0.5
m = (np.array(v)*s).astype(np.int) m = (np.array(v)*s).astype(np.int)
return m//reduce(np.gcd,m) return m//reduce(np.gcd,m)
class return_message(): class return_message():
"""Object with formatted return message.""" """Object with formatted return message."""
def __init__(self,message):
"""
Sets return message.
def __init__(self,message): Parameters
""" ----------
Sets return message. message : str or list of str
message for output to screen
Parameters """
---------- self.message = message
message : str or list of str
message for output to screen def __repr__(self):
"""Return message suitable for interactive shells."""
return srepr(self.message)
"""
self.message = message
def __repr__(self):
"""Return message suitable for interactive shells."""
return srepr(self.message)
class ThreadPool:
"""Pool of threads consuming tasks from a queue."""
class Worker(Thread):
"""Thread executing tasks from a given tasks queue."""
def __init__(self, tasks):
"""Worker for tasks."""
Thread.__init__(self)
self.tasks = tasks
self.daemon = True
self.start()
def run(self):
while True:
func, args, kargs = self.tasks.get()
try:
func(*args, **kargs)
except Exception as e:
# An exception happened in this thread
print(e)
finally:
# Mark this task as done, whether an exception happened or not
self.tasks.task_done()
def __init__(self, num_threads):
"""
Thread pool.
Parameters
----------
num_threads : int
number of threads
"""
self.tasks = Queue(num_threads)
for _ in range(num_threads):
self.Worker(self.tasks)
def add_task(self, func, *args, **kargs):
"""Add a task to the queue."""
self.tasks.put((func, args, kargs))
def map(self, func, args_list):
"""Add a list of tasks to the queue."""
for args in args_list:
self.add_task(func, args)
def wait_completion(self):
"""Wait for completion of all the tasks in the queue."""
self.tasks.join()

0
python/tests/reference/DADF5/12grains6x7x8.geom → python/tests/reference/Result/12grains6x7x8.geom
View File

0
python/tests/reference/DADF5/12grains6x7x8_tensionY.hdf5 → python/tests/reference/Result/12grains6x7x8_tensionY.hdf5
View File

0
python/tests/reference/DADF5/material.config → python/tests/reference/Result/material.config
View File

0
python/tests/reference/DADF5/tensionY.load → python/tests/reference/Result/tensionY.load
View File

65
python/tests/test_Orientation.py Normal file
View File

@ -0,0 +1,65 @@
import os
from itertools import permutations
import pytest
import numpy as np
import damask
from damask import Rotation
from damask import Orientation
from damask import Lattice
n = 1000
@pytest.fixture
def default():
"""A set of n random rotations."""
return [Rotation.fromRandom() for r in range(n)]
@pytest.fixture
def reference_dir(reference_dir_base):
"""Directory containing reference results."""
return os.path.join(reference_dir_base,'Rotation')
class TestOrientation:
@pytest.mark.parametrize('color',[{'label':'red', 'RGB':[1,0,0],'direction':[0,0,1]},
{'label':'green','RGB':[0,1,0],'direction':[0,1,1]},
{'label':'blue', 'RGB':[0,0,1],'direction':[1,1,1]}])
@pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_IPF_cubic(self,default,color,lattice):
cube = damask.Orientation(damask.Rotation(),lattice)
for direction in set(permutations(np.array(color['direction']))):
assert np.allclose(cube.IPFcolor(direction),np.array(color['RGB']))
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_IPF(self,lattice):
direction = np.random.random(3)*2.0-1
for rot in [Rotation.fromRandom() for r in range(n//100)]:
R = damask.Orientation(rot,lattice)
color = R.IPFcolor(direction)
for equivalent in R.equivalentOrientations():
assert np.allclose(color,R.IPFcolor(direction))
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_relationship_forward_backward(self,model,lattice):
ori = Orientation(Rotation.fromRandom(),lattice)
for i,r in enumerate(ori.relatedOrientations(model)):
ori2 = r.relatedOrientations(model)[i]
misorientation = ori.rotation.misorientation(ori2.rotation)
assert misorientation.asAxisAngle(degrees=True)[3]<1.0e-5
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_relationship_reference(self,update,reference_dir,model,lattice):
reference = os.path.join(reference_dir,'{}_{}.txt'.format(lattice,model))
ori = Orientation(Rotation(),lattice)
eu = np.array([o.rotation.asEulers(degrees=True) for o in ori.relatedOrientations(model)])
if update:
coords = np.array([(1,i+1) for i,x in enumerate(eu)])
table = damask.Table(eu,{'Eulers':(3,)})
table.add('pos',coords)
table.to_ASCII(reference)
assert np.allclose(eu,damask.Table.from_ASCII(reference).get('Eulers'))

10
python/tests/test_DADF5.py → python/tests/test_Result.py
View File

@ -4,25 +4,25 @@ import os
import pytest import pytest
import numpy as np import numpy as np
from damask import DADF5 from damask import Result
from damask import mechanics from damask import mechanics
@pytest.fixture @pytest.fixture
def default(tmp_path,reference_dir): def default(tmp_path,reference_dir):
"""Small DADF5 file in temp location for modification.""" """Small Result file in temp location for modification."""
fname = '12grains6x7x8_tensionY.hdf5' fname = '12grains6x7x8_tensionY.hdf5'
shutil.copy(os.path.join(reference_dir,fname),tmp_path) shutil.copy(os.path.join(reference_dir,fname),tmp_path)
f = DADF5(os.path.join(tmp_path,fname)) f = Result(os.path.join(tmp_path,fname))
f.set_by_time(20.0,20.0) f.set_by_time(20.0,20.0)
return f return f
@pytest.fixture @pytest.fixture
def reference_dir(reference_dir_base): def reference_dir(reference_dir_base):
"""Directory containing reference results.""" """Directory containing reference results."""
return os.path.join(reference_dir_base,'DADF5') return os.path.join(reference_dir_base,'Result')
class TestDADF5: class TestResult:
def test_time_increments(self,default): def test_time_increments(self,default):
shape = default.read_dataset(default.get_dataset_location('F'),0).shape shape = default.read_dataset(default.get_dataset_location('F'),0).shape

45
python/tests/test_Rotation.py
View File

@ -1,13 +1,9 @@
import os import os
from itertools import permutations
import pytest import pytest
import numpy as np import numpy as np
import damask
from damask import Rotation from damask import Rotation
from damask import Orientation
from damask import Lattice
n = 1000 n = 1000
@ -58,44 +54,3 @@ class TestRotation:
for rot in default: for rot in default:
assert np.allclose(rot.asCubochoric(), assert np.allclose(rot.asCubochoric(),
Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric()) Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric())
@pytest.mark.parametrize('color',[{'label':'red', 'RGB':[1,0,0],'direction':[0,0,1]},
{'label':'green','RGB':[0,1,0],'direction':[0,1,1]},
{'label':'blue', 'RGB':[0,0,1],'direction':[1,1,1]}])
@pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_IPF_cubic(self,default,color,lattice):
cube = damask.Orientation(damask.Rotation(),lattice)
for direction in set(permutations(np.array(color['direction']))):
assert np.allclose(cube.IPFcolor(direction),np.array(color['RGB']))
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_IPF(self,lattice):
direction = np.random.random(3)*2.0-1
for rot in [Rotation.fromRandom() for r in range(n//100)]:
R = damask.Orientation(rot,lattice)
color = R.IPFcolor(direction)
for equivalent in R.equivalentOrientations():
assert np.allclose(color,R.IPFcolor(direction))
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_relationship_forward_backward(self,model,lattice):
ori = Orientation(Rotation.fromRandom(),lattice)
for i,r in enumerate(ori.relatedOrientations(model)):
ori2 = r.relatedOrientations(model)[i]
misorientation = ori.rotation.misorientation(ori2.rotation)
assert misorientation.asAxisAngle(degrees=True)[3]<1.0e-5
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_relationship_reference(self,update,reference_dir,model,lattice):
reference = os.path.join(reference_dir,'{}_{}.txt'.format(lattice,model))
ori = Orientation(Rotation(),lattice)
eu = np.array([o.rotation.asEulers(degrees=True) for o in ori.relatedOrientations(model)])
if update:
coords = np.array([(1,i+1) for i,x in enumerate(eu)])
table = damask.Table(eu,{'Eulers':(3,)})
table.add('pos',coords)
table.to_ASCII(reference)
assert np.allclose(eu,damask.Table.from_ASCII(reference).get('Eulers'))

91
src/CPFEM.f90
View File

@ -69,31 +69,32 @@ contains
!> @brief call (thread safe) all module initializations !> @brief call (thread safe) all module initializations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_initAll(el,ip) subroutine CPFEM_initAll(el,ip)
integer(pInt), intent(in) :: el, & !< FE el number
ip !< FE integration point number integer(pInt), intent(in) :: el, & !< FE el number
ip !< FE integration point number
!$OMP CRITICAL (init) !$OMP CRITICAL(init)
if (.not. CPFEM_init_done) then if (.not. CPFEM_init_done) then
call DAMASK_interface_init ! Spectral and FEM interface to commandline call DAMASK_interface_init
call prec_init call prec_init
call IO_init call IO_init
call numerics_init call numerics_init
call debug_init call debug_init
call config_init call config_init
call math_init call math_init
call rotations_init call rotations_init
call HDF5_utilities_init call HDF5_utilities_init
call results_init call results_init
call mesh_init(ip, el) call mesh_init(ip, el)
call lattice_init call lattice_init
call material_init call material_init
call constitutive_init call constitutive_init
call crystallite_init call crystallite_init
call homogenization_init call homogenization_init
call CPFEM_init call CPFEM_init
CPFEM_init_done = .true. CPFEM_init_done = .true.
endif endif
!$OMP END CRITICAL (init) !$OMP END CRITICAL(init)
end subroutine CPFEM_initAll end subroutine CPFEM_initAll
@ -174,35 +175,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
CPFEM_dcsde = CPFEM_dcsde_knownGood CPFEM_dcsde = CPFEM_dcsde_knownGood
!*** age results !*** age results
if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) then if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward
crystallite_F0 = crystallite_partionedF ! crystallite deformation
crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation
crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity
crystallite_Fi0 = crystallite_Fi ! crystallite intermediate deformation
crystallite_Li0 = crystallite_Li ! crystallite intermediate velocity
crystallite_S0 = crystallite_S ! crystallite 2nd Piola Kirchhoff stress
forall (i = 1:size(plasticState)) plasticState(i)%state0 = plasticState(i)%state
do i = 1, size(sourceState)
do mySource = 1,phase_Nsources(i)
sourceState(i)%p(mySource)%state0 = sourceState(i)%p(mySource)%state
enddo; enddo
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) then
write(6,'(a)') '<< CPFEM >> aging states'
if (debug_e <= discretization_nElem .and. debug_i <=discretization_nIP) then
write(6,'(a,1x,i8,1x,i2,1x,i4,/,(12x,6(e20.8,1x)),/)') &
'<< CPFEM >> aged state of elFE ip grain',debug_e, debug_i, 1, &
plasticState(material_phaseAt(1,debug_e))%state(:,material_phasememberAt(1,debug_i,debug_e))
endif
endif
do homog = 1_pInt, material_Nhomogenization
homogState (homog)%state0 = homogState (homog)%state
thermalState (homog)%state0 = thermalState (homog)%state
damageState (homog)%state0 = damageState (homog)%state
enddo
endif
!*** collection of FEM input with returning of randomize odd stress and jacobian !*** collection of FEM input with returning of randomize odd stress and jacobian
@ -360,7 +333,17 @@ end subroutine CPFEM_general
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief triggers writing of the results !> @brief Forward data for new time increment.
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_forward
call crystallite_forward
end subroutine CPFEM_forward
!--------------------------------------------------------------------------------------------------
!> @brief Trigger writing of results.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_results(inc,time) subroutine CPFEM_results(inc,time)

118
src/CPFEM2.f90
View File

@ -64,131 +64,45 @@ end subroutine CPFEM_initAll
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocate the arrays defined in module CPFEM and initialize them !> @brief Read restart information if needed.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_init subroutine CPFEM_init
integer :: i
integer(HID_T) :: fileHandle, groupHandle
character(len=pStringLen) :: fileName, datasetName
write(6,'(/,a)') ' <<<+- CPFEM init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- CPFEM init -+>>>'; flush(6)
if (interface_restartInc > 0) then if (interface_restartInc > 0) call crystallite_restartRead
write(6,'(/,a,i0,a)') ' reading restart information of increment ', interface_restartInc, ' from file'
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName)
call HDF5_read(fileHandle,crystallite_F0, 'F')
call HDF5_read(fileHandle,crystallite_Fp0,'Fp')
call HDF5_read(fileHandle,crystallite_Fi0,'Fi')
call HDF5_read(fileHandle,crystallite_Lp0,'Lp')
call HDF5_read(fileHandle,crystallite_Li0,'Li')
call HDF5_read(fileHandle,crystallite_S0, 'S')
groupHandle = HDF5_openGroup(fileHandle,'constituent')
do i = 1,size(phase_plasticity)
write(datasetName,'(i0,a)') i,'_omega_plastic'
call HDF5_read(groupHandle,plasticState(i)%state0,datasetName)
enddo
call HDF5_closeGroup(groupHandle)
groupHandle = HDF5_openGroup(fileHandle,'materialpoint')
do i = 1, material_Nhomogenization
write(datasetName,'(i0,a)') i,'_omega_homogenization'
call HDF5_read(groupHandle,homogState(i)%state0,datasetName)
enddo
call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle)
endif
end subroutine CPFEM_init end subroutine CPFEM_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Forward data after successful increment. !> @brief Write restart information.
! ToDo: Any guessing for the current states possible?
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_forward
integer :: i, j
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) &
write(6,'(a)') '<< CPFEM >> aging states'
crystallite_F0 = crystallite_partionedF
crystallite_Fp0 = crystallite_Fp
crystallite_Lp0 = crystallite_Lp
crystallite_Fi0 = crystallite_Fi
crystallite_Li0 = crystallite_Li
crystallite_S0 = crystallite_S
do i = 1, size(plasticState)
plasticState(i)%state0 = plasticState(i)%state
enddo
do i = 1, size(sourceState)
do j = 1,phase_Nsources(i)
sourceState(i)%p(j)%state0 = sourceState(i)%p(j)%state
enddo; enddo
do i = 1, material_Nhomogenization
homogState (i)%state0 = homogState (i)%state
thermalState(i)%state0 = thermalState(i)%state
damageState (i)%state0 = damageState (i)%state
enddo
end subroutine CPFEM_forward
!--------------------------------------------------------------------------------------------------
!> @brief Write current restart information (Field and constitutive data) to file.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_restartWrite subroutine CPFEM_restartWrite
integer :: i call crystallite_restartWrite
integer(HID_T) :: fileHandle, groupHandle
character(len=pStringLen) :: fileName, datasetName
write(6,'(a)') ' writing field and constitutive data required for restart to file';flush(6)
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName,'a')
call HDF5_write(fileHandle,crystallite_partionedF,'F')
call HDF5_write(fileHandle,crystallite_Fp, 'Fp')
call HDF5_write(fileHandle,crystallite_Fi, 'Fi')
call HDF5_write(fileHandle,crystallite_Lp, 'Lp')
call HDF5_write(fileHandle,crystallite_Li, 'Li')
call HDF5_write(fileHandle,crystallite_S, 'S')
groupHandle = HDF5_addGroup(fileHandle,'constituent')
do i = 1,size(phase_plasticity)
write(datasetName,'(i0,a)') i,'_omega_plastic'
call HDF5_write(groupHandle,plasticState(i)%state,datasetName)
enddo
call HDF5_closeGroup(groupHandle)
groupHandle = HDF5_addGroup(fileHandle,'materialpoint')
do i = 1, material_Nhomogenization
write(datasetName,'(i0,a)') i,'_omega_homogenization'
call HDF5_write(groupHandle,homogState(i)%state,datasetName)
enddo
call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle)
end subroutine CPFEM_restartWrite end subroutine CPFEM_restartWrite
!--------------------------------------------------------------------------------------------------
!> @brief Forward data for new time increment.
!--------------------------------------------------------------------------------------------------
subroutine CPFEM_forward
call crystallite_forward
end subroutine CPFEM_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Trigger writing of results. !> @brief Trigger writing of results.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_results(inc,time) subroutine CPFEM_results(inc,time)
integer, intent(in) :: inc integer, intent(in) :: inc
real(pReal), intent(in) :: time real(pReal), intent(in) :: time
call results_openJobFile call results_openJobFile
call results_addIncrement(inc,time) call results_addIncrement(inc,time)
call constitutive_results call constitutive_results

18
src/DAMASK_marc.f90
View File

@ -40,9 +40,9 @@ module DAMASK_interface
implicit none implicit none
private private
logical, protected, public :: symmetricSolver logical, protected, public :: symmetricSolver
character(len=*), parameter, public :: INPUTFILEEXTENSION = '.dat' character(len=*), parameter, public :: INPUTFILEEXTENSION = '.dat'
public :: & public :: &
DAMASK_interface_init, & DAMASK_interface_init, &
getSolverJobName getSolverJobName
@ -57,14 +57,14 @@ subroutine DAMASK_interface_init
integer, dimension(8) :: dateAndTime integer, dimension(8) :: dateAndTime
integer :: ierr integer :: ierr
character(len=pPathLen) :: wd character(len=pPathLen) :: wd
write(6,'(/,a)') ' <<<+- DAMASK_marc init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_marc init -+>>>'
write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2019' write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2019'
write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2018.04.030' write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2018.04.030'
write(6,'(/,a)') ' Version: '//DAMASKVERSION write(6,'(/,a)') ' Version: '//DAMASKVERSION
! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md ! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md
#if __INTEL_COMPILER >= 1800 #if __INTEL_COMPILER >= 1800
write(6,'(/,a)') ' Compiled with: '//compiler_version() write(6,'(/,a)') ' Compiled with: '//compiler_version()
@ -73,13 +73,13 @@ subroutine DAMASK_interface_init
write(6,'(/,a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,& write(6,'(/,a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,&
', build date :', __INTEL_COMPILER_BUILD_DATE ', build date :', __INTEL_COMPILER_BUILD_DATE
#endif #endif
write(6,'(/,a)') ' Compiled on: '//__DATE__//' at '//__TIME__ write(6,'(/,a)') ' Compiled on: '//__DATE__//' at '//__TIME__
call date_and_time(values = dateAndTime) call date_and_time(values = dateAndTime)
write(6,'(/,a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',dateAndTime(2),'/', dateAndTime(1) write(6,'(/,a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',dateAndTime(2),'/', dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':', dateAndTime(6),':', dateAndTime(7) write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':', dateAndTime(6),':', dateAndTime(7)
inquire(5, name=wd) inquire(5, name=wd)
wd = wd(1:scan(wd,'/',back=.true.)) wd = wd(1:scan(wd,'/',back=.true.))
ierr = CHDIR(wd) ierr = CHDIR(wd)

188
src/crystallite.f90
View File

@ -11,6 +11,8 @@
module crystallite module crystallite
use prec use prec
use IO use IO
use HDF5_utilities
use DAMASK_interface
use config use config
use debug use debug
use numerics use numerics
@ -36,25 +38,25 @@ module crystallite
crystallite_orientation !< current orientation crystallite_orientation !< current orientation
real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: &
crystallite_Fe, & !< current "elastic" def grad (end of converged time step) crystallite_Fe, & !< current "elastic" def grad (end of converged time step)
crystallite_P !< 1st Piola-Kirchhoff stress per grain crystallite_P, & !< 1st Piola-Kirchhoff stress per grain
crystallite_S0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc
crystallite_Fp0, & !< plastic def grad at start of FE inc
crystallite_Fi0, & !< intermediate def grad at start of FE inc
crystallite_F0, & !< def grad at start of FE inc
crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc
crystallite_Li0 !< intermediate velocitiy grad at start of FE inc
real(pReal), dimension(:,:,:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:,:,:), allocatable, public :: &
crystallite_S, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step) crystallite_S, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step)
crystallite_S0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc
crystallite_partionedS0, & !< 2nd Piola-Kirchhoff stress vector at start of homog inc crystallite_partionedS0, & !< 2nd Piola-Kirchhoff stress vector at start of homog inc
crystallite_Fp, & !< current plastic def grad (end of converged time step) crystallite_Fp, & !< current plastic def grad (end of converged time step)
crystallite_Fp0, & !< plastic def grad at start of FE inc
crystallite_partionedFp0,& !< plastic def grad at start of homog inc crystallite_partionedFp0,& !< plastic def grad at start of homog inc
crystallite_Fi, & !< current intermediate def grad (end of converged time step) crystallite_Fi, & !< current intermediate def grad (end of converged time step)
crystallite_Fi0, & !< intermediate def grad at start of FE inc
crystallite_partionedFi0,& !< intermediate def grad at start of homog inc crystallite_partionedFi0,& !< intermediate def grad at start of homog inc
crystallite_F0, & !< def grad at start of FE inc
crystallite_partionedF, & !< def grad to be reached at end of homog inc crystallite_partionedF, & !< def grad to be reached at end of homog inc
crystallite_partionedF0, & !< def grad at start of homog inc crystallite_partionedF0, & !< def grad at start of homog inc
crystallite_Lp, & !< current plastic velocitiy grad (end of converged time step) crystallite_Lp, & !< current plastic velocitiy grad (end of converged time step)
crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc
crystallite_partionedLp0, & !< plastic velocity grad at start of homog inc crystallite_partionedLp0, & !< plastic velocity grad at start of homog inc
crystallite_Li, & !< current intermediate velocitiy grad (end of converged time step) crystallite_Li, & !< current intermediate velocitiy grad (end of converged time step)
crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc
crystallite_partionedLi0 !< intermediate velocity grad at start of homog inc crystallite_partionedLi0 !< intermediate velocity grad at start of homog inc
real(pReal), dimension(:,:,:,:,:), allocatable :: & real(pReal), dimension(:,:,:,:,:), allocatable :: &
crystallite_subFp0,& !< plastic def grad at start of crystallite inc crystallite_subFp0,& !< plastic def grad at start of crystallite inc
@ -104,7 +106,10 @@ module crystallite
crystallite_stressTangent, & crystallite_stressTangent, &
crystallite_orientations, & crystallite_orientations, &
crystallite_push33ToRef, & crystallite_push33ToRef, &
crystallite_results crystallite_results, &
crystallite_restartWrite, &
crystallite_restartRead, &
crystallite_forward
contains contains
@ -130,38 +135,30 @@ subroutine crystallite_init
iMax = discretization_nIP iMax = discretization_nIP
eMax = discretization_nElem eMax = discretization_nElem
allocate(crystallite_S0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_partionedF(3,3,cMax,iMax,eMax),source=0.0_pReal)
allocate(crystallite_partionedS0(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_S(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_S0, &
allocate(crystallite_P(3,3,cMax,iMax,eMax), source=0.0_pReal) crystallite_F0, crystallite_Fi0,crystallite_Fp0, &
allocate(crystallite_F0(3,3,cMax,iMax,eMax), source=0.0_pReal) crystallite_Li0,crystallite_Lp0, &
allocate(crystallite_partionedF0(3,3,cMax,iMax,eMax), source=0.0_pReal) crystallite_partionedS0, &
allocate(crystallite_partionedF(3,3,cMax,iMax,eMax), source=0.0_pReal) crystallite_partionedF0,crystallite_partionedFp0,crystallite_partionedFi0, &
allocate(crystallite_subF0(3,3,cMax,iMax,eMax), source=0.0_pReal) crystallite_partionedLp0,crystallite_partionedLi0, &
allocate(crystallite_subF(3,3,cMax,iMax,eMax), source=0.0_pReal) crystallite_S,crystallite_P, &
allocate(crystallite_Fp0(3,3,cMax,iMax,eMax), source=0.0_pReal) crystallite_Fe,crystallite_Fi,crystallite_Fp, &
allocate(crystallite_partionedFp0(3,3,cMax,iMax,eMax), source=0.0_pReal) crystallite_Li,crystallite_Lp, &
allocate(crystallite_subFp0(3,3,cMax,iMax,eMax), source=0.0_pReal) crystallite_subF,crystallite_subF0, &
allocate(crystallite_Fp(3,3,cMax,iMax,eMax), source=0.0_pReal) crystallite_subFp0,crystallite_subFi0, &
allocate(crystallite_Fi0(3,3,cMax,iMax,eMax), source=0.0_pReal) crystallite_subLi0,crystallite_subLp0, &
allocate(crystallite_partionedFi0(3,3,cMax,iMax,eMax), source=0.0_pReal) source = crystallite_partionedF)
allocate(crystallite_subFi0(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Fi(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_dPdF(3,3,3,3,cMax,iMax,eMax),source=0.0_pReal)
allocate(crystallite_Fe(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Lp0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_dt(cMax,iMax,eMax),source=0.0_pReal)
allocate(crystallite_partionedLp0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_subdt,crystallite_subFrac,crystallite_subStep, &
allocate(crystallite_subLp0(3,3,cMax,iMax,eMax), source=0.0_pReal) source = crystallite_dt)
allocate(crystallite_Lp(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Li0(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_partionedLi0(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_subLi0(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_Li(3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_dPdF(3,3,3,3,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_dt(cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_subdt(cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_subFrac(cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_subStep(cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_orientation(cMax,iMax,eMax)) allocate(crystallite_orientation(cMax,iMax,eMax))
allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.) allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.)
allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) allocate(crystallite_requested(cMax,iMax,eMax), source=.false.)
allocate(crystallite_todo(cMax,iMax,eMax), source=.false.) allocate(crystallite_todo(cMax,iMax,eMax), source=.false.)
@ -1844,4 +1841,117 @@ logical function stateJump(ipc,ip,el)
end function stateJump end function stateJump
!--------------------------------------------------------------------------------------------------
!> @brief Write current restart information (Field and constitutive data) to file.
! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively
!--------------------------------------------------------------------------------------------------
subroutine crystallite_restartWrite
integer :: i
integer(HID_T) :: fileHandle, groupHandle
character(len=pStringLen) :: fileName, datasetName
write(6,'(a)') ' writing field and constitutive data required for restart to file';flush(6)
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName,'a')
call HDF5_write(fileHandle,crystallite_partionedF,'F')
call HDF5_write(fileHandle,crystallite_Fp, 'Fp')
call HDF5_write(fileHandle,crystallite_Fi, 'Fi')
call HDF5_write(fileHandle,crystallite_Lp, 'Lp')
call HDF5_write(fileHandle,crystallite_Li, 'Li')
call HDF5_write(fileHandle,crystallite_S, 'S')
groupHandle = HDF5_addGroup(fileHandle,'constituent')
do i = 1,size(phase_plasticity)
write(datasetName,'(i0,a)') i,'_omega_plastic'
call HDF5_write(groupHandle,plasticState(i)%state,datasetName)
enddo
call HDF5_closeGroup(groupHandle)
groupHandle = HDF5_addGroup(fileHandle,'materialpoint')
do i = 1, material_Nhomogenization
write(datasetName,'(i0,a)') i,'_omega_homogenization'
call HDF5_write(groupHandle,homogState(i)%state,datasetName)
enddo
call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle)
end subroutine crystallite_restartWrite
!--------------------------------------------------------------------------------------------------
!> @brief Read data for restart
! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively
!--------------------------------------------------------------------------------------------------
subroutine crystallite_restartRead
integer :: i
integer(HID_T) :: fileHandle, groupHandle
character(len=pStringLen) :: fileName, datasetName
write(6,'(/,a,i0,a)') ' reading restart information of increment from file'
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName)
call HDF5_read(fileHandle,crystallite_F0, 'F')
call HDF5_read(fileHandle,crystallite_Fp0,'Fp')
call HDF5_read(fileHandle,crystallite_Fi0,'Fi')
call HDF5_read(fileHandle,crystallite_Lp0,'Lp')
call HDF5_read(fileHandle,crystallite_Li0,'Li')
call HDF5_read(fileHandle,crystallite_S0, 'S')
groupHandle = HDF5_openGroup(fileHandle,'constituent')
do i = 1,size(phase_plasticity)
write(datasetName,'(i0,a)') i,'_omega_plastic'
call HDF5_read(groupHandle,plasticState(i)%state0,datasetName)
enddo
call HDF5_closeGroup(groupHandle)
groupHandle = HDF5_openGroup(fileHandle,'materialpoint')
do i = 1, material_Nhomogenization
write(datasetName,'(i0,a)') i,'_omega_homogenization'
call HDF5_read(groupHandle,homogState(i)%state0,datasetName)
enddo
call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle)
end subroutine crystallite_restartRead
!--------------------------------------------------------------------------------------------------
!> @brief Forward data after successful increment.
! ToDo: Any guessing for the current states possible?
!--------------------------------------------------------------------------------------------------
subroutine crystallite_forward
integer :: i, j
crystallite_F0 = crystallite_partionedF
crystallite_Fp0 = crystallite_Fp
crystallite_Lp0 = crystallite_Lp
crystallite_Fi0 = crystallite_Fi
crystallite_Li0 = crystallite_Li
crystallite_S0 = crystallite_S
do i = 1, size(plasticState)
plasticState(i)%state0 = plasticState(i)%state
enddo
do i = 1, size(sourceState)
do j = 1,phase_Nsources(i)
sourceState(i)%p(j)%state0 = sourceState(i)%p(j)%state
enddo; enddo
do i = 1, material_Nhomogenization
homogState (i)%state0 = homogState (i)%state
thermalState(i)%state0 = thermalState(i)%state
damageState (i)%state0 = damageState (i)%state
enddo
end subroutine crystallite_forward
end module crystallite end module crystallite

4
src/homogenization_mech_RGC.f90
View File

@ -669,7 +669,7 @@ module procedure mech_RGC_updateState
nDef = 0.0_pReal nDef = 0.0_pReal
do i = 1,3; do j = 1,3 do i = 1,3; do j = 1,3
do k = 1,3; do l = 1,3 do k = 1,3; do l = 1,3
nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_civita(j,k,l) ! compute the interface mismatch tensor from the jump of deformation gradient nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_LeviCivita(j,k,l) ! compute the interface mismatch tensor from the jump of deformation gradient
enddo; enddo enddo; enddo
nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal ! compute the norm of the mismatch tensor nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal ! compute the norm of the mismatch tensor
enddo; enddo enddo; enddo
@ -689,7 +689,7 @@ module procedure mech_RGC_updateState
rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xiAlpha & rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xiAlpha &
*surfCorr(abs(intFace(1)))/prm%dAlpha(abs(intFace(1))) & *surfCorr(abs(intFace(1)))/prm%dAlpha(abs(intFace(1))) &
*cosh(prm%ciAlpha*nDefNorm) & *cosh(prm%ciAlpha*nDefNorm) &
*0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_civita(k,l,j) & *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
*tanh(nDefNorm/xSmoo_RGC) *tanh(nDefNorm/xSmoo_RGC)
enddo; enddo;enddo; enddo enddo; enddo;enddo; enddo
enddo interfaceLoop enddo interfaceLoop

186
src/math.f90
View File

@ -72,7 +72,12 @@ module math
3,2, & 3,2, &
3,3 & 3,3 &
],[2,9]) !< arrangement in Plain notation ],[2,9]) !< arrangement in Plain notation
interface math_mul33xx33
module procedure math_tensordot
end interface math_mul33xx33
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
private :: & private :: &
unitTest unitTest
@ -88,7 +93,7 @@ subroutine math_init
real(pReal), dimension(4) :: randTest real(pReal), dimension(4) :: randTest
integer :: randSize integer :: randSize
integer, dimension(:), allocatable :: randInit integer, dimension(:), allocatable :: randInit
write(6,'(/,a)') ' <<<+- math init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- math init -+>>>'; flush(6)
call random_seed(size=randSize) call random_seed(size=randSize)
@ -140,7 +145,7 @@ recursive subroutine math_sort(a, istart, iend, sortDim)
else else
e = ubound(a,2) e = ubound(a,2)
endif endif
if(present(sortDim)) then if(present(sortDim)) then
d = sortDim d = sortDim
else else
@ -160,12 +165,12 @@ recursive subroutine math_sort(a, istart, iend, sortDim)
!> @brief Partitioning required for quicksort !> @brief Partitioning required for quicksort
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
integer function qsort_partition(a, istart, iend, sort) integer function qsort_partition(a, istart, iend, sort)
integer, dimension(:,:), intent(inout) :: a integer, dimension(:,:), intent(inout) :: a
integer, intent(in) :: istart,iend,sort integer, intent(in) :: istart,iend,sort
integer, dimension(size(a,1)) :: tmp integer, dimension(size(a,1)) :: tmp
integer :: i,j integer :: i,j
do do
! find the first element on the right side less than or equal to the pivot point ! find the first element on the right side less than or equal to the pivot point
do j = iend, istart, -1 do j = iend, istart, -1
@ -187,7 +192,7 @@ recursive subroutine math_sort(a, istart, iend, sortDim)
a(:,j) = tmp a(:,j) = tmp
endif cross endif cross
enddo enddo
end function qsort_partition end function qsort_partition
end subroutine math_sort end subroutine math_sort
@ -196,7 +201,7 @@ end subroutine math_sort
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief vector expansion !> @brief vector expansion
!> @details takes a set of numbers (a,b,c,...) and corresponding multiples (x,y,z,...) !> @details takes a set of numbers (a,b,c,...) and corresponding multiples (x,y,z,...)
!> to return a vector of x times a, y times b, z times c, ... !> to return a vector of x times a, y times b, z times c, ...
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_expand(what,how) pure function math_expand(what,how)
@ -204,9 +209,9 @@ pure function math_expand(what,how)
integer, dimension(:), intent(in) :: how integer, dimension(:), intent(in) :: how
real(pReal), dimension(sum(how)) :: math_expand real(pReal), dimension(sum(how)) :: math_expand
integer :: i integer :: i
if (sum(how) == 0) return if (sum(how) == 0) return
do i = 1, size(how) do i = 1, size(how)
math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1,size(what))+1) math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1,size(what))+1)
enddo enddo
@ -266,31 +271,30 @@ end function math_identity4th
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief permutation tensor e_ijk used for computing cross product of two tensors !> @brief permutation tensor e_ijk
! e_ijk = 1 if even permutation of ijk ! e_ijk = 1 if even permutation of ijk
! e_ijk = -1 if odd permutation of ijk ! e_ijk = -1 if odd permutation of ijk
! e_ijk = 0 otherwise ! e_ijk = 0 otherwise
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
real(pReal) pure function math_civita(i,j,k) real(pReal) pure function math_LeviCivita(i,j,k)
integer, intent(in) :: i,j,k integer, intent(in) :: i,j,k
math_civita = 0.0_pReal if (all([i,j,k] == [1,2,3]) .or. all([i,j,k] == [2,3,1]) .or. all([i,j,k] == [3,1,2])) then
if (((i == 1).and.(j == 2).and.(k == 3)) .or. & math_LeviCivita = +1.0_pReal
((i == 2).and.(j == 3).and.(k == 1)) .or. & elseif (all([i,j,k] == [3,2,1]) .or. all([i,j,k] == [2,1,3]) .or. all([i,j,k] == [1,3,2])) then
((i == 3).and.(j == 1).and.(k == 2))) math_civita = 1.0_pReal math_LeviCivita = -1.0_pReal
if (((i == 1).and.(j == 3).and.(k == 2)) .or. & else
((i == 2).and.(j == 1).and.(k == 3)) .or. & math_LeviCivita = 0.0_pReal
((i == 3).and.(j == 2).and.(k == 1))) math_civita = -1.0_pReal endif
end function math_civita end function math_LeviCivita
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief kronecker delta function d_ij !> @brief kronecker delta function d_ij
! d_ij = 1 if i = j ! d_ij = 1 if i = j
! d_ij = 0 otherwise ! d_ij = 0 otherwise
! inspired by http://fortraninacworld.blogspot.de/2012/12/ternary-operator.html
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
real(pReal) pure function math_delta(i,j) real(pReal) pure function math_delta(i,j)
@ -317,7 +321,7 @@ end function math_cross
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief outer product A \otimes B of arbitrary sized vectors A and B !> @brief outer product of arbitrary sized vectors (A ⊗ B / i,j)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_outer(A,B) pure function math_outer(A,B)
@ -333,7 +337,7 @@ end function math_outer
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief outer product A \otimes B of arbitrary sized vectors A and B !> @brief inner product of arbitrary sized vectors (A · B / i,i)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
real(pReal) pure function math_inner(A,B) real(pReal) pure function math_inner(A,B)
@ -346,31 +350,26 @@ end function math_inner
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief matrix multiplication 33xx33 = 1 (double contraction --> ij * ij) !> @brief double contraction of 3x3 matrices (A : B / ij,ij)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
real(pReal) pure function math_mul33xx33(A,B) real(pReal) pure function math_tensordot(A,B)
real(pReal), dimension(3,3), intent(in) :: A,B real(pReal), dimension(3,3), intent(in) :: A,B
integer :: i,j
real(pReal), dimension(3,3) :: C
do i=1,3; do j=1,3 math_tensordot = sum(A*B)
C(i,j) = A(i,j) * B(i,j)
enddo; enddo
math_mul33xx33 = sum(C)
end function math_mul33xx33 end function math_tensordot
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief matrix multiplication 3333x33 = 33 (double contraction --> ijkl *kl = ij) !> @brief matrix double contraction 3333x33 = 33 (ijkl,kl)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_mul3333xx33(A,B) pure function math_mul3333xx33(A,B)
real(pReal), dimension(3,3,3,3), intent(in) :: A real(pReal), dimension(3,3,3,3), intent(in) :: A
real(pReal), dimension(3,3), intent(in) :: B real(pReal), dimension(3,3), intent(in) :: B
real(pReal), dimension(3,3) :: math_mul3333xx33 real(pReal), dimension(3,3) :: math_mul3333xx33
integer :: i,j integer :: i,j
do i=1,3; do j=1,3 do i=1,3; do j=1,3
math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
@ -380,7 +379,7 @@ end function math_mul3333xx33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief matrix multiplication 3333x3333 = 3333 (ijkl *klmn = ijmn) !> @brief matrix multiplication 3333x3333 = 3333 (ijkl,klmn)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_mul3333xx3333(A,B) pure function math_mul3333xx3333(A,B)
@ -402,8 +401,9 @@ end function math_mul3333xx3333
pure function math_exp33(A,n) pure function math_exp33(A,n)
real(pReal), dimension(3,3), intent(in) :: A real(pReal), dimension(3,3), intent(in) :: A
integer, intent(in), optional :: n integer, intent(in), optional :: n
real(pReal), dimension(3,3) :: B, math_exp33 real(pReal), dimension(3,3) :: B, math_exp33
real(pReal) :: invFac real(pReal) :: invFac
integer :: n_,i integer :: n_,i
@ -428,15 +428,16 @@ end function math_exp33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Cramer inversion of 33 matrix (function) !> @brief Cramer inversion of 33 matrix (function)
!> @details Direct Cramer inversion of matrix A. Returns all zeroes if not possible, i.e. !> @details Direct Cramer inversion of matrix A. Returns all zeroes if not possible, i.e.
! if determinant is close to zero ! if determinant is close to zero
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_inv33(A) pure function math_inv33(A)
real(pReal),dimension(3,3),intent(in) :: A real(pReal), dimension(3,3), intent(in) :: A
real(pReal) :: DetA real(pReal), dimension(3,3) :: math_inv33
real(pReal),dimension(3,3) :: math_inv33
logical :: error real(pReal) :: DetA
logical :: error
call math_invert33(math_inv33,DetA,error,A) call math_invert33(math_inv33,DetA,error,A)
if(error) math_inv33 = 0.0_pReal if(error) math_inv33 = 0.0_pReal
@ -451,10 +452,10 @@ end function math_inv33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine math_invert33(InvA, DetA, error, A) pure subroutine math_invert33(InvA, DetA, error, A)
logical, intent(out) :: error real(pReal), dimension(3,3), intent(out) :: InvA
real(pReal),dimension(3,3),intent(in) :: A real(pReal), intent(out) :: DetA
real(pReal),dimension(3,3),intent(out) :: InvA logical, intent(out) :: error
real(pReal), intent(out) :: DetA real(pReal), dimension(3,3), intent(in) :: A
InvA(1,1) = A(2,2) * A(3,3) - A(2,3) * A(3,2) InvA(1,1) = A(2,2) * A(3,3) - A(2,3) * A(3,2)
InvA(2,1) = -A(2,1) * A(3,3) + A(2,3) * A(3,1) InvA(2,1) = -A(2,1) * A(3,3) + A(2,3) * A(3,1)
@ -529,12 +530,12 @@ subroutine math_invert(InvA, error, A)
dgetrf, & dgetrf, &
dgetri dgetri
invA = A invA = A
call dgetrf(size(A,1),size(A,1),invA,size(A,1),ipiv,ierr) call dgetrf(size(A,1),size(A,1),invA,size(A,1),ipiv,ierr)
error = (ierr /= 0) error = (ierr /= 0)
call dgetri(size(A,1),InvA,size(A,1),ipiv,work,size(work,1),ierr) call dgetri(size(A,1),InvA,size(A,1),ipiv,work,size(work,1),ierr)
error = error .or. (ierr /= 0) error = error .or. (ierr /= 0)
end subroutine math_invert end subroutine math_invert
@ -621,7 +622,7 @@ end function math_trace33
real(pReal) pure function math_det33(m) real(pReal) pure function math_det33(m)
real(pReal), dimension(3,3), intent(in) :: m real(pReal), dimension(3,3), intent(in) :: m
math_det33 = m(1,1)* (m(2,2)*m(3,3)-m(2,3)*m(3,2)) & math_det33 = m(1,1)* (m(2,2)*m(3,3)-m(2,3)*m(3,2)) &
- m(1,2)* (m(2,1)*m(3,3)-m(2,3)*m(3,1)) & - m(1,2)* (m(2,1)*m(3,3)-m(2,3)*m(3,1)) &
+ m(1,3)* (m(2,1)*m(3,2)-m(2,2)*m(3,1)) + m(1,3)* (m(2,1)*m(3,2)-m(2,2)*m(3,1))
@ -635,7 +636,7 @@ end function math_det33
real(pReal) pure function math_detSym33(m) real(pReal) pure function math_detSym33(m)
real(pReal), dimension(3,3), intent(in) :: m real(pReal), dimension(3,3), intent(in) :: m
math_detSym33 = -(m(1,1)*m(2,3)**2 + m(2,2)*m(1,3)**2 + m(3,3)*m(1,2)**2) & math_detSym33 = -(m(1,1)*m(2,3)**2 + m(2,2)*m(1,3)**2 + m(3,3)*m(1,2)**2) &
+ m(1,1)*m(2,2)*m(3,3) + 2.0_pReal * m(1,2)*m(1,3)*m(2,3) + m(1,1)*m(2,2)*m(3,3) + 2.0_pReal * m(1,2)*m(1,3)*m(2,3)
@ -649,9 +650,9 @@ pure function math_33to9(m33)
real(pReal), dimension(9) :: math_33to9 real(pReal), dimension(9) :: math_33to9
real(pReal), dimension(3,3), intent(in) :: m33 real(pReal), dimension(3,3), intent(in) :: m33
integer :: i integer :: i
do i = 1, 9 do i = 1, 9
math_33to9(i) = m33(MAPPLAIN(1,i),MAPPLAIN(2,i)) math_33to9(i) = m33(MAPPLAIN(1,i),MAPPLAIN(2,i))
enddo enddo
@ -666,7 +667,7 @@ pure function math_9to33(v9)
real(pReal), dimension(3,3) :: math_9to33 real(pReal), dimension(3,3) :: math_9to33
real(pReal), dimension(9), intent(in) :: v9 real(pReal), dimension(9), intent(in) :: v9
integer :: i integer :: i
do i = 1, 9 do i = 1, 9
@ -687,10 +688,10 @@ pure function math_sym33to6(m33,weighted)
real(pReal), dimension(6) :: math_sym33to6 real(pReal), dimension(6) :: math_sym33to6
real(pReal), dimension(3,3), intent(in) :: m33 !< symmetric matrix (no internal check) real(pReal), dimension(3,3), intent(in) :: m33 !< symmetric matrix (no internal check)
logical, optional, intent(in) :: weighted !< weight according to Mandel (.true. by default) logical, optional, intent(in) :: weighted !< weight according to Mandel (.true. by default)
real(pReal), dimension(6) :: w real(pReal), dimension(6) :: w
integer :: i integer :: i
if(present(weighted)) then if(present(weighted)) then
w = merge(NRMMANDEL,1.0_pReal,weighted) w = merge(NRMMANDEL,1.0_pReal,weighted)
else else
@ -699,7 +700,7 @@ pure function math_sym33to6(m33,weighted)
do i = 1, 6 do i = 1, 6
math_sym33to6(i) = w(i)*m33(MAPNYE(1,i),MAPNYE(2,i)) math_sym33to6(i) = w(i)*m33(MAPNYE(1,i),MAPNYE(2,i))
enddo enddo
end function math_sym33to6 end function math_sym33to6
@ -718,7 +719,7 @@ pure function math_6toSym33(v6,weighted)
real(pReal), dimension(6) :: w real(pReal), dimension(6) :: w
integer :: i integer :: i
if(present(weighted)) then if(present(weighted)) then
w = merge(INVNRMMANDEL,1.0_pReal,weighted) w = merge(INVNRMMANDEL,1.0_pReal,weighted)
else else
@ -781,7 +782,7 @@ pure function math_sym3333to66(m3333,weighted)
real(pReal), dimension(6) :: w real(pReal), dimension(6) :: w
integer :: i,j integer :: i,j
if(present(weighted)) then if(present(weighted)) then
w = merge(NRMMANDEL,1.0_pReal,weighted) w = merge(NRMMANDEL,1.0_pReal,weighted)
else else
@ -806,10 +807,10 @@ pure function math_66toSym3333(m66,weighted)
real(pReal), dimension(3,3,3,3) :: math_66toSym3333 real(pReal), dimension(3,3,3,3) :: math_66toSym3333
real(pReal), dimension(6,6), intent(in) :: m66 real(pReal), dimension(6,6), intent(in) :: m66
logical, optional, intent(in) :: weighted !< weight according to Mandel (.true. by default) logical, optional, intent(in) :: weighted !< weight according to Mandel (.true. by default)
real(pReal), dimension(6) :: w real(pReal), dimension(6) :: w
integer :: i,j integer :: i,j
if(present(weighted)) then if(present(weighted)) then
w = merge(INVNRMMANDEL,1.0_pReal,weighted) w = merge(INVNRMMANDEL,1.0_pReal,weighted)
else else
@ -909,48 +910,48 @@ end subroutine math_eigenValuesVectorsSym
! ToDo: has wrong oder of arguments ! ToDo: has wrong oder of arguments
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine math_eigenValuesVectorsSym33(m,values,vectors) subroutine math_eigenValuesVectorsSym33(m,values,vectors)
real(pReal), dimension(3,3),intent(in) :: m real(pReal), dimension(3,3),intent(in) :: m
real(pReal), dimension(3), intent(out) :: values real(pReal), dimension(3), intent(out) :: values
real(pReal), dimension(3,3),intent(out) :: vectors real(pReal), dimension(3,3),intent(out) :: vectors
real(pReal) :: T, U, norm, threshold real(pReal) :: T, U, norm, threshold
logical :: error logical :: error
values = math_eigenvaluesSym33(m) values = math_eigenvaluesSym33(m)
vectors(1:3,2) = [ m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2), & vectors(1:3,2) = [ m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2), &
m(1, 3) * m(1, 2) - m(2, 3) * m(1, 1), & m(1, 3) * m(1, 2) - m(2, 3) * m(1, 1), &
m(1, 2)**2] m(1, 2)**2]
T = maxval(abs(values)) T = maxval(abs(values))
U = max(T, T**2) U = max(T, T**2)
threshold = sqrt(5.68e-14_pReal * U**2) threshold = sqrt(5.68e-14_pReal * U**2)
! Calculate first eigenvector by the formula v[0] = (m - lambda[0]).e1 x (m - lambda[0]).e2 ! Calculate first eigenvector by the formula v[0] = (m - lambda[0]).e1 x (m - lambda[0]).e2
vectors(1:3,1) = [ vectors(1,2) + m(1, 3) * values(1), & vectors(1:3,1) = [ vectors(1,2) + m(1, 3) * values(1), &
vectors(2,2) + m(2, 3) * values(1), & vectors(2,2) + m(2, 3) * values(1), &
(m(1,1) - values(1)) * (m(2,2) - values(1)) - vectors(3,2)] (m(1,1) - values(1)) * (m(2,2) - values(1)) - vectors(3,2)]
norm = norm2(vectors(1:3, 1)) norm = norm2(vectors(1:3, 1))
fallback1: if(norm < threshold) then fallback1: if(norm < threshold) then
call math_eigenValuesVectorsSym(m,values,vectors,error) call math_eigenValuesVectorsSym(m,values,vectors,error)
return return
endif fallback1 endif fallback1
vectors(1:3,1) = vectors(1:3, 1) / norm vectors(1:3,1) = vectors(1:3, 1) / norm
! Calculate second eigenvector by the formula v[1] = (m - lambda[1]).e1 x (m - lambda[1]).e2 ! Calculate second eigenvector by the formula v[1] = (m - lambda[1]).e1 x (m - lambda[1]).e2
vectors(1:3,2) = [ vectors(1,2) + m(1, 3) * values(2), & vectors(1:3,2) = [ vectors(1,2) + m(1, 3) * values(2), &
vectors(2,2) + m(2, 3) * values(2), & vectors(2,2) + m(2, 3) * values(2), &
(m(1,1) - values(2)) * (m(2,2) - values(2)) - vectors(3,2)] (m(1,1) - values(2)) * (m(2,2) - values(2)) - vectors(3,2)]
norm = norm2(vectors(1:3, 2)) norm = norm2(vectors(1:3, 2))
fallback2: if(norm < threshold) then fallback2: if(norm < threshold) then
call math_eigenValuesVectorsSym(m,values,vectors,error) call math_eigenValuesVectorsSym(m,values,vectors,error)
return return
endif fallback2 endif fallback2
vectors(1:3,2) = vectors(1:3, 2) / norm vectors(1:3,2) = vectors(1:3, 2) / norm
! Calculate third eigenvector according to v[2] = v[0] x v[1] ! Calculate third eigenvector according to v[2] = v[0] x v[1]
vectors(1:3,3) = math_cross(vectors(1:3,1),vectors(1:3,2)) vectors(1:3,3) = math_cross(vectors(1:3,1),vectors(1:3,2))
@ -968,11 +969,11 @@ function math_eigenvectorBasisSym(m)
real(pReal), dimension(size(m,1),size(m,1)) :: math_eigenvectorBasisSym real(pReal), dimension(size(m,1),size(m,1)) :: math_eigenvectorBasisSym
logical :: error logical :: error
integer :: i integer :: i
math_eigenvectorBasisSym = 0.0_pReal math_eigenvectorBasisSym = 0.0_pReal
call math_eigenValuesVectorsSym(m,values,vectors,error) call math_eigenValuesVectorsSym(m,values,vectors,error)
if(error) return if(error) return
do i=1, size(m,1) do i=1, size(m,1)
math_eigenvectorBasisSym = math_eigenvectorBasisSym & math_eigenvectorBasisSym = math_eigenvectorBasisSym &
+ sqrt(values(i)) * math_outer(vectors(:,i),vectors(:,i)) + sqrt(values(i)) * math_outer(vectors(:,i),vectors(:,i))
@ -992,13 +993,13 @@ pure function math_eigenvectorBasisSym33(m)
real(pReal) :: P, Q, rho, phi real(pReal) :: P, Q, rho, phi
real(pReal), parameter :: TOL=1.e-14_pReal real(pReal), parameter :: TOL=1.e-14_pReal
real(pReal), dimension(3,3,3) :: N, EB real(pReal), dimension(3,3,3) :: N, EB
invariants = math_invariantsSym33(m) invariants = math_invariantsSym33(m)
EB = 0.0_pReal EB = 0.0_pReal
P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal
Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3) Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3)
threeSimilarEigenvalues: if(all(abs([P,Q]) < TOL)) then threeSimilarEigenvalues: if(all(abs([P,Q]) < TOL)) then
values = invariants(1)/3.0_pReal values = invariants(1)/3.0_pReal
! this is not really correct, but at least the basis is correct ! this is not really correct, but at least the basis is correct
@ -1016,7 +1017,7 @@ pure function math_eigenvectorBasisSym33(m)
N(1:3,1:3,1) = m-values(1)*math_I3 N(1:3,1:3,1) = m-values(1)*math_I3
N(1:3,1:3,2) = m-values(2)*math_I3 N(1:3,1:3,2) = m-values(2)*math_I3
N(1:3,1:3,3) = m-values(3)*math_I3 N(1:3,1:3,3) = m-values(3)*math_I3
twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then
EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ & EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ &
((values(3)-values(1))*(values(3)-values(2))) ((values(3)-values(1))*(values(3)-values(2)))
EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3) EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3)
@ -1024,7 +1025,7 @@ pure function math_eigenvectorBasisSym33(m)
EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ & EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ &
((values(1)-values(2))*(values(1)-values(3))) ((values(1)-values(2))*(values(1)-values(3)))
EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1) EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1)
elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues
EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ & EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ &
((values(2)-values(1))*(values(2)-values(3))) ((values(2)-values(1))*(values(2)-values(3)))
EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2) EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2)
@ -1080,7 +1081,7 @@ pure function math_eigenvectorBasisSym33_log(m)
N(1:3,1:3,1) = m-values(1)*math_I3 N(1:3,1:3,1) = m-values(1)*math_I3
N(1:3,1:3,2) = m-values(2)*math_I3 N(1:3,1:3,2) = m-values(2)*math_I3
N(1:3,1:3,3) = m-values(3)*math_I3 N(1:3,1:3,3) = m-values(3)*math_I3
twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then
EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ & EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ &
((values(3)-values(1))*(values(3)-values(2))) ((values(3)-values(1))*(values(3)-values(2)))
EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3) EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3)
@ -1088,7 +1089,7 @@ pure function math_eigenvectorBasisSym33_log(m)
EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ & EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ &
((values(1)-values(2))*(values(1)-values(3))) ((values(1)-values(2))*(values(1)-values(3)))
EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1) EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1)
elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues
EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ & EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ &
((values(2)-values(1))*(values(2)-values(3))) ((values(2)-values(1))*(values(2)-values(3)))
EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2) EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2)
@ -1136,7 +1137,7 @@ end function math_rotationalPart33
! will return NaN on error ! will return NaN on error
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function math_eigenvaluesSym(m) function math_eigenvaluesSym(m)
real(pReal), dimension(:,:), intent(in) :: m real(pReal), dimension(:,:), intent(in) :: m
real(pReal), dimension(size(m,1)) :: math_eigenvaluesSym real(pReal), dimension(size(m,1)) :: math_eigenvaluesSym
real(pReal), dimension(size(m,1),size(m,1)) :: m_ real(pReal), dimension(size(m,1),size(m,1)) :: m_
@ -1221,7 +1222,7 @@ integer pure function math_binomial(n,k)
integer, intent(in) :: n, k integer, intent(in) :: n, k
integer :: i, k_, n_ integer :: i, k_, n_
k_ = min(k,n-k) k_ = min(k,n-k)
n_ = n n_ = n
math_binomial = merge(1,0,k_>-1) ! handling special cases k < 0 or k > n math_binomial = merge(1,0,k_>-1) ! handling special cases k < 0 or k > n
@ -1244,7 +1245,7 @@ integer pure function math_multinomial(alpha)
math_multinomial = 1 math_multinomial = 1
do i = 1, size(alpha) do i = 1, size(alpha)
math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i)) math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i))
enddo enddo
end function math_multinomial end function math_multinomial
@ -1307,6 +1308,7 @@ subroutine unitTest
sort_out_ = reshape([-1,-1, +1,+5, +5,+6, +3,-2],[2,4]) sort_out_ = reshape([-1,-1, +1,+5, +5,+6, +3,-2],[2,4])
integer, dimension(5) :: range_out_ = [1,2,3,4,5] integer, dimension(5) :: range_out_ = [1,2,3,4,5]
integer, dimension(3) :: ijk
real(pReal) :: det real(pReal) :: det
real(pReal), dimension(3) :: v3_1,v3_2,v3_3,v3_4 real(pReal), dimension(3) :: v3_1,v3_2,v3_3,v3_4
@ -1349,7 +1351,7 @@ subroutine unitTest
call random_number(v9) call random_number(v9)
if(any(dNeq(math_33to9(math_9to33(v9)),v9))) & if(any(dNeq(math_33to9(math_9to33(v9)),v9))) &
call IO_error(0,ext_msg='math_33to9/math_9to33') call IO_error(0,ext_msg='math_33to9/math_9to33')
call random_number(t99) call random_number(t99)
if(any(dNeq(math_3333to99(math_99to3333(t99)),t99))) & if(any(dNeq(math_3333to99(math_99to3333(t99)),t99))) &
call IO_error(0,ext_msg='math_3333to99/math_99to3333') call IO_error(0,ext_msg='math_3333to99/math_99to3333')
@ -1378,7 +1380,7 @@ subroutine unitTest
call random_number(t33) call random_number(t33)
if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) & if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) &
call IO_error(0,ext_msg='math_det33/math_detSym33') call IO_error(0,ext_msg='math_det33/math_detSym33')
if(any(dNeq0(math_identity2nd(3),math_inv33(math_I3)))) & if(any(dNeq0(math_identity2nd(3),math_inv33(math_I3)))) &
call IO_error(0,ext_msg='math_inv33(math_I3)') call IO_error(0,ext_msg='math_inv33(math_I3)')
@ -1410,18 +1412,28 @@ subroutine unitTest
if(any(dNeq0(txx_2,txx) .or. e)) & if(any(dNeq0(txx_2,txx) .or. e)) &
call IO_error(0,ext_msg='math_invert(txx)/math_identity2nd') call IO_error(0,ext_msg='math_invert(txx)/math_identity2nd')
call math_invert(t99_2,e,t99) ! not sure how likely it is that we get a singular matrix call math_invert(t99_2,e,t99) ! not sure how likely it is that we get a singular matrix
if(any(dNeq0(matmul(t99_2,t99)-math_identity2nd(9),tol=1.0e-9_pReal)) .or. e) & if(any(dNeq0(matmul(t99_2,t99)-math_identity2nd(9),tol=1.0e-9_pReal)) .or. e) &
call IO_error(0,ext_msg='math_invert(t99)') call IO_error(0,ext_msg='math_invert(t99)')
if(any(dNeq(math_clip([4.0_pReal,9.0_pReal],5.0_pReal,6.5_pReal),[5.0_pReal,6.5_pReal]))) & if(any(dNeq(math_clip([4.0_pReal,9.0_pReal],5.0_pReal,6.5_pReal),[5.0_pReal,6.5_pReal]))) &
call IO_error(0,ext_msg='math_clip') call IO_error(0,ext_msg='math_clip')
if(math_factorial(10) /= 3628800) & if(math_factorial(10) /= 3628800) &
call IO_error(0,ext_msg='math_factorial') call IO_error(0,ext_msg='math_factorial')
if(math_binomial(49,6) /= 13983816) & if(math_binomial(49,6) /= 13983816) &
call IO_error(0,ext_msg='math_binomial') call IO_error(0,ext_msg='math_binomial')
ijk = cshift([1,2,3],int(r*1.0e2_pReal))
if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) &
call IO_error(0,ext_msg='math_LeviCivita(even)')
ijk = cshift([3,2,1],int(r*2.0e2_pReal))
if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),-1.0_pReal)) &
call IO_error(0,ext_msg='math_LeviCivita(odd)')
ijk = cshift([2,2,1],int(r*2.0e2_pReal))
if(dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3))))&
call IO_error(0,ext_msg='math_LeviCivita')
end subroutine unitTest end subroutine unitTest