bassi_rebay_mesh.cc

#include "libmesh/node.h"
#include "libmesh/elem.h"
#include "libmesh/mesh.h"
#include "libmesh/face_tri3.h"
#include "libmesh/getpot.h"

using namespace libMesh;

// Commands to approximately generate the Bassi & Rebay meshes:
// In my opinion, the meshes look best when nodes_per_ring = 2 * num_rings
// ./bassi_rebay_mesh-opt --inner-radius 1 --outer-radius 40 --nodes-per-ring 16 --num-rings 8
// ./bassi_rebay_mesh-opt --inner-radius 1 --outer-radius 40 --nodes-per-ring 32 --num-rings 16
// ./bassi_rebay_mesh-opt --inner-radius 1 --outer-radius 40 --nodes-per-ring 64 --num-rings 32
// ./bassi_rebay_mesh-opt --inner-radius 1 --outer-radius 40 --nodes-per-ring 128 --num-rings 64
// ./bassi_rebay_mesh-opt --inner-radius 1 --outer-radius 40 --nodes-per-ring 256 --num-rings 128
int main (int argc, char ** argv)
{
  LibMeshInit init (argc, argv);

  // Create a GetPot object to parse the command line
  GetPot command_line (argc, argv);

  // Table of values to approximately reproduce Bassi and Rebay's meshes.
  // inner_radius & outer_radius & nodes_per_ring & num_rings   & Total N. Elem & radial_bias (computed)
  // 1            & 40           & 16             & 8           & 256           & 1.70329
  // 1            & 40           & 32             & 16          & 1024          & 1.28961
  // 1            & 40           & 64             & 32          & 4096          & 1.13253
  // 1            & 40           & 128            & 64          & 16384         & 1.06351
  // 1            & 40           & 256            & 128         & 65536         & 1.0311

  // Size of the cylinder whose exterior is being meshed.
  Real inner_radius = 1.0;
  if (command_line.search(1, "--inner-radius"))
    inner_radius = command_line.next(inner_radius);

  // How far away is the "far field". In Bassi & Rebay, pg. 258, it states
  // that "all grids extend about 20 diameters away from the circle."
  Real outer_radius = 40.0;
  if (command_line.search(1, "--outer-radius"))
    outer_radius = command_line.next(outer_radius);

  // Factor to increase initial_delta_r for each ring.
  // Uniform grid: radial_bias = 1.0
  Real radial_bias = 1.2;
  bool user_set_radial_bias = false;
  if (command_line.search(1, "--radial-bias"))
    {
      radial_bias = command_line.next(radial_bias);
      user_set_radial_bias = true;
    }

  // Number of nodes on each ring.
  int nodes_per_ring = 16;
  if (command_line.search(1, "--nodes-per-ring"))
    nodes_per_ring = command_line.next(nodes_per_ring);

  // Generate mesh of TRI6 elements instead of TRI3 elements.
  bool second_order = false;
  if (command_line.search(1, "--second-order"))
    second_order = true;

  // Read the number of rings from the command line. You can't specify
  // both the number of rings and the radial bias if you want to match
  // a specified outer radius exactly... you have to leave one of
  // those parameters free so that it can be determined.
  int num_rings = 18;
  bool user_set_num_rings = false;
  if (command_line.search(1, "--num-rings"))
    {
      num_rings = command_line.next(num_rings);
      user_set_num_rings = true;
    }

  // Input check: can't have 0 rings.
  if (num_rings == 0)
    libmesh_error_msg("Cannot create mesh with num_rings == 0.");

  // Width of the initial layer of elements around the cylinder.  This
  // _could_ be user-selectable, but we set it to to ensure that the
  // initial layer of elements are nearly equilateral triangles.
  Real initial_delta_r = 2 * libMesh::pi * inner_radius / nodes_per_ring;

  // Input check: if the user set both the radial bias and the number
  // of rings, they may not be consistent. We let the number of rings
  // take precendence and compute the required radial bias. If it does
  // not match the user's selected radial bias, print a warning message.

  // The boundary id to use for the cylinder
  boundary_id_type
    cylinder_bid = 1,
    exterior_bid = 2;

  {
    // Compute the radial bias given:
    // .) the inner radius
    // .) the outer radius
    // .) the initial_delta_r
    // .) the desired number of intervals
    // Note: the exponent n used in the formula is one less than the
    // number of rings the user requests.
    Real alpha = 1.1;
    int n = num_rings - 1;

    // lambda used to compute the residual and Jacobian for the Newton iterations.
    // We capture parameters which don't need to change from the current scope at
    // the time this lambda is declared. The values are not updated later, so we
    // can't use this for e.g. f, df, and alpha.
    auto newton = [outer_radius, inner_radius, initial_delta_r, n](Real & f, Real & df, const Real & alpha)
      {
        f  = (1. - std::pow(alpha, n + 1)) / (1. - alpha) - (outer_radius - inner_radius) / initial_delta_r;
        df = (-(n + 1) * (1 - alpha) * std::pow(alpha, n) + (1. - std::pow(alpha, n + 1))) / (1. - alpha) / (1. - alpha);
      };

    Real f, df;
    newton(f, df, alpha);
    // std::cout << "Initial |residual| = " << std::abs(f) << std::endl;
    // std::cout << "Initial derivative = " << df << std::endl;

    while (std::abs(f) > 1.e-9)
      {
        // Compute and apply update.
        Real dx = -f/df;
        alpha += dx;
        newton(f, df, alpha);
        // std::cout << "dx=" << dx << std::endl;
        // std::cout << "abs(f)=" << std::abs(f) << std::endl;
      }
    std::cout << "Computed radial bias = " << alpha << std::endl;

    if (user_set_radial_bias && std::abs(radial_bias - alpha) > TOLERANCE)
      std::cout << "Warning: user-selected radial bias " << radial_bias << " will not hit the outer_radius exactly. Using " << alpha << " instead." << std::endl;

    // Set radial basis to the value of alpha that we computed with Newton.
    radial_bias = alpha;
  }

  // The number of rings specified by the user does not include the ring at
  // the surface of the cylinder itself, so we increment it by one now.
  num_rings += 1;

  // Mesh we are eventually going to create.
  Mesh mesh(init.comm());

  // Data structure that holds pointers to the Nodes of each ring.
  std::vector<std::vector<Node *>> ring_nodes(num_rings);

  // Initialize radius and delta_r variables.
  Real radius = inner_radius;
  Real delta_r = initial_delta_r;

  // Node id counter.
  unsigned int current_node_id = 0;

  for (std::size_t r=0; r<num_rings; ++r)
    {
      ring_nodes[r].resize(nodes_per_ring);

      // Debugging
      // libMesh::out << "Adding nodes for ring " << r << " at radius " << radius << std::endl;

      // Add nodes starting from either theta=0 or theta=pi/nodes_per_ring
      Real theta = r%2 == 0 ? 0 : (libMesh::pi / nodes_per_ring);
      for (std::size_t n=0; n<nodes_per_ring; ++n)
        {
          // Debugging:
          // libMesh::out << "Add node at theta = " << 180 * theta / libMesh::pi << std::endl;

          ring_nodes[r][n] = mesh.add_point(Point(radius * std::cos(theta),
                                                  radius * std::sin(theta)),
                                            current_node_id++);
          // Update angle
          theta += 2 * libMesh::pi / nodes_per_ring;

          // Debugging
          // libMesh::out << *ring_nodes[r][n] << std::endl;
        }

      // Go to next ring
      radius += delta_r;
      delta_r *= radial_bias;
    }

  // Add elements
  for (std::size_t r=0; r<num_rings - 1; ++r)
    {
      // even -> odd ring
      if (r%2 == 0)
        {
          // Inner ring (n, n*, n+1)
          // Starred indices refer to nodes on the "outer" ring of this pair.
          for (std::size_t n=0; n<nodes_per_ring; ++n)
            {
              // Wrap around
              unsigned int np1 = (n == nodes_per_ring - 1) ? 0 : n + 1;
              Elem * elem = mesh.add_elem(new Tri3);
              elem->set_node(0) = ring_nodes[r][n];
              elem->set_node(1) = ring_nodes[r+1][n];
              elem->set_node(2) = ring_nodes[r][np1];

              // Add interior faces to 'cylinder' sideset if we are on ring 0.
              if (r == 0)
                mesh.boundary_info->add_side(elem->id(), /*side=*/2, cylinder_bid);
            }

          // Outer ring (n*, n+1*, n+1)
          for (std::size_t n=0; n<nodes_per_ring; ++n)
            {
              // Wrap around
              unsigned int np1 = (n == nodes_per_ring - 1) ? 0 : n + 1;
              Elem * elem = mesh.add_elem(new Tri3);
              elem->set_node(0) = ring_nodes[r+1][n];
              elem->set_node(1) = ring_nodes[r+1][np1];
              elem->set_node(2) = ring_nodes[r][np1];

              // Add exterior faces to 'exterior' sideset if we're on the last ring.
              // Note: this code appears in two places since we could end on either an even or odd ring.
              if (r == num_rings - 2)
                mesh.boundary_info->add_side(elem->id(), /*side=*/0, exterior_bid);
            }
        }
      else
        {
          // odd -> even ring
          // Inner ring (n, n+1*, n+1)
          for (std::size_t n=0; n<nodes_per_ring; ++n)
            {
              // Wrap around
              unsigned int np1 = (n == nodes_per_ring - 1) ? 0 : n + 1;
              Elem * elem = mesh.add_elem(new Tri3);
              elem->set_node(0) = ring_nodes[r][n];
              elem->set_node(1) = ring_nodes[r+1][np1];
              elem->set_node(2) = ring_nodes[r][np1];
            }

          // Outer ring (n*, n+1*, n)
          for (std::size_t n=0; n<nodes_per_ring; ++n)
            {
              // Wrap around
              unsigned int np1 = (n == nodes_per_ring - 1) ? 0 : n + 1;
              Elem * elem = mesh.add_elem(new Tri3);
              elem->set_node(0) = ring_nodes[r+1][n];
              elem->set_node(1) = ring_nodes[r+1][np1];
              elem->set_node(2) = ring_nodes[r][n];

              // Add exterior faces to 'exterior' sideset if we're on the last ring.
              if (r == num_rings - 2)
                mesh.boundary_info->add_side(elem->id(), /*side=*/0, exterior_bid);
            }
        }
    }

  // Sanity check: make sure all elements have positive area. Note: we
  // can't use elem->volume() for this, as that always returns a
  // positive area regardless of the node ordering.
  // We compute (p1-p0) \cross (p2-p0) and check that the z-component is positive.
  for (const auto & elem : mesh.element_ptr_range())
    {
      Point cp = (elem->point(1) - elem->point(0)).cross(elem->point(2) - elem->point(0));
      if (cp(2) < 0.)
        libmesh_error_msg("Invalid elem " << elem->id() << " found with negative area = " << cp(2));
    }

  // Create sideset names.
  mesh.boundary_info->sideset_name(cylinder_bid) = "cylinder";
  mesh.boundary_info->sideset_name(exterior_bid) = "exterior";

  // Find neighbors, etc.
  mesh.prepare_for_use();

  if (second_order)
    {
      mesh.all_second_order(/*full_ordered=*/true);
      std::vector<unsigned int> nos;

      // Loop over the elements, moving mid-edge nodes onto the
      // nearest radius as applicable. For each element, exactly one
      // edge should lie on the same radius, so we move only that
      // mid-edge node.
      for (const auto & elem : mesh.element_ptr_range())
        {
          // Make sure we are dealing only with triangles
          libmesh_assert(elem->n_vertices() == 3);

          // Compute vertex radii
          Real radii[3] = {elem->point(0).norm(), elem->point(1).norm(), elem->point(2).norm()};

          // Compute absolute differences between radii so we can determine which two are on the same circular arc.
          Real dr[3] = {std::abs(radii[0] - radii[1]), std::abs(radii[1] - radii[2]), std::abs(radii[2] - radii[0])};

          // Compute index of minimum dr.
          auto index = std::distance(std::begin(dr), std::min_element(std::begin(dr), std::end(dr)));

          // Make sure that the minimum found is also (almost) zero.
          if (dr[index] > TOLERANCE)
            libmesh_error_msg("Error: element " << elem->id() << " had no sides with nodes on same radius.");

          // Get list of all local node ids on this side. The first
          // two entries in nos correspond to the vertices, the last
          // entry corresponds to the mid-edge node.
          nos = elem->nodes_on_side(index);

          // Compute the angles associated with nodes nos[0] and nos[1].
          Real
            theta0 = std::atan2(elem->point(nos[0])(1), elem->point(nos[0])(0)),
            theta1 = std::atan2(elem->point(nos[1])(1), elem->point(nos[1])(0));

          // atan2 returns values in the range (-pi, pi).  If theta0
          // and theta1 have the same sign, we can simply average them
          // to get half of the acute angle between them. On the other
          // hand, if theta0 and theta1 are of opposite sign _and_ both
          // are larger than pi/2, we need to add 2*pi when averaging,
          // otherwise we will get half of the _obtuse_ angle between
          // them, and the point will flip to the other side of the
          // circle (see below).
          Real new_theta = 0.5 * (theta0 + theta1);

          // It should not be possible for both:
          // 1.) |theta0| > pi/2, and
          // 2.) |theta1| < pi/2
          // as this would not be a well-formed element.
          if ((theta0*theta1 < 0) &&
              (std::abs(theta0) > 0.5*libMesh::pi) &&
              (std::abs(theta1) > 0.5*libMesh::pi))
            new_theta = 0.5 * (theta0 + theta1 + 2*libMesh::pi);

          // The new radius will be the radius of point nos[0] or nos[1] (they are the same!).
          Real new_r = elem->point(nos[0]).norm();

          // Finally, move the point to its new location.
          elem->point(nos[2]) = Point(new_r * std::cos(new_theta), new_r * std::sin(new_theta), 0.);
        }
    }

  mesh.write("bassi_rebay_mesh.e");

  // Print summary
  libMesh::out << "Number of elements: " << mesh.n_elem() << std::endl;

  return 0;
}