how to delete an element on one MPI processus?

[fdrui]

Hello everyone,
I am trying to modify a distributed mesh (that I have built manually and that is already "prepared for use") by deleting one cell on one specific MPI proc (MPI proc 0 in my test). Thus I used the delete_elem() method on MPI rank 0, but then, my test crashes when calling prepare_for_use().
Has anyone an idea about this problem and how to delete (or create) cells properly?
Thank you in advance!
F.D.

[roystgnr]

If you're manually modifying a distributed mesh, you'll have to keep it in sync manually too. Deleting an element would have to be done on every rank where it's ghosted, not just on the owner (and definitely not just on rank 0 if that's not the owner).

How to best do that is a different question. Assuming you don't have any a priori way to tell which elements need deletion on every processor, you may be stuck doing a pull_parallel_vector_data, querying all ghost ids and returning which ones are in need of deletion. This would need to be preceded by a push_parallel_vector_data if there's a chance that even an element's owner might not be the rank that can tell it should be deleted.

[fdrui]

Thank you for your answer. However, in my test the cell I want to delete is only accessible on proc 0. Here is the test:

#include "libmesh/elem.h"
#include "libmesh/mesh.h"
#include "libmesh/distributed_mesh.h"
#include "libmesh/dof_map.h"
#include "libmesh/gmsh_io.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_communication.h"
#include <iostream>
#include <stdio.h>
#include <type_traits>

int main(int argc, char **argv)
{
  // Initialize libMesh.
  LibMeshInit init (argc, argv);
  
  int MPI_rank;
  MPI_Comm_rank (libMesh::GLOBAL_COMM_WORLD, &MPI_rank);

  // create a mesh from scratch
  libMesh::DistributedMesh mesh0(init.comm());
  
  // 4x4 mesh distributed over procs
  // 20 -- 21 -- 22 -- 23 -- 24
  //  |  12 |     |     | 15  |
  // 15 -- 16 -- 17 -- 18 -- 19
  //  |  8  |     |     | 11  |
  // 10 -- 11 -- 12 -- 13 -- 14
  //  |  4  |     |     |  7  |
  // 05 -- 06 -- 07 -- 08 -- 09
  //  |  0  |  1  |  2  |  3  |
  //  0 -- 01 -- 02 -- 03 -- 04

  // add nodes
  int node_sum = 0;
  for (int j=0; j<3; j++){
    double yy = 2*(MPI_rank/2)+j;
    for (int i=0; i<3; i++){
      double xx = 2*(MPI_rank%2)+i;
      unsigned int node_id = 5*floor(yy)+floor(xx); //3*j+i;
      Node * newnode = mesh0.add_point (Point(xx, yy, 0), node_id); //4*node_id+MPI_rank);
      newnode->set_unique_id() = node_id;
      node_sum++;
    }
  }

  int maxnode = 25;

  int elem_sum = 0;
  for (int j=0; j<2; j++){
    for (int i=0; i<2; i++){
      unsigned int elem_id = 4*(2*(MPI_rank/2)+j) + 2*(MPI_rank%2)+i;
      unsigned int first_node_id = elem_id+2*(MPI_rank/2)+j;
      auto elem = Elem::build_with_id(QUAD4, elem_id); 
      elem->processor_id(MPI_rank);
      elem->set_node(0) = mesh0.node_ptr(first_node_id); 
      elem->set_node(1) = mesh0.node_ptr(first_node_id+1); 
      elem->set_node(2) = mesh0.node_ptr(first_node_id+6);
      elem->set_node(3) = mesh0.node_ptr(first_node_id+5); 
      elem->set_unique_id() = maxnode + elem_id;
      mesh0.add_elem(std::move(elem));
      elem_sum++;
    }
  }

  MeshCommunication().gather_neighboring_elements(cast_ref<DistributedMesh &>(mesh0));
  mesh0.prepare_for_use (true);

  MPI_Barrier(libMesh::GLOBAL_COMM_WORLD);

  // mesh0.set_isnt_prepared();
  if (MPI_rank == 0){
  auto delelem = mesh0.elem_ptr(0);
  mesh0.delete_elem(delelem);
  }

  mesh0.prepare_for_use();
  return 0;
}

