nodal displacements for grid solver

currently for non-MPI only

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)

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:

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
 

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