Merge branch 'development' into state-integration-cleaning

This commit is contained in:
Martin Diehl 2020-04-10 13:20:38 +02:00
commit 6082f3c1aa
20 changed files with 314 additions and 119 deletions

View File

@ -148,7 +148,7 @@ else ()
set(OPENMP "${OPENMP}")
endif ()
# syntax check only (mainly for pre-receive hook, works only with gfortran)
# syntax check only (mainly for pre-receive hook)
if (CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
set (BUILDCMD_POST "${BUILDCMD_POST} -fsyntax-only")
endif ()

View File

@ -1 +1 @@
v2.0.3-2194-g369682aa
v2.0.3-2243-gbb03483b

View File

@ -21,36 +21,42 @@ def findClosestSeed(seeds, weights, point):
return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights)
def Laguerre_tessellation(grid, seeds, grains, size, periodic, weights, cpus):
def Laguerre_tessellation(grid, size, seeds, weights, origin = np.zeros(3), periodic = True, cpus = 2):
if periodic:
weights_p = np.tile(weights,27).flatten(order='F') # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
coords = damask.grid_filters.cell_coord0(grid*3,size*3,-origin-size).reshape(-1,3,order='F')
else:
weights_p = weights.flatten()
seeds_p = seeds
coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F')
if cpus > 1:
default_args = partial(findClosestSeed,seeds_p,weights_p)
pool = multiprocessing.Pool(processes = cpus) # initialize workers
result = pool.map_async(default_args, [point for point in grid]) # evaluate function in parallel
pool = multiprocessing.Pool(processes = cpus)
result = pool.map_async(partial(findClosestSeed,seeds_p,weights_p), [coord for coord in coords])
pool.close()
pool.join()
closestSeeds = np.array(result.get()).flatten()
closest_seed = np.array(result.get())
else:
closestSeeds= np.array([findClosestSeed(seeds_p,weights_p,point) for point in grid])
closest_seed= np.array([findClosestSeed(seeds_p,weights_p,coord) for coord in coords])
return grains[closestSeeds%seeds.shape[0]]
if periodic:
closest_seed = closest_seed.reshape(grid*3)
return closest_seed[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0]
else:
return closest_seed
def Voronoi_tessellation(grid, seeds, grains, size, periodic = True):
def Voronoi_tessellation(grid, size, seeds, origin = np.zeros(3), periodic = True):
coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F')
KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
devNull,closestSeeds = KDTree.query(grid)
devNull,closest_seed = KDTree.query(coords)
return grains[closestSeeds]
return closest_seed
# --------------------------------------------------------------------
@ -186,18 +192,15 @@ for name in filenames:
grains = table.get(options.microstructure) if options.microstructure in table.labels else np.arange(len(seeds))+1
grainIDs = np.unique(grains).astype('i')
NgrainIDs = len(grainIDs)
if options.eulers in table.labels:
eulers = table.get(options.eulers)
coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F')
if options.laguerre:
indices = Laguerre_tessellation(coords,seeds,grains,size,options.periodic,
table.get(options.weight),options.cpus)
indices = grains[Laguerre_tessellation(grid,size,seeds,table.get(options.weight),origin,
options.periodic,options.cpus)]
else:
indices = Voronoi_tessellation (coords,seeds,grains,size,options.periodic)
indices = grains[Voronoi_tessellation (grid,size,seeds,origin,options.periodic)]
config_header = []
if options.config:

View File

@ -52,7 +52,7 @@ for name in filenames:
geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
substituted = geom.get_microstructure()
for old,new in zip(sub[0::2],sub[1::2]): substituted[substituted==old] = new # substitute microstructure indices
for old,new in zip(sub[0::2],sub[1::2]): substituted[geom.microstructure==old] = new # substitute microstructure indices
substituted += options.microstructure # constant shift
damask.util.croak(geom.update(substituted,origin=geom.get_origin()+options.origin))

2
python/.coveragerc Normal file
View File

