Merge branch 'misc-improvements' into 'development'

polishing

See merge request damask/DAMASK!497

.gitlab-ci.yml

@@ -45,7 +45,7 @@ variables:
   MPI_INTEL:      "MPI/Intel/2022.0.1/IntelMPI/2021.5.0"
   # ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   PETSC_GNU:       "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1"
-  PETSC_INTELLLVM: "Libraries/PETSc/3.16.2/oneAPI-2022.0.1-IntelMPI-2021.5.0"
+  PETSC_INTELLLVM: "Libraries/PETSc/3.16.3/oneAPI-2022.0.1-IntelMPI-2021.5.0"
   PETSC_INTEL:    "Libraries/PETSc/3.16.2/Intel-2022.0.1-IntelMPI-2021.5.0"
   # ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   MSC:            "FEM/MSC/2021.3.1"

CMakeLists.txt

@@ -10,7 +10,7 @@ endif()
 # Dummy project to determine compiler names and version
 project(Prerequisites LANGUAGES)
 set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig")
-pkg_check_modules(PETSC REQUIRED PETSc>=3.12.0 PETSc<3.17.0)
+pkg_search_module(PETSC REQUIRED PETSc>=3.12.0 PETSc<3.17.0)
 pkg_get_variable(CMAKE_Fortran_COMPILER PETSc fcompiler)
 pkg_get_variable(CMAKE_C_COMPILER PETSc ccompiler)
 
@@ -88,16 +88,12 @@ else()
     message(FATAL_ERROR "Compiler type(CMAKE_Fortran_COMPILER_ID) not recognized")
 endif()
 
-file(STRINGS "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/petsc/conf/petscvariables" PETSC_EXTERNAL_LIB REGEX "PETSC_WITH_EXTERNAL_LIB = .*$?")
+file(STRINGS "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/petsc/conf/petscvariables" PETSC_EXTERNAL_LIB REGEX "PETSC_EXTERNAL_LIB_BASIC = .*$?")
-string(REGEX MATCHALL "-[lLW]([^\" ]+)" PETSC_EXTERNAL_LIB "${PETSC_EXTERNAL_LIB}")
+string(REPLACE "PETSC_EXTERNAL_LIB_BASIC = " "" PETSC_EXTERNAL_LIB "${PETSC_EXTERNAL_LIB}")
-list(REMOVE_DUPLICATES PETSC_EXTERNAL_LIB)
-string(REPLACE ";" " " PETSC_EXTERNAL_LIB "${PETSC_EXTERNAL_LIB}")
 message("PETSC_EXTERNAL_LIB:\n${PETSC_EXTERNAL_LIB}\n")
 
 file(STRINGS "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/petsc/conf/petscvariables" PETSC_INCLUDES REGEX "PETSC_FC_INCLUDES = .*$?")
-string(REGEX MATCHALL "-I([^\" ]+)" PETSC_INCLUDES "${PETSC_INCLUDES}")
+string(REPLACE "PETSC_FC_INCLUDES = " "" PETSC_INCLUDES "${PETSC_INCLUDES}")
-list(REMOVE_DUPLICATES PETSC_INCLUDES)
-string(REPLACE ";" " " PETSC_INCLUDES "${PETSC_INCLUDES}")
 message("PETSC_INCLUDES:\n${PETSC_INCLUDES}\n")
 
 set(CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE} "${BUILDCMD_PRE} ${OPENMP_FLAGS} ${STANDARD_CHECK} ${OPTIMIZATION_FLAGS} ${COMPILE_FLAGS} ${PRECISION_FLAGS}")
@@ -109,7 +105,7 @@ if(CMAKE_BUILD_TYPE STREQUAL "DEBUG")
 endif()
 
 set(CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}   "${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}} ${PETSC_INCLUDES} ${BUILDCMD_POST}")
-set(CMAKE_Fortran_LINK_EXECUTABLE             "${CMAKE_Fortran_LINK_EXECUTABLE} <OBJECTS> -o <TARGET> <LINK_LIBRARIES> ${PETSC_EXTERNAL_LIB} -lz ${BUILDCMD_POST}")
+set(CMAKE_Fortran_LINK_EXECUTABLE             "${CMAKE_Fortran_LINK_EXECUTABLE} <OBJECTS> -o <TARGET> <LINK_LIBRARIES> -L${PETSC_LIBRARY_DIRS} -lpetsc ${PETSC_EXTERNAL_LIB} -lz ${BUILDCMD_POST}")
 
 message("Fortran Compiler Flags:\n${CMAKE_Fortran_FLAGS_${CMAKE_BUILD_TYPE}}\n")
 message("C Compiler Flags:\n${CMAKE_C_FLAGS_${CMAKE_BUILD_TYPE}}\n")

