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

Merge branch 'development' into vectorize_rotation

This commit is contained in:
Martin Diehl 2020-04-21 19:29:30 +02:00
parent 97a5880d76 f1afc159ec
commit 7efe14be35
16 changed files with 492 additions and 446 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
VERSION
View File

@ -1 +1 @@
v2.0.3-2243-gbb03483b v2.0.3-2303-g2a6132b7

2
python/damask/__init__.py
View File

@ -6,7 +6,7 @@ name = 'damask'
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 = _re.sub(r'^v','',_f.readline().strip()) version = _re.sub(r'^v','',_f.readline().strip())
# classes # make classes directly accessible as damask.Class
from ._environment import Environment # noqa from ._environment import Environment # noqa
from ._table import Table # noqa from ._table import Table # noqa
from ._vtk import VTK # noqa from ._vtk import VTK # noqa

4
python/damask/_geom.py
View File

@ -291,7 +291,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.split('#')[0].lower().strip().split()
key = items[0] if items 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']])
@ -307,7 +307,7 @@ class Geom:
microstructure = np.empty(grid.prod()) # initialize as flat array microstructure = np.empty(grid.prod()) # initialize as flat array
i = 0 i = 0
for line in content[header_length:]: for line in content[header_length:]:
items = line.split() items = line.split('#')[0].split()
if len(items) == 3: if len(items) == 3:
if items[1].lower() == 'of': if items[1].lower() == 'of':
items = np.ones(int(items[0]))*float(items[2]) items = np.ones(int(items[0]))*float(items[2])

18
python/damask/_result.py
View File

@ -68,12 +68,12 @@ class Result:
self.con_physics = [] self.con_physics = []
for c in self.constituents: for c in self.constituents:
self.con_physics += f['/'.join([self.increments[0],'constituent',c])].keys() self.con_physics += f['/'.join([self.increments[0],'constituent',c])].keys()
self.con_physics = list(set(self.con_physics)) # make unique self.con_physics = list(set(self.con_physics)) # make unique
self.mat_physics = [] self.mat_physics = []
for m in self.materialpoints: for m in self.materialpoints:
self.mat_physics += f['/'.join([self.increments[0],'materialpoint',m])].keys() self.mat_physics += f['/'.join([self.increments[0],'materialpoint',m])].keys()
self.mat_physics = list(set(self.mat_physics)) # make unique self.mat_physics = list(set(self.mat_physics)) # make unique
self.selection = {'increments': self.increments, self.selection = {'increments': self.increments,
'constituents': self.constituents,'materialpoints': self.materialpoints, 'constituents': self.constituents,'materialpoints': self.materialpoints,
@ -86,13 +86,19 @@ class Result:
def __repr__(self): def __repr__(self):
"""Show selected data.""" """Show selected data."""
all_selected_increments = self.selection['increments'] all_selected_increments = self.selection['increments']
self.pick('increments',all_selected_increments[0:1]) self.pick('increments',all_selected_increments[0:1])
first = self.list_data() first = self.list_data()
self.pick('increments',all_selected_increments[-1:]) self.pick('increments',all_selected_increments[-1:])
last = self.list_data() last = '' if len(all_selected_increments) < 2 else self.list_data()
self.pick('increments',all_selected_increments) self.pick('increments',all_selected_increments)
in_between = ''.join(['\n{}\n ...\n'.format(inc) for inc in all_selected_increments[1:-2]])
return util.srepr(first+ in_between + last) in_between = '' if len(all_selected_increments) < 3 else \
''.join(['\n{}\n ...\n'.format(inc) for inc in all_selected_increments[1:-2]])
return util.srepr(first + in_between + last)
def _manage_selection(self,action,what,datasets): def _manage_selection(self,action,what,datasets):
@ -1009,7 +1015,7 @@ class Result:
continue continue
lock.acquire() lock.acquire()
with h5py.File(self.fname, 'a') as f: with h5py.File(self.fname, 'a') as f:
try: try: # ToDo: Replace if exists?
dataset = f[result[0]].create_dataset(result[1]['label'],data=result[1]['data']) dataset = f[result[0]].create_dataset(result[1]['label'],data=result[1]['data'])
for l,v in result[1]['meta'].items(): for l,v in result[1]['meta'].items():
dataset.attrs[l]=v.encode() dataset.attrs[l]=v.encode()

142
python/damask/grid_filters.py
View File

@ -1,5 +1,5 @@
from scipy import spatial from scipy import spatial as _spatial
import numpy as np import numpy as _np
def _ks(size,grid,first_order=False): def _ks(size,grid,first_order=False):
""" """
@ -11,16 +11,16 @@ def _ks(size,grid,first_order=False):
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]
if grid[0]%2 == 0 and first_order: k_sk[grid[0]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) if grid[0]%2 == 0 and first_order: k_sk[grid[0]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/size[1] k_sj = _np.where(_np.arange(grid[1])>grid[1]//2,_np.arange(grid[1])-grid[1],_np.arange(grid[1]))/size[1]
if grid[1]%2 == 0 and first_order: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) if grid[1]%2 == 0 and first_order: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
k_si = np.arange(grid[2]//2+1)/size[2] k_si = _np.arange(grid[2]//2+1)/size[2]
kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij') kk, kj, ki = _np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij')
return np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3) return _np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3)
def curl(size,field): def curl(size,field):
@ -33,18 +33,18 @@ def curl(size,field):
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))
curl_ = (np.einsum('slm,ijkl,ijkm ->ijks', e,k_s,field_fourier)*2.0j*np.pi if n == 3 else # vector, 3 -> 3 curl_ = (_np.einsum('slm,ijkl,ijkm ->ijks', e,k_s,field_fourier)*2.0j*_np.pi if n == 3 else # vector, 3 -> 3
np.einsum('slm,ijkl,ijknm->ijksn',e,k_s,field_fourier)*2.0j*np.pi) # tensor, 3x3 -> 3x3 _np.einsum('slm,ijkl,ijknm->ijksn',e,k_s,field_fourier)*2.0j*_np.pi) # tensor, 3x3 -> 3x3
return np.fft.irfftn(curl_,axes=(0,1,2),s=field.shape[:3]) return _np.fft.irfftn(curl_,axes=(0,1,2),s=field.shape[:3])
def divergence(size,field): def divergence(size,field):
@ -57,14 +57,14 @@ def divergence(size,field):
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))
div_ = (np.einsum('ijkl,ijkl ->ijk', k_s,field_fourier)*2.0j*np.pi if n == 3 else # vector, 3 -> 1 div_ = (_np.einsum('ijkl,ijkl ->ijk', k_s,field_fourier)*2.0j*_np.pi if n == 3 else # vector, 3 -> 1
np.einsum('ijkm,ijklm->ijkl',k_s,field_fourier)*2.0j*np.pi) # tensor, 3x3 -> 3 _np.einsum('ijkm,ijklm->ijkl',k_s,field_fourier)*2.0j*_np.pi) # tensor, 3x3 -> 3
return np.fft.irfftn(div_,axes=(0,1,2),s=field.shape[:3]) return _np.fft.irfftn(div_,axes=(0,1,2),s=field.shape[:3])
def gradient(size,field): def gradient(size,field):
@ -77,17 +77,17 @@ def gradient(size,field):
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))
grad_ = (np.einsum('ijkl,ijkm->ijkm', field_fourier,k_s)*2.0j*np.pi if n == 1 else # scalar, 1 -> 3 grad_ = (_np.einsum('ijkl,ijkm->ijkm', field_fourier,k_s)*2.0j*_np.pi if n == 1 else # scalar, 1 -> 3
np.einsum('ijkl,ijkm->ijklm',field_fourier,k_s)*2.0j*np.pi) # vector, 3 -> 3x3 _np.einsum('ijkl,ijkm->ijklm',field_fourier,k_s)*2.0j*_np.pi) # vector, 3 -> 3x3
return np.fft.irfftn(grad_,axes=(0,1,2),s=field.shape[:3]) return _np.fft.irfftn(grad_,axes=(0,1,2),s=field.shape[:3])
def cell_coord0(grid,size,origin=np.zeros(3)): def cell_coord0(grid,size,origin=_np.zeros(3)):
""" """
Cell center positions (undeformed). Cell center positions (undeformed).
@ -103,7 +103,7 @@ def cell_coord0(grid,size,origin=np.zeros(3)):
""" """
start = origin + size/grid*.5 start = origin + size/grid*.5
end = origin + size - size/grid*.5 end = origin + size - size/grid*.5
return np.mgrid[start[0]:end[0]:grid[0]*1j,start[1]:end[1]:grid[1]*1j,start[2]:end[2]:grid[2]*1j].T return _np.mgrid[start[0]:end[0]:grid[0]*1j,start[1]:end[1]:grid[1]*1j,start[2]:end[2]:grid[2]*1j].T
def cell_displacement_fluct(size,F): def cell_displacement_fluct(size,F):
@ -118,19 +118,19 @@ def cell_displacement_fluct(size,F):
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
displacement = -np.einsum('ijkml,ijkl,l->ijkm', displacement = -_np.einsum('ijkml,ijkl,l->ijkm',
np.fft.rfftn(F,axes=(0,1,2)), _np.fft.rfftn(F,axes=(0,1,2)),
k_s, k_s,
integrator, integrator,
) / k_s_squared[...,np.newaxis] ) / k_s_squared[...,_np.newaxis]
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):
@ -145,8 +145,8 @@ def cell_displacement_avg(size,F):
deformation gradient field. deformation gradient field.
""" """
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):
@ -164,7 +164,7 @@ def cell_displacement(size,F):
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.
@ -193,17 +193,17 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True):
expect coord0 data to be ordered (x fast, z slow). expect coord0 data to be ordered (x fast, z slow).
""" """
coords = [np.unique(coord0[:,i]) for i in range(3)] coords = [_np.unique(coord0[:,i]) for i in range(3)]
mincorner = np.array(list(map(min,coords))) mincorner = _np.array(list(map(min,coords)))
maxcorner = np.array(list(map(max,coords))) maxcorner = _np.array(list(map(max,coords)))
grid = np.array(list(map(len,coords)),'i') grid = _np.array(list(map(len,coords)),'i')
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))
@ -211,12 +211,12 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True):
start = origin + delta*.5 start = origin + delta*.5
end = origin - delta*.5 + size end = origin - delta*.5 + size
if not np.allclose(coords[0],np.linspace(start[0],end[0],grid[0])) and \ if not _np.allclose(coords[0],_np.linspace(start[0],end[0],grid[0])) and \
np.allclose(coords[1],np.linspace(start[1],end[1],grid[1])) and \ _np.allclose(coords[1],_np.linspace(start[1],end[1],grid[1])) and \
np.allclose(coords[2],np.linspace(start[2],end[2],grid[2])): _np.allclose(coords[2],_np.linspace(start[2],end[2],grid[2])):
raise ValueError('Regular grid spacing violated.') raise ValueError('Regular grid spacing violated.')
if ordered and not np.allclose(coord0.reshape(tuple(grid[::-1])+(3,)),cell_coord0(grid,size,origin)): if ordered and not _np.allclose(coord0.reshape(tuple(grid[::-1])+(3,)),cell_coord0(grid,size,origin)):
raise ValueError('Input data is not a regular grid.') raise ValueError('Input data is not a regular grid.')
return (grid,size,origin) return (grid,size,origin)
@ -235,7 +235,7 @@ 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).
@ -249,9 +249,9 @@ def node_coord0(grid,size,origin=np.zeros(3)):
physical origin of the periodic field. Defaults to [0.0,0.0,0.0]. physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
""" """
return np.mgrid[origin[0]:size[0]+origin[0]:(grid[0]+1)*1j, return _np.mgrid[origin[0]:size[0]+origin[0]:(grid[0]+1)*1j,
origin[1]:size[1]+origin[1]:(grid[1]+1)*1j, origin[1]:size[1]+origin[1]:(grid[1]+1)*1j,
origin[2]:size[2]+origin[2]:(grid[2]+1)*1j].T origin[2]:size[2]+origin[2]:(grid[2]+1)*1j].T
def node_displacement_fluct(size,F): def node_displacement_fluct(size,F):
@ -281,8 +281,8 @@ def node_displacement_avg(size,F):
deformation gradient field. deformation gradient field.
""" """
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):
@ -300,7 +300,7 @@ def node_displacement(size,F):
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.
@ -319,18 +319,18 @@ def node_coord(size,F,origin=np.zeros(3)):
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]
@ -347,22 +347,22 @@ def node_coord0_gridSizeOrigin(coord0,ordered=False):
expect coord0 data to be ordered (x fast, z slow). expect coord0 data to be ordered (x fast, z slow).
""" """
coords = [np.unique(coord0[:,i]) for i in range(3)] coords = [_np.unique(coord0[:,i]) for i in range(3)]
mincorner = np.array(list(map(min,coords))) mincorner = _np.array(list(map(min,coords)))
maxcorner = np.array(list(map(max,coords))) maxcorner = _np.array(list(map(max,coords)))
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))
if not np.allclose(coords[0],np.linspace(mincorner[0],maxcorner[0],grid[0]+1)) and \ if not _np.allclose(coords[0],_np.linspace(mincorner[0],maxcorner[0],grid[0]+1)) and \
np.allclose(coords[1],np.linspace(mincorner[1],maxcorner[1],grid[1]+1)) and \ _np.allclose(coords[1],_np.linspace(mincorner[1],maxcorner[1],grid[1]+1)) and \
np.allclose(coords[2],np.linspace(mincorner[2],maxcorner[2],grid[2]+1)): _np.allclose(coords[2],_np.linspace(mincorner[2],maxcorner[2],grid[2]+1)):
raise ValueError('Regular grid spacing violated.') raise ValueError('Regular grid spacing violated.')
if ordered and not np.allclose(coord0.reshape(tuple((grid+1)[::-1])+(3,)),node_coord0(grid,size,origin)): if ordered and not _np.allclose(coord0.reshape(tuple((grid+1)[::-1])+(3,)),node_coord0(grid,size,origin)):
raise ValueError('Input data is not a regular grid.') raise ValueError('Input data is not a regular grid.')
return (grid,size,origin) return (grid,size,origin)
@ -386,10 +386,10 @@ def regrid(size,F,new_grid):
+ cell_displacement_avg(size,F) \ + cell_displacement_avg(size,F) \
+ cell_displacement_fluct(size,F) + cell_displacement_fluct(size,F)
outer = np.dot(np.average(F,axis=(0,1,2)),size) outer = _np.dot(_np.average(F,axis=(0,1,2)),size)
for d in range(3): for d in range(3):
c[np.where(c[:,:,:,d]<0)] += outer[d] c[_np.where(c[:,:,:,d]<0)] += outer[d]
c[np.where(c[:,:,:,d]>outer[d])] -= outer[d] c[_np.where(c[:,:,:,d]>outer[d])] -= outer[d]
tree = spatial.cKDTree(c.reshape(-1,3),boxsize=outer) tree = _spatial.cKDTree(c.reshape(-1,3),boxsize=outer)
return tree.query(cell_coord0(new_grid,outer))[1].flatten() return tree.query(cell_coord0(new_grid,outer))[1].flatten()

