connected_components.C

Go to the documentation of this file.

#include <libmesh/libmesh.h>
#include <libmesh/elem.h>
#include <libmesh/mesh_generation.h>
#include <libmesh/mesh_tools.h>
#include <libmesh/mesh.h>
#include <libmesh/node.h>
#include <libmesh/sparse_matrix.h>
#include <libmesh/system.h>
 
#include <timpi/parallel_implementation.h>
 
#include "test_comm.h"
#include "libmesh_cppunit.h"
 
 
using namespace libMesh;
 
class ConnectedComponentsTest : public CppUnit::TestCase
{
public:
  LIBMESH_CPPUNIT_TEST_SUITE( ConnectedComponentsTest );
 
#ifdef LIBMESH_HAVE_SOLVER
  CPPUNIT_TEST( testEdge2 );
  CPPUNIT_TEST( testEdge3 );
  CPPUNIT_TEST( testEdge4 );
#endif
 
  CPPUNIT_TEST_SUITE_END();
 
protected:
 
  void testEdge(unsigned n_elem,
                ElemType elem_type)
  {
    Mesh mesh(*TestCommWorld);
 
    // We don't yet support n_connected_components() on distributed
    // meshes
    mesh.allow_remote_element_removal(false);
 
    MeshTools::Generation::build_line(mesh,
                                      n_elem,
                                      /*xmin=*/0.,
                                      /*xmax=*/1.,
                                      elem_type);
 
    const dof_id_type n_components_orig =
      MeshTools::n_connected_components(mesh);
 
    CPPUNIT_ASSERT_EQUAL(n_components_orig, dof_id_type(1));
 
    // We can't disconnect a too-small mesh
    if (n_elem < 3)
      return;
 
    for (const auto & elem : mesh.element_ptr_range())
      {
        const Point c = elem->vertex_average();
 
        if (std::abs(c(0) - 0.5) < 0.2)
          mesh.delete_elem(elem);
      }
 
    mesh.prepare_for_use();
 
    const dof_id_type n_components_broken =
      MeshTools::n_connected_components(mesh);
 
    CPPUNIT_ASSERT_EQUAL(n_components_broken, dof_id_type(2));
 
    // Add a constraint connecting the original end nodes
    auto matrix = SparseMatrix<Number>::build (mesh.comm());
 
    dof_id_type left_node = 0, right_node = 0;
    for (const auto & elem : mesh.element_ptr_range())
      for (const auto & node : elem->node_ref_range())
        {
          if (node(0) == Real(0))
            left_node = node.id();
          if (node(0) == Real(1))
            right_node = node.id();
        }
 
    mesh.comm().max(left_node);
    mesh.comm().max(right_node);
 
    const dof_id_type n_nodes = mesh.n_nodes();
    CPPUNIT_ASSERT_EQUAL(n_nodes, mesh.max_node_id());
 
    processor_id_type my_rank = mesh.comm().rank();
    matrix->init(n_nodes, n_nodes-1,
                 (my_rank ? 0 : n_nodes),
                 (my_rank ? 0 : n_nodes-1),
                 (my_rank ? 0 : n_nodes-1),
                 0);
 
    if (!my_rank)
      {
        for (auto i : make_range(n_nodes))
          {
            if (i < right_node)
              matrix->set(i,i,1);
            if (i > right_node)
              matrix->set(i,i-1,1);
          }
        if (left_node > right_node)
          matrix->set(right_node, left_node-1, 1);
        else
          matrix->set(right_node, left_node, 1);
      }
 
    matrix->close();
 
    mesh.copy_constraint_rows(*matrix);
 
    const dof_id_type n_components_constrained =
      MeshTools::n_connected_components(mesh);
 
    CPPUNIT_ASSERT_EQUAL(n_components_constrained, dof_id_type(1));
  }
 
public:
  void setUp() {}
 
  void tearDown() {}
 
  void testEdge2()
  {
    LOG_UNIT_TEST;
 
    testEdge(/*n_elem=*/10, EDGE2);
  }
 
  void testEdge3()
  {
    LOG_UNIT_TEST;
 
    testEdge(/*n_elem=*/10, EDGE3);
  }
 
  void testEdge4()
  {
    LOG_UNIT_TEST;
 
    testEdge(/*n_elem=*/10, EDGE4);
  }
};
 
 
CPPUNIT_TEST_SUITE_REGISTRATION( ConnectedComponentsTest );