MeRelaxer.cpp

Go to the documentation of this file.

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------

//----- Included files ---------------------------------------------------------

// 1. Precompiled header

// 2. My own header
#include <xmsmesh/meshing/detail/MeRelaxer.h>

// 3. Standard library headers
#include <cfloat>

// 4. External library headers

// 5. Shared code headers
#include <xmscore/math/math.h>
#include <xmscore/misc/StringUtil.h>
#include <xmscore/misc/XmConst.h>
#include <xmscore/misc/XmError.h>
#include <xmscore/misc/carray.h>
#include <xmscore/misc/xmstype.h> // for XM_PI
#include <xmscore/points/pt.h>
#include <xmscore/stl/vector.h>
#include <xmsinterp/geometry/geoms.h>
#include <xmsinterp/interpolate/InterpBase.h>
#include <xmsinterp/triangulate/TrTin.h>
#include <xmsinterp/triangulate/triangles.h>
#include <xmsmesh/meshing/MePolyRedistributePts.h>

// 6. Non-shared code headers

//----- Forward declarations ---------------------------------------------------

//----- External globals -------------------------------------------------------

//----- Namespace declaration --------------------------------------------------

namespace xms
{
//----- Constants / Enumerations -----------------------------------------------

//----- Classes / Structs ------------------------------------------------------

class MeRelaxerImpl : public MeRelaxer
{
public:
  enum RelaxFlagEnum { RELAX_RELAX = 1 };
  enum RelaxTypeEnum { RELAXTYPE_AREA, RELAXTYPE_ANGLE, RELAXTYPE_SPRING };

  MeRelaxerImpl();
  virtual ~MeRelaxerImpl();

  // protected:
  //  MeRelaxerImpl(const MeRelaxerImpl&);

  XM_DISALLOW_COPY_AND_ASSIGN(MeRelaxerImpl)

  virtual void Relax(/*MeshPolyEnum a_meshPolyEnum,*/
                     const VecInt& a_fixedPoints,
                     BSHP<TrTin> a_tin) override;
  virtual bool SetRelaxationMethod(const std::string& a_relaxMethod) override;
  //------------------------------------------------------------------------------
  //------------------------------------------------------------------------------
  virtual void SetPointSizer(BSHP<MePolyRedistributePts> a_sizer) override { m_sizer = a_sizer; }

  void ComputeCentroids();
  void RelaxMarkedPoints(RelaxTypeEnum a_relaxType, int a_iteration, int a_numiterations);
  void AreaRelax(int a_point, Pt3d& a_newLocation);
  void AngleRelax(int a_point, Pt3d& a_newLocation);
  void SpringRelaxSinglePoint(int a_point, Pt3d& a_newLocation);
  void SetupNeighbors();
  void SetupPointSizes();
  bool NewLocationIsValid(size_t a_idx, Pt3d& a_newLocation);
  bool AllTrianglesHavePositiveArea(BSHP<TrTin> a_tin);