Makefile

@@ -10,14 +10,12 @@ all: grid mesh
 .PHONY: grid
 grid:
 	@cmake -B build/grid -DDAMASK_SOLVER=grid -DCMAKE_INSTALL_PREFIX=${PWD} -DCMAKE_BUILD_TYPE=${BUILD_TYPE} -DBUILDCMD_POST=${BUILDCMD_POST} -DBUILDCMD_PRE=${BUILDCMD_PRE} -DOPTIMIZATION=${OPTIMIZATION} -DOPENMP=${OPENMP}
-	@cmake --build build/grid --parallel
+	@cmake --build build/grid --parallel --target install
-	@cmake --install build/grid
 
 .PHONY: mesh
 mesh:
 	@cmake -B build/mesh -DDAMASK_SOLVER=mesh -DCMAKE_INSTALL_PREFIX=${PWD} -DCMAKE_BUILD_TYPE=${BUILD_TYPE} -DBUILDCMD_POST=${BUILDCMD_POST} -DBUILDCMD_PRE=${BUILDCMD_PRE} -DOPTIMIZATION=${OPTIMIZATION} -DOPENMP=${OPENMP}
-	@cmake --build build/mesh --parallel
+	@cmake --build build/mesh --parallel --target install
-	@cmake --install build/mesh
 
 .PHONY: clean
 clean:

install/MarcMentat/apply_DAMASK_modifications.py

@@ -67,9 +67,7 @@ os.system(f'xvfb-run -a {executable} -compile {menu_file}')
 
 print('setting file access rights...')
 
-files = (glob.glob(str(marc_root/f'marc{marc_version}/tools/*_damask*')) +
+for file in (glob.glob(str(marc_root/f'marc{marc_version}/tools/*_damask*')) +
-         glob.glob(str(marc_root/f'mentat{marc_version}/bin/kill[4-6]')) +
+             glob.glob(str(marc_root/f'mentat{marc_version}/bin/kill[4-6]')) +
-         glob.glob(str(marc_root/f'mentat{marc_version}/bin/submit[4-6]')))
+             glob.glob(str(marc_root/f'mentat{marc_version}/bin/submit[4-6]'))):
-
-for file in files:
     os.chmod(file , 0o755)

processing/legacy/addDisplacement.py

@@ -1,71 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(usage='%prog options [ASCIItable(s)]', description = """
-Add displacments resulting from deformation gradient field.
-Operates on periodic three-dimensional x,y,z-ordered data sets.
-Outputs at cell centers or cell nodes (into separate file).
-
-""", version = scriptID)
-
-parser.add_option('-f',
-                  '--defgrad',
-                  dest    = 'f',
-                  metavar = 'string',
-                  help    = 'label of deformation gradient [%default]')
-parser.add_option('-p',
-                  '--pos', '--position',
-                  dest    = 'pos',
-                  metavar = 'string',
-                  help    = 'label of coordinates [%default]')
-parser.add_option('--nodal',
-                  dest    = 'nodal',
-                  action  = 'store_true',
-                  help    = 'output nodal (instead of cell-centered) displacements')
-
-parser.set_defaults(f   = 'f',
-                    pos = 'pos',
-                   )
-
-(options,filenames) = parser.parse_args()
-
-for name in filenames:
-    damask.util.report(scriptName,name)
-
-    table = damask.Table.load(StringIO(''.join(sys.stdin.read())) if name is None else name)
-    grid,size,origin = damask.grid_filters.cellsSizeOrigin_coordinates0_point(table.get(options.pos))
-
-    F = table.get(options.f).reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+(3,3))
-    if options.nodal:
-        damask.Table(damask.grid_filters.coordinates0_node(grid,size).reshape(-1,3,order='F'),
-                     {'pos':(3,)})\
-              .add('avg({}).{}'.format(options.f,options.pos),
-                   damask.grid_filters.displacement_avg_node(size,F).reshape(-1,3,order='F'),
-                   scriptID+' '+' '.join(sys.argv[1:]))\
-              .add('fluct({}).{}'.format(options.f,options.pos),
-                   damask.grid_filters.displacement_fluct_node(size,F).reshape(-1,3,order='F'),
-                    scriptID+' '+' '.join(sys.argv[1:]))\
-               .save((sys.stdout if name is None else os.path.splitext(name)[0]+'_nodal.txt'))
-    else:
-        table.add('avg({}).{}'.format(options.f,options.pos),
-                  damask.grid_filters.displacement_avg_point(size,F).reshape(-1,3,order='F'),
-                  scriptID+' '+' '.join(sys.argv[1:]))\
-             .add('fluct({}).{}'.format(options.f,options.pos),
-                  damask.grid_filters.displacement_fluct_point(size,F).reshape(-1,3,order='F'),
-                  scriptID+' '+' '.join(sys.argv[1:]))\
-             .save((sys.stdout if name is None else name))

