level-set/src/internal_mesh.cpp

/*
    GNU Octave level-set package.
    Copyright (C) 2013-2014  Daniel Kraft <d@domob.eu>

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

/* Internal C++ implementation for mesh building.  Inside and outside elements
   are simply split in two triangles, while boundary elements are split
   depending on how they are intersected by the boundary.  For them, the
   geom.internal data about "boundary segments" is used.  */

#include <octave/oct.h>
#include <octave/Cell.h>
#include <octave/ov-struct.h>

#include <stdexcept>
#include <vector>

#include "Utils.hpp"

/**
 * Type used to hold maps of indices.  It needs to allow negative values, too,
 * since they are used as "indices" of intersection points in the geometry.
 * Furthermore, -1 is used to denote not-yet-set indices (since NA is not
 * available as an integer value).
 */
typedef std::vector<int> indexArr;

/** Type used for coordinates in space.  */
typedef double realT;

/* ************************************************************************** */
/* Containers for the datas of the mesh.  */

/**
 * Represent a vertex in the list of vertices of a mesh (the mesh.p result).
 */
class Vertex
{

private:

  /** This vertex' x coordinate in space.  */
  realT x;
  /** This vertex' y coordinate in space.  */
  realT y;

public:

  /**
   * Construct with given coordinates.
   * @param mx x coordinate.
   * @param my y coordinate.
   */
  inline Vertex (const realT mx, const realT my)
    : x(mx), y(my)
  {
    // Nothing else to do.
  }

  /* Default copying.  */
  Vertex () = default;
  Vertex (const Vertex& o) = default;
  Vertex& operator= (const Vertex& o) = default;

  /**
   * Convert to row in the vertex matrix to be returned to Octave.
   * @param i Index for the row to write in.
   * @param m Matrix to write to.
   */
  inline void
  toMatrix (const unsigned i, Matrix& m) const
  {
    m(0, i) = x;
    m(1, i) = y;
  }

};

/**
 * Represent a triangle in the list of triangles of a mesh (mesh.t field).
 */
class Triangle
{

private:

  /** Indices of the vertices.  */
  unsigned vertices[3];

public:

  /**
   * Construct with given vertices.
   * @param a First vertex.
   * @param b Second vertex.
   * @param c Third vertex.
   */
  inline Triangle (const unsigned a, const unsigned b, const unsigned c)
    : vertices{a, b, c}
  {
    // Nothing else to do.
  }

  /* Default copying.  */
  Triangle () = default;
  Triangle (const Triangle& o) = default;
  Triangle& operator= (const Triangle& o) = default;

  /**
   * Convert to row in the triangle matrix to be returned to Octave.
   * @param i Index for the row to write in.
   * @param m Matrix to write to.
   */
  inline void
  toMatrix (const unsigned i, Matrix& m) const
  {
    for (unsigned j = 0; j < 3; ++j)
      m(j, i) = vertices[j];
  }

};

/* ************************************************************************** */
/* Mesh class itself.  */

/**
 * This class represents the mesh that is being built up.  It keeps track of
 * the vertices and triangles as they are added.  In addition,
 * in order to do that, it also keeps a map between original vertex indices
 * and indices into the already set-up list of vertices.
 */
class Mesh
{

private:

  /** The list of vertices.  */
  std::vector<Vertex> vertices;
  /** The list of triangles.  */
  std::vector<Triangle> triangles;

  /** Map ordinary grid point index to that in list of vertices.  */
  indexArr innerVertexInd;
  /** Map intersection point index to that in list of vertices.  */
  indexArr bdryVertexInd;
  /**
   * Map index of vertices on the mesh to the vertex code on the geometry.
   * I. e., may be positive or negative depending on grid point or ispt.
   */
  indexArr meshVertexInd;

public:

  /**
   * Construct the yet empty mesh.
   * @param n Number of inner grid points.
   * @param nBdry Number of intersection points.
   */
  Mesh (const unsigned n, const unsigned nBdry);

  /* No copying.  */
  Mesh () = delete;
  Mesh (const Mesh& o) = delete;
  Mesh& operator= (const Mesh& o) = delete;

  /**
   * Get internal vertex index for a given grid point.  Insert it into
   * the list of vertices if necessary.
   * @param ind Grid point index.
   * @param coord Full coordinate array, so we know the coordinates.
   * @return Index into internal vertex list.
   */
  unsigned getInnerVertex (const unsigned ind, const Matrix& coord);

