Merge branch 'development' into HMS

.gitignore

@@ -1,4 +1,3 @@
-.noH5py
 *.pyc
 *.mod
 *.o
@@ -8,3 +7,4 @@
 bin
 PRIVATE
 build
+system_report.txt

CMakeLists.txt

@@ -124,6 +124,7 @@ endif ()
 
 # Predefined sets for OPTIMIZATION/OPENMP based on BUILD_TYPE
 if ("${CMAKE_BUILD_TYPE}" STREQUAL "DEBUG" OR "${CMAKE_BUILD_TYPE}" STREQUAL "SYNTAXONLY" )
+    set (DEBUG_FLAGS "${DEBUG_FLAGS} -DDEBUG")
     set (PARALLEL "OFF")
     set (OPTI "OFF")
 elseif ("${CMAKE_BUILD_TYPE}" STREQUAL "RELEASE")

DAMASK_prerequisites.sh

@@ -1,9 +1,17 @@
 #!/usr/bin/env bash
 
-OUTFILE="system_report.txt"
-echo generating $OUTFILE
+#==================================================================================================
+# Execute this script (type './DAMASK_prerequisites.sh') 
+# and send system_report.txt to damask@mpie.de for support
+#==================================================================================================
+
+OUTFILE="system_report.txt"
+echo ===========================================
+echo +  Generating $OUTFILE                    
+echo +  Send to damask@mpie.de for support
+echo ===========================================
+
 
-echo date +"%m-%d-%y" >$OUTFILE
 
 # redirect STDOUT and STDERR to logfile
 # https://stackoverflow.com/questions/11229385/redirect-all-output-in-a-bash-script-when-using-set-x^
@@ -13,6 +21,10 @@ exec > $OUTFILE 2>&1
 # https://stackoverflow.com/questions/59895/getting-the-source-directory-of-a-bash-script-from-within
 DAMASK_ROOT="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )"
 
+echo XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
+echo System report for \'$(hostname)\' created on $(date '+%Y-%m-%d %H:%M:%S') by \'$(whoami)\'
+echo XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
+echo
 echo ==============================================================================================
 echo DAMASK settings
 echo ==============================================================================================
@@ -30,19 +42,24 @@ echo System
 echo ==============================================================================================
 uname -a
 echo
+echo PATH: $PATH
+echo LD_LIBRARY_PATH: $LD_LIBRARY_PATH
+echo PYTHONPATH: $PYTHONPATH
+echo SHELL: $SHELL
+echo 
 echo ==============================================================================================
 echo Python
 echo ==============================================================================================
 
 DEFAULT_PYTHON=python2.7
 for executable in python python2 python3 python2.7; do
-  if [[ "$(which $executable)x" != "x" ]]; then
+  if which $executable &> /dev/null; then
     echo $executable version: $($executable --version 2>&1)
   else
     echo $executable does not exist
   fi
 done
-echo Location of $DEFAULT_PYTHON: $(ls -la $(which $DEFAULT_PYTHON))
+echo Details on $DEFAULT_PYTHON: $(ls -la $(which $DEFAULT_PYTHON))
 echo
 for module in numpy scipy;do
   echo ----------------------------------------------------------------------------------------------
@@ -69,7 +86,7 @@ echo ===========================================================================
 echo GCC
 echo ==============================================================================================
 for executable in gcc g++ gfortran ;do
-  if [[ "$(which $executable)x" != "x" ]]; then
+  if which $executable &> /dev/null; then
     echo $(which $executable) version: $($executable --version 2>&1)
   else
      echo $executable does not exist
@@ -80,7 +97,7 @@ echo ===========================================================================
 echo Intel Compiler Suite
 echo ==============================================================================================
 for executable in icc icpc ifort ;do
-  if [[ "$(which $executable)x" != "x" ]]; then
+  if which $executable &> /dev/null; then
     echo $(which $executable) version: $($executable --version 2>&1)
   else
     echo $executable does not exist
@@ -91,7 +108,7 @@ echo ===========================================================================
 echo MPI Wrappers
 echo ==============================================================================================
 for executable in mpicc mpiCC mpicxx mpicxx mpifort mpif90 mpif77; do
-  if [[ "$(which $executable)x" != "x" ]]; then
+  if which $executable &> /dev/null; then
     echo $(which $executable) version: $($executable --show 2>&1)
   else
      echo $executable does not exist

LICENSE

@@ -1,6 +1,6 @@
 Copyright 2011-17 Max-Planck-Institut für Eisenforschung GmbH
 
-This program is free software: you can redistribute it and/or modify
+DAMASK is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.

PRIVATE

@@ -1 +1 @@
-Subproject commit 55a263fc30c40c16ef337be050f8901dd2747390
+Subproject commit 5fc3188c86ea1f4159db87529ac3e3169fb56e5d

VERSION

@@ -1 +1 @@
-v2.0.1-960-geddc2a6
+v2.0.1-1035-gd80a255

examples/AbaqusStandard/SX_PX_compression.jnl

@@ -1 +1,23 @@
-# Save by abaqususer on Thu May 12 10:22:10 2011
+# Save by m.diehl on 2017_12_06-18.38.26; build 2017 2016_09_27-23.54.59 126836
+from abaqus import *
+upgradeMdb(
+    '/nethome/storage/raid4/m.diehl/DAMASK/examples/AbaqusStandard/SX_PX_compression-6.9-1.cae'
+    , 
+    '/nethome/storage/raid4/m.diehl/DAMASK/examples/AbaqusStandard/SX_PX_compression.cae')
+# Save by m.diehl on 2017_12_06-18.38.26; build 2017 2016_09_27-23.54.59 126836
+from part import *
+from material import *
+from section import *
+from assembly import *
+from step import *
+from interaction import *
+from load import *
+from mesh import *
+from optimization import *
+from job import *
+from sketch import *
+from visualization import *
+from connectorBehavior import *
+mdb.jobs['Job_sx-px'].setValues(description='compression', userSubroutine=
+    '$HOME/DAMASK/src/DAMASK_abaqus_std.f')
+# Save by m.diehl on 2017_12_06-18.39.44; build 2017 2016_09_27-23.54.59 126836

examples/ConfigFiles/numerics.config

@@ -49,7 +49,7 @@ maxVolDiscrepancy_RGC     1.0e-5       # maximum allowable relative volume discr
 volDiscrepancyMod_RGC     1.0e+12
 discrepancyPower_RGC      5.0
 
-fixed_seed                0            # put any number larger than zero, integer, if you want to have a pseudo random distribution
+random_seed               0            # any integer larger than zero seeds the random generator, otherwise random seeding
 
 ## spectral parameters ##
 err_div_tolAbs            1.0e-3       # absolute tolerance for fulfillment of stress equilibrium

img/COPYING

@@ -0,0 +1,94 @@
+Creative Commons Attribution-NoDerivatives 4.0 International Public License
+
+By exercising the Licensed Rights (defined below), You accept and agree to be bound by the terms and conditions of this Creative Commons Attribution-NoDerivatives 4.0 International Public License ("Public License"). To the extent this Public License may be interpreted as a contract, You are granted the Licensed Rights in consideration of Your acceptance of these terms and conditions, and the Licensor grants You such rights in consideration of benefits the Licensor receives from making the Licensed Material available under these terms and conditions.
+
+Section 1 – Definitions.
+
+    Adapted Material means material subject to Copyright and Similar Rights that is derived from or based upon the Licensed Material and in which the Licensed Material is translated, altered, arranged, transformed, or otherwise modified in a manner requiring permission under the Copyright and Similar Rights held by the Licensor. For purposes of this Public License, where the Licensed Material is a musical work, performance, or sound recording, Adapted Material is always produced where the Licensed Material is synched in timed relation with a moving image.
+    Copyright and Similar Rights means copyright and/or similar rights closely related to copyright including, without limitation, performance, broadcast, sound recording, and Sui Generis Database Rights, without regard to how the rights are labeled or categorized. For purposes of this Public License, the rights specified in Section 2(b)(1)-(2) are not Copyright and Similar Rights.
+    Effective Technological Measures means those measures that, in the absence of proper authority, may not be circumvented under laws fulfilling obligations under Article 11 of the WIPO Copyright Treaty adopted on December 20, 1996, and/or similar international agreements.
+    Exceptions and Limitations means fair use, fair dealing, and/or any other exception or limitation to Copyright and Similar Rights that applies to Your use of the Licensed Material.
+    Licensed Material means the artistic or literary work, database, or other material to which the Licensor applied this Public License.
+    Licensed Rights means the rights granted to You subject to the terms and conditions of this Public License, which are limited to all Copyright and Similar Rights that apply to Your use of the Licensed Material and that the Licensor has authority to license.
+    Licensor means the individual(s) or entity(ies) granting rights under this Public License.
+    Share means to provide material to the public by any means or process that requires permission under the Licensed Rights, such as reproduction, public display, public performance, distribution, dissemination, communication, or importation, and to make material available to the public including in ways that members of the public may access the material from a place and at a time individually chosen by them.
+    Sui Generis Database Rights means rights other than copyright resulting from Directive 96/9/EC of the European Parliament and of the Council of 11 March 1996 on the legal protection of databases, as amended and/or succeeded, as well as other essentially equivalent rights anywhere in the world.
+    You means the individual or entity exercising the Licensed Rights under this Public License. Your has a corresponding meaning.
+
+Section 2 – Scope.
+
+    License grant.
+        Subject to the terms and conditions of this Public License, the Licensor hereby grants You a worldwide, royalty-free, non-sublicensable, non-exclusive, irrevocable license to exercise the Licensed Rights in the Licensed Material to:
+            reproduce and Share the Licensed Material, in whole or in part; and
+            produce and reproduce, but not Share, Adapted Material.
+        Exceptions and Limitations. For the avoidance of doubt, where Exceptions and Limitations apply to Your use, this Public License does not apply, and You do not need to comply with its terms and conditions.
+        Term. The term of this Public License is specified in Section 6(a).
+        Media and formats; technical modifications allowed. The Licensor authorizes You to exercise the Licensed Rights in all media and formats whether now known or hereafter created, and to make technical modifications necessary to do so. The Licensor waives and/or agrees not to assert any right or authority to forbid You from making technical modifications necessary to exercise the Licensed Rights, including technical modifications necessary to circumvent Effective Technological Measures. For purposes of this Public License, simply making modifications authorized by this Section 2(a)(4) never produces Adapted Material.
+        Downstream recipients.
+            Offer from the Licensor – Licensed Material. Every recipient of the Licensed Material automatically receives an offer from the Licensor to exercise the Licensed Rights under the terms and conditions of this Public License.
+            No downstream restrictions. You may not offer or impose any additional or different terms or conditions on, or apply any Effective Technological Measures to, the Licensed Material if doing so restricts exercise of the Licensed Rights by any recipient of the Licensed Material.
+        No endorsement. Nothing in this Public License constitutes or may be construed as permission to assert or imply that You are, or that Your use of the Licensed Material is, connected with, or sponsored, endorsed, or granted official status by, the Licensor or others designated to receive attribution as provided in Section 3(a)(1)(A)(i).
+
+    Other rights.
+        Moral rights, such as the right of integrity, are not licensed under this Public License, nor are publicity, privacy, and/or other similar personality rights; however, to the extent possible, the Licensor waives and/or agrees not to assert any such rights held by the Licensor to the limited extent necessary to allow You to exercise the Licensed Rights, but not otherwise.
+        Patent and trademark rights are not licensed under this Public License.
+        To the extent possible, the Licensor waives any right to collect royalties from You for the exercise of the Licensed Rights, whether directly or through a collecting society under any voluntary or waivable statutory or compulsory licensing scheme. In all other cases the Licensor expressly reserves any right to collect such royalties.
+
+Section 3 – License Conditions.
+
+Your exercise of the Licensed Rights is expressly made subject to the following conditions.
+
+    Attribution.
+
+        If You Share the Licensed Material, You must:
+            retain the following if it is supplied by the Licensor with the Licensed Material:
+                identification of the creator(s) of the Licensed Material and any others designated to receive attribution, in any reasonable manner requested by the Licensor (including by pseudonym if designated);
+                a copyright notice;
+                a notice that refers to this Public License;
+                a notice that refers to the disclaimer of warranties;
+                a URI or hyperlink to the Licensed Material to the extent reasonably practicable;
+            indicate if You modified the Licensed Material and retain an indication of any previous modifications; and
+            indicate the Licensed Material is licensed under this Public License, and include the text of, or the URI or hyperlink to, this Public License.
+        For the avoidance of doubt, You do not have permission under this Public License to Share Adapted Material.
+        You may satisfy the conditions in Section 3(a)(1) in any reasonable manner based on the medium, means, and context in which You Share the Licensed Material. For example, it may be reasonable to satisfy the conditions by providing a URI or hyperlink to a resource that includes the required information.
+        If requested by the Licensor, You must remove any of the information required by Section 3(a)(1)(A) to the extent reasonably practicable.
+
+Section 4 – Sui Generis Database Rights.
+
+Where the Licensed Rights include Sui Generis Database Rights that apply to Your use of the Licensed Material:
+
+    for the avoidance of doubt, Section 2(a)(1) grants You the right to extract, reuse, reproduce, and Share all or a substantial portion of the contents of the database, provided You do not Share Adapted Material;
+    if You include all or a substantial portion of the database contents in a database in which You have Sui Generis Database Rights, then the database in which You have Sui Generis Database Rights (but not its individual contents) is Adapted Material; and
+    You must comply with the conditions in Section 3(a) if You Share all or a substantial portion of the contents of the database.
+
+For the avoidance of doubt, this Section 4 supplements and does not replace Your obligations under this Public License where the Licensed Rights include other Copyright and Similar Rights.
+
+Section 5 – Disclaimer of Warranties and Limitation of Liability.
+
+    Unless otherwise separately undertaken by the Licensor, to the extent possible, the Licensor offers the Licensed Material as-is and as-available, and makes no representations or warranties of any kind concerning the Licensed Material, whether express, implied, statutory, or other. This includes, without limitation, warranties of title, merchantability, fitness for a particular purpose, non-infringement, absence of latent or other defects, accuracy, or the presence or absence of errors, whether or not known or discoverable. Where disclaimers of warranties are not allowed in full or in part, this disclaimer may not apply to You.
+    To the extent possible, in no event will the Licensor be liable to You on any legal theory (including, without limitation, negligence) or otherwise for any direct, special, indirect, incidental, consequential, punitive, exemplary, or other losses, costs, expenses, or damages arising out of this Public License or use of the Licensed Material, even if the Licensor has been advised of the possibility of such losses, costs, expenses, or damages. Where a limitation of liability is not allowed in full or in part, this limitation may not apply to You.
+
+    The disclaimer of warranties and limitation of liability provided above shall be interpreted in a manner that, to the extent possible, most closely approximates an absolute disclaimer and waiver of all liability.
+
+Section 6 – Term and Termination.
+
+    This Public License applies for the term of the Copyright and Similar Rights licensed here. However, if You fail to comply with this Public License, then Your rights under this Public License terminate automatically.
+
+    Where Your right to use the Licensed Material has terminated under Section 6(a), it reinstates:
+        automatically as of the date the violation is cured, provided it is cured within 30 days of Your discovery of the violation; or
+        upon express reinstatement by the Licensor.
+    For the avoidance of doubt, this Section 6(b) does not affect any right the Licensor may have to seek remedies for Your violations of this Public License.
+    For the avoidance of doubt, the Licensor may also offer the Licensed Material under separate terms or conditions or stop distributing the Licensed Material at any time; however, doing so will not terminate this Public License.
+    Sections 1, 5, 6, 7, and 8 survive termination of this Public License.
+
+Section 7 – Other Terms and Conditions.
+
+    The Licensor shall not be bound by any additional or different terms or conditions communicated by You unless expressly agreed.
+    Any arrangements, understandings, or agreements regarding the Licensed Material not stated herein are separate from and independent of the terms and conditions of this Public License.
+
+Section 8 – Interpretation.
+
+    For the avoidance of doubt, this Public License does not, and shall not be interpreted to, reduce, limit, restrict, or impose conditions on any use of the Licensed Material that could lawfully be made without permission under this Public License.
+    To the extent possible, if any provision of this Public License is deemed unenforceable, it shall be automatically reformed to the minimum extent necessary to make it enforceable. If the provision cannot be reformed, it shall be severed from this Public License without affecting the enforceability of the remaining terms and conditions.
+    No term or condition of this Public License will be waived and no failure to comply consented to unless expressly agreed to by the Licensor.
+    Nothing in this Public License constitutes or may be interpreted as a limitation upon, or waiver of, any privileges and immunities that apply to the Licensor or You, including from the legal processes of any jurisdiction or authority.

img/DAMASK_QR-CodeBW.svg

@@ -5,11 +5,37 @@
    xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
    xmlns:svg="http://www.w3.org/2000/svg"
    xmlns="http://www.w3.org/2000/svg"
+   xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
+   xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
    viewBox="0 0 104.4 104.4"
    height="104.4"
    width="104.4"
    version="1.1"
-   id="svg4314">
+   id="svg4314"
+   sodipodi:docname="DAMASK_QR-CodeBW.svg"
+   inkscape:version="0.92.1 r15371">
+  <title
+     id="title4483">DAMASK QR code</title>
+  <sodipodi:namedview
+     pagecolor="#ffffff"
+     bordercolor="#666666"
+     borderopacity="1"
+     objecttolerance="10"
+     gridtolerance="10"
+     guidetolerance="10"
+     inkscape:pageopacity="0"
+     inkscape:pageshadow="2"
+     inkscape:window-width="1920"
+     inkscape:window-height="1001"
+     id="namedview6"
+     showgrid="false"
+     inkscape:zoom="2.2605364"
+     inkscape:cx="52.200001"
+     inkscape:cy="52.200001"
+     inkscape:window-x="-9"
+     inkscape:window-y="-9"
+     inkscape:window-maximized="1"
+     inkscape:current-layer="svg4314" />
   <metadata
      id="metadata4320">
     <rdf:RDF>
@@ -18,8 +44,32 @@
         <dc:format>image/svg+xml</dc:format>
         <dc:type
            rdf:resource="http://purl.org/dc/dcmitype/StillImage" />
-        <dc:title></dc:title>
+        <dc:title>DAMASK QR code</dc:title>
+        <cc:license
+           rdf:resource="http://creativecommons.org/licenses/by-nd/4.0/" />
+        <dc:creator>
+          <cc:Agent>
+            <dc:title>Franz Roters</dc:title>
+          </cc:Agent>
+        </dc:creator>
+        <dc:subject>
+          <rdf:Bag>
+            <rdf:li>DAMASK; Crystal Plasticity; Multi-Physics</rdf:li>
+          </rdf:Bag>
+        </dc:subject>
+        <dc:description>link to https://damask.mpie.de</dc:description>
       </cc:Work>
+      <cc:License
+         rdf:about="http://creativecommons.org/licenses/by-nd/4.0/">
+        <cc:permits
+           rdf:resource="http://creativecommons.org/ns#Reproduction" />
+        <cc:permits
+           rdf:resource="http://creativecommons.org/ns#Distribution" />
+        <cc:requires
+           rdf:resource="http://creativecommons.org/ns#Notice" />
+        <cc:requires
+           rdf:resource="http://creativecommons.org/ns#Attribution" />
+      </cc:License>
     </rdf:RDF>
   </metadata>
   <defs

img/LICENSE

@@ -0,0 +1,32 @@
+DAMASK logo:    Copyright 2011 Philip Eisenlohr, 
+                Max-Planck-Institut für Eisenforschung GmbH
+DAMASK QR code: Copyright 2016 Franz Roters,
+                Max-Planck-Institut für Eisenforschung GmbH
+
+These images are licensed according to the Attribution-NoDerivatives 4.0 International 
+(CC BY-ND 4.0) license.
+
+You should have received a copy of the Attribution-NoDerivatives 4.0 International license 
+along with these images. If not, see <http://creativecommons.org/licenses/by-nd/4.0/>
+
+The following is a human-readable summary of (and not a substitute for) the license:
+
+You are free to:
+
+    Share — copy and redistribute the material in any medium or format for any purpose,
+    even commercially.
+
+    The licensor cannot revoke these freedoms as long as you follow the license terms.
+
+Under the following terms:
+
+    Attribution — You must give appropriate credit, provide a link to the license, and 
+    indicate if changes were made. You may do so in any reasonable manner, but not in 
+    any way that suggests the licensor endorses you or your use.
+
+    NoDerivatives — If you remix, transform, or build upon the material, you may not
+    distribute the modified material.
+
+    No additional restrictions — You may not apply legal terms or technological measures
+    that legally restrict others from doing anything the license permits.
+

installation/mods_Abaqus/abaqus_v6.env

@@ -44,16 +44,6 @@ double_precision=BOTH
 # The user will not be asked whether old job files of the same name should be deleted.
 ask_delete=OFF
 
-# You can compile DAMASK into a library to be used with abaqus
-# it saves you from compiling the subroutine for each job
-# in this case you do not have to specify a usersubroutine file
-# however if you still do, the compiled version will override that in the library
-# Procedure:
-# 1. create a library directory, e.g. abqlib, in your prefered location
-# 2. build the library replacing your_prefered_location/abqlib with the correct path to the directory created in 1.:
-#    abaqus make -l DAMASK_abaqus_std.f -dir your_prefered_location/abqlib
-#    abaqus make -l DAMASK_abaqus_exp.f -dir your_prefered_location/abqlib
-# 3. uncomment the next line after replacing your_prefered_location/abqlib with the correct path to the directory created in 1.
 # usub_lib_dir='your_prefered_location/abqlib'
 
 # Remove the temporary names from the namespace

installation/mods_Abaqus/abaqus_v6_serial.env

@@ -44,16 +44,6 @@ double_precision=BOTH
 # The user will not be asked whether old job files of the same name should be deleted.
 ask_delete=OFF
 
-# You can compile DAMASK into a library to be used with abaqus
-# it saves you from compiling the subroutine for each job
-# in this case you do not have to specify a usersubroutine file
-# however if you still do, the compiled version will override that in the library
-# Procedure:
-# 1. create a library directory, e.g. abqlib, in your prefered location
-# 2. build the library replacing your_prefered_location/abqlib with the correct path to the directory created in 1.:
-#    abaqus make -l DAMASK_abaqus_std.f -dir your_prefered_location/abqlib
-#    abaqus make -l DAMASK_abaqus_exp.f -dir your_prefered_location/abqlib
-# 3. uncomment the next line after replacing your_prefered_location/abqlib with the correct path to the directory created in 1.
 # usub_lib_dir='your_prefered_location/abqlib'
 
 # Remove the temporary names from the namespace

installation/patch/README.md

@@ -1,6 +1,6 @@
 # DAMASK patching
 
-This folder contains patches that modify the functionality of the current version of DAMASK prior to the corresponding inclusion in the official release.
+This folder contains patches that modify the functionality of the current development version of DAMASK ahead of the corresponding adoption in the official release.
 
 ## Usage
 
@@ -13,3 +13,13 @@ patch -p1 < installation/patch/nameOfPatch
 
   * **fwbw_derivative** switches the default spatial derivative from continuous to forward/backward difference.  
     This generally reduces spurious oscillations in the result as the spatial accuracy of the derivative is then compatible with the underlying solution grid.
+
+  * **PETSc-3.8** adjusts all includes and calls to PETSc to follow the 3.8.x API.
+    This allows to use the most recent version of PETSc.
+
+## Create patch
+commit your changes
+
+```bash
+git format-patch PATH_TO_COMPARE --stdout >
+```

lib/damask/DS_HDF5.xml