python/damask/_crystal.py

@@ -114,12 +114,13 @@ class Crystal():
 
     def __repr__(self):
         """Represent."""
-        return '\n'.join([f'Crystal family {self.family}']
+        family = f'Crystal family: {self.family}'
-                       + ([] if self.lattice is None else [f'Bravais lattice {self.lattice}']+
+        return family if self.lattice is None else \
-                                                           list(map(lambda x:f'{x[0]}: {x[1]:.5g}',
+               '\n'.join([family,
-                                                                    zip(['a','b','c','α','β','γ',],
+                          f'Bravais lattice: {self.lattice}',
-                                                                         self.parameters))))
+                          'a={:.5g}m, b={:.5g}m, c={:.5g}m'.format(*self.parameters[:3]),
-                        )
+                          'α={:.5g}°, β={:.5g}°, γ={:.5g}°'.format(*np.degrees(self.parameters[3:]))])
+
 
     def __eq__(self,other):
         """
@@ -378,7 +379,7 @@ class Crystal():
         """
         _kinematics = {
             'cF': {
-                'slip' :[np.array([
+                'slip': [np.array([
                            [+0,+1,-1, +1,+1,+1],
                            [-1,+0,+1, +1,+1,+1],
                            [+1,-1,+0, +1,+1,+1],
@@ -398,7 +399,7 @@ class Crystal():
                            [+1,+0,-1, +1,+0,+1],
                            [+0,+1,+1, +0,+1,-1],
                            [+0,+1,-1, +0,+1,+1]])],
-                'twin' :[np.array([
+                'twin': [np.array([
                            [-2, 1, 1,  1, 1, 1],
                            [ 1,-2, 1,  1, 1, 1],
                            [ 1, 1,-2,  1, 1, 1],
@@ -413,7 +414,7 @@ class Crystal():
                            [-1, 1, 2, -1, 1,-1]])]
             },
             'cI': {
-                'slip' :[np.array([
+                'slip': [np.array([
                            [+1,-1,+1, +0,+1,+1],
                            [-1,-1,+1, +0,+1,+1],
                            [+1,+1,+1, +0,-1,+1],
@@ -464,7 +465,7 @@ class Crystal():
                            [+1,+1,+1, -3,+2,+1],
                            [+1,+1,-1, +3,-2,+1],
                            [+1,-1,+1, +3,+2,-1]])],
-                'twin' :[np.array([
+                'twin': [np.array([
                            [-1, 1, 1,  2, 1, 1],
                            [ 1, 1, 1, -2, 1, 1],
                            [ 1, 1,-1,  2,-1, 1],
@@ -479,7 +480,7 @@ class Crystal():
                            [ 1, 1, 1,  1, 1,-2]])]
             },
             'hP': {
-                'slip' :[np.array([
+                'slip': [np.array([
                            [+2,-1,-1,+0, +0,+0,+0,+1],
                            [-1,+2,-1,+0, +0,+0,+0,+1],
                            [-1,-1,+2,+0, +0,+0,+0,+1]]),
@@ -514,7 +515,7 @@ class Crystal():
                            [+1,+1,-2,+3, -1,-1,+2,+2],
                            [-1,+2,-1,+3, +1,-2,+1,+2],
                            [-2,+1,+1,+3, +2,-1,-1,+2]])],
-                'twin' :[np.array([
+                'twin': [np.array([
                            [-1, 0, 1, 1,  1, 0,-1, 2],   # shear = (3-(c/a)^2)/(sqrt(3) c/a) <-10.1>{10.2}
                            [ 0,-1, 1, 1,  0, 1,-1, 2],
                            [ 1,-1, 0, 1, -1, 1, 0, 2],
@@ -542,7 +543,74 @@ class Crystal():
                            [-1,-1, 2,-3, -1,-1, 2, 2],
                            [ 1,-2, 1,-3,  1,-2, 1, 2],
                            [ 2,-1,-1,-3,  2,-1,-1, 2]])]
-                },
+            },
+            'tI': {
+                'slip': [np.array([
+                           [+0,+0,+1, +1,+0,+0],
+                           [+0,+0,+1, +0,+1,+0]]),
+                         np.array([
+                           [+0,+0,+1, +1,+1,+0],
+                           [+0,+0,+1, -1,+1,+0]]),
+                         np.array([
+                           [+0,+1,+0, +1,+0,+0],
+                           [+1,+0,+0, +0,+1,+0]]),
+                         np.array([
+                           [+1,-1,+1, +1,+1,+0],
+                           [+1,-1,-1, +1,+1,+0],
+                           [-1,-1,-1, -1,+1,+0],
+                           [-1,-1,+1, -1,+1,+0]]),
+                         np.array([
+                           [+1,-1,+0, +1,+1,+0],
+                           [+1,+1,+0, +1,-1,+0]]),
+                         np.array([
+                           [+0,+1,+1, +1,+0,+0],
+                           [+0,-1,+1, +1,+0,+0],
+                           [-1,+0,+1, +0,+1,+0],
+                           [+1,+0,+1, +0,+1,+0]]),
+                         np.array([
+                           [+0,+1,+0, +0,+0,+1],
+                           [+1,+0,+0, +0,+0,+1]]),
+                         np.array([
+                           [+1,+1,+0, +0,+0,+1],
+                           [-1,+1,+0, +0,+0,+1]]),
+                         np.array([
+                           [+0,+1,-1, +0,+1,+1],
+                           [+0,-1,-1, +0,-1,+1],
+                           [-1,+0,-1, -1,+0,+1],
+                           [+1,+0,-1, +1,+0,+1]]),
+                         np.array([
+                           [+1,-1,+1, +0,+1,+1],
+                           [+1,+1,-1, +0,+1,+1],
+                           [+1,+1,+1, +0,+1,-1],
+                           [-1,+1,+1, +0,+1,-1],
+                           [+1,-1,-1, +1,+0,+1],
+                           [-1,-1,+1, +1,+0,+1],
+                           [+1,+1,+1, +1,+0,-1],
+                           [+1,-1,+1, +1,+0,-1]]),
+                         np.array([
+                           [+1,+0,+0, +0,+1,+1],
+                           [+1,+0,+0, +0,+1,-1],
+                           [+0,+1,+0, +1,+0,+1],
+                           [+0,+1,+0, +1,+0,-1]]),
+                         np.array([
+                           [+0,+1,-1, +2,+1,+1],
+                           [+0,-1,-1, +2,-1,+1],
+                           [+1,+0,-1, +1,+2,+1],
+                           [-1,+0,-1, -1,+2,+1],
+                           [+0,+1,-1, -2,+1,+1],
+                           [+0,-1,-1, -2,-1,+1],
+                           [-1,+0,-1, -1,-2,+1],
+                           [+1,+0,-1, +1,-2,+1]]),
+                         np.array([
+                           [-1,+1,+1, +2,+1,+1],
+                           [-1,-1,+1, +2,-1,+1],
+                           [+1,-1,+1, +1,+2,+1],
+                           [-1,-1,+1, -1,+2,+1],
+                           [+1,+1,+1, -2,+1,+1],
+                           [+1,-1,+1, -2,-1,+1],
+                           [-1,+1,+1, -1,-2,+1],
+                           [+1,+1,+1, +1,-2,+1]])]
+                }
         }
         master = _kinematics[self.lattice][mode]
         if self.lattice == 'hP':

python/damask/_orientation.py

@@ -514,6 +514,17 @@ class Orientation(Rotation,Crystal):
          [ 0.07359167 -0.36505797  0.92807163]]
         Bunge Eulers / deg: (11.40, 21.86, 0.60)
 
+        Plot a sample from the Mackenzie distribution.
+
+        >>> import matplotlib.pyplot as plt
+        >>> import damask
+        >>> N = 10000
+        >>> a = damask.Orientation.from_random(shape=N,family='cubic')
+        >>> b = damask.Orientation.from_random(shape=N,family='cubic')
+        >>> d = a.disorientation(b).as_axis_angle(degrees=True,pair=True)[1]
+        >>> plt.hist(d,25)
+        >>> plt.show()
+
         """
         if self.family != other.family:
             raise NotImplementedError('disorientation between different crystal families')

python/tests/test_Crystal.py

@@ -79,3 +79,23 @@ class TestCrystal:
                     a=a,b=b,c=c,
                     alpha=alpha,beta=beta,gamma=gamma)
         assert np.allclose(points,c.lattice_points)
+
+    @pytest.mark.parametrize('crystal,length',
+                            [(Crystal(lattice='cF'),[12,6]),
+                             (Crystal(lattice='cI'),[12,12,24]),
+                             (Crystal(lattice='hP'),[3,3,6,12,6]),
+                             (Crystal(lattice='tI',c=1.2),[2,2,2,4,2,4,2,2,4,8,4,8,8])
+                            ])
+    def test_N_slip(self,crystal,length):
+        assert [len(s) for s in crystal.kinematics('slip')['direction']] == length
+        assert [len(s) for s in crystal.kinematics('slip')['plane']] == length
+
+    @pytest.mark.parametrize('crystal,length',
+                            [(Crystal(lattice='cF'),[12]),
+                             (Crystal(lattice='cI'),[12]),
+                             (Crystal(lattice='hP'),[6,6,6,6]),
+                            ])
+    def test_N_twin(self,crystal,length):
+        assert [len(s) for s in crystal.kinematics('twin')['direction']] == length
+        assert [len(s) for s in crystal.kinematics('twin')['plane']] == length
+

python/tests/test_grid_filters.py

@@ -2,6 +2,8 @@ import pytest
 import numpy as np
 
 from damask import grid_filters
+from damask import Grid
+from damask import seeds
 
 class TestGridFilters:
 
@@ -139,12 +141,19 @@ class TestGridFilters:
         else:
             function(unordered,mode)
 
-    def test_regrid(self):
+    def test_regrid_identity(self):
          size = np.random.random(3)
          cells = np.random.randint(8,32,(3))
-         F    = np.broadcast_to(np.eye(3), tuple(cells)+(3,3))
+         F = np.broadcast_to(np.eye(3), tuple(cells)+(3,3))
          assert all(grid_filters.regrid(size,F,cells) == np.arange(cells.prod()))
 
+    def test_regrid_double_cells(self):
+         size = np.random.random(3)
+         cells = np.random.randint(8,32,(3))
+         g = Grid.from_Voronoi_tessellation(cells,size,seeds.from_random(size,10))
+         F = np.broadcast_to(np.eye(3), tuple(cells)+(3,3))
+         assert all(g.scale(cells*2).material.flatten() ==
+                    g.material.flatten()[grid_filters.regrid(size,F,cells*2)])
 
     @pytest.mark.parametrize('differential_operator',[grid_filters.curl,
                                                       grid_filters.divergence,

src/phase_mechanical_elastic.f90

@@ -30,7 +30,7 @@ module subroutine elastic_init(phases)
     phase, &
     mech, &
     elastic
-  logical :: thermal_active
+
 
   print'(/,1x,a)', '<<<+-  phase:mechanical:elastic init  -+>>>'
   print'(/,1x,a)', '<<<+-  phase:mechanical:elastic:Hooke init  -+>>>'

src/rotations.f90

@@ -372,7 +372,7 @@ end function rotTensor4
 
 
 !---------------------------------------------------------------------------------------------------
-!> @brief Rotate a rank-4 tensor in Voigt 6x6 notation passively (default) or actively.
+!> @brief Rotate a rank-4 stiffness tensor in Voigt 6x6 notation passively (default) or actively.
 !> @details: https://scicomp.stackexchange.com/questions/35600
 !! ToDo: Need to check active/passive !!!
 !---------------------------------------------------------------------------------------------------
@@ -393,11 +393,11 @@ pure function rotStiffness(self,C,active) result(cRot)
     R = self%asMatrix()
   endif
 
-  M = reshape([R(1,1)**2.0_pReal,           R(2,1)**2.0_pReal,           R(3,1)**2.0_pReal, &
+  M = reshape([R(1,1)**2,                   R(2,1)**2,                   R(3,1)**2, &
                R(2,1)*R(3,1),               R(1,1)*R(3,1),               R(1,1)*R(2,1), &
-               R(1,2)**2.0_pReal,           R(2,2)**2.0_pReal,           R(3,2)**2.0_pReal, &
+               R(1,2)**2,                   R(2,2)**2,                   R(3,2)**2, &
                R(2,2)*R(3,2),               R(1,2)*R(3,2),               R(1,2)*R(2,2), &
-               R(1,3)**2.0_pReal,           R(2,3)**2.0_pReal,           R(3,3)**2.0_pReal, &
+               R(1,3)**2,                   R(2,3)**2,                   R(3,3)**2, &
                R(2,3)*R(3,3),               R(1,3)*R(3,3),               R(1,3)*R(2,3), &
                2.0_pReal*R(1,2)*R(1,3),     2.0_pReal*R(2,2)*R(2,3),     2.0_pReal*R(3,2)*R(3,3), &
                R(2,2)*R(3,3)+R(2,3)*R(3,2), R(1,2)*R(3,3)+R(1,3)*R(3,2), R(1,2)*R(2,3)+R(1,3)*R(2,2), &