  /**
   * Get internal vertex index for a given intersection point.  Insert it into
   * the list of vertices if necessary.
   * @param ind Grid point index.
   * @param coord Array of coordinates of boundary points.
   * @return Index into internal vertex list.
   */
  unsigned getBoundaryVertex (const unsigned ind, const Matrix& coord);

  /**
   * Add a triangle with the given vertices.
   * @param a First vertex.
   * @param b Second vertex.
   * @param c Third vertex.
   */
  inline void
  addTriangle (const unsigned a, const unsigned b, const unsigned c)
  {
    triangles.push_back (Triangle(a, b, c));
  }

  /**
   * Wrap the mesh's content up into the struct expected by Octave.
   * @param res Store the msh structure here.
   * @param vertexMap Store the vertex map here.
   */
  void toStruct (octave_scalar_map& res, octave_scalar_map& vertexMap) const;

};

Mesh::Mesh (const unsigned n, const unsigned nBdry)
  : vertices(), triangles(),
    innerVertexInd(n, -1), bdryVertexInd(nBdry, -1), meshVertexInd()
{
  // Nothing else to do.
}

unsigned
Mesh::getInnerVertex (const unsigned ind, const Matrix& coord)
{
  if (innerVertexInd[ind] < 0)
    {
      const unsigned next = vertices.size ();
      assert (next == meshVertexInd.size ());

      vertices.push_back (Vertex(coord(ind, 0), coord(ind, 1)));
      meshVertexInd.push_back (ind + 1);

      innerVertexInd[ind] = next + 1;
    }

  assert (innerVertexInd[ind] >= 0);
  return innerVertexInd[ind];
}

unsigned
Mesh::getBoundaryVertex (const unsigned ind, const Matrix& coord)
{
  if (bdryVertexInd[ind] < 0)
    {
      const unsigned next = vertices.size ();
      assert (next == meshVertexInd.size ());

      vertices.push_back (Vertex(coord(ind, 0), coord(ind, 1)));
      /* FIXME: Can get rid of cast?  */
      meshVertexInd.push_back (-static_cast<int> (ind) - 1);

      bdryVertexInd[ind] = next + 1;
    }

  assert (bdryVertexInd[ind] >= 0);
  return bdryVertexInd[ind];
}

void
Mesh::toStruct (octave_scalar_map& res, octave_scalar_map& vertexMap) const
{
  Matrix vert(2, vertices.size ());
  Matrix tri(3, triangles.size ());

  for (unsigned i = 0; i < vertices.size (); ++i)
    vertices[i].toMatrix (i, vert);
  for (unsigned i = 0; i < triangles.size (); ++i)
    triangles[i].toMatrix (i, tri);

  vertexMap = octave_scalar_map ();
  vertexMap.assign ("grid", convertToColumnVector (innerVertexInd));
  vertexMap.assign ("ispt", convertToColumnVector (bdryVertexInd));
  vertexMap.assign ("mesh", convertToColumnVector (meshVertexInd));

  res = octave_scalar_map ();
  res.assign ("p", vert);
  res.assign ("t", tri);
  /* e is filled out from geom.bedges in the m-file itself.  */
}

/* ************************************************************************** */

/**
 * Extract the relevant information for building up the triangles from
 * the boundary element segments information provided by ls_find_geometry.
 * It basically just extracts everything of interest and reorders the points.
 * @param segs Boundary element segment info from ls_find_geometry.
 * @param bdryPoints Return boundary points here.
 * @param innerPts Add inner points here.
 */
void
getInnerSegment (const octave_scalar_map& segs,
                 unsigned bdryPts[2], indexArr& innerPts)
{
  bdryPts[0] = segs.contents ("endPt").int_value () - 1;
  bdryPts[1] = segs.contents ("startPt").int_value () - 1;

  assert (innerPts.empty ());
  const ColumnVector inners = segs.contents ("inners").column_vector_value ();
  const unsigned nInners = inners.numel ();
  for (unsigned i = 0; i < nInners; ++i)
    innerPts.push_back (inners(nInners - i - 1) - 1);
}