@@ -1,198 +0,0 @@
-<?xml version="1.0"?>
-<h5ds>
-    <!--Top level attributes-->
-    <history>
-        <type>attr</type>
-        <h5path>/</h5path>
-        <category></category>
-        <description>store cmd history</description>
-    </history>
-
-    <inc>
-        <type>attr</type>
-        <h5path>/</h5path>
-        <category></category>
-        <description></description>
-    </inc>
-
-    <!--Geometry section-->
-    <Vx>
-        <type>Scalar</type>
-        <h5path>/Geometry/Vx</h5path>
-        <category>Geometry</category>
-        <description>Vector along x define the spectral mesh</description>
-    </Vx>
-
-    <Vy>
-        <type>Scalar</type>
-        <h5path>/Geometry/Vy</h5path>
-        <category>Geometry</category>
-        <description>Vector along y defines the spectral mesh</description>
-    </Vy>
-
-    <Vz>
-        <type>Scalar</type>
-        <h5path>/Geometry/Vz</h5path>
-        <category>Geometry</category>
-        <description>Vector along z defines the spectral mesh</description>
-    </Vz>
-
-    <ip>
-        <type>Scalar</type>
-        <h5path>/Geometry/ip</h5path>
-        <category>Geometry</category>
-        <description></description>
-    </ip>
-
-    <node>
-        <type>Scalar</type>
-        <h5path>/Geometry/node</h5path>
-        <category>Geometry</category>
-        <description></description>
-    </node>
-
-    <grain>
-        <type>Scalar</type>
-        <h5path>/Geometry/grain</h5path>
-        <category>Geometry</category>
-        <description></description>
-    </grain>
-
-    <pos>
-        <type>Vector</type>
-        <h5path>/Geometry/pos</h5path>
-        <category>Geometry</category>
-        <description></description>
-    </pos>
-
-    <elem>
-        <type>Scalar</type>
-        <h5path>/Geometry/elem</h5path>
-        <category>Geometry</category>
-        <description></description>
-    </elem>
-
-    <!--Crystallite section-->
-    <phase>
-        <type>Scalar</type>
-        <h5path>/Crystallite/phase</h5path>
-        <category>Crystallite</category>
-        <description></description>
-    </phase>
-
-    <texture>
-        <type>Scalar</type>
-        <h5path>/Crystallite/texture</h5path>
-        <category>Crystallite</category>
-        <description></description>
-    </texture>
-
-    <volume>
-        <type>Scalar</type>
-        <h5path>/Crystallite/volume</h5path>
-        <category>Crystallite</category>
-        <description></description>
-    </volume>
-
-    <orientation>
-        <type>Vector</type>
-        <h5path>/Crystallite/orientation</h5path>
-        <category>Crystallite</category>
-        <description></description>
-    </orientation>
-
-    <eulerangles>
-        <type>Vector</type>
-        <h5path>/Crystallite/eulerangles</h5path>
-        <category>Crystallite</category>
-        <description>Bunnge Euler angles in degrees</description>
-    </eulerangles>
-
-    <grainrotation>
-        <type>Vector</type>
-        <h5path>/Crystallite/grainrotation</h5path>
-        <category>Crystallite</category>
-        <description></description>
-    </grainrotation>
-
-    <f>
-        <type>Tensor</type>
-        <h5path>/Crystallite/f</h5path>
-        <category>Crystallite</category>
-        <description>deformation gradient (F)</description>
-    </f>
-
-    <p>
-        <type>Tensor</type>
-        <h5path>/Crystallite/p</h5path>
-        <category>Crystallite</category>
-        <description>Pikola Kirkhoff stress</description>
-    </p>
-
-    <Cauchy>
-        <type>Tensor</type>
-        <h5path>/Crystallite/Cauchy</h5path>
-        <category>Crystallite</category>
-        <description>Cauchy stress tensor</description>
-    </Cauchy>
-
-    <lnV>
-        <type>Tensor</type>
-        <h5path>/Crystallite/lnV</h5path>
-        <category>Crystallite</category>
-        <description></description>
-    </lnV>
-
-    <MisesCauchy>
-        <type>Scalar</type>
-        <h5path>/Crystallite/MisesCauchy</h5path>
-        <category>Crystallite</category>
-        <description>von Mises equivalent of Cauchy stress</description>
-    </MisesCauchy>
-
-    <MiseslnV>
-        <type>Scalar</type>
-        <h5path>/Crystallite/MiseslnV</h5path>
-        <category>Crystallite</category>
-        <description>left von Mises strain</description>
-    </MiseslnV>
-
-    <!--Constitutive section-->
-    <resistance_slip>
-        <type>Vector</type>
-        <h5path>/Constitutive/resistance_slip</h5path>
-        <category>Constitutive</category>
-        <description></description>
-    </resistance_slip>
-
-    <shearrate_slip>
-        <type>Vector</type>
-        <h5path>/Constitutive/shearrate_slip</h5path>
-        <category>Constitutive</category>
-        <description></description>
-    </shearrate_slip>
-
-    <resolvedstress_slip>
-        <type>Vector</type>
-        <h5path>/Constitutive/resolvedstress_slip</h5path>
-        <category>Constitutive</category>
-        <description></description>
-    </resolvedstress_slip>
-
-    <totalshear>
-        <type>Scalar</type>
-        <h5path>/Constitutive/totalshear</h5path>
-        <category>Constitutive</category>
-        <description></description>
-    </totalshear>
-
-    <accumulatedshear_slip>
-        <type>Matrix</type>
-        <h5path>/Constitutive/accumulatedshear_slip</h5path>
-        <category>Constitutive</category>
-        <description>vector contains accumulated shear per slip system</description>
-    </accumulatedshear_slip>
-
-    <!--Derived section-->
-
-</h5ds>

lib/damask/__init__.py

@@ -1,31 +1,13 @@
 # -*- coding: UTF-8 no BOM -*-
 
 """Main aggregator"""
-import os,sys,time
-
-h5py_flag = os.path.join(os.path.dirname(__file__),'../../.noH5py')
-h5py_grace = 7200                         # only complain once every 7200 sec (2 hours)
-h5py_msg = "h5py module not found."
-
-now = time.time()
+import os
 
 with open(os.path.join(os.path.dirname(__file__),'../../VERSION')) as f:
   version = f.readline()[:-1]
 
 from .environment import Environment      # noqa
 from .asciitable  import ASCIItable       # noqa
-try:
-  from .h5table   import H5Table          # noqa
-  if os.path.exists(h5py_flag): os.remove(h5py_flag)  # delete flagging file on success
-except ImportError:
-  if os.path.exists(h5py_flag):
-    if now - os.path.getmtime(h5py_flag) > h5py_grace: # complain (again) every so-and-so often
-      sys.stderr.write(h5py_msg+'\n')
-      with open(h5py_flag, 'a'):
-        os.utime(h5py_flag,(now,now))     # update flag modification time to "now"
-  else:
-    open(h5py_flag, 'a').close()          # create flagging file
-    sys.stderr.write(h5py_msg+'\n')       # complain for the first time
     
 from .config      import Material         # noqa
 from .colormaps   import Colormap, Color  # noqa

lib/damask/h5table.py

@@ -1,146 +0,0 @@
-# -*- coding: UTF-8 no BOM -*-
-
-# ----------------------------------------------------------- #
-# Ideally the h5py should be enough to serve as the data      #
-# interface for future DAMASK, but since we are still not     #
-# sure when this major shift will happen, it seems to be a    #
-# good idea to provide a interface class that help user ease  #
-# into using HDF5 as the new daily storage driver.            #
-# ----------------------------------------------------------- #
-
-import os
-import h5py
-import numpy as np
-import xml.etree.cElementTree as ET
-
-# ---------------------------------------------------------------- #
-# python 3 has no unicode object, this ensures that the code works #
-# on Python 2&3                                                    #
-# ---------------------------------------------------------------- #
-try:
-    test = isinstance('test', unicode)
-except(NameError):
-    unicode = str
-
-
-def lables_to_path(label, dsXMLPath=None):
-    """Read the XML definition file and return the path."""
-    if dsXMLPath is None:
-        # use the default storage layout in DS_HDF5.xml
-        if "h5table.pyc" in __file__:
-            dsXMLPath = os.path.abspath(__file__).replace("h5table.pyc",
-                                                          "DS_HDF5.xml")
-        else:
-            dsXMLPath = os.path.abspath(__file__).replace("h5table.py",
-                                                          "DS_HDF5.xml")
-    # This current implementation requires that all variables
-    # stay under the root node, the nesting is defined through the
-    # h5path.
-    # Allow new derived data to be put under the root
-    tree = ET.parse(dsXMLPath)
-    try:
-        dataType = tree.find('{}/type'.format(label)).text
-        h5path = tree.find('{}/h5path'.format(label)).text
-    except:
-        dataType = "Scalar"
-        h5path = "/{}".format(label)  # just put it under root
-    return (dataType, h5path)
-
-
-class H5Table(object):
-    """
-    Lightweight interface class for h5py
-
-    DESCRIPTION
-    -----------
-        Interface/wrapper class for manipulating data in HDF5 with DAMASK
-        specialized data structure.
-        --> try to maintain a minimal API design.
-    PARAMETERS
-    ----------
-    h5f_path: str
-        Absolute path of the HDF5 file
-    METHOD
-    ------
-    del_entry()  -- Force delete attributes/group/datasets (dangerous)
-    get_attr()   -- Return attributes if possible
-    add_attr()   -- Add NEW attributes to dataset/group (no force overwrite)
-    get_data()   -- Retrieve data in numpy.ndarray
-    add_data()   -- Add dataset to H5 file
-    get_cmdlog() -- Return the command used to generate the data if possible
-    NOTE
-    ----
-        1. As an interface class, it uses the lazy evaluation design
-           that reads the data only when it is absolutely necessary.
-        2. The command line used to generate each new feature is stored with
-           each dataset as dataset attribute.
-
-    """
-
-    def __init__(self, h5f_path, new_file=False, dsXMLFile=None):
-        self.h5f_path = h5f_path
-        self.dsXMLFile = dsXMLFile
-        msg = 'Created by H5Talbe from DAMASK'
-        mode = 'w' if new_file else 'a'
-        with h5py.File(self.h5f_path, mode) as h5f:
-            h5f['/'].attrs['description'] = msg
-
-    def del_entry(self, feature_name):
-        """Delete entry in HDF5 table"""
-        dataType, h5f_path = lables_to_path(feature_name,
-                                            dsXMLPath=self.dsXMLFile)
-        with h5py.File(self.h5f_path, 'a') as h5f:
-            del h5f[h5f_path]
-
-    def get_attr(self, attr_name):
-        dataType, h5f_path = lables_to_path(attr_name,
-                                            dsXMLPath=self.dsXMLFile)
-        with h5py.File(self.h5f_path, 'a') as h5f:
-            rst_attr = h5f[h5f_path].attrs[attr_name]
-        return rst_attr
-
-    def add_attr(self, attr_name, attr_data):
-        dataType, h5f_path = lables_to_path(attr_name,
-                                            dsXMLPath=self.dsXMLFile)
-        with h5py.File(self.h5f_path, 'a') as h5f:
-            h5f[h5f_path].attrs[attr_name] = attr_data
-            h5f.flush()
-
-    def get_data(self, feature_name=None):
-        """Extract dataset from HDF5 table and return it in a numpy array"""
-        dataType, h5f_path = lables_to_path(feature_name,
-                                            dsXMLPath=self.dsXMLFile)
-        with h5py.File(self.h5f_path, 'a') as h5f:
-            h5f_dst = h5f[h5f_path]  # get the handle for target dataset(table)
-            rst_data = np.zeros(h5f_dst.shape)
-            h5f_dst.read_direct(rst_data)
-        return rst_data
-
-    def add_data(self, feature_name, dataset, cmd_log=None):
-        """Adding new feature into existing HDF5 file"""
-        dataType, h5f_path = lables_to_path(feature_name,
-                                            dsXMLPath=self.dsXMLFile)
-        with h5py.File(self.h5f_path, 'a') as h5f:
-            # NOTE:
-            # --> If dataset exists, delete the old one so as to write
-            #     a new one. For brand new dataset. For brand new one,
-            #     record its state as fresh in the cmd log.
-            try:
-                del h5f[h5f_path]
-                print("***deleting old {} from {}".format(feature_name,self.h5f_path))
-            except:
-                # if no cmd log, None will used
-                cmd_log = str(cmd_log) + " [FRESH]"
-            h5f.create_dataset(h5f_path, data=dataset)
-            # store the cmd in log is possible
-            if cmd_log is not None:
-                h5f[h5f_path].attrs['log'] = str(cmd_log)
-            h5f.flush()
-
-    def get_cmdlog(self, feature_name):
-        """Get cmd history used to generate the feature"""
-        dataType, h5f_path = lables_to_path(feature_name,
-                                            dsXMLPath=self.dsXMLFile)
-        with h5py.File(self.h5f_path, 'a') as h5f:
-            cmd_logs = h5f[h5f_path].attrs['log']
-        return cmd_logs

lib/damask/test/test.py

@@ -49,7 +49,8 @@ class Test():
 
     self.dirBase = os.path.dirname(os.path.realpath(sys.modules[self.__class__.__module__].__file__))
 
-    self.parser = OptionParser(description = '{} (Test class version: {})'.format(self.description,damask.version),
+    self.parser = OptionParser(option_class=damask.extendableOption,
+                               description = '{} (Test class version: {})'.format(self.description,damask.version),
                                usage = './test.py [options]')
     self.parser.add_option("-k", "--keep",
                            action = "store_true",
@@ -65,7 +66,8 @@ class Test():
                            help   = "show all test variants without actual calculation")
     self.parser.add_option("-s",   "--select",
                            dest   = "select",
-                           help   = "run test of given name only")
+                           action = 'extend', metavar = '<string LIST>',
+                           help   = "run test(s) of given name only")
     self.parser.set_defaults(keep   = self.keep,
                              accept = self.accept,
                              update = self.updateRequest,
@@ -90,7 +92,7 @@ class Test():
       if self.options.show:
         logging.critical('{}: {}'.format(variant+1,name))
       elif self.options.select is not None \
-           and not (name == self.options.select or str(variant+1) == self.options.select):
+           and not (name in self.options.select or str(variant+1) in self.options.select):
         pass
       else:
         try:
@@ -106,7 +108,7 @@ class Test():
             return variant+1                                                                          # return culprit
 
         except Exception as e:
-          logging.critical('exception during variant execution: "{}"'.format(e.message))
+          logging.critical('exception during variant execution: "{}"'.format(str(e)))
           return variant+1                                                                            # return culprit
     return 0
   
@@ -320,7 +322,9 @@ class Test():
     cur1Name = self.fileInCurrent(cur1)
     return self.compare_Array(cur0Name,cur1Name)
 
-  def compare_Table(self,headings0,file0,headings1,file1,normHeadings='',normType=None,
+  def compare_Table(self,headings0,file0,
+                         headings1,file1,
+                         normHeadings='',normType=None,
                          absoluteTolerance=False,perLine=False,skipLines=[]):
 
     import numpy as np
@@ -366,11 +370,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 1. table...\n'.format(key0))
+        raise Exception('column "{}" not found in first table...\n'.format(key0))
       elif key1 not in table1.labels(raw = True):
-        raise Exception('column {} not found in 2. table...\n'.format(key1))
+        raise Exception('column "{}" not found in second table...\n'.format(key1))
       elif normKey not in table0.labels(raw = True):
-        raise Exception('column {} not found in 1. table...\n'.format(normKey))
+        raise Exception('column "{}" not found in first table...\n'.format(normKey))
       else:
         column[0][i]  = table0.label_index(key0)
         column[1][i]  = table1.label_index(key1)
@@ -398,9 +402,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 1. table is 0.0, using absolute tolerance'.format(headings0[i]['label']))
+          logging.warning('At least one norm of "{}" in first table is 0.0, using absolute tolerance'.format(headings0[i]['label']))
         else:
-          logging.warning('Maximum norm of {} in 1. table is 0.0, using absolute tolerance'.format(headings0[i]['label']))
+          logging.warning('Maximum norm of "{}" in first table is 0.0, using absolute tolerance'.format(headings0[i]['label']))
 
     line1 = 0
     while table1.data_read():                                                  # read next data line of ASCII table
@@ -412,7 +416,7 @@ class Test():
                                                                                    norm[i][line1-len(skipLines)])
       line1 +=1
 
-    if (line0 != line1): raise Exception('found {} lines in 1. table but {} in 2. table'.format(line0,line1))
+    if (line0 != line1): raise Exception('found {} lines in first table but {} in second table'.format(line0,line1))
 
     logging.info(' ********')
     for i in range(dataLength):
@@ -559,25 +563,28 @@ class Test():
     return allclose
 
 
-  def compare_TableRefCur(self,headingsRef,ref,headingsCur='',cur='',normHeadings='',normType=None,\
+  def compare_TableRefCur(self,headingsRef,ref,headingsCur='',cur='',
+                               normHeadings='',normType=None,
                                absoluteTolerance=False,perLine=False,skipLines=[]):
 
-    if cur == '': cur = ref
-    if headingsCur == '': headingsCur = headingsRef
-    refName = self.fileInReference(ref)
-    curName = self.fileInCurrent(cur)
-    return self.compare_Table(headingsRef,refName,headingsCur,curName,normHeadings,normType,
+    return self.compare_Table(headingsRef,
+                              self.fileInReference(ref),
+                              headingsRef if headingsCur == '' else headingsCur,
+                              self.fileInCurrent(ref if cur == '' else cur),
+                              normHeadings,normType,
                               absoluteTolerance,perLine,skipLines)
 
 
-  def compare_TableCurCur(self,headingsCur0,Cur0,Cur1,headingsCur1='',normHeadings='',normType=None,\
+  def compare_TableCurCur(self,headingsCur0,Cur0,Cur1,
+                               headingsCur1='',
+                               normHeadings='',normType=None,
                                absoluteTolerance=False,perLine=False,skipLines=[]):
 
-    if headingsCur1 == '': headingsCur1 = headingsCur0
-    cur0Name = self.fileInCurrent(Cur0)
-    cur1Name = self.fileInCurrent(Cur1)
-    return self.compare_Table(headingsCur0,cur0Name,headingsCur1,cur1Name,normHeadings,normType,
-                                                               absoluteTolerance,perLine,skipLines)
+    return self.compare_Table(headingsCur0,
+                              self.fileInCurrent(Cur0),
+                              headingsCur0 if headingsCur1 == '' else headingsCur1,
+                              self.fileInCurrent(Cur1),
+                              normHeadings,normType,absoluteTolerance,perLine,skipLines)
 
 
   def report_Success(self,culprit):
@@ -585,13 +592,13 @@ class Test():
     ret = culprit
 
     if culprit == 0:
-      msg = 'The test passed.' if (self.options.select is not None or len(self.variants) == 1) \
-       else 'All {} tests passed.'.format(len(self.variants))
+      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)
     elif culprit == -1:
-      msg = 'Warning: Could not start test...'
+      msg = 'Warning: could not start test...'
       ret = 0
     else:
-      msg = ' * Test "{}" failed.'.format(self.variants[culprit-1])
+      msg = 'Test "{}" failed.'.format(self.variantName(culprit-1))
 
     logging.critical('\n'.join(['*'*40,msg,'*'*40]) + '\n')
     return ret

lib/damask/util.py

@@ -100,6 +100,18 @@ def execute(cmd,
   if process.returncode != 0: raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode))
   return out,error
 
+
+def coordGridAndSize(coordinates):
+  """Determines grid count and overall physical size along each dimension of an ordered array of coordinates"""
+  dim    = coordinates.shape[1]
+  coords = [np.unique(coordinates[:,i]) for i in range(dim)]
+  mincorner = np.array(map(min,coords))
+  maxcorner = np.array(map(max,coords))
+  grid   = np.array(map(len,coords),'i')
+  size   = grid/np.maximum(np.ones(dim,'d'), grid-1.0) * (maxcorner-mincorner)                      # size from edge to edge = dim * n/(n-1) 
+  size   = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1]))                             # spacing for grid==1 equal to smallest among other ones
+  return grid,size
+  
 # -----------------------------
 class extendableOption(Option):
   """
@@ -130,7 +142,7 @@ class backgroundMessage(threading.Thread):
              'hexagon':  ['⬢', '⬣'],
              'square':   ['▖', '▘', '▝', '▗'],
              'triangle': ['ᐊ', 'ᐊ', 'ᐃ', 'ᐅ', 'ᐅ', 'ᐃ'],
-             'amoeba':   ['▖', '▏', '▘', '▔', '▝', '▕', '▗', '▂'],
+             'amoeba':   ['▖', '▏', '▘', '▔', '▝', '▕', '▗', '▁'],
              'beat':     ['▁', '▂', '▃', '▅', '▆', '▇', '▇', '▆', '▅', '▃', '▂'],
              'prison':   ['ᚋ', 'ᚌ', 'ᚍ', 'ᚏ', 'ᚎ', 'ᚍ', 'ᚌ', 'ᚋ'],
              'breath':   ['ᚐ', 'ᚑ', 'ᚒ', 'ᚓ', 'ᚔ', 'ᚓ', 'ᚒ', 'ᚑ', 'ᚐ'],

processing/post/addCauchy.py

@@ -72,6 +72,7 @@ for name in filenames:
   table.head_write()
 
 # ------------------------------------------ process data ------------------------------------------
+
   outputAlive = True
   while outputAlive and table.data_read():                                                          # read next data line of ASCII table
     F = np.array(map(float,table.data[column[options.defgrad]:column[options.defgrad]+9]),'d').reshape(3,3)

processing/post/addCumulative.py

@@ -73,22 +73,17 @@ for name in filenames:
   table.head_write()
 
 # ------------------------------------------ process data ------------------------------------------ 
-
-  table.data_readArray()
-
   mask = []
   for col,dim in zip(columns,dims): mask += range(col,col+dim)                                      # isolate data columns to cumulate
+  cumulated = np.zeros(len(mask),dtype=float)                                                       # prepare output field
 
-  cumulated = np.zeros((len(table.data),len(mask)))                                                 # prepare output field
+  outputAlive = True
+  while outputAlive and table.data_read():                                                          # read next data line of ASCII table
+    for i,col in enumerate(mask):
+      cumulated[i] += float(table.data[col])                                                        # cumulate values
+    table.data_append(cumulated)
 
-  for i,values in enumerate(table.data[:,mask]):
-    cumulated[i,:] = cumulated[max(0,i-1),:] + values                                               # cumulate values
-  
-  table.data = np.hstack((table.data,cumulated))
-
-# ------------------------------------------ output result -----------------------------------------
-
-  table.data_writeArray()
+    outputAlive = table.data_write()                                                                # output processed line
 
 # ------------------------------------------ output finalization -----------------------------------  
 

processing/post/addCurl.py

@@ -9,25 +9,34 @@ import damask
 scriptName = os.path.splitext(os.path.basename(__file__))[0]
 scriptID   = ' '.join([scriptName,damask.version])
 
+def merge_dicts(*dict_args):
+  """Given any number of dicts, shallow copy and merge into a new dict, with precedence going to key value pairs in latter dicts."""
+  result = {}
+  for dictionary in dict_args:
+      result.update(dictionary)
+  return result
+
 def curlFFT(geomdim,field):
+  """Calculate curl of a vector or tensor field by transforming into Fourier space."""
   shapeFFT    = np.array(np.shape(field))[0:3]
   grid = np.array(np.shape(field)[2::-1])
   N = grid.prod()                                                                                    # field size
   n = np.array(np.shape(field)[3:]).prod()                                                           # data size
 
- if   n == 3:   dataType = 'vector'
- elif n == 9:   dataType = 'tensor'
-
   field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT)
   curl_fourier  = np.empty(field_fourier.shape,'c16')
 
   # differentiation in Fourier space
   TWOPIIMG = 2.0j*math.pi
+  einsums = { 
+              3:'slm,ijkl,ijkm->ijks',                                                               # vector, 3 -> 3
+              9:'slm,ijkl,ijknm->ijksn',                                                             # tensor, 3x3 -> 3x3
+            }
   k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0]
- if grid[2]%2 == 0: k_sk[grid[2]//2] = 0                                                            # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
+  if grid[2]%2 == 0: k_sk[grid[2]//2] = 0                                                            # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
 
   k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1]
- if grid[1]%2 == 0: k_sj[grid[1]//2] = 0                                                            # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
+  if grid[1]%2 == 0: k_sj[grid[1]//2] = 0                                                            # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
 
   k_si = np.arange(grid[0]//2+1)/geomdim[2]
 
@@ -38,10 +47,7 @@ def curlFFT(geomdim,field):
   e[0, 1, 2] = e[1, 2, 0] = e[2, 0, 1] = 1.0                                                         # Levi-Civita symbols 
   e[0, 2, 1] = e[2, 1, 0] = e[1, 0, 2] = -1.0
 
- if dataType == 'tensor':                                                                           # tensor, 3x3 -> 3x3 
-   curl_fourier = np.einsum('slm,ijkl,ijknm->ijksn',e,k_s,field_fourier)*TWOPIIMG
- elif dataType == 'vector':                                                                         # vector, 3 -> 3
-   curl_fourier = np.einsum('slm,ijkl,ijkm->ijks',e,k_s,field_fourier)*TWOPIIMG
+  curl_fourier = np.einsum(einsums[n],e,k_s,field_fourier)*TWOPIIMG
 
   return np.fft.irfftn(curl_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n])
 
@@ -52,31 +58,37 @@ def curlFFT(geomdim,field):
 
 parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """
 Add column(s) containing curl of requested column(s).
