Mixed-dimension refinement test that hits an assert in debug mode

test_mixed_dimension_refinement.cc

// libMesh includes
#include "libmesh/libmesh.h"
#include "libmesh/elem.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/replicated_mesh.h"
#include "libmesh/mesh_refinement.h"

using namespace libMesh;

int main (int argc, char** argv)
{
  LibMeshInit init(argc, argv);

  // This code should eventually be made to work for both Replicated
  // and Distributed meshes, but for now I'm just focusing on
  // ReplicatedMesh.
  ReplicatedMesh mesh(init.comm());

  // Build 3x3 square mesh of QUAD4 elements. Note that the same error
  // can also be triggered using QUADSHELL4 elements.
  MeshTools::Generation::build_square (mesh,
                                       /*nx*/3,
                                       /*ny*/3,
                                       /*xmin*/0., /*xmax*/1.,
                                       /*ymin*/0., /*ymax*/1.,
                                       QUAD4);

  // Get pointer to middle Elem
  Elem * middle_elem = mesh.elem_ptr(4);

  // Construct bottom side of middle Elem (for connectivity only)
  std::unique_ptr<Elem> bottom_side = middle_elem->build_side_ptr(0);

  // Attach Edge3 to vertices of bottom_side
  Elem * beam_elem = mesh.add_elem(Elem::build(EDGE3));

  // Set vertices of beam_elem to match bottom_side
  for (unsigned int i=0; i<2; ++i)
    beam_elem->set_node(i) = bottom_side->node_ptr(i);

  // Add node at midpoint of bottom_side to mesh
  Node * midpoint_node = mesh.add_point(bottom_side->vertex_average());

  // Connect midpoint node
  beam_elem->set_node(2) = midpoint_node;

  // New element has different geometric type; put it in a different
  // subdomain so that we can use the Exodus writer with it.
  beam_elem->subdomain_id() = 1;

  // prepare_for_use(): includes calling find_neighbors() and detect_interior_parents()
  mesh.prepare_for_use();

  // Check that the beam_elem has an interior parent set by
  // prepare_for_use() Note: when there is more than one candidate for
  // an interior parent (i.e. for an internal side) the candidate with
  // the lowest elem id (here, 1) takes precedence.
  if (beam_elem->interior_parent())
    libMesh::out << "beam_elem has interior parent " << beam_elem->interior_parent()->id() << std::endl;

  // Now simulate the refinement pattern that we do in the real application
  MeshRefinement mesh_refinement(mesh);

  unsigned int max_refinement_level = 2;
  for (unsigned int r=0; r<max_refinement_level; ++r)
  {
    // Flag active 1D elements for refinement
    for (auto & elem : mesh.active_element_ptr_range())
    {
      if (elem->dim() == 1)
        elem->set_refinement_flag(Elem::REFINE);
      else
        elem->set_refinement_flag(Elem::DO_NOTHING);
    }
    mesh_refinement.refine_elements();
  }

  // Write out the refined mesh
  mesh.write("test_mixed_dimension_refinement.exo");

  return 0;
}