/* Octave routine for mesh building.  */
DEFUN_DLD (__levelset_internal_mesh, args, nargout,
  "  [MESH, VMAP] = internal_mesh (L, M, INSIDE, COORDS, ISPTCOORDS,\n"
  "                                NODELIST, INELEMS, OUTELEMS,\n"
  "                                BDRYELEMS, BDRYELSEGS, GLOBAL)\n\n"
  "Internal, fast implementation of mesh-building.  Returned are already\n"
  "the mesh structure (except the 'e' field) and vertex map returned\n"
  "from ls_build_mesh, but the arguments are already picked apart by the\n"
  "wrapper routine.\n\n"
  "L, M are the size of the grid, COORDS give the coordinates of grid points\n"
  "in space and INSIDE flags which of those points are inside the domain.\n"
  "ISPTCOORDS are the coordinates of all intersection points indexed by\n"
  "intersection point number.  GLOBAL signals that we want a global mesh\n"
  "instead of a list of meshes for only the inner parts.\n\n"
  "The remaining arguments contain the relevant information from the\n"
  "geometry struct picked apart.\n")
{
  try
    {
      if (args.length () != 11 || nargout > 2)
        throw std::runtime_error ("invalid argument counts");

      /* Get arguments in.  */
      const unsigned L = args(0).int_value ();
      const unsigned M = args(1).int_value ();
      const boolNDArray inside = args(2).bool_array_value ();
      const Matrix coords = args(3).matrix_value ();
      const Matrix isptCoords = args(4).matrix_value ();
      const Matrix nodelist = args(5).matrix_value ();
      const ColumnVector inElems = args(6).column_vector_value ();
      const ColumnVector outElems = args(7).column_vector_value ();
      const ColumnVector bdryElems = args(8).column_vector_value ();
      const Cell bdryelSegs = args(9).cell_value ();
      const bool global = args(10).bool_value ();

      /* Check and extract dimensions.  */
      getDimension (inside, L * M, 1);
      getDimension (coords, L * M, 2);
      const unsigned numIntsec = getDimension (isptCoords, -1, 2);
      getDimension (nodelist, (L - 1) * (M - 1), 4);
      const unsigned numInElems = getDimension (inElems, -1, 1);
      const unsigned numOutElems = getDimension (outElems, -1, 1);
      const unsigned numBdryElems = getDimension (bdryElems, -1, 1);

      /* Construct empty mesh.  */
      Mesh mesh(L * M, numIntsec);

      /* Construct triangles for inner elements.  */
      for (unsigned i = 0; i < numInElems; ++i)
        {
          const unsigned cur = inElems(i) - 1;

          unsigned vertInds[4];
          for (unsigned j = 0; j < 4; ++j)
            vertInds[j] = mesh.getInnerVertex (nodelist(cur, j) - 1, coords);

          mesh.addTriangle (vertInds[0], vertInds[1], vertInds[2]);
          mesh.addTriangle (vertInds[0], vertInds[2], vertInds[3]);
        }

      /* In global mode, do that also for outer elements.  */
      if (global)
        for (unsigned i = 0; i < numOutElems; ++i)
          {
            const unsigned cur = outElems(i) - 1;

            unsigned vertInds[4];
            for (unsigned j = 0; j < 4; ++j)
              vertInds[j] = mesh.getInnerVertex (nodelist(cur, j) - 1, coords);

            mesh.addTriangle (vertInds[0], vertInds[1], vertInds[2]);
            mesh.addTriangle (vertInds[0], vertInds[2], vertInds[3]);
          }

      /* Go over boundary elements and handle them.  We use the information
         provided already by ls_find_geometry internally about the segments
         of each boundary element that belong to the domain.  These need
         to be split into triangles.  */
      for (unsigned i = 0; i < numBdryElems; ++i)
        {
          const unsigned cur = bdryElems(i) - 1;
          const Cell cellSegs = bdryelSegs(i).cell_value ();
          const unsigned nSegs = cellSegs.numel ();

          std::vector<octave_scalar_map> segs;
          indexArr endEdges;
          for (unsigned j = 0; j < nSegs; ++j)
            {
              segs.push_back (cellSegs(j).scalar_map_value ());
              int edge = segs.back ().contents ("endEdge").int_value () - 1;
              endEdges.push_back (edge);
            }

          /* Special case is that of *two* starting points, which
             indicates that we have a narrow pair.  */
          if (segs.size () == 2)
            {
              assert ((endEdges[0] + 2) % 4 == endEdges[1]);

              unsigned bdryPts[4];
              indexArr innerPts[2];
              for (unsigned j = 0; j < 2; ++j)
                {
                  getInnerSegment (segs[j], bdryPts + 2 * j, innerPts[j]);
                  assert (innerPts[j].size () == 1);

                  for (unsigned k = 0; k < 2; ++k)
                    {
                      const unsigned ind = 2 * j + k;
                      bdryPts[ind] = mesh.getBoundaryVertex (bdryPts[ind],
                                                             isptCoords);
                    }
                  for (auto& p : innerPts[j])
                    p = mesh.getInnerVertex (p, coords);

                  mesh.addTriangle (bdryPts[2 * j], innerPts[j][0],
                                    bdryPts[2 * j + 1]);
                }

              if (global)
                {
                  int outerPts[] = {endEdges[0], endEdges[1]};
                  for (auto& p : outerPts)
                    {
                      p = nodelist(cur, p) - 1;
                      assert (!inside (p));
                      p = mesh.getInnerVertex (p, coords);
                    }

                  mesh.addTriangle (bdryPts[1], outerPts[1], bdryPts[2]);
                  mesh.addTriangle (bdryPts[0], bdryPts[1], bdryPts[2]);
                  mesh.addTriangle (bdryPts[0], bdryPts[2], bdryPts[3]);
                  mesh.addTriangle (bdryPts[0], bdryPts[3], outerPts[0]);
                }
            }

          /* "Ordinary" case.  */
          else
            {
              assert (segs.size () == 1);

              /* Find whole list of points forming the segment.  */
              unsigned bdryPts[2];
              indexArr innerPts;
              getInnerSegment (segs[0], bdryPts, innerPts);
              assert (!innerPts.empty ());

              indexArr outerPts;
              for (unsigned p = endEdges[0] % 4; ; p = (p + 3) % 4)
                {
                  const unsigned pInd = nodelist (cur, p) - 1;
                  if (static_cast<int> (pInd) == innerPts.back ())
                    break;
                  assert (!inside (pInd));
                  outerPts.push_back (pInd);
                }
              assert (innerPts.size () + outerPts.size () == 4);

              /* Add to list of vertices in the mesh and overwrite the
                 indices on the way.  */
              for (auto& p : bdryPts)
                p = mesh.getBoundaryVertex (p, isptCoords);
              for (auto& p : innerPts)
                p = mesh.getInnerVertex (p, coords);
              if (global)
                for (auto& p : outerPts)
                  p = mesh.getInnerVertex (p, coords);

              /* Build the triangles accordingly and also the list of
                 boundary edges in this case.  */
              switch (innerPts.size ())
                {
                  case 1:
                    mesh.addTriangle (bdryPts[0], innerPts[0], bdryPts[1]);
                    if (global)
                      {
                        mesh.addTriangle (bdryPts[0], outerPts[1], outerPts[0]);
                        mesh.addTriangle (bdryPts[0], bdryPts[1], outerPts[1]);
                        mesh.addTriangle (bdryPts[1], outerPts[2], outerPts[1]);
                      }
                    break;
                  case 2:
                    mesh.addTriangle (bdryPts[0], innerPts[0], innerPts[1]);
                    mesh.addTriangle (innerPts[1], bdryPts[1], bdryPts[0]);
                    if (global)
                      {
                        mesh.addTriangle (bdryPts[1], outerPts[1], outerPts[0]);
                        mesh.addTriangle (outerPts[0], bdryPts[0], bdryPts[1]);
                      }
                    break;
                  case 3:
                    mesh.addTriangle (bdryPts[0], innerPts[0], innerPts[1]);
                    mesh.addTriangle (innerPts[1], bdryPts[1], bdryPts[0]);
                    mesh.addTriangle (innerPts[1], innerPts[2], bdryPts[1]);
                    if (global)
                      mesh.addTriangle (bdryPts[1], outerPts[0], bdryPts[0]);
                    break;
                  default:
                    assert (false);
                }
            }
        }

      /* Build return results.  */
      octave_scalar_map msh, vertexMap;
      mesh.toStruct (msh, vertexMap);
      octave_value_list res;
      res(0) = msh;
      res(1) = vertexMap;

      return res;
    }
  catch (const std::runtime_error& exc)
    {
      error (exc.what ());
      return octave_value_list ();
    }
}