-Operates on periodic ordered three-dimensional data sets.
-Deals with both vector- and tensor fields.
-
+Operates on periodic ordered three-dimensional data sets
+of vector and tensor fields.
 """, version = scriptID)
 
 parser.add_option('-p','--pos','--periodiccellcenter',
                   dest = 'pos',
                   type = 'string', metavar = 'string',
                   help = 'label of coordinates [%default]')
-parser.add_option('-v','--vector',
-                  dest = 'vector',
+parser.add_option('-l','--label',
+                  dest = 'data',
                   action = 'extend', metavar = '<string LIST>',
-                  help = 'label(s) of vector field values')
-parser.add_option('-t','--tensor',
-                  dest = 'tensor',
-                  action = 'extend', metavar = '<string LIST>',
-                  help = 'label(s) of tensor field values')
+                  help = 'label(s) of field values')
 
 parser.set_defaults(pos = 'pos',
                    )
 
+
 (options,filenames) = parser.parse_args()
 
-if options.vector is None and options.tensor is None:
-  parser.error('no data column specified.')
+if options.data is None: parser.error('no data column specified.')
+
+# --- define possible data types -------------------------------------------------------------------
+
+datatypes = {
+              3: {'name': 'vector',
+                  'shape': [3],
+                 },
+              9: {'name': 'tensor',
+                  'shape': [3,3],
+                 },
+            }
 
 # --- loop over input files ------------------------------------------------------------------------
 
@@ -87,30 +99,27 @@ for name in filenames:
   except: continue
   damask.util.report(scriptName,name)
 
-# ------------------------------------------ read header ------------------------------------------
+# --- interpret header ----------------------------------------------------------------------------
 
   table.head_read()
 
-# ------------------------------------------ sanity checks ----------------------------------------
-
-  items = {
-            'tensor': {'dim': 9, 'shape': [3,3], 'labels':options.tensor, 'active':[], 'column': []},
-            'vector': {'dim': 3, 'shape': [3],   'labels':options.vector, 'active':[], 'column': []},
-          }
-  errors  = []
   remarks = []
-  column = {}
+  errors  = []
+  active = []
 
-  if table.label_dimension(options.pos) != 3: errors.append('coordinates {} are not a vector.'.format(options.pos))
-  else: colCoord = table.label_index(options.pos)
+  coordDim = table.label_dimension(options.pos)
+  if coordDim != 3:
+    errors.append('coordinates "{}" must be three-dimensional.'.format(options.pos))
+  else: coordCol = table.label_index(options.pos)
 
-  for type, data in items.iteritems():
-    for what in (data['labels'] if data['labels'] is not None else []):
-      dim = table.label_dimension(what)
-      if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type))
+  for me in options.data:
+    dim = table.label_dimension(me)
+    if dim in datatypes:
+      active.append(merge_dicts({'label':me},datatypes[dim]))
+      remarks.append('differentiating {} "{}"...'.format(datatypes[dim]['name'],me))
     else:
-        items[type]['active'].append(what)
-        items[type]['column'].append(table.label_index(what))
+      remarks.append('skipping "{}" of dimension {}...'.format(me,dim) if dim != -1 else \
+                     '"{}" not found...'.format(me) )
 
   if remarks != []: damask.util.croak(remarks)
   if errors  != []:
@@ -121,30 +130,24 @@ for name in filenames:
 # ------------------------------------------ assemble header --------------------------------------
 
   table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
-  for type, data in items.iteritems():
-    for label in data['active']:
-      table.labels_append(['{}_curlFFT({})'.format(i+1,label) for i in range(data['dim'])])         # extend ASCII header with new labels
+  for data in active:
+    table.labels_append(['{}_curlFFT({})'.format(i+1,data['label']) 
+                        for i in range(np.prod(np.array(data['shape'])))])                         # extend ASCII header with new labels
   table.head_write()
 
 # --------------- figure out size and grid ---------------------------------------------------------
 
   table.data_readArray()
 
-  coords = [np.unique(table.data[:,colCoord+i]) for i in range(3)]
-  mincorner = np.array(map(min,coords))
-  maxcorner = np.array(map(max,coords))
-  grid   = np.array(map(len,coords),'i')
-  size   = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner)                        # size from edge to edge = dim * n/(n-1) 
-  size   = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1]))                             # spacing for grid==1 equal to smallest among other ones
+  grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)])
 
 # ------------------------------------------ process value field -----------------------------------
 
   stack = [table.data]
-  for type, data in items.iteritems():
-    for i,label in enumerate(data['active']):
+  for data in active:
     # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation
     stack.append(curlFFT(size[::-1],
-                           table.data[:,data['column'][i]:data['column'][i]+data['dim']].
+                         table.data[:,table.label_indexrange(data['label'])].
                          reshape(grid[::-1].tolist()+data['shape'])))
 
 # ------------------------------------------ output result -----------------------------------------

processing/post/addDerivative.py

@@ -0,0 +1,121 @@
+#!/usr/bin/env python2.7
+# -*- coding: UTF-8 no BOM -*-
+
+import os,sys
+import numpy as np
+from optparse import OptionParser
+import damask
+
+scriptName = os.path.splitext(os.path.basename(__file__))[0]
+scriptID   = ' '.join([scriptName,damask.version])
+
+def derivative(coordinates,what):
+  
+  result = np.empty_like(what)
+  
+  # use differentiation by interpolation
+  # as described in http://www2.math.umd.edu/~dlevy/classes/amsc466/lecture-notes/differentiation-chap.pdf
+
+  result[1:-1,:] = + what[1:-1,:] * (2.*coordinates[1:-1]-coordinates[:-2]-coordinates[2:]) / \
+                     ((coordinates[1:-1]-coordinates[:-2])*(coordinates[1:-1]-coordinates[2:])) \
+                   + what[2:,:] * (coordinates[1:-1]-coordinates[:-2]) / \
+                     ((coordinates[2:]-coordinates[1:-1])*(coordinates[2:]-coordinates[:-2])) \
+                   + what[:-2,:] * (coordinates[1:-1]-coordinates[2:]) / \
+                     ((coordinates[:-2]-coordinates[1:-1])*(coordinates[:-2]-coordinates[2:])) \
+
+  result[0,:]    = (what[0,:] - what[1,:]) / \
+                   (coordinates[0] - coordinates[1])
+  result[-1,:]   = (what[-1,:] - what[-2,:]) / \
+                   (coordinates[-1] - coordinates[-2])
+
+  return result
+  
+# --------------------------------------------------------------------
+#                                MAIN
+# --------------------------------------------------------------------
+
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+Add column(s) containing numerical derivative of requested column(s) with respect to given coordinates.
+
+""", version = scriptID)
+
+parser.add_option('-c','--coordinates',
+                  dest = 'coordinates',
+                  type = 'string', metavar='string',
+                  help = 'heading of coordinate column')
+parser.add_option('-l','--label',
+                  dest = 'label',
+                  action = 'extend', metavar = '<string LIST>',
+                  help = 'heading of column(s) to differentiate')
+
+
+(options,filenames) = parser.parse_args()
+
+if options.coordinates is None:
+  parser.error('no coordinate column specified.')
+if options.label is None:
+  parser.error('no data column specified.')
+
+# --- loop over input files -------------------------------------------------------------------------
+
+if filenames == []: filenames = [None]
+
+for name in filenames:
+  try:    table = damask.ASCIItable(name = name,
+                                    buffered = False)
+  except: continue
+  damask.util.report(scriptName,name)
+
+# ------------------------------------------ read header ------------------------------------------
+
+  table.head_read()
+
+# ------------------------------------------ sanity checks ----------------------------------------
+
+  errors  = []
+  remarks = []
+  columns = []
+  dims    = []
+  
+  if table.label_dimension(options.coordinates) != 1:
+    errors.append('coordinate column {} is not scalar.'.format(options.coordinates))
+
+  for what in options.label:
+    dim = table.label_dimension(what)
+    if dim < 0: remarks.append('column {} not found...'.format(what))
+    else:
+      dims.append(dim)
+      columns.append(table.label_index(what))
+      table.labels_append('d({})/d({})'.format(what,options.coordinates) if dim == 1 else
+                         ['{}_d({})/d({})'.format(i+1,what,options.coordinates) for i in range(dim)]  )  # extend ASCII header with new labels
+
+  if remarks != []: damask.util.croak(remarks)
+  if errors  != []:
+    damask.util.croak(errors)
+    table.close(dismiss = True)
+    continue
+
+# ------------------------------------------ assemble header --------------------------------------
+
+  table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
+  table.head_write()
+
+# ------------------------------------------ process data ------------------------------------------
+
+  table.data_readArray()
+
+  mask = []
+  for col,dim in zip(columns,dims): mask += range(col,col+dim)                                      # isolate data columns to differentiate
+  
+  differentiated = derivative(table.data[:,table.label_index(options.coordinates)].reshape((len(table.data),1)),
+                              table.data[:,mask])                                                   # calculate numerical derivative
+  
+  table.data = np.hstack((table.data,differentiated))
+
+# ------------------------------------------ output result -----------------------------------------
+
+  table.data_writeArray()
+
+# ------------------------------------------ output finalization -----------------------------------  
+
+  table.close()                                                                                     # close ASCII tables

processing/post/addDivergence.py

@@ -9,34 +9,41 @@ import damask
 scriptName = os.path.splitext(os.path.basename(__file__))[0]
 scriptID   = ' '.join([scriptName,damask.version])
 
+def merge_dicts(*dict_args):
+  """Given any number of dicts, shallow copy and merge into a new dict, with precedence going to key value pairs in latter dicts."""
+  result = {}
+  for dictionary in dict_args:
+      result.update(dictionary)
+  return result
+
 def divFFT(geomdim,field):
+  """Calculate divergence of a vector or tensor field by transforming into Fourier space."""
   shapeFFT    = np.array(np.shape(field))[0:3]
   grid = np.array(np.shape(field)[2::-1])
   N = grid.prod()                                                                                    # field size
   n = np.array(np.shape(field)[3:]).prod()                                                           # data size
 
- if   n == 3:   dataType = 'vector'
- elif n == 9:   dataType = 'tensor'
-
   field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT)
   div_fourier   = np.empty(field_fourier.shape[0:len(np.shape(field))-1],'c16')
 
   # differentiation in Fourier space
   TWOPIIMG = 2.0j*math.pi
+  einsums = { 
+              3:'ijkl,ijkl->ijk',                                                               # vector, 3 -> 1
+              9:'ijkm,ijklm->ijkl',                                                             # tensor, 3x3 -> 3
+            }
   k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0]
- if grid[2]%2 == 0: k_sk[grid[2]//2] = 0                                                            # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
+  if grid[2]%2 == 0: k_sk[grid[2]//2] = 0                                                            # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
 
   k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1]
- if grid[1]%2 == 0: k_sj[grid[1]//2] = 0                                                            # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
+  if grid[1]%2 == 0: k_sj[grid[1]//2] = 0                                                            # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
 
   k_si = np.arange(grid[0]//2+1)/geomdim[2]
 
   kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij')
   k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')                           
- if dataType == 'tensor':                                                                           # tensor, 3x3 -> 3
-   div_fourier = np.einsum('ijklm,ijkm->ijkl',field_fourier,k_s)*TWOPIIMG
- elif dataType == 'vector':                                                                         # vector, 3 -> 1
-   div_fourier = np.einsum('ijkl,ijkl->ijk',field_fourier,k_s)*TWOPIIMG
+
+  div_fourier = np.einsum(einsums[n],k_s,field_fourier)*TWOPIIMG
 
   return np.fft.irfftn(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n/3])
 
@@ -46,32 +53,38 @@ def divFFT(geomdim,field):
 # --------------------------------------------------------------------
 
 parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """
-Add column(s) containing divergence of requested column(s).
-Operates on periodic ordered three-dimensional data sets.
-Deals with both vector- and tensor-valued fields.
-
+Add column(s) containing curl of requested column(s).
+Operates on periodic ordered three-dimensional data sets
+of vector and tensor fields.
 """, version = scriptID)
 
 parser.add_option('-p','--pos','--periodiccellcenter',
                   dest = 'pos',
                   type = 'string', metavar = 'string',
                   help = 'label of coordinates [%default]')
-parser.add_option('-v','--vector',
-                  dest = 'vector',
+parser.add_option('-l','--label',
+                  dest = 'data',
                   action = 'extend', metavar = '<string LIST>',
-                  help = 'label(s) of vector field values')
-parser.add_option('-t','--tensor',
-                  dest = 'tensor',
-                  action = 'extend', metavar = '<string LIST>',
-                  help = 'label(s) of tensor field values')
+                  help = 'label(s) of field values')
 
 parser.set_defaults(pos = 'pos',
                    )
 
+
 (options,filenames) = parser.parse_args()
 
-if options.vector is None and options.tensor is None:
-  parser.error('no data column specified.')
+if options.data is None: parser.error('no data column specified.')
+
+# --- define possible data types -------------------------------------------------------------------
+
+datatypes = {
+              3: {'name': 'vector',
+                  'shape': [3],
+                 },
+              9: {'name': 'tensor',
+                  'shape': [3,3],
+                 },
+            }
 
 # --- loop over input files ------------------------------------------------------------------------
 
@@ -82,30 +95,27 @@ for name in filenames:
   except: continue
   damask.util.report(scriptName,name)
 
-# ------------------------------------------ read header ------------------------------------------
+# --- interpret header ----------------------------------------------------------------------------
 
   table.head_read()
 
-# ------------------------------------------ sanity checks ----------------------------------------
-
-  items = {
-            'tensor': {'dim': 9, 'shape': [3,3], 'labels':options.tensor, 'active':[], 'column': []},
-            'vector': {'dim': 3, 'shape': [3],   'labels':options.vector, 'active':[], 'column': []},
-          }
-  errors  = []
   remarks = []
-  column = {}
+  errors  = []
+  active = []
 
-  if table.label_dimension(options.pos) != 3: errors.append('coordinates {} are not a vector.'.format(options.pos))
-  else: colCoord = table.label_index(options.pos)
+  coordDim = table.label_dimension(options.pos)
+  if coordDim != 3:
+    errors.append('coordinates "{}" must be three-dimensional.'.format(options.pos))
+  else: coordCol = table.label_index(options.pos)
 
-  for type, data in items.iteritems():
-    for what in (data['labels'] if data['labels'] is not None else []):
-      dim = table.label_dimension(what)
-      if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type))
+  for me in options.data:
+    dim = table.label_dimension(me)
+    if dim in datatypes:
+      active.append(merge_dicts({'label':me},datatypes[dim]))
+      remarks.append('differentiating {} "{}"...'.format(datatypes[dim]['name'],me))
     else:
-        items[type]['active'].append(what)
-        items[type]['column'].append(table.label_index(what))
+      remarks.append('skipping "{}" of dimension {}...'.format(me,dim) if dim != -1 else \
+                     '"{}" not found...'.format(me) )
 
   if remarks != []: damask.util.croak(remarks)
   if errors  != []:
@@ -116,31 +126,25 @@ for name in filenames:
 # ------------------------------------------ assemble header --------------------------------------
 
   table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
-  for type, data in items.iteritems():
-    for label in data['active']:
-      table.labels_append(['divFFT({})'.format(label) if type == 'vector' else
-                           '{}_divFFT({})'.format(i+1,label) for i in range(data['dim']//3)])       # extend ASCII header with new labels
+  for data in active:
+    table.labels_append(['divFFT({})'.format(data['label']) if data['shape'] == [3] \
+                 else '{}_divFFT({})'.format(i+1,data['label']) 
+                 for i in range(np.prod(np.array(data['shape']))//3)])                              # extend ASCII header with new labels
   table.head_write()
 
 # --------------- figure out size and grid ---------------------------------------------------------
 
   table.data_readArray()
 
-  coords = [np.unique(table.data[:,colCoord+i]) for i in range(3)]
-  mincorner = np.array(map(min,coords))
-  maxcorner = np.array(map(max,coords))
-  grid   = np.array(map(len,coords),'i')
-  size   = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner)                        # size from edge to edge = dim * n/(n-1) 
-  size   = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1]))                             # spacing for grid==1 equal to smallest among other ones
+  grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)])
 
 # ------------------------------------------ process value field -----------------------------------
 
   stack = [table.data]
-  for type, data in items.iteritems():
-    for i,label in enumerate(data['active']):
+  for data in active:
     # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation
     stack.append(divFFT(size[::-1],
-                          table.data[:,data['column'][i]:data['column'][i]+data['dim']].
+                        table.data[:,table.label_indexrange(data['label'])].
                         reshape(grid[::-1].tolist()+data['shape'])))
 
 # ------------------------------------------ output result -----------------------------------------

processing/post/addGaussian.py

@@ -18,7 +18,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """
 Add column(s) containing Gaussian filtered values of requested column(s).
 Operates on periodic and non-periodic ordered three-dimensional data sets.
-For Details see scipy.ndimage documentation.
+For details see scipy.ndimage documentation.
 
 """, version = scriptID)
 
@@ -43,15 +43,14 @@ parser.add_option('--sigma',
 parser.add_option('--periodic',
                   dest = 'periodic',
                   action = 'store_true',
-                  help = 'assume periodic grain structure'
-                 )
+                  help = 'assume periodic grain structure')
 
 
 
 parser.set_defaults(pos = 'pos',
                     order = 0,
                     sigma = 1,
-                    periodic = False
+                    periodic = False,
                    )
 
 (options,filenames) = parser.parse_args()
@@ -110,12 +109,7 @@ for name in filenames:
 
   table.data_readArray()
 
-  coords = [np.unique(table.data[:,colCoord+i]) for i in range(3)]
-  mincorner = np.array(map(min,coords))
-  maxcorner = np.array(map(max,coords))
-  grid   = np.array(map(len,coords),'i')
-  size   = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner)                        # size from edge to edge = dim * n/(n-1) 
-  size   = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1]))                             # spacing for grid==1 equal to smallest among other ones
+  grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)])
 
 # ------------------------------------------ process value field -----------------------------------
 

processing/post/addGradient.py

@@ -9,34 +9,41 @@ import damask
 scriptName = os.path.splitext(os.path.basename(__file__))[0]
 scriptID   = ' '.join([scriptName,damask.version])
 
+def merge_dicts(*dict_args):
+  """Given any number of dicts, shallow copy and merge into a new dict, with precedence going to key value pairs in latter dicts."""
+  result = {}
+  for dictionary in dict_args:
+      result.update(dictionary)
+  return result
+
 def gradFFT(geomdim,field):
+  """Calculate gradient of a vector or scalar field by transforming into Fourier space."""
   shapeFFT = np.array(np.shape(field))[0:3]
   grid     = np.array(np.shape(field)[2::-1])
   N = grid.prod()                                                                                    # field size
   n = np.array(np.shape(field)[3:]).prod()                                                           # data size
 
- if   n == 3:   dataType = 'vector'
- elif n == 1:   dataType = 'scalar'
-
   field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT)
   grad_fourier  = np.empty(field_fourier.shape+(3,),'c16')
 
   # differentiation in Fourier space
   TWOPIIMG = 2.0j*math.pi
+  einsums = { 
+              1:'ijkl,ijkm->ijkm',                                                                   # scalar, 1 -> 3
+              3:'ijkl,ijkm->ijklm',                                                                  # vector, 3 -> 3x3
+            }
+
   k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0]
- if grid[2]%2 == 0: k_sk[grid[2]//2] = 0                                                            # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
+  if grid[2]%2 == 0: k_sk[grid[2]//2] = 0                                                            # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
 
   k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1]
- if grid[1]%2 == 0: k_sj[grid[1]//2] = 0                                                            # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
+  if grid[1]%2 == 0: k_sj[grid[1]//2] = 0                                                            # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
 
   k_si = np.arange(grid[0]//2+1)/geomdim[2]
 
   kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij')
   k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')                           
- if dataType == 'vector':                                                                           # vector, 3 -> 3x3
-   grad_fourier = np.einsum('ijkl,ijkm->ijklm',field_fourier,k_s)*TWOPIIMG
- elif dataType == 'scalar':                                                                         # scalar, 1 -> 3
-   grad_fourier = np.einsum('ijkl,ijkl->ijkl',field_fourier,k_s)*TWOPIIMG
+  grad_fourier = np.einsum(einsums[n],field_fourier,k_s)*TWOPIIMG
 
   return np.fft.irfftn(grad_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,3*n])
 
@@ -47,8 +54,8 @@ def gradFFT(geomdim,field):
 
 parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """
 Add column(s) containing gradient of requested column(s).
-Operates on periodic ordered three-dimensional data sets.
-Deals with both vector- and scalar fields.
+Operates on periodic ordered three-dimensional data sets
+of vector and scalar fields.
 
 """, version = scriptID)
 
@@ -56,22 +63,28 @@ parser.add_option('-p','--pos','--periodiccellcenter',
                   dest = 'pos',
                   type = 'string', metavar = 'string',
                   help = 'label of coordinates [%default]')
-parser.add_option('-v','--vector',
-                  dest = 'vector',
+parser.add_option('-l','--label',
+                  dest = 'data',
                   action = 'extend', metavar = '<string LIST>',
-                  help = 'label(s) of vector field values')
-parser.add_option('-s','--scalar',
-                  dest = 'scalar',
-                  action = 'extend', metavar = '<string LIST>',
-                  help = 'label(s) of scalar field values')
+                  help = 'label(s) of field values')
 
 parser.set_defaults(pos = 'pos',
                    )
 
 (options,filenames) = parser.parse_args()
 
-if options.vector is None and options.scalar is None:
-  parser.error('no data column specified.')
+if options.data is None: parser.error('no data column specified.')
+
+# --- define possible data types -------------------------------------------------------------------
+
+datatypes = {
+              1: {'name': 'scalar',
+                  'shape': [1],
+                 },
+              3: {'name': 'vector',
+                  'shape': [3],
+                 },
+            }
 
 # --- loop over input files ------------------------------------------------------------------------
 
@@ -82,30 +95,27 @@ for name in filenames:
   except: continue
   damask.util.report(scriptName,name)
 
-# ------------------------------------------ read header ------------------------------------------
+# --- interpret header ----------------------------------------------------------------------------
 
   table.head_read()
 
-# ------------------------------------------ sanity checks ----------------------------------------
-
-  items = {
-            'scalar': {'dim': 1, 'shape': [1], 'labels':options.scalar, 'active':[], 'column': []},
-            'vector': {'dim': 3, 'shape': [3], 'labels':options.vector, 'active':[], 'column': []},
-          }
-  errors  = []
   remarks = []
-  column = {}
+  errors  = []
+  active = []
 
-  if table.label_dimension(options.pos) != 3: errors.append('coordinates {} are not a vector.'.format(options.pos))
-  else: colCoord = table.label_index(options.pos)
+  coordDim = table.label_dimension(options.pos)
+  if coordDim != 3:
+    errors.append('coordinates "{}" must be three-dimensional.'.format(options.pos))
+  else: coordCol = table.label_index(options.pos)
 
-  for type, data in items.iteritems():
-    for what in (data['labels'] if data['labels'] is not None else []):
-      dim = table.label_dimension(what)
-      if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type))
+  for me in options.data:
+    dim = table.label_dimension(me)
+    if dim in datatypes:
+      active.append(merge_dicts({'label':me},datatypes[dim]))
+      remarks.append('differentiating {} "{}"...'.format(datatypes[dim]['name'],me))
     else:
-        items[type]['active'].append(what)
-        items[type]['column'].append(table.label_index(what))
+      remarks.append('skipping "{}" of dimension {}...'.format(me,dim) if dim != -1 else \
+                     '"{}" not found...'.format(me) )
 
   if remarks != []: damask.util.croak(remarks)
   if errors  != []:
@@ -116,30 +126,24 @@ for name in filenames:
 # ------------------------------------------ assemble header --------------------------------------
 
   table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
-  for type, data in items.iteritems():
-    for label in data['active']:
-      table.labels_append(['{}_gradFFT({})'.format(i+1,label) for i in range(3 * data['dim'])])     # extend ASCII header with new labels
+  for data in active:
+    table.labels_append(['{}_gradFFT({})'.format(i+1,data['label']) 
+                        for i in range(coordDim*np.prod(np.array(data['shape'])))])                 # extend ASCII header with new labels
   table.head_write()
 
 # --------------- figure out size and grid ---------------------------------------------------------
 
   table.data_readArray()
 
-  coords = [np.unique(table.data[:,colCoord+i]) for i in range(3)]
-  mincorner = np.array(map(min,coords))
-  maxcorner = np.array(map(max,coords))
-  grid   = np.array(map(len,coords),'i')
-  size   = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner)                        # size from edge to edge = dim * n/(n-1) 
-  size   = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1]))                             # spacing for grid==1 equal to smallest among other ones
+  grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)])
 
 # ------------------------------------------ process value field -----------------------------------
 
   stack = [table.data]
-  for type, data in items.iteritems():
-    for i,label in enumerate(data['active']):
+  for data in active:
     # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation
     stack.append(gradFFT(size[::-1],
-                           table.data[:,data['column'][i]:data['column'][i]+data['dim']].
+                         table.data[:,table.label_indexrange(data['label'])].
                          reshape(grid[::-1].tolist()+data['shape'])))
 
 # ------------------------------------------ output result -----------------------------------------

processing/post/addSpectralDecomposition.py

@@ -49,7 +49,8 @@ for name in filenames:
 
   table.head_read()
 
