delaunay_refinement.cpp

#ifdef _OPENMP
#include <omp.h>
#endif
#include <stack>
#include <set>
#include <vector>
#include <algorithm>
#include <math.h>
#include "GmshMessage.h"
#include "OS.h"
#include "SPoint3.h"
#include "delaunay3d_private.h"
#include "delaunay3d.h"
#include "rtree.h"
#include "MVertex.h"
#include "MTetrahedron.h"
#include "MTriangle.h"
#include "GRegion.h"
#include "GFace.h"

typedef  std::set< Edge > edgeContainer2;

long int AVGSEARCH;

struct edgeContainer
{
  std::set< Edge > _hash2;
  std::vector<std::vector<Edge> > _hash;
  edgeContainer (unsigned int N = 1000000) {
    _hash.resize(N);
  }
  bool addNewEdge2 (const Edge &e) {
    std::set< Edge >::iterator it = _hash2.find(e);
    if (it != _hash2.end())return false;
    _hash2.insert(e);
    return true;
  }
  bool isNewEdge (const Edge &e) {
    size_t h = (size_t) e.first >> 3;
    std::vector<Edge> &v = _hash[h %_hash.size()];
    AVGSEARCH+=v.size();
    for (unsigned int i=0; i< v.size();i++)if (e == v[i]) {return false;}
    return true;
  }

  bool addNewEdge (const Edge &e)
  {
    size_t h = (size_t) e.first >> 3;
    std::vector<Edge> &v = _hash[h %_hash.size()];
    AVGSEARCH+=v.size();
    for (unsigned int i=0; i< v.size();i++)if (e == v[i]) {return false;}
    v.push_back(e);
    return true;
  }
};

struct IPT {
  double _x1,_x2,_x3,_x4;
  IPT(double x1, double x2, double x3, double x4) :
    _x1(x1),_x2(x2),_x3(x3),_x4(x4){};
};

double adaptiveTrapezoidalRule (SPoint3 p1 , SPoint3 p2 ,
				double lc1 , double lc2 ,
				double (*f)(const SPoint3 &p, void *),
				void *data, std::vector< IPT > & _result,
				double &dl, double epsilon = 1.e-6)
{
  std::stack< IPT > _stack;
  _result.clear();