diff --git a/src/vtkhdf_h5.F90.fypp b/src/vtkhdf_h5.F90.fypp
index 5194347..dd9b5c4 100644
--- a/src/vtkhdf_h5.F90.fypp
+++ b/src/vtkhdf_h5.F90.fypp
@@ -221,21 +221,24 @@ contains
   !! dataset itself is a collective operation.
 #:for s in num_specs
 
-  subroutine write_dataset_${s["kind"]}$(ctx, loc_id, name, array, root)
+  subroutine write_dataset_${s["kind"]}$(ctx, loc_id, name, array, root, force_leading_dims)
 
     type(vtkhdf_ctx), intent(in) :: ctx
     integer(hid_t), intent(in) :: loc_id
     character(*), intent(in) :: name
     ${s["decl"]}$, intent(in) :: array(..)
     integer, intent(in), optional :: root
+    integer(hsize_t), intent(in), optional :: force_leading_dims(:)
 
     integer(hid_t) :: type_id, dset_id, mem_space_id, data_space_id, dxpl_id
     integer :: ierr
+    ${s["decl"]}$ :: fill_value
 
     type_id = ${s["h5type"]}$
+    fill_value = 0
 
     call write_dataset_common(ctx, loc_id, name, type_id, shape(array,hsize_t), &
-        dset_id, mem_space_id, data_space_id, root)
+        dset_id, mem_space_id, data_space_id, root, force_leading_dims, fill_value)
 
     call ctx%get_dxpl(dxpl_id) ! collective data transfer in MPI build
 
@@ -270,7 +273,7 @@ contains
   !! in the subsequent collective H5Dwrite call.
 
   subroutine write_dataset_common(ctx, loc_id, name, type_id, dims, &
-      dset_id, mem_space_id, data_space_id, root)
+      dset_id, mem_space_id, data_space_id, root, force_leading_dims, fill_value)
 
     type(vtkhdf_ctx), intent(in) :: ctx
     integer(hid_t), intent(in) :: loc_id, type_id
@@ -278,11 +281,14 @@ contains
     integer(hsize_t), intent(in) :: dims(:)
     integer(hid_t), intent(out) :: dset_id, mem_space_id, data_space_id
     integer, intent(in), optional :: root
+    integer(hsize_t), intent(in), optional :: force_leading_dims(:)
+    type(*), intent(in), optional :: fill_value
 
     integer :: ierr, rank
     integer(hsize_t), allocatable :: mem_dims(:), data_dims(:), start(:), count(:)
     integer(hsize_t) :: n
     logical :: rank_writes
+    integer(hid_t) :: dcpl_id
 
     !! Convert DIMS to its equivalent C form (MEM_DIMS)
     if (size(dims) == 0) then ! scalar data
@@ -308,10 +314,26 @@ contains
     INSIST(n > 0) ! writing an empy dataset not allowed
     data_dims = mem_dims
     data_dims(1) = n
+
+    if (present(force_leading_dims)) then
+      INSIST(size(force_leading_dims) == size(dims)-1)
+      INSIST(all(force_leading_dims >= dims(1:size(dims)-1)))
+      data_dims(2:) = force_leading_dims(size(force_leading_dims):1:-1)
+    end if
+
     data_space_id = H5Screate(data_dims)
     INSIST(data_space_id > 0)
-    dset_id = H5Dcreate(loc_id, name, type_id, data_space_id)
+
+    dcpl_id = H5Pcreate(H5P_DATASET_CREATE)
+    if (present(force_leading_dims)) then
+      if (present(fill_value)) then
+        ierr = H5Pset_fill_value(dcpl_id, type_id, fill_value)
+        INSIST(ierr >= 0)
+      end if
+    end if
+    dset_id = H5Dcreate(loc_id, name, type_id, data_space_id, dcpl_id=dcpl_id)
     INSIST(ctx%global_all(dset_id >= 0))
+    ierr = H5Pclose(dcpl_id)
 
     !! Starting index for the dataset hyperslab for this rank
     allocate(start, mold=mem_dims)
diff --git a/src/vtkhdf_h5_c_binding.F90 b/src/vtkhdf_h5_c_binding.F90
index d75dff0..0a37aee 100644
--- a/src/vtkhdf_h5_c_binding.F90
+++ b/src/vtkhdf_h5_c_binding.F90
@@ -361,6 +361,13 @@ function H5Pset_dxpl_mpio(dxpl_id, xfer_mode) &
       integer(c_int) :: hdferr
     end function
 
+    function H5Pset_fill_value(plist_id, type_id, value) result(h5err) bind(c,name='H5Pset_fill_value')
+      import :: hid_t, c_int
+      integer(hid_t), value :: plist_id, type_id
+      type(*), intent(in) :: value(*)
+      integer(c_int) :: h5err
+    end function
+
     function H5Pclose(prp_id) result(h5err) bind(c,name='H5Pclose')
       import :: hid_t, c_int
       integer(hid_t), value :: prp_id