  BSHP<TrTin> m_tin;           
  const VecInt* m_fixedPoints; 
  VecPt3d m_centroids;         
  double m_slideangle;         
  VecInt m_flags;              
  RelaxTypeEnum m_relaxType;   
  BSHP<MePolyRedistributePts> m_sizer; 
  VecDbl m_pointSizes;                 
  VecInt2d m_pointNeighbors;           
  VecInt m_pointsToDelete;             
};                                     // class MeRelaxerImpl

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
MeRelaxerImpl::MeRelaxerImpl()
: m_tin()
/*, m_meshPolyEnum*/
, m_fixedPoints(nullptr)
, m_centroids()
, m_slideangle(5.0)
, m_flags()
, m_relaxType(RELAXTYPE_AREA)
, m_sizer()
, m_pointSizes()
, m_pointsToDelete()
{
} // MeRelaxerImpl::MeRelaxerImpl
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
MeRelaxerImpl::~MeRelaxerImpl()
{
} // MeRelaxerImpl::~MeRelaxerImpl
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeRelaxerImpl::Relax(/*MeshPolyEnum a_meshPolyEnum,*/
                          const VecInt& a_fixedPoints,
                          BSHP<TrTin> a_tin)
{
  XM_ENSURE_TRUE_VOID_NO_ASSERT(a_tin);
  XM_ENSURE_TRUE_VOID(!a_tin->TrisAdjToPts().empty());
  XM_ENSURE_TRUE_VOID(AllTrianglesHavePositiveArea(a_tin));
  // Set up member variables
  // m_meshPolyEnum = a_meshPolyEnum;
  m_tin = a_tin;
  m_fixedPoints = &a_fixedPoints;

  // Mark fixed points to not be relaxed (breaklines)
  m_flags.assign(m_tin->Points().size(), RELAX_RELAX);
  for (size_t i = 0; i < a_fixedPoints.size(); ++i)
  {
    XM_ASSERT(a_fixedPoints[i] >= 0 && a_fixedPoints[i] < (int)m_flags.size());
    m_flags[a_fixedPoints[i]] = 0;
  }

  ComputeCentroids();

  // Set up iterations
#ifdef _DEBUG
  VecPt3d& pts(m_tin->Points());
#endif
  int numiterations(3); // default to 3 for area relaxation
  RelaxTypeEnum relaxtype(m_relaxType);
  if (relaxtype == RELAXTYPE_SPRING)
  {
    if (!m_sizer)
    {
      std::string msg =
        "No size function specified with spring relaxation "
        "method. Relaxation method has been set to AREA "
        "relaxation.";
      XM_LOG(xmlog::warning, msg);
      relaxtype = RELAXTYPE_AREA;
    }
    else
    {
      if (m_pointSizes.empty())
      {
        SetupPointSizes();
      }
      SetupNeighbors();
      XM_ASSERT(pts.size() == m_pointNeighbors.size());
      XM_ASSERT(pts.size() == m_pointSizes.size());
    }
  }
  XM_ASSERT(pts.size() == m_flags.size());

  int iteration = 0;

  // Perform relaxation and swap
  for (int i = 0; i < numiterations; ++i)
  {
    ++iteration;
    RelaxMarkedPoints(relaxtype, iteration, numiterations);
    // delete points and triangles
    if (!m_pointsToDelete.empty())
    {
      m_tin->Triangles().clear();
      m_tin->TrisAdjToPts().clear();
      VecPt3d& pts(m_tin->Points());
      VecDbl& szs(m_pointSizes);
      for (auto it = m_pointsToDelete.rbegin(); it != m_pointsToDelete.rend(); ++it)
      {
        int idx = *it;
        pts[idx] = pts.back();
        pts.pop_back();
        szs[idx] = szs.back();
        szs.pop_back();
        m_flags[idx] = m_flags.back();
        m_flags.pop_back();
      }
      m_pointsToDelete.clear();
      break;
    }
    bool trianglesChanged(false);
    trianglesChanged = m_tin->OptimizeTriangulation();

    if (trianglesChanged && RELAXTYPE_SPRING == relaxtype)
    { // set up spring relax if the triangles have changed
      SetupNeighbors();
    }
  }

} // MeRelaxerImpl::Relax
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool MeRelaxerImpl::SetRelaxationMethod(const std::string& a_relaxType)
{
  bool rval(false);
  std::string type = stToLowerCopy(a_relaxType);
  if (type == "spring_relaxation")
  {
    m_relaxType = RELAXTYPE_SPRING;
    rval = true;
  }
  return rval;
} // MeRelaxerImpl::SetRelaxationMethod
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeRelaxerImpl::ComputeCentroids()
{
  size_t numTri = m_tin->NumTriangles();
  m_centroids.assign(numTri, Pt3d());
  for (size_t t = 0; t < numTri; ++t)
  {
    m_centroids[t] = m_tin->TriangleCentroid((int)t);
  }
} // MeRelaxerImpl::ComputeCentroids
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeRelaxerImpl::RelaxMarkedPoints(RelaxTypeEnum a_relaxType,
                                      int a_iteration,
                                      int a_numiterations)
{
  XM_ASSERT(a_iteration > 0 && a_numiterations > 0 && a_iteration <= a_numiterations);

  // Relax the marked nodes
  VecPt3d& points = m_tin->Points(); // for convenience
  size_t nPoints = points.size();
  Pt3d origlocation, newlocation;
  for (size_t p = 0; p < nPoints; ++p)
  {
    if (m_flags[p])
    {
      origlocation = points[p];
      // do general point relax
      switch (a_relaxType)
      {
      case RELAXTYPE_AREA:
        AreaRelax((int)p, newlocation);
        break;
      case RELAXTYPE_ANGLE:
        AngleRelax((int)p, newlocation);
        break;
      case RELAXTYPE_SPRING:
        SpringRelaxSinglePoint((int)p, newlocation);
        break;
      default:
        XM_ASSERT(false);
        break;
      }
      // preserve the elevation
      newlocation.z = origlocation.z;

      // Change the points location and check for bad triangles
      if (NewLocationIsValid(p, newlocation))
      {
        points[p] = newlocation;
        if (m_sizer)
        {
          // point moved and we must update its target size when doing
          // spring_relax
          m_pointSizes[p] = m_sizer->SizeFromLocation(points[p]);
        }
      }
      else if (RELAXTYPE_SPRING == a_relaxType)
      { // mark point for deletion
        m_pointsToDelete.push_back(static_cast<int>(p));
      }

    } // if (m_flags[p])

  } // for (int p = 0; p < nPoints; ++p)
} // MeRelaxerImpl::RelaxMarkedPoints
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeRelaxerImpl::AreaRelax(int a_point, Pt3d& a_newLocation)
{
  double sumx(0.0), sumy(0.0), area, sumarea(0.0);
  const VecInt& adjTris = m_tin->TrisAdjToPts()[a_point];
  for (size_t i = 0; i < adjTris.size(); ++i)
  {
    area = m_tin->TriangleArea(adjTris[i]);
    sumx += area * m_centroids[adjTris[i]].x;
    sumy += area * m_centroids[adjTris[i]].y;
    sumarea += area;
  }

  a_newLocation.x = sumx / sumarea;
  a_newLocation.y = sumy / sumarea;
} // MeRelaxerImpl::AreaRelax
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeRelaxerImpl::AngleRelax(int a_point, Pt3d& a_newLocation)
{
  int id;
  double targetangle, mag, a, c, d;
  double da, db, dx, dy, sx, sy, radius;
  Pt2d sum, pt, p1, p2, p3, centroid;
  int point1, point2;

  // get the target angle
  const VecInt& adjTris = m_tin->TrisAdjToPts()[a_point];
  VecPt3d& points = m_tin->Points();
  size_t num = adjTris.size();
  XM_ASSERT(num > 0);

  targetangle = 2.0 * XM_PI / (double)num;
  // loop through adjacent elems and find the
  // circle with the two adjacent points and
  // the new point
  sum.x = sum.y = 0.0;
  for (size_t i = 0; i < adjTris.size(); ++i)
  {
    int tri = adjTris[i];
    id = m_tin->LocalIndex(tri, a_point);
    point1 = m_tin->GlobalIndex(tri, trIncrementIndex(id));
    point2 = m_tin->GlobalIndex(tri, trDecrementIndex(id));
    // set up points
    pt.x = points[a_point].x;
    pt.y = points[a_point].y;
    p1.x = points[point1].x;
    p1.y = points[point1].y;
    p2.x = points[point2].x;
    p2.y = points[point2].y;
    // find length of line bisecting p1 and p2
    dx = p1.x - p2.x;
    dy = p1.y - p2.y;
    c = sqrt(sqr(dx) + sqr(dy));
    sx = p1.x + dx / 2.0;
    sy = p1.y + dy / 2.0;
    // compute length of triangle sides
    a = (c / 2.0) / sin(targetangle / 2.0);
    d = a * cos(targetangle / 2.0);
    radius = (a * a) / (2.0 * d);
    // find candidate point p3 and circle centroid
    dx /= c;
    dy /= c;
    p3.x = sx + d * dy;
    p3.y = sy - d * dx;
    da = sqr(p3.x - pt.x) + sqr(p3.y - pt.y);
    p3.x = sx - d * dy;
    p3.y = sy + d * dx;
    db = sqr(p3.x - pt.x) + sqr(p3.y - pt.y);
    if (da < db)
    {
      centroid.x = sx + (d - radius) * dy;
      centroid.y = sy - (d - radius) * dx;
    }
    else
    {
      centroid.x = sx - (d - radius) * dy;
      centroid.y = sy + (d - radius) * dx;
    }
    // find p3 along circle, radius away from centroid
    dx = pt.x - centroid.x;
    dy = pt.y - centroid.y;
    mag = sqrt(dx * dx + dy * dy);
    dx /= mag;
    dy /= mag;
    // find new point
    p3.x = centroid.x + radius * dx;
    p3.y = centroid.y + radius * dy;
    // sum new point's x and y
    sum.x += p3.x;
    sum.y += p3.y;
  }
  a_newLocation.x = (sum.x / (double)num);
  a_newLocation.y = (sum.y / (double)num);
} // MeRelaxerImpl::AngleRelax
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeRelaxerImpl::SpringRelaxSinglePoint(int a_point, Pt3d& a_newLocation)
{
  double fx(0.0);
  double fy(0.0);

  double springStiffness(1.0);
  VecPt3d& points = m_tin->Points();
  VecInt& neighbors(m_pointNeighbors[a_point]);
  double numPtsFactor(1.7 / (double)neighbors.size());
  Pt3d& p0(points[a_point]);
  double size0(m_pointSizes[a_point]);
  for (size_t i = 0; i < neighbors.size(); ++i)
  {
    int neighborIdx = neighbors[i];
    Pt3d& p1(points[neighborIdx]);
    double size1(m_pointSizes[neighborIdx]);
    Pt3d vec = gmCreateVector(p0, p1);
    double dist = Mdist(p0.x, p0.y, p1.x, p1.y);
    double targetDist = 0.5 * (size0 + size1); 
    double force(springStiffness * (dist - targetDist));
    double _fx = force * vec.x / dist;
    double _fy = force * vec.y / dist;
    fx += _fx;
    fy += _fy;
  }
  fx *= numPtsFactor;
  fy *= numPtsFactor;
  a_newLocation.x = p0.x + fx;
  a_newLocation.y = p0.y + fy;
} // MeRelaxerImpl::SpringRelax
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeRelaxerImpl::SetupNeighbors()
{
  m_pointNeighbors.clear();
  VecPt3d& points = m_tin->Points();
  size_t nPts(points.size());
  if (m_flags.empty())
  {
    m_flags.assign(nPts, RELAX_RELAX);
  }
  m_pointNeighbors.assign(nPts, VecInt());
  VecInt& tris = m_tin->Triangles();
  VecInt2d& adjTris2d = m_tin->TrisAdjToPts();
  for (size_t i = 0; i < points.size(); ++i)
  {
    if (m_flags[i] == 0) // We don't relax fixed points, so we don't need its neighbors.
    {
      continue;
    }
    const VecInt& adjTris = adjTris2d[i];
    size_t numTri = adjTris.size();
    // find the other points that make up the edges that a_point is connected to
    std::set<int> setPoints;
    for (size_t j = 0; j < numTri; ++j)
    {
      int idx = adjTris[j] * 3;
      for (int t = 0; t < 3; ++t)
      {
#ifdef _DEBUG
        int ptIdx = tris[idx + t];
#endif
        XM_ASSERT(ptIdx < nPts);
        setPoints.insert(tris[idx + t]);
      }
    }
    setPoints.erase(static_cast<int>(i));
    m_pointNeighbors[i] = VecInt(setPoints.begin(), setPoints.end());
  }
} // MeRelaxerImpl::SetupNeighbors()
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeRelaxerImpl::SetupPointSizes()
{
  m_pointSizes.clear();
  VecPt3d& points = m_tin->Points();
  size_t nPts(points.size());
  m_pointSizes.reserve(nPts);
  for (size_t i = 0; i < points.size(); ++i)
  {
    // compute point sizes prior to relaxing below
    Pt3d& p(m_tin->Points()[i]);
    m_pointSizes.push_back(m_sizer->SizeFromLocation(p));
  }
} // MeRelaxerImpl::SetupPointSizes
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool MeRelaxerImpl::NewLocationIsValid(size_t a_idx, Pt3d& a_newLocation)
{
  Pt3d originalLocation = m_tin->Points()[a_idx];
  m_tin->Points()[a_idx] = a_newLocation;
  const VecInt& adjTris = m_tin->TrisAdjToPts()[a_idx];
  bool rval(true);
  for (size_t i = 0; rval && i < adjTris.size(); ++i)
  {
    // Is triangle ill-formed? See that it has some area
    if (!GT_EPS(m_tin->TriangleArea(adjTris[i]), 0.0, FLT_EPSILON))
    {
      rval = false;
    }
  }
  m_tin->Points()[a_idx] = originalLocation;
  return rval;
} // MeRelaxerImpl::NewLocationIsValid
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool MeRelaxerImpl::AllTrianglesHavePositiveArea(BSHP<TrTin> a_tin)
{
  for (size_t i = 0; i < a_tin->NumTriangles(); ++i)
  {
    if (a_tin->TriangleArea(static_cast<int>(i)) < 0.0)
      return false;
  }
  return true;
} // MeRelaxerImpl::AllTrianglesHavePositiveArea

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
BSHP<MeRelaxer> MeRelaxer::New()
{
  BSHP<MeRelaxer> polyMesher(new MeRelaxerImpl);
  // MeRelaxerImpl->SetBSHP(polyMesher); // If MeRelaxerImpl needs ptr to
  // MeRelaxer
  return polyMesher;
} // MeRelaxer::New
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
MeRelaxer::MeRelaxer()
{
} // MeRelaxer::MeRelaxer
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
MeRelaxer::~MeRelaxer()
{
} // MeRelaxer::~MeRelaxer

} // namespace xms