76
python/damask/mechanics.py
View File

@ -1,4 +1,4 @@
import numpy as np import numpy as _np
def Cauchy(P,F): def Cauchy(P,F):
""" """
@ -14,10 +14,10 @@ def Cauchy(P,F):
First Piola-Kirchhoff stress. First Piola-Kirchhoff stress.
""" """
if np.shape(F) == np.shape(P) == (3,3): if _np.shape(F) == _np.shape(P) == (3,3):
sigma = 1.0/np.linalg.det(F) * np.dot(P,F.T) sigma = 1.0/_np.linalg.det(F) * _np.dot(P,F.T)
else: else:
sigma = np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F),P,F) sigma = _np.einsum('i,ijk,ilk->ijl',1.0/_np.linalg.det(F),P,F)
return symmetric(sigma) return symmetric(sigma)
@ -31,8 +31,8 @@ def deviatoric_part(T):
Tensor of which the deviatoric part is computed. Tensor of which the deviatoric part is computed.
""" """
return T - np.eye(3)*spherical_part(T) if np.shape(T) == (3,3) else \ return T - _np.eye(3)*spherical_part(T) if _np.shape(T) == (3,3) else \
T - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[T.shape[0],3,3]),spherical_part(T)) T - _np.einsum('ijk,i->ijk',_np.broadcast_to(_np.eye(3),[T.shape[0],3,3]),spherical_part(T))
def eigenvalues(T_sym): def eigenvalues(T_sym):
@ -48,7 +48,7 @@ def eigenvalues(T_sym):
Symmetric tensor of which the eigenvalues are computed. Symmetric tensor of which the eigenvalues are computed.
""" """
return np.linalg.eigvalsh(symmetric(T_sym)) return _np.linalg.eigvalsh(symmetric(T_sym))
def eigenvectors(T_sym,RHS=False): def eigenvectors(T_sym,RHS=False):
@ -65,13 +65,13 @@ def eigenvectors(T_sym,RHS=False):
Enforce right-handed coordinate system. Default is False. Enforce right-handed coordinate system. Default is False.
""" """
(u,v) = np.linalg.eigh(symmetric(T_sym)) (u,v) = _np.linalg.eigh(symmetric(T_sym))
if RHS: if RHS:
if np.shape(T_sym) == (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
@ -99,7 +99,7 @@ def maximum_shear(T_sym):
""" """
w = eigenvalues(T_sym) w = eigenvalues(T_sym)
return (w[0] - w[2])*0.5 if np.shape(T_sym) == (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
@ -141,10 +141,10 @@ def PK2(P,F):
Deformation gradient. Deformation gradient.
""" """
if np.shape(F) == np.shape(P) == (3,3): if _np.shape(F) == _np.shape(P) == (3,3):
S = np.dot(np.linalg.inv(F),P) S = _np.dot(_np.linalg.inv(F),P)
else: else:
S = np.einsum('...jk,...kl->...jl',np.linalg.inv(F),P) S = _np.einsum('...jk,...kl->...jl',_np.linalg.inv(F),P)
return symmetric(S) return symmetric(S)
@ -187,14 +187,14 @@ def spherical_part(T,tensor=False):
""" """
if T.shape == (3,3): if T.shape == (3,3):
sph = np.trace(T)/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(T,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),(T.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):
@ -216,22 +216,22 @@ def strain_tensor(F,t,m):
""" """
F_ = F.reshape(1,3,3) if F.shape == (3,3) else F F_ = F.reshape(1,3,3) if F.shape == (3,3) else F
if t == 'V': if t == 'V':
B = np.matmul(F_,transpose(F_)) B = _np.matmul(F_,transpose(F_))
w,n = np.linalg.eigh(B) w,n = _np.linalg.eigh(B)
elif t == 'U': elif t == 'U':
C = np.matmul(transpose(F_),F_) C = _np.matmul(transpose(F_),F_)
w,n = np.linalg.eigh(C) w,n = _np.linalg.eigh(C)
if m > 0.0: if m > 0.0:
eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n)) eps = 1.0/(2.0*abs(m)) * (+ _np.matmul(n,_np.einsum('ij,ikj->ijk',w**m,n))
- np.broadcast_to(np.eye(3),[F_.shape[0],3,3])) - _np.broadcast_to(_np.eye(3),[F_.shape[0],3,3]))
elif m < 0.0: elif m < 0.0:
eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n)) eps = 1.0/(2.0*abs(m)) * (- _np.matmul(n,_np.einsum('ij,ikj->ijk',w**m,n))
+ np.broadcast_to(np.eye(3),[F_.shape[0],3,3])) + _np.broadcast_to(_np.eye(3),[F_.shape[0],3,3]))
else: else:
eps = np.matmul(n,np.einsum('ij,ikj->ijk',0.5*np.log(w),n)) eps = _np.matmul(n,_np.einsum('ij,ikj->ijk',0.5*_np.log(w),n))
return eps.reshape(3,3) if np.shape(F) == (3,3) else \ return eps.reshape(3,3) if _np.shape(F) == (3,3) else \
eps eps
@ -258,8 +258,8 @@ def transpose(T):
Tensor of which the transpose is computed. Tensor of which the transpose is computed.
""" """
return T.T if np.shape(T) == (3,3) else \ return T.T if _np.shape(T) == (3,3) else \
np.swapaxes(T,axis2=-2,axis1=-1) _np.swapaxes(T,axis2=-2,axis1=-1)
def _polar_decomposition(T,requested): def _polar_decomposition(T,requested):
@ -275,17 +275,17 @@ def _polar_decomposition(T,requested):
‘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(T) u, s, vh = _np.linalg.svd(T)
R = np.dot(u,vh) if np.shape(T) == (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(T,R.T) if np.shape(T) == (3,3) else np.einsum('ijk,ilk->ijl',T,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,T) if np.shape(T) == (3,3) else np.einsum('ikj,ikl->ijl',R,T)) 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)
@ -303,5 +303,5 @@ def _Mises(T_sym,s):
""" """
d = deviatoric_part(T_sym) d = deviatoric_part(T_sym)
return np.sqrt(s*(np.sum(d**2.0))) if np.shape(T_sym) == (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))

