Skip to content
Snippets Groups Projects
Commit 11e8ab0b authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

#ifdef levelset calls

parent 9633baa6
No related branches found
No related tags found
No related merge requests found
......@@ -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]);
}
......
......@@ -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
}
......@@ -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
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment