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

further commenting, removing not used variables and code (marc2vtk)

This commit is contained in:
Martin Diehl 2016-03-04 18:50:13 +01:00
parent d1e8f69857
commit a4bbdd5ecb
8 changed files with 194 additions and 273 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
lib/damask/__init__.py
View File

@ -1,27 +1,27 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
# $Id$ """Main aggregator"""
import sys, os import sys, os
with open(os.path.join(os.path.dirname(__file__),'../../VERSION')) as f: with open(os.path.join(os.path.dirname(__file__),'../../VERSION')) as f:
version = f.readline()[:-1] version = f.readline()[:-1]
from .environment import Environment # only one class from .environment import Environment # noqa
from .asciitable import ASCIItable # only one class from .asciitable import ASCIItable # noqa
from .config import Material # will be extended to debug and numerics from .config import Material # noqa
from .colormaps import Colormap, Color from .colormaps import Colormap, Color # noqa
from .orientation import Quaternion, Rodrigues, Symmetry, Orientation from .orientation import Quaternion, Rodrigues, Symmetry, Orientation # noqa
# try: # try:
# from .corientation import Quaternion, Rodrigues, Symmetry, Orientation # from .corientation import Quaternion, Rodrigues, Symmetry, Orientation
# print "Import Cython version of Orientation module" # print "Import Cython version of Orientation module"
# except: # except:
# from .orientation import Quaternion, Rodrigues, Symmetry, Orientation # from .orientation import Quaternion, Rodrigues, Symmetry, Orientation
#from .block import Block # only one class #from .block import Block # only one class
from .result import Result # only one class from .result import Result # noqa
from .geometry import Geometry # one class with subclasses from .geometry import Geometry # noqa
from .solver import Solver # one class with subclasses from .solver import Solver # noqa
from .test import Test from .test import Test # noqa
from .util import extendableOption from .util import extendableOption # noqa
try: try:
from . import core from . import core

128
lib/damask/asciitable.py
View File