153
python/damask/util.py
View File

@ -9,37 +9,22 @@ from optparse import Option
import numpy as np import numpy as np
class bcolors: # limit visibility
""" __all__=[
ASCII Colors. 'srepr',
'croak',
https://svn.blender.org/svnroot/bf-blender/trunk/blender/build_files/scons/tools/bcolors.py 'report',
https://stackoverflow.com/questions/287871 'emph','deemph','delete','strikeout',
""" 'execute',
'show_progress',
HEADER = '\033[95m' 'scale_to_coprime',
OKBLUE = '\033[94m' 'return_message',
OKGREEN = '\033[92m' 'extendableOption',
WARNING = '\033[93m' ]
FAIL = '\033[91m'
ENDC = '\033[0m'
BOLD = '\033[1m'
DIM = '\033[2m'
UNDERLINE = '\033[4m'
CROSSOUT = '\033[9m'
def disable(self):
self.HEADER = ''
self.OKBLUE = ''
self.OKGREEN = ''
self.WARNING = ''
self.FAIL = ''
self.ENDC = ''
self.BOLD = ''
self.UNDERLINE = ''
self.CROSSOUT = ''
####################################################################################################
# Functions
####################################################################################################
def srepr(arg,glue = '\n'): def srepr(arg,glue = '\n'):
r""" r"""
Join arguments as individual lines. Join arguments as individual lines.
@ -144,6 +129,52 @@ def execute(cmd,
return out,error return out,error
def show_progress(iterable,N_iter=None,prefix='',bar_length=50):
"""
Decorate a loop with a status bar.
Use similar like enumerate.
Parameters
----------
iterable : iterable/function with yield statement
Iterable (or function with yield statement) to be decorated.
N_iter : int
Total # of iterations. Needed if number of iterations can not be obtained as len(iterable).
prefix : str, optional.
Prefix string.
bar_length : int, optional
Character length of bar. Defaults to 50.
"""
status = _ProgressBar(N_iter if N_iter else len(iterable),prefix,bar_length)
for i,item in enumerate(iterable):
yield item
status.update(i)
def scale_to_coprime(v):
"""Scale vector to co-prime (relatively prime) integers."""
MAX_DENOMINATOR = 1000
def get_square_denominator(x):
"""Denominator of the square of a number."""
return fractions.Fraction(x ** 2).limit_denominator(MAX_DENOMINATOR).denominator
def lcm(a, b):
"""Least common multiple."""
return a * b // np.gcd(a, b)
denominators = [int(get_square_denominator(i)) for i in v]
s = reduce(lcm, denominators) ** 0.5
m = (np.array(v)*s).astype(np.int)
return m//reduce(np.gcd,m)
####################################################################################################
# Classes
####################################################################################################
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.
@ -215,47 +246,36 @@ class _ProgressBar:
sys.stderr.write('\n') sys.stderr.write('\n')
sys.stderr.flush() sys.stderr.flush()
def show_progress(iterable,N_iter=None,prefix='',bar_length=50):
class bcolors:
""" """
Decorate a loop with a status bar. ASCII Colors.
Use similar like enumerate.
Parameters
----------
iterable : iterable/function with yield statement
Iterable (or function with yield statement) to be decorated.
N_iter : int
Total # of iterations. Needed if number of iterations can not be obtained as len(iterable).
prefix : str, optional.
Prefix string.
bar_length : int, optional
Character length of bar. Defaults to 50.
https://svn.blender.org/svnroot/bf-blender/trunk/blender/build_files/scons/tools/bcolors.py
https://stackoverflow.com/questions/287871
""" """
status = _ProgressBar(N_iter if N_iter else len(iterable),prefix,bar_length)
for i,item in enumerate(iterable): HEADER = '\033[95m'
yield item OKBLUE = '\033[94m'
status.update(i) OKGREEN = '\033[92m'
WARNING = '\033[93m'
FAIL = '\033[91m'
ENDC = '\033[0m'
BOLD = '\033[1m'
DIM = '\033[2m'
UNDERLINE = '\033[4m'
CROSSOUT = '\033[9m'
def disable(self):
def scale_to_coprime(v): self.HEADER = ''
"""Scale vector to co-prime (relatively prime) integers.""" self.OKBLUE = ''
MAX_DENOMINATOR = 1000 self.OKGREEN = ''
self.WARNING = ''
def get_square_denominator(x): self.FAIL = ''
"""Denominator of the square of a number.""" self.ENDC = ''
return fractions.Fraction(x ** 2).limit_denominator(MAX_DENOMINATOR).denominator self.BOLD = ''
self.UNDERLINE = ''
def lcm(a, b): self.CROSSOUT = ''
"""Least common multiple."""
return a * b // np.gcd(a, b)
denominators = [int(get_square_denominator(i)) for i in v]
s = reduce(lcm, denominators) ** 0.5
m = (np.array(v)*s).astype(np.int)
return m//reduce(np.gcd,m)
class return_message: class return_message:
@ -276,4 +296,3 @@ class return_message:
def __repr__(self): def __repr__(self):
"""Return message suitable for interactive shells.""" """Return message suitable for interactive shells."""
return srepr(self.message) return srepr(self.message)

4
python/tests/test_Result.py
View File

@ -24,6 +24,10 @@ def reference_dir(reference_dir_base):
class TestResult: class TestResult:
def test_self_report(self,default):
print(default)
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
default.set_by_time(0.0,20.0) default.set_by_time(0.0,20.0)

24
python/tests/test_grid_filters.py
View File

@ -24,7 +24,7 @@ class TestGridFilters:
n = grid_filters.node_coord0(grid,size) + size/grid*.5 n = grid_filters.node_coord0(grid,size) + size/grid*.5
assert np.allclose(c,n) assert np.allclose(c,n)
@pytest.mark.parametrize('mode',[('cell'),('node')]) @pytest.mark.parametrize('mode',['cell','node'])
def test_grid_DNA(self,mode): def test_grid_DNA(self,mode):
"""Ensure that xx_coord0_gridSizeOrigin is the inverse of xx_coord0.""" """Ensure that xx_coord0_gridSizeOrigin is the inverse of xx_coord0."""
grid = np.random.randint(8,32,(3)) grid = np.random.randint(8,32,(3))
@ -49,7 +49,7 @@ class TestGridFilters:
assert np.allclose(grid_filters.node_coord(size,F) [1:-1,1:-1,1:-1],grid_filters.cell_2_node( assert np.allclose(grid_filters.node_coord(size,F) [1:-1,1:-1,1:-1],grid_filters.cell_2_node(
grid_filters.cell_coord(size,F))[1:-1,1:-1,1:-1]) grid_filters.cell_coord(size,F))[1:-1,1:-1,1:-1])
@pytest.mark.parametrize('mode',[('cell'),('node')]) @pytest.mark.parametrize('mode',['cell','node'])
def test_coord0_origin(self,mode): def test_coord0_origin(self,mode):
origin= np.random.random(3) origin= np.random.random(3)
size = np.random.random(3) # noqa size = np.random.random(3) # noqa
@ -61,22 +61,24 @@ class TestGridFilters:
elif mode == 'node': elif mode == 'node':
assert np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid[::-1]+1)+(3,))) assert np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid[::-1]+1)+(3,)))
@pytest.mark.parametrize('mode',[('cell'),('node')]) @pytest.mark.parametrize('function',[grid_filters.cell_displacement_avg,
def test_displacement_avg_vanishes(self,mode): grid_filters.node_displacement_avg])
def test_displacement_avg_vanishes(self,function):
"""Ensure that random fluctuations in F do not result in average displacement.""" """Ensure that random fluctuations in F do not result in average displacement."""
size = np.random.random(3) # noqa size = np.random.random(3)
grid = np.random.randint(8,32,(3)) grid = np.random.randint(8,32,(3))
F = np.random.random(tuple(grid)+(3,3)) F = np.random.random(tuple(grid)+(3,3))
F += np.eye(3) - np.average(F,axis=(0,1,2)) F += np.eye(3) - np.average(F,axis=(0,1,2))
assert np.allclose(eval('grid_filters.{}_displacement_avg(size,F)'.format(mode)),0.0) assert np.allclose(function(size,F),0.0)
@pytest.mark.parametrize('mode',[('cell'),('node')]) @pytest.mark.parametrize('function',[grid_filters.cell_displacement_fluct,
def test_displacement_fluct_vanishes(self,mode): grid_filters.node_displacement_fluct])
def test_displacement_fluct_vanishes(self,function):
"""Ensure that constant F does not result in fluctuating displacement.""" """Ensure that constant F does not result in fluctuating displacement."""
size = np.random.random(3) # noqa size = np.random.random(3)
grid = np.random.randint(8,32,(3)) grid = np.random.randint(8,32,(3))
F = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3)) # noqa F = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3))
assert np.allclose(eval('grid_filters.{}_displacement_fluct(size,F)'.format(mode)),0.0) assert np.allclose(function(size,F),0.0)
def test_regrid(self): def test_regrid(self):
size = np.random.random(3) size = np.random.random(3)

