Skip to content
Snippets Groups Projects
Select Git revision
  • gmsh_4_12_2
  • master default protected
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • alphashapes
  • bl
  • relaying
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 3115-issue-fix
  • 3023-Fillet2D-Update
  • convert_fdivs
  • tmp_jcjc24
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
40 results

t14.geo

Blame
  • DivideAndConquer.cpp 23.21 KiB
    // Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
    //
    // See the LICENSE.txt file for license information. Please report all
    // bugs and problems to <gmsh@geuz.org>.
    
    // Triangulation using a divide and conquer algorithm
    //
    // The main idea is to be able to merge two Delaunay triangulations
    // into a single one (see the 'Merge' function). Points are then
    // separated into two groups, then each group into two, ... until
    // having 1, 2 or 3 points (the triangulation is then trivial).
    //
    // This is used to contruct the initial mesh.
    //
    // Warning: point positions must be PERTURBED by a small random
    // value to avoid 3 aligned points or 4 cocyclical points!
    
    #include "GmshMessage.h"
    #include "DivideAndConquer.h"
    #include "Numeric.h"
    #include "robustPredicates.h"
    
    // FIXME: should not introduce dependencies to Geo/ code in Numeric
    #include "GPoint.h"
    #include "GFace.h"
    #include "MLine.h"
    
    #define Pred(x) ((x)->prev)
    #define Succ(x) ((x)->next)
    
    PointNumero DocRecord::Predecessor(PointNumero a, PointNumero b)
    {
      DListPeek p = points[a].adjacent;
      if(p == NULL)
        return -1;
    
      do {
        if(p->point_num == b)
          return Pred(p)->point_num;
        p = Pred(p);
      } while(p != points[a].adjacent);
    
      return -1;
    }
    
    PointNumero DocRecord::Successor(PointNumero a, PointNumero b)
    {
      DListPeek p = points[a].adjacent;
      if(p == NULL)
        return -1;
    
      do {
        if(p->point_num == b)
          return Succ(p)->point_num;
        p = Succ(p);
      } while(p != points[a].adjacent);
    
      return -1;
    }
    
    int DocRecord::FixFirst(PointNumero x, PointNumero f)
    {
      DListPeek p = points[x].adjacent;
      if(p == NULL)
        return 0;
    
      int out = 0;
      DListPeek copy = p;
      do {
        if(p->point_num == f) {
          points[x].adjacent = p;
          out = 1;
        }
        else
          p = p->next;
      } while((p != copy) && !out);
      return out;
    }
    
    PointNumero DocRecord::First(PointNumero x)
    {
      return (points[x].adjacent)->point_num;
    }
    
    int DocRecord::IsLeftOf(PointNumero x, PointNumero y, PointNumero check)
    {
      double pa[2] = {(double)points[x].where.h, (double)points[x].where.v};
      double pb[2] = {(double)points[y].where.h, (double)points[y].where.v};
      double pc[2] = {(double)points[check].where.h, (double)points[check].where.v};
    
      // we use robust predicates here
      double result = robustPredicates::orient2d(pa, pb, pc);
    
      return result > 0;
    }
    
    int DocRecord::IsRightOf(PointNumero x, PointNumero y, PointNumero check)
    {
      return IsLeftOf(y, x, check);
    }
    
    Segment DocRecord::LowerCommonTangent(DT vl, DT vr)
    {
      PointNumero x, y, z, z1, z2, temp;
      Segment s;
    
      x = vl.end;   // vu le tri, c'est le point le + a droite
      y = vr.begin; // idem, le + a gauche
      z = First(y);
      z1 = First(x);
      z2 = Predecessor(x, z1);
      for(;;) {
        if(IsRightOf(x, y, z)) {
          temp = z;
          z = Successor(z, y);
          y = temp;
        }
        else if(IsRightOf(x, y, z2)) {
          temp = z2;
          z2 = Predecessor(z2, x);
          x = temp;
        }
        else {
          s.from = x;
          s.to = y;
          return s;
        }
      }
    }
    
    Segment DocRecord::UpperCommonTangent(DT vl, DT vr)
    {
      PointNumero x, y, z, z1, z2, temp;
      Segment s;
    
      x = vl.end;   // vu le tri, c'est le point le + a droite
      y = vr.begin; // idem, le + a gauche
      z = First(y);
      z1 = First(x);
      z2 = Predecessor(y, z);
      for(;;) {
        if(IsLeftOf(x, y, z2)) {
          temp = z2;
          z2 = Predecessor(z2, y);
          y = temp;
        }
        else if(IsLeftOf(x, y, z1)) {
          temp = z1;
          z1 = Successor(z1, x);
          x = temp;
        }
        else {
          s.from = x;
          s.to = y;
          return s;
        }
      }
    }
    
    // return 1 if the point k is NOT in the circumcircle of triangle  hij
    int DocRecord::Qtest(PointNumero h, PointNumero i, PointNumero j, PointNumero k)
    {
      if((h == i) && (h == j) && (h == k)) {
        Msg::Error("Identical points in triangulation: increase element size "
                   "or Mesh.RandomFactor");
        return 0;
      }
      
      double pa[2] = {(double)points[h].where.h, (double)points[h].where.v};
      double pb[2] = {(double)points[i].where.h, (double)points[i].where.v};
      double pc[2] = {(double)points[j].where.h, (double)points[j].where.v};
      double pd[2] = {(double)points[k].where.h, (double)points[k].where.v};
    
      // we use robust predicates here  
      double result = robustPredicates::incircle(pa, pb, pc, pd) * 
        robustPredicates::orient2d(pa, pb, pc);
      
      return (result < 0) ? 1 : 0;
    }
    
    int DocRecord::Merge(DT vl, DT vr)
    {
      Segment bt, ut;
      int a, b, out;
      PointNumero r, r1, r2, l, l1, l2;
    
      bt = LowerCommonTangent(vl, vr);
      ut = UpperCommonTangent(vl, vr);
      l = bt.from;  // left endpoint of BT
      r = bt.to;    // right endpoint of BT
    
      while((l != ut.from) || (r != ut.to)) {
        a = b = 0;
        if(!Insert(l, r))
          return 0;
    
        r1 = Predecessor(r, l);
        if(r1 == -1)
          return 0;
        if(IsRightOf(l, r, r1))
          a = 1;
        else {
          out = 0;
          while(!out) {
            r2 = Predecessor(r, r1);
            if(r2 == -1)
              return 0;
            if(r2 < vr.begin)
              out = 1;
            else if(Qtest(l, r, r1, r2))
              out = 1;
            else {
              if(!Delete(r, r1))
                return 0;
              r1 = r2;
              if(IsRightOf(l, r, r1))
                out = a = 1;
            }
          }
        }
    
        l1 = Successor(l, r);
        if(l1 == -1)
          return 0;
        if(IsLeftOf(r, l, l1))
          b = 1;
        else {
          out = 0;
          while(!out) {
            l2 = Successor(l, l1);
            if(l2 == -1)
              return 0;
            if(l2 > vl.end)
              out = 1;
            else if(Qtest(r, l, l1, l2))
              out = 1;
            else {
              if(!Delete(l, l1))
                return 0;
              l1 = l2;
              if(IsLeftOf(r, l, l1))
                out = b = 1;
            }
          }
        }
    
        if(a)
          l = l1;
        else if(b)
          r = r1;
        else {
          if(Qtest(l, r, r1, l1))
            r = r1;
          else
            l = l1;
        }
      }
      if(!Insert(l, r))
        return 0;
    
      if(!FixFirst(ut.to, ut.from))
        return 0;
      if(!FixFirst(bt.from, bt.to))
        return 0;
      return 1;
    }
    
    DT DocRecord::RecurTrig(PointNumero left, PointNumero right)
    {
      int n, m;
      DT dt;
    
      dt.begin = left;
      dt.end = right;
      n = right - left + 1; // nombre de points a triangulariser
      switch (n) {
      case 0:
      case 1:
        // 0 ou 1 points -> rien a faire
        break;
    
      case 2: // deux points : cas trivial
        Insert(left, right);
        FixFirst(left, right);
        FixFirst(right, left);
        break;
    
      case 3: // trois points : cas trivial
        Insert(left, right);
        Insert(left, left + 1);
        Insert(left + 1, right);
        if(IsRightOf(left, right, left + 1)) {
          FixFirst(left, left + 1);
          FixFirst(left + 1, right);
          FixFirst(right, left);
        }
        else {
          FixFirst(left, right);
          FixFirst(left + 1, left);
          FixFirst(right, left + 1);
        }
        break;
    
      default: // plus de trois points : cas recursif
        m = (left + right) >> 1;
        if(!Merge(RecurTrig(left, m), RecurTrig(m + 1, right)))
          break;
      }
      return dt;
    }
    
    static int comparePoints(const void *i, const void *j)
    {
      double x, y;
    
      x = ((PointRecord *) i)->where.h - ((PointRecord *) j)->where.h;
      if(x == 0.) {
        y = ((PointRecord *) i)->where.v - ((PointRecord *) j)->where.v;
        return (y < 0.) ? -1 : 1;
      }
      else
        return (x < 0.) ? -1 : 1;
    }
    
    // this fonction builds the delaunay triangulation for a window
    int DocRecord::BuildDelaunay()
    {
      qsort(points, numPoints, sizeof(PointRecord), comparePoints);
      RecurTrig(0, numPoints - 1);
      return 1;
    }
    
    // This routine insert the point 'newPoint' in the list dlist,
    // respecting the clock-wise orientation
    int DocRecord::DListInsert(DListRecord **dlist, DPoint center, PointNumero newPoint)
    {
      DListRecord *p, *newp;
      double alpha1, alpha, beta, xx, yy;
      int first;
    
      newp = new DListRecord;
      newp->point_num = newPoint;
    
      if(*dlist == NULL) {
        *dlist = newp;
        Pred(*dlist) = newp;
        Succ(*dlist) = newp;
        return 1;
      }
      if(Succ(*dlist) == *dlist) {
        Pred(*dlist) = newp;
        Succ(*dlist) = newp;
        Pred(newp) = *dlist;
        Succ(newp) = *dlist;
        return 1;
      }
    
      // If we are here, the double-linked circular list has 2 or more
      // elements, so we have to calculate where to put the new one
    
      p = *dlist;
      first = p->point_num;
    
      // first, compute polar coord. of the first point
      yy = (double)(points[first].where.v - center.v);
      xx = (double)(points[first].where.h - center.h);
      alpha1 = atan2(yy, xx);
    
      // compute polar coord of the point to insert
      yy = (double)(points[newPoint].where.v - center.v);
      xx = (double)(points[newPoint].where.h - center.h);
      beta = atan2(yy, xx) - alpha1;
      if(beta <= 0)
        beta += 2. * M_PI;
    
      do {
        yy = (double)(points[Succ(p)->point_num].where.v - center.v);
        xx = (double)(points[Succ(p)->point_num].where.h - center.h);
        alpha = atan2(yy, xx) - alpha1;
        if(alpha <= 1.e-15)
          alpha += 2. * M_PI;
        if(alpha >= beta) {
          Succ(newp) = Succ(p);
          Succ(p) = newp;
          Pred(newp) = p;
          Pred(Succ(newp)) = newp;
          return 1;
        }
        p = Succ(p);
      } while(p != *dlist);
    
      // never here!
      return 0;
    }
    
    // This routine inserts the point 'a' in the adjency list of 'b' and
    // the point 'b' in the adjency list of 'a'
    int DocRecord::Insert(PointNumero a, PointNumero b)
    {
      int rslt = DListInsert(&points[a].adjacent, points[a].where, b);
      rslt &= DListInsert(&points[b].adjacent, points[b].where, a);
      return rslt;
    }
    
    int DocRecord::DListDelete(DListPeek *dlist, PointNumero oldPoint)
    {
      DListPeek p;
    
      if(*dlist == NULL)
        return 0;
      if(Succ(*dlist) == *dlist) {
        if((*dlist)->point_num == oldPoint) {
          delete *dlist;
          *dlist = NULL;
          return 1;
        }
        else
          return 0;
      }
      p = *dlist;
      do {
        if(p->point_num == oldPoint) {
          Succ(Pred(p)) = Succ(p);
          Pred(Succ(p)) = Pred(p);
          if(p == *dlist) {
            *dlist = Succ(p);
          }
          delete p;
          return 1;
        }
        p = Succ(p);
      } while(p != *dlist);
    
      return 0;
    }
    
    // This routine removes the point 'a' in the adjency list of 'b' and
    // the point 'b' in the adjency list of 'a'
    int DocRecord::Delete(PointNumero a, PointNumero b)
    {
      int rslt = DListDelete(&points[a].adjacent, b);
      rslt &= DListDelete(&points[b].adjacent, a);
      return rslt;
    }
    
    // compte les points sur le polygone convexe
    int DocRecord::CountPointsOnHull()
    {
      PointNumero p, p2, temp;
      int i, n=numPoints;
    
      if(points[0].adjacent == NULL)
        return 0;
      i = 1;
      p = 0;
      p2 = First(0);
      while((p2 != 0) && (i < n)) {
        i++;
        temp = p2;
        p2 = Successor(p2, p);
        p = temp;
      }
      return (i <= n) ? i : -1;
    }
    
    // compte les points sur le polygone convexe
    void DocRecord::ConvexHull()
    {
      PointNumero p, p2, temp;
    
      if(points[0].adjacent == NULL)
        return;
      int count = 0;
      p = 0;
      _hull[count++] = p;
      p2 = First(0);
      while((p2 != 0) && (count < numPoints)) {
        temp = p2;
        p2 = Successor(p2, p);
        p = temp;
        _hull[count++] = p;
      }
    }
    
    PointNumero *DocRecord::ConvertDlistToArray(DListPeek *dlist, int *n)
    {
      DListPeek p, temp;
      int i, max = 0;
      PointNumero *ptr;
    
      p = *dlist;
      do {
        max++;
        p = Pred(p);
      } while(p != *dlist);
      ptr = new PointNumero[max + 1];
      if(ptr == NULL)
        return NULL;
      p = *dlist;
      for(i = 0; i < max; i++) {
        ptr[i] = p->point_num;
        temp = p;
        p = Pred(p);
        delete temp;
      }
      ptr[max] = ptr[0];
      *dlist = NULL;
      *n = max;
      return ptr;
    }
    
    // build ready to use voronoi datas
    void DocRecord::ConvertDListToVoronoiData()
    {
      if (_adjacencies){
        for(int i = 0; i < numPoints; i++)
          if (_adjacencies[i].t) 
            delete [] _adjacencies[i].t;
        delete [] _adjacencies;
      }  
      if (_hull)delete [] _hull;
      _hullSize = CountPointsOnHull ();
      _hull = new PointNumero[_hullSize];
      ConvexHull();
      std::sort(_hull, _hull+_hullSize);
      _adjacencies = new STriangle[numPoints];
      
      for(PointNumero i = 0; i < numPoints; i++) 
        _adjacencies[i].t = ConvertDlistToArray(&points[i].adjacent,
                                                &_adjacencies[i].t_length);    
    }
    
    void DocRecord::voronoiCell(PointNumero pt, std::vector<SPoint2> &pts) const
    {
      if (!_adjacencies){
        Msg::Error("No adjacencies were created");
      }
      const int n = _adjacencies[pt].t_length;
      for(int j = 0; j < n; j++) {
        PointNumero a = _adjacencies[pt].t[j];
        PointNumero b = _adjacencies[pt].t[(j+1) %n];
        double pa[2] = {(double)points[a].where.h, (double)points[a].where.v};
        double pb[2] = {(double)points[b].where.h, (double)points[b].where.v};
        double pc[2] = {(double)points[pt].where.h, (double)points[pt].where.v};
        double center[2];
        circumCenterXY (pa,pb,pc,center);
        pts.push_back(SPoint2(center[0],center[1]));
      }
    }
    
    /*
      Consider N points {X_1, \dots, X_N}
      We want to find the point X_P that verifies
      
      min max | X_i - X_P | , i=1,\dots,N
    
      L2 norm
    
      min \int_V || X - X_P||^2 dv 
    
      =>  2 \int_V || X - X_P|| dv = 0 => X_P is the centroid 
    
      min \int_V || X - X_P||^{2m} dv 
    
      =>  2 \int_V || X - X_P||^{2m-1} dv = 0 => X_P is somewhere ...
    
      now in infinite norm, how to find X_P ?
    
    */
    
    void DocRecord::makePosView(std::string fileName, GFace *gf)
    {
      FILE *f = fopen(fileName.c_str(),"w");
       if (_adjacencies){
        fprintf(f,"View \"voronoi\" {\n");
        for(PointNumero i = 0; i < numPoints; i++) {
          std::vector<SPoint2> pts;
          double pc[2] = {(double)points[i].where.h, (double)points[i].where.v};
          if (!onHull(i)){
    	GPoint p0(pc[0], pc[1], 0.0);
    	//if (gf) p0 = gf->point(pc[0], pc[1]);
            fprintf(f,"SP(%g,%g,%g)  {%g};\n",p0.x(),p0.y(),p0.z(),(double)i);
            voronoiCell (i,pts);
            for (unsigned int j = 0; j < pts.size(); j++){
    	  SPoint2 pp1 = pts[j];
    	  SPoint2 pp2 = pts[(j+1)%pts.size()];
    	  GPoint p1(pp1.x(), pp1.y(), 0.0);
    	  GPoint p2(pp2.x(), pp2.y(), 0.0);
    	  // if (gf) {
    	  //    p1 = gf->point(p1.x(), p1.y());
    	  //    p2 = gf->point(p2.x(), p2.y());
    	  // }
              fprintf(f,"SL(%g,%g,%g,%g,%g,%g)  {%g,%g};\n",
                      p1.x(),p1.y(),p1.z(),p2.x(),p2.y(),p2.z(),
                      (double)i,(double)i);
            }
           }
          else {
            fprintf(f,"SP(%g,%g,%g)  {%g};\n",pc[0],pc[1],0.0,-(double)i);
          }
        }
        fprintf(f,"};\n");    
      }
      fclose(f);
    }
    
    void DocRecord::printMedialAxis(Octree *_octree, std::string fileName, GFace *gf, GEdge *ge)
    {
      
      FILE *f = fopen(fileName.c_str(),"w");
       if (_adjacencies){
        fprintf(f,"View \"medial axis\" {\n");
        for(PointNumero i = 0; i < numPoints; i++) {
          std::vector<SPoint2> pts;
          SPoint2 pc((double)points[i].where.h, (double)points[i].where.v);
          if (!onHull(i)){
    	GPoint p0(pc[0], pc[1], 0.0);
    	if (gf) p0 = gf->point(pc[0], pc[1]);
            fprintf(f,"SP(%g,%g,%g)  {%g};\n",p0.x(),p0.y(),p0.z(),(double)i);
            voronoiCell (i,pts);
            for (unsigned int j = 0; j < pts.size(); j++){
    	  SVector3 pp1(pts[j].x(), pts[j].y(), 0.0);
    	  SVector3 pp2(pts[(j+1)%pts.size()].x(),pts[(j+1)%pts.size()].y(), 0.0);
    	  SVector3 v1(pp1.x()-pc.x(),pp1.y()-pc.y(),0.0);
    	  SVector3 v2(pp2.x()-pc.x(),pp2.y()-pc.y(),0.0);
    	  GPoint p1(pp1.x(), pp1.y(), 0.0);
    	  GPoint p2(pp2.x(), pp2.y(), 0.0);
    	  if (gf) {
    	     p1 = gf->point(p1.x(), p1.y());
    	     p2 = gf->point(p2.x(), p2.y());
    	  }
    	  double P1[3] = {p1.x(), p1.y(), p1.z()};
    	  double P2[3] = {p2.x(), p2.y(), p2.z()};
    	  MElement *m1 = (MElement*)Octree_Search(P1, _octree);
    	  MElement *m2 = (MElement*)Octree_Search(P2, _octree);
    	  if (m1 && m2){
    	    MVertex *v0 = new MVertex(p1.x(), p1.y(), p1.z());
    	    MVertex *v1 = new MVertex(p2.x(), p2.y(), p2.z());
    	    ge->lines.push_back(new MLine(v0, v1));
    	    ge->mesh_vertices.push_back(v0);
    	    ge->mesh_vertices.push_back(v1);
    	    fprintf(f,"SL(%g,%g,%g,%g,%g,%g)  {%g,%g};\n",
    		    p1.x(),p1.y(),p1.z(),
    		    p2.x(),p2.y(),p2.z(),
    		    (double)i,(double)i);
    	 }
            }
           }
        }
        fprintf(f,"};\n");    
      }
      fclose(f);
    
    }
    
    void centroidOfOrientedBox(std::vector<SPoint2> &pts, const double &angle,
                               double &xc, double &yc, double &inertia, double &area)
    {
      const int N = pts.size();
      
      double sina = sin(angle);
      double cosa = cos(angle);
      
      double xmin = cosa* pts[0].x()+ sina*pts[0].y();
      double xmax = cosa* pts[0].x()+ sina*pts[0].y();
      double ymin = -sina* pts[0].x()+ cosa*pts[0].y();
      double ymax = -sina* pts[0].x()+ cosa*pts[0].y();
      for (int j=1;j<N;j++){
        xmin = std::min(cosa* pts[j].x()+ sina*pts[j].y(),xmin);
        ymin = std::min(-sina* pts[j].x()+ cosa*pts[j].y(),ymin);
        xmax = std::max(cosa* pts[j].x()+ sina*pts[j].y(),xmax);
        ymax = std::max(-sina* pts[j].x()+ cosa*pts[j].y(),ymax);
      }
      double XC = 0.5*(xmax+xmin);
      double YC = 0.5*(ymax+ymin);
      xc = XC*cosa - YC *sina;
      yc = XC*sina + YC *cosa;
      inertia = std::max(xmax-xmin,ymax-ymin);
      area = (xmax-xmin) * (ymax-ymin);
    }
    
    void centroidOfPolygon(SPoint2 &pc, std::vector<SPoint2> &pts,
                           double &xc, double &yc, double &inertia, double &areaCell,
                           simpleFunction<double> *bgm)
    {
      double area_tot = 0;
      areaCell = 0.0;
      SPoint2 center(0,0);
      for (unsigned int j = 0; j < pts.size(); j++){
        SPoint2 &pa = pts[j];
        SPoint2 &pb = pts[(j+1)%pts.size()];
        const double area  = triangle_area2d(pa,pb,pc);         
        const double lc = bgm ? (*bgm)((pa.x()+pb.x()+pc.x())/3.0,
                                       (pa.y()+pb.y()+pc.y())/3.0,
                                       0.0) : 1.0;
        const double fact = 1./(lc*lc*lc*lc);//rho = 1/lc^4 (emi)
        areaCell += area;
        area_tot += area*fact;
        center += ((pa+pb+pc) * (area*fact/3.0));
      }
      SPoint2 x = center * (1.0/area_tot);
      inertia = 0;
      for (unsigned int j = 0; j < pts.size(); j++){
        SPoint2 &pa = pts[j];
        SPoint2 &pb = pts[(j+1)%pts.size()];
        const double area  = triangle_area2d(pa,pb,pc);     
        
        const double b = sqrt (  (pa.x()-pb.x())*(pa.x()-pb.x()) + 
                                 (pa.y()-pb.y())*(pa.y()-pb.y()) );
        const double h = 2.0 * area / b;
        const double a = fabs((pb.x()-pa.x())*(pc.x()-pa.x())*
                              (pb.y()-pa.y())*(pc.y()-pa.y()))/b;
        
        const double j2 = (h*b*b*b + h*a*b*b + h*a*a*b + b*h*h*h) / 12.0;
        
        center = (pa+pb+pc) * (1.0/3.0);
        const SPoint2 dx = x - center;
        inertia += j2 + area*area*(dx.x()+dx.x()+dx.y()*dx.y());
      }
      xc = x.x();
      yc = x.y();
    }
    
    double DocRecord::Lloyd(int type)
    {
      fullMatrix<double> cgs(numPoints,2);
      double inertia_tot = 0.;
      for(PointNumero i = 0; i < numPoints; i++) {
        PointRecord &pt = points[i];
        std::vector<SPoint2> pts;
        voronoiCell (i,pts); 
        double E, A;
        
        if (!points[i].data){
          SPoint2 p (pt.where.h,pt.where.v);
          if (type == 0)    
            centroidOfPolygon (p,pts, cgs(i,0), cgs(i,1),E, A);              
          else 
            centroidOfOrientedBox (pts, 0.0, cgs(i,0),cgs(i,1),E, A);    
        }
        inertia_tot += E;
      } 
      
      for(PointNumero i = 0; i < numPoints; i++) {
        if (!points[i].data){
          points[i].where.h = cgs(i,0);
          points[i].where.v = cgs(i,1);
        }
      }
      return inertia_tot;
    }
    
    // Convertir les listes d'adjacence en triangles
    int DocRecord::ConvertDListToTriangles()
    {
      // on suppose que n >= 3. points est suppose OK.
    
      int n = numPoints, i, j;
      int count = 0, count2 = 0;
      PointNumero aa, bb, cc;
    
      STriangle *striangle = new STriangle[n];
    
      count2 = CountPointsOnHull();
    
      // nombre de triangles que l'on doit obtenir
      count2 = 2 * (n - 1) - count2;
    
      // FIXME ::: WHY 2 * ???????????????????
      triangles = new Triangle[2 * count2];
    
      for(i = 0; i < n; i++) {
        // on cree une liste de points connectes au point i (t) + nombre
        // de points (t_length)
        striangle[i].t = ConvertDlistToArray(&points[i].adjacent,
                                             &striangle[i].t_length);
      }
    
      // on balaye les noeuds de gauche a droite -> on cree les triangles
      count = 0;
      for(i = 0; i < n; i++) {
        for(j = 0; j < striangle[i].t_length; j++) {
          if((striangle[i].t[j] > i) && (striangle[i].t[j + 1] > i) &&
             (IsRightOf(i, striangle[i].t[j], striangle[i].t[j + 1]))) {
            aa = i;
            bb = striangle[i].t[j];
            cc = striangle[i].t[j + 1];
            triangles[count].a = aa;
            triangles[count].b = bb;
            triangles[count].c = cc;
            count++;
          }
        }
      }
      numTriangles = count2;
    
      for(int i = 0; i < n; i++)
        delete [] striangle[i].t;
      delete []striangle;
    
      return 1;
    }
    
    //  Cette routine efface toutes les listes d'adjacence de points
    void DocRecord::RemoveAllDList()
    {
      int i;
      DListPeek p, temp;
    
      for(i = 0; i < numPoints; i++)
        if(points[i].adjacent != NULL) {
          p = points[i].adjacent;
          do {
            temp = p;
            p = Pred(p);
            delete temp;
          } while(p != points[i].adjacent);
          points[i].adjacent = NULL;
        }
    }
    
    DocRecord::DocRecord(int n) 
      : _hullSize(0), _hull(NULL), _adjacencies(NULL),
        numPoints(n), points(NULL), numTriangles(0), triangles(NULL)
    {
      if(numPoints)
        points = new PointRecord[numPoints+1000];
    }
    
    DocRecord::~DocRecord()
    {
      if(points) delete [] points;
      if(triangles) delete []triangles;
      if(_hull) delete [] _hull;
      if (_adjacencies){
        for(int i = 0; i < numPoints; i++)
          delete [] _adjacencies[i].t;
        delete _adjacencies;
      }
    }
    
    void DocRecord::MakeMeshWithPoints()
    {
      if(numPoints < 3) return;
      BuildDelaunay();
      ConvertDListToTriangles();
      RemoveAllDList();
    }
    
    void DocRecord::Voronoi()
    {
      if(numPoints < 3) return;
      BuildDelaunay();
      ConvertDListToVoronoiData();
    }
    
    void DocRecord::setPoints(fullMatrix<double> *p)
    { 
      if (numPoints != p->size1())throw;
      for (int i = 0; i < p->size1(); i++){
        x(i) = (*p)(i, 0);
        y(i) = (*p)(i, 1);
        data(i) = (*p)(i, 2) < 0  ? (void *) 1 : NULL;
      }
    } 
    
    void DocRecord::initialize(){
      int i;
      for(i=0;i<numPoints;i++){
    	points[i].flag = 0;
      }
    }
    
    bool DocRecord::remove_point(int index){
      if(points[index].flag == 0){
        points[index].flag = 1;
        return 1;
      }
      else return 0;
    }
    
    void DocRecord::remove_all(){
      int i;
      int index;
      int numPoints2;
      PointRecord* points2;
      numPoints2 = 0;
      for(i=0;i<numPoints;i++){
        if(points[i].flag==0){
          numPoints2++;
        }
      }
      points2 = new PointRecord[numPoints2];
      index = 0;
      for(i=0;i<numPoints;i++){
    	if(points[i].flag==0){
    	  points2[index].where.h = points[i].where.h;
    	  points2[index].where.v = points[i].where.v;
    	  points2[index].data = points[i].data;
    	  points2[index].flag = points[i].flag;
    	  index++;
    	}
      }
      delete [] points;
      points = points2;
      numPoints = numPoints2;
    }
    
    #include "Bindings.h"
    
    void DocRecord::registerBindings(binding *b)
    {
      classBinding *cb = b->addClass<DocRecord>("Triangulator");
      cb->setDescription("A class that does 2D delaunay triangulation "
                         "(JF's SANDBOX for the moment)");
      methodBinding *cm;
    
      cm = cb->addMethod("setPoints", &DocRecord::setPoints);
      cm->setDescription("Set the NumPoints points of the triangulation (x,y,fixed)");
      cm->setArgNames("points",NULL);
      cm = cb->addMethod("Triangulate", &DocRecord::MakeMeshWithPoints);
      cm->setDescription("Compute the Delaunay triangulation");
      cm = cb->addMethod("Voronoi", &DocRecord::Voronoi);
      cm->setDescription("Compute the Voronoi cells");
      cm = cb->addMethod("hullSize", &DocRecord::hullSize);
      cm->setDescription("returns the size of the hull");
      cm = cb->addMethod("makePosView", &DocRecord::makePosView);
      cm->setDescription("save a .pos file with the voronoi in GFace");
      cm->setArgNames("FileName", "GFace", NULL);
      cm = cb->addMethod("Lloyd", &DocRecord::Lloyd);
      cm->setDescription("do one iteration of Lloyd's algorithm");
      cm->setArgNames("Type",NULL);
      cm = cb->setConstructor<DocRecord,int>();
      cm->setDescription ("A Triangulator");
      cm->setArgNames("NumPoints",NULL);
    }