#if CXX_TEST
// UNIT TESTS

#include <xmsmesh/meshing/detail/MeRelaxer.t.h>

#include <xmscore/testing/TestTools.h>
#include <xmsinterp/triangulate/TrTriangulatorPoints.h>

//----- Namespace declaration --------------------------------------------------

// namespace xms {
using namespace xms;

//------------------------------------------------------------------------------
//  40-   20-----21-----22-----23-----24
//    |   |\2   1|\2   1|\0   2|\1   0|
//    |   |0\ 12 |1\ 26 |2\ 25 |1\ 28 |
//    |   |  \   |  \   |  \   |  \   |
//    |   | 8 \  |   \  |   \  |   \  |
//    |   |    \0| 11 \0| 23 \1| 27 \2|
//    |   |1   2\|2   0\|0   1\|2   0\|
//  30-   15-----16-----17-----18-----19
//    |   |\2   1|\2   1|0   2/|\1   0|
//    |   |0\    |1\ 31 | 24 /2|2\ 29 |
//    |   |  \ 9 |  \   |   /  |  \   |
//    |   | 7 \  |   \  |  /   |   \  |
//    |   |    \0| 10 \0|1/ 22 | 30 \2|
//    |   |1   2\|2   0\|/0   1|0   1\|
//  20-   10-----11-----12-----13-----14
//    |   |0   2/|\1   0|\0   2|\0   2|
//    |   | 13 /2|2\ 14 |2\ 18 |2\ 21 |
//    |   |   /  |  \   |  \   |  \   |
//    |   |  / 6 | 5 \  |   \  |   \  |
//    |   |1/    |    \2| 15 \1| 19 \1|
//    |   |/0   1|0   1\|0   1\|0   1\|
//  10-   5------6------7------8------9
//    |   |\0   2|\0   2|\1   0|\1   0|
//    |   |2\    |1\    |0\ 16 |0\ 20 |
//    |   |  \ 1 |  \ 2 |  \   |  \   |
//    |   | 0 \  | 3 \  | 4 \  |   \  |
//    |   |    \1|    \1|    \2| 17 \2|
//    |   |0   1\|2   0\|1   2\|1   2\|
//  0 -   0------1------2------3------4
//
//        |------|------|------|------|
//        0     10     20     30     40
//------------------------------------------------------------------------------
void MeRelaxerUnitTests::testRelaxWhileMeshing()
{
  // Create tin and get some convenience variables
  BSHP<TrTin> tin = TrTin::New();
  VecInt2d& trisAdjToPts = tin->TrisAdjToPts();

  // Set up to our example tin points
  tin->Points().reserve(5 * 5);
  for (size_t i = 0; i < 5; ++i)
  {
    for (size_t j = 0; j < 5; ++j)
    {
      tin->Points().push_back(Pt3d(j * 10.0, i * 10.0, static_cast<double>(j)));
    }
  }

  // Triangulate the points
  TrTriangulatorPoints client(tin->Points(), tin->Triangles(), &trisAdjToPts);
  client.Triangulate();

  // Set up fixed points (outer boundary)
  int outPolyA[] = {0, 1, 2, 3, 4, 9, 14, 19, 24, 23, 22, 21, 20, 15, 10, 5};
  VecInt outPoly(&outPolyA[0], &outPolyA[XM_COUNTOF(outPolyA)]);

  BSHP<MeRelaxer> relaxer = MeRelaxer::New();
  relaxer->Relax(outPoly, tin);
  const double kDelta = 1e-2;
  VecPt3d& tinPts = tin->Points();
  VecPt3d expectedPts{
    {0.00, 0.00, 0.00},   {10.00, 0.00, 1.00},  {20.00, 0.00, 2.00},  {30.00, 0.00, 3.00},
    {40.00, 0.00, 4.00},  {0.00, 10.00, 0.00},  {11.10, 8.88, 1.00},  {20.02, 10.20, 2.00},
    {29.94, 9.92, 3.00},  {40.00, 10.00, 4.00}, {0.00, 20.00, 0.00},  {9.33, 19.09, 1.00},
    {20.68, 20.86, 2.00}, {31.15, 18.86, 3.00}, {40.00, 20.00, 4.00}, {0.00, 30.00, 0.00},
    {9.94, 29.75, 1.00},  {18.83, 31.19, 2.00}, {29.14, 29.13, 3.00}, {40.00, 30.00, 4.00},
    {0.00, 40.00, 0.00},  {10.00, 40.00, 1.00}, {20.00, 40.00, 2.00}, {30.00, 40.00, 3.00},
    {40.00, 40.00, 4.00}};

  TS_ASSERT_DELTA_VECPT3D(expectedPts, tinPts, kDelta);

} // MeRelaxerUnitTests::testRelaxWhileMeshing
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeRelaxerUnitTests::testSpringRelaxSetup()
{
  VecPt3d points = {{-10, -10, 0}, {0, -10, 0},  {10, -10, 0}, {-10, 0, 0}, {0, 0, 0},
                    {10, 0, 0},    {-10, 10, 0}, {0, 10, 0},   {10, 10, 0}, {-5, 10, 0}};
  VecInt tris = {0, 4, 3, 0, 1, 4, 1, 2, 4, 2, 5, 4, 3, 7, 9, 4, 8, 7, 4, 5, 8, 3, 9, 6, 3, 4, 7};
  // create a test tin
  BSHP<TrTin> tin = TrTin::New();
  tin->Points() = points;
  tin->Triangles() = tris;
  tin->BuildTrisAdjToPts();
  // create a sizer that returns a constant value
  BSHP<MePolyRedistributePts> sizer = MePolyRedistributePts::New();
  sizer->SetConstantSizeFunc(3.0);

  MeRelaxerImpl r;
  r.m_tin = tin;
  r.m_sizer = sizer;
  r.SetupPointSizes();
  r.SetupNeighbors();

  VecDbl expectedPointSizes(points.size(), 3);
  TS_ASSERT_EQUALS_VEC(expectedPointSizes, r.m_pointSizes);
  VecInt2d expectedPointNeighbors = {
    {1, 3, 4}, {0, 2, 4}, {1, 4, 5},    {0, 4, 6, 7, 9}, {0, 1, 2, 3, 5, 7, 8},
    {2, 4, 8}, {3, 9},    {3, 4, 8, 9}, {4, 5, 7},       {3, 6, 7}};
  TS_ASSERT_EQUALS_VEC2D(expectedPointNeighbors, r.m_pointNeighbors);
} // MeRelaxerUnitTests::testSpringRelaxSetup
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeRelaxerUnitTests::testSpringRelaxSetup2()
{
  VecPt3d points = {{-10, -10, 0}, {0, -10, 0},  {10, -10, 0}, {-10, 0, 0}, {0, 0, 0},
                    {10, 0, 0},    {-10, 10, 0}, {0, 10, 0},   {10, 10, 0}, {-5, 10, 0}};
  VecInt tris = {0, 4, 3, 0, 1, 4, 1, 2, 4, 2, 5, 4, 4, 8, 7, 4, 5, 8, 3, 9, 6};
  // create a test tin
  BSHP<TrTin> tin = TrTin::New();
  tin->Points() = points;
  tin->Triangles() = tris;
  tin->BuildTrisAdjToPts();
  // create a sizer that returns a constant value
  BSHP<MePolyRedistributePts> sizer = MePolyRedistributePts::New();
  sizer->SetConstantSizeFunc(3.0);

  MeRelaxerImpl r;
  r.m_tin = tin;
  r.m_sizer = sizer;
  r.SetupPointSizes();
  r.SetupNeighbors();

  VecDbl expectedPointSizes(points.size(), 3);
  TS_ASSERT_EQUALS_VEC(expectedPointSizes, r.m_pointSizes);
  VecInt2d expectedPointNeighbors = {
    {1, 3, 4}, {0, 2, 4}, {1, 4, 5}, {0, 4, 6, 9}, {0, 1, 2, 3, 5, 7, 8},
    {2, 4, 8}, {3, 9},    {4, 8},    {4, 5, 7},    {3, 6}};
  TS_ASSERT_EQUALS_VEC2D(expectedPointNeighbors, r.m_pointNeighbors);
} // MeRelaxerUnitTests::testSpringRelaxSetup
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeRelaxerUnitTests::testSpringRelaxSinglePoint()
{
  VecPt3d points = {{-10, -10, 1}, {0, -10, 2},  {10, -10, 3}, {-10, 0, 4}, {8, 7, 5},
                    {10, 0, 6},    {-10, 10, 7}, {0, 10, 8},   {10, 10, 9}};
  VecInt tris = {0, 4, 3, 0, 1, 4, 1, 2, 4, 2, 5, 4, 3, 4, 6, 4, 7, 6, 4, 5, 8, 4, 8, 7};
  // create a test tin
  BSHP<TrTin> tin = TrTin::New();
  tin->Points() = points;
  tin->Triangles() = tris;
  tin->BuildTrisAdjToPts();
  // create a sizer that returns a constant value
  BSHP<MePolyRedistributePts> sizer = MePolyRedistributePts::New();
  sizer->SetConstantSizeFunc(10.0);

  MeRelaxerImpl r;
  r.m_tin = tin;
  r.m_sizer = sizer;
  r.SetupPointSizes();
  r.SetupNeighbors();

  Pt3d newloc, startPt(8, 7, 0), ptZero;
  double startDist(Mdist(0.0, 0.0, startPt.x, startPt.y));
  r.SpringRelaxSinglePoint(4, newloc);
  double newDist(Mdist(0.0, 0.0, newloc.x, newloc.y));
  TS_ASSERT(newDist < startDist);
  TS_ASSERT(r.NewLocationIsValid(4, newloc));

  Pt3d expectedNewLoc(0.905, 0.542, 0);
  TS_ASSERT_DELTA_PT3D(expectedNewLoc, newloc, 1e-2);
  TS_ASSERT(r.NewLocationIsValid(4, newloc));

  tin->Points()[4] = newloc;
  startDist = newDist;
  r.SpringRelaxSinglePoint(4, newloc);
  newDist = Mdist(0.0, 0.0, newloc.x, newloc.y);
  TS_ASSERT(newDist < startDist);
  expectedNewLoc = Pt3d(0.023, 0.016, 0);
  TS_ASSERT_DELTA_PT3D(expectedNewLoc, newloc, 1e-2);

  tin->Points()[4] = newloc;
  startDist = newDist;
  r.SpringRelaxSinglePoint(4, newloc);
  newDist = Mdist(0.0, 0.0, newloc.x, newloc.y);
  TS_ASSERT(newDist < startDist);
  expectedNewLoc = Pt3d(0, 0, 0);
  TS_ASSERT_DELTA_PT3D(expectedNewLoc, newloc, 1e-2);
  TS_ASSERT(r.NewLocationIsValid(4, newloc));
} // MeRelaxerUnitTests::testSpringRelaxSinglePoint
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeRelaxerUnitTests::testSpringRelaxSinglePoint2()
{
  VecPt3d points = {{-10, -10, 0}, {0, -10, 0},  {10, -10, 0}, {-10, 0, 0}, {8, 7, 0},
                    {15, 0, 0},    {-10, 10, 0}, {0, 10, 0},   {10, 10, 0}};
  VecInt tris = {0, 4, 3, 0, 1, 4, 1, 2, 4, 2, 5, 4, 3, 4, 6, 4, 7, 6, 4, 5, 8, 4, 8, 7};
  // create a test tin
  BSHP<TrTin> tin = TrTin::New();
  tin->Points() = points;
  tin->Triangles() = tris;
  tin->BuildTrisAdjToPts();
  // create a sizer that returns a constant value
  BSHP<MePolyRedistributePts> sizer = MePolyRedistributePts::New();
  sizer->SetConstantSizeFunc(10.0);

  MeRelaxerImpl r;
  r.m_tin = tin;
  r.m_sizer = sizer;
  r.SetupPointSizes();
  r.SetupNeighbors();
  r.m_pointSizes[5] = 15;

  Pt3d newloc, startPt(8, 7, 0), ptZero;
  double startDist(Mdist(0.0, 0.0, startPt.x, startPt.y));
  r.SpringRelaxSinglePoint(4, newloc);
  double newDist(Mdist(0.0, 0.0, newloc.x, newloc.y));
  TS_ASSERT(newDist < startDist);
  TS_ASSERT(r.NewLocationIsValid(4, newloc));

  Pt3d expectedNewLoc(0.673, 0.377, 0);
  TS_ASSERT_DELTA_PT3D(expectedNewLoc, newloc, 1e-2);
  TS_ASSERT(r.NewLocationIsValid(4, newloc));

  tin->Points()[4] = newloc;
  startDist = newDist;
  r.SpringRelaxSinglePoint(4, newloc);
  newDist = Mdist(0.0, 0.0, newloc.x, newloc.y);
  TS_ASSERT(newDist < startDist);
  expectedNewLoc = Pt3d(0.547, -0.005, 0);
  TS_ASSERT_DELTA_PT3D(expectedNewLoc, newloc, 1e-2);
  TS_ASSERT(r.NewLocationIsValid(4, newloc));

  tin->Points()[4] = newloc;
  startDist = newDist;
  r.SpringRelaxSinglePoint(4, newloc);
  newDist = Mdist(0.0, 0.0, newloc.x, newloc.y);
  TS_ASSERT(newDist < startDist);
  expectedNewLoc = Pt3d(0.545, 0, 0);
  TS_ASSERT_DELTA_PT3D(expectedNewLoc, newloc, 1e-2);
  TS_ASSERT(r.NewLocationIsValid(4, newloc));

} // MeRelaxerUnitTests::testSpringRelaxSinglePoint2
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeRelaxerUnitTests::testSpringRelaxSinglePoint3()
{
  VecPt3d points = {{-10, -10, 0}, {0, -10, 0},  {10, -10, 0}, {-10, 0, 0}, {8, 7, 0},
                    {10, 0, 0},    {-10, 10, 0}, {0, 10, 0},   {10, 10, 0}};
  VecInt tris = {0, 1, 3, 1, 4, 3, 1, 5, 4, 1, 2, 5, 3, 4, 6, 4, 7, 6, 4, 5, 8, 4, 8, 7};
  // create a test tin
  BSHP<TrTin> tin = TrTin::New();
  tin->Points() = points;
  tin->Triangles() = tris;
  tin->BuildTrisAdjToPts();
  // create a sizer that returns a constant value
  BSHP<MePolyRedistributePts> sizer = MePolyRedistributePts::New();
  sizer->SetConstantSizeFunc(100.0);

  MeRelaxerImpl r;
  r.m_tin = tin;
  r.m_sizer = sizer;
  r.SetupPointSizes();
  r.SetupNeighbors();

  Pt3d newloc, startPt(8, 7, 0), ptZero;
  r.SpringRelaxSinglePoint(4, newloc);
  TS_ASSERT(!r.NewLocationIsValid(4, newloc));
} // MeRelaxerUnitTests::testSpringRelaxSinglePoint3
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeRelaxerUnitTests::testNewLocationIsValid()
{
  VecPt3d points = {{-10, -10, 0}, {0, -10, 0},  {10, -10, 0}, {-10, 0, 0}, {0, 0, 0},
                    {5, 0, 0},     {-10, 10, 0}, {0, 10, 0},   {10, 10, 0}};
  VecInt tris = {0, 4, 3, 0, 1, 4, 1, 2, 4, 2, 5, 4, 4, 5, 7, 5, 8, 7, 3, 4, 7, 3, 7, 6};
  // create a test tin
  BSHP<TrTin> tin = TrTin::New();
  tin->Points() = points;
  tin->Triangles() = tris;
  tin->BuildTrisAdjToPts();
  // create a sizer that returns a constant value
  BSHP<MePolyRedistributePts> sizer = MePolyRedistributePts::New();
  sizer->SetConstantSizeFunc(3.0);

  MeRelaxerImpl r;
  r.m_tin = tin;
  r.m_sizer = sizer;
  r.SetupPointSizes();
  r.SetupNeighbors();

  Pt3d newloc(-2, -2, 0), ptZero;
  TS_ASSERT(r.NewLocationIsValid(4, newloc));

  tin->Points()[4] = ptZero;
  newloc = Pt3d(-2, 2, 0);
  TS_ASSERT(r.NewLocationIsValid(4, newloc));

  tin->Points()[4] = ptZero;
  newloc = Pt3d(-7, 7, 0);
  TS_ASSERT(!r.NewLocationIsValid(4, newloc));

  tin->Points()[4] = ptZero;
  newloc = Pt3d(7, 0, 0);
  TS_ASSERT(!r.NewLocationIsValid(4, newloc));

  tin->Points()[4] = ptZero;
  newloc = Pt3d(5, 5, 0);
  TS_ASSERT(!r.NewLocationIsValid(4, newloc));

  tin->Points()[4] = ptZero;
  newloc = Pt3d(-10, 0, 0);
  TS_ASSERT(!r.NewLocationIsValid(4, newloc));

  tin->Points()[4] = ptZero;
  newloc = Pt3d(-11, 0, 0);
  TS_ASSERT(!r.NewLocationIsValid(4, newloc));

  tin->Points()[4] = ptZero;
  newloc = Pt3d(-2, -11, 0);
  TS_ASSERT(!r.NewLocationIsValid(4, newloc));

} // MeRelaxerUnitTests::testNewLocationIsValid
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeRelaxerUnitTests::testAllTrianglesHavePositiveArea()
{
  VecPt3d points = {{-10, -10, 0}, {0, -10, 0},  {10, -10, 0}, {-10, 0, 0}, {8, 7, 0},
                    {10, 0, 0},    {-10, 10, 0}, {0, 10, 0},   {10, 10, 0}};
  VecInt tris = {0, 4, 3, 0, 1, 4, 1, 2, 4, 2, 5, 4, 3, 4, 6, 4, 7, 6, 4, 5, 8, 4, 8, 7};
  // create a test tin
  BSHP<TrTin> tin = TrTin::New();
  tin->Points() = points;
  tin->Triangles() = tris;
  tin->BuildTrisAdjToPts();

  MeRelaxerImpl r;
  TS_ASSERT(r.AllTrianglesHavePositiveArea(tin));
  // create an invalid triangle
  tris[1] = 3;
  tris[2] = 4;
  tin->Triangles() = tris;
  TS_ASSERT(!r.AllTrianglesHavePositiveArea(tin));
} // MeRelaxerUnitTests::testAllTrianglesHavePositiveArea

//} // namespace xms
#endif // CXX_TEST