SliverChecker.h

Go to the documentation of this file.

#ifndef SliverChecker_HEADER
#define SliverChecker_HEADER

#ifdef HAVE_CONFIG_H
# include <config.h>
#endif

#include "MeshTypes.h"
#include "ProtectedBall.h"
#include "RefineVertex.h"
#include <boost/static_assert.hpp>
#include <geometry/Circles.h>
#include <geometry/predicates.h>
#include <limits>
#include <refine-common/RefineConstants.h>
#include <utilities/SmallList.h>

template <unsigned ambient> struct SliverChecker;

template <>
struct SliverChecker<3> {
  static const unsigned ambient = 3;

  typedef MeshTypes<ambient>::Vertex Vertex;
  typedef MeshTypes<ambient>::Point Point;
  typedef MeshTypes<ambient>::Ball Ball;

  struct Triangle {
    Vertex *const a_, *const b_, *const c_;

    /* The goal is to return 'sigma', which is volume / e^3, where e is the
     * shortest edge.  We compute the edge length in the ctor, and orient3d
     * returns six times the volume, so I store one-sixth the edge cubed. */
    double sixth_e3_;

    Triangle(Vertex *a, Vertex *b, Vertex *c): a_(a), b_(b), c_(c) {
      double ab2 = a->toPoint().dist2(b->toPoint());
      double ac2 = a->toPoint().dist2(c->toPoint());
      double bc2 = b->toPoint().dist2(c->toPoint());
      double min2 = std::min(ab2, std::min(ac2, bc2));
      double min = sqrt(min2);
      sixth_e3_ = min * min * min / 6.0;
    }

    double sigma(const Point& p) const {
      /* Return sigma = volume / e^3, a dimensionless quantity.  Larger is better. */
      double six_volume = - orient3d(a_->x(), b_->x(), c_->x(), p.data());
      if (six_volume < 0.0) {
        // this tet won't be created, so just return that it's perfect quality
        return std::numeric_limits<double>::max();
      } else {
        return six_volume / sixth_e3_;
      }
    }

    /* Circumradius of the new sliver.  Only called if it is, in fact, a sliver;
       therefore, we really shouldn't need to recompute it since it should be
       constant. */
    double radius2(const Point& p) const {
      return Geometry::circumcenter_tetrahedron<ambient>(
          a_->toPoint(), b_->toPoint(), c_->toPoint(), p).radius2;
    }
  };

  template <class ball_iterator>
  SliverChecker(const ball_iterator& begin, const ball_iterator& end, const
      RefineConstants & consts, double radius2):
  consts_(consts), radius2_(radius2) {
    BOUNDED_FOREACH(Ball *b, begin, end) {
      // ignore uppers that aren't up enough to form slivers at all
      if (b->topological() < 3) { continue; }

      // since we're in ambient dimension 3... TODO
      assert(b->topological() == 3);

      // Copy all triangles that this protected simplex can make.
      typename Vertex::VertexUniqueList verts;
      b->gatherVertices(verts);
      assert(verts.size() == 4);
      boost::array<Vertex*,4> varray;
      std::copy(verts.begin(), verts.end(), varray.begin());
      for(int i = 0; i < 4; ++i) {
        for(int j = i+1; j < 4; ++j) {
          for(int k = j+1; k < 4; ++k) {
            Triangle t(varray[i], varray[j], varray[k]);
            triangles_.push_front(t);
          }
        }
      }
    }
  }

  bool formsSliver(const Point& p) const {
    double b2 = consts_.sliverGrowth * consts_.sliverGrowth;
    BOOST_FOREACH(const Triangle& t, triangles_) {
      if (t.sigma(p) < consts_.sigma) {
        // this forms a sliver.  Is it a small sliver?
        if (t.radius2(p) < b2 * radius2_) {
          // it's too small => report!
          return true;
        }
      }
    }
    // Nothing got reported?  Then this is a good point.
    return false;
  }

  private:
  const RefineConstants consts_;
  hudson::SmallList<Triangle> triangles_;
  double radius2_;
};

template <>
struct SliverChecker<2> {
  template <class ball_iterator>
  SliverChecker(const ball_iterator&, const ball_iterator&,
      const RefineConstants &, double) { }

  /* We never form slivers. */
  bool formsSliver(const Geometry::Point<2>&) const {
    return false;
  }
};

template <>
struct SliverChecker<1> {
  template <class ball_iterator>
  SliverChecker(const ball_iterator&, const ball_iterator&, 
      const RefineConstants &, double) { }

  /* We never form slivers. */
  bool formsSliver(const Geometry::Point<1>&) const {
    return false;
  }
};

#endif

---