boundary_ids_test.cc

#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/ucd_io.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/elem.h"

using namespace libMesh;

// Uncomment this or compile with -DOLD_STYLE_BOUNDARY_IDS to test
// old-style boundary_ids() function.
// #define OLD_STYLE_BOUNDARY_IDS 1

// Build a mesh, add lots of boundary IDs to different sides and
// nodes, time how long it takes to call the old and new versions of
// boundary_ids() (both pre- and post-patch).
//
// Timing results
// n. boundary ids  |  std::set (new style)                        | std::vector (old style)
// -----------------------------------------------------------------------------------------
// 10               | 0.010026,0.012408,0.012863,0.009972,0.010144 | 0.005616,0.004844,0.005025,0.005055,0.004839
// 20               | 0.024317,0.020987,0.021466,0.020348,0.023336 | 0.008200,0.008237,0.008183,0.007636,0.008194
// 40               | 0.040803,0.041244,0.040899,0.039981,0.040076 | 0.012971,0.011933,0.011961,0.012046,0.013681
// 80               | 0.084363,0.083282,0.084440,0.083495,0.083140 | 0.022848,0.024616,0.024637,0.024383,0.024176
// 160              | 0.174268,0.166348,0.175957,0.178958,0.171316 | 0.039649,0.048735,0.038687,0.046747,0.048566
//
// Unfortunately, it appears the old return-vector-by-value version is faster
// over a wide range of boundary ids.
int main (int argc, char** argv)
{
  LibMeshInit init(argc, argv);

  Mesh mesh(init.comm());
  MeshTools::Generation::build_cube (mesh,
                                     /*nx=*/10, /*ny=*/10, /*nz=*/10,
                                     /*xmin=*/0, /*xmax=*/1.,
                                     /*ymin=*/0, /*ymax=*/1.,
                                     /*zmin=*/0, /*zmax=*/1.,
                                     HEX8);

  {
    // Number of boundary IDs to add to each side in the Mesh.  The
    // performance of the different containers may depend on the
    // number of boundary IDs present for a given side.
    unsigned n_boundary_ids=160;

    // Add lots of boundary_id info to every side in the Mesh.  We could
    // instead randomize this, etc...
    MeshBase::element_iterator           el = mesh.level_elements_begin(0);
    const MeshBase::element_iterator end_el = mesh.level_elements_end(0);
    for (; el != end_el; ++el)
      {
        Elem *elem = *el;

        unsigned int n_sides = elem->n_sides();
        for (unsigned short s=0; s != n_sides; ++s)
          {
            // Add n_boundary_ids IDs to each side in the mesh.  Start
            // with 6, since build_cube presumably uses 0-5 already.
            for (boundary_id_type bid=6; bid<=6+n_boundary_ids; ++bid)
              mesh.get_boundary_info().add_side(elem, s, bid);
          }
      }
  }

  // Verify that worked as expected.
  // mesh.get_boundary_info().print_info();

  // We will accumulate ids into this total to try and prevent unused
  // containers from being optimized away.
  unsigned long running_sum = 0;

  PerfLog perf_log ("Boundary Info");

  perf_log.start_event("Get all boundary ids");

#ifdef OLD_STYLE_BOUNDARY_IDS

  // Here's a typical "old style" code snippet that loops over all the
  // sides of a level-0 element and requests boundary IDs.  In this
  // code snippet, we are also replacing new ids with old ones, but
  // that is not strictly required for the test. However, we probably
  // need to do *something* with the ids that are returned, to prevent
  // the container from being optimized away...
  MeshBase::element_iterator           el = mesh.level_elements_begin(0);
  const MeshBase::element_iterator end_el = mesh.level_elements_end(0);
  for (; el != end_el; ++el)
    {
      Elem *elem = *el;

      unsigned int n_sides = elem->n_sides();
      for (unsigned short s=0; s != n_sides; ++s)
        {
          std::vector<boundary_id_type> ids =
            mesh.get_boundary_info().boundary_ids(elem, s);

          // Accumulate into running_sum to avoid optimizing away an
          // unused container. Don't time this, though.
          perf_log.pause_event("Get all boundary ids");
          for (unsigned i=0; i<ids.size(); ++i)
            running_sum += ids[i];
          perf_log.restart_event("Get all boundary ids");
        }
    }

#else // OLD_STYLE_BOUNDARY_IDS==0

  // Here's a new style code snippet that does the same thing, but
  // with the std::set interface.  We probably don't need to do the
  // same we did with the vector to ensure that the container is not
  // optimized away, but we need to make sure both cases do the same
  // amount of work.
  MeshBase::element_iterator           el = mesh.level_elements_begin(0);
  const MeshBase::element_iterator end_el = mesh.level_elements_end(0);

  // Container used to hold boundary ids for different geometric entities.
  std::set<boundary_id_type> bndry_ids;

  for (; el != end_el; ++el)
    {
      Elem *elem = *el;

      unsigned int n_sides = elem->n_sides();
      for (unsigned short s=0; s != n_sides; ++s)
        {
          mesh.get_boundary_info().boundary_ids(elem, s, bndry_ids);

          // Accumulate a running sum to be consistent with the std::vector
          // case. Don't time this, though.
          perf_log.pause_event("Get all boundary ids");
          for (std::set<boundary_id_type>::iterator i=bndry_ids.begin();
               i != bndry_ids.end(); ++i)
            running_sum += *i;
          perf_log.restart_event("Get all boundary ids");
        }
    }

#endif

  perf_log.stop_event("Get all boundary ids");

  libMesh::out << "running_sum=" << running_sum << std::endl;

  return 0;
}