@ -4,12 +4,9 @@
import os,sys import os,sys
import numpy as np import numpy as np
import util
class ASCIItable(): class ASCIItable():
''' """Read and write to ASCII tables"""
There should be a doc string here :)
'''
__slots__ = ['__IO__', __slots__ = ['__IO__',
'info', 'info',
@ -58,8 +55,8 @@ class ASCIItable():
self.data = [] self.data = []
self.line = '' self.line = ''
if self.__IO__['in'] == None \ if self.__IO__['in'] is None \
or self.__IO__['out'] == 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 _transliterateToFloat(self, def _transliterateToFloat(self,
@ -86,9 +83,7 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def output_write(self, def output_write(self,
what): what):
''' """aggregate a single row (string) or list of (possibly containing further lists of) rows into output"""
aggregate a single row (string) or list of (possibly containing further lists of) rows into output
'''
if not isinstance(what, (str, unicode)): if not isinstance(what, (str, unicode)):
try: try:
for item in what: self.output_write(item) for item in what: self.output_write(item)
@ -104,7 +99,7 @@ class ASCIItable():
clear = True): clear = True):
try: try:
self.__IO__['output'] == [] or self.__IO__['out'].write('\n'.join(self.__IO__['output']) + '\n') self.__IO__['output'] == [] or self.__IO__['out'].write('\n'.join(self.__IO__['output']) + '\n')
except IOError as e: except IOError:
return False return False
if clear: self.output_clear() if clear: self.output_clear()
return True return True
@ -127,11 +122,12 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def head_read(self): def head_read(self):
''' """
get column labels by either reading get column labels by either reading
the first row or, if keyword "head[*]" is present,
the last line of the header the first row or, if keyword "head[*]" is present,
''' the last line of the header
"""
import re import re
try: try:
@ -180,10 +176,7 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def head_write(self, def head_write(self,
header = True): header = True):
''' """write current header information (info + labels)"""
write current header information (info + labels)
'''
head = ['{}\theader'.format(len(self.info)+self.__IO__['labeled'])] if header else [] head = ['{}\theader'.format(len(self.info)+self.__IO__['labeled'])] if header else []
head.append(self.info) head.append(self.info)
if self.__IO__['labeled']: head.append('\t'.join(self.labels)) if self.__IO__['labeled']: head.append('\t'.join(self.labels))
@ -192,9 +185,7 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def head_getGeom(self): def head_getGeom(self):
''' """interpret geom header"""
interpret geom header
'''
identifiers = { identifiers = {
'grid': ['a','b','c'], 'grid': ['a','b','c'],
'size': ['x','y','z'], 'size': ['x','y','z'],
@ -234,9 +225,7 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def head_putGeom(self,info): def head_putGeom(self,info):
''' """translate geometry description to header"""
translate geometry description to header
'''
self.info_append([ self.info_append([
"grid\ta {}\tb {}\tc {}".format(*info['grid']), "grid\ta {}\tb {}\tc {}".format(*info['grid']),
"size\tx {}\ty {}\tz {}".format(*info['size']), "size\tx {}\ty {}\tz {}".format(*info['size']),
@ -249,9 +238,7 @@ class ASCIItable():
def labels_append(self, def labels_append(self,
what, what,
reset = False): reset = False):
''' """add item or list to existing set of labels (and switch on labeling)"""
add item or list to existing set of labels (and switch on labeling)
'''
if not isinstance(what, (str, unicode)): if not isinstance(what, (str, unicode)):
try: try:
for item in what: self.labels_append(item) for item in what: self.labels_append(item)
@ -265,26 +252,25 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def labels_clear(self): def labels_clear(self):
''' """delete existing labels and switch to no labeling"""
delete existing labels and switch to no labeling
'''
self.labels = [] self.labels = []
self.__IO__['labeled'] = False self.__IO__['labeled'] = False
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def label_index(self, def label_index(self,
labels): labels):
''' """
tell index of column label(s). tell index of column label(s).
return numpy array if asked for list of labels.
transparently deals with label positions implicitly given as numbers or their headings given as strings. return numpy array if asked for list of labels.
''' transparently deals with label positions implicitly given as numbers or their headings given as strings.
"""
from collections import Iterable from collections import Iterable
if isinstance(labels, Iterable) and not isinstance(labels, str): # check whether list of labels is requested if isinstance(labels, Iterable) and not isinstance(labels, str): # check whether list of labels is requested
idx = [] idx = []
for label in labels: for label in labels:
if label != None: if label is not None:
try: try:
idx.append(int(label)) # column given as integer number? idx.append(int(label)) # column given as integer number?
except ValueError: except ValueError:
@ -305,25 +291,25 @@ class ASCIItable():
try: try:
idx = self.labels.index('1_'+labels) # locate '1_'+string in label list idx = self.labels.index('1_'+labels) # locate '1_'+string in label list
except ValueError: except ValueError:
idx = None if labels == None else -1 idx = None if labels is None else -1
return np.array(idx) if isinstance(idx,list) else idx return np.array(idx) if isinstance(idx,list) else idx
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def label_dimension(self, def label_dimension(self,
labels): labels):
''' """
tell dimension (length) of column label(s). tell dimension (length) of column label(s).
return numpy array if asked for list of labels.
transparently deals with label positions implicitly given as numbers or their headings given as strings.
'''
return numpy array if asked for list of labels.
transparently deals with label positions implicitly given as numbers or their headings given as strings.
"""
from collections import Iterable from collections import Iterable
if isinstance(labels, Iterable) and not isinstance(labels, str): # check whether list of labels is requested if isinstance(labels, Iterable) and not isinstance(labels, str): # check whether list of labels is requested
dim = [] dim = []
for label in labels: for label in labels:
if label != None: if label is not None:
myDim = -1 myDim = -1
try: # column given as number? try: # column given as number?
idx = int(label) idx = int(label)
@ -364,12 +350,12 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def label_indexrange(self, def label_indexrange(self,
labels): labels):
''' """
tell index range for given label(s). tell index range for given label(s).
return numpy array if asked for list of labels.
transparently deals with label positions implicitly given as numbers or their headings given as strings.
'''
return numpy array if asked for list of labels.
transparently deals with label positions implicitly given as numbers or their headings given as strings.
"""
from collections import Iterable from collections import Iterable
start = self.label_index(labels) start = self.label_index(labels)
@ -381,9 +367,7 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def info_append(self, def info_append(self,
what): what):
''' """add item or list to existing set of infos"""
add item or list to existing set of infos
'''
if not isinstance(what, (str, unicode)): if not isinstance(what, (str, unicode)):
try: try:
for item in what: self.info_append(item) for item in what: self.info_append(item)
@ -394,9 +378,7 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def info_clear(self): def info_clear(self):
''' """delete any info block"""
delete any info block
'''
self.info = [] self.info = []
# ------------------------------------------------------------------ # ------------------------------------------------------------------
@ -409,9 +391,7 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def data_skipLines(self, def data_skipLines(self,
count): count):
''' """wind forward by count number of lines"""
wind forward by count number of lines
'''
for i in xrange(count): for i in xrange(count):
alive = self.data_read() alive = self.data_read()
@ -421,9 +401,7 @@ class ASCIItable():
def data_read(self, def data_read(self,
advance = True, advance = True,
respectLabels = True): respectLabels = True):
''' """read next line (possibly buffered) and parse it into data array"""
read next line (possibly buffered) and parse it into data array
'''
self.line = self.__IO__['readBuffer'].pop(0) if len(self.__IO__['readBuffer']) > 0 \ self.line = self.__IO__['readBuffer'].pop(0) if len(self.__IO__['readBuffer']) > 0 \
else self.__IO__['in'].readline().strip() # take buffered content or get next data row from file else self.__IO__['in'].readline().strip() # take buffered content or get next data row from file
@ -434,7 +412,7 @@ class ASCIItable():
if self.__IO__['labeled'] and respectLabels: # if table has labels if self.__IO__['labeled'] and respectLabels: # if table has labels
items = self.line.split()[:len(self.__IO__['labels'])] # use up to label count (from original file info) items = self.line.split()[:len(self.__IO__['labels'])] # use up to label count (from original file info)
self.data = items if len(items) == len(self.__IO__['labels']) else [] # take entries if correct number, i.e. not too few compared to label count self.data = items if len(items) == len(self.__IO__['labels']) else [] # take entries if label count matches
else: else:
self.data = self.line.split() # otherwise take all self.data = self.line.split() # otherwise take all
@ -443,9 +421,7 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def data_readArray(self, def data_readArray(self,
labels = []): labels = []):
''' """read whole data of all (given) labels as numpy array"""
read whole data of all (given) labels as numpy array
'''
from collections import Iterable from collections import Iterable
try: try:
@ -453,7 +429,7 @@ class ASCIItable():
except: except:
pass # assume/hope we are at data start already... pass # assume/hope we are at data start already...
if labels == None or labels == []: if labels is None or labels == []:
use = None # use all columns (and keep labels intact) use = None # use all columns (and keep labels intact)
labels_missing = [] labels_missing = []
else: else:
@ -467,9 +443,10 @@ 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
columns += range(c,c + \ columns += range(c,c + \
(d if str(c) != str(labels[present[i]]) else \ (d if str(c) != str(labels[present[i]]) else \
1)) # ... transparently add all components unless column referenced by number or with explicit dimension 1))
use = np.array(columns) use = np.array(columns)
self.labels = list(np.array(self.labels)[use]) # update labels with valid subset self.labels = list(np.array(self.labels)[use]) # update labels with valid subset
@ -481,9 +458,7 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def data_write(self, def data_write(self,
delimiter = '\t'): delimiter = '\t'):
''' """write current data array and report alive output back"""
write current data array and report alive output back
'''
if len(self.data) == 0: return True if len(self.data) == 0: return True
if isinstance(self.data[0],list): if isinstance(self.data[0],list):
@ -495,9 +470,7 @@ class ASCIItable():
def data_writeArray(self, def data_writeArray(self,
fmt = None, fmt = None,
delimiter = '\t'): delimiter = '\t'):
''' """write whole numpy array data"""
write whole numpy array data
'''
for row in self.data: for row in self.data:
try: try:
output = [fmt % value for value in row] if fmt else map(repr,row) output = [fmt % value for value in row] if fmt else map(repr,row)
@ -520,9 +493,7 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def data_set(self, def data_set(self,
what, where): what, where):
''' """update data entry in column "where". grows data array if needed."""
update data entry in column "where". grows data array if needed.
'''
idx = -1 idx = -1
try: try:
idx = self.label_index(where) idx = self.label_index(where)
@ -547,10 +518,7 @@ class ASCIItable():
# ------------------------------------------------------------------ # ------------------------------------------------------------------
def microstructure_read(self, def microstructure_read(self,
grid): grid):
''' """read microstructure data (from .geom format)"""
read microstructure data (from .geom format)
'''
N = grid.prod() # expected number of microstructure indices in data N = grid.prod() # expected number of microstructure indices in data
microstructure = np.zeros(N,'i') # initialize as flat array microstructure = np.zeros(N,'i') # initialize as flat array

2
lib/damask/config/__init__.py
View File

@ -2,4 +2,4 @@
"""Aggregator for configuration file handling""" """Aggregator for configuration file handling"""
from .material import Material #noqa from .material import Material # noqa

3
lib/damask/config/material.py
View File

@ -243,7 +243,8 @@ class Material():
except AttributeError: except AttributeError:
pass pass
for (phase,texture,fraction,crystallite) in zip(components['phase'],components['texture'],components['fraction'],components['crystallite']): for (phase,texture,fraction,crystallite) in zip(components['phase'],components['texture'],
components['fraction'],components['crystallite']):
microstructure.add_multiKey('constituent','phase %i\ttexture %i\tfraction %g\ncrystallite %i'%( microstructure.add_multiKey('constituent','phase %i\ttexture %i\tfraction %g\ncrystallite %i'%(
self.data['phase']['__order__'].index(phase)+1, self.data['phase']['__order__'].index(phase)+1,
self.data['texture']['__order__'].index(texture)+1, self.data['texture']['__order__'].index(texture)+1,

144
lib/damask/orientation.py
View File

@ -9,7 +9,6 @@ import numpy as np
# ****************************************************************************************** # ******************************************************************************************
class Rodrigues: class Rodrigues:
# ******************************************************************************************
def __init__(self, vector = np.zeros(3)): def __init__(self, vector = np.zeros(3)):
self.vector = vector self.vector = vector
@ -28,20 +27,22 @@ class Rodrigues:
# ****************************************************************************************** # ******************************************************************************************
class Quaternion: class Quaternion:
# ****************************************************************************************** """
# All methods and naming conventions based off Orientation represented as unit quaternion
# http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions
# w is the real part, (x, y, z) are the imaginary parts All methods and naming conventions based on http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions
# Representation of rotation is in ACTIVE form! w is the real part, (x, y, z) are the imaginary parts
# (derived directly or through angleAxis, Euler angles, or active matrix) Representation of rotation is in ACTIVE form!
# vector "a" (defined in coordinate system "A") is actively rotated to new coordinates "b" (derived directly or through angleAxis, Euler angles, or active matrix)
# b = Q * a vector "a" (defined in coordinate system "A") is actively rotated to new coordinates "b"
# b = np.dot(Q.asMatrix(),a) b = Q * a
b = np.dot(Q.asMatrix(),a)
"""
def __init__(self, def __init__(self,
quatArray = [1.0,0.0,0.0,0.0]): quatArray = [1.0,0.0,0.0,0.0]):
"""initializes to identity if not given"""
self.w, \ self.w, \
self.x, \ self.x, \
self.y, \ self.y, \
@ -49,19 +50,23 @@ class Quaternion:
self.homomorph() self.homomorph()
def __iter__(self): def __iter__(self):
"""components"""
return iter([self.w,self.x,self.y,self.z]) return iter([self.w,self.x,self.y,self.z])
def __copy__(self): def __copy__(self):
"""create copy"""
Q = Quaternion([self.w,self.x,self.y,self.z]) Q = Quaternion([self.w,self.x,self.y,self.z])
return Q return Q
copy = __copy__ copy = __copy__
def __repr__(self): def __repr__(self):
"""readbable string"""
return 'Quaternion(real=%+.6f, imag=<%+.6f, %+.6f, %+.6f>)' % \ return 'Quaternion(real=%+.6f, imag=<%+.6f, %+.6f, %+.6f>)' % \
(self.w, self.x, self.y, self.z) (self.w, self.x, self.y, self.z)
def __pow__(self, exponent): def __pow__(self, exponent):
"""power"""
omega = math.acos(self.w) omega = math.acos(self.w)
vRescale = math.sin(exponent*omega)/math.sin(omega) vRescale = math.sin(exponent*omega)/math.sin(omega)
Q = Quaternion() Q = Quaternion()
@ -72,6 +77,7 @@ class Quaternion:
return Q return Q
def __ipow__(self, exponent): def __ipow__(self, exponent):
"""in place power"""
omega = math.acos(self.w) omega = math.acos(self.w)
vRescale = math.sin(exponent*omega)/math.sin(omega) vRescale = math.sin(exponent*omega)/math.sin(omega)
self.w = np.cos(exponent*omega) self.w = np.cos(exponent*omega)
@ -81,6 +87,7 @@ class Quaternion:
return self return self
def __mul__(self, other): def __mul__(self, other):
"""multiplication"""
try: # quaternion try: # quaternion
Aw = self.w Aw = self.w
Ax = self.x Ax = self.x
@ -128,6 +135,7 @@ class Quaternion:
return self.copy() return self.copy()
def __imul__(self, other): def __imul__(self, other):
"""in place multiplication"""
try: # Quaternion try: # Quaternion
Ax = self.x Ax = self.x
Ay = self.y Ay = self.y
@ -145,6 +153,7 @@ class Quaternion:
return self return self
def __div__(self, other): def __div__(self, other):
"""division"""
if isinstance(other, (int,float,long)): if isinstance(other, (int,float,long)):
w = self.w / other w = self.w / other
x = self.x / other x = self.x / other
@ -155,6 +164,7 @@ class Quaternion:
return NotImplemented return NotImplemented
def __idiv__(self, other): def __idiv__(self, other):
"""in place division"""
if isinstance(other, (int,float,long)): if isinstance(other, (int,float,long)):
self.w /= other self.w /= other
self.x /= other self.x /= other
@ -163,6 +173,7 @@ class Quaternion:
return self return self
def __add__(self, other): def __add__(self, other):
"""addition"""
if isinstance(other, Quaternion): if isinstance(other, Quaternion):
w = self.w + other.w w = self.w + other.w
x = self.x + other.x x = self.x + other.x
@ -173,6 +184,7 @@ class Quaternion:
return NotImplemented return NotImplemented
def __iadd__(self, other): def __iadd__(self, other):
"""in place division"""
if isinstance(other, Quaternion): if isinstance(other, Quaternion):
self.w += other.w self.w += other.w
self.x += other.x self.x += other.x
@ -181,6 +193,7 @@ class Quaternion:
return self return self
def __sub__(self, other): def __sub__(self, other):
"""subtraction"""
if isinstance(other, Quaternion): if isinstance(other, Quaternion):
Q = self.copy() Q = self.copy()
Q.w -= other.w Q.w -= other.w
@ -192,6 +205,7 @@ class Quaternion:
return self.copy() return self.copy()
def __isub__(self, other): def __isub__(self, other):
"""in place subtraction"""
if isinstance(other, Quaternion): if isinstance(other, Quaternion):
self.w -= other.w self.w -= other.w
self.x -= other.x self.x -= other.x
@ -200,6 +214,7 @@ class Quaternion:
return self return self
def __neg__(self): def __neg__(self):
"""additive inverse"""
self.w = -self.w self.w = -self.w
self.x = -self.x self.x = -self.x
self.y = -self.y self.y = -self.y
@ -207,6 +222,7 @@ class Quaternion:
return self return self
def __abs__(self): def __abs__(self):
"""norm"""
return math.sqrt(self.w ** 2 + \ return math.sqrt(self.w ** 2 + \
self.x ** 2 + \ self.x ** 2 + \
self.y ** 2 + \ self.y ** 2 + \
@ -215,6 +231,7 @@ class Quaternion:
magnitude = __abs__ magnitude = __abs__
def __eq__(self,other): def __eq__(self,other):
"""equal at e-8 precision"""
return (abs(self.w-other.w) < 1e-8 and \ return (abs(self.w-other.w) < 1e-8 and \
abs(self.x-other.x) < 1e-8 and \ abs(self.x-other.x) < 1e-8 and \
abs(self.y-other.y) < 1e-8 and \ abs(self.y-other.y) < 1e-8 and \
@ -226,9 +243,11 @@ class Quaternion:
abs(-self.z-other.z) < 1e-8) abs(-self.z-other.z) < 1e-8)
def __ne__(self,other): def __ne__(self,other):
"""not equal at e-8 precision"""
return not self.__eq__(self,other) return not self.__eq__(self,other)
def __cmp__(self,other): def __cmp__(self,other):
"""linear ordering"""
return cmp(self.Rodrigues(),other.Rodrigues()) return cmp(self.Rodrigues(),other.Rodrigues())
def magnitude_squared(self): def magnitude_squared(self):
@ -290,9 +309,10 @@ class Quaternion:
return np.outer([i for i in self],[i for i in self]) return np.outer([i for i in self],[i for i in self])
def asMatrix(self): def asMatrix(self):
return np.array([[1.0-2.0*(self.y*self.y+self.z*self.z), 2.0*(self.x*self.y-self.z*self.w), 2.0*(self.x*self.z+self.y*self.w)], return np.array(
[ 2.0*(self.x*self.y+self.z*self.w), 1.0-2.0*(self.x*self.x+self.z*self.z), 2.0*(self.y*self.z-self.x*self.w)], [[1.0-2.0*(self.y*self.y+self.z*self.z), 2.0*(self.x*self.y-self.z*self.w), 2.0*(self.x*self.z+self.y*self.w)],
[ 2.0*(self.x*self.z-self.y*self.w), 2.0*(self.x*self.w+self.y*self.z), 1.0-2.0*(self.x*self.x+self.y*self.y)]]) [ 2.0*(self.x*self.y+self.z*self.w), 1.0-2.0*(self.x*self.x+self.z*self.z), 2.0*(self.y*self.z-self.x*self.w)],
[ 2.0*(self.x*self.z-self.y*self.w), 2.0*(self.x*self.w+self.y*self.z), 1.0-2.0*(self.x*self.x+self.y*self.y)]])
def asAngleAxis(self, def asAngleAxis(self,
degrees = False): degrees = False):
@ -315,15 +335,17 @@ class Quaternion:
return np.inf*np.ones(3) if self.w == 0.0 else np.array([self.x, self.y, self.z])/self.w return np.inf*np.ones(3) if self.w == 0.0 else np.array([self.x, self.y, self.z])/self.w
def asEulers(self, def asEulers(self,
type = 'bunge', type = "bunge",
degrees = False, degrees = False,
standardRange = False): standardRange = False):
''' u"""
Orientation as Bunge-Euler angles
conversion of ACTIVE rotation to Euler angles taken from: conversion of ACTIVE rotation to Euler angles taken from:
Melcher, A.; Unser, A.; Reichhardt, M.; Nestler, B.; Pötschke, M.; Selzer, M. Melcher, A.; Unser, A.; Reichhardt, M.; Nestler, B.; Pötschke, M.; Selzer, M.
Conversion of EBSD data by a quaternion based algorithm to be used for grain structure simulations Conversion of EBSD data by a quaternion based algorithm to be used for grain structure simulations
Technische Mechanik 30 (2010) pp 401--413 Technische Mechanik 30 (2010) pp 401--413
''' """
angles = [0.0,0.0,0.0] angles = [0.0,0.0,0.0]
if type.lower() == 'bunge' or type.lower() == 'zxz': if type.lower() == 'bunge' or type.lower() == 'zxz':
@ -369,7 +391,7 @@ class Quaternion:
@classmethod @classmethod
def fromRandom(cls,randomSeed = None): def fromRandom(cls,randomSeed = None):
if randomSeed == None: if randomSeed is None:
randomSeed = int(os.urandom(4).encode('hex'), 16) randomSeed = int(os.urandom(4).encode('hex'), 16)
np.random.seed(randomSeed) np.random.seed(randomSeed)
r = np.random.random(3) r = np.random.random(3)
@ -420,7 +442,6 @@ class Quaternion:
y = - c1 * s2 * s3 + s1 * s2 * c3 y = - c1 * s2 * s3 + s1 * s2 * c3
z = c1 * c2 * s3 + s1 * c2 * c3 z = c1 * c2 * s3 + s1 * c2 * c3
else: else:
# print 'unknown Euler convention'
w = c1 * c2 * c3 - s1 * s2 * s3 w = c1 * c2 * c3 - s1 * s2 * s3
x = s1 * s2 * c3 + c1 * c2 * s3 x = s1 * s2 * c3 + c1 * c2 * s3
y = s1 * c2 * c3 + c1 * s2 * s3 y = s1 * c2 * c3 + c1 * s2 * s3
@ -428,7 +449,8 @@ class Quaternion:
return cls([w,x,y,z]) return cls([w,x,y,z])
## Modified Method to calculate Quaternion from Orientation Matrix, Source: http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/ # Modified Method to calculate Quaternion from Orientation Matrix,
# Source: http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/
@classmethod @classmethod
def fromMatrix(cls, m): def fromMatrix(cls, m):
@ -482,8 +504,12 @@ class Quaternion:
@classmethod @classmethod
def new_interpolate(cls, q1, q2, t): def new_interpolate(cls, q1, q2, t):
# see http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20070017872_2007014421.pdf for (another?) way to interpolate quaternions """
interpolation
see http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20070017872_2007014421.pdf
for (another?) way to interpolate quaternions
"""
assert isinstance(q1, Quaternion) and isinstance(q2, Quaternion) assert isinstance(q1, Quaternion) and isinstance(q2, Quaternion)
Q = cls() Q = cls()
@ -522,11 +548,11 @@ class Quaternion:
# ****************************************************************************************** # ******************************************************************************************
class Symmetry: class Symmetry:
# ******************************************************************************************
lattices = [None,'orthorhombic','tetragonal','hexagonal','cubic',] lattices = [None,'orthorhombic','tetragonal','hexagonal','cubic',]
def __init__(self, symmetry = None): def __init__(self, symmetry = None):
"""lattice with given symmetry, defaults to None"""
if isinstance(symmetry, basestring) and symmetry.lower() in Symmetry.lattices: if isinstance(symmetry, basestring) and symmetry.lower() in Symmetry.lattices:
self.lattice = symmetry.lower() self.lattice = symmetry.lower()
else: else:
@ -534,29 +560,31 @@ class Symmetry:
def __copy__(self): def __copy__(self):
"""copy"""
return self.__class__(self.lattice) return self.__class__(self.lattice)
copy = __copy__ copy = __copy__
def __repr__(self): def __repr__(self):
"""readbable string"""
return '%s' % (self.lattice) return '%s' % (self.lattice)
def __eq__(self, other): def __eq__(self, other):
"""equal"""
return self.lattice == other.lattice return self.lattice == other.lattice
def __neq__(self, other): def __neq__(self, other):
"""not equal"""
return not self.__eq__(other) return not self.__eq__(other)
def __cmp__(self,other): def __cmp__(self,other):
"""linear ordering"""
return cmp(Symmetry.lattices.index(self.lattice),Symmetry.lattices.index(other.lattice)) return cmp(Symmetry.lattices.index(self.lattice),Symmetry.lattices.index(other.lattice))
def symmetryQuats(self,who = []): def symmetryQuats(self,who = []):
''' """List of symmetry operations as quaternions."""
List of symmetry operations as quaternions.
'''
if self.lattice == 'cubic': if self.lattice == 'cubic':
symQuats = [ symQuats = [
[ 1.0, 0.0, 0.0, 0.0 ], [ 1.0, 0.0, 0.0, 0.0 ],
@ -629,18 +657,15 @@ class Symmetry:
def equivalentQuaternions(self, def equivalentQuaternions(self,
quaternion, quaternion,
who = []): who = []):
''' """List of symmetrically equivalent quaternions based on own symmetry."""
List of symmetrically equivalent quaternions based on own symmetry.
'''
return [quaternion*q for q in self.symmetryQuats(who)] return [quaternion*q for q in self.symmetryQuats(who)]
def inFZ(self,R): def inFZ(self,R):
''' """Check whether given Rodrigues vector falls into fundamental zone of own symmetry."""
Check whether given Rodrigues vector falls into fundamental zone of own symmetry.
'''
if isinstance(R, Quaternion): R = R.asRodrigues() # translate accidentially passed quaternion if isinstance(R, Quaternion): R = R.asRodrigues() # translate accidentially passed quaternion
R = abs(R) # fundamental zone in Rodrigues space is point symmetric around origin # fundamental zone in Rodrigues space is point symmetric around origin
R = abs(R)
if self.lattice == 'cubic': if self.lattice == 'cubic':
return math.sqrt(2.0)-1.0 >= R[0] \ return math.sqrt(2.0)-1.0 >= R[0] \
and math.sqrt(2.0)-1.0 >= R[1] \ and math.sqrt(2.0)-1.0 >= R[1] \
@ -662,12 +687,13 @@ class Symmetry:
def inDisorientationSST(self,R): def inDisorientationSST(self,R):
''' """
Check whether given Rodrigues vector (of misorientation) falls into standard stereographic triangle of own symmetry. Check whether given Rodrigues vector (of misorientation) falls into standard stereographic triangle of own symmetry.
Determination of disorientations follow the work of A. Heinz and P. Neumann: Determination of disorientations follow the work of A. Heinz and P. Neumann:
Representation of Orientation and Disorientation Data for Cubic, Hexagonal, Tetragonal and Orthorhombic Crystals Representation of Orientation and Disorientation Data for Cubic, Hexagonal, Tetragonal and Orthorhombic Crystals
Acta Cryst. (1991). A47, 780-789 Acta Cryst. (1991). A47, 780-789
''' """
if isinstance(R, Quaternion): R = R.asRodrigues() # translate accidentially passed quaternion if isinstance(R, Quaternion): R = R.asRodrigues() # translate accidentially passed quaternion
epsilon = 0.0 epsilon = 0.0
@ -691,11 +717,12 @@ class Symmetry:
vector, vector,
proper = False, proper = False,
color = False): color = False):
''' """
Check whether given vector falls into standard stereographic triangle of own symmetry. Check whether given vector falls into standard stereographic triangle of own symmetry.
proper considers only vectors with z >= 0, hence uses two neighboring SSTs. proper considers only vectors with z >= 0, hence uses two neighboring SSTs.
Return inverse pole figure color if requested. Return inverse pole figure color if requested.
''' """
# basis = {'cubic' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red # basis = {'cubic' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
# [1.,0.,1.]/np.sqrt(2.), # direction of green # [1.,0.,1.]/np.sqrt(2.), # direction of green
# [1.,1.,1.]/np.sqrt(3.)]).transpose()), # direction of blue # [1.,1.,1.]/np.sqrt(3.)]).transpose()), # direction of blue
@ -752,15 +779,15 @@ class Symmetry:
inSST = np.all(theComponents >= 0.0) inSST = np.all(theComponents >= 0.0)
else: else:
v = np.array(vector,dtype = float) v = np.array(vector,dtype = float)
if proper: # check both improper ... if proper: # check both improper ...
theComponents = np.dot(basis['improper'],v) theComponents = np.dot(basis['improper'],v)
inSST = np.all(theComponents >= 0.0) inSST = np.all(theComponents >= 0.0)
if not inSST: # ... and proper SST if not inSST: # ... and proper SST
theComponents = np.dot(basis['proper'],v) theComponents = np.dot(basis['proper'],v)
inSST = np.all(theComponents >= 0.0) inSST = np.all(theComponents >= 0.0)
else: else:
v[2] = abs(v[2]) # z component projects identical for positive and negative values v[2] = abs(v[2]) # z component projects identical
theComponents = np.dot(basis['improper'],v) theComponents = np.dot(basis['improper'],v) # for positive and negative values
inSST = np.all(theComponents >= 0.0) inSST = np.all(theComponents >= 0.0)
if color: # have to return color array if color: # have to return color array
@ -781,7 +808,6 @@ class Symmetry:
# ****************************************************************************************** # ******************************************************************************************
class Orientation: class Orientation:
# ******************************************************************************************
__slots__ = ['quaternion','symmetry'] __slots__ = ['quaternion','symmetry']
@ -791,7 +817,7 @@ class Orientation:
angleAxis = None, angleAxis = None,
matrix = None, matrix = None,
Eulers = None, Eulers = None,
random = False, # put any integer to have a fixed seed or True for real random random = False, # integer to have a fixed seed or True for real random
symmetry = None, symmetry = None,
): ):
if random: # produce random orientation if random: # produce random orientation
@ -815,12 +841,14 @@ class Orientation:
self.symmetry = Symmetry(symmetry) self.symmetry = Symmetry(symmetry)
def __copy__(self): def __copy__(self):
"""copy"""
return self.__class__(quaternion=self.quaternion,symmetry=self.symmetry.lattice) return self.__class__(quaternion=self.quaternion,symmetry=self.symmetry.lattice)
copy = __copy__ copy = __copy__
def __repr__(self): def __repr__(self):
"""value as all implemented representations"""
return 'Symmetry: %s\n' % (self.symmetry) + \ return 'Symmetry: %s\n' % (self.symmetry) + \
'Quaternion: %s\n' % (self.quaternion) + \ 'Quaternion: %s\n' % (self.quaternion) + \
'Matrix:\n%s\n' % ( '\n'.join(['\t'.join(map(str,self.asMatrix()[i,:])) for i in range(3)]) ) + \ 'Matrix:\n%s\n' % ( '\n'.join(['\t'.join(map(str,self.asMatrix()[i,:])) for i in range(3)]) ) + \
@ -863,10 +891,7 @@ class Orientation:
self.equivalentQuaternions(who)) self.equivalentQuaternions(who))
def reduced(self): def reduced(self):
''' """Transform orientation to fall into fundamental zone according to symmetry"""
Transform orientation to fall into fundamental zone according to symmetry
'''
for me in self.symmetry.equivalentQuaternions(self.quaternion): for me in self.symmetry.equivalentQuaternions(self.quaternion):
if self.symmetry.inFZ(me.asRodrigues()): break if self.symmetry.inFZ(me.asRodrigues()): break
@ -876,13 +901,13 @@ class Orientation:
def disorientation(self, def disorientation(self,
other, other,
SST = True): SST = True):
''' """
Disorientation between myself and given other orientation. Disorientation between myself and given other orientation.
Rotation axis falls into SST if SST == True. Rotation axis falls into SST if SST == True.
(Currently requires same symmetry for both orientations. (Currently requires same symmetry for both orientations.
Look into A. Heinz and P. Neumann 1991 for cases with differing sym.) Look into A. Heinz and P. Neumann 1991 for cases with differing sym.)
''' """
if self.symmetry != other.symmetry: raise TypeError('disorientation between different symmetry classes not supported yet.') if self.symmetry != other.symmetry: raise TypeError('disorientation between different symmetry classes not supported yet.')
misQ = self.quaternion.conjugated()*other.quaternion misQ = self.quaternion.conjugated()*other.quaternion
@ -900,32 +925,27 @@ class Orientation:
if breaker: break if breaker: break
if breaker: break if breaker: break
# disorientation, own sym, other sym, self-->other: True, self<--other: False
return (Orientation(quaternion = theQ,symmetry = self.symmetry.lattice), return (Orientation(quaternion = theQ,symmetry = self.symmetry.lattice),
i,j,k == 1) # disorientation, own sym, other sym, self-->other: True, self<--other: False i,j,k == 1)
def inversePole(self, def inversePole(self,
axis, axis,
proper = False, proper = False,
SST = True): SST = True):
''' """axis rotated according to orientation (using crystal symmetry to ensure location falls into SST)"""
axis rotated according to orientation (using crystal symmetry to ensure location falls into SST)
'''
if SST: # pole requested to be within SST if SST: # pole requested to be within SST
for i,q in enumerate(self.symmetry.equivalentQuaternions(self.quaternion)): # test all symmetric equivalent quaternions for i,q in enumerate(self.symmetry.equivalentQuaternions(self.quaternion)): # test all symmetric equivalent quaternions
pole = q.conjugated()*axis # align crystal direction to axis pole = q.conjugated()*axis # align crystal direction to axis
if self.symmetry.inSST(pole,proper): break # found SST version if self.symmetry.inSST(pole,proper): break # found SST version
else: else:
pole = self.quaternion.conjugated()*axis # align crystal direction to axis pole = self.quaternion.conjugated()*axis # align crystal direction to axis
return (pole,i if SST else 0) return (pole,i if SST else 0)
def IPFcolor(self,axis): def IPFcolor(self,axis):
''' """TSL color of inverse pole figure for given axis"""
TSL color of inverse pole figure for given axis
'''
color = np.zeros(3,'d') color = np.zeros(3,'d')
for q in self.symmetry.equivalentQuaternions(self.quaternion): for q in self.symmetry.equivalentQuaternions(self.quaternion):
@ -939,7 +959,9 @@ class Orientation:
def average(cls, def average(cls,
orientations, orientations,
multiplicity = []): multiplicity = []):
"""RETURN THE AVERAGE ORIENTATION """
average orientation
ref: F. Landis Markley, Yang Cheng, John Lucas Crassidis, and Yaakov Oshman. ref: F. Landis Markley, Yang Cheng, John Lucas Crassidis, and Yaakov Oshman.
Averaging Quaternions, Averaging Quaternions,
Journal of Guidance, Control, and Dynamics, Vol. 30, No. 4 (2007), pp. 1193-1197. Journal of Guidance, Control, and Dynamics, Vol. 30, No. 4 (2007), pp. 1193-1197.
@ -949,7 +971,6 @@ class Orientation:
b = Orientation(Eulers=np.radians([20, 0, 0]), symmetry='hexagonal') b = Orientation(Eulers=np.radians([20, 0, 0]), symmetry='hexagonal')
avg = Orientation.average([a,b]) avg = Orientation.average([a,b])
""" """
if not all(isinstance(item, Orientation) for item in orientations): if not all(isinstance(item, Orientation) for item in orientations):
raise TypeError("Only instances of Orientation can be averaged.") raise TypeError("Only instances of Orientation can be averaged.")
@ -960,8 +981,7 @@ class Orientation:
reference = orientations[0] # take first as reference reference = orientations[0] # take first as reference
for i,(o,n) in enumerate(zip(orientations,multiplicity)): for i,(o,n) in enumerate(zip(orientations,multiplicity)):
closest = o.equivalentOrientations(reference.disorientation(o,SST = False)[2])[0] # select sym orientation with lowest misorientation closest = o.equivalentOrientations(reference.disorientation(o,SST = False)[2])[0] # select sym orientation with lowest misorientation
M = closest.quaternion.asM() * n if i == 0 else M + closest.quaternion.asM() * n # add (multiples) of this orientation to average M = closest.quaternion.asM() * n if i == 0 else M + closest.quaternion.asM() * n # noqa add (multiples) of this orientation to average noqa
eig, vec = np.linalg.eig(M/N) eig, vec = np.linalg.eig(M/N)
return Orientation(quaternion = Quaternion(quatArray = np.real(vec.T[eig.argmax()])), return Orientation(quaternion = Quaternion(quatArray = np.real(vec.T[eig.argmax()])),

9
lib/damask/result.py
View File

@ -11,10 +11,11 @@ except (ImportError) as e:
sys.stderr.write('\nREMARK: h5py module not available \n\n') sys.stderr.write('\nREMARK: h5py module not available \n\n')
class Result(): class Result():
''' """
General class for result parsing. General class for result parsing.
Needs h5py to be installed
''' Needs h5py to be installed
"""
def __init__(self,resultsFile): def __init__(self,resultsFile):
self.data=h5py.File(resultsFile,"r") self.data=h5py.File(resultsFile,"r")