59
src/LAPACK_interface.f90 Normal file
View File

@ -0,0 +1,59 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Fortran interfaces for LAPACK routines
!> @details https://www.netlib.org/lapack/
!--------------------------------------------------------------------------------------------------
module LAPACK_interface
interface
subroutine dgeev(jobvl,jobvr,n,a,lda,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info)
use prec
character, intent(in) :: jobvl,jobvr
integer, intent(in) :: n,lda,ldvl,ldvr,lwork
real(pReal), intent(inout), dimension(lda,n) :: a
real(pReal), intent(out), dimension(n) :: wr,wi
real(pReal), intent(out), dimension(ldvl,n) :: vl
real(pReal), intent(out), dimension(ldvr,n) :: vr
real(pReal), intent(out), dimension(max(1,lwork)) :: work
integer, intent(out) :: info
end subroutine dgeev
subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info)
use prec
integer, intent(in) :: n,nrhs,lda,ldb
real(pReal), intent(inout), dimension(lda,n) :: a
integer, intent(out), dimension(n) :: ipiv
real(pReal), intent(out), dimension(ldb,nrhs) :: b
integer, intent(out) :: info
end subroutine dgesv
subroutine dgetrf(m,n,a,lda,ipiv,info)
use prec
integer, intent(in) :: m,n,lda
real(pReal), intent(inout), dimension(lda,n) :: a
integer, intent(out), dimension(min(m,n)) :: ipiv
integer, intent(out) :: info
end subroutine dgetrf
subroutine dgetri(n,a,lda,ipiv,work,lwork,info)
use prec
integer, intent(in) :: n,lda,lwork
real(pReal), intent(inout), dimension(lda,n) :: a
integer, intent(out), dimension(n) :: ipiv
real(pReal), intent(out), dimension(max(1,lwork)) :: work
integer, intent(out) :: info
end subroutine dgetri
subroutine dsyev(jobz,uplo,n,a,lda,w,work,lwork,info)
use prec
character, intent(in) :: jobz,uplo
integer, intent(in) :: n,lda,lwork
real(pReal), intent(inout), dimension(lda,n) :: a
real(pReal), intent(out), dimension(n) :: w
real(pReal), intent(out), dimension(max(1,lwork)) :: work
integer, intent(out) :: info
end subroutine dsyev
end interface
end module LAPACK_interface

1
src/commercialFEM_fileList.f90
View File

@ -9,6 +9,7 @@
#include "list.f90" #include "list.f90"
#include "future.f90" #include "future.f90"
#include "config.f90" #include "config.f90"
#include "LAPACK_interface.f90"
#include "math.f90" #include "math.f90"
#include "quaternions.f90" #include "quaternions.f90"
#include "Lambert.f90" #include "Lambert.f90"

2
src/crystallite.f90
View File

@ -835,8 +835,6 @@ logical function integrateStress(ipc,ip,el,timeFraction)
jacoCounterLp, & jacoCounterLp, &
jacoCounterLi ! counters to check for Jacobian update jacoCounterLi ! counters to check for Jacobian update
logical :: error logical :: error
external :: &
dgesv
integrateStress = .false. integrateStress = .false.

343
src/homogenization.f90
View File

