Bugfix: writevtk with empty parts

NEWS.md

@@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Fixed
 
 - Fixed bug in `num_cells` in the case where a `DistributedTriangulation` contained ghost cells. Since PR[#160](https://github.com/gridap/GridapDistributed.jl/pull/160).
+- Fixed bug in writevtk when dealing with empty processors. Since PR[#158](https://github.com/gridap/GridapDistributed.jl/pull/158).
 
 ## [0.4.4] 2024-08-14
 

src/Visualization.jl

@@ -154,7 +154,11 @@ function Visualization.write_vtk_file(
     parts,grid,filebase;celldata=celldata,nodaldata=nodaldata,
     compress=compress,append=append,ascii=ascii,vtkversion=vtkversion
   )
-  map(vtk_save,pvtk)
+  map(pvtk) do pvtk
+    if !isnothing(pvtk)
+      vtk_save(pvtk)
+    end
+  end
 end
 
 function Visualization.create_vtk_file(
@@ -164,15 +168,27 @@ function Visualization.create_vtk_file(
   celldata, nodaldata,
   compress=false,append=true,ascii=false,vtkversion=:default
 )
-  nparts = length(parts)
-  map(parts,grid,celldata,nodaldata) do part,g,c,n
-    Visualization.create_pvtk_file(
-      g,filebase;
-      part=part,nparts=nparts,
-      celldata=c,nodaldata=n,
-      compress=compress,append=append,ascii=ascii,vtkversion=vtkversion
-    )
+  nparts, new_parts = filter_empty_parts(parts,grid)
+  map(new_parts,grid,celldata,nodaldata) do part,g,c,n
+    if part > 0
+      Visualization.create_pvtk_file(
+        g,filebase;
+        part=part,nparts=nparts,
+        celldata=c,nodaldata=n,
+        compress=compress,append=append,ascii=ascii,vtkversion=vtkversion
+      )
+    end
+  end
+end
+
+function filter_empty_parts(parts,grid)
+  notempty = map(g -> Int(num_cells(g) > 0), grid)
+  new_parts = scan(+,notempty,type=:inclusive,init=0)
+  new_parts = map(new_parts,notempty) do part, e
+    iszero(e) ? -1 : part
   end
+  nparts = reduce(+,notempty)
+  return nparts, new_parts
 end
 
 const DistributedModelOrTriangulation = Union{DistributedDiscreteModel,DistributedTriangulation}
@@ -181,7 +197,7 @@ function Visualization.writevtk(
   arg::DistributedModelOrTriangulation,args...;
   compress=false,append=true,ascii=false,vtkversion=:default,kwargs...
 )
-  parts=get_parts(arg)
+  parts = get_parts(arg)
   map(visualization_data(arg,args...;kwargs...)) do visdata
     write_vtk_file(
       parts,visdata.grid,visdata.filebase,celldata=visdata.celldata,nodaldata=visdata.nodaldata,
@@ -195,7 +211,7 @@ function Visualization.createvtk(
   compress=false,append=true,ascii=false,vtkversion=:default,kwargs...
 )
   v = visualization_data(arg,args...;kwargs...)
-  parts=get_parts(arg)
+  parts = get_parts(arg)
   @notimplementedif length(v) != 1
   visdata = first(v)
   Visualization.create_vtk_file(

test/TestApp/src/TestApp.jl

@@ -17,4 +17,5 @@ module TestApp
   include("../../AdaptivityUnstructuredTests.jl")
   include("../../AdaptivityMultiFieldTests.jl")
   include("../../BlockSparseMatrixAssemblersTests.jl")
+  include("../../VisualizationTests.jl")
 end

test/VisualizationTests.jl

@@ -0,0 +1,31 @@
+module VisualizationTests
+
+using Gridap, GridapDistributed, PartitionedArrays
+
+function half_empty_trian(ranks,model)
+  cell_ids = get_cell_gids(model)
+  trians = map(ranks,local_views(model),partition(cell_ids)) do rank, model, ids
+    cell_mask = zeros(Bool, num_cells(model))
+    if rank ∈ (3,4)
+      cell_mask[own_to_local(ids)] .= true
+    end
+    Triangulation(model,cell_mask)
+  end
+  GridapDistributed.DistributedTriangulation(trians,model)
+end
+
+function main(distribute,parts)
+  ranks = distribute(LinearIndices((prod(parts),)))
+
+  model = CartesianDiscreteModel(ranks,parts,(0,1,0,1),(8,8))
+  V = FESpace(model, ReferenceFE(lagrangian,Float64,1))
+  uh = interpolate(x -> x[1]+x[2], V)
+
+  t1 = Triangulation(model)
+  writevtk(t1,"output/t1",cellfields=["uh" => uh]) 
+
+  t2 = half_empty_trian(ranks,model)
+  writevtk(t2,"output/t2",cellfields=["uh" => uh])
+end
+
+end # module

test/mpi/runtests_np4_body.jl

@@ -54,5 +54,10 @@ function all_tests(distribute,parts)
   TestApp.BlockSparseMatrixAssemblersTests.main(distribute,parts)
   PArrays.toc!(t,"BlockSparseMatrixAssemblers")
 
+  if prod(parts) == 4
+    TestApp.VisualizationTests.main(distribute,parts)
+    PArrays.toc!(t,"Visualization")
+  end
+
   display(t)
 end