From e7ca371b14a0fea1a64bf57cd0fdc9d682dcc01e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 May 2019 19:54:57 +0200 Subject: [PATCH 01/44] routines to do pointwise operations --- python/damask/dadf5.py | 129 +++++++++++++++++++++++++++++++++++++++++ python/damask/util.py | 46 +++++++++++++++ 2 files changed, 175 insertions(+) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index f5fc81b16..28cfea4b5 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -2,6 +2,8 @@ import h5py import re import numpy as np +from queue import Queue +from . import util # ------------------------------------------------------------------ class DADF5(): @@ -15,6 +17,7 @@ class DADF5(): if mode not in ['a','r']: print('Invalid file access mode') + else: with h5py.File(filename,mode): pass @@ -64,8 +67,33 @@ class DADF5(): self.filename = filename self.mode = mode + + def get_candidates(self,l): + groups = [] + if type(l) is not list: + print('mist') + with h5py.File(self.filename,'r') as f: + for g in self.get_active_groups(): + if set(l).issubset(f[g].keys()): groups.append(g) + return groups + def get_active_groups(self): + groups = [] + for i,x in enumerate(self.active['increments']): + group_inc = 'inc{:05}'.format(self.active['increments'][i]['inc']) + for c in self.active['constituents']: + group_constituent = group_inc+'/constituent/'+c + for t in self.active['c_output_types']: + group_output_types = group_constituent+'/'+t + groups.append(group_output_types) + for m in self.active['materialpoints']: + group_materialpoint = group_inc+'/materialpoint/'+m + for t in self.active['m_output_types']: + group_output_types = group_materialpoint+'/'+t + groups.append(group_output_types) + return groups + def list_data(self): """Shows information on all datasets in the file""" with h5py.File(self.filename,'r') as f: @@ -153,5 +181,106 @@ class DADF5(): print('unable to read materialpoint: '+ str(e)) return dataset + + + def add_Cauchy(self,PK2='P',F='F'): + + def Cauchy(F,P): + return 1.0/np.linalg.det(F)*np.dot(P,F.T) + + args = [{'label':F, 'shape':[3,3],'unit':'-'}, + {'label':PK2,'shape':[3,3],'unit':'Pa'} ] + result = {'label':'Cauchy','unit':'Pa'} + + self.add_generic_pointwise(Cauchy,args,result) + + + def add_Mises_stress(self,stress='Cauchy'): + + def Mises_stress(stress): + dev = stress - np.trace(stress)/3.0*np.eye(3) + symdev = 0.5*(dev+dev.T) + return np.sqrt(np.sum(symdev*symdev.T)*3.0/2.0) + + args = [{'label':stress,'shape':[3,3],'unit':'Pa'}] + result = {'label':'Mises({})'.format(stress),'unit':'Pa'} + + self.add_generic_pointwise(Mises_stress,args,result) + + def get_fitting(self,data): + groups = [] + if type(data) is not list: + print('mist') + with h5py.File(self.filename,'r') as f: + for g in self.get_candidates([l['label'] for l in data]): + print(g) + fits = True + for d in data: + fits = fits and np.all(np.array(f[g+'/'+d['label']].shape[1:]) == np.array(d['shape'])) # ToDo: allow here shape none and check for unit + if fits: groups.append(g) + return groups + + + def add_generic_pointwise(self,func,args,result): + """ + Ggeneral function to add pointwise data + + function 'func' first needs to have data arguments before other arguments + """ + groups = self.get_fitting(args) + + def job(args): + out = args['out'] + datasets_in = args['dat'] + func = args['fun'] + # calling the function per point might be performance-wise not optimal + # could be worth to investigate the performance for vectorized add_XXX functions that do the + # loops internally + for i in range(out.shape[0]): + arg = tuple([d[i,] for d in datasets_in]) + out[i,] = func(*arg) + args['results'].put({'out':out,'group':args['group']}) + + Nthreads = 4 # ToDo: should be a parameter + results = Queue(Nthreads+1) + + todo = [] + + for g in groups: + with h5py.File(self.filename,'r') as f: + datasets_in = [f[g+'/'+u['label']][()] for u in args] + + # figure out dimension of results + testArg = tuple([d[0,] for d in datasets_in]) # to call function with first point + out = np.empty([datasets_in[0].shape[0]] + list(func(*testArg).shape)) # shape is Npoints x shape of the results for one point + todo.append({'dat':datasets_in,'fun':func,'out':out,'group':g,'results':results}) + + # Instantiate a thread pool with worker threads + pool = util.ThreadPool(Nthreads) + missingResults = len(todo) + + + # Add the jobs in bulk to the thread pool. Alternatively you could use + # `pool.add_task` to add single jobs. The code will block here, which + # makes it possible to cancel the thread pool with an exception when + # the currently running batch of workers is finished. + + pool.map(job, todo[:Nthreads+1]) + i = 0 + while missingResults > 0: + r=results.get() # noqa + print(r['group']) + with h5py.File(self.filename,'r+') as f: + dataset_out = f[r['group']].create_dataset(result['label'],data=r['out']) + dataset_out.attrs['unit'] = result['unit'] + missingResults-=1 + try: + pool.add_task(job,todo[Nthreads+1+i]) + except: + pass + i+=1 + + pool.wait_completion() + diff --git a/python/damask/util.py b/python/damask/util.py index 02cb4a2c6..e8bb4b3ec 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -2,6 +2,9 @@ import sys,time,os,subprocess,shlex import numpy as np from optparse import Option +from queue import Queue +from threading import Thread + class bcolors: """ @@ -425,3 +428,46 @@ def curve_fit_bound(f, xdata, ydata, p0=None, sigma=None, bounds=None, **kw): pcov = np.inf return (popt, pcov, infodict, errmsg, ier) if return_full else (popt, pcov) + +class Worker(Thread): + """Thread executing tasks from a given tasks queue""" + + def __init__(self, tasks): + Thread.__init__(self) + self.tasks = tasks + self.daemon = True + self.start() + + def run(self): + while True: + func, args, kargs = self.tasks.get() + try: + func(*args, **kargs) + except Exception as e: + # An exception happened in this thread + print(e) + finally: + # Mark this task as done, whether an exception happened or not + self.tasks.task_done() + + +class ThreadPool: + """Pool of threads consuming tasks from a queue""" + + def __init__(self, num_threads): + self.tasks = Queue(num_threads) + for _ in range(num_threads): + Worker(self.tasks) + + def add_task(self, func, *args, **kargs): + """Add a task to the queue""" + self.tasks.put((func, args, kargs)) + + def map(self, func, args_list): + """Add a list of tasks to the queue""" + for args in args_list: + self.add_task(func, args) + + def wait_completion(self): + """Wait for completion of all the tasks in the queue""" + self.tasks.join() From f6dda99bfb39cc209e88d16486096edc83a299cc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 23 May 2019 08:54:20 +0200 Subject: [PATCH 02/44] more post processing functionality --- python/damask/dadf5.py | 43 +++++++++++++++++++++++++++++++++--------- 1 file changed, 34 insertions(+), 9 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 28cfea4b5..984a81576 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -183,14 +183,16 @@ class DADF5(): return dataset - def add_Cauchy(self,PK2='P',F='F'): + def add_Cauchy(self,P='P',F='F'): def Cauchy(F,P): return 1.0/np.linalg.det(F)*np.dot(P,F.T) - args = [{'label':F, 'shape':[3,3],'unit':'-'}, - {'label':PK2,'shape':[3,3],'unit':'Pa'} ] - result = {'label':'Cauchy','unit':'Pa'} + args = [{'label':F,'shape':[3,3],'unit':'-'}, + {'label':P,'shape':[3,3],'unit':'Pa'} ] + result = {'label':'sigma', + 'unit':'Pa', + 'Description': 'Cauchy stress calculated from 1st Piola-Kirchhoff stress and deformation gradient'} self.add_generic_pointwise(Cauchy,args,result) @@ -203,9 +205,29 @@ class DADF5(): return np.sqrt(np.sum(symdev*symdev.T)*3.0/2.0) args = [{'label':stress,'shape':[3,3],'unit':'Pa'}] - result = {'label':'Mises({})'.format(stress),'unit':'Pa'} + result = {'label':'Mises({})'.format(stress), + 'unit':'Pa', + 'Description': 'Equivalent Mises stress'} self.add_generic_pointwise(Mises_stress,args,result) + + def add_norm(self,x,ord=None): + # ToDo: The output unit should be the input unit + args = [{'label':x,'shape':None,'unit':None}] + result = {'label':'norm_{}({})'.format(str(ord),x), + 'unit':'n/a', + 'Description': 'Norm of vector or tensor or magnitude of a scalar. See numpy.linalg.norm manual for details'} + + self.add_generic_pointwise(np.linalg.norm,args,result) + + def add_determinant(self,a): + # ToDo: The output unit should be the input unit + args = [{'label':a,'shape':[3,3],'unit':None}] + result = {'label':'det({})'.format(a), + 'unit':'n/a', + 'Description': 'Determinan of a tensor'} + + self.add_generic_pointwise(np.linalg.det,args,result) def get_fitting(self,data): groups = [] @@ -215,8 +237,9 @@ class DADF5(): for g in self.get_candidates([l['label'] for l in data]): print(g) fits = True - for d in data: - fits = fits and np.all(np.array(f[g+'/'+d['label']].shape[1:]) == np.array(d['shape'])) # ToDo: allow here shape none and check for unit + for d in data: # ToDo: check for unit + if d['shape'] is not None: + fits = fits and np.all(np.array(f[g+'/'+d['label']].shape[1:]) == np.array(d['shape'])) if fits: groups.append(g) return groups @@ -263,7 +286,7 @@ class DADF5(): # Add the jobs in bulk to the thread pool. Alternatively you could use # `pool.add_task` to add single jobs. The code will block here, which # makes it possible to cancel the thread pool with an exception when - # the currently running batch of workers is finished. + # the currently running batch of workers is finishnumpy.linalg.normed. pool.map(job, todo[:Nthreads+1]) i = 0 @@ -272,7 +295,9 @@ class DADF5(): print(r['group']) with h5py.File(self.filename,'r+') as f: dataset_out = f[r['group']].create_dataset(result['label'],data=r['out']) - dataset_out.attrs['unit'] = result['unit'] + dataset_out.attrs['Unit'] = result['unit'] + dataset_out.attrs['Description'] = result['Description'] + dataset_out.attrs['Creator'] = 'dadf5.py v{}'.format('n/a') missingResults-=1 try: pool.add_task(job,todo[Nthreads+1+i]) From 508f69691b60846d94f39c4eac91298f2df1edf3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 23 May 2019 10:10:40 +0200 Subject: [PATCH 03/44] relax tolerance for nonlocal --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 3a2f89547..d31da38cf 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 3a2f89547c264044a7bfab9d33aee78eec495a76 +Subproject commit d31da38cf25734a91e994a3d5d33bb048eb2f44f From a280f9a4a20ef0b0c9efa1845c94392d681d5a00 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 23 May 2019 17:22:57 +0200 Subject: [PATCH 04/44] polishing --- python/damask/dadf5.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 984a81576..716747288 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -184,7 +184,7 @@ class DADF5(): def add_Cauchy(self,P='P',F='F'): - + """Adds Cauchy stress calculated from 1st Piola-Kirchhoff stress and deformation gradient""" def Cauchy(F,P): return 1.0/np.linalg.det(F)*np.dot(P,F.T) @@ -197,8 +197,8 @@ class DADF5(): self.add_generic_pointwise(Cauchy,args,result) - def add_Mises_stress(self,stress='Cauchy'): - + def add_Mises_stress(self,stress='sigma'): + """Adds equivalent von Mises stress""" def Mises_stress(stress): dev = stress - np.trace(stress)/3.0*np.eye(3) symdev = 0.5*(dev+dev.T) @@ -211,7 +211,9 @@ class DADF5(): self.add_generic_pointwise(Mises_stress,args,result) + def add_norm(self,x,ord=None): + """Adds norm of vector or tensor or magnitude of a scalar.""" # ToDo: The output unit should be the input unit args = [{'label':x,'shape':None,'unit':None}] result = {'label':'norm_{}({})'.format(str(ord),x), @@ -219,15 +221,18 @@ class DADF5(): 'Description': 'Norm of vector or tensor or magnitude of a scalar. See numpy.linalg.norm manual for details'} self.add_generic_pointwise(np.linalg.norm,args,result) - + + def add_determinant(self,a): + """Adds the determinan of a tensor""" # ToDo: The output unit should be the input unit args = [{'label':a,'shape':[3,3],'unit':None}] result = {'label':'det({})'.format(a), 'unit':'n/a', - 'Description': 'Determinan of a tensor'} + 'Description': 'Determinant of a tensor'} self.add_generic_pointwise(np.linalg.det,args,result) + def get_fitting(self,data): groups = [] From dbfaf66472bbabe2f26b1befe14ab47ca0c1dc72 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 31 May 2019 20:14:35 +0200 Subject: [PATCH 05/44] pInt not needed --- src/mesh_grid.f90 | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index d873e3542..df4f200f7 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -13,18 +13,18 @@ module mesh implicit none private - integer(pInt), public, protected :: & + integer, public, protected :: & mesh_Nnodes - integer(pInt), dimension(:), allocatable, private :: & + integer, dimension(:), allocatable, private :: & microGlobal - integer(pInt), dimension(:), allocatable, private :: & + integer, dimension(:), allocatable, private :: & mesh_homogenizationAt - integer(pInt), dimension(:,:), allocatable, public, protected :: & + integer, dimension(:,:), allocatable, public, protected :: & mesh_element !< entryCount and list of elements containing node - integer(pInt), dimension(:,:,:,:), allocatable, public, protected :: & + integer, dimension(:,:,:,:), allocatable, public, protected :: & mesh_ipNeighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] real(pReal), public, protected :: & @@ -51,9 +51,9 @@ module mesh ! grid specific - integer(pInt), dimension(3), public, protected :: & + integer, dimension(3), public, protected :: & grid !< (global) grid - integer(pInt), public, protected :: & + integer, public, protected :: & mesh_NcpElemsGlobal, & !< total number of CP elements in global mesh grid3, & !< (local) grid in 3rd direction grid3Offset !< (local) grid offset in 3rd direction @@ -76,9 +76,9 @@ module mesh type, public, extends(tMesh) :: tMesh_grid - integer(pInt), dimension(3), public :: & + integer, dimension(3), public :: & grid !< (global) grid - integer(pInt), public :: & + integer, public :: & mesh_NcpElemsGlobal, & !< total number of CP elements in global mesh grid3, & !< (local) grid in 3rd direction grid3Offset !< (local) grid offset in 3rd direction @@ -135,7 +135,7 @@ subroutine mesh_init(ip,el) include 'fftw3-mpi.f03' integer(C_INTPTR_T) :: devNull, local_K, local_K_offset integer :: ierr, worldsize, j - integer(pInt), intent(in), optional :: el, ip + integer, intent(in), optional :: el, ip logical :: myDebug write(6,'(/,a)') ' <<<+- mesh init -+>>>' @@ -233,9 +233,9 @@ subroutine mesh_spectral_read_grid() implicit none character(len=:), allocatable :: rawData character(len=65536) :: line - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: h =- 1_pInt - integer(pInt) :: & + integer, allocatable, dimension(:) :: chunkPos + integer :: h =- 1_pInt + integer :: & headerLength = -1_pInt, & !< length of header (in lines) fileLength, & !< length of the geom file (in characters) fileUnit, & @@ -416,7 +416,7 @@ end function mesh_build_ipCoordinates !! Allocates global array 'mesh_element' !-------------------------------------------------------------------------------------------------- subroutine mesh_spectral_build_elements() - integer(pInt) :: & + integer :: & e, & elemOffset @@ -450,7 +450,7 @@ end subroutine mesh_spectral_build_elements subroutine mesh_spectral_build_ipNeighborhood implicit none - integer(pInt) :: & + integer :: & x,y,z, & e allocate(mesh_ipNeighborhood(3,6,1,theMesh%nElems),source=0_pInt) @@ -518,16 +518,16 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes) real(pReal), dimension(3,size(centres,2)+2,size(centres,3)+2,size(centres,4)+2) :: & wrappedCentres - integer(pInt) :: & + integer :: & i,j,k,n - integer(pInt), dimension(3), parameter :: & + integer, dimension(3), parameter :: & diag = 1_pInt - integer(pInt), dimension(3) :: & + integer, dimension(3) :: & shift = 0_pInt, & lookup = 0_pInt, & me = 0_pInt, & iRes = 0_pInt - integer(pInt), dimension(3,8) :: & + integer, dimension(3,8) :: & neighbor = reshape([ & 0_pInt, 0_pInt, 0_pInt, & 1_pInt, 0_pInt, 0_pInt, & From 95d93ce58c8c802b25c3764afe20fff705c0f419 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 31 May 2019 20:17:21 +0200 Subject: [PATCH 06/44] one implicit none is enough --- src/mesh_grid.f90 | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index df4f200f7..cb08a6482 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -13,6 +13,7 @@ module mesh implicit none private + integer, public, protected :: & mesh_Nnodes @@ -99,7 +100,6 @@ contains subroutine tMesh_grid_init(self,nodes) - implicit none class(tMesh_grid) :: self real(pReal), dimension(:,:), intent(in) :: nodes @@ -131,7 +131,6 @@ subroutine mesh_init(ip,el) FEsolving_execElem, & FEsolving_execIP - implicit none include 'fftw3-mpi.f03' integer(C_INTPTR_T) :: devNull, local_K, local_K_offset integer :: ierr, worldsize, j @@ -230,7 +229,6 @@ subroutine mesh_spectral_read_grid() use DAMASK_interface, only: & geometryFile - implicit none character(len=:), allocatable :: rawData character(len=65536) :: line integer, allocatable, dimension(:) :: chunkPos @@ -449,7 +447,6 @@ end subroutine mesh_spectral_build_elements !-------------------------------------------------------------------------------------------------- subroutine mesh_spectral_build_ipNeighborhood - implicit none integer :: & x,y,z, & e From c045bd6355716e3e03489f3bcdd8346877541e19 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 31 May 2019 20:48:16 +0200 Subject: [PATCH 07/44] cleaning --- src/mesh_grid.f90 | 649 ++++++++++++++++++++++------------------------ 1 file changed, 309 insertions(+), 340 deletions(-) diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index cb08a6482..17eb677d2 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -6,95 +6,104 @@ !> @brief Sets up the mesh for the solvers MSC.Marc, Abaqus and the spectral solver !-------------------------------------------------------------------------------------------------- module mesh - use, intrinsic :: iso_c_binding - use prec - use geometry_plastic_nonlocal - use mesh_base + use, intrinsic :: iso_c_binding + use prec + use geometry_plastic_nonlocal + use mesh_base +#include + use PETScsys + use DAMASK_interface + use IO + use debug + use numerics + use FEsolving - implicit none - private - - integer, public, protected :: & - mesh_Nnodes - - integer, dimension(:), allocatable, private :: & - microGlobal - integer, dimension(:), allocatable, private :: & - mesh_homogenizationAt - - integer, dimension(:,:), allocatable, public, protected :: & - mesh_element !< entryCount and list of elements containing node - - integer, dimension(:,:,:,:), allocatable, public, protected :: & - mesh_ipNeighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] - - real(pReal), public, protected :: & - mesh_unitlength !< physical length of one unit in mesh - - real(pReal), dimension(:,:), allocatable, private :: & - mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) - - - real(pReal), dimension(:,:), allocatable, public, protected :: & - mesh_ipVolume, & !< volume associated with IP (initially!) - mesh_node0 !< node x,y,z coordinates (initially!) - - real(pReal), dimension(:,:,:), allocatable, public, protected :: & - mesh_ipArea !< area of interface to neighboring IP (initially!) - - real(pReal), dimension(:,:,:), allocatable, public :: & - mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) - - real(pReal),dimension(:,:,:,:), allocatable, public, protected :: & - mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) - - logical, dimension(3), public, parameter :: mesh_periodicSurface = .true. !< flag indicating periodic outer surfaces (used for fluxes) - - -! grid specific - integer, dimension(3), public, protected :: & - grid !< (global) grid - integer, public, protected :: & - mesh_NcpElemsGlobal, & !< total number of CP elements in global mesh - grid3, & !< (local) grid in 3rd direction - grid3Offset !< (local) grid offset in 3rd direction - real(pReal), dimension(3), public, protected :: & - geomSize - real(pReal), public, protected :: & - size3, & !< (local) size in 3rd direction - size3offset !< (local) size offset in 3rd direction - - public :: & - mesh_init - - private :: & - mesh_build_ipAreas, & - mesh_build_ipNormals, & - mesh_spectral_build_nodes, & - mesh_spectral_build_elements, & - mesh_spectral_build_ipNeighborhood, & - mesh_build_ipCoordinates - - type, public, extends(tMesh) :: tMesh_grid - integer, dimension(3), public :: & - grid !< (global) grid - integer, public :: & - mesh_NcpElemsGlobal, & !< total number of CP elements in global mesh - grid3, & !< (local) grid in 3rd direction - grid3Offset !< (local) grid offset in 3rd direction - real(pReal), dimension(3), public :: & - geomSize - real(pReal), public :: & - size3, & !< (local) size in 3rd direction - size3offset - - contains - procedure, pass(self) :: tMesh_grid_init - generic, public :: init => tMesh_grid_init - end type tMesh_grid + implicit none + private - type(tMesh_grid), public, protected :: theMesh + include 'fftw3-mpi.f03' + integer, public, protected :: & + mesh_Nnodes + + integer, dimension(:), allocatable, private :: & + microGlobal + integer, dimension(:), allocatable, private :: & + mesh_homogenizationAt + + integer, dimension(:,:), allocatable, public, protected :: & + mesh_element !< entryCount and list of elements containing node + + integer, dimension(:,:,:,:), allocatable, public, protected :: & + mesh_ipNeighborhood !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me] + + real(pReal), public, protected :: & + mesh_unitlength !< physical length of one unit in mesh + + real(pReal), dimension(:,:), allocatable, private :: & + mesh_node !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!) + + + real(pReal), dimension(:,:), allocatable, public, protected :: & + mesh_ipVolume, & !< volume associated with IP (initially!) + mesh_node0 !< node x,y,z coordinates (initially!) + + real(pReal), dimension(:,:,:), allocatable, public, protected :: & + mesh_ipArea !< area of interface to neighboring IP (initially!) + + real(pReal), dimension(:,:,:), allocatable, public :: & + mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) + + real(pReal),dimension(:,:,:,:), allocatable, public, protected :: & + mesh_ipAreaNormal !< area normal of interface to neighboring IP (initially!) + + logical, dimension(3), public, parameter :: mesh_periodicSurface = .true. !< flag indicating periodic outer surfaces (used for fluxes) + + + ! grid specific + integer, dimension(3), public, protected :: & + grid !< (global) grid + integer, public, protected :: & + mesh_NcpElemsGlobal, & !< total number of CP elements in global mesh + grid3, & !< (local) grid in 3rd direction + grid3Offset !< (local) grid offset in 3rd direction + real(pReal), dimension(3), public, protected :: & + geomSize + real(pReal), public, protected :: & + size3, & !< (local) size in 3rd direction + size3offset !< (local) size offset in 3rd direction + + public :: & + mesh_init + + private :: & + mesh_build_ipAreas, & + mesh_build_ipNormals, & + mesh_spectral_build_nodes, & + mesh_spectral_build_elements, & + mesh_spectral_build_ipNeighborhood, & + mesh_build_ipCoordinates + + type, public, extends(tMesh) :: tMesh_grid + + integer, dimension(3), public :: & + grid !< (global) grid + integer, public :: & + mesh_NcpElemsGlobal, & !< total number of CP elements in global mesh + grid3, & !< (local) grid in 3rd direction + grid3Offset !< (local) grid offset in 3rd direction + real(pReal), dimension(3), public :: & + geomSize + real(pReal), public :: & + size3, & !< (local) size in 3rd direction + size3offset + + contains + procedure, pass(self) :: tMesh_grid_init + generic, public :: init => tMesh_grid_init + end type tMesh_grid + + type(tMesh_grid), public, protected :: theMesh contains @@ -103,7 +112,7 @@ subroutine tMesh_grid_init(self,nodes) class(tMesh_grid) :: self real(pReal), dimension(:,:), intent(in) :: nodes - call self%tMesh%init('grid',10_pInt,nodes) + call self%tMesh%init('grid',10,nodes) end subroutine tMesh_grid_init @@ -113,101 +122,82 @@ end subroutine tMesh_grid_init !-------------------------------------------------------------------------------------------------- subroutine mesh_init(ip,el) -#include - use PETScsys + integer(C_INTPTR_T) :: devNull, local_K, local_K_offset + integer :: ierr, worldsize, j + integer, intent(in), optional :: el, ip + logical :: myDebug - use DAMASK_interface - use IO, only: & - IO_error - use debug, only: & - debug_e, & - debug_i, & - debug_level, & - debug_mesh, & - debug_levelBasic - use numerics, only: & - numerics_unitlength - use FEsolving, only: & - FEsolving_execElem, & - FEsolving_execIP + write(6,'(/,a)') ' <<<+- mesh init -+>>>' - include 'fftw3-mpi.f03' - integer(C_INTPTR_T) :: devNull, local_K, local_K_offset - integer :: ierr, worldsize, j - integer, intent(in), optional :: el, ip - logical :: myDebug + mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh - write(6,'(/,a)') ' <<<+- mesh init -+>>>' + myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0) - mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh - - myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) - - call fftw_mpi_init() - call mesh_spectral_read_grid() + call fftw_mpi_init() + call mesh_spectral_read_grid() - call MPI_comm_size(PETSC_COMM_WORLD, worldsize, ierr) - if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_comm_size') - if(worldsize>grid(3)) call IO_error(894_pInt, ext_msg='number of processes exceeds grid(3)') + call MPI_comm_size(PETSC_COMM_WORLD, worldsize, ierr) + if(ierr /=0) call IO_error(894, ext_msg='MPI_comm_size') + if(worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)') - devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T), & - int(grid(2),C_INTPTR_T), & - int(grid(1),C_INTPTR_T)/2+1, & - PETSC_COMM_WORLD, & - local_K, & ! domain grid size along z - local_K_offset) ! domain grid offset along z - grid3 = int(local_K,pInt) - grid3Offset = int(local_K_offset,pInt) - size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal) - size3Offset = geomSize(3)*real(grid3Offset,pReal)/real(grid(3),pReal) + devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T), & + int(grid(2),C_INTPTR_T), & + int(grid(1),C_INTPTR_T)/2+1, & + PETSC_COMM_WORLD, & + local_K, & ! domain grid size along z + local_K_offset) ! domain grid offset along z + grid3 = int(local_K,pInt) + grid3Offset = int(local_K_offset,pInt) + size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal) + size3Offset = geomSize(3)*real(grid3Offset,pReal)/real(grid(3),pReal) - mesh_NcpElemsGlobal = product(grid) + mesh_NcpElemsGlobal = product(grid) - mesh_Nnodes = product(grid(1:2) + 1_pInt)*(grid3 + 1_pInt) + mesh_Nnodes = product(grid(1:2) + 1)*(grid3 + 1) - mesh_node0 = mesh_spectral_build_nodes() - mesh_node = mesh_node0 - if (myDebug) write(6,'(a)') ' Built nodes'; flush(6) + mesh_node0 = mesh_spectral_build_nodes() + mesh_node = mesh_node0 + if (myDebug) write(6,'(a)') ' Built nodes'; flush(6) - call theMesh%init(mesh_node) - call theMesh%setNelems(product(grid(1:2))*grid3) - call mesh_spectral_build_elements() - mesh_homogenizationAt = mesh_homogenizationAt(product(grid(1:2))*grid3Offset+1: & - product(grid(1:2))*(grid3Offset+grid3)) ! reallocate/shrink in case of MPI - - if (myDebug) write(6,'(a)') ' Built elements'; flush(6) - - + call theMesh%init(mesh_node) + call theMesh%setNelems(product(grid(1:2))*grid3) + call mesh_spectral_build_elements() + mesh_homogenizationAt = mesh_homogenizationAt(product(grid(1:2))*grid3Offset+1: & + product(grid(1:2))*(grid3Offset+grid3)) ! reallocate/shrink in case of MPI + + if (myDebug) write(6,'(a)') ' Built elements'; flush(6) + + - if (myDebug) write(6,'(a)') ' Built cell nodes'; flush(6) - mesh_ipCoordinates = mesh_build_ipCoordinates() - if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) - allocate(mesh_ipVolume(1,theMesh%nElems),source=product([geomSize(1:2),size3]/real([grid(1:2),grid3]))) - if (myDebug) write(6,'(a)') ' Built IP volumes'; flush(6) - mesh_ipArea = mesh_build_ipAreas() - mesh_ipAreaNormal = mesh_build_ipNormals() - if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) + if (myDebug) write(6,'(a)') ' Built cell nodes'; flush(6) + mesh_ipCoordinates = mesh_build_ipCoordinates() + if (myDebug) write(6,'(a)') ' Built IP coordinates'; flush(6) + allocate(mesh_ipVolume(1,theMesh%nElems),source=product([geomSize(1:2),size3]/real([grid(1:2),grid3]))) + if (myDebug) write(6,'(a)') ' Built IP volumes'; flush(6) + mesh_ipArea = mesh_build_ipAreas() + mesh_ipAreaNormal = mesh_build_ipNormals() + if (myDebug) write(6,'(a)') ' Built IP areas'; flush(6) - call mesh_spectral_build_ipNeighborhood - call geometry_plastic_nonlocal_set_IPneighborhood(mesh_ipNeighborhood) + call mesh_spectral_build_ipNeighborhood + call geometry_plastic_nonlocal_set_IPneighborhood(mesh_ipNeighborhood) - if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6) + if (myDebug) write(6,'(a)') ' Built IP neighborhood'; flush(6) - if (debug_e < 1 .or. debug_e > theMesh%nElems) & - call IO_error(602_pInt,ext_msg='element') ! selected element does not exist - if (debug_i < 1 .or. debug_i > theMesh%elem%nIPs) & - call IO_error(602_pInt,ext_msg='IP') ! selected element does not have requested IP + if (debug_e < 1 .or. debug_e > theMesh%nElems) & + call IO_error(602,ext_msg='element') ! selected element does not exist + if (debug_i < 1 .or. debug_i > theMesh%elem%nIPs) & + call IO_error(602,ext_msg='IP') ! selected element does not have requested IP - FEsolving_execElem = [ 1_pInt,theMesh%nElems ] ! parallel loop bounds set to comprise all DAMASK elements - allocate(FEsolving_execIP(2_pInt,theMesh%nElems), source=1_pInt) ! parallel loop bounds set to comprise from first IP... - forall (j = 1_pInt:theMesh%nElems) FEsolving_execIP(2,j) = theMesh%elem%nIPs ! ...up to own IP count for each element + FEsolving_execElem = [ 1,theMesh%nElems ] ! parallel loop bounds set to comprise all DAMASK elements + allocate(FEsolving_execIP(2,theMesh%nElems), source=1) ! parallel loop bounds set to comprise from first IP... + forall (j = 1:theMesh%nElems) FEsolving_execIP(2,j) = theMesh%elem%nIPs ! ...up to own IP count for each element !!!! COMPATIBILITY HACK !!!! - theMesh%homogenizationAt = mesh_element(3,:) - theMesh%microstructureAt = mesh_element(4,:) + theMesh%homogenizationAt = mesh_element(3,:) + theMesh%microstructureAt = mesh_element(4,:) !!!!!!!!!!!!!!!!!!!!!!!! end subroutine mesh_init @@ -219,22 +209,13 @@ end subroutine mesh_init ! supposed to be called only once! !-------------------------------------------------------------------------------------------------- subroutine mesh_spectral_read_grid() - use IO, only: & - IO_stringPos, & - IO_lc, & - IO_stringValue, & - IO_intValue, & - IO_floatValue, & - IO_error - use DAMASK_interface, only: & - geometryFile character(len=:), allocatable :: rawData character(len=65536) :: line integer, allocatable, dimension(:) :: chunkPos - integer :: h =- 1_pInt + integer :: h =- 1 integer :: & - headerLength = -1_pInt, & !< length of header (in lines) + headerLength = -1, & !< length of header (in lines) fileLength, & !< length of the geom file (in characters) fileUnit, & startPos, endPos, & @@ -245,7 +226,7 @@ subroutine mesh_spectral_read_grid() e, & !< "element", i.e. spectral collocation point i, j - grid = -1_pInt + grid = -1 geomSize = -1.0_pReal !-------------------------------------------------------------------------------------------------- @@ -253,7 +234,7 @@ subroutine mesh_spectral_read_grid() inquire(file = trim(geometryFile), size=fileLength) open(newunit=fileUnit, file=trim(geometryFile), access='stream',& status='old', position='rewind', action='read',iostat=myStat) - if(myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=trim(geometryFile)) + if(myStat /= 0) call IO_error(100,ext_msg=trim(geometryFile)) allocate(character(len=fileLength)::rawData) read(fileUnit) rawData close(fileUnit) @@ -263,106 +244,106 @@ subroutine mesh_spectral_read_grid() endPos = index(rawData,new_line('')) if(endPos <= index(rawData,'head')) then startPos = len(rawData) - call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_read_grid') + call IO_error(error_ID=841, ext_msg='mesh_spectral_read_grid') else chunkPos = IO_stringPos(rawData(1:endPos)) - if (chunkPos(1) < 2_pInt) call IO_error(error_ID=841_pInt, ext_msg='mesh_spectral_read_grid') - headerLength = IO_intValue(rawData(1:endPos),chunkPos,1_pInt) - startPos = endPos + 1_pInt + if (chunkPos(1) < 2) call IO_error(error_ID=841, ext_msg='mesh_spectral_read_grid') + headerLength = IO_intValue(rawData(1:endPos),chunkPos,1) + startPos = endPos + 1 endif !-------------------------------------------------------------------------------------------------- ! read and interprete header l = 0 do while (l < headerLength .and. startPos < len(rawData)) - endPos = startPos + index(rawData(startPos:),new_line('')) - 1_pInt + endPos = startPos + index(rawData(startPos:),new_line('')) - 1 if (endPos < startPos) endPos = len(rawData) ! end of file without new line line = rawData(startPos:endPos) - startPos = endPos + 1_pInt - l = l + 1_pInt + startPos = endPos + 1 + l = l + 1 chunkPos = IO_stringPos(trim(line)) if (chunkPos(1) < 2) cycle ! need at least one keyword value pair - select case ( IO_lc(IO_StringValue(trim(line),chunkPos,1_pInt,.true.)) ) + select case ( IO_lc(IO_StringValue(trim(line),chunkPos,1,.true.)) ) case ('grid') if (chunkPos(1) > 6) then - do j = 2_pInt,6_pInt,2_pInt + do j = 2,6,2 select case (IO_lc(IO_stringValue(line,chunkPos,j))) case('a') - grid(1) = IO_intValue(line,chunkPos,j+1_pInt) + grid(1) = IO_intValue(line,chunkPos,j+1) case('b') - grid(2) = IO_intValue(line,chunkPos,j+1_pInt) + grid(2) = IO_intValue(line,chunkPos,j+1) case('c') - grid(3) = IO_intValue(line,chunkPos,j+1_pInt) + grid(3) = IO_intValue(line,chunkPos,j+1) end select enddo endif case ('size') if (chunkPos(1) > 6) then - do j = 2_pInt,6_pInt,2_pInt + do j = 2,6,2 select case (IO_lc(IO_stringValue(line,chunkPos,j))) case('x') - geomSize(1) = IO_floatValue(line,chunkPos,j+1_pInt) + geomSize(1) = IO_floatValue(line,chunkPos,j+1) case('y') - geomSize(2) = IO_floatValue(line,chunkPos,j+1_pInt) + geomSize(2) = IO_floatValue(line,chunkPos,j+1) case('z') - geomSize(3) = IO_floatValue(line,chunkPos,j+1_pInt) + geomSize(3) = IO_floatValue(line,chunkPos,j+1) end select enddo endif case ('homogenization') - if (chunkPos(1) > 1) h = IO_intValue(line,chunkPos,2_pInt) + if (chunkPos(1) > 1) h = IO_intValue(line,chunkPos,2) end select enddo !-------------------------------------------------------------------------------------------------- ! sanity checks - if(h < 1_pInt) & - call IO_error(error_ID = 842_pInt, ext_msg='homogenization (mesh_spectral_read_grid)') - if(any(grid < 1_pInt)) & - call IO_error(error_ID = 842_pInt, ext_msg='grid (mesh_spectral_read_grid)') + if(h < 1) & + call IO_error(error_ID = 842, ext_msg='homogenization (mesh_spectral_read_grid)') + if(any(grid < 1)) & + call IO_error(error_ID = 842, ext_msg='grid (mesh_spectral_read_grid)') if(any(geomSize < 0.0_pReal)) & - call IO_error(error_ID = 842_pInt, ext_msg='size (mesh_spectral_read_grid)') + call IO_error(error_ID = 842, ext_msg='size (mesh_spectral_read_grid)') - allocate(microGlobal(product(grid)), source = -1_pInt) + allocate(microGlobal(product(grid)), source = -1) allocate(mesh_homogenizationAt(product(grid)), source = h) ! too large in case of MPI (shrink later, not very elegant) !-------------------------------------------------------------------------------------------------- ! read and interpret content - e = 1_pInt + e = 1 do while (startPos < len(rawData)) - endPos = startPos + index(rawData(startPos:),new_line('')) - 1_pInt + endPos = startPos + index(rawData(startPos:),new_line('')) - 1 if (endPos < startPos) endPos = len(rawData) ! end of file without new line line = rawData(startPos:endPos) - startPos = endPos + 1_pInt - l = l + 1_pInt + startPos = endPos + 1 + l = l + 1 chunkPos = IO_stringPos(trim(line)) noCompression: if (chunkPos(1) /= 3) then c = chunkPos(1) - microGlobal(e:e+c-1_pInt) = [(IO_intValue(line,chunkPos,i+1_pInt), i=0_pInt, c-1_pInt)] + microGlobal(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)] else noCompression compression: if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'of') then c = IO_intValue(line,chunkPos,1) - microGlobal(e:e+c-1_pInt) = [(IO_intValue(line,chunkPos,3),i = 1_pInt,IO_intValue(line,chunkPos,1))] + microGlobal(e:e+c-1) = [(IO_intValue(line,chunkPos,3),i = 1,IO_intValue(line,chunkPos,1))] else if (IO_lc(IO_stringValue(line,chunkPos,2)) == 'to') then compression - c = abs(IO_intValue(line,chunkPos,3) - IO_intValue(line,chunkPos,1)) + 1_pInt - o = merge(+1_pInt, -1_pInt, IO_intValue(line,chunkPos,3) > IO_intValue(line,chunkPos,1)) - microGlobal(e:e+c-1_pInt) = [(i, i = IO_intValue(line,chunkPos,1),IO_intValue(line,chunkPos,3),o)] + c = abs(IO_intValue(line,chunkPos,3) - IO_intValue(line,chunkPos,1)) + 1 + o = merge(+1, -1, IO_intValue(line,chunkPos,3) > IO_intValue(line,chunkPos,1)) + microGlobal(e:e+c-1) = [(i, i = IO_intValue(line,chunkPos,1),IO_intValue(line,chunkPos,3),o)] else compression c = chunkPos(1) - microGlobal(e:e+c-1_pInt) = [(IO_intValue(line,chunkPos,i+1_pInt), i=0_pInt, c-1_pInt)] + microGlobal(e:e+c-1) = [(IO_intValue(line,chunkPos,i+1), i=0, c-1)] endif compression endif noCompression e = e+c end do - if (e-1 /= product(grid)) call IO_error(error_ID = 843_pInt, el=e) + if (e-1 /= product(grid)) call IO_error(error_ID = 843, el=e) end subroutine mesh_spectral_read_grid @@ -413,30 +394,30 @@ end function mesh_build_ipCoordinates !> @brief Store FEid, type, material, texture, and node list per element. !! Allocates global array 'mesh_element' !-------------------------------------------------------------------------------------------------- -subroutine mesh_spectral_build_elements() - integer :: & - e, & - elemOffset +subroutine mesh_spectral_build_elements + integer :: & + e, & + elemOffset - allocate(mesh_element (4_pInt+8_pInt,theMesh%nElems), source = 0_pInt) + allocate(mesh_element (4+8,theMesh%nElems), source = 0) - elemOffset = product(grid(1:2))*grid3Offset - do e=1, theMesh%nElems - mesh_element( 1,e) = -1_pInt ! DEPRECATED - mesh_element( 2,e) = -1_pInt ! DEPRECATED - mesh_element( 3,e) = mesh_homogenizationAt(e) - mesh_element( 4,e) = microGlobal(e+elemOffset) ! microstructure - mesh_element( 5,e) = e + (e-1_pInt)/grid(1) + & - ((e-1_pInt)/(grid(1)*grid(2)))*(grid(1)+1_pInt) ! base node - mesh_element( 6,e) = mesh_element(5,e) + 1_pInt - mesh_element( 7,e) = mesh_element(5,e) + grid(1) + 2_pInt - mesh_element( 8,e) = mesh_element(5,e) + grid(1) + 1_pInt - mesh_element( 9,e) = mesh_element(5,e) +(grid(1) + 1_pInt) * (grid(2) + 1_pInt) ! second floor base node - mesh_element(10,e) = mesh_element(9,e) + 1_pInt - mesh_element(11,e) = mesh_element(9,e) + grid(1) + 2_pInt - mesh_element(12,e) = mesh_element(9,e) + grid(1) + 1_pInt - enddo + elemOffset = product(grid(1:2))*grid3Offset + do e=1, theMesh%nElems + mesh_element( 1,e) = -1 ! DEPRECATED + mesh_element( 2,e) = -1 ! DEPRECATED + mesh_element( 3,e) = mesh_homogenizationAt(e) + mesh_element( 4,e) = microGlobal(e+elemOffset) ! microstructure + mesh_element( 5,e) = e + (e-1)/grid(1) + & + ((e-1)/(grid(1)*grid(2)))*(grid(1)+1) ! base node + mesh_element( 6,e) = mesh_element(5,e) + 1 + mesh_element( 7,e) = mesh_element(5,e) + grid(1) + 2 + mesh_element( 8,e) = mesh_element(5,e) + grid(1) + 1 + mesh_element( 9,e) = mesh_element(5,e) +(grid(1) + 1) * (grid(2) + 1) ! second floor base node + mesh_element(10,e) = mesh_element(9,e) + 1 + mesh_element(11,e) = mesh_element(9,e) + grid(1) + 2 + mesh_element(12,e) = mesh_element(9,e) + grid(1) + 1 + enddo end subroutine mesh_spectral_build_elements @@ -447,50 +428,50 @@ end subroutine mesh_spectral_build_elements !-------------------------------------------------------------------------------------------------- subroutine mesh_spectral_build_ipNeighborhood - integer :: & - x,y,z, & - e - allocate(mesh_ipNeighborhood(3,6,1,theMesh%nElems),source=0_pInt) + integer :: & + x,y,z, & + e + allocate(mesh_ipNeighborhood(3,6,1,theMesh%nElems),source=0) - e = 0_pInt - do z = 0_pInt,grid3-1_pInt - do y = 0_pInt,grid(2)-1_pInt - do x = 0_pInt,grid(1)-1_pInt - e = e + 1_pInt - mesh_ipNeighborhood(1,1,1,e) = z * grid(1) * grid(2) & - + y * grid(1) & - + modulo(x+1_pInt,grid(1)) & - + 1_pInt - mesh_ipNeighborhood(1,2,1,e) = z * grid(1) * grid(2) & - + y * grid(1) & - + modulo(x-1_pInt,grid(1)) & - + 1_pInt - mesh_ipNeighborhood(1,3,1,e) = z * grid(1) * grid(2) & - + modulo(y+1_pInt,grid(2)) * grid(1) & - + x & - + 1_pInt - mesh_ipNeighborhood(1,4,1,e) = z * grid(1) * grid(2) & - + modulo(y-1_pInt,grid(2)) * grid(1) & - + x & - + 1_pInt - mesh_ipNeighborhood(1,5,1,e) = modulo(z+1_pInt,grid3) * grid(1) * grid(2) & - + y * grid(1) & - + x & - + 1_pInt - mesh_ipNeighborhood(1,6,1,e) = modulo(z-1_pInt,grid3) * grid(1) * grid(2) & - + y * grid(1) & - + x & - + 1_pInt - mesh_ipNeighborhood(2,1:6,1,e) = 1_pInt - mesh_ipNeighborhood(3,1,1,e) = 2_pInt - mesh_ipNeighborhood(3,2,1,e) = 1_pInt - mesh_ipNeighborhood(3,3,1,e) = 4_pInt - mesh_ipNeighborhood(3,4,1,e) = 3_pInt - mesh_ipNeighborhood(3,5,1,e) = 6_pInt - mesh_ipNeighborhood(3,6,1,e) = 5_pInt - enddo - enddo - enddo + e = 0 + do z = 0,grid3-1 + do y = 0,grid(2)-1 + do x = 0,grid(1)-1 + e = e + 1 + mesh_ipNeighborhood(1,1,1,e) = z * grid(1) * grid(2) & + + y * grid(1) & + + modulo(x+1,grid(1)) & + + 1 + mesh_ipNeighborhood(1,2,1,e) = z * grid(1) * grid(2) & + + y * grid(1) & + + modulo(x-1,grid(1)) & + + 1 + mesh_ipNeighborhood(1,3,1,e) = z * grid(1) * grid(2) & + + modulo(y+1,grid(2)) * grid(1) & + + x & + + 1 + mesh_ipNeighborhood(1,4,1,e) = z * grid(1) * grid(2) & + + modulo(y-1,grid(2)) * grid(1) & + + x & + + 1 + mesh_ipNeighborhood(1,5,1,e) = modulo(z+1,grid3) * grid(1) * grid(2) & + + y * grid(1) & + + x & + + 1 + mesh_ipNeighborhood(1,6,1,e) = modulo(z-1,grid3) * grid(1) * grid(2) & + + y * grid(1) & + + x & + + 1 + mesh_ipNeighborhood(2,1:6,1,e) = 1 + mesh_ipNeighborhood(3,1,1,e) = 2 + mesh_ipNeighborhood(3,2,1,e) = 1 + mesh_ipNeighborhood(3,3,1,e) = 4 + mesh_ipNeighborhood(3,4,1,e) = 3 + mesh_ipNeighborhood(3,5,1,e) = 6 + mesh_ipNeighborhood(3,6,1,e) = 5 + enddo + enddo + enddo end subroutine mesh_spectral_build_ipNeighborhood @@ -499,41 +480,37 @@ end subroutine mesh_spectral_build_ipNeighborhood !> @brief builds mesh of (distorted) cubes for given coordinates (= center of the cubes) !-------------------------------------------------------------------------------------------------- function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes) - use debug, only: & - debug_mesh, & - debug_level, & - debug_levelBasic - - real(pReal), intent(in), dimension(:,:,:,:) :: & - centres - real(pReal), dimension(3,size(centres,2)+1,size(centres,3)+1,size(centres,4)+1) :: & - nodes - real(pReal), intent(in), dimension(3) :: & - gDim - real(pReal), intent(in), dimension(3,3) :: & - Favg - real(pReal), dimension(3,size(centres,2)+2,size(centres,3)+2,size(centres,4)+2) :: & - wrappedCentres - - integer :: & - i,j,k,n - integer, dimension(3), parameter :: & - diag = 1_pInt - integer, dimension(3) :: & - shift = 0_pInt, & - lookup = 0_pInt, & - me = 0_pInt, & - iRes = 0_pInt - integer, dimension(3,8) :: & - neighbor = reshape([ & - 0_pInt, 0_pInt, 0_pInt, & - 1_pInt, 0_pInt, 0_pInt, & - 1_pInt, 1_pInt, 0_pInt, & - 0_pInt, 1_pInt, 0_pInt, & - 0_pInt, 0_pInt, 1_pInt, & - 1_pInt, 0_pInt, 1_pInt, & - 1_pInt, 1_pInt, 1_pInt, & - 0_pInt, 1_pInt, 1_pInt ], [3,8]) + + real(pReal), intent(in), dimension(:,:,:,:) :: & + centres + real(pReal), dimension(3,size(centres,2)+1,size(centres,3)+1,size(centres,4)+1) :: & + nodes + real(pReal), intent(in), dimension(3) :: & + gDim + real(pReal), intent(in), dimension(3,3) :: & + Favg + real(pReal), dimension(3,size(centres,2)+2,size(centres,3)+2,size(centres,4)+2) :: & + wrappedCentres + + integer :: & + i,j,k,n + integer, dimension(3), parameter :: & + diag = 1 + integer, dimension(3) :: & + shift = 0, & + lookup = 0, & + me = 0, & + iRes = 0 + integer, dimension(3,8) :: & + neighbor = reshape([ & + 0, 0, 0, & + 1, 0, 0, & + 1, 1, 0, & + 0, 1, 0, & + 0, 0, 1, & + 1, 0, 1, & + 1, 1, 1, & + 0, 1, 1 ], [3,8]) !-------------------------------------------------------------------------------------------------- ! initializing variables @@ -541,43 +518,35 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes) nodes = 0.0_pReal wrappedCentres = 0.0_pReal -!-------------------------------------------------------------------------------------------------- -! report - if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then - write(6,'(a)') ' Meshing cubes around centroids' - write(6,'(a,3(e12.5))') ' Dimension: ', gDim - write(6,'(a,3(i5))') ' Resolution:', iRes - endif - !-------------------------------------------------------------------------------------------------- ! building wrappedCentres = centroids + ghosts - wrappedCentres(1:3,2_pInt:iRes(1)+1_pInt,2_pInt:iRes(2)+1_pInt,2_pInt:iRes(3)+1_pInt) = centres - do k = 0_pInt,iRes(3)+1_pInt - do j = 0_pInt,iRes(2)+1_pInt - do i = 0_pInt,iRes(1)+1_pInt - if (k==0_pInt .or. k==iRes(3)+1_pInt .or. & ! z skin - j==0_pInt .or. j==iRes(2)+1_pInt .or. & ! y skin - i==0_pInt .or. i==iRes(1)+1_pInt ) then ! x skin - me = [i,j,k] ! me on skin - shift = sign(abs(iRes+diag-2_pInt*me)/(iRes+diag),iRes+diag-2_pInt*me) - lookup = me-diag+shift*iRes - wrappedCentres(1:3,i+1_pInt, j+1_pInt, k+1_pInt) = & - centres(1:3,lookup(1)+1_pInt,lookup(2)+1_pInt,lookup(3)+1_pInt) & - - matmul(Favg, real(shift,pReal)*gDim) - endif - enddo; enddo; enddo + wrappedCentres(1:3,2:iRes(1)+1,2:iRes(2)+1,2:iRes(3)+1) = centres + do k = 0,iRes(3)+1 + do j = 0,iRes(2)+1 + do i = 0,iRes(1)+1 + if (k==0 .or. k==iRes(3)+1 .or. & ! z skin + j==0 .or. j==iRes(2)+1 .or. & ! y skin + i==0 .or. i==iRes(1)+1 ) then ! x skin + me = [i,j,k] ! me on skin + shift = sign(abs(iRes+diag-2*me)/(iRes+diag),iRes+diag-2*me) + lookup = me-diag+shift*iRes + wrappedCentres(1:3,i+1, j+1, k+1) = & + centres(1:3,lookup(1)+1,lookup(2)+1,lookup(3)+1) & + - matmul(Favg, real(shift,pReal)*gDim) + endif + enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! averaging - do k = 0_pInt,iRes(3); do j = 0_pInt,iRes(2); do i = 0_pInt,iRes(1) - do n = 1_pInt,8_pInt - nodes(1:3,i+1_pInt,j+1_pInt,k+1_pInt) = & - nodes(1:3,i+1_pInt,j+1_pInt,k+1_pInt) + wrappedCentres(1:3,i+1_pInt+neighbor(1,n), & - j+1_pInt+neighbor(2,n), & - k+1_pInt+neighbor(3,n) ) - enddo - enddo; enddo; enddo - nodes = nodes/8.0_pReal + do k = 0,iRes(3); do j = 0,iRes(2); do i = 0,iRes(1) + do n = 1,8 + nodes(1:3,i+1,j+1,k+1) = & + nodes(1:3,i+1,j+1,k+1) + wrappedCentres(1:3,i+1+neighbor(1,n), & + j+1+neighbor(2,n), & + k+1+neighbor(3,n) ) + enddo + enddo; enddo; enddo + nodes = nodes/8.0_pReal end function mesh_nodesAroundCentres From 9612dd1d829a1d2be37c91ec979e7655d472b98e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 4 Jun 2019 22:38:01 +0200 Subject: [PATCH 08/44] simple data structure needs explanation --- src/mesh_grid.f90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index 17eb677d2..fce7d4b82 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -438,6 +438,7 @@ subroutine mesh_spectral_build_ipNeighborhood do y = 0,grid(2)-1 do x = 0,grid(1)-1 e = e + 1 + ! neigboring element mesh_ipNeighborhood(1,1,1,e) = z * grid(1) * grid(2) & + y * grid(1) & + modulo(x+1,grid(1)) & @@ -462,7 +463,9 @@ subroutine mesh_spectral_build_ipNeighborhood + y * grid(1) & + x & + 1 + ! neigboring IP mesh_ipNeighborhood(2,1:6,1,e) = 1 + ! neigboring face mesh_ipNeighborhood(3,1,1,e) = 2 mesh_ipNeighborhood(3,2,1,e) = 1 mesh_ipNeighborhood(3,3,1,e) = 4 From b6bb8a465009b414a00927a3a5b4a5f6e51485e9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 6 Jul 2019 21:41:29 -0700 Subject: [PATCH 09/44] simple strain calculation for DADF5 further enhancement requires to give optional arguments to add_genericpontwise --- processing/post/addStrainTensors.py | 6 +++--- python/damask/dadf5.py | 29 ++++++++++++++++++++++++++++- 2 files changed, 31 insertions(+), 4 deletions(-) diff --git a/processing/post/addStrainTensors.py b/processing/post/addStrainTensors.py index 2d62c31ae..7addc450a 100755 --- a/processing/post/addStrainTensors.py +++ b/processing/post/addStrainTensors.py @@ -136,9 +136,9 @@ for name in filenames: for column in items['tensor']['column']: # loop over all requested defgrads F = np.array(list(map(float,table.data[column:column+items['tensor']['dim']])),'d').reshape(items['tensor']['shape']) (U,S,Vh) = np.linalg.svd(F) # singular value decomposition - R = np.dot(U,Vh) # rotation of polar decomposition - stretch['U'] = np.dot(np.linalg.inv(R),F) # F = RU - stretch['V'] = np.dot(F,np.linalg.inv(R)) # F = VR + R_inv = np.linalg.inv(np.dot(U,Vh)) # inverse rotation of polar decomposition + stretch['U'] = np.dot(R_inv,F) # F = RU + stretch['V'] = np.dot(F,R_inv) # F = VR for theStretch in stretches: stretch[theStretch] = np.where(abs(stretch[theStretch]) < 1e-12, 0, stretch[theStretch]) # kill nasty noisy data diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 716747288..a96f4cfdb 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -224,7 +224,7 @@ class DADF5(): def add_determinant(self,a): - """Adds the determinan of a tensor""" + """Adds the determinant of a tensor""" # ToDo: The output unit should be the input unit args = [{'label':a,'shape':[3,3],'unit':None}] result = {'label':'det({})'.format(a), @@ -233,6 +233,33 @@ class DADF5(): self.add_generic_pointwise(np.linalg.det,args,result) + + def add_strain_tensors(self,defgrad='F'): + """Adds a strain definition""" + def strain(defgrad): + (U,S,Vh) = np.linalg.svd(defgrad) # singular value decomposition + R_inv = np.linalg.inv(np.dot(U,Vh)) # inverse rotation of polar decomposition + U = np.dot(R_inv,defgrad) # F = RU + U = np.where(abs(U) < 1e-12, 0, U) # kill nasty noisy data + (D,V) = np.linalg.eig(U) # eigen decomposition (of symmetric matrix) + neg = np.where(D < 0.0) # find negative eigenvalues ... + D[neg] *= -1. # ... flip value ... + V[:,neg] *= -1. # ... and vector + for i,eigval in enumerate(D): + if np.dot(V[:,i],V[:,(i+1)%3]) != 0.0: # check each vector for orthogonality + V[:,(i+1)%3] = np.cross(V[:,(i+2)%3],V[:,i]) # correct next vector + V[:,(i+1)%3] /= np.sqrt(np.dot(V[:,(i+1)%3],V[:,(i+1)%3].conj())) # and renormalize (hyperphobic?) + d = np.log(D) # operate on eigenvalues of U o r V + return np.dot(V,np.dot(np.diag(d),V.T)).real # build tensor back from eigenvalue/vector basis + + # ToDo: The output unit should be the input unit + args = [{'label':defgrad,'shape':[3,3],'unit':None}] + result = {'label':'strain({})'.format(defgrad), + 'unit':'-', + 'Description': 'strain (ln(V)) of a deformation gradient'} + + self.add_generic_pointwise(strain,args,result) + def get_fitting(self,data): groups = [] From fa6e88970cf93f219dc939f51d209d0b37895b36 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 15 Jul 2019 13:53:34 -0700 Subject: [PATCH 10/44] avoid empty entries --- src/config.f90 | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/config.f90 b/src/config.f90 index 85efcb82f..90233057c 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -97,11 +97,16 @@ subroutine config_init enddo - if (size(config_homogenization) < 1) call IO_error(160,ext_msg='') - if (size(config_microstructure) < 1) call IO_error(160,ext_msg='') - if (size(config_crystallite) < 1) call IO_error(160,ext_msg='') - if (size(config_phase) < 1) call IO_error(160,ext_msg='') - if (size(config_texture) < 1) call IO_error(160,ext_msg='') + if (.not. allocated(config_homogenization) .or. size(config_homogenization) < 1) & + call IO_error(160,ext_msg='') + if (.not. allocated(config_microstructure) .or. size(config_microstructure) < 1) & + call IO_error(160,ext_msg='') + if (.not. allocated(config_crystallite) .or. size(config_crystallite) < 1) & + call IO_error(160,ext_msg='') + if (.not. allocated(config_phase) .or. size(config_phase) < 1) & + call IO_error(160,ext_msg='') + if (.not. allocated(config_texture) .or. size(config_texture) < 1) & + call IO_error(160,ext_msg='') inquire(file='numerics.config', exist=fileExists) From c51cf8d506f4767e71b3f8373324acaebfc6da89 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 15 Jul 2019 13:55:14 -0700 Subject: [PATCH 11/44] transferring post processing capabilities --- python/damask/dadf5.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index a96f4cfdb..fcfa9ccd1 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -232,6 +232,35 @@ class DADF5(): 'Description': 'Determinant of a tensor'} self.add_generic_pointwise(np.linalg.det,args,result) + + + def add_spherical(self,a): + """Adds the spherical component of a tensor""" + def spherical(m): + return (m[0,0]+m[1,1]+m[2,2])/3.0 + + # ToDo: The output unit should be the input unit + args = [{'label':a,'shape':[3,3],'unit':None}] + result = {'label':'sph({})'.format(a), + 'unit':'n/a', + 'Description': 'Spherical component of a tensor'} + + self.add_generic_pointwise(spherical,args,result) + + + def add_deviator(self,a): + """Adds the deviator of a tensor""" + def deviator(m): + return m - np.eye(3)*(m[0,0]+m[1,1]+m[2,2])/3.0 + + # ToDo: The output unit should be the input unit + args = [{'label':a,'shape':[3,3],'unit':'Pa'}] + result = {'label':'dev({})'.format(a), + 'unit':'n/a', + 'Description': 'Deviatoric component of a tensor'} + + self.add_generic_pointwise(deviator,args,result) + def add_strain_tensors(self,defgrad='F'): From d446248d75404413c60022ebc9b8ea84c816b103 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 15 Jul 2019 17:07:04 -0700 Subject: [PATCH 12/44] also ouptut materialpoint results if requested --- processing/post/DADF5_vtk_cells.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py index ef27a973c..75301386c 100755 --- a/processing/post/DADF5_vtk_cells.py +++ b/processing/post/DADF5_vtk_cells.py @@ -60,6 +60,7 @@ for filename in options.filenames: vtk_data = [] results.active['increments'] = [inc] + results.active['m_output_types'] = [] for label in options.con: for o in results.c_output_types: @@ -85,6 +86,32 @@ for filename in options.filenames: vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) vtk_data[-1].SetName('1_'+x[0].split('/')[1]+'/generic/'+label) rGrid.GetCellData().AddArray(vtk_data[-1]) + + results.active['c_output_types'] = [] + for label in options.mat: + for o in results.m_output_types: + results.active['m_output_types'] = [o] + if o != 'generic': + for m in results.materialpoints: + results.active['materialpoints'] = [m] + x = results.get_dataset_location(label) + if len(x) == 0: + continue + array = results.read_dataset(x,0) + shape = [array.shape[0],np.product(array.shape[1:])] + vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) + vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) + rGrid.GetCellData().AddArray(vtk_data[-1]) + else: + results.active['materialpoints'] = results.materialpoints + x = results.get_dataset_location(label) + if len(x) == 0: + continue + array = results.read_dataset(x,0) + shape = [array.shape[0],np.product(array.shape[1:])] + vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) + vtk_data[-1].SetName('1_'+x[0].split('/')[1]+'/generic/'+label) + rGrid.GetCellData().AddArray(vtk_data[-1]) if results.structured: writer = vtk.vtkXMLRectilinearGridWriter() From 7ccc097406ed59532821d738b66fc0bdc4eefa15 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 15 Jul 2019 17:08:18 -0700 Subject: [PATCH 13/44] for testing output of materialpoint results --- src/homogenization.f90 | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 2ee921d10..64019b880 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -800,6 +800,8 @@ subroutine homogenization_results integer :: p character(len=256) :: group + + real(pReal), dimension(:,:,:), allocatable :: temp do p=1,size(config_name_homogenization) group = trim('current/materialpoint')//'/'//trim(config_name_homogenization(p)) @@ -812,7 +814,17 @@ subroutine homogenization_results case(HOMOGENIZATION_rgc_ID) call mech_RGC_results(homogenization_typeInstance(p),group) end select - + + group = trim('current/materialpoint')//'/'//trim(config_name_homogenization(p))//'/generic' + call HDF5_closeGroup(results_addGroup(group)) + + !temp = reshape(materialpoint_F,[3,3,discretization_nIP*discretization_nElem]) + !call results_writeDataset(group,temp,'F',& + ! 'deformation gradient','1') + !temp = reshape(materialpoint_P,[3,3,discretization_nIP*discretization_nElem]) + !call results_writeDataset(group,temp,'P',& + ! '1st Piola-Kirchoff stress','Pa') + enddo #endif end subroutine homogenization_results From 2fc66cff5b90f1c81158f424f926d57ce619b03b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 4 Sep 2019 14:30:26 -0700 Subject: [PATCH 14/44] better readable --- src/lattice.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 126d76d4d..7a47aa055 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -747,13 +747,13 @@ pure function lattice_symmetrizeC66(struct,C66) case (LATTICE_iso_ID) do k=1,3 forall(j=1:3) lattice_symmetrizeC66(k,j) = C66(1,2) - lattice_symmetrizeC66(k,k) = C66(1,1) + lattice_symmetrizeC66(k,k) = C66(1,1) lattice_symmetrizeC66(k+3,k+3) = 0.5_pReal*(C66(1,1)-C66(1,2)) enddo case (LATTICE_fcc_ID,LATTICE_bcc_ID) do k=1,3 - forall(j=1:3) lattice_symmetrizeC66(k,j) = C66(1,2) - lattice_symmetrizeC66(k,k) = C66(1,1) + forall(j=1:3) lattice_symmetrizeC66(k,j) = C66(1,2) + lattice_symmetrizeC66(k,k) = C66(1,1) lattice_symmetrizeC66(k+3,k+3) = C66(4,4) enddo case (LATTICE_hex_ID) From bc893762eca13b0a0928c944a476839ecd0bfc24 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 11 Sep 2019 07:29:34 -0700 Subject: [PATCH 15/44] no need to inverse a rotation, transpose is faster --- python/damask/dadf5.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index fcfa9ccd1..2ec1da286 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -267,7 +267,7 @@ class DADF5(): """Adds a strain definition""" def strain(defgrad): (U,S,Vh) = np.linalg.svd(defgrad) # singular value decomposition - R_inv = np.linalg.inv(np.dot(U,Vh)) # inverse rotation of polar decomposition + R_inv = np.dot(U,Vh).T # inverse rotation of polar decomposition U = np.dot(R_inv,defgrad) # F = RU U = np.where(abs(U) < 1e-12, 0, U) # kill nasty noisy data (D,V) = np.linalg.eig(U) # eigen decomposition (of symmetric matrix) From a6c69a744b29a748ef9756fbeafaaa27c91a45dc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 11 Sep 2019 17:57:24 -0700 Subject: [PATCH 16/44] do operations vectorized --- python/damask/dadf5.py | 89 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 77 insertions(+), 12 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 2ec1da286..fd3e24b6c 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -184,9 +184,15 @@ class DADF5(): def add_Cauchy(self,P='P',F='F'): - """Adds Cauchy stress calculated from 1st Piola-Kirchhoff stress and deformation gradient""" + """ + Adds Cauchy stress calculated from 1st Piola-Kirchhoff stress and deformation gradient. + + Todo + ---- + The einsum formula is completely untested! + """ def Cauchy(F,P): - return 1.0/np.linalg.det(F)*np.dot(P,F.T) + return np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F),F,P) args = [{'label':F,'shape':[3,3],'unit':'-'}, {'label':P,'shape':[3,3],'unit':'Pa'} ] @@ -194,7 +200,7 @@ class DADF5(): 'unit':'Pa', 'Description': 'Cauchy stress calculated from 1st Piola-Kirchhoff stress and deformation gradient'} - self.add_generic_pointwise(Cauchy,args,result) + self.add_generic_pointwise_vectorized(Cauchy,args,result) def add_Mises_stress(self,stress='sigma'): @@ -220,7 +226,7 @@ class DADF5(): 'unit':'n/a', 'Description': 'Norm of vector or tensor or magnitude of a scalar. See numpy.linalg.norm manual for details'} - self.add_generic_pointwise(np.linalg.norm,args,result) + self.add_generic_pointwise_vectorized(np.linalg.norm,args,result) def add_determinant(self,a): @@ -231,7 +237,7 @@ class DADF5(): 'unit':'n/a', 'Description': 'Determinant of a tensor'} - self.add_generic_pointwise(np.linalg.det,args,result) + self.add_generic_pointwise_vectorized(np.linalg.det,args,result) def add_spherical(self,a): @@ -307,9 +313,10 @@ class DADF5(): def add_generic_pointwise(self,func,args,result): """ - Ggeneral function to add pointwise data + Ggeneral function to add pointwise data. function 'func' first needs to have data arguments before other arguments + Works for functions that are pointwise defined. """ groups = self.get_fitting(args) @@ -317,9 +324,6 @@ class DADF5(): out = args['out'] datasets_in = args['dat'] func = args['fun'] - # calling the function per point might be performance-wise not optimal - # could be worth to investigate the performance for vectorized add_XXX functions that do the - # loops internally for i in range(out.shape[0]): arg = tuple([d[i,] for d in datasets_in]) out[i,] = func(*arg) @@ -347,7 +351,7 @@ class DADF5(): # Add the jobs in bulk to the thread pool. Alternatively you could use # `pool.add_task` to add single jobs. The code will block here, which # makes it possible to cancel the thread pool with an exception when - # the currently running batch of workers is finishnumpy.linalg.normed. + # the currently running batch of workers is finished pool.map(job, todo[:Nthreads+1]) i = 0 @@ -367,6 +371,67 @@ class DADF5(): i+=1 pool.wait_completion() - + + + def add_generic_pointwise_vectorized(self,func,args,result): + """ + Ggeneral function to add pointwise data + + function 'func' first needs to have data arguments before other arguments + Works for vectorized functions. + """ + groups = self.get_fitting(args) + + def job(args): + out = args['out'] + datasets_in = args['dat'] + func = args['fun'] + out = func(*datasets_in) + args['results'].put({'out':out,'group':args['group']}) + + Nthreads = 4 # ToDo: should be a parameter + results = Queue(Nthreads+1) + + todo = [] + + for g in groups: + with h5py.File(self.filename,'r') as f: + datasets_in = [f[g+'/'+u['label']][()] for u in args] - + # figure out dimension of results + for d in datasets_in: + print('input shape',d.shape) + testArg = tuple([d[0:1,] for d in datasets_in]) # to call function with first point + #print('testArg',testArg) + out = np.empty([datasets_in[0].shape[0]] + list(func(*testArg).shape[1:])) # shape is Npoints x shape of the results for one point + print('estimated output shape',out.shape) + todo.append({'dat':datasets_in,'fun':func,'out':out,'group':g,'results':results}) + + # Instantiate a thread pool with worker threads + pool = util.ThreadPool(Nthreads) + missingResults = len(todo) + + + # Add the jobs in bulk to the thread pool. Alternatively you could use + # `pool.add_task` to add single jobs. The code will block here, which + # makes it possible to cancel the thread pool with an exception when + # the currently running batch of workers is finished + + pool.map(job, todo[:Nthreads+1]) + i = 0 + while missingResults > 0: + r=results.get() # noqa + print(r['group']) + with h5py.File(self.filename,'r+') as f: + dataset_out = f[r['group']].create_dataset(result['label'],data=r['out']) + dataset_out.attrs['Unit'] = result['unit'] + dataset_out.attrs['Description'] = result['Description'] + dataset_out.attrs['Creator'] = 'dadf5.py v{}'.format('n/a') + missingResults-=1 + try: + pool.add_task(job,todo[Nthreads+1+i]) + except: + pass + i+=1 + + pool.wait_completion() From 953ba5321163531ab259a9f1998bae93ad8c264c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 11 Sep 2019 18:03:19 -0700 Subject: [PATCH 17/44] adjusting for strict prospector checking --- python/damask/dadf5.py | 66 +++++++++++++++++++++++++----------------- 1 file changed, 39 insertions(+), 27 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index fd3e24b6c..771b35936 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -1,20 +1,35 @@ -# -*- coding: UTF-8 no BOM -*- -import h5py -import re -import numpy as np from queue import Queue +import re + +import h5py +import numpy as np + from . import util # ------------------------------------------------------------------ class DADF5(): - """Read and write to DADF5 files""" + """ + Read and write to DADF5 files. + + DADF5 files contain DAMASK results. + """ # ------------------------------------------------------------------ def __init__(self, filename, mode = 'r', ): + """ + Opens an existing DADF5 file. + Parameters + ---------- + filename : str + name of the DADF5 file to be openend. + mode : str, optional + filemode for opening, either 'r' or 'a'. + + """ if mode not in ['a','r']: print('Invalid file access mode') else: @@ -95,7 +110,7 @@ class DADF5(): return groups def list_data(self): - """Shows information on all datasets in the file""" + """Shows information on all datasets in the file.""" with h5py.File(self.filename,'r') as f: group_inc = 'inc{:05}'.format(self.active['increments'][0]['inc']) for c in self.active['constituents']: @@ -107,7 +122,7 @@ class DADF5(): try: for x in f[group_output_types].keys(): print(' {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode())) - except: + except KeyError: pass for m in self.active['materialpoints']: group_materialpoint = group_inc+'/materialpoint/'+m @@ -117,12 +132,12 @@ class DADF5(): try: for x in f[group_output_types].keys(): print(' {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode())) - except: + except KeyError: pass def get_dataset_location(self,label): - """Returns the location of all active datasets with given label""" + """Returns the location of all active datasets with given label.""" path = [] with h5py.File(self.filename,'r') as f: for i in self.active['increments']: @@ -134,7 +149,7 @@ class DADF5(): try: f[group_constituent+'/'+t+'/'+label] path.append(group_constituent+'/'+t+'/'+label) - except Exception as e: + except KeyError as e: print('unable to locate constituents dataset: '+ str(e)) for m in self.active['materialpoints']: @@ -143,7 +158,7 @@ class DADF5(): try: f[group_materialpoint+'/'+t+'/'+label] path.append(group_materialpoint+'/'+t+'/'+label) - except Exception as e: + except KeyError as e: print('unable to locate materialpoints dataset: '+ str(e)) return path @@ -151,7 +166,7 @@ class DADF5(): def read_dataset(self,path,c): """ - Dataset for all points/cells + Dataset for all points/cells. If more than one path is given, the dataset is composed of the individual contributions """ @@ -168,7 +183,7 @@ class DADF5(): if len(a.shape) == 1: a=a.reshape([a.shape[0],1]) dataset[p,:] = a[u,:] - except Exception as e: + except KeyError as e: print('unable to read constituent: '+ str(e)) try: p = np.where(f['mapping/cellResults/materialpoint']['Name'] == str.encode(label))[0] @@ -177,7 +192,7 @@ class DADF5(): if len(a.shape) == 1: a=a.reshape([a.shape[0],1]) dataset[p,:] = a[u,:] - except Exception as e: + except KeyError as e: print('unable to read materialpoint: '+ str(e)) return dataset @@ -190,6 +205,7 @@ class DADF5(): Todo ---- The einsum formula is completely untested! + """ def Cauchy(F,P): return np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F),F,P) @@ -204,7 +220,7 @@ class DADF5(): def add_Mises_stress(self,stress='sigma'): - """Adds equivalent von Mises stress""" + """Adds equivalent von Mises stress.""" def Mises_stress(stress): dev = stress - np.trace(stress)/3.0*np.eye(3) symdev = 0.5*(dev+dev.T) @@ -230,7 +246,7 @@ class DADF5(): def add_determinant(self,a): - """Adds the determinant of a tensor""" + """Adds the determinant of a tensor.""" # ToDo: The output unit should be the input unit args = [{'label':a,'shape':[3,3],'unit':None}] result = {'label':'det({})'.format(a), @@ -241,7 +257,7 @@ class DADF5(): def add_spherical(self,a): - """Adds the spherical component of a tensor""" + """Adds the spherical component of a tensor.""" def spherical(m): return (m[0,0]+m[1,1]+m[2,2])/3.0 @@ -255,7 +271,7 @@ class DADF5(): def add_deviator(self,a): - """Adds the deviator of a tensor""" + """Adds the deviator of a tensor.""" def deviator(m): return m - np.eye(3)*(m[0,0]+m[1,1]+m[2,2])/3.0 @@ -270,7 +286,7 @@ class DADF5(): def add_strain_tensors(self,defgrad='F'): - """Adds a strain definition""" + """Adds a strain definition.""" def strain(defgrad): (U,S,Vh) = np.linalg.svd(defgrad) # singular value decomposition R_inv = np.dot(U,Vh).T # inverse rotation of polar decomposition @@ -313,7 +329,7 @@ class DADF5(): def add_generic_pointwise(self,func,args,result): """ - Ggeneral function to add pointwise data. + General function to add pointwise data. function 'func' first needs to have data arguments before other arguments Works for functions that are pointwise defined. @@ -366,7 +382,7 @@ class DADF5(): missingResults-=1 try: pool.add_task(job,todo[Nthreads+1+i]) - except: + except IndexError: pass i+=1 @@ -375,7 +391,7 @@ class DADF5(): def add_generic_pointwise_vectorized(self,func,args,result): """ - Ggeneral function to add pointwise data + General function to add pointwise data. function 'func' first needs to have data arguments before other arguments Works for vectorized functions. @@ -399,12 +415,8 @@ class DADF5(): datasets_in = [f[g+'/'+u['label']][()] for u in args] # figure out dimension of results - for d in datasets_in: - print('input shape',d.shape) testArg = tuple([d[0:1,] for d in datasets_in]) # to call function with first point - #print('testArg',testArg) out = np.empty([datasets_in[0].shape[0]] + list(func(*testArg).shape[1:])) # shape is Npoints x shape of the results for one point - print('estimated output shape',out.shape) todo.append({'dat':datasets_in,'fun':func,'out':out,'group':g,'results':results}) # Instantiate a thread pool with worker threads @@ -430,7 +442,7 @@ class DADF5(): missingResults-=1 try: pool.add_task(job,todo[Nthreads+1+i]) - except: + except IndexError: pass i+=1 From de313269d95c8a382443535f63d8c81dc660fdcc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 11 Sep 2019 18:54:26 -0700 Subject: [PATCH 18/44] bugfix --- python/damask/dadf5.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 771b35936..3e919a77c 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -235,14 +235,22 @@ class DADF5(): def add_norm(self,x,ord=None): - """Adds norm of vector or tensor or magnitude of a scalar.""" - # ToDo: The output unit should be the input unit + """ + Adds norm of vector or tensor or magnitude of a scalar. + + Todo + ---- + The output unit should be the input unit. + The ord parameter should be taken into account. + The whole thing should be vectorized. This requires to parse optional arguments to func. + + """ args = [{'label':x,'shape':None,'unit':None}] result = {'label':'norm_{}({})'.format(str(ord),x), 'unit':'n/a', 'Description': 'Norm of vector or tensor or magnitude of a scalar. See numpy.linalg.norm manual for details'} - self.add_generic_pointwise_vectorized(np.linalg.norm,args,result) + self.add_generic_pointwise(np.linalg.norm,args,result) def add_determinant(self,a): From 3db3e9e7628d6d512c0031b1187ed2e1c907d7e1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 11 Sep 2019 21:20:14 -0700 Subject: [PATCH 19/44] preparing for use of optional arguments to function --- python/damask/dadf5.py | 48 ++++++++++++++++++++++++++++-------------- 1 file changed, 32 insertions(+), 16 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 3e919a77c..27a93097e 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -82,11 +82,21 @@ class DADF5(): self.filename = filename self.mode = mode - + + def get_candidates(self,l): + """ + Get groups that contain all requested datasets. + + Parameters + ---------- + l : list of str + Names of datasets that need to be located in the group. + + """ groups = [] if type(l) is not list: - print('mist') + raise TypeError('Candidates should be given as a list') with h5py.File(self.filename,'r') as f: for g in self.get_active_groups(): if set(l).issubset(f[g].keys()): groups.append(g) @@ -94,6 +104,9 @@ class DADF5(): def get_active_groups(self): + """ + Get groups that are currently considered for evaluation. + """ groups = [] for i,x in enumerate(self.active['increments']): group_inc = 'inc{:05}'.format(self.active['increments'][i]['inc']) @@ -107,10 +120,11 @@ class DADF5(): for t in self.active['m_output_types']: group_output_types = group_materialpoint+'/'+t groups.append(group_output_types) - return groups + return groups + def list_data(self): - """Shows information on all datasets in the file.""" + """Shows information on all active datasets in the file.""" with h5py.File(self.filename,'r') as f: group_inc = 'inc{:05}'.format(self.active['increments'][0]['inc']) for c in self.active['constituents']: @@ -216,7 +230,7 @@ class DADF5(): 'unit':'Pa', 'Description': 'Cauchy stress calculated from 1st Piola-Kirchhoff stress and deformation gradient'} - self.add_generic_pointwise_vectorized(Cauchy,args,result) + self.add_generic_pointwise_vectorized(Cauchy,args,None,result) def add_Mises_stress(self,stress='sigma'): @@ -397,7 +411,7 @@ class DADF5(): pool.wait_completion() - def add_generic_pointwise_vectorized(self,func,args,result): + def add_generic_pointwise_vectorized(self,func,args,args2=None,result=None): """ General function to add pointwise data. @@ -407,9 +421,16 @@ class DADF5(): groups = self.get_fitting(args) def job(args): + """ + A job. It has different args! + """ + print('args for job',args) out = args['out'] datasets_in = args['dat'] func = args['fun'] +# try: + # out = func(*datasets_in,*args['fun_args']) + # except: out = func(*datasets_in) args['results'].put({'out':out,'group':args['group']}) @@ -421,22 +442,17 @@ class DADF5(): for g in groups: with h5py.File(self.filename,'r') as f: datasets_in = [f[g+'/'+u['label']][()] for u in args] - - # figure out dimension of results - testArg = tuple([d[0:1,] for d in datasets_in]) # to call function with first point - out = np.empty([datasets_in[0].shape[0]] + list(func(*testArg).shape[1:])) # shape is Npoints x shape of the results for one point - todo.append({'dat':datasets_in,'fun':func,'out':out,'group':g,'results':results}) + + if args2 is not None: + todo.append({'dat':datasets_in,'fun':func,'group':g,'results':results,'func_args':args,'out':None}) + else: + todo.append({'dat':datasets_in,'fun':func,'group':g,'results':results,'out':None}) # Instantiate a thread pool with worker threads pool = util.ThreadPool(Nthreads) missingResults = len(todo) - # Add the jobs in bulk to the thread pool. Alternatively you could use - # `pool.add_task` to add single jobs. The code will block here, which - # makes it possible to cancel the thread pool with an exception when - # the currently running batch of workers is finished - pool.map(job, todo[:Nthreads+1]) i = 0 while missingResults > 0: From 6f008c5d5f157288079652a41308055913647656 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 13 Sep 2019 06:02:42 -0700 Subject: [PATCH 20/44] rewrite pointwise add function - all vectorized: Much faster - passing in all relevant information allows to do sanity checks and add useful meta data in HDF5 - improved readability --- python/damask/dadf5.py | 402 +++++++++++++++++++---------------------- 1 file changed, 187 insertions(+), 215 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 27a93097e..6edad9c37 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -84,7 +84,7 @@ class DADF5(): self.mode = mode - def get_candidates(self,l): + def get_groups(self,l): """ Get groups that contain all requested datasets. @@ -210,8 +210,8 @@ class DADF5(): print('unable to read materialpoint: '+ str(e)) return dataset - - + + def add_Cauchy(self,P='P',F='F'): """ Adds Cauchy stress calculated from 1st Piola-Kirchhoff stress and deformation gradient. @@ -222,252 +222,224 @@ class DADF5(): """ def Cauchy(F,P): - return np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F),F,P) + return { + 'data' : np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F['data']),F['data'],P['data']), + 'label' : 'sigma', + 'meta' : { + 'Unit' : P['meta']['Unit'], + 'Description' : 'Cauchy stress calculated from {} ({}) '.format(P['label'],P['meta']['Description'])+\ + 'and deformation gradient {} ({})'.format(F['label'],P['meta']['Description']), + 'Creator' : 'dadf5.py:add_Cauchy vXXXXX' + } + } - args = [{'label':F,'shape':[3,3],'unit':'-'}, - {'label':P,'shape':[3,3],'unit':'Pa'} ] - result = {'label':'sigma', - 'unit':'Pa', - 'Description': 'Cauchy stress calculated from 1st Piola-Kirchhoff stress and deformation gradient'} + requested = [{'label':F,'arg':'F'}, + {'label':P,'arg':'P'} ] - self.add_generic_pointwise_vectorized(Cauchy,args,None,result) + self.__add_generic_pointwise(Cauchy,requested) - - def add_Mises_stress(self,stress='sigma'): - """Adds equivalent von Mises stress.""" - def Mises_stress(stress): - dev = stress - np.trace(stress)/3.0*np.eye(3) - symdev = 0.5*(dev+dev.T) - return np.sqrt(np.sum(symdev*symdev.T)*3.0/2.0) + + def add_Mises(self,x): + """Adds the equivalent Mises stres of a tensor.""" + def deviator(x): - args = [{'label':stress,'shape':[3,3],'unit':'Pa'}] - result = {'label':'Mises({})'.format(stress), - 'unit':'Pa', - 'Description': 'Equivalent Mises stress'} - - self.add_generic_pointwise(Mises_stress,args,result) - + if x['meta']['Unit'] == 'Pa': + factor = 3.0/2.0 + elif x['meta']['Unit'] == '-': + factor = 2.0/3.0 + else: + ValueError + + d = x['data'] + dev = d - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[d.shape[0],3,3]),np.trace(d,axis1=1,axis2=2)/3.0) + dev_sym = (dev + np.einsum('ikj',dev))*0.5 + + return { + 'data' : np.sqrt(np.einsum('ijk->i',dev_sym**2)*factor), + 'label' : 'dev({})'.format(x['label']), + 'meta' : { + 'Unit' : x['meta']['Unit'], + 'Description' : 'Mises equivalent stress of {} ({})'.format(x['label'],x['meta']['Description']), + 'Creator' : 'dadf5.py:add_Mises_stress vXXXXX' + } + } + + requested = [{'label':x,'arg':'x'}] + self.__add_generic_pointwise(deviator,requested) + + def add_norm(self,x,ord=None): """ Adds norm of vector or tensor or magnitude of a scalar. - - Todo - ---- - The output unit should be the input unit. - The ord parameter should be taken into account. - The whole thing should be vectorized. This requires to parse optional arguments to func. - + See numpy.linalg.norm manual for details. """ - args = [{'label':x,'shape':None,'unit':None}] - result = {'label':'norm_{}({})'.format(str(ord),x), - 'unit':'n/a', - 'Description': 'Norm of vector or tensor or magnitude of a scalar. See numpy.linalg.norm manual for details'} + def norm(x,ord): + + if len(x['data'].shape) == 1: + axis = 0 + t = 'scalar' + elif len(x['data'].shape) == 2: + axis = 1 + t = 'vector' + elif len(x['data'].shape) == 3: + axis = (1,2) + t = 'tensor' + else: + raise ValueError + + return { + 'data' : np.linalg.norm(x['data'],ord=ord,axis=axis,keepdims=True), + 'label' : 'norm({})'.format(x['label']), + 'meta' : { + 'Unit' : x['meta']['Unit'], + 'Description' : 'Norm of {} {} ({})'.format(t,x['label'],x['meta']['Description']), + 'Creator' : 'dadf5.py:add_norm vXXXXX' + } + } + + requested = [{'label':x,'arg':'x'}] - self.add_generic_pointwise(np.linalg.norm,args,result) - - - def add_determinant(self,a): - """Adds the determinant of a tensor.""" - # ToDo: The output unit should be the input unit - args = [{'label':a,'shape':[3,3],'unit':None}] - result = {'label':'det({})'.format(a), - 'unit':'n/a', - 'Description': 'Determinant of a tensor'} - - self.add_generic_pointwise_vectorized(np.linalg.det,args,result) + self.__add_generic_pointwise(norm,requested,{'ord':ord}) - def add_spherical(self,a): + def add_determinant(self,x): + """Adds the determinant component of a tensor.""" + def determinant(x): + + return { + 'data' : np.linalg.det(x['data']), + 'label' : 'det({})'.format(x['label']), + 'meta' : { + 'Unit' : x['meta']['Unit'], + 'Description' : 'Determinant of tensor {} ({})'.format(x['label'],x['meta']['Description']), + 'Creator' : 'dadf5.py:add_determinant vXXXXX' + } + } + + requested = [{'label':x,'arg':'x'}] + + self.__add_generic_pointwise(determinant,requested) + + + def add_spherical(self,x): """Adds the spherical component of a tensor.""" - def spherical(m): - return (m[0,0]+m[1,1]+m[2,2])/3.0 + def spherical(x): + + if not np.all(np.array(x['data'].shape[1:]) == np.array([3,3])): + raise ValueError - # ToDo: The output unit should be the input unit - args = [{'label':a,'shape':[3,3],'unit':None}] - result = {'label':'sph({})'.format(a), - 'unit':'n/a', - 'Description': 'Spherical component of a tensor'} + return { + 'data' : np.trace(x['data'],axis1=1,axis2=2)/3.0, + 'label' : 'sph({})'.format(x['label']), + 'meta' : { + 'Unit' : x['meta']['Unit'], + 'Description' : 'Spherical component of tensor {} ({})'.format(x['label'],x['meta']['Description']), + 'Creator' : 'dadf5.py:add_spherical vXXXXX' + } + } - self.add_generic_pointwise(spherical,args,result) - - - def add_deviator(self,a): + requested = [{'label':x,'arg':'x'}] + + self.__add_generic_pointwise(spherical,requested) + + + def add_deviator(self,x): """Adds the deviator of a tensor.""" - def deviator(m): - return m - np.eye(3)*(m[0,0]+m[1,1]+m[2,2])/3.0 - - # ToDo: The output unit should be the input unit - args = [{'label':a,'shape':[3,3],'unit':'Pa'}] - result = {'label':'dev({})'.format(a), - 'unit':'n/a', - 'Description': 'Deviatoric component of a tensor'} - - self.add_generic_pointwise(deviator,args,result) - + def deviator(x): + d = x['data'] + return { + 'data' : d - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[d.shape[0],3,3]),np.trace(d,axis1=1,axis2=2)/3.0), + 'label' : 'dev({})'.format(x['label']), + 'meta' : { + 'Unit' : x['meta']['Unit'], + 'Description' : 'Deviator of tensor {} ({})'.format(x['label'],x['meta']['Description']), + 'Creator' : 'dadf5.py:add_deviator vXXXXX' + } + } + + requested = [{'label':x,'arg':'x'}] - - def add_strain_tensors(self,defgrad='F'): - """Adds a strain definition.""" - def strain(defgrad): - (U,S,Vh) = np.linalg.svd(defgrad) # singular value decomposition - R_inv = np.dot(U,Vh).T # inverse rotation of polar decomposition - U = np.dot(R_inv,defgrad) # F = RU - U = np.where(abs(U) < 1e-12, 0, U) # kill nasty noisy data - (D,V) = np.linalg.eig(U) # eigen decomposition (of symmetric matrix) - neg = np.where(D < 0.0) # find negative eigenvalues ... - D[neg] *= -1. # ... flip value ... - V[:,neg] *= -1. # ... and vector - for i,eigval in enumerate(D): - if np.dot(V[:,i],V[:,(i+1)%3]) != 0.0: # check each vector for orthogonality - V[:,(i+1)%3] = np.cross(V[:,(i+2)%3],V[:,i]) # correct next vector - V[:,(i+1)%3] /= np.sqrt(np.dot(V[:,(i+1)%3],V[:,(i+1)%3].conj())) # and renormalize (hyperphobic?) - d = np.log(D) # operate on eigenvalues of U o r V - return np.dot(V,np.dot(np.diag(d),V.T)).real # build tensor back from eigenvalue/vector basis + self.__add_generic_pointwise(deviator,requested) - # ToDo: The output unit should be the input unit - args = [{'label':defgrad,'shape':[3,3],'unit':None}] - result = {'label':'strain({})'.format(defgrad), - 'unit':'-', - 'Description': 'strain (ln(V)) of a deformation gradient'} - - self.add_generic_pointwise(strain,args,result) - - def get_fitting(self,data): - groups = [] - if type(data) is not list: - print('mist') - with h5py.File(self.filename,'r') as f: - for g in self.get_candidates([l['label'] for l in data]): - print(g) - fits = True - for d in data: # ToDo: check for unit - if d['shape'] is not None: - fits = fits and np.all(np.array(f[g+'/'+d['label']].shape[1:]) == np.array(d['shape'])) - if fits: groups.append(g) - return groups + def add_strain_tensor(self,t,ord,defgrad='F'): + """Adds the a strain tensor.""" + def strain_tensor(defgrad,t,ord): + (U,S,Vh) = np.linalg.svd(defgrad['data']) # singular value decomposition + R_inv = np.einsum('ijk->ikj',np.matmul(U,Vh)) # inverse rotation of polar decomposition + U = np.matmul(R_inv,defgrad['data']) # F = RU + (D,V) = np.linalg.eigh((U+np.einsum('ikj',U))*.5) # eigen decomposition (of symmetric(ed) matrix) + neg = np.where(D < 0.0) # find negative eigenvalues ... + D[neg[0],neg[1]] = D[neg[0],neg[1]]* -1 # ... flip value ... + V[neg[0],:,neg[1]] = V[neg[0],:,neg[1]]* -1 # ... and vector + + d = np.log(D) + a = np.matmul(V,np.einsum('ij,ikj->ijk',d,V)) # this is wrong ... + for j in range(V.shape[0]): # but this is slow ... + a[j,:,:] = np.dot(V[j,:,:],np.dot(np.diag(d[j,:]),V[j,:,:].T)) + print(np.max(a)) + + return { + 'data' : a, + 'label' : 'lnV({})'.format(defgrad['label']), + 'meta' : { + 'Unit' : defgrad['meta']['Unit'], + 'Description' : 'Strain tensor {} ({})'.format(defgrad['label'],defgrad['meta']['Description']), + 'Creator' : 'dadf5.py:add_deviator vXXXXX' + } + } + + requested = [{'label':defgrad,'arg':'defgrad'}] - def add_generic_pointwise(self,func,args,result): + self.__add_generic_pointwise(strain_tensor,requested,{'t':t,'ord':ord}) + + + def __add_generic_pointwise(self,func,datasets_requested,extra_args={}): """ General function to add pointwise data. - - function 'func' first needs to have data arguments before other arguments - Works for functions that are pointwise defined. """ - groups = self.get_fitting(args) - - def job(args): - out = args['out'] - datasets_in = args['dat'] - func = args['fun'] - for i in range(out.shape[0]): - arg = tuple([d[i,] for d in datasets_in]) - out[i,] = func(*arg) - args['results'].put({'out':out,'group':args['group']}) - - Nthreads = 4 # ToDo: should be a parameter - results = Queue(Nthreads+1) - - todo = [] - - for g in groups: - with h5py.File(self.filename,'r') as f: - datasets_in = [f[g+'/'+u['label']][()] for u in args] - - # figure out dimension of results - testArg = tuple([d[0,] for d in datasets_in]) # to call function with first point - out = np.empty([datasets_in[0].shape[0]] + list(func(*testArg).shape)) # shape is Npoints x shape of the results for one point - todo.append({'dat':datasets_in,'fun':func,'out':out,'group':g,'results':results}) - - # Instantiate a thread pool with worker threads - pool = util.ThreadPool(Nthreads) - missingResults = len(todo) - - # Add the jobs in bulk to the thread pool. Alternatively you could use - # `pool.add_task` to add single jobs. The code will block here, which - # makes it possible to cancel the thread pool with an exception when - # the currently running batch of workers is finished - - pool.map(job, todo[:Nthreads+1]) - i = 0 - while missingResults > 0: - r=results.get() # noqa - print(r['group']) - with h5py.File(self.filename,'r+') as f: - dataset_out = f[r['group']].create_dataset(result['label'],data=r['out']) - dataset_out.attrs['Unit'] = result['unit'] - dataset_out.attrs['Description'] = result['Description'] - dataset_out.attrs['Creator'] = 'dadf5.py v{}'.format('n/a') - missingResults-=1 - try: - pool.add_task(job,todo[Nthreads+1+i]) - except IndexError: - pass - i+=1 - - pool.wait_completion() - - - def add_generic_pointwise_vectorized(self,func,args,args2=None,result=None): - """ - General function to add pointwise data. - - function 'func' first needs to have data arguments before other arguments - Works for vectorized functions. - """ - groups = self.get_fitting(args) - def job(args): """ - A job. It has different args! + Call function with input data + extra arguments, returns results + group. """ - print('args for job',args) - out = args['out'] - datasets_in = args['dat'] - func = args['fun'] -# try: - # out = func(*datasets_in,*args['fun_args']) - # except: - out = func(*datasets_in) - args['results'].put({'out':out,'group':args['group']}) + args['results'].put({**args['func'](**args['in']),'group':args['group']}) - Nthreads = 4 # ToDo: should be a parameter - results = Queue(Nthreads+1) + + N_threads = 1 # ToDo: should be a parameter + + results = Queue(N_threads) + pool = util.ThreadPool(N_threads) + N_added = N_threads + 1 todo = [] - - for g in groups: + # ToDo: It would be more memory efficient to read only from file when required, i.e. do to it in pool.add_task + for group in self.get_groups([d['label'] for d in datasets_requested]): with h5py.File(self.filename,'r') as f: - datasets_in = [f[g+'/'+u['label']][()] for u in args] + datasets_in = {} + for d in datasets_requested: + loc = f[group+'/'+d['label']] + data = loc[()] + meta = {k:loc.attrs[k] for k in loc.attrs.keys()} + datasets_in[d['arg']] = {'data': data, 'meta' : meta, 'label' : d['label']} - if args2 is not None: - todo.append({'dat':datasets_in,'fun':func,'group':g,'results':results,'func_args':args,'out':None}) - else: - todo.append({'dat':datasets_in,'fun':func,'group':g,'results':results,'out':None}) + todo.append({'in':{**datasets_in,**extra_args},'func':func,'group':group,'results':results}) - # Instantiate a thread pool with worker threads - pool = util.ThreadPool(Nthreads) - missingResults = len(todo) + pool.map(job, todo[:N_added]) # initialize - - pool.map(job, todo[:Nthreads+1]) - i = 0 - while missingResults > 0: - r=results.get() # noqa - print(r['group']) - with h5py.File(self.filename,'r+') as f: - dataset_out = f[r['group']].create_dataset(result['label'],data=r['out']) - dataset_out.attrs['Unit'] = result['unit'] - dataset_out.attrs['Description'] = result['Description'] - dataset_out.attrs['Creator'] = 'dadf5.py v{}'.format('n/a') - missingResults-=1 - try: - pool.add_task(job,todo[Nthreads+1+i]) - except IndexError: - pass - i+=1 + N_not_calculated = len(todo) + while N_not_calculated > 0: + result = results.get() + with h5py.File(self.filename,self.mode) as f: # write to file + dataset_out = f[result['group']].create_dataset(result['label'],data=result['data']) + for k in result['meta'].keys(): + dataset_out.attrs[k] = result['meta'][k] + N_not_calculated-=1 + + if N_added < len(todo): # add more jobs + pool.add_task(job,todo[N_added]) + N_added +=1 pool.wait_completion() From f83a167414fb8917a8708c85f13501763de21230 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 13 Sep 2019 09:48:39 -0700 Subject: [PATCH 21/44] file mode for high level funtion not useful --- python/damask/dadf5.py | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 6edad9c37..688cbba08 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -17,7 +17,7 @@ class DADF5(): # ------------------------------------------------------------------ def __init__(self, filename, - mode = 'r', + mode = 'a', ): """ Opens an existing DADF5 file. @@ -26,16 +26,9 @@ class DADF5(): ---------- filename : str name of the DADF5 file to be openend. - mode : str, optional - filemode for opening, either 'r' or 'a'. """ - if mode not in ['a','r']: - print('Invalid file access mode') - else: - with h5py.File(filename,mode): - pass - + with h5py.File(filename,'r') as f: if f.attrs['DADF5-major'] != 0 or f.attrs['DADF5-minor'] != 2: From b2b625af3eff14e0e77772c25315b1b525312fa6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 13 Sep 2019 15:17:46 -0700 Subject: [PATCH 22/44] notes from discussion with Philip --- python/damask/dadf5.py | 39 +++++++++++++++++++-------------------- 1 file changed, 19 insertions(+), 20 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 688cbba08..44238572f 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -45,15 +45,11 @@ class DADF5(): 'time': round(f[u].attrs['time/s'],12), } for u in f.keys() if r.match(u)] - self.constituents = np.unique(f['mapping/cellResults/constituent']['Name']).tolist() # ToDo: I am not to happy with the name - self.constituents = [c.decode() for c in self.constituents] - - self.materialpoints = np.unique(f['mapping/cellResults/materialpoint']['Name']).tolist() # ToDo: I am not to happy with the name - self.materialpoints = [m.decode() for m in self.materialpoints] - - self.Nconstituents = [i for i in range(np.shape(f['mapping/cellResults/constituent'])[1])] - self.Nmaterialpoints = np.shape(f['mapping/cellResults/constituent'])[0] - + self.Nmaterialpoints, self.Nconstituents = np.shape(f['mapping/cellResults/constituent']) + self.materialpoints = [m.decode() for m in np.unique(f['mapping/cellResults/materialpoint']['Name'])] + self.constituents = [c.decode() for c in np.unique(f['mapping/cellResults/constituent'] ['Name'])] + + self.c_output_types = [] for c in self.constituents: for o in f['inc{:05}/constituent/{}'.format(self.increments[0]['inc'],c)].keys(): @@ -65,19 +61,22 @@ class DADF5(): for o in f['inc{:05}/materialpoint/{}'.format(self.increments[0]['inc'],m)].keys(): self.m_output_types.append(o) self.m_output_types = list(set(self.m_output_types)) # make unique - - self.active= {'increments': self.increments, + + #self.on_air + self.active= {'increments': self.increments, # ToDo:simplify, activity only positions that translate into (no complex types) 'constituents': self.constituents, 'materialpoints': self.materialpoints, - 'constituent': self.Nconstituents, + 'constituent': range(self.Nconstituents), # ToDo: stupid naming 'c_output_types': self.c_output_types, 'm_output_types': self.m_output_types} +# ToDo: store increments, select icrements (trivial), position, and time + self.filename = filename self.mode = mode - def get_groups(self,l): + def get_groups(self,l): #group_with_data(datasets) """ Get groups that contain all requested datasets. @@ -96,13 +95,13 @@ class DADF5(): return groups - def get_active_groups(self): + def get_active_groups(self): # rename: get_groups needed? merge with datasets and have [] and ['*'] """ Get groups that are currently considered for evaluation. """ groups = [] for i,x in enumerate(self.active['increments']): - group_inc = 'inc{:05}'.format(self.active['increments'][i]['inc']) + group_inc = 'inc{:05}'.format(self.active['increments'][i]['inc']) #ToDo: Merge path only once at the end '/'.join(listE) for c in self.active['constituents']: group_constituent = group_inc+'/constituent/'+c for t in self.active['c_output_types']: @@ -116,16 +115,16 @@ class DADF5(): return groups - def list_data(self): + def list_data(self): # print_datasets and have [] and ['*'], loop over all increment, soll auf anderen basieren (get groups with sternchen) """Shows information on all active datasets in the file.""" with h5py.File(self.filename,'r') as f: - group_inc = 'inc{:05}'.format(self.active['increments'][0]['inc']) + group_inc = 'inc{:05}'.format(self.active['increments'][0]['inc']) #ToDo: Merge path only once at the end '/'.join(listE) for c in self.active['constituents']: print('\n'+c) group_constituent = group_inc+'/constituent/'+c for t in self.active['c_output_types']: print(' {}'.format(t)) - group_output_types = group_constituent+'/'+t + group_output_types = group_constituent+'/'+t try: for x in f[group_output_types].keys(): print(' {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode())) @@ -143,8 +142,8 @@ class DADF5(): pass - def get_dataset_location(self,label): - """Returns the location of all active datasets with given label.""" + def get_dataset_location(self,label): # names + """Returns the location of all active datasets with given label.""" #ToDo: Merge path only once at the end '/'.join(listE) path = [] with h5py.File(self.filename,'r') as f: for i in self.active['increments']: From 002383afc2902b95238a9e520e44b1d47118d6b7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 13 Sep 2019 16:01:30 -0700 Subject: [PATCH 23/44] solved problem with postprocessing - to not 'try' with h5py library, it might have another 'try'. Check explicitly for empty argument also some polishing --- processing/post/DADF5_postResults.py | 5 ++-- processing/post/DADF5_vtk_cells.py | 10 ++++--- python/damask/dadf5.py | 40 ++++++++++++---------------- 3 files changed, 26 insertions(+), 29 deletions(-) diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py index fa47805bb..136824282 100755 --- a/processing/post/DADF5_postResults.py +++ b/processing/post/DADF5_postResults.py @@ -1,9 +1,10 @@ #!/usr/bin/env python3 -# -*- coding: UTF-8 no BOM -*- import os -import numpy as np import argparse + +import numpy as np + import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py index 75301386c..aaf8eff26 100755 --- a/processing/post/DADF5_vtk_cells.py +++ b/processing/post/DADF5_vtk_cells.py @@ -1,12 +1,14 @@ #!/usr/bin/env python3 -# -*- coding: UTF-8 no BOM -*- -import os,vtk -import numpy as np +import os import argparse -import damask + +import numpy as np +import vtk from vtk.util import numpy_support +import damask + scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 44238572f..d428cfb59 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -100,18 +100,14 @@ class DADF5(): Get groups that are currently considered for evaluation. """ groups = [] - for i,x in enumerate(self.active['increments']): - group_inc = 'inc{:05}'.format(self.active['increments'][i]['inc']) #ToDo: Merge path only once at the end '/'.join(listE) + for i in self.active['increments']: + group_inc = 'inc{:05}'.format(i['inc']) #ToDo: Merge path only once at the end '/'.join(listE) for c in self.active['constituents']: - group_constituent = group_inc+'/constituent/'+c for t in self.active['c_output_types']: - group_output_types = group_constituent+'/'+t - groups.append(group_output_types) + groups.append('/'.join([group_inc,'constituent',c,t])) for m in self.active['materialpoints']: - group_materialpoint = group_inc+'/materialpoint/'+m for t in self.active['m_output_types']: - group_output_types = group_materialpoint+'/'+t - groups.append(group_output_types) + groups.append('/'.join([group_inc,'materialpoint',m,t])) return groups @@ -150,20 +146,20 @@ class DADF5(): group_inc = 'inc{:05}'.format(i['inc']) for c in self.active['constituents']: - group_constituent = group_inc+'/constituent/'+c for t in self.active['c_output_types']: try: - f[group_constituent+'/'+t+'/'+label] - path.append(group_constituent+'/'+t+'/'+label) + p = '/'.join([group_inc,'constituent',c,t,label]) + f[p] + path.append(p) except KeyError as e: print('unable to locate constituents dataset: '+ str(e)) for m in self.active['materialpoints']: - group_materialpoint = group_inc+'/materialpoint/'+m for t in self.active['m_output_types']: try: - f[group_materialpoint+'/'+t+'/'+label] - path.append(group_materialpoint+'/'+t+'/'+label) + p = '/'.join([group_inc,'materialpoint',m,t,label]) + f[p] + path.append(p) except KeyError as e: print('unable to locate materialpoints dataset: '+ str(e)) @@ -182,24 +178,22 @@ class DADF5(): dataset = np.full(shape,np.nan) for pa in path: label = pa.split('/')[2] - try: - p = np.where(f['mapping/cellResults/constituent'][:,c]['Name'] == str.encode(label))[0] + + p = np.where(f['mapping/cellResults/constituent'][:,c]['Name'] == str.encode(label))[0] + if len(p)>0: u = (f['mapping/cellResults/constituent'][p,c]['Position']) a = np.array(f[pa]) if len(a.shape) == 1: a=a.reshape([a.shape[0],1]) dataset[p,:] = a[u,:] - except KeyError as e: - print('unable to read constituent: '+ str(e)) - try: - p = np.where(f['mapping/cellResults/materialpoint']['Name'] == str.encode(label))[0] + + p = np.where(f['mapping/cellResults/materialpoint']['Name'] == str.encode(label))[0] + if len(p)>0: u = (f['mapping/cellResults/materialpoint'][p.tolist()]['Position']) a = np.array(f[pa]) if len(a.shape) == 1: a=a.reshape([a.shape[0],1]) dataset[p,:] = a[u,:] - except KeyError as e: - print('unable to read materialpoint: '+ str(e)) return dataset @@ -424,7 +418,7 @@ class DADF5(): N_not_calculated = len(todo) while N_not_calculated > 0: result = results.get() - with h5py.File(self.filename,self.mode) as f: # write to file + with h5py.File(self.filename,'a') as f: # write to file dataset_out = f[result['group']].create_dataset(result['label'],data=result['data']) for k in result['meta'].keys(): dataset_out.attrs[k] = result['meta'][k] From c5006e264b92be816a4e2e10614eaf15ce9f46ef Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 13 Sep 2019 16:06:47 -0700 Subject: [PATCH 24/44] handling prospector complaints --- python/damask/dadf5.py | 27 ++++++++++++--------------- 1 file changed, 12 insertions(+), 15 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index d428cfb59..20ea08137 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -15,10 +15,7 @@ class DADF5(): """ # ------------------------------------------------------------------ - def __init__(self, - filename, - mode = 'a', - ): + def __init__(self,filename): """ Opens an existing DADF5 file. @@ -28,7 +25,6 @@ class DADF5(): name of the DADF5 file to be openend. """ - with h5py.File(filename,'r') as f: if f.attrs['DADF5-major'] != 0 or f.attrs['DADF5-minor'] != 2: @@ -63,7 +59,7 @@ class DADF5(): self.m_output_types = list(set(self.m_output_types)) # make unique #self.on_air - self.active= {'increments': self.increments, # ToDo:simplify, activity only positions that translate into (no complex types) + self.active= {'increments': self.increments, # ToDo:simplify, activity only positions that translate into (no complex types) 'constituents': self.constituents, 'materialpoints': self.materialpoints, 'constituent': range(self.Nconstituents), # ToDo: stupid naming @@ -73,8 +69,6 @@ class DADF5(): # ToDo: store increments, select icrements (trivial), position, and time self.filename = filename - self.mode = mode - def get_groups(self,l): #group_with_data(datasets) """ @@ -96,9 +90,7 @@ class DADF5(): def get_active_groups(self): # rename: get_groups needed? merge with datasets and have [] and ['*'] - """ - Get groups that are currently considered for evaluation. - """ + """Get groups that are currently considered for evaluation.""" groups = [] for i in self.active['increments']: group_inc = 'inc{:05}'.format(i['inc']) #ToDo: Merge path only once at the end '/'.join(listE) @@ -258,6 +250,7 @@ class DADF5(): def add_norm(self,x,ord=None): """ Adds norm of vector or tensor or magnitude of a scalar. + See numpy.linalg.norm manual for details. """ def norm(x,ord): @@ -385,12 +378,16 @@ class DADF5(): def __add_generic_pointwise(self,func,datasets_requested,extra_args={}): """ General function to add pointwise data. - """ + + Parameters + ---------- + func : function + datasets_requested : list of dictionaries + extra_args : dictitionary + """ def job(args): - """ - Call function with input data + extra arguments, returns results + group. - """ + """Call function with input data + extra arguments, returns results + group.""" args['results'].put({**args['func'](**args['in']),'group':args['group']}) From 38f6609ad779a06522aad3c656affc4f1e54d875 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 13 Sep 2019 18:36:06 -0700 Subject: [PATCH 25/44] high level functions for selecting output tested Cauchy stress calculation (comparison with addCauchy.py) --- python/damask/dadf5.py | 88 +++++++++++++++++++++++++++++++++++++++--- 1 file changed, 82 insertions(+), 6 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 20ea08137..8a697e82a 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -65,10 +65,87 @@ class DADF5(): 'constituent': range(self.Nconstituents), # ToDo: stupid naming 'c_output_types': self.c_output_types, 'm_output_types': self.m_output_types} + + self.filename = filename + + def __on_air_set(self,output,t,p): + valid = set(p) + choice = [output] if isinstance(output,str) else output + self.active[t] = list(valid.intersection(choice)) + + + def __on_air_add(self,output,t,p): + choice = [output] if isinstance(output,str) else output + valid = set(p).intersection(choice) + existing = set(self.active[t]) + self.active[t] = list(existing.add(valid)) + + + def __on_air_del(self,output): + choice = [output] if isinstance(output,str) else output + existing = set(self.active[t]) + self.active[t] = list(existing.remove(choice)) + + + def constitutive_output_set(self,output): + self.__on_air_set(output,'c_output_types',self.c_output_types) + + + def constitutive_output_add(self,output): + self.__on_air_add(output,'c_output_types',self.c_output_types) + + + def constitutive_output_del(self,output): + self.__on_air_del(output,'c_output_types') + + + def materialpoint_output_set(self,output): + self.__on_air_set(output,'m_output_types',self.m_output_types) + + + def materialpoint_output_add(self,output): + self.__on_air_add(output,'m_output_types',self.m_output_types) + + + def materialpoint_output_del(self,output): + self.__on_air_del(output,'m_output_types') + + + def constitutive_set(self,output): + self.__on_air_set(output,'constituents',self.constituents) + + + def constitutive_add(self,output): + self.__on_air_add(output,'constituents',self.constituents) + + + def constitutive_del(self,output): + self.__on_air_del(output,'constituents') + + + def materialpoint_set(self,output): + self.__on_air_set(output,'materialpoints',self.materialpoints) + + + def materialpoint_add(self,output): + self.__on_air_add(output,'materialpoints',self.materialpoints) + + + def materialpoint_del(self,output): + self.__on_air_del(output,'materialpoints') + + def c_it(self): + a = self.active['constituents']#.copy() + for i in a: + self.constitutive_set(i) + yield i + self.constitutive_set(a) + + # ToDo: store increments, select icrements (trivial), position, and time - self.filename = filename + def get_groups(self,l): #group_with_data(datasets) """ @@ -194,14 +271,13 @@ class DADF5(): """ Adds Cauchy stress calculated from 1st Piola-Kirchhoff stress and deformation gradient. - Todo - ---- - The einsum formula is completely untested! - + Resulting tensor is symmetrized as the Cauchy stress should be symmetric. """ def Cauchy(F,P): + sigma = np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F['data']),P['data'],F['data']) + sigma = (sigma + np.einsum('ijk->ikj',sigma))*0.5 # enforce symmetry return { - 'data' : np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F['data']),F['data'],P['data']), + 'data' : sigma, 'label' : 'sigma', 'meta' : { 'Unit' : P['meta']['Unit'], From 69462f81908b7ce2bb4adc0e34e470359388f3ba Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 13 Sep 2019 19:18:41 -0700 Subject: [PATCH 26/44] polished strain calculation agrees up to 1e-4 with results from addStrainTensors. Not too exciting, but ok --- python/damask/dadf5.py | 33 +++++++++++++++++---------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 8a697e82a..3d1a0b567 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -81,7 +81,7 @@ class DADF5(): self.active[t] = list(existing.add(valid)) - def __on_air_del(self,output): + def __on_air_del(self,output,t): choice = [output] if isinstance(output,str) else output existing = set(self.active[t]) self.active[t] = list(existing.remove(choice)) @@ -275,7 +275,7 @@ class DADF5(): """ def Cauchy(F,P): sigma = np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F['data']),P['data'],F['data']) - sigma = (sigma + np.einsum('ijk->ikj',sigma))*0.5 # enforce symmetry + sigma = (sigma + np.einsum('ikj',sigma))*0.5 # enforce symmetry return { 'data' : sigma, 'label' : 'sigma', @@ -294,10 +294,10 @@ class DADF5(): def add_Mises(self,x): - """Adds the equivalent Mises stres of a tensor.""" + """Adds the equivalent Mises stress or strain of a tensor.""" def deviator(x): - if x['meta']['Unit'] == 'Pa': + if x['meta']['Unit'] == 'Pa': #ToDo: Should we use this? Then add_Cauchy and add_strain_tensors also should perform sanity checks factor = 3.0/2.0 elif x['meta']['Unit'] == '-': factor = 2.0/3.0 @@ -306,10 +306,10 @@ class DADF5(): d = x['data'] dev = d - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[d.shape[0],3,3]),np.trace(d,axis1=1,axis2=2)/3.0) - dev_sym = (dev + np.einsum('ikj',dev))*0.5 + #dev_sym = (dev + np.einsum('ikj',dev))*0.5 # ToDo: this is not needed (only if the input is not symmetric, but then the whole concept breaks down) return { - 'data' : np.sqrt(np.einsum('ijk->i',dev_sym**2)*factor), + 'data' : np.sqrt(np.einsum('ijk->i',dev**2)*factor), 'label' : 'dev({})'.format(x['label']), 'meta' : { 'Unit' : x['meta']['Unit'], @@ -348,7 +348,7 @@ class DADF5(): 'label' : 'norm({})'.format(x['label']), 'meta' : { 'Unit' : x['meta']['Unit'], - 'Description' : 'Norm of {} {} ({})'.format(t,x['label'],x['meta']['Description']), + 'Description' : 'Norm of {} {} ({})'.format(t,x['label'],x['meta']['Description']), # ToDo: add details about norm used 'Creator' : 'dadf5.py:add_norm vXXXXX' } } @@ -403,6 +403,10 @@ class DADF5(): """Adds the deviator of a tensor.""" def deviator(x): d = x['data'] + + if not np.all(np.array(d.shape[1:]) == np.array([3,3])): + raise ValueError + return { 'data' : d - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[d.shape[0],3,3]),np.trace(d,axis1=1,axis2=2)/3.0), 'label' : 'dev({})'.format(x['label']), @@ -418,11 +422,11 @@ class DADF5(): self.__add_generic_pointwise(deviator,requested) - def add_strain_tensor(self,t,ord,defgrad='F'): + def add_strain_tensor(self,t,ord,defgrad='F'): #ToDo: Use t and ord """Adds the a strain tensor.""" def strain_tensor(defgrad,t,ord): (U,S,Vh) = np.linalg.svd(defgrad['data']) # singular value decomposition - R_inv = np.einsum('ijk->ikj',np.matmul(U,Vh)) # inverse rotation of polar decomposition + R_inv = np.einsum('ikj',np.matmul(U,Vh)) # inverse rotation of polar decomposition U = np.matmul(R_inv,defgrad['data']) # F = RU (D,V) = np.linalg.eigh((U+np.einsum('ikj',U))*.5) # eigen decomposition (of symmetric(ed) matrix) @@ -431,17 +435,14 @@ class DADF5(): V[neg[0],:,neg[1]] = V[neg[0],:,neg[1]]* -1 # ... and vector d = np.log(D) - a = np.matmul(V,np.einsum('ij,ikj->ijk',d,V)) # this is wrong ... - for j in range(V.shape[0]): # but this is slow ... - a[j,:,:] = np.dot(V[j,:,:],np.dot(np.diag(d[j,:]),V[j,:,:].T)) - print(np.max(a)) - + a = np.matmul(V,np.einsum('ij,ikj->ijk',d,V)) + return { 'data' : a, - 'label' : 'lnV({})'.format(defgrad['label']), + 'label' : 'ln(V)({})'.format(defgrad['label']), 'meta' : { 'Unit' : defgrad['meta']['Unit'], - 'Description' : 'Strain tensor {} ({})'.format(defgrad['label'],defgrad['meta']['Description']), + 'Description' : 'Strain tensor ln(V){} ({})'.format(defgrad['label'],defgrad['meta']['Description']), 'Creator' : 'dadf5.py:add_deviator vXXXXX' } } From 898f53295dbe95bf738655b42002bc9d815ad37b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 13 Sep 2019 21:14:52 -0700 Subject: [PATCH 27/44] iter functions for groups restore original selection after iterating (i.e. temp setting of active) polishing: correct names and more details for HDF5 meta data --- python/damask/dadf5.py | 81 +++++++++++++++++++++++++++++++----------- 1 file changed, 61 insertions(+), 20 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 3d1a0b567..b14754b7c 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -87,18 +87,39 @@ class DADF5(): self.active[t] = list(existing.remove(choice)) - def constitutive_output_set(self,output): + def __on_air_iter(self,t): + a = self.active[t] + last_a = a.copy() + for i in a: + if last_a != self.active[t]: + self.__on_air_set(a,t,a) + raise Exception + self.__on_air_set(i,t,a) + last_a = self.active[t] + yield i + self.__on_air_set(a,t,a) + + + def constituent_output_iter(self): + return self.__on_air_iter('c_output_types') + + + def constituent_output_set(self,output): self.__on_air_set(output,'c_output_types',self.c_output_types) - def constitutive_output_add(self,output): + def constituent_output_add(self,output): self.__on_air_add(output,'c_output_types',self.c_output_types) - def constitutive_output_del(self,output): + def constituent_output_del(self,output): self.__on_air_del(output,'c_output_types') - + + def materialpoint_output_iter(self): + return self.__on_air_iter('m_output_types') + + def materialpoint_output_set(self,output): self.__on_air_set(output,'m_output_types',self.m_output_types) @@ -111,18 +132,26 @@ class DADF5(): self.__on_air_del(output,'m_output_types') - def constitutive_set(self,output): + def constituent_iter(self): + return self.__on_air_iter('constituents') + + + def constituent_set(self,output): self.__on_air_set(output,'constituents',self.constituents) - def constitutive_add(self,output): + def constituent_add(self,output): self.__on_air_add(output,'constituents',self.constituents) - def constitutive_del(self,output): + def constituent_del(self,output): self.__on_air_del(output,'constituents') - + + def materialpoint_iter(self): + return self.__on_air_iter('materialpoints') + + def materialpoint_set(self,output): self.__on_air_set(output,'materialpoints',self.materialpoints) @@ -134,13 +163,6 @@ class DADF5(): def materialpoint_del(self,output): self.__on_air_del(output,'materialpoints') - def c_it(self): - a = self.active['constituents']#.copy() - for i in a: - self.constitutive_set(i) - yield i - self.constitutive_set(a) - # ToDo: store increments, select icrements (trivial), position, and time @@ -310,7 +332,7 @@ class DADF5(): return { 'data' : np.sqrt(np.einsum('ijk->i',dev**2)*factor), - 'label' : 'dev({})'.format(x['label']), + 'label' : 'Mises({})'.format(x['label']), 'meta' : { 'Unit' : x['meta']['Unit'], 'Description' : 'Mises equivalent stress of {} ({})'.format(x['label'],x['meta']['Description']), @@ -329,26 +351,32 @@ class DADF5(): See numpy.linalg.norm manual for details. """ + + def norm(x,ord): + o = ord if len(x['data'].shape) == 1: axis = 0 t = 'scalar' + if o is None: o = 2 elif len(x['data'].shape) == 2: axis = 1 t = 'vector' + if o is None: o = 2 elif len(x['data'].shape) == 3: axis = (1,2) t = 'tensor' + if o is None: o = 'fro' else: raise ValueError return { - 'data' : np.linalg.norm(x['data'],ord=ord,axis=axis,keepdims=True), - 'label' : 'norm({})'.format(x['label']), + 'data' : np.linalg.norm(x['data'],ord=o,axis=axis,keepdims=True), + 'label' : '|{}|_{}'.format(x['label'],ord), 'meta' : { 'Unit' : x['meta']['Unit'], - 'Description' : 'Norm of {} {} ({})'.format(t,x['label'],x['meta']['Description']), # ToDo: add details about norm used + 'Description' : '{}-Norm of {} {} ({})'.format(ord,t,x['label'],x['meta']['Description']), 'Creator' : 'dadf5.py:add_norm vXXXXX' } } @@ -425,6 +453,16 @@ class DADF5(): def add_strain_tensor(self,t,ord,defgrad='F'): #ToDo: Use t and ord """Adds the a strain tensor.""" def strain_tensor(defgrad,t,ord): + # def operator(stretch,strain,eigenvalues): + #"""Albrecht Bertram: Elasticity and Plasticity of Large Deformations An Introduction (3rd Edition, 2012), p. 102""" + #return { + # 'V#ln': np.log(eigenvalues) , + # 'U#ln': np.log(eigenvalues) , + # 'V#Biot': ( np.ones(3,'d') - 1.0/eigenvalues ) , + # 'U#Biot': ( eigenvalues - np.ones(3,'d') ) , + # 'V#Green': ( np.ones(3,'d') - 1.0/eigenvalues/eigenvalues) *0.5, + # 'U#Green': ( eigenvalues*eigenvalues - np.ones(3,'d')) *0.5, + # }[stretch+'#'+strain] (U,S,Vh) = np.linalg.svd(defgrad['data']) # singular value decomposition R_inv = np.einsum('ikj',np.matmul(U,Vh)) # inverse rotation of polar decomposition U = np.matmul(R_inv,defgrad['data']) # F = RU @@ -459,8 +497,11 @@ class DADF5(): Parameters ---------- func : function + Function that calculates a new dataset from one or more datasets per HDF5 group. datasets_requested : list of dictionaries - extra_args : dictitionary + Details of the datasets to be used: label (in HDF5 file) and arg (argument to which the data is parsed in func). + extra_args : dictionary, optional + Any extra arguments parsed to func. """ def job(args): From a6567e0cc6729188403f2fb775a1861cc4d59659 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 13 Sep 2019 21:41:35 -0700 Subject: [PATCH 28/44] safer to use functions don't mess with attributes of the object, they might be renamed --- processing/post/DADF5_postResults.py | 13 +++++-------- processing/post/DADF5_vtk_cells.py | 20 ++++++++------------ 2 files changed, 13 insertions(+), 20 deletions(-) diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py index 136824282..d4d8678ec 100755 --- a/processing/post/DADF5_postResults.py +++ b/processing/post/DADF5_postResults.py @@ -60,11 +60,10 @@ for filename in options.filenames: header+=' 1_pos 2_pos 3_pos' results.active['increments'] = [inc] + for label in options.con: - for o in results.c_output_types: - results.active['c_output_types'] = [o] - for c in results.constituents: - results.active['constituents'] = [c] + for o in results.constituent_output_iter(): + for c in results.constituent_iter(): x = results.get_dataset_location(label) if len(x) == 0: continue @@ -80,10 +79,8 @@ for filename in options.filenames: header+=' '+label for label in options.mat: - for o in results.m_output_types: - results.active['m_output_types'] = [o] - for m in results.materialpoints: - results.active['materialpoints'] = [m] + for o in results.materialpoint_output_iter(): + for m in results.materialpoint_iter(): x = results.get_dataset_location(label) if len(x) == 0: continue diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py index aaf8eff26..533ea6920 100755 --- a/processing/post/DADF5_vtk_cells.py +++ b/processing/post/DADF5_vtk_cells.py @@ -62,14 +62,13 @@ for filename in options.filenames: vtk_data = [] results.active['increments'] = [inc] - results.active['m_output_types'] = [] + results.materialpoint_set([]) + results.constituent_set(results.constituents) for label in options.con: - for o in results.c_output_types: - results.active['c_output_types'] = [o] + for o in results.constituent_output_iter(): if o != 'generic': - for c in results.constituents: - results.active['constituents'] = [c] + for c in results.constituent_iter(): x = results.get_dataset_location(label) if len(x) == 0: continue @@ -79,7 +78,6 @@ for filename in options.filenames: vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) rGrid.GetCellData().AddArray(vtk_data[-1]) else: - results.active['constituents'] = results.constituents x = results.get_dataset_location(label) if len(x) == 0: continue @@ -89,13 +87,12 @@ for filename in options.filenames: vtk_data[-1].SetName('1_'+x[0].split('/')[1]+'/generic/'+label) rGrid.GetCellData().AddArray(vtk_data[-1]) - results.active['c_output_types'] = [] + results.constituent_set([]) + results.materialpoint_set(results.materialpoints) for label in options.mat: - for o in results.m_output_types: - results.active['m_output_types'] = [o] + for o in results.materialpoint_output_iter(): if o != 'generic': - for m in results.materialpoints: - results.active['materialpoints'] = [m] + for m in results.materialpoint_iter(): x = results.get_dataset_location(label) if len(x) == 0: continue @@ -105,7 +102,6 @@ for filename in options.filenames: vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) rGrid.GetCellData().AddArray(vtk_data[-1]) else: - results.active['materialpoints'] = results.materialpoints x = results.get_dataset_location(label) if len(x) == 0: continue From c13db4b3cafd22c2f20a6d0bef38506ffde44e51 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 13 Sep 2019 21:49:22 -0700 Subject: [PATCH 29/44] renaming on_air/active are not clear. visible seems to be the most appropriate name --- processing/post/DADF5_postResults.py | 2 +- processing/post/DADF5_vtk_cells.py | 2 +- python/damask/dadf5.py | 95 ++++++++++++++-------------- 3 files changed, 49 insertions(+), 50 deletions(-) diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py index d4d8678ec..4330690f9 100755 --- a/processing/post/DADF5_postResults.py +++ b/processing/post/DADF5_postResults.py @@ -59,7 +59,7 @@ for filename in options.filenames: data = np.concatenate((data,coords),1) header+=' 1_pos 2_pos 3_pos' - results.active['increments'] = [inc] + results.visible['increments'] = [inc] for label in options.con: for o in results.constituent_output_iter(): diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py index 533ea6920..fbd89180b 100755 --- a/processing/post/DADF5_vtk_cells.py +++ b/processing/post/DADF5_vtk_cells.py @@ -60,7 +60,7 @@ for filename in options.filenames: for i,inc in enumerate(results.increments): print('Output step {}/{}'.format(i+1,len(results.increments))) vtk_data = [] - results.active['increments'] = [inc] + results.visible['increments'] = [inc] results.materialpoint_set([]) results.constituent_set(results.constituents) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index b14754b7c..15b25dde8 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -58,8 +58,7 @@ class DADF5(): self.m_output_types.append(o) self.m_output_types = list(set(self.m_output_types)) # make unique - #self.on_air - self.active= {'increments': self.increments, # ToDo:simplify, activity only positions that translate into (no complex types) + self.visible= {'increments': self.increments, # ToDo:simplify, activity only positions that translate into (no complex types) 'constituents': self.constituents, 'materialpoints': self.materialpoints, 'constituent': range(self.Nconstituents), # ToDo: stupid naming @@ -68,100 +67,100 @@ class DADF5(): self.filename = filename - def __on_air_set(self,output,t,p): + def __visible_set(self,output,t,p): valid = set(p) choice = [output] if isinstance(output,str) else output - self.active[t] = list(valid.intersection(choice)) + self.visible[t] = list(valid.intersection(choice)) - def __on_air_add(self,output,t,p): + def __visible_add(self,output,t,p): choice = [output] if isinstance(output,str) else output valid = set(p).intersection(choice) - existing = set(self.active[t]) - self.active[t] = list(existing.add(valid)) + existing = set(self.visible[t]) + self.visible[t] = list(existing.add(valid)) - def __on_air_del(self,output,t): + def __visible_del(self,output,t): choice = [output] if isinstance(output,str) else output - existing = set(self.active[t]) - self.active[t] = list(existing.remove(choice)) + existing = set(self.visible[t]) + self.visible[t] = list(existing.remove(choice)) - def __on_air_iter(self,t): - a = self.active[t] + def __visible_iter(self,t): + a = self.visible[t] last_a = a.copy() for i in a: - if last_a != self.active[t]: - self.__on_air_set(a,t,a) + if last_a != self.visible[t]: + self.__visible_set(a,t,a) raise Exception - self.__on_air_set(i,t,a) - last_a = self.active[t] + self.__visible_set(i,t,a) + last_a = self.visible[t] yield i - self.__on_air_set(a,t,a) + self.__visible_set(a,t,a) def constituent_output_iter(self): - return self.__on_air_iter('c_output_types') + return self.__visible_iter('c_output_types') def constituent_output_set(self,output): - self.__on_air_set(output,'c_output_types',self.c_output_types) + self.__visible_set(output,'c_output_types',self.c_output_types) def constituent_output_add(self,output): - self.__on_air_add(output,'c_output_types',self.c_output_types) + self.__visible_add(output,'c_output_types',self.c_output_types) def constituent_output_del(self,output): - self.__on_air_del(output,'c_output_types') + self.__visible_del(output,'c_output_types') def materialpoint_output_iter(self): - return self.__on_air_iter('m_output_types') + return self.__visible_iter('m_output_types') def materialpoint_output_set(self,output): - self.__on_air_set(output,'m_output_types',self.m_output_types) + self.__visible_set(output,'m_output_types',self.m_output_types) def materialpoint_output_add(self,output): - self.__on_air_add(output,'m_output_types',self.m_output_types) + self.__visible_add(output,'m_output_types',self.m_output_types) def materialpoint_output_del(self,output): - self.__on_air_del(output,'m_output_types') + self.__visible_del(output,'m_output_types') def constituent_iter(self): - return self.__on_air_iter('constituents') + return self.__visible_iter('constituents') def constituent_set(self,output): - self.__on_air_set(output,'constituents',self.constituents) + self.__visible_set(output,'constituents',self.constituents) def constituent_add(self,output): - self.__on_air_add(output,'constituents',self.constituents) + self.__visible_add(output,'constituents',self.constituents) def constituent_del(self,output): - self.__on_air_del(output,'constituents') + self.__visible_del(output,'constituents') def materialpoint_iter(self): - return self.__on_air_iter('materialpoints') + return self.__visible_iter('materialpoints') def materialpoint_set(self,output): - self.__on_air_set(output,'materialpoints',self.materialpoints) + self.__visible_set(output,'materialpoints',self.materialpoints) def materialpoint_add(self,output): - self.__on_air_add(output,'materialpoints',self.materialpoints) + self.__visible_add(output,'materialpoints',self.materialpoints) def materialpoint_del(self,output): - self.__on_air_del(output,'materialpoints') + self.__visible_del(output,'materialpoints') @@ -191,13 +190,13 @@ class DADF5(): def get_active_groups(self): # rename: get_groups needed? merge with datasets and have [] and ['*'] """Get groups that are currently considered for evaluation.""" groups = [] - for i in self.active['increments']: + for i in self.visible['increments']: group_inc = 'inc{:05}'.format(i['inc']) #ToDo: Merge path only once at the end '/'.join(listE) - for c in self.active['constituents']: - for t in self.active['c_output_types']: + for c in self.visible['constituents']: + for t in self.visible['c_output_types']: groups.append('/'.join([group_inc,'constituent',c,t])) - for m in self.active['materialpoints']: - for t in self.active['m_output_types']: + for m in self.visible['materialpoints']: + for t in self.visible['m_output_types']: groups.append('/'.join([group_inc,'materialpoint',m,t])) return groups @@ -205,11 +204,11 @@ class DADF5(): def list_data(self): # print_datasets and have [] and ['*'], loop over all increment, soll auf anderen basieren (get groups with sternchen) """Shows information on all active datasets in the file.""" with h5py.File(self.filename,'r') as f: - group_inc = 'inc{:05}'.format(self.active['increments'][0]['inc']) #ToDo: Merge path only once at the end '/'.join(listE) - for c in self.active['constituents']: + group_inc = 'inc{:05}'.format(self.visible['increments'][0]['inc']) #ToDo: Merge path only once at the end '/'.join(listE) + for c in self.visible['constituents']: print('\n'+c) group_constituent = group_inc+'/constituent/'+c - for t in self.active['c_output_types']: + for t in self.visible['c_output_types']: print(' {}'.format(t)) group_output_types = group_constituent+'/'+t try: @@ -217,9 +216,9 @@ class DADF5(): print(' {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode())) except KeyError: pass - for m in self.active['materialpoints']: + for m in self.visible['materialpoints']: group_materialpoint = group_inc+'/materialpoint/'+m - for t in self.active['m_output_types']: + for t in self.visible['m_output_types']: print(' {}'.format(t)) group_output_types = group_materialpoint+'/'+t try: @@ -233,11 +232,11 @@ class DADF5(): """Returns the location of all active datasets with given label.""" #ToDo: Merge path only once at the end '/'.join(listE) path = [] with h5py.File(self.filename,'r') as f: - for i in self.active['increments']: + for i in self.visible['increments']: group_inc = 'inc{:05}'.format(i['inc']) - for c in self.active['constituents']: - for t in self.active['c_output_types']: + for c in self.visible['constituents']: + for t in self.visible['c_output_types']: try: p = '/'.join([group_inc,'constituent',c,t,label]) f[p] @@ -245,8 +244,8 @@ class DADF5(): except KeyError as e: print('unable to locate constituents dataset: '+ str(e)) - for m in self.active['materialpoints']: - for t in self.active['m_output_types']: + for m in self.visible['materialpoints']: + for t in self.visible['m_output_types']: try: p = '/'.join([group_inc,'materialpoint',m,t,label]) f[p] From b85ac11c490a7e452f4f2b6346c4a85d7023a5c5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 14 Sep 2019 08:52:07 -0700 Subject: [PATCH 30/44] implemented ideas from discussion with Philip group matching unified and with wildcard support time step handling should become more convenient (WIP) add_norm can not compute abs of scalar, added function for that general polishing here and there --- processing/post/DADF5_postResults.py | 12 +-- python/damask/dadf5.py | 155 ++++++++++++++++++--------- 2 files changed, 110 insertions(+), 57 deletions(-) diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py index 4330690f9..d1553f025 100755 --- a/processing/post/DADF5_postResults.py +++ b/processing/post/DADF5_postResults.py @@ -24,9 +24,9 @@ parser.add_argument('filenames', nargs='+', parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string', help='name of subdirectory to hold output') parser.add_argument('--mat', nargs='+', - help='labels for materialpoint/homogenization',dest='mat') + help='labels for materialpoint',dest='mat') parser.add_argument('--con', nargs='+', - help='labels for constituent/crystallite/constitutive',dest='con') + help='labels for constituent',dest='con') options = parser.parse_args() @@ -67,11 +67,9 @@ for filename in options.filenames: x = results.get_dataset_location(label) if len(x) == 0: continue - label = x[0].split('/')[-1] array = results.read_dataset(x,0) d = int(np.product(np.shape(array)[1:])) - array = np.reshape(array,[np.product(results.grid),d]) - data = np.concatenate((data,array),1) + data = np.concatenate((data,np.reshape(array,[np.product(results.grid),d])),1) if d>1: header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)]) @@ -84,11 +82,9 @@ for filename in options.filenames: x = results.get_dataset_location(label) if len(x) == 0: continue - label = x[0].split('/')[-1] array = results.read_dataset(x,0) d = int(np.product(np.shape(array)[1:])) - array = np.reshape(array,[np.product(results.grid),d]) - data = np.concatenate((data,array),1) + data = np.concatenate((data,np.reshape(array,[np.product(results.grid),d])),1) if d>1: header+= ''.join([' {}_{}'.format(j+1,label) for j in range(d)]) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 15b25dde8..256135c66 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -1,5 +1,6 @@ from queue import Queue import re +import glob import h5py import numpy as np @@ -37,9 +38,11 @@ class DADF5(): self.size = f['geometry'].attrs['size'] r=re.compile('inc[0-9]+') - self.increments = [{'inc': int(u[3:]), + self.time_information = [{'inc': int(u[3:]), 'time': round(f[u].attrs['time/s'],12), } for u in f.keys() if r.match(u)] + + self.increments = self.time_information.copy() # unify later self.Nmaterialpoints, self.Nconstituents = np.shape(f['mapping/cellResults/constituent']) self.materialpoints = [m.decode() for m in np.unique(f['mapping/cellResults/materialpoint']['Name'])] @@ -58,7 +61,7 @@ class DADF5(): self.m_output_types.append(o) self.m_output_types = list(set(self.m_output_types)) # make unique - self.visible= {'increments': self.increments, # ToDo:simplify, activity only positions that translate into (no complex types) + self.visible= {'increments': self.increments, # ToDo:simplify, activity only positions that translate into (no complex types) 'constituents': self.constituents, 'materialpoints': self.materialpoints, 'constituent': range(self.Nconstituents), # ToDo: stupid naming @@ -97,8 +100,25 @@ class DADF5(): last_a = self.visible[t] yield i self.__visible_set(a,t,a) - - + + + def increment_set_by_time(self,start,end): + for t in self.time_information: + if start<= t['time']< end: + print(t) + + + def increment_set_by_position(self,start,end): + for t in self.time_information[start:end]: + print(t) + + + def increment_set(self,start,end): + for t in self.time_information: + if start<= t['inc']< end: + print(t) + + def constituent_output_iter(self): return self.__visible_iter('c_output_types') @@ -167,39 +187,60 @@ class DADF5(): # ToDo: store increments, select icrements (trivial), position, and time - - def get_groups(self,l): #group_with_data(datasets) + def groups_with_datasets(self,datasets): """ Get groups that contain all requested datasets. + Only groups within inc?????/constituent/*_*/* inc?????/materialpoint/*_*/* + are considered as they contain the data. + Single strings will be treated as list with one entry. + + Wild card matching is allowed, but the number of arguments need to fit. + Parameters ---------- - l : list of str - Names of datasets that need to be located in the group. - - """ - groups = [] - if type(l) is not list: - raise TypeError('Candidates should be given as a list') - with h5py.File(self.filename,'r') as f: - for g in self.get_active_groups(): - if set(l).issubset(f[g].keys()): groups.append(g) - return groups - - - def get_active_groups(self): # rename: get_groups needed? merge with datasets and have [] and ['*'] - """Get groups that are currently considered for evaluation.""" - groups = [] - for i in self.visible['increments']: - group_inc = 'inc{:05}'.format(i['inc']) #ToDo: Merge path only once at the end '/'.join(listE) - for c in self.visible['constituents']: - for t in self.visible['c_output_types']: - groups.append('/'.join([group_inc,'constituent',c,t])) - for m in self.visible['materialpoints']: - for t in self.visible['m_output_types']: - groups.append('/'.join([group_inc,'materialpoint',m,t])) - return groups + datasets : iterable or str or boolean + Examples + -------- + datasets = False matches no group + datasets = True matches all groups + datasets = ['F','P'] matches a group with ['F','P','sigma'] + datasets = ['*','P'] matches a group with ['F','P','sigma'] + datasets = ['*'] does not matche a group with ['F','P','sigma'] + datasets = ['*','*'] does not matche a group with ['F','P','sigma'] + datasets = ['*','*','*'] matches a group with ['F','P','sigma'] + + """ + if datasets is False: return [] + if isinstance(datasets,str): + s = [datasets] + else: + s = datasets + + groups = [] + + with h5py.File(self.filename,'r') as f: + for i in self.visible['increments']: + group_inc = 'inc{:05}'.format(i['inc']) #ToDo: Merge path only once at the end '/'.join(listE) + for c in self.constituent_iter(): + for t in self.constituent_output_iter(): + group = '/'.join([group_inc,'constituent',c,t]) + if datasets is True: + groups.append(group) + else: + match = [e for e_ in [glob.fnmatch.filter(f[group].keys(),s) for s in datasets] for e in e_] + if len(set(match)) == len(s) : groups.append(group) + for m in self.materialpoint_iter(): + for t in self.materialpoint_output_iter(): + group = '/'.join([group_inc,'materialpoint',m,t]) + if datasets is True: + groups.append(group) + else: + match = [e for e_ in [glob.fnmatch.filter(f[group].keys(),s) for s in datasets] for e in e_] + if len(set(match)) == len(s) : groups.append(group) + return groups + def list_data(self): # print_datasets and have [] and ['*'], loop over all increment, soll auf anderen basieren (get groups with sternchen) """Shows information on all active datasets in the file.""" @@ -316,14 +357,17 @@ class DADF5(): def add_Mises(self,x): """Adds the equivalent Mises stress or strain of a tensor.""" - def deviator(x): + def Mises(x): - if x['meta']['Unit'] == 'Pa': #ToDo: Should we use this? Then add_Cauchy and add_strain_tensors also should perform sanity checks + if x['meta']['Unit'] == b'Pa': #ToDo: Should we use this? Then add_Cauchy and add_strain_tensors also should perform sanity checks factor = 3.0/2.0 - elif x['meta']['Unit'] == '-': + t = 'stress' + elif x['meta']['Unit'] == b'1': factor = 2.0/3.0 + t = 'strain' else: - ValueError + print(x['meta']['Unit']) + raise ValueError d = x['data'] dev = d - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[d.shape[0],3,3]),np.trace(d,axis1=1,axis2=2)/3.0) @@ -331,35 +375,29 @@ class DADF5(): return { 'data' : np.sqrt(np.einsum('ijk->i',dev**2)*factor), - 'label' : 'Mises({})'.format(x['label']), + 'label' : '{}_vM'.format(x['label']), 'meta' : { 'Unit' : x['meta']['Unit'], - 'Description' : 'Mises equivalent stress of {} ({})'.format(x['label'],x['meta']['Description']), + 'Description' : 'Mises equivalent {} of {} ({})'.format(t,x['label'],x['meta']['Description']), 'Creator' : 'dadf5.py:add_Mises_stress vXXXXX' } } requested = [{'label':x,'arg':'x'}] - self.__add_generic_pointwise(deviator,requested) + self.__add_generic_pointwise(Mises,requested) def add_norm(self,x,ord=None): """ - Adds norm of vector or tensor or magnitude of a scalar. + Adds norm of vector or tensor. See numpy.linalg.norm manual for details. """ - - def norm(x,ord): o = ord - if len(x['data'].shape) == 1: - axis = 0 - t = 'scalar' - if o is None: o = 2 - elif len(x['data'].shape) == 2: + if len(x['data'].shape) == 2: axis = 1 t = 'vector' if o is None: o = 2 @@ -372,7 +410,7 @@ class DADF5(): return { 'data' : np.linalg.norm(x['data'],ord=o,axis=axis,keepdims=True), - 'label' : '|{}|_{}'.format(x['label'],ord), + 'label' : '|{}|_{}'.format(x['label'],o), 'meta' : { 'Unit' : x['meta']['Unit'], 'Description' : '{}-Norm of {} {} ({})'.format(ord,t,x['label'],x['meta']['Description']), @@ -381,8 +419,27 @@ class DADF5(): } requested = [{'label':x,'arg':'x'}] - + self.__add_generic_pointwise(norm,requested,{'ord':ord}) + + + def add_absolute(self,x): + """Adds absolute value.""" + def absolute(x): + + return { + 'data' : np.abs(x['data']), + 'label' : '|{}|'.format(x['label']), + 'meta' : { + 'Unit' : x['meta']['Unit'], + 'Description' : 'Absolute value of {} ({})'.format(x['label'],x['meta']['Description']), + 'Creator' : 'dadf5.py:add_abs vXXXXX' + } + } + + requested = [{'label':x,'arg':'x'}] + + self.__add_generic_pointwise(absolute,requested) def add_determinant(self,x): @@ -516,7 +573,7 @@ class DADF5(): todo = [] # ToDo: It would be more memory efficient to read only from file when required, i.e. do to it in pool.add_task - for group in self.get_groups([d['label'] for d in datasets_requested]): + for group in self.groups_with_datasets([d['label'] for d in datasets_requested]): with h5py.File(self.filename,'r') as f: datasets_in = {} for d in datasets_requested: From e4e9c5f55804411fdbc8cb23806153f73b19ab6f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 14 Sep 2019 09:38:04 -0700 Subject: [PATCH 31/44] support glob matching more flexibility in selecting active datasets and groups --- python/damask/dadf5.py | 44 +++++++++++++++++++++++++++++++----------- 1 file changed, 33 insertions(+), 11 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 256135c66..25987704b 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -71,22 +71,47 @@ class DADF5(): self.filename = filename def __visible_set(self,output,t,p): - valid = set(p) + """Sets visible.""" + # allow True/False and string arguments + if output is True: + output = ['*'] + elif output is False: + output = [] choice = [output] if isinstance(output,str) else output - self.visible[t] = list(valid.intersection(choice)) + + valid = [e for e_ in [glob.fnmatch.filter(p,s) for s in choice] for e in e_] + + self.visible[t] = valid def __visible_add(self,output,t,p): - choice = [output] if isinstance(output,str) else output - valid = set(p).intersection(choice) + """Adds from visible.""" + # allow True/False and string arguments + if output is True: + output = ['*'] + elif output is False: + output = [] + choice = [output] if isinstance(output,str) else output + existing = set(self.visible[t]) - self.visible[t] = list(existing.add(valid)) + valid = [e for e_ in [glob.fnmatch.filter(p,s) for s in choice] for e in e_] + + self.visible[t] = list(existing.union(valid)) def __visible_del(self,output,t): - choice = [output] if isinstance(output,str) else output + """Deletes from visible.""" + # allow True/False and string arguments + if output is True: + output = ['*'] + elif output is False: + output = [] + choice = [output] if isinstance(output,str) else output + existing = set(self.visible[t]) - self.visible[t] = list(existing.remove(choice)) + valid = [e for e_ in [glob.fnmatch.filter(existing,s) for s in choice] for e in e_] + + self.visible[t] = list(existing.difference_update(valid)) def __visible_iter(self,t): @@ -102,6 +127,7 @@ class DADF5(): self.__visible_set(a,t,a) +# ToDo: store increments, select icrements (trivial), position, and time def increment_set_by_time(self,start,end): for t in self.time_information: if start<= t['time']< end: @@ -181,10 +207,6 @@ class DADF5(): def materialpoint_del(self,output): self.__visible_del(output,'materialpoints') - - - -# ToDo: store increments, select icrements (trivial), position, and time def groups_with_datasets(self,datasets): From 8251725bceb9530b343043900cfaa3ec574a3d44 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 14 Sep 2019 10:53:33 -0700 Subject: [PATCH 32/44] WIP: different norm types --- python/damask/dadf5.py | 30 +++++++++++++++++------------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 25987704b..8354015f5 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -85,7 +85,7 @@ class DADF5(): def __visible_add(self,output,t,p): - """Adds from visible.""" + """Adds to visible.""" # allow True/False and string arguments if output is True: output = ['*'] @@ -529,18 +529,22 @@ class DADF5(): def add_strain_tensor(self,t,ord,defgrad='F'): #ToDo: Use t and ord - """Adds the a strain tensor.""" + """ + Adds the a strain tensor. + + Albrecht Bertram: Elasticity and Plasticity of Large Deformations An Introduction (3rd Edition, 2012), p. 102. + """ def strain_tensor(defgrad,t,ord): - # def operator(stretch,strain,eigenvalues): - #"""Albrecht Bertram: Elasticity and Plasticity of Large Deformations An Introduction (3rd Edition, 2012), p. 102""" - #return { - # 'V#ln': np.log(eigenvalues) , - # 'U#ln': np.log(eigenvalues) , - # 'V#Biot': ( np.ones(3,'d') - 1.0/eigenvalues ) , - # 'U#Biot': ( eigenvalues - np.ones(3,'d') ) , - # 'V#Green': ( np.ones(3,'d') - 1.0/eigenvalues/eigenvalues) *0.5, - # 'U#Green': ( eigenvalues*eigenvalues - np.ones(3,'d')) *0.5, - # }[stretch+'#'+strain] + + operator = { + 'V#ln': lambda V: np.log(V), + 'U#ln': lambda U: np.log(U), + 'V#Biot': lambda V: np.broadcast_to(np.ones(3),[V.shape[0],3]) - 1.0/V, + 'U#Biot': lambda U: U - np.broadcast_to(np.ones(3),[U.shape[0],3]), + 'V#Green':lambda V: np.broadcast_to(np.ones(3),[V.shape[0],3]) - 1.0/V**2, + 'U#Biot': lambda U: U**2 - np.broadcast_to(np.ones(3),[U.shape[0],3]), + } + (U,S,Vh) = np.linalg.svd(defgrad['data']) # singular value decomposition R_inv = np.einsum('ikj',np.matmul(U,Vh)) # inverse rotation of polar decomposition U = np.matmul(R_inv,defgrad['data']) # F = RU @@ -550,7 +554,7 @@ class DADF5(): D[neg[0],neg[1]] = D[neg[0],neg[1]]* -1 # ... flip value ... V[neg[0],:,neg[1]] = V[neg[0],:,neg[1]]* -1 # ... and vector - d = np.log(D) + d = operator['V#ln'](D) a = np.matmul(V,np.einsum('ij,ikj->ijk',d,V)) return { From d413aef7c3ee29a776bddc756293567979ffade5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 14 Sep 2019 12:00:23 -0700 Subject: [PATCH 33/44] functionality for general calculations on datasets currently limited to vectorized expressions. --- python/damask/dadf5.py | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 8354015f5..ecfd04374 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -526,6 +526,38 @@ class DADF5(): requested = [{'label':x,'arg':'x'}] self.__add_generic_pointwise(deviator,requested) + + + def add_calculation(self,formula,label,unit='n/a',description=None,vectorized=True): + """ + General formula. + + Works currently only for vectorized expressions + + """ + if vectorized is not True: + raise NotImplementedError + + def calculation(**kwargs): + + formula = kwargs['formula'] + for d in re.findall(r'#(.*?)#',formula): + formula = re.sub('#{}#'.format(d),"kwargs['{}']['data']".format(d),formula) + + return { + 'data' : eval(formula), + 'label' : kwargs['label'], + 'meta' : { + 'Unit' : kwargs['unit'], + 'Description' : '{}'.format(kwargs['description']), + 'Creator' : 'dadf5.py:add_calculation vXXXXX' + } + } + + requested = [{'label':d,'arg':d} for d in re.findall(r'#(.*?)#',formula)] # datasets used in the formula + pass_through = {'formula':formula,'label':label,'unit':unit,'description':description} + + self.__add_generic_pointwise(calculation,requested,pass_through) def add_strain_tensor(self,t,ord,defgrad='F'): #ToDo: Use t and ord From cc20fd0809b4379e89aef4ea2f633bc89c690b9c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 14 Sep 2019 12:48:57 -0700 Subject: [PATCH 34/44] tests use new HDF5 based post processing --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 5c5adbd8c..b17290e28 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 5c5adbd8ccc0210fd6507431db8ec82ecec75352 +Subproject commit b17290e28fadeee9e868e551b2df7e5109514cc7 From 04aefa84bc14f3b8cca62babcd512622a97c9ea1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 15 Sep 2019 15:02:16 -0700 Subject: [PATCH 35/44] more appropriate names --- processing/post/DADF5_vtk_cells.py | 4 +-- python/damask/dadf5.py | 40 +++++++++++++++--------------- 2 files changed, 22 insertions(+), 22 deletions(-) diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py index fbd89180b..f35c51361 100755 --- a/processing/post/DADF5_vtk_cells.py +++ b/processing/post/DADF5_vtk_cells.py @@ -26,9 +26,9 @@ parser.add_argument('filenames', nargs='+', parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string', help='name of subdirectory to hold output') parser.add_argument('--mat', nargs='+', - help='labels for materialpoint/homogenization',dest='mat') + help='labels for materialpoint',dest='mat') parser.add_argument('--con', nargs='+', - help='labels for constituent/crystallite/constitutive',dest='con') + help='labels for constituent',dest='con') options = parser.parse_args() diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index ecfd04374..3f27dd611 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -49,24 +49,24 @@ class DADF5(): self.constituents = [c.decode() for c in np.unique(f['mapping/cellResults/constituent'] ['Name'])] - self.c_output_types = [] + self.con_physics = [] for c in self.constituents: for o in f['inc{:05}/constituent/{}'.format(self.increments[0]['inc'],c)].keys(): - self.c_output_types.append(o) - self.c_output_types = list(set(self.c_output_types)) # make unique + self.con_physics.append(o) + self.con_physics = list(set(self.con_physics)) # make unique - self.m_output_types = [] + self.mat_physics = [] for m in self.materialpoints: for o in f['inc{:05}/materialpoint/{}'.format(self.increments[0]['inc'],m)].keys(): - self.m_output_types.append(o) - self.m_output_types = list(set(self.m_output_types)) # make unique + self.mat_physics.append(o) + self.mat_physics = list(set(self.mat_physics)) # make unique self.visible= {'increments': self.increments, # ToDo:simplify, activity only positions that translate into (no complex types) 'constituents': self.constituents, 'materialpoints': self.materialpoints, 'constituent': range(self.Nconstituents), # ToDo: stupid naming - 'c_output_types': self.c_output_types, - 'm_output_types': self.m_output_types} + 'con_physics': self.con_physics, + 'mat_physics': self.mat_physics} self.filename = filename @@ -146,35 +146,35 @@ class DADF5(): def constituent_output_iter(self): - return self.__visible_iter('c_output_types') + return self.__visible_iter('con_physics') def constituent_output_set(self,output): - self.__visible_set(output,'c_output_types',self.c_output_types) + self.__visible_set(output,'con_physics',self.con_physics) def constituent_output_add(self,output): - self.__visible_add(output,'c_output_types',self.c_output_types) + self.__visible_add(output,'con_physics',self.con_physics) def constituent_output_del(self,output): - self.__visible_del(output,'c_output_types') + self.__visible_del(output,'con_physics') def materialpoint_output_iter(self): - return self.__visible_iter('m_output_types') + return self.__visible_iter('mat_physics') def materialpoint_output_set(self,output): - self.__visible_set(output,'m_output_types',self.m_output_types) + self.__visible_set(output,'mat_physics',self.mat_physics) def materialpoint_output_add(self,output): - self.__visible_add(output,'m_output_types',self.m_output_types) + self.__visible_add(output,'mat_physics',self.mat_physics) def materialpoint_output_del(self,output): - self.__visible_del(output,'m_output_types') + self.__visible_del(output,'mat_physics') def constituent_iter(self): @@ -271,7 +271,7 @@ class DADF5(): for c in self.visible['constituents']: print('\n'+c) group_constituent = group_inc+'/constituent/'+c - for t in self.visible['c_output_types']: + for t in self.visible['con_physics']: print(' {}'.format(t)) group_output_types = group_constituent+'/'+t try: @@ -281,7 +281,7 @@ class DADF5(): pass for m in self.visible['materialpoints']: group_materialpoint = group_inc+'/materialpoint/'+m - for t in self.visible['m_output_types']: + for t in self.visible['mat_physics']: print(' {}'.format(t)) group_output_types = group_materialpoint+'/'+t try: @@ -299,7 +299,7 @@ class DADF5(): group_inc = 'inc{:05}'.format(i['inc']) for c in self.visible['constituents']: - for t in self.visible['c_output_types']: + for t in self.visible['con_physics']: try: p = '/'.join([group_inc,'constituent',c,t,label]) f[p] @@ -308,7 +308,7 @@ class DADF5(): print('unable to locate constituents dataset: '+ str(e)) for m in self.visible['materialpoints']: - for t in self.visible['m_output_types']: + for t in self.visible['mat_physics']: try: p = '/'.join([group_inc,'materialpoint',m,t,label]) f[p] From 88eba27cecac66ff2567c0f43ca98d6f2c161743 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 15 Sep 2019 15:10:32 -0700 Subject: [PATCH 36/44] avoid name duplication + polishing --- python/damask/dadf5.py | 31 +++++++++++++------------------ 1 file changed, 13 insertions(+), 18 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 3f27dd611..60ef12dcb 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -51,14 +51,12 @@ class DADF5(): self.con_physics = [] for c in self.constituents: - for o in f['inc{:05}/constituent/{}'.format(self.increments[0]['inc'],c)].keys(): - self.con_physics.append(o) + self.con_physics += f['inc{:05}/constituent/{}'.format(self.increments[0]['inc'],c)].keys() self.con_physics = list(set(self.con_physics)) # make unique self.mat_physics = [] for m in self.materialpoints: - for o in f['inc{:05}/materialpoint/{}'.format(self.increments[0]['inc'],m)].keys(): - self.mat_physics.append(o) + self.mat_physics += f['inc{:05}/materialpoint/{}'.format(self.increments[0]['inc'],m)].keys() self.mat_physics = list(set(self.mat_physics)) # make unique self.visible= {'increments': self.increments, # ToDo:simplify, activity only positions that translate into (no complex types) @@ -228,17 +226,14 @@ class DADF5(): datasets = False matches no group datasets = True matches all groups datasets = ['F','P'] matches a group with ['F','P','sigma'] - datasets = ['*','P'] matches a group with ['F','P','sigma'] - datasets = ['*'] does not matche a group with ['F','P','sigma'] - datasets = ['*','*'] does not matche a group with ['F','P','sigma'] + datasets = ['*','P'] matches a group with ['F','P'] + datasets = ['*'] does not match a group with ['F','P','sigma'] + datasets = ['*','*'] does not match a group with ['F','P','sigma'] datasets = ['*','*','*'] matches a group with ['F','P','sigma'] """ if datasets is False: return [] - if isinstance(datasets,str): - s = [datasets] - else: - s = datasets + sets = [datasets] if isinstance(datasets,str) else datasets groups = [] @@ -248,19 +243,19 @@ class DADF5(): for c in self.constituent_iter(): for t in self.constituent_output_iter(): group = '/'.join([group_inc,'constituent',c,t]) - if datasets is True: + if sets is True: groups.append(group) else: - match = [e for e_ in [glob.fnmatch.filter(f[group].keys(),s) for s in datasets] for e in e_] - if len(set(match)) == len(s) : groups.append(group) + match = [e for e_ in [glob.fnmatch.filter(f[group].keys(),s) for s in sets] for e in e_] + if len(set(match)) == len(sets) : groups.append(group) for m in self.materialpoint_iter(): for t in self.materialpoint_output_iter(): group = '/'.join([group_inc,'materialpoint',m,t]) - if datasets is True: + if sets is True: groups.append(group) else: - match = [e for e_ in [glob.fnmatch.filter(f[group].keys(),s) for s in datasets] for e in e_] - if len(set(match)) == len(s) : groups.append(group) + match = [e for e_ in [glob.fnmatch.filter(f[group].keys(),s) for s in sets] for e in e_] + if len(set(match)) == len(sets) : groups.append(group) return groups @@ -366,7 +361,7 @@ class DADF5(): 'meta' : { 'Unit' : P['meta']['Unit'], 'Description' : 'Cauchy stress calculated from {} ({}) '.format(P['label'],P['meta']['Description'])+\ - 'and deformation gradient {} ({})'.format(F['label'],P['meta']['Description']), + 'and deformation gradient {} ({})'.format(F['label'],F['meta']['Description']), 'Creator' : 'dadf5.py:add_Cauchy vXXXXX' } } From c76d4d3f87a4fcc75a8ffe45c31992b81f5ffeb4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 15 Sep 2019 15:24:23 -0700 Subject: [PATCH 37/44] avoid unneeded arguments --- python/damask/dadf5.py | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 60ef12dcb..db549b0f0 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -68,7 +68,7 @@ class DADF5(): self.filename = filename - def __visible_set(self,output,t,p): + def __visible_set(self,output,t): """Sets visible.""" # allow True/False and string arguments if output is True: @@ -77,12 +77,12 @@ class DADF5(): output = [] choice = [output] if isinstance(output,str) else output - valid = [e for e_ in [glob.fnmatch.filter(p,s) for s in choice] for e in e_] + valid = [e for e_ in [glob.fnmatch.filter(getattr(self,t),s) for s in choice] for e in e_] self.visible[t] = valid - def __visible_add(self,output,t,p): + def __visible_add(self,output,t): """Adds to visible.""" # allow True/False and string arguments if output is True: @@ -92,7 +92,7 @@ class DADF5(): choice = [output] if isinstance(output,str) else output existing = set(self.visible[t]) - valid = [e for e_ in [glob.fnmatch.filter(p,s) for s in choice] for e in e_] + valid = [e for e_ in [glob.fnmatch.filter(getattr(self,t),s) for s in choice] for e in e_] self.visible[t] = list(existing.union(valid)) @@ -117,12 +117,12 @@ class DADF5(): last_a = a.copy() for i in a: if last_a != self.visible[t]: - self.__visible_set(a,t,a) + self.__visible_set(a,t) raise Exception - self.__visible_set(i,t,a) + self.__visible_set(i,t) last_a = self.visible[t] yield i - self.__visible_set(a,t,a) + self.__visible_set(a,t) # ToDo: store increments, select icrements (trivial), position, and time @@ -148,11 +148,11 @@ class DADF5(): def constituent_output_set(self,output): - self.__visible_set(output,'con_physics',self.con_physics) + self.__visible_set(output,'con_physics') def constituent_output_add(self,output): - self.__visible_add(output,'con_physics',self.con_physics) + self.__visible_add(output,'con_physics') def constituent_output_del(self,output): @@ -164,11 +164,11 @@ class DADF5(): def materialpoint_output_set(self,output): - self.__visible_set(output,'mat_physics',self.mat_physics) + self.__visible_set(output,'mat_physics') def materialpoint_output_add(self,output): - self.__visible_add(output,'mat_physics',self.mat_physics) + self.__visible_add(output,'mat_physics') def materialpoint_output_del(self,output): @@ -180,11 +180,11 @@ class DADF5(): def constituent_set(self,output): - self.__visible_set(output,'constituents',self.constituents) + self.__visible_set(output,'constituents') def constituent_add(self,output): - self.__visible_add(output,'constituents',self.constituents) + self.__visible_add(output,'constituents') def constituent_del(self,output): @@ -196,11 +196,11 @@ class DADF5(): def materialpoint_set(self,output): - self.__visible_set(output,'materialpoints',self.materialpoints) + self.__visible_set(output,'materialpoints') def materialpoint_add(self,output): - self.__visible_add(output,'materialpoints',self.materialpoints) + self.__visible_add(output,'materialpoints') def materialpoint_del(self,output): From f6ac8c995f76a7a3245f83f43fe7ca406241a6f3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 15 Sep 2019 16:00:19 -0700 Subject: [PATCH 38/44] simplified and selected better names --- processing/post/DADF5_postResults.py | 8 +- processing/post/DADF5_vtk_cells.py | 20 ++-- python/damask/dadf5.py | 161 +++++++-------------------- 3 files changed, 57 insertions(+), 132 deletions(-) diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py index d1553f025..2d2749e33 100755 --- a/processing/post/DADF5_postResults.py +++ b/processing/post/DADF5_postResults.py @@ -62,8 +62,8 @@ for filename in options.filenames: results.visible['increments'] = [inc] for label in options.con: - for o in results.constituent_output_iter(): - for c in results.constituent_iter(): + for o in results.iter_visible('con_physics'): + for c in results.iter_visible('constituents'): x = results.get_dataset_location(label) if len(x) == 0: continue @@ -77,8 +77,8 @@ for filename in options.filenames: header+=' '+label for label in options.mat: - for o in results.materialpoint_output_iter(): - for m in results.materialpoint_iter(): + for o in results.iter_visible('mat_physics'): + for m in results.iter_visible('materialpoints'): x = results.get_dataset_location(label) if len(x) == 0: continue diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py index f35c51361..47d78918f 100755 --- a/processing/post/DADF5_vtk_cells.py +++ b/processing/post/DADF5_vtk_cells.py @@ -62,13 +62,13 @@ for filename in options.filenames: vtk_data = [] results.visible['increments'] = [inc] - results.materialpoint_set([]) - results.constituent_set(results.constituents) + results.set_visible('materialpoints',False) + results.set_visible('constituents', True) for label in options.con: - for o in results.constituent_output_iter(): + for o in results.iter_visible('con_physics'): if o != 'generic': - for c in results.constituent_iter(): + for c in results.iter_visible('constituents'): x = results.get_dataset_location(label) if len(x) == 0: continue @@ -84,15 +84,15 @@ for filename in options.filenames: array = results.read_dataset(x,0) shape = [array.shape[0],np.product(array.shape[1:])] vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) - vtk_data[-1].SetName('1_'+x[0].split('/')[1]+'/generic/'+label) + vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) rGrid.GetCellData().AddArray(vtk_data[-1]) - results.constituent_set([]) - results.materialpoint_set(results.materialpoints) + results.set_visible('constituents', False) + results.set_visible('materialpoints',True) for label in options.mat: - for o in results.materialpoint_output_iter(): + for o in results.iter_visible('mat_physics'): if o != 'generic': - for m in results.materialpoint_iter(): + for m in results.iter_visible('materialpoints'): x = results.get_dataset_location(label) if len(x) == 0: continue @@ -108,7 +108,7 @@ for filename in options.filenames: array = results.read_dataset(x,0) shape = [array.shape[0],np.product(array.shape[1:])] vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE)) - vtk_data[-1].SetName('1_'+x[0].split('/')[1]+'/generic/'+label) + vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) rGrid.GetCellData().AddArray(vtk_data[-1]) if results.structured: diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index db549b0f0..9eea34ca4 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -68,61 +68,38 @@ class DADF5(): self.filename = filename - def __visible_set(self,output,t): + + def __manage_visible(self,datasets,what,action): """Sets visible.""" # allow True/False and string arguments - if output is True: - output = ['*'] - elif output is False: - output = [] - choice = [output] if isinstance(output,str) else output + if datasets is True: + datasets = ['*'] + elif datasets is False: + datasets = [] + choice = [datasets] if isinstance(datasets,str) else datasets - valid = [e for e_ in [glob.fnmatch.filter(getattr(self,t),s) for s in choice] for e in e_] + valid = [e for e_ in [glob.fnmatch.filter(getattr(self,what) ,s) for s in choice] for e in e_] + existing = set(self.visible[what]) - self.visible[t] = valid + if action == 'set': + self.visible[what] = valid + elif action == 'add': + self.visible[what] = list(existing.union(valid)) + elif action == 'del': + self.visible[what] = list(existing.difference_update(valid)) - def __visible_add(self,output,t): - """Adds to visible.""" - # allow True/False and string arguments - if output is True: - output = ['*'] - elif output is False: - output = [] - choice = [output] if isinstance(output,str) else output - - existing = set(self.visible[t]) - valid = [e for e_ in [glob.fnmatch.filter(getattr(self,t),s) for s in choice] for e in e_] - - self.visible[t] = list(existing.union(valid)) - - - def __visible_del(self,output,t): - """Deletes from visible.""" - # allow True/False and string arguments - if output is True: - output = ['*'] - elif output is False: - output = [] - choice = [output] if isinstance(output,str) else output - - existing = set(self.visible[t]) - valid = [e for e_ in [glob.fnmatch.filter(existing,s) for s in choice] for e in e_] - - self.visible[t] = list(existing.difference_update(valid)) - - - def __visible_iter(self,t): - a = self.visible[t] - last_a = a.copy() - for i in a: - if last_a != self.visible[t]: - self.__visible_set(a,t) + def iter_visible(self,what): + datasets = self.visible[what] + last_datasets = datasets.copy() + for dataset in datasets: + if last_datasets != self.visible[what]: + self.__manage_visible(datasets,what,'set') raise Exception - self.__visible_set(i,t) - last_a = self.visible[t] - yield i - self.__visible_set(a,t) + self.__manage_visible(dataset,what,'set') + last_datasets = self.visible[what] + yield dataset + self.__manage_visible(datasets,what,'set') # ToDo: store increments, select icrements (trivial), position, and time @@ -141,70 +118,18 @@ class DADF5(): for t in self.time_information: if start<= t['inc']< end: print(t) - - - def constituent_output_iter(self): - return self.__visible_iter('con_physics') - - def constituent_output_set(self,output): - self.__visible_set(output,'con_physics') - - - def constituent_output_add(self,output): - self.__visible_add(output,'con_physics') - - - def constituent_output_del(self,output): - self.__visible_del(output,'con_physics') - - def materialpoint_output_iter(self): - return self.__visible_iter('mat_physics') - - - def materialpoint_output_set(self,output): - self.__visible_set(output,'mat_physics') + def set_visible(self,what,datasets): + self.__manage_visible(datasets,what,'set') - def materialpoint_output_add(self,output): - self.__visible_add(output,'mat_physics') + def add_visible(self,what,datasets): + self.__manage_visible(datasets,what,'add') - def materialpoint_output_del(self,output): - self.__visible_del(output,'mat_physics') - - - def constituent_iter(self): - return self.__visible_iter('constituents') - - - def constituent_set(self,output): - self.__visible_set(output,'constituents') - - - def constituent_add(self,output): - self.__visible_add(output,'constituents') - - - def constituent_del(self,output): - self.__visible_del(output,'constituents') - - - def materialpoint_iter(self): - return self.__visible_iter('materialpoints') - - - def materialpoint_set(self,output): - self.__visible_set(output,'materialpoints') - - - def materialpoint_add(self,output): - self.__visible_add(output,'materialpoints') - - - def materialpoint_del(self,output): - self.__visible_del(output,'materialpoints') + def del_visible(self,what,datasets): + self.__manage_visible(datasets,what,'del') def groups_with_datasets(self,datasets): @@ -240,16 +165,16 @@ class DADF5(): with h5py.File(self.filename,'r') as f: for i in self.visible['increments']: group_inc = 'inc{:05}'.format(i['inc']) #ToDo: Merge path only once at the end '/'.join(listE) - for c in self.constituent_iter(): - for t in self.constituent_output_iter(): + for c in self.iter_visible('constituents'): + for t in self.iter_visible('con_physics'): group = '/'.join([group_inc,'constituent',c,t]) if sets is True: groups.append(group) else: match = [e for e_ in [glob.fnmatch.filter(f[group].keys(),s) for s in sets] for e in e_] if len(set(match)) == len(sets) : groups.append(group) - for m in self.materialpoint_iter(): - for t in self.materialpoint_output_iter(): + for m in self.iter_visible('materialpoints'): + for t in self.iter_visible('mat_physics'): group = '/'.join([group_inc,'materialpoint',m,t]) if sets is True: groups.append(group) @@ -263,10 +188,10 @@ class DADF5(): """Shows information on all active datasets in the file.""" with h5py.File(self.filename,'r') as f: group_inc = 'inc{:05}'.format(self.visible['increments'][0]['inc']) #ToDo: Merge path only once at the end '/'.join(listE) - for c in self.visible['constituents']: + for c in self.iter_visible('constituents'): print('\n'+c) group_constituent = group_inc+'/constituent/'+c - for t in self.visible['con_physics']: + for t in self.iter_visible('con_physics'): print(' {}'.format(t)) group_output_types = group_constituent+'/'+t try: @@ -274,9 +199,9 @@ class DADF5(): print(' {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode())) except KeyError: pass - for m in self.visible['materialpoints']: + for m in self.iter_visible('materialpoints'): group_materialpoint = group_inc+'/materialpoint/'+m - for t in self.visible['mat_physics']: + for t in self.iter_visible('mat_physics'): print(' {}'.format(t)) group_output_types = group_materialpoint+'/'+t try: @@ -293,8 +218,8 @@ class DADF5(): for i in self.visible['increments']: group_inc = 'inc{:05}'.format(i['inc']) - for c in self.visible['constituents']: - for t in self.visible['con_physics']: + for c in self.iter_visible('constituents'): + for t in self.iter_visible('con_physics'): try: p = '/'.join([group_inc,'constituent',c,t,label]) f[p] @@ -302,8 +227,8 @@ class DADF5(): except KeyError as e: print('unable to locate constituents dataset: '+ str(e)) - for m in self.visible['materialpoints']: - for t in self.visible['mat_physics']: + for m in self.iter_visible('materialpoints'): + for t in self.iter_visible('mat_physics'): try: p = '/'.join([group_inc,'materialpoint',m,t,label]) f[p] From 4cedcee0b4fc1e79c29de31227ea2d1d53232638 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 15 Sep 2019 17:08:03 -0700 Subject: [PATCH 39/44] handling of increments follows handling of other 'visible'-items --- processing/post/DADF5_postResults.py | 12 ++--- processing/post/DADF5_vtk_cells.py | 13 +++-- python/damask/dadf5.py | 73 +++++++++------------------- 3 files changed, 35 insertions(+), 63 deletions(-) diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py index 2d2749e33..4ef55a19e 100755 --- a/processing/post/DADF5_postResults.py +++ b/processing/post/DADF5_postResults.py @@ -47,22 +47,20 @@ for filename in options.filenames: coords = np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) - for i,inc in enumerate(results.increments): + for i,inc in enumerate(results.iter_visible('increments')): print('Output step {}/{}'.format(i+1,len(results.increments))) header = '1 header\n' - data = np.array([inc['inc'] for j in range(np.product(results.grid))]).reshape([np.product(results.grid),1]) + data = np.array([int(inc[3:]) for j in range(np.product(results.grid))]).reshape([np.product(results.grid),1]) header+= 'inc' coords = coords.reshape([np.product(results.grid),3]) data = np.concatenate((data,coords),1) header+=' 1_pos 2_pos 3_pos' - - results.visible['increments'] = [inc] for label in options.con: - for o in results.iter_visible('con_physics'): + for p in results.iter_visible('con_physics'): for c in results.iter_visible('constituents'): x = results.get_dataset_location(label) if len(x) == 0: @@ -77,7 +75,7 @@ for filename in options.filenames: header+=' '+label for label in options.mat: - for o in results.iter_visible('mat_physics'): + for p in results.iter_visible('mat_physics'): for m in results.iter_visible('materialpoints'): x = results.get_dataset_location(label) if len(x) == 0: @@ -96,5 +94,5 @@ for filename in options.filenames: os.mkdir(dirname) except FileExistsError: pass - file_out = '{}_inc{:04d}.txt'.format(filename.split('.')[0],inc['inc']) + file_out = '{}_{}.txt'.format(filename.split('.')[0],inc) np.savetxt(os.path.join(dirname,file_out),data,header=header,comments='') diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py index 47d78918f..8c5bb4838 100755 --- a/processing/post/DADF5_vtk_cells.py +++ b/processing/post/DADF5_vtk_cells.py @@ -57,17 +57,16 @@ for filename in options.filenames: rGrid.SetZCoordinates(coordArray[2]) - for i,inc in enumerate(results.increments): + for i,inc in enumerate(results.iter_visible('increments')): print('Output step {}/{}'.format(i+1,len(results.increments))) vtk_data = [] - results.visible['increments'] = [inc] results.set_visible('materialpoints',False) results.set_visible('constituents', True) for label in options.con: - for o in results.iter_visible('con_physics'): - if o != 'generic': + for p in results.iter_visible('con_physics'): + if p != 'generic': for c in results.iter_visible('constituents'): x = results.get_dataset_location(label) if len(x) == 0: @@ -90,8 +89,8 @@ for filename in options.filenames: results.set_visible('constituents', False) results.set_visible('materialpoints',True) for label in options.mat: - for o in results.iter_visible('mat_physics'): - if o != 'generic': + for p in results.iter_visible('mat_physics'): + if p != 'generic': for m in results.iter_visible('materialpoints'): x = results.get_dataset_location(label) if len(x) == 0: @@ -120,7 +119,7 @@ for filename in options.filenames: os.mkdir(dirname) except FileExistsError: pass - file_out = '{}_inc{:04d}.{}'.format(filename.split('.')[0],inc['inc'],writer.GetDefaultFileExtension()) + file_out = '{}_{}.{}'.format(filename.split('.')[0],inc,writer.GetDefaultFileExtension()) writer.SetCompressorTypeToZLib() writer.SetDataModeToBinary() diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 9eea34ca4..ff067be23 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -38,11 +38,7 @@ class DADF5(): self.size = f['geometry'].attrs['size'] r=re.compile('inc[0-9]+') - self.time_information = [{'inc': int(u[3:]), - 'time': round(f[u].attrs['time/s'],12), - } for u in f.keys() if r.match(u)] - - self.increments = self.time_information.copy() # unify later + self.increments = [u for u in f.keys() if r.match(u)] self.Nmaterialpoints, self.Nconstituents = np.shape(f['mapping/cellResults/constituent']) self.materialpoints = [m.decode() for m in np.unique(f['mapping/cellResults/materialpoint']['Name'])] @@ -51,20 +47,20 @@ class DADF5(): self.con_physics = [] for c in self.constituents: - self.con_physics += f['inc{:05}/constituent/{}'.format(self.increments[0]['inc'],c)].keys() + self.con_physics += f['inc00000/constituent/{}'.format(c)].keys() self.con_physics = list(set(self.con_physics)) # make unique self.mat_physics = [] for m in self.materialpoints: - self.mat_physics += f['inc{:05}/materialpoint/{}'.format(self.increments[0]['inc'],m)].keys() + self.mat_physics += f['inc00000/materialpoint/{}'.format(m)].keys() self.mat_physics = list(set(self.mat_physics)) # make unique self.visible= {'increments': self.increments, # ToDo:simplify, activity only positions that translate into (no complex types) 'constituents': self.constituents, 'materialpoints': self.materialpoints, 'constituent': range(self.Nconstituents), # ToDo: stupid naming - 'con_physics': self.con_physics, - 'mat_physics': self.mat_physics} + 'con_physics': self.con_physics, + 'mat_physics': self.mat_physics} self.filename = filename @@ -100,24 +96,6 @@ class DADF5(): last_datasets = self.visible[what] yield dataset self.__manage_visible(datasets,what,'set') - - -# ToDo: store increments, select icrements (trivial), position, and time - def increment_set_by_time(self,start,end): - for t in self.time_information: - if start<= t['time']< end: - print(t) - - - def increment_set_by_position(self,start,end): - for t in self.time_information[start:end]: - print(t) - - - def increment_set(self,start,end): - for t in self.time_information: - if start<= t['inc']< end: - print(t) def set_visible(self,what,datasets): @@ -163,19 +141,18 @@ class DADF5(): groups = [] with h5py.File(self.filename,'r') as f: - for i in self.visible['increments']: - group_inc = 'inc{:05}'.format(i['inc']) #ToDo: Merge path only once at the end '/'.join(listE) + for i in self.iter_visible('increments'): #ToDo: Merge path only once at the end '/'.join(listE) for c in self.iter_visible('constituents'): - for t in self.iter_visible('con_physics'): - group = '/'.join([group_inc,'constituent',c,t]) + for p in self.iter_visible('con_physics'): + group = '/'.join([i,'constituent',c,p]) if sets is True: groups.append(group) else: match = [e for e_ in [glob.fnmatch.filter(f[group].keys(),s) for s in sets] for e in e_] if len(set(match)) == len(sets) : groups.append(group) for m in self.iter_visible('materialpoints'): - for t in self.iter_visible('mat_physics'): - group = '/'.join([group_inc,'materialpoint',m,t]) + for p in self.iter_visible('mat_physics'): + group = '/'.join([i,'materialpoint',m,p]) if sets is True: groups.append(group) else: @@ -187,23 +164,23 @@ class DADF5(): def list_data(self): # print_datasets and have [] and ['*'], loop over all increment, soll auf anderen basieren (get groups with sternchen) """Shows information on all active datasets in the file.""" with h5py.File(self.filename,'r') as f: - group_inc = 'inc{:05}'.format(self.visible['increments'][0]['inc']) #ToDo: Merge path only once at the end '/'.join(listE) + i = 'inc{:05}'.format(0) #ToDo: Merge path only once at the end '/'.join(listE) for c in self.iter_visible('constituents'): print('\n'+c) - group_constituent = group_inc+'/constituent/'+c - for t in self.iter_visible('con_physics'): + group_constituent = i+'/constituent/'+c + for p in self.iter_visible('con_physics'): print(' {}'.format(t)) - group_output_types = group_constituent+'/'+t + group_output_types = group_constituent+'/'+p try: for x in f[group_output_types].keys(): print(' {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode())) except KeyError: pass for m in self.iter_visible('materialpoints'): - group_materialpoint = group_inc+'/materialpoint/'+m - for t in self.iter_visible('mat_physics'): + group_materialpoint = i+'/materialpoint/'+m + for p in self.iter_visible('mat_physics'): print(' {}'.format(t)) - group_output_types = group_materialpoint+'/'+t + group_output_types = group_materialpoint+'/'+p try: for x in f[group_output_types].keys(): print(' {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode())) @@ -215,24 +192,22 @@ class DADF5(): """Returns the location of all active datasets with given label.""" #ToDo: Merge path only once at the end '/'.join(listE) path = [] with h5py.File(self.filename,'r') as f: - for i in self.visible['increments']: - group_inc = 'inc{:05}'.format(i['inc']) - + for i in self.iter_visible('increments'): for c in self.iter_visible('constituents'): for t in self.iter_visible('con_physics'): try: - p = '/'.join([group_inc,'constituent',c,t,label]) - f[p] - path.append(p) + k = '/'.join([i,'constituent',c,t,label]) + f[k] + path.append(k) except KeyError as e: print('unable to locate constituents dataset: '+ str(e)) for m in self.iter_visible('materialpoints'): for t in self.iter_visible('mat_physics'): try: - p = '/'.join([group_inc,'materialpoint',m,t,label]) - f[p] - path.append(p) + k = '/'.join([i,'materialpoint',m,t,label]) + f[k] + path.append(k) except KeyError as e: print('unable to locate materialpoints dataset: '+ str(e)) From b3b710c84878f1f72db1914b00da0774ba52a785 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 15 Sep 2019 19:56:07 -0700 Subject: [PATCH 40/44] polishing --- python/damask/dadf5.py | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index ff067be23..65f63f0ed 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -66,7 +66,7 @@ class DADF5(): def __manage_visible(self,datasets,what,action): - """Sets visible.""" + """Manages the visibility of the groups.""" # allow True/False and string arguments if datasets is True: datasets = ['*'] @@ -86,6 +86,7 @@ class DADF5(): def iter_visible(self,what): + """Iterates over visible items by setting each one visible.""" datasets = self.visible[what] last_datasets = datasets.copy() for dataset in datasets: @@ -161,51 +162,50 @@ class DADF5(): return groups - def list_data(self): # print_datasets and have [] and ['*'], loop over all increment, soll auf anderen basieren (get groups with sternchen) + def list_data(self): """Shows information on all active datasets in the file.""" with h5py.File(self.filename,'r') as f: - i = 'inc{:05}'.format(0) #ToDo: Merge path only once at the end '/'.join(listE) + i = 'inc{:05}'.format(0) for c in self.iter_visible('constituents'): - print('\n'+c) - group_constituent = i+'/constituent/'+c + print('{}'.format(c)) for p in self.iter_visible('con_physics'): - print(' {}'.format(t)) - group_output_types = group_constituent+'/'+p + print(' {}'.format(p)) try: - for x in f[group_output_types].keys(): - print(' {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode())) + k = '/'.join([i,'constituent',c,p]) + for d in f[k].keys(): + print(' {} ({})'.format(d,f[k+'/'+d].attrs['Description'].decode())) except KeyError: pass for m in self.iter_visible('materialpoints'): - group_materialpoint = i+'/materialpoint/'+m + print('{}'.format(m)) for p in self.iter_visible('mat_physics'): - print(' {}'.format(t)) - group_output_types = group_materialpoint+'/'+p + print(' {}'.format(p)) try: - for x in f[group_output_types].keys(): - print(' {} ({})'.format(x,f[group_output_types+'/'+x].attrs['Description'].decode())) + k = '/'.join([i,'materialpoint',m,p]) + for d in f[k].keys(): + print(' {} ({})'.format(d,f[k+'/'+d].attrs['Description'].decode())) except KeyError: pass - def get_dataset_location(self,label): # names - """Returns the location of all active datasets with given label.""" #ToDo: Merge path only once at the end '/'.join(listE) + def get_dataset_location(self,label): + """Returns the location of all active datasets with given label.""" path = [] with h5py.File(self.filename,'r') as f: for i in self.iter_visible('increments'): for c in self.iter_visible('constituents'): - for t in self.iter_visible('con_physics'): + for p in self.iter_visible('con_physics'): try: - k = '/'.join([i,'constituent',c,t,label]) + k = '/'.join([i,'constituent',c,p,label]) f[k] path.append(k) except KeyError as e: print('unable to locate constituents dataset: '+ str(e)) for m in self.iter_visible('materialpoints'): - for t in self.iter_visible('mat_physics'): + for p in self.iter_visible('mat_physics'): try: - k = '/'.join([i,'materialpoint',m,t,label]) + k = '/'.join([i,'materialpoint',m,p,label]) f[k] path.append(k) except KeyError as e: From 86fb0a79420a41a8b85789983b8ac268d478a31d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 15 Sep 2019 20:04:52 -0700 Subject: [PATCH 41/44] time info needed for filtering of increments --- python/damask/dadf5.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 65f63f0ed..091f1fe23 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -38,8 +38,9 @@ class DADF5(): self.size = f['geometry'].attrs['size'] r=re.compile('inc[0-9]+') - self.increments = [u for u in f.keys() if r.match(u)] - + self.increments = [i for i in f.keys() if r.match(i)] + self.times = [round(f[i].attrs['time/s'],12) for i in self.increments] + self.Nmaterialpoints, self.Nconstituents = np.shape(f['mapping/cellResults/constituent']) self.materialpoints = [m.decode() for m in np.unique(f['mapping/cellResults/materialpoint']['Name'])] self.constituents = [c.decode() for c in np.unique(f['mapping/cellResults/constituent'] ['Name'])] From 980c02b0c3aad60a298dcbd55014ca6991dd6a40 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 15 Sep 2019 20:19:14 -0700 Subject: [PATCH 42/44] selected increments by simulation time --- python/damask/dadf5.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 091f1fe23..586ce193f 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -84,6 +84,26 @@ class DADF5(): self.visible[what] = list(existing.union(valid)) elif action == 'del': self.visible[what] = list(existing.difference_update(valid)) + + + def __time_to_inc(self,start,end): + selected = [] + for i,time in enumerate(self.times): + if start <= time < end: + selected.append(self.increments[i]) + return selected + + + def set_by_time(self,start,end): + self.__manage_visible(self.__time_to_inc(start,end),'increments','set') + + + def add_by_time(self,start,end): + self.__manage_visible(self.__time_to_inc(start,end),'increments','add') + + + def del_by_time(self,start,end): + self.__manage_visible(self.__time_to_inc(start,end),'increments','del') def iter_visible(self,what): From b2cab84efde72e02140894b324c8cfafddba33c6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 15 Sep 2019 21:57:25 -0700 Subject: [PATCH 43/44] updated tests required to reflect changes in dadf5.py --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index b17290e28..7a989197f 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit b17290e28fadeee9e868e551b2df7e5109514cc7 +Subproject commit 7a989197fc4ed1b557f9554618cfc65c9fc4eacf From ac2f366a7248ff14e52fc3ead5c2fe8975e9b67d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 15 Sep 2019 22:25:35 -0700 Subject: [PATCH 44/44] test failed --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 7a989197f..93564092d 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 7a989197fc4ed1b557f9554618cfc65c9fc4eacf +Subproject commit 93564092d20e0c9553245418874ddc3b4484f3dd