-# ------------------------------------------ sanity checks ----------------------------------------
+
+# ------------------------------------------ assemble header 1 ------------------------------------
 
   items = {
             'tensor': {'dim': 9, 'shape': [3,3], 'labels':options.tensor, 'column': []},
@@ -74,12 +75,12 @@ for name in filenames:
     table.close(dismiss = True)
     continue
 
-# ------------------------------------------ assemble header --------------------------------------
+# ------------------------------------------ assemble header 2 ------------------------------------
 
   table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
   table.head_write()
 
-# ------------------------------------------ process data ------------------------------------------
+# ------------------------------------------ process data -----------------------------------------
 
   outputAlive = True
   while outputAlive and table.data_read():                                                          # read next data line of ASCII table

processing/post/h5_addCalculation.py

@@ -1,72 +0,0 @@
-#!/usr/bin/env python2.7
-# -*- coding: UTF-8 no BOM -*-
-
-import os
-# import re
-# import sys
-import collections
-# import math
-import damask
-# import numpy as np
-from optparse import OptionParser
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID = ' '.join([scriptName, damask.version])
-
-
-# ----- Helper functions ----- #
-def listify(x):
-    return x if isinstance(x, collections.Iterable) else [x]
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-usageEx = """
-usage_in_details:
-    Column labels are tagged by '#label#' in formulas.
-    Use ';' for ',' in functions. Numpy is available as 'np'.
-    Special variables: #_row_# -- row index
-
-    Examples:
-    (1) magnitude of vector -- "np.linalg.norm(#vec#)"
-    (2) rounded root of row number -- "round(math.sqrt(#_row_#);3)"
-"""
-desp = "Add or alter column(s) with derived values according to "
-desp += "user-defined arithmetic operation between column(s)."
-
-parser = OptionParser(option_class=damask.extendableOption,
-                      usage='%prog options [file[s]]' + usageEx,
-                      description=desp,
-                      version=scriptID)
-parser.add_option('-l', '--label',
-                  dest='labels',
-                  action='extend', metavar='<string LIST>',
-                  help='(list of) new column labels')
-parser.add_option('-f', '--formula',
-                  dest='formulas',
-                  action='extend', metavar='<string LIST>',
-                  help='(list of) formulas corresponding to labels')
-parser.add_option('-c', '--condition',
-                  dest='condition', metavar='string',
-                  help='condition to filter rows')
-
-parser.set_defaults(condition=None)
-
-(options, filenames) = parser.parse_args()
-
-# ----- parse formulas ----- #
-for i in range(len(options.formulas)):
-    options.formulas[i] = options.formulas[i].replace(';', ',')
-
-# ----- loop over input files ----- #
-for name in filenames:
-    try:
-        h5f = damask.H5Table(name, new_file=False)
-    except:
-        print("!!!Cannot process {}".format(name))
-        continue
-    damask.util.report(scriptName, name)
-
-# Note:
-# --> not immediately needed, come back later

processing/post/h5_addCauchy.py

@@ -1,61 +0,0 @@
-#!/usr/bin/env python2.7
-# -*- coding: UTF-8 no BOM -*-
-
-import os
-import damask
-import numpy as np
-from optparse import OptionParser
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID = ' '.join([scriptName, damask.version])
-
-
-def getCauchy(f, p):
-    """Return Cauchy stress for given f and p"""
-    # [Cauchy] = (1/det(F)) * [P].[F_transpose]
-    f = f.reshape((3, 3))
-    p = p.reshape((3, 3))
-    return 1.0/np.linalg.det(f)*np.dot(p, f.T).reshape(9)
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-desp = "Add column(s) containing Cauchy stress based on given column(s)"
-desp += "of deformation gradient and first Piola--Kirchhoff stress."
-parser = OptionParser(option_class=damask.extendableOption,
-                      usage='%prog options [file[s]]',
-                      description=desp,
-                      version=scriptID)
-parser.add_option('-f', '--defgrad',
-                  dest='defgrad',
-                  type='string', metavar='string',
-                  help='heading for deformation gradient [%default]')
-parser.add_option('-p', '--stress',
-                  dest='stress',
-                  type='string', metavar='string',
-                  help='heading for first Piola--Kirchhoff stress [%default]')
-
-parser.set_defaults(defgrad='f',
-                    stress='p')
-
-(options, filenames) = parser.parse_args()
-
-# ----- loop over input H5 files ----- #
-for name in filenames:
-    try:
-        h5f = damask.H5Table(name, new_file=False)
-    except:
-        continue
-    damask.util.report(scriptName, name)
-
-    # ----- read in data ----- #
-    f = h5f.get_data("f")
-    p = h5f.get_data("p")
-
-    # ----- calculate Cauchy stress ----- #
-    cauchy = [getCauchy(f_i, p_i) for f_i, p_i in zip(f, p)]
-
-    # ----- write to HDF5 file ----- #
-    cmd_log = " ".join([scriptID, name])
-    h5f.add_data('Cauchy', np.array(cauchy), cmd_log=cmd_log)

processing/post/h5_addIPFcolor.py

@@ -1,145 +0,0 @@
-#!/usr/bin/env python2.7
-# -*- coding: UTF-8 no BOM -*-
-
-import os
-import sys
-import math
-import damask
-import numpy as np
-from optparse import OptionParser
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID = ' '.join([scriptName, damask.version])
-
-# TODO
-#   This implementation will have to iterate through the array one
-#   element at a time, maybe there are some other ways to make this
-#   faster.
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-desp = "Add RGB color value corresponding to TSL-OIM scheme for IPF."
-parser = OptionParser(option_class=damask.extendableOption,
-                      usage='%prog options [file[s]]',
-                      description=desp,
-                      version=scriptID)
-parser.add_option('-p', '--pole',
-                  dest='pole',
-                  type='float', nargs=3, metavar='float float float',
-                  help='lab frame direction for IPF [%default]')
-msg = ', '.join(damask.Symmetry.lattices[1:])
-parser.add_option('-s', '--symmetry',
-                  dest='symmetry',
-                  type='choice', choices=damask.Symmetry.lattices[1:],
-                  metavar='string',
-                  help='crystal symmetry [%default] {{{}}} '.format(msg))
-parser.add_option('-e', '--eulers',
-                  dest='eulers',
-                  type='string', metavar='string',
-                  help='Euler angles label')
-parser.add_option('-d', '--degrees',
-                  dest='degrees',
-                  action='store_true',
-                  help='Euler angles are given in degrees [%default]')
-parser.add_option('-m', '--matrix',
-                  dest='matrix',
-                  type='string', metavar='string',
-                  help='orientation matrix label')
-parser.add_option('-a',
-                  dest='a',
-                  type='string', metavar='string',
-                  help='crystal frame a vector label')
-parser.add_option('-b',
-                  dest='b',
-                  type='string', metavar='string',
-                  help='crystal frame b vector label')
-parser.add_option('-c',
-                  dest='c',
-                  type='string', metavar='string',
-                  help='crystal frame c vector label')
-parser.add_option('-q', '--quaternion',
-                  dest='quaternion',
-                  type='string', metavar='string',
-                  help='quaternion label')
-
-parser.set_defaults(pole=(0.0, 0.0, 1.0),
-                    symmetry=damask.Symmetry.lattices[-1],
-                    degrees=False)
-
-(options, filenames) = parser.parse_args()
-
-# safe guarding to have only one orientation representation
-# use dynamic typing to group a,b,c into frame
-options.frame = [options.a, options.b, options.c]
-input = [options.eulers is not None,
-         all(options.frame),
-         options.matrix is not None,
-         options.quaternion is not None]
-
-if np.sum(input) != 1:
-    parser.error('needs exactly one input format.')
-
-# select input label that was requested (active)
-label_active = np.where(input)[0][0]
-(label, dim, inputtype) = [(options.eulers, 3, 'eulers'),
-                           (options.frame, [3, 3, 3], 'frame'),
-                           (options.matrix, 9, 'matrix'),
-                           (options.quaternion, 4, 'quaternion')][label_active]
-
-# rescale degrees to radians
-toRadians = math.pi/180.0 if options.degrees else 1.0
-
-# only use normalized pole
-pole = np.array(options.pole)
-pole /= np.linalg.norm(pole)
-
-# ----- Loop over input files ----- #
-for name in filenames:
-    try:
-        h5f = damask.H5Table(name, new_file=False)
-    except:
-        continue
-    damask.util.report(scriptName, name)
-
-    # extract data from HDF5 file
-    if inputtype == 'eulers':
-        orieData = h5f.get_data(label)
-    elif inputtype == 'matrix':
-        orieData = h5f.get_data(label)
-        orieData = orieData.reshape(orieData.shape[0], 3, 3)
-    elif inputtype == 'frame':
-        vctr_a = h5f.get_data(label[0])
-        vctr_b = h5f.get_data(label[1])
-        vctr_c = h5f.get_data(label[2])
-        frame = np.column_stack((vctr_a, vctr_b, vctr_c))
-        orieData = frame.reshape(frame.shape[0], 3, 3)
-    elif inputtype == 'quaternion':
-        orieData = h5f.get_data(label)
-
-    # calculate the IPF color
-    rgbArrays = np.zeros((orieData.shape[0], 3))
-    for ci in range(rgbArrays.shape[0]):
-        if inputtype == 'eulers':
-            o = damask.Orientation(Eulers=np.array(orieData[ci, :])*toRadians,
-                                   symmetry=options.symmetry).reduced()
-        elif inputtype == 'matrix':
-            o = damask.Orientation(matrix=orieData[ci, :, :].transpose(),
-                                   symmetry=options.symmetry).reduced()
-        elif inputtype == 'frame':
-            o = damask.Orientation(matrix=orieData[ci, :, :],
-                                   symmetry=options.symmetry).reduced()
-        elif inputtype == 'quaternion':
-            o = damask.Orientation(quaternion=orieData[ci, :],
-                                   symmetry=options.symmetry).reduced()
-        rgbArrays[ci, :] = o.IPFcolor(pole)
-
-    # compose labels/headers for IPF color (RGB)
-    labelIPF = 'IPF_{:g}{:g}{:g}_{sym}'.format(*options.pole,
-                                               sym=options.symmetry.lower())
-
-    # compose cmd history (go with dataset)
-    cmd_log = scriptID + '\t' + ' '.join(sys.argv[1:])
-
-    # write data to HDF5 file
-    h5f.add_data(labelIPF, rgbArrays, cmd_log=cmd_log)

processing/post/h5_addMises.py

@@ -1,85 +0,0 @@
-#!/usr/bin/env python2.7
-# -*- coding: UTF-8 no BOM -*-
-
-import os
-import sys
-import math
-import damask
-import numpy as np
-from optparse import OptionParser
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID = ' '.join([scriptName, damask.version])
-
-
-# ----- Helper functions ----- #
-def calcMises(what, tensor):
-    """Calculate von Mises equivalent"""
-    dev = tensor - np.trace(tensor)/3.0*np.eye(3)
-    symdev = 0.5*(dev+dev.T)
-    return math.sqrt(np.sum(symdev*symdev.T) *
-                     {
-                     'stress': 3.0/2.0,
-                     'strain': 2.0/3.0,
-                     }[what.lower()])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-desp = "Add von Mises equivalent values for symmetric part of requested"
-parser = OptionParser(option_class=damask.extendableOption,
-                      usage='%prog options [file[s]]',
-                      description=desp,
-                      version=scriptID)
-parser.add_option('-e', '--strain',
-                  dest='strain',
-                  metavar='string',
-                  help='name of dataset containing strain tensors')
-parser.add_option('-s', '--stress',
-                  dest='stress',
-                  metavar='string',
-                  help='name of dataset containing stress tensors')
-
-parser.set_defaults(strain=None, stress=None)
-
-(options, filenames) = parser.parse_args()
-
-# ----- Loop over input files ----- #
-for name in filenames:
-    try:
-        h5f = damask.H5Table(name, new_file=False)
-    except:
-        continue
-    damask.util.report(scriptName, name)
-
-    # TODO:
-    #   Could use some refactoring here
-    if options.stress is not None:
-        # extract stress tensor from HDF5
-        tnsr = h5f.get_data(options.stress)
-
-        # calculate von Mises equivalent row by row
-        vmStress = np.zeros(tnsr.shape[0])
-        for ri in range(tnsr.shape[0]):
-            stressTnsr = tnsr[ri, :].reshape(3, 3)
-            vmStress[ri] = calcMises('stress', stressTnsr)
-
-        # compose label
-        label = "Mises{}".format(options.stress)
-
-        # prepare log info
-        cmd_log = scriptID + '\t' + ' '.join(sys.argv[1:])
-
-        # write data to HDF5 file
-        h5f.add_data(label, vmStress, cmd_log=cmd_log)
-
-    if options.strain is not None:
-        tnsr = h5f.get_data(options.strain)
-        vmStrain = np.zeros(tnsr.shape[0])
-        for ri in range(tnsr.shape[0]):
-            strainTnsr = tnsr[ri, :].reshape(3, 3)
-            vmStrain[ri] = calcMises('strain', strainTnsr)
-        label = "Mises{}".format(options.strain)
-        cmd_log = scriptID + '\t' + ' '.join(sys.argv[1:])
-        h5f.add_data(label, vmStrain, cmd_log=cmd_log)

processing/post/h5_addStrainTensors.py

@@ -1,156 +0,0 @@
-#!/usr/bin/env python2.7
-# -*- coding: UTF-8 no BOM -*-
-
-import os
-import sys
-import damask
-import numpy as np
-from optparse import OptionParser
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID = ' '.join([scriptName, damask.version])
-
-
-# ----- Helper functions ----- #
-def operator(stretch, strain, eigenvalues):
-    # Albrecht Bertram: Elasticity and Plasticity of Large Deformations
-    # An Introduction (3rd Edition, 2012), p. 102
-    return {'V#ln': np.log(eigenvalues),
-            'U#ln': np.log(eigenvalues),
-            'V#Biot': (np.ones(3, 'd') - 1.0/eigenvalues),
-            'U#Biot': (eigenvalues - np.ones(3, 'd')),
-            'V#Green': (np.ones(3, 'd') - 1.0/eigenvalues/eigenvalues)*0.5,
-            'U#Green': (eigenvalues*eigenvalues - np.ones(3, 'd'))*0.5,
-            }[stretch+'#'+strain]
-
-
-def calcEPS(defgrads, stretchType, strainType):
-    """Calculate specific type of strain tensor"""
-    eps = np.zeros(defgrads.shape)  # initialize container
-
-    # TODO:
-    # this loop can use some performance boost
-    # (multi-threading?)
-    for ri in range(defgrads.shape[0]):
-        f = defgrads[ri, :].reshape(3, 3)
-        U, S, Vh = np.linalg.svd(f)
-        R = np.dot(U, Vh)  # rotation of polar decomposition
-        if stretchType == 'U':
-            stretch = np.dot(np.linalg.inv(R), f)  # F = RU
-        elif stretchType == 'V':
-            stretch = np.dot(f, np.linalg.inv(R))  # F = VR
-
-        # kill nasty noisy data
-        stretch = np.where(abs(stretch) < 1e-12, 0, stretch)
-
-        (D, V) = np.linalg.eig(stretch)
-        # flip principal component with negative Eigen values
-        neg = np.where(D < 0.0)
-        D[neg] *= -1.
-        V[:, neg] *= -1.
-
-        # check each vector for orthogonality
-        # --> brutal force enforcing orthogonal base
-        #     and re-normalize
-        for i, eigval in enumerate(D):
-            if np.dot(V[:, i], V[:, (i+1) % 3]) != 0.0:
-                V[:, (i+1) % 3] = np.cross(V[:, (i+2) % 3], V[:, i])
-                V[:, (i+1) % 3] /= np.sqrt(np.dot(V[:, (i+1) % 3],
-                                                  V[:, (i+1) % 3].conj()))
-
-        # calculate requested version of strain tensor
-        d = operator(stretchType, strainType, D)
-        eps[ri] = (np.dot(V, np.dot(np.diag(d), V.T)).real).reshape(9)
-
-    return eps
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-desp = "Add column(s) containing given strains based on given stretches"
-parser = OptionParser(option_class=damask.extendableOption,
-                      usage='%prog options [file[s]]',
-                      description=desp,
-                      version=scriptID)
-msg = 'material strains based on right Cauchy-Green deformation, i.e., C and U'
-parser.add_option('-u', '--right',
-                  dest='right',
-                  action='store_true',
-                  help=msg)
-msg = 'spatial strains based on left Cauchy--Green deformation, i.e., B and V'
-parser.add_option('-v', '--left',
-                  dest='left',
-                  action='store_true',
-                  help=msg)
-parser.add_option('-0', '--logarithmic',
-                  dest='logarithmic',
-                  action='store_true',
-                  help='calculate logarithmic strain tensor')
-parser.add_option('-1', '--biot',
-                  dest='biot',
-                  action='store_true',
-                  help='calculate biot strain tensor')
-parser.add_option('-2', '--green',
-                  dest='green',
-                  action='store_true',
-                  help='calculate green strain tensor')
-# NOTE:
-#   It might be easier to just calculate one type of deformation gradient
-#   at a time.
-msg = 'heading(s) of columns containing deformation tensor values'
-parser.add_option('-f', '--defgrad',
-                  dest='defgrad',
-                  action='extend',
-                  metavar='<string LIST>',
-                  help=msg)
-
-parser.set_defaults(right=False, left=False,
-                    logarithmic=False, biot=False, green=False,
-                    defgrad='f')
-
-(options, filenames) = parser.parse_args()
-
-stretches = []
-strains = []
-
-if options.right:
-    stretches.append('U')
-if options.left:
-    stretches.append('V')
-
-if options.logarithmic:
-    strains.append('ln')
-if options.biot:
-    strains.append('Biot')
-if options.green:
-    strains.append('Green')
-
-if options.defgrad is None:
-    parser.error('no data column specified.')
-
-# ----- Loop over input files ----- #
-for name in filenames:
-    try:
-        h5f = damask.H5Table(name, new_file=False)
-    except:
-        continue
-    damask.util.report(scriptName, name)
-
-    # extract defgrads from HDF5 storage
-    F = h5f.get_data(options.defgrad)
-
-    # allow calculate multiple types of strain within the
-    # same cmd call
-    for stretchType in stretches:
-        for strainType in strains:
-            # calculate strain tensor for this type
-            eps = calcEPS(F, stretchType, strainType)
-
-            # compose labels/headers for this strain tensor
-            labelsStrain = strainType + stretchType
-
-            # prepare log info
-            cmd_log = scriptID + '\t' + ' '.join(sys.argv[1:])
-
-            # write data to HDF5 file
-            h5f.add_data(labelsStrain, eps, cmd_log=cmd_log)

processing/post/h5_addXdmfWapper.py

@@ -1,130 +0,0 @@
-#!/usr/bin/env python2.7
-# -*- coding: UTF-8 no BOM -*-
-
-# ------------------------------------------------------------------- #
-# NOTE:                                                               #
-#   1. Current Xdmf rendering in Paraview has some memory issue where #
-#      large number of polyvertices will cause segmentation fault. By #
-#      default, paraview output a cell based xdmf description, which  #
-#      is working for small and medium mesh (<10,000) points. Hence a #
-#      rectangular mesh is used as the de facto Geometry description  #
-#      here.                                                          #
-#   2. Due to the unstable state Xdmf, it is safer to use port data   #
-#      to VTR rather than using xdmf as interpretive layer for data   #
-#      visualization.                                                 #
-# ------------------------------------------------------------------- #
-
-
-import os
-import damask
-import h5py
-import xml.etree.cElementTree as ET
-from optparse import OptionParser
-from xml.dom import minidom
-from damask.h5table import lables_to_path
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# ----- HELPER FUNCTIONS -----#
-def addTopLvlCmt(xmlstr, topLevelCmt):
-    """Add top level comment to string from ET"""
-    # a quick hack to add the top level comment to XML file
-    # --> somehow Elementtree does not provide this functionality
-    # --> by default
-    strList = xmlstr.split("\n")
-    strList[0] += "\n"+topLevelCmt
-    return "\n".join(strList)
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-msg = 'Generate Xdmf wrapper for HDF5 file.'
-parser = OptionParser(option_class=damask.extendableOption,
-                      usage='%prog options [file[s]]',
-                      description = msg,
-                      version = scriptID)
-
-(options, filenames) = parser.parse_args()
-
-h5f = filenames[0]
-h5f_base = h5f.split("/")[-1]
-
-# ----- parse HDF5 file ----- #
-h5f_dataDim = {}
-h5f_dataPath = {}
-h5f_dataType = {}
-with h5py.File(h5f, 'a') as f:
-    labels = f.keys()
-    labels += f['/Constitutive'].keys()
-    labels += f['/Crystallite'].keys()
-    labels += ['Vx', 'Vy', "Vz"]
-    # remove group names as they do not contain real data
-    # TODO: use h5py/H5table API to detect dataset name to
-    #       avoid necessary name space pruning.
-    labels.remove('Constitutive')
-    labels.remove('Crystallite')
-    labels.remove('Geometry')
-    # loop through remaining labels
-    for label in labels:
-        dataType, h5Path = lables_to_path(label)
-        h5f_dataType[label] = dataType
-        h5f_dataDim[label] = " ".join(map(str,f[h5Path].shape))
-        h5f_dataPath[label] = h5Path
-
-# ----- constructing xdmf elements ----- #
-root = ET.Element("Xdmf", version='3.3')
-root.set('xmlns:xi', "http://www.w3.org/2001/XInclude")
-root.append(ET.Comment('Generated Xdmf wapper for DAMASH H5 output'))
-
-# usually there should only be ONE domain
-domain = ET.SubElement(root, 'Domain',
-                       Name=h5f_base.split(".")[0])
-
-# use global topology through reference
-grid = ET.SubElement(domain, 'Grid', GridType="Uniform")
-# geometry section
-geometry = ET.SubElement(grid, 'Geometry', GeometryType="VXVYVZ")
-for vector in ["Vz", "Vy", "Vx"]:
-    dataitem = ET.SubElement(geometry, "DataItem",
-                             DataType="Float",
-                             Dimensions=h5f_dataDim[vector],
-                             Name=vector,
-                             Format="HDF")
-    dataitem.text = h5f_base.split("/")[-1] + ":{}".format(h5f_dataPath[vector])
-# topology section
-# TODO: support for other format based on given option
-meshDim = [h5f_dataDim["Vz"], h5f_dataDim["Vy"], h5f_dataDim["Vx"]]
-topology = ET.SubElement(grid, 'Topology',
-                         TopologyType="3DRectMesh",
-                         Dimensions=" ".join(map(str, meshDim)))
-
-# attributes section
-# Question: how to properly handle data mapping for multiphase situations
-labelsProcessed = ['Vx', 'Vy', 'Vz']
-# walk through each attributes
-for label in labels:
-    if label in labelsProcessed: continue
-    print("adding {}...".format(label))
-    attr = ET.SubElement(grid, 'Attribute',
-                         Name=label,
-                         Type="None",
-                         Center="Cell")
-    dataitem = ET.SubElement(attr, 'DataItem',
-                             Name=label,
-                             Format='HDF',
-                             Dimensions=h5f_dataDim[label])
-    dataitem.text = h5f_base + ":" + h5f_dataPath[label]
-    # update progress list
-    labelsProcessed.append(label)
-
-
-# pretty print the xdmf(xml) file content
-xmlstr = minidom.parseString(ET.tostring(root)).toprettyxml(indent="\t")
-xmlstr = addTopLvlCmt(xmlstr, '<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>')
-# write str to file through native python API
-with open(h5f.replace(".h5", ".xmf"), 'w') as f:
-    f.write(xmlstr)

processing/post/h5_vtkAddRectilinearGridData.py

@@ -1,191 +0,0 @@
-#!/usr/bin/env python2.7
-# -*- coding: UTF-8 no BOM -*-
-
-import os
-import vtk
-import damask
-from vtk.util import numpy_support
-from optparse import OptionParser
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID = ' '.join([scriptName, damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-msg = "Add scalars, vectors, and/or an RGB tuple from"
-msg += "an HDF5 to existing VTK rectilinear grid (.vtr/.vtk)."
-parser = OptionParser(option_class=damask.extendableOption,
-                      usage='%prog options [file[s]]',
-                      description=msg,
-                      version=scriptID)
-parser.add_option('--vtk',
-                  dest='vtk',
-                  type='string', metavar='string',
-                  help='VTK file name')
-parser.add_option('--inplace',
-                  dest='inplace',
-                  action='store_true',
-                  help='modify VTK file in-place')
-parser.add_option('-r', '--render',
-                  dest='render',
-                  action='store_true',
-                  help='open output in VTK render window')
-parser.add_option('-d', '--data',
-                  dest='data',
-                  action='extend', metavar='<string LIST>',
-                  help='scalar/vector value(s) label(s)')
-parser.add_option('-t', '--tensor',
-                  dest='tensor',
-                  action='extend', metavar='<string LIST>',
-                  help='tensor (3x3) value label(s)')
-parser.add_option('-c', '--color',
-                  dest='color',
-                  action='extend', metavar='<string LIST>',
-                  help='RGB color tuple label')
-parser.add_option('-m',
-                  '--mode',
-                  dest='mode',
-                  metavar='string',
-                  type='choice', choices=['cell', 'point'],
-                  help='cell-centered or point-centered coordinates')
-
-parser.set_defaults(data=[],
-                    tensor=[],
-                    color=[],
-                    mode='cell',
-                    inplace=False,
-                    render=False)
-
-(options, filenames) = parser.parse_args()
-
-# ----- Legacy VTK format support ----- #
-if os.path.splitext(options.vtk)[1] == '.vtr':
-    reader = vtk.vtkXMLRectilinearGridReader()
-    reader.SetFileName(options.vtk)
-    reader.Update()
-    rGrid = reader.GetOutput()
-elif os.path.splitext(options.vtk)[1] == '.vtk':
-    reader = vtk.vtkGenericDataObjectReader()
-    reader.SetFileName(options.vtk)
-    reader.Update()
-    rGrid = reader.GetRectilinearGridOutput()
-else:
-    parser.error('Unsupported VTK file type extension.')
-
-Npoints = rGrid.GetNumberOfPoints()
-Ncells = rGrid.GetNumberOfCells()
-
-# ----- Summary output (Sanity Check) ----- #
-msg = '{}: {} points and {} cells...'.format(options.vtk,
-                                             Npoints,
-                                             Ncells)
-damask.util.croak(msg)
-
-# ----- Read HDF5 file ----- #
-# NOTE:
-# --> It is possible in the future we are trying to add data
-#     from different increment into the same VTK file, but
-#     this feature is not supported for the moment.
-# --> Let it fail, if the HDF5 is invalid, python interpretor
-# --> should be able to catch this error.
-h5f = damask.H5Table(filenames[0], new_file=False)
-
-# ----- Process data ----- #
-featureToAdd = {'data': options.data,
-                'tensor': options.tensor,
-                'color': options.color}
-VTKarray = {}  # store all vtkData in dict, then ship them to file
-for dataType in featureToAdd.keys():
-    featureNames = featureToAdd[dataType]
-    for featureName in featureNames:
-        VTKtype = vtk.VTK_DOUBLE
-        VTKdata = h5f.get_data(featureName)
-        if dataType == 'color':
-            VTKtype = vtk.VTK_UNSIGNED_CHAR
-            VTKdata = (VTKdata*255).astype(int)
-        elif dataType == 'tensor':
-            # Force symmetries tensor type data
-            VTKdata[:, 1] = VTKdata[:, 3] = 0.5*(VTKdata[:, 1]+VTKdata[:, 3])
-            VTKdata[:, 2] = VTKdata[:, 6] = 0.5*(VTKdata[:, 2]+VTKdata[:, 6])
-            VTKdata[:, 5] = VTKdata[:, 7] = 0.5*(VTKdata[:, 5]+VTKdata[:, 7])
-        # use vtk build-in numpy support to add data (much faster)
-        # NOTE:
-        # --> deep copy is necessary here, otherwise memory leak could occur
-        VTKarray[featureName] = numpy_support.numpy_to_vtk(num_array=VTKdata,
-                                                           deep=True,
-                                                           array_type=VTKtype)
-        VTKarray[featureName].SetName(featureName)
-
-# ----- ship data to vtkGrid ----- #
-mode = options.mode
-damask.util.croak('{} mode...'.format(mode))
-
-# NOTE:
-# --> For unknown reason, Paraview only recognize one
-#     tensor attributes per cell, thus it would be safe
-#     to only add one attributes as tensor.
-for dataType in featureToAdd.keys():
-    featureNames = featureToAdd[dataType]
-    for featureName in featureNames:
-        if dataType == 'color':
-            if mode == 'cell':
-                rGrid.GetCellData().SetScalars(VTKarray[featureName])
-            elif mode == 'point':
-                rGrid.GetPointData().SetScalars(VTKarray[featureName])
-        elif dataType == 'tensor':
-            if mode == 'cell':
-                rGrid.GetCellData().SetTensors(VTKarray[featureName])
-            elif mode == 'point':
-                rGrid.GetPointData().SetTensors(VTKarray[featureName])
-        else:
-            if mode == 'cell':
-                rGrid.GetCellData().AddArray(VTKarray[featureName])
-            elif mode == 'point':
-                rGrid.GetPointData().AddArray(VTKarray[featureName])
-
-rGrid.Modified()
-if vtk.VTK_MAJOR_VERSION <= 5:
-    rGrid.Update()
-
-# ----- write Grid to VTK file ----- #
-writer = vtk.vtkXMLRectilinearGridWriter()
-writer.SetDataModeToBinary()
-writer.SetCompressorTypeToZLib()
-vtkFileN = os.path.splitext(options.vtk)[0]
-vtkExtsn = '.vtr' if options.inplace else '_added.vtr'
-writer.SetFileName(vtkFileN+vtkExtsn)
-if vtk.VTK_MAJOR_VERSION <= 5:
-    writer.SetInput(rGrid)
-else:
-    writer.SetInputData(rGrid)
-writer.Write()
-
-# ----- render results from script ----- #
-if options.render:
-    mapper = vtk.vtkDataSetMapper()
-    mapper.SetInputData(rGrid)
-    actor = vtk.vtkActor()
-    actor.SetMapper(mapper)
-
-    # Create the graphics structure. The renderer renders into the
-    # render window. The render window interactor captures mouse events
-    # and will perform appropriate camera or actor manipulation
-    # depending on the nature of the events.
-
-    ren = vtk.vtkRenderer()
-
-    renWin = vtk.vtkRenderWindow()
-    renWin.AddRenderer(ren)
-
-    ren.AddActor(actor)
-    ren.SetBackground(1, 1, 1)
-    renWin.SetSize(200, 200)
-
-    iren = vtk.vtkRenderWindowInteractor()
-    iren.SetRenderWindow(renWin)
-
-    iren.Initialize()
-    renWin.Render()
-    iren.Start()

processing/post/h5_vtkRectilinearGrid.py

@@ -1,135 +0,0 @@
-#!/usr/bin/env python2.7
-# -*- coding: UTF-8 no BOM -*-
-
-# ------------------------------------------------------------------ #
-# NOTE:                                                              #
-#   1. It might be a good idea to separate IO and calculation.       #
-#   2. Some of the calculation could be useful in other situations,  #
-#      why not build a math_util, or math_sup module that contains   #
-#      all the useful functions.                                     #
-# ------------------------------------------------------------------ #
-
-import os
-import vtk
-import numpy as np
-import damask
-from optparse import OptionParser
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID = ' '.join([scriptName, damask.version])
-
-
-# ----- HELPER FUNCTION ----- #
-def getMeshFromXYZ(xyzArray, mode):
-    """Calc Vx,Vy,Vz vectors for vtk rectangular mesh"""
-    # NOTE:
-    # --> np.unique will automatically sort the list
-    # --> although not exactly n(1), but since mesh dimension should
-    #     small anyway, so this is still light weight.
-    dim = xyzArray.shape[1]  # 2D:2, 3D:3
-    coords = [np.unique(xyzArray[:, i]) for i in range(dim)]
-
-    if mode == 'cell':
-        # since x, y, z might now have the same number of elements,
-        # we have to deal with them individually
-        for ri in range(dim):
-            vctr_pt = coords[ri]
-            vctr_cell = np.empty(len(vctr_pt)+1)
-            # calculate first and last end point
-            vctr_cell[0] = vctr_pt[0] - 0.5*abs(vctr_pt[1] - vctr_pt[0])
-            vctr_cell[-1] = vctr_pt[-1] + 0.5*abs(vctr_pt[-2] - vctr_pt[-1])
-            for cj in range(1, len(vctr_cell)-1):
-                vctr_cell[cj] = 0.5*(vctr_pt[cj-1] + vctr_pt[cj])
-            # update the coords
-            coords[ri] = vctr_cell
-
-    if dim < 3:
-        coords.append([0])  # expand to a 3D with 0 for z
-
-    # auxiliary description
-    grid = np.array(map(len, coords), 'i')
-    N = grid.prod() if mode == 'point' else (grid-1).prod()
-    return coords, grid, N
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-msg = "Create regular voxel grid from points in an ASCIItable."
-parser = OptionParser(option_class=damask.extendableOption,
-                      usage='%prog options [file[s]]',
-                      description=msg,
-                      version=scriptID)
-
-parser.add_option('-m',
-                  '--mode',
-                  dest='mode',
-                  metavar='string',
-                  type='choice', choices=['cell', 'point'],
-                  help='cell-centered or point-centered coordinates')
-parser.add_option('-p',
-                  '--pos', '--position',
-                  dest='pos',
-                  type='string', metavar='string',
-                  help='label of coordinates [%default]')
-
-parser.set_defaults(mode='cell',
-                    pos='pos')
-
-(options, filenames) = parser.parse_args()
-
-# ----- loop over input files ----- #
-for name in filenames:
-    try:
-        h5f = damask.H5Table(name, new_file=False)
-    except:
-        continue
-    damask.util.report(scriptName, name)
-
-    # ----- read xyzArray from HDF5 file ----- #
-    xyzArray = h5f.get_data(options.pos)
-
-    # ----- figure out size and grid ----- #
-    coords, grid, N = getMeshFromXYZ(xyzArray, options.mode)
-
-    # ----- process data ----- #
-    rGrid = vtk.vtkRectilinearGrid()
-    # WARNING: list expansion does not work here as these are
-    #          just pointers for a vtk instance. Simply put,
-    #          DON't USE
-    #            [<VTK_CONSTRUCTOR>] * <NUM_OF_ELEMENTS>
-    coordArray = [vtk.vtkDoubleArray(),
-                  vtk.vtkDoubleArray(),
-                  vtk.vtkDoubleArray()]
-
-    rGrid.SetDimensions(*grid)
-    for i, points in enumerate(coords):
-        for point in points:
-            coordArray[i].InsertNextValue(point)
-
-    rGrid.SetXCoordinates(coordArray[0])
-    rGrid.SetYCoordinates(coordArray[1])
-    rGrid.SetZCoordinates(coordArray[2])
-
-    # ----- output result ----- #
-    dirPath = os.path.split(name)[0]
-    if name:
-        writer = vtk.vtkXMLRectilinearGridWriter()
-        writer.SetCompressorTypeToZLib()
-        writer.SetDataModeToBinary()
-        # getting the name is a little bit tricky
-        vtkFileName = os.path.splitext(os.path.split(name)[1])[0]
-        vtkFileName += '_{}({})'.format(options.pos, options.mode)
-        vtkFileName += '.' + writer.GetDefaultFileExtension()
-        writer.SetFileName(os.path.join(dirPath, vtkFileName))
-    else:
-        writer = vtk.vtkDataSetWriter()
-        writer.SetHeader('# powered by '+scriptID)
-        writer.WriteToOutputStringOn()
-
-    if vtk.VTK_MAJOR_VERSION <= 5:
-        writer.SetInput(rGrid)
-    else:
-        writer.SetInputData(rGrid)
-
-    writer.Write()

processing/post/marc_to_vtk.py

@@ -0,0 +1,199 @@
+#!/usr/bin/env python2.7
+# -*- coding: UTF-8 no BOM -*-
+
+import os,sys,re
+import argparse
+import damask
+import vtk, numpy as np
+
+scriptName = os.path.splitext(os.path.basename(__file__))[0]
+scriptID   = ' '.join([scriptName, damask.version])
+
+parser = argparse.ArgumentParser(description='Convert from Marc input file format (.dat) to VTK format (.vtu)', version = scriptID)
+parser.add_argument('filename', type=str, help='file to convert')
+parser.add_argument('-t', '--table',  type=str, help='ASCIItable file containing nodal data to subdivide and interpolate')
+
+args = parser.parse_args()
+
+with open(args.filename, 'r') as marcfile:
+  marctext = marcfile.read();
+  
+# Load table (if any)
+if args.table is not None:
+  try:
+    table = damask.ASCIItable(
+      name=args.table,
+      outname='subdivided_{}'.format(args.table),
+      buffered=True
+    )
+    
+    table.head_read()
+    table.data_readArray()
+    
+    # Python list is faster for appending
+    nodal_data = list(table.data)
+  except: args.table = None
+
+# Extract connectivity chunk from file...
+connectivity_text = re.findall(r'connectivity[\n\r]+(.*?)[\n\r]+[a-zA-Z]', marctext, flags=(re.MULTILINE | re.DOTALL))[0]
+connectivity_lines = re.split(r'[\n\r]+', connectivity_text, flags=(re.MULTILINE | re.DOTALL))
+connectivity_header = connectivity_lines[0]
+connectivity_lines = connectivity_lines[1:]
+
+# Construct element map
+elements = dict(map(lambda line: 
+  (
+    int(line[0:10]), # index
+    { 
+      'type':  int(line[10:20]), 
+      'verts': list(map(int, re.split(r' +', line[20:].strip())))
+    }
+  ), connectivity_lines))
+
+# Extract coordinate chunk from file
+coordinates_text = re.findall(r'coordinates[\n\r]+(.*?)[\n\r]+[a-zA-Z]', marctext, flags=(re.MULTILINE | re.DOTALL))[0]
+coordinates_lines = re.split(r'[\n\r]+', coordinates_text, flags=(re.MULTILINE | re.DOTALL))
+coordinates_header = coordinates_lines[0]
+coordinates_lines = coordinates_lines[1:]
+
+# marc input file does not use "e" in scientific notation, this adds it and converts
+fl_format = lambda string: float(re.sub(r'(\d)([\+\-])', r'\1e\2', string))
+# Construct coordinate map
+coordinates = dict(map(lambda line:
+  (
+    int(line[0:10]),
+    np.array([
+      fl_format(line[10:30]), 
+      fl_format(line[30:50]), 
+      fl_format(line[50:70])
+    ])
+  ), coordinates_lines))
+
+# Subdivide volumes
+grid = vtk.vtkUnstructuredGrid()
+vertex_count = len(coordinates)
+edge_to_vert = dict() # when edges are subdivided, a new vertex in the middle is produced and placed in here
+ordered_pair = lambda a, b: (a, b) if a < b else (b, a) # edges are bidirectional
+
+def subdivide_edge(vert1, vert2):
+  edge = ordered_pair(vert1, vert2)
+  
+  if edge in edge_to_vert:
+    return edge_to_vert[edge]
+  
+  # Vertex does not exist, create it
+  newvert = len(coordinates) + 1
+  coordinates[newvert] = 0.5 * (coordinates[vert1] + coordinates[vert2]) # Average
+  edge_to_vert[edge] = newvert;
+  
+  # Interpolate nodal data
+  if args.table is not None:
+    nodal_data.append(0.5 * (nodal_data[vert1 - 1] + nodal_data[vert2 - 1]))
+  return newvert;
+
+for el_id in range(1, len(elements) + 1): # Marc starts counting at 1
+  el = elements[el_id]
+  if el['type'] == 7:
+    # Hexahedron, subdivided
+    
+    # There may be a better way to iterate over these, but this is consistent
+    # with the ordering scheme provided at https://damask.mpie.de/pub/Documentation/ElementType
+    
+    subverts = np.zeros((3,3,3), dtype=int)
+    # Get corners
+    subverts[0, 0, 0] = el['verts'][0]
+    subverts[2, 0, 0] = el['verts'][1]
+    subverts[2, 2, 0] = el['verts'][2]
+    subverts[0, 2, 0] = el['verts'][3]
+    subverts[0, 0, 2] = el['verts'][4]
+    subverts[2, 0, 2] = el['verts'][5]
+    subverts[2, 2, 2] = el['verts'][6]
+    subverts[0, 2, 2] = el['verts'][7]
+    
+    # lower edges
+    subverts[1, 0, 0] = subdivide_edge(subverts[0, 0, 0], subverts[2, 0, 0])
+    subverts[2, 1, 0] = subdivide_edge(subverts[2, 0, 0], subverts[2, 2, 0])
+    subverts[1, 2, 0] = subdivide_edge(subverts[2, 2, 0], subverts[0, 2, 0])
+    subverts[0, 1, 0] = subdivide_edge(subverts[0, 2, 0], subverts[0, 0, 0])
+    
+    # middle edges
+    subverts[0, 0, 1] = subdivide_edge(subverts[0, 0, 0], subverts[0, 0, 2])
+    subverts[2, 0, 1] = subdivide_edge(subverts[2, 0, 0], subverts[2, 0, 2])
+    subverts[2, 2, 1] = subdivide_edge(subverts[2, 2, 0], subverts[2, 2, 2])
+    subverts[0, 2, 1] = subdivide_edge(subverts[0, 2, 0], subverts[0, 2, 2])
+    
+    # top edges
+    subverts[1, 0, 2] = subdivide_edge(subverts[0, 0, 2], subverts[2, 0, 2])
+    subverts[2, 1, 2] = subdivide_edge(subverts[2, 0, 2], subverts[2, 2, 2])
+    subverts[1, 2, 2] = subdivide_edge(subverts[2, 2, 2], subverts[0, 2, 2])
+    subverts[0, 1, 2] = subdivide_edge(subverts[0, 2, 2], subverts[0, 0, 2])
+    
+    # then faces...  The edge_to_vert addition is due to there being two ways
+    # to calculate a face vertex, depending on which opposite vertices are used to subdivide.
+    # This way, we avoid creating duplicate vertices.
+    subverts[1, 1, 0] = subdivide_edge(subverts[1, 0, 0], subverts[1, 2, 0])
+    edge_to_vert[ordered_pair(subverts[0, 1, 0], subverts[2, 1, 0])] = subverts[1, 1, 0]
+    
+    subverts[1, 0, 1] = subdivide_edge(subverts[1, 0, 0], subverts[1, 0, 2])
+    edge_to_vert[ordered_pair(subverts[0, 0, 1], subverts[2, 0, 1])] = subverts[1, 0, 1]
+    
+    subverts[2, 1, 1] = subdivide_edge(subverts[2, 1, 0], subverts[2, 1, 2])
+    edge_to_vert[ordered_pair(subverts[2, 0, 1], subverts[2, 2, 1])] = subverts[2, 1, 1]
+    
+    subverts[1, 2, 1] = subdivide_edge(subverts[1, 2, 0], subverts[1, 2, 2])
+    edge_to_vert[ordered_pair(subverts[0, 2, 1], subverts[2, 2, 1])] = subverts[1, 2, 1]
+    
+    subverts[0, 1, 1] = subdivide_edge(subverts[0, 1, 0], subverts[0, 1, 2])
+    edge_to_vert[ordered_pair(subverts[0, 0, 1], subverts[0, 2, 1])] = subverts[0, 1, 1]
+    
+    subverts[1, 1, 2] = subdivide_edge(subverts[1, 0, 2], subverts[1, 2, 2])
+    edge_to_vert[ordered_pair(subverts[0, 1, 2], subverts[2, 1, 2])] = subverts[1, 1, 2]
+    
+    # and finally the center.  There are three ways to calculate, but elements should
+    # not intersect, so the edge_to_vert part isn't needed here.
+    subverts[1, 1, 1] = subdivide_edge(subverts[1, 1, 0], subverts[1, 1, 2])
+    
+    
+    # Now make the hexahedron subelements
+    # order in which vtk expects vertices for a hexahedron
+    order = np.array([(0,0,0),(1,0,0),(1,1,0),(0,1,0),(0,0,1),(1,0,1),(1,1,1),(0,1,1)])
+    for z in range(2):
+      for y in range(2):
+        for x in range(2):
+          hex_ = vtk.vtkHexahedron()
+          for vert_id in range(8):
+            coord = order[vert_id] + (x, y, z)
+            # minus one, since vtk starts at zero but marc starts at one
+            hex_.GetPointIds().SetId(vert_id, subverts[coord[0], coord[1], coord[2]] - 1) 
+          grid.InsertNextCell(hex_.GetCellType(), hex_.GetPointIds())
+    
+  else:
+    damask.util.croak('Unsupported Marc element type: {} (skipping)'.format(el['type']))
+
+# Load all points
+points = vtk.vtkPoints()
+for point in range(1, len(coordinates) + 1): # marc indices start at 1
+  points.InsertNextPoint(coordinates[point].tolist())
+
+grid.SetPoints(points)
+
+# grid now contains the elements from the given marc file
+writer = vtk.vtkXMLUnstructuredGridWriter()
+writer.SetFileName(re.sub(r'\..+', ".vtu", args.filename))  # *.vtk extension does not work in paraview
+
+if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(grid)
+else:                          writer.SetInputData(grid)
+writer.Write()
+
+if args.table is not None:  
+  table.info_append([
+    scriptID + ' ' + ' '.join(sys.argv[1:]),
+  ])
+  table.head_write()
+  table.output_flush()
+  
+  table.data = np.array(nodal_data)
+  
+  table.data_writeArray()
+
+table.close()

processing/post/vtk_addGridData.py

@@ -0,0 +1,206 @@
+#!/usr/bin/env python2.7
+# -*- coding: UTF-8 no BOM -*-
+
+import os,vtk
+import damask
+from vtk.util import numpy_support
+from collections import defaultdict
+from optparse import OptionParser
+
+scriptName = os.path.splitext(os.path.basename(__file__))[0]
+scriptID   = ' '.join([scriptName,damask.version])
+
+# --------------------------------------------------------------------
+#                                MAIN
+# --------------------------------------------------------------------
+
+msg = "Add scalars, vectors, and/or an RGB tuple from"
+msg += "an ASCIItable to existing VTK grid (.vtr/.vtk/.vtu)."
+parser = OptionParser(option_class=damask.extendableOption,
+                      usage='%prog options [file[s]]',
+                      description = msg,
+                      version = scriptID)
+
+parser.add_option(      '--vtk',
+                  dest = 'vtk',
+                  type = 'string', metavar = 'string',
+                  help = 'VTK file name')
+parser.add_option(      '--inplace',
+                  dest = 'inplace',
+                  action = 'store_true',
+                  help = 'modify VTK file in-place')
+parser.add_option('-r', '--render',
+                  dest = 'render',
+                  action = 'store_true',
+                  help = 'open output in VTK render window')
+parser.add_option('-d', '--data',
+                  dest = 'data',
+                  action = 'extend', metavar = '<string LIST>',
+                  help = 'scalar/vector value(s) label(s)')
+parser.add_option('-t', '--tensor',
+                  dest = 'tensor',
+                  action = 'extend', metavar = '<string LIST>',
+                  help = 'tensor (3x3) value label(s)')
+parser.add_option('-c', '--color',
+                  dest = 'color',
+                  action = 'extend', metavar = '<string LIST>',
+                  help = 'RGB color tuple label')
+
+parser.set_defaults(data = [],
+                    tensor = [],
+                    color = [],
+                    inplace = False,
+                    render = False,
+)
+
+(options, filenames) = parser.parse_args()
+
+if not options.vtk:                 parser.error('No VTK file specified.')
+if not os.path.exists(options.vtk): parser.error('VTK file does not exist.')
+
+if os.path.splitext(options.vtk)[1] == '.vtr':
+  reader = vtk.vtkXMLRectilinearGridReader()
+  reader.SetFileName(options.vtk)
+  reader.Update()
+  rGrid = reader.GetOutput()
+  writer = vtk.vtkXMLRectilinearGridWriter()
+  writer.SetFileName(os.path.splitext(options.vtk)[0]+('.vtr' if options.inplace else '_added.vtr'))
+elif os.path.splitext(options.vtk)[1] == '.vtk':
+  reader = vtk.vtkGenericDataObjectReader()
+  reader.SetFileName(options.vtk)
+  reader.Update()
+  rGrid = reader.GetRectilinearGridOutput()
+  writer = vtk.vtkXMLRectilinearGridWriter()
+  writer.SetFileName(os.path.splitext(options.vtk)[0]+('.vtr' if options.inplace else '_added.vtr'))
+elif os.path.splitext(options.vtk)[1] == '.vtu':
+  reader = vtk.vtkXMLUnstructuredGridReader()
+  reader.SetFileName(options.vtk)
+  reader.Update()
+  rGrid = reader.GetOutput()
+  writer = vtk.vtkXMLUnstructuredGridWriter()
+  writer.SetFileName(os.path.splitext(options.vtk)[0]+('.vtu' if options.inplace else '_added.vtu'))
+else:
+  parser.error('Unsupported VTK file type extension.')
+
+Npoints = rGrid.GetNumberOfPoints()
+Ncells  = rGrid.GetNumberOfCells()
+
+damask.util.croak('{}: {} points and {} cells...'.format(options.vtk,Npoints,Ncells))
+
+# --- loop over input files -------------------------------------------------------------------------
+
+if filenames == []: filenames = [None]
+
+for name in filenames:
+  try:    table = damask.ASCIItable(name = name,
+                                    buffered = False,
+                                    readonly = True)
+  except: continue
+  damask.util.report(scriptName, name)
+
+# --- interpret header ----------------------------------------------------------------------------
+
+  table.head_read()
+
+  remarks = []
+  errors  = []
+  VTKarray = {}
+  active = defaultdict(list)
+
+  for datatype,dimension,label in [['data',99,options.data],
+                                   ['tensor',9,options.tensor],
+                                   ['color' ,3,options.color],
+                                   ]:
+    for i,dim in enumerate(table.label_dimension(label)):
+      me = label[i]
+      if dim == -1:         remarks.append('{} "{}" not found...'.format(datatype,me))
+      elif dim > dimension: remarks.append('"{}" not of dimension {}...'.format(me,dimension))
+      else:
+        remarks.append('adding {} "{}"...'.format(datatype,me))
+        active[datatype].append(me)
+
+  if remarks != []: damask.util.croak(remarks)
+  if errors  != []:
+    damask.util.croak(errors)
+    table.close(dismiss = True)
+    continue
+
+# ------------------------------------------ process data ---------------------------------------
+
+  table.data_readArray([item for sublist in active.values() for item in sublist])                 # read all requested data
+
+  for datatype,labels in active.items():                                                          # loop over scalar,color
+    for me in labels:                                                                             # loop over all requested items
+      VTKtype = vtk.VTK_DOUBLE
+      VTKdata = table.data[:, table.label_indexrange(me)].copy()                                  # copy to force contiguous layout
+
+      if datatype == 'color':
+        VTKtype = vtk.VTK_UNSIGNED_CHAR
+        VTKdata = (VTKdata*255).astype(int)                                                       # translate to 0..255 UCHAR
+      elif datatype == 'tensor':
+        VTKdata[:,1] = VTKdata[:,3] = 0.5*(VTKdata[:,1]+VTKdata[:,3])
+        VTKdata[:,2] = VTKdata[:,6] = 0.5*(VTKdata[:,2]+VTKdata[:,6])
+        VTKdata[:,5] = VTKdata[:,7] = 0.5*(VTKdata[:,5]+VTKdata[:,7])
+
+      VTKarray[me] = numpy_support.numpy_to_vtk(num_array=VTKdata,deep=True,array_type=VTKtype)
+      VTKarray[me].SetName(me)
+
+  table.close()                                                                                     # close input ASCII table
+
+# ------------------------------------------ add data ---------------------------------------
+
+  if   len(table.data) == Npoints:  mode = 'point'
+  elif len(table.data) == Ncells:   mode = 'cell'
+  else:
+    damask.util.croak('Data count is incompatible with grid...')
+    continue
+
+  damask.util.croak('{} mode...'.format(mode))
+
+  for datatype,labels in active.items():                                                            # loop over scalar,color
+    if datatype == 'color':
+      if   mode == 'cell':  rGrid.GetCellData().SetScalars(VTKarray[active['color'][0]])
+      elif mode == 'point': rGrid.GetPointData().SetScalars(VTKarray[active['color'][0]])
+    for me in labels:                                                                               # loop over all requested items
+      if   mode == 'cell':  rGrid.GetCellData().AddArray(VTKarray[me])
+      elif mode == 'point': rGrid.GetPointData().AddArray(VTKarray[me])
+
+  rGrid.Modified()
+  if vtk.VTK_MAJOR_VERSION <= 5: rGrid.Update()
+
+# ------------------------------------------ output result ---------------------------------------
+
+  writer.SetDataModeToBinary()
+  writer.SetCompressorTypeToZLib()
+  if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(rGrid)
+  else:                          writer.SetInputData(rGrid)
+  writer.Write()
+
+# ------------------------------------------ render result ---------------------------------------
+
+if options.render:
+  mapper = vtk.vtkDataSetMapper()
+  mapper.SetInputData(rGrid)
+  actor = vtk.vtkActor()
+  actor.SetMapper(mapper)
+
+# Create the graphics structure. The renderer renders into the
+# render window. The render window interactor captures mouse events
+# and will perform appropriate camera or actor manipulation
+# depending on the nature of the events.
+
+  ren = vtk.vtkRenderer()
+
+  renWin = vtk.vtkRenderWindow()
+  renWin.AddRenderer(ren)
+
+  ren.AddActor(actor)
+  ren.SetBackground(1, 1, 1)
+  renWin.SetSize(200, 200)
+
+  iren = vtk.vtkRenderWindowInteractor()
+  iren.SetRenderWindow(renWin)
+
+  iren.Initialize()
+  renWin.Render()
+  iren.Start()

src/CMakeLists.txt

@@ -74,10 +74,8 @@ add_library (PLASTIC OBJECT
   "plastic_disloUCLA.f90"
   "plastic_isotropic.f90"
   "plastic_phenopowerlaw.f90"
-  "plastic_titanmod.f90"
   "plastic_nonlocal.f90"
-  "plastic_none.f90"
-  "plastic_phenoplus.f90")
+  "plastic_none.f90")
 add_dependencies(PLASTIC DAMASK_HELPERS)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:PLASTIC>)
 

src/DAMASK_spectral.f90

@@ -485,8 +485,8 @@ program DAMASK_spectral
                              min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
      call MPI_file_write(resUnit, &
                          reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)), &
