From c5d0c7e52ea63eb1e0d15b04897cd15171b6b2af Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 10:43:07 +0200 Subject: [PATCH 01/46] easier to read, more flexible --- python/damask/__init__.py | 4 ++-- python/damask/_result.py | 10 ++++------ python/damask/_vtk.py | 12 ++++++------ 3 files changed, 12 insertions(+), 14 deletions(-) diff --git a/python/damask/__init__.py b/python/damask/__init__.py index 0f343a581..6d65f7cd1 100644 --- a/python/damask/__init__.py +++ b/python/damask/__init__.py @@ -1,9 +1,9 @@ """Tools for pre and post processing of DAMASK simulations.""" -import os as _os +from pathlib import Path as _Path import re as _re name = 'damask' -with open(_os.path.join(_os.path.dirname(__file__),'VERSION')) as _f: +with open(_Path(__file__).parent/_Path('VERSION')) as _f: version = _re.sub(r'^v','',_f.readline().strip()) # make classes directly accessible as damask.Class diff --git a/python/damask/_result.py b/python/damask/_result.py index 73a038c4c..b6ac3782e 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -6,6 +6,7 @@ import os import datetime import xml.etree.ElementTree as ET import xml.dom.minidom +from pathlib import Path from functools import partial import h5py @@ -88,7 +89,7 @@ class Result: 'con_physics': self.con_physics, 'mat_physics': self.mat_physics } - self.fname = os.path.abspath(fname) + self.fname = Path(fname).absolute() self._allow_overwrite = False @@ -1163,7 +1164,7 @@ class Result: 'Dimensions': '{} {} {} {}'.format(*self.grid,np.prod(shape))} data_items[-1].text='{}:{}'.format(os.path.split(self.fname)[1],name) - with open(os.path.splitext(self.fname)[0]+'.xdmf','w') as f: + with open(self.fname.with_suffix('.xdmf').name,'w') as f: f.write(xml.dom.minidom.parseString(ET.tostring(xdmf).decode()).toprettyxml()) @@ -1239,7 +1240,4 @@ class Result: u = self.read_dataset(self.get_dataset_location('u_n' if mode.lower() == 'cell' else 'u_p')) v.add(u,'u') - file_out = '{}_inc{}'.format(os.path.splitext(os.path.basename(self.fname))[0], - inc[3:].zfill(N_digits)) - - v.write(file_out) + v.write('{}_inc{}'.format(self.fname.stem,inc[3:].zfill(N_digits))) diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py index 58b43dc1f..72e224bda 100644 --- a/python/damask/_vtk.py +++ b/python/damask/_vtk.py @@ -1,4 +1,4 @@ -import os +from pathlib import Path import pandas as pd import numpy as np @@ -126,8 +126,8 @@ class VTK: vtkUnstructuredGrid, and vtkPolyData. """ - ext = os.path.splitext(fname)[1] - if ext == '.vtk': + ext = Path(fname).suffix + if ext == '.vtk' or dataset_type: reader = vtk.vtkGenericDataObjectReader() reader.SetFileName(fname) reader.Update() @@ -176,10 +176,10 @@ class VTK: writer = vtk.vtkXMLPolyDataWriter() default_ext = writer.GetDefaultFileExtension() - name, ext = os.path.splitext(fname) + ext = Path(fname).suffix if ext and ext != '.'+default_ext: - raise ValueError('Given extension {} is not .{}'.format(ext,default_ext)) - writer.SetFileName('{}.{}'.format(name,default_ext)) + raise ValueError('Given extension {} does not match default .{}'.format(ext,default_ext)) + writer.SetFileName(str(Path(fname).with_suffix('.'+default_ext))) writer.SetCompressorTypeToZLib() writer.SetDataModeToBinary() writer.SetInputData(self.geom) From c67fbacfc7fb828bf9ac8e462c3b0e3e87510b71 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 11:03:31 +0200 Subject: [PATCH 02/46] higher test coverage - invalid operations - legacy output --- python/tests/test_VTK.py | 34 +++++++++++++++++++++++++--------- 1 file changed, 25 insertions(+), 9 deletions(-) diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py index 627ea4007..8795e7161 100644 --- a/python/tests/test_VTK.py +++ b/python/tests/test_VTK.py @@ -17,18 +17,24 @@ class TestVTK: size = np.random.random(3) + 1.0 origin = np.random.random(3) v = VTK.from_rectilinearGrid(grid,size,origin) - s = v.__repr__() + string = v.__repr__() v.write(os.path.join(tmp_path,'rectilinearGrid')) - v = VTK.from_file(os.path.join(tmp_path,'rectilinearGrid.vtr')) - assert(s == v.__repr__()) + vtr = VTK.from_file(os.path.join(tmp_path,'rectilinearGrid.vtr')) + with open(os.path.join(tmp_path,'rectilinearGrid.vtk'),'w') as f: + f.write(string) + vtk = VTK.from_file(os.path.join(tmp_path,'rectilinearGrid.vtk'),'VTK_rectilinearGrid') + assert(string == vtr.__repr__() == vtk.__repr__()) def test_polyData(self,tmp_path): points = np.random.rand(3,100) v = VTK.from_polyData(points) - s = v.__repr__() + string = v.__repr__() v.write(os.path.join(tmp_path,'polyData')) - v = VTK.from_file(os.path.join(tmp_path,'polyData.vtp')) - assert(s == v.__repr__()) + vtp = VTK.from_file(os.path.join(tmp_path,'polyData.vtp')) + with open(os.path.join(tmp_path,'polyData.vtk'),'w') as f: + f.write(string) + vtk = VTK.from_file(os.path.join(tmp_path,'polyData.vtk'),'polyData') + assert(string == vtp.__repr__() == vtk.__repr__()) @pytest.mark.parametrize('cell_type,n',[ ('VTK_hexahedron',8), @@ -41,7 +47,17 @@ class TestVTK: nodes = np.random.rand(n,3) connectivity = np.random.choice(np.arange(n),n,False).reshape(-1,n) v = VTK.from_unstructuredGrid(nodes,connectivity,cell_type) - s = v.__repr__() + string = v.__repr__() v.write(os.path.join(tmp_path,'unstructuredGrid')) - v = VTK.from_file(os.path.join(tmp_path,'unstructuredGrid.vtu')) - assert(s == v.__repr__()) + vtu = VTK.from_file(os.path.join(tmp_path,'unstructuredGrid.vtu')) + with open(os.path.join(tmp_path,'unstructuredGrid.vtk'),'w') as f: + f.write(string) + vtk = VTK.from_file(os.path.join(tmp_path,'unstructuredGrid.vtk'),'unstructuredgrid') + assert(string == vtu.__repr__() == vtk.__repr__()) + + @pytest.mark.parametrize('name,dataset_type',[('this_file_does_not_exist.vtk',None), + ('this_file_does_not_exist.vtk','vtk'), + ('this_file_does_not_exist.vtx', None)]) + def test_invalid_dataset_type(self,dataset_type,name): + with pytest.raises(TypeError): + VTK.from_file('this_file_does_not_exist.vtk',dataset_type) From 3be1a33820c5ab276c4e487b0276b689732849c0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 13:32:47 +0200 Subject: [PATCH 03/46] easier to read --- installation/symlink_Processing.py | 56 +++++++-------------------- processing/pre/mentat_pbcOnBoxMesh.py | 2 +- processing/pre/mentat_spectralBox.py | 2 +- python/damask/_environment.py | 25 ++++++------ python/damask/solver/_marc.py | 27 ++++++------- 5 files changed, 44 insertions(+), 68 deletions(-) diff --git a/installation/symlink_Processing.py b/installation/symlink_Processing.py index 5972a7f6a..9ee6db99d 100755 --- a/installation/symlink_Processing.py +++ b/installation/symlink_Processing.py @@ -1,56 +1,30 @@ #!/usr/bin/env python3 # Makes postprocessing routines accessible from everywhere. -import os import sys +from pathlib import Path import damask -damaskEnv = damask.Environment() -baseDir = damaskEnv.relPath('processing/') -binDir = damaskEnv.relPath('bin/') +env = damask.Environment() +bin_dir = env.root_dir/Path('bin') -if not os.path.isdir(binDir): - os.mkdir(binDir) +if not bin_dir.exists(): + bin_dir.mkdir() -#define ToDo list -processing_subDirs = ['pre', - 'post', - ] sys.stdout.write('\nsymbolic linking...\n') +for sub_dir in ['pre','post']: + the_dir = env.root_dir/Path('processing')/Path(sub_dir) -for subDir in processing_subDirs: - theDir = os.path.abspath(os.path.join(baseDir,subDir)) - - sys.stdout.write('\n'+binDir+' ->\n'+theDir+damask.util.deemph(' ...')+'\n') - - for theFile in os.listdir(theDir): - theName,theExt = os.path.splitext(theFile) - if theExt in ['.py']: - - src = os.path.abspath(os.path.join(theDir,theFile)) - sym_link = os.path.abspath(os.path.join(binDir,theName)) - - if os.path.lexists(sym_link): - os.remove(sym_link) - output = theName+damask.util.deemph(theExt) - else: - output = damask.util.emph(theName)+damask.util.deemph(theExt) - - sys.stdout.write(damask.util.deemph('... ')+output+'\n') - os.symlink(src,sym_link) + for the_file in the_dir.glob('*.py'): + src = the_dir/the_file + dst = bin_dir/Path(the_file.with_suffix('').name) + dst.unlink(True) + dst.symlink_to(src) sys.stdout.write('\npruning broken links...\n') - -brokenLinks = 0 - -for filename in os.listdir(binDir): - path = os.path.join(binDir,filename) - if os.path.islink(path) and not os.path.exists(path): - sys.stdout.write(' '+damask.util.delete(path)+'\n') - os.remove(path) - brokenLinks += 1 - -sys.stdout.write(('none.' if brokenLinks == 0 else '')+'\n') +for filename in bin_dir.glob('*'): + if not filename.is_file(): + filename.unlink diff --git a/processing/pre/mentat_pbcOnBoxMesh.py b/processing/pre/mentat_pbcOnBoxMesh.py index c171c6ccd..defcbb09b 100755 --- a/processing/pre/mentat_pbcOnBoxMesh.py +++ b/processing/pre/mentat_pbcOnBoxMesh.py @@ -6,7 +6,7 @@ import numpy as np from optparse import OptionParser import damask -sys.path.append(damask.solver.Marc().libraryPath()) +sys.path.append(str(damask.solver.Marc().library_path)) scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) diff --git a/processing/pre/mentat_spectralBox.py b/processing/pre/mentat_spectralBox.py index a61bef57a..e6d138952 100755 --- a/processing/pre/mentat_spectralBox.py +++ b/processing/pre/mentat_spectralBox.py @@ -9,7 +9,7 @@ import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) -sys.path.append(damask.solver.Marc().libraryPath()) +sys.path.append(str(damask.solver.Marc().library_path)) #------------------------------------------------------------------------------------------------- def outMentat(cmd,locals): diff --git a/python/damask/_environment.py b/python/damask/_environment.py index df690314c..7d12050aa 100644 --- a/python/damask/_environment.py +++ b/python/damask/_environment.py @@ -1,11 +1,11 @@ import os +from pathlib import Path import tkinter class Environment: def __init__(self): """Read and provide values of DAMASK configuration.""" - self.options = self._get_options() try: tk = tkinter.Tk() self.screen_width = tk.winfo_screenwidth() @@ -15,17 +15,8 @@ class Environment: self.screen_width = 1024 self.screen_height = 768 - def relPath(self,relative = '.'): - """Return absolute path from path relative to DAMASK root.""" - return os.path.join(self.rootDir(),relative) - - - def rootDir(self): - """Return DAMASK root path.""" - return os.path.normpath(os.path.join(os.path.realpath(__file__),'../../../')) - - - def _get_options(self): + @property + def options(self): options = {} for item in ['DAMASK_NUM_THREADS', 'MSC_ROOT', @@ -34,3 +25,13 @@ class Environment: options[item] = os.environ[item] if item in os.environ else None return options + + @property + def root_dir(self): + """Return DAMASK root path.""" + return Path(__file__).parents[2] + + + # for compatibility + def rootDir(self): + return str(self.root_dir) diff --git a/python/damask/solver/_marc.py b/python/damask/solver/_marc.py index 22e9842f5..9d8b30406 100644 --- a/python/damask/solver/_marc.py +++ b/python/damask/solver/_marc.py @@ -2,6 +2,7 @@ import os import subprocess import shlex import string +from pathlib import Path from .._environment import Environment @@ -25,22 +26,22 @@ class Marc: self.version = -1 -#-------------------------- - def libraryPath(self): + @property + def library_path(self): path_MSC = Environment().options['MSC_ROOT'] - path_lib = '{}/mentat{}/shlib/linux64'.format(path_MSC,self.version) + path_lib = Path('{}/mentat{}/shlib/linux64'.format(path_MSC,self.version)) - return path_lib if os.path.exists(path_lib) else '' + return path_lib if path_lib.is_file() else None -#-------------------------- - def toolsPath(self): + @property + def tools_path(self): path_MSC = Environment().options['MSC_ROOT'] path_tools = '{}/marc{}/tools'.format(path_MSC,self.version) - return path_tools if os.path.exists(path_tools) else '' + return path_tools if path_tools.is_file() else None #-------------------------- @@ -53,21 +54,21 @@ class Marc: ): - damaskEnv = Environment() + env = Environment() - user = os.path.join(damaskEnv.relPath('src'),'DAMASK_marc{}.{}'.format(self.version,'f90' if compile else 'marc')) - if not os.path.isfile(user): + user = env.root_dir/Path('src')/Path('DAMASK_marc{}'.format(version)).with_suffix('.f90' if compile else '.marc') + if not user.is_file(): raise FileNotFoundError("DAMASK4Marc ({}) '{}' not found".format(('source' if compile else 'binary'),user)) # Define options [see Marc Installation and Operation Guide, pp 23] script = 'run_damask_{}mp'.format(optimization) - cmd = os.path.join(self.toolsPath(),script) + \ + cmd = str(self.tools_path/Path(script)) + \ ' -jid ' + model + '_' + job + \ ' -nprocd 1 -autorst 0 -ci n -cr n -dcoup 0 -b no -v no' - if compile: cmd += ' -u ' + user + ' -save y' - else: cmd += ' -prog ' + os.path.splitext(user)[0] + if compile: cmd += ' -u ' + str(user) + ' -save y' + else: cmd += ' -prog ' + str(user.with_suffix('')) print('job submission {} compilation: {}'.format('with' if compile else 'without',user)) if logfile: log = open(logfile, 'w') From 33a8aa4db72ee8bb16ba872599eb620dd6dba8ef Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 13:46:53 +0200 Subject: [PATCH 04/46] don't catch all exceptions and give meaningful meassages --- processing/pre/mentat_pbcOnBoxMesh.py | 23 ++++++++++------------- 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/processing/pre/mentat_pbcOnBoxMesh.py b/processing/pre/mentat_pbcOnBoxMesh.py index defcbb09b..adf0b1aff 100755 --- a/processing/pre/mentat_pbcOnBoxMesh.py +++ b/processing/pre/mentat_pbcOnBoxMesh.py @@ -197,14 +197,15 @@ def add_servoLinks(mfd_data,active=[True,True,True]): # directions on which to if mfd_data[i]['uid'] == 1705: del mfd_data[i] mfd_data.insert(i, links) - + + #-------------------------------------------------------------------------------------------------- # MAIN #-------------------------------------------------------------------------------------------------- + parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ Set up servo linking to achieve periodic boundary conditions for a regular hexahedral mesh. Use *py_connection to operate on model presently opened in MSC.Mentat. - """, version = scriptID) parser.add_option('-p', '--port', @@ -229,10 +230,7 @@ if remote and filenames != []: if filenames == []: filenames = [None] if remote: - try: import py_mentat - except: - damask.util.croak('no valid Mentat release found.') - sys.exit(-1) + import py_mentat damask.util.report(scriptName, 'waiting to connect...') filenames = [os.path.join(tempfile._get_default_tempdir(), next(tempfile._get_candidate_names()) + '.mfd')] @@ -240,14 +238,14 @@ if remote: py_mentat.py_connect('',options.port) py_mentat.py_send('*set_save_formatted on') py_mentat.py_send('*save_as_model "{}" yes'.format(filenames[0])) - py_mentat.py_get_int("nnodes()") # hopefully blocks until file is written - except: - damask.util.croak('failed. try setting Tools/Python/"Run as Separate Process" & "Initiate".') - sys.exit() + py_mentat.py_get_int("nnodes()") + except py_mentat.InputError as err: + damask.util.croak('{}. Try Tools/Python/"Run as Separate Process" & "Initiate".'.format(err)) + sys.exit(-1) damask.util.croak( 'connected...') for name in filenames: - while remote and not os.path.exists(name): time.sleep(0.5) # wait for Mentat to write MFD file + while remote and not os.path.exists(name): time.sleep(0.5) with open( name,'r') if name is not None else sys.stdin as fileIn: damask.util.report(scriptName, name) mfd = parseMFD(fileIn) @@ -257,5 +255,4 @@ for name in filenames: fileOut.write(asMFD(mfd)) if remote: - try: py_mentat.py_send('*open_model "{}"'.format(filenames[0])) - except: damask.util.croak('lost connection on sending open command for "{}".'.format(filenames[0])) + py_mentat.py_send('*open_model "{}"'.format(filenames[0])) From 9b04a45bbd23a6e46307c0d30dbc8f0dfbc5f39d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 13:53:00 +0200 Subject: [PATCH 05/46] bugfix (wrong variable name) --- python/damask/solver/_marc.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/python/damask/solver/_marc.py b/python/damask/solver/_marc.py index 9d8b30406..8bee9efc0 100644 --- a/python/damask/solver/_marc.py +++ b/python/damask/solver/_marc.py @@ -1,4 +1,3 @@ -import os import subprocess import shlex import string @@ -56,7 +55,7 @@ class Marc: env = Environment() - user = env.root_dir/Path('src')/Path('DAMASK_marc{}'.format(version)).with_suffix('.f90' if compile else '.marc') + user = env.root_dir/Path('src/DAMASK_marc{}'.format(self.version)).with_suffix('.f90' if compile else '.marc') if not user.is_file(): raise FileNotFoundError("DAMASK4Marc ({}) '{}' not found".format(('source' if compile else 'binary'),user)) From 18849c3a69b411bbe5b178509f1979a2c805c959 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 14:03:39 +0200 Subject: [PATCH 06/46] support for older python --- installation/symlink_Processing.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/installation/symlink_Processing.py b/installation/symlink_Processing.py index 9ee6db99d..441e9a9ea 100755 --- a/installation/symlink_Processing.py +++ b/installation/symlink_Processing.py @@ -20,7 +20,7 @@ for sub_dir in ['pre','post']: for the_file in the_dir.glob('*.py'): src = the_dir/the_file dst = bin_dir/Path(the_file.with_suffix('').name) - dst.unlink(True) + if dst.is_file(): dst.unlink() # dst.unlink(True) for Python >3.8 dst.symlink_to(src) From 5da1aa49bc1778c8f38e6f555d9f1891fc977198 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 15:17:00 +0200 Subject: [PATCH 07/46] string should be a Path object --- python/damask/solver/_marc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/solver/_marc.py b/python/damask/solver/_marc.py index 8bee9efc0..ce32b17e5 100644 --- a/python/damask/solver/_marc.py +++ b/python/damask/solver/_marc.py @@ -38,7 +38,7 @@ class Marc: def tools_path(self): path_MSC = Environment().options['MSC_ROOT'] - path_tools = '{}/marc{}/tools'.format(path_MSC,self.version) + path_tools = Path('{}/marc{}/tools'.format(path_MSC,self.version)) return path_tools if path_tools.is_file() else None From bb5485927ed74814d4f6c1151e757e7d32a1564f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 15:18:46 +0200 Subject: [PATCH 08/46] names like 2020.2 cannot be converted to int --- python/damask/solver/_marc.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/python/damask/solver/_marc.py b/python/damask/solver/_marc.py index ce32b17e5..0484ba7e6 100644 --- a/python/damask/solver/_marc.py +++ b/python/damask/solver/_marc.py @@ -19,11 +19,7 @@ class Marc: """ self.solver = 'Marc' - try: - self.version = int(version) - except TypeError: - self.version = -1 - + self.version = version @property def library_path(self): From bda555fd1c55fcf46aea4288cd523321bb12b84c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 16:47:03 +0200 Subject: [PATCH 09/46] we are looking for a path, not for a file --- python/damask/solver/_marc.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/damask/solver/_marc.py b/python/damask/solver/_marc.py index 0484ba7e6..dfe45a46d 100644 --- a/python/damask/solver/_marc.py +++ b/python/damask/solver/_marc.py @@ -27,7 +27,7 @@ class Marc: path_MSC = Environment().options['MSC_ROOT'] path_lib = Path('{}/mentat{}/shlib/linux64'.format(path_MSC,self.version)) - return path_lib if path_lib.is_file() else None + return path_lib if path_lib.is_dir() else None @property @@ -36,7 +36,7 @@ class Marc: path_MSC = Environment().options['MSC_ROOT'] path_tools = Path('{}/marc{}/tools'.format(path_MSC,self.version)) - return path_tools if path_tools.is_file() else None + return path_tools if path_tools.is_dir() else None #-------------------------- From a9c61ede69e21c094f304828afbe47b31bdfd29d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 20:24:18 +0200 Subject: [PATCH 10/46] bugfix: should also work if DAMASK_NUM_THREADS is not set --- python/damask/_result.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index b6ac3782e..f68dfbea0 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -1030,14 +1030,17 @@ class Result: Parameters ---------- func : function - Callback function that calculates a new dataset from one or more datasets per HDF5 group. + Callback function that calculates a new dataset from one or + more datasets per HDF5 group. datasets : dictionary - Details of the datasets to be used: label (in HDF5 file) and arg (argument to which the data is parsed in func). + Details of the datasets to be used: label (in HDF5 file) and + arg (argument to which the data is parsed in func). args : dictionary, optional Arguments parsed to func. """ - pool = multiprocessing.Pool(int(Environment().options['DAMASK_NUM_THREADS'])) + num_threads = Environment().options['DAMASK_NUM_THREADS'] + pool = multiprocessing.Pool(int(num_threads) if num_threads is not None else None) lock = multiprocessing.Manager().Lock() groups = self.groups_with_datasets(datasets.values()) From 82326ed8126cb6a7bc6e91a1fe77d6ef1e446044 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 11 Jun 2020 08:21:11 +0200 Subject: [PATCH 11/46] drop support for ping-pong scheme --- src/DAMASK_marc.f90 | 33 +-------------------------------- src/crystallite.f90 | 2 +- src/numerics.f90 | 5 ----- 3 files changed, 2 insertions(+), 38 deletions(-) diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index fa0711b02..dae77fb5d 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -305,36 +305,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & endif ! convergence treatment end - if (usePingPong) then - calcMode(nn,cp_en) = .not. calcMode(nn,cp_en) ! ping pong (calc <--> collect) - if (calcMode(nn,cp_en)) then ! now --- CALC --- - computationMode = CPFEM_CALCRESULTS - if (lastLovl /= lovl) then ! first after ping pong - call debug_reset() ! resets debugging - outdatedFFN1 = .false. - cycleCounter = cycleCounter + 1 - !mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates - !call mesh_build_ipCoordinates() ! update ip coordinates - endif - if (outdatedByNewInc) then - computationMode = ior(computationMode,CPFEM_AGERESULTS) ! calc and age results - outdatedByNewInc = .false. ! reset flag - endif - else ! now --- COLLECT --- - computationMode = CPFEM_COLLECT ! plain collect - if (lastLovl /= lovl .and. & .not. terminallyIll) & - call debug_info() ! first after ping pong reports (meaningful) debugging - if (lastIncConverged) then - computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) ! collect and backup Jacobian after convergence - lastIncConverged = .false. ! reset flag - endif - !do node = 1,theMesh%elem%nNodes - !CPnodeID = mesh_element(4+node,cp_en) - !mesh_node(1:ndeg,CPnodeID) = mesh_node0(1:ndeg,CPnodeID) + numerics_unitlength * dispt(1:ndeg,node) - !enddo - endif - - else ! --- PLAIN MODE --- computationMode = CPFEM_CALCRESULTS ! always calc if (lastLovl /= lovl) then if (.not. terminallyIll) & @@ -353,7 +323,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) ! backup Jacobian after convergence lastIncConverged = .false. ! reset flag endif - endif theTime = cptim ! record current starting time theDelta = timinc ! record current time increment @@ -362,7 +331,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & endif lastLovl = lovl ! record lovl - call CPFEM_general(computationMode,usePingPong,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) + call CPFEM_general(computationMode,.false.,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) d = ddsdde(1:ngens,1:ngens) s = stress(1:ndi+nshear) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 929fec862..c9fef3a51 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -245,7 +245,7 @@ subroutine crystallite_init enddo !$OMP END PARALLEL DO - if(any(plasticState%nonlocal) .and. .not. usePingPong) call IO_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal? + !if(any(plasticState%nonlocal) .and. .not. usePingPong) call IO_error(601) crystallite_partionedFp0 = crystallite_Fp0 crystallite_partionedFi0 = crystallite_Fi0 diff --git a/src/numerics.f90 b/src/numerics.f90 index b0163aee3..77c3397de 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -28,8 +28,6 @@ module numerics numerics_unitlength = 1.0_pReal, & !< determines the physical length of one computational length unit charLength = 1.0_pReal, & !< characteristic length scale for gradient problems residualStiffness = 1.0e-6_pReal !< non-zero residual damage - logical, protected, public :: & - usePingPong = .true. !-------------------------------------------------------------------------------------------------- ! field parameters: @@ -133,8 +131,6 @@ subroutine numerics_init defgradTolerance = IO_floatValue(line,chunkPos,2) case ('ijacostiffness') iJacoStiffness = IO_intValue(line,chunkPos,2) - case ('usepingpong') - usepingpong = IO_intValue(line,chunkPos,2) > 0 case ('unitlength') numerics_unitlength = IO_floatValue(line,chunkPos,2) @@ -221,7 +217,6 @@ subroutine numerics_init ! writing parameters to output write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness - write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength !-------------------------------------------------------------------------------------------------- From b353129ba88000e3c21542097e4380af2bfc947a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 11 Jun 2020 08:36:21 +0200 Subject: [PATCH 12/46] cleaning --- src/CPFEM.f90 | 74 +++------------------------------------------ src/DAMASK_marc.f90 | 4 +-- 2 files changed, 6 insertions(+), 72 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 7791d2b12..7ab6a2ece 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -124,7 +124,7 @@ end subroutine CPFEM_init !-------------------------------------------------------------------------------------------------- !> @brief perform initialization at first call, update variables and call the actual material model !-------------------------------------------------------------------------------------------------- -subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyStress, jacobian) +subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyStress, jacobian) integer(pInt), intent(in) :: elFE, & !< FE element number ip !< integration point number @@ -133,17 +133,14 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt ffn1 !< deformation gradient for t=t1 integer(pInt), intent(in) :: mode !< computation mode 1: regular computation plus aging of results real(pReal), intent(in) :: temperature_inp !< temperature - logical, intent(in) :: parallelExecution !< flag indicating parallel computation of requested IPs real(pReal), dimension(6), intent(out) :: cauchyStress !< stress as 6 vector real(pReal), dimension(6,6), intent(out) :: jacobian !< jacobian as 66 tensor (Consistent tangent dcs/dE) real(pReal) J_inverse, & ! inverse of Jacobian rnd - real(pReal), dimension (3,3) :: Kirchhoff, & ! Piola-Kirchhoff stress - cauchyStress33 ! stress vector + real(pReal), dimension (3,3) :: Kirchhoff ! Piola-Kirchhoff stress real(pReal), dimension (3,3,3,3) :: H_sym, & - H, & - jacobian3333 ! jacobian in Matrix notation + H integer(pInt) elCP, & ! crystal plasticity element number i, j, k, l, m, n, ph, homog, mySource @@ -178,9 +175,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward - !*** collection of FEM input with returning of randomize odd stress and jacobian - !* If no parallel execution is required, there is no need to collect FEM input - if (.not. parallelExecution) then chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP))) case (THERMAL_conduction_ID) chosenThermal1 temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = & @@ -189,38 +183,11 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt materialpoint_F0(1:3,1:3,ip,elCP) = ffn materialpoint_F(1:3,1:3,ip,elCP) = ffn1 - elseif (iand(mode, CPFEM_COLLECT) /= 0_pInt) then - call random_number(rnd) - if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal - CPFEM_cs(1:6,ip,elCP) = rnd * ODD_STRESS - CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6) - chosenThermal2: select case (thermal_type(material_homogenizationAt(elCP))) - case (THERMAL_conduction_ID) chosenThermal2 - temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = & - temperature_inp - end select chosenThermal2 - materialpoint_F0(1:3,1:3,ip,elCP) = ffn - materialpoint_F(1:3,1:3,ip,elCP) = ffn1 - CPFEM_calc_done = .false. - endif - - !*** calculation of stress and jacobian if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then !*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one - validCalculation: if (terminallyIll & - .or. outdatedFFN1 & - .or. any(abs(ffn1 - materialpoint_F(1:3,1:3,ip,elCP)) > defgradTolerance)) then - if (any(abs(ffn1 - materialpoint_F(1:3,1:3,ip,elCP)) > defgradTolerance)) then - if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) then - write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at elFE ip',elFE,ip - write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',& - transpose(materialpoint_F(1:3,1:3,ip,elCP)) - write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',transpose(ffn1) - endif - outdatedFFN1 = .true. - endif + validCalculation: if (terminallyIll .or. outdatedFFN1) then call random_number(rnd) if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd @@ -229,23 +196,12 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt !*** deformation gradient is not outdated else validCalculation updateJaco = mod(cycleCounter,iJacoStiffness) == 0 - !* no parallel computation, so we use just one single elFE and ip for computation - if (.not. parallelExecution) then FEsolving_execElem = elCP FEsolving_execIP = ip if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip call materialpoint_stressAndItsTangent(updateJaco, dt) - !* parallel computation and calulation not yet done - elseif (.not. CPFEM_calc_done) then - if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & - write(6,'(a,i8,a,i8)') '<< CPFEM >> calculation for elements ',FEsolving_execElem(1),& - ' to ',FEsolving_execElem(2) - call materialpoint_stressAndItsTangent(updateJaco, dt) - CPFEM_calc_done = .true. - endif - !* map stress and stiffness (or return odd values if terminally ill) terminalIllness: if (terminallyIll) then @@ -300,28 +256,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt cauchyStress = CPFEM_cs (1:6, ip,elCP) jacobian = CPFEM_dcsdE(1:6,1:6,ip,elCP) - - !*** remember extreme values of stress ... - cauchyStress33 = math_6toSym33(CPFEM_cs(1:6,ip,elCP),weighted=.false.) - if (maxval(cauchyStress33) > debug_stressMax) then - debug_stressMaxLocation = [elCP, ip] - debug_stressMax = maxval(cauchyStress33) - endif - if (minval(cauchyStress33) < debug_stressMin) then - debug_stressMinLocation = [elCP, ip] - debug_stressMin = minval(cauchyStress33) - endif - !*** ... and Jacobian - jacobian3333 = math_66toSym3333(CPFEM_dcsdE(1:6,1:6,ip,elCP),weighted=.false.) - if (maxval(jacobian3333) > debug_jacobianMax) then - debug_jacobianMaxLocation = [elCP, ip] - debug_jacobianMax = maxval(jacobian3333) - endif - if (minval(jacobian3333) < debug_jacobianMin) then - debug_jacobianMinLocation = [elCP, ip] - debug_jacobianMin = minval(jacobian3333) - endif - end subroutine CPFEM_general diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index dae77fb5d..fbff127e9 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -261,7 +261,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & endif !$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc - !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS + !$ call omp_set_num_threads(1) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS if (.not. CPFEM_init_done) call CPFEM_initAll(m(1),nn) @@ -331,7 +331,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & endif lastLovl = lovl ! record lovl - call CPFEM_general(computationMode,.false.,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) + call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) d = ddsdde(1:ngens,1:ngens) s = stress(1:ndi+nshear) From 579ced6a52aab319c3b4630459648438be12500b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 11 Jun 2020 08:44:24 +0200 Subject: [PATCH 13/46] removed global public variables --- src/DAMASK_marc.f90 | 3 --- src/crystallite.f90 | 3 --- src/debug.f90 | 55 +-------------------------------------------- 3 files changed, 1 insertion(+), 60 deletions(-) diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index fbff127e9..a653bcc00 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -307,9 +307,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & computationMode = CPFEM_CALCRESULTS ! always calc if (lastLovl /= lovl) then - if (.not. terminallyIll) & - call debug_info() ! first reports (meaningful) debugging - call debug_reset() ! and resets debugging outdatedFFN1 = .false. cycleCounter = cycleCounter + 1 !mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates diff --git a/src/crystallite.f90 b/src/crystallite.f90 index c9fef3a51..14b624b0a 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -276,9 +276,6 @@ subroutine crystallite_init write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax flush(6) endif - - call debug_info - call debug_reset #endif end subroutine crystallite_init diff --git a/src/debug.f90 b/src/debug.f90 index a0947c410..7564c3037 100644 --- a/src/debug.f90 +++ b/src/debug.f90 @@ -49,26 +49,11 @@ module debug debug_i = 1, & debug_g = 1 - integer, dimension(2), public :: & - debug_stressMaxLocation = 0, & - debug_stressMinLocation = 0, & - debug_jacobianMaxLocation = 0, & - debug_jacobianMinLocation = 0 - - - real(pReal), public :: & - debug_stressMax = -huge(1.0_pReal), & - debug_stressMin = huge(1.0_pReal), & - debug_jacobianMax = -huge(1.0_pReal), & - debug_jacobianMin = huge(1.0_pReal) - #ifdef PETSc character(len=1024), parameter, public :: & PETSCDEBUG = ' -snes_view -snes_monitor ' #endif - public :: debug_init, & - debug_reset, & - debug_info + public :: debug_init contains @@ -230,42 +215,4 @@ subroutine debug_init end subroutine debug_init - -!-------------------------------------------------------------------------------------------------- -!> @brief resets all debug values -!-------------------------------------------------------------------------------------------------- -subroutine debug_reset - - debug_stressMaxLocation = 0 - debug_stressMinLocation = 0 - debug_jacobianMaxLocation = 0 - debug_jacobianMinLocation = 0 - debug_stressMax = -huge(1.0_pReal) - debug_stressMin = huge(1.0_pReal) - debug_jacobianMax = -huge(1.0_pReal) - debug_jacobianMin = huge(1.0_pReal) - -end subroutine debug_reset - - -!-------------------------------------------------------------------------------------------------- -!> @brief writes debug statements to standard out -!-------------------------------------------------------------------------------------------------- -subroutine debug_info - - !$OMP CRITICAL (write2out) - debugOutputCPFEM: if (iand(debug_level(debug_CPFEM),debug_LEVELBASIC) /= 0 & - .and. any(debug_stressMinLocation /= 0) & - .and. any(debug_stressMaxLocation /= 0) ) then - write(6,'(2/,a,/)') ' Extreme values of returned stress and Jacobian' - write(6,'(a39)') ' value el ip' - write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' stress min :', debug_stressMin, debug_stressMinLocation - write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' max :', debug_stressMax, debug_stressMaxLocation - write(6,'(a14,1x,e12.3,1x,i8,1x,i4)') ' Jacobian min :', debug_jacobianMin, debug_jacobianMinLocation - write(6,'(a14,1x,e12.3,1x,i8,1x,i4,/)') ' max :', debug_jacobianMax, debug_jacobianMaxLocation - endif debugOutputCPFEM - !$OMP END CRITICAL (write2out) - -end subroutine debug_info - end module debug From e5c9380bac2be818b6886f01a006def560b7a7bc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 11 Jun 2020 08:51:17 +0200 Subject: [PATCH 14/46] cleaning --- src/CPFEM.f90 | 12 ++++-------- src/DAMASK_marc.f90 | 10 +--------- src/marc/discretization_marc.f90 | 10 ++-------- 3 files changed, 7 insertions(+), 25 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 7ab6a2ece..5acfa9aa1 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -43,7 +43,6 @@ module CPFEM theTime = 0.0_pReal, & !< needs description theDelta = 0.0_pReal logical, public :: & - outdatedFFN1 = .false., & !< needs description lastIncConverged = .false., & !< needs description outdatedByNewInc = .false. !< needs description @@ -68,12 +67,9 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief call (thread safe) all module initializations +!> @brief call all module initializations !-------------------------------------------------------------------------------------------------- -subroutine CPFEM_initAll(el,ip) - - integer(pInt), intent(in) :: el, & !< FE el number - ip !< FE integration point number +subroutine CPFEM_initAll CPFEM_init_done = .true. call DAMASK_interface_init @@ -88,7 +84,7 @@ subroutine CPFEM_initAll(el,ip) call YAML_init call HDF5_utilities_init call results_init(.false.) - call discretization_marc_init(ip, el) + call discretization_marc_init call lattice_init call material_init(.false.) call constitutive_init @@ -187,7 +183,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then !*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one - validCalculation: if (terminallyIll .or. outdatedFFN1) then + validCalculation: if (terminallyIll) then call random_number(rnd) if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index a653bcc00..6de01273d 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -43,9 +43,6 @@ module DAMASK_interface logical, protected, public :: symmetricSolver character(len=*), parameter, public :: INPUTFILEEXTENSION = '.dat' - logical, dimension(:,:), public, allocatable :: & - calcMode !< calculate or collect (ping pong scheme) - public :: & DAMASK_interface_init, & getSolverJobName @@ -263,7 +260,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & !$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc !$ call omp_set_num_threads(1) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS - if (.not. CPFEM_init_done) call CPFEM_initAll(m(1),nn) + if (.not. CPFEM_init_done) call CPFEM_initAll computationMode = 0 ! save initialization value, since it does not result in any calculation if (lovl == 4 ) then ! jacobian requested by marc @@ -277,20 +274,17 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & if (inc == 0) then ! >> start of analysis << lastIncConverged = .false. ! no Jacobian backup outdatedByNewInc = .false. ! no aging of state - calcMode = .false. ! pretend last step was collection lastLovl = lovl ! pretend that this is NOT the first after a lovl change write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> start of analysis..! ',m(1),nn flush(6) else if (inc - theInc > 1) then ! >> restart of broken analysis << lastIncConverged = .false. ! no Jacobian backup outdatedByNewInc = .false. ! no aging of state - calcMode = .true. ! pretend last step was calculation write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> restart of analysis..! ',m(1),nn flush(6) else ! >> just the next inc << lastIncConverged = .true. ! request Jacobian backup outdatedByNewInc = .true. ! request aging of state - calcMode = .true. ! assure last step was calculation write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> new increment..! ',m(1),nn flush(6) endif @@ -299,7 +293,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & outdatedByNewInc = .false. ! no aging of state terminallyIll = .false. cycleCounter = -1 ! first calc step increments this to cycle = 0 - calcMode = .true. ! pretend last step was calculation write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> cutback detected..! ',m(1),nn flush(6) endif ! convergence treatment end @@ -307,7 +300,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & computationMode = CPFEM_CALCRESULTS ! always calc if (lastLovl /= lovl) then - outdatedFFN1 = .false. cycleCounter = cycleCounter + 1 !mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates !call mesh_build_ipCoordinates() ! update ip coordinates diff --git a/src/marc/discretization_marc.f90 b/src/marc/discretization_marc.f90 index 0221ad693..d1026d568 100644 --- a/src/marc/discretization_marc.f90 +++ b/src/marc/discretization_marc.f90 @@ -45,10 +45,8 @@ contains !> @brief initializes the mesh by calling all necessary private routines the mesh module !! Order and routines strongly depend on type of solver !-------------------------------------------------------------------------------------------------- -subroutine discretization_marc_init(ip,el) +subroutine discretization_marc_init - integer, intent(in) :: el, ip - real(pReal), dimension(:,:), allocatable :: & node0_elem, & !< node x,y,z coordinates (initially!) node0_cell @@ -70,7 +68,7 @@ subroutine discretization_marc_init(ip,el) real(pReal), dimension(:,:,:,:),allocatable :: & unscaledNormals - write(6,'(/,a)') ' <<<+- mesh init -+>>>'; flush(6) + write(6,'(/,a)') ' <<<+- discretization_marc init -+>>>'; flush(6) mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh @@ -83,10 +81,6 @@ subroutine discretization_marc_init(ip,el) FEsolving_execElem = [1,nElems] FEsolving_execIP = [1,elem%nIPs] - allocate(calcMode(elem%nIPs,nElems),source=.false.) ! pretend to have collected what first call is asking (F = I) - calcMode(ip,mesh_FEM2DAMASK_elem(el)) = .true. ! first ip,el needs to be already pingponged to "calc" - - allocate(cellNodeDefinition(elem%nNodes-1)) allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems)) call buildCells(connectivity_cell,cellNodeDefinition,& From e952ab71278602a9215b6afe7bf2f3b405aa0ed4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 15 Jun 2020 23:12:49 +0200 Subject: [PATCH 15/46] bugfix do not access unitinialized memory --- src/homogenization.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 62150899e..4fdded94e 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -249,7 +249,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) NiterationHomog = 0 cutBackLooping: do while (.not. terminallyIll .and. & - any(subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > num%subStepMinHomog)) + any(subStep(FEsolving_execIP(1):FEsolving_execIP(2),& + FEsolving_execElem(1):FEsolving_execElem(2)) > num%subStepMinHomog)) !$OMP PARALLEL DO PRIVATE(myNgrains) elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) From 0a9902818c8a3731a81647a152493c57850a7998 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 06:34:12 +0200 Subject: [PATCH 16/46] polishing --- src/DAMASK_marc.f90 | 272 ++++++++++++++++++++++---------------------- 1 file changed, 134 insertions(+), 138 deletions(-) diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 6de01273d..5349608c9 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -172,62 +172,62 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & dispt,coord,ffn,frotn,strechn,eigvn,ffn1,frotn1, & strechn1,eigvn1,ncrd,itel,ndeg,ndm,nnode, & jtype,lclass,ifr,ifu) - use prec - use DAMASK_interface - use numerics - use FEsolving - use debug - use discretization_marc - use homogenization - use CPFEM + use prec + use DAMASK_interface + use numerics + use FEsolving + use debug + use discretization_marc + use homogenization + use CPFEM - implicit none -!$ include "omp_lib.h" ! the openMP function library - integer, intent(in) :: & ! according to MSC.Marc 2012 Manual D - ngens, & !< size of stress-strain law - nn, & !< integration point number - ndi, & !< number of direct components - nshear, & !< number of shear components - ncrd, & !< number of coordinates - itel, & !< dimension of F and R, either 2 or 3 - ndeg, & !< number of degrees of freedom - ndm, & !< not specified in MSC.Marc 2012 Manual D - nnode, & !< number of nodes per element - jtype, & !< element type - ifr, & !< set to 1 if R has been calculated - ifu !< set to 1 if stretch has been calculated - integer, dimension(2), intent(in) :: & ! according to MSC.Marc 2012 Manual D - m, & !< (1) user element number, (2) internal element number - matus, & !< (1) user material identification number, (2) internal material identification number - kcus, & !< (1) layer number, (2) internal layer number - lclass !< (1) element class, (2) 0: displacement, 1: low order Herrmann, 2: high order Herrmann - real(pReal), dimension(*), intent(in) :: & ! has dimension(1) according to MSC.Marc 2012 Manual D, but according to example hypela2.f dimension(*) - e, & !< total elastic strain - de, & !< increment of strain - dt !< increment of state variables - real(pReal), dimension(itel), intent(in) :: & ! according to MSC.Marc 2012 Manual D - strechn, & !< square of principal stretch ratios, lambda(i) at t=n - strechn1 !< square of principal stretch ratios, lambda(i) at t=n+1 - real(pReal), dimension(3,3), intent(in) :: & ! has dimension(itel,*) according to MSC.Marc 2012 Manual D, but we alway assume dimension(3,3) - ffn, & !< deformation gradient at t=n - ffn1 !< deformation gradient at t=n+1 - real(pReal), dimension(itel,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D - frotn, & !< rotation tensor at t=n - eigvn, & !< i principal direction components for j eigenvalues at t=n - frotn1, & !< rotation tensor at t=n+1 - eigvn1 !< i principal direction components for j eigenvalues at t=n+1 - real(pReal), dimension(ndeg,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D - disp, & !< incremental displacements - dispt !< displacements at t=n (at assembly, lovl=4) and displacements at t=n+1 (at stress recovery, lovl=6) - real(pReal), dimension(ncrd,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D - coord !< coordinates - real(pReal), dimension(*), intent(inout) :: & ! according to MSC.Marc 2012 Manual D - t !< state variables (comes in at t=n, must be updated to have state variables at t=n+1) - real(pReal), dimension(ndi+nshear), intent(out) :: & ! has dimension(*) according to MSC.Marc 2012 Manual D, but we need to loop over it - s, & !< stress - should be updated by user - g !< change in stress due to temperature effects - real(pReal), dimension(ngens,ngens), intent(out) :: & ! according to MSC.Marc 2012 Manual D, but according to example hypela2.f dimension(ngens,*) - d !< stress-strain law to be formed + implicit none + include "omp_lib.h" ! the openMP function library + integer, intent(in) :: & ! according to MSC.Marc 2012 Manual D + ngens, & !< size of stress-strain law + nn, & !< integration point number + ndi, & !< number of direct components + nshear, & !< number of shear components + ncrd, & !< number of coordinates + itel, & !< dimension of F and R, either 2 or 3 + ndeg, & !< number of degrees of freedom + ndm, & !< not specified in MSC.Marc 2012 Manual D + nnode, & !< number of nodes per element + jtype, & !< element type + ifr, & !< set to 1 if R has been calculated + ifu !< set to 1 if stretch has been calculated + integer, dimension(2), intent(in) :: & ! according to MSC.Marc 2012 Manual D + m, & !< (1) user element number, (2) internal element number + matus, & !< (1) user material identification number, (2) internal material identification number + kcus, & !< (1) layer number, (2) internal layer number + lclass !< (1) element class, (2) 0: displacement, 1: low order Herrmann, 2: high order Herrmann + real(pReal), dimension(*), intent(in) :: & ! has dimension(1) according to MSC.Marc 2012 Manual D, but according to example hypela2.f dimension(*) + e, & !< total elastic strain + de, & !< increment of strain + dt !< increment of state variables + real(pReal), dimension(itel), intent(in) :: & ! according to MSC.Marc 2012 Manual D + strechn, & !< square of principal stretch ratios, lambda(i) at t=n + strechn1 !< square of principal stretch ratios, lambda(i) at t=n+1 + real(pReal), dimension(3,3), intent(in) :: & ! has dimension(itel,*) according to MSC.Marc 2012 Manual D, but we alway assume dimension(3,3) + ffn, & !< deformation gradient at t=n + ffn1 !< deformation gradient at t=n+1 + real(pReal), dimension(itel,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D + frotn, & !< rotation tensor at t=n + eigvn, & !< i principal direction components for j eigenvalues at t=n + frotn1, & !< rotation tensor at t=n+1 + eigvn1 !< i principal direction components for j eigenvalues at t=n+1 + real(pReal), dimension(ndeg,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D + disp, & !< incremental displacements + dispt !< displacements at t=n (at assembly, lovl=4) and displacements at t=n+1 (at stress recovery, lovl=6) + real(pReal), dimension(ncrd,*), intent(in) :: & ! according to MSC.Marc 2012 Manual D + coord !< coordinates + real(pReal), dimension(*), intent(inout) :: & ! according to MSC.Marc 2012 Manual D + t !< state variables (comes in at t=n, must be updated to have state variables at t=n+1) + real(pReal), dimension(ndi+nshear), intent(out) :: & ! has dimension(*) according to MSC.Marc 2012 Manual D, but we need to loop over it + s, & !< stress - should be updated by user + g !< change in stress due to temperature effects + real(pReal), dimension(ngens,ngens), intent(out) :: & ! according to MSC.Marc 2012 Manual D, but according to example hypela2.f dimension(ngens,*) + d !< stress-strain law to be formed !-------------------------------------------------------------------------------------------------- ! Marc common blocks are in fixed format so they have to be reformated to free format (f90) @@ -236,98 +236,94 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & #include QUOTE(PASTE(./marc/include/concom,Marc4DAMASK)) ! concom is needed for inc, lovl #include QUOTE(PASTE(./marc/include/creeps,Marc4DAMASK)) ! creeps is needed for timinc (time increment) - logical :: cutBack - real(pReal), dimension(6) :: stress - real(pReal), dimension(6,6) :: ddsdde - integer :: computationMode, i, cp_en, node, CPnodeID - !$ integer(4) :: defaultNumThreadsInt !< default value set by Marc + logical :: cutBack + real(pReal), dimension(6) :: stress + real(pReal), dimension(6,6) :: ddsdde + integer :: computationMode, i, cp_en, node, CPnodeID + integer(4) :: defaultNumThreadsInt !< default value set by Marc + + if (iand(debug_level(debug_MARC),debug_LEVELBASIC) /= 0) then + write(6,'(a,/,i8,i8,i2)') ' MSC.MARC information on shape of element(2), IP:', m, nn + write(6,'(a,2(i1))') ' Jacobian: ', ngens,ngens + write(6,'(a,i1)') ' Direct stress: ', ndi + write(6,'(a,i1)') ' Shear stress: ', nshear + write(6,'(a,i2)') ' DoF: ', ndeg + write(6,'(a,i2)') ' Coordinates: ', ncrd + write(6,'(a,i12)') ' Nodes: ', nnode + write(6,'(a,i1)') ' Deformation gradient: ', itel + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n:', & + transpose(ffn) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', & + transpose(ffn1) + endif - if (iand(debug_level(debug_MARC),debug_LEVELBASIC) /= 0) then - write(6,'(a,/,i8,i8,i2)') ' MSC.MARC information on shape of element(2), IP:', m, nn - write(6,'(a,2(i1))') ' Jacobian: ', ngens,ngens - write(6,'(a,i1)') ' Direct stress: ', ndi - write(6,'(a,i1)') ' Shear stress: ', nshear - write(6,'(a,i2)') ' DoF: ', ndeg - write(6,'(a,i2)') ' Coordinates: ', ncrd - write(6,'(a,i12)') ' Nodes: ', nnode - write(6,'(a,i1)') ' Deformation gradient: ', itel - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n:', & - transpose(ffn) - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', & - transpose(ffn1) - endif + defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc + call omp_set_num_threads(1) ! no openMP - !$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc - !$ call omp_set_num_threads(1) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS + if (.not. CPFEM_init_done) call CPFEM_initAll - if (.not. CPFEM_init_done) call CPFEM_initAll + computationMode = 0 ! save initialization value, since it does not result in any calculation + if (lovl == 4 ) then ! jacobian requested by marc + if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) & ! first after cutback + computationMode = CPFEM_RESTOREJACOBIAN + elseif (lovl == 6) then ! stress requested by marc + computationMode = CPFEM_CALCRESULTS ! always calc + cp_en = mesh_FEM2DAMASK_elem(m(1)) + if (cptim > theTime .or. inc /= theInc) then ! reached "convergence" + terminallyIll = .false. + cycleCounter = -1 ! first calc step increments this to cycle = 0 + if (inc == 0) then ! >> start of analysis << + lastIncConverged = .false. + outdatedByNewInc = .false. + lastLovl = lovl ! pretend that this is NOT the first after a lovl change + write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> start of analysis..! ',m(1),nn + else if (inc - theInc > 1) then ! >> restart of broken analysis << + lastIncConverged = .false. + outdatedByNewInc = .false. + write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> restart of analysis..! ',m(1),nn + else ! >> just the next inc << + lastIncConverged = .true. + outdatedByNewInc = .true. + write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> new increment..! ',m(1),nn + endif + else if ( timinc < theDelta ) then ! >> cutBack << + lastIncConverged = .false. + outdatedByNewInc = .false. + terminallyIll = .false. + cycleCounter = -1 ! first calc step increments this to cycle = 0 + write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> cutback detected..! ',m(1),nn + endif ! convergence treatment end + flush(6) - computationMode = 0 ! save initialization value, since it does not result in any calculation - if (lovl == 4 ) then ! jacobian requested by marc - if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) & ! first after cutback - computationMode = CPFEM_RESTOREJACOBIAN - elseif (lovl == 6) then ! stress requested by marc - cp_en = mesh_FEM2DAMASK_elem(m(1)) - if (cptim > theTime .or. inc /= theInc) then ! reached "convergence" - terminallyIll = .false. - cycleCounter = -1 ! first calc step increments this to cycle = 0 - if (inc == 0) then ! >> start of analysis << - lastIncConverged = .false. ! no Jacobian backup - outdatedByNewInc = .false. ! no aging of state - lastLovl = lovl ! pretend that this is NOT the first after a lovl change - write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> start of analysis..! ',m(1),nn - flush(6) - else if (inc - theInc > 1) then ! >> restart of broken analysis << - lastIncConverged = .false. ! no Jacobian backup - outdatedByNewInc = .false. ! no aging of state - write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> restart of analysis..! ',m(1),nn - flush(6) - else ! >> just the next inc << - lastIncConverged = .true. ! request Jacobian backup - outdatedByNewInc = .true. ! request aging of state - write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> new increment..! ',m(1),nn - flush(6) - endif - else if ( timinc < theDelta ) then ! >> cutBack << - lastIncConverged = .false. ! no Jacobian backup - outdatedByNewInc = .false. ! no aging of state - terminallyIll = .false. - cycleCounter = -1 ! first calc step increments this to cycle = 0 - write(6,'(a,i6,1x,i2)') '<< HYPELA2 >> cutback detected..! ',m(1),nn - flush(6) - endif ! convergence treatment end + if (lastLovl /= lovl) then + cycleCounter = cycleCounter + 1 + !mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates + !call mesh_build_ipCoordinates() ! update ip coordinates + endif + if (outdatedByNewInc) then + computationMode = ior(computationMode,CPFEM_AGERESULTS) + outdatedByNewInc = .false. + endif + if (lastIncConverged) then + computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) + lastIncConverged = .false. + endif + theTime = cptim + theDelta = timinc + theInc = inc - computationMode = CPFEM_CALCRESULTS ! always calc - if (lastLovl /= lovl) then - cycleCounter = cycleCounter + 1 - !mesh_cellnode = mesh_build_cellnodes() ! update cell node coordinates - !call mesh_build_ipCoordinates() ! update ip coordinates - endif - if (outdatedByNewInc) then - computationMode = ior(computationMode,CPFEM_AGERESULTS) - outdatedByNewInc = .false. ! reset flag - endif - if (lastIncConverged) then - computationMode = ior(computationMode,CPFEM_BACKUPJACOBIAN) ! backup Jacobian after convergence - lastIncConverged = .false. ! reset flag - endif + endif + lastLovl = lovl - theTime = cptim ! record current starting time - theDelta = timinc ! record current time increment - theInc = inc ! record current increment number + call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) - endif - lastLovl = lovl ! record lovl + d = ddsdde(1:ngens,1:ngens) + s = stress(1:ndi+nshear) + g = 0.0_pReal + if(symmetricSolver) d = 0.5_pReal*(d+transpose(d)) - call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,m(1),nn,stress,ddsdde) - - d = ddsdde(1:ngens,1:ngens) - s = stress(1:ndi+nshear) - g = 0.0_pReal - if(symmetricSolver) d = 0.5_pReal*(d+transpose(d)) - - !$ call omp_set_num_threads(defaultNumThreadsInt) ! reset number of threads to stored default value + call omp_set_num_threads(defaultNumThreadsInt) ! reset number of threads to stored default value end subroutine hypela2 From 54aa5a67ff57ca7db298d06441226dcb3a08810f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 06:41:53 +0200 Subject: [PATCH 17/46] polishing --- src/CPFEM.f90 | 38 +++++++++++++------------------------- 1 file changed, 13 insertions(+), 25 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 5acfa9aa1..8c53eb9a4 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -48,15 +48,12 @@ module CPFEM logical, public, protected :: & CPFEM_init_done = .false. !< remember whether init has been done already - logical, private :: & - CPFEM_calc_done = .false. !< remember whether first ip has already calced the results integer(pInt), parameter, public :: & - CPFEM_COLLECT = 2_pInt**0_pInt, & - CPFEM_CALCRESULTS = 2_pInt**1_pInt, & - CPFEM_AGERESULTS = 2_pInt**2_pInt, & - CPFEM_BACKUPJACOBIAN = 2_pInt**3_pInt, & - CPFEM_RESTOREJACOBIAN = 2_pInt**4_pInt + CPFEM_CALCRESULTS = 2_pInt**0_pInt, & + CPFEM_AGERESULTS = 2_pInt**1_pInt, & + CPFEM_BACKUPJACOBIAN = 2_pInt**2_pInt, & + CPFEM_RESTOREJACOBIAN = 2_pInt**3_pInt public :: & CPFEM_general, & @@ -134,7 +131,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS real(pReal) J_inverse, & ! inverse of Jacobian rnd - real(pReal), dimension (3,3) :: Kirchhoff ! Piola-Kirchhoff stress + real(pReal), dimension (3,3) :: Kirchhoff ! Piola-Kirchhoff stress real(pReal), dimension (3,3,3,3) :: H_sym, & H @@ -142,8 +139,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS i, j, k, l, m, n, ph, homog, mySource logical updateJaco ! flag indicating if Jacobian has to be updated - real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress in case of ping pong dummy cycle - ODD_JACOBIAN = 1e50_pReal !< return value for jacobian in case of ping pong dummy cycle + real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress if terminallyIll + ODD_JACOBIAN = 1e50_pReal !< return value for jacobian if terminallyIll elCP = mesh_FEM2DAMASK_elem(elFE) @@ -167,10 +164,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS if (iand(mode, CPFEM_RESTOREJACOBIAN) /= 0_pInt) & CPFEM_dcsde = CPFEM_dcsde_knownGood - !*** age results if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward - chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP))) case (THERMAL_conduction_ID) chosenThermal1 temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = & @@ -179,26 +174,22 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS materialpoint_F0(1:3,1:3,ip,elCP) = ffn materialpoint_F(1:3,1:3,ip,elCP) = ffn1 - !*** calculation of stress and jacobian if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then - !*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one validCalculation: if (terminallyIll) then call random_number(rnd) if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6) - !*** deformation gradient is not outdated else validCalculation updateJaco = mod(cycleCounter,iJacoStiffness) == 0 - FEsolving_execElem = elCP - FEsolving_execIP = ip - if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & - write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip - call materialpoint_stressAndItsTangent(updateJaco, dt) + FEsolving_execElem = elCP + FEsolving_execIP = ip + if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & + write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip + call materialpoint_stressAndItsTangent(updateJaco, dt) - !* map stress and stiffness (or return odd values if terminally ill) terminalIllness: if (terminallyIll) then call random_number(rnd) @@ -208,7 +199,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS else terminalIllness - ! translate from P to CS + ! translate from P to sigma Kirchhoff = matmul(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP))) J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP)) CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.) @@ -232,7 +223,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS endif terminalIllness endif validCalculation - !* report stress and stiffness if ((iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & .and. ((debug_e == elCP .and. debug_i == ip) & .or. .not. iand(debug_level(debug_CPFEM), debug_levelSelective) /= 0_pInt)) then @@ -245,10 +235,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS endif - !*** warn if stiffness close to zero if (all(abs(CPFEM_dcsdE(1:6,1:6,ip,elCP)) < 1e-10_pReal)) call IO_warning(601,elCP,ip) - !*** copy to output if using commercial FEM solver cauchyStress = CPFEM_cs (1:6, ip,elCP) jacobian = CPFEM_dcsdE(1:6,1:6,ip,elCP) From 06f6e15123a674efce1b5b1c9d2f3adfa1c145f4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 07:02:29 +0200 Subject: [PATCH 18/46] avoid public variables --- src/CPFEM.f90 | 18 ++---------------- src/DAMASK_marc.f90 | 19 +++++++++++++++++-- 2 files changed, 19 insertions(+), 18 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 8c53eb9a4..d285a1e01 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -35,19 +35,9 @@ module CPFEM CPFEM_dcsdE !< Cauchy stress tangent real(pReal), dimension (:,:,:,:), allocatable, private :: & CPFEM_dcsdE_knownGood !< known good tangent - integer(pInt), public :: & - cycleCounter = 0_pInt, & !< needs description - theInc = -1_pInt, & !< needs description - lastLovl = 0_pInt !< lovl in previous call to marc hypela2 - real(pReal), public :: & - theTime = 0.0_pReal, & !< needs description - theDelta = 0.0_pReal - logical, public :: & - lastIncConverged = .false., & !< needs description - outdatedByNewInc = .false. !< needs description - logical, public, protected :: & - CPFEM_init_done = .false. !< remember whether init has been done already + integer(pInt), public :: & + cycleCounter = 0_pInt !< needs description integer(pInt), parameter, public :: & CPFEM_CALCRESULTS = 2_pInt**0_pInt, & @@ -68,7 +58,6 @@ contains !-------------------------------------------------------------------------------------------------- subroutine CPFEM_initAll - CPFEM_init_done = .true. call DAMASK_interface_init call prec_init call IO_init @@ -149,9 +138,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS write(6,'(/,a)') '#############################################' write(6,'(a1,a22,1x,i8,a13)') '#','element', elCP, '#' write(6,'(a1,a22,1x,i8,a13)') '#','ip', ip, '#' - write(6,'(a1,a22,1x,f15.7,a6)') '#','theTime', theTime, '#' - write(6,'(a1,a22,1x,f15.7,a6)') '#','theDelta', theDelta, '#' - write(6,'(a1,a22,1x,i8,a13)') '#','theInc', theInc, '#' write(6,'(a1,a22,1x,i8,a13)') '#','cycleCounter', cycleCounter, '#' write(6,'(a1,a22,1x,i8,a13)') '#','computationMode',mode, '#' if (terminallyIll) & diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 5349608c9..efa054dbf 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -42,6 +42,7 @@ module DAMASK_interface logical, protected, public :: symmetricSolver character(len=*), parameter, public :: INPUTFILEEXTENSION = '.dat' + public :: & DAMASK_interface_init, & @@ -241,6 +242,17 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & real(pReal), dimension(6,6) :: ddsdde integer :: computationMode, i, cp_en, node, CPnodeID integer(4) :: defaultNumThreadsInt !< default value set by Marc + + integer(pInt), save :: & + theInc = -1_pInt, & !< needs description + lastLovl = 0_pInt !< lovl in previous call to marc hypela2 + real(pReal), save :: & + theTime = 0.0_pReal, & !< needs description + theDelta = 0.0_pReal + logical, save :: & + lastIncConverged = .false., & !< needs description + outdatedByNewInc = .false., & !< needs description + CPFEM_init_done = .false. !< remember whether init has been done already if (iand(debug_level(debug_MARC),debug_LEVELBASIC) /= 0) then write(6,'(a,/,i8,i8,i2)') ' MSC.MARC information on shape of element(2), IP:', m, nn @@ -260,14 +272,17 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, & defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc call omp_set_num_threads(1) ! no openMP - if (.not. CPFEM_init_done) call CPFEM_initAll + if (.not. CPFEM_init_done) then + CPFEM_init_done = .true. + call CPFEM_initAll + endif computationMode = 0 ! save initialization value, since it does not result in any calculation if (lovl == 4 ) then ! jacobian requested by marc if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) & ! first after cutback computationMode = CPFEM_RESTOREJACOBIAN elseif (lovl == 6) then ! stress requested by marc - computationMode = CPFEM_CALCRESULTS ! always calc + computationMode = CPFEM_CALCRESULTS cp_en = mesh_FEM2DAMASK_elem(m(1)) if (cptim > theTime .or. inc /= theInc) then ! reached "convergence" terminallyIll = .false. From 753fbb70fd7401745f9a3259b5c124affd5e0df6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 13:27:21 +0200 Subject: [PATCH 19/46] cleaning --- src/discretization.f90 | 2 +- src/mesh/discretization_mesh.f90 | 146 ++++++++++++++----------------- 2 files changed, 67 insertions(+), 81 deletions(-) diff --git a/src/discretization.f90 b/src/discretization.f90 index e9a987912..16f76b158 100644 --- a/src/discretization.f90 +++ b/src/discretization.f90 @@ -49,7 +49,7 @@ subroutine discretization_init(homogenizationAt,microstructureAt,& IPcoords0, & NodeCoords0 integer, optional, intent(in) :: & - sharedNodesBegin + sharedNodesBegin ! index of first node shared among different processes (MPI) write(6,'(/,a)') ' <<<+- discretization init -+>>>'; flush(6) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 0880de115..329205f77 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -8,53 +8,50 @@ module discretization_mesh #include #include #include - use PETScdmplex - use PETScdmda - use PETScis + use PETScdmplex + use PETScdmda + use PETScis - use DAMASK_interface - use IO - use debug - use discretization - use numerics - use FEsolving - use FEM_quadrature - use prec - - implicit none - private - - integer, public, protected :: & - mesh_Nboundaries, & - mesh_NcpElemsGlobal + use DAMASK_interface + use IO + use debug + use discretization + use numerics + use FEsolving + use FEM_quadrature + use prec - integer :: & - mesh_NcpElems !< total number of CP elements in mesh + implicit none + private + + integer, public, protected :: & + mesh_Nboundaries, & + mesh_NcpElemsGlobal + + integer :: & + mesh_NcpElems !< total number of CP elements in mesh !!!! BEGIN DEPRECATED !!!!! - integer, public, protected :: & - mesh_maxNips !< max number of IPs in any CP element + integer, public, protected :: & + mesh_maxNips !< max number of IPs in any CP element !!!! BEGIN DEPRECATED !!!!! - integer, dimension(:,:), allocatable :: & - mesh_element !DEPRECATED + real(pReal), dimension(:,:), allocatable :: & + mesh_ipVolume, & !< volume associated with IP (initially!) + mesh_node0 !< node x,y,z coordinates (initially!) - real(pReal), dimension(:,:), allocatable :: & - mesh_ipVolume, & !< volume associated with IP (initially!) - mesh_node0 !< node x,y,z coordinates (initially!) - - real(pReal), dimension(:,:,:), allocatable :: & - mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) + real(pReal), dimension(:,:,:), allocatable :: & + mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) - DM, public :: geomMesh - - PetscInt, dimension(:), allocatable, public, protected :: & - mesh_boundaries + DM, public :: geomMesh - public :: & - discretization_mesh_init, & - mesh_FEM_build_ipVolumes, & - mesh_FEM_build_ipCoordinates + PetscInt, dimension(:), allocatable, public, protected :: & + mesh_boundaries + + public :: & + discretization_mesh_init, & + mesh_FEM_build_ipVolumes, & + mesh_FEM_build_ipCoordinates contains @@ -67,39 +64,32 @@ subroutine discretization_mesh_init(restart) logical, intent(in) :: restart - integer, dimension(1), parameter:: FE_geomtype = [1] !< geometry type of particular element type - integer, dimension(1) :: FE_Nips !< number of IPs in a specific type of element - integer, allocatable, dimension(:) :: chunkPos integer :: dimPlex, & mesh_Nnodes, & !< total number of nodes in mesh j, l - integer, parameter :: & - mesh_ElemType=1 !< Element type of the mesh (only support homogeneous meshes) PetscSF :: sf DM :: globalMesh PetscInt :: nFaceSets PetscInt, pointer, dimension(:) :: pFaceSets character(len=pStringLen), dimension(:), allocatable :: fileContent - IS :: faceSetIS + IS :: faceSetIS PetscErrorCode :: ierr + integer, dimension(:,:), allocatable :: & + mesh_element + - write(6,'(/,a)') ' <<<+- mesh init -+>>>' - ! read in file call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr) CHKERRQ(ierr) - ! get spatial dimension (2 or 3?) call DMGetDimension(globalMesh,dimPlex,ierr) CHKERRQ(ierr) - write(6,*) 'dimension',dimPlex;flush(6) call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr) CHKERRQ(ierr) ! get number of IDs in face sets (for boundary conditions?) call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr) CHKERRQ(ierr) - write(6,*) 'number of "Face Sets"',mesh_Nboundaries;flush(6) call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) call MPI_Bcast(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr) @@ -142,48 +132,44 @@ subroutine discretization_mesh_init(restart) enddo call DMClone(globalMesh,geomMesh,ierr) CHKERRQ(ierr) - else + else call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr) CHKERRQ(ierr) - endif + endif call DMDestroy(globalMesh,ierr); CHKERRQ(ierr) - + call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,ierr) CHKERRQ(ierr) call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr) CHKERRQ(ierr) - FE_Nips(FE_geomtype(1)) = FEM_nQuadrature(dimPlex,integrationOrder) - mesh_maxNips = FE_Nips(1) - - write(6,*) 'mesh_maxNips',mesh_maxNips + mesh_maxNips = FEM_nQuadrature(dimPlex,integrationOrder) + call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p) call mesh_FEM_build_ipVolumes(dimPlex) - + allocate (mesh_element (4,mesh_NcpElems)); mesh_element = 0 do j = 1, mesh_NcpElems - mesh_element( 1,j) = -1 ! DEPRECATED - mesh_element( 2,j) = mesh_elemType ! elem type + mesh_element( 1,j) = -1 + mesh_element( 2,j) = -1 mesh_element( 3,j) = 1 ! homogenization call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr) CHKERRQ(ierr) - end do + end do + + if (debug_e < 1 .or. debug_e > mesh_NcpElems) call IO_error(602,ext_msg='element') + if (debug_i < 1 .or. debug_i > mesh_maxNips) call IO_error(602,ext_msg='IP') + + FEsolving_execElem = [1,mesh_NcpElems] ! parallel loop bounds set to comprise all DAMASK elements + FEsolving_execIP = [1,mesh_maxNips] - if (debug_e < 1 .or. debug_e > mesh_NcpElems) & - call IO_error(602,ext_msg='element') ! selected element does not exist - if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2,debug_e)))) & - call IO_error(602,ext_msg='IP') ! selected element does not have requested IP - - FEsolving_execElem = [ 1,mesh_NcpElems ] ! parallel loop bounds set to comprise all DAMASK elements - FEsolving_execIP = [1,FE_Nips(FE_geomtype(mesh_element(2,1)))] - allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) call discretization_init(mesh_element(3,:),mesh_element(4,:),& reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & mesh_node0) - + end subroutine discretization_mesh_init @@ -191,23 +177,23 @@ end subroutine discretization_mesh_init !> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume' !-------------------------------------------------------------------------------------------------- subroutine mesh_FEM_build_ipVolumes(dimPlex) - + PetscInt,intent(in):: dimPlex PetscReal :: vol PetscReal, pointer,dimension(:) :: pCent, pNorm PetscInt :: cellStart, cellEnd, cell PetscErrorCode :: ierr - + allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems),source=0.0_pReal) - + call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) allocate(pCent(dimPlex)) allocate(pNorm(dimPlex)) do cell = cellStart, cellEnd-1 - call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr) + call DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr) CHKERRQ(ierr) mesh_ipVolume(:,cell+1) = vol/real(mesh_maxNips,pReal) - enddo + enddo end subroutine mesh_FEM_build_ipVolumes @@ -219,20 +205,20 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints) PetscInt, intent(in) :: dimPlex PetscReal, intent(in) :: qPoints(mesh_maxNips*dimPlex) - + PetscReal, pointer,dimension(:) :: pV0, pCellJ, pInvcellJ PetscReal :: detJ PetscInt :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset PetscErrorCode :: ierr - - + + allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal) - + allocate(pV0(dimPlex)) allocatE(pCellJ(dimPlex**2)) allocatE(pinvCellJ(dimPlex**2)) call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) - do cell = cellStart, cellEnd-1 !< loop over all elements + do cell = cellStart, cellEnd-1 !< loop over all elements call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) CHKERRQ(ierr) qOffset = 0 @@ -246,7 +232,7 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints) enddo qOffset = qOffset + dimPlex enddo - enddo + enddo end subroutine mesh_FEM_build_ipCoordinates From e0d4ee44a340a7adf327732c3fb9e921ee917b2a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 13:59:59 +0200 Subject: [PATCH 20/46] better name --- src/mesh/discretization_mesh.f90 | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 329205f77..f462f7196 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -75,8 +75,9 @@ subroutine discretization_mesh_init(restart) character(len=pStringLen), dimension(:), allocatable :: fileContent IS :: faceSetIS PetscErrorCode :: ierr - integer, dimension(:,:), allocatable :: & - mesh_element + integer, dimension(:), allocatable :: & + homogenizationAt, & + microstructureAt write(6,'(/,a)') ' <<<+- mesh init -+>>>' @@ -149,12 +150,10 @@ subroutine discretization_mesh_init(restart) call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p) call mesh_FEM_build_ipVolumes(dimPlex) - allocate (mesh_element (4,mesh_NcpElems)); mesh_element = 0 + allocate(microstructureAt(mesh_NcpElems)) + allocate(homogenizationAt(mesh_NcpElems),source=1) do j = 1, mesh_NcpElems - mesh_element( 1,j) = -1 - mesh_element( 2,j) = -1 - mesh_element( 3,j) = 1 ! homogenization - call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr) + call DMGetLabelValue(geomMesh,'material',j-1,microstructureAt(j),ierr) CHKERRQ(ierr) end do @@ -166,7 +165,7 @@ subroutine discretization_mesh_init(restart) allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) - call discretization_init(mesh_element(3,:),mesh_element(4,:),& + call discretization_init(microstructureAt,homogenizationAt,& reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & mesh_node0) From d0d9245707a665952042db5a2d2090af8a574939 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Jun 2020 15:02:39 +0200 Subject: [PATCH 21/46] clearer intention --- python/tests/test_Result.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index d7946e5e0..38c0d84cc 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -302,8 +302,8 @@ class TestResult: @pytest.mark.parametrize('allowed',['off','on']) def test_rename(self,default,allowed): - F = default.read_dataset(default.get_dataset_location('F')) if allowed == 'on': + F = default.read_dataset(default.get_dataset_location('F')) default.allow_modification() default.rename('F','new_name') assert np.all(F == default.read_dataset(default.get_dataset_location('new_name'))) From 829896390ca3089d2f9e2e1d7d9d3b2659f2878e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 23:37:24 +0200 Subject: [PATCH 22/46] hopefully not needed any more --- src/crystallite.f90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 929fec862..e85f9e517 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -287,11 +287,9 @@ end subroutine crystallite_init !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- -function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) +function crystallite_stress() logical, dimension(discretization_nIP,discretization_nElem) :: crystallite_stress - real(pReal), intent(in), optional :: & - dummyArgumentToPreventInternalCompilerErrorWithGCC real(pReal) :: & formerSubStep integer :: & From 62384b583677347de327158d431511f5ee6347da Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Jun 2020 23:43:15 +0200 Subject: [PATCH 23/46] bugfix: invalid description/unit --- src/constitutive_plastic_dislotwin.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index 7c7d24ab8..dc40d269e 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -830,7 +830,7 @@ module subroutine plastic_dislotwin_results(instance,group) 'mobile dislocation density','1/m²') case('rho_dip') if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,'rho_dip',& - 'dislocation dipole density''1/m²') + 'dislocation dipole density','1/m²') case('gamma_sl') if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,'gamma_sl',& 'plastic shear','1') From 8fa68101b8e8c13774c4b28be67b3f49dcefd65f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 17 Jun 2020 22:28:29 +0200 Subject: [PATCH 24/46] not needed --- examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config b/examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config index 6920aee2b..2cdd5c58f 100644 --- a/examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config +++ b/examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config @@ -27,7 +27,7 @@ SolidSolutionStrength 0.0 # Strength due to elements in solid solution ### Dislocation glide parameters ### #per family -Nslip 12 0 +Nslip 12 slipburgers 2.72e-10 # Burgers vector of slip system [m] rhoedge0 1.0e12 # Initial edge dislocation density [m/m**3] rhoedgedip0 1.0 # Initial edged dipole dislocation density [m/m**3] From 01f21a9e3baa5aca8fe8eec0421bdd76197e5675 Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 18 Jun 2020 18:04:51 +0200 Subject: [PATCH 25/46] [skip ci] updated version information after successful test of v2.0.3-2664-ge959aaab --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 4a0b49cc6..b25957f75 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2654-g5c544a6e +v2.0.3-2664-ge959aaab From 055fa64f5f0270feb4f41e4384586b6a70928f57 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 19 Jun 2020 12:25:46 +0200 Subject: [PATCH 26/46] better readable --- python/damask/_rotation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index b1b1cd5ad..7e9a03c3a 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -530,7 +530,7 @@ class Rotation: np.zeros(qu.shape[:-1]+(2,))]), np.where(np.abs(q03_s) < 1.0e-8, np.block([np.arctan2( 2.0*qu[...,1:2]*qu[...,2:3],qu[...,1:2]**2-qu[...,2:3]**2), - np.broadcast_to(np.pi,qu.shape[:-1]+(1,)), + np.broadcast_to(np.pi,qu[...,0:1].shape), np.zeros(qu.shape[:-1]+(1,))]), np.block([np.arctan2((-_P*q02+q13)*chi, (-_P*q01-q23)*chi), np.arctan2( 2.0*chi, q03_s-q12_s ), @@ -553,7 +553,7 @@ class Rotation: s = np.sign(qu[...,0:1])/np.sqrt(qu[...,1:2]**2+qu[...,2:3]**2+qu[...,3:4]**2) omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0)) ax = np.where(np.broadcast_to(qu[...,0:1] < 1.0e-8,qu.shape), - np.block([qu[...,1:4],np.broadcast_to(np.pi,qu.shape[:-1]+(1,))]), + np.block([qu[...,1:4],np.broadcast_to(np.pi,qu[...,0:1].shape)]), np.block([qu[...,1:4]*s,omega])) ax[np.isclose(qu[...,0],1.,rtol=0.0)] = [0.0, 0.0, 1.0, 0.0] return ax @@ -564,7 +564,7 @@ class Rotation: with np.errstate(invalid='ignore',divide='ignore'): s = np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True) ro = np.where(np.broadcast_to(np.abs(qu[...,0:1]) < 1.0e-12,qu.shape), - np.block([qu[...,1:2], qu[...,2:3], qu[...,3:4], np.broadcast_to(np.inf,qu.shape[:-1]+(1,))]), + np.block([qu[...,1:2], qu[...,2:3], qu[...,3:4], np.broadcast_to(np.inf,qu[...,0:1].shape)]), np.block([qu[...,1:2]/s,qu[...,2:3]/s,qu[...,3:4]/s, np.tan(np.arccos(np.clip(qu[...,0:1],-1.0,1.0))) ]) From d93ed2bc5c9dfa4a48c4539c74d4dee48087c6db Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Jun 2020 12:20:43 +0200 Subject: [PATCH 27/46] several improvements: - vectorized from_directions - more tests (96% coverage, only random functionality is untested) - updated documentation, folloing numpy standard - inverse operator '~' introduced --- python/damask/_rotation.py | 293 ++++++++++++++++++++++++++-------- python/tests/test_Rotation.py | 41 ++++- 2 files changed, 263 insertions(+), 71 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 7e9a03c3a..a33b4bc51 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -20,8 +20,8 @@ class Rotation: when viewing from the end point of the rotation axis towards the origin. - rotations will be interpreted in the passive sense. - Euler angle triplets are implemented using the Bunge convention, - with the angular ranges as [0, 2π],[0, π],[0, 2π]. - - the rotation angle ω is limited to the interval [0, π]. + with the angular ranges as [0,2π], [0,π], [0,2π]. + - the rotation angle ω is limited to the interval [0,π]. - the real part of a quaternion is positive, Re(q) > 0 - P = -1 (as default). @@ -49,7 +49,8 @@ class Rotation: Parameters ---------- quaternion : numpy.ndarray, optional - Unit quaternion that follows the conventions. Use .from_quaternion to perform a sanity check. + Unit quaternion in positive real hemisphere. + Use .from_quaternion to perform a sanity check. """ self.quaternion = quaternion.copy() @@ -73,7 +74,7 @@ class Rotation: raise NotImplementedError('Support for multiple rotations missing') return '\n'.join([ 'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)), - 'Matrix:\n{}'.format(self.as_matrix()), + 'Matrix:\n{}'.format(np.round(self.as_matrix(),8)), 'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.as_Eulers(degrees=True)), ]) @@ -87,10 +88,6 @@ class Rotation: other : numpy.ndarray or Rotation Vector, second or fourth order tensor, or rotation object that is rotated. - Todo - ---- - Check rotation of 4th order tensor - """ if isinstance(other, Rotation): q_m = self.quaternion[...,0:1] @@ -99,7 +96,7 @@ class Rotation: p_o = other.quaternion[...,1:] q = (q_m*q_o - np.einsum('...i,...i',p_m,p_o).reshape(self.shape+(1,))) p = q_m*p_o + q_o*p_m + _P * np.cross(p_m,p_o) - return self.__class__(np.block([q,p])).standardize() + return self.__class__(np.block([q,p]))._standardize() elif isinstance(other,np.ndarray): if self.shape + (3,) == other.shape: @@ -124,24 +121,24 @@ class Rotation: else: raise TypeError('Cannot rotate {}'.format(type(other))) - def inverse(self): - """In-place inverse rotation/backward rotation.""" - self.quaternion[...,1:] *= -1 - return self - def inversed(self): - """Inverse rotation/backward rotation.""" - return self.copy().inverse() - - - def standardize(self): - """In-place quaternion representation with positive real part.""" + def _standardize(self): + """Standardize (ensure positive real hemisphere).""" self.quaternion[self.quaternion[...,0] < 0.0] *= -1 return self - def standardized(self): - """Quaternion representation with positive real part.""" - return self.copy().standardize() + def inverse(self): + """In-place inverse rotation (backward rotation).""" + self.quaternion[...,1:] *= -1 + return self + + def __invert__(self): + """Inverse rotation (backward rotation).""" + return self.copy().inverse() + + def inversed(self): + """Inverse rotation (backward rotation).""" + return ~ self def misorientation(self,other): @@ -154,7 +151,7 @@ class Rotation: Rotation to which the misorientation is computed. """ - return other*self.inversed() + return other@~self def broadcast_to(self,shape): @@ -169,7 +166,7 @@ class Rotation: return self.__class__(q) - def average(self,other): + def average(self,other): #ToDo: discuss calling for vectors """ Calculate the average rotation. @@ -189,25 +186,31 @@ class Rotation: def as_quaternion(self): """ - Unit quaternion [q, p_1, p_2, p_3]. + Represent as unit quaternion. - Parameters - ---------- - quaternion : bool, optional - return quaternion as DAMASK object. + Returns + ------- + q : numpy.ndarray of shape (...,4) + Unit quaternion in positive real hemisphere: (q_0, q_1, q_2, q_3), |q|=1, q_0 ≥ 0. """ - return self.quaternion + return self.quaternion.copy() def as_Eulers(self, degrees = False): """ - Bunge-Euler angles: (φ_1, ϕ, φ_2). + Represent as Bunge-Euler angles. Parameters ---------- degrees : bool, optional - return angles in degrees. + Return angles in degrees. + + Returns + ------- + phi : numpy.ndarray of shape (...,3) + Bunge-Euler angles: (φ_1, ϕ, φ_2), φ_1 ∈ [0,2π], ϕ ∈ [0,π], φ_2 ∈ [0,2π] + unless degrees == True: φ_1 ∈ [0,360], ϕ ∈ [0,180], φ_2 ∈ [0,360] """ eu = Rotation._qu2eu(self.quaternion) @@ -218,14 +221,21 @@ class Rotation: degrees = False, pair = False): """ - Axis angle representation [n_1, n_2, n_3, ω] unless pair == True: ([n_1, n_2, n_3], ω). + Represent as axis angle pair. Parameters ---------- degrees : bool, optional - return rotation angle in degrees. + Return rotation angle in degrees. Defaults to False. pair : bool, optional - return tuple of axis and angle. + Return tuple of axis and angle. Defaults to False. + + Returns + ------- + axis_angle : numpy.ndarray of shape (...,4) unless pair == True: + tuple containing numpy.ndarray of shapes (...,3) and (...) + Axis angle pair: (n_1, n_2, n_3, ω), |n| = 1 and ω ∈ [0,π] + unless degrees = True: ω ∈ [0,180]. """ ax = Rotation._qu2ax(self.quaternion) @@ -233,29 +243,60 @@ class Rotation: return (ax[...,:3],ax[...,3]) if pair else ax def as_matrix(self): - """Rotation matrix.""" + """ + Represent as rotation matrix. + + Returns + ------- + R : numpy.ndarray of shape (...,3,3) + Rotation matrix R, det(R) = 1, R.T∙R=I. + + """ return Rotation._qu2om(self.quaternion) def as_Rodrigues(self, vector = False): """ - Rodrigues-Frank vector representation [n_1, n_2, n_3, tan(ω/2)] unless vector == True: [n_1, n_2, n_3] * tan(ω/2). + Represent as Rodrigues-Frank vector with separated axis and angle argument. Parameters ---------- vector : bool, optional - return as actual Rodrigues--Frank vector, i.e. rotation axis scaled by tan(ω/2). + Return as actual Rodrigues-Frank vector, i.e. axis + and angle argument are not separated. + + Returns + ------- + rho : numpy.ndarray of shape (...,4) unless vector == True: + numpy.ndarray of shape (...,3) + Rodrigues-Frank vector: [n_1, n_2, n_3, tan(ω/2)], |n| = 1 and ω ∈ [0,π]. """ ro = Rotation._qu2ro(self.quaternion) return ro[...,:3]*ro[...,3] if vector else ro def as_homochoric(self): - """Homochoric vector: (h_1, h_2, h_3).""" + """ + Represent as homochoric vector. + + Returns + ------- + h : numpy.ndarray of shape (...,3) + Homochoric vector: (h_1, h_2, h_3), |h| < 1/2*π^(2/3). + + """ return Rotation._qu2ho(self.quaternion) def as_cubochoric(self): - """Cubochoric vector: (c_1, c_2, c_3).""" + """ + Represent as cubocoric vector. + + Returns + ------- + c : numpy.ndarray of shape (...,3) + Cubochoric vector: (c_1, c_2, c_3), max(c_i) < 1/2*π^(2/3). + + """ return Rotation._qu2cu(self.quaternion) def M(self): # ToDo not sure about the name: as_M or M? we do not have a from_M @@ -275,18 +316,34 @@ class Rotation: # Static constructors. The input data needs to follow the conventions, options allow to # relax the conventions. @staticmethod - def from_quaternion(quaternion, + def from_quaternion(q, accept_homomorph = False, P = -1, - acceptHomomorph = None): + acceptHomomorph = None): # old name (for compatibility) + """ + Initialize from quaternion. + Parameters + ---------- + q: numpy.ndarray of shape (...,4) + Unit quaternion in positive real hemisphere: (q_0, q_1, q_2, q_3), + |q|=1, q_0 ≥ 0. + accept_homomorph: boolean, optional + Allow homomorphic variants, i.e. q_0 < 0 (negative real hemisphere). + Defaults to False. + P: integer ∈ {-1,1}, optional + Convention used. Defaults to -1. + + """ if acceptHomomorph is not None: - accept_homomorph = acceptHomomorph - qu = np.array(quaternion,dtype=float) + accept_homomorph = acceptHomomorph # for compatibility + qu = np.array(q,dtype=float) if qu.shape[:-2:-1] != (4,): raise ValueError('Invalid shape.') + if abs(P) != 1: + raise ValueError('P ∉ {-1,1}') - if P > 0: qu[...,1:4] *= -1 # convert from P=1 to P=-1 + if P == 1: qu[...,1:4] *= -1 if accept_homomorph: qu[qu[...,0] < 0.0] *= -1 else: @@ -298,10 +355,21 @@ class Rotation: return Rotation(qu) @staticmethod - def from_Eulers(eulers, + def from_Eulers(phi, degrees = False): + """ + Initialize from Bunge-Euler angles. - eu = np.array(eulers,dtype=float) + Parameters + ---------- + phi: numpy.ndarray of shape (...,3) + Bunge-Euler angles: (φ_1, ϕ, φ_2), φ_1 ∈ [0,2π], ϕ ∈ [0,π], φ_2 ∈ [0,2π] + unless degrees == True: φ_1 ∈ [0,360], ϕ ∈ [0,180], φ_2 ∈ [0,360]. + degrees: boolean, optional + Bunge-Euler angles are given in degrees. Defaults to False. + + """ + eu = np.array(phi,dtype=float) if eu.shape[:-2:-1] != (3,): raise ValueError('Invalid shape.') @@ -316,12 +384,29 @@ class Rotation: degrees = False, normalise = False, P = -1): + """ + Initialize from Axis angle pair. + Parameters + ---------- + axis_angle: numpy.ndarray of shape (...,4) + Axis angle pair: [n_1, n_2, n_3, ω], |n| = 1 and ω ∈ [0,π] + unless degrees = True: ω ∈ [0,180]. + degrees: boolean, optional + Angle ω is given in degrees. Defaults to False. + normalize: boolean, optional + Allow |n| ≠ 1. Defaults to False. + P: integer ∈ {-1,1}, optional + Convention used. Defaults to -1. + + """ ax = np.array(axis_angle,dtype=float) if ax.shape[:-2:-1] != (4,): raise ValueError('Invalid shape.') + if abs(P) != 1: + raise ValueError('P ∉ {-1,1}') - if P > 0: ax[...,0:3] *= -1 # convert from P=1 to P=-1 + if P == 1: ax[...,0:3] *= -1 if degrees: ax[..., 3] = np.radians(ax[...,3]) if normalise: ax[...,0:3] /= np.linalg.norm(ax[...,0:3],axis=-1) if np.any(ax[...,3] < 0.0) or np.any(ax[...,3] > np.pi): @@ -335,7 +420,19 @@ class Rotation: def from_basis(basis, orthonormal = True, reciprocal = False): + """ + Initialize from tbd. + Parameters + ---------- + basis: numpy.ndarray of shape (...,3,3) + tbd + orthonormal: boolean, optional + tbd. Defaults to True. + reciprocal: boolean, optional + tbd. Defaults to False. + + """ om = np.array(basis,dtype=float) if om.shape[:-3:-1] != (3,3): raise ValueError('Invalid shape.') @@ -356,20 +453,64 @@ class Rotation: return Rotation(Rotation._om2qu(om)) @staticmethod - def from_matrix(om): + def from_directions(hkl,uvw): + """ + Initialize from pair of directions/planes. - return Rotation.from_basis(om) + Parameters + ---------- + hkl: numpy.ndarray of shape (...,3) + Direction parallel to z direction, i.e. (h k l) || (0,0,1). + uvw: numpy.ndarray of shape (...,3) + Direction parallel to x direction, i.e. || (1,0,0). + + """ + hkl_ = hkl/np.linalg.norm(hkl,axis=-1,keepdims=True) + uvw_ = uvw/np.linalg.norm(uvw,axis=-1,keepdims=True) + v_1 = np.block([uvw_,np.cross(hkl_,uvw_),hkl_]).reshape(hkl_.shape+(3,)) + v_2 = np.block([uvw_,np.cross(uvw_,hkl_),hkl_]).reshape(hkl_.shape+(3,)) + R = np.where(np.broadcast_to(np.expand_dims(np.expand_dims(np.linalg.det(v_1)>0,-1),-1),v_1.shape), + v_1,v_2) + return Rotation.from_basis(np.swapaxes(R,axis2=-2,axis1=-1)) @staticmethod - def from_Rodrigues(rodrigues, + def from_matrix(R): + """ + Initialize from rotation matrix. + + Parameters + ---------- + R: numpy.ndarray of shape (...,3,3) + Rotation matrix: det(R) = 1, R.T∙R=I. + + """ + return Rotation.from_basis(R) + + @staticmethod + def from_Rodrigues(rho, normalise = False, P = -1): + """ + Initialize from Rodrigues-Frank vector. - ro = np.array(rodrigues,dtype=float) + Parameters + ---------- + rho: numpy.ndarray of shape (...,4) + Rodrigues-Frank vector (angle separated from axis). + (n_1, n_2, n_3, tan(ω/2)), |n| = 1 and ω ∈ [0,π]. + normalize: boolean, optional + Allow |n| ≠ 1. Defaults to False. + P: integer ∈ {-1,1}, optional + Convention used. Defaults to -1. + + """ + ro = np.array(rho,dtype=float) if ro.shape[:-2:-1] != (4,): raise ValueError('Invalid shape.') + if abs(P) != 1: + raise ValueError('P ∉ {-1,1}') - if P > 0: ro[...,0:3] *= -1 # convert from P=1 to P=-1 + if P == 1: ro[...,0:3] *= -1 if normalise: ro[...,0:3] /= np.linalg.norm(ro[...,0:3],axis=-1) if np.any(ro[...,3] < 0.0): raise ValueError('Rodrigues vector rotation angle not positive.') @@ -379,14 +520,26 @@ class Rotation: return Rotation(Rotation._ro2qu(ro)) @staticmethod - def from_homochoric(homochoric, + def from_homochoric(h, P = -1): + """ + Initialize from homochoric vector. - ho = np.array(homochoric,dtype=float) + Parameters + ---------- + h: numpy.ndarray of shape (...,3) + Homochoric vector: (h_1, h_2, h_3), |h| < (3/4*π)^(1/3). + P: integer ∈ {-1,1}, optional + Convention used. Defaults to -1. + + """ + ho = np.array(h,dtype=float) if ho.shape[:-2:-1] != (3,): raise ValueError('Invalid shape.') + if abs(P) != 1: + raise ValueError('P ∉ {-1,1}') - if P > 0: ho *= -1 # convert from P=1 to P=-1 + if P == 1: ho *= -1 if np.any(np.linalg.norm(ho,axis=-1) >_R1+1e-9): raise ValueError('Homochoric coordinate outside of the sphere.') @@ -394,18 +547,30 @@ class Rotation: return Rotation(Rotation._ho2qu(ho)) @staticmethod - def from_cubochoric(cubochoric, - P = -1): + def from_cubochoric(c, + P = -1): + """ + Initialize from cubochoric vector. - cu = np.array(cubochoric,dtype=float) + Parameters + ---------- + c: numpy.ndarray of shape (...,3) + Cubochoric vector: (c_1, c_2, c_3), max(c_i) < 1/2*π^(2/3). + P: integer ∈ {-1,1}, optional + Convention used. Defaults to -1. + + """ + cu = np.array(c,dtype=float) if cu.shape[:-2:-1] != (3,): raise ValueError('Invalid shape.') + if abs(P) != 1: + raise ValueError('P ∉ {-1,1}') - if np.abs(np.max(cu))>np.pi**(2./3.) * 0.5+1e-9: - raise ValueError('Cubochoric coordinate outside of the cube: {} {} {}.'.format(*cu)) + if np.abs(np.max(cu)) > np.pi**(2./3.) * 0.5+1e-9: + raise ValueError('Cubochoric coordinate outside of the cube.') ho = Rotation._cu2ho(cu) - if P > 0: ho *= -1 # convert from P=1 to P=-1 + if P == 1: ho *= -1 return Rotation(Rotation._ho2qu(ho)) @@ -458,7 +623,7 @@ class Rotation: np.cos(2.0*np.pi*r[...,1])*B, np.sin(2.0*np.pi*r[...,0])*A],axis=-1) - return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q).standardize() + return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q)._standardize() # for compatibility (old names do not follow convention) diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 5a1cd145d..a6b4fea69 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -14,7 +14,7 @@ scatter=1.e-2 @pytest.fixture def default(): - """A set of n random rotations.""" + """A set of n rotations (corner cases and random).""" specials = np.array([ [1.0, 0.0, 0.0, 0.0], #---------------------- @@ -170,7 +170,7 @@ def qu2ax(qu): Modified version of the original formulation, should be numerically more stable """ - if np.isclose(qu[0],1.,rtol=0.0): # set axis to [001] if the angle is 0/360 + if np.isclose(qu[0],1.,rtol=0.0): # set axis to [001] if the angle is 0/360 ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) elif qu[0] > 1.e-8: s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) @@ -534,7 +534,7 @@ def mul(me, other): other_p = other.quaternion[1:] R = me.__class__(np.append(me_q*other_q - np.dot(me_p,other_p), me_q*other_p + other_q*me_p + _P * np.cross(me_p,other_p))) - return R.standardize() + return R._standardize() elif isinstance(other, np.ndarray): if other.shape == (3,): A = me.quaternion[0]**2.0 - np.dot(me.quaternion[1:],me.quaternion[1:]) @@ -554,7 +554,7 @@ def mul(me, other): axes = ((0, 2, 4, 6), (0, 1, 2, 3)) return np.tensordot(RRRR, other, axes) else: - raise ValueError('Can only rotate vectors, 2nd order ternsors, and 4th order tensors') + raise ValueError('Can only rotate vectors, 2nd order tensors, and 4th order tensors') else: raise TypeError('Cannot rotate {}'.format(type(other))) @@ -854,6 +854,15 @@ class TestRotation: rot = Rotation.from_basis(om,False,reciprocal=reciprocal) assert np.isclose(np.linalg.det(rot.as_matrix()),1.0) + def test_directions(self): + hkl = np.array([0.,0.,1.]) + uvw = np.array([1.,0.,0.]) + assert np.allclose(Rotation.from_directions(hkl,uvw).as_matrix(),np.eye(3)) + + @pytest.mark.parametrize('shape',[None,1,(4,4)]) + def test_random(self,shape): + Rotation.from_random(shape) + @pytest.mark.parametrize('function',[Rotation.from_quaternion, Rotation.from_Eulers, Rotation.from_axis_angle, @@ -866,6 +875,16 @@ class TestRotation: with pytest.raises(ValueError): function(invalid_shape) + @pytest.mark.parametrize('fr,to',[(Rotation.from_quaternion,'as_quaternion'), + (Rotation.from_axis_angle,'as_axis_angle'), + (Rotation.from_Rodrigues, 'as_Rodrigues'), + (Rotation.from_homochoric,'as_homochoric'), + (Rotation.from_cubochoric,'as_cubochoric')]) + def test_invalid_P(self,fr,to): + R = Rotation.from_random(np.random.randint(8,32,(3))) # noqa + with pytest.raises(ValueError): + fr(eval('R.{}()'.format(to)),P=-30) + @pytest.mark.parametrize('shape',[None,(3,),(4,2)]) def test_broadcast(self,shape): rot = Rotation.from_random(shape) @@ -932,14 +951,18 @@ class TestRotation: phi_2 = 2*np.pi - phi_1 R_1 = Rotation.from_Eulers(np.array([phi_1,0.,0.])) R_2 = Rotation.from_Eulers(np.array([0.,0.,phi_2])) - assert np.allclose(data,R_2*(R_1*data)) + assert np.allclose(data,R_2@(R_1@data)) + + def test_rotate_inverse(self): + R = Rotation.from_random() + assert np.allclose(np.eye(3),(R.inversed()@R).as_matrix()) @pytest.mark.parametrize('data',[np.random.rand(3), np.random.rand(3,3), np.random.rand(3,3,3,3)]) - def test_rotate_inverse(self,data): + def test_rotate_inverse_array(self,data): R = Rotation.from_random() - assert np.allclose(data,R.inversed()*(R*data)) + assert np.allclose(data,R.inversed()@(R@data)) @pytest.mark.parametrize('data',[np.random.rand(4), np.random.rand(3,2), @@ -956,3 +979,7 @@ class TestRotation: R = Rotation.from_random() with pytest.raises(TypeError): R*data + + def test_misorientation(self): + R = Rotation.from_random() + assert np.allclose(R.misorientation(R).as_matrix(),np.eye(3)) From 3f63a4fdbc56f96dae3f87952b91c9c7d8a667db Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Jun 2020 18:13:34 +0200 Subject: [PATCH 28/46] [skip ci] typo --- python/damask/_rotation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index a33b4bc51..09aeab256 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -289,7 +289,7 @@ class Rotation: def as_cubochoric(self): """ - Represent as cubocoric vector. + Represent as cubochoric vector. Returns ------- From 15b43bcebf9ae932d54e5a623b82177351214569 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Jun 2020 19:57:49 +0200 Subject: [PATCH 29/46] from_directions is not general, removed polishing --- python/damask/_rotation.py | 67 ++++++++++++----------------------- python/tests/test_Rotation.py | 5 --- 2 files changed, 23 insertions(+), 49 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 09aeab256..7c469a2b5 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -325,13 +325,13 @@ class Rotation: Parameters ---------- - q: numpy.ndarray of shape (...,4) + q : numpy.ndarray of shape (...,4) Unit quaternion in positive real hemisphere: (q_0, q_1, q_2, q_3), |q|=1, q_0 ≥ 0. - accept_homomorph: boolean, optional + accept_homomorph : boolean, optional Allow homomorphic variants, i.e. q_0 < 0 (negative real hemisphere). Defaults to False. - P: integer ∈ {-1,1}, optional + P : integer ∈ {-1,1}, optional Convention used. Defaults to -1. """ @@ -362,10 +362,10 @@ class Rotation: Parameters ---------- - phi: numpy.ndarray of shape (...,3) + phi : numpy.ndarray of shape (...,3) Bunge-Euler angles: (φ_1, ϕ, φ_2), φ_1 ∈ [0,2π], ϕ ∈ [0,π], φ_2 ∈ [0,2π] unless degrees == True: φ_1 ∈ [0,360], ϕ ∈ [0,180], φ_2 ∈ [0,360]. - degrees: boolean, optional + degrees : boolean, optional Bunge-Euler angles are given in degrees. Defaults to False. """ @@ -389,14 +389,14 @@ class Rotation: Parameters ---------- - axis_angle: numpy.ndarray of shape (...,4) + axis_angle : numpy.ndarray of shape (...,4) Axis angle pair: [n_1, n_2, n_3, ω], |n| = 1 and ω ∈ [0,π] unless degrees = True: ω ∈ [0,180]. - degrees: boolean, optional + degrees : boolean, optional Angle ω is given in degrees. Defaults to False. normalize: boolean, optional Allow |n| ≠ 1. Defaults to False. - P: integer ∈ {-1,1}, optional + P : integer ∈ {-1,1}, optional Convention used. Defaults to -1. """ @@ -421,16 +421,16 @@ class Rotation: orthonormal = True, reciprocal = False): """ - Initialize from tbd. + Initialize from lattice basis vectors. Parameters ---------- - basis: numpy.ndarray of shape (...,3,3) - tbd - orthonormal: boolean, optional - tbd. Defaults to True. - reciprocal: boolean, optional - tbd. Defaults to False. + basis : numpy.ndarray of shape (...,3,3) + Three lattice basis vectors in three dimensions. + orthonormal : boolean, optional + Basis is strictly orthonormal, i.e. is free of stretch components. Defaults to True. + reciprocal : boolean, optional + Basis vectors are given in reciprocal (instead of real) space. Defaults to False. """ om = np.array(basis,dtype=float) @@ -452,27 +452,6 @@ class Rotation: return Rotation(Rotation._om2qu(om)) - @staticmethod - def from_directions(hkl,uvw): - """ - Initialize from pair of directions/planes. - - Parameters - ---------- - hkl: numpy.ndarray of shape (...,3) - Direction parallel to z direction, i.e. (h k l) || (0,0,1). - uvw: numpy.ndarray of shape (...,3) - Direction parallel to x direction, i.e. || (1,0,0). - - """ - hkl_ = hkl/np.linalg.norm(hkl,axis=-1,keepdims=True) - uvw_ = uvw/np.linalg.norm(uvw,axis=-1,keepdims=True) - v_1 = np.block([uvw_,np.cross(hkl_,uvw_),hkl_]).reshape(hkl_.shape+(3,)) - v_2 = np.block([uvw_,np.cross(uvw_,hkl_),hkl_]).reshape(hkl_.shape+(3,)) - R = np.where(np.broadcast_to(np.expand_dims(np.expand_dims(np.linalg.det(v_1)>0,-1),-1),v_1.shape), - v_1,v_2) - return Rotation.from_basis(np.swapaxes(R,axis2=-2,axis1=-1)) - @staticmethod def from_matrix(R): """ @@ -480,7 +459,7 @@ class Rotation: Parameters ---------- - R: numpy.ndarray of shape (...,3,3) + R : numpy.ndarray of shape (...,3,3) Rotation matrix: det(R) = 1, R.T∙R=I. """ @@ -495,12 +474,12 @@ class Rotation: Parameters ---------- - rho: numpy.ndarray of shape (...,4) + rho : numpy.ndarray of shape (...,4) Rodrigues-Frank vector (angle separated from axis). (n_1, n_2, n_3, tan(ω/2)), |n| = 1 and ω ∈ [0,π]. - normalize: boolean, optional + normalize : boolean, optional Allow |n| ≠ 1. Defaults to False. - P: integer ∈ {-1,1}, optional + P : integer ∈ {-1,1}, optional Convention used. Defaults to -1. """ @@ -527,9 +506,9 @@ class Rotation: Parameters ---------- - h: numpy.ndarray of shape (...,3) + h : numpy.ndarray of shape (...,3) Homochoric vector: (h_1, h_2, h_3), |h| < (3/4*π)^(1/3). - P: integer ∈ {-1,1}, optional + P : integer ∈ {-1,1}, optional Convention used. Defaults to -1. """ @@ -554,9 +533,9 @@ class Rotation: Parameters ---------- - c: numpy.ndarray of shape (...,3) + c : numpy.ndarray of shape (...,3) Cubochoric vector: (c_1, c_2, c_3), max(c_i) < 1/2*π^(2/3). - P: integer ∈ {-1,1}, optional + P : integer ∈ {-1,1}, optional Convention used. Defaults to -1. """ diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index a6b4fea69..933719d6f 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -854,11 +854,6 @@ class TestRotation: rot = Rotation.from_basis(om,False,reciprocal=reciprocal) assert np.isclose(np.linalg.det(rot.as_matrix()),1.0) - def test_directions(self): - hkl = np.array([0.,0.,1.]) - uvw = np.array([1.,0.,0.]) - assert np.allclose(Rotation.from_directions(hkl,uvw).as_matrix(),np.eye(3)) - @pytest.mark.parametrize('shape',[None,1,(4,4)]) def test_random(self,shape): Rotation.from_random(shape) From 4c5939ef233ff4021d0f836f836df3d21ae25c91 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Jun 2020 22:51:00 +0200 Subject: [PATCH 30/46] small polishing --- python/damask/util.py | 14 +++++++------- python/tests/test_Rotation.py | 4 ++-- python/tests/test_util.py | 11 +++++++++++ src/discretization.f90 | 2 +- 4 files changed, 21 insertions(+), 10 deletions(-) create mode 100644 python/tests/test_util.py diff --git a/python/damask/util.py b/python/damask/util.py index cb1d8757d..a3fc71d3a 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -93,7 +93,7 @@ def strikeout(what): def execute(cmd, - streamIn = None, + stream_in = None, wd = './', env = None): """ @@ -103,7 +103,7 @@ def execute(cmd, ---------- cmd : str Command to be executed. - streanIn :, optional + stream_in : file object, optional Input (via pipe) for executed process. wd : str, optional Working directory of process. Defaults to ./ . @@ -119,14 +119,14 @@ def execute(cmd, stderr = subprocess.PIPE, stdin = subprocess.PIPE, env = myEnv) - out,error = [i for i in (process.communicate() if streamIn is None - else process.communicate(streamIn.read().encode('utf-8')))] + stdout, stderr = [i for i in (process.communicate() if stream_in is None + else process.communicate(stream_in.read().encode('utf-8')))] os.chdir(initialPath) - out = out.decode('utf-8').replace('\x08','') - error = error.decode('utf-8').replace('\x08','') + stdout = stdout.decode('utf-8').replace('\x08','') + stderr = stderr.decode('utf-8').replace('\x08','') if process.returncode != 0: raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode)) - return out,error + return stdout, stderr def show_progress(iterable,N_iter=None,prefix='',bar_length=50): diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 933719d6f..86e5fa2a7 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -950,14 +950,14 @@ class TestRotation: def test_rotate_inverse(self): R = Rotation.from_random() - assert np.allclose(np.eye(3),(R.inversed()@R).as_matrix()) + assert np.allclose(np.eye(3),(~R@R).as_matrix()) @pytest.mark.parametrize('data',[np.random.rand(3), np.random.rand(3,3), np.random.rand(3,3,3,3)]) def test_rotate_inverse_array(self,data): R = Rotation.from_random() - assert np.allclose(data,R.inversed()@(R@data)) + assert np.allclose(data,~R@(R@data)) @pytest.mark.parametrize('data',[np.random.rand(4), np.random.rand(3,2), diff --git a/python/tests/test_util.py b/python/tests/test_util.py new file mode 100644 index 000000000..e67478d82 --- /dev/null +++ b/python/tests/test_util.py @@ -0,0 +1,11 @@ +from damask import util + +class TestUtil: + + def test_execute_direct(self): + out,err = util.execute('echo test') + assert out=='test\n' and err=='' + + def test_execute_env(self): + out,err = util.execute('sh -c "echo $test_for_execute"',env={'test_for_execute':'test'}) + assert out=='test\n' and err=='' diff --git a/src/discretization.f90 b/src/discretization.f90 index 16f76b158..e52784c6d 100644 --- a/src/discretization.f90 +++ b/src/discretization.f90 @@ -49,7 +49,7 @@ subroutine discretization_init(homogenizationAt,microstructureAt,& IPcoords0, & NodeCoords0 integer, optional, intent(in) :: & - sharedNodesBegin ! index of first node shared among different processes (MPI) + sharedNodesBegin !< index of first node shared among different processes (MPI) write(6,'(/,a)') ' <<<+- discretization init -+>>>'; flush(6) From d4efadb333f864b4ebef645373cdbeccb2356e83 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 21 Jun 2020 10:03:52 +0200 Subject: [PATCH 31/46] should be availabe outside of this module --- src/quaternions.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index f5f5276e8..c337bd681 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -101,6 +101,8 @@ module quaternions assignment(=), & conjg, aimag, & log, exp, & + abs, dot_product, & + inverse, & real contains From c6a5bb8a3b873a66ba968339d76c754bc63c0d4d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 21 Jun 2020 10:04:45 +0200 Subject: [PATCH 32/46] is 2020 --- python/damask/_rotation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 7c469a2b5..f433b2b37 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -615,7 +615,7 @@ class Rotation: #################################################################################################### # Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations #################################################################################################### -# Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH +# Copyright (c) 2017-2020, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH # Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University # All rights reserved. # From deaf055f8876a55c6bdfee0e0089cb8bfefbadb0 Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 22 Jun 2020 15:59:31 +0200 Subject: [PATCH 33/46] [skip ci] updated version information after successful test of v2.0.3-2680-g17212b27 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index b25957f75..d041f794f 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2664-ge959aaab +v2.0.3-2680-g17212b27 From 4314ec1f3778e5a1acd979761dee384682425bfb Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 24 Jun 2020 11:02:15 -0400 Subject: [PATCH 34/46] [skip ci] started to replace .format() with f-strings --- python/damask/_asciitable.py | 4 ++-- python/damask/grid_filters.py | 4 ++-- python/damask/util.py | 15 +++++++-------- 3 files changed, 11 insertions(+), 12 deletions(-) diff --git a/python/damask/_asciitable.py b/python/damask/_asciitable.py index 8410b23a1..9d762369a 100644 --- a/python/damask/_asciitable.py +++ b/python/damask/_asciitable.py @@ -157,7 +157,7 @@ class ASCIItable(): def head_write(self, header = True): """Write current header information (info + labels).""" - head = ['{}\theader'.format(len(self.info)+self.__IO__['labeled'])] if header else [] + head = [f"{len(self.info)+self.__IO__['labeled']}\theader"] if header else [] head.append(self.info) if self.__IO__['labeled']: head.append('\t'.join(map(self._quote,self.tags))) @@ -209,7 +209,7 @@ class ASCIItable(): labelList.append(tags[id]) else: label = tags[id][2:] # get label - while id < len(tags) and tags[id] == '{}_{}'.format(dim,label): # check successors + while id < len(tags) and tags[id] == f'{dim}_{label}': # check successors id += 1 # next label... dim += 1 # ...should be one higher dimension labelList.append(label) # reached end --> store diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 76ab20872..6c3f24cff 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -233,7 +233,7 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True): origin[_np.where(grid==1)] = 0.0 if grid.prod() != len(coord0): - raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid)) + raise ValueError('Data count {len(coord0)} does not match grid {grid}.') start = origin + delta*.5 end = origin - delta*.5 + size @@ -384,7 +384,7 @@ def node_coord0_gridSizeOrigin(coord0,ordered=True): origin = mincorner if (grid+1).prod() != len(coord0): - raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid)) + raise ValueError('Data count {len(coord0)} does not match grid {grid}.') atol = _np.max(size)*5e-2 if not (_np.allclose(coords[0],_np.linspace(mincorner[0],maxcorner[0],grid[0]+1),atol=atol) and \ diff --git a/python/damask/util.py b/python/damask/util.py index a3fc71d3a..f67df461a 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -125,7 +125,7 @@ def execute(cmd, stdout = stdout.decode('utf-8').replace('\x08','') stderr = stderr.decode('utf-8').replace('\x08','') if process.returncode != 0: - raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode)) + raise RuntimeError(f'{cmd} failed with returncode {process.returncode}') return stdout, stderr @@ -223,21 +223,20 @@ class _ProgressBar: self.start_time = datetime.datetime.now() self.last_fraction = 0.0 - sys.stderr.write('{} {} 0% ETA n/a'.format(self.prefix, '░'*self.bar_length)) + sys.stderr.write(f"{self.prefix} {'░'*self.bar_length} 0% ETA n/a") sys.stderr.flush() def update(self,iteration): fraction = (iteration+1) / self.total + filled_length = int(self.bar_length * fraction) - if int(self.bar_length * fraction) > int(self.bar_length * self.last_fraction): + if filled_length > int(self.bar_length * self.last_fraction): + bar = '█' * filled_length + '░' * (self.bar_length - filled_length) delta_time = datetime.datetime.now() - self.start_time remaining_time = (self.total - (iteration+1)) * delta_time / (iteration+1) remaining_time -= datetime.timedelta(microseconds=remaining_time.microseconds) # remove μs - - filled_length = int(self.bar_length * fraction) - bar = '█' * filled_length + '░' * (self.bar_length - filled_length) - sys.stderr.write('\r{} {} {:>4.0%} ETA {}'.format(self.prefix, bar, fraction, remaining_time)) + sys.stderr.write('\r{self.prefix} {bar} {fraction:>4.0%} ETA {remaining_time}') sys.stderr.flush() self.last_fraction = fraction @@ -249,7 +248,7 @@ class _ProgressBar: class bcolors: """ - ASCII Colors. + ASCII colors. https://svn.blender.org/svnroot/bf-blender/trunk/blender/build_files/scons/tools/bcolors.py https://stackoverflow.com/questions/287871 From a71b8aea5caffa0fed3b01dd950c25dc72cc058d Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 24 Jun 2020 17:29:13 +0200 Subject: [PATCH 35/46] [skip ci] updated version information after successful test of v2.0.3-2692-g3d93a5ff --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index d041f794f..d94a70ab1 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2680-g17212b27 +v2.0.3-2692-g3d93a5ff From ff858fd4c82188ca8b8e070030ac4885e9f62151 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 24 Jun 2020 20:13:09 +0200 Subject: [PATCH 36/46] [skip ci] corrected help string for "srepr" --- python/damask/util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/util.py b/python/damask/util.py index a3fc71d3a..6c8fabf90 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -27,7 +27,7 @@ __all__=[ #################################################################################################### def srepr(arg,glue = '\n'): r""" - Join arguments as individual lines. + Join arguments with glue string. Parameters ---------- From f8f433e826eab368b11f9af2e228a4afcaab682d Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 24 Jun 2020 12:05:12 -0400 Subject: [PATCH 37/46] introduced f-strings --- python/damask/_colormaps.py | 23 ++++++------ python/damask/_geom.py | 73 +++++++++++++++++++------------------ python/damask/_rotation.py | 2 +- python/damask/_vtk.py | 10 ++--- 4 files changed, 54 insertions(+), 54 deletions(-) diff --git a/python/damask/_colormaps.py b/python/damask/_colormaps.py index 2117b0b57..94416ffcd 100644 --- a/python/damask/_colormaps.py +++ b/python/damask/_colormaps.py @@ -450,7 +450,7 @@ class Colormap: def rad_diff(a,b): return abs(a[2]-b[2]) - + def adjust_hue(Msh_sat, Msh_unsat): """If saturation of one of the two colors is too less than the other, hue of the less.""" if Msh_sat[0] >= Msh_unsat[0]: @@ -502,7 +502,7 @@ class Colormap: [RGB] colormap for use in paraview or gmsh, or as raw string, or array. Arguments: name, format, steps, crop. - Format is one of (paraview, gmsh, raw, list). + Format is one of (paraview, gmsh, gom, raw, list). Crop selects a (sub)range in [-1.0,1.0]. Generates sequential map if one limiting color is either white or black, diverging map otherwise. @@ -511,23 +511,22 @@ class Colormap: frac = 0.5*(np.array(crop) + 1.0) # rescale crop range to fractions colors = [self.color(float(i)/(steps-1)*(frac[1]-frac[0])+frac[0]).express_as(model).color for i in range(steps)] if format == 'paraview': - colormap = ['[\n {{\n "ColorSpace": "RGB", "Name": "{}", "DefaultMap": true,\n "RGBPoints" : ['.format(name)] \ - + [' {:4d},{:8.6f},{:8.6f},{:8.6f},'.format(i,color[0],color[1],color[2],) \ - for i,color in enumerate(colors[:-1])] \ - + [' {:4d},{:8.6f},{:8.6f},{:8.6f} '.format(len(colors),colors[-1][0],colors[-1][1],colors[-1][2],)] \ + colormap = [f'[\n {{\n "ColorSpace": "RGB", "Name": "{name}", "DefaultMap": true,\n "RGBPoints" : ['] \ + + [f' {i:4d},{color[0]:8.6f},{color[1]:8.6f},{color[2]:8.6f}{"," if i+1 Date: Wed, 24 Jun 2020 14:18:06 -0400 Subject: [PATCH 38/46] shapes init copes with integers instead of strict tuples; introduced f-strings --- python/damask/_table.py | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/python/damask/_table.py b/python/damask/_table.py index d180e5914..c2d304138 100644 --- a/python/damask/_table.py +++ b/python/damask/_table.py @@ -4,6 +4,7 @@ import pandas as pd import numpy as np from . import version +from . import util class Table: """Store spreadsheet-like data.""" @@ -24,7 +25,7 @@ class Table: """ self.comments = [] if comments is None else [c for c in comments] self.data = pd.DataFrame(data=data) - self.shapes = shapes + self.shapes = { k:(v,) if isinstance(v,(np.integer,int)) else v for k,v in shapes.items() } self._label_condensed() @@ -33,7 +34,7 @@ class Table: labels = [] for label,shape in self.shapes.items(): size = int(np.prod(shape)) - labels += ['{}{}'.format('' if size == 1 else '{}_'.format(i+1),label) for i in range(size)] + labels += [('' if size == 1 else f'{i+1}_')+label for i in range(size)] self.data.columns = labels @@ -47,8 +48,7 @@ class Table: def _add_comment(self,label,shape,info): if info is not None: - c = '{}{}: {}'.format(label,' '+str(shape) if np.prod(shape,dtype=int) > 1 else '',info) - self.comments.append(c) + self.comments.append(f'{label}{" "+str(shape) if np.prod(shape,dtype=int) > 1 else ""}: {info}') @staticmethod @@ -136,7 +136,7 @@ class Table: content = f.readlines() - comments = ['table.py:from_ang v {}'.format(version)] + comments = [f'table.py:from_ang v {version}'] for line in content: if line.startswith('#'): comments.append(line.strip()) @@ -145,7 +145,7 @@ class Table: data = np.loadtxt(content) for c in range(data.shape[1]-10): - shapes['n/a_{}'.format(c+1)] = (1,) + shapes[f'n/a_{c+1}'] = (1,) return Table(data,shapes,comments) @@ -251,8 +251,7 @@ class Table: """ self.data.rename(columns={label_old:label_new},inplace=True) - c = '{} => {}{}'.format(label_old,label_new,'' if info is None else ': {}'.format(info)) - self.comments.append(c) + self.comments.append(f'{label_old} => {label_new}'+('' if info is None else f': {info}')) self.shapes = {(label if label != label_old else label_new):self.shapes[label] for label in self.shapes} @@ -271,7 +270,7 @@ class Table: self._label_flat() self.data.sort_values(labels,axis=0,inplace=True,ascending=ascending) self._label_condensed() - self.comments.append('sorted by [{}]'.format(', '.join(labels))) + self.comments.append(f'sorted by [{", ".join(labels)}]') def append(self,other): @@ -328,18 +327,18 @@ class Table: labels = [] for l in [x for x in self.data.columns if not (x in seen or seen.add(x))]: if self.shapes[l] == (1,): - labels.append('{}'.format(l)) + labels.append(f'{l}') elif len(self.shapes[l]) == 1: - labels += ['{}_{}'.format(i+1,l) \ + labels += [f'{i+1}_{l}' \ for i in range(self.shapes[l][0])] else: - labels += ['{}:{}_{}'.format('x'.join([str(d) for d in self.shapes[l]]),i+1,l) \ + labels += [f'{util.srepr(self.shapes[l],"x")}:{i+1}_{l}' \ for i in range(np.prod(self.shapes[l]))] if new_style: - header = ['# {}'.format(comment) for comment in self.comments] + header = [f'# {comment}' for comment in self.comments] else: - header = ['{} header'.format(len(self.comments)+1)] \ + header = [f'{len(self.comments)+1} header'] \ + self.comments \ try: From 14d3b7e66d456d2e63a3c57048fb6d8e2a7c77f9 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 24 Jun 2020 15:34:51 -0400 Subject: [PATCH 39/46] more f-stringing --- python/damask/_colormaps.py | 5 +- python/damask/_geom.py | 13 +--- python/damask/_lattice.py | 12 ++-- python/damask/_result.py | 113 +++++++++++++++--------------- python/damask/_table.py | 2 +- python/damask/_test.py | 78 +++++++++------------ python/damask/config/material.py | 91 ++++++++++++------------ python/damask/solver/_marc.py | 18 ++--- python/damask/util.py | 2 +- python/tests/test_Geom.py | 17 ++--- python/tests/test_Lattice.py | 2 +- python/tests/test_Orientation.py | 2 +- python/tests/test_Result.py | 8 +-- python/tests/test_Rotation.py | 4 +- python/tests/test_grid_filters.py | 8 +-- 15 files changed, 180 insertions(+), 195 deletions(-) diff --git a/python/damask/_colormaps.py b/python/damask/_colormaps.py index 94416ffcd..6c6f82604 100644 --- a/python/damask/_colormaps.py +++ b/python/damask/_colormaps.py @@ -1,4 +1,5 @@ import numpy as np +from . import util class Color: """Color representation in and conversion between different color-spaces.""" @@ -525,8 +526,8 @@ class Colormap: colormap = [ f'1 1 {name}' + f' 9 {name}' + ' 0 1 0 3 0 0 -1 9 \\ 0 0 0 255 255 255 0 0 255 ' - + '30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 ' + str(len(colors)) - + ' '.join([' 0 {} 255 1'.format(' '.join([str(int(x*255.0)) for x in color])) for color in reversed(colors)])] + + f'30 NO_UNIT 1 1 64 64 64 255 1 0 0 0 0 0 0 3 0 {len(colors)}' + + ' '.join([f' 0 {util.srepr((255*np.array(c)).astype(int)," ")} 255 1' for c in reversed(colors)])] elif format == 'raw': colormap = ['\t'.join(map(str,color)) for color in colors] diff --git a/python/damask/_geom.py b/python/damask/_geom.py index 7deedee5e..5a3e6d47f 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -257,11 +257,11 @@ class Geom: def get_header(self): """Return the full header (grid, size, origin, homogenization, comments).""" - header = ['{} header'.format(len(self.comments)+4)] + self.comments + header = [f'{len(self.comments)+4} header'] + self.comments header.append('grid a {} b {} c {}'.format(*self.get_grid())) header.append('size x {} y {} z {}'.format(*self.get_size())) header.append('origin x {} y {} z {}'.format(*self.get_origin())) - header.append('homogenization {}'.format(self.get_homogenization())) + header.append(f'homogenization {self.get_homogenization()}') return header @@ -374,7 +374,6 @@ class Geom: else: microstructure = microstructure.reshape(grid) - #comments = 'geom.py:from_Laguerre_tessellation v{}'.format(version) return Geom(microstructure+1,size,homogenization=1) @@ -399,7 +398,6 @@ class Geom: KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds) devNull,microstructure = KDTree.query(coords) - #comments = 'geom.py:from_Voronoi_tessellation v{}'.format(version) return Geom(microstructure.reshape(grid)+1,size,homogenization=1) @@ -524,7 +522,6 @@ class Geom: if 'x' in directions: ms = np.concatenate([ms,ms[limits[0]:limits[1]:-1,:,:]],0) - #self.add_comments('geom.py:mirror v{}'.format(version) return self.update(ms,rescale=True) @@ -538,7 +535,6 @@ class Geom: number of grid points in x,y,z direction. """ - #self.add_comments('geom.py:scale v{}'.format(version) return self.update( ndimage.interpolation.zoom( self.microstructure, @@ -565,7 +561,6 @@ class Geom: unique, inverse = np.unique(arr, return_inverse=True) return unique[np.argmax(np.bincount(inverse))] - #self.add_comments('geom.py:clean v{}'.format(version) return self.update(ndimage.filters.generic_filter( self.microstructure, mostFrequent, @@ -580,7 +575,6 @@ class Geom: for i, oldID in enumerate(np.unique(self.microstructure)): renumbered = np.where(self.microstructure == oldID, i+1, renumbered) - #self.add_comments('geom.py:renumber v{}'.format(version) return self.update(renumbered) @@ -615,7 +609,6 @@ class Geom: origin = self.origin-(np.asarray(microstructure_in.shape)-self.grid)*.5 * self.size/self.grid - #self.add_comments('geom.py:rotate v{}'.format(version) return self.update(microstructure_in,origin=origin,rescale=True) @@ -647,7 +640,6 @@ class Geom: canvas[l[0]:r[0],l[1]:r[1],l[2]:r[2]] = self.microstructure[L[0]:R[0],L[1]:R[1],L[2]:R[2]] - #self.add_comments('geom.py:canvas v{}'.format(version) return self.update(canvas,origin=self.origin+offset*self.size/self.grid,rescale=True) @@ -667,5 +659,4 @@ class Geom: for from_ms,to_ms in zip(from_microstructure,to_microstructure): substituted[self.microstructure==from_ms] = to_ms - #self.add_comments('geom.py:substitute v{}'.format(version) return self.update(substituted) diff --git a/python/damask/_lattice.py b/python/damask/_lattice.py index 0f6cffd04..2a9c51b9b 100644 --- a/python/damask/_lattice.py +++ b/python/damask/_lattice.py @@ -26,7 +26,7 @@ class Symmetry: """ if symmetry is not None and symmetry.lower() not in Symmetry.lattices: - raise KeyError('Symmetry/crystal system "{}" is unknown'.format(symmetry)) + raise KeyError(f'Symmetry/crystal system "{symmetry}" is unknown') self.lattice = symmetry.lower() if isinstance(symmetry,str) else symmetry @@ -40,7 +40,7 @@ class Symmetry: def __repr__(self): """Readable string.""" - return '{}'.format(self.lattice) + return f'{self.lattice}' def __eq__(self, other): @@ -348,7 +348,7 @@ class Lattice: def __repr__(self): """Report basic lattice information.""" - return 'Bravais lattice {} ({} symmetry)'.format(self.lattice,self.symmetry) + return f'Bravais lattice {self.lattice} ({self.symmetry} symmetry)' # Kurdjomov--Sachs orientation relationship for fcc <-> bcc transformation @@ -613,17 +613,17 @@ class Lattice: try: relationship = models[model] except KeyError : - raise KeyError('Orientation relationship "{}" is unknown'.format(model)) + raise KeyError(f'Orientation relationship "{model}" is unknown') if self.lattice not in relationship['mapping']: - raise ValueError('Relationship "{}" not supported for lattice "{}"'.format(model,self.lattice)) + raise ValueError(f'Relationship "{model}" not supported for lattice "{self.lattice}"') r = {'lattice':Lattice((set(relationship['mapping'])-{self.lattice}).pop()), # target lattice 'rotations':[] } myPlane_id = relationship['mapping'][self.lattice] otherPlane_id = (myPlane_id+1)%2 - myDir_id = myPlane_id +2 + myDir_id = myPlane_id +2 otherDir_id = otherPlane_id +2 for miller in np.hstack((relationship['planes'],relationship['directions'])): diff --git a/python/damask/_result.py b/python/damask/_result.py index e1846cba9..354c1c3f7 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -50,8 +50,7 @@ class Result: self.version_minor = f.attrs['DADF5-minor'] if self.version_major != 0 or not 2 <= self.version_minor <= 6: - raise TypeError('Unsupported DADF5 version {}.{} '.format(self.version_major, - self.version_minor)) + raise TypeError(f'Unsupported DADF5 version {self.version_major}.{self.version_minor}') self.structured = 'grid' in f['geometry'].attrs.keys() @@ -107,7 +106,7 @@ class Result: self.pick('increments',all_selected_increments) in_between = '' if len(all_selected_increments) < 3 else \ - ''.join(['\n{}\n ...\n'.format(inc) for inc in all_selected_increments[1:-2]]) + ''.join([f'\n{inc}\n ...\n' for inc in all_selected_increments[1:-2]]) return util.srepr(first + in_between + last) @@ -137,7 +136,7 @@ class Result: if what == 'increments': choice = [c if isinstance(c,str) and c.startswith('inc') else - 'inc{}'.format(c) for c in choice] + f'inc{c}' for c in choice] elif what == 'times': what = 'increments' if choice == ['*']: @@ -268,7 +267,7 @@ class Result: def rename(self,name_old,name_new): """ Rename datasets. - + Parameters ---------- name_old : str @@ -412,21 +411,19 @@ class Result: message = '' with h5py.File(self.fname,'r') as f: for i in self.iterate('increments'): - message += '\n{} ({}s)\n'.format(i,self.times[self.increments.index(i)]) + message += f'\n{i} ({self.times[self.increments.index(i)]}s)\n' for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']): for oo in self.iterate(o): - message += ' {}\n'.format(oo) + message += f' {oo}\n' for pp in self.iterate(p): - message += ' {}\n'.format(pp) + message += f' {pp}\n' group = '/'.join([i,o[:-1],oo,pp]) # o[:-1]: plural/singular issue for d in f[group].keys(): try: dataset = f['/'.join([group,d])] + unit = f" / {dataset.attrs['Unit'].decode()}" if 'Unit' in dataset.attrs else '' description = dataset.attrs['Description'].decode() - if 'Unit' in dataset.attrs: - message += ' {} / ({}): {}\n'.format(d,dataset.attrs['Unit'].decode(),description) - else: - message += ' {}: {}\n'.format(d,description) + message += f' {d}{unit}: {description}\n' except KeyError: pass return message @@ -528,10 +525,10 @@ class Result: def _add_absolute(x): return { 'data': np.abs(x['data']), - 'label': '|{}|'.format(x['label']), + 'label': f'|{x["label"]}|', 'meta': { 'Unit': x['meta']['Unit'], - 'Description': 'Absolute value of {} ({})'.format(x['label'],x['meta']['Description']), + 'Description': f"Absolute value of {x['label']} ({x['meta']['Description']})", 'Creator': inspect.stack()[0][3][1:] } } @@ -552,14 +549,14 @@ class Result: def _add_calculation(**kwargs): formula = kwargs['formula'] for d in re.findall(r'#(.*?)#',formula): - formula = formula.replace('#{}#'.format(d),"kwargs['{}']['data']".format(d)) + formula = formula.replace(f'#{d}#',f"kwargs['{d}']['data']") return { 'data': eval(formula), 'label': kwargs['label'], 'meta': { 'Unit': kwargs['unit'], - 'Description': '{} (formula: {})'.format(kwargs['description'],kwargs['formula']), + 'Description': f"{kwargs['description']} (formula: {kwargs['formula']})", 'Creator': inspect.stack()[0][3][1:] } } @@ -596,8 +593,9 @@ class Result: 'label': 'sigma', 'meta': { 'Unit': P['meta']['Unit'], - 'Description': 'Cauchy stress calculated from {} ({}) and {} ({})'\ - .format(P['label'],P['meta']['Description'],F['label'],F['meta']['Description']), + 'Description': "Cauchy stress calculated " + f"from {P['label']} ({P['meta']['Description']})" + f" and {F['label']} ({F['meta']['Description']})", 'Creator': inspect.stack()[0][3][1:] } } @@ -620,10 +618,10 @@ class Result: def _add_determinant(T): return { 'data': np.linalg.det(T['data']), - 'label': 'det({})'.format(T['label']), + 'label': f"det({T['label']})", 'meta': { 'Unit': T['meta']['Unit'], - 'Description': 'Determinant of tensor {} ({})'.format(T['label'],T['meta']['Description']), + 'Description': f"Determinant of tensor {T['label']} ({T['meta']['Description']})", 'Creator': inspect.stack()[0][3][1:] } } @@ -644,10 +642,10 @@ class Result: def _add_deviator(T): return { 'data': mechanics.deviatoric_part(T['data']), - 'label': 's_{}'.format(T['label']), + 'label': f"s_{T['label']}", 'meta': { 'Unit': T['meta']['Unit'], - 'Description': 'Deviator of tensor {} ({})'.format(T['label'],T['meta']['Description']), + 'Description': f"Deviator of tensor {T['label']} ({T['meta']['Description']})", 'Creator': inspect.stack()[0][3][1:] } } @@ -675,10 +673,10 @@ class Result: return { 'data': mechanics.eigenvalues(T_sym['data'])[:,p], - 'label': 'lambda_{}({})'.format(eigenvalue,T_sym['label']), + 'label': f"lambda_{eigenvalue}({T_sym['label']})", 'meta' : { 'Unit': T_sym['meta']['Unit'], - 'Description': '{} eigenvalue of {} ({})'.format(label,T_sym['label'],T_sym['meta']['Description']), + 'Description': f"{label} eigenvalue of {T_sym['label']} ({T_sym['meta']['Description']})", 'Creator': inspect.stack()[0][3][1:] } } @@ -708,14 +706,14 @@ class Result: print('p',eigenvalue) return { 'data': mechanics.eigenvectors(T_sym['data'])[:,p], - 'label': 'v_{}({})'.format(eigenvalue,T_sym['label']), + 'label': f"v_{eigenvalue}({T_sym['label']})", 'meta' : { 'Unit': '1', - 'Description': 'Eigenvector corresponding to {} eigenvalue of {} ({})'\ - .format(label,T_sym['label'],T_sym['meta']['Description']), + 'Description': f"Eigenvector corresponding to {label} eigenvalue" + f" of {T_sym['label']} ({T_sym['meta']['Description']})", 'Creator': inspect.stack()[0][3][1:] } - } + } def add_eigenvector(self,T_sym,eigenvalue='max'): """ Add eigenvector of symmetric tensor. @@ -774,10 +772,10 @@ class Result: def _add_maximum_shear(T_sym): return { 'data': mechanics.maximum_shear(T_sym['data']), - 'label': 'max_shear({})'.format(T_sym['label']), + 'label': f"max_shear({T_sym['label']})", 'meta': { 'Unit': T_sym['meta']['Unit'], - 'Description': 'Maximum shear component of {} ({})'.format(T_sym['label'],T_sym['meta']['Description']), + 'Description': f"Maximum shear component of {T_sym['label']} ({T_sym['meta']['Description']})", 'Creator': inspect.stack()[0][3][1:] } } @@ -800,11 +798,11 @@ class Result: 'stress' return { - 'data': mechanics.Mises_strain(T_sym['data']) if t=='strain' else mechanics.Mises_stress(T_sym['data']), - 'label': '{}_vM'.format(T_sym['label']), + 'data': (mechanics.Mises_strain if t=='strain' else mechanics.Mises_stress)(T_sym['data']), + 'label': f"{T_sym['label']}_vM", 'meta': { 'Unit': T_sym['meta']['Unit'], - 'Description': 'Mises equivalent {} of {} ({})'.format(t,T_sym['label'],T_sym['meta']['Description']), + 'Description': f"Mises equivalent {t} of {T_sym['label']} ({T_sym['meta']['Description']})", 'Creator': inspect.stack()[0][3][1:] } } @@ -837,10 +835,10 @@ class Result: return { 'data': np.linalg.norm(x['data'],ord=o,axis=axis,keepdims=True), - 'label': '|{}|_{}'.format(x['label'],o), + 'label': f"|{x['label']}|_{o}", 'meta': { 'Unit': x['meta']['Unit'], - 'Description': '{}-norm of {} {} ({})'.format(o,t,x['label'],x['meta']['Description']), + 'Description': f"{o}-norm of {t} {x['label']} ({x['meta']['Description']})", 'Creator': inspect.stack()[0][3][1:] } } @@ -866,19 +864,20 @@ class Result: 'label': 'S', 'meta': { 'Unit': P['meta']['Unit'], - 'Description': '2. Piola-Kirchhoff stress calculated from {} ({}) and {} ({})'\ - .format(P['label'],P['meta']['Description'],F['label'],F['meta']['Description']), + 'Description': "2. Piola-Kirchhoff stress calculated " + f"from {P['label']} ({P['meta']['Description']})" + f" and {F['label']} ({F['meta']['Description']})", 'Creator': inspect.stack()[0][3][1:] } } def add_PK2(self,P='P',F='F'): """ - Add second Piola-Kirchhoff calculated from first Piola-Kirchhoff stress and deformation gradient. + Add second Piola-Kirchhoff stress calculated from first Piola-Kirchhoff stress and deformation gradient. Parameters ---------- P : str, optional - Label first Piola-Kirchhoff stress dataset. Defaults to ‘P’. + Label of first Piola-Kirchhoff stress dataset. Defaults to ‘P’. F : str, optional Label of deformation gradient dataset. Defaults to ‘F’. @@ -928,10 +927,10 @@ class Result: def _add_rotational_part(F): return { 'data': mechanics.rotational_part(F['data']), - 'label': 'R({})'.format(F['label']), + 'label': f"R({F['label']})", 'meta': { 'Unit': F['meta']['Unit'], - 'Description': 'Rotational part of {} ({})'.format(F['label'],F['meta']['Description']), + 'Description': f"Rotational part of {F['label']} ({F['meta']['Description']})", 'Creator': inspect.stack()[0][3][1:] } } @@ -952,10 +951,10 @@ class Result: def _add_spherical(T): return { 'data': mechanics.spherical_part(T['data']), - 'label': 'p_{}'.format(T['label']), + 'label': f"p_{T['label']}", 'meta': { 'Unit': T['meta']['Unit'], - 'Description': 'Spherical component of tensor {} ({})'.format(T['label'],T['meta']['Description']), + 'Description': f"Spherical component of tensor {T['label']} ({T['meta']['Description']})", 'Creator': inspect.stack()[0][3][1:] } } @@ -976,10 +975,10 @@ class Result: def _add_strain_tensor(F,t,m): return { 'data': mechanics.strain_tensor(F['data'],t,m), - 'label': 'epsilon_{}^{}({})'.format(t,m,F['label']), + 'label': f"epsilon_{t}^{m}({F['label']})", 'meta': { 'Unit': F['meta']['Unit'], - 'Description': 'Strain tensor of {} ({})'.format(F['label'],F['meta']['Description']), + 'Description': f"Strain tensor of {F['label']} ({F['meta']['Description']})", 'Creator': inspect.stack()[0][3][1:] } } @@ -1006,11 +1005,11 @@ class Result: @staticmethod def _add_stretch_tensor(F,t): return { - 'data': mechanics.left_stretch(F['data']) if t == 'V' else mechanics.right_stretch(F['data']), - 'label': '{}({})'.format(t,F['label']), + 'data': (mechanics.left_stretch if t.upper() == 'V' else mechanics.right_stretch)(F['data']), + 'label': f"{t}({F['label']})", 'meta': { 'Unit': F['meta']['Unit'], - 'Description': '{} stretch tensor of {} ({})'.format('Left' if t == 'V' else 'Right', + 'Description': '{} stretch tensor of {} ({})'.format('Left' if t.upper() == 'V' else 'Right', F['label'],F['meta']['Description']), 'Creator': inspect.stack()[0][3][1:] } @@ -1046,7 +1045,7 @@ class Result: r = func(**datasets_in,**args) return [group,r] except Exception as err: - print('Error during calculation: {}.'.format(err)) + print(f'Error during calculation: {err}.') return None @@ -1091,11 +1090,11 @@ class Result: for l,v in result[1]['meta'].items(): dataset.attrs[l]=v.encode() - creator = 'damask.Result.{} v{}'.format(dataset.attrs['Creator'].decode(),version) + creator = f"damask.Result.{dataset.attrs['Creator'].decode()} v{version}" dataset.attrs['Creator'] = creator.encode() except (OSError,RuntimeError) as err: - print('Could not add dataset: {}.'.format(err)) + print(f'Could not add dataset: {err}.') lock.release() pool.close() @@ -1128,7 +1127,7 @@ class Result: time_data = ET.SubElement(time, 'DataItem') time_data.attrib={'Format': 'XML', 'NumberType': 'Float', - 'Dimensions': '{}'.format(len(self.times))} + 'Dimensions': f'{len(self.times)}'} time_data.text = ' '.join(map(str,self.times)) attributes = [] @@ -1169,7 +1168,7 @@ class Result: data_items[-1].attrib={'Format': 'HDF', 'Precision': '8', 'Dimensions': '{} {} {} 3'.format(*(self.grid+1))} - data_items[-1].text='{}:/{}/geometry/u_n'.format(os.path.split(self.fname)[1],inc) + data_items[-1].text=f'{os.path.split(self.fname)[1]}:/{inc}/geometry/u_n' for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']): for oo in getattr(self,o): @@ -1184,15 +1183,15 @@ class Result: if (shape not in [(1,), (3,), (3,3)]) or dtype != np.float64: continue attributes.append(ET.SubElement(grid, 'Attribute')) - attributes[-1].attrib={'Name': '{}'.format(name.split('/',2)[2]), + attributes[-1].attrib={'Name': name.split('/',2)[2], 'Center': 'Cell', 'AttributeType': 'Tensor'} data_items.append(ET.SubElement(attributes[-1], 'DataItem')) data_items[-1].attrib={'Format': 'HDF', 'NumberType': 'Float', - 'Precision': '{}'.format(prec), + 'Precision': f'{prec}', 'Dimensions': '{} {} {} {}'.format(*self.grid,np.prod(shape))} - data_items[-1].text='{}:{}'.format(os.path.split(self.fname)[1],name) + data_items[-1].text=f'{os.path.split(self.fname)[1]}:{name}' with open(self.fname.with_suffix('.xdmf').name,'w') as f: f.write(xml.dom.minidom.parseString(ET.tostring(xdmf).decode()).toprettyxml()) @@ -1270,4 +1269,4 @@ class Result: u = self.read_dataset(self.get_dataset_location('u_n' if mode.lower() == 'cell' else 'u_p')) v.add(u,'u') - v.write('{}_inc{}'.format(self.fname.stem,inc[3:].zfill(N_digits))) + v.write(f'{self.fname.stem}_inc{inc[3:].zfill(N_digits)}') diff --git a/python/damask/_table.py b/python/damask/_table.py index c2d304138..8c09e2b78 100644 --- a/python/damask/_table.py +++ b/python/damask/_table.py @@ -136,7 +136,7 @@ class Table: content = f.readlines() - comments = [f'table.py:from_ang v {version}'] + comments = [f'table.py:from_ang v{version}'] for line in content: if line.startswith('#'): comments.append(line.strip()) diff --git a/python/damask/_test.py b/python/damask/_test.py index a46eedac8..9c4c7ad3d 100644 --- a/python/damask/_test.py +++ b/python/damask/_test.py @@ -53,7 +53,7 @@ class Test: self.dirBase = os.path.dirname(os.path.realpath(sys.modules[self.__class__.__module__].__file__)) self.parser = OptionParser(option_class=damask.extendableOption, - description = '{} (Test class version: {})'.format(self.description,damask.version), + description = f'{self.description} (Test class version: {damask.version})', usage = './test.py [options]') self.parser.add_option("-k", "--keep", action = "store_true", @@ -93,7 +93,7 @@ class Test: for variant,object in enumerate(self.variants): name = self.variantName(variant) if self.options.show: - logging.critical('{}: {}'.format(variant+1,name)) + logging.critical(f'{variant+1}: {name}') elif self.options.select is not None \ and not (name in self.options.select or str(variant+1) in self.options.select): pass @@ -106,12 +106,12 @@ class Test: self.postprocess(variant) if self.options.update: - if self.update(variant) != 0: logging.critical('update for "{}" failed.'.format(name)) + if self.update(variant) != 0: logging.critical(f'update for "{name}" failed.') elif not (self.options.accept or self.compare(variant)): # no update, do comparison return variant+1 # return culprit except Exception as e: - logging.critical('exception during variant execution: "{}"'.format(str(e))) + logging.critical(f'exception during variant execution: "{e}"') return variant+1 # return culprit return 0 @@ -124,13 +124,13 @@ class Test: try: shutil.rmtree(self.dirCurrent()) except FileNotFoundError: - logging.warning('removal of directory "{}" not possible...'.format(self.dirCurrent())) + logging.warning(f'removal of directory "{self.dirCurrent()}" not possible...') try: os.mkdir(self.dirCurrent()) return True except FileExistsError: - logging.critical('creation of directory "{}" failed.'.format(self.dirCurrent())) + logging.critical(f'creation of directory "{self.dirCurrent()}" failed.') return False def prepareAll(self): @@ -211,7 +211,7 @@ class Test: try: shutil.copy2(source,target) except FileNotFoundError: - logging.critical('error copying {} to {}'.format(source,target)) + logging.critical(f'error copying {source} to {target}') raise FileNotFoundError @@ -222,7 +222,7 @@ class Test: try: shutil.copy2(self.fileInReference(f),self.fileInCurrent(targetfiles[i])) except FileNotFoundError: - logging.critical('Reference2Current: Unable to copy file "{}"'.format(f)) + logging.critical(f'Reference2Current: Unable to copy file "{f}"') raise FileNotFoundError @@ -235,7 +235,7 @@ class Test: shutil.copy2(os.path.join(source,f),self.fileInCurrent(targetfiles[i])) except FileNotFoundError: logging.error(os.path.join(source,f)) - logging.critical('Base2Current: Unable to copy file "{}"'.format(f)) + logging.critical(f'Base2Current: Unable to copy file "{f}"') raise FileNotFoundError @@ -246,7 +246,7 @@ class Test: try: shutil.copy2(self.fileInCurrent(f),self.fileInReference(targetfiles[i])) except FileNotFoundError: - logging.critical('Current2Reference: Unable to copy file "{}"'.format(f)) + logging.critical(f'Current2Reference: Unable to copy file "{f}"') raise FileNotFoundError @@ -257,7 +257,7 @@ class Test: try: shutil.copy2(self.fileInProof(f),self.fileInCurrent(targetfiles[i])) except FileNotFoundError: - logging.critical('Proof2Current: Unable to copy file "{}"'.format(f)) + logging.critical(f'Proof2Current: Unable to copy file "{f}"') raise FileNotFoundError @@ -267,7 +267,7 @@ class Test: try: shutil.copy2(self.fileInReference(f),self.fileInCurrent(targetfiles[i])) except FileNotFoundError: - logging.critical('Current2Current: Unable to copy file "{}"'.format(f)) + logging.critical(f'Current2Current: Unable to copy file "{f}"') raise FileNotFoundError @@ -302,9 +302,7 @@ class Test: max_loc=np.argmax(abs(refArrayNonZero[curArray.nonzero()]/curArray[curArray.nonzero()]-1.)) refArrayNonZero = refArrayNonZero[curArray.nonzero()] curArray = curArray[curArray.nonzero()] - print(' ********\n * maximum relative error {} between {} and {}\n ********'.format(max_err, - refArrayNonZero[max_loc], - curArray[max_loc])) + print(f' ********\n * maximum relative error {max_err} between {refArrayNonZero[max_loc]} and {curArray[max_loc]}\n ********') return max_err else: raise Exception('mismatch in array size to compare') @@ -350,7 +348,7 @@ class Test: for i in range(dataLength): if headings0[i]['shape'] != headings1[i]['shape']: - raise Exception('shape mismatch between {} and {} '.format(headings0[i]['label'],headings1[i]['label'])) + raise Exception(f"shape mismatch between {headings0[i]['label']} and {headings1[i]['label']}") shape[i] = headings0[i]['shape'] for j in range(np.shape(shape[i])[0]): length[i] *= shape[i][j] @@ -358,9 +356,7 @@ class Test: for j in range(np.shape(normShape[i])[0]): normLength[i] *= normShape[i][j] else: - raise Exception('trying to compare {} with {} normed by {} data sets'.format(len(headings0), - len(headings1), - len(normHeadings))) + raise Exception(f'trying to compare {len(headings0)} with {len(headings1)} normed by {len(normHeadings)} data sets') table0 = damask.ASCIItable(name=file0,readonly=True) table0.head_read() @@ -372,11 +368,11 @@ class Test: key1 = ('1_' if length[i]>1 else '') + headings1[i]['label'] normKey = ('1_' if normLength[i]>1 else '') + normHeadings[i]['label'] if key0 not in table0.labels(raw = True): - raise Exception('column "{}" not found in first table...\n'.format(key0)) + raise Exception(f'column "{key0}" not found in first table...') elif key1 not in table1.labels(raw = True): - raise Exception('column "{}" not found in second table...\n'.format(key1)) + raise Exception(f'column "{key1}" not found in second table...') elif normKey not in table0.labels(raw = True): - raise Exception('column "{}" not found in first table...\n'.format(normKey)) + raise Exception(f'column "{normKey}" not found in first table...') else: column[0][i] = table0.label_index(key0) column[1][i] = table1.label_index(key1) @@ -404,9 +400,9 @@ class Test: norm[i] = [1.0 for j in range(line0-len(skipLines))] absTol[i] = True if perLine: - logging.warning('At least one norm of "{}" in first table is 0.0, using absolute tolerance'.format(headings0[i]['label'])) + logging.warning(f"At least one norm of \"{headings0[i]['label']}\" in first table is 0.0, using absolute tolerance") else: - logging.warning('Maximum norm of "{}" in first table is 0.0, using absolute tolerance'.format(headings0[i]['label'])) + logging.warning(f"Maximum norm of \"{headings0[i]['label']}\" in first table is 0.0, using absolute tolerance") line1 = 0 while table1.data_read(): # read next data line of ASCII table @@ -418,18 +414,14 @@ class Test: norm[i][line1-len(skipLines)]) line1 +=1 - if (line0 != line1): raise Exception('found {} lines in first table but {} in second table'.format(line0,line1)) + if (line0 != line1): raise Exception(f'found {line0} lines in first table but {line1} in second table') logging.info(' ********') for i in range(dataLength): if absTol[i]: - logging.info(' * maximum absolute error {} between {} and {}'.format(maxError[i], - headings0[i]['label'], - headings1[i]['label'])) + logging.info(f" * maximum absolute error {maxError[i]} between {headings0[i]['label']} and {headings1[i]['label']}") else: - logging.info(' * maximum relative error {} between {} and {}'.format(maxError[i], - headings0[i]['label'], - headings1[i]['label'])) + logging.info(f" * maximum relative error {maxError[i]} between {headings0[i]['label']} and {headings1[i]['label']}") logging.info(' ********') return maxError @@ -480,8 +472,8 @@ class Test: normedDelta = np.where(normBy>preFilter,delta/normBy,0.0) mean = np.amax(np.abs(np.mean(normedDelta,0))) std = np.amax(np.std(normedDelta,0)) - logging.info('mean: {:f}'.format(mean)) - logging.info('std: {:f}'.format(std)) + logging.info(f'mean: {mean:f}') + logging.info(f'std: {std:f}') return (mean atol + rtol*np.absolute(data[1]) - logging.info('shape of violators: {}'.format(violators.shape)) + logging.info(f'shape of violators: {violators.shape}') for j,culprits in enumerate(violators): goodguys = np.logical_not(culprits) if culprits.any(): - logging.info('{} has {}'.format(j,np.sum(culprits))) - logging.info('deviation: {}'.format(np.absolute(data[0][j]-data[1][j])[culprits])) - logging.info('data : {}'.format(np.absolute(data[1][j])[culprits])) - logging.info('deviation: {}'.format(np.absolute(data[0][j]-data[1][j])[goodguys])) - logging.info('data : {}'.format(np.absolute(data[1][j])[goodguys])) + logging.info(f'{j} has {np.sum(culprits)}') + logging.info(f'deviation: {np.absolute(data[0][j]-data[1][j])[culprits]}') + logging.info(f'data : {np.absolute(data[1][j])[culprits]}') + logging.info(f'deviation: {np.absolute(data[0][j]-data[1][j])[goodguys]}') + logging.info(f'data : {np.absolute(data[1][j])[goodguys]}') allclose = True # start optimistic for i in range(1,len(data)): @@ -588,12 +580,12 @@ class Test: if culprit == 0: count = len(self.variants) if self.options.select is None else len(self.options.select) - msg = 'Test passed.' if count == 1 else 'All {} tests passed.'.format(count) + msg = 'Test passed.' if count == 1 else f'All {count} tests passed.' elif culprit == -1: msg = 'Warning: could not start test...' ret = 0 else: - msg = 'Test "{}" failed.'.format(self.variantName(culprit-1)) + msg = f'Test "{self.variantName(culprit-1)}" failed.' logging.critical('\n'.join(['*'*40,msg,'*'*40]) + '\n') return ret diff --git a/python/damask/config/material.py b/python/damask/config/material.py index a30db1242..433c70d03 100644 --- a/python/damask/config/material.py +++ b/python/damask/config/material.py @@ -1,6 +1,7 @@ import re import os +from damask import util class Section(): def __init__(self,data = {'__order__':[]},part = ''): @@ -27,7 +28,7 @@ class Section(): if multiKey not in self.parameters: self.parameters[multiKey] = [] if multiKey not in self.parameters['__order__']: self.parameters['__order__'] += [multiKey] self.parameters[multiKey] += [[item] for item in data] if isinstance(data, list) else [[data]] - + def data(self): return self.parameters @@ -42,27 +43,27 @@ class Crystallite(Section): def __init__(self,data = {'__order__':[]}): """New material.config section.""" Section.__init__(self,data) - + class Phase(Section): def __init__(self,data = {'__order__':[]}): """New material.config section.""" Section.__init__(self,data) - + class Microstructure(Section): def __init__(self,data = {'__order__':[]}): """New material.config section.""" Section.__init__(self,data) - - + + class Texture(Section): def __init__(self,data = {'__order__':[]}): """New material.config section.""" Section.__init__(self,data) - + def add_component(self,theType,properties): - + scatter = properties['scatter'] if 'scatter' in list(map(str.lower,list(properties.keys()))) else 0.0 fraction = properties['fraction'] if 'fraction' in list(map(str.lower,list(properties.keys()))) else 1.0 @@ -81,45 +82,45 @@ class Texture(Section): ) ) - + class Material(): """Read, manipulate, and write material.config files.""" def __init__(self,verbose=True): """Generates ordered list of parts.""" self.parts = [ - 'homogenization', - 'crystallite', - 'phase', - 'texture', - 'microstructure', + 'homogenization', + 'crystallite', + 'phase', + 'texture', + 'microstructure', ] - self.data = {\ - 'homogenization': {'__order__': []}, - 'microstructure': {'__order__': []}, - 'crystallite': {'__order__': []}, - 'phase': {'__order__': []}, - 'texture': {'__order__': []}, - } + self.data = { + 'homogenization': {'__order__': []}, + 'microstructure': {'__order__': []}, + 'crystallite': {'__order__': []}, + 'phase': {'__order__': []}, + 'texture': {'__order__': []}, + } self.verbose = verbose - + def __repr__(self): """Returns current data structure in material.config format.""" me = [] for part in self.parts: - if self.verbose: print('processing <{}>'.format(part)) + if self.verbose: print(f'processing <{part}>') me += ['', '#'*100, - '<{}>'.format(part), + f'<{part}>', '#'*100, ] for section in self.data[part]['__order__']: - me += ['[{}] {}'.format(section,'#'+'-'*max(0,96-len(section)))] + me += [f'[{section}] {"#"+"-"*max(0,96-len(section))}'] for key in self.data[part][section]['__order__']: if key.startswith('(') and key.endswith(')'): # multiple (key) - me += ['{}\t{}'.format(key,' '.join(values)) for values in self.data[part][section][key]] + me += [f'{key}\t{" ".join(values)}' for values in self.data[part][section][key]] else: # plain key - me += ['{}\t{}'.format(key,' '.join(map(str,self.data[part][section][key])))] + me += [f'{key}\t{util.srepr(self.data[part][section][key]," ")}'] return '\n'.join(me) + '\n' def parse(self, part=None, sections=[], content=None): @@ -158,7 +159,7 @@ class Material(): self.data[part][name_section][items[0]].append(items[1:]) else: # plain key self.data[part][name_section][items[0]] = items[1:] - + def read(self,filename=None): @@ -174,20 +175,20 @@ class Material(): recursiveRead(match.group(1) if match.group(1).startswith('/') else os.path.normpath(os.path.join(os.path.dirname(filename),match.group(1)))) return result - + c = recursiveRead(filename) for p in self.parts: self.parse(part=p, content=c) - + def write(self,filename='material.config', overwrite=False): """Write to material.config.""" i = 0 outname = filename while os.path.exists(outname) and not overwrite: i += 1 - outname = '{}_{}'.format(filename,i) + outname = f'{filename}_{i}' - if self.verbose: print('Writing material data to {}'.format(outname)) + if self.verbose: print(f'Writing material data to {outname}') with open(outname,'w') as f: f.write(str(self)) return outname @@ -196,11 +197,11 @@ class Material(): """Add Update.""" part = part.lower() section = section.lower() - if part not in self.parts: raise Exception('invalid part {}'.format(part)) + if part not in self.parts: raise Exception(f'invalid part {part}') if not isinstance(initialData, dict): initialData = initialData.data() - + if section not in self.data[part]: self.data[part]['__order__'] += [section] if section in self.data[part] and merge: for existing in self.data[part][section]['__order__']: # replace existing @@ -215,9 +216,9 @@ class Material(): self.data[part][section]['__order__'] += [new] else: self.data[part][section] = initialData - - - + + + def add_microstructure(self, section='', components={}, # dict of phase,texture, and fraction lists @@ -226,7 +227,7 @@ class Material(): microstructure = Microstructure() # make keys lower case (http://stackoverflow.com/questions/764235/dictionary-to-lowercase-in-python) components=dict((k.lower(), v) for k,v in components.items()) - + for key in ['phase','texture','fraction','crystallite']: if isinstance(components[key], list): for i, x in enumerate(components[key]): @@ -239,7 +240,7 @@ class Material(): components[key] = [components[key].lower()] except AttributeError: components[key] = [components[key]] - + for (phase,texture,fraction,crystallite) in zip(components['phase'],components['texture'], components['fraction'],components['crystallite']): microstructure.add_multiKey('constituent','phase %i\ttexture %i\tfraction %g\ncrystallite %i'%( @@ -247,19 +248,19 @@ class Material(): self.data['texture']['__order__'].index(texture)+1, fraction, self.data['crystallite']['__order__'].index(crystallite)+1)) - + self.add_section('microstructure',section,microstructure) - - def change_value(self, part=None, - section=None, - key=None, + + def change_value(self, part=None, + section=None, + key=None, value=None): - if not isinstance(value,list): + if not isinstance(value,list): if not isinstance(value,str): value = '%s'%value value = [value] - newlen = len(value) + newlen = len(value) oldval = self.data[part.lower()][section.lower()][key.lower()] oldlen = len(oldval) print('changing %s:%s:%s from %s to %s '%(part.lower(),section.lower(),key.lower(),oldval,value)) diff --git a/python/damask/solver/_marc.py b/python/damask/solver/_marc.py index dfe45a46d..23ae5639d 100644 --- a/python/damask/solver/_marc.py +++ b/python/damask/solver/_marc.py @@ -25,7 +25,7 @@ class Marc: def library_path(self): path_MSC = Environment().options['MSC_ROOT'] - path_lib = Path('{}/mentat{}/shlib/linux64'.format(path_MSC,self.version)) + path_lib = Path(f'{path_MSC}/mentat{self.version}/shlib/linux64') return path_lib if path_lib.is_dir() else None @@ -34,7 +34,7 @@ class Marc: def tools_path(self): path_MSC = Environment().options['MSC_ROOT'] - path_tools = Path('{}/marc{}/tools'.format(path_MSC,self.version)) + path_tools = Path(f'{path_MSC}/marc{self.version}/tools') return path_tools if path_tools.is_dir() else None @@ -51,21 +51,21 @@ class Marc: env = Environment() - user = env.root_dir/Path('src/DAMASK_marc{}'.format(self.version)).with_suffix('.f90' if compile else '.marc') - if not user.is_file(): - raise FileNotFoundError("DAMASK4Marc ({}) '{}' not found".format(('source' if compile else 'binary'),user)) + usersub = env.root_dir/Path(f'src/DAMASK_marc{self.version}').with_suffix('.f90' if compile else '.marc') + if not usersub.is_file(): + raise FileNotFoundError("DAMASK4Marc ({}) '{}' not found".format(('source' if compile else 'binary'),usersub)) # Define options [see Marc Installation and Operation Guide, pp 23] - script = 'run_damask_{}mp'.format(optimization) + script = f'run_damask_{optimization}mp' cmd = str(self.tools_path/Path(script)) + \ ' -jid ' + model + '_' + job + \ ' -nprocd 1 -autorst 0 -ci n -cr n -dcoup 0 -b no -v no' - if compile: cmd += ' -u ' + str(user) + ' -save y' - else: cmd += ' -prog ' + str(user.with_suffix('')) + if compile: cmd += ' -u ' + str(usersub) + ' -save y' + else: cmd += ' -prog ' + str(usersub.with_suffix('')) - print('job submission {} compilation: {}'.format('with' if compile else 'without',user)) + print('job submission {} compilation: {}'.format(('with' if compile else 'without'),usersub)) if logfile: log = open(logfile, 'w') print(cmd) process = subprocess.Popen(shlex.split(cmd),stdout = log,stderr = subprocess.STDOUT) diff --git a/python/damask/util.py b/python/damask/util.py index f67df461a..471a5069b 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -236,7 +236,7 @@ class _ProgressBar: delta_time = datetime.datetime.now() - self.start_time remaining_time = (self.total - (iteration+1)) * delta_time / (iteration+1) remaining_time -= datetime.timedelta(microseconds=remaining_time.microseconds) # remove μs - sys.stderr.write('\r{self.prefix} {bar} {fraction:>4.0%} ETA {remaining_time}') + sys.stderr.write(f'\r{self.prefix} {bar} {fraction:>4.0%} ETA {remaining_time}') sys.stderr.flush() self.last_fraction = fraction diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index 7a3823ab3..28c9d14b0 100644 --- a/python/tests/test_Geom.py +++ b/python/tests/test_Geom.py @@ -6,6 +6,7 @@ import numpy as np from damask import Geom from damask import Rotation +from damask import util def geom_equal(a,b): @@ -85,8 +86,8 @@ class TestGeom: def test_mirror(self,default,update,reference_dir,directions,reflect): modified = copy.deepcopy(default) modified.mirror(directions,reflect) - tag = 'directions={}_reflect={}'.format('-'.join(directions),reflect) - reference = os.path.join(reference_dir,'mirror_{}.geom'.format(tag)) + tag = f'directions={"-".join(directions)}_reflect={reflect}' + reference = os.path.join(reference_dir,f'mirror_{tag}.geom') if update: modified.to_file(reference) assert geom_equal(modified,Geom.from_file(reference)) @@ -94,8 +95,8 @@ class TestGeom: def test_clean(self,default,update,reference_dir,stencil): modified = copy.deepcopy(default) modified.clean(stencil) - tag = 'stencil={}'.format(stencil) - reference = os.path.join(reference_dir,'clean_{}.geom'.format(tag)) + tag = f'stencil={stencil}' + reference = os.path.join(reference_dir,f'clean_{tag}.geom') if update: modified.to_file(reference) assert geom_equal(modified,Geom.from_file(reference)) @@ -111,8 +112,8 @@ class TestGeom: def test_scale(self,default,update,reference_dir,grid): modified = copy.deepcopy(default) modified.scale(grid) - tag = 'grid={}'.format('-'.join([str(x) for x in grid])) - reference = os.path.join(reference_dir,'scale_{}.geom'.format(tag)) + tag = f'grid={util.srepr(grid,"-")}' + reference = os.path.join(reference_dir,f'scale_{tag}.geom') if update: modified.to_file(reference) assert geom_equal(modified,Geom.from_file(reference)) @@ -150,8 +151,8 @@ class TestGeom: def test_rotate(self,default,update,reference_dir,Eulers): modified = copy.deepcopy(default) modified.rotate(Rotation.from_Eulers(Eulers,degrees=True)) - tag = 'Eulers={}-{}-{}'.format(*Eulers) - reference = os.path.join(reference_dir,'rotate_{}.geom'.format(tag)) + tag = f'Eulers={util.srepr(Eulers,"-")}' + reference = os.path.join(reference_dir,f'rotate_{tag}.geom') if update: modified.to_file(reference) assert geom_equal(modified,Geom.from_file(reference)) diff --git a/python/tests/test_Lattice.py b/python/tests/test_Lattice.py index fb8ce8ea5..01d1aac0c 100644 --- a/python/tests/test_Lattice.py +++ b/python/tests/test_Lattice.py @@ -38,4 +38,4 @@ class TestSymmetry: def test_invalid_argument(self,function): s = Symmetry() # noqa with pytest.raises(ValueError): - eval('s.{}(np.ones(4))'.format(function)) + eval(f's.{function}(np.ones(4))') diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index a8b8afdac..e811b2043 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -49,7 +49,7 @@ class TestOrientation: @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) @pytest.mark.parametrize('lattice',['fcc','bcc']) def test_relationship_reference(self,update,reference_dir,model,lattice): - reference = os.path.join(reference_dir,'{}_{}.txt'.format(lattice,model)) + reference = os.path.join(reference_dir,f'{lattice}_{model}.txt') ori = Orientation(Rotation(),lattice) eu = np.array([o.rotation.as_Eulers(degrees=True) for o in ori.relatedOrientations(model)]) if update: diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index 38c0d84cc..6940d573d 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -137,7 +137,7 @@ class TestResult: default.add_Cauchy('P','F') default.add_eigenvalue('sigma',eigenvalue) loc = {'sigma' :default.get_dataset_location('sigma'), - 'lambda':default.get_dataset_location('lambda_{}(sigma)'.format(eigenvalue))} + 'lambda':default.get_dataset_location(f'lambda_{eigenvalue}(sigma)')} in_memory = function(mechanics.eigenvalues(default.read_dataset(loc['sigma'],0)),axis=1,keepdims=True) in_file = default.read_dataset(loc['lambda'],0) assert np.allclose(in_memory,in_file) @@ -147,7 +147,7 @@ class TestResult: default.add_Cauchy('P','F') default.add_eigenvector('sigma',eigenvalue) loc = {'sigma' :default.get_dataset_location('sigma'), - 'v(sigma)':default.get_dataset_location('v_{}(sigma)'.format(eigenvalue))} + 'v(sigma)':default.get_dataset_location(f'v_{eigenvalue}(sigma)')} in_memory = mechanics.eigenvectors(default.read_dataset(loc['sigma'],0))[:,idx] in_file = default.read_dataset(loc['v(sigma)'],0) assert np.allclose(in_memory,in_file) @@ -179,7 +179,7 @@ class TestResult: t = ['V','U'][np.random.randint(0,2)] m = np.random.random()*2.0 - 1.0 default.add_strain_tensor('F',t,m) - label = 'epsilon_{}^{}(F)'.format(t,m) + label = f'epsilon_{t}^{m}(F)' default.add_Mises(label) loc = {label :default.get_dataset_location(label), label+'_vM':default.get_dataset_location(label+'_vM')} @@ -248,7 +248,7 @@ class TestResult: t = ['V','U'][np.random.randint(0,2)] m = np.random.random()*2.0 - 1.0 default.add_strain_tensor('F',t,m) - label = 'epsilon_{}^{}(F)'.format(t,m) + label = f'epsilon_{t}^{m}(F)' loc = {'F': default.get_dataset_location('F'), label: default.get_dataset_location(label)} in_memory = mechanics.strain_tensor(default.read_dataset(loc['F'],0),t,m) diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 86e5fa2a7..89173bb5a 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -556,7 +556,7 @@ def mul(me, other): else: raise ValueError('Can only rotate vectors, 2nd order tensors, and 4th order tensors') else: - raise TypeError('Cannot rotate {}'.format(type(other))) + raise TypeError(f'Cannot rotate {type(other)}') class TestRotation: @@ -878,7 +878,7 @@ class TestRotation: def test_invalid_P(self,fr,to): R = Rotation.from_random(np.random.randint(8,32,(3))) # noqa with pytest.raises(ValueError): - fr(eval('R.{}()'.format(to)),P=-30) + fr(eval(f'R.{to}()'),P=-30) @pytest.mark.parametrize('shape',[None,(3,),(4,2)]) def test_broadcast(self,shape): diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py index 8a343e26b..3acb554a7 100644 --- a/python/tests/test_grid_filters.py +++ b/python/tests/test_grid_filters.py @@ -30,8 +30,8 @@ class TestGridFilters: grid = np.random.randint(8,32,(3)) size = np.random.random(3) origin = np.random.random(3) - coord0 = eval('grid_filters.{}_coord0(grid,size,origin)'.format(mode)) # noqa - _grid,_size,_origin = eval('grid_filters.{}_coord0_gridSizeOrigin(coord0.reshape(-1,3,order="F"))'.format(mode)) + coord0 = eval(f'grid_filters.{mode}_coord0(grid,size,origin)') # noqa + _grid,_size,_origin = eval(f'grid_filters.{mode}_coord0_gridSizeOrigin(coord0.reshape(-1,3,order="F"))') assert np.allclose(grid,_grid) and np.allclose(size,_size) and np.allclose(origin,_origin) def test_displacement_fluct_equivalence(self): @@ -67,8 +67,8 @@ class TestGridFilters: origin= np.random.random(3) size = np.random.random(3) # noqa grid = np.random.randint(8,32,(3)) - shifted = eval('grid_filters.{}_coord0(grid,size,origin)'.format(mode)) - unshifted = eval('grid_filters.{}_coord0(grid,size)'.format(mode)) + shifted = eval(f'grid_filters.{mode}_coord0(grid,size,origin)') + unshifted = eval(f'grid_filters.{mode}_coord0(grid,size)') if mode == 'cell': assert np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid) +(3,))) elif mode == 'node': From 99995602479bdab432cb095513af1d515ecc303e Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 24 Jun 2020 18:36:43 -0400 Subject: [PATCH 40/46] easier understanding of from_ang data layout interpretation --- python/damask/_table.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/python/damask/_table.py b/python/damask/_table.py index 8c09e2b78..435b7afa9 100644 --- a/python/damask/_table.py +++ b/python/damask/_table.py @@ -25,7 +25,7 @@ class Table: """ self.comments = [] if comments is None else [c for c in comments] self.data = pd.DataFrame(data=data) - self.shapes = { k:(v,) if isinstance(v,(np.integer,int)) else v for k,v in shapes.items() } + self.shapes = { k:(v,) if isinstance(v,(np.int,int)) else v for k,v in shapes.items() } self._label_condensed() @@ -126,8 +126,6 @@ class Table: Filename or file for reading. """ - shapes = {'eu':(3,), 'pos':(2,), - 'IQ':(1,), 'CI':(1,), 'ID':(1,), 'intensity':(1,), 'fit':(1,)} try: f = open(fname) except TypeError: @@ -144,8 +142,11 @@ class Table: break data = np.loadtxt(content) - for c in range(data.shape[1]-10): - shapes[f'n/a_{c+1}'] = (1,) + + shapes = {'eu':3, 'pos':2, 'IQ':1, 'CI':1, 'ID':1, 'intensity':1, 'fit':1} + remainder = data.shape[1]-sum(shapes.values()) + if remainder > 0: # 3.8 can do: if (remainder := data.shape[1]-sum(shapes.values())) > 0 + shapes['unknown'] = remainder return Table(data,shapes,comments) @@ -234,7 +235,6 @@ class Table: """ self.data.drop(columns=label,inplace=True) - del self.shapes[label] From e64d353865337d48ba5123be583a71e08983f7af Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 24 Jun 2020 18:37:33 -0400 Subject: [PATCH 41/46] condensed scale_to_coprime; added test of scale_to_coprime --- python/damask/util.py | 11 +++++++---- python/tests/test_util.py | 23 ++++++++++++++++++++++- 2 files changed, 29 insertions(+), 5 deletions(-) diff --git a/python/damask/util.py b/python/damask/util.py index 471a5069b..cf99dc1b9 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -166,10 +166,13 @@ def scale_to_coprime(v): """Least common multiple.""" return a * b // np.gcd(a, b) - denominators = [int(get_square_denominator(i)) for i in v] - s = reduce(lcm, denominators) ** 0.5 - m = (np.array(v)*s).astype(np.int) - return m//reduce(np.gcd,m) + m = (np.array(v) * reduce(lcm, map(lambda x: int(get_square_denominator(x)),v)) ** 0.5).astype(np.int) + m = m//reduce(np.gcd,m) + + if not np.allclose(v/m,v[0]/m[0]): + raise ValueError(f'Invalid result {m} for input {v}. Insufficient precision?') + + return m #################################################################################################### diff --git a/python/tests/test_util.py b/python/tests/test_util.py index e67478d82..de6926e8c 100644 --- a/python/tests/test_util.py +++ b/python/tests/test_util.py @@ -1,11 +1,32 @@ +import pytest +import numpy as np from damask import util + class TestUtil: def test_execute_direct(self): out,err = util.execute('echo test') assert out=='test\n' and err=='' - + def test_execute_env(self): out,err = util.execute('sh -c "echo $test_for_execute"',env={'test_for_execute':'test'}) assert out=='test\n' and err=='' + + def test_croak(self): + util.croak('Burp!') + + @pytest.mark.parametrize('input,output', + [ + ([0.5,0.5],[1,1]), + ([1./2.,1./3.],[3,2]), + ([2./3.,1./2.,1./3.],[4,3,2]), + ]) + + def test_scale2coprime(self,input,output): + assert np.allclose(util.scale_to_coprime(np.array(input)), + np.array(output).astype(int)) + + def test_lackofprecision(self): + with pytest.raises(ValueError): + util.scale_to_coprime(np.array([0.66,0.5,0.33])) From c9829f0f1f0ba834d16906d06b229c848b0ca262 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 25 Jun 2020 08:18:39 +0200 Subject: [PATCH 42/46] only Chuck Norris can divide by zero --- python/damask/util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/util.py b/python/damask/util.py index 3c1f29c85..8668a9060 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -169,7 +169,7 @@ def scale_to_coprime(v): m = (np.array(v) * reduce(lcm, map(lambda x: int(get_square_denominator(x)),v)) ** 0.5).astype(np.int) m = m//reduce(np.gcd,m) - if not np.allclose(v/m,v[0]/m[0]): + if not np.allclose(v[v.nonzero()]/m[v.nonzero()],v[v.nonzero()][0]/m[m.nonzero()][0]): raise ValueError(f'Invalid result {m} for input {v}. Insufficient precision?') return m From 5d7213b062a6da6771e3ed1353398af7455cdf47 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 25 Jun 2020 08:19:07 +0200 Subject: [PATCH 43/46] restore reminder to document changes causes conflicts with outdated tests and is not urgent: postpone until release. --- python/damask/_geom.py | 11 ++++++++++- python/damask/_vtk.py | 2 +- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/python/damask/_geom.py b/python/damask/_geom.py index 5a3e6d47f..bbdcaf9f8 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -230,7 +230,7 @@ class Geom: return np.copy(self.microstructure) - def get_size(self,): + def get_size(self): """Return the physical size in meter.""" return np.copy(self.size) @@ -374,6 +374,7 @@ class Geom: else: microstructure = microstructure.reshape(grid) + #ToDo: comments = 'geom.py:from_Laguerre_tessellation v{}'.format(version) return Geom(microstructure+1,size,homogenization=1) @@ -398,6 +399,7 @@ class Geom: KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds) devNull,microstructure = KDTree.query(coords) + #ToDo: comments = 'geom.py:from_Voronoi_tessellation v{}'.format(version) return Geom(microstructure.reshape(grid)+1,size,homogenization=1) @@ -522,6 +524,7 @@ class Geom: if 'x' in directions: ms = np.concatenate([ms,ms[limits[0]:limits[1]:-1,:,:]],0) + #ToDo: self.add_comments('geom.py:mirror v{}'.format(version) return self.update(ms,rescale=True) @@ -535,6 +538,7 @@ class Geom: number of grid points in x,y,z direction. """ + #ToDo: self.add_comments('geom.py:scale v{}'.format(version) return self.update( ndimage.interpolation.zoom( self.microstructure, @@ -561,6 +565,7 @@ class Geom: unique, inverse = np.unique(arr, return_inverse=True) return unique[np.argmax(np.bincount(inverse))] + #ToDo: self.add_comments('geom.py:clean v{}'.format(version) return self.update(ndimage.filters.generic_filter( self.microstructure, mostFrequent, @@ -575,6 +580,7 @@ class Geom: for i, oldID in enumerate(np.unique(self.microstructure)): renumbered = np.where(self.microstructure == oldID, i+1, renumbered) + #ToDo: self.add_comments('geom.py:renumber v{}'.format(version) return self.update(renumbered) @@ -609,6 +615,7 @@ class Geom: origin = self.origin-(np.asarray(microstructure_in.shape)-self.grid)*.5 * self.size/self.grid + #ToDo: self.add_comments('geom.py:rotate v{}'.format(version) return self.update(microstructure_in,origin=origin,rescale=True) @@ -640,6 +647,7 @@ class Geom: canvas[l[0]:r[0],l[1]:r[1],l[2]:r[2]] = self.microstructure[L[0]:R[0],L[1]:R[1],L[2]:R[2]] + #ToDo: self.add_comments('geom.py:canvas v{}'.format(version) return self.update(canvas,origin=self.origin+offset*self.size/self.grid,rescale=True) @@ -659,4 +667,5 @@ class Geom: for from_ms,to_ms in zip(from_microstructure,to_microstructure): substituted[self.microstructure==from_ms] = to_ms + #ToDo: self.add_comments('geom.py:substitute v{}'.format(version) return self.update(substituted) diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py index 55b84445c..f4855820e 100644 --- a/python/damask/_vtk.py +++ b/python/damask/_vtk.py @@ -178,7 +178,7 @@ class VTK: default_ext = writer.GetDefaultFileExtension() ext = Path(fname).suffix if ext and ext != '.'+default_ext: - raise ValueError(f'Given extension {ext} does not match default {default_ext}') + raise ValueError(f'Given extension {ext} does not match default .{default_ext}') writer.SetFileName(str(Path(fname).with_suffix('.'+default_ext))) writer.SetCompressorTypeToZLib() writer.SetDataModeToBinary() From 3290e2c585aab9edad7e4de0b852ef334ca1f83c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 25 Jun 2020 08:29:36 +0200 Subject: [PATCH 44/46] handle even obscure directions --- python/damask/util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/util.py b/python/damask/util.py index 8668a9060..733fd0010 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -156,7 +156,7 @@ def show_progress(iterable,N_iter=None,prefix='',bar_length=50): def scale_to_coprime(v): """Scale vector to co-prime (relatively prime) integers.""" - MAX_DENOMINATOR = 1000 + MAX_DENOMINATOR = 1000000 def get_square_denominator(x): """Denominator of the square of a number.""" From 27220a03bcf63b6268c91b7eee7890545b92af95 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 25 Jun 2020 09:23:43 +0200 Subject: [PATCH 45/46] not invalid anymore --- python/tests/test_util.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/tests/test_util.py b/python/tests/test_util.py index de6926e8c..c786e0de1 100644 --- a/python/tests/test_util.py +++ b/python/tests/test_util.py @@ -18,6 +18,7 @@ class TestUtil: @pytest.mark.parametrize('input,output', [ + ([2,0],[1,0]), ([0.5,0.5],[1,1]), ([1./2.,1./3.],[3,2]), ([2./3.,1./2.,1./3.],[4,3,2]), @@ -29,4 +30,4 @@ class TestUtil: def test_lackofprecision(self): with pytest.raises(ValueError): - util.scale_to_coprime(np.array([0.66,0.5,0.33])) + util.scale_to_coprime(np.array([1/3333,1,1])) From 9c89d537d25a72c1afe48018a25099a21096dd1d Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 25 Jun 2020 20:17:00 +0200 Subject: [PATCH 46/46] [skip ci] updated version information after successful test of v2.0.3-2706-g8ce5da55 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index d94a70ab1..4246a9b8c 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2692-g3d93a5ff +v2.0.3-2706-g8ce5da55