Add #include <vector> for compilation on linux

CMakeLists.txt

@@ -1,14 +1,25 @@
 # Set  the minimum  required version  of cmake  for a  project.
-cmake_minimum_required(VERSION 3.5)
+cmake_minimum_required(VERSION 3.6)
+
+include_directories(./)
+
+set(INC
+  tetgen.h
+)
+
+set(SRC
+  tetgen.cxx
+  predicates.cxx
+)
 
 # Add an executable to the project using the specified source files.
-add_executable(tetgen tetgen.cxx predicates.cxx)
+add_executable(tetgen ${INC} ${SRC})
 
 #Add a library to the project using the specified source files. 
 # In Linux/Unix, it will creates the libtet.a
-add_library(tet STATIC tetgen.cxx predicates.cxx)
+#add_library(tet STATIC tetgen.cxx predicates.cxx)
 
 #Set properties on a target. 
 #We use this here to set -DTETLIBRARY for when compiling the
 #library
-set_target_properties(tet PROPERTIES "COMPILE_DEFINITIONS" TETLIBRARY)
+#set_target_properties(tet PROPERTIES "COMPILE_DEFINITIONS" TETLIBRARY)

example.poly

@@ -1,84 +0,0 @@
-28  3  0  1
-1  0  0  0  1
-2  2  0  0  1
-3  2  2  0  1
-4  0  2  0  1
-5  0  0  4  9
-6  2  0  4  9
-7  2  2  3  9
-8  0  2  3  9
-9  0  0  5  2
-10  2  0  5  2
-11  2  2  5  2
-12  0  2  5  2
-13  0.25  0.25  0.5  4
-14  1.75  0.25  0.5  4
-15  1.75  1.5  0.5  4
-16  0.25  1.5  0.5  4
-17  0.25  0.25  1  4
-18  1.75  0.25  1  4
-19  1.75  1.5  1  4
-20  0.25  1.5  1  4
-21  0.25  0  2  4
-22  1.75  0  2  4
-23  1.75  1.5  2  4
-24  0.25  1.5  2  4
-25  0.25  0  2.5  4
-26  1.75  0  2.5  4
-27  1.75  1.5  2.5  4
-28  0.25  1.5  2.5  4
-23  1
-1  0  1  # 1
-4    1  2  3  4
-1  0  9  # 2
-4    5  6  7  8
-2  1  3  # 3
-4    1  2  6  5
-4    21  22  26  25
-1  1  0  2.25
-1  0  3  # 4
-4    2  3  7  6
-1  0  3  # 5
-4    3  4  8  7
-1  0  3  # 6
-4    4  1  5  8
-1  0  2  # 7
-4    9  10  11  12
-1  0  3  # 8
-4    9  10  6  5
-1  0  3  # 9
-4    10  11  7  6
-1  0  3  # 10
-4    11  12  8  7
-1  0  3  # 11
-4    12  9  5  8
-1  0  4  # 12
-4    13  14  15  16
-1  0  4  # 13
-4    17  18  19  20
-1  0  4  # 14
-4    13  14  18  17
-1  0  4  # 15
-4    14  15  19  18
-1  0  4  # 16
-4    15  16  20  19
-1  0  4  # 17
-4    16  13  17  20
-1  0  4  # 18
-4    21  22  23  24
-1  0  4  # 19
-4    25  26  27  28
-1  0  4  # 20
-4    21  22  26  25
-1  0  4  # 21
-4    22  23  27  26
-1  0  4  # 22
-4    23  24  28  27
-1  0  4  # 23
-4    24  21  25  28
-2
-1  1  0.4  2.25
-2  1  0.4  0.75
-2
-1  1  0.25  0.1  10  0.001
-2  1  0.5  4  20  0.01

predicates.cxx

@@ -2190,6 +2190,9 @@ REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
 
 REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
 {
+//#ifdef USING_GMP
+//  return orient3d_gmp(pa, pb, pc, pd);
+//#endif
   REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
   REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
   REAL det;
@@ -4703,7 +4706,623 @@ REAL orient4d(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL* pe,
 
 
 
+//==============================================================================
 
+static REAL det3x3(REAL adx, REAL ady, REAL adz,
+                   REAL bdx, REAL bdy, REAL bdz,
+                   REAL cdx, REAL cdy, REAL cdz)
+{
+  return adx * (bdy * cdz - bdz * cdy)
+       + bdx * (cdy * adz - cdz * ady)
+       + cdx * (ady * bdz - adz * bdy);
+}
+
+void tetgenmesh::pre_calculate_insphere(point pa, point pb, point pc, point pd,
+                                        REAL* dets)
+{
+  if (pd != dummypoint) {
+    REAL ba_x = pb[0] - pa[0];
+    REAL ba_y = pb[1] - pa[1];
+    REAL ba_z = pb[2] - pa[2];
+    REAL ba_norm = ba_x*ba_x + ba_y*ba_y + ba_z*ba_z;
+    
+    REAL ca_x = pc[0] - pa[0];
+    REAL ca_y = pc[1] - pa[1];
+    REAL ca_z = pc[2] - pa[2];
+    REAL ca_norm = ca_x*ca_x + ca_y*ca_y + ca_z*ca_z;
+    
+    REAL da_x = pd[0] - pa[0];
+    REAL da_y = pd[1] - pa[1];
+    REAL da_z = pd[2] - pa[2];
+    REAL da_norm = da_x*da_x + da_y*da_y + da_z*da_z;
+
+    dets[0] = det3x3(ba_y, ba_z, ba_norm,
+                     ca_y, ca_z, ca_norm,
+                     da_y, da_z, da_norm);
+    
+    dets[1] = det3x3(ba_x, ba_z, ba_norm,
+                     ca_x, ca_z, ca_norm,
+                     da_x, da_z, da_norm);
+
+    dets[2] = det3x3(ba_x, ba_y, ba_norm,
+                     ca_x, ca_y, ca_norm,
+                     da_x, da_y, da_norm);
+    
+    dets[3] = det3x3(ba_x, ba_y, ba_z,
+                     ca_x, ca_y, ca_z,
+                     da_x, da_y, da_z);
+  } else {
+    double ab[4],ac[4];
+    double* a = pa; // mesh->vertices[Node[0]].coord;
+    double* b = pb; //mesh->vertices[Node[1]].coord;
+    double* c = pc; //mesh->vertices[Node[2]].coord;
+    unsigned i;
+    for (i=0; i<3; i++)
+    {
+      ab[i]=b[i]-a[i]; //AB
+      ac[i]=c[i]-a[i]; //AC
+    }
+    dets[0] = ac[1]*ab[2] - ac[2]*ab[1];
+    dets[1] = ac[2]*ab[0] - ac[0]*ab[2];
+    dets[2] = ac[0]*ab[1] - ac[1]*ab[0];
+    dets[3] = dets[0]*dets[0] + dets[1]*dets[1] + dets[2]*dets[2];
+  }
+}
+
+REAL tetgenmesh::insphere_use_subdets(tetrahedron *tet, REAL* pe)
+{
+  REAL *dets = get_polar(tet);
+
+  if (dets[3] == 0) { // Only calculate once.
+    point *pts = (point *) tet;
+    pre_calculate_insphere(pts[4], pts[5], pts[6], pts[7], dets);
+  }
+
+  point pa = (point) tet[4];
+  
+  if (((point) tet[7]) == dummypoint) {
+    double aex = pe[0] - pa[0];
+    double aey = pe[1] - pa[1];
+    double aez = pe[2] - pa[2];
+    double det = aex*dets[0]+aey*dets[1]+aez*dets[2];
+    if(fabs(det) > o3dstaticfilter) return det;
+    point *pts = (point *) tet;
+    det = orient3d(pts[4],pts[4],pts[6],pe);
+    return det;
+  }
+  
+  REAL ea_x = pe[0] - pa[0];
+  REAL ea_y = pe[1] - pa[1];
+  REAL ea_z = pe[2] - pa[2];
+  REAL ea_norm = ea_x * ea_x + ea_y * ea_y + ea_z * ea_z;
+
+  REAL det = -ea_x * dets[0]
+           +  ea_y * dets[1]
+           -  ea_z * dets[2]
+           +  ea_norm * dets[3];
+
+  if (fabs(det) < ispstaticfilter) {
+    point *pts = (point *) tet;
+    det = insphere_s(pts[4], pts[5], pts[6], pts[7], pe);
+  }
+  return det;
+}
+
+#ifdef USING_GMP
+
+#include <gmp.h>
+
+//============================================================================//
+//                                                                            //
+// orient3d_gmp()                                                             //
+//                                                                            //
+//============================================================================//
+
+REAL orient3d_gmp(REAL *aa, REAL *bb, REAL *cc, REAL *dd)
+{
+  mpz_t pa[3], pb[3], pc[3], pd[3];
+  
+  for (int i = 0; i < 3; i++) {
+    mpz_init_set_d(pa[i], aa[i]*1.e+15);
+    mpz_init_set_d(pb[i], bb[i]*1.e+15);
+    mpz_init_set_d(pc[i], cc[i]*1.e+15);
+    mpz_init_set_d(pd[i], dd[i]*1.e+15);
+  }
+  
+  mpz_t adx, bdx, cdx;
+  mpz_t ady, bdy, cdy;
+  mpz_t adz, bdz, cdz;
+  
+  //adx = pa[0] - pd[0];
+  mpz_init(adx);
+  mpz_sub(adx, pa[0], pd[0]);
+  //bdx = pb[0] - pd[0];
+  mpz_init(bdx);
+  mpz_sub(bdx, pb[0], pd[0]);
+  //cdx = pc[0] - pd[0];
+  mpz_init(cdx);
+  mpz_sub(cdx, pc[0], pd[0]);
+  //ady = pa[1] - pd[1];
+  mpz_init(ady);
+  mpz_sub(ady, pa[1], pd[1]);
+  //bdy = pb[1] - pd[1];
+  mpz_init(bdy);
+  mpz_sub(bdy, pb[1], pd[1]);
+  //cdy = pc[1] - pd[1];
+  mpz_init(cdy);
+  mpz_sub(cdy, pc[1], pd[1]);
+  //adz = pa[2] - pd[2];
+  mpz_init(adz);
+  mpz_sub(adz, pa[2], pd[2]);
+  //bdz = pb[2] - pd[2];
+  mpz_init(bdz);
+  mpz_sub(bdz, pb[2], pd[2]);
+  //cdz = pc[2] - pd[2];
+  mpz_init(cdz);
+  mpz_sub(cdz, pc[2], pd[2]);
+
+  /*
+  return adx * (bdy * cdz - bdz * cdy)
+       + bdx * (cdy * adz - cdz * ady)
+       + cdx * (ady * bdz - adz * bdy);
+   */
+  mpz_t bdy_cdz;
+  mpz_init(bdy_cdz);
+  mpz_mul(bdy_cdz, bdy, cdz);
+  mpz_t bdz_cdy;
+  mpz_init(bdz_cdy);
+  mpz_mul(bdz_cdy, bdz, cdy);
+  
+  mpz_t cdy_adz;
+  mpz_init(cdy_adz);
+  mpz_mul(cdy_adz, cdy, adz);
+  mpz_t cdz_ady;
+  mpz_init(cdz_ady);
+  mpz_mul(cdz_ady, cdz, ady);
+  
+  mpz_t ady_bdz;
+  mpz_init(ady_bdz);
+  mpz_mul(ady_bdz, ady, bdz);
+  mpz_t adz_bdy;
+  mpz_init(adz_bdy);
+  mpz_mul(adz_bdy, adz, bdy);
+  
+  mpz_t bdy_cdz_bdz_cdy;
+  mpz_init(bdy_cdz_bdz_cdy);
+  mpz_sub(bdy_cdz_bdz_cdy, bdy_cdz, bdz_cdy);
+  
+  mpz_t cdy_adz_cdz_ady;
+  mpz_init(cdy_adz_cdz_ady);
+  mpz_sub(cdy_adz_cdz_ady, cdy_adz, cdz_ady);
+  
+  mpz_t ady_bdz_adz_bdy;
+  mpz_init(ady_bdz_adz_bdy);
+  mpz_sub(ady_bdz_adz_bdy, ady_bdz, adz_bdy);
+  
+  mpz_t adx_;
+  mpz_init(adx_);
+  mpz_mul(adx_, adx, bdy_cdz_bdz_cdy);
+  mpz_t bdx_;
+  mpz_init(bdx_);
+  mpz_mul(bdx_, bdx, cdy_adz_cdz_ady);
+  mpz_t cdx_;
+  mpz_init(cdx_);
+  mpz_mul(cdx_, cdx, ady_bdz_adz_bdy);
+  
+  mpz_t det;
+  mpz_init(det);
+  mpz_add(det, adx_, bdx_);
+  mpz_add(det, det, cdx_);
+  
+  bool debug_flag = false;
+  if (debug_flag) { // Debug only
+    char str[1024];
+    mpz_get_str(str, 10, det);
+    printf("\ndet_str = %s\n", str);
+    double detd = mpz_get_d(det);
+    printf ("\ndetd = %.17g\n", detd);
+  }
+  
+  int sign = mpz_sgn(det);
+  
+  for (int i = 0; i < 3; i++) {
+    mpz_clear(pa[i]);
+    mpz_clear(pb[i]);
+    mpz_clear(pc[i]);
+    mpz_clear(pd[i]);
+  }
+  
+  mpz_clear(adx); mpz_clear(bdx); mpz_clear(cdx);
+  mpz_clear(ady); mpz_clear(bdy); mpz_clear(cdy);
+  mpz_clear(adz); mpz_clear(bdz); mpz_clear(cdz);
+  
+  mpz_clear(bdy_cdz);
+  mpz_clear(bdz_cdy);
+  mpz_clear(cdy_adz);
+  mpz_clear(cdz_ady);
+  mpz_clear(ady_bdz);
+  mpz_clear(adz_bdy);
+  mpz_clear(bdy_cdz_bdz_cdy);
+  mpz_clear(cdy_adz_cdz_ady);
+  mpz_clear(ady_bdz_adz_bdy);
+  mpz_clear(adx_);
+  mpz_clear(bdx_);
+  mpz_clear(cdx_);
+  mpz_clear(det);
+  
+  return double(sign);
+}
+
+//============================================================================//
+//                                                                            //
+// insphere_gmp()                                                             //
+//                                                                            //
+//============================================================================//
+
+REAL insphere_gmp(REAL *aa, REAL *bb, REAL *cc, REAL *dd, REAL *ee)
+{
+  mpz_t pa[3], pb[3], pc[3], pd[3], pe[3];
+
+  for (int i = 0; i < 3; i++) {
+    mpz_init_set_d(pa[i], aa[i]*1.e+15);
+    mpz_init_set_d(pb[i], bb[i]*1.e+15);
+    mpz_init_set_d(pc[i], cc[i]*1.e+15);
+    mpz_init_set_d(pd[i], dd[i]*1.e+15);
+    mpz_init_set_d(pe[i], ee[i]*1.e+15);
+  }
+
+  mpz_t aex, bex, cex, dex;
+  mpz_t aey, bey, cey, dey;
+  mpz_t aez, bez, cez, dez;
+
+  //aex = pa[0] - pe[0];
+  mpz_init(aex);
+  mpz_sub(aex, pa[0], pe[0]);
+  //bex = pb[0] - pe[0];
+  mpz_init(bex);
+  mpz_sub(bex, pb[0], pe[0]);
+  //cex = pc[0] - pe[0];
+  mpz_init(cex);
+  mpz_sub(cex, pc[0], pe[0]);
+  //dex = pd[0] - pe[0];
+  mpz_init(dex);
+  mpz_sub(dex, pd[0], pe[0]);
+  //aey = pa[1] - pe[1];
+  mpz_init(aey);
+  mpz_sub(aey, pa[1], pe[1]);
+  //bey = pb[1] - pe[1];
+  mpz_init(bey);
+  mpz_sub(bey, pb[1], pe[1]);
+  //cey = pc[1] - pe[1];
+  mpz_init(cey);
+  mpz_sub(cey, pc[1], pe[1]);
+  //dey = pd[1] - pe[1];
+  mpz_init(dey);
+  mpz_sub(dey, pd[1], pe[1]);
+  //aez = pa[2] - pe[2];
+  mpz_init(aez);
+  mpz_sub(aez, pa[2], pe[2]);
+  //bez = pb[2] - pe[2];
+  mpz_init(bez);
+  mpz_sub(bez, pb[2], pe[2]);
+  //cez = pc[2] - pe[2];
+  mpz_init(cez);
+  mpz_sub(cez, pc[2], pe[2]);
+  //dez = pd[2] - pe[2];
+  mpz_init(dez);
+  mpz_sub(dez, pd[2], pe[2]);
+
+  mpz_t aexbey, bexaey, bexcey, cexbey, cexdey, dexcey, dexaey, aexdey;
+  mpz_t aexcey, cexaey, bexdey, dexbey;
+  mpz_t alift, blift, clift, dlift;
+  mpz_t ab, bc, cd, da, ac, bd;
+  mpz_t abc, bcd, cda, dab;
+  mpz_t det;
+
+  //aexbey = aex * bey;
+  mpz_init(aexbey);
+  mpz_mul(aexbey, aex, bey);
+  
+  //bexaey = bex * aey;
+  mpz_init(bexaey);
+  mpz_mul(bexaey, bex, aey);
+  
+  //ab = aexbey - bexaey;
+  mpz_init(ab);
+  mpz_sub(ab, aexbey, bexaey);
+  
+  //bexcey = bex * cey;
+  mpz_init(bexcey);
+  mpz_mul(bexcey, bex, cey);
+  
+  //cexbey = cex * bey;
+  mpz_init(cexbey);
+  mpz_mul(cexbey, cex, bey);
+  
+  //bc = bexcey - cexbey;
+  mpz_init(bc);
+  mpz_sub(bc, bexcey, cexbey);
+  
+  //cexdey = cex * dey;
+  mpz_init(cexdey);
+  mpz_mul(cexdey, cex, dey);
+  
+  //dexcey = dex * cey;
+  mpz_init(dexcey);
+  mpz_mul(dexcey, dex, cey);
+  
+  //cd = cexdey - dexcey;
+  mpz_init(cd);
+  mpz_sub(cd, cexdey, dexcey);
+  
+  //dexaey = dex * aey;
+  mpz_init(dexaey);
+  mpz_mul(dexaey, dex, aey);
+  
+  //aexdey = aex * dey;
+  mpz_init(aexdey);
+  mpz_mul(aexdey, aex, dey);
+  
+  //da = dexaey - aexdey;
+  mpz_init(da);
+  mpz_sub(da, dexaey, aexdey);
+  
+  //aexcey = aex * cey;
+  mpz_init(aexcey);
+  mpz_mul(aexcey, aex, cey);
+  
+  //cexaey = cex * aey;
+  mpz_init(cexaey);
+  mpz_mul(cexaey, cex, aey);
+  
+  //ac = aexcey - cexaey;
+  mpz_init(ac);
+  mpz_sub(ac, aexcey, cexaey);
+  
+  //bexdey = bex * dey;
+  mpz_init(bexdey);
+  mpz_mul(bexdey, bex, dey);
+  
+  //dexbey = dex * bey;
+  mpz_init(dexbey);
+  mpz_mul(dexbey, dex, bey);
+  
+  //bd = bexdey - dexbey;
+  mpz_init(bd);
+  mpz_sub(bd, bexdey, dexbey);
+  
+  //abc = aez * bc - bez * ac + cez * ab;
+  mpz_t cez_ab;
+  mpz_init(cez_ab);
+  mpz_mul(cez_ab, cez, ab);
+  mpz_t bez_ac;
+  mpz_init(bez_ac);
+  mpz_mul(bez_ac, bez, ac);
+  mpz_t aez_bc;
+  mpz_init(aez_bc);
+  mpz_mul(aez_bc, aez, bc);
+  mpz_init_set(abc, aez_bc);
+  mpz_sub(abc, abc, bez_ac);
+  mpz_add(abc, abc, cez_ab);
+  
+  //bcd = bez * cd - cez * bd + dez * bc;
+  mpz_t bez_cd;
+  mpz_init(bez_cd);
+  mpz_mul(bez_cd, bez, cd);
+  mpz_t cez_bd;
+  mpz_init(cez_bd);
+  mpz_mul(cez_bd, cez, bd);
+  mpz_t dez_bc;
+  mpz_init(dez_bc);
+  mpz_mul(dez_bc, dez, bc);
+  mpz_init_set(bcd, bez_cd);
+  mpz_sub(bcd, bcd, cez_bd);
+  mpz_add(bcd, bcd, dez_bc);
+  
+  //cda = cez * da + dez * ac + aez * cd;
+  mpz_t cez_da;
+  mpz_init(cez_da);
+  mpz_mul(cez_da, cez, da);
+  mpz_t dez_ac;
+  mpz_init(dez_ac);
+  mpz_mul(dez_ac, dez, ac);
+  mpz_t aez_cd;
+  mpz_init(aez_cd);
+  mpz_mul(aez_cd, aez, cd);
+  mpz_init_set(cda, cez_da);
+  mpz_add(cda, cda, dez_ac);
+  mpz_add(cda, cda, aez_cd);
+  
+  //dab = dez * ab + aez * bd + bez * da;
+  mpz_t dez_ab;
+  mpz_init(dez_ab);
+  mpz_mul(dez_ab, dez, ab);
+  mpz_t aez_bd;
+  mpz_init(aez_bd);
+  mpz_mul(aez_bd, aez, bd);
+  mpz_t bez_da;
+  mpz_init(bez_da);
+  mpz_mul(bez_da, bez, da);
+  mpz_init_set(dab, dez_ab);
+  mpz_add(dab, dab, aez_bd);
+  mpz_add(dab, dab, bez_da);
+  
+  //alift = aex * aex + aey * aey + aez * aez;
+  mpz_t aex_aex;
+  mpz_init(aex_aex);
+  mpz_mul(aex_aex, aex, aex);
+  mpz_t aey_aey;
+  mpz_init(aey_aey);
+  mpz_mul(aey_aey, aey, aey);
+  mpz_t aez_aez;
+  mpz_init(aez_aez);
+  mpz_mul(aez_aez, aez, aez);
+  mpz_init_set(alift, aex_aex);
+  mpz_add(alift, alift, aey_aey);
+  mpz_add(alift, alift, aez_aez);
+  
+  //blift = bex * bex + bey * bey + bez * bez;
+  mpz_t bex_bex;
+  mpz_init(bex_bex);
+  mpz_mul(bex_bex, bex, bex);
+  mpz_t bey_bey;
+  mpz_init(bey_bey);
+  mpz_mul(bey_bey, bey, bey);
+  mpz_t bez_bez;
+  mpz_init(bez_bez);
+  mpz_mul(bez_bez, bez, bez);
+  mpz_init_set(blift, bex_bex);
+  mpz_add(blift, blift, bey_bey);
+  mpz_add(blift, blift, bez_bez);
+  
+  //clift = cex * cex + cey * cey + cez * cez;
+  mpz_t cex_cex;
+  mpz_init(cex_cex);
+  mpz_mul(cex_cex, cex, cex);
+  mpz_t cey_cey;
+  mpz_init(cey_cey);
+  mpz_mul(cey_cey, cey, cey);
+  mpz_t cez_cez;
+  mpz_init(cez_cez);
+  mpz_mul(cez_cez, cez, cez);
+  mpz_init_set(clift, cex_cex);
+  mpz_add(clift, clift, cey_cey);
+  mpz_add(clift, clift, cez_cez);
+  
+  //dlift = dex * dex + dey * dey + dez * dez;
+  mpz_t dex_dex;
+  mpz_init(dex_dex);
+  mpz_mul(dex_dex, dex, dex);
+  mpz_t dey_dey;
+  mpz_init(dey_dey);
+  mpz_mul(dey_dey, dey, dey);
+  mpz_t dez_dez;
+  mpz_init(dez_dez);
+  mpz_mul(dez_dez, dez, dez);
+  mpz_init_set(dlift, dex_dex);
+  mpz_add(dlift, dlift, dey_dey);
+  mpz_add(dlift, dlift, dez_dez);
+  
+  //det = (dlift * abc - clift * dab) + (blift * cda - alift * bcd);
+  mpz_t dlift_abc;
+  mpz_init(dlift_abc);
+  mpz_mul(dlift_abc, dlift, abc);
+  mpz_t clift_dab;
+  mpz_init(clift_dab);
+  mpz_mul(clift_dab, clift, dab);
+  mpz_t blift_cda;
+  mpz_init(blift_cda);
+  mpz_mul(blift_cda, blift, cda);
+  mpz_t alift_bcd;
+  mpz_init(alift_bcd);
+  mpz_mul(alift_bcd, alift, bcd);
+  mpz_init_set(det, dlift_abc);
+  mpz_sub(det, det, clift_dab);
+  mpz_add(det, det, blift_cda);
+  mpz_sub(det, det, alift_bcd);
+  
+  bool debug_flag = false;
+  if (debug_flag) { // Debug only
+    char str[1024];
+    mpz_get_str(str, 10, det);
+    printf("\ndet_str = %s\n", str);
+    double detd = mpz_get_d(det);
+    printf ("\ndetd = %.17g\n", detd);
+  }
+  
+  int sign = mpz_sgn(det);
+  
+  // Clear memory
+  for (int i = 0; i < 3; i++) {
+    mpz_clear(pa[i]);
+    mpz_clear(pb[i]);
+    mpz_clear(pc[i]);
+    mpz_clear(pd[i]);
+    mpz_clear(pe[i]);
+  }
+  
+  mpz_clear(aex);
+  mpz_clear(bex);
+  mpz_clear(cex);
+  mpz_clear(dex);
+  mpz_clear(aey);
+  mpz_clear(bey);
+  mpz_clear(cey);
+  mpz_clear(dey);
+  mpz_clear(aez);
+  mpz_clear(bez);
+  mpz_clear(cez);
+  mpz_clear(dez);
+  
+  mpz_clear(aexbey);
+  mpz_clear(bexaey);
+  mpz_clear(bexcey);
+  mpz_clear(cexbey);
+  mpz_clear(cexdey);
+  mpz_clear(dexcey);
+  mpz_clear(dexaey);
+  mpz_clear(aexdey);
+  mpz_clear(aexcey);
+  mpz_clear(cexaey);
+  mpz_clear(bexdey);
+  mpz_clear(dexbey);
+  mpz_clear(alift);
+  mpz_clear(blift);
+  mpz_clear(clift);
+  mpz_clear(dlift);
+  mpz_clear(ab);
+  mpz_clear(bc);
+  mpz_clear(cd);
+  mpz_clear(da);
+  mpz_clear(ac);
+  mpz_clear(bd);
+  mpz_clear(abc);
+  mpz_clear(bcd);
+  mpz_clear(cda);
+  mpz_clear(dab);
+  mpz_clear(det);
+  
+  mpz_clear(cez_ab);
+  mpz_clear(bez_ac);
+  mpz_clear(aez_bc);
+  
+  mpz_clear(bez_cd);
+  mpz_clear(cez_bd);
+  mpz_clear(dez_bc);
+  
+  mpz_clear(cez_da);
+  mpz_clear(dez_ac);
+  mpz_clear(aez_cd);
+  
+  mpz_clear(dez_ab);
+  mpz_clear(aez_bd);
+  mpz_clear(bez_da);
+  
+  mpz_clear(aex_aex);
+  mpz_clear(aey_aey);
+  mpz_clear(aez_aez);
+  
+  mpz_clear(bex_bex);
+  mpz_clear(bey_bey);
+  mpz_clear(bez_bez);
+  
+  mpz_clear(cex_cex);
+  mpz_clear(cey_cey);
+  mpz_clear(cez_cez);
+  
+  mpz_clear(dex_dex);
+  mpz_clear(dey_dey);
+  mpz_clear(dez_dez);
+  
+  mpz_clear(dlift_abc);
+  mpz_clear(clift_dab);
+  mpz_clear(blift_cda);
+  mpz_clear(alift_bcd);
+  
+  return double(sign);
+}
+
+#endif
 
 //==============================================================================
 

tetgen.h

@@ -1,19 +1,21 @@
 //============================================================================//
-//                                                                            //
-// TetGen                                                                     //
+//  _____     _    _____                                                      //
+// |_   _|__ | |_ / ____|___ _ __                                             //
+//   | |/ _ \| __| |  _ / _ \ '_ \                                            //
+//   | |  __/| |_| |_| |  __/ | | |                                           //
+//   |_|\___/ \__|\_____\___|_| |_|                                           //
 //                                                                            //
 // A Quality Tetrahedral Mesh Generator and A 3D Delaunay Triangulator        //
 //                                                                            //
-// Version 1.6.0                                                              //
-// August 31, 2020                                                            //
+// Version 1.6.1                                                              //
+// xxx xx, 2025                                                               //
 //                                                                            //
-// Copyright (C) 2002--2020                                                   //
+// Copyright (C) 2002--2025                                                   //
 //                                                                            //
 // Hang Si                                                                    //
-// Research Group: Numerical Mathematics and Scientific Computing             //
-// Weierstrass Institute for Applied Analysis and Stochastics (WIAS)          //
-// Mohrenstr. 39, 10117 Berlin, Germany                                       //
-// si@wias-berlin.de                                                          //
+// hangsi@dlut.edu.cn                                                         //
+// http://www.tetgen.org                                                      //
+// https://codeberg.org/TetGen/TetGen                                         //
 //                                                                            //
 // TetGen is a tetrahedral mesh generator. It creates 3d triangulations of    //
 //   polyhedral domains. It generates meshes with well-shaped elements whose  //
@@ -78,6 +80,15 @@
 #include <string.h>
 #include <math.h>
 #include <time.h>
+#include <vector>
+#include <iomanip>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+#ifdef USING_GMP
+#include <gmpxx.h>
+#endif
 
 // The types 'intptr_t' and 'uintptr_t' are signed and unsigned integer types,
 //   respectively. They are guaranteed to be the same width as a pointer.
@@ -155,6 +166,13 @@ public:
   } vorofacet;
 
 
+  typedef struct
+  {
+    int v1, v2, v3;
+    int tag;
+    double maxvol;
+  } Subdomain_Facet;
+
   // Additional parameters associated with an input (or mesh) vertex.
   //   These informations are provided by CAD libraries. 
   typedef struct {
@@ -248,7 +266,10 @@ public:
   REAL *regionlist;
   int numberofregions;
 
-  // 'refine_elem_list': An array of tetrahedra to be refined.  The first 
+  // a list of regions (subdomains) defined by facets.
+  std::vector<Subdomain_Facet> Subdomain_Facets;
+
+  // 'refine_elem_list': An array of tetrahedra to be refined.  The first
   //   tetrahedron's first corner is at index [0], followed by its other 
   //   corners. Four integers per element.
   // 'refine_elem_vol_list':  An array of constraints, i.e. tetrahedron's
@@ -339,6 +360,7 @@ public:
   bool load_edge(char*);
   bool load_face(char*);
   bool load_tet(char*);
+  bool load_region(char*);
   bool load_vol(char*);
   bool load_var(char*);
   bool load_mtr(char*);
@@ -349,7 +371,6 @@ public:
   bool load_stl(char*);
   bool load_vtk(char*);
   bool load_medit(char*, int);
-  bool load_neumesh(char*, int);
   bool load_plc(char*, int);
   bool load_tetmesh(char*, int);
   void save_nodes(const char*);
@@ -389,7 +410,7 @@ public:
     pointattributelist = (REAL *) NULL;
     pointmtrlist = (REAL *) NULL;
     pointmarkerlist = (int *) NULL;
-	point2tetlist = (int *) NULL;
+	  point2tetlist = (int *) NULL;
     pointparamlist = (pointparam *) NULL;
     numberofpoints = 0;
     numberofpointattributes = 0;
@@ -399,8 +420,8 @@ public:
     tetrahedronattributelist = (REAL *) NULL;
     tetrahedronvolumelist = (REAL *) NULL;
     neighborlist = (int *) NULL;
-	tet2facelist = (int *) NULL;
-	tet2edgelist = (int *) NULL;
+	  tet2facelist = (int *) NULL;
+	  tet2edgelist = (int *) NULL;
     numberoftetrahedra = 0;
     numberofcorners = 4; 
     numberoftetrahedronattributes = 0;
@@ -409,8 +430,8 @@ public:
     trifacemarkerlist = (int *) NULL;
     o2facelist = (int *) NULL;
     face2tetlist = (int *) NULL;
-	face2edgelist = (int *) NULL;
-    numberoftrifaces = 0; 
+	  face2edgelist = (int *) NULL;
+    numberoftrifaces = 0;
 
     edgelist = (int *) NULL;
     edgemarkerlist = (int *) NULL;
@@ -476,7 +497,7 @@ public:
     if (pointmarkerlist != (int *) NULL) {
       delete [] pointmarkerlist;
     }
-	if (point2tetlist != (int *) NULL) {
+	  if (point2tetlist != (int *) NULL) {
       delete [] point2tetlist;
     }
     if (pointparamlist != (pointparam *) NULL) {
@@ -496,11 +517,11 @@ public:
       delete [] neighborlist;
     }
     if (tet2facelist != (int *) NULL) {
-	  delete [] tet2facelist;
-	}
-	if (tet2edgelist != (int *) NULL) {
-	  delete [] tet2edgelist;
-	}
+	    delete [] tet2facelist;
+	  }
+	  if (tet2edgelist != (int *) NULL) {
+	    delete [] tet2edgelist;
+	  }
 
     if (trifacelist != (int *) NULL) {
       delete [] trifacelist;
@@ -514,7 +535,7 @@ public:
     if (face2tetlist != (int *) NULL) {
       delete [] face2tetlist;
     }
-	if (face2edgelist != (int *) NULL) {
+	  if (face2edgelist != (int *) NULL) {
       delete [] face2edgelist;
     }
 
@@ -599,13 +620,9 @@ public:
 //                                                                            //
 // tetgenbehavior                                                             //
 //                                                                            //
-// A structure for maintaining the switches and parameters used by TetGen's   //
-// internal data structure and algorithms.                                    //
+// A structure for maintaining TetGen's options and parameters.               //
 //                                                                            //
-// All switches and parameters are initialized with default values. They are  //
-// set by the command line arguments (argc, argv).                            //
-//                                                                            //
-// NOTE: Some switches are incompatible with others. While some may depend    //
+// NOTE: Some options are incompatible with others. While some may depend     //
 // on other switches. The routine parse_commandline() sets the switches from  //
 // the command line (a list of strings) and checks the consistency of the     //
 // applied switches.                                                          //
@@ -616,7 +633,7 @@ class tetgenbehavior {
 
 public:
 
-  // Switches of TetGen. 
+  // Options
   int plc;                                                         // '-p', 0.
   int psc;                                                         // '-s', 0.
   int refine;                                                      // '-r', 0.
@@ -638,13 +655,14 @@ public:
   int nomergefacet;                                                // '-M', 0.
   int nomergevertex;                                               // '-M', 0.
   int noexact;                                                     // '-X', 0.
-  int nostaticfilter;                                              // '-X', 0.
+  int nostaticfilter;                                             // '-X1', 0.
   int zeroindex;                                                   // '-z', 0.
   int facesout;                                                    // '-f', 0.
   int edgesout;                                                    // '-e', 0.
   int neighout;                                                    // '-n', 0.
   int voroout;                                                     // '-v', 0.
   int meditview;                                                   // '-g', 0.
+  int out_mesh_to_nas;                                            // '-g2', 0.
   int vtkview;                                                     // '-k', 0.
   int vtksurfview;                                                 // '-k', 0.
   int nobound;                                                     // '-B', 0.
@@ -658,7 +676,7 @@ public:
   int nowarning;                                                   // '-W', 0.
   int verbose;                                                     // '-V', 0.
 
-  // Parameters of TetGen.
+  // Parameters
   int vertexperblock;                                           // '-x', 4092.
   int tetrahedraperblock;                                       // '-x', 8188.
   int shellfaceperblock;                                        // '-x', 2044.
@@ -669,6 +687,7 @@ public:
   int fliplinklevel;                                                    // -1.
   int flipstarsize;                                                     // -1.
   int fliplinklevelinc;                                                 //  1.
+  int flip_valid_angles;                                                //  0.
   int opt_max_flip_level;                                          // '-O', 3.
   int opt_scheme;                                                // '-O/#', 7.
   int opt_iterations;                                             // -O//#, 3.
@@ -690,16 +709,26 @@ public:
   REAL facet_small_ang_tol;                                   // '-p//', 15.0.
   REAL maxvolume;                                               // '-a', -1.0.
   REAL maxvolume_length;                                        // '-a', -1.0.
-  REAL minratio;                                                 // '-q', 0.0.
+  REAL minratio;                                                // '-q#', 2.0.
+  REAL mindihedral;                                            // '-q/#', 0.0.
+  REAL segment_enc_angle;                             // '-q//#', 90.0 degree.
+  REAL facet_enc_dihed_angle;                        // '-q///#', 90.0 degree.
+  REAL min_edge_length;                                         // '-L#', 0.0.
+  REAL growth_ratio;                                            //  '-L/#', 0.
   REAL opt_max_asp_ratio;                                           // 1000.0.
   REAL opt_max_edge_ratio;                                           // 100.0.
-  REAL mindihedral;                                              // '-q', 5.0.
   REAL optmaxdihedral;                                          // -o/# 177.0.
   REAL metric_scale;                                              // -m#, 1.0.
   REAL smooth_alpha;                                             // '-s', 0.3.
   REAL coarsen_percent;                                         // -R1/#, 1.0.
   REAL elem_growth_ratio;             // Growth ratio of # elements, -r#, 0.0.
-  REAL refine_progress_ratio;                                   // -r/#, 0.333.
+  REAL refine_progress_ratio;                                  // -r/#, 0.333.
+
+  int convert_interior_subdomains_to_holes;                         // -A2, 0.
+
+  int use_volume_size_map;                         // -rS use volume size map.
+  int debug_out_Steiner_tags;  // -GN use out_nodes_Steiner() to output nodes.
+  int debug_volume_size_map;            // -GS Output volume size map to .vtk.
 
   // Strings of command line arguments and input/output file names.
   char commandline[1024];
@@ -708,11 +737,7 @@ public:
   char addinfilename[1024];
   char bgmeshfilename[1024];
 
-  // Read an additional tetrahedral mesh and treat it as holes [2018-07-30].
-  int hole_mesh;                                                   // '-H', 0.
-  char hole_mesh_filename[1024];
-
-  // The input object of TetGen. They are recognized by either the input 
+  // The input object of TetGen. They are recognized by either the input
   //   file extensions or by the specified options. 
   // Currently the following objects are supported:
   //   - NODES, a list of nodes (.node); 
@@ -724,7 +749,7 @@ public:
   //   - MESH, a tetrahedral mesh (.ele).
   // If no extension is available, the imposed command line switch
   //   (-p or -r) implies the object. 
-  enum objecttype {NODES, POLY, OFF, PLY, STL, MEDIT, VTK, MESH, NEU_MESH} object;
+  enum objecttype {NODES, POLY, OFF, PLY, STL, MEDIT, VTK, MESH} object;
 
 
   void syntax();
@@ -765,6 +790,7 @@ public:
     neighout = 0;
     voroout = 0;
     meditview = 0;
+    out_mesh_to_nas = 0;
     vtkview = 0;
     vtksurfview = 0;
     nobound = 0;
@@ -790,6 +816,7 @@ public:
     fliplinklevel = -1; 
     flipstarsize = -1;  
     fliplinklevelinc = 1;
+    flip_valid_angles = 0;
     opt_scheme = 7;
     opt_max_flip_level = 3;
     opt_iterations = 3;
@@ -809,9 +836,13 @@ public:
     maxvolume = -1.0;
     maxvolume_length = -1.0;
     minratio = 2.0;
+    mindihedral = 3.5;
+    segment_enc_angle = 90.0;
+    facet_enc_dihed_angle = 90.0;
+    min_edge_length = 0.0;
+    growth_ratio = 0.;
     opt_max_asp_ratio = 1000.;
     opt_max_edge_ratio = 100.;
-    mindihedral = 3.5;
     optmaxdihedral = 177.00;
     epsilon = 1.0e-8;
     coarsen_percent = 1.0;
@@ -819,7 +850,13 @@ public:
     elem_growth_ratio = 0.0; // -r#
     refine_progress_ratio = 0.333; // -r/#
     object = NODES;
-    
+
+    int convert_interior_subdomains_to_holes = 0; // -A2
+
+    use_volume_size_map = 0; // -rS
+    debug_out_Steiner_tags = 0; // -GN
+    debug_volume_size_map = 0; // -GS
+
     smooth_cirterion = 3; // -s# default smooth surface and volume vertices.
     smooth_maxiter = 7;   // set by -s#/7
     smooth_alpha = 0.3;   // relax parameter, set by -s#/#/0.3
@@ -830,13 +867,11 @@ public:
     addinfilename[0] = '\0';
     bgmeshfilename[0] = '\0';
 
-    hole_mesh = 0;
-    hole_mesh_filename[0] = '\0';
-
   }
 
 }; // class tetgenbehavior
 
+#ifndef USING_GMP
 //============================================================================//
 //                                                                            //
 // Robust Geometric predicates                                                //
@@ -866,16 +901,17 @@ REAL orient3dexact(REAL *pa, REAL *pb, REAL *pc, REAL *pd);
 REAL orient4dexact(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL* pe,
                    REAL ah, REAL bh, REAL ch, REAL dh, REAL eh);
 
+#endif
 
 //============================================================================//
 //                                                                            //
-// tetgenmesh    TetGen's internal mesh data structure.                       //
+// tetgenmesh    TetGen's mesh data structure.                                //
 //                                                                            //
-// It uses a tetrahedron-based mesh data structure. It implements elementary  //
-// flip operations to locally modify the mesh. It implements basic meshing    //
-// algorithms to create Delaunay tetrahedraliations, to perform boundary      //
-// recovery, to place Steiner points in the mesh domain, and to optimize the  //
-// quality of the mesh.                                                       //
+// It implements a tetrahedron-based mesh data structure. It implements algo- //
+// rithms to perform elementary flip operations to locally modify the mesh,   //
+// to create Delaunay tetrahedraliations, to perform boundary recovery, and   //
+// to place Steiner points in the mesh domain, and to optimize the quality of //
+// the mesh.                                                                  //
 //                                                                            //
 //============================================================================//
 
@@ -947,6 +983,7 @@ public:
   //   - a metric tensor (optional, for -q or -m switch);
   //   - a pointer to an adjacent tetrahedron;
   //   - a pointer to a parent (or a duplicate) point;
+  //   - a pointer to the Cpoint (to access additional attributes);
   //   - a pointer to an adjacent subface or segment (optional, -p switch);
   //   - a pointer to a tet in background mesh (optional, for -m switch);
   //   - an integer for boundary marker (point index);
@@ -956,6 +993,23 @@ public:
   //   determined at the runtime.
 
   typedef REAL *point;
+  
+  // A compound point (which includes rational coordinates, etc)
+  class Cpoint { public:
+    REAL *pt;
+    char p[3][64];
+    char q[3][64];
+    Cpoint *next;
+    
+    void init() {
+      pt = nullptr;
+      p[0][0] = p[1][0] = p[2][0] = '\0';
+      q[0][0] = q[1][0] = q[2][0] = '\0';
+      next = nullptr;
+    }
+    
+    Cpoint() {init();}
+  };
 
 //============================================================================//
 //                                                                            //
@@ -1189,6 +1243,9 @@ public:
     REAL smlen; // for useinsertradius.
     point parentpt;
 
+    int check_gr_ratio;
+    double gr_len;
+
     void init() {
       iloc = bowywat = lawson = 0;
       splitbdflag = validflag = respectbdflag = 0;
@@ -1205,6 +1262,9 @@ public:
       smlenflag = 0;
       smlen = 0.0;
       parentpt = NULL;
+      
+      check_gr_ratio = 0;
+      gr_len = 0.;
     }
 
     insertvertexflags() {
@@ -1289,7 +1349,7 @@ public:
 
     // The one of goals of optimization.
     int max_min_volume;      // Maximize the minimum volume.
-	int min_max_aspectratio; // Minimize the maximum aspect ratio. 
+	  int min_max_aspectratio; // Minimize the maximum aspect ratio.
     int min_max_dihedangle;  // Minimize the maximum dihedral angle.
 
     // The initial and improved value.
@@ -1336,12 +1396,12 @@ public:
   // Labels that signify the result of point location.
   enum locateresult {UNKNOWN, OUTSIDE, INTETRAHEDRON, ONFACE, ONEDGE, ONVERTEX,
                      ENCVERTEX, ENCSEGMENT, ENCSUBFACE, NEARVERTEX, NONREGULAR,
-                     INSTAR, BADELEMENT, NULLCAVITY, SHARPCORNER, FENSEDIN,
-                     NONCOPLANAR, SELF_ENCROACH};
+                     INSTAR, BADELEMENT, MINEDGELENGTH, NULLCAVITY, SHARPCORNER,
+                     FENSEDIN, NONCOPLANAR, SELF_ENCROACH};
 
 //============================================================================//
 //                                                                            //
-// Variables of TetGen                                                        //
+// TetGen's internally-used variables                                         //
 //                                                                            //
 //============================================================================//
 
@@ -1358,6 +1418,9 @@ public:
   //   segments) and extra pointers between tetrahedra, subfaces, and segments.
   memorypool *tetrahedrons, *subfaces, *subsegs, *points;
   memorypool *tet2subpool, *tet2segpool;
+  
+  // Memorypool to store compound points.
+  memorypool *Cpointpool;
 
   // Memorypools to store bad-quality (or encroached) elements.
   memorypool *badtetrahedrons, *badsubfacs, *badsubsegs;
@@ -1413,7 +1476,7 @@ public:
   int    segmentendpointslist_length;
   point  *segmentendpointslist;
   double *segment_info_list;
-  int    *idx_segment_ridge_vertex_list; // are two ridge vertices form a segment?
+  int    *idx_segment_ridge_vertex_list; // store two endpoints of a segment
   point  *segment_ridge_vertex_list;
 
   // The infinite vertex.
@@ -1422,6 +1485,11 @@ public:
   triface recenttet;
   face recentsh;
 
+  // Pre-calculated above points.
+  //   required by the tri_edge_2d() test, and incircle3d() test.
+  point _above_points[3];
+  point _above_point;
+  
   // PI is the ratio of a circle's circumference to its diameter.
   static REAL PI;
 
@@ -1429,7 +1497,7 @@ public:
   int subdomains;                                    // Number of subdomains.
   int *subdomain_markers;
 
-  // Various variables.
+  // Variables.
   int numpointattrib;                          // Number of point attributes.
   int numelemattrib;                     // Number of tetrahedron attributes.
   int sizeoftensor;                     // Number of REALs per metric tensor.
@@ -1458,12 +1526,18 @@ public:
   REAL cosslidihed;      // The cosine value of the max dihedral of a sliver.
   REAL cos_large_dihed;   // The cosine value of large dihedral (135 degree).
   REAL opt_max_sliver_asp_ratio;              // = 10 x b->opt_max_asp_ratio.
-  REAL minfaceang, minfacetdihed;     // The minimum input (dihedral) angles.
+  REAL minratio_ang;                 // The angle corresponds to b->minratio.
+  REAL min_input_facet_dihed;            // The minimum input dihedral angle.
   REAL cos_facet_separate_ang_tol;
   REAL cos_collinear_ang_tol;
+  REAL cos_flip_valid_ang_tol;
+  REAL cos_flip_valid_dihed_ang_tol;
+  REAL cos_segment_enc_angle;                  // -q//# b->segment_enc_angle.
+  REAL cos_facet_enc_dihed_angle;         // -q///# b->facet_enc_dihed_angle.
   REAL tetprism_vol_sum;   // The total volume of tetrahedral-prisms (in 4D).
   REAL longest;                          // The longest possible edge length.
   REAL minedgelength;                               // = longest * b->epsion.
+  REAL total_ele_volume;                        // -r sum of all tet volumes.
   REAL xmax, xmin, ymax, ymin, zmax, zmin;         // Bounding box of points.
 
   // Options for mesh refinement.
@@ -1553,7 +1627,6 @@ public:
   inline REAL elemattribute(tetrahedron* ptr, int attnum);
   inline void setelemattribute(tetrahedron* ptr, int attnum, REAL value);
   inline REAL* get_polar(tetrahedron* ptr);
-  inline REAL get_volume(tetrahedron* ptr);
   inline REAL volumebound(tetrahedron* ptr);
   inline void setvolumebound(tetrahedron* ptr, REAL value);
   inline int  elemindex(tetrahedron* ptr);
@@ -1671,6 +1744,8 @@ public:
   inline void setpoint2sh(point pt, shellface value);
   inline point point2ppt(point pt);
   inline void setpoint2ppt(point pt, point value);
+  inline Cpoint* point2Cp(point pt);
+  inline void setpoint2Cp(point pt, Cpoint* value);
   inline tetrahedron point2bgmtet(point pt);
   inline void setpoint2bgmtet(point pt, tetrahedron value);
   inline void setpointinsradius(point pt, REAL value);
@@ -1725,14 +1800,68 @@ public:
 //                                                                            //
 //============================================================================//
 
-  // Symbolic perturbations (robust)
+#ifdef USING_GMP
+  void print_mpq(mpq_class *pa); // debug
+  
+  std::pair<std::string, std::string> double_to_strs(double d, int precision = 15);
+  mpq_class double_to_mpq(double d, int precision = 15);
+  void store_mpq(Cpoint *cp, int i, std::pair<std::string, std::string>&);
+  void store_mpq(Cpoint *cp, int i, double d, int precision = 15);
+  void store_mpq(Cpoint *cp, int i, mpq_class *pa);
+  mpq_class restore_mpq(Cpoint *cp, int i);
+  
+  void clear_mpq_point(point pt); // Recycle the space of compound points.
+  void get_mpq_from_point(point pt, mpq_class *m_pt); // mpq_class m_pt[3], m_pt = pt
+  void set_mpq_to_point(mpq_class *m_pt, point pt); // mpq_class m_pt[3], pt = m_pt;
+  
+  void create_mpq_point(point pt, int precision = 15); // used in transfernodes()
+  void create_mpq_point(double *d, point pt, int precision = 15);
+  void get_d3_from_mpq3(mpq_class *m_pt, double *d); // d[] = m_pt[].get_d()
+  void get_mpq3_from_d3(double *d, mpq_class *m_pt); // m_pt[] = double_to_mpq(d[]).
+  
+  mpq_class orient2d_mpq(mpq_class *pa, mpq_class *pb, mpq_class *pc);
+  mpq_class orient3d_mpq(mpq_class *pa, mpq_class *pb, mpq_class *pc, mpq_class *pd);
+  mpq_class orient3d_mpq(point pa, point pb, point pc, point pd); // Delaunay
+  mpq_class orient3d(point pa, point pb, point pc, mpq_class *m_pd);
+  
+  mpq_class insphere_mpq(mpq_class *pa, mpq_class *pb, mpq_class *pc, mpq_class *pd, mpq_class *pe);
+  mpq_class insphere_mpq(point pa, point pb, point pc, point pd, point pe);
+  
+  mpq_class orient4d_mpq(mpq_class *pa, mpq_class *pb, mpq_class *pc, mpq_class *pd, mpq_class *pe, mpq_class &ah, mpq_class &bh, mpq_class &ch, mpq_class &dh, mpq_class &eh);
+  mpq_class orient4d_mpq(point pa, point pb, point pc, point pd, point pe,
+                         REAL ah, REAL bh, REAL ch, REAL dh, REAL eh);
+  
+  void facenormal_mpq(mpq_class *pa, mpq_class *pb, mpq_class *pc, mpq_class *n);
+  bool planelineint_mpq(mpq_class*, mpq_class*, mpq_class*, mpq_class*, mpq_class*, mpq_class*, mpq_class&);
+  void projpt2edge_mpq(mpq_class *p, mpq_class *e1, mpq_class *e2, mpq_class* prj);
+  void projpt2face_mpq(mpq_class *p, mpq_class *f1, mpq_class *f2, mpq_class *f3,
+                       mpq_class *prj);
+  
+  bool get_subface_ccent_mpq(face *chkfac, mpq_class *m_ccent);
+#endif
+  
+  // Above/below plane predicate (robust)
+  REAL orient3d(point pa, point pb, point pc, point pd);
+  REAL orient4d(point pa, point pb, point pc, point pd, point pe,
+                REAL, REAL, REAL, REAL, REAL);
+  
+  // In-Sphere predicate with symbolic perturbation (robust)
   REAL insphere_s(REAL*, REAL*, REAL*, REAL*, REAL*);
   REAL orient4d_s(REAL*, REAL*, REAL*, REAL*, REAL*, 
                   REAL, REAL, REAL, REAL, REAL);
 
-  // An embedded 2-dimensional geometric predicate (non-robust)
-  REAL incircle3d(point pa, point pb, point pc, point pd);
+  // Degenerated predicates (robust)
+  REAL incircle3d(point pa, point pb, point pc, point pd, point R = NULL);
 
+  // predicates.cxx
+  void pre_calculate_insphere(point, point, point, point, REAL *subdets);
+  REAL insphere_use_subdets(tetrahedron*, REAL*);
+
+  // Set an above point for the plane passing through a set of (co-planar) points.
+  //   Required by the tri_edge_2d() test, and incircle3d() test.
+  int set_above_point3(point, point, point);
+  int set_above_point4(point, point, point, point);
+  
   // Triangle-edge intersection test (robust)
   int tri_edge_2d(point, point, point, point, point, point, int, int*, int*);
   int tri_edge_tail(point,point,point,point,point,point,REAL,REAL,int,int*,int*);
@@ -1755,6 +1884,7 @@ public:
   inline REAL distance2(REAL* p1, REAL* p2);
   void facenormal(point pa, point pb, point pc, REAL *n, int pivot, REAL *lav);
   REAL facedihedral(REAL* pa, REAL* pb, REAL* pc1, REAL* pc2);
+  REAL cos_facedihedral(REAL* pa, REAL* pb, REAL* pc1, REAL* pc2);
   REAL triarea(REAL* pa, REAL* pb, REAL* pc);
   REAL interiorangle(REAL* o, REAL* p1, REAL* p2, REAL* n);
   REAL cos_interiorangle(REAL* o, REAL* p1, REAL* p2);
@@ -1765,8 +1895,7 @@ public:
   void planelineint(REAL*, REAL*, REAL*, REAL*, REAL*, REAL*, REAL*);
   int  linelineint(REAL*, REAL*, REAL*, REAL*, REAL*, REAL*, REAL*, REAL*);
   REAL tetprismvol(REAL* pa, REAL* pb, REAL* pc, REAL* pd);
-  bool calculateabovepoint(arraypool*, point*, point*, point*);
-  void calculateabovepoint4(point, point, point, point);
+
 
 //============================================================================//
 //                                                                            //
@@ -1913,6 +2042,7 @@ public:
   enum locateresult slocate(point, face*, int, int, int);
   enum interresult sscoutsegment(face*, point, int, int, int);
   void scarveholes(int, REAL*);
+  int get_first_triangle(arraypool*, point*, point*, point*);
   int triangulate(int, arraypool*, arraypool*, int, REAL*);
 
   void unifysegments();
@@ -2027,28 +2157,29 @@ public:
 
   bool is_collinear_at(point mid, point left, point right);
   bool is_segment(point p1, point p2);
-  bool valid_constrained_f23(triface&, point pd, point pe);
+  bool valid_constrained_f23(triface& flipface, point pd, point pe, int f44);
   bool valid_constrained_f32(triface*, point pa, point pb);
-  
   int checkflipeligibility(int fliptype, point, point, point, point, point,
                            int level, int edgepivot, flipconstraints* fc);
 
   int removeedgebyflips(triface*, flipconstraints*);
   int removefacebyflips(triface*, flipconstraints*);
 
-  int recoveredgebyflips(point, point, face*, triface*, int fullsearch, int& idir);
+  int recoveredgebyflips(point, point, face*, triface*, int fullsearch,
+                         flipconstraints &fc, int& idir);
   int add_steinerpt_in_schoenhardtpoly(triface*, int, int, int chkencflag);
   int add_steinerpt_in_segment(face*, int searchlevel, int& idir); 
   int add_steinerpt_to_recover_edge(point, point, face*, int, int, int& idir);
   int recoversegments(arraypool*, int fullsearch, int steinerflag);
 
-  int recoverfacebyflips(point,point,point,face*,triface*,int&,point*,point*);
+  int recoverfacebyflips(point, point, point, face*, triface*,
+                        flipconstraints&, int&,point*,point*);
   int recoversubfaces(arraypool*, int steinerflag);
 
   int getvertexstar(int, point searchpt, arraypool*, arraypool*, arraypool*);
   int getedge(point, point, triface*);
-  int reduceedgesatvertex(point startpt, arraypool* endptlist);
-  int removevertexbyflips(point steinerpt);
+  int reduceedgesatvertex(point startpt, arraypool* endptlist, flipconstraints &fc);
+  int removevertexbyflips(point steinerpt, flipconstraints &fc);
 
   int smoothpoint(point smtpt, arraypool*, int ccw, optparameters *opm);
   int suppressbdrysteinerpoint(point steinerpt);
@@ -2114,6 +2245,9 @@ public:
 //                                                                            //
 //============================================================================//
 
+  bool compute_field_point(triface *chktet, REAL* offcent);
+  int  insert_field_points();
+
   void makesegmentendpointsmap();
   REAL set_ridge_vertex_protecting_ball(point);
   REAL get_min_angle_at_ridge_vertex(face* seg);
@@ -2134,8 +2268,7 @@ public:
 
   void enqueuesubface(memorypool*, face*);
   void enqueuetetrahedron(triface*);
-
-  bool check_encroachment(point pa, point pb, point checkpt);
+  
   bool check_enc_segment(face *chkseg, point *pencpt);  
   bool get_steiner_on_segment(face* seg, point encpt, point newpt);
   bool split_segment(face *splitseg, point encpt, REAL *param, int qflag, int, int*);
@@ -2151,9 +2284,10 @@ public:
   enum locateresult locate_on_surface(point searchpt, face* searchsh);
   bool split_subface(face *splitfac, point encpt, REAL *ccent, REAL*, int, int, int*);
   void repairencfacs(REAL *param, int qflag, int chkencflag);
-
-  bool check_tetrahedron(triface *chktet, REAL* param, int& qflag);
+  
   bool checktet4split(triface *chktet, REAL* param, int& qflag);
+  REAL get_longest_edge(triface *chktet, triface *longedge);
+  REAL search_terminal_edge(triface *chktet, triface *termedge, face *termsh, face *termseg, int maxlevel);
   enum locateresult locate_point_walk(point searchpt, triface*, int chkencflag);
   bool split_tetrahedron(triface*, REAL*, int, int, insertvertexflags &ivf);
   void repairbadtets(REAL queratio, int chkencflag);
@@ -2167,7 +2301,7 @@ public:
 //============================================================================//
 
   long lawsonflip3d(flipconstraints *fc);
-  void recoverdelaunay();
+  void recoverdelaunay(flipconstraints &fc);
 
   int  get_seg_laplacian_center(point mesh_vert, REAL target[3]);
   int  get_surf_laplacian_center(point mesh_vert, REAL target[3]);
@@ -2182,10 +2316,17 @@ public:
   badface* top_badtet();
   void dequeue_badtet();
 
-  bool add_steinerpt_to_repair(badface *bf, bool bSmooth);
-  bool flip_edge_to_improve(triface *sliver_edge, REAL& improved_cosmaxd);
-  bool repair_tet(badface *bf, bool bFlips, bool bSmooth, bool bSteiners);
-  long repair_badqual_tets(bool bFlips, bool bSmooth, bool bSteiners);
+  bool construct_point_from_edge(REAL* e1, REAL* e2, REAL* mov_vec, REAL t, REAL angle, REAL* newpt);
+  bool move_vertex_to_improve(point mesh_vert, REAL* target_vert, REAL* move_dir, REAL in_asp, REAL in_cosmaxd);
+  bool add_steinerpt_to_remove_edge(triface *sliver_tet, triface *short_edge, REAL in_asp, REAL in_cosmaxd);
+  bool smooth_flat_S_tet(triface *sliver_edge, REAL in_asp, REAL in_cosmaxd);
+  bool is_edge_collapsible(triface *check_edge, REAL* lambda);
+  bool collapse_edge_to_improve(triface *short_edge, REAL in_asp, REAL in_cosmaxd);
+  bool flip_edge_to_improve(triface *sliver_edge, REAL in_cosmaxd);
+  bool flip_face_to_improve(triface *flip_face, REAL in_cosmaxd);
+  void get_flat_T_tet_shape(badface *bf, REAL, REAL, triface*, triface*, triface*);
+  bool repair_tet(badface *bf, bool bFlips, bool bCollapse, bool bSteiners, bool bSmooth);
+  long repair_badqual_tets(bool bFlips, bool bCollapse, bool bSteiners, bool bSmooth, bool, REAL, REAL);
   void improve_mesh();
 
 //============================================================================//
@@ -2230,9 +2371,11 @@ public:
   void outvoronoi(tetgenio*);
   void outsmesh(char*);
   void outmesh2medit(char*);
+  void out_mesh_to_nas();
   void outmesh2vtk(char*, int);
   void out_surfmesh_vtk(char*, int);
   void out_intersected_facets();
+  void out_nodes_Steiner_tags();
 
 
 
@@ -2252,6 +2395,9 @@ public:
     tetrahedrons = subfaces = subsegs = points = NULL;
     tet2segpool = tet2subpool = NULL;
     dummypoint = NULL;
+    
+    _above_points[0] = _above_points[1] = _above_points[2] = NULL;
+    _above_point = NULL;
 
     badtetrahedrons = badsubfacs = badsubsegs = NULL;
     split_segments_pool = split_subfaces_pool = NULL;
@@ -2259,6 +2405,8 @@ public:
     check_tets_list = NULL;
     badqual_tets_pool = NULL;
 
+    Cpointpool = NULL;
+    
     stack_enc_segments = stack_enc_subfaces = NULL;
   
     flippool = NULL;
@@ -2314,11 +2462,17 @@ public:
     useinsertradius = 0;
     samples = 0l;
     randomseed = 1l;
-    minfaceang = minfacetdihed = PI;
+    minratio_ang = 0;
+    min_input_facet_dihed = PI;
     cos_facet_separate_ang_tol = cos(179.9/180.*PI);
     cos_collinear_ang_tol = cos(179.9/180.*PI);
+    cos_flip_valid_ang_tol = cos(179.97/180.*PI);
+    cos_flip_valid_dihed_ang_tol = cos(179.97/180.*PI);
+    cos_segment_enc_angle = cos(90.0/180.0*PI);
+    cos_facet_enc_dihed_angle = cos(90.0/180.0*PI);
     tetprism_vol_sum = 0.0;
     longest = minedgelength = 0.0;
+    total_ele_volume = 0.;
     xmax = xmin = ymax = ymin = zmax = zmin = 0.0;
 
     smallest_insradius = 1.e+30;
@@ -2357,6 +2511,9 @@ public:
     if (points != (memorypool *) NULL) {
       delete points;
       delete [] dummypoint;
+      for (int i = 0; i < 3; i++) {
+        delete [] _above_points[i];
+      }
     }
     if (tetrahedrons != (memorypool *) NULL) {
       delete tetrahedrons;
@@ -2438,7 +2595,10 @@ public:
     if (subdomain_markers != NULL) {
       delete [] subdomain_markers;
     }
-
+#ifdef USING_GMP
+    delete Cpointpool;
+#endif
+    
     initializetetgenmesh();
   }
 
@@ -2494,7 +2654,7 @@ inline void terminatetetgen(tetgenmesh *m, int x)
     printf("Error:  Out of memory.\n"); 
     break;
   case 2: // Encounter an internal error.
-    printf("Please report this bug to Hang.Si@wias-berlin.de. Include\n");
+    printf("Please report this bug to hangsi@dlut.edu.cn. Include\n");
     printf("  the message above, your input data set, and the exact\n");
     printf("  command line you used to run this program, thank you.\n");
     break;
@@ -2731,10 +2891,7 @@ inline REAL* tetgenmesh::get_polar(tetrahedron* ptr)
 {
   return &(((REAL *) (ptr))[polarindex]);
 }
-inline REAL tetgenmesh::get_volume(tetrahedron* ptr)
-{
-  return ((REAL *) (ptr))[polarindex + 4];
-}
+
 
 // Check or set a tetrahedron's attributes.
 
@@ -3472,21 +3629,29 @@ inline void tetgenmesh::setpoint2ppt(point pt, point value) {
   ((tetrahedron *) (pt))[point2simindex + 1] = (tetrahedron) value;
 }
 
+inline tetgenmesh::Cpoint* tetgenmesh::point2Cp(point pt) {
+  return (tetgenmesh::Cpoint*) ((tetrahedron *) (pt))[point2simindex + 2];
+}
+
+inline void tetgenmesh::setpoint2Cp(point pt, Cpoint* value) {
+  ((tetrahedron *) (pt))[point2simindex + 2] = (tetrahedron) value;
+}
+
 inline tetgenmesh::shellface tetgenmesh::point2sh(point pt) {
-  return (shellface) ((tetrahedron *) (pt))[point2simindex + 2];
+  return (shellface) ((tetrahedron *) (pt))[point2simindex + 3];
 }
 
 inline void tetgenmesh::setpoint2sh(point pt, shellface value) {
-  ((tetrahedron *) (pt))[point2simindex + 2] = (tetrahedron) value;
+  ((tetrahedron *) (pt))[point2simindex + 3] = (tetrahedron) value;
 }
 
 
 inline tetgenmesh::tetrahedron tetgenmesh::point2bgmtet(point pt) {
-  return ((tetrahedron *) (pt))[point2simindex + 3];
+  return ((tetrahedron *) (pt))[point2simindex + 4];
 }
 
 inline void tetgenmesh::setpoint2bgmtet(point pt, tetrahedron value) {
-  ((tetrahedron *) (pt))[point2simindex + 3] = value;
+  ((tetrahedron *) (pt))[point2simindex + 4] = value;
 }