@ -27,33 +27,22 @@ module homogenization
implicit none implicit none
private private
!--------------------------------------------------------------------------------------------------
! General variables for the homogenization at a material point
logical, public :: & logical, public :: &
terminallyIll = .false. !< at least one material point is terminally ill terminallyIll = .false. !< at least one material point is terminally ill
real(pReal), dimension(:,:,:,:), allocatable, public :: &
materialpoint_F0, & !< def grad of IP at start of FE increment
materialpoint_F, & !< def grad of IP to be reached at end of FE increment
materialpoint_P !< first P--K stress of IP
real(pReal), dimension(:,:,:,:,:,:), allocatable, public :: &
materialpoint_dPdF !< tangent of first P--K stress at IP
real(pReal), dimension(:,:,:,:), allocatable :: & !--------------------------------------------------------------------------------------------------
materialpoint_subF0, & !< def grad of IP at beginning of homogenization increment ! General variables for the homogenization at a material point
materialpoint_subF !< def grad of IP to be reached at end of homog inc real(pReal), dimension(:,:,:,:), allocatable, public :: &
real(pReal), dimension(:,:), allocatable :: & materialpoint_F0, & !< def grad of IP at start of FE increment
materialpoint_subFrac, & materialpoint_F !< def grad of IP to be reached at end of FE increment
materialpoint_subStep, & real(pReal), dimension(:,:,:,:), allocatable, public, protected :: &
materialpoint_subdt materialpoint_P !< first P--K stress of IP
logical, dimension(:,:), allocatable :: & real(pReal), dimension(:,:,:,:,:,:), allocatable, public, protected :: &
materialpoint_requested, & materialpoint_dPdF !< tangent of first P--K stress at IP
materialpoint_converged
logical, dimension(:,:,:), allocatable :: &
materialpoint_doneAndHappy
type :: tNumerics type :: tNumerics
integer :: & integer :: &
nMPstate !< materialpoint state loop limit nMPstate !< materialpoint state loop limit
real(pReal) :: & real(pReal) :: &
subStepMinHomog, & !< minimum (relative) size of sub-step allowed during cutback in homogenization subStepMinHomog, & !< minimum (relative) size of sub-step allowed during cutback in homogenization
subStepSizeHomog, & !< size of first substep when cutback in homogenization subStepSizeHomog, & !< size of first substep when cutback in homogenization
@ -161,15 +150,7 @@ subroutine homogenization_init
allocate(materialpoint_dPdF(3,3,3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) allocate(materialpoint_dPdF(3,3,3,3,discretization_nIP,discretization_nElem), source=0.0_pReal)
materialpoint_F0 = spread(spread(math_I3,3,discretization_nIP),4,discretization_nElem) ! initialize to identity materialpoint_F0 = spread(spread(math_I3,3,discretization_nIP),4,discretization_nElem) ! initialize to identity
materialpoint_F = materialpoint_F0 ! initialize to identity materialpoint_F = materialpoint_F0 ! initialize to identity
allocate(materialpoint_subF0(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal)
allocate(materialpoint_subF(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal)
allocate(materialpoint_P(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) allocate(materialpoint_P(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal)
allocate(materialpoint_subFrac(discretization_nIP,discretization_nElem), source=0.0_pReal)
allocate(materialpoint_subStep(discretization_nIP,discretization_nElem), source=0.0_pReal)
allocate(materialpoint_subdt(discretization_nIP,discretization_nElem), source=0.0_pReal)
allocate(materialpoint_requested(discretization_nIP,discretization_nElem), source=.false.)
allocate(materialpoint_converged(discretization_nIP,discretization_nElem), source=.true.)
allocate(materialpoint_doneAndHappy(2,discretization_nIP,discretization_nElem), source=.true.)
write(6,'(/,a)') ' <<<+- homogenization init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- homogenization init -+>>>'; flush(6)
@ -203,6 +184,14 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
e, & !< element number e, & !< element number
mySource, & mySource, &
myNgrains myNgrains
real(pReal), dimension(discretization_nIP,discretization_nElem) :: &
subFrac, &
subStep
logical, dimension(discretization_nIP,discretization_nElem) :: &
requested, &
converged
logical, dimension(2,discretization_nIP,discretization_nElem) :: &
doneAndHappy
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then
@ -216,7 +205,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
#endif #endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize restoration points of ... ! initialize restoration points
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(material_homogenizationAt(e)) myNgrains = homogenization_Ngrains(material_homogenizationAt(e))
do i = FEsolving_execIP(1),FEsolving_execIP(2); do i = FEsolving_execIP(1),FEsolving_execIP(2);
@ -238,74 +227,60 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
enddo enddo
subFrac(i,e) = 0.0_pReal
materialpoint_subF0(1:3,1:3,i,e) = materialpoint_F0(1:3,1:3,i,e) converged(i,e) = .false. ! pretend failed step ...
materialpoint_subFrac(i,e) = 0.0_pReal subStep(i,e) = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation
materialpoint_subStep(i,e) = 1.0_pReal/num%subStepSizeHomog ! <<added to adopt flexibility in cutback size>> requested(i,e) = .true. ! everybody requires calculation
materialpoint_converged(i,e) = .false. ! pretend failed step of twice the required size
materialpoint_requested(i,e) = .true. ! everybody requires calculation
if (homogState(material_homogenizationAt(e))%sizeState > 0) & if (homogState(material_homogenizationAt(e))%sizeState > 0) &
homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
homogState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) ! ...internal homogenization state homogState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e))
if (thermalState(material_homogenizationAt(e))%sizeState > 0) & if (thermalState(material_homogenizationAt(e))%sizeState > 0) &
thermalState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & thermalState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
thermalState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) ! ...internal thermal state thermalState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e))
if (damageState(material_homogenizationAt(e))%sizeState > 0) & if (damageState(material_homogenizationAt(e))%sizeState > 0) &
damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
damageState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) ! ...internal damage state damageState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e))
enddo enddo
enddo enddo
NiterationHomog = 0 NiterationHomog = 0
cutBackLooping: do while (.not. terminallyIll .and. & cutBackLooping: do while (.not. terminallyIll .and. &
any(materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > num%subStepMinHomog)) any(subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > num%subStepMinHomog))
!$OMP PARALLEL DO PRIVATE(myNgrains) !$OMP PARALLEL DO PRIVATE(myNgrains)
elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(material_homogenizationAt(e)) myNgrains = homogenization_Ngrains(material_homogenizationAt(e))
IpLooping1: do i = FEsolving_execIP(1),FEsolving_execIP(2) IpLooping1: do i = FEsolving_execIP(1),FEsolving_execIP(2)
converged: if (materialpoint_converged(i,e)) then if (converged(i,e)) then
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0 & if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0 &
.and. ((e == debug_e .and. i == debug_i) & .and. ((e == debug_e .and. i == debug_i) &
.or. .not. iand(debug_level(debug_homogenization),debug_levelSelective) /= 0)) then .or. .not. iand(debug_level(debug_homogenization),debug_levelSelective) /= 0)) then
write(6,'(a,1x,f12.8,1x,a,1x,f12.8,1x,a,i8,1x,i2/)') '<< HOMOG >> winding forward from', & write(6,'(a,1x,f12.8,1x,a,1x,f12.8,1x,a,i8,1x,i2/)') '<< HOMOG >> winding forward from', &
materialpoint_subFrac(i,e), 'to current materialpoint_subFrac', & subFrac(i,e), 'to current subFrac', &
materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),'in materialpoint_stressAndItsTangent at el ip',e,i subFrac(i,e)+subStep(i,e),'in materialpoint_stressAndItsTangent at el ip',e,i
endif endif
#endif #endif
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
! calculate new subStep and new subFrac ! calculate new subStep and new subFrac
materialpoint_subFrac(i,e) = materialpoint_subFrac(i,e) + materialpoint_subStep(i,e) subFrac(i,e) = subFrac(i,e) + subStep(i,e)
materialpoint_subStep(i,e) = min(1.0_pReal-materialpoint_subFrac(i,e), & subStep(i,e) = min(1.0_pReal-subFrac(i,e),num%stepIncreaseHomog*subStep(i,e)) ! introduce flexibility for step increase/acceleration
num%stepIncreaseHomog*materialpoint_subStep(i,e)) ! introduce flexibility for step increase/acceleration
steppingNeeded: if (materialpoint_subStep(i,e) > num%subStepMinHomog) then steppingNeeded: if (subStep(i,e) > num%subStepMinHomog) then
! wind forward grain starting point of... ! wind forward grain starting point
crystallite_partionedF0 (1:3,1:3,1:myNgrains,i,e) = & crystallite_partionedF0 (1:3,1:3,1:myNgrains,i,e) = crystallite_partionedF(1:3,1:3,1:myNgrains,i,e)
crystallite_partionedF(1:3,1:3,1:myNgrains,i,e) crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Fp (1:3,1:3,1:myNgrains,i,e)
crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Lp (1:3,1:3,1:myNgrains,i,e)
crystallite_partionedFp0 (1:3,1:3,1:myNgrains,i,e) = & crystallite_partionedFi0(1:3,1:3,1:myNgrains,i,e) = crystallite_Fi (1:3,1:3,1:myNgrains,i,e)
crystallite_Fp (1:3,1:3,1:myNgrains,i,e) crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e) = crystallite_Li (1:3,1:3,1:myNgrains,i,e)
crystallite_partionedS0 (1:3,1:3,1:myNgrains,i,e) = crystallite_S (1:3,1:3,1:myNgrains,i,e)
crystallite_partionedLp0 (1:3,1:3,1:myNgrains,i,e) = &
crystallite_Lp (1:3,1:3,1:myNgrains,i,e)
crystallite_partionedFi0 (1:3,1:3,1:myNgrains,i,e) = &
crystallite_Fi (1:3,1:3,1:myNgrains,i,e)
crystallite_partionedLi0 (1:3,1:3,1:myNgrains,i,e) = &
crystallite_Li (1:3,1:3,1:myNgrains,i,e)
crystallite_partionedS0 (1:3,1:3,1:myNgrains,i,e) = &
crystallite_S (1:3,1:3,1:myNgrains,i,e)
do g = 1,myNgrains do g = 1,myNgrains
plasticState (material_phaseAt(g,e))%partionedState0(:,material_phasememberAt(g,i,e)) = & plasticState (material_phaseAt(g,e))%partionedState0(:,material_phasememberAt(g,i,e)) = &
@ -326,15 +301,12 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
damageState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e)) damageState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e))
materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e)
endif steppingNeeded endif steppingNeeded
else converged else
if ( (myNgrains == 1 .and. materialpoint_subStep(i,e) <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite if ( (myNgrains == 1 .and. subStep(i,e) <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite
num%subStepSizeHomog * materialpoint_subStep(i,e) <= num%subStepMinHomog ) then ! would require too small subStep num%subStepSizeHomog * subStep(i,e) <= num%subStepMinHomog ) then ! would require too small subStep
! cutback makes no sense ! cutback makes no sense
!$OMP FLUSH(terminallyIll)
if (.not. terminallyIll) then ! so first signals terminally ill... if (.not. terminallyIll) then ! so first signals terminally ill...
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) 'Integration point ', i,' at element ', e, ' terminally ill' write(6,*) 'Integration point ', i,' at element ', e, ' terminally ill'
@ -342,32 +314,27 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
endif endif
terminallyIll = .true. ! ...and kills all others terminallyIll = .true. ! ...and kills all others
else ! cutback makes sense else ! cutback makes sense
materialpoint_subStep(i,e) = num%subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback subStep(i,e) = num%subStepSizeHomog * subStep(i,e) ! crystallite had severe trouble, so do a significant cutback
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0 & if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0 &
.and. ((e == debug_e .and. i == debug_i) & .and. ((e == debug_e .and. i == debug_i) &
.or. .not. iand(debug_level(debug_homogenization), debug_levelSelective) /= 0)) then .or. .not. iand(debug_level(debug_homogenization), debug_levelSelective) /= 0)) then
write(6,'(a,1x,f12.8,a,i8,1x,i2/)') & write(6,'(a,1x,f12.8,a,i8,1x,i2/)') &
'<< HOMOG >> cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep:',& '<< HOMOG >> cutback step in materialpoint_stressAndItsTangent with new subStep:',&
materialpoint_subStep(i,e),' at el ip',e,i subStep(i,e),' at el ip',e,i
endif endif
#endif #endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! restore... ! restore
if (materialpoint_subStep(i,e) < 1.0_pReal) then ! protect against fake cutback from \Delta t = 2 to 1. Maybe that "trick" is not necessary anymore at all? I.e. start with \Delta t = 1 if (subStep(i,e) < 1.0_pReal) then ! protect against fake cutback from \Delta t = 2 to 1. Maybe that "trick" is not necessary anymore at all? I.e. start with \Delta t = 1
crystallite_Lp(1:3,1:3,1:myNgrains,i,e) = & crystallite_Lp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e)
crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) crystallite_Li(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e)
crystallite_Li(1:3,1:3,1:myNgrains,i,e) = &
crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e)
endif ! maybe protecting everything from overwriting (not only L) makes even more sense endif ! maybe protecting everything from overwriting (not only L) makes even more sense
crystallite_Fp(1:3,1:3,1:myNgrains,i,e) = & crystallite_Fp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e)
crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) crystallite_Fi(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedFi0(1:3,1:3,1:myNgrains,i,e)
crystallite_Fi(1:3,1:3,1:myNgrains,i,e) = & crystallite_S (1:3,1:3,1:myNgrains,i,e) = crystallite_partionedS0 (1:3,1:3,1:myNgrains,i,e)
crystallite_partionedFi0(1:3,1:3,1:myNgrains,i,e)
crystallite_S(1:3,1:3,1:myNgrains,i,e) = &
crystallite_partionedS0(1:3,1:3,1:myNgrains,i,e)
do g = 1, myNgrains do g = 1, myNgrains
plasticState (material_phaseAt(g,e))%state( :,material_phasememberAt(g,i,e)) = & plasticState (material_phaseAt(g,e))%state( :,material_phasememberAt(g,i,e)) = &
plasticState (material_phaseAt(g,e))%partionedState0(:,material_phasememberAt(g,i,e)) plasticState (material_phaseAt(g,e))%partionedState0(:,material_phasememberAt(g,i,e))
@ -386,15 +353,11 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
damageState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = & damageState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = &
damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e))
endif endif
endif converged endif
if (materialpoint_subStep(i,e) > num%subStepMinHomog) then if (subStep(i,e) > num%subStepMinHomog) then
materialpoint_requested(i,e) = .true. requested(i,e) = .true.
materialpoint_subF(1:3,1:3,i,e) = materialpoint_subF0(1:3,1:3,i,e) & doneAndHappy(1:2,i,e) = [.false.,.true.]
+ materialpoint_subStep(i,e) * (materialpoint_F(1:3,1:3,i,e) &
- materialpoint_F0(1:3,1:3,i,e))
materialpoint_subdt(i,e) = materialpoint_subStep(i,e) * dt
materialpoint_doneAndHappy(1:2,i,e) = [.false.,.true.]
endif endif
enddo IpLooping1 enddo IpLooping1
enddo elementLooping1 enddo elementLooping1
@ -403,24 +366,24 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
NiterationMPstate = 0 NiterationMPstate = 0
convergenceLooping: do while (.not. terminallyIll .and. & convergenceLooping: do while (.not. terminallyIll .and. &
any( materialpoint_requested(:,FEsolving_execELem(1):FEsolving_execElem(2)) & any( requested(:,FEsolving_execELem(1):FEsolving_execElem(2)) &
.and. .not. materialpoint_doneAndHappy(1,:,FEsolving_execELem(1):FEsolving_execElem(2)) & .and. .not. doneAndHappy(1,:,FEsolving_execELem(1):FEsolving_execElem(2)) &
) .and. & ) .and. &
NiterationMPstate < num%nMPstate) NiterationMPstate < num%nMPstate)
NiterationMPstate = NiterationMPstate + 1 NiterationMPstate = NiterationMPstate + 1
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! deformation partitioning ! deformation partitioning
! based on materialpoint_subF0,.._subF,crystallite_partionedF0, and homogenization_state,
! results in crystallite_partionedF
!$OMP PARALLEL DO PRIVATE(myNgrains) !$OMP PARALLEL DO PRIVATE(myNgrains)
elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(material_homogenizationAt(e)) myNgrains = homogenization_Ngrains(material_homogenizationAt(e))
IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2) IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2)
if ( materialpoint_requested(i,e) .and. & ! process requested but... if(requested(i,e) .and. .not. doneAndHappy(1,i,e)) then ! requested but not yet done
.not. materialpoint_doneAndHappy(1,i,e)) then ! ...not yet done material points call partitionDeformation(materialpoint_F0(1:3,1:3,i,e) &
call partitionDeformation(i,e) ! partition deformation onto constituents + (materialpoint_F(1:3,1:3,i,e)-materialpoint_F0(1:3,1:3,i,e))&
crystallite_dt(1:myNgrains,i,e) = materialpoint_subdt(i,e) ! propagate materialpoint dt to grains *(subStep(i,e)+subFrac(i,e)), &
i,e)
crystallite_dt(1:myNgrains,i,e) = dt*subStep(i,e) ! propagate materialpoint dt to grains
crystallite_requested(1:myNgrains,i,e) = .true. ! request calculation for constituents crystallite_requested(1:myNgrains,i,e) = .true. ! request calculation for constituents
else else
crystallite_requested(1:myNgrains,i,e) = .false. ! calculation for constituents not required anymore crystallite_requested(1:myNgrains,i,e) = .false. ! calculation for constituents not required anymore
@ -431,23 +394,23 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! crystallite integration ! crystallite integration
! based on crystallite_partionedF0,.._partionedF converged = crystallite_stress() !ToDo: MD not sure if that is the best logic
! incrementing by crystallite_dt
materialpoint_converged = crystallite_stress() !ToDo: MD not sure if that is the best logic
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state update ! state update
!$OMP PARALLEL DO !$OMP PARALLEL DO
elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2)
IpLooping3: do i = FEsolving_execIP(1),FEsolving_execIP(2) IpLooping3: do i = FEsolving_execIP(1),FEsolving_execIP(2)
if ( materialpoint_requested(i,e) .and. & if (requested(i,e) .and. .not. doneAndHappy(1,i,e)) then
.not. materialpoint_doneAndHappy(1,i,e)) then if (.not. converged(i,e)) then
if (.not. materialpoint_converged(i,e)) then doneAndHappy(1:2,i,e) = [.true.,.false.]
materialpoint_doneAndHappy(1:2,i,e) = [.true.,.false.]
else else
materialpoint_doneAndHappy(1:2,i,e) = updateState(i,e) doneAndHappy(1:2,i,e) = updateState(dt*subStep(i,e), &
materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(1:2,i,e)) ! converged if done and happy materialpoint_F0(1:3,1:3,i,e) &
+ (materialpoint_F(1:3,1:3,i,e)-materialpoint_F0(1:3,1:3,i,e)) &
*(subStep(i,e)+subFrac(i,e)), &
i,e)
converged(i,e) = all(doneAndHappy(1:2,i,e)) ! converged if done and happy
endif endif
endif endif
enddo IpLooping3 enddo IpLooping3
@ -481,29 +444,31 @@ end subroutine materialpoint_stressAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief partition material point def grad onto constituents !> @brief partition material point def grad onto constituents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine partitionDeformation(ip,el) subroutine partitionDeformation(subF,ip,el)
integer, intent(in) :: & real(pReal), intent(in), dimension(3,3) :: &
ip, & !< integration point subF
el !< element number integer, intent(in) :: &
ip, & !< integration point
el !< element number
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
crystallite_partionedF(1:3,1:3,1,ip,el) = materialpoint_subF(1:3,1:3,ip,el) crystallite_partionedF(1:3,1:3,1,ip,el) = subF
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call mech_isostrain_partitionDeformation(& call mech_isostrain_partitionDeformation(&
crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
subF)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call mech_RGC_partitionDeformation(&
crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
materialpoint_subF(1:3,1:3,ip,el)) subF,&
ip, &
case (HOMOGENIZATION_RGC_ID) chosenHomogenization el)
call mech_RGC_partitionDeformation(& end select chosenHomogenization
crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
materialpoint_subF(1:3,1:3,ip,el),&
ip, &
el)
end select chosenHomogenization
end subroutine partitionDeformation end subroutine partitionDeformation
@ -512,45 +477,49 @@ end subroutine partitionDeformation
!> @brief update the internal state of the homogenization scheme and tell whether "done" and !> @brief update the internal state of the homogenization scheme and tell whether "done" and
!> "happy" with result !> "happy" with result
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function updateState(ip,el) function updateState(subdt,subF,ip,el)
integer, intent(in) :: & real(pReal), intent(in) :: &
ip, & !< integration point subdt !< current time step
el !< element number real(pReal), intent(in), dimension(3,3) :: &
logical, dimension(2) :: updateState subF
integer, intent(in) :: &
ip, & !< integration point
el !< element number
logical, dimension(2) :: updateState
updateState = .true. updateState = .true.
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
updateState = & updateState = &
updateState .and. & updateState .and. &
mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el),& crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el),&
materialpoint_subF(1:3,1:3,ip,el),& subF,&
materialpoint_subdt(ip,el), & subdt, &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
ip, & ip, &
el) el)
end select chosenHomogenization end select chosenHomogenization
chosenThermal: select case (thermal_type(material_homogenizationAt(el))) chosenThermal: select case (thermal_type(material_homogenizationAt(el)))
case (THERMAL_adiabatic_ID) chosenThermal case (THERMAL_adiabatic_ID) chosenThermal
updateState = & updateState = &
updateState .and. & updateState .and. &
thermal_adiabatic_updateState(materialpoint_subdt(ip,el), & thermal_adiabatic_updateState(subdt, &
ip, & ip, &
el) el)
end select chosenThermal end select chosenThermal
chosenDamage: select case (damage_type(material_homogenizationAt(el))) chosenDamage: select case (damage_type(material_homogenizationAt(el)))
case (DAMAGE_local_ID) chosenDamage case (DAMAGE_local_ID) chosenDamage
updateState = & updateState = &
updateState .and. & updateState .and. &
damage_local_updateState(materialpoint_subdt(ip,el), & damage_local_updateState(subdt, &
ip, & ip, &
el) el)
end select chosenDamage end select chosenDamage
end function updateState end function updateState
@ -560,31 +529,31 @@ end function updateState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine averageStressAndItsTangent(ip,el) subroutine averageStressAndItsTangent(ip,el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point
el !< element number el !< element number
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
materialpoint_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el) materialpoint_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el)
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_dPdF(1:3,1:3,1:3,1:3,1,ip,el) materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_dPdF(1:3,1:3,1:3,1:3,1,ip,el)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call mech_isostrain_averageStressAndItsTangent(& call mech_isostrain_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), & materialpoint_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
homogenization_typeInstance(material_homogenizationAt(el))) homogenization_typeInstance(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call mech_RGC_averageStressAndItsTangent(& call mech_RGC_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), & materialpoint_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
homogenization_typeInstance(material_homogenizationAt(el))) homogenization_typeInstance(material_homogenizationAt(el)))
end select chosenHomogenization end select chosenHomogenization
end subroutine averageStressAndItsTangent end subroutine averageStressAndItsTangent

