diff --git a/processing/post/DADF5_vtk_cells.py b/processing/post/DADF5_vtk_cells.py index 6d1e8e340..0b8fd07de 100755 --- a/processing/post/DADF5_vtk_cells.py +++ b/processing/post/DADF5_vtk_cells.py @@ -88,7 +88,7 @@ for filename in options.filenames: results.set_visible('constituents', False) results.set_visible('materialpoints',True) - for label in options.mat: + for label in options.mat: for p in results.iter_visible('mat_physics'): if p != 'generic': for m in results.iter_visible('materialpoints'): @@ -113,7 +113,13 @@ for filename in options.filenames: if results.structured: writer = vtk.vtkXMLRectilinearGridWriter() - + results.set_visible('constituents', False) + results.set_visible('materialpoints',False) + x = results.get_dataset_location('u_n') + vtk_data.append(numpy_support.numpy_to_vtk(num_array=results.read_dataset(x,0),deep=True,array_type=vtk.VTK_DOUBLE)) + vtk_data[-1].SetName('u') + rGrid.GetPointData().AddArray(vtk_data[-1]) + dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) if not os.path.isdir(dirname): os.mkdir(dirname,0o755) diff --git a/python/damask/dadf5.py b/python/damask/dadf5.py index 278a499e7..f7b065271 100644 --- a/python/damask/dadf5.py +++ b/python/damask/dadf5.py @@ -289,7 +289,13 @@ class DADF5(): """Returns the location of all active datasets with given label.""" path = [] with h5py.File(self.filename,'r') as f: - for i in self.iter_visible('increments'): + for i in self.iter_visible('increments'): + k = '/'.join([i,'geometry',label]) + try: + f[k] + path.append(k) + except KeyError as e: + print('unable to locate geometry dataset: {}'.format(str(e))) for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']): for oo in self.iter_visible(o): for pp in self.iter_visible(p): @@ -313,7 +319,11 @@ class DADF5(): if len(shape) == 1: shape = shape +(1,) dataset = np.full(shape,np.nan) for pa in path: - label = pa.split('/')[2] + label = pa.split('/')[2] + + if (pa.split('/')[1] == 'geometry'): + dataset = np.array(f[pa]) + continue p = np.where(f['mapping/cellResults/constituent'][:,c]['Name'] == str.encode(label))[0] if len(p)>0: diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index f19434415..6a2dcfe6f 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -24,7 +24,6 @@ module spectral_utilities include 'fftw3-mpi.f03' logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill - integer, public, parameter :: maxPhaseFields = 2 integer, public :: nActiveFields = 0 !-------------------------------------------------------------------------------------------------- @@ -1068,6 +1067,63 @@ subroutine utilities_updateCoords(F) enddo; enddo; enddo call discretization_setIPcoords(reshape(IPcoords,[3,grid(1)*grid(2)*grid3])) + call discretization_setNodeCoords(reshape(NodeCoords(IPcoords,Favg),[3,(grid(1)+1)*(grid(2)+1)*(grid3+1)])) + + contains + + function nodeCoords(IPcoords,Favg) + real(pReal), dimension(:,:,:,:), intent(in) :: IPcoords + real(pReal), dimension(3,3), intent(in) :: Favg + + real(pReal), dimension(3,size(IPcoords,2)+1,size(IPcoords,3)+1,size(IPcoords,3)+1) :: nodeCoords + real(pReal), dimension(3,size(IPcoords,2)+2,size(IPcoords,3)+2,size(IPcoords,3)+2) :: IPCoords_wrapped + + integer :: & + i,j,k,n + integer, dimension(3) :: & + shift, & + lookup, & + me + integer, dimension(3,8) :: & + neighbor = reshape([ & + 0, 0, 0, & + 1, 0, 0, & + 1, 1, 0, & + 0, 1, 0, & + 0, 0, 1, & + 1, 0, 1, & + 1, 1, 1, & + 0, 1, 1 ], [3,8]) + + do k = grid3offset,grid3offset+grid3+1 + do j = 0,grid(2)+1 + do i = 0,grid(1)+1 + if (k==0 .or. k==grid(3)+1 .or. & ! z skin + j==0 .or. j==grid(2)+1 .or. & ! y skin + i==0 .or. i==grid(1)+1 ) then ! x skin + me = [i,j,k] ! me on skin + shift = sign(abs(grid+1-2*me)/(grid+1),grid+1-2*me) + lookup = me-1+shift*grid + IPCoords_wrapped(1:3,i+1,j+1,k+1) = IPcoords(1:3,lookup(1)+1,lookup(2)+1,lookup(3)+1) & + - matmul(Favg, real(shift,pReal)*geomSize) + else + IPCoords_wrapped(1:3,i+1,j+1,k+1) = IPCoords(1:3,i,j,k) + endif + enddo; enddo; enddo + + + nodeCoords = 0.0_pReal + do k = 0,grid3; do j = 0,grid(2); do i = 0,grid(1) + do n = 1,8 + nodeCoords(1:3,i+1,j+1,k+1) = nodeCoords(1:3,i+1,j+1,k+1) & + + IPCoords_wrapped(1:3,i+1+neighbor(1,n), & + j+1+neighbor(2,n), & + k+1+neighbor(3,n) ) + enddo + enddo; enddo; enddo + nodeCoords = nodeCoords/8.0_pReal + + end function nodeCoords end subroutine utilities_updateCoords diff --git a/src/mesh_grid.f90 b/src/mesh_grid.f90 index 26e017aa5..a69c9b807 100644 --- a/src/mesh_grid.f90 +++ b/src/mesh_grid.f90 @@ -399,78 +399,4 @@ pure function IPneighborhood(grid) end function IPneighborhood -!!-------------------------------------------------------------------------------------------------- -!!> @brief builds mesh of (distorted) cubes for given coordinates (= center of the cubes) -!!-------------------------------------------------------------------------------------------------- -!function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes) -! -! real(pReal), intent(in), dimension(:,:,:,:) :: & -! centres -! real(pReal), dimension(3,size(centres,2)+1,size(centres,3)+1,size(centres,4)+1) :: & -! nodes -! real(pReal), intent(in), dimension(3) :: & -! gDim -! real(pReal), intent(in), dimension(3,3) :: & -! Favg -! real(pReal), dimension(3,size(centres,2)+2,size(centres,3)+2,size(centres,4)+2) :: & -! wrappedCentres -! -! integer :: & -! i,j,k,n -! integer, dimension(3), parameter :: & -! diag = 1 -! integer, dimension(3) :: & -! shift = 0, & -! lookup = 0, & -! me = 0, & -! iRes = 0 -! integer, dimension(3,8) :: & -! neighbor = reshape([ & -! 0, 0, 0, & -! 1, 0, 0, & -! 1, 1, 0, & -! 0, 1, 0, & -! 0, 0, 1, & -! 1, 0, 1, & -! 1, 1, 1, & -! 0, 1, 1 ], [3,8]) -! -!!-------------------------------------------------------------------------------------------------- -!! initializing variables -! iRes = [size(centres,2),size(centres,3),size(centres,4)] -! nodes = 0.0_pReal -! wrappedCentres = 0.0_pReal -! -!!-------------------------------------------------------------------------------------------------- -!! building wrappedCentres = centroids + ghosts -! wrappedCentres(1:3,2:iRes(1)+1,2:iRes(2)+1,2:iRes(3)+1) = centres -! do k = 0,iRes(3)+1 -! do j = 0,iRes(2)+1 -! do i = 0,iRes(1)+1 -! if (k==0 .or. k==iRes(3)+1 .or. & ! z skin -! j==0 .or. j==iRes(2)+1 .or. & ! y skin -! i==0 .or. i==iRes(1)+1 ) then ! x skin -! me = [i,j,k] ! me on skin -! shift = sign(abs(iRes+diag-2*me)/(iRes+diag),iRes+diag-2*me) -! lookup = me-diag+shift*iRes -! wrappedCentres(1:3,i+1, j+1, k+1) = & -! centres(1:3,lookup(1)+1,lookup(2)+1,lookup(3)+1) & -! - matmul(Favg, real(shift,pReal)*gDim) -! endif -! enddo; enddo; enddo -! -!!-------------------------------------------------------------------------------------------------- -!! averaging -! do k = 0,iRes(3); do j = 0,iRes(2); do i = 0,iRes(1) -! do n = 1,8 -! nodes(1:3,i+1,j+1,k+1) = & -! nodes(1:3,i+1,j+1,k+1) + wrappedCentres(1:3,i+1+neighbor(1,n), & -! j+1+neighbor(2,n), & -! k+1+neighbor(3,n) ) -! enddo -! enddo; enddo; enddo -! nodes = nodes/8.0_pReal -! -!end function mesh_nodesAroundCentres - end module mesh_grid