debugging mixed grids

examples/test_parallel_grid.rs

@@ -280,7 +280,7 @@ fn test_parallel_assembly_flat_triangle_grid<C: Communicator>(
     let rank = comm.rank();
     let size = comm.size();
 
-    let n = 5;
+    let n = 10;
     let grid = example_flat_triangle_grid(comm, n);
     let element = LagrangeElementFamily::<f64>::new(degree, cont);
     let space = ParallelFunctionSpace::new(&grid, &element);
@@ -460,7 +460,7 @@ fn test_parallel_assembly_mixed_grid<C: Communicator>(comm: &C, degree: usize, c
     let rank = comm.rank();
     let size = comm.size();
 
-    let n = 10;
+    let n = 5;
     let grid = example_mixed_grid(comm, n);
     let element = LagrangeElementFamily::<f64>::new(degree, cont);
     let space = ParallelFunctionSpace::new(&grid, &element);
@@ -522,12 +522,11 @@ fn test_parallel_assembly_mixed_grid<C: Communicator>(comm: &C, degree: usize, c
         println!("== full_matrix ==");
         println!("{:?}", full_matrix.indptr());
         println!("{:?}", full_matrix.indices());
-        println!("{:?}", full_matrix.data());
-
+        //println!("{:?}", full_matrix.data());
         println!("== serial_matrix ==");
         println!("{:?}", serial_matrix.indptr());
         println!("{:?}", serial_matrix.indices());
-        println!("{:?}", serial_matrix.data());
+        //println!("{:?}", serial_matrix.data());
 
         for (i, j) in full_matrix.indices().iter().zip(serial_matrix.indices()) {
             assert_eq!(i, j);
@@ -559,11 +558,8 @@ fn main() {
     let world = universe.world();
     let rank = world.rank();
 
-    test_parallel_assembly_flat_triangle_grid(&world, 1, Continuity::Continuous);
-
-    if rank < 100 {
-        return;
-    }
+    test_parallel_assembly_mixed_grid(&world, 2, Continuity::Continuous);
+    if rank < 100 { return; }
 
     if rank == 0 {
         println!("Testing FlatTriangleGrid in parallel.");
@@ -589,11 +585,11 @@ fn main() {
         }
         test_parallel_assembly_flat_triangle_grid(&world, degree, Continuity::Continuous);
         if rank == 0 {
-            println!("Testing assembly with DP{degree} using SingleElementGrid in parallel.");
+            println!("Testing assembly with P{degree} using SingleElementGrid in parallel.");
         }
         test_parallel_assembly_single_element_grid(&world, degree, Continuity::Continuous);
         if rank == 0 {
-            println!("Testing assembly with DP{degree} using MixedGrid in parallel.");
+            println!("Testing assembly with P{degree} using MixedGrid in parallel.");
         }
         test_parallel_assembly_mixed_grid(&world, degree, Continuity::Continuous);
     }

src/assembly/batched.rs

@@ -210,7 +210,7 @@ fn assemble_batch_singular<
                         .unwrap();
                     }
                 }
-                if trial_space.ownership(*trial_dof) == Ownership::Owned {
+                if trial_space.grid().cell_from_index(trial_space.grid().cell_index_from_id(*trial_cell)).ownership() == Ownership::Owned {
                     output.rows.push(test_space.global_dof_index(*test_dof));
                     output.cols.push(trial_space.global_dof_index(*trial_dof));
                     output.data.push(sum);
@@ -472,7 +472,7 @@ fn assemble_batch_singular_correction<
                         };
                     }
                 }
-                if trial_space.ownership(*trial_dof) == Ownership::Owned {
+                if trial_space.grid().cell_from_index(trial_space.grid().cell_index_from_id(*trial_cell)).ownership() == Ownership::Owned {
                     output.rows.push(test_space.global_dof_index(*test_dof));
                     output.cols.push(trial_space.global_dof_index(*trial_dof));
                     output.data.push(sum);

src/function/function_space/parallel.rs

@@ -86,6 +86,7 @@ impl<'a, T: RlstScalar, GridImpl: ParallelGridType + GridType<T = T::Real>>
         let (cell_dofs, entity_dofs, dofmap_size, owner_data) =
             assign_dofs(rank as usize, grid, e_family);
 
+
         let mut elements = HashMap::new();
         for cell in grid.cell_types() {
             elements.insert(*cell, e_family.element(*cell));
@@ -186,6 +187,8 @@ impl<'a, T: RlstScalar, GridImpl: ParallelGridType + GridType<T = T::Real>>
             }
         }
 
+        println!("{rank} {:?}", cell_dofs.iter().map(|cell| cell.iter().map(|i| global_dof_numbers[*i]).collect::<Vec<_>>()).collect::<Vec<_>>());
+
         let serial_space = SerialFunctionSpace {
             grid: grid.local_grid(),
             elements,

src/grid/mixed_grid/parallel.rs

@@ -280,6 +280,9 @@ where
             });
         }
 
+        println!("{rank} {:?}", cells.iter().map(|i| point_ids[*i]).collect::<Vec<_>>());
+        println!("{rank} {cell_ownership:?}");
+
         let serial_grid = MixedGrid::new(
             coordinates,
             cells,