30
src/math.f90
View File

@ -9,6 +9,7 @@ module math
use prec use prec
use IO use IO
use numerics use numerics
use LAPACK_interface
implicit none implicit none
public public
@ -485,18 +486,14 @@ function math_invSym3333(A)
real(pReal),dimension(3,3,3,3),intent(in) :: A real(pReal),dimension(3,3,3,3),intent(in) :: A
integer, dimension(6) :: ipiv6 integer, dimension(6) :: ipiv6
real(pReal), dimension(6,6) :: temp66 real(pReal), dimension(6,6) :: temp66
real(pReal), dimension(6*(64+2)) :: work real(pReal), dimension(6*6) :: work
integer :: ierr_i, ierr_f integer :: ierr_i, ierr_f
external :: &
dgetrf, &
dgetri
temp66 = math_sym3333to66(A) temp66 = math_sym3333to66(A)
call dgetrf(6,6,temp66,6,ipiv6,ierr_i) call dgetrf(6,6,temp66,6,ipiv6,ierr_i)
call dgetri(6,temp66,6,ipiv6,work,size(work,1),ierr_f) call dgetri(6,temp66,6,ipiv6,work,size(work,1),ierr_f)
if (ierr_i /= 0 .or. ierr_f /= 0) then if (ierr_i /= 0 .or. ierr_f /= 0) then
call IO_error(400, ext_msg = 'math_invSym3333') call IO_error(400, ext_msg = 'math_invSym3333')
else else
@ -515,12 +512,9 @@ subroutine math_invert(InvA, error, A)
real(pReal), dimension(size(A,1),size(A,1)), intent(out) :: invA real(pReal), dimension(size(A,1),size(A,1)), intent(out) :: invA
logical, intent(out) :: error logical, intent(out) :: error
integer, dimension(size(A,1)) :: ipiv integer, dimension(size(A,1)) :: ipiv
real(pReal), dimension(size(A,1)*(64+2)) :: work real(pReal), dimension(size(A,1)**2) :: work
integer :: ierr integer :: ierr
external :: &
dgetrf, &
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)
@ -884,9 +878,7 @@ subroutine math_eigh(m,w,v,error)
logical, intent(out) :: error logical, intent(out) :: error
integer :: ierr integer :: ierr
real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f real(pReal), dimension(size(m,1)**2) :: work
external :: &
dsyev
v = m ! copy matrix to input (doubles as output) array v = m ! copy matrix to input (doubles as output) array
call dsyev('V','U',size(m,1),v,size(m,1),w,work,size(work,1),ierr) call dsyev('V','U',size(m,1),v,size(m,1),w,work,size(work,1),ierr)
@ -1041,9 +1033,7 @@ function math_eigvalsh(m)
real(pReal), dimension(size(m,1),size(m,1)) :: m_ real(pReal), dimension(size(m,1),size(m,1)) :: m_
integer :: ierr integer :: ierr
real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f real(pReal), dimension(size(m,1)**2) :: work
external :: &
dsyev
m_= m ! copy matrix to input (will be destroyed) m_= m ! copy matrix to input (will be destroyed)
call dsyev('N','U',size(m,1),m_,size(m,1),math_eigvalsh,work,size(work,1),ierr) call dsyev('N','U',size(m,1),m_,size(m,1),math_eigvalsh,work,size(work,1),ierr)