@ -0,0 +1,2 @@
[run]
omit = tests/*

1
python/.gitignore vendored
View File

@ -1,3 +1,4 @@
*.pyc
dist
damask.egg-info
.coverage

View File

@ -1,11 +1,15 @@
import sys
from io import StringIO
import multiprocessing
from functools import partial
import numpy as np
from scipy import ndimage
from scipy import ndimage,spatial
from . import VTK
from . import util
from . import Environment
from . import grid_filters
class Geom:
@ -50,7 +54,7 @@ class Geom:
def update(self,microstructure=None,size=None,origin=None,rescale=False):
"""
Updates microstructure and size.
Update microstructure and size.
Parameters
----------
@ -113,7 +117,7 @@ class Geom:
def set_comments(self,comments):
"""
Replaces all existing comments.
Replace all existing comments.
Parameters
----------
@ -127,7 +131,7 @@ class Geom:
def add_comments(self,comments):
"""
Appends comments to existing comments.
Append comments to existing comments.
Parameters
----------
@ -140,7 +144,7 @@ class Geom:
def set_microstructure(self,microstructure):
"""
Replaces the existing microstructure representation.
Replace the existing microstructure representation.
Parameters
----------
@ -159,7 +163,7 @@ class Geom:
def set_size(self,size):
"""
Replaces the existing size information.
Replace the existing size information.
Parameters
----------
@ -179,7 +183,7 @@ class Geom:
def set_origin(self,origin):
"""
Replaces the existing origin information.
Replace the existing origin information.
Parameters
----------
@ -196,7 +200,7 @@ class Geom:
def set_homogenization(self,homogenization):
"""
Replaces the existing homogenization index.
Replace the existing homogenization index.
Parameters
----------
@ -264,7 +268,7 @@ class Geom:
@staticmethod
def from_file(fname):
"""
Reads a geom file.
Read a geom file.
Parameters
----------
@ -325,6 +329,81 @@ class Geom:
return Geom(microstructure.reshape(grid),size,origin,homogenization,comments)
@staticmethod
def _find_closest_seed(seeds, weights, point):
return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights)
@staticmethod
def from_Laguerre_tessellation(grid,size,seeds,weights,periodic=True):
"""
Generate geometry from Laguerre tessellation.
Parameters
----------
grid : numpy.ndarray of shape (3)
number of grid points in x,y,z direction.
size : list or numpy.ndarray of shape (3)
physical size of the microstructure in meter.
seeds : numpy.ndarray of shape (:,3)
position of the seed points in meter. All points need to lay within the box.
weights : numpy.ndarray of shape (seeds.shape[0])
weights of the seeds. Setting all weights to 1.0 gives a standard Voronoi tessellation.
periodic : Boolean, optional
perform a periodic tessellation. Defaults to True.
"""
if periodic:
weights_p = np.tile(weights,27).flatten(order='F') # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
coords = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3,order='F')
else:
weights_p = weights.flatten()
seeds_p = seeds
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
pool = multiprocessing.Pool(processes = int(Environment().options['DAMASK_NUM_THREADS']))
result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords])
pool.close()
pool.join()
microstructure = np.array(result.get())
if periodic:
microstructure = microstructure.reshape(grid*3)
microstructure = microstructure[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0]
else:
microstructure = microstructure.reshape(grid)
#comments = 'geom.py:from_Laguerre_tessellation v{}'.format(version)
return Geom(microstructure+1,size,homogenization=1)
@staticmethod
def from_Voronoi_tessellation(grid,size,seeds,periodic=True):
"""
Generate geometry from Voronoi tessellation.
Parameters
----------
grid : numpy.ndarray of shape (3)
number of grid points in x,y,z direction.
size : list or numpy.ndarray of shape (3)
physical size of the microstructure in meter.
seeds : numpy.ndarray of shape (:,3)
position of the seed points in meter. All points need to lay within the box.
periodic : Boolean, optional
perform a periodic tessellation. Defaults to True.
"""
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
devNull,microstructure = KDTree.query(coords)
#comments = 'geom.py:from_Voronoi_tessellation v{}'.format(version)
return Geom(microstructure.reshape(grid)+1,size,homogenization=1)
def to_file(self,fname,pack=None):
"""
Writes a geom file.

View File

@ -72,7 +72,7 @@ class VTK:
connectivity : numpy.ndarray of np.dtype = int
Cell connectivity (0-based), first dimension determines #Cells, second dimension determines #Nodes/Cell.
cell_type : str
Name of the vtk.vtkCell subclass. Tested for TRIANGLE, QUAD, and HEXAHEDRON.
Name of the vtk.vtkCell subclass. Tested for TRIANGLE, QUAD, TETRA, and HEXAHEDRON.
"""
vtk_nodes = vtk.vtkPoints()

View File

@ -28,7 +28,7 @@ def reference_dir(reference_dir_base):
class TestGeom:
def test_update(self,default):
modified = copy.deepcopy(default)
modified.update(
@ -72,7 +72,7 @@ class TestGeom:
if update: modified.to_file(reference)
assert geom_equal(modified,Geom.from_file(reference))
@pytest.mark.parametrize('stencil',[(1),(2),(3),(4)])
@pytest.mark.parametrize('stencil',[1,2,3,4])
def test_clean(self,default,update,reference_dir,stencil):
modified = copy.deepcopy(default)
modified.clean(stencil)
@ -82,12 +82,12 @@ class TestGeom:
assert geom_equal(modified,Geom.from_file(reference))
@pytest.mark.parametrize('grid',[
((10,11,10)),
([10,13,10]),
(np.array((10,10,10))),
(np.array((8, 10,12))),
(np.array((5, 4, 20))),
(np.array((10,20,2)) )
(10,11,10),
[10,13,10],
np.array((10,10,10)),
np.array((8, 10,12)),
np.array((5, 4, 20)),
np.array((10,20,2))
]
)
def test_scale(self,default,update,reference_dir,grid):
@ -97,3 +97,37 @@ class TestGeom:
reference = os.path.join(reference_dir,'scale_{}.geom'.format(tag))
if update: modified.to_file(reference)
assert geom_equal(modified,Geom.from_file(reference))
@pytest.mark.parametrize('periodic',[True,False])
def test_tessellation_approaches(self,periodic):
grid = np.random.randint(10,20,3)
size = np.random.random(3) + 1.0
N_seeds= np.random.randint(10,30)
seeds = np.random.rand(N_seeds,3) * np.broadcast_to(size,(N_seeds,3))
Voronoi = Geom.from_Voronoi_tessellation( grid,size,seeds, periodic)
Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(N_seeds),periodic)
assert geom_equal(Laguerre,Voronoi)
def test_Laguerre_weights(self):
grid = np.random.randint(10,20,3)
size = np.random.random(3) + 1.0
N_seeds= np.random.randint(10,30)
seeds = np.random.rand(N_seeds,3) * np.broadcast_to(size,(N_seeds,3))
weights= np.full((N_seeds),-np.inf)
ms = np.random.randint(1, N_seeds+1)
weights[ms-1] = np.random.random()
Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,np.random.random()>0.5)
assert np.all(Laguerre.microstructure == ms)
@pytest.mark.parametrize('approach',['Laguerre','Voronoi'])
def test_tessellate_bicrystal(self,approach):
grid = np.random.randint(5,10,3)*2
size = grid.astype(np.float)
seeds = np.vstack((size*np.array([0.5,0.25,0.5]),size*np.array([0.5,0.75,0.5])))
microstructure = np.ones(grid)
microstructure[:,grid[1]//2:,:] = 2
if approach == 'Laguerre':
geom = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(2),np.random.random()>0.5)
elif approach == 'Voronoi':
geom = Geom.from_Voronoi_tessellation(grid,size,seeds, np.random.random()>0.5)
assert np.all(geom.microstructure == microstructure)

47
python/tests/test_VTK.py Normal file
View File

@ -0,0 +1,47 @@
import os
import pytest
import numpy as np
from damask import VTK
@pytest.fixture
def reference_dir(reference_dir_base):
"""Directory containing reference results."""
return os.path.join(reference_dir_base,'Result')
class TestVTK:
def test_rectilinearGrid(self,tmp_path):
grid = np.random.randint(5,10,3)*2
size = np.random.random(3) + 1.0
origin = np.random.random(3)
v = VTK.from_rectilinearGrid(grid,size,origin)
s = v.__repr__()
v.write(os.path.join(tmp_path,'rectilinearGrid'))
v = VTK.from_file(os.path.join(tmp_path,'rectilinearGrid.vtr'))
assert(s == v.__repr__())
def test_polyData(self,tmp_path):
points = np.random.rand(3,100)
v = VTK.from_polyData(points)
s = v.__repr__()
v.write(os.path.join(tmp_path,'polyData'))
v = VTK.from_file(os.path.join(tmp_path,'polyData.vtp'))
assert(s == v.__repr__())
@pytest.mark.parametrize('cell_type,n',[
('VTK_hexahedron',8),
('TETRA',4),
('quad',4),
('VTK_TRIANGLE',3)
]
)
def test_unstructuredGrid(self,tmp_path,cell_type,n):
nodes = np.random.rand(n,3)
connectivity = np.random.choice(np.arange(n),n,False).reshape(-1,n)
v = VTK.from_unstructuredGrid(nodes,connectivity,cell_type)
s = v.__repr__()
v.write(os.path.join(tmp_path,'unstructuredGrid'))
v = VTK.from_file(os.path.join(tmp_path,'unstructuredGrid.vtu'))
assert(s == v.__repr__())

View File

@ -13,7 +13,7 @@
#define INTEL_MIN 1700
#define PETSC_MAJOR 3
#define PETSC_MINOR_MIN 10
#define PETSC_MINOR_MAX 12
#define PETSC_MINOR_MAX 13
module DAMASK_interface
use, intrinsic :: iso_fortran_env
@ -25,7 +25,7 @@ module DAMASK_interface
implicit none
private
logical, public, protected :: &
logical, volatile, public, protected :: &
SIGTERM, & !< termination signal
SIGUSR1, & !< 1. user-defined signal
SIGUSR2 !< 2. user-defined signal

View File

@ -43,6 +43,9 @@ module DAMASK_interface
logical, protected, public :: symmetricSolver
character(len=*), parameter, public :: INPUTFILEEXTENSION = '.dat'
logical, dimension(:,:), public, allocatable :: &
calcMode !< calculate or collect (ping pong scheme)
public :: &
DAMASK_interface_init, &
getSolverJobName
@ -102,7 +105,6 @@ function getSolverJobName()
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash
integer :: extPos
inputName=''
inquire(5, name=inputName) ! determine inputfile
extPos = len_trim(inputName)-4
getSolverJobName=inputName(scan(inputName,pathSep,back=.true.)+1:extPos)
@ -179,6 +181,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
use FEsolving
use debug
use discretization_marc
use homogenization
use CPFEM
implicit none

View File

@ -4,20 +4,12 @@
!> @brief global variables for flow control
!--------------------------------------------------------------------------------------------------
module FEsolving
implicit none
public
logical :: &
terminallyIll = .false. !< at least one material point is terminally ill
integer, dimension(2) :: &
FEsolving_execElem, & !< for ping-pong scheme always whole range, otherwise one specific element
FEsolving_execIP !< for ping-pong scheme always range to max IP, otherwise one specific IP
#if defined(Marc4DAMASK)
logical, dimension(:,:), allocatable :: &
calcMode !< do calculation or simply collect when using ping pong scheme
#endif
end module FEsolving

View File

@ -568,7 +568,22 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
+ ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
enddo slipContribution
!ToDo: Why do this before shear banding?
call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,dot_gamma_twin,ddot_gamma_dtau_twin)
twinContibution: do i = 1, prm%sum_N_tw
Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
enddo twinContibution
call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,dot_gamma_tr,ddot_gamma_dtau_trans)
transContibution: do i = 1, prm%sum_N_tr
Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau_trans(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i)
enddo transContibution
Lp = Lp * f_unrotated
dLp_dMp = dLp_dMp * f_unrotated
@ -598,23 +613,6 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
endif shearBandingContribution
call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,dot_gamma_twin,ddot_gamma_dtau_twin)
twinContibution: do i = 1, prm%sum_N_tw
Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i) * f_unrotated
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) * f_unrotated
enddo twinContibution
call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,dot_gamma_tr,ddot_gamma_dtau_trans)
transContibution: do i = 1, prm%sum_N_tr
Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) * f_unrotated
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ ddot_gamma_dtau_trans(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) * f_unrotated
enddo transContibution
end associate
end subroutine plastic_dislotwin_LpAndItsTangent

View File

@ -16,15 +16,15 @@ submodule(constitutive) plastic_nonlocal
kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin
! storage order of dislocation types
integer, dimension(8), parameter :: &
integer, dimension(*), parameter :: &
sgl = [1,2,3,4,5,6,7,8] !< signed (single)
integer, dimension(5), parameter :: &
integer, dimension(*), parameter :: &
edg = [1,2,5,6,9], & !< edge
scr = [3,4,7,8,10] !< screw
integer, dimension(4), parameter :: &
integer, dimension(*), parameter :: &
mob = [1,2,3,4], & !< mobile
imm = [5,6,7,8] !< immobile (blocked)
integer, dimension(2), parameter :: &
integer, dimension(*), parameter :: &
dip = [9,10], & !< dipole
imm_edg = imm(1:2), & !< immobile edge
imm_scr = imm(3:4) !< immobile screw
@ -1612,7 +1612,7 @@ end subroutine stateInit
!--------------------------------------------------------------------------------------------------
!> @brief calculates kinetics
!--------------------------------------------------------------------------------------------------
subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Temperature, instance)
pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, Temperature, instance)
integer, intent(in) :: &
c, & !< dislocation character (1:edge, 2:screw)
@ -1727,7 +1727,7 @@ end subroutine kinetics
!> @brief returns copy of current dislocation densities from state
!> @details raw values is rectified
!--------------------------------------------------------------------------------------------------
function getRho(instance,of,ip,el)
pure function getRho(instance,of,ip,el)
integer, intent(in) :: instance, of,ip,el
real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho
@ -1752,7 +1752,7 @@ end function getRho
!> @brief returns copy of current dislocation densities from state
!> @details raw values is rectified
!--------------------------------------------------------------------------------------------------
function getRho0(instance,of,ip,el)
pure function getRho0(instance,of,ip,el)
integer, intent(in) :: instance, of,ip,el
real(pReal), dimension(param(instance)%sum_N_sl,10) :: getRho0

View File

@ -301,7 +301,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
e, & !< counter in element loop
startIP, endIP, &
s
logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false in hase of different Ngrains
logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false in has of different Ngrains
#ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 &
@ -1626,6 +1626,59 @@ logical pure function converged(residuum,state,atol)
end function converged
!--------------------------------------------------------------------------------------------------
!> @brief calculates a jump in the state according to the current state and the current stress
!> returns true, if state jump was successfull or not needed. false indicates NaN in delta state
!--------------------------------------------------------------------------------------------------
logical function stateJump(ipc,ip,el)
integer, intent(in):: &
el, & ! element index
ip, & ! integration point index
ipc ! grain index
integer :: &
c, &
p, &
mySource, &
myOffset, &
mySize
c = material_phaseMemberAt(ipc,ip,el)
p = material_phaseAt(ipc,el)
call constitutive_collectDeltaState(crystallite_S(1:3,1:3,ipc,ip,el), &
crystallite_Fe(1:3,1:3,ipc,ip,el), &
crystallite_Fi(1:3,1:3,ipc,ip,el), &
ipc,ip,el)
myOffset = plasticState(p)%offsetDeltaState
mySize = plasticState(p)%sizeDeltaState
if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c)))) then
stateJump = .false.
return
endif
plasticState(p)%state(myOffset + 1:myOffset + mySize,c) = &
plasticState(p)%state(myOffset + 1:myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c)
do mySource = 1, phase_Nsources(p)
myOffset = sourceState(p)%p(mySource)%offsetDeltaState
mySize = sourceState(p)%p(mySource)%sizeDeltaState
if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySize,c)))) then
stateJump = .false.
return
endif
sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) = &
sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) + sourceState(p)%p(mySource)%deltaState(1:mySize,c)
enddo
stateJump = .true.
end function stateJump
!--------------------------------------------------------------------------------------------------
!> @brief Write current restart information (Field and constitutive data) to file.
! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively

View File

@ -1056,14 +1056,15 @@ subroutine utilities_updateCoords(F)
call MPI_Irecv(IPfluct_padded(:,:,:,grid3+2),c,MPI_DOUBLE,rank_t,0,PETSC_COMM_WORLD,r,ierr)
if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Irecv')
call MPI_Wait(r,s,ierr)
if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Wait')
! send top layer to process above
if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Wait')
call MPI_Isend(IPfluct_padded(:,:,:,grid3+1),c,MPI_DOUBLE,rank_t,0,PETSC_COMM_WORLD,r,ierr)
if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Isend')
call MPI_Irecv(IPfluct_padded(:,:,:,1), c,MPI_DOUBLE,rank_b,0,PETSC_COMM_WORLD,r,ierr)
if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Irecv')
call MPI_Wait(r,s,ierr)
if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords/MPI_Wait')
!--------------------------------------------------------------------------------------------------
! calculate nodal displacements

View File

@ -29,6 +29,8 @@ module homogenization
!--------------------------------------------------------------------------------------------------
! General variables for the homogenization at a material point
logical, public :: &
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
@ -171,22 +173,6 @@ subroutine homogenization_init
write(6,'(/,a)') ' <<<+- homogenization init -+>>>'; flush(6)
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F0: ', shape(materialpoint_F0)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F: ', shape(materialpoint_F)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subF0: ', shape(materialpoint_subF0)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subF: ', shape(materialpoint_subF)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_P: ', shape(materialpoint_P)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subFrac: ', shape(materialpoint_subFrac)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subStep: ', shape(materialpoint_subStep)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subdt: ', shape(materialpoint_subdt)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_requested: ', shape(materialpoint_requested)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_converged: ', shape(materialpoint_converged)
write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy)
endif
flush(6)
if (debug_g < 1 .or. debug_g > homogenization_Ngrains(material_homogenizationAt(debug_e))) &
call IO_error(602,ext_msg='constituent', el=debug_e, g=debug_g)

View File

@ -18,16 +18,16 @@ module lattice
!--------------------------------------------------------------------------------------------------
! face centered cubic
integer, dimension(2), parameter :: &
integer, dimension(*), parameter :: &
FCC_NSLIPSYSTEM = [12, 6] !< # of slip systems per family for fcc
integer, dimension(1), parameter :: &
integer, dimension(*), parameter :: &
FCC_NTWINSYSTEM = [12] !< # of twin systems per family for fcc
integer, dimension(1), parameter :: &
integer, dimension(*), parameter :: &
FCC_NTRANSSYSTEM = [12] !< # of transformation systems per family for fcc
integer, dimension(1), parameter :: &
integer, dimension(*), parameter :: &
FCC_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for fcc
integer, parameter :: &
@ -109,13 +109,13 @@ module lattice
!--------------------------------------------------------------------------------------------------
! body centered cubic
integer, dimension(2), parameter :: &
integer, dimension(*), parameter :: &
BCC_NSLIPSYSTEM = [12, 12] !< # of slip systems per family for bcc
integer, dimension(1), parameter :: &
integer, dimension(*), parameter :: &
BCC_NTWINSYSTEM = [12] !< # of twin systems per family for bcc
integer, dimension(1), parameter :: &
integer, dimension(*), parameter :: &
BCC_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for bcc
integer, parameter :: &
@ -187,10 +187,10 @@ module lattice
!--------------------------------------------------------------------------------------------------
! hexagonal
integer, dimension(6), parameter :: &
integer, dimension(*), parameter :: &
HEX_NSLIPSYSTEM = [3, 3, 3, 6, 12, 6] !< # of slip systems per family for hex
integer, dimension(4), parameter :: &
integer, dimension(*), parameter :: &
HEX_NTWINSYSTEM = [6, 6, 6, 6] !< # of slip systems per family for hex
integer, parameter :: &
@ -280,7 +280,7 @@ module lattice
!--------------------------------------------------------------------------------------------------
! body centered tetragonal
integer, dimension(13), parameter :: &
integer, dimension(*), parameter :: &
BCT_NSLIPSYSTEM = [2, 2, 2, 4, 2, 4, 2, 2, 4, 8, 4, 8, 8 ] !< # of slip systems per family for bct (Sn) Bieler J. Electr Mater 2009
integer, parameter :: &
@ -362,7 +362,7 @@ module lattice
!--------------------------------------------------------------------------------------------------
! orthorhombic
integer, dimension(3), parameter :: &
integer, dimension(*), parameter :: &
ORT_NCLEAVAGESYSTEM = [1, 1, 1] !< # of cleavage systems per family for ortho
integer, parameter :: &

View File

@ -30,14 +30,12 @@ module math
1.0_pReal,0.0_pReal,0.0_pReal, &
0.0_pReal,1.0_pReal,0.0_pReal, &
0.0_pReal,0.0_pReal,1.0_pReal &
],[3,3]) !< 3x3 Identity
],shape(math_I3)) !< 3x3 Identity
real(pReal), dimension(6), parameter, private :: &
NRMMANDEL = [&
1.0_pReal, 1.0_pReal, 1.0_pReal, &
sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal) ] !< forward weighting for Mandel notation
real(pReal), dimension(*), parameter, private :: &
NRMMANDEL = [1.0_pReal, 1.0_pReal,1.0_pReal, sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal)] !< forward weighting for Mandel notation
real(pReal), dimension(6), parameter, private :: &
real(pReal), dimension(*), parameter, private :: &
INVNRMMANDEL = 1.0_pReal/NRMMANDEL !< backward weighting for Mandel notation
integer, dimension (2,6), parameter, private :: &
@ -48,7 +46,7 @@ module math
1,2, &
2,3, &
1,3 &
],[2,6]) !< arrangement in Nye notation.
],shape(MAPNYE)) !< arrangement in Nye notation.
integer, dimension (2,6), parameter, private :: &
MAPVOIGT = reshape([&
@ -58,7 +56,7 @@ module math
2,3, &
1,3, &
1,2 &
],[2,6]) !< arrangement in Voigt notation
],shape(MAPVOIGT)) !< arrangement in Voigt notation
integer, dimension (2,9), parameter, private :: &
MAPPLAIN = reshape([&
@ -71,7 +69,7 @@ module math
3,1, &
3,2, &
3,3 &
],[2,9]) !< arrangement in Plain notation
],shape(MAPPLAIN)) !< arrangement in Plain notation
interface math_eye
module procedure math_identity2nd
@ -487,21 +485,19 @@ function math_invSym3333(A)
real(pReal),dimension(3,3,3,3),intent(in) :: A
integer :: ierr
integer, dimension(6) :: ipiv6
real(pReal), dimension(6,6) :: temp66
real(pReal), dimension(6*(64+2)) :: work
logical :: error
integer :: ierr_i, ierr_f
external :: &
dgetrf, &
dgetri
temp66 = math_sym3333to66(A)
call dgetrf(6,6,temp66,6,ipiv6,ierr)
error = (ierr /= 0)
call dgetri(6,temp66,6,ipiv6,work,size(work,1),ierr)
error = error .or. (ierr /= 0)
if (error) then
call dgetrf(6,6,temp66,6,ipiv6,ierr_i)
call dgetri(6,temp66,6,ipiv6,work,size(work,1),ierr_f)
if (ierr_i /= 0 .or. ierr_f /= 0) then
call IO_error(400, ext_msg = 'math_invSym3333')
else
math_invSym3333 = math_66toSym3333(temp66)