[roystgnr]

However, in my test the cell I want to delete is only accessible on proc 0

I wouldn't be so confident, without checking that in the code. When you call that first prepare_for_use() you haven't disabled repartitioning (so depending on ParMETIS's whims, elem 0 might now be owned by another proc and ghosted on any combination of non-owners), and you haven't disabled renumbering (so ids in general, albeit probably not 0 in particular, might now refer to different elements in your mesh than the ones you expect).

Also, I'm now guessing you're running an opt mode build of this? You want to be trying out new code in dbg mode for a test case this small, or devel mode at the worst.

[fdrui]

Hi !
I used to develop my code with a version of the master branch (of LibMesh) that dated from the end of year 2022. A few days ago, I updated to a recent state of master and my tests now fail due to an assert at the end of the libmesh_assert_valid_boundary_ids() function:

libmesh_assert(mesh.comm().semiverify
                     (elem ? &all_bcids : nullptr));

(l.1698 of mesh_tools.C), called by prepare for use in debug mode.
In my code, I am building a mesh in parallel, where each process creates its own cells. I don't define any boundary condition nor boundary element.
Then, for instance on proc 0, cells nr 0 exists, while it does not exist on proc 1.
On the one hand, on proc 0, although there is no boundary condition, all_bcids will be non null because of "separators" made of BoundaryInfo::invalid_id.
On the other hand, all_bcids will be null on proc 1, since elem does not exist.
I believe that this difference makes the assert raise an exception here, am I right? and in the previous version of libMesh I was using, this test was not present here, instead, asserts were made at different stages sooner on the bcids vectors...

What could I do to avoid this problem?
Thank you very much for your help!

Best regards,

Florence

[fdrui]

Ok, unfortunately, I won't be able to check that before thursday... I will let you know there!
As I remember, it is the really first element (i=0) whatever its type (edge, quad, ...).

[roystgnr]

Also, if you can replicate the problem with a test case as compact as the code above, consider submitting a PR to add it to our unit test suite? Then if there's any proposed library changes that would cause a regression in behavior we'll find the regression as it happens, and either fix the bug in the changes (if that's what it is) or fix the test (if the test turned out to be relying on undefined behavior) before merging the changes. In the former cases you'd never see a regression in your own applications; in the latter you'd see a regression but you'd also see an example of how to fix it.

I personally never write new test coverage just because I like writing test coverage; sometimes I do it to defend myself against collaborators and coworkers screwing things up, and most of the time I do it to defend myself against future-Roy screwing things up.

[fdrui]

Hi!
What I can observe in a debugger is:

1. the arguments of semiverify as defined at l. 2107 of parallel_implementation.h are:

• proc 0: r non null pointer (to std::vector)
• other procs: r is a null pointer

2. at line 2107 of parallel_implementation.h:

• proc 0: psize is a non null pointer to an integer (whose value is 24 in my case)
• other procs: psize is a null pointer

3. in the semiverify function defined l. 2071 of parallel_implementation.h, at line2082:

• proc 0: tempmin = tempmax = 24
• other procs: tempmin = 18446744073709551615, tempmax = 0

4. strangely after line 2083 on all procs (after the mpi_allreduce): tempmin=18446744073709551615 (tempmax = 24)
5. consequently, the variable invalid is true on proc 0 first, and then on all other procs... making the assertions fail...

[fdrui]

My compiler is gcc 11.4.0, with OpenMPI v4.1.2

[roystgnr]

Oh! This may be an OpenMPI bug: open-mpi/ompi#10525

It should be fixed in v4.1.5.