diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 0ec42e3d6a9d37b40150d19c2797b95a9e7dccc0..558cc50cd5f41253c8b86c5d18bc6ee27c5ee136 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -116,8 +116,6 @@ class ChainComplex{
virtual void matrixTest();
};
-#endif
-
// A class representing a chain.
// Used to visualize generators of the homology groups.
class Chain{
@@ -165,3 +163,5 @@ class Chain{
};
#endif
+
+#endif
diff --git a/Geo/Makefile b/Geo/Makefile
index dbb7036a8223db67b97ca26761effcb3805fed5d..b8c6c75bca84136ad0407212be324ceafb0c31d0 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -40,9 +40,7 @@ SRC = GEntity.cpp STensor3.cpp\
MLine.cpp MTriangle.cpp MQuadrangle.cpp MTetrahedron.cpp\
MHexahedron.cpp MPrism.cpp MPyramid.cpp\
MZone.cpp MZoneBoundary.cpp\
- CellComplex.cpp ChainComplex.cpp Homology.cpp\
- STensor3.cpp
-
+ CellComplex.cpp ChainComplex.cpp Homology.cpp
OBJ = ${SRC:.cpp=${OBJEXT}}
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 627a8776c5790d6272067c3c92ca607af858adb2..5f7557cb8d726b61b3e2da63c4f9cc7209f00ca6 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -158,7 +158,10 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
double guess[2] = {0,0};
int Count = 0;
for(int j = 0; j < 3; j++){
- if(v[j]->onWhat()->dim() == 2){
+ if(!v[j]->onWhat()){
+ Msg::Error("Uncategorized vertex %d", v[j]->getNum());
+ }
+ else if(v[j]->onWhat()->dim() == 2){
double uu,vv;
v[j]->getParameter(0,uu);
v[j]->getParameter(1,vv);
diff --git a/contrib/Tetgen/Makefile b/contrib/Tetgen/Makefile
index 525e7e56bee6fb55606c09325815d1291517a779..e2296d4dfa15f7073e43534eeab75ed58c7ae377 100644
--- a/contrib/Tetgen/Makefile
+++ b/contrib/Tetgen/Makefile
@@ -21,6 +21,7 @@ SRC = behavior.cxx\
meshio.cxx\
meshstat.cxx\
predicates.cxx\
+ refine.cxx\
surface.cxx
OBJ = ${SRC:.cxx=${OBJEXT}}
diff --git a/contrib/Tetgen/behavior.cxx b/contrib/Tetgen/behavior.cxx
index db89a9a53edc7e2e7632a967c282b70c4488b2d1..e25f7e772cc7286f3daf0d8dd1cfeb4449e96ae7 100644
--- a/contrib/Tetgen/behavior.cxx
+++ b/contrib/Tetgen/behavior.cxx
@@ -50,7 +50,7 @@ tetgenbehavior::tetgenbehavior()
nofacewritten = 0;
noiterationnum = 0;
nobisect = 0;
- steiner = -1;
+ steinerleft = -1l;
nomerge = 0;
docheck = 0;
quiet = 0;
@@ -360,7 +360,7 @@ bool tetgenbehavior::parse_commandline(int argc, char **argv)
k++;
}
workstring[k] = '\0';
- steiner = (int) strtol(workstring, (char **) NULL, 0);
+ steinerleft = (int) strtol(workstring, (char **) NULL, 0);
}
} else if (argv[i][j] == 'D') {
conformdel++;
diff --git a/contrib/Tetgen/constrain.cxx b/contrib/Tetgen/constrain.cxx
index 377ded1e02394c3173a112c45816ce797db28f00..b2e175540f0d656e3cd9f10ff8a88fdeec2a0e69 100644
--- a/contrib/Tetgen/constrain.cxx
+++ b/contrib/Tetgen/constrain.cxx
@@ -357,16 +357,14 @@ enum tetgenmesh::intersection tetgenmesh::finddirection(triface* searchtet,
// scoutsegment() Look for a given segment in the tetrahedralization T. //
// //
// Search an edge in the tetrahedralization that matches the given segmment. //
-// If such an edge is found, the segment is 'locked' at the edge. //
-// //
-// If 'searchtet' != NULL, it's origin must be the origin of 'sseg'. It is //
-// used as the starting tet for searching the edge. //
+// If such an edge exists, the segment is 'locked' at the edge. 'searchtet' //
+// returns this (constrained) edge. Otherwise, the segment is missing. //
// //
// The returned value indicates one of the following cases: //
-// - SHAREVERT, the segment exists and is inserted in T; //
-// - ACROSSVERT, a vertex ('refpt') lies on the segment; //
-// - ACROSSEDGE, the segment is missing; //
-// - ACROSSFACE, the segment is missing; //
+// - SHAREEDGE, the segment exists and is inserted in T; //
+// - ACROSSVERT, the segment intersects a vertex ('refpt'). //
+// - ACROSSEDGE, the segment intersects an edge (in 'searchtet'). //
+// - ACROSSFACE, the segment crosses a face (in 'searchtet'). //
// //
// If the returned value is ACROSSEDGE or ACROSSFACE, i.e., the segment is //
// missing, 'refpt' returns the reference point for splitting thus segment, //
@@ -385,44 +383,22 @@ enum tetgenmesh::intersection tetgenmesh::scoutsegment(face* sseg,
enum intersection dir;
REAL angmax, ang;
long facecount;
- bool orgflag;
int types[2], poss[4];
- int shver, pos, i;
+ int pos, i;
tetrahedron ptr;
shellface sptr;
int *iptr, tver;
+ startpt = sorg(*sseg);
+ endpt = sdest(*sseg);
+
// Is 'searchtet' a valid handle?
if (searchtet->tet == NULL) {
- orgflag = false;
- // Search a tet whose origin is one of the endpoints of 'sseg'.
- for (shver = 0; shver < 2 && !orgflag; shver++) {
- startpt = (point) sseg->sh[shver + 3];
- decode(point2tet(startpt), *searchtet);
- if ((searchtet->tet != NULL) && (searchtet->tet[4] != NULL)) {
- // Check if this tet contains pa.
- for (i = 4; i < 8 && !orgflag; i++) {
- if ((point) searchtet->tet[i] == startpt) {
- // Found. Set pa as its origin.
- switch (i) {
- case 4: searchtet->loc = 0; searchtet->ver = 0; break;
- case 5: searchtet->loc = 0; searchtet->ver = 2; break;
- case 6: searchtet->loc = 0; searchtet->ver = 4; break;
- case 7: searchtet->loc = 1; searchtet->ver = 2; break;
- }
- sseg->shver = shver;
- orgflag = true;
- }
- }
- }
- }
- assert(orgflag); // SELF_CHECK
+ point2tetorg(startpt, *searchtet);
} else {
- startpt = sorg(*sseg);
assert(org(*searchtet) == startpt); // SELF_CHECK
}
- endpt = sdest(*sseg);
if (b->verbose > 1) {
printf(" Scout seg (%d, %d).\n", pointmark(startpt), pointmark(endpt));
@@ -436,6 +412,8 @@ enum tetgenmesh::intersection tetgenmesh::scoutsegment(face* sseg,
// Found! Insert the segment.
tsspivot(*searchtet, checkseg); // SELF_CHECK
if (checkseg.sh == NULL) {
+ // Let the segment remember an adjacent tet.
+ sstbond(*sseg, *searchtet);
neightet = *searchtet;
do {
tssbond1(neightet, *sseg);
@@ -447,7 +425,7 @@ enum tetgenmesh::intersection tetgenmesh::scoutsegment(face* sseg,
assert(checkseg.sh == sseg->sh); // SELF_CHECK
}
// The job is done.
- return SHAREVERT;
+ return SHAREEDGE;
} else {
// A point is on the path.
*refpt = pd;
@@ -773,6 +751,216 @@ void tetgenmesh::getsegmentsplitpoint(face* sseg, point refpt, REAL* vt)
}
}
+// This function only use Rule-1 to split the segment. refpt may be NULL.
+// Added 2009-05-20.
+
+void tetgenmesh::getsegmentsplitpoint2(face* sseg, point refpt, REAL* vt)
+{
+ point ei, ej;
+ REAL split, L, d, d1, d2;
+ int i;
+
+ ei = sorg(*sseg);
+ ej = sdest(*sseg);
+
+ if (b->verbose > 1) {
+ printf(" Split seg(%d, %d) ref(%d).\n", pointmark(ei), pointmark(ej),
+ refpt ? pointmark(refpt) : -1);
+ }
+
+ if (refpt != NULL) {
+ L = DIST(ei, ej);
+ d1 = DIST(ei, refpt);
+ d2 = DIST(ej, refpt);
+ if (d1 < d2) {
+ // Choose ei as center.
+ if (d1 < 0.5 * L) {
+ split = d1 / L;
+ // Adjust split if it is close to middle. (2009-02-01)
+ // if ((split > 0.4) || (split < 0.6)) split = 0.5;
+ } else {
+ split = 0.5;
+ }
+ for (i = 0; i < 3; i++) {
+ vt[i] = ei[i] + split * (ej[i] - ei[i]);
+ }
+ } else {
+ // Choose ej as center.
+ if (d2 < 0.5 * L) {
+ split = d2 / L;
+ // Adjust split if it is close to middle. (2009-02-01)
+ // if ((split > 0.4) || (split < 0.6)) split = 0.5;
+ } else {
+ split = 0.5;
+ }
+ for (i = 0; i < 3; i++) {
+ vt[i] = ej[i] + split * (ei[i] - ej[i]);
+ }
+ }
+ } else {
+ split = 0.5;
+ for (i = 0; i < 3; i++) {
+ vt[i] = ei[i] + split * (ej[i] - ei[i]);
+ }
+ }
+ r1count++;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// getsegmentsplitpoint3() Calculate a split point in the given segment.
+// This routine does not check far origin and far dest.
+
+void tetgenmesh::getsegmentsplitpoint3(face* sseg, point refpt, REAL* vt)
+{
+ point ei, ej, ek;
+ REAL split, L, d, d1, d2, d3;
+ int stype, sign;
+ int i;
+
+ // Decide the type of this segment.
+ sign = 1;
+ ei = sorg(*sseg);
+ ej = sdest(*sseg);
+ ek = NULL;
+
+ // ei can be an ACUTEVERTEX, or a RIDGEVERTEX, or a STEINERVERTEX.
+ if (getpointtype(ei) == ACUTEVERTEX) {
+ if (getpointtype(ej) == ACUTEVERTEX) {
+ // ej is an ACUTEVERTEX.
+ stype = 0;
+ } else {
+ // ej is either a RIDGEVERTEX or a STEINERVERTEX.
+ stype = 1;
+ ek = ei;
+ }
+ } else {
+ if (getpointtype(ei) == RIDGEVERTEX) {
+ if (getpointtype(ej) == ACUTEVERTEX) {
+ // ej is an ACUTEVERTEX.
+ stype = 1; sign = -1;
+ ek = ej;
+ } else {
+ if (getpointtype(ej) == RIDGEVERTEX) {
+ // ej is a RIDGEVERTEX.
+ stype = 0;
+ } else {
+ // ej is a STEINERVETEX.
+ stype = 0;
+ /* ek = farsdest(*sseg);
+ if (getpointtype(ek) == ACUTEVERTEX) {
+ stype = 1; sign = -1;
+ } else {
+ stype = 0;
+ }*/
+ }
+ }
+ } else {
+ // ei is a STEINERVERTEX.
+ if (getpointtype(ej) == ACUTEVERTEX) {
+ stype = 1; sign = -1;
+ ek = ej;
+ } else {
+ // ej is either a RIDGEVERTEX or STEINERVERTEX.
+ stype = 0;
+ /*ek = farsorg(*sseg);
+ if (getpointtype(ej) == RIDGEVERTEX) {
+ if (getpointtype(ek) == ACUTEVERTEX) {
+ stype = 1;
+ } else {
+ stype = 0;
+ }
+ } else {
+ // Both ei and ej are STEINERVETEXs. ei has priority.
+ if (getpointtype(ek) == ACUTEVERTEX) {
+ stype = 1;
+ } else {
+ ek = farsdest(*sseg);
+ if (getpointtype(ek) == ACUTEVERTEX) {
+ stype = 1; sign = -1;
+ } else {
+ stype = 0;
+ }
+ }
+ }*/
+ }
+ }
+ }
+
+ // Adjust the endpoints: ei, ej.
+ if (sign == -1) {
+ sesymself(*sseg);
+ ei = sorg(*sseg);
+ ej = sdest(*sseg);
+ }
+
+ if (b->verbose > 1) {
+ printf(" Split a type-%d seg(%d, %d) ref(%d).\n", stype,
+ pointmark(ei), pointmark(ej), pointmark(refpt));
+ }
+
+ // Calculate the split point.
+ if (stype == 0) {
+ // Use rule-1.
+ L = DIST(ei, ej);
+ d1 = DIST(ei, refpt);
+ d2 = DIST(ej, refpt);
+ if (d1 < d2) {
+ // Choose ei as center.
+ if (d1 < 0.5 * L) {
+ split = d1 / L;
+ // Adjust split if it is close to middle. (2009-02-01)
+ if ((split > 0.4) || (split < 0.6)) split = 0.5;
+ } else {
+ split = 0.5;
+ }
+ for (i = 0; i < 3; i++) {
+ vt[i] = ei[i] + split * (ej[i] - ei[i]);
+ }
+ } else {
+ // Choose ej as center.
+ if (d2 < 0.5 * L) {
+ split = d2 / L;
+ // Adjust split if it is close to middle. (2009-02-01)
+ if ((split > 0.4) || (split < 0.6)) split = 0.5;
+ } else {
+ split = 0.5;
+ }
+ for (i = 0; i < 3; i++) {
+ vt[i] = ej[i] + split * (ei[i] - ej[i]);
+ }
+ }
+ r1count++;
+ } else {
+ // Use rule-2.
+ // ek = farsorg(*sseg);
+ L = DIST(ek, ej);
+ d = DIST(ek, refpt);
+ split = d / L;
+ for (i = 0; i < 3; i++) {
+ vt[i] = ek[i] + split * (ej[i] - ek[i]);
+ }
+ d1 = DIST(vt, refpt);
+ d2 = DIST(vt, ej);
+ if (d1 > d2) {
+ // Use rule-3.
+ d3 = DIST(ei, refpt);
+ if (d1 < 0.5 * d3) {
+ split = (d - d1) / L;
+ } else {
+ split = (d - 0.5 * d3) / L;
+ }
+ for (i = 0; i < 3; i++) {
+ vt[i] = ek[i] + split * (ej[i] - ek[i]);
+ }
+ }
+ d1 > d2 ? r3count++ : r2count++;
+ }
+
+ if (b->verbose > 1) {
+ printf(" split (%g), vt (%g, %g, %g).\n", split, vt[0], vt[1], vt[2]);
+ }
+}
+
///////////////////////////////////////////////////////////////////////////////
// //
// delaunizesegments() Recover segments in a Delaunay tetrahedralization. //
@@ -788,10 +976,6 @@ void tetgenmesh::delaunizesegments()
enum intersection dir;
bool visflag;
- if (b->verbose) {
- printf(" Delaunizing segments.\n");
- }
-
// Loop until 'subsegstack' is empty.
while (subsegstack->objects > 0l) {
// seglist is used as a stack.
@@ -806,7 +990,7 @@ void tetgenmesh::delaunizesegments()
searchtet.tet = NULL;
dir = scoutsegment(&sseg, &searchtet, &refpt);
- if (dir != SHAREVERT) {
+ if (dir != SHAREEDGE) {
// The segment is missing, split it.
spivot(sseg, splitsh);
if (dir != ACROSSVERT) {
@@ -814,6 +998,10 @@ void tetgenmesh::delaunizesegments()
makepoint(&newpt);
getsegmentsplitpoint(&sseg, refpt, newpt);
setpointtype(newpt, STEINERVERTEX);
+ // DEBUG
+ if (pointmark(newpt) == 4157) {
+ printf("debug point\n");
+ }
// Split the segment by newpt.
sinsertvertex(newpt, &splitsh, &sseg, true, false);
// Insert newpt into the DT. If 'checksubfaces == 1' the current
@@ -836,10 +1024,6 @@ void tetgenmesh::delaunizesegments()
}
}
}
-
- if (b->verbose) {
- printf(" %d protecting points.\n", r1count + r2count + r3count);
- }
}
///////////////////////////////////////////////////////////////////////////////
@@ -1618,6 +1802,7 @@ bool tetgenmesh::delaunizecavity(arraypool *cavpoints, arraypool *cavfaces,
{
triface *parytet, searchtet, neightet, spintet, *parytet1;
face checksh, tmpsh, *parysh;
+ face checkseg;
point pa, pb, pc, pd, pt[3], *parypt;
// badface *newflipface;
enum intersection dir;
@@ -1673,18 +1858,18 @@ bool tetgenmesh::delaunizecavity(arraypool *cavpoints, arraypool *cavfaces,
pt[0] = * (point *) fastlookup(cavpoints, i);
assert(pt[0] != dummypoint); // SELF_CHECK
if (!pinfected(pt[0])) {
- pinfect(pt[0]); // Mark it as inserted.
+ // pinfect(pt[0]); // Mark it as inserted.
searchtet = recenttet;
insertvertex(pt[0], &searchtet, true, false, false, false);
} else {
- // puninfect(pt[0]); // It is already inserted.
+ puninfect(pt[0]); // It is already inserted.
}
}
- // Comment: All vertices of the cavity are marked.
+ // Comment: All vertices of the cavity are NOT marked.
while (1) {
- // Indentify boundary faces. Mark interior tets. Save missing faces.
+ // Identify boundary faces. Mark interior tets. Save missing faces.
for (i = 0; i < cavfaces->objects; i++) {
parytet = (triface *) fastlookup(cavfaces, i);
// Skip an interior face (due to the enlargement of the cavity).
@@ -1713,10 +1898,11 @@ bool tetgenmesh::delaunizecavity(arraypool *cavpoints, arraypool *cavfaces,
tsbond(neightet, tmpsh);
}
assert(dest(neightet) == pt[0]); // SELF_CHECK
- // Mark neightet as interior.
- if (!infected(neightet)) {
- infect(neightet);
- }
+ // The following step is now done in fillcavity(), 2009-04-24.
+ // // Mark neightet as interior.
+ // if (!infected(neightet)) {
+ // infect(neightet);
+ // }
} else if (dir == COLLISIONFACE) {
// A subface is already inserted (see fig/dum-cavity-case6).
assert(oppo(*parytet) == NULL); // It must be a faked tet.
@@ -1777,22 +1963,17 @@ bool tetgenmesh::delaunizecavity(arraypool *cavpoints, arraypool *cavfaces,
shellfacedealloc(subfacepool, parysh->sh);
}
cavshells->restart();
+
+ // Infect the points which are of the cavity.
+ for (i = 0; i < cavpoints->objects; i++) {
+ pt[0] = * (point *) fastlookup(cavpoints, i);
+ pinfect(pt[0]); // Mark it as inserted.
+ }
+
// Enlarge the cavity.
for (i = 0; i < misfaces->objects; i++) {
// Get a missing face.
parytet = (triface *) fastlookup(misfaces, i);
- // Check for a missing subface.
- tspivot(*parytet, checksh);
- if (checksh.sh != NULL) {
- if (b->verbose > 1) {
- printf(" Queue a subface x%lx (%d, %d, %d).\n",
- (unsigned long) checksh.sh, pointmark(sorg(checksh)),
- pointmark(sdest(checksh)), pointmark(sapex(checksh)));
- }
- stdissolve(checksh);
- subfacstack->newindex((void **) &parysh);
- *parysh = checksh;
- }
if (!infected(*parytet)) {
// Put it into crossing tet list.
infect(*parytet);
@@ -1810,6 +1991,18 @@ bool tetgenmesh::delaunizecavity(arraypool *cavpoints, arraypool *cavfaces,
searchtet = recenttet;
insertvertex(pd, &searchtet, true, false, false, false);
}
+ // Check for a missing subface.
+ tspivot(*parytet, checksh);
+ if (checksh.sh != NULL) {
+ if (b->verbose > 1) {
+ printf(" Queue a subface x%lx (%d, %d, %d).\n",
+ (unsigned long) checksh.sh, pointmark(sorg(checksh)),
+ pointmark(sdest(checksh)), pointmark(sapex(checksh)));
+ }
+ stdissolve(checksh);
+ subfacstack->newindex((void **) &parysh);
+ *parysh = checksh;
+ }
// Add three opposite faces into the boundary list.
for (j = 0; j < 3; j++) {
enext0fnext(*parytet, neightet);
@@ -1840,6 +2033,13 @@ bool tetgenmesh::delaunizecavity(arraypool *cavpoints, arraypool *cavfaces,
} // j
} // if (!infected(parytet))
}
+
+ // Uninfect the points which are of the cavity.
+ for (i = 0; i < cavpoints->objects; i++) {
+ pt[0] = * (point *) fastlookup(cavpoints, i);
+ puninfect(pt[0]);
+ }
+
misfaces->restart();
cavityexpcount++;
continue;
@@ -1865,11 +2065,6 @@ bool tetgenmesh::delaunizecavity(arraypool *cavpoints, arraypool *cavfaces,
}
}
- // Uninfect all points of the DT.
- for (i = 0; i < cavpoints->objects; i++) {
- parypt = (point *) fastlookup(cavpoints, i);
- puninfect(*parypt);
- }
cavpoints->restart();
// Comment: Now no vertex is marked.
cavfaces->restart();
@@ -1909,50 +2104,59 @@ bool tetgenmesh::fillcavity(arraypool* topshells, arraypool* botshells,
// Connect newtets to tets outside the cavity.
for (k = 0; k < 2; k++) {
cavshells = (k == 0 ? topshells : botshells);
- for (i = 0; i < cavshells->objects; i++) {
- // Get a temp subface.
- parysh = (face *) fastlookup(cavshells, i);
- // Get the boundary tet outsode the cavity.
- decode(parysh->sh[0], bdrytet);
- pa = org(bdrytet);
- pb = dest(bdrytet);
- pc = apex(bdrytet);
- // Get the adjacent new tet.
- stpivot(*parysh, neightet);
- assert(org(neightet) == pb); // SELF_CHECK
- assert(dest(neightet) == pa); // SELF_CHECK
- if (oppo(bdrytet) != NULL) {
- // Bond the two tets.
- bond(bdrytet, neightet); // Also cleared the pointer to tmpsh.
- }
- // Bond a subface (if it exists).
- tspivot(bdrytet, checksh);
- if (checksh.sh != NULL) {
- tsbond(neightet, checksh); // Also cleared the pointer to tmpsh.
- } else {
- tsdissolve(neightet); // No subface, clear the pointer to tmpsh.
- }
- // Update the point-to-tets map.
- point2tet(pa) = encode(neightet);
- point2tet(pb) = encode(neightet);
- point2tet(pc) = encode(neightet);
- // Delete the temp subface.
- // shellfacedealloc(subfacepool, parysh->sh);
- if (oppo(bdrytet) == NULL) {
- // Delete a faked tet.
- tetrahedrondealloc(bdrytet.tet);
+ if (cavshells != NULL) {
+ for (i = 0; i < cavshells->objects; i++) {
+ // Get a temp subface.
+ parysh = (face *) fastlookup(cavshells, i);
+ // Get the boundary tet outsode the cavity.
+ decode(parysh->sh[0], bdrytet);
+ pa = org(bdrytet);
+ pb = dest(bdrytet);
+ pc = apex(bdrytet);
+ // Get the adjacent new tet.
+ stpivot(*parysh, neightet);
+ assert(org(neightet) == pb); // SELF_CHECK
+ assert(dest(neightet) == pa); // SELF_CHECK
+ // Mark neightet as an interior tet of this cavity, 2009-04-24.
+ // Comment: We know neightet is an interior tet.
+ if (!infected(neightet)) {
+ infect(neightet);
+ }
+ if (oppo(bdrytet) != NULL) {
+ // Bond the two tets.
+ bond(bdrytet, neightet); // Also cleared the pointer to tmpsh.
+ }
+ // Bond a subface (if it exists).
+ tspivot(bdrytet, checksh);
+ if (checksh.sh != NULL) {
+ tsbond(neightet, checksh); // Also cleared the pointer to tmpsh.
+ } else {
+ tsdissolve(neightet); // No subface, clear the pointer to tmpsh.
+ }
+ // Update the point-to-tets map.
+ point2tet(pa) = encode(neightet);
+ point2tet(pb) = encode(neightet);
+ point2tet(pc) = encode(neightet);
+ // Delete the temp subface.
+ // shellfacedealloc(subfacepool, parysh->sh);
+ if (oppo(bdrytet) == NULL) {
+ // Delete a faked tet.
+ tetrahedrondealloc(bdrytet.tet);
+ }
}
- }
+ } // if (cavshells != NULL)
}
+ mflag = true; // Initialize it.
+
+ if (midfaces != NULL) {
+
// Mark all facet vertices for finding middle subfaces.
for (i = 0; i < facpoints->objects; i++) {
pf = * (point *) fastlookup(facpoints, i);
pinfect(pf);
}
- mflag = true; // Initialize it.
-
// The first pair of top and bottom tets share the same edge [a, b].
// toptet = * (triface *) fastlookup(topfaces, 0);
if (infected(firsttopface)) {
@@ -2086,81 +2290,90 @@ bool tetgenmesh::fillcavity(arraypool* topshells, arraypool* botshells,
} // j
} // i
+ } // if (midfaces != NULL)
if (mflag) {
- if (b->verbose > 1) {
- printf(" Found %ld middle subfaces.\n", midfaces->objects);
- }
- if (midfaces->objects > maxregionsize) {
- maxregionsize = midfaces->objects;
- }
- // Unmark middle faces.
- for (i = 0; i < midfaces->objects; i++) {
- // Get a matched middle face [a, b, c]
- midface = * (triface *) fastlookup(midfaces, i);
- assert(facemarked(midface)); // SELF_CHECK
- unmarkface(midface);
+ if (midfaces != NULL) {
+ if (b->verbose > 1) {
+ printf(" Found %ld middle subfaces.\n", midfaces->objects);
+ }
+ if (midfaces->objects > maxregionsize) {
+ maxregionsize = midfaces->objects;
+ }
+ // Unmark middle faces.
+ for (i = 0; i < midfaces->objects; i++) {
+ // Get a matched middle face [a, b, c]
+ midface = * (triface *) fastlookup(midfaces, i);
+ assert(facemarked(midface)); // SELF_CHECK
+ unmarkface(midface);
+ }
}
// Bond subsegments to new tets.
// Comment: *** The following code does redundant job. Should be
// re-placed in the future.
for (k = 0; k < 2; k++) {
cavshells = (k == 0 ? topshells : botshells);
- for (i = 0; i < cavshells->objects; i++) {
- parysh = (face *) fastlookup(cavshells, i);
- decode(parysh->sh[0], bdrytet);
- if (bdrytet.tet[4] != NULL) {
- // Not a faked tet. Bond a subsegment (if it exists).
- for (j = 0; j < 3; j++) {
- tsspivot(bdrytet, checkseg);
- if (checkseg.sh != NULL) {
- symedge(bdrytet, neightet);
- assert(marktested(neightet)); // SELF_CHECK
- while (1) {
- tssbond1(neightet, checkseg);
- fnextself(neightet);
- if (!marktested(neightet)) break;
+ if (cavshells != NULL) {
+ for (i = 0; i < cavshells->objects; i++) {
+ parysh = (face *) fastlookup(cavshells, i);
+ decode(parysh->sh[0], bdrytet);
+ if (bdrytet.tet[4] != NULL) {
+ // Not a faked tet. Bond a subsegment (if it exists).
+ for (j = 0; j < 3; j++) {
+ tsspivot(bdrytet, checkseg);
+ if (checkseg.sh != NULL) {
+ symedge(bdrytet, neightet);
+ assert(marktested(neightet)); // SELF_CHECK
+ // Let the segment remember an adjacent tet.
+ sstbond(checkseg, neightet);
+ while (1) {
+ tssbond1(neightet, checkseg);
+ fnextself(neightet);
+ if (!marktested(neightet)) break;
+ }
}
+ enextself(bdrytet);
}
- enextself(bdrytet);
- }
- } else {
- // A faked tet. There is an interior subface. Use it.
- // See fig/dump-cavity-case19.
- stpivot(*parysh, neightet);
- assert(marktested(neightet)); // SELF_CHECK
- tspivot(neightet, checksh);
- assert(checksh.sh != NULL); // SELF_CHECK
- assert(checksh.sh != parysh->sh); // // SELF_CHECK
- // Align them at the same directed edge.
- pa = org(neightet);
- pb = dest(neightet);
- for (j = 0; j < 3; j++) {
- if (sorg(checksh) == pa) break;
- senextself(checksh);
- }
- assert(j < 3); // SELF_CHECK
- if (sdest(checksh) != pb) {
- senext2self(checksh);
- sesymself(checksh);
- }
- assert(sdest(checksh) == pb); // SELF_CHECK
- // Bond a subsegment (if it exists).
- for (j = 0; j < 3; j++) {
- sspivot(checksh, checkseg);
- if (checkseg.sh != NULL) {
- toptet = neightet;
- while (1) {
- tssbond1(toptet, checkseg);
- fnextself(toptet);
- if (apex(toptet) == apex(neightet)) break;
+ } else {
+ // A faked tet. There is an interior subface. Use it.
+ // See fig/dump-cavity-case19.
+ stpivot(*parysh, neightet);
+ assert(marktested(neightet)); // SELF_CHECK
+ tspivot(neightet, checksh);
+ assert(checksh.sh != NULL); // SELF_CHECK
+ assert(checksh.sh != parysh->sh); // // SELF_CHECK
+ // Align them at the same directed edge.
+ pa = org(neightet);
+ pb = dest(neightet);
+ for (j = 0; j < 3; j++) {
+ if (sorg(checksh) == pa) break;
+ senextself(checksh);
+ }
+ assert(j < 3); // SELF_CHECK
+ if (sdest(checksh) != pb) {
+ senext2self(checksh);
+ sesymself(checksh);
+ }
+ assert(sdest(checksh) == pb); // SELF_CHECK
+ // Bond a subsegment (if it exists).
+ for (j = 0; j < 3; j++) {
+ sspivot(checksh, checkseg);
+ if (checkseg.sh != NULL) {
+ // Let the segment remember an adjacent tet.
+ sstbond(checkseg, neightet);
+ toptet = neightet;
+ while (1) {
+ tssbond1(toptet, checkseg);
+ fnextself(toptet);
+ if (apex(toptet) == apex(neightet)) break;
+ }
}
+ senextself(checksh);
+ enextself(neightet);
}
- senextself(checksh);
- enextself(neightet);
}
}
- }
+ } // if (cavshells != NULL)
}
} else {
// Faces at top and bottom are not matched. There exists non-Delaunay
@@ -2224,22 +2437,33 @@ bool tetgenmesh::fillcavity(arraypool* topshells, arraypool* botshells,
// dummypoint[0] = dummypoint[1] = dummypoint[2] = 0;
ndelaunayedgecount++;
}
- // Unmark all facet vertices.
- for (i = 0; i < facpoints->objects; i++) {
- pf = * (point *) fastlookup(facpoints, i);
- puninfect(pf);
+
+ if (facpoints != NULL) {
+ // Unmark all facet vertices.
+ for (i = 0; i < facpoints->objects; i++) {
+ pf = * (point *) fastlookup(facpoints, i);
+ puninfect(pf);
+ }
}
+
// Delete the temp subfaces.
for (k = 0; k < 2; k++) {
cavshells = (k == 0 ? topshells : botshells);
- for (i = 0; i < cavshells->objects; i++) {
- parysh = (face *) fastlookup(cavshells, i);
- shellfacedealloc(subfacepool, parysh->sh);
+ if (cavshells != NULL) {
+ for (i = 0; i < cavshells->objects; i++) {
+ parysh = (face *) fastlookup(cavshells, i);
+ shellfacedealloc(subfacepool, parysh->sh);
+ }
}
}
+
topshells->restart();
- botshells->restart();
- midfaces->restart();
+ if (botshells != NULL) {
+ botshells->restart();
+ }
+ if (midfaces != NULL) {
+ midfaces->restart();
+ }
// Comment: Now no vertex is marked.
return mflag;
@@ -2261,6 +2485,8 @@ void tetgenmesh::carvecavity(arraypool *crosstets, arraypool *topnewtets,
// NOTE: Some subsegments may contained inside the cavity. They must be
// queued for recovery. See fig/dump-cavity-case20.
+ // Comment: This check should be avoided in the future. Do the check in
+ // routine delaunizecavity() is NOT enough. (2009-04-24).
for (i = 0; i < crosstets->objects; i++) {
parytet = (triface *) fastlookup(crosstets, i);
assert(infected(*parytet)); // SELF_CHECK
@@ -2283,6 +2509,8 @@ void tetgenmesh::carvecavity(arraypool *crosstets, arraypool *topnewtets,
printf(" Queue a missing segment (%d, %d).\n",
pointmark(sorg(checkseg)), pointmark(sdest(checkseg)));
}
+ // Clean the seg-to-tet pointer.
+ stdissolve(checkseg);
sinfect(checkseg);
subsegstack->newindex((void **) &parysh);
*parysh = checkseg;
@@ -2303,17 +2531,22 @@ void tetgenmesh::carvecavity(arraypool *crosstets, arraypool *topnewtets,
// Collect infected new tets in cavity.
for (k = 0; k < 2; k++) {
newtets = (k == 0 ? topnewtets : botnewtets);
- for (i = 0; i < newtets->objects; i++) {
- parytet = (triface *) fastlookup(newtets, i);
- if (infected(*parytet)) {
- crosstets->newindex((void **) &pnewtet);
- *pnewtet = *parytet;
+ if (newtets != NULL) {
+ for (i = 0; i < newtets->objects; i++) {
+ parytet = (triface *) fastlookup(newtets, i);
+ if (infected(*parytet)) {
+ crosstets->newindex((void **) &pnewtet);
+ *pnewtet = *parytet;
+ }
}
}
}
// Collect all new tets in cavity.
for (i = 0; i < crosstets->objects; i++) {
parytet = (triface *) fastlookup(crosstets, i);
+ if (i == 0) {
+ recenttet = *parytet; // Remember a live handle.
+ }
for (j = 0; j < 4; j++) {
decode(parytet->tet[j], neightet);
if (marktested(neightet)) { // Is it a new tet?
@@ -2331,22 +2564,26 @@ void tetgenmesh::carvecavity(arraypool *crosstets, arraypool *topnewtets,
// Delete outer new tets.
for (k = 0; k < 2; k++) {
newtets = (k == 0 ? topnewtets : botnewtets);
- for (i = 0; i < newtets->objects; i++) {
- parytet = (triface *) fastlookup(newtets, i);
- if (infected(*parytet)) {
- // This is an interior tet.
- uninfect(*parytet);
- unmarktest(*parytet);
- } else {
- // An outer tet. Delete it.
- tetrahedrondealloc(parytet->tet);
+ if (newtets != NULL) {
+ for (i = 0; i < newtets->objects; i++) {
+ parytet = (triface *) fastlookup(newtets, i);
+ if (infected(*parytet)) {
+ // This is an interior tet.
+ uninfect(*parytet);
+ unmarktest(*parytet);
+ } else {
+ // An outer tet. Delete it.
+ tetrahedrondealloc(parytet->tet);
+ }
}
}
}
crosstets->restart();
topnewtets->restart();
- botnewtets->restart();
+ if (botnewtets != NULL) {
+ botnewtets->restart();
+ }
}
///////////////////////////////////////////////////////////////////////////////
@@ -2367,6 +2604,9 @@ void tetgenmesh::restorecavity(arraypool *crosstets, arraypool *topnewtets,
for (i = 0; i < crosstets->objects; i++) {
parytet = (triface *) fastlookup(crosstets, i);
assert(infected(*parytet)); // SELF_CHECK
+ if (i == 0) {
+ recenttet = *parytet; // Remember a live handle.
+ }
parytet->ver = 0;
for (parytet->loc = 0; parytet->loc < 4; parytet->loc++) {
symedge(*parytet, neightet);
@@ -2398,14 +2638,18 @@ void tetgenmesh::restorecavity(arraypool *crosstets, arraypool *topnewtets,
tetrahedrondealloc(parytet->tet);
}
- for (i = 0; i < botnewtets->objects; i++) {
- parytet = (triface *) fastlookup(botnewtets, i);
- tetrahedrondealloc(parytet->tet);
+ if (botnewtets != NULL) {
+ for (i = 0; i < botnewtets->objects; i++) {
+ parytet = (triface *) fastlookup(botnewtets, i);
+ tetrahedrondealloc(parytet->tet);
+ }
}
crosstets->restart();
topnewtets->restart();
- botnewtets->restart();
+ if (botnewtets != NULL) {
+ botnewtets->restart();
+ }
}
///////////////////////////////////////////////////////////////////////////////
@@ -2478,27 +2722,23 @@ void tetgenmesh::splitsubedge(point newpt, face *searchsh, arraypool *facfaces,
void tetgenmesh::constrainedfacets()
{
+ /*
arraypool *crosstets, *topnewtets, *botnewtets;
arraypool *topfaces, *botfaces, *midfaces;
arraypool *topshells, *botshells, *facfaces;
arraypool *toppoints, *botpoints, *facpoints;
+ */
triface *parytet, searchtet, neightet;
face *pssub, ssub, neighsh;
face checkseg;
point *ppt, pt, newpt;
enum intersection dir;
bool success, delaunayflag;
- long bakflip22count;
- long cavitycount;
int facetcount;
int bakhullsize;
int s, i, j;
- if (b->verbose) {
- printf(" Constraining facets.\n");
- }
-
- // Initialize arrays.
+ /* // Initialize arrays.
crosstets = new arraypool(sizeof(triface), 10);
topnewtets = new arraypool(sizeof(triface), 10);
botnewtets = new arraypool(sizeof(triface), 10);
@@ -2511,9 +2751,8 @@ void tetgenmesh::constrainedfacets()
facfaces = new arraypool(sizeof(face), 10);
topshells = new arraypool(sizeof(face), 10);
botshells = new arraypool(sizeof(face), 10);
+ */
- bakflip22count = flip22count;
- cavitycount = 0;
facetcount = 0;
// Loop until 'subfacstack' is empty.
@@ -2528,11 +2767,11 @@ void tetgenmesh::constrainedfacets()
if (neightet.tet == NULL) {
// Find an unrecovered subface.
smarktest(ssub);
- facfaces->newindex((void **) &pssub);
+ tg_facfaces->newindex((void **) &pssub);
*pssub = ssub;
// Get all subfaces and vertices of the same facet.
- for (i = 0; i < facfaces->objects; i++) {
- ssub = * (face *) fastlookup(facfaces, i);
+ for (i = 0; i < tg_facfaces->objects; i++) {
+ ssub = * (face *) fastlookup(tg_facfaces, i);
for (j = 0; j < 3; j++) {
sspivot(ssub, checkseg);
if (checkseg.sh == NULL) {
@@ -2543,7 +2782,7 @@ void tetgenmesh::constrainedfacets()
stpivot(neighsh, neightet);
if (neightet.tet == NULL) {
smarktest(neighsh);
- facfaces->newindex((void **) &pssub);
+ tg_facfaces->newindex((void **) &pssub);
*pssub = neighsh;
}
}
@@ -2551,32 +2790,32 @@ void tetgenmesh::constrainedfacets()
pt = sorg(ssub);
if (!pinfected(pt)) {
pinfect(pt);
- facpoints->newindex((void **) &ppt);
+ tg_facpoints->newindex((void **) &ppt);
*ppt = pt;
}
senextself(ssub);
} // j
} // i
// Have found all facet subfaces (vertices). Uninfect them.
- for (i = 0; i < facfaces->objects; i++) {
- pssub = (face *) fastlookup(facfaces, i);
+ for (i = 0; i < tg_facfaces->objects; i++) {
+ pssub = (face *) fastlookup(tg_facfaces, i);
sunmarktest(*pssub);
}
- for (i = 0; i < facpoints->objects; i++) {
- ppt = (point *) fastlookup(facpoints, i);
+ for (i = 0; i < tg_facpoints->objects; i++) {
+ ppt = (point *) fastlookup(tg_facpoints, i);
puninfect(*ppt);
}
if (b->verbose > 1) {
printf(" Recover facet #%d: %ld subfaces, %ld vertices.\n",
- facetcount + 1, facfaces->objects, facpoints->objects);
+ facetcount + 1, tg_facfaces->objects, tg_facpoints->objects);
}
facetcount++;
- // Loop until 'facfaces' is empty.
- while (facfaces->objects > 0l) {
+ // Loop until 'tg_facfaces' is empty.
+ while (tg_facfaces->objects > 0l) {
// Get the last subface of this array.
- facfaces->objects--;
- pssub = (face *) fastlookup(facfaces, facfaces->objects);
+ tg_facfaces->objects--;
+ pssub = (face *) fastlookup(tg_facfaces, tg_facfaces->objects);
ssub = *pssub;
stpivot(ssub, neightet);
@@ -2589,40 +2828,41 @@ void tetgenmesh::constrainedfacets()
assert(dir != COLLISIONFACE); // SELF_CHECK
// Not exist. Push the subface back into stack.
- s = randomnation(facfaces->objects + 1);
- facfaces->newindex((void **) &pssub);
- *pssub = * (face *) fastlookup(facfaces, s);
- * (face *) fastlookup(facfaces, s) = ssub;
+ s = randomnation(tg_facfaces->objects + 1);
+ tg_facfaces->newindex((void **) &pssub);
+ *pssub = * (face *) fastlookup(tg_facfaces, s);
+ * (face *) fastlookup(tg_facfaces, s) = ssub;
if (dir == EDGETRIINT) continue; // All three edges are missing.
// Search for a crossing tet.
- dir = scoutcrosstet(&ssub, &searchtet, facpoints);
+ dir = scoutcrosstet(&ssub, &searchtet, tg_facpoints);
if (dir == ACROSSTET) {
// Recover subfaces by local retetrahedralization.
cavitycount++;
bakhullsize = hullsize;
checksubsegs = checksubfaces = 0;
- crosstets->newindex((void **) &parytet);
+ tg_crosstets->newindex((void **) &parytet);
*parytet = searchtet;
// Form a cavity of crossing tets.
- formcavity(&ssub, crosstets, topfaces, botfaces, toppoints,
- botpoints, facpoints);
+ formcavity(&ssub, tg_crosstets, tg_topfaces, tg_botfaces,
+ tg_toppoints, tg_botpoints, tg_facpoints);
delaunayflag = true;
- // Tetrahedralize the top part. Re-use 'midfaces'.
- success = delaunizecavity(toppoints, topfaces, topshells,
- topnewtets, crosstets, midfaces);
+ // Tetrahedralize the top part. Re-use 'tg_midfaces'.
+ success = delaunizecavity(tg_toppoints, tg_topfaces, tg_topshells,
+ tg_topnewtets, tg_crosstets, tg_midfaces);
if (success) {
- // Tetrahedralize the bottom part. Re-use 'midfaces'.
- success = delaunizecavity(botpoints, botfaces, botshells,
- botnewtets, crosstets, midfaces);
+ // Tetrahedralize the bottom part. Re-use 'tg_midfaces'.
+ success = delaunizecavity(tg_botpoints, tg_botfaces, tg_botshells,
+ tg_botnewtets, tg_crosstets, tg_midfaces);
if (success) {
// Fill the cavity with new tets.
- success = fillcavity(topshells, botshells, midfaces, facpoints);
+ success = fillcavity(tg_topshells, tg_botshells, tg_midfaces,
+ tg_facpoints);
if (success) {
// Delete old tets and outer new tets.
- carvecavity(crosstets, topnewtets, botnewtets);
+ carvecavity(tg_crosstets, tg_topnewtets, tg_botnewtets);
}
} else {
delaunayflag = false;
@@ -2632,16 +2872,16 @@ void tetgenmesh::constrainedfacets()
}
if (!success) {
// Restore old tets and delete new tets.
- restorecavity(crosstets, topnewtets, botnewtets);
+ restorecavity(tg_crosstets, tg_topnewtets, tg_botnewtets);
}
/*if (!delaunayflag) {
dump_facetof(&ssub, "facet1.lua");
while (futureflip != NULL) {
- formedgecavity(futureflip->forg, futureflip->fdest, crosstets,
- topfaces, toppoints);
- crosstets->restart();
- topfaces->restart();
- toppoints->restart();
+ formedgecavity(futureflip->forg, futureflip->fdest, tg_crosstets,
+ tg_topfaces, tg_toppoints);
+ tg_crosstets->restart();
+ tg_topfaces->restart();
+ tg_toppoints->restart();
futureflip = futureflip->nextitem;
}
flippool->restart();
@@ -2654,7 +2894,7 @@ void tetgenmesh::constrainedfacets()
checksubsegs = checksubfaces = 1;
} else if (dir == ACROSSFACE) {
// Recover subfaces by flipping edges in surface mesh.
- recoversubfacebyflips(&ssub, &searchtet, facfaces);
+ recoversubfacebyflips(&ssub, &searchtet, tg_facfaces);
success = true;
} else { // dir == TOUCHFACE
assert(0);
@@ -2662,7 +2902,7 @@ void tetgenmesh::constrainedfacets()
if (!success) break;
} // while
- if (facfaces->objects > 0l) {
+ if (tg_facfaces->objects > 0l) {
// Found a non-Delaunay edge, split it (or a segment close to it).
// Create a new point at the middle of this edge, its coordinates
// were saved in dummypoint in 'fillcavity()'.
@@ -2671,28 +2911,21 @@ void tetgenmesh::constrainedfacets()
setpointtype(newpt, STEINERVERTEX);
dummypoint[0] = dummypoint[1] = dummypoint[2] = 0;
// Insert the new point. Starting search it from 'ssub'.
- splitsubedge(newpt, &ssub, facfaces, facpoints);
- facfaces->restart();
+ splitsubedge(newpt, &ssub, tg_facfaces, tg_facpoints);
+ tg_facfaces->restart();
}
// Clear the list of facet vertices.
- facpoints->restart();
+ tg_facpoints->restart();
// Some subsegments may be queued, recover them.
if (subsegstack->objects > 0l) {
- b->verbose--; // Suppress the message output.
delaunizesegments();
- b->verbose++;
}
// Now the mesh should be constrained Delaunay.
} // if (neightet.tet == NULL)
}
- if (b->verbose) {
- printf(" %ld subedge flips.\n", flip22count - bakflip22count);
- printf(" %ld cavities remeshed.\n", cavitycount);
- }
-
- // Delete arrays.
+ /* // Delete arrays.
delete crosstets;
delete topnewtets;
delete botnewtets;
@@ -2705,26 +2938,37 @@ void tetgenmesh::constrainedfacets()
delete facfaces;
delete topshells;
delete botshells;
+ */
}
///////////////////////////////////////////////////////////////////////////////
// //
-// scoutedge() Search an edge in current tetrahedralization. //
+// scoutsegment2() Search a segment in tetrahedralization. //
+// //
+// Search an edge in tetrahedralization that matches the given segmment. If //
+// such an edge exists, the segment is 'locked' at that edge. 'searchtet' //
+// returns this (constrained) edge. Otherwise, the segment is missing. //
// //
-// The edge [a, b] will be searched. If the edge exists, return a handle in //
-// the tetrahedralization ('searchtet') referring to this edge. If the edge //
-// is missing, and the array 'crosstets' is not a NULL, return the cavity //
-// of all corssing tets of this edge. //
+// If the segment is missing, and the array 'crosstets' is given, it returns //
+// all crossing tetrahedra by this segment. //
+// //
+// The returned value indicates one of the following cases: //
+// - SHAREEDGE, the segment exists and is inserted in T; //
+// - ACROSSVERT, the segment passes a vertex (the origin of 'searchtet'). //
+// - ACROSSEDGE, the segment intersects an edge (in 'searchtet'). //
+// - ACROSSFACE, the segment crosses a face (in 'searchtet'). //
+// - ACROSSSUBSEG, the segment intersects a segment (in 'searchtet'). //
+// - ACROSSSUBFACE, the segment crosses a subface (in 'searchtet'). //
+// - ACROSSTET, the segment is missing and the cavity has formed. //
// //
///////////////////////////////////////////////////////////////////////////////
-enum tetgenmesh::intersection tetgenmesh::scoutedge(point pa, point pb,
- triface* searchtet, arraypool* crosstets, arraypool* cavfaces,
- arraypool* cavpoints)
+enum tetgenmesh::intersection tetgenmesh::scoutsegment2(face* sseg,
+ triface* searchtet, arraypool* crosstets)
{
triface neightet, spintet, *parytet;
face checksh, checkseg;
- point pc, pd, pe, pf, *ppt, *parypt;
+ point pa, pb, pc, pd, pe, pf;// *ppt, *parypt;
enum intersection dir;
int types[2], poss[4];
int pos, i;
@@ -2732,8 +2976,11 @@ enum tetgenmesh::intersection tetgenmesh::scoutedge(point pa, point pb,
tetrahedron ptr;
int *iptr, tver;
+ pa = sorg(*sseg);
+ pb = sdest(*sseg);
+
if (b->verbose > 1) {
- printf(" Search edge (%d, %d).\n", pointmark(pa), pointmark(pb));
+ printf(" Search edge (%d, %d).\n", pointmark(pa), pointmark(pb));
}
// Search a tet whose origin is pa.
@@ -2743,15 +2990,32 @@ enum tetgenmesh::intersection tetgenmesh::scoutedge(point pa, point pb,
dir = finddirection(searchtet, pb);
if (dir == ACROSSVERT) {
if (dest(*searchtet) == pb) {
- // Found! Does there already exist a segment.
+ // Found! Insert the segment.
+ tsspivot(*searchtet, checkseg); // SELF_CHECK
+ if (checkseg.sh == NULL) {
+ // Let the segment remember an adjacent tet.
+ sstbond(*sseg, *searchtet);
+ neightet = *searchtet;
+ do {
+ tssbond1(neightet, *sseg);
+ fnextself(neightet);
+ } while (neightet.tet != searchtet->tet);
+ } else {
+ // Collision! This can happy during facet recovery.
+ // See fig/dump-cavity-case19, -case20.
+ assert(checkseg.sh == sseg->sh); // SELF_CHECK
+ }
return SHAREEDGE; // The edge is not missing.
} else {
+ enextself(*searchtet);
return ACROSSVERT; // The edge intersects a vertex.
}
}
// The edge is missing, shall we form the edge cavity?
- if ((crosstets == NULL) && (cavfaces == NULL) && (cavpoints == NULL)) {
+ if (crosstets == NULL) {
+ // Go to the face/edge it crosses.
+ enextfnextself(*searchtet);
return dir; // ACROSSFACE and ACROSSEDGE
}
@@ -2762,13 +3026,13 @@ enum tetgenmesh::intersection tetgenmesh::scoutedge(point pa, point pb,
infect(*searchtet);
crosstets->newindex((void **) &parytet);
*parytet = *searchtet;
- // Add all vertices of this tet into list.
+ /*// Add all vertices of this tet into list.
ppt = (point *) &(searchtet->tet[4]);
for (i = 0; i < 4; i++) {
pinfect(ppt[i]);
cavpoints->newindex((void **) &parypt);
*parypt = ppt[i];
- }
+ }*/
// Collect all crossing tets of the edge [a, b].
while (1) {
@@ -2783,11 +3047,11 @@ enum tetgenmesh::intersection tetgenmesh::scoutedge(point pa, point pb,
crosstets->newindex((void **) &parytet);
*parytet = *searchtet;
}
- if (!pinfected(pf)) { // Add the opposite point into list.
+ /*if (!pinfected(pf)) { // Add the opposite point into list.
pinfect(pf);
cavpoints->newindex((void **) &parypt);
*parypt = pf;
- }
+ }*/
tspivot(*searchtet, checksh); // Check if a subface is crossed.
if (checksh.sh != NULL) {
// The edge intersect a subface.
@@ -2806,16 +3070,17 @@ enum tetgenmesh::intersection tetgenmesh::scoutedge(point pa, point pb,
*parytet = spintet;
}
pd = oppo(spintet);
- if (!pinfected(pd)) { // Add the opposite point into list.
+ /*if (!pinfected(pd)) { // Add the opposite point into list.
pinfect(pd);
cavpoints->newindex((void **) &parypt);
*parypt = pd;
- }
+ }*/
if (apex(spintet) == pc) break;
}
tsspivot(spintet, checkseg); // Check if a segment is crossed.
if (checkseg.sh != NULL) {
// The edge intersects a subsegment.
+ *searchtet = spintet;
dir = ACROSSSUBSEG;
break;
}
@@ -2900,16 +3165,16 @@ enum tetgenmesh::intersection tetgenmesh::scoutedge(point pa, point pb,
parytet = (triface *) fastlookup(crosstets, i);
uninfect(*parytet);
}
- for (i = 0; i < cavpoints->objects; i++) {
+ /*for (i = 0; i < cavpoints->objects; i++) {
parypt = (point *) fastlookup(cavpoints, i);
puninfect(*parypt);
- }
+ }*/
crosstets->restart();
- cavpoints->restart();
+ // cavpoints->restart();
return dir;
}
- //We can form the edge cavity.
+ /*//We can form the edge cavity.
for (i = 0; i < crosstets->objects; i++) {
parytet = (triface *) fastlookup(crosstets, i);
*searchtet = *parytet;
@@ -2922,11 +3187,55 @@ enum tetgenmesh::intersection tetgenmesh::scoutedge(point pa, point pb,
}
}
}
+ // Uninfect the vertices.
+ for (i = 0; i < cavpoints->objects; i++) {
+ parypt = (point *) fastlookup(cavpoints, i);
+ puninfect(*parypt);
+ }
+ // Comment: All crossing tets are infected.
+ */
+
+ if (b->verbose > 1) {
+ printf(" Formed edge cavity: %ld tets.\n", crosstets->objects);
+ }
+
+ return ACROSSTET;
+}
- // Uninfect the collected crossing tets and vertices.
+///////////////////////////////////////////////////////////////////////////////
+// //
+// tetrasegcavity() Tetrahedralize a cavity for recovering a segment. //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+bool tetgenmesh::tetrasegcavity(face* sseg, arraypool* crosstets,
+ arraypool* cavpoints, arraypool* cavfaces, arraypool* cavshells,
+ arraypool* newtets, arraypool* misfaces, arraypool* missegs,
+ arraypool* fixedseglist)
+{
+ triface searchtet, neightet, spintet, *parytet, *parytet1;
+ face tmpsh, *parysh;
+ face checkseg, *paryseg;
+ point pt[4], pswap, *parypt;
+ enum intersection dir;
+ REAL ori;
+ bool success;
+ int i, j;
+
+ tetrahedron ptr;
+ int *iptr, tver;
+
+ // cavpoints are collected.
for (i = 0; i < crosstets->objects; i++) {
parytet = (triface *) fastlookup(crosstets, i);
- uninfect(*parytet);
+ for (j = 0; j < 4; j++) {
+ pt[0] = (point) parytet->tet[4 + j];
+ if (!pinfected(pt[0])) {
+ pinfect(pt[0]);
+ cavpoints->newindex((void **) &parypt);
+ *parypt = pt[0];
+ }
+ }
}
for (i = 0; i < cavpoints->objects; i++) {
parypt = (point *) fastlookup(cavpoints, i);
@@ -2934,11 +3243,533 @@ enum tetgenmesh::intersection tetgenmesh::scoutedge(point pa, point pb,
}
if (b->verbose > 1) {
- printf(" Formed edge cavity: %ld tets, %ld faces, %ld vertices.\n",
- crosstets->objects, cavfaces->objects, cavpoints->objects);
+ printf(" Tetrahedralizing segment cavity: %ld points.\n",
+ cavpoints->objects, cavfaces->objects);
}
- return ACROSSFACE; // or ACROSSEDGE
+ // Form an initial constrained tetrahedralization (CT) which has only one
+ // tetrahedron containing the given segment.
+
+ // The first two points are the endpoints of the segment.
+ pt[0] = sorg(*sseg);
+ pt[1] = sdest(*sseg);
+ pinfect(pt[0]);
+ pinfect(pt[1]);
+
+ // Get the third and fourth points.
+ for (i = 0; i < cavpoints->objects; i++) {
+ parypt = (point *) fastlookup(cavpoints, i);
+ if (!pinfected(*parypt)) {
+ pt[2] = *parypt;
+ i++;
+ for (; i < cavpoints->objects; i++) {
+ parypt = (point *) fastlookup(cavpoints, i);
+ if (!pinfected(*parypt)) {
+ ori = orient3d(pt[0], pt[1], pt[2], *parypt);
+ if (ori != 0) {
+ pt[3] = *parypt;
+ if (ori > 0) { // Swap pa and pb.
+ pswap = pt[0]; pt[0] = pt[1]; pt[1] = pswap;
+ }
+ break;
+ }
+ }
+ } // i
+ break;
+ }
+ } // i
+ assert(i < cavpoints->objects); // SELF_CHECK
+ pinfect(pt[2]);
+ pinfect(pt[3]);
+
+ // Create an initial DT.
+ initialDT(pt[0], pt[1], pt[2], pt[3]);
+
+ // Insert the segment (turns a DT into a CT).
+ scoutsegment2(sseg, &searchtet, NULL);
+
+ // Insert the other points into the CT (the segment is respected).
+ for (i = 0; i < cavpoints->objects; i++) {
+ parypt = (point *) fastlookup(cavpoints, i);
+ if (!pinfected(*parypt)) {
+ // pinfect(*parypt);
+ searchtet = recenttet; // No random samples.
+ // Insert the point by flips. Set 'flipflag' = 3, do not flip a segment.
+ flipinsertvertex(*parypt, &searchtet, 3);
+ } else {
+ puninfect(*parypt);
+ }
+ }
+ // Comment: All vertices of the cavity are NOT marked.
+
+ success = true;
+
+ while (1) {
+
+ // cavfaces are collected.
+ for (i = 0; i < crosstets->objects; i++) {
+ parytet = (triface *) fastlookup(crosstets, i);
+ searchtet = *parytet;
+ for (searchtet.loc = 0; searchtet.loc < 4; searchtet.loc++) {
+ sym(searchtet, neightet);
+ if (!infected(neightet)) {
+ // A cavity bounday face.
+ cavfaces->newindex((void **) &parytet);
+ *parytet = neightet;
+ }
+ }
+ }
+
+ // Identify boundary faces. Mark interior tets. Save missing faces.
+ for (i = 0; i < cavfaces->objects; i++) {
+ parytet = (triface *) fastlookup(cavfaces, i);
+ // The tet of this face is exteriorly adjacent to the cavity.
+ assert(!infected(*parytet)); // SELF_CHECK
+ // This face may contain dummypoint (See fig/dump-cavity-case2).
+ // If so, dummypoint must be its apex.
+ parytet->ver = 4;
+ pt[0] = org(*parytet);
+ pt[1] = dest(*parytet);
+ pt[2] = apex(*parytet);
+ // Create a temp subface.
+ makeshellface(subfacepool, &tmpsh);
+ setshvertices(tmpsh, pt[0], pt[1], pt[2]);
+ // Insert tmpsh in CT.
+ searchtet.tet = NULL;
+ dir = scoutsubface(&tmpsh, &searchtet);
+ if (dir == SHAREFACE) {
+ // Identify the inter and outer tets at tempsh.
+ stpivot(tmpsh, neightet);
+ // neightet and tmpsh refer to the same edge [pt[0], pt[1]].
+ // Moreover, neightet is in 0th edge ring (see decode()).
+ if (org(neightet) != pt[1]) {
+ symedgeself(neightet);
+ assert(org(neightet) == pt[1]); // SELF_CHECK
+ // Make sure that tmpsh is connected with an interior tet.
+ tsbond(neightet, tmpsh);
+ }
+ assert(dest(neightet) == pt[0]); // SELF_CHECK
+ } else if (dir == COLLISIONFACE) {
+ // A subface is already inserted.
+ assert(0); // Not handled yet.
+ } else {
+ if (b->verbose > 1) {
+ printf(" p:draw_subface(%d, %d, %d) -- %d is missing\n",
+ pointmark(pt[0]), pointmark(pt[1]), pointmark(pt[2]), i);
+ }
+ shellfacedealloc(subfacepool, tmpsh.sh);
+ // Save this face in list.
+ misfaces->newindex((void **) &parytet1);
+ *parytet1 = *parytet;
+ continue;
+ }
+ // Remember tmpsh (use the adjacent tet slot).
+ // parytet->tet[parytet->loc] = (tetrahedron) sencode(tmpsh);
+ tmpsh.sh[0] = (shellface) encode(*parytet);
+ // Save this subface.
+ cavshells->newindex((void **) &parysh);
+ *parysh = tmpsh;
+ } // for (i)
+
+ if (misfaces->objects > 0) {
+ // Removing tempoaray subfaces.
+ for (i = 0; i < cavshells->objects; i++) {
+ parysh = (face *) fastlookup(cavshells, i);
+ stpivot(*parysh, neightet);
+ uninfect(neightet);
+ tsdissolve(neightet); // Detach it from adj. tets.
+ symself(neightet);
+ tsdissolve(neightet);
+ shellfacedealloc(subfacepool, parysh->sh);
+ }
+ cavshells->restart();
+
+ for (i = 0; i < cavpoints->objects; i++) {
+ parypt = (point *) fastlookup(cavpoints, i);
+ pinfect(*parypt);
+ }
+
+ // Enlarge the cavity.
+ for (i = 0; i < misfaces->objects; i++) {
+ // Get a missing face.
+ parytet = (triface *) fastlookup(misfaces, i);
+ // For this routine we do not check subface(s).
+ if (!infected(*parytet)) {
+ // This face should not be on the hull.
+ assert((point) parytet->tet[7] != dummypoint); // SELF_CHECK
+ // Check if we can enclose this tet into our cavity.
+ infect(*parytet);
+ // Check its six edges for enclosed segments.
+ neightet.tet = parytet->tet;
+ for (j = 0; j < 6; j++) {
+ neightet.loc = edge2locver[j][0];
+ neightet.ver = edge2locver[j][1];
+ tsspivot(neightet, checkseg);
+ if (checkseg.sh != NULL) {
+ // Check if this segment is inside our cavity.
+ spintet = neightet;
+ while (1) {
+ fnextself(spintet);
+ if (!infected(spintet)) break; // Not inside.
+ if (apex(spintet) == apex(neightet)) {
+ if (b->verbose > 1) {
+ printf(" p:draw_subseg(%d, %d) -- is inside.\n",
+ pointmark(sorg(checkseg)), pointmark(sdest(checkseg)));
+ }
+ if (!smarktested(checkseg)) {
+ missegs->newindex((void **) &parysh);
+ *parysh = checkseg;
+ } else {
+ if (b->verbose > 1) {
+ printf(" !! A fixed segment.\n");
+ }
+ success = false;
+ }
+ break;
+ }
+ } // while (1)
+ if (!success) break;
+ }
+ } // j
+ if (success) {
+ // We can enlarge the cavity.
+ if (b->verbose > 1) {
+ printf(" Add a crosstet (%d, %d, %d, %d).\n",
+ pointmark(org(*parytet)), pointmark(dest(*parytet)),
+ pointmark(apex(*parytet)), pointmark(oppo(*parytet)));
+ }
+ crosstets->newindex((void **) &parytet1);
+ *parytet1 = *parytet;
+ // Insert the opposite point if it is not in CT.
+ pt[0] = oppo(*parytet);
+ if (!pinfected(pt[0])) {
+ if (b->verbose > 1) {
+ printf(" Insert oppo-point %d.\n", pointmark(pt[0]));
+ }
+ pinfect(pt[0]);
+ cavpoints->newindex((void **) &parypt);
+ *parypt = pt[0];
+ searchtet = recenttet;
+ flipinsertvertex(pt[0], &searchtet, 3);
+ }
+ /*// Add three opposite faces into the boundary list.
+ for (j = 0; j < 3; j++) {
+ enext0fnext(*parytet, neightet);
+ symself(neightet);
+ if (!infected(neightet)) {
+ if (b->verbose > 1) {
+ printf(" Add a cavface (%d, %d, %d).\n",
+ pointmark(org(neightet)), pointmark(dest(neightet)),
+ pointmark(apex(neightet)));
+ }
+ cavfaces->newindex((void **) &parytet1);
+ *parytet1 = neightet;
+ } else {
+ // It is an interior face, we do not check subface
+ // for this routine.
+ }
+ enextself(*parytet);
+ } // j */
+ } else {
+ // Do not enlarge it due to an existing segment.
+ uninfect(*parytet);
+ }
+ } // if (!infected(*parytet))
+ if (!success) break;
+ } // i
+
+ for (i = 0; i < cavpoints->objects; i++) {
+ parypt = (point *) fastlookup(cavpoints, i);
+ puninfect(*parypt);
+ }
+
+ misfaces->restart();
+ cavfaces->restart();
+
+ if (success && (missegs->objects > 0l)) {
+ // The cavity has been enlarged, but some segments are enclosed.
+ if (!smarktested(*sseg)) {
+ smarktest(*sseg);
+ fixedseglist->newindex((void **) &paryseg);
+ *paryseg = *sseg;
+ }
+ // SELF_CHECK
+ assert(cavshells->objects == 0l);
+ assert(newtets->objects == 0l);
+ // Recover the missing segments inside the enlarged cavity.
+ success = recoversegments(fixedseglist, missegs, cavfaces, cavshells,
+ misfaces, newtets);
+ // SELF_CHECK
+ assert(missegs->objects == 0l);
+ assert(cavshells->objects == 0l);
+ assert(newtets->objects == 0l);
+ assert(cavfaces->objects == 0l);
+ }
+
+ if (success) {
+ // Cavity has enlarged. Continue to recover the segment.
+ cavityexpcount++;
+ continue;
+ }
+ } // if (misfaces->objects > 0)
+
+ break; // Leave the loop.
+
+ } // while (1)
+
+ // Collect all tets of the DT.
+ marktest(recenttet);
+ newtets->newindex((void **) &parytet);
+ *parytet = recenttet;
+ for (i = 0; i < newtets->objects; i++) {
+ searchtet = * (triface *) fastlookup(newtets, i);
+ for (searchtet.loc = 0; searchtet.loc < 4; searchtet.loc++) {
+ sym(searchtet, neightet);
+ if (!marktested(neightet)) {
+ marktest(neightet);
+ newtets->newindex((void **) &parytet);
+ *parytet = neightet;
+ }
+ }
+ }
+ // Comment: All new tets are marktested.
+
+ cavpoints->restart(); // Comment: Now no vertex is marked.
+ cavfaces->restart();
+
+ return success;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// //
+// recoversegments() Recover segments by local remeshing. //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+bool tetgenmesh::recoversegments(arraypool* fixedseglist, arraypool *recosegs,
+ arraypool* cavfaces, arraypool* cavshells, arraypool* misfaces,
+ arraypool* newtets)
+{
+ arraypool *crosstets;
+ arraypool *missegs;
+ arraypool *cavpoints;
+ triface searchtet;
+ face *psseg, sseg;
+ enum intersection dir;
+ bool success;
+ long bakhullsize;
+ long cavitycount;
+ int i;
+
+ if (b->verbose > 1) {
+ printf(" Recovering %ld segments (%ld fixed segments).\n",
+ recosegs->objects, fixedseglist->objects);
+ }
+
+ // Initialize arrays.
+ crosstets = new arraypool(sizeof(triface), 10);
+ missegs = new arraypool(sizeof(face), 8);
+ cavpoints = new arraypool(sizeof(point), 8);
+
+ success = true;
+
+ // Loop until 'recosegs' is empty.
+ while (recosegs->objects > 0l) {
+ // seglist is used as a stack.
+ recosegs->objects--;
+ psseg = (face *) fastlookup(recosegs, recosegs->objects);
+ sseg = *psseg;
+
+ // It is not a global missing segment.
+ assert(!sinfected(sseg));
+
+ // Form the segment cavity.
+ searchtet.tet = NULL;
+ dir = scoutsegment2(&sseg, &searchtet, crosstets);
+
+ // The segment is missing.
+ assert(dir != SHAREEDGE); // SELF_CHECK
+
+ if (dir == ACROSSTET) {
+ // Recover the segment by local re-meshing.
+ bakhullsize = hullsize;
+ success = tetrasegcavity(&sseg, crosstets, cavpoints, cavfaces,
+ cavshells, newtets, misfaces, missegs, fixedseglist);
+ // Do not clean 'fixedseglist'.
+ hullsize = bakhullsize;
+ if (success) {
+ // The segment is recovered.
+ fillcavity(cavshells, NULL, NULL, NULL);
+ carvecavity(crosstets, newtets, NULL);
+ } else {
+ // Unable to recover the segment.
+ restorecavity(crosstets, newtets, NULL);
+ break;
+ }
+ } else {
+ // Unexpected return type.
+ assert(0); // Not handled yet.
+ }
+ }
+
+ if (recosegs->objects > 0l) {
+ recosegs->restart();
+ }
+
+ delete crosstets;
+ delete missegs;
+ delete cavpoints;
+
+ return success;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// //
+// constrainedsegments() Recover segments in a constrained //
+// tetrahedralization. //
+// //
+// 'recoverseglist' contains a list of recovering segments, 'fixedseglist' //
+// is an accumulated list of segments which are inside current tetrahedrali- //
+// zation and must "fix" at their places (do not remove them). All segments //
+// in 'fixedseglist' are smarktested. //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+void tetgenmesh::constrainedsegments()
+{
+ arraypool *crosstets, *newtets, *misfaces;
+ arraypool *cavfaces;
+ arraypool *cavpoints;
+ arraypool *cavshells;
+ arraypool *fixedseglist, *missegs;
+ triface searchtet;
+ face splitsh;
+ face *psseg, sseg;
+ point newpt, pa, pb;
+ enum intersection dir;
+ bool success;
+ long bakhullsize;
+ long cavitycount;
+ long steinptcount;
+ int i;
+
+ if (b->verbose) {
+ printf(" Recovering %ld segments.\n", subsegstack->objects);
+ }
+
+ // Initialize arrays.
+ crosstets = new arraypool(sizeof(triface), 10);
+ newtets = new arraypool(sizeof(triface), 10);
+ misfaces = new arraypool(sizeof(triface), 10);
+ cavfaces = new arraypool(sizeof(triface), 10);
+ cavshells = new arraypool(sizeof(face), 10);
+ cavpoints = new arraypool(sizeof(point), 8);
+ fixedseglist = new arraypool(sizeof(face), 8);
+ missegs = new arraypool(sizeof(face), 8);
+
+ cavitycount = 0l;
+ steinptcount = 0l;
+
+ // Loop until 'subsegstack' is empty.
+ while (subsegstack->objects > 0l) {
+ // seglist is used as a stack.
+ subsegstack->objects--;
+ psseg = (face *) fastlookup(subsegstack, subsegstack->objects);
+ sseg = *psseg;
+
+ if (!sinfected(sseg)) continue; // Not a missing segment.
+ suninfect(sseg);
+
+ // Insert the segment.
+ searchtet.tet = NULL;
+ dir = scoutsegment2(&sseg, &searchtet, crosstets);
+
+ if (dir != SHAREEDGE) {
+ // The segment is missing.
+ if (dir == ACROSSTET) {
+ // Recover the segment by local re-meshing.
+ bakhullsize = hullsize;
+ success = tetrasegcavity(&sseg, crosstets, cavpoints, cavfaces,
+ cavshells, newtets, misfaces, missegs, fixedseglist);
+ hullsize = bakhullsize;
+ // Unmark the testmarked segments.
+ for (i = 0; i < fixedseglist->objects; i++) {
+ psseg = (face *) fastlookup(fixedseglist, i);
+ assert(smarktested(*psseg));
+ sunmarktest(*psseg);
+ }
+ fixedseglist->restart(); // Clear this list.
+ if (success) {
+ // The segment is recovered.
+ fillcavity(cavshells, NULL, NULL, NULL);
+ carvecavity(crosstets, newtets, NULL);
+ cavitycount++;
+ } else {
+ // Unable to recover the segment.
+ restorecavity(crosstets, newtets, NULL);
+ /*// Split the segment at its middle.
+ makepoint(&newpt);
+ pa = sorg(sseg);
+ pb = sdest(sseg);
+ for (i = 0; i < 3; i++) {
+ newpt[i] = 0.5 * (pa[i] + pb[i]);
+ }
+ setpointtype(newpt, STEINERVERTEX);
+ // Insert the point into the surface mesh.
+ spivot(sseg, splitsh);
+ // Two subsegments are queued in 'subsegstack' for recovery.
+ sinsertvertex(newpt, &splitsh, &sseg, true, false);
+ // Insert the point into the CT.
+ point2tetorg(pa, searchtet);
+ // Set 'flipflag' = 3, do not flip a segment.
+ flipinsertvertex(newpt, &searchtet, 3);
+ */
+ steinptcount++;
+ }
+ } else if (dir == ACROSSVERT) {
+ printf("Error: Invalid PLC! A point and a segment intersect.\n");
+ pa = farsorg(sseg);
+ pb = farsdest(sseg);
+ printf(" Point: %d. Segment: (%d, %d).\n", pointmark(org(searchtet)),
+ pointmark(pa), pointmark(pb));
+ terminatetetgen(1);
+ } else if (dir == ACROSSSUBSEG) {
+ printf("Error: Invalid PLC! Two segments intersect.\n");
+ pa = farsorg(sseg);
+ pb = farsdest(sseg);
+ printf(" 1st: (%d, %d)", pointmark(pa), pointmark(pb));
+ tsspivot(searchtet, sseg);
+ assert(sseg.sh != NULL);
+ pa = farsorg(sseg);
+ pb = farsdest(sseg);
+ printf(" 2nd: (%d, %d).\n", pointmark(pa), pointmark(pb));
+ terminatetetgen(1);
+ } else if (dir == ACROSSSUBFACE) {
+ printf("Error: Invalid PLC! A segment and a subface intersect.\n");
+ pa = farsorg(sseg);
+ pb = farsdest(sseg);
+ printf(" Segment: (%d, %d)", pointmark(pa), pointmark(pb));
+ tspivot(searchtet, splitsh);
+ printf(" Subface: (%d, %d, %d)", pointmark(sorg(splitsh)),
+ pointmark(sdest(splitsh)), pointmark(sapex(splitsh)));
+ terminatetetgen(1);
+ }
+ }
+ }
+
+ if (b->verbose) {
+ printf(" %ld cavities remeshed.\n", cavitycount);
+ printf(" %ld Steiner points inserted.\n", steinptcount);
+ }
+
+ delete crosstets;
+ delete newtets;
+ delete misfaces;
+ delete cavfaces;
+ delete cavshells;
+ delete cavpoints;
+ delete fixedseglist;
+ delete missegs;
}
///////////////////////////////////////////////////////////////////////////////
@@ -2953,6 +3784,8 @@ void tetgenmesh::formskeleton()
{
face *pssub, ssub;
REAL bakeps;
+ long bakflip22count;
+ long bakcavitycount;
int s, i;
if (!b->quiet) {
@@ -2991,8 +3824,22 @@ void tetgenmesh::formskeleton()
// Segments will be introduced.
checksubsegs = 1;
+
// Recover segments.
- delaunizesegments();
+ if (b->nobisect == 0) {
+ if (b->verbose) {
+ printf(" Delaunizing segments.\n");
+ }
+
+ delaunizesegments();
+
+ if (b->verbose) {
+ printf(" %d protecting points.\n", r1count + r2count + r3count);
+ }
+ } else {
+ // -Y option, constrained recover.
+ constrainedsegments();
+ }
// Randomly order the subfaces.
subfacepool->traversalinit();
@@ -3009,9 +3856,49 @@ void tetgenmesh::formskeleton()
// Subfaces will be introduced.
checksubfaces = 1;
+ bakflip22count = flip22count;
+ bakcavitycount = cavitycount;
+
+ if (b->verbose) {
+ printf(" Constraining facets.\n");
+ }
+
+ // Initialize arrays.
+ tg_crosstets = new arraypool(sizeof(triface), 10);
+ tg_topnewtets = new arraypool(sizeof(triface), 10);
+ tg_botnewtets = new arraypool(sizeof(triface), 10);
+ tg_topfaces = new arraypool(sizeof(triface), 10);
+ tg_botfaces = new arraypool(sizeof(triface), 10);
+ tg_midfaces = new arraypool(sizeof(triface), 10);
+ tg_toppoints = new arraypool(sizeof(point), 8);
+ tg_botpoints = new arraypool(sizeof(point), 8);
+ tg_facpoints = new arraypool(sizeof(point), 8);
+ tg_facfaces = new arraypool(sizeof(face), 10);
+ tg_topshells = new arraypool(sizeof(face), 10);
+ tg_botshells = new arraypool(sizeof(face), 10);
+
// Recover facets.
constrainedfacets();
+ // Delete arrays.
+ delete tg_crosstets;
+ delete tg_topnewtets;
+ delete tg_botnewtets;
+ delete tg_topfaces;
+ delete tg_botfaces;
+ delete tg_midfaces;
+ delete tg_toppoints;
+ delete tg_botpoints;
+ delete tg_facpoints;
+ delete tg_facfaces;
+ delete tg_topshells;
+ delete tg_botshells;
+
+ if (b->verbose) {
+ printf(" %ld subedge flips.\n", flip22count - bakflip22count);
+ printf(" %ld cavities remeshed.\n", cavitycount - bakcavitycount);
+ }
+
// checksubsegs = 0;
b->epsilon = bakeps;
}
@@ -3108,7 +3995,7 @@ void tetgenmesh::carveholes()
}
}
- // Find and infect all exterior tets.
+ // Find and infect all exterior tets (in concave place and in holes).
for (i = 0; i < tetarray->objects; i++) {
parytet = (triface *) fastlookup(tetarray, i);
tetloop = *parytet;
@@ -3177,7 +4064,8 @@ void tetgenmesh::carveholes()
tetarray->restart(); // Re-use it for new hull tets.
- // Create new hull faces and update the point-to-tet map.
+ // Create new hull tets.
+ // Update point-to-tet map, segment-to-tet map, and subface-to-tet map.
tetrahedronpool->traversalinit();
tetloop.ver = 0;
tetloop.tet = tetrahedrontraverse();
@@ -3193,7 +4081,18 @@ void tetgenmesh::carveholes()
pc = apex(tetloop);
setvertices(hulltet, pb, pa, pc, dummypoint);
bond(tetloop, hulltet);
+ // Update subface-to-tet map.
tsbond(hulltet, checksh);
+ // Update segment-to-tet map.
+ for (i = 0; i < 3; i++) {
+ tsspivot(tetloop, checkseg);
+ if (checkseg.sh != NULL) {
+ tssbond1(hulltet, checkseg);
+ sstbond(checkseg, hulltet);
+ }
+ enextself(tetloop);
+ enext2self(hulltet);
+ }
// Save this hull tet in list.
tetarray->newindex((void **) &parytet);
*parytet = hulltet;
@@ -3234,6 +4133,10 @@ void tetgenmesh::carveholes()
//////////////////////////////////////////////////////////////////////
// Peel off "flat" tetrahedra at boundary.
+ //
+ // A tet is flat if it contains two subfaces of the same facet.
+ // Flat tets are possible when a facet is defined by non-exactly
+ // coplanar vertices.
tetarray->restart(); // Re-use this array.
flatcount = 0;
diff --git a/contrib/Tetgen/delaunay.cxx b/contrib/Tetgen/delaunay.cxx
index d97a0101e8b9387e356a2789ce69d7bdfeba71d6..a023c94b201a1c2bc4c9e64930f56ab40ed8475d 100644
--- a/contrib/Tetgen/delaunay.cxx
+++ b/contrib/Tetgen/delaunay.cxx
@@ -517,9 +517,9 @@ void tetgenmesh::initialDT(point pa, point pb, point pc, point pd)
// If 'visflag' is TRUE, force to check the visibility of the boundary faces //
// of cavity. This is needed when T is not Delaunay. //
// //
-// If 'noencflag' is TRUE, only insert the new point p if it does not cause //
-// any existing (sub)segment be non-Delaunay. This option only is checked //
-// when the global variable 'checksubsegs' is set. //
+// If 'noencsegflag' is TRUE, only insert the point if it does not encroach //
+// on any existing segment of the mesh. Otherwise, do not insert the point, //
+// and all encroaching segments are returned in subsegstack. //
// //
///////////////////////////////////////////////////////////////////////////////
@@ -763,8 +763,9 @@ enum tetgenmesh::location tetgenmesh::insertvertex(point insertpt,
maxbowatcavsize = cavebdrylist->objects;
}
- if (checksubsegs || noencsegflag) {
- // Check if some (sub)segments are inside the cavity.
+ if (checksubsegs) {
+ // Check if some (sub)segments are inside the cavity. Such segments
+ // are queued in 'subsegstack'.
for (i = 0; i < caveoldtetlist->objects; i++) {
cavetet = (triface *) fastlookup(caveoldtetlist, i);
for (j = 0; j < 6; j++) {
@@ -788,6 +789,9 @@ enum tetgenmesh::location tetgenmesh::insertvertex(point insertpt,
printf(" Queue a missing segment (%d, %d).\n",
pointmark(sorg(sseg)), pointmark(sdest(sseg)));
}
+ // All tets containing this segment will be dead, clean the
+ // seg-to-tet pointer.
+ stdissolve(sseg);
sinfect(sseg); // Only save it once.
subsegstack->newindex((void **) &psseg);
*psseg = sseg;
@@ -795,6 +799,35 @@ enum tetgenmesh::location tetgenmesh::insertvertex(point insertpt,
}
}
}
+ if (noencsegflag) {
+ // Check for encroaching segment on the boundary of the cavity.
+ // Encroached segments are queued in 'subsegstack'.
+ for (i = 0; i < cavebdrylist->objects; i++) {
+ cavetet = (triface *) fastlookup(cavebdrylist, i);
+ // 'cavetet' is an exterior tet adjacent to the cavity.
+ assert(cavetet->ver == 4); // SELF_CHECK
+ for (j = 0; j < 3; j++) {
+ tsspivot(*cavetet, sseg);
+ if (sseg.sh != NULL) {
+ // Found a segment. Check it if it is not queued yet.
+ if (!sinfected(sseg)) {
+ if (checkedge4encroach(sseg, insertpt, 0)) {
+ if (b->verbose > 1) {
+ printf(" Queue an encroaching segment (%d, %d).\n",
+ pointmark(sorg(sseg)), pointmark(sdest(sseg)));
+ }
+ // This segment will still be connected to a tet after the
+ // insertion.
+ sinfect(sseg); // Only save it once.
+ subsegstack->newindex((void **) &psseg);
+ *psseg = sseg;
+ }
+ }
+ }
+ enextself(*cavetet);
+ }
+ }
+ }
}
if (noencsegflag && (subsegstack->objects > 0)) {
@@ -818,8 +851,9 @@ enum tetgenmesh::location tetgenmesh::insertvertex(point insertpt,
return ENCSEGMENT;
}
- if (checksubfaces || noencsubflag) {
- // Check if some subfaces are inside the cavity.
+ if (checksubfaces) {
+ // Check if some subfaces are inside the cavity. Such subfaces
+ // are queued in 'subfacstack'.
for (i = 0; i < caveoldtetlist->objects; i++) {
cavetet = (triface *) fastlookup(caveoldtetlist, i);
neightet.tet = cavetet->tet;
@@ -843,6 +877,19 @@ enum tetgenmesh::location tetgenmesh::insertvertex(point insertpt,
}
}
}
+ if (noencsubflag) {
+ // Check for encroaching subface on the boundary of the cavity.
+ // Encroached subfaces are queued in 'subfacstack'.
+ for (i = 0; i < cavebdrylist->objects; i++) {
+ cavetet = (triface *) fastlookup(cavebdrylist, i);
+ // 'cavetet' is an exterior tet adjacent to the cavity.
+ assert(cavetet->ver == 4); // SELF_CHECK
+ tspivot(*cavetet, checksh);
+ if (checksh.sh != NULL) {
+ // checkface4encroach();
+ }
+ }
+ }
}
if (noencsubflag && (subfacstack->objects > 0)) {
@@ -869,6 +916,8 @@ enum tetgenmesh::location tetgenmesh::insertvertex(point insertpt,
if (visflag) {
// If T is not a Delaunay triangulation, the formed cavity may not be
// star-shaped (fig/dump-cavity-case8). Validation is needed.
+ // Comment: The validation is done by removing tets from the cavity
+ // until the cavity is star-shaped.
cavetetlist->restart(); // Re-use it.
for (i = 0; i < cavebdrylist->objects; i++) {
cavetet = (triface *) fastlookup(cavebdrylist, i);
@@ -881,7 +930,7 @@ enum tetgenmesh::location tetgenmesh::insertvertex(point insertpt,
pb = dest(*cavetet);
pc = apex(*cavetet);
ori = orient3d(pa, pb, pc, insertpt); orient3dcount++;
- assert(ori != 0.0); // SELF_CHECK
+ // assert(ori != 0.0); // SELF_CHECK
enqflag = (ori > 0.0);
} else {
enqflag = true; // A hull face.
@@ -996,6 +1045,12 @@ enum tetgenmesh::location tetgenmesh::insertvertex(point insertpt,
suninfect(sseg); // Dequeue a non-missing segment.
}
tssbond1(newtet, sseg);
+ sstbond(sseg, newtet);
+ // Do we need to care about encroached segments?
+ if ((badsegpool != NULL) && !smarktested(sseg)) {
+ // Queue the subsegment if it is encroached by the new point.
+ checkedge4encroach(sseg, insertpt, 1);
+ }
}
enextself(neightet);
enext2self(newtet);
@@ -1005,6 +1060,11 @@ enum tetgenmesh::location tetgenmesh::insertvertex(point insertpt,
tspivot(neightet, checksh);
if (checksh.sh != NULL) {
tsbond(newtet, checksh); // Also disconnect the old bond.
+ // Do we need to care about encroached segments?
+ if ((badsubpool != NULL) && !smarktested(checksh)) {
+ // Queue the subface if it is encroached by the new point.
+ // checkface4encroach(checksh, insertpt, 1);
+ }
}
}
if (updatecount > 0l) {
@@ -1165,7 +1225,13 @@ enum tetgenmesh::location tetgenmesh::insertvertex(point insertpt,
futureflip = flippush(futureflip, parytet, pa);
}
// Recover Delaunay faces.
- // Set 'flipflag' = 2, s.t. all faces are checked for flipping.
+ // Set 'flipflag' = 2, s.t. all non-flippable faces are not ignoring,
+ // they are queued in the flip queue for the future flips. One key
+ // is that we also flip non-Delaunay segments and subfaces. This way
+ // we are able to recover all non-Delaunay edges/faces (no proof yet),
+ // i.e., the flip will terminate.
+ // The flipped segments and subfaces are queued in 'subsegstack' and
+ // 'subfacstack' for recovery.
lawsonflip3d(2);
}
diff --git a/contrib/Tetgen/flip.cxx b/contrib/Tetgen/flip.cxx
index 12b0ba76c15c49c4df0e35bf9c0018049feadc47..d32dd83fc9a3d961d3470c536a13e086ac060a29 100644
--- a/contrib/Tetgen/flip.cxx
+++ b/contrib/Tetgen/flip.cxx
@@ -573,18 +573,22 @@ void tetgenmesh::flip14(point newpt, triface* splittet, int flipflag)
if (checkseg.sh != NULL) {
tssbond1(*splittet, checkseg);
tssbond1(newface, checkseg);
+ // Let the segment remember an adjacent tet.
+ sstbond(checkseg, *splittet);
}
enextself(newface); // [b,d], [c,d], [a,d]
enext(castets[i], casface);
tsspivot(casface, checkseg);
if (checkseg.sh != NULL) {
tssbond1(newface, checkseg);
+ sstbond(checkseg, newface);
}
enextself(newface); // [d,a], [d,b], [d,c]
enext2(castets[i], casface);
tsspivot(casface, checkseg);
if (checkseg.sh != NULL) {
tssbond1(newface, checkseg);
+ sstbond(checkseg, newface);
}
enextself(*splittet);
}
@@ -769,6 +773,7 @@ void tetgenmesh::flip26(point newpt, triface* splitface, int flipflag)
tsspivot(casface, checkseg);
if (checkseg.sh != NULL) {
tssbond1(newface, checkseg);
+ sstbond(checkseg, newface);
}
enextself(casface);
enextself(newface);
@@ -786,6 +791,7 @@ void tetgenmesh::flip26(point newpt, triface* splitface, int flipflag)
tsspivot(casface, checkseg);
if (checkseg.sh != NULL) {
tssbond1(newface, checkseg);
+ sstbond(checkseg, newface);
}
enextself(casface);
enextself(newface);
@@ -993,6 +999,7 @@ void tetgenmesh::flipn2n(point newpt, triface* splitedge, int flipflag)
tsspivot(casface, checkseg);
if (checkseg.sh != NULL) {
tssbond1(newface, checkseg);
+ sstbond(checkseg, newface);
}
enextself(casface);
enextself(newface);
@@ -1005,6 +1012,7 @@ void tetgenmesh::flipn2n(point newpt, triface* splitedge, int flipflag)
tsspivot(casface, checkseg);
if (checkseg.sh != NULL) {
tssbond1(newface, checkseg);
+ sstbond(checkseg, newface);
}
enextself(casface);
enextself(newface);
@@ -1257,6 +1265,7 @@ void tetgenmesh::flip23(triface* fliptets, int hullflag, int flipflag)
enextfnext(fliptets[i], newface);
enextself(newface);
tssbond1(newface, checkseg);
+ sstbond(checkseg, newface);
}
}
// The top three: da, db, dc. Each edge belongs to two tets.
@@ -1266,9 +1275,11 @@ void tetgenmesh::flip23(triface* fliptets, int hullflag, int flipflag)
if (checkseg.sh != NULL) {
enext(fliptets[i], newface);
tssbond1(newface, checkseg);
+ sstbond(checkseg, newface);
enext0fnext(fliptets[(i + 2) % 3], newface);
enextself(newface);
tssbond1(newface, checkseg);
+ sstbond(checkseg, newface);
}
}
// The bot three: ae, be, ce. Each edge belongs to two tets.
@@ -1278,9 +1289,11 @@ void tetgenmesh::flip23(triface* fliptets, int hullflag, int flipflag)
if (checkseg.sh != NULL) {
enext2(fliptets[i], newface);
tssbond1(newface, checkseg);
+ sstbond(checkseg, newface);
enext0fnext(fliptets[(i + 2) % 3], newface);
enext2self(newface);
tssbond1(newface, checkseg);
+ sstbond(checkseg, newface);
}
}
}
@@ -1455,15 +1468,19 @@ void tetgenmesh::flip32(triface* fliptets, int hullflag, int flipflag)
if (checksubsegs) {
// Check if the flip edge is subsegment.
tsspivot(fliptets[0], checkseg);
- if ((checkseg.sh != NULL) && !sinfected(checkseg)) {
- // This subsegment will be flipped. Queue it.
- if (b->verbose > 1) {
- printf(" Queue a flipped segment (%d, %d).\n",
- pointmark(sorg(checkseg)), pointmark(sdest(checkseg)));
+ if ((checkseg.sh != NULL)) {
+ if (!sinfected(checkseg)) {
+ // This subsegment will be flipped. Queue it.
+ if (b->verbose > 1) {
+ printf(" Queue a flipped segment (%d, %d).\n",
+ pointmark(sorg(checkseg)), pointmark(sdest(checkseg)));
+ }
+ sinfect(checkseg); // Only save it once.
+ subsegstack->newindex((void **) &pssub);
+ *pssub = checkseg;
}
- sinfect(checkseg); // Only save it once.
- subsegstack->newindex((void **) &pssub);
- *pssub = checkseg;
+ // Clean the seg-to-tet pointer.
+ stdissolve(checkseg);
}
}
@@ -1566,6 +1583,7 @@ void tetgenmesh::flip32(triface* fliptets, int hullflag, int flipflag)
tsspivot(topcastets[i], checkseg);
if (checkseg.sh != NULL) {
tssbond1(fliptets[0], checkseg);
+ sstbond(checkseg, fliptets[0]);
}
enextself(fliptets[0]);
}
@@ -1577,6 +1595,7 @@ void tetgenmesh::flip32(triface* fliptets, int hullflag, int flipflag)
tsspivot(casface, checkseg);
if (checkseg.sh != NULL) {
tssbond1(newface, checkseg);
+ sstbond(checkseg, newface);
}
enextself(fliptets[0]);
}
@@ -1585,6 +1604,7 @@ void tetgenmesh::flip32(triface* fliptets, int hullflag, int flipflag)
tsspivot(botcastets[i], checkseg);
if (checkseg.sh != NULL) {
tssbond1(fliptets[1], checkseg);
+ sstbond(checkseg, fliptets[1]);
}
enext2self(fliptets[1]);
}
@@ -1596,6 +1616,7 @@ void tetgenmesh::flip32(triface* fliptets, int hullflag, int flipflag)
tsspivot(casface, checkseg);
if (checkseg.sh != NULL) {
tssbond1(newface, checkseg);
+ sstbond(checkseg, newface);
}
enext2self(fliptets[1]);
}
@@ -1873,7 +1894,7 @@ void tetgenmesh::lawsonflip3d(int flipflag)
point *pt, pd, pe;
REAL sign, ori;
long flipcount;
- // bool bflag; // for flipflag = 3.
+ bool bflag; // for flipflag = 3.
int n, i;
int *iptr;
@@ -1963,29 +1984,29 @@ void tetgenmesh::lawsonflip3d(int flipflag)
}
if (i == 3) {
// A 2-to-3 flip is found.
- // if (flipflag == 3) {
- // // Do not flip a subface.
- // tspivot(fliptet, checksh);
- // bflag = (checksh.sh == NULL);
- // } else {
- // bflag = true;
- // }
- // if (bflag) {
+ if (flipflag == 3) {
+ // Do not flip a subface.
+ tspivot(fliptet, checksh);
+ bflag = (checksh.sh == NULL);
+ } else {
+ bflag = true;
+ }
+ if (bflag) {
fliptets[0] = fliptet; // tet abcd, d is the new vertex.
symedge(fliptets[0], fliptets[1]); // tet bace.
flip23(fliptets, 0, flipflag);
recenttet = fliptets[0]; // for point location.
- // }
+ }
} else {
// A 3-to-2 or 4-to-4 may possible.
- // if (flipflag == 3) {
- // // Do not flip a subsegment.
- // tsspivot(fliptet, checkseg);
- // bflag = (checkseg.sh == NULL);
- // } else {
- // bflag = true;
- // }
- // if (bflag) {
+ if (flipflag == 3) {
+ // Do not flip a subsegment.
+ tsspivot(fliptet, checkseg);
+ bflag = (checkseg.sh == NULL);
+ } else {
+ bflag = true;
+ }
+ if (bflag) {
enext0fnext(fliptet, fliptets[0]);
esymself(fliptets[0]); // tet badc, d is the new vertex.
n = 0;
@@ -2015,19 +2036,29 @@ void tetgenmesh::lawsonflip3d(int flipflag)
recenttet = fliptets[0]; // for point location.
} else {
// An unflipable face. Will be flipped later.
- if (flipflag > 1) {
+ if (flipflag == 2) { // if (flipflag > 1) {
// Queue all other faces at the edge for flipping.
+ if (b->verbose > 1) {
+ printf(" Queue an N32 edge (%d, %d)",
+ pointmark(org(fliptets[0])), pointmark(dest(fliptets[0])));
+ }
pe = apex(fliptets[0]);
fliptets[1] = fliptets[0];
while (1) {
+ if (b->verbose > 1) {
+ printf(" %d", pointmark(apex(fliptets[1])));
+ }
pd = oppo(fliptets[1]);
futureflip = flippush(futureflip, &fliptets[1], pd);
fnextself(fliptets[1]);
if (apex(fliptets[1]) == pe) break;
}
+ if (b->verbose > 1) {
+ printf("\n");
+ }
}
}
- // } // bflag
+ } // bflag
} // if (i == 3)
} // if (sign < 0)
}
diff --git a/contrib/Tetgen/main.cxx b/contrib/Tetgen/main.cxx
index dafc3a88779712980127b355ae6e461a29329e9b..82edf9afd5dbe80b560e2f0d3c051697d15835ab 100644
--- a/contrib/Tetgen/main.cxx
+++ b/contrib/Tetgen/main.cxx
@@ -187,20 +187,32 @@ void tetrahedralize(tetgenbehavior *b, tetgenio *in, tetgenio *out,
printf("seconds: %g\n", (tv[4] - tv[3]) / (REAL) CLOCKS_PER_SEC);
}
}
+
+ if (b->quality) {
+ m.enforcequality();
+ }
+
+ tv[5] = clock();
+
+ if (!b->quiet) {
+ if (b->quality) {
+ printf("Quality seconds: %g\n", (tv[5] - tv[4])/(REAL) CLOCKS_PER_SEC);
+ }
+ }
- if (b->plc) {
+ if (b->plc) {
if (b->convexity == 0) { // if has no -c option.
m.carveholes();
}
}
- tv[5] = clock();
+ tv[6] = clock();
if (!b->quiet) {
if (b->plc) {
if (b->convexity == 0) { // if has no -c option.
printf("Holes and region seconds: %g\n",
- (tv[5] - tv[4]) / (REAL) CLOCKS_PER_SEC);
+ (tv[6] - tv[5]) / (REAL) CLOCKS_PER_SEC);
}
}
}
@@ -284,7 +296,8 @@ void tetrahedralize(tetgenbehavior *b, tetgenio *in, tetgenio *out,
if (b->docheck) {
if (b->plc) {
- m.checkshells(1);
+ if (m.checkshells(1) > 0) assert(0);
+ if (m.checksegments() > 0) assert(0);
}
if (m.checkdelaunay(b->plc) > 0) assert(0);
}
@@ -348,7 +361,7 @@ void tetrahedralize(char *switches, tetgenio *in, tetgenio *out,
}
// FOR DEBUG -S1 option.
- if (b.steiner > 0) {
+ if (b.steinerleft > 0) {
test_tri_tri(&b, &in);
terminatetetgen(0);
}
diff --git a/contrib/Tetgen/meshstat.cxx b/contrib/Tetgen/meshstat.cxx
index 8c4e7ec30492c9bd11a3358b269b6ba20013d6da..c5c0c992c1ee91c73a2b9e6bac8a8c0288ddb249 100644
--- a/contrib/Tetgen/meshstat.cxx
+++ b/contrib/Tetgen/meshstat.cxx
@@ -104,13 +104,13 @@ int tetgenmesh::mi1mo3[3] = {2, 0, 1};
// //
///////////////////////////////////////////////////////////////////////////////
-void tetgenmesh::checkmesh()
+int tetgenmesh::checkmesh()
{
triface tetloop;
triface neightet, symtet;
point pa, pb, pc, pd;
REAL ori;
- int horrors;
+ int horrors, i;
if (!b->quiet) {
printf(" Checking consistency of mesh...\n");
@@ -138,11 +138,13 @@ void tetgenmesh::checkmesh()
horrors++;
}
}
- if (infected(tetloop)) {
+ if (infected(tetloop)) {
+ // This may be a bug. Report it.
printf(" !! (%d, %d, %d, %d) is infected.\n", pointmark(pa),
pointmark(pb), pointmark(pc), pointmark(pd));
}
if (marktested(tetloop)) {
+ // This may be a bug. Report it.
printf(" !! (%d, %d, %d, %d) is marked.\n", pointmark(pa),
pointmark(pb), pointmark(pc), pointmark(pd));
}
@@ -178,6 +180,16 @@ void tetgenmesh::checkmesh()
pointmark(oppo(neightet)));
horrors++;
}
+ // Check if they have the same apex.
+ if (apex(neightet) != pc) {
+ printf(" !! !! Wrong face-face bond:\n");
+ printf(" First: (%d, %d, %d, %d)\n", pointmark(pa),
+ pointmark(pb), pointmark(pc), pointmark(pd));
+ printf(" Second: (%d, %d, %d, %d)\n", pointmark(org(neightet)),
+ pointmark(dest(neightet)), pointmark(apex(neightet)),
+ pointmark(oppo(neightet)));
+ horrors++;
+ }
// Check if they have the same opposite.
if (oppo(neightet) == pd) {
printf(" !! !! Two identical tetra:\n");
@@ -189,6 +201,22 @@ void tetgenmesh::checkmesh()
horrors++;
}
}
+ if (facemarked(tetloop)) {
+ // This may be a bug. Report it.
+ printf(" !! tetface (%d, %d, %d) %d is marked.\n", pointmark(pa),
+ pointmark(pb), pointmark(pc), pointmark(pd));
+ }
+ }
+ // Check the six edges of this tet.
+ for (i = 0; i < 6; i++) {
+ tetloop.loc = edge2locver[i][0];
+ tetloop.ver = edge2locver[i][1];
+ if (edgemarked(tetloop)) {
+ // This may be a bug. Report it.
+ printf(" !! tetedge (%d, %d) %d, %d is marked.\n",
+ pointmark(org(tetloop)), pointmark(dest(tetloop)),
+ pointmark(apex(tetloop)), pointmark(oppo(tetloop)));
+ }
}
tetloop.tet = tetrahedrontraverse();
}
@@ -200,6 +228,8 @@ void tetgenmesh::checkmesh()
printf(" !! !! !! !! %d %s witnessed.\n", horrors,
horrors > 1 ? "abnormity" : "abnormities");
}
+
+ return horrors;
}
///////////////////////////////////////////////////////////////////////////////
@@ -215,7 +245,8 @@ int tetgenmesh::checkshells(int sub2tet)
triface neightet, symtet;
face shloop, spinsh, nextsh;
face checkseg;
- point pa, pb;
+ point pa, pb, *ppt;
+ int count;
int horrors, i;
tetrahedron ptr;
@@ -273,9 +304,33 @@ int tetgenmesh::checkshells(int sub2tet)
// Check the tet-subface connections.
stpivot(shloop, neightet);
if (neightet.tet != NULL) {
+ // Check if the tet and subface have the same vertices.
+ ppt = (point *) shloop.sh;
+ for (i = 0; i < 3; i++) {
+ pinfect(ppt[3 + i]); // Infect the subface vertices.
+ }
+ ppt = (point *) neightet.tet;
+ count = 0; // Count the infected vertices of this tet.
+ for (i = 0; i < 4; i++) {
+ if (pinfected(ppt[4 + i])) count++;
+ }
+ ppt = (point *) shloop.sh;
+ for (i = 0; i < 3; i++) {
+ puninfect(ppt[3 + i]); // Uninfect the subface vertices.
+ }
+ if (count != 3) {
+ printf(" !! !! Wrong tet-sub connection (face does not match).\n");
+ printf(" Sub: x%lx (%d, %d, %d).\n", (unsigned long) shloop.sh,
+ pmark(sorg(shloop)), pmark(sdest(shloop)), pmark(sapex(shloop)));
+ printf(" Tet: x%lx (%d, %d, %d, %d).\n",
+ (unsigned long) neightet.tet, pmark(org(neightet)),
+ pmark(dest(neightet)), pmark(apex(neightet)),
+ pmark(oppo(neightet)));
+ horrors++;
+ }
tspivot(neightet, spinsh);
if (spinsh.sh != shloop.sh) {
- printf(" !! !! Wrong connection betwee tet and subface.\n");
+ printf(" !! !! Wrong tet-sub connection (wrong subfaces).\n");
printf(" Sub: x%lx (%d, %d, %d).\n", (unsigned long) shloop.sh,
pmark(sorg(shloop)), pmark(sdest(shloop)), pmark(sapex(shloop)));
printf(" Tet: x%lx (%d, %d, %d, %d).\n",
@@ -287,7 +342,7 @@ int tetgenmesh::checkshells(int sub2tet)
symself(neightet);
tspivot(neightet, spinsh);
if (spinsh.sh != shloop.sh) {
- printf(" !! !! Wrong connection betwee tet and subface.\n");
+ printf(" !! !! Wrong tet-sub connection (wrong subfaces).\n");
printf(" Sub: x%lx (%d, %d, %d).\n", (unsigned long) shloop.sh,
pmark(sorg(shloop)), pmark(sdest(shloop)), pmark(sapex(shloop)));
printf(" Tet: x%lx (%d, %d, %d, %d).\n",
@@ -304,6 +359,16 @@ int tetgenmesh::checkshells(int sub2tet)
// horrors++;
}
}
+ if (sinfected(shloop)) {
+ // This may be a bug. report it.
+ printf(" !! A infected subface: (%d, %d, %d).\n", pmark(sorg(shloop)),
+ pmark(sdest(shloop)), pmark(sapex(shloop)));
+ }
+ if (smarktested(shloop)) {
+ // This may be a bug. report it.
+ printf(" !! A marked subface: (%d, %d, %d).\n", pmark(sorg(shloop)),
+ pmark(sdest(shloop)), pmark(sapex(shloop)));
+ }
shloop.sh = shellfacetraverse(subfacepool);
}
@@ -440,11 +505,11 @@ int tetgenmesh::checksegments()
do {
tsspivot(neightet, checkseg);
if (checkseg.sh != sseg.sh) {
- printf(" !! Wrong tet-seg connection.\n");
+ printf(" !! Wrong tet->seg connection.\n");
printf(" Tet: x%lx (%d, %d, %d, %d) - ",
- (unsigned long) tetloop.tet, pointmark(org(tetloop)),
- pointmark(dest(tetloop)), pointmark(apex(tetloop)),
- pointmark(oppo(tetloop)));
+ (unsigned long) neightet.tet, pointmark(org(neightet)),
+ pointmark(dest(neightet)), pointmark(apex(neightet)),
+ pointmark(oppo(neightet)));
if (checkseg.sh != NULL) {
printf("Seg x%lx (%d, %d).\n", (unsigned long) checkseg.sh,
pointmark(sorg(checkseg)), pointmark(sdest(checkseg)));
@@ -456,6 +521,23 @@ int tetgenmesh::checksegments()
fnextself(neightet);
} while (neightet.tet != tetloop.tet);
}
+ // Check the seg->tet pointer.
+ stpivot(sseg, neightet);
+ if (neightet.tet == NULL) {
+ printf(" !! Wrong seg->tet connection (A NULL tet).\n");
+ horrors++;
+ } else {
+ if (!(((org(neightet) == pa) && (dest(neightet) == pb)) ||
+ ((org(neightet) == pb) && (dest(neightet) == pa)))) {
+ printf(" !! Wrong seg->tet connection (Wrong edge).\n");
+ printf(" Tet: x%lx (%d, %d, %d, %d) - Seg: x%lx (%d, %d).\n",
+ (unsigned long) neightet.tet, pointmark(org(neightet)),
+ pointmark(dest(neightet)), pointmark(apex(neightet)),
+ pointmark(oppo(neightet)), (unsigned long) sseg.sh,
+ pointmark(pa), pointmark(pb));
+ horrors++;
+ }
+ }
}
}
}
@@ -471,6 +553,42 @@ int tetgenmesh::checksegments()
return horrors;
}
+///////////////////////////////////////////////////////////////////////////////
+// //
+// checkconforming() Check if the mesh is boundary conforming Delaunay. //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+int tetgenmesh::checkconforming()
+{
+ face shloop;
+ int horrors;
+
+ horrors = 0;
+
+ // Find all encroached segments.
+ subsegpool->traversalinit();
+ shloop.shver = 0;
+ shloop.sh = shellfacetraverse(subsegpool);
+ while (shloop.sh != NULL) {
+ if (checkedge4encroach(shloop, NULL, 0)) {
+ printf(" !! Segment x%lx (%d, %d) is encroached.\n",
+ (unsigned long) shloop.sh, pointmark(sorg(shloop)),
+ pointmark(sdest(shloop)));
+ horrors++;
+ }
+ shloop.sh = shellfacetraverse(subsegpool);
+ }
+
+ if (horrors == 0) {
+ printf(" Segments are boundary conforming Delaunay.\n");
+ } else {
+ printf(" !! !! !! !! Found %d encroached segments.\n", horrors);
+ }
+
+ return horrors;
+}
+
///////////////////////////////////////////////////////////////////////////////
// //
// algorithmstatistics() Report algorithmic performances. //
@@ -592,8 +710,8 @@ void tetgenmesh::statistics()
if (b->verbose > 0) {
printf(" Euler characteristic of mesh domain: %ld\n", pointpool->items
- meshedges + facenumber - tetnumber);
- printf(" Euler characteristic of boundary: %ld\n", pointpool->items
- - meshsubedges + subfacepool->items);
+ //printf(" Euler characteristic of boundary: %ld\n", pointpool->items
+ // - meshsubedges + subfacepool->items);
printf("\n");
}
@@ -871,13 +989,16 @@ void tetgenmesh::psh(face *s)
} else {
printf(" [8] = x%lx %d\n", (unsigned long) prtsh.sh, prtsh.shver);
}
+ }
- decode(s->sh[9], prttet);
- if (prttet.tet == NULL) {
+ // Print the adjacent tet of this subface or segment.
+ decode(s->sh[9], prttet);
+ if (prttet.tet == NULL) {
printf(" [9] = Outer space\n");
- } else {
- printf(" [9] = x%lx %d\n",(unsigned long) prttet.tet, prttet.loc);
- }
+ } else {
+ printf(" [9] = x%lx (%d, %d, %d, %d)\n",(unsigned long) prttet.tet,
+ pointmark(org(prttet)), pointmark(dest(prttet)),
+ pointmark(apex(prttet)), pointmark(oppo(prttet)));
}
}
@@ -1129,9 +1250,11 @@ void tetgenmesh::psubface(int i, int j, int k)
///////////////////////////////////////////////////////////////////////////////
// Print the information of the subsegment (i, j).
+// Return 1 if seg-to-tet pointer is broken.
-void tetgenmesh::psubseg(int i, int j)
+int tetgenmesh::psubseg(int i, int j)
{
+ triface t;
face s; //, s1;
point forg, fdest;
bool bflag;
@@ -1157,6 +1280,26 @@ void tetgenmesh::psubseg(int i, int j)
fdest = farsdest(s);
printf(" seg x%lx (%d, %d) < (%d, %d)\n", (unsigned long) s.sh, i, j,
pointmark(forg), pointmark(fdest));
+ // Chekc seg-to-tet pointer
+ stpivot(s, t);
+ if (t.tet != NULL) {
+ if (t.tet[4] != NULL) {
+ forg = org(t);
+ fdest = dest(t);
+ if (((pointmark(forg) == i) && (pointmark(fdest) == j)) ||
+ ((pointmark(forg) == j) && (pointmark(fdest) == i))) {
+ printf(" adj tet x%lx (%d, %d, %d, %d)\n", (unsigned long) t.tet,
+ pointmark(forg), pointmark(fdest), pointmark(apex(t)),
+ pointmark(oppo(t)));
+ } else {
+ printf(" !! Wrong seg-to-tet pointer.\n");
+ return 1;
+ }
+ } else {
+ printf(" !! A DEAD seg-to-tet pointer.\n");
+ return 1;
+ }
+ }
/*// Print the adjacent subsegments at i and j.
senext2(s, s1);
spivotself(s1);
@@ -1184,6 +1327,7 @@ void tetgenmesh::psubseg(int i, int j)
if (!bflag) {
printf(" !! Not exist.\n");
}
+ return 0;
}
///////////////////////////////////////////////////////////////////////////////
@@ -1571,4 +1715,50 @@ void tetgenmesh::dump_facetof(face *pssub, char *filename)
delete tmpfaces;
}
+///////////////////////////////////////////////////////////////////////////////
+// Dump the new tets of a cavity (in draw_tet() lua commands)
+
+void tetgenmesh::dump_cavitynewtets()
+{
+ arraypool *newtets;
+ triface searchtet, neightet, *parytet;
+ point *ppt;
+ int i;
+
+ newtets = new arraypool(sizeof(triface), 8);
+
+ // Collect all tets of the DT.
+ marktest(recenttet);
+ newtets->newindex((void **) &parytet);
+ *parytet = recenttet;
+ for (i = 0; i < newtets->objects; i++) {
+ searchtet = * (triface *) fastlookup(newtets, i);
+ for (searchtet.loc = 0; searchtet.loc < 4; searchtet.loc++) {
+ sym(searchtet, neightet);
+ if (!marktested(neightet)) {
+ marktest(neightet);
+ newtets->newindex((void **) &parytet);
+ *parytet = neightet;
+ }
+ }
+ }
+ // Comment: All new tets are marktested.
+
+ // Output the new tets.
+ for (i = 0; i < newtets->objects; i++) {
+ parytet = (triface *) fastlookup(newtets, i);
+ ppt = (point *) parytet->tet;
+ if (ppt[7] != dummypoint) {
+ printf("p:draw_tet(%d, %d, %d, %d) -- %i\n", pointmark(ppt[4]),
+ pointmark(ppt[5]), pointmark(ppt[6]), pointmark(ppt[7]), i);
+ } else {
+ printf("-- p:draw_tet(%d, %d, %d, -1) -- %i\n", pointmark(ppt[4]),
+ pointmark(ppt[5]), pointmark(ppt[6]), i);
+ }
+ unmarktest(*parytet); // Unmarktest it.
+ }
+
+ delete newtets;
+}
+
#endif // #ifndef meshstatCXX
\ No newline at end of file
diff --git a/contrib/Tetgen/surface.cxx b/contrib/Tetgen/surface.cxx
index b0305591abb2247454d9e24ef258673532bb8f1e..e9a72b232d3f59cb1a2872968275a0cf63d7125d 100644
--- a/contrib/Tetgen/surface.cxx
+++ b/contrib/Tetgen/surface.cxx
@@ -260,7 +260,7 @@ enum tetgenmesh::location tetgenmesh::sinsertvertex(point insertpt,
face neighsh, newsh, casout, casin;
face aseg, bseg, aoutseg, boutseg;
face checkseg;
- triface neightet;
+ triface neightet, spintet;
point pa, pb, pc;
enum location loc;
REAL sign, ori, area;
@@ -598,6 +598,22 @@ enum tetgenmesh::location tetgenmesh::sinsertvertex(point insertpt,
printf(" Split seg (%d, %d) by %d.\n", pointmark(pa), pointmark(pb),
pointmark(insertpt));
}
+ // Detach the adjacent tets from this segment.
+ stpivot(aseg, neightet);
+ if (neightet.tet != NULL) {
+ // It should not be a dead tet.
+ assert(neightet.tet[4] != NULL); // SELF_CHECK
+ assert(((org(neightet) == pa) && (dest(neightet) == pb)) ||
+ ((org(neightet) == pb) && (dest(neightet) == pa))); // SELF_CHECK
+ spintet = neightet;
+ while (1) {
+ tssdissolve(spintet);
+ fnextself(spintet);
+ if (spintet.tet == neightet.tet) break;
+ }
+ // Clean the seg-to-tet pointer.
+ stdissolve(aseg);
+ }
// Insert the new point p.
makeshellface(subsegpool, &bseg);
setshvertices(bseg, insertpt, pb, NULL);
diff --git a/contrib/Tetgen/tetgen.h b/contrib/Tetgen/tetgen.h
index b44a68876f4033b6c4503dc76c28ee9bd8247e00..b7e971f039fc8b72e3d06494d3e3b60989664b53 100644
--- a/contrib/Tetgen/tetgen.h
+++ b/contrib/Tetgen/tetgen.h
@@ -442,11 +442,11 @@ class tetgenbehavior { // Begin of class tetgenbehavior
int nomerge; // -M
int nobisect; // -Y
int nojettison; // -J
- int steiner; // -S with a number
int docheck; // -C
int quiet; // -Q
int verbose; // -V
int useshelles; // -p, -r, -q, -d, or -R
+ long steinerleft; // -S with a number
REAL minratio; // number after -q
REAL goodratio; // number calculated from 'minratio'
REAL minangle; // minimum angle bound
@@ -1373,9 +1373,13 @@ void tssbond1(triface& t, face& s) {
}
// Bond the segment.
((shellface *) (t).tet[8])[locver2edge[(t).loc][(t).ver]] = sencode((s));
-}
+}
+
+// sstbond() -- bond a tet edge to a segment (only at the segment side).
-// tssdissolve() -- dissolve a tet-seg bond at the tet edge.
+#define sstbond(s, t) (s).sh[9] = (shellface) encode(t)
+
+// tssdissolve() -- dissolve a tet-seg bond at the tet edge side.
void tssdissolve(triface& t) {
if ((t).tet[8] != NULL) {
@@ -1383,6 +1387,9 @@ void tssdissolve(triface& t) {
}
}
+// sstdissolve() -- dissolve a tet-seg bond at the segment side.
+// Use stdissolve().
+
///////////////////////////////////////////////////////////////////////////////
// //
// Primitives for points //
@@ -1473,6 +1480,7 @@ memorypool *subsegpool, *subfacepool;
memorypool *tet2segpool, *tet2subpool;
memorypool *pointpool;
memorypool *flippool;
+memorypool *badsegpool, *badsubpool, *badtetpool;
// A dummy point at infinity (see [Guibas & Stolfi 85]). All the hull
// tetrahedra having this point as a vertex.
@@ -1488,6 +1496,12 @@ arraypool *caveshlist, *caveshbdlist;
// Stacks used by the boundary recovery algorithm.
arraypool *subsegstack, *subfacstack;
+// Arrays used for facet recovery.
+arraypool *tg_crosstets, *tg_topnewtets, *tg_botnewtets;
+arraypool *tg_topfaces, *tg_botfaces, *tg_midfaces;
+arraypool *tg_topshells, *tg_botshells, *tg_facfaces;
+arraypool *tg_toppoints, *tg_botpoints, *tg_facpoints;
+
// Variables for accessing data fields (initialized in initializepools()).
int point2tetindex, pointmarkindex;
int elemmarkerindex;
@@ -1529,7 +1543,7 @@ long triedgcount, triedgcopcount, trivercopcount;
long across_face_count, across_edge_count, across_max_count;
long r1count, r2count, r3count;
long maxcavsize, maxregionsize;
-long ndelaunayedgecount, cavityexpcount;
+long cavitycount, ndelaunayedgecount, cavityexpcount;
///////////////////////////////////////////////////////////////////////////////
// //
@@ -1681,6 +1695,8 @@ void meshsurface();
enum intersection finddirection(triface* searchtet, point endpt);
enum intersection scoutsegment(face* sseg, triface* searchtet, point* refpt);
void getsegmentsplitpoint(face* sseg, point refpt, REAL* vt);
+void getsegmentsplitpoint2(face* sseg, point refpt, REAL* vt);
+void getsegmentsplitpoint3(face* sseg, point refpt, REAL* vt);
void delaunizesegments();
enum intersection scoutsubface(face* ssub, triface* searchtet);
@@ -1696,12 +1712,27 @@ void restorecavity(arraypool*, arraypool*, arraypool*);
void splitsubedge(point, face*, arraypool*, arraypool*);
void constrainedfacets();
-enum intersection scoutedge(point, point, triface*, arraypool*, arraypool*,
- arraypool*);
+enum intersection scoutsegment2(face*, triface*, arraypool*);
+bool tetrasegcavity(face*, arraypool*, arraypool*, arraypool*, arraypool*,
+ arraypool*, arraypool*, arraypool*, arraypool*);
+bool recoversegments(arraypool*, arraypool*, arraypool*, arraypool*,
+ arraypool*, arraypool*);
+void constrainedsegments();
void formskeleton();
void carveholes();
+///////////////////////////////////////////////////////////////////////////////
+// //
+// Mesh refinement functions. //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+bool checkedge4encroach(face& seg, point testpt, int enqflag);
+void repairencsegs();
+
+void enforcequality();
+
///////////////////////////////////////////////////////////////////////////////
// //
// Mesh input & output functions. //
@@ -1730,11 +1761,11 @@ void outvoronoi(tetgenio* out);
// //
///////////////////////////////////////////////////////////////////////////////
-void checkmesh();
+int checkmesh();
int checkshells(int);
int checkdelaunay(int);
int checksegments();
-void checkconforming();
+int checkconforming();
void algorithmstatistics();
void qualitystatistics();
void statistics();
@@ -1754,11 +1785,16 @@ void initialize()
tet2segpool = tet2subpool = (memorypool *) NULL;
pointpool = (memorypool *) NULL;
flippool = (memorypool *) NULL;
+ badsegpool = badsubpool = badtetpool = (memorypool *) NULL;
dummypoint = (point) NULL;
futureflip = (badface *) NULL;
cavetetlist = cavebdrylist = caveoldtetlist = (arraypool *) NULL;
caveshlist = caveshbdlist = (arraypool *) NULL;
subsegstack = subfacstack = (arraypool *) NULL;
+ arraypool *tg_crosstets, *tg_topnewtets, *tg_botnewtets;
+ tg_topfaces = tg_botfaces = tg_midfaces = NULL;
+ tg_topshells = tg_botshells = tg_facfaces = NULL;
+ tg_toppoints = tg_botpoints = tg_facpoints = NULL;
point2tetindex = pointmarkindex = 0;
elemmarkerindex = 0;
elemattribindex = volumeboundindex = highorderindex = 0;
@@ -1781,7 +1817,7 @@ void initialize()
across_face_count = across_edge_count = across_max_count = 0l;
r1count = r2count = r3count = 0l;
maxcavsize = maxregionsize = 0l;
- ndelaunayedgecount = cavityexpcount = 0l;
+ cavitycount = ndelaunayedgecount = cavityexpcount = 0l;
}
void deinitialize()
@@ -1834,7 +1870,7 @@ void pteti(int i, int j, int k, int l);
void pface(int i, int j, int k);
bool pedge(int i, int j);
void psubface(int i, int j, int k);
-void psubseg(int i, int j);
+int psubseg(int i, int j);
int pmark(point p);
void pvert(point p);
int pverti(int i);
@@ -1848,6 +1884,7 @@ void print_facearray(arraypool* facearray);
void print_subfacearray(arraypool* subfacearray);
void dump_cavity(arraypool *topfaces, arraypool *botfaces);
void dump_facetof(face* pssub, char* filename);
+void dump_cavitynewtets();
///////////////////////////////////////////////////////////////////////////////
}; // End of class tetgenmesh;