76
src/quaternions.f90
View File

@ -63,27 +63,27 @@ module quaternions
module procedure assign_quat__ module procedure assign_quat__
module procedure assign_vec__ module procedure assign_vec__
end interface assignment (=) end interface assignment (=)
interface quaternion interface quaternion
module procedure init__ module procedure init__
end interface quaternion end interface quaternion
interface abs interface abs
procedure abs__ procedure abs__
end interface abs end interface abs
interface dot_product interface dot_product
procedure dot_product__ procedure dot_product__
end interface dot_product end interface dot_product
interface conjg interface conjg
module procedure conjg__ module procedure conjg__
end interface conjg end interface conjg
interface exp interface exp
module procedure exp__ module procedure exp__
end interface exp end interface exp
interface log interface log
module procedure log__ module procedure log__
end interface log end interface log
@ -95,7 +95,7 @@ module quaternions
interface aimag interface aimag
module procedure aimag__ module procedure aimag__
end interface aimag end interface aimag
public :: & public :: &
quaternions_init, & quaternions_init, &
assignment(=), & assignment(=), &
@ -118,7 +118,7 @@ end subroutine quaternions_init
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> construct a quaternion from a 4-vector !> @brief construct a quaternion from a 4-vector
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) pure function init__(array) type(quaternion) pure function init__(array)
@ -133,7 +133,7 @@ end function init__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> assign a quaternion !> @brief assign a quaternion
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
elemental pure subroutine assign_quat__(self,other) elemental pure subroutine assign_quat__(self,other)
@ -141,12 +141,12 @@ elemental pure subroutine assign_quat__(self,other)
type(quaternion), intent(in) :: other type(quaternion), intent(in) :: other
self = [other%w,other%x,other%y,other%z] self = [other%w,other%x,other%y,other%z]
end subroutine assign_quat__ end subroutine assign_quat__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> assign a 4-vector !> @brief assign a 4-vector
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
pure subroutine assign_vec__(self,other) pure subroutine assign_vec__(self,other)
@ -162,7 +162,7 @@ end subroutine assign_vec__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> add a quaternion !> @brief add a quaternion
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function add__(self,other) type(quaternion) elemental pure function add__(self,other)
@ -170,24 +170,24 @@ type(quaternion) elemental pure function add__(self,other)
add__ = [ self%w, self%x, self%y ,self%z] & add__ = [ self%w, self%x, self%y ,self%z] &
+ [other%w, other%x, other%y,other%z] + [other%w, other%x, other%y,other%z]
end function add__ end function add__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> return (unary positive operator) !> @brief return (unary positive operator)
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function pos__(self) type(quaternion) elemental pure function pos__(self)
class(quaternion), intent(in) :: self class(quaternion), intent(in) :: self
pos__ = self * (+1.0_pReal) pos__ = self * (+1.0_pReal)
end function pos__ end function pos__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> subtract a quaternion !> @brief subtract a quaternion
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function sub__(self,other) type(quaternion) elemental pure function sub__(self,other)
@ -195,24 +195,24 @@ type(quaternion) elemental pure function sub__(self,other)
sub__ = [ self%w, self%x, self%y ,self%z] & sub__ = [ self%w, self%x, self%y ,self%z] &
- [other%w, other%x, other%y,other%z] - [other%w, other%x, other%y,other%z]
end function sub__ end function sub__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> negate (unary negative operator) !> @brief negate (unary negative operator)
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function neg__(self) type(quaternion) elemental pure function neg__(self)
class(quaternion), intent(in) :: self class(quaternion), intent(in) :: self
neg__ = self * (-1.0_pReal) neg__ = self * (-1.0_pReal)
end function neg__ end function neg__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> multiply with a quaternion !> @brief multiply with a quaternion
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function mul_quat__(self,other) type(quaternion) elemental pure function mul_quat__(self,other)
@ -227,7 +227,7 @@ end function mul_quat__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> multiply with a scalar !> @brief multiply with a scalar
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function mul_scal__(self,scal) type(quaternion) elemental pure function mul_scal__(self,scal)
@ -235,12 +235,12 @@ type(quaternion) elemental pure function mul_scal__(self,scal)
real(pReal), intent(in) :: scal real(pReal), intent(in) :: scal
mul_scal__ = [self%w,self%x,self%y,self%z]*scal mul_scal__ = [self%w,self%x,self%y,self%z]*scal
end function mul_scal__ end function mul_scal__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> divide by a quaternion !> @brief divide by a quaternion
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function div_quat__(self,other) type(quaternion) elemental pure function div_quat__(self,other)
@ -252,7 +252,7 @@ end function div_quat__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> divide by a scalar !> @brief divide by a scalar
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function div_scal__(self,scal) type(quaternion) elemental pure function div_scal__(self,scal)
@ -265,7 +265,7 @@ end function div_scal__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> test equality !> @brief test equality
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
logical elemental pure function eq__(self,other) logical elemental pure function eq__(self,other)
@ -278,7 +278,7 @@ end function eq__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> test inequality !> @brief test inequality
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
logical elemental pure function neq__(self,other) logical elemental pure function neq__(self,other)
@ -290,7 +290,7 @@ end function neq__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> raise to the power of a quaternion !> @brief raise to the power of a quaternion
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function pow_quat__(self,expon) type(quaternion) elemental pure function pow_quat__(self,expon)
@ -303,7 +303,7 @@ end function pow_quat__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> raise to the power of a scalar !> @brief raise to the power of a scalar
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function pow_scal__(self,expon) type(quaternion) elemental pure function pow_scal__(self,expon)
@ -316,7 +316,7 @@ end function pow_scal__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> take exponential !> @brief take exponential
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function exp__(a) type(quaternion) elemental pure function exp__(a)
@ -336,7 +336,7 @@ end function exp__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> take logarithm !> @brief take logarithm
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function log__(a) type(quaternion) elemental pure function log__(a)
@ -356,7 +356,7 @@ end function log__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> return norm !> @brief return norm
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
real(pReal) elemental pure function abs__(self) real(pReal) elemental pure function abs__(self)
@ -368,7 +368,7 @@ end function abs__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> calculate dot product !> @brief calculate dot product
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
real(pReal) elemental pure function dot_product__(a,b) real(pReal) elemental pure function dot_product__(a,b)
@ -380,7 +380,7 @@ end function dot_product__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> take conjugate complex !> @brief take conjugate complex
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function conjg__(self) type(quaternion) elemental pure function conjg__(self)
@ -392,7 +392,7 @@ end function conjg__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> homomorph !> @brief homomorph
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function homomorphed(self) type(quaternion) elemental pure function homomorphed(self)
@ -404,7 +404,7 @@ end function homomorphed
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> return as plain array !> @brief return as plain array
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
pure function asArray(self) pure function asArray(self)
@ -417,7 +417,7 @@ end function asArray
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> real part (scalar) !> @brief real part (scalar)
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
pure function real__(self) pure function real__(self)
@ -430,7 +430,7 @@ end function real__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> imaginary part (3-vector) !> @brief imaginary part (3-vector)
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
pure function aimag__(self) pure function aimag__(self)
@ -443,7 +443,7 @@ end function aimag__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> inverse !> @brief inverse
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
type(quaternion) elemental pure function inverse(self) type(quaternion) elemental pure function inverse(self)

2
src/rotations.f90
View File

@ -595,8 +595,6 @@ function om2ax(om) result(ax)
real(pReal), dimension(3,3) :: VR, devNull, om_ real(pReal), dimension(3,3) :: VR, devNull, om_
integer :: ierr, i integer :: ierr, i
external :: dgeev
om_ = om om_ = om
! first get the rotation angle ! first get the rotation angle