From 718880f17722ffc48577548d3d12f1164ef1d685 Mon Sep 17 00:00:00 2001 From: Abisheik Panneerselvam Date: Mon, 31 May 2021 14:09:11 +0200 Subject: [PATCH 1/4] Initial nodal and ip coordinates added to the HDF5 output file --- src/mesh/discretization_mesh.f90 | 40 ++++++++++++++++++++++++++++---- 1 file changed, 36 insertions(+), 4 deletions(-) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 4d8546b6a..349e66ff2 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/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,ierr) + CHKERRQ(ierr) + allocate(mesh_node0_temp(dimPlex*mesh_Nnodes)) + call VecGetArrayF90(coords, 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,14 @@ 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]) + deallocate(mesh_node0_temp) call discretization_init(materialAt,& reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & mesh_node0) - call results_openJobFile - call results_closeGroup(results_addGroup('geometry')) - call results_closeJobFile + call writeGeometry(reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]),mesh_node0) end subroutine discretization_mesh_init @@ -231,4 +241,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 From 650fbd75099abcdd63916bbbad0f549947292f7d Mon Sep 17 00:00:00 2001 From: Abisheik Panneerselvam Date: Mon, 31 May 2021 15:01:42 +0200 Subject: [PATCH 2/4] Initial nodal and ip coordinates added to the HDF5 output file --- src/mesh/discretization_mesh.f90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 349e66ff2..af3021be0 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -144,10 +144,10 @@ subroutine discretization_mesh_init(restart) CHKERRQ(ierr) ! Get initial nodal coordinates - call DMGetCoordinates(geomMesh,coords,ierr) + call DMGetCoordinates(geomMesh,coords_node0,ierr) CHKERRQ(ierr) allocate(mesh_node0_temp(dimPlex*mesh_Nnodes)) - call VecGetArrayF90(coords, mesh_node0_temp,ierr) + call VecGetArrayF90(coords_node0, mesh_node0_temp,ierr) CHKERRQ(ierr) mesh_maxNips = FEM_nQuadrature(dimPlex,integrationOrder) @@ -167,7 +167,6 @@ subroutine discretization_mesh_init(restart) allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) mesh_node0 = reshape(mesh_node0_temp,[dimPlex,mesh_Nnodes]) - deallocate(mesh_node0_temp) call discretization_init(materialAt,& reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & From 91015fd07c2406d0c6c11d2f3e5f8e67817a1bdf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 31 May 2021 19:41:22 +0200 Subject: [PATCH 3/4] write out undeformed configuration --- src/mesh/DAMASK_mesh.f90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 979845484..ad0ec9794 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/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 From bd10ee033e5aeb644aaaa749c259cf33241832f7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 1 Jun 2021 06:49:14 +0200 Subject: [PATCH 4/4] dealing with user errors/incomplete files --- python/damask/_result.py | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index 8dd915d1b..1d927ac00 100644 --- a/python/damask/_result.py +++ b/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 @@ -111,6 +110,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']) @@ -727,8 +728,8 @@ class Result: ---------- T_sym : str Name of symmetric tensor dataset. - eigenvalue : str, optional - Eigenvalue. Select from 'max', 'mid', 'min'. Defaults to 'max'. + eigenvalue : {'max', 'mid', 'min'} + Eigenvalue. Defaults to 'max'. """ self._add_generic_pointwise(self._add_eigenvalue,{'T_sym':T_sym},{'eigenvalue':eigenvalue}) @@ -760,9 +761,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}) @@ -771,14 +772,8 @@ class Result: @staticmethod def _add_IPF_color(l,q): m = util.scale_to_coprime(np.array(l)) - try: - lattice = {'fcc':'cF','bcc':'cI','hex':'hP'}[q['meta']['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) + lattice = q['meta']['lattice'] + o = Orientation(rotation = q['data'],lattice=lattice) return { 'data': np.uint8(o.IPF_color(l)*255), @@ -901,7 +896,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), @@ -1407,6 +1402,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