Merge branch 'HDFgeometryupdate' into 'development'

HDF geometry update

See merge request damask/DAMASK!395

python/damask/_result.py

@@ -14,7 +14,6 @@ from collections.abc import Iterable
 import h5py
 import numpy as np
 import numpy.ma as ma
-from numpy.lib import recfunctions as rfn
 
 import damask
 from . import VTK
@@ -115,6 +114,8 @@ class Result:
             r=re.compile('increment_[0-9]+')
             self.increments = sorted([i for i in f.keys() if r.match(i)],key=util.natural_sort)
             self.times      = [round(f[i].attrs['t/s'],12) for i in self.increments]
+            if len(self.increments) == 0:
+                raise ValueError('incomplete DADF5 file')
 
             self.N_materialpoints, self.N_constituents = np.shape(f['cell_to/phase'])
 
@@ -733,8 +734,8 @@ class Result:
         ----------
         T_sym : str
             Name of symmetric tensor dataset.
-        eigenvalue : str, optional
+        eigenvalue : {'max', 'mid', 'min'}
-            Eigenvalue. Select from 'max', 'mid', 'min'. Defaults to 'max'.
+            Eigenvalue. Defaults to 'max'.
 
         """
         self._add_generic_pointwise(self._add_eigenvalue,{'T_sym':T_sym},{'eigenvalue':eigenvalue})
@@ -766,9 +767,9 @@ class Result:
         ----------
         T_sym : str
             Name of symmetric tensor dataset.
-        eigenvalue : str, optional
+        eigenvalue : {'max', 'mid', 'min'}
             Eigenvalue to which the eigenvector corresponds.
-            Select from 'max', 'mid', 'min'. Defaults to 'max'.
+            Defaults to 'max'.
 
         """
         self._add_generic_pointwise(self._add_eigenvector,{'T_sym':T_sym},{'eigenvalue':eigenvalue})
@@ -777,14 +778,8 @@ class Result:
     @staticmethod
     def _add_IPF_color(l,q):
         m = util.scale_to_coprime(np.array(l))
-        try:
+        lattice =  q['meta']['lattice']
-            lattice = {'fcc':'cF','bcc':'cI','hex':'hP'}[q['meta']['lattice']]
+        o = Orientation(rotation = q['data'],lattice=lattice)
-        except KeyError:
-            lattice =  q['meta']['lattice']
-        try:
-            o = Orientation(rotation = (rfn.structured_to_unstructured(q['data'])),lattice=lattice)
-        except ValueError:
-            o = Orientation(rotation = q['data'],lattice=lattice)
 
         return {
                 'data': np.uint8(o.IPF_color(l)*255),
@@ -907,7 +902,7 @@ class Result:
             t = 'tensor'
             if o is None: o = 'fro'
         else:
-            raise ValueError
+            raise ValueError(f'invalid norm order {ord}')
 
         return {
                 'data':  np.linalg.norm(x['data'],ord=o,axis=axis,keepdims=True),
@@ -1552,6 +1547,9 @@ class Result:
             v = self.geometry0
         elif mode.lower()=='point':
             v = VTK.from_poly_data(self.coordinates0_point)
+        else:
+            raise ValueError(f'invalid mode {mode}')
+
         v.set_comments(util.execution_stamp('Result','save_VTK'))
 
         N_digits = int(np.floor(np.log10(max(1,int(self.increments[-1][10:])))))+1

src/mesh/DAMASK_mesh.f90

@@ -252,6 +252,10 @@ program DAMASK_mesh
     write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded'                  ! statistics file
   endif
 
+  print'(/,a)', ' ... writing initial configuration to file ........................'
+  flush(IO_STDOUT)
+  call CPFEM_results(0,0.0_pReal)
+
   loadCaseLooping: do currentLoadCase = 1, size(loadCases)
     time0 = time                                                                                    ! load case start time
     guess = loadCases(currentLoadCase)%followFormerTrajectory                                       ! change of load case? homogeneous guess for the first inc

src/mesh/discretization_mesh.f90

@@ -41,6 +41,9 @@ module discretization_mesh
     mesh_ipVolume, &                                                                                !< volume associated with IP (initially!)
     mesh_node0                                                                                      !< node x,y,z coordinates (initially!)
 
+  real(pReal), pointer, dimension(:) :: &
+    mesh_node0_temp
+  
   real(pReal), dimension(:,:,:), allocatable :: &
     mesh_ipCoordinates                                                                              !< IP x,y,z coordinates (after deformation!)
 
@@ -82,7 +85,7 @@ subroutine discretization_mesh_init(restart)
   class(tNode), pointer :: &
     num_mesh
   integer :: integrationOrder                                                                       !< order of quadrature rule required
-
+  type(tvec)  :: coords_node0
 
   print'(/,a)',   ' <<<+-  discretization_mesh init  -+>>>'
 
@@ -140,6 +143,13 @@ subroutine discretization_mesh_init(restart)
   call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr)
   CHKERRQ(ierr)
 
+! Get initial nodal coordinates
+  call DMGetCoordinates(geomMesh,coords_node0,ierr)
+  CHKERRQ(ierr)
+  allocate(mesh_node0_temp(dimPlex*mesh_Nnodes))
+  call VecGetArrayF90(coords_node0, mesh_node0_temp,ierr)
+  CHKERRQ(ierr)
+  
   mesh_maxNips = FEM_nQuadrature(dimPlex,integrationOrder)
 
   call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p)
@@ -156,14 +166,13 @@ subroutine discretization_mesh_init(restart)
   if (debug_ip < 1 .or. debug_ip > mesh_maxNips)            call IO_error(602,ext_msg='IP')
 
   allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal)
+  mesh_node0 = reshape(mesh_node0_temp,[dimPlex,mesh_Nnodes])
 
   call discretization_init(materialAt,&
                            reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), &
                            mesh_node0)
 
-  call results_openJobFile
+  call writeGeometry(reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]),mesh_node0)  
-  call results_closeGroup(results_addGroup('geometry'))
-  call results_closeJobFile
 
 end subroutine discretization_mesh_init
 
@@ -231,4 +240,26 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
 
 end subroutine mesh_FEM_build_ipCoordinates
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Write all information needed for the DADF5 geometry
+!--------------------------------------------------------------------------------------------------
+subroutine writeGeometry(coordinates_points,coordinates_nodes)
+
+  real(pReal), dimension(:,:), intent(in) :: &
+  coordinates_nodes, &
+  coordinates_points
+  
+  call results_openJobFile
+  call results_closeGroup(results_addGroup('geometry'))
+  
+  call results_writeDataset('geometry',coordinates_nodes,'x_n', &
+        'initial coordinates of the nodes','m')
+  
+  call results_writeDataset('geometry',coordinates_points,'x_p', &
+        'initial coordinates of the materialpoints (cell centers)','m')
+  
+  call results_closeJobFile
+  
+  end subroutine writeGeometry
+
 end module discretization_mesh