149
lib/damask/result/marc2vtk.py
View File

@ -1,6 +1,5 @@
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
# $Id$
# This tool converts a msc.marc result file into the vtk format that # This tool converts a msc.marc result file into the vtk format that
# can be viewed by Paraview software (Kitware), or MayaVi (needs xml-vtk, or ... # can be viewed by Paraview software (Kitware), or MayaVi (needs xml-vtk, or ...
# #
@ -8,13 +7,8 @@
# Some example vtk files: http://people.sc.fsu.edu/~jburkardt/data/vtk/vtk.html # Some example vtk files: http://people.sc.fsu.edu/~jburkardt/data/vtk/vtk.html
# www.paraview.org # www.paraview.org
import os,sys,math,time,re import os,sys,re
# python external import numpy as np
try:
import numpy as N
import numpy
except:
print('Could not import numpy.')
import py_post # MSC closed source module to access marc result files import py_post # MSC closed source module to access marc result files
@ -27,7 +21,7 @@ class MARC_POST():
self.fpath=os.path.join(self.projdir,self.postname) self.fpath=os.path.join(self.projdir,self.postname)
print('Trying to open ',self.fpath,' ...') print('Trying to open ',self.fpath,' ...')
self.p=py_post.post_open(self.fpath) self.p=py_post.post_open(self.fpath)
if self.p==None: if self.p is None:
print('Could not open %s.'%self.postname); #return 'err'#; sys.exit(1) print('Could not open %s.'%self.postname); #return 'err'#; sys.exit(1)
raise Exception('Could not open t16') raise Exception('Could not open t16')
print('Postfile %s%s is open ...'%(self.projdir,self.postname)) print('Postfile %s%s is open ...'%(self.projdir,self.postname))
@ -105,7 +99,6 @@ class MARC_POST():
def writeNodes2VTK(self, fobj): def writeNodes2VTK(self, fobj):
self.points=[] self.points=[]
self.VTKcnt=200 # number of values per line in vtk file self.VTKcnt=200 # number of values per line in vtk file
ndCnt=1
fobj.write('POINTS %i'%self.p.nodes()+' float\n') fobj.write('POINTS %i'%self.p.nodes()+' float\n')
self.nodes_dict={} # store the node IDs in case of holes in the numbering self.nodes_dict={} # store the node IDs in case of holes in the numbering
for iNd in self.nodes: for iNd in self.nodes:
@ -126,8 +119,6 @@ class MARC_POST():
el=self.p.element(iEl) el=self.p.element(iEl)
cell_nodes=[] # for pyvtk cell_nodes=[] # for pyvtk
ndlist=el.items ndlist=el.items
#for k in [0, 1, 3, 2, 4, 5, 7, 6]: # FOR CELL TPYE VTK_VOXEL
#for k in [0, 4, 3, 1, 5, 7, 6, 2]:
for k in [0, 1, 2, 3, 4, 5, 6, 7]: # FOR CELL TYPE VTK_HEXAHEDRON for k in [0, 1, 2, 3, 4, 5, 6, 7]: # FOR CELL TYPE VTK_HEXAHEDRON
node=ndlist[k]-1 node=ndlist[k]-1
cell_nodes.append(self.nodes_dict[node]) cell_nodes.append(self.nodes_dict[node])
@ -147,7 +138,6 @@ class MARC_POST():
fobj.write('\n');cnt=0 fobj.write('\n');cnt=0
fobj.write('\n') fobj.write('\n')
print('Elements written to VTK: %i'%self.p.elements()) print('Elements written to VTK: %i'%self.p.elements())
#print('Nr of nodes: ',self.nodes)
def writeElScalars2NodesVTK(self,fobj): def writeElScalars2NodesVTK(self,fobj):
fobj.write('\nPOINT_DATA %i\n'%self.p.nodes()) fobj.write('\nPOINT_DATA %i\n'%self.p.nodes())
@ -157,7 +147,6 @@ class MARC_POST():
fobj.write('LOOKUP_TABLE default\n') fobj.write('LOOKUP_TABLE default\n')
idxScal=self.nscal_list.index('Displacement Z') idxScal=self.nscal_list.index('Displacement Z')
for iNd in self.nodes: for iNd in self.nodes:
#fobj.write('%f %f '%(self.p.node_scalar(iNd,idxScal), N.random.rand()))
fobj.write('%f '%(self.p.node_scalar(iNd,idxScal))) fobj.write('%f '%(self.p.node_scalar(iNd,idxScal)))
for iEl in range(0,self.nel): for iEl in range(0,self.nel):
el=self.p.element(iEl) el=self.p.element(iEl)
@ -173,8 +162,6 @@ class MARC_POST():
def writeNodeScalars2VTK(self,fobj): def writeNodeScalars2VTK(self,fobj):
fobj.write('\nPOINT_DATA %i\n'%self.p.nodes()) fobj.write('\nPOINT_DATA %i\n'%self.p.nodes())
nNdDat=self.nscals
nComponents=1+nNdDat
self.pointDataScalars=[] self.pointDataScalars=[]
for idxNdScal in range(-3,self.nscals): #now include node x,y,z for idxNdScal in range(-3,self.nscals): #now include node x,y,z
if idxNdScal>=0: if idxNdScal>=0:
@ -209,8 +196,6 @@ class MARC_POST():
idx_sig_vMises=self.getLabelNr('Equivalent Von Mises Stress') idx_sig_vMises=self.getLabelNr('Equivalent Von Mises Stress')
idx_sig33=self.getLabelNr('Comp 33 of Cauchy Stress') idx_sig33=self.getLabelNr('Comp 33 of Cauchy Stress')
fobj.write('\nCELL_DATA %i\n'%self.p.elements()) fobj.write('\nCELL_DATA %i\n'%self.p.elements())
nElDat=self.elscals
nComponents=1+nElDat
for idxElScal in range(0,self.elscals): for idxElScal in range(0,self.elscals):
datalabel=self.elscal_list[idxElScal] datalabel=self.elscal_list[idxElScal]
datalabel=re.sub("\s",'_',datalabel) datalabel=re.sub("\s",'_',datalabel)
@ -251,18 +236,15 @@ class MARC_POST():
return result return result
def writeUniaxiality2VTK(self,fobj): def writeUniaxiality2VTK(self,fobj):
#fobj.write('\nCELL_DATA %i\n'%self.p.elements())
datalabel='uniaxiality_sig_vMises_durch_sig33' datalabel='uniaxiality_sig_vMises_durch_sig33'
fobj.write('SCALARS %s float %i\n'%(datalabel,1)) fobj.write('SCALARS %s float %i\n'%(datalabel,1))
fobj.write('LOOKUP_TABLE default\n') fobj.write('LOOKUP_TABLE default\n')
cnt=0 cnt=0
for iEl in range(0,self.nel): for iEl in range(0,self.nel):
cnt=cnt+1 cnt=cnt+1
#if abs(self.sig33[iEl])<1e-5:
if abs(self.sig_vMises[iEl])<1e-5: if abs(self.sig_vMises[iEl])<1e-5:
datum=0. datum=0.
else: else:
#datum=self.sig_vMises[iEl]/self.sig33[iEl]
datum=self.sig33[iEl]/self.sig_vMises[iEl] datum=self.sig33[iEl]/self.sig_vMises[iEl]
fobj.write('%E '%(datum)) fobj.write('%E '%(datum))
if cnt>self.VTKcnt: if cnt>self.VTKcnt:
@ -283,8 +265,8 @@ class MARC_POST():
self.mean_stress.append(self.meanStress(sig)) self.mean_stress.append(self.meanStress(sig))
def triaxiality_per_element(self): def triaxiality_per_element(self):
# classical triaxiality # classical triaxiality
# 1/3 : uniax tension # 1/3 : uniax tension
self.triaxiality=[] self.triaxiality=[]
for iEl in range(0,self.nel): for iEl in range(0,self.nel):
t=self.mean_stress[iEl]/self.sig_vMises[iEl] t=self.mean_stress[iEl]/self.sig_vMises[iEl]
@ -303,10 +285,6 @@ class MARC_POST():
fobj.write('\n') fobj.write('\n')
def calc_lode_parameter(self): def calc_lode_parameter(self):
# [-1 ... +1] see e.g. Wippler & Boehlke triaxiality measures doi:10.1002/pamm.201010061
# +1 : uniax tensile?
# 0 : shear
# -1 : uniax compr ?
self.lode=[] self.lode=[]
try: try:
self.stress self.stress
@ -328,10 +306,11 @@ class MARC_POST():
def princStress(self, stress): def princStress(self, stress):
""" """
Function to compute 3D principal stresses and sort them. Function to compute 3D principal stresses and sort them.
from: http://geodynamics.org/svn/cig/short/3D/PyLith/trunk/playpen/postproc/vtkcff.py from: http://geodynamics.org/svn/cig/short/3D/PyLith/trunk/playpen/postproc/vtkcff.py
""" """
stressMat=N.array(stress) stressMat=np.array(stress)
(princStress, princAxes) = numpy.linalg.eigh(stressMat) (princStress, princAxes) = np.linalg.eigh(stressMat)
idx = princStress.argsort() idx = princStress.argsort()
princStressOrdered = princStress[idx] princStressOrdered = princStress[idx]
princAxesOrdered = princAxes[:,idx] princAxesOrdered = princAxes[:,idx]
@ -339,36 +318,28 @@ class MARC_POST():
def avg_elten(self, def avg_elten(self,
idxElTen, mat=0, elID=None): idxElTen, mat=0, elID=None):
tensum=N.zeros((3,3)); tensum=np.zeros((3,3));
T=N.zeros((3,3)); T=np.zeros((3,3));
pts=0; pts=0;
avg=N.zeros((3,3)); avg=np.zeros((3,3));
#print 'Element Scalars'
#print self.p.element_scalar_label(elscal2) if elID is None:
if elID==None:
averaged_elements=range(0,self.nel) averaged_elements=range(0,self.nel)
else: else:
averaged_elements=[elID] averaged_elements=[elID]
#for i in range (0,self.nel):
for i in averaged_elements: for i in averaged_elements:
if mat==0 or int(self.p.element_scalar(i,4)[0].value)==mat: if mat==0 or int(self.p.element_scalar(i,4)[0].value)==mat:
eldata=self.p.element(i)
T=self.p.element_tensor(i,idxElTen) T=self.p.element_tensor(i,idxElTen)
for k in range (0,8): for k in range (0,8):
tensum[0][0] = tensum[0][0] + T[k].t11 tensum[0][0] = tensum[0][0] + T[k].t11
tensum[0][1] = tensum[0][1] + T[k].t12 tensum[0][1] = tensum[0][1] + T[k].t12
tensum[0][2] = tensum[0][2] + T[k].t13 tensum[0][2] = tensum[0][2] + T[k].t13
#tensum1[1][0] = tensum1[1][0] + T1[k].t21
tensum[1][1] = tensum[1][1] + T[k].t22 tensum[1][1] = tensum[1][1] + T[k].t22
tensum[1][2] = tensum[1][2] + T[k].t23 tensum[1][2] = tensum[1][2] + T[k].t23
#tensum1[2][0] = tensum1[2][0] + T1[k].t31
#tensum1[2][1] = tensum1[2][1] + T1[k].t32
tensum[2][2] = tensum[2][2] + T[k].t33 tensum[2][2] = tensum[2][2] + T[k].t33
pts=pts+1 pts=pts+1
avg=tensum/pts avg=tensum/pts
#print avg
avg=self.fillComponents(avg) avg=self.fillComponents(avg)
#print avg
del [T] del [T]
return (avg,tensum,pts) return (avg,tensum,pts)
@ -384,7 +355,7 @@ class MARC_POST():
t=tensor33 t=tensor33
s=(t[0,0]-t[1,1])**2+(t[1,1]-t[2,2])**2+(t[0,0]-t[2,2])**2+\ s=(t[0,0]-t[1,1])**2+(t[1,1]-t[2,2])**2+(t[0,0]-t[2,2])**2+\
6*(t[0,1]**2+t[1,2]**2+t[2,0]**2) 6*(t[0,1]**2+t[1,2]**2+t[2,0]**2)
vM=N.sqrt(s/2.) vM=np.sqrt(s/2.)
return vM return vM
def meanStress(self,tensor33): def meanStress(self,tensor33):
@ -398,7 +369,6 @@ class MARC_POST():
I1=t[0,0]+t[1,1]+t[2,2] I1=t[0,0]+t[1,1]+t[2,2]
I2=t[0,0]*t[1,1]+t[1,1]*t[2,2]+t[0,0]*t[2,2]-\ I2=t[0,0]*t[1,1]+t[1,1]*t[2,2]+t[0,0]*t[2,2]-\
t[0,1]**2-t[1,2]**2-t[0,2]**2 t[0,1]**2-t[1,2]**2-t[0,2]**2
# I3 = det(t)
I3=t[0,0]*t[1,1]*t[2,2]+\ I3=t[0,0]*t[1,1]*t[2,2]+\
2*t[0,1]*t[1,2]*t[2,0]-\ 2*t[0,1]*t[1,2]*t[2,0]-\
t[2,2]*t[0,1]**2-t[0,0]*t[1,2]**2-t[1,1]*t[0,2]**2 t[2,2]*t[0,1]**2-t[0,0]*t[1,2]**2-t[1,1]*t[0,2]**2
@ -406,17 +376,18 @@ class MARC_POST():
class VTK_WRITER(): class VTK_WRITER():
''' """
The resulting vtk-file can be imported in Paraview 3.12 The resulting vtk-file can be imported in Paraview 3.12
Then use Filters: Cell Data to Point Data + Contour
to plot semi-transparent iso-surfaces. Then use Filters: Cell Data to Point Data + Contour
''' to plot semi-transparent iso-surfaces.
"""
import re import re
def __init__(self): def __init__(self):
self.p=MARC_POST() # self.p self.p=MARC_POST() # self.p
def openFile(self, filename='test.vtp'): def openFile(self, filename='test.vtp'):
#if not self.f:#==None:
self.f=open(filename,'w+') self.f=open(filename,'w+')
self.fname=filename self.fname=filename
@ -427,7 +398,7 @@ class VTK_WRITER():
dformat='ASCII', # BINARY | [ASCII] dformat='ASCII', # BINARY | [ASCII]
dtype='UNSTRUCTURED_GRID' # UNSTRUCTURED GRID dtype='UNSTRUCTURED_GRID' # UNSTRUCTURED GRID
): ):
if vtkFile==None: if vtkFile is None:
vtkFile=self.f vtkFile=self.f
# First Line contains Data format version # First Line contains Data format version
self.versionVTK=version self.versionVTK=version
@ -440,7 +411,6 @@ class VTK_WRITER():
def marc2vtkBatch(self): def marc2vtkBatch(self):
for iori in range(1,63): for iori in range(1,63):
#self.p=msc_post.MSC_POST()
self.p.postname='indent_fric0.3_R2.70_cA146.0_h0.320_ori%03i_OST_h19d.t16'%(iori) self.p.postname='indent_fric0.3_R2.70_cA146.0_h0.320_ori%03i_OST_h19d.t16'%(iori)
if os.path.exists(self.p.postname): if os.path.exists(self.p.postname):
self.marc2vtk(mode='fast', batchMode=1) self.marc2vtk(mode='fast', batchMode=1)
@ -496,14 +466,14 @@ class VTK_WRITER():
def scaleBar(self, length=1.0, posXYZ=[0., 0., 0.]): def scaleBar(self, length=1.0, posXYZ=[0., 0., 0.]):
self.fsb=open('micronbar_l%.1f.vtp'%length,'w+') self.fsb=open('micronbar_l%.1f.vtp'%length,'w+')
self.writeFirstLines(self.fsb, comment='micronbar') self.writeFirstLines(self.fsb, comment='micronbar')
pts=N.array([]) pts=np.array([])
width=length*1. width=length*1.
height=length*1. height=length*1.
wVec=N.array([0., width, 0.]) wVec=np.array([0., width, 0.])
lVec=N.array([length,0.,0.]) lVec=np.array([length,0.,0.])
hVec=N.array([0.,0.,height]) hVec=np.array([0.,0.,height])
posXYZ=posXYZ-0.5*wVec-0.5*lVec#-0.5*hVec # CENTERING Y/N posXYZ=posXYZ-0.5*wVec-0.5*lVec#-0.5*hVec # CENTERING Y/N
posXYZ=N.array(posXYZ) posXYZ=np.array(posXYZ)
pts=[posXYZ, posXYZ+lVec, pts=[posXYZ, posXYZ+lVec,
posXYZ+wVec, posXYZ+wVec,
posXYZ+wVec+lVec] posXYZ+wVec+lVec]
@ -514,34 +484,22 @@ class VTK_WRITER():
self.fsb.write('%f %f %f\n'%(pts[npts][0], pts[npts][1], pts[npts][2])) self.fsb.write('%f %f %f\n'%(pts[npts][0], pts[npts][1], pts[npts][2]))
if 1: #Triad if 1: #Triad
nCells=3 nCells=3
#nCells=1 #One Line
ptsPerCell=2 # Lines (Type=3) ptsPerCell=2 # Lines (Type=3)
#ptsPerCell=4 # Quads (Type=9)
#ptsPerCell=8 # Hexahedron (Type=12)
cellSize=(ptsPerCell+1)*nCells cellSize=(ptsPerCell+1)*nCells
self.fsb.write('CELLS %i %i\n'%(nCells,cellSize)) self.fsb.write('CELLS %i %i\n'%(nCells,cellSize))
self.fsb.write('2 0 1\n') #X-Line self.fsb.write('2 0 1\n') #X-Line
self.fsb.write('2 0 2\n') #Y-Line self.fsb.write('2 0 2\n') #Y-Line
self.fsb.write('2 0 4\n') #Z-Line self.fsb.write('2 0 4\n') #Z-Line
#self.fsb.write('4 0 1 3 2\n') #Quad
#self.fsb.write('%i 0 1 3 2 4 5 7 6\n'%ptsPerCell) #Hexahedron
self.fsb.write('CELL_TYPES %i\n'%(nCells)) self.fsb.write('CELL_TYPES %i\n'%(nCells))
self.fsb.write('3\n3\n3\n')#Line self.fsb.write('3\n3\n3\n')#Line
#self.fsb.write('12\n')#Hexahedron
else: # Cube, change posXYZ else: # Cube, change posXYZ
nCells=1 nCells=1
ptsPerCell=2 # Lines (Type=3) ptsPerCell=2 # Lines (Type=3)
#ptsPerCell=4 # Quads (Type=9)
#ptsPerCell=8 # Hexahedron (Type=12)
cellSize=(ptsPerCell+1)*nCells cellSize=(ptsPerCell+1)*nCells
self.fsb.write('CELLS %i %i\n'%(nCells,cellSize)) self.fsb.write('CELLS %i %i\n'%(nCells,cellSize))
self.fsb.write('2 0 1\n') #Line self.fsb.write('2 0 1\n') #Line
#self.fsb.write('4 0 1 3 2\n') #Quad
#self.fsb.write('%i 0 1 3 2 4 5 7 6\n'%ptsPerCell) #Hexahedron
self.fsb.write('CELL_TYPES %i\n'%(nCells)) self.fsb.write('CELL_TYPES %i\n'%(nCells))
self.fsb.write('3\n')#Line self.fsb.write('3\n')#Line
#self.fsb.write('12\n')#Hexahedron
self.fsb.write('\n') self.fsb.write('\n')
self.fsb.close() self.fsb.close()
@ -549,8 +507,7 @@ class VTK_WRITER():
def example_unstructured(self): def example_unstructured(self):
self.openFile(filename='example_unstructured_grid.vtk') self.openFile(filename='example_unstructured_grid.vtk')
#self.writeFirstLines() self.f.write("""
self.f.write('''
# vtk DataFile Version 2.0 # vtk DataFile Version 2.0
example_unstruct_grid example_unstruct_grid
ASCII ASCII
@ -590,61 +547,40 @@ LOOKUP_TABLE default
1.02 1.02
1.50 1.50
0.00 0.00
3 5 6 23423423423423423423.23423423''') 3 5 6 23423423423423423423.23423423""")
self.f.close() self.f.close()
def writeNodes2VTK(self, fobj): def writeNodes2VTK(self, fobj):
self.VTKcnt=200 # how many numbers per line in vtk file self.VTKcnt=200 # how many numbers per line in vtk file
#self.VTKcnt=6
ndCnt=1
#self.nodes=range(0,10)
fobj.write('POINTS %i'%self.p.nodes()+' float\n') fobj.write('POINTS %i'%self.p.nodes()+' float\n')
for iNd in self.nodes: for iNd in self.nodes:
nd=self.p.node(iNd) nd=self.p.node(iNd)
disp=self.p.node_displacement(iNd) disp=self.p.node_displacement(iNd)
#contact=self.p.node_scalar(iNd,contactNr)
#ndCnt=ndCnt+1
fobj.write('%f %f %f \n'% fobj.write('%f %f %f \n'%
#(nd.x, nd.y, nd.z))
(nd.x+disp[0], nd.y+disp[1], nd.z+disp[2])) (nd.x+disp[0], nd.y+disp[1], nd.z+disp[2]))
#if ndCnt>6:
# fobj.write('\n')
# ndCnt=1
fobj.write('\n') fobj.write('\n')
print('Nodes written to VTK: %i'%self.p.nodes()) print('Nodes written to VTK: %i'%self.p.nodes())
#print('Nr of nodes: ',self.nodes)
def writeElements2VTK(self, fobj): def writeElements2VTK(self, fobj):
fobj.write('\nCELLS %i %i'%(self.p.elements(),self.p.elements()*9)+'\n') fobj.write('\nCELLS %i %i'%(self.p.elements(),self.p.elements()*9)+'\n')
for iEl in range(0,self.nel): for iEl in range(0,self.nel):
el=self.p.element(iEl) el=self.p.element(iEl)
#disp=self.p.node_displacement(iNd)
#contact=self.p.node_scalar(iNd,contactNr)
#ndCnt=ndCnt+1
fobj.write('8 ') fobj.write('8 ')
ndlist=el.items ndlist=el.items
#for k in [0, 1, 3, 2, 4, 5, 7, 6]: # FOR CELL TPYE VTK_VOXEL
#for k in [0, 4, 3, 1, 5, 7, 6, 2]:
for k in [0, 1, 2, 3, 4, 5, 6, 7]: # FOR CELL TYPE VTK_HEXAHEDRON for k in [0, 1, 2, 3, 4, 5, 6, 7]: # FOR CELL TYPE VTK_HEXAHEDRON
fobj.write('%6i '%(ndlist[k]-1)) fobj.write('%6i '%(ndlist[k]-1))
fobj.write('\n') fobj.write('\n')
#if ndCnt>6:
# fobj.write('\n')
# ndCnt=1
fobj.write('\nCELL_TYPES %i'%self.p.elements()+'\n') fobj.write('\nCELL_TYPES %i'%self.p.elements()+'\n')
cnt=0 cnt=0
for iEl in range(0,self.nel): for iEl in range(0,self.nel):
cnt=cnt+1 cnt=cnt+1
#fobj.write('11\n') #VTK_VOXEL
fobj.write('12 ') #VTK_HEXAHEDRON fobj.write('12 ') #VTK_HEXAHEDRON
if cnt>self.VTKcnt: if cnt>self.VTKcnt:
fobj.write('\n');cnt=0 fobj.write('\n');cnt=0
fobj.write('\n') fobj.write('\n')
print('Elements written to VTK: %i'%self.p.elements()) print('Elements written to VTK: %i'%self.p.elements())
#print('Nr of nodes: ',self.nodes)
def writeElScalars2NodesVTK(self,fobj): def writeElScalars2NodesVTK(self,fobj):
fobj.write('\nPOINT_DATA %i\n'%self.p.nodes()) fobj.write('\nPOINT_DATA %i\n'%self.p.nodes())
@ -668,10 +604,7 @@ LOOKUP_TABLE default
fobj.write('\n') fobj.write('\n')
def writeNodeScalars2VTK(self,fobj): def writeNodeScalars2VTK(self,fobj):
#print('writeElementData2VTK')
fobj.write('\nPOINT_DATA %i\n'%self.p.nodes()) fobj.write('\nPOINT_DATA %i\n'%self.p.nodes())
nNdDat=self.nscals
nComponents=1+nNdDat
for idxNdScal in range(-3,self.nscals): # include node x,y,z for idxNdScal in range(-3,self.nscals): # include node x,y,z
if idxNdScal>=0: if idxNdScal>=0:
datalabel=self.nscal_list[idxNdScal] datalabel=self.nscal_list[idxNdScal]
@ -700,10 +633,7 @@ LOOKUP_TABLE default
fobj.write('\n') fobj.write('\n')
def writeElementData2VTK(self,fobj): def writeElementData2VTK(self,fobj):
#print('writeElementData2VTK')
fobj.write('\nCELL_DATA %i\n'%self.p.elements()) fobj.write('\nCELL_DATA %i\n'%self.p.elements())
nElDat=self.elscals
nComponents=1+nElDat
for idxElScal in range(0,self.elscals): for idxElScal in range(0,self.elscals):
datalabel=self.elscal_list[idxElScal] datalabel=self.elscal_list[idxElScal]
datalabel=re.sub("\s",'_',datalabel) datalabel=re.sub("\s",'_',datalabel)
@ -730,7 +660,7 @@ LOOKUP_TABLE default
def example1(self): def example1(self):
self.openFile() self.openFile()
self.writeFirstLines() self.writeFirstLines()
self.f.write('''DATASET POLYDATA self.f.write("""DATASET POLYDATA
POINTS 8 float POINTS 8 float
0.0 0.0 0.0 0.0 0.0 0.0
1.0 0.0 0.0 1.0 0.0 0.0
@ -789,18 +719,20 @@ LOOKUP_TABLE my_table 8
0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0
1.0 0.0 1.0 1.0 1.0 0.0 1.0 1.0
0.0 1.0 1.0 1.0 0.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0''') 1.0 1.0 1.0 1.0""")
self.f.close() self.f.close()
import pyvtk import pyvtk
class marc_to_vtk(): class marc_to_vtk():
''' """
Anybody wants to implement it with pyvtk? Anybody wants to implement it with pyvtk?
The advantage would be that pyvtk can also wirte the
<xml>-VTK format and binary. The advantage would be that pyvtk can also wirte the
These can be plotted with mayavi. <xml>-VTK format and binary.
''' These can be plotted with mayavi.
"""
def __init__(self): def __init__(self):
self.p=[]#MARC_POST() # self.p self.p=[]#MARC_POST() # self.p
@ -810,5 +742,4 @@ class marc_to_vtk():
hexahedron=self.p.cells), hexahedron=self.p.cells),
'm2v output') 'm2v output')
vtk.tofile('m2v_file') vtk.tofile('m2v_file')
#vtk.tofile('example3b','binary')
#VtkData('example3')

2
lib/damask/test/__init__.py
View File

@ -2,4 +2,4 @@
"""Test functionality""" """Test functionality"""
from .test import Test "noqa from .test import Test #noqa