-                                [(outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults]), &
-                         (outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults, &
+                                [(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)]), &
+                         (outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt), &
                          MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
      if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_write')
    enddo
@@ -683,8 +683,8 @@ program DAMASK_spectral
            outputIndex=int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, &
                       min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
            call MPI_file_write(resUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)),&
-                                         [(outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults]), &
-                               (outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults,&
+                                         [(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)]), &
+                               (outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt),&
                                MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
            if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_write')
          enddo

src/IO.f90

@@ -50,6 +50,7 @@ module IO
    IO_skipChunks, &
    IO_extractValue, &
    IO_countDataLines, &
+   IO_countNumericalDataLines, &
    IO_countContinuousIntValues, &
    IO_continuousIntValues, &
    IO_error, &
@@ -1247,6 +1248,37 @@ integer(pInt) function IO_countDataLines(fileUnit)
 end function IO_countDataLines
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief count lines containig data up to next *keyword
+!--------------------------------------------------------------------------------------------------
+integer(pInt) function IO_countNumericalDataLines(fileUnit)
+
+ implicit none
+ integer(pInt), intent(in)                :: fileUnit                                               !< file handle
+
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ character(len=65536)                     :: line, &
+                                             tmp
+
+ IO_countNumericalDataLines = 0_pInt
+ line = ''
+
+ do while (trim(line) /= IO_EOF)
+   line = IO_read(fileUnit)
+   chunkPos = IO_stringPos(line)
+   tmp = IO_lc(IO_stringValue(line,chunkPos,1_pInt))
+   if (verify(trim(tmp) ,"0123456789")/=0) then                                                            ! found keyword
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit
+   else
+     IO_countNumericalDataLines = IO_countNumericalDataLines + 1_pInt
+   endif
+ enddo
+ backspace(fileUnit)
+
+end function IO_countNumericalDataLines
+
 !--------------------------------------------------------------------------------------------------
 !> @brief count items in consecutive lines depending on lines
 !> @details Marc:      ints concatenated by "c" as last char or range of values a "to" b