@@ -369,7 +376,8 @@ function H5Pclose(prp_id) result(h5err) bind(c,name='H5Pclose')
   end interface
 
   public :: H5Pcreate, H5Pclose, H5Pset_chunk, H5Pset_link_creation_order, &
-            H5Pset_all_coll_metadata_ops, H5Pset_coll_metadata_write, H5Pset_dxpl_mpio
+            H5Pset_all_coll_metadata_ops, H5Pset_coll_metadata_write, H5Pset_dxpl_mpio, &
+            H5Pset_fill_value
 
   interface ! to wrapper functions that take Fortran comm instead of C comm
     function H5Pset_fapl_mpio(fapl_id, comm) &
diff --git a/src/vtkhdf_ug_type.F90.fypp b/src/vtkhdf_ug_type.F90.fypp
index 37c0d1b..59faca7 100644
--- a/src/vtkhdf_ug_type.F90.fypp
+++ b/src/vtkhdf_ug_type.F90.fypp
@@ -276,7 +276,9 @@ contains
 #endif
 
     call h5_write_dataset(this%ctx, this%root_id, 'Types', types)
-    call h5_write_dataset(this%ctx, this%root_id, 'Points', x)
+
+    INSIST(size(x,1) <= 3)
+    call h5_write_dataset(this%ctx, this%root_id, 'Points', x, force_leading_dims=[3_hsize_t])
 
   end subroutine write_mesh_${r["kind"]}$_${i["kind"]}$
 
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index fb1339e..b831912 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -24,6 +24,15 @@ else()
   set_tests_properties(vtkhdf_mb_file_type PROPERTIES PROCESSORS 1)
 endif()
 
+add_executable(vtkhdf_2d_test vtkhdf_2d_test.F90)
+target_link_libraries(vtkhdf_2d_test PRIVATE fvtkhdf)
+if(ENABLE_MPI)
+  add_test(vtkhdf_2d_test ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 1 ./vtkhdf_2d_test ${MPIEXEC_POSTFLAGS})
+else()
+  add_test(vtkhdf_2d_test ./vtkhdf_2d_test)
+endif()
+set_tests_properties(vtkhdf_2d_test PROPERTIES PROCESSORS 1)
+
 add_executable(test_leak test_leak.F90)
 target_link_libraries(test_leak PRIVATE fvtkhdf)
 add_test(NAME test_leak COMMAND test_leak)
diff --git a/test/vtkhdf_2d_test.F90 b/test/vtkhdf_2d_test.F90
new file mode 100644
index 0000000..0e57d17
--- /dev/null
+++ b/test/vtkhdf_2d_test.F90
@@ -0,0 +1,53 @@
+program vtkhdf_2d_test
+
+  use,intrinsic :: iso_fortran_env, only: r8 => real64, int8
+  use vtkhdf_ug_file_type
+  use vtkhdf_vtk_cell_types
+#ifdef USE_MPI
+  use mpi_f08
+#endif
+  implicit none
+
+  real(r8), allocatable :: x(:,:), scalar_cell_data(:), vector_cell_data(:,:)
+  real(r8), allocatable :: scalar_point_data(:), vector_point_data(:,:)
+  integer, allocatable :: cnode(:), xcnode(:)
+  integer(int8), allocatable :: types(:)
+  character(:), allocatable :: errmsg
+  integer :: stat
+#ifdef USE_MPI
+  integer :: istat
+#endif
+
+  type(vtkhdf_ug_file) :: vizfile
+
+#ifdef USE_MPI
+  call MPI_Init(istat)
+  call vizfile%create('ug_2d_test.vtkhdf', MPI_COMM_WORLD, stat, errmsg)
+#else
+  call vizfile%create('ug_2d_test.vtkhdf', stat, errmsg)
+#endif
+  if (stat /= 0) error stop errmsg
+
+  ! A 2D mesh (two triangles forming a square)
+  x = reshape([0.0_r8,0.0_r8, 1.0_r8,0.0_r8, 1.0_r8,1.0_r8, 0.0_r8,1.0_r8], shape=[2,4])
+  cnode = [1,2,3, 1,3,4]
+  xcnode = [1,4,7]
+  types = [VTK_TRIANGLE, VTK_TRIANGLE]
+
+  call vizfile%write_mesh(x, cnode, xcnode, types)
+
+  !! Add some 2D point data to make sure it still works
+  allocate(scalar_point_data(4))
+  scalar_point_data = [1.0_r8, 2.0_r8, 3.0_r8, 4.0_r8]
+  call vizfile%write_point_data('scalar-data', scalar_point_data)
+
+  allocate(vector_point_data(2,4))
+  vector_point_data = x
+  call vizfile%write_point_data('vector-data', vector_point_data)
+
+  call vizfile%close
+#ifdef USE_MPI
+  call MPI_Finalize(istat)
+#endif
+
+end program