Merge pull request #28 from LEoPart-project/nate/transfer-fixes

Add particle deficiency test

Author: nate-sime

demo/demo_rayleigh-taylor-instability/demo_rayleigh-taylor-instability.py

@@ -157,6 +157,10 @@ def velocity(t):
     pyleopart.rk(mesh._cpp_object, ptcls, tableau,
                  velocity, t, dt)
 
+    deficient_cells = pyleopart.find_deficient_cells(phi._cpp_object, ptcls)
+    if len(deficient_cells) > 0:
+        pprint(f"Particle deficient cells found: {deficient_cells}",
+               rank=mesh.comm.rank)
     pyleopart.transfer_to_function(phi._cpp_object, ptcls, ptcls.field("phi"))
     ptcl_file.write_particles(ptcls, t, ["phi"])
     phi_file.write_function(phi, t=t)

leopart/cpp/transfer.h

@@ -49,6 +49,50 @@ void transfer_to_particles(const Particles<T>& particles, Field<T>& field,
   f->eval(x, xshape, p2c, u, ushape);
 }
 
+/// Transfer the provided particle field data to the finite element
+/// function using local l_2 projection.
+/// We solve the problem: find \f(u_h \in V\f) such that
+///
+/// \f[
+///    u_h(x_p) v(x_p) = u_p v(x_p) \quad \forall v \in V, \; p = 1,\ldots,n_p.
+/// \f]
+///
+/// Here \f(u_p\f) is the \f(p\f)th particle's data, \f(x_p\f) is the \f(p\f)th
+/// particle's position, \f(n_p\f) is the total number of particles
+/// and \f(V\f) is the function space to which the provided finite element
+/// function belongs.
+///
+/// @tparam T The function scalar type
+/// @tparam U The function geometry type
+/// @param f The finite element function
+/// @param pax The particles collection
+/// @param field The field data to be transferred
+template <dolfinx::scalar T, std::floating_point U>
+std::vector<std::int32_t> find_deficient_cells(
+  std::shared_ptr<dolfinx::fem::Function<T>> f,
+  const Particles<T>& pax)
+{
+  std::shared_ptr<const dolfinx::mesh::Mesh<U>> mesh
+    = f->function_space()->mesh();
+  std::shared_ptr<const dolfinx::fem::FunctionSpace<T>> V
+    = f->function_space();
+  const int tdim = mesh->topology()->dim();
+  std::int32_t ncells = mesh->topology()->index_map(tdim)->size_local();
+  std::shared_ptr<const dolfinx::fem::FiniteElement<T>> element
+    = f->function_space()->element();
+
+  const std::vector<std::vector<std::size_t>>& c2p
+    = pax.cell_to_particle();
+
+  std::vector<std::int32_t> deficient_cells;
+  for (std::int32_t c = 0; c < ncells; ++c)
+  {
+    const std::size_t ncdof = V->dofmap()->cell_dofs(c).size();
+    if (c2p[c].size() < ncdof)
+      deficient_cells.push_back(c);
+  }
+  return deficient_cells;
+}
 
 /// Transfer the provided particle field data to the finite element
 /// function using local l_2 projection. The solve_callback function
@@ -185,6 +229,7 @@ template <dolfinx::scalar T, std::floating_point U>
 void transfer_to_function(std::shared_ptr<dolfinx::fem::Function<T>> f,
                           const Particles<T>& pax, const Field<T>& field)
 {
+  // Simply solve the particle mass matrix / rhs system
   std::function<const std::vector<T>(mdspan_t<T, 2>, mdspan_t<T, 2>)>
    solve_function = [](mdspan_t<T, 2> QT_Q, mdspan_t<T, 2> QT_L)
   {
@@ -238,6 +283,7 @@ void transfer_to_function_constrained(
     ci0[i + space_dimension] = u;
   }
 
+  // Solve the  mass matrix / rhs optimisation problem with quadprog
   std::function<std::vector<T>(mdspan_t<T, 2>, mdspan_t<T, 2>)>
     solve_function = [&CE, &ce0, &CI, &ci0](
       mdspan_t<T, 2> QT_Q, mdspan_t<T, 2> QT_L)

leopart/cpp/wrapper.cpp

@@ -168,6 +168,11 @@ PYBIND11_MODULE(cpp, m)
         }, py::return_value_policy::move);
 
   // Transfer functions
+  m.def("find_deficient_cells",
+        [](std::shared_ptr<dolfinx::fem::Function<dtype>> f,
+        const Particles<dtype>& pax){
+          return leopart::transfer::find_deficient_cells<dtype, dtype_geom>(f, pax);
+        });
   m.def("transfer_to_particles", &leopart::transfer::transfer_to_particles<dtype>);
   m.def("transfer_to_function",
         [](std::shared_ptr<dolfinx::fem::Function<dtype>> f,

test/test_interpolate_project.py

@@ -24,6 +24,43 @@ def create_mesh(cell_type, dtype, n):
                              cell_type=cell_type, dtype=dtype)
 
 
+@pytest.mark.parametrize("k", [1, 2])
+@pytest.mark.parametrize("dtype", [np.float64])
+@pytest.mark.parametrize("cell_type", [dolfinx.mesh.CellType.triangle,
+                                       dolfinx.mesh.CellType.tetrahedron,
+                                       dolfinx.mesh.CellType.quadrilateral,
+                                       dolfinx.mesh.CellType.hexahedron])
+@pytest.mark.parametrize("shape", [(1,), (2,)])
+def test_find_deficient_cells(k, dtype, cell_type, shape):
+    mesh = create_mesh(cell_type, dtype, 4)
+    Q = dolfinx.fem.FunctionSpace(mesh, ("DG", k, shape))
+
+    npart = Q.dofmap.dof_layout.num_entity_closure_dofs(mesh.topology.dim)
+    x, c = pyleopart.mesh_fill(mesh._cpp_object, npart, seed=1)
+    if mesh.geometry.dim == 2:
+        x = np.c_[x, np.zeros_like(x[:, 0])]
+    ptcls = pyleopart.Particles(x, c)
+
+    num_cells = mesh.topology.index_map(mesh.topology.dim).size_local
+    cells_to_cull = [0, num_cells - 1]
+    assert cells_to_cull[0] != cells_to_cull[1]
+
+    num_to_delete = 1
+    for c in cells_to_cull:
+        pidxs = ptcls.cell_to_particle()[c]
+        for _ in range(num_to_delete):
+            ptcls.delete_particle(c, 0)
+        assert len(ptcls.cell_to_particle()[c]) == len(pidxs) - num_to_delete
+
+    u = dolfinx.fem.Function(Q)
+    deficient_cells = pyleopart.find_deficient_cells(u._cpp_object, ptcls)
+
+    cells_to_cull = np.sort(np.array(cells_to_cull, dtype=np.int32))
+    deficient_cells = np.sort(np.array(deficient_cells, dtype=np.int32))
+
+    assert np.all(cells_to_cull == deficient_cells)
+
+
 @pytest.mark.parametrize("k", [1, 2])
 @pytest.mark.parametrize("dtype", [np.float64])
 @pytest.mark.parametrize("cell_type", [dolfinx.mesh.CellType.triangle,