From 11e8ab0bb06721d173679cc2fe04ce1a8bd15223 Mon Sep 17 00:00:00 2001 From: Christophe Geuzaine <cgeuzaine@ulg.ac.be> Date: Thu, 25 Jun 2009 05:40:02 +0000 Subject: [PATCH] #ifdef levelset calls --- Common/OpenFile.cpp | 4 +-- Geo/MElementCut.cpp | 14 +++++++++++ contrib/Tetgen/flip.cxx | 54 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 70 insertions(+), 2 deletions(-) diff --git a/Common/OpenFile.cpp b/Common/OpenFile.cpp index 059fd69b05..8f44c1df37 100644 --- a/Common/OpenFile.cpp +++ b/Common/OpenFile.cpp @@ -80,8 +80,8 @@ void SetBoundingBox(double xmin, double xmax, CTX::instance()->min[2] = zmin; CTX::instance()->max[2] = zmax; FinishUpBoundingBox(); CTX::instance()->lc = sqrt(SQU(CTX::instance()->max[0] - CTX::instance()->min[0]) + - SQU(CTX::instance()->max[1] - CTX::instance()->min[1]) + - SQU(CTX::instance()->max[2] - CTX::instance()->min[2])); + SQU(CTX::instance()->max[1] - CTX::instance()->min[1]) + + SQU(CTX::instance()->max[2] - CTX::instance()->min[2])); for(int i = 0; i < 3; i++) CTX::instance()->cg[i] = 0.5 * (CTX::instance()->min[i] + CTX::instance()->max[i]); } diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp index 5eb387210d..bdf8e614ed 100644 --- a/Geo/MElementCut.cpp +++ b/Geo/MElementCut.cpp @@ -6,8 +6,13 @@ #include "GModel.h" #include "MElement.h" #include "MElementCut.h" + +//#define HAVE_LEVELSET + +#if defined(HAVE_LEVELSET) #include "DILevelset.h" #include "Integration3D.h" +#endif //---------------------------------------- MPolyhedron ---------------------------- @@ -175,6 +180,8 @@ void MPolygon::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const //---------------------------------------- CutMesh ---------------------------- +#if defined(HAVE_LEVELSET) + int getElementVertexNum (DI_Point p, MElement *e) { for(int i = 0; i < e->getNumVertices(); i++) @@ -573,11 +580,14 @@ static void elementCutMesh (MElement *e, gLevelset *ls, GEntity *ge, GModel *GM, } } +#endif + GModel *buildCutMesh(GModel *gm, gLevelset *ls, std::map<int, std::vector<MElement*> > elements[10], std::map<int, MVertex*> &vertexMap, std::map<int, std::map<int, std::string> > physicals[4]) { +#if defined(HAVE_LEVELSET) GModel *cutGM = new GModel; std::map<int, std::vector<MElement*> > border[2]; @@ -635,5 +645,9 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls, }printf("numbering borders finished \n"); return cutGM; +#else + Msg::Error("Gmsh need to be compiled with levelset support to cut mesh"); + return 0; +#endif } diff --git a/contrib/Tetgen/flip.cxx b/contrib/Tetgen/flip.cxx index d32dd83fc9..fbaf010acf 100644 --- a/contrib/Tetgen/flip.cxx +++ b/contrib/Tetgen/flip.cxx @@ -1899,11 +1899,17 @@ void tetgenmesh::lawsonflip3d(int flipflag) int *iptr; + // for flipflag = 2. + point pa, pb, ppa, ppb; + int ecount; + if (b->verbose > 1) { printf(" Lawson flip %ld faces.\n", flippool->items); flipcount = flip23count + flip32count; } + ecount = 0; // Initialize the counter. + while (futureflip != (badface *) NULL) { // Pop a face from the stack. @@ -2018,6 +2024,9 @@ void tetgenmesh::lawsonflip3d(int flipflag) // Found a 3-to-2 flip. flip32(fliptets, 0, flipflag); recenttet = fliptets[0]; // for point location. + if (flipflag == 2) { + ecount = 0; // Clear the counter. + } } else if ((n == 4) && (ori == 0)) { // Find a 4-to-4 flip. flipnmcount++; @@ -2034,6 +2043,9 @@ void tetgenmesh::lawsonflip3d(int flipflag) fliptets[2] = baktets[1]; flip32(fliptets, 1, flipflag); // hull tet may involve. recenttet = fliptets[0]; // for point location. + if (flipflag == 2) { + ecount = 0; // Clear the counter. + } } else { // An unflipable face. Will be flipped later. if (flipflag == 2) { // if (flipflag > 1) { @@ -2056,6 +2068,48 @@ void tetgenmesh::lawsonflip3d(int flipflag) if (b->verbose > 1) { printf("\n"); } + ecount++; // Increase the counter. + if (ecount >= 1000) { + // assert(0); // A flip deadlock. + // Dump the dead lock case. + pa = org(fliptets[0]); + pb = dest(fliptets[0]); + if (ecount == 1000) { + printf("-- Flip failed after inserting point %ld.\n", + pointpool->items); + ppa = pa; // Bakup the startiing loop edge. + ppb = pb; + } else { + if (((ppa == pa) && (ppb == pb)) || + ((ppa == pb) && (ppb == pa))) { + // Dump the current mesh and exit. + outnodes(0); + outelements(0); + exit(1); + } + } + printf("p:draw_subseg(%d, %d)\n",pointmark(pa),pointmark(pb)); + printf("p:draw_subface(%d, %d, %d)\n", pointmark(org(fliptet)), + pointmark(dest(fliptet)), pointmark(apex(fliptet))); + printf("p:draw_tet(%d, %d, %d, %d)\n", pointmark(org(fliptet)), + pointmark(dest(fliptet)), pointmark(apex(fliptet)), + pointmark(oppo(fliptet))); + symedge(fliptet, fliptets[1]); + printf("p:draw_tet(%d, %d, %d, %d)\n", + pointmark(org(fliptets[1])), + pointmark(dest(fliptets[1])), + pointmark(apex(fliptets[1])), + pointmark(oppo(fliptets[1]))); + pe = apex(fliptets[0]); + fliptets[1] = fliptets[0]; + while (1) { + pd = oppo(fliptets[1]); + printf(" -- p:draw_tet(%d, %d, %d, %d)\n", pointmark(pa), + pointmark(pb),pointmark(apex(fliptets[1])),pointmark(pd)); + fnextself(fliptets[1]); + if (apex(fliptets[1]) == pe) break; + } + } // if (ecount >= 1000) } } } // bflag -- GitLab