@@ -1313,7 +1345,7 @@ end function IO_countContinuousIntValues
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief return integer list corrsponding to items in consecutive lines.
+!> @brief return integer list corresponding to items in consecutive lines.
 !! First integer in array is counter
 !> @details Marc:      ints concatenated by "c" as last char, range of a "to" b, or named set
 !! Abaqus:    triplet of start,stop,inc or named set
@@ -1807,12 +1839,12 @@ integer(pInt) function IO_verifyIntValue (string,validChars,myName)
  invalidWhere = verify(string,validChars)
  if (invalidWhere == 0_pInt) then
    read(UNIT=string,iostat=readStatus,FMT=*) IO_verifyIntValue                                        ! no offending chars found
-   if (readStatus /= 0_pInt) &                                                                        ! error during string to float conversion
+   if (readStatus /= 0_pInt) &                                                                        ! error during string to integer conversion
      call IO_warning(203_pInt,ext_msg=myName//'"'//string//'"')
  else
    call IO_warning(202_pInt,ext_msg=myName//'"'//string//'"')                                         ! complain about offending characters
    read(UNIT=string(1_pInt:invalidWhere-1_pInt),iostat=readStatus,FMT=*) IO_verifyIntValue            ! interpret remaining string
-   if (readStatus /= 0_pInt) &                                                                        ! error during string to float conversion
+   if (readStatus /= 0_pInt) &                                                                        ! error during string to integer conversion
      call IO_warning(203_pInt,ext_msg=myName//'"'//string(1_pInt:invalidWhere-1_pInt)//'"')
  endif
 

src/commercialFEM_fileList.f90

@@ -28,8 +28,6 @@
 #include "plastic_none.f90"
 #include "plastic_isotropic.f90"
 #include "plastic_phenopowerlaw.f90"
-#include "plastic_phenoplus.f90"
-#include "plastic_titanmod.f90"
 #include "plastic_dislotwin.f90"
 #include "plastic_disloUCLA.f90"
 #include "plastic_nonlocal.f90"

src/constitutive.f90

@@ -74,10 +74,8 @@ subroutine constitutive_init()
    PLASTICITY_none_ID, &
    PLASTICITY_isotropic_ID, &
    PLASTICITY_phenopowerlaw_ID, &
-   PLASTICITY_phenoplus_ID, &
    PLASTICITY_dislotwin_ID, &
    PLASTICITY_disloucla_ID, &
-   PLASTICITY_titanmod_ID, &
    PLASTICITY_nonlocal_ID ,&
    SOURCE_thermal_dissipation_ID, &
    SOURCE_thermal_externalheat_ID, &
@@ -97,10 +95,8 @@ subroutine constitutive_init()
    PLASTICITY_NONE_label, &
    PLASTICITY_ISOTROPIC_label, &
    PLASTICITY_PHENOPOWERLAW_label, &
-   PLASTICITY_PHENOPLUS_label, &
    PLASTICITY_DISLOTWIN_label, &
    PLASTICITY_DISLOUCLA_label, &
-   PLASTICITY_TITANMOD_label, &
    PLASTICITY_NONLOCAL_label, &
    SOURCE_thermal_dissipation_label, &
    SOURCE_thermal_externalheat_label, &
@@ -117,10 +113,8 @@ subroutine constitutive_init()
  use plastic_none
  use plastic_isotropic
  use plastic_phenopowerlaw
- use plastic_phenoplus
  use plastic_dislotwin
  use plastic_disloucla
- use plastic_titanmod
  use plastic_nonlocal
  use source_thermal_dissipation
  use source_thermal_externalheat
@@ -162,10 +156,8 @@ subroutine constitutive_init()
  if (any(phase_plasticity == PLASTICITY_NONE_ID))          call plastic_none_init
  if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID))     call plastic_isotropic_init(FILEUNIT)
  if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init(FILEUNIT)
- if (any(phase_plasticity == PLASTICITY_PHENOPLUS_ID))     call plastic_phenoplus_init(FILEUNIT)
  if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID))     call plastic_dislotwin_init(FILEUNIT)
  if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID))     call plastic_disloucla_init(FILEUNIT)
- if (any(phase_plasticity == PLASTICITY_TITANMOD_ID))      call plastic_titanmod_init(FILEUNIT)
  if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then
   call plastic_nonlocal_init(FILEUNIT)
   call plastic_nonlocal_stateInit()
@@ -194,11 +186,11 @@ subroutine constitutive_init()
  if (any(phase_kinematics == KINEMATICS_hydrogen_strain_ID))   call kinematics_hydrogen_strain_init(FILEUNIT)
  close(FILEUNIT)
 
- mainProcess: if (worldrank == 0) then
  write(6,'(/,a)')   ' <<<+-  constitutive init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
 
+ mainProcess: if (worldrank == 0) then
 !--------------------------------------------------------------------------------------------------
 ! write description file for constitutive output
    call IO_write_jobFile(FILEUNIT,'outputConstitutive')
@@ -222,11 +214,6 @@ subroutine constitutive_init()
            thisNoutput => plastic_phenopowerlaw_Noutput
            thisOutput => plastic_phenopowerlaw_output
            thisSize   => plastic_phenopowerlaw_sizePostResult
-         case (PLASTICITY_PHENOPLUS_ID) plasticityType
-           outputName = PLASTICITY_PHENOPLUS_label
-           thisNoutput => plastic_phenoplus_Noutput
-           thisOutput => plastic_phenoplus_output
-           thisSize   => plastic_phenoplus_sizePostResult
          case (PLASTICITY_DISLOTWIN_ID) plasticityType
            outputName = PLASTICITY_DISLOTWIN_label
            thisNoutput => plastic_dislotwin_Noutput
@@ -237,11 +224,6 @@ subroutine constitutive_init()
            thisNoutput => plastic_disloucla_Noutput
            thisOutput => plastic_disloucla_output
            thisSize   => plastic_disloucla_sizePostResult
-         case (PLASTICITY_TITANMOD_ID) plasticityType
-           outputName = PLASTICITY_TITANMOD_label
-           thisNoutput => plastic_titanmod_Noutput
-           thisOutput => plastic_titanmod_output
-           thisSize   => plastic_titanmod_sizePostResult
          case (PLASTICITY_NONLOCAL_ID) plasticityType
            outputName = PLASTICITY_NONLOCAL_label
            thisNoutput => plastic_nonlocal_Noutput
@@ -396,11 +378,8 @@ function constitutive_homogenizedC(ipc,ip,el)
  use material, only: &
    phase_plasticity, &
    material_phase, &
-   PLASTICITY_TITANMOD_ID, &
    PLASTICITY_DISLOTWIN_ID, &
    PLASTICITY_DISLOUCLA_ID
- use plastic_titanmod, only: &
-   plastic_titanmod_homogenizedC
  use plastic_dislotwin, only: &
    plastic_dislotwin_homogenizedC
  use lattice, only: &
@@ -416,8 +395,6 @@ function constitutive_homogenizedC(ipc,ip,el)
  plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
    case (PLASTICITY_DISLOTWIN_ID) plasticityType
      constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el)
-   case (PLASTICITY_TITANMOD_ID) plasticityType
-     constitutive_homogenizedC = plastic_titanmod_homogenizedC (ipc,ip,el)
    case default plasticityType
      constitutive_homogenizedC = lattice_C66(1:6,1:6,material_phase (ipc,ip,el))
  end select plasticityType
@@ -438,19 +415,13 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
    thermalMapping, &
    PLASTICITY_dislotwin_ID, &
    PLASTICITY_disloucla_ID, &
-   PLASTICITY_titanmod_ID, &
-   PLASTICITY_nonlocal_ID, &
-   PLASTICITY_phenoplus_ID
- use plastic_titanmod, only: &
-   plastic_titanmod_microstructure
+   PLASTICITY_nonlocal_ID
  use plastic_nonlocal, only: &
    plastic_nonlocal_microstructure
  use plastic_dislotwin, only: &
    plastic_dislotwin_microstructure
  use plastic_disloucla, only: &
    plastic_disloucla_microstructure
- use plastic_phenoplus, only: &
-   plastic_phenoplus_microstructure
 
  implicit none
  integer(pInt), intent(in) :: &
@@ -474,12 +445,8 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
      call plastic_dislotwin_microstructure(temperature(ho)%p(tme),ipc,ip,el)
    case (PLASTICITY_DISLOUCLA_ID) plasticityType
      call plastic_disloucla_microstructure(temperature(ho)%p(tme),ipc,ip,el)
-   case (PLASTICITY_TITANMOD_ID) plasticityType
-     call plastic_titanmod_microstructure (temperature(ho)%p(tme),ipc,ip,el)
    case (PLASTICITY_NONLOCAL_ID) plasticityType
      call plastic_nonlocal_microstructure (Fe,Fp,ip,el)
-   case (PLASTICITY_PHENOPLUS_ID) plasticityType
-     call plastic_phenoplus_microstructure(orientations,ipc,ip,el)
  end select plasticityType
 
 end subroutine constitutive_microstructure
@@ -505,23 +472,17 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
    PLASTICITY_NONE_ID, &
    PLASTICITY_ISOTROPIC_ID, &
    PLASTICITY_PHENOPOWERLAW_ID, &
-   PLASTICITY_PHENOPLUS_ID, &
    PLASTICITY_DISLOTWIN_ID, &
    PLASTICITY_DISLOUCLA_ID, &
-   PLASTICITY_TITANMOD_ID, &
    PLASTICITY_NONLOCAL_ID
  use plastic_isotropic, only: &
    plastic_isotropic_LpAndItsTangent
  use plastic_phenopowerlaw, only: &
    plastic_phenopowerlaw_LpAndItsTangent
- use plastic_phenoplus, only: &
-   plastic_phenoplus_LpAndItsTangent
  use plastic_dislotwin, only: &
    plastic_dislotwin_LpAndItsTangent
  use plastic_disloucla, only: &
    plastic_disloucla_LpAndItsTangent
- use plastic_titanmod, only: &
-   plastic_titanmod_LpAndItsTangent
  use plastic_nonlocal, only: &
    plastic_nonlocal_LpAndItsTangent
 
@@ -564,8 +525,6 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
      call plastic_isotropic_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
    case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
      call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
-   case (PLASTICITY_PHENOPLUS_ID) plasticityType
-     call plastic_phenoplus_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
    case (PLASTICITY_NONLOCAL_ID) plasticityType
      call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v, &
                                            temperature(ho)%p(tme),ip,el)
@@ -575,9 +534,6 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
    case (PLASTICITY_DISLOUCLA_ID) plasticityType
      call plastic_disloucla_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v, &
                                             temperature(ho)%p(tme), ipc,ip,el)
-   case (PLASTICITY_TITANMOD_ID) plasticityType
-     call plastic_titanmod_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v, &
-                                           temperature(ho)%p(tme), ipc,ip,el)
  end select plasticityType
 
  dLp_dTstar3333 = math_Plain99to3333(dLp_dMstar)
@@ -888,10 +844,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
    PLASTICITY_none_ID, &
    PLASTICITY_isotropic_ID, &
    PLASTICITY_phenopowerlaw_ID, &
-   PLASTICITY_phenoplus_ID, &
    PLASTICITY_dislotwin_ID, &
    PLASTICITY_disloucla_ID, &
-   PLASTICITY_titanmod_ID, &
    PLASTICITY_nonlocal_ID, &
    SOURCE_damage_isoDuctile_ID, &
    SOURCE_damage_anisoBrittle_ID, &
@@ -901,14 +855,10 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
    plastic_isotropic_dotState
  use plastic_phenopowerlaw, only: &
    plastic_phenopowerlaw_dotState
- use plastic_phenoplus, only: &
-   plastic_phenoplus_dotState
  use plastic_dislotwin, only: &
    plastic_dislotwin_dotState
  use plastic_disloucla, only: &
    plastic_disloucla_dotState
- use plastic_titanmod, only: &
-   plastic_titanmod_dotState
  use plastic_nonlocal, only: &
    plastic_nonlocal_dotState
  use source_damage_isoDuctile, only: &
@@ -954,17 +904,12 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
      call plastic_isotropic_dotState    (Tstar_v,ipc,ip,el)
    case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
      call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
-   case (PLASTICITY_PHENOPLUS_ID) plasticityType
-     call plastic_phenoplus_dotState    (Tstar_v,ipc,ip,el)
    case (PLASTICITY_DISLOTWIN_ID) plasticityType
      call plastic_dislotwin_dotState    (Tstar_v,temperature(ho)%p(tme), &
                                          ipc,ip,el)
    case (PLASTICITY_DISLOUCLA_ID) plasticityType
      call plastic_disloucla_dotState    (Tstar_v,temperature(ho)%p(tme), &
                                          ipc,ip,el)
-   case (PLASTICITY_TITANMOD_ID) plasticityType
-     call plastic_titanmod_dotState     (Tstar_v,temperature(ho)%p(tme), &
-                                         ipc,ip,el)
    case (PLASTICITY_NONLOCAL_ID) plasticityType
      call plastic_nonlocal_dotState     (Tstar_v,FeArray,FpArray,temperature(ho)%p(tme), &
                                          subdt,subfracArray,ip,el)
@@ -1097,10 +1042,8 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
    PLASTICITY_NONE_ID, &
    PLASTICITY_ISOTROPIC_ID, &
    PLASTICITY_PHENOPOWERLAW_ID, &
-   PLASTICITY_PHENOPLUS_ID, &
    PLASTICITY_DISLOTWIN_ID, &
    PLASTICITY_DISLOUCLA_ID, &
-   PLASTICITY_TITANMOD_ID, &
    PLASTICITY_NONLOCAL_ID, &
    SOURCE_damage_isoBrittle_ID, &
    SOURCE_damage_isoDuctile_ID, &
@@ -1110,14 +1053,10 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
    plastic_isotropic_postResults
  use plastic_phenopowerlaw, only: &
    plastic_phenopowerlaw_postResults
- use plastic_phenoplus, only: &
-   plastic_phenoplus_postResults
  use plastic_dislotwin, only: &
    plastic_dislotwin_postResults
  use plastic_disloucla, only: &
    plastic_disloucla_postResults
- use plastic_titanmod, only: &
-   plastic_titanmod_postResults
  use plastic_nonlocal, only: &
    plastic_nonlocal_postResults
  use source_damage_isoBrittle, only: &
@@ -1157,16 +1096,11 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
  endPos = plasticState(material_phase(ipc,ip,el))%sizePostResults
 
  plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
-   case (PLASTICITY_TITANMOD_ID) plasticityType
-     constitutive_postResults(startPos:endPos) = plastic_titanmod_postResults(ipc,ip,el)
    case (PLASTICITY_ISOTROPIC_ID) plasticityType
      constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
    case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
      constitutive_postResults(startPos:endPos) = &
        plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)
-   case (PLASTICITY_PHENOPLUS_ID) plasticityType
-     constitutive_postResults(startPos:endPos) = &
-       plastic_phenoplus_postResults(Tstar_v,ipc,ip,el)
    case (PLASTICITY_DISLOTWIN_ID) plasticityType
      constitutive_postResults(startPos:endPos) = &
        plastic_dislotwin_postResults(Tstar_v,temperature(ho)%p(tme),ipc,ip,el)

src/crystallite.f90

@@ -554,7 +554,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
    FEsolving_execIP
  use mesh, only: &
    mesh_element, &
-   mesh_NcpElems, &
    mesh_maxNips, &
    mesh_ipNeighborhood, &
    FE_NipNeighbors, &
@@ -565,8 +564,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
    plasticState, &
    sourceState, &
    phase_Nsources, &
-   phaseAt, phasememberAt, &
-   homogenization_maxNgrains
+   phaseAt, phasememberAt
  use constitutive, only:  &
    constitutive_TandItsTangent, &
    constitutive_LpAndItsTangent, &
@@ -794,7 +792,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
                        if (.not. crystallite_localPlasticity(1,neighboring_i,neighboring_e) &
                            .and. .not. crystallite_converged(1,neighboring_i,neighboring_e)) then
                          crystallite_neighborEnforcedCutback(i,e) = .true.
-#ifndef _OPENMP
+#ifdef DEBUG
                          if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) &
                             write(6,'(a12,i5,1x,i2,a,i5,1x,i2)') '<< CRYST >> ', neighboring_e,neighboring_i, &
                                                                  ' enforced cutback at ',e,i
@@ -829,7 +827,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
                      if (.not. crystallite_localPlasticity(1,neighboring_i,neighboring_e) &
                          .and. .not. crystallite_converged(1,neighboring_i,neighboring_e)) then
                        crystallite_syncSubFrac(i,e) = .true.
-#ifndef _OPENMP
+#ifdef DEBUG
                        if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) &
                          write(6,'(a12,i5,1x,i2,a,i5,1x,i2)') '<< CRYST >> ',neighboring_e,neighboring_i, &
                                                              ' enforced time synchronization at ',e,i
@@ -937,7 +935,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
                  crystallite_todo(c,i,e) = .true.
                endif
                !$OMP FLUSH(crystallite_todo)
-#ifndef _OPENMP
+#ifdef DEBUG
                if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt &
                    .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) &
                           .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) &
@@ -987,7 +985,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
                                                                                                      ! cant restore dotState here, since not yet calculated in first cutback after initialization
              crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > subStepMinCryst                  ! still on track or already done (beyond repair)
              !$OMP FLUSH(crystallite_todo)
-#ifndef _OPENMP
+#ifdef DEBUG
              if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt) then
                if (crystallite_todo(c,i,e)) then
                  write(6,'(a,f12.8,a,i8,1x,i2,1x,i3,/)') '<< CRYST >> cutback step in crystallite_stressAndItsTangent &
@@ -1393,7 +1391,7 @@ subroutine crystallite_integrateStateRK4()
          * crystallite_subdt(g,i,e) * timeStepFraction(n)
          enddo
 
-#ifndef _OPENMP
+#ifdef DEBUG
          if (n == 4 &
              .and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
              .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
@@ -1784,7 +1782,7 @@ subroutine crystallite_integrateStateRKCK45()
 
 
    ! --- dot state and RK dot state---
-#ifndef _OPENMP
+#ifdef DEBUG
    if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
      write(6,'(a,1x,i1)') '<< CRYST >> Runge--Kutta step',stage+1_pInt
 #endif
@@ -1933,7 +1931,7 @@ subroutine crystallite_integrateStateRKCK45()
                                        sourceState(p)%p(mySource)%aTolState(1:mySizeSourceDotState))
        enddo
 
-#ifndef _OPENMP
+#ifdef DEBUG
        if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt&
            .and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
                   .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
@@ -2317,7 +2315,7 @@ subroutine crystallite_integrateStateAdaptiveEuler()
          !$OMP FLUSH(relPlasticStateResiduum)
          !$OMP FLUSH(relSourceStateResiduum)
 
-#ifndef _OPENMP
+#ifdef DEBUG
 
          if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
              .and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
@@ -2513,7 +2511,7 @@ eIter = FEsolving_execElem(1:2)
          * crystallite_subdt(g,i,e)
          enddo
 
-#ifndef _OPENMP
+#ifdef DEBUG
          if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
              .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
                      .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
@@ -2962,7 +2960,7 @@ subroutine crystallite_integrateStateFPI()
            * (1.0_pReal - sourceStateDamper)
          enddo
 
-#ifndef _OPENMP
+#ifdef DEBUG
          if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
              .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
                     .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
@@ -3134,7 +3132,7 @@ logical function crystallite_stateJump(ipc,ip,el)
      sourceState(p)%p(mySource)%deltaState(1:mySizeSourceDeltaState,c)
  enddo
 
-#ifndef _OPENMP
+#ifdef DEBUG
  if (any(dNeq0(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c))) &
      .and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
      .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
@@ -3309,7 +3307,7 @@ logical function crystallite_integrateStress(&
 
  !* be pessimistic
  crystallite_integrateStress = .false.
-#ifndef _OPENMP
+#ifdef DEBUG
  if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
      .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
             .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) &
@@ -3342,9 +3340,9 @@ logical function crystallite_integrateStress(&
 
  invFp_current = math_inv33(Fp_current)
  failedInversionFp: if (all(dEq0(invFp_current))) then
-#ifndef _OPENMP
+#ifdef DEBUG
    if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
-     write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fp_current at el (elFE) ip g ',&
+     write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fp_current at el (elFE) ip ipc ',&
        el,'(',mesh_element(1,el),')',ip,ipc
      if (iand(debug_level(debug_crystallite), debug_levelExtensive) > 0_pInt) &
        write(6,'(/,a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp_current',math_transpose33(Fp_current(1:3,1:3))
@@ -3358,7 +3356,7 @@ logical function crystallite_integrateStress(&
 
  invFi_current = math_inv33(Fi_current)
  failedInversionFi: if (all(dEq0(invFi_current))) then
-#ifndef _OPENMP
+#ifdef DEBUG
    if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
      write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fi_current at el (elFE) ip ipc ',&
        el,'(',mesh_element(1,el),')',ip,ipc
@@ -3379,10 +3377,10 @@ logical function crystallite_integrateStress(&
  LiLoop: do
    NiterationStressLi = NiterationStressLi + 1_pInt
    IloopsExeced: if (NiterationStressLi > nStress) then
-#ifndef _OPENMP
+#ifdef DEBUG
      if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
        write(6,'(a,i3,a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> integrateStress reached inelastic loop limit',nStress, &
-       ' at el (elFE) ip ipc ', el,mesh_element(1,el),ip,ipc
+       ' at el (elFE) ip ipc ', el,'(',mesh_element(1,el),')',ip,ipc
 #endif
      return
    endif IloopsExeced
@@ -3400,7 +3398,7 @@ logical function crystallite_integrateStress(&
    LpLoop: do                                    ! inner stress integration loop for consistency with Fi
      NiterationStressLp = NiterationStressLp + 1_pInt
      loopsExeced: if (NiterationStressLp > nStress) then
-#ifndef _OPENMP
+#ifdef DEBUG
      if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
          write(6,'(a,i3,a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> integrateStress reached loop limit',nStress, &
          ' at el (elFE) ip ipc ', el,mesh_element(1,el),ip,ipc
@@ -3433,7 +3431,7 @@ logical function crystallite_integrateStress(&
        !$OMP END CRITICAL (debugTimingLpTangent)
      endif
 
-#ifndef _OPENMP
+#ifdef DEBUG
      if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
          .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
                 .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
@@ -3451,10 +3449,10 @@ logical function crystallite_integrateStress(&
      residuumLp = Lpguess - Lp_constitutive
 
      if (any(IEEE_is_NaN(residuumLp))) then                                                          ! NaN in residuum...
-#ifndef _OPENMP
+#ifdef DEBUG
        if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
          write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el (elFE) ip ipc ', &
-           el,mesh_element(1,el),ip,ipc, &
+           el,'(',mesh_element(1,el),')',ip,ipc, &
                                                         ' ; iteration ', NiterationStressLp,&
                                                         ' >> returning..!'
 #endif
@@ -3486,10 +3484,10 @@ logical function crystallite_integrateStress(&
        work = math_plain33to9(residuumLp)
        call dgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr)                                                   ! solve dRLp/dLp * delta Lp = -res for delta Lp
        if (ierr /= 0_pInt) then
-#ifndef _OPENMP
+#ifdef DEBUG
          if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
-           write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLp inversion at el ip ipc ', &
-             el,mesh_element(1,el),ip,ipc
+           write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on dR/dLp inversion at el (elFE) ip ipc ', &
+             el,'(',mesh_element(1,el),')',ip,ipc
            if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
                .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g)&
                       .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
@@ -3527,7 +3525,7 @@ logical function crystallite_integrateStress(&
    call constitutive_LiAndItsTangent(Li_constitutive, dLi_dT3333, dLi_dFi3333, &
                                      Tstar_v, Fi_new, ipc, ip, el)
 
-#ifndef _OPENMP
+#ifdef DEBUG
      if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
          .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
                 .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
@@ -3575,10 +3573,10 @@ logical function crystallite_integrateStress(&
      work = math_plain33to9(residuumLi)
      call dgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr)                                                    ! solve dRLi/dLp * delta Li = -res for delta Li
      if (ierr /= 0_pInt) then
-#ifndef _OPENMP
+#ifdef DEBUG
        if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
-         write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLi inversion at el ip ipc ', &
-               el,mesh_element(1,el),ip,ipc
+         write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on dR/dLi inversion at el (elFE) ip ipc ', &
+               el,'(',mesh_element(1,el),')',ip,ipc
          if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
              .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g)&
                     .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
@@ -3615,10 +3613,10 @@ logical function crystallite_integrateStress(&
  invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal)                               ! regularize by det
  Fp_new = math_inv33(invFp_new)
  failedInversionInvFp: if (all(dEq0(Fp_new))) then
-#ifndef _OPENMP
+#ifdef DEBUG
    if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
-     write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip ipc ',&
-       el,mesh_element(1,el),ip,ipc, ' ; iteration ', NiterationStressLp
+     write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el (elFE) ip ipc ',&
+       el,'(',mesh_element(1,el),')',ip,ipc, ' ; iteration ', NiterationStressLp
      if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
          .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
                 .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) &
@@ -3649,7 +3647,7 @@ logical function crystallite_integrateStress(&
  !* set return flag to true
 
  crystallite_integrateStress = .true.
-#ifndef _OPENMP
+#ifdef DEBUG
  if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt &
      .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
              .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then

src/damage_local.f90

@@ -72,8 +72,6 @@ subroutine damage_local_init(fileUnit)
    damage, &
    damage_initialPhi, &
    material_partHomogenization
- use numerics,only: &
-   worldrank
  
  implicit none
  integer(pInt), intent(in) :: fileUnit
@@ -86,11 +84,9 @@ subroutine damage_local_init(fileUnit)
    tag  = '', &
    line = ''
 
- mainProcess: if (worldrank == 0) then 
  write(6,'(/,a)')   ' <<<+-  damage_'//DAMAGE_local_label//' init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
- endif mainProcess
 
  maxNinstance = int(count(damage_type == DAMAGE_local_ID),pInt)
  if (maxNinstance == 0_pInt) return

src/damage_none.f90

@@ -26,19 +26,15 @@ subroutine damage_none_init()
  use IO, only: &
    IO_timeStamp
  use material
- use numerics, only: &
-   worldrank
  
  implicit none
  integer(pInt) :: &
    homog, &
    NofMyHomog
 
- mainProcess: if (worldrank == 0) then 
  write(6,'(/,a)')   ' <<<+-  damage_'//DAMAGE_none_label//' init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
- endif mainProcess
 
   initializeInstances: do homog = 1_pInt, material_Nhomogenization
    

src/damage_nonlocal.f90

@@ -77,8 +77,6 @@ subroutine damage_nonlocal_init(fileUnit)
    damage, &
    damage_initialPhi, &
    material_partHomogenization
- use numerics,only: &
-   worldrank
 
  implicit none
  integer(pInt), intent(in) :: fileUnit
@@ -91,11 +89,9 @@ subroutine damage_nonlocal_init(fileUnit)
    tag  = '', &
    line = ''
 
- mainProcess: if (worldrank == 0) then 
  write(6,'(/,a)')   ' <<<+-  damage_'//DAMAGE_nonlocal_label//' init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
- endif mainProcess
  
  maxNinstance = int(count(damage_type == DAMAGE_nonlocal_ID),pInt)
  if (maxNinstance == 0_pInt) return

src/homogenization.f90

@@ -542,6 +542,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
    debug_level, &
    debug_homogenization, &
    debug_levelBasic, &
+   debug_levelExtensive, &
    debug_levelSelective, &
    debug_e, &
    debug_i, &
@@ -638,8 +639,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
      IpLooping1: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
 
        converged: if ( materialpoint_converged(i,e) ) then
-#ifndef _OPENMP
-         if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt &
+#ifdef DEBUG
+         if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt &
             .and. ((e == debug_e .and. i == debug_i) &
                    .or. .not. iand(debug_level(debug_homogenization),debug_levelSelective) /= 0_pInt)) then
            write(6,'(a,1x,f12.8,1x,a,1x,f12.8,1x,a,i8,1x,i2/)') '<< HOMOG >> winding forward from', &
@@ -741,8 +742,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
            materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e)               ! crystallite had severe trouble, so do a significant cutback
            !$OMP FLUSH(materialpoint_subStep)
 
-#ifndef _OPENMP
-           if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt &
+#ifdef DEBUG
+           if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt &
               .and. ((e == debug_e .and. i == debug_i) &
                     .or. .not. iand(debug_level(debug_homogenization), debug_levelSelective) /= 0_pInt)) then
              write(6,'(a,1x,f12.8,a,i8,1x,i2/)') &

src/homogenization_RGC.f90

@@ -100,8 +100,6 @@ subroutine homogenization_RGC_init(fileUnit)
    FE_geomtype
  use IO
  use material
- use numerics, only: &
-   worldrank
 
  implicit none
  integer(pInt), intent(in) :: fileUnit                                                              !< file pointer to material configuration
@@ -117,11 +115,9 @@ subroutine homogenization_RGC_init(fileUnit)
    tag = '', &
    line = ''
  
- mainProcess: if (worldrank == 0) then 
  write(6,'(/,a)')   ' <<<+-  homogenization_'//HOMOGENIZATION_RGC_label//' init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
- endif mainProcess
 
  maxNinstance = int(count(homogenization_type == HOMOGENIZATION_RGC_ID),pInt)
  if (maxNinstance == 0_pInt) return

src/homogenization_isostrain.f90

@@ -62,8 +62,6 @@ subroutine homogenization_isostrain_init(fileUnit)
    debug_levelBasic
  use IO
  use material
- use numerics, only: &
-   worldrank
  
  implicit none
  integer(pInt),                                      intent(in) :: fileUnit
@@ -80,11 +78,9 @@ subroutine homogenization_isostrain_init(fileUnit)
    tag  = '', &
    line = ''
  
- mainProcess: if (worldrank == 0) then 
  write(6,'(/,a)')   ' <<<+-  homogenization_'//HOMOGENIZATION_ISOSTRAIN_label//' init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
- endif mainProcess
 
  maxNinstance = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
  if (maxNinstance == 0) return

src/homogenization_none.f90

@@ -29,19 +29,15 @@ subroutine homogenization_none_init()
  use IO, only: &
    IO_timeStamp
  use material
- use numerics, only: &
-   worldrank
  
  implicit none
  integer(pInt) :: &
    homog, &
    NofMyHomog
 
- mainProcess: if (worldrank == 0) then 
  write(6,'(/,a)')   ' <<<+-  homogenization_'//HOMOGENIZATION_NONE_label//' init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
- endif mainProcess
 
  initializeInstances: do homog = 1_pInt, material_Nhomogenization
    

src/hydrogenflux_cahnhilliard.f90

@@ -84,8 +84,6 @@ subroutine hydrogenflux_cahnhilliard_init(fileUnit)
    hydrogenflux_initialCh, &
    material_partHomogenization, &
    material_partPhase
- use numerics,only: &
-   worldrank
 
  implicit none
  integer(pInt), intent(in) :: fileUnit
@@ -98,11 +96,9 @@ subroutine hydrogenflux_cahnhilliard_init(fileUnit)
    tag  = '', &
    line = ''
 
- mainProcess: if (worldrank == 0) then 
  write(6,'(/,a)')   ' <<<+-  hydrogenflux_'//HYDROGENFLUX_cahnhilliard_label//' init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
- endif mainProcess
  
  maxNinstance = int(count(hydrogenflux_type == HYDROGENFLUX_cahnhilliard_ID),pInt)
  if (maxNinstance == 0_pInt) return

src/hydrogenflux_isoconc.f90

@@ -27,19 +27,15 @@ subroutine hydrogenflux_isoconc_init()
  use IO, only: &
    IO_timeStamp
  use material
- use numerics, only: &
-   worldrank
  
  implicit none
  integer(pInt) :: &
    homog, &
    NofMyHomog
 
- mainProcess: if (worldrank == 0) then 
  write(6,'(/,a)')   ' <<<+-  hydrogenflux_'//HYDROGENFLUX_isoconc_label//' init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
- endif mainProcess
 
  initializeInstances: do homog = 1_pInt, material_Nhomogenization
    

src/kinematics_cleavage_opening.f90

@@ -81,8 +81,6 @@ subroutine kinematics_cleavage_opening_init(fileUnit)
    KINEMATICS_cleavage_opening_ID, &
    material_Nphase, &
    MATERIAL_partPhase
- use numerics,only: &
-   worldrank
  use lattice, only: &
    lattice_maxNcleavageFamily, &
    lattice_NcleavageSystem
@@ -97,11 +95,9 @@ subroutine kinematics_cleavage_opening_init(fileUnit)
    tag  = '', &
    line = ''
 
- mainProcess: if (worldrank == 0) then 
  write(6,'(/,a)')   ' <<<+-  kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
- endif mainProcess
 
  maxNinstance = int(count(phase_kinematics == KINEMATICS_cleavage_opening_ID),pInt)
  if (maxNinstance == 0_pInt) return

src/kinematics_slipplane_opening.f90

@@ -81,8 +81,6 @@ subroutine kinematics_slipplane_opening_init(fileUnit)
    KINEMATICS_slipplane_opening_ID, &
    material_Nphase, &
    MATERIAL_partPhase
- use numerics,only: &
-   worldrank
  use lattice, only: &
    lattice_maxNslipFamily, &
    lattice_NslipSystem
@@ -97,11 +95,9 @@ subroutine kinematics_slipplane_opening_init(fileUnit)
    tag  = '', &
    line = ''
 
- mainProcess: if (worldrank == 0) then 
  write(6,'(/,a)')   ' <<<+-  kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
- endif mainProcess
 
  maxNinstance = int(count(phase_kinematics == KINEMATICS_slipplane_opening_ID),pInt)
  if (maxNinstance == 0_pInt) return

src/kinematics_thermal_expansion.f90

@@ -71,8 +71,6 @@ subroutine kinematics_thermal_expansion_init(fileUnit)
    KINEMATICS_thermal_expansion_ID, &
    material_Nphase, &
    MATERIAL_partPhase
- use numerics,only: &
-   worldrank
 
  implicit none
  integer(pInt), intent(in) :: fileUnit
@@ -83,11 +81,9 @@ subroutine kinematics_thermal_expansion_init(fileUnit)
    tag     = '', &
    line    = ''
 
- mainProcess: if (worldrank == 0) then 
  write(6,'(/,a)')   ' <<<+-  kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
- endif mainProcess
 
  maxNinstance = int(count(phase_kinematics == KINEMATICS_thermal_expansion_ID),pInt)
  if (maxNinstance == 0_pInt) return

src/kinematics_vacancy_strain.f90

@@ -71,8 +71,6 @@ subroutine kinematics_vacancy_strain_init(fileUnit)
    KINEMATICS_vacancy_strain_ID, &
    material_Nphase, &
    MATERIAL_partPhase
- use numerics,only: &
-   worldrank
 
  implicit none
  integer(pInt), intent(in) :: fileUnit
@@ -83,11 +81,9 @@ subroutine kinematics_vacancy_strain_init(fileUnit)
    tag  = '', &
    line = ''
 
- mainProcess: if (worldrank == 0) then 
  write(6,'(/,a)')   ' <<<+-  kinematics_'//KINEMATICS_vacancy_strain_LABEL//' init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
- endif mainProcess
 
  maxNinstance = int(count(phase_kinematics == KINEMATICS_vacancy_strain_ID),pInt)
  if (maxNinstance == 0_pInt) return

src/lattice.f90

@@ -96,19 +96,19 @@ module lattice
 
  real(pReal), dimension(3+3,LATTICE_fcc_Nslip), parameter, private :: &
    LATTICE_fcc_systemSlip = reshape(real([&
-    ! Slip direction     Plane normal
-      0, 1,-1,     1, 1, 1, &
-     -1, 0, 1,     1, 1, 1, &
-      1,-1, 0,     1, 1, 1, &
-      0,-1,-1,    -1,-1, 1, &
-      1, 0, 1,    -1,-1, 1, &
-     -1, 1, 0,    -1,-1, 1, &
-      0,-1, 1,     1,-1,-1, &
-     -1, 0,-1,     1,-1,-1, &
-      1, 1, 0,     1,-1,-1, &
-      0, 1, 1,    -1, 1,-1, &
-      1, 0,-1,    -1, 1,-1, &
-     -1,-1, 0,    -1, 1,-1  &
+    ! Slip direction     Plane normal                                                               ! SCHMID-BOAS notation
+      0, 1,-1,     1, 1, 1, &                                                                       ! B2
+     -1, 0, 1,     1, 1, 1, &                                                                       ! B4
+      1,-1, 0,     1, 1, 1, &                                                                       ! B5
+      0,-1,-1,    -1,-1, 1, &                                                                       ! C1
+      1, 0, 1,    -1,-1, 1, &                                                                       ! C3
+     -1, 1, 0,    -1,-1, 1, &                                                                       ! C5
+      0,-1, 1,     1,-1,-1, &                                                                       ! A2
+     -1, 0,-1,     1,-1,-1, &                                                                       ! A3
+      1, 1, 0,     1,-1,-1, &                                                                       ! A6
+      0, 1, 1,    -1, 1,-1, &                                                                       ! D1
+      1, 0,-1,    -1, 1,-1, &                                                                       ! D4
+     -1,-1, 0,    -1, 1,-1  &                                                                       ! D6
      ],pReal),[ 3_pInt + 3_pInt,LATTICE_fcc_Nslip])                                                 !< Slip system <110>{111} directions. Sorted according to Eisenlohr & Hantcherli
 
  real(pReal), dimension(3+3,LATTICE_fcc_Ntwin), parameter, private :: &

src/material.f90

@@ -25,10 +25,8 @@ module material
    PLASTICITY_none_label                = 'none', &
    PLASTICITY_isotropic_label           = 'isotropic', &
    PLASTICITY_phenopowerlaw_label       = 'phenopowerlaw', &
-   PLASTICITY_phenoplus_label           = 'phenoplus', &
    PLASTICITY_dislotwin_label           = 'dislotwin', &
    PLASTICITY_disloucla_label           = 'disloucla', &
-   PLASTICITY_titanmod_label            = 'titanmod', &
    PLASTICITY_nonlocal_label            = 'nonlocal', &
    SOURCE_thermal_dissipation_label     = 'thermal_dissipation', &
    SOURCE_thermal_externalheat_label    = 'thermal_externalheat', &
@@ -74,10 +72,8 @@ module material
                  PLASTICITY_none_ID, &
                  PLASTICITY_isotropic_ID, &
                  PLASTICITY_phenopowerlaw_ID, &
-                 PLASTICITY_phenoplus_ID, &
                  PLASTICITY_dislotwin_ID, &
                  PLASTICITY_disloucla_ID, &
-                 PLASTICITY_titanmod_ID, &
                  PLASTICITY_nonlocal_ID
  end enum
 
@@ -312,10 +308,8 @@ module material
    PLASTICITY_none_ID, &
    PLASTICITY_isotropic_ID, &
    PLASTICITY_phenopowerlaw_ID, &
-   PLASTICITY_phenoplus_ID, &
    PLASTICITY_dislotwin_ID, &
    PLASTICITY_disloucla_ID, &
-   PLASTICITY_titanmod_ID, &
    PLASTICITY_nonlocal_ID, &
    SOURCE_thermal_dissipation_ID, &
    SOURCE_thermal_externalheat_ID, &
@@ -989,14 +983,10 @@ subroutine material_parsePhase(fileUnit,myPart)
              phase_plasticity(section) = PLASTICITY_ISOTROPIC_ID
            case (PLASTICITY_PHENOPOWERLAW_label)
              phase_plasticity(section) = PLASTICITY_PHENOPOWERLAW_ID
-           case (PLASTICITY_PHENOPLUS_label)
-             phase_plasticity(section) = PLASTICITY_PHENOPLUS_ID
            case (PLASTICITY_DISLOTWIN_label)
              phase_plasticity(section) = PLASTICITY_DISLOTWIN_ID
            case (PLASTICITY_DISLOUCLA_label)
              phase_plasticity(section) = PLASTICITY_DISLOUCLA_ID
-           case (PLASTICITY_TITANMOD_label)
-             phase_plasticity(section) = PLASTICITY_TITANMOD_ID
            case (PLASTICITY_NONLOCAL_label)
              phase_plasticity(section) = PLASTICITY_NONLOCAL_ID
            case default

src/math.f90

@@ -178,7 +178,7 @@ subroutine math_init
    compiler_version, &
    compiler_options
 #endif
- use numerics, only: fixedSeed
+ use numerics, only: randomSeed
  use IO,       only: IO_timeStamp
 
  implicit none
@@ -195,8 +195,8 @@ subroutine math_init
  call random_seed(size=randSize)
  if (allocated(randInit)) deallocate(randInit)
  allocate(randInit(randSize))
- if (fixedSeed > 0_pInt) then
-   randInit(1:randSize) = int(fixedSeed)                                                            ! fixedSeed is of type pInt, randInit not
+ if (randomSeed > 0_pInt) then
+   randInit(1:randSize) = int(randomSeed)                                                           ! randomSeed is of type pInt, randInit not
    call random_seed(put=randInit)
  else
    call random_seed()
@@ -1440,34 +1440,36 @@ end function math_RtoQ
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief rotation matrix from Euler angles (in radians)
-!> @details rotation matrix is meant to represent a PASSIVE rotation, 
-!> @details composed of INTRINSIC rotations around the axes of the 
-!> @details rotating reference frame 
-!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
+!> @brief rotation matrix from Bunge-Euler (3-1-3) angles (in radians)
+!> @details rotation matrix is meant to represent a PASSIVE rotation, composed of INTRINSIC
+!> @details rotations around the axes of the details rotating reference frame.
+!> @details similar to eu2om from "D Rowenhorst et al. Consistent representations of and conversions
+!> @details between 3D rotations, Model. Simul. Mater. Sci. Eng. 23-8 (2015)", but R is transposed
 !--------------------------------------------------------------------------------------------------
 pure function math_EulerToR(Euler)
 
  implicit none
  real(pReal), dimension(3), intent(in) :: Euler
  real(pReal), dimension(3,3) :: math_EulerToR
- real(pReal) c1, c, c2, s1, s, s2
+ real(pReal) :: c1, C, c2, s1, S, s2
 
- C1 = cos(Euler(1))
+ c1 = cos(Euler(1))
  C  = cos(Euler(2))
- C2 = cos(Euler(3))
- S1 = sin(Euler(1))
+ c2 = cos(Euler(3))
+ s1 = sin(Euler(1))
  S  = sin(Euler(2))
- S2 = sin(Euler(3))
+ s2 = sin(Euler(3))
 
- math_EulerToR(1,1)=C1*C2-S1*S2*C
- math_EulerToR(1,2)=-C1*S2-S1*C2*C
- math_EulerToR(1,3)=S1*S
- math_EulerToR(2,1)=S1*C2+C1*S2*C
- math_EulerToR(2,2)=-S1*S2+C1*C2*C
- math_EulerToR(2,3)=-C1*S
- math_EulerToR(3,1)=S2*S
- math_EulerToR(3,2)=C2*S
+ math_EulerToR(1,1) =  c1*c2 -s1*C*s2
+ math_EulerToR(1,2) = -c1*s2 -s1*C*c2
+ math_EulerToR(1,3) =  s1*S
+
+ math_EulerToR(2,1) =  s1*c2 +c1*C*s2
+ math_EulerToR(2,2) = -s1*s2 +c1*C*c2
+ math_EulerToR(2,3) = -c1*S
+
+ math_EulerToR(3,1) =  S*s2
+ math_EulerToR(3,2) =  S*c2
  math_EulerToR(3,3) =  C
  
  math_EulerToR = transpose(math_EulerToR)  ! convert to passive rotation
@@ -1476,29 +1478,29 @@ end function math_EulerToR
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief quaternion (w+ix+jy+kz) from 3-1-3 Euler angles (in radians)
-!> @details quaternion is meant to represent a PASSIVE rotation, 
-!> @details composed of INTRINSIC rotations around the axes of the 
-!> @details rotating reference frame 
-!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
+!> @brief quaternion (w+ix+jy+kz) from Bunge-Euler (3-1-3) angles (in radians)
+!> @details rotation matrix is meant to represent a PASSIVE rotation, composed of INTRINSIC
+!> @details rotations around the axes of the details rotating reference frame.
+!> @details similar to eu2qu from "D Rowenhorst et al. Consistent representations of and
+!> @details conversions between 3D rotations, Model. Simul. Mater. Sci. Eng. 23-8 (2015)", but
+!> @details Q is conjucated and Q is not reversed for Q(0) < 0.
 !--------------------------------------------------------------------------------------------------
 pure function math_EulerToQ(eulerangles)
 
  implicit none
  real(pReal), dimension(3), intent(in) :: eulerangles
  real(pReal), dimension(4) :: math_EulerToQ
- real(pReal), dimension(3) :: halfangles
- real(pReal) :: c, s
+ real(pReal) :: c, s, sigma, delta
 
- halfangles = 0.5_pReal * eulerangles
+ c = cos(0.5_pReal * eulerangles(2))
+ s = sin(0.5_pReal * eulerangles(2))
+ sigma = 0.5_pReal * (eulerangles(1)+eulerangles(3))
+ delta = 0.5_pReal * (eulerangles(1)-eulerangles(3))
  
- c = cos(halfangles(2))
- s = sin(halfangles(2))
-
- math_EulerToQ= [cos(halfangles(1)+halfangles(3)) * c, &
-                 cos(halfangles(1)-halfangles(3)) * s, &
-                 sin(halfangles(1)-halfangles(3)) * s, &
-                 sin(halfangles(1)+halfangles(3)) * c ]
+ math_EulerToQ= [c * cos(sigma), &
+                 s * cos(delta), &
+                 s * sin(delta), &
+                 c * sin(sigma) ]
  math_EulerToQ = math_qConj(math_EulerToQ)                                                          ! convert to passive rotation
 
 end function math_EulerToQ
@@ -1509,6 +1511,8 @@ end function math_EulerToQ
 !> @details rotation matrix is meant to represent a ACTIVE rotation
 !> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
 !> @details formula for active rotation taken from http://mathworld.wolfram.com/RodriguesRotationFormula.html
+!> @details equivalent to eu2om (P=-1) from "D Rowenhorst et al. Consistent representations of and
+!> @details conversions between 3D rotations, Model. Simul. Mater. Sci. Eng. 23-8 (2015)"
 !--------------------------------------------------------------------------------------------------
 pure function math_axisAngleToR(axis,omega)
 
@@ -1516,31 +1520,31 @@ pure function math_axisAngleToR(axis,omega)
  real(pReal), dimension(3,3) :: math_axisAngleToR
  real(pReal), dimension(3), intent(in) :: axis
  real(pReal), intent(in) :: omega
- real(pReal), dimension(3) :: axisNrm
+ real(pReal), dimension(3) :: n
  real(pReal) :: norm,s,c,c1
 
  norm = norm2(axis)
- if (norm > 1.0e-8_pReal) then                             ! non-zero rotation
-   axisNrm = axis/norm                                     ! normalize axis to be sure
+ wellDefined: if (norm > 1.0e-8_pReal) then
+   n = axis/norm                                           ! normalize axis to be sure
 
    s = sin(omega)
    c = cos(omega)
    c1 = 1.0_pReal - c
 
-   math_axisAngleToR(1,1) =  c + c1*axisNrm(1)**2.0_pReal
-   math_axisAngleToR(1,2) = -s*axisNrm(3) + c1*axisNrm(1)*axisNrm(2)
-   math_axisAngleToR(1,3) =  s*axisNrm(2) + c1*axisNrm(1)*axisNrm(3)
+   math_axisAngleToR(1,1) =  c + c1*n(1)**2.0_pReal
+   math_axisAngleToR(1,2) =  c1*n(1)*n(2) - s*n(3)
+   math_axisAngleToR(1,3) =  c1*n(1)*n(3) + s*n(2)
                              
-   math_axisAngleToR(2,1) =  s*axisNrm(3) + c1*axisNrm(2)*axisNrm(1)
-   math_axisAngleToR(2,2) =  c + c1*axisNrm(2)**2.0_pReal
-   math_axisAngleToR(2,3) = -s*axisNrm(1) + c1*axisNrm(2)*axisNrm(3)
+   math_axisAngleToR(2,1) =  c1*n(1)*n(2) + s*n(3)
+   math_axisAngleToR(2,2) =  c + c1*n(2)**2.0_pReal
+   math_axisAngleToR(2,3) =  c1*n(2)*n(3) - s*n(1)
                              
-   math_axisAngleToR(3,1) = -s*axisNrm(2) + c1*axisNrm(3)*axisNrm(1)
-   math_axisAngleToR(3,2) =  s*axisNrm(1) + c1*axisNrm(3)*axisNrm(2)
-   math_axisAngleToR(3,3) =  c + c1*axisNrm(3)**2.0_pReal
- else
+   math_axisAngleToR(3,1) =  c1*n(1)*n(3) - s*n(2)
+   math_axisAngleToR(3,2) =  c1*n(2)*n(3) + s*n(1)
+   math_axisAngleToR(3,3) =  c + c1*n(3)**2.0_pReal
+ else wellDefined
    math_axisAngleToR = math_I3
- endif
+ endif wellDefined
  
 end function math_axisAngleToR
 
@@ -1549,6 +1553,8 @@ end function math_axisAngleToR
 !> @brief rotation matrix from axis and angle (in radians)
 !> @details rotation matrix is meant to represent a PASSIVE rotation
 !> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
+!> @details eq-uivalent to eu2qu (P=+1) from "D Rowenhorst et al. Consistent representations of and
+!> @details conversions between 3D rotations, Model. Simul. Mater. Sci. Eng. 23-8 (2015)"
 !--------------------------------------------------------------------------------------------------
 pure function math_EulerAxisAngleToR(axis,omega)
 
@@ -1587,6 +1593,8 @@ end function math_EulerAxisAngleToQ
 !> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
 !> @details formula for active rotation taken from
 !> @details http://en.wikipedia.org/wiki/Rotation_representation_%28mathematics%29#Rodrigues_parameters
+!> @details equivalent to eu2qu (P=+1) from "D Rowenhorst et al. Consistent representations of and
+!> @details conversions between 3D rotations, Model. Simul. Mater. Sci. Eng. 23-8 (2015)"
 !--------------------------------------------------------------------------------------------------
 pure function math_axisAngleToQ(axis,omega)
 
@@ -1597,13 +1605,13 @@ pure function math_axisAngleToQ(axis,omega)
  real(pReal), dimension(3) :: axisNrm
  real(pReal) :: norm
 
- norm = sqrt(math_mul3x3(axis,axis))
- rotation: if (norm > 1.0e-8_pReal) then
+ norm = norm2(axis)
+ wellDefined: if (norm > 1.0e-8_pReal) then
    axisNrm = axis/norm                                                                              ! normalize axis to be sure
    math_axisAngleToQ = [cos(0.5_pReal*omega), sin(0.5_pReal*omega) * axisNrm(1:3)]
- else rotation
+ else wellDefined
    math_axisAngleToQ = [1.0_pReal,0.0_pReal,0.0_pReal,0.0_pReal]
- endif rotation
+ endif wellDefined
 
 end function math_axisAngleToQ
 

src/mesh.f90

@@ -73,8 +73,11 @@ module mesh
 
 #ifdef Marc4DAMASK
  integer(pInt), private :: &
+   MarcVersion, &                                                                                   !< Version of input file format (Marc only)
    hypoelasticTableStyle, &                                                                         !< Table style (Marc only)
    initialcondTableStyle                                                                            !< Table style (Marc only)
+ integer(pInt), dimension(:), allocatable, private :: &
+   Marc_matNumber                                                                                   !< array of material numbers for hypoelastic material (Marc only)
 #endif
 
  integer(pInt), dimension(2), private :: &
@@ -428,7 +431,9 @@ module mesh
    mesh_spectral_build_elements, &
    mesh_spectral_build_ipNeighborhood, &
 #elif defined Marc4DAMASK
+   mesh_marc_get_fileFormat, &
    mesh_marc_get_tableStyles, &
+   mesh_marc_get_matNumber, &
    mesh_marc_count_nodesAndElements, &
    mesh_marc_count_elementSets, &
    mesh_marc_map_elementSets, &
@@ -519,6 +524,8 @@ subroutine mesh_init(ip,el)
  integer(pInt) :: j
  logical :: myDebug
 
+ external :: MPI_comm_size
+
  write(6,'(/,a)')   ' <<<+-  mesh init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
@@ -577,8 +584,14 @@ subroutine mesh_init(ip,el)
 #elif defined Marc4DAMASK
  call IO_open_inputFile(FILEUNIT,modelName)                                                         ! parse info from input file...
  if (myDebug) write(6,'(a)') ' Opened input file'; flush(6)
+ call mesh_marc_get_fileFormat(FILEUNIT)
+ if (myDebug) write(6,'(a)') ' Got input file format'; flush(6)
  call mesh_marc_get_tableStyles(FILEUNIT)
  if (myDebug) write(6,'(a)') ' Got table styles'; flush(6)
+ if (MarcVersion > 12) then
+   call mesh_marc_get_matNumber(FILEUNIT)
+   if (myDebug) write(6,'(a)') ' Got hypoleastic material number'; flush(6)
+ endif
  call mesh_marc_count_nodesAndElements(FILEUNIT)
  if (myDebug) write(6,'(a)') ' Counted nodes/elements'; flush(6)
  call mesh_marc_count_elementSets(FILEUNIT)
@@ -1532,6 +1545,37 @@ end function mesh_nodesAroundCentres
 #endif
 
 #ifdef Marc4DAMASK
+!--------------------------------------------------------------------------------------------------
+!> @brief Figures out version of Marc input file format and stores ist as MarcVersion
+!--------------------------------------------------------------------------------------------------
+subroutine mesh_marc_get_fileFormat(fileUnit)
+ use IO, only: &
+   IO_lc, &
+   IO_intValue, &
+   IO_stringValue, &
+   IO_stringPos
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ character(len=300) line
+
+610 FORMAT(A300)
+
+ rewind(fileUnit)
+ do
+   read (fileUnit,610,END=620) line
+   chunkPos = IO_stringPos(line)
+   if ( IO_lc(IO_stringValue(line,chunkPos,1_pInt)) == 'version') then
+     MarcVersion = IO_intValue(line,chunkPos,2_pInt)
+     exit
+   endif
+ enddo
+
+620 end subroutine mesh_marc_get_fileFormat
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief Figures out table styles (Marc only) and stores to 'initialcondTableStyle' and
 !! 'hypoelasticTableStyle'
@@ -1568,6 +1612,52 @@ subroutine mesh_marc_get_tableStyles(fileUnit)
 
 620 end subroutine mesh_marc_get_tableStyles
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Figures out material number of hypoelastic material and stores it in Marc_matNumber array
+!--------------------------------------------------------------------------------------------------
+subroutine mesh_marc_get_matNumber(fileUnit)
+ use IO, only: &
+   IO_lc, &
+   IO_intValue, &
+   IO_stringValue, &
+   IO_stringPos
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: i, j, data_blocks
+ character(len=300) line
+
+610 FORMAT(A300)
+
+ rewind(fileUnit)
+
+ data_blocks = 1_pInt
+ do
+  read (fileUnit,610,END=620) line
+  chunkPos = IO_stringPos(line)
+  if ( IO_lc(IO_stringValue(line,chunkPos,1_pInt)) == 'hypoelastic') then
+    read (fileUnit,610,END=620) line
+    if (len(trim(line))/=0_pInt) then
+      chunkPos = IO_stringPos(line)
+      data_blocks = IO_intValue(line,chunkPos,1_pInt)
+    endif
+    allocate(Marc_matNumber(data_blocks))
+    do i=1_pInt,data_blocks                                                                        ! read all data blocks
+      read (fileUnit,610,END=620) line
+      chunkPos = IO_stringPos(line)
+      Marc_matNumber(i) = IO_intValue(line,chunkPos,1_pInt)
+      do j=1_pint,2_pInt + hypoelasticTableStyle                                                   ! read 2 or 3 remaining lines of data block
+        read (fileUnit,610,END=620) line
+      enddo
+    enddo
+    exit
+  endif
+ enddo
+
+620 end subroutine mesh_marc_get_matNumber
+
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Count overall number of nodes and elements in mesh and stores the numbers in
@@ -1697,13 +1787,14 @@ subroutine mesh_marc_count_cpElements(fileUnit)
                  IO_stringPos, &
                  IO_countContinuousIntValues, &
                  IO_error, &
-                 IO_intValue
+                 IO_intValue, &
+                 IO_countNumericalDataLines
 
  implicit none
  integer(pInt), intent(in) :: fileUnit
 
  integer(pInt), allocatable, dimension(:) :: chunkPos
- integer(pInt) :: i, version
+ integer(pInt) :: i
  character(len=300):: line
 
  mesh_NcpElems = 0_pInt
@@ -1711,13 +1802,7 @@ subroutine mesh_marc_count_cpElements(fileUnit)
 610 FORMAT(A300)
 
  rewind(fileUnit)
- do 
-   read (fileUnit,610,END=620) line
-   chunkPos = IO_stringPos(line)
-   if ( IO_lc(IO_stringValue(line,chunkPos,1_pInt)) == 'version') then
-     version = IO_intValue(line,chunkPos,2_pInt)
-     if (version < 13) then                                                                             ! Marc 2016 or earlier
-       rewind(fileUnit)
+ if (MarcVersion < 13) then                                                                             ! Marc 2016 or earlier
    do
      read (fileUnit,610,END=620) line
      chunkPos = IO_stringPos(line)
@@ -1725,15 +1810,23 @@ subroutine mesh_marc_count_cpElements(fileUnit)
        do i=1_pInt,3_pInt+hypoelasticTableStyle  ! Skip 3 or 4 lines
          read (fileUnit,610,END=620) line
        enddo
-           mesh_NcpElems = mesh_NcpElems + IO_countContinuousIntValues(fileUnit)                        ! why not simply mesh_NcpElems = IO_countContinuousIntValues(fileUnit)? keyword hypoelastic might appear several times
+       mesh_NcpElems = mesh_NcpElems + IO_countContinuousIntValues(fileUnit)                            ! why not simply mesh_NcpElems = IO_countContinuousIntValues(fileUnit)? not fully correct as hypoelastic can have multiple data fields, needs update
        exit
      endif
    enddo
  else                                                                                                   ! Marc2017 and later
-       call IO_error(error_ID=701_pInt)
+   do
+     read (fileUnit,610,END=620) line
+     chunkPos = IO_stringPos(line)
+     if ( IO_lc(IO_stringValue(line,chunkPos,1_pInt)) == 'connectivity') then
+       read (fileUnit,610,END=620) line
+       chunkPos = IO_stringPos(line)
+       if (any(Marc_matNumber==IO_intValue(line,chunkPos,6_pInt))) then
+         mesh_NcpElems = mesh_NcpElems + IO_countNumericalDataLines(fileUnit)
        endif
      endif
    enddo
+ end if
 
 620 end subroutine mesh_marc_count_cpElements
 
@@ -1746,6 +1839,7 @@ subroutine mesh_marc_map_elements(fileUnit)
 
  use math, only: math_qsort
  use IO,   only: IO_lc, &
+                 IO_intValue, &
                  IO_stringValue, &
                  IO_stringPos, &
                  IO_continuousIntValues
@@ -1754,7 +1848,8 @@ subroutine mesh_marc_map_elements(fileUnit)
  integer(pInt), intent(in) :: fileUnit
 
  integer(pInt), allocatable, dimension(:) :: chunkPos
- character(len=300) :: line
+ character(len=300) :: line, &
+                       tmp
 
  integer(pInt), dimension (1_pInt+mesh_NcpElems) :: contInts
  integer(pInt) :: i,cpElem = 0_pInt
@@ -1763,25 +1858,47 @@ subroutine mesh_marc_map_elements(fileUnit)
 
 610 FORMAT(A300)
 
+ contInts = 0_pInt
  rewind(fileUnit)
  do
    read (fileUnit,610,END=660) line
    chunkPos = IO_stringPos(line)
+   if (MarcVersion < 13) then                                                                             ! Marc 2016 or earlier
      if( IO_lc(IO_stringValue(line,chunkPos,1_pInt)) == 'hypoelastic' ) then
        do i=1_pInt,3_pInt+hypoelasticTableStyle                                                           ! skip three (or four if new table style!) lines
          read (fileUnit,610,END=660) line
        enddo
        contInts = IO_continuousIntValues(fileUnit,mesh_NcpElems,mesh_nameElemSet,&
                                               mesh_mapElemSet,mesh_NelemSets)
-     do i = 1_pInt,contInts(1)
+       exit
+     endif  
+   else                                                                                                   ! Marc2017 and later
+     if ( IO_lc(IO_stringValue(line,chunkPos,1_pInt)) == 'connectivity') then
+       read (fileUnit,610,END=660) line
+       chunkPos = IO_stringPos(line)
+       if(any(Marc_matNumber==IO_intValue(line,chunkPos,6_pInt))) then
+         do 
+           read (fileUnit,610,END=660) line
+           chunkPos = IO_stringPos(line)
+           tmp = IO_lc(IO_stringValue(line,chunkPos,1_pInt))
+           if (verify(trim(tmp),"0123456789")/=0) then                                                          ! found keyword
+             exit
+           else
+             contInts(1) = contInts(1) + 1_pInt  
+             read (tmp,*) contInts(contInts(1)+1)     
+           endif
+         enddo
+       endif  
+     endif
+   endif    
+ enddo    
+660 do i = 1_pInt,contInts(1)
       cpElem = cpElem+1_pInt
       mesh_mapFEtoCPelem(1,cpElem) = contInts(1_pInt+i)
       mesh_mapFEtoCPelem(2,cpElem) = cpElem
     enddo
-   endif
- enddo
  
-660 call math_qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt))            ! should be mesh_NcpElems
+call math_qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt))                  ! should be mesh_NcpElems
 
 end subroutine mesh_marc_map_elements
 
@@ -2936,7 +3053,6 @@ subroutine mesh_build_sharedElems
 
  allocate(node_seen(maxval(FE_NmatchingNodes)))
 
- 
  node_count = 0_pInt
 
  do e = 1_pInt,mesh_NcpElems

src/numerics.f90

@@ -25,7 +25,7 @@ module numerics
    nState                     = 10_pInt, &                                                          !< state loop limit
    nStress                    = 40_pInt, &                                                          !< stress loop limit
    pert_method                =  1_pInt, &                                                          !< method used in perturbation technique for tangent
-   fixedSeed                  =  0_pInt, &                                                          !< fixed seeding for pseudo-random number generator, Default 0: use random seed
+   randomSeed                 =  0_pInt, &                                                          !< fixed seeding for pseudo-random number generator, Default 0: use random seed
    worldrank                  =  0_pInt, &                                                          !< MPI worldrank (/=0 for MPI simulations only)
    worldsize                  =  0_pInt                                                             !< MPI worldsize (/=0 for MPI simulations only)
  integer(4), protected, public :: &
@@ -359,8 +359,8 @@ subroutine numerics_init
 
 !--------------------------------------------------------------------------------------------------
 ! random seeding parameter
-       case ('fixed_seed')
-         fixedSeed = IO_intValue(line,chunkPos,2_pInt)
+       case ('random_seed','fixed_seed')
+         randomSeed = IO_intValue(line,chunkPos,2_pInt)
 
 !--------------------------------------------------------------------------------------------------
 ! gradient parameter
@@ -560,9 +560,9 @@ subroutine numerics_init
 
 !--------------------------------------------------------------------------------------------------
 ! Random seeding parameter
- write(6,'(a24,1x,i16,/)')    ' fixed_seed:            ',fixedSeed
- if (fixedSeed <= 0_pInt) &
-   write(6,'(a,/)') ' No fixed Seed: Random is random!'
+ write(6,'(a24,1x,i16,/)')    ' random_seed:            ',randomSeed
+ if (randomSeed <= 0_pInt) &
+   write(6,'(a,/)')           ' random seed will be generated!'
 
 !--------------------------------------------------------------------------------------------------
 ! gradient parameter

src/plastic_disloUCLA.f90

@@ -1178,7 +1178,7 @@ end subroutine plastic_disloUCLA_dotState
 function plastic_disloUCLA_postResults(Tstar_v,Temperature,ipc,ip,el)
  use prec, only: &
    tol_math_check, &
-   dEq
+   dEq, dNeq0
  use math, only: &
    pi
  use material, only: &
@@ -1445,9 +1445,13 @@ function plastic_disloUCLA_postResults(Tstar_v,Temperature,ipc,ip,el)
            index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph))                                 ! at which index starts my family
            slipSystems2: do i = 1_pInt,plastic_disloUCLA_Nslip(f,instance)
               j = j + 1_pInt
+              if (dNeq0(abs(dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph))))) then
               plastic_disloUCLA_postResults(c+j) = &
                 (3.0_pReal*lattice_mu(ph)*plastic_disloUCLA_burgersPerSlipSystem(j,instance))/&
                 (16.0_pReal*pi*abs(dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph))))
+              else
+              plastic_disloUCLA_postResults(c+j) = huge(1.0_pReal)
+              endif
               plastic_disloUCLA_postResults(c+j)=min(plastic_disloUCLA_postResults(c+j),&
                                                             state(instance)%mfp_slip(j,of))
         enddo slipSystems2; enddo slipFamilies2

src/plastic_dislotwin.f90

@@ -1029,7 +1029,7 @@ subroutine plastic_dislotwin_init(fileUnit)
            do p = 1_pInt,3_pInt; do q = 1_pInt,3_pInt; do r = 1_pInt,3_pInt; do s = 1_pInt,3_pInt
              plastic_dislotwin_Ctwin3333(l,m,n,o,index_myFamily+j,instance) = &
              plastic_dislotwin_Ctwin3333(l,m,n,o,index_myFamily+j,instance) + &
-               lattice_C3333(p,q,r,s,instance) * &
+               lattice_C3333(p,q,r,s,phase) * &
                lattice_Qtwin(l,p,index_otherFamily+j,phase) * &
                lattice_Qtwin(m,q,index_otherFamily+j,phase) * &
                lattice_Qtwin(n,r,index_otherFamily+j,phase) * &
@@ -1087,7 +1087,7 @@ subroutine plastic_dislotwin_init(fileUnit)
            do p = 1_pInt,3_pInt; do q = 1_pInt,3_pInt; do r = 1_pInt,3_pInt; do s = 1_pInt,3_pInt
              plastic_dislotwin_Ctrans3333(l,m,n,o,index_myFamily+j,instance) = &
              plastic_dislotwin_Ctrans3333(l,m,n,o,index_myFamily+j,instance) + &
-               lattice_trans_C3333(p,q,r,s,instance) * &
+               lattice_trans_C3333(p,q,r,s,phase) * &
                lattice_Qtrans(l,p,index_otherFamily+j,phase) * &
                lattice_Qtrans(m,q,index_otherFamily+j,phase) * &
                lattice_Qtrans(n,r,index_otherFamily+j,phase) * &

src/plastic_nonlocal.f90

@@ -1823,14 +1823,14 @@ plasticState(ph)%state(iRhoF(1:ns,instance),of) = rhoForest
 plasticState(ph)%state(iTauF(1:ns,instance),of) = tauThreshold
 plasticState(ph)%state(iTauB(1:ns,instance),of) = tauBack
 
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
       .and. ((debug_e == el .and. debug_i == ip)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
     write(6,'(/,a,i8,1x,i2,1x,i1,/)') '<< CONST >> nonlocal_microstructure at el ip ',el,ip
     write(6,'(a,/,12x,12(e10.3,1x))') '<< CONST >> rhoForest', rhoForest
-    write(6,'(a,/,12x,12(f10.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold/1e6
-    write(6,'(a,/,12x,12(f10.5,1x),/)') '<< CONST >> tauBack / MPa', tauBack/1e6
+    write(6,'(a,/,12x,12(f10.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold*1e-6
+    write(6,'(a,/,12x,12(f10.5,1x),/)') '<< CONST >> tauBack / MPa', tauBack*1e-6
   endif
 #endif
 
@@ -1978,15 +1978,15 @@ if (Temperature > 0.0_pReal) then
 endif
     
 
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
       .and. ((debug_e == el .and. debug_i == ip)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
     write(6,'(/,a,i8,1x,i2,1x,i1,/)') '<< CONST >> nonlocal_kinetics at el ip',el,ip
-    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold / 1e6_pReal
-    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tau / MPa', tau / 1e6_pReal
-    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauNS / MPa', tauNS / 1e6_pReal
-    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> v / 1e-3m/s', v * 1e3
+    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold * 1e-6_pReal
+    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tau / MPa', tau * 1e-6_pReal
+    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauNS / MPa', tauNS * 1e-6_pReal
+    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> v / mm/s', v * 1e3
     write(6,'(a,/,12x,12(e12.5,1x))') '<< CONST >> dv_dtau', dv_dtau
     write(6,'(a,/,12x,12(e12.5,1x))') '<< CONST >> dv_dtauNS', dv_dtauNS
   endif
@@ -2176,12 +2176,12 @@ enddo
 dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
 
 
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
       .and. ((debug_e == el .and. debug_i == ip)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
     write(6,'(/,a,i8,1x,i2,1x,i1,/)') '<< CONST >> nonlocal_LpandItsTangent at el ip',el,ip
-    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> gdot total / 1e-3',gdotTotal*1e3_pReal
+    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> gdot total',gdotTotal
     write(6,'(a,/,3(12x,3(f12.7,1x),/))') '<< CONST >> Lp',transpose(Lp)
   endif
 #endif
@@ -2248,7 +2248,7 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,e
                                             dUpperOld, &                  ! old maximum stable dipole distance for edges and screws
                                             deltaDUpper                   ! change in maximum stable dipole distance for edges and screws
 
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt &
       .and. ((debug_e == el .and. debug_i == ip)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) &
@@ -2361,7 +2361,7 @@ forall (s = 1:ns, c = 1_pInt:2_pInt) &
   plasticState(ph)%deltaState(iRhoD(s,c,instance),of) = deltaRho(s,c+8_pInt)
 
 
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
       .and. ((debug_e == el .and. debug_i == ip)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
@@ -2522,11 +2522,11 @@ logical                                     considerEnteringFlux, &
 
 
 
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt &
       .and. ((debug_e == el .and. debug_i == ip)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) &
-    write(6,'(/,a,i8,1x,i2,1x,i1,/)') '<< CONST >> nonlocal_dotState at el ip ',el,ip
+    write(6,'(/,a,i8,1x,i2,/)') '<< CONST >> nonlocal_dotState at el ip ',el,ip
 #endif
 
 ph = material_phase(1_pInt,ip,el)
@@ -2589,7 +2589,7 @@ endif
 forall (t = 1_pInt:4_pInt) &
   gdot(1_pInt:ns,t) = rhoSgl(1_pInt:ns,t) * burgers(1:ns,instance) * v(1:ns,t)
 
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt &
       .and. ((debug_e == el .and. debug_i == ip)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
@@ -2663,7 +2663,7 @@ else
           / burgers(s,instance) * sqrt(rhoForest(s)) / lambda0(s,instance)
       endif
     enddo
-#ifndef _OPENMP
+#ifdef DEBUG
     if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
         .and. ((debug_e == el .and. debug_i == ip)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) &
@@ -2690,7 +2690,7 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then
   if (any( abs(gdot) > 0.0_pReal &                                                                  ! any active slip system ...
           .and. CFLfactor(instance) * abs(v) * timestep &
               > mesh_ipVolume(ip,el) / maxval(mesh_ipArea(:,ip,el)))) then                          ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
-#ifndef _OPENMP
+#ifdef DEBUG
     if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt) then 
       write(6,'(a,i5,a,i2)') '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
       write(6,'(a,e10.3,a,e10.3)') '<< CONST >> velocity is at  ', &
@@ -2952,7 +2952,7 @@ if (numerics_integrationMode == 1_pInt) then
 endif
 
 
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
       .and. ((debug_e == el .and. debug_i == ip .and. debug_g == 1_pInt)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
@@ -2978,7 +2978,7 @@ endif
 
 if (    any(rhoSglOriginal(1:ns,1:4) + rhoDot(1:ns,1:4) * timestep < -aTolRho(instance)) &
    .or. any(rhoDipOriginal(1:ns,1:2) + rhoDot(1:ns,9:10) * timestep < -aTolRho(instance))) then
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt) then 
     write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip
     write(6,'(a)') '<< CONST >> enforcing cutback !!!'

src/prec.f90

@@ -137,6 +137,7 @@ end subroutine prec_init
 !> @brief equality comparison for float with double precision
 ! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq
 ! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
+! AlmostEqualRelative
 !--------------------------------------------------------------------------------------------------
 logical elemental pure function dEq(a,b,tol)
 
@@ -153,6 +154,7 @@ end function dEq
 !> @brief inequality comparison for float with double precision
 ! replaces "!=" but for certain (relative) tolerance. Counterpart to dEq
 ! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
+! AlmostEqualRelative NOT
 !--------------------------------------------------------------------------------------------------
 logical elemental pure function dNeq(a,b,tol)
 
@@ -167,33 +169,35 @@ end function dNeq
 
 !--------------------------------------------------------------------------------------------------
 !> @brief equality to 0 comparison for float with double precision
-! replaces "==0" but for certain (absolute) tolerance. Counterpart to dNeq0
-! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
+! replaces "==0" but everything not representable as a normal number is treated as 0. Counterpart to dNeq0
+! https://de.mathworks.com/help/matlab/ref/realmin.html
+! https://docs.oracle.com/cd/E19957-01/806-3568/ncg_math.html
 !--------------------------------------------------------------------------------------------------
 logical elemental pure function dEq0(a,tol)
 
  implicit none
  real(pReal), intent(in)           :: a
  real(pReal), intent(in), optional :: tol
- real(pReal), parameter            :: eps = 2.220446049250313E-16                                   ! DBL_EPSILON in C
+ real(pReal), parameter            :: eps = 2.2250738585072014E-308                                 ! smallest non-denormalized number
 
- dEq0 = merge(.True., .False.,abs(a) <= merge(tol,eps,present(tol))*10.0_pReal)
+ dEq0 = merge(.True., .False.,abs(a) <= merge(tol,eps,present(tol)))
 end function dEq0
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief inequality to 0 comparison for float with double precision
-! replaces "!=0" but for certain (absolute) tolerance. Counterpart to dEq0
-! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
+! replaces "!=0" but everything not representable as a normal number is treated as 0. Counterpart to dEq0
+! https://de.mathworks.com/help/matlab/ref/realmin.html
+! https://docs.oracle.com/cd/E19957-01/806-3568/ncg_math.html
 !--------------------------------------------------------------------------------------------------
 logical elemental pure function dNeq0(a,tol)
 
  implicit none
  real(pReal), intent(in)           :: a
  real(pReal), intent(in), optional :: tol
- real(pReal), parameter            :: eps = 2.220446049250313E-16                                   ! DBL_EPSILON in C
+ real(pReal), parameter            :: eps = 2.2250738585072014E-308                                 ! smallest non-denormalized number
 
- dNeq0 = merge(.False., .True.,abs(a) <= merge(tol,eps,present(tol))*10.0_pReal)
+ dNeq0 = merge(.False., .True.,abs(a) <= merge(tol,eps,present(tol)))
 end function dNeq0
 
 

src/spectral_mech_Basic.f90

@@ -310,7 +310,6 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
    debug_spectral, &
    debug_spectralRotation
  use spectral_utilities, only: &
-   wgt, &
    tensorField_real, &
    utilities_FFTtensorForward, &
    utilities_fourierGammaConvolution, &

src/spectral_utilities.f90

@@ -971,6 +971,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
  real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF
  real(pReal)   :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet
 
+ external :: &
+  MPI_Allreduce
+
  write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
  flush(6)
  age = .False.