diff --git a/CMakeLists.txt b/CMakeLists.txt index bf44985a111f724d490149c7671c0cd1ee6dfff4..63c7263b358e83e5a07918df9dba3a04e56dbaf8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,7 +19,7 @@ project(gmsh CXX C) # this variable controls the default value of the "ENABLE_XXX" options which are # normally set to ON (useful if you want to configure a minimal version of Gmsh: # e.g. "cmake -DDEFAULT=0 -DENABLE_POST=1 -DENABLE_PARSER=1") -set(DEFAULT ON CACHE INTERNAL "Default value for options") +set(DEFAULT ON CACHE INTERNAL "Default value for enabled-by-default options") option(ENABLE_ACIS "Enable ACIS geometrical models" ${DEFAULT}) option(ENABLE_ANN "Enable ANN to compute Approximate Nearest Neighbors" ${DEFAULT}) @@ -61,7 +61,7 @@ option(ENABLE_SLEPC "Enable SLEPc eigensolvers" ${DEFAULT}) option(ENABLE_SOLVER "Enable solver components" ${DEFAULT}) option(ENABLE_TAUCS "Enable Taucs linear algebra solver" ${DEFAULT}) option(ENABLE_TETGEN "Enable Tetgen mesh generator" ${DEFAULT}) -option(ENABLE_TETGEN_NEW "Enable experimental version of Tetgen" OFF) +option(ENABLE_TETGEN_OLD "Use old version of Tetgen" OFF) option(ENABLE_VORO3D "Enable Voro3D" ${DEFAULT}) option(ENABLE_WRAP_JAVA "Build Java wrappers" OFF) option(ENABLE_WRAP_PYTHON "Build Python wrappers" ${DEFAULT}) @@ -95,7 +95,8 @@ set(GMSH_API Geo/SOrientedBoundingBox.h Geo/CellComplex.h Geo/ChainComplex.h Geo/Cell.h Geo/Homology.h Geo/partitionEdge.h Geo/CGNSOptions.h Geo/gmshLevelset.h Mesh/meshGEdge.h Mesh/meshGFace.h Mesh/meshGFaceOptimize.h Mesh/meshPartition.h - Mesh/meshGFaceDelaunayInsertion.h Mesh/simple3D.h Mesh/meshPartitionOptions.h Mesh/Voronoi3D.h Mesh/Levy3D.h Mesh/periodical.h + Mesh/meshGFaceDelaunayInsertion.h Mesh/simple3D.h Mesh/meshPartitionOptions.h + Mesh/Voronoi3D.h Mesh/Levy3D.h Mesh/periodical.h Numeric/mathEvaluator.h Solver/dofManager.h Solver/femTerm.h Solver/laplaceTerm.h Solver/elasticityTerm.h Solver/crossConfTerm.h Solver/orthogonalTerm.h @@ -141,6 +142,7 @@ endif(APPLE) include(CheckTypeSize) include(CheckFunctionExists) include(CheckIncludeFile) +include(CheckCXXCompilerFlag) if(MSVC) # remove annoying warning about bool/int cast performance @@ -594,17 +596,15 @@ if(HAVE_MESH) set_config_option(HAVE_MMG3D "Mmg3d") endif(ENABLE_MMG3D) - if(ENABLE_TETGEN_NEW AND EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/contrib/TetgenNew/tetgen.h) - add_subdirectory(contrib/TetgenNew) - include_directories(contrib/TetgenNew) - set_config_option(HAVE_TETGEN "Tetgen(New)") + if(ENABLE_TETGEN AND EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/contrib/Tetgen1.5/tetgen.h) + add_subdirectory(contrib/Tetgen1.5) + include_directories(contrib/Tetgen1.5) + set_config_option(HAVE_TETGEN "Tetgen(1.5)") add_definitions(-DTETLIBRARY) - message("WARNING: You are including an experimental version of Tetgen " - "that is KNOWN TO BE BUGGY on 64 bits archs and on WIN32/MSVC.") - elseif(ENABLE_TETGEN AND EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/contrib/Tetgen/tetgen.h) - add_subdirectory(contrib/Tetgen) - include_directories(contrib/Tetgen) - set_config_option(HAVE_TETGEN "Tetgen") + elseif(ENABLE_TETGEN_OLD AND EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/contrib/Tetgen1.4/tetgen.h) + add_subdirectory(contrib/Tetgen1.4) + include_directories(contrib/Tetgen1.4) + set_config_option(HAVE_TETGEN "Tetgen(1.4)") add_definitions(-DTETLIBRARY) elseif(ENABLE_TETGEN) find_library(TETGEN_LIB tet PATH_SUFFIXES lib) @@ -614,7 +614,7 @@ if(HAVE_MESH) list(APPEND EXTERNAL_INCLUDES ${TETGEN_INC}) set_config_option(HAVE_TETGEN "Tetgen") endif(TETGEN_LIB AND TETGEN_INC) - endif(ENABLE_TETGEN_NEW AND EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/contrib/TetgenNew/tetgen.h) + endif(ENABLE_TETGEN AND EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/contrib/Tetgen1.5/tetgen.h) if(HAVE_TETGEN) message("WARNING: By including Tetgen you have to comply with Tetgen's " "special licensing requirements stated in contrib/Tetgen/LICENSE.") @@ -909,13 +909,6 @@ if(DLFCN_H) list(APPEND EXTERNAL_LIBRARIES ${CMAKE_DL_LIBS}) endif(DLFCN_H) -if(UNIX) - # do not optimize some files on Unix - file(GLOB_RECURSE NON_OPTIMIZED_SRC Common/GmshPredicates.cpp Mesh/BDS.cpp - Parser/Gmsh.tab.cpp contrib/Tetgen*/*.cxx) - set_source_files_properties(${NON_OPTIMIZED_SRC} COMPILE_FLAGS "-O0") -endif(UNIX) - if(MSVC) add_definitions(-D_USE_MATH_DEFINES -DNOMINMAX -D_CRT_SECURE_NO_DEPRECATE -D_SCL_SECURE_NO_DEPRECATE) @@ -926,6 +919,13 @@ if(WIN32 OR CYGWIN) list(APPEND EXTERNAL_LIBRARIES wsock32 ws2_32) endif(WIN32 OR CYGWIN) +check_cxx_compiler_flag("-O0" NOOPT) +if(NOOPT) + file(GLOB_RECURSE NON_OPTIMIZED_SRC Numeric/robustPredicates.cpp Mesh/BDS.cpp + Parser/Gmsh.tab.cpp contrib/Tetgen*/*.cxx) + set_source_files_properties(${NON_OPTIMIZED_SRC} COMPILE_FLAGS "-O0") +endif(NOOPT) + list(SORT CONFIG_OPTIONS) set(GMSH_CONFIG_OPTIONS "") foreach(OPT ${CONFIG_OPTIONS}) @@ -1054,7 +1054,6 @@ elseif(MSVC) endif(WIN32 AND NOT MSVC OR CYGWIN) # force full warnings to encourage everybody to write clean(er) code -include(CheckCXXCompilerFlag) check_cxx_compiler_flag("-Wall" WALL) if(WALL) set_target_properties(gmsh lib shared PROPERTIES COMPILE_FLAGS "-Wall") @@ -1235,7 +1234,7 @@ set(CPACK_SOURCE_GENERATOR TGZ) set(CPACK_SOURCE_IGNORE_FILES "${CMAKE_CURRENT_BINARY_DIR}" "/CVS/" "/.svn" "~$" "DS_Store$" "GmshConfig.h$" "GmshVersion.h$" "/benchmarks/" "/tmp/" "/bin/" "/lib/" "/nightly/" "GPATH" "GRTAGS" "GSYMS" "GTAGS" "/HTML/" - "/projects/" "/TetgenNew/") #"/Tetgen.*/.*(cxx|h)") + "/projects/") if(UNIX) # make sure we remove previous installs before doing the next one diff --git a/Numeric/robustPredicates.cpp b/Numeric/robustPredicates.cpp index a861a3094e15308f253965bf9916370534dd96db..dbf750bf557f6c4fbef09be1c377cef542d169cf 100644 --- a/Numeric/robustPredicates.cpp +++ b/Numeric/robustPredicates.cpp @@ -138,7 +138,7 @@ namespace robustPredicates #define INEXACT /* Nothing */ /* #define INEXACT volatile */ -#define REAL double /* float or double */ +#define REAL double /* float or double */ #define REALPRINT doubleprint #define REALRAND doublerand #define NARROWRAND narrowdoublerand @@ -1870,6 +1870,17 @@ REAL orient3dadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL permanent) REAL fin1[192], fin2[192]; int finlength; + //////////////////////////////////////////////////////// + // To avoid uninitialized warnings reported by valgrind. + int i; + for (i = 0; i < 8; i++) { + adet[i] = bdet[i] = cdet[i] = 0.0; + } + for (i = 0; i < 16; i++) { + abdet[i] = 0.0; + } + //////////////////////////////////////////////////////// + REAL adxtail, bdxtail, cdxtail; REAL adytail, bdytail, cdytail; REAL adztail, bdztail, cdztail; @@ -2253,6 +2264,31 @@ REAL orient3dadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL permanent) return finnow[finlength - 1]; } +#ifdef INEXACT_GEOM_PRED + +REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd) +{ + REAL adx, bdx, cdx; + REAL ady, bdy, cdy; + REAL adz, bdz, cdz; + + adx = pa[0] - pd[0]; + bdx = pb[0] - pd[0]; + cdx = pc[0] - pd[0]; + ady = pa[1] - pd[1]; + bdy = pb[1] - pd[1]; + cdy = pc[1] - pd[1]; + adz = pa[2] - pd[2]; + bdz = pb[2] - pd[2]; + cdz = pc[2] - pd[2]; + + return adx * (bdy * cdz - bdz * cdy) + + bdx * (cdy * adz - cdz * ady) + + cdx * (ady * bdz - adz * bdy); +} + +#else + REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd) { REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz; @@ -2279,7 +2315,7 @@ REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd) adxbdy = adx * bdy; bdxady = bdx * ady; - det = adz * (bdxcdy - cdxbdy) + det = adz * (bdxcdy - cdxbdy) + bdz * (cdxady - adxcdy) + cdz * (adxbdy - bdxady); @@ -2294,6 +2330,8 @@ REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd) return orient3dadapt(pa, pb, pc, pd, permanent); } +#endif // ifdef INEXACT_GEOM_PRED + /*****************************************************************************/ /* */ /* incirclefast() Approximate 2D incircle test. Nonrobust. */ @@ -2638,7 +2676,7 @@ REAL incircleadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL permanent) REAL cxtaa[8], cxtbb[8], cytaa[8], cytbb[8]; int cxtaalen, cxtbblen, cytaalen, cytbblen; REAL axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8]; - int axtbclen=0, aytbclen=0, bxtcalen=0, bytcalen=0, cxtablen=0, cytablen=0; + int axtbclen, aytbclen, bxtcalen, bytcalen, cxtablen, cytablen; REAL axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16]; int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen; REAL axtbctt[8], aytbctt[8], bxtcatt[8]; @@ -2660,6 +2698,9 @@ REAL incircleadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL permanent) INEXACT REAL _i, _j; REAL _0; + // Avoid compiler warnings. H. Si, 2012-02-16. + axtbclen = aytbclen = bxtcalen = bytcalen = cxtablen = cytablen = 0; + adx = (REAL) (pa[0] - pd[0]); bdx = (REAL) (pb[0] - pd[0]); cdx = (REAL) (pc[0] - pd[0]); @@ -4074,6 +4115,52 @@ REAL insphereadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe, return insphereexact(pa, pb, pc, pd, pe); } +#ifdef INEXACT_GEOM_PRED + +REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe) +{ + REAL aex, bex, cex, dex; + REAL aey, bey, cey, dey; + REAL aez, bez, cez, dez; + REAL alift, blift, clift, dlift; + REAL ab, bc, cd, da, ac, bd; + REAL abc, bcd, cda, dab; + + aex = pa[0] - pe[0]; + bex = pb[0] - pe[0]; + cex = pc[0] - pe[0]; + dex = pd[0] - pe[0]; + aey = pa[1] - pe[1]; + bey = pb[1] - pe[1]; + cey = pc[1] - pe[1]; + dey = pd[1] - pe[1]; + aez = pa[2] - pe[2]; + bez = pb[2] - pe[2]; + cez = pc[2] - pe[2]; + dez = pd[2] - pe[2]; + + ab = aex * bey - bex * aey; + bc = bex * cey - cex * bey; + cd = cex * dey - dex * cey; + da = dex * aey - aex * dey; + + ac = aex * cey - cex * aey; + bd = bex * dey - dex * bey; + + abc = aez * bc - bez * ac + cez * ab; + bcd = bez * cd - cez * bd + dez * bc; + cda = cez * da + dez * ac + aez * cd; + dab = dez * ab + aez * bd + bez * da; + + alift = aex * aex + aey * aey + aez * aez; + blift = bex * bex + bey * bey + bez * bez; + clift = cex * cex + cey * cey + cez * cez; + dlift = dex * dex + dey * dey + dez * dez; + + return (dlift * abc - clift * dab) + (blift * cda - alift * bcd); +} +#else + REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe) { REAL aex, bex, cex, dex; @@ -4176,4 +4263,547 @@ REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe) return insphereadapt(pa, pb, pc, pd, pe, permanent); } -} // namespace robustPredicates +#endif // #ifdef INEXACT_GEOM_PRED + +/*****************************************************************************/ +/* */ +/* orient4d() Return a positive value if the point pe lies above the */ +/* hyperplane passing through pa, pb, pc, and pd; "above" is */ +/* defined in a manner best found by trial-and-error. Returns */ +/* a negative value if pe lies below the hyperplane. Returns */ +/* zero if the points are co-hyperplanar (not affinely */ +/* independent). The result is also a rough approximation of */ +/* 24 times the signed volume of the 4-simplex defined by the */ +/* five points. */ +/* */ +/* Uses exact arithmetic if necessary to ensure a correct answer. The */ +/* result returned is the determinant of a matrix. This determinant is */ +/* computed adaptively, in the sense that exact arithmetic is used only to */ +/* the degree it is needed to ensure that the returned value has the */ +/* correct sign. Hence, orient4d() is usually quite fast, but will run */ +/* more slowly when the input points are hyper-coplanar or nearly so. */ +/* */ +/* See my Robust Predicates paper for details. */ +/* */ +/*****************************************************************************/ + +REAL orient4dexact(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL* pe, + REAL aheight, REAL bheight, REAL cheight, REAL dheight, + REAL eheight) +{ + INEXACT REAL axby1, bxcy1, cxdy1, dxey1, exay1; + INEXACT REAL bxay1, cxby1, dxcy1, exdy1, axey1; + INEXACT REAL axcy1, bxdy1, cxey1, dxay1, exby1; + INEXACT REAL cxay1, dxby1, excy1, axdy1, bxey1; + REAL axby0, bxcy0, cxdy0, dxey0, exay0; + REAL bxay0, cxby0, dxcy0, exdy0, axey0; + REAL axcy0, bxdy0, cxey0, dxay0, exby0; + REAL cxay0, dxby0, excy0, axdy0, bxey0; + REAL ab[4], bc[4], cd[4], de[4], ea[4]; + REAL ac[4], bd[4], ce[4], da[4], eb[4]; + REAL temp8a[8], temp8b[8], temp16[16]; + int temp8alen, temp8blen, temp16len; + REAL abc[24], bcd[24], cde[24], dea[24], eab[24]; + REAL abd[24], bce[24], cda[24], deb[24], eac[24]; + int abclen, bcdlen, cdelen, dealen, eablen; + int abdlen, bcelen, cdalen, deblen, eaclen; + REAL temp48a[48], temp48b[48]; + int temp48alen, temp48blen; + REAL abcd[96], bcde[96], cdea[96], deab[96], eabc[96]; + int abcdlen, bcdelen, cdealen, deablen, eabclen; + REAL adet[192], bdet[192], cdet[192], ddet[192], edet[192]; + int alen, blen, clen, dlen, elen; + REAL abdet[384], cddet[384], cdedet[576]; + int ablen, cdlen; + REAL deter[960]; + int deterlen; + int i; + + INEXACT REAL bvirt; + REAL avirt, bround, around; + INEXACT REAL c; + INEXACT REAL abig; + REAL ahi, alo, bhi, blo; + REAL err1, err2, err3; + INEXACT REAL _i, _j; + REAL _0; + + Two_Product(pa[0], pb[1], axby1, axby0); + Two_Product(pb[0], pa[1], bxay1, bxay0); + Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]); + + Two_Product(pb[0], pc[1], bxcy1, bxcy0); + Two_Product(pc[0], pb[1], cxby1, cxby0); + Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]); + + Two_Product(pc[0], pd[1], cxdy1, cxdy0); + Two_Product(pd[0], pc[1], dxcy1, dxcy0); + Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]); + + Two_Product(pd[0], pe[1], dxey1, dxey0); + Two_Product(pe[0], pd[1], exdy1, exdy0); + Two_Two_Diff(dxey1, dxey0, exdy1, exdy0, de[3], de[2], de[1], de[0]); + + Two_Product(pe[0], pa[1], exay1, exay0); + Two_Product(pa[0], pe[1], axey1, axey0); + Two_Two_Diff(exay1, exay0, axey1, axey0, ea[3], ea[2], ea[1], ea[0]); + + Two_Product(pa[0], pc[1], axcy1, axcy0); + Two_Product(pc[0], pa[1], cxay1, cxay0); + Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]); + + Two_Product(pb[0], pd[1], bxdy1, bxdy0); + Two_Product(pd[0], pb[1], dxby1, dxby0); + Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]); + + Two_Product(pc[0], pe[1], cxey1, cxey0); + Two_Product(pe[0], pc[1], excy1, excy0); + Two_Two_Diff(cxey1, cxey0, excy1, excy0, ce[3], ce[2], ce[1], ce[0]); + + Two_Product(pd[0], pa[1], dxay1, dxay0); + Two_Product(pa[0], pd[1], axdy1, axdy0); + Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]); + + Two_Product(pe[0], pb[1], exby1, exby0); + Two_Product(pb[0], pe[1], bxey1, bxey0); + Two_Two_Diff(exby1, exby0, bxey1, bxey0, eb[3], eb[2], eb[1], eb[0]); + + temp8alen = scale_expansion_zeroelim(4, bc, pa[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, ac, -pb[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, ab, pc[2], temp8a); + abclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + abc); + + temp8alen = scale_expansion_zeroelim(4, cd, pb[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, bd, -pc[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, bc, pd[2], temp8a); + bcdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + bcd); + + temp8alen = scale_expansion_zeroelim(4, de, pc[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, ce, -pd[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, cd, pe[2], temp8a); + cdelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + cde); + + temp8alen = scale_expansion_zeroelim(4, ea, pd[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, da, -pe[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, de, pa[2], temp8a); + dealen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + dea); + + temp8alen = scale_expansion_zeroelim(4, ab, pe[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, eb, -pa[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, ea, pb[2], temp8a); + eablen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + eab); + + temp8alen = scale_expansion_zeroelim(4, bd, pa[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, da, pb[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, ab, pd[2], temp8a); + abdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + abd); + + temp8alen = scale_expansion_zeroelim(4, ce, pb[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, eb, pc[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, bc, pe[2], temp8a); + bcelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + bce); + + temp8alen = scale_expansion_zeroelim(4, da, pc[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, ac, pd[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, cd, pa[2], temp8a); + cdalen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + cda); + + temp8alen = scale_expansion_zeroelim(4, eb, pd[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, bd, pe[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, de, pb[2], temp8a); + deblen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + deb); + + temp8alen = scale_expansion_zeroelim(4, ac, pe[2], temp8a); + temp8blen = scale_expansion_zeroelim(4, ce, pa[2], temp8b); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, + temp16); + temp8alen = scale_expansion_zeroelim(4, ea, pc[2], temp8a); + eaclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, + eac); + + temp48alen = fast_expansion_sum_zeroelim(cdelen, cde, bcelen, bce, temp48a); + temp48blen = fast_expansion_sum_zeroelim(deblen, deb, bcdlen, bcd, temp48b); + for (i = 0; i < temp48blen; i++) { + temp48b[i] = -temp48b[i]; + } + bcdelen = fast_expansion_sum_zeroelim(temp48alen, temp48a, + temp48blen, temp48b, bcde); + alen = scale_expansion_zeroelim(bcdelen, bcde, aheight, adet); + + temp48alen = fast_expansion_sum_zeroelim(dealen, dea, cdalen, cda, temp48a); + temp48blen = fast_expansion_sum_zeroelim(eaclen, eac, cdelen, cde, temp48b); + for (i = 0; i < temp48blen; i++) { + temp48b[i] = -temp48b[i]; + } + cdealen = fast_expansion_sum_zeroelim(temp48alen, temp48a, + temp48blen, temp48b, cdea); + blen = scale_expansion_zeroelim(cdealen, cdea, bheight, bdet); + + temp48alen = fast_expansion_sum_zeroelim(eablen, eab, deblen, deb, temp48a); + temp48blen = fast_expansion_sum_zeroelim(abdlen, abd, dealen, dea, temp48b); + for (i = 0; i < temp48blen; i++) { + temp48b[i] = -temp48b[i]; + } + deablen = fast_expansion_sum_zeroelim(temp48alen, temp48a, + temp48blen, temp48b, deab); + clen = scale_expansion_zeroelim(deablen, deab, cheight, cdet); + + temp48alen = fast_expansion_sum_zeroelim(abclen, abc, eaclen, eac, temp48a); + temp48blen = fast_expansion_sum_zeroelim(bcelen, bce, eablen, eab, temp48b); + for (i = 0; i < temp48blen; i++) { + temp48b[i] = -temp48b[i]; + } + eabclen = fast_expansion_sum_zeroelim(temp48alen, temp48a, + temp48blen, temp48b, eabc); + dlen = scale_expansion_zeroelim(eabclen, eabc, dheight, ddet); + + temp48alen = fast_expansion_sum_zeroelim(bcdlen, bcd, abdlen, abd, temp48a); + temp48blen = fast_expansion_sum_zeroelim(cdalen, cda, abclen, abc, temp48b); + for (i = 0; i < temp48blen; i++) { + temp48b[i] = -temp48b[i]; + } + abcdlen = fast_expansion_sum_zeroelim(temp48alen, temp48a, + temp48blen, temp48b, abcd); + elen = scale_expansion_zeroelim(abcdlen, abcd, eheight, edet); + + ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); + cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet); + cdelen = fast_expansion_sum_zeroelim(cdlen, cddet, elen, edet, cdedet); + deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdelen, cdedet, deter); + + return deter[deterlen - 1]; +} + +REAL orient4dadapt(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL* pe, + REAL aheight, REAL bheight, REAL cheight, REAL dheight, + REAL eheight, REAL permanent) +{ + INEXACT REAL aex, bex, cex, dex, aey, bey, cey, dey, aez, bez, cez, dez; + INEXACT REAL aeheight, beheight, ceheight, deheight; + REAL det, errbound; + + INEXACT REAL aexbey1, bexaey1, bexcey1, cexbey1; + INEXACT REAL cexdey1, dexcey1, dexaey1, aexdey1; + INEXACT REAL aexcey1, cexaey1, bexdey1, dexbey1; + REAL aexbey0, bexaey0, bexcey0, cexbey0; + REAL cexdey0, dexcey0, dexaey0, aexdey0; + REAL aexcey0, cexaey0, bexdey0, dexbey0; + REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4]; + INEXACT REAL ab3, bc3, cd3, da3, ac3, bd3; + REAL abeps, bceps, cdeps, daeps, aceps, bdeps; + REAL temp8a[8], temp8b[8], temp8c[8], temp16[16], temp24[24]; + int temp8alen, temp8blen, temp8clen, temp16len, temp24len; + REAL adet[48], bdet[48], cdet[48], ddet[48]; + int alen, blen, clen, dlen; + REAL abdet[96], cddet[96]; + int ablen, cdlen; + REAL fin1[192]; + int finlength; + + REAL aextail, bextail, cextail, dextail; + REAL aeytail, beytail, ceytail, deytail; + REAL aeztail, beztail, ceztail, deztail; + REAL aeheighttail, beheighttail, ceheighttail, deheighttail; + + INEXACT REAL bvirt; + REAL avirt, bround, around; + INEXACT REAL c; + INEXACT REAL abig; + REAL ahi, alo, bhi, blo; + REAL err1, err2, err3; + INEXACT REAL _i, _j; + REAL _0; + + aex = (REAL) (pa[0] - pe[0]); + bex = (REAL) (pb[0] - pe[0]); + cex = (REAL) (pc[0] - pe[0]); + dex = (REAL) (pd[0] - pe[0]); + aey = (REAL) (pa[1] - pe[1]); + bey = (REAL) (pb[1] - pe[1]); + cey = (REAL) (pc[1] - pe[1]); + dey = (REAL) (pd[1] - pe[1]); + aez = (REAL) (pa[2] - pe[2]); + bez = (REAL) (pb[2] - pe[2]); + cez = (REAL) (pc[2] - pe[2]); + dez = (REAL) (pd[2] - pe[2]); + aeheight = (REAL) (aheight - eheight); + beheight = (REAL) (bheight - eheight); + ceheight = (REAL) (cheight - eheight); + deheight = (REAL) (dheight - eheight); + + Two_Product(aex, bey, aexbey1, aexbey0); + Two_Product(bex, aey, bexaey1, bexaey0); + Two_Two_Diff(aexbey1, aexbey0, bexaey1, bexaey0, ab3, ab[2], ab[1], ab[0]); + ab[3] = ab3; + + Two_Product(bex, cey, bexcey1, bexcey0); + Two_Product(cex, bey, cexbey1, cexbey0); + Two_Two_Diff(bexcey1, bexcey0, cexbey1, cexbey0, bc3, bc[2], bc[1], bc[0]); + bc[3] = bc3; + + Two_Product(cex, dey, cexdey1, cexdey0); + Two_Product(dex, cey, dexcey1, dexcey0); + Two_Two_Diff(cexdey1, cexdey0, dexcey1, dexcey0, cd3, cd[2], cd[1], cd[0]); + cd[3] = cd3; + + Two_Product(dex, aey, dexaey1, dexaey0); + Two_Product(aex, dey, aexdey1, aexdey0); + Two_Two_Diff(dexaey1, dexaey0, aexdey1, aexdey0, da3, da[2], da[1], da[0]); + da[3] = da3; + + Two_Product(aex, cey, aexcey1, aexcey0); + Two_Product(cex, aey, cexaey1, cexaey0); + Two_Two_Diff(aexcey1, aexcey0, cexaey1, cexaey0, ac3, ac[2], ac[1], ac[0]); + ac[3] = ac3; + + Two_Product(bex, dey, bexdey1, bexdey0); + Two_Product(dex, bey, dexbey1, dexbey0); + Two_Two_Diff(bexdey1, bexdey0, dexbey1, dexbey0, bd3, bd[2], bd[1], bd[0]); + bd[3] = bd3; + + temp8alen = scale_expansion_zeroelim(4, cd, bez, temp8a); + temp8blen = scale_expansion_zeroelim(4, bd, -cez, temp8b); + temp8clen = scale_expansion_zeroelim(4, bc, dez, temp8c); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, + temp8blen, temp8b, temp16); + temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, + temp16len, temp16, temp24); + alen = scale_expansion_zeroelim(temp24len, temp24, -aeheight, adet); + + temp8alen = scale_expansion_zeroelim(4, da, cez, temp8a); + temp8blen = scale_expansion_zeroelim(4, ac, dez, temp8b); + temp8clen = scale_expansion_zeroelim(4, cd, aez, temp8c); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, + temp8blen, temp8b, temp16); + temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, + temp16len, temp16, temp24); + blen = scale_expansion_zeroelim(temp24len, temp24, beheight, bdet); + + temp8alen = scale_expansion_zeroelim(4, ab, dez, temp8a); + temp8blen = scale_expansion_zeroelim(4, bd, aez, temp8b); + temp8clen = scale_expansion_zeroelim(4, da, bez, temp8c); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, + temp8blen, temp8b, temp16); + temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, + temp16len, temp16, temp24); + clen = scale_expansion_zeroelim(temp24len, temp24, -ceheight, cdet); + + temp8alen = scale_expansion_zeroelim(4, bc, aez, temp8a); + temp8blen = scale_expansion_zeroelim(4, ac, -bez, temp8b); + temp8clen = scale_expansion_zeroelim(4, ab, cez, temp8c); + temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, + temp8blen, temp8b, temp16); + temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, + temp16len, temp16, temp24); + dlen = scale_expansion_zeroelim(temp24len, temp24, deheight, ddet); + + ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet); + cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet); + finlength = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, fin1); + + det = estimate(finlength, fin1); + errbound = isperrboundB * permanent; + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + Two_Diff_Tail(pa[0], pe[0], aex, aextail); + Two_Diff_Tail(pa[1], pe[1], aey, aeytail); + Two_Diff_Tail(pa[2], pe[2], aez, aeztail); + Two_Diff_Tail(aheight, eheight, aeheight, aeheighttail); + Two_Diff_Tail(pb[0], pe[0], bex, bextail); + Two_Diff_Tail(pb[1], pe[1], bey, beytail); + Two_Diff_Tail(pb[2], pe[2], bez, beztail); + Two_Diff_Tail(bheight, eheight, beheight, beheighttail); + Two_Diff_Tail(pc[0], pe[0], cex, cextail); + Two_Diff_Tail(pc[1], pe[1], cey, ceytail); + Two_Diff_Tail(pc[2], pe[2], cez, ceztail); + Two_Diff_Tail(cheight, eheight, ceheight, ceheighttail); + Two_Diff_Tail(pd[0], pe[0], dex, dextail); + Two_Diff_Tail(pd[1], pe[1], dey, deytail); + Two_Diff_Tail(pd[2], pe[2], dez, deztail); + Two_Diff_Tail(dheight, eheight, deheight, deheighttail); + if ((aextail == 0.0) && (aeytail == 0.0) && (aeztail == 0.0) + && (bextail == 0.0) && (beytail == 0.0) && (beztail == 0.0) + && (cextail == 0.0) && (ceytail == 0.0) && (ceztail == 0.0) + && (dextail == 0.0) && (deytail == 0.0) && (deztail == 0.0) + && (aeheighttail == 0.0) && (beheighttail == 0.0) + && (ceheighttail == 0.0) && (deheighttail == 0.0)) { + return det; + } + + errbound = isperrboundC * permanent + resulterrbound * Absolute(det); + abeps = (aex * beytail + bey * aextail) + - (aey * bextail + bex * aeytail); + bceps = (bex * ceytail + cey * bextail) + - (bey * cextail + cex * beytail); + cdeps = (cex * deytail + dey * cextail) + - (cey * dextail + dex * ceytail); + daeps = (dex * aeytail + aey * dextail) + - (dey * aextail + aex * deytail); + aceps = (aex * ceytail + cey * aextail) + - (aey * cextail + cex * aeytail); + bdeps = (bex * deytail + dey * bextail) + - (bey * dextail + dex * beytail); + det += ((beheight + * ((cez * daeps + dez * aceps + aez * cdeps) + + (ceztail * da3 + deztail * ac3 + aeztail * cd3)) + + deheight + * ((aez * bceps - bez * aceps + cez * abeps) + + (aeztail * bc3 - beztail * ac3 + ceztail * ab3))) + - (aeheight + * ((bez * cdeps - cez * bdeps + dez * bceps) + + (beztail * cd3 - ceztail * bd3 + deztail * bc3)) + + ceheight + * ((dez * abeps + aez * bdeps + bez * daeps) + + (deztail * ab3 + aeztail * bd3 + beztail * da3)))) + + ((beheighttail * (cez * da3 + dez * ac3 + aez * cd3) + + deheighttail * (aez * bc3 - bez * ac3 + cez * ab3)) + - (aeheighttail * (bez * cd3 - cez * bd3 + dez * bc3) + + ceheighttail * (dez * ab3 + aez * bd3 + bez * da3))); + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + return orient4dexact(pa, pb, pc, pd, pe, + aheight, bheight, cheight, dheight, eheight); +} + +REAL orient4d(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL* pe, + REAL aheight, REAL bheight, REAL cheight, REAL dheight, + REAL eheight) +{ + REAL aex, bex, cex, dex; + REAL aey, bey, cey, dey; + REAL aez, bez, cez, dez; + REAL aexbey, bexaey, bexcey, cexbey, cexdey, dexcey, dexaey, aexdey; + REAL aexcey, cexaey, bexdey, dexbey; + REAL aeheight, beheight, ceheight, deheight; + REAL ab, bc, cd, da, ac, bd; + REAL abc, bcd, cda, dab; + REAL aezplus, bezplus, cezplus, dezplus; + REAL aexbeyplus, bexaeyplus, bexceyplus, cexbeyplus; + REAL cexdeyplus, dexceyplus, dexaeyplus, aexdeyplus; + REAL aexceyplus, cexaeyplus, bexdeyplus, dexbeyplus; + REAL det; + REAL permanent, errbound; + + //orient4dcount++; + + aex = pa[0] - pe[0]; + bex = pb[0] - pe[0]; + cex = pc[0] - pe[0]; + dex = pd[0] - pe[0]; + aey = pa[1] - pe[1]; + bey = pb[1] - pe[1]; + cey = pc[1] - pe[1]; + dey = pd[1] - pe[1]; + aez = pa[2] - pe[2]; + bez = pb[2] - pe[2]; + cez = pc[2] - pe[2]; + dez = pd[2] - pe[2]; + aeheight = aheight - eheight; + beheight = bheight - eheight; + ceheight = cheight - eheight; + deheight = dheight - eheight; + + aexbey = aex * bey; + bexaey = bex * aey; + ab = aexbey - bexaey; + bexcey = bex * cey; + cexbey = cex * bey; + bc = bexcey - cexbey; + cexdey = cex * dey; + dexcey = dex * cey; + cd = cexdey - dexcey; + dexaey = dex * aey; + aexdey = aex * dey; + da = dexaey - aexdey; + + aexcey = aex * cey; + cexaey = cex * aey; + ac = aexcey - cexaey; + bexdey = bex * dey; + dexbey = dex * bey; + bd = bexdey - dexbey; + + abc = aez * bc - bez * ac + cez * ab; + bcd = bez * cd - cez * bd + dez * bc; + cda = cez * da + dez * ac + aez * cd; + dab = dez * ab + aez * bd + bez * da; + + det = (deheight * abc - ceheight * dab) + (beheight * cda - aeheight * bcd); + + if (0) { //if (noexact) { + return det; + } + + aezplus = Absolute(aez); + bezplus = Absolute(bez); + cezplus = Absolute(cez); + dezplus = Absolute(dez); + aexbeyplus = Absolute(aexbey); + bexaeyplus = Absolute(bexaey); + bexceyplus = Absolute(bexcey); + cexbeyplus = Absolute(cexbey); + cexdeyplus = Absolute(cexdey); + dexceyplus = Absolute(dexcey); + dexaeyplus = Absolute(dexaey); + aexdeyplus = Absolute(aexdey); + aexceyplus = Absolute(aexcey); + cexaeyplus = Absolute(cexaey); + bexdeyplus = Absolute(bexdey); + dexbeyplus = Absolute(dexbey); + permanent = ((cexdeyplus + dexceyplus) * bezplus + + (dexbeyplus + bexdeyplus) * cezplus + + (bexceyplus + cexbeyplus) * dezplus) + * aeheight + + ((dexaeyplus + aexdeyplus) * cezplus + + (aexceyplus + cexaeyplus) * dezplus + + (cexdeyplus + dexceyplus) * aezplus) + * beheight + + ((aexbeyplus + bexaeyplus) * dezplus + + (bexdeyplus + dexbeyplus) * aezplus + + (dexaeyplus + aexdeyplus) * bezplus) + * ceheight + + ((bexceyplus + cexbeyplus) * aezplus + + (cexaeyplus + aexceyplus) * bezplus + + (aexbeyplus + bexaeyplus) * cezplus) + * deheight; + errbound = isperrboundA * permanent; + if ((det > errbound) || (-det > errbound)) { + return det; + } + + return orient4dadapt(pa, pb, pc, pd, pe, + aheight, bheight, cheight, dheight, eheight, permanent); +} + +} // end namespace diff --git a/contrib/Tetgen/CMakeLists.txt b/contrib/Tetgen1.4/CMakeLists.txt similarity index 84% rename from contrib/Tetgen/CMakeLists.txt rename to contrib/Tetgen1.4/CMakeLists.txt index e56b3669232e262f0c77366a24affd308306e15b..a8435d0935fc0ecf8802549de397dfad1fa19559 100644 --- a/contrib/Tetgen/CMakeLists.txt +++ b/contrib/Tetgen1.4/CMakeLists.txt @@ -9,4 +9,4 @@ set(SRC ) file(GLOB_RECURSE HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h) -append_gmsh_src(contrib/Tetgen "${SRC};${HDR}") +append_gmsh_src(contrib/Tetgen1.4 "${SRC};${HDR}") diff --git a/contrib/Tetgen/LICENSE b/contrib/Tetgen1.4/LICENSE similarity index 100% rename from contrib/Tetgen/LICENSE rename to contrib/Tetgen1.4/LICENSE diff --git a/contrib/Tetgen/README.txt b/contrib/Tetgen1.4/README.txt similarity index 100% rename from contrib/Tetgen/README.txt rename to contrib/Tetgen1.4/README.txt diff --git a/contrib/Tetgen/predicates.cxx b/contrib/Tetgen1.4/predicates.cxx similarity index 100% rename from contrib/Tetgen/predicates.cxx rename to contrib/Tetgen1.4/predicates.cxx diff --git a/contrib/Tetgen/tetgen.cxx b/contrib/Tetgen1.4/tetgen.cxx similarity index 100% rename from contrib/Tetgen/tetgen.cxx rename to contrib/Tetgen1.4/tetgen.cxx diff --git a/contrib/Tetgen/tetgen.h b/contrib/Tetgen1.4/tetgen.h similarity index 100% rename from contrib/Tetgen/tetgen.h rename to contrib/Tetgen1.4/tetgen.h diff --git a/contrib/TetgenNew/CMakeLists.txt b/contrib/Tetgen1.5/CMakeLists.txt similarity index 84% rename from contrib/TetgenNew/CMakeLists.txt rename to contrib/Tetgen1.5/CMakeLists.txt index f4a0965a4c93631cf27eefc5c9df8402d1156a8c..c228a2f0d95cdb8c512832d609b92bae0e508af7 100644 --- a/contrib/TetgenNew/CMakeLists.txt +++ b/contrib/Tetgen1.5/CMakeLists.txt @@ -9,4 +9,4 @@ set(SRC ) file(GLOB_RECURSE HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h) -append_gmsh_src(contrib/TetgenNew "${SRC};${HDR}") +append_gmsh_src(contrib/Tetgen1.5 "${SRC};${HDR}") diff --git a/contrib/TetgenNew/LICENSE b/contrib/Tetgen1.5/LICENSE similarity index 100% rename from contrib/TetgenNew/LICENSE rename to contrib/Tetgen1.5/LICENSE diff --git a/contrib/TetgenNew/README b/contrib/Tetgen1.5/README similarity index 100% rename from contrib/TetgenNew/README rename to contrib/Tetgen1.5/README diff --git a/contrib/TetgenNew/example.poly b/contrib/Tetgen1.5/example.poly similarity index 100% rename from contrib/TetgenNew/example.poly rename to contrib/Tetgen1.5/example.poly diff --git a/contrib/TetgenNew/makefile b/contrib/Tetgen1.5/makefile similarity index 100% rename from contrib/TetgenNew/makefile rename to contrib/Tetgen1.5/makefile diff --git a/contrib/TetgenNew/predicates.cxx b/contrib/Tetgen1.5/predicates.cxx similarity index 100% rename from contrib/TetgenNew/predicates.cxx rename to contrib/Tetgen1.5/predicates.cxx diff --git a/contrib/TetgenNew/tetgen.cxx b/contrib/Tetgen1.5/tetgen.cxx similarity index 100% rename from contrib/TetgenNew/tetgen.cxx rename to contrib/Tetgen1.5/tetgen.cxx diff --git a/contrib/TetgenNew/tetgen.h b/contrib/Tetgen1.5/tetgen.h similarity index 100% rename from contrib/TetgenNew/tetgen.h rename to contrib/Tetgen1.5/tetgen.h