3D图形编程

开发平台：

C/C++

triangle.c：源码内容

/* A live triangle. The vertex survives. */
killorg = 0;
}
/* Walk clockwise about the vertex. */
oprevself(neighbor);
}
}
if (killorg) {
if (b->verbose > 1) {
printf(" Deleting vertex (%.12g, %.12g)n",
testvertex[0], testvertex[1]);
}
setvertextype(testvertex, UNDEADVERTEX);
m->undeads++;
}
}
}
/* Record changes in the number of boundary edges, and disconnect */
/* dead triangles from their neighbors. */
for (testtri.orient = 0; testtri.orient < 3; testtri.orient++) {
sym(testtri, neighbor);
if (neighbor.tri == m->dummytri) {
/* There is no neighboring triangle on this edge, so this edge */
/* is a boundary edge. This triangle is being deleted, so this */
/* boundary edge is deleted. */
m->hullsize--;
} else {
/* Disconnect the triangle from its neighbor. */
dissolve(neighbor);
/* There is a neighboring triangle on this edge, so this edge */
/* becomes a boundary edge when this triangle is deleted. */
m->hullsize++;
}
}
/* Return the dead triangle to the pool of triangles. */
triangledealloc(m, testtri.tri);
virusloop = (triangle **) traverse(&m->viri);
}
/* Empty the virus pool. */
poolrestart(&m->viri);
}
/*****************************************************************************/
/* */
/* regionplague() Spread regional attributes and/or area constraints */
/* (from a .poly file) throughout the mesh. */
/* */
/* This procedure operates in two phases. The first phase spreads an */
/* attribute and/or an area constraint through a (segment-bounded) region. */
/* The triangles are marked to ensure that each triangle is added to the */
/* virus pool only once, so the procedure will terminate. */
/* */
/* The second phase uninfects all infected triangles, returning them to */
/* normal. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void regionplague(struct mesh *m, struct behavior *b,
REAL attribute, REAL area)
#else /* not ANSI_DECLARATORS */
void regionplague(m, b, attribute, area)
struct mesh *m;
struct behavior *b;
REAL attribute;
REAL area;
#endif /* not ANSI_DECLARATORS */
{
struct otri testtri;
struct otri neighbor;
triangle **virusloop;
triangle **regiontri;
struct osub neighborsubseg;
vertex regionorg, regiondest, regionapex;
triangle ptr; /* Temporary variable used by sym() and onext(). */
subseg sptr; /* Temporary variable used by tspivot(). */
if (b->verbose > 1) {
printf(" Marking neighbors of marked triangles.n");
}
/* Loop through all the infected triangles, spreading the attribute */
/* and/or area constraint to their neighbors, then to their neighbors' */
/* neighbors. */
traversalinit(&m->viri);
virusloop = (triangle **) traverse(&m->viri);
while (virusloop != (triangle **) NULL) {
testtri.tri = *virusloop;
/* A triangle is marked as infected by messing with one of its pointers */
/* to subsegments, setting it to an illegal value. Hence, we have to */
/* temporarily uninfect this triangle so that we can examine its */
/* adjacent subsegments. */
uninfect(testtri);
if (b->regionattrib) {
/* Set an attribute. */
setelemattribute(testtri, m->eextras, attribute);
}
if (b->vararea) {
/* Set an area constraint. */
setareabound(testtri, area);
}
if (b->verbose > 2) {
/* Assign the triangle an orientation for convenience in */
/* checking its vertices. */
testtri.orient = 0;
org(testtri, regionorg);
dest(testtri, regiondest);
apex(testtri, regionapex);
printf(" Checking (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)n",
regionorg[0], regionorg[1], regiondest[0], regiondest[1],
regionapex[0], regionapex[1]);
}
/* Check each of the triangle's three neighbors. */
for (testtri.orient = 0; testtri.orient < 3; testtri.orient++) {
/* Find the neighbor. */
sym(testtri, neighbor);
/* Check for a subsegment between the triangle and its neighbor. */
tspivot(testtri, neighborsubseg);
/* Make sure the neighbor exists, is not already infected, and */
/* isn't protected by a subsegment. */
if ((neighbor.tri != m->dummytri) && !infected(neighbor)
&& (neighborsubseg.ss == m->dummysub)) {
if (b->verbose > 2) {
org(neighbor, regionorg);
dest(neighbor, regiondest);
apex(neighbor, regionapex);
printf(" Marking (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)n",
regionorg[0], regionorg[1], regiondest[0], regiondest[1],
regionapex[0], regionapex[1]);
}
/* Infect the neighbor. */
infect(neighbor);
/* Ensure that the neighbor's neighbors will be infected. */
regiontri = (triangle **) poolalloc(&m->viri);
*regiontri = neighbor.tri;
}
}
/* Remark the triangle as infected, so it doesn't get added to the */
/* virus pool again. */
infect(testtri);
virusloop = (triangle **) traverse(&m->viri);
}
/* Uninfect all triangles. */
if (b->verbose > 1) {
printf(" Unmarking marked triangles.n");
}
traversalinit(&m->viri);
virusloop = (triangle **) traverse(&m->viri);
while (virusloop != (triangle **) NULL) {
testtri.tri = *virusloop;
uninfect(testtri);
virusloop = (triangle **) traverse(&m->viri);
}
/* Empty the virus pool. */
poolrestart(&m->viri);
}
/*****************************************************************************/
/* */
/* carveholes() Find the holes and infect them. Find the area */
/* constraints and infect them. Infect the convex hull. */
/* Spread the infection and kill triangles. Spread the */
/* area constraints. */
/* */
/* This routine mainly calls other routines to carry out all these */
/* functions. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void carveholes(struct mesh *m, struct behavior *b, REAL *holelist, int holes,
REAL *regionlist, int regions)
#else /* not ANSI_DECLARATORS */
void carveholes(m, b, holelist, holes, regionlist, regions)
struct mesh *m;
struct behavior *b;
REAL *holelist;
int holes;
REAL *regionlist;
int regions;
#endif /* not ANSI_DECLARATORS */
{
struct otri searchtri;
struct otri triangleloop;
struct otri *regiontris;
triangle **holetri;
triangle **regiontri;
vertex searchorg, searchdest;
enum locateresult intersect;
int i;
triangle ptr; /* Temporary variable used by sym(). */
if (!(b->quiet || (b->noholes && b->convex))) {
printf("Removing unwanted triangles.n");
if (b->verbose && (holes > 0)) {
printf(" Marking holes for elimination.n");
}
}
if (regions > 0) {
/* Allocate storage for the triangles in which region points fall. */
regiontris = (struct otri *) trimalloc(regions *
(int) sizeof(struct otri));
} else {
regiontris = (struct otri *) NULL;
}
if (((holes > 0) && !b->noholes) || !b->convex || (regions > 0)) {
/* Initialize a pool of viri to be used for holes, concavities, */
/* regional attributes, and/or regional area constraints. */
poolinit(&m->viri, sizeof(triangle *), VIRUSPERBLOCK, VIRUSPERBLOCK, 0);
}
if (!b->convex) {
/* Mark as infected any unprotected triangles on the boundary. */
/* This is one way by which concavities are created. */
infecthull(m, b);
}
if ((holes > 0) && !b->noholes) {
/* Infect each triangle in which a hole lies. */
for (i = 0; i < 2 * holes; i += 2) {
/* Ignore holes that aren't within the bounds of the mesh. */
if ((holelist[i] >= m->xmin) && (holelist[i] <= m->xmax)
&& (holelist[i + 1] >= m->ymin) && (holelist[i + 1] <= m->ymax)) {
/* Start searching from some triangle on the outer boundary. */
searchtri.tri = m->dummytri;
searchtri.orient = 0;
symself(searchtri);
/* Ensure that the hole is to the left of this boundary edge; */
/* otherwise, locate() will falsely report that the hole */
/* falls within the starting triangle. */
org(searchtri, searchorg);
dest(searchtri, searchdest);
if (counterclockwise(m, b, searchorg, searchdest, &holelist[i]) >
0.0) {
/* Find a triangle that contains the hole. */
intersect = locate(m, b, &holelist[i], &searchtri);
if ((intersect != OUTSIDE) && (!infected(searchtri))) {
/* Infect the triangle. This is done by marking the triangle */
/* as infected and including the triangle in the virus pool. */
infect(searchtri);
holetri = (triangle **) poolalloc(&m->viri);
*holetri = searchtri.tri;
}
}
}
}
}
/* Now, we have to find all the regions BEFORE we carve the holes, because */
/* locate() won't work when the triangulation is no longer convex. */
/* (Incidentally, this is the reason why regional attributes and area */
/* constraints can't be used when refining a preexisting mesh, which */
/* might not be convex; they can only be used with a freshly */
/* triangulated PSLG.) */
if (regions > 0) {
/* Find the starting triangle for each region. */
for (i = 0; i < regions; i++) {
regiontris[i].tri = m->dummytri;
/* Ignore region points that aren't within the bounds of the mesh. */
if ((regionlist[4 * i] >= m->xmin) && (regionlist[4 * i] <= m->xmax) &&
(regionlist[4 * i + 1] >= m->ymin) &&
(regionlist[4 * i + 1] <= m->ymax)) {
/* Start searching from some triangle on the outer boundary. */
searchtri.tri = m->dummytri;
searchtri.orient = 0;
symself(searchtri);
/* Ensure that the region point is to the left of this boundary */
/* edge; otherwise, locate() will falsely report that the */
/* region point falls within the starting triangle. */
org(searchtri, searchorg);
dest(searchtri, searchdest);
if (counterclockwise(m, b, searchorg, searchdest, &regionlist[4 * i]) >
0.0) {
/* Find a triangle that contains the region point. */
intersect = locate(m, b, &regionlist[4 * i], &searchtri);
if ((intersect != OUTSIDE) && (!infected(searchtri))) {
/* Record the triangle for processing after the */
/* holes have been carved. */
otricopy(searchtri, regiontris[i]);
}
}
}
}
}
if (m->viri.items > 0) {
/* Carve the holes and concavities. */
plague(m, b);
}
/* The virus pool should be empty now. */
if (regions > 0) {
if (!b->quiet) {
if (b->regionattrib) {
if (b->vararea) {
printf("Spreading regional attributes and area constraints.n");
} else {
printf("Spreading regional attributes.n");
}
} else {
printf("Spreading regional area constraints.n");
}
}
if (b->regionattrib && !b->refine) {
/* Assign every triangle a regional attribute of zero. */
traversalinit(&m->triangles);
triangleloop.orient = 0;
triangleloop.tri = triangletraverse(m);
while (triangleloop.tri != (triangle *) NULL) {
setelemattribute(triangleloop, m->eextras, 0.0);
triangleloop.tri = triangletraverse(m);
}
}
for (i = 0; i < regions; i++) {
if (regiontris[i].tri != m->dummytri) {
/* Make sure the triangle under consideration still exists. */
/* It may have been eaten by the virus. */
if (!deadtri(regiontris[i].tri)) {
/* Put one triangle in the virus pool. */
infect(regiontris[i]);
regiontri = (triangle **) poolalloc(&m->viri);
*regiontri = regiontris[i].tri;
/* Apply one region's attribute and/or area constraint. */
regionplague(m, b, regionlist[4 * i + 2], regionlist[4 * i + 3]);
/* The virus pool should be empty now. */
}
}
}
if (b->regionattrib && !b->refine) {
/* Note the fact that each triangle has an additional attribute. */
m->eextras++;
}
}
/* Free up memory. */
if (((holes > 0) && !b->noholes) || !b->convex || (regions > 0)) {
pooldeinit(&m->viri);
}
if (regions > 0) {
trifree((VOID *) regiontris);
}
}
/** **/
/** **/
/********* Carving out holes and concavities ends here *********/
/********* Mesh quality maintenance begins here *********/
/** **/
/** **/
/*****************************************************************************/
/* */
/* tallyencs() Traverse the entire list of subsegments, and check each */
/* to see if it is encroached. If so, add it to the list. */
/* */
/*****************************************************************************/
#ifndef CDT_ONLY
#ifdef ANSI_DECLARATORS
void tallyencs(struct mesh *m, struct behavior *b)
#else /* not ANSI_DECLARATORS */
void tallyencs(m, b)
struct mesh *m;
struct behavior *b;
#endif /* not ANSI_DECLARATORS */
{
struct osub subsegloop;
int dummy;
traversalinit(&m->subsegs);
subsegloop.ssorient = 0;
subsegloop.ss = subsegtraverse(m);
while (subsegloop.ss != (subseg *) NULL) {
/* If the segment is encroached, add it to the list. */
dummy = checkseg4encroach(m, b, &subsegloop);
subsegloop.ss = subsegtraverse(m);
}
}
#endif /* not CDT_ONLY */
/*****************************************************************************/
/* */
/* precisionerror() Print an error message for precision problems. */
/* */
/*****************************************************************************/
#ifndef CDT_ONLY
void precisionerror()
{
printf("Try increasing the area criterion and/or reducing the minimumn");
printf(" allowable angle so that tiny triangles are not created.n");
#ifdef SINGLE
printf("Alternatively, try recompiling me with double precisionn");
printf(" arithmetic (by removing "#define SINGLE" from then");
printf(" source file or "-DSINGLE" from the makefile).n");
#endif /* SINGLE */
}
#endif /* not CDT_ONLY */
/*****************************************************************************/
/* */
/* splitencsegs() Split all the encroached subsegments. */
/* */
/* Each encroached subsegment is repaired by splitting it - inserting a */
/* vertex at or near its midpoint. Newly inserted vertices may encroach */
/* upon other subsegments; these are also repaired. */
/* */
/* `triflaws' is a flag that specifies whether one should take note of new */
/* bad triangles that result from inserting vertices to repair encroached */
/* subsegments. */
/* */
/*****************************************************************************/
#ifndef CDT_ONLY
#ifdef ANSI_DECLARATORS
void splitencsegs(struct mesh *m, struct behavior *b, int triflaws)
#else /* not ANSI_DECLARATORS */
void splitencsegs(m, b, triflaws)
struct mesh *m;
struct behavior *b;
int triflaws;
#endif /* not ANSI_DECLARATORS */
{
struct otri enctri;
struct otri testtri;
struct osub testsh;
struct osub currentenc;
struct badsubseg *encloop;
vertex eorg, edest, eapex;
vertex newvertex;
enum insertvertexresult success;
REAL segmentlength, nearestpoweroftwo;
REAL split;
REAL multiplier, divisor;
int acuteorg, acuteorg2, acutedest, acutedest2;
int dummy;
int i;
triangle ptr; /* Temporary variable used by stpivot(). */
subseg sptr; /* Temporary variable used by snext(). */
/* Note that steinerleft == -1 if an unlimited number */
/* of Steiner points is allowed. */
while ((m->badsubsegs.items > 0) && (m->steinerleft != 0)) {
traversalinit(&m->badsubsegs);
encloop = badsubsegtraverse(m);
while ((encloop != (struct badsubseg *) NULL) && (m->steinerleft != 0)) {
sdecode(encloop->encsubseg, currentenc);
sorg(currentenc, eorg);
sdest(currentenc, edest);
/* Make sure that this segment is still the same segment it was */
/* when it was determined to be encroached. If the segment was */
/* enqueued multiple times (because several newly inserted */
/* vertices encroached it), it may have already been split. */
if (!deadsubseg(currentenc.ss) &&
(eorg == encloop->subsegorg) && (edest == encloop->subsegdest)) {
/* To decide where to split a segment, we need to know if the */
/* segment shares an endpoint with an adjacent segment. */
/* The concern is that, if we simply split every encroached */
/* segment in its center, two adjacent segments with a small */
/* angle between them might lead to an infinite loop; each */
/* vertex added to split one segment will encroach upon the */
/* other segment, which must then be split with a vertex that */
/* will encroach upon the first segment, and so on forever. */
/* To avoid this, imagine a set of concentric circles, whose */
/* radii are powers of two, about each segment endpoint. */
/* These concentric circles determine where the segment is */
/* split. (If both endpoints are shared with adjacent */
/* segments, split the segment in the middle, and apply the */
/* concentric circles for later splittings.) */
/* Is the origin shared with another segment? */
stpivot(currentenc, enctri);
lnext(enctri, testtri);
tspivot(testtri, testsh);
acuteorg = testsh.ss != m->dummysub;
/* Is the destination shared with another segment? */
lnextself(testtri);
tspivot(testtri, testsh);
acutedest = testsh.ss != m->dummysub;
/* If we're using Chew's algorithm (rather than Ruppert's) */
/* to define encroachment, delete free vertices from the */
/* subsegment's diametral circle. */
if (!b->conformdel && !acuteorg && !acutedest) {
apex(enctri, eapex);
while ((vertextype(eapex) == FREEVERTEX) &&
((eorg[0] - eapex[0]) * (edest[0] - eapex[0]) +
(eorg[1] - eapex[1]) * (edest[1] - eapex[1]) < 0.0)) {
deletevertex(m, b, &testtri);
stpivot(currentenc, enctri);
apex(enctri, eapex);
lprev(enctri, testtri);
}
}
/* Now, check the other side of the segment, if there's a triangle */
/* there. */
sym(enctri, testtri);
if (testtri.tri != m->dummytri) {
/* Is the destination shared with another segment? */
lnextself(testtri);
tspivot(testtri, testsh);
acutedest2 = testsh.ss != m->dummysub;
acutedest = acutedest || acutedest2;
/* Is the origin shared with another segment? */
lnextself(testtri);
tspivot(testtri, testsh);
acuteorg2 = testsh.ss != m->dummysub;
acuteorg = acuteorg || acuteorg2;
/* Delete free vertices from the subsegment's diametral circle. */
if (!b->conformdel && !acuteorg2 && !acutedest2) {
org(testtri, eapex);
while ((vertextype(eapex) == FREEVERTEX) &&
((eorg[0] - eapex[0]) * (edest[0] - eapex[0]) +
(eorg[1] - eapex[1]) * (edest[1] - eapex[1]) < 0.0)) {
deletevertex(m, b, &testtri);
sym(enctri, testtri);
apex(testtri, eapex);
lprevself(testtri);
}
}
}
/* Use the concentric circles if exactly one endpoint is shared */
/* with another adjacent segment. */
if (acuteorg || acutedest) {
segmentlength = sqrt((edest[0] - eorg[0]) * (edest[0] - eorg[0]) +
(edest[1] - eorg[1]) * (edest[1] - eorg[1]));
/* Find the power of two that most evenly splits the segment. */
/* The worst case is a 2:1 ratio between subsegment lengths. */
nearestpoweroftwo = 1.0;
while (segmentlength > 3.0 * nearestpoweroftwo) {
nearestpoweroftwo *= 2.0;
}
while (segmentlength < 1.5 * nearestpoweroftwo) {
nearestpoweroftwo *= 0.5;
}
/* Where do we split the segment? */
split = nearestpoweroftwo / segmentlength;
if (acutedest) {
split = 1.0 - split;
}
} else {
/* If we're not worried about adjacent segments, split */
/* this segment in the middle. */
split = 0.5;
}
/* Create the new vertex. */
newvertex = (vertex) poolalloc(&m->vertices);
/* Interpolate its coordinate and attributes. */
for (i = 0; i < 2 + m->nextras; i++) {
newvertex[i] = eorg[i] + split * (edest[i] - eorg[i]);
}
if (!b->noexact) {
/* Roundoff in the above calculation may yield a `newvertex' */
/* that is not precisely collinear with `eorg' and `edest'. */
/* Improve collinearity by one step of iterative refinement. */
multiplier = counterclockwise(m, b, eorg, edest, newvertex);
divisor = ((eorg[0] - edest[0]) * (eorg[0] - edest[0]) +
(eorg[1] - edest[1]) * (eorg[1] - edest[1]));
if ((multiplier != 0.0) && (divisor != 0.0)) {
multiplier = multiplier / divisor;
/* Watch out for NANs. */
if (multiplier == multiplier) {
newvertex[0] += multiplier * (edest[1] - eorg[1]);
newvertex[1] += multiplier * (eorg[0] - edest[0]);
}
}
}
setvertexmark(newvertex, mark(currentenc));
setvertextype(newvertex, SEGMENTVERTEX);
if (b->verbose > 1) {
printf(
" Splitting subsegment (%.12g, %.12g) (%.12g, %.12g) at (%.12g, %.12g).n",
eorg[0], eorg[1], edest[0], edest[1],
newvertex[0], newvertex[1]);
}
/* Check whether the new vertex lies on an endpoint. */
if (((newvertex[0] == eorg[0]) && (newvertex[1] == eorg[1])) ||
((newvertex[0] == edest[0]) && (newvertex[1] == edest[1]))) {
printf("Error: Ran out of precision at (%.12g, %.12g).n",
newvertex[0], newvertex[1]);
printf("I attempted to split a segment to a smaller size thann");
printf(" can be accommodated by the finite precision ofn");
printf(" floating point arithmetic.n");
precisionerror();
triexit(1);
}
/* Insert the splitting vertex. This should always succeed. */
success = insertvertex(m, b, newvertex, &enctri, &currentenc,
1, triflaws);
if ((success != SUCCESSFULVERTEX) && (success != ENCROACHINGVERTEX)) {
printf("Internal error in splitencsegs():n");
printf(" Failure to split a segment.n");
internalerror();
}
if (m->steinerleft > 0) {
m->steinerleft--;
}
/* Check the two new subsegments to see if they're encroached. */
dummy = checkseg4encroach(m, b, &currentenc);
snextself(currentenc);
dummy = checkseg4encroach(m, b, &currentenc);
}
badsubsegdealloc(m, encloop);
encloop = badsubsegtraverse(m);
}
}
}
#endif /* not CDT_ONLY */
/*****************************************************************************/
/* */
/* tallyfaces() Test every triangle in the mesh for quality measures. */
/* */
/*****************************************************************************/
#ifndef CDT_ONLY
#ifdef ANSI_DECLARATORS
void tallyfaces(struct mesh *m, struct behavior *b)
#else /* not ANSI_DECLARATORS */
void tallyfaces(m, b)
struct mesh *m;
struct behavior *b;
#endif /* not ANSI_DECLARATORS */
{
struct otri triangleloop;
if (b->verbose) {
printf(" Making a list of bad triangles.n");
}
traversalinit(&m->triangles);
triangleloop.orient = 0;
triangleloop.tri = triangletraverse(m);
while (triangleloop.tri != (triangle *) NULL) {
/* If the triangle is bad, enqueue it. */
testtriangle(m, b, &triangleloop);
triangleloop.tri = triangletraverse(m);
}
}
#endif /* not CDT_ONLY */
/*****************************************************************************/
/* */
/* splittriangle() Inserts a vertex at the circumcenter of a triangle. */
/* Deletes the newly inserted vertex if it encroaches */
/* upon a segment. */
/* */
/*****************************************************************************/
#ifndef CDT_ONLY
#ifdef ANSI_DECLARATORS
void splittriangle(struct mesh *m, struct behavior *b,
struct badtriang *badtri)
#else /* not ANSI_DECLARATORS */
void splittriangle(m, b, badtri)
struct mesh *m;
struct behavior *b;
struct badtriang *badtri;
#endif /* not ANSI_DECLARATORS */
{
struct otri badotri;
vertex borg, bdest, bapex;
vertex newvertex;
REAL xi, eta;
enum insertvertexresult success;
int errorflag;
int i;
decode(badtri->poortri, badotri);
org(badotri, borg);
dest(badotri, bdest);
apex(badotri, bapex);
/* Make sure that this triangle is still the same triangle it was */
/* when it was tested and determined to be of bad quality. */
/* Subsequent transformations may have made it a different triangle. */
if (!deadtri(badotri.tri) && (borg == badtri->triangorg) &&
(bdest == badtri->triangdest) && (bapex == badtri->triangapex)) {
if (b->verbose > 1) {
printf(" Splitting this triangle at its circumcenter:n");
printf(" (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)n", borg[0],
borg[1], bdest[0], bdest[1], bapex[0], bapex[1]);
}
errorflag = 0;
/* Create a new vertex at the triangle's circumcenter. */
newvertex = (vertex) poolalloc(&m->vertices);
findcircumcenter(m, b, borg, bdest, bapex, newvertex, &xi, &eta, 1);
/* Check whether the new vertex lies on a triangle vertex. */
if (((newvertex[0] == borg[0]) && (newvertex[1] == borg[1])) ||
((newvertex[0] == bdest[0]) && (newvertex[1] == bdest[1])) ||
((newvertex[0] == bapex[0]) && (newvertex[1] == bapex[1]))) {
if (!b->quiet) {
printf(
"Warning: New vertex (%.12g, %.12g) falls on existing vertex.n",
newvertex[0], newvertex[1]);
errorflag = 1;
}
vertexdealloc(m, newvertex);
} else {
for (i = 2; i < 2 + m->nextras; i++) {
/* Interpolate the vertex attributes at the circumcenter. */
newvertex[i] = borg[i] + xi * (bdest[i] - borg[i])
+ eta * (bapex[i] - borg[i]);
}
/* The new vertex must be in the interior, and therefore is a */
/* free vertex with a marker of zero. */
setvertexmark(newvertex, 0);
setvertextype(newvertex, FREEVERTEX);
/* Ensure that the handle `badotri' does not represent the longest */
/* edge of the triangle. This ensures that the circumcenter must */
/* fall to the left of this edge, so point location will work. */
/* (If the angle org-apex-dest exceeds 90 degrees, then the */
/* circumcenter lies outside the org-dest edge, and eta is */
/* negative. Roundoff error might prevent eta from being */
/* negative when it should be, so I test eta against xi.) */
if (eta < xi) {
lprevself(badotri);
}
/* Insert the circumcenter, searching from the edge of the triangle, */
/* and maintain the Delaunay property of the triangulation. */
success = insertvertex(m, b, newvertex, &badotri, (struct osub *) NULL,
1, 1);
if (success == SUCCESSFULVERTEX) {
if (m->steinerleft > 0) {
m->steinerleft--;
}
} else if (success == ENCROACHINGVERTEX) {
/* If the newly inserted vertex encroaches upon a subsegment, */
/* delete the new vertex. */
undovertex(m, b);
if (b->verbose > 1) {
printf(" Rejecting (%.12g, %.12g).n", newvertex[0], newvertex[1]);
}
vertexdealloc(m, newvertex);
} else if (success == VIOLATINGVERTEX) {
/* Failed to insert the new vertex, but some subsegment was */
/* marked as being encroached. */
vertexdealloc(m, newvertex);
} else { /* success == DUPLICATEVERTEX */
/* Couldn't insert the new vertex because a vertex is already there. */
if (!b->quiet) {
printf(
"Warning: New vertex (%.12g, %.12g) falls on existing vertex.n",
newvertex[0], newvertex[1]);
errorflag = 1;
}
vertexdealloc(m, newvertex);
}
}
if (errorflag) {
if (b->verbose) {
printf(" The new vertex is at the circumcenter of trianglen");
printf(" (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)n",
borg[0], borg[1], bdest[0], bdest[1], bapex[0], bapex[1]);
}
printf("This probably means that I am trying to refine trianglesn");
printf(" to a smaller size than can be accommodated by the finiten");
printf(" precision of floating point arithmetic. (You can ben");
printf(" sure of this if I fail to terminate.)n");
precisionerror();
}
}
}
#endif /* not CDT_ONLY */
/*****************************************************************************/
/* */
/* enforcequality() Remove all the encroached subsegments and bad */
/* triangles from the triangulation. */
/* */
/*****************************************************************************/
#ifndef CDT_ONLY
#ifdef ANSI_DECLARATORS
void enforcequality(struct mesh *m, struct behavior *b)
#else /* not ANSI_DECLARATORS */
void enforcequality(m, b)
struct mesh *m;
struct behavior *b;
#endif /* not ANSI_DECLARATORS */
{
struct badtriang *badtri;
int i;
if (!b->quiet) {
printf("Adding Steiner points to enforce quality.n");
}
/* Initialize the pool of encroached subsegments. */
poolinit(&m->badsubsegs, sizeof(struct badsubseg), BADSUBSEGPERBLOCK,
BADSUBSEGPERBLOCK, 0);
if (b->verbose) {
printf(" Looking for encroached subsegments.n");
}
/* Test all segments to see if they're encroached. */
tallyencs(m, b);
if (b->verbose && (m->badsubsegs.items > 0)) {
printf(" Splitting encroached subsegments.n");
}
/* Fix encroached subsegments without noting bad triangles. */
splitencsegs(m, b, 0);
/* At this point, if we haven't run out of Steiner points, the */
/* triangulation should be (conforming) Delaunay. */
/* Next, we worry about enforcing triangle quality. */
if ((b->minangle > 0.0) || b->vararea || b->fixedarea || b->usertest) {
/* Initialize the pool of bad triangles. */
poolinit(&m->badtriangles, sizeof(struct badtriang), BADTRIPERBLOCK,
BADTRIPERBLOCK, 0);
/* Initialize the queues of bad triangles. */
for (i = 0; i < 4096; i++) {
m->queuefront[i] = (struct badtriang *) NULL;
}
m->firstnonemptyq = -1;
/* Test all triangles to see if they're bad. */
tallyfaces(m, b);
/* Initialize the pool of recently flipped triangles. */
poolinit(&m->flipstackers, sizeof(struct flipstacker), FLIPSTACKERPERBLOCK,
FLIPSTACKERPERBLOCK, 0);
m->checkquality = 1;
if (b->verbose) {
printf(" Splitting bad triangles.n");
}
while ((m->badtriangles.items > 0) && (m->steinerleft != 0)) {
/* Fix one bad triangle by inserting a vertex at its circumcenter. */
badtri = dequeuebadtriang(m);
splittriangle(m, b, badtri);
if (m->badsubsegs.items > 0) {
/* Put bad triangle back in queue for another try later. */
enqueuebadtriang(m, b, badtri);
/* Fix any encroached subsegments that resulted. */
/* Record any new bad triangles that result. */
splitencsegs(m, b, 1);
} else {
/* Return the bad triangle to the pool. */
pooldealloc(&m->badtriangles, (VOID *) badtri);
}
}
}
/* At this point, if the "-D" switch was selected and we haven't run out */
/* of Steiner points, the triangulation should be (conforming) Delaunay */
/* and have no low-quality triangles. */
/* Might we have run out of Steiner points too soon? */
if (!b->quiet && b->conformdel && (m->badsubsegs.items > 0) &&
(m->steinerleft == 0)) {
printf("nWarning: I ran out of Steiner points, but the mesh hasn");
if (m->badsubsegs.items == 1) {
printf(" one encroached subsegment, and therefore might not be trulyn"
);
} else {
printf(" %ld encroached subsegments, and therefore might not be trulyn"
, m->badsubsegs.items);
}
printf(" Delaunay. If the Delaunay property is important to you,n");
printf(" try increasing the number of Steiner points (controlled byn");
printf(" the -S switch) slightly and try again.nn");
}
}
#endif /* not CDT_ONLY */
/** **/
/** **/
/********* Mesh quality maintenance ends here *********/
/*****************************************************************************/
/* */
/* highorder() Create extra nodes for quadratic subparametric elements. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void highorder(struct mesh *m, struct behavior *b)
#else /* not ANSI_DECLARATORS */
void highorder(m, b)
struct mesh *m;
struct behavior *b;
#endif /* not ANSI_DECLARATORS */
{
struct otri triangleloop, trisym;
struct osub checkmark;
vertex newvertex;
vertex torg, tdest;
int i;
triangle ptr; /* Temporary variable used by sym(). */
subseg sptr; /* Temporary variable used by tspivot(). */
if (!b->quiet) {
printf("Adding vertices for second-order triangles.n");
}
/* The following line ensures that dead items in the pool of nodes */
/* cannot be allocated for the extra nodes associated with high */
/* order elements. This ensures that the primary nodes (at the */
/* corners of elements) will occur earlier in the output files, and */
/* have lower indices, than the extra nodes. */
m->vertices.deaditemstack = (VOID *) NULL;
traversalinit(&m->triangles);
triangleloop.tri = triangletraverse(m);
/* To loop over the set of edges, loop over all triangles, and look at */
/* the three edges of each triangle. If there isn't another triangle */
/* adjacent to the edge, operate on the edge. If there is another */
/* adjacent triangle, operate on the edge only if the current triangle */
/* has a smaller pointer than its neighbor. This way, each edge is */
/* considered only once. */
while (triangleloop.tri != (triangle *) NULL) {
for (triangleloop.orient = 0; triangleloop.orient < 3;
triangleloop.orient++) {
sym(triangleloop, trisym);
if ((triangleloop.tri < trisym.tri) || (trisym.tri == m->dummytri)) {
org(triangleloop, torg);
dest(triangleloop, tdest);
/* Create a new node in the middle of the edge. Interpolate */
/* its attributes. */
newvertex = (vertex) poolalloc(&m->vertices);
for (i = 0; i < 2 + m->nextras; i++) {
newvertex[i] = 0.5 * (torg[i] + tdest[i]);
}
/* Set the new node's marker to zero or one, depending on */
/* whether it lies on a boundary. */
setvertexmark(newvertex, trisym.tri == m->dummytri);
setvertextype(newvertex,
trisym.tri == m->dummytri ? FREEVERTEX : SEGMENTVERTEX);
if (b->usesegments) {
tspivot(triangleloop, checkmark);
/* If this edge is a segment, transfer the marker to the new node. */
if (checkmark.ss != m->dummysub) {
setvertexmark(newvertex, mark(checkmark));
setvertextype(newvertex, SEGMENTVERTEX);
}
}
if (b->verbose > 1) {
printf(" Creating (%.12g, %.12g).n", newvertex[0], newvertex[1]);
}
/* Record the new node in the (one or two) adjacent elements. */
triangleloop.tri[m->highorderindex + triangleloop.orient] =
(triangle) newvertex;
if (trisym.tri != m->dummytri) {
trisym.tri[m->highorderindex + trisym.orient] = (triangle) newvertex;
}
}
}
triangleloop.tri = triangletraverse(m);
}
}
/********* File I/O routines begin here *********/
/** **/
/** **/
/*****************************************************************************/
/* */
/* readline() Read a nonempty line from a file. */
/* */
/* A line is considered "nonempty" if it contains something that looks like */
/* a number. Comments (prefaced by `#') are ignored. */
/* */
/*****************************************************************************/
#ifndef TRILIBRARY
#ifdef ANSI_DECLARATORS
char *readline(char *string, FILE *infile, char *infilename)
#else /* not ANSI_DECLARATORS */
char *readline(string, infile, infilename)
char *string;
FILE *infile;
char *infilename;
#endif /* not ANSI_DECLARATORS */
{
char *result;
/* Search for something that looks like a number. */
do {
result = fgets(string, INPUTLINESIZE, infile);
if (result == (char *) NULL) {
printf(" Error: Unexpected end of file in %s.n", infilename);
triexit(1);
}
/* Skip anything that doesn't look like a number, a comment, */
/* or the end of a line. */
while ((*result != '') && (*result != '#')
&& (*result != '.') && (*result != '+') && (*result != '-')
&& ((*result < '0') || (*result > '9'))) {
result++;
}
/* If it's a comment or end of line, read another line and try again. */
} while ((*result == '#') || (*result == ''));
return result;
}
#endif /* not TRILIBRARY */
/*****************************************************************************/
/* */
/* findfield() Find the next field of a string. */
/* */
/* Jumps past the current field by searching for whitespace, then jumps */
/* past the whitespace to find the next field. */
/* */
/*****************************************************************************/
#ifndef TRILIBRARY
#ifdef ANSI_DECLARATORS
char *findfield(char *string)
#else /* not ANSI_DECLARATORS */
char *findfield(string)
char *string;
#endif /* not ANSI_DECLARATORS */
{
char *result;
result = string;
/* Skip the current field. Stop upon reaching whitespace. */
while ((*result != '') && (*result != '#')
&& (*result != ' ') && (*result != 't')) {
result++;
}
/* Now skip the whitespace and anything else that doesn't look like a */
/* number, a comment, or the end of a line. */
while ((*result != '') && (*result != '#')
&& (*result != '.') && (*result != '+') && (*result != '-')
&& ((*result < '0') || (*result > '9'))) {
result++;
}
/* Check for a comment (prefixed with `#'). */
if (*result == '#') {
*result = '';
}
return result;
}
#endif /* not TRILIBRARY */
/*****************************************************************************/
/* */
/* readnodes() Read the vertices from a file, which may be a .node or */
/* .poly file. */
/* */
/*****************************************************************************/
#ifndef TRILIBRARY
#ifdef ANSI_DECLARATORS
void readnodes(struct mesh *m, struct behavior *b, char *nodefilename,
char *polyfilename, FILE **polyfile)
#else /* not ANSI_DECLARATORS */
void readnodes(m, b, nodefilename, polyfilename, polyfile)
struct mesh *m;
struct behavior *b;
char *nodefilename;
char *polyfilename;
FILE **polyfile;
#endif /* not ANSI_DECLARATORS */
{
FILE *infile;
vertex vertexloop;
char inputline[INPUTLINESIZE];
char *stringptr;
char *infilename;
REAL x, y;
int firstnode;
int nodemarkers;
int currentmarker;
int i, j;
if (b->poly) {
/* Read the vertices from a .poly file. */
if (!b->quiet) {
printf("Opening %s.n", polyfilename);
}
*polyfile = fopen(polyfilename, "r");
if (*polyfile == (FILE *) NULL) {
printf(" Error: Cannot access file %s.n", polyfilename);
triexit(1);
}
/* Read number of vertices, number of dimensions, number of vertex */
/* attributes, and number of boundary markers. */
stringptr = readline(inputline, *polyfile, polyfilename);
m->invertices = (int) strtol(stringptr, &stringptr, 0);
stringptr = findfield(stringptr);
if (*stringptr == '') {
m->mesh_dim = 2;
} else {
m->mesh_dim = (int) strtol(stringptr, &stringptr, 0);
}
stringptr = findfield(stringptr);
if (*stringptr == '') {
m->nextras = 0;
} else {
m->nextras = (int) strtol(stringptr, &stringptr, 0);
}
stringptr = findfield(stringptr);
if (*stringptr == '') {
nodemarkers = 0;
} else {
nodemarkers = (int) strtol(stringptr, &stringptr, 0);
}
if (m->invertices > 0) {
infile = *polyfile;
infilename = polyfilename;
m->readnodefile = 0;
} else {
/* If the .poly file claims there are zero vertices, that means that */
/* the vertices should be read from a separate .node file. */
m->readnodefile = 1;
infilename = nodefilename;
}
} else {
m->readnodefile = 1;
infilename = nodefilename;
*polyfile = (FILE *) NULL;
}
if (m->readnodefile) {
/* Read the vertices from a .node file. */
if (!b->quiet) {
printf("Opening %s.n", nodefilename);
}
infile = fopen(nodefilename, "r");
if (infile == (FILE *) NULL) {
printf(" Error: Cannot access file %s.n", nodefilename);
triexit(1);
}
/* Read number of vertices, number of dimensions, number of vertex */
/* attributes, and number of boundary markers. */
stringptr = readline(inputline, infile, nodefilename);
m->invertices = (int) strtol(stringptr, &stringptr, 0);
stringptr = findfield(stringptr);
if (*stringptr == '') {
m->mesh_dim = 2;
} else {
m->mesh_dim = (int) strtol(stringptr, &stringptr, 0);
}
stringptr = findfield(stringptr);
if (*stringptr == '') {
m->nextras = 0;
} else {
m->nextras = (int) strtol(stringptr, &stringptr, 0);
}
stringptr = findfield(stringptr);
if (*stringptr == '') {
nodemarkers = 0;
} else {
nodemarkers = (int) strtol(stringptr, &stringptr, 0);
}
}
if (m->invertices < 3) {
printf("Error: Input must have at least three input vertices.n");
triexit(1);
}
if (m->mesh_dim != 2) {
printf("Error: Triangle only works with two-dimensional meshes.n");
triexit(1);
}
if (m->nextras == 0) {
b->weighted = 0;
}
initializevertexpool(m, b);
/* Read the vertices. */
for (i = 0; i < m->invertices; i++) {
vertexloop = (vertex) poolalloc(&m->vertices);
stringptr = readline(inputline, infile, infilename);
if (i == 0) {
firstnode = (int) strtol(stringptr, &stringptr, 0);
if ((firstnode == 0) || (firstnode == 1)) {
b->firstnumber = firstnode;
}
}
stringptr = findfield(stringptr);
if (*stringptr == '') {
printf("Error: Vertex %d has no x coordinate.n", b->firstnumber + i);
triexit(1);
}
x = (REAL) strtod(stringptr, &stringptr);
stringptr = findfield(stringptr);
if (*stringptr == '') {
printf("Error: Vertex %d has no y coordinate.n", b->firstnumber + i);
triexit(1);
}
y = (REAL) strtod(stringptr, &stringptr);
vertexloop[0] = x;
vertexloop[1] = y;
/* Read the vertex attributes. */
for (j = 2; j < 2 + m->nextras; j++) {
stringptr = findfield(stringptr);
if (*stringptr == '') {
vertexloop[j] = 0.0;
} else {
vertexloop[j] = (REAL) strtod(stringptr, &stringptr);
}
}
if (nodemarkers) {
/* Read a vertex marker. */
stringptr = findfield(stringptr);
if (*stringptr == '') {
setvertexmark(vertexloop, 0);
} else {
currentmarker = (int) strtol(stringptr, &stringptr, 0);
setvertexmark(vertexloop, currentmarker);
}
} else {
/* If no markers are specified in the file, they default to zero. */
setvertexmark(vertexloop, 0);
}
setvertextype(vertexloop, INPUTVERTEX);
/* Determine the smallest and largest x and y coordinates. */
if (i == 0) {
m->xmin = m->xmax = x;
m->ymin = m->ymax = y;
} else {
m->xmin = (x < m->xmin) ? x : m->xmin;
m->xmax = (x > m->xmax) ? x : m->xmax;
m->ymin = (y < m->ymin) ? y : m->ymin;
m->ymax = (y > m->ymax) ? y : m->ymax;
}
}
if (m->readnodefile) {
fclose(infile);
}
/* Nonexistent x value used as a flag to mark circle events in sweepline */
/* Delaunay algorithm. */
m->xminextreme = 10 * m->xmin - 9 * m->xmax;
}
#endif /* not TRILIBRARY */
/*****************************************************************************/
/* */
/* transfernodes() Read the vertices from memory. */
/* */
/*****************************************************************************/
#ifdef TRILIBRARY
#ifdef ANSI_DECLARATORS
void transfernodes(struct mesh *m, struct behavior *b, REAL *pointlist,
REAL *pointattriblist, int *pointmarkerlist,
int numberofpoints, int numberofpointattribs)
#else /* not ANSI_DECLARATORS */
void transfernodes(m, b, pointlist, pointattriblist, pointmarkerlist,
numberofpoints, numberofpointattribs)
struct mesh *m;
struct behavior *b;
REAL *pointlist;
REAL *pointattriblist;
int *pointmarkerlist;
int numberofpoints;
int numberofpointattribs;
#endif /* not ANSI_DECLARATORS */
{
vertex vertexloop;
REAL x, y;
int i, j;
int coordindex;
int attribindex;
m->invertices = numberofpoints;
m->mesh_dim = 2;
m->nextras = numberofpointattribs;
m->readnodefile = 0;
if (m->invertices < 3) {
printf("Error: Input must have at least three input vertices.n");
triexit(1);
}
if (m->nextras == 0) {
b->weighted = 0;
}
initializevertexpool(m, b);
/* Read the vertices. */
coordindex = 0;
attribindex = 0;
for (i = 0; i < m->invertices; i++) {
vertexloop = (vertex) poolalloc(&m->vertices);
/* Read the vertex coordinates. */
x = vertexloop[0] = pointlist[coordindex++];
y = vertexloop[1] = pointlist[coordindex++];
/* Read the vertex attributes. */
for (j = 0; j < numberofpointattribs; j++) {
vertexloop[2 + j] = pointattriblist[attribindex++];
}
if (pointmarkerlist != (int *) NULL) {
/* Read a vertex marker. */
setvertexmark(vertexloop, pointmarkerlist[i]);
} else {
/* If no markers are specified, they default to zero. */
setvertexmark(vertexloop, 0);
}
setvertextype(vertexloop, INPUTVERTEX);
/* Determine the smallest and largest x and y coordinates. */
if (i == 0) {
m->xmin = m->xmax = x;
m->ymin = m->ymax = y;
} else {
m->xmin = (x < m->xmin) ? x : m->xmin;
m->xmax = (x > m->xmax) ? x : m->xmax;
m->ymin = (y < m->ymin) ? y : m->ymin;
m->ymax = (y > m->ymax) ? y : m->ymax;
}
}
/* Nonexistent x value used as a flag to mark circle events in sweepline */
/* Delaunay algorithm. */
m->xminextreme = 10 * m->xmin - 9 * m->xmax;
}
#endif /* TRILIBRARY */
/*****************************************************************************/
/* */
/* readholes() Read the holes, and possibly regional attributes and area */
/* constraints, from a .poly file. */
/* */
/*****************************************************************************/
#ifndef TRILIBRARY
#ifdef ANSI_DECLARATORS
void readholes(struct mesh *m, struct behavior *b,
FILE *polyfile, char *polyfilename, REAL **hlist, int *holes,
REAL **rlist, int *regions)
#else /* not ANSI_DECLARATORS */
void readholes(m, b, polyfile, polyfilename, hlist, holes, rlist, regions)
struct mesh *m;
struct behavior *b;
FILE *polyfile;
char *polyfilename;
REAL **hlist;
int *holes;
REAL **rlist;
int *regions;
#endif /* not ANSI_DECLARATORS */
{
REAL *holelist;
REAL *regionlist;
char inputline[INPUTLINESIZE];
char *stringptr;
int index;
int i;
/* Read the holes. */
stringptr = readline(inputline, polyfile, polyfilename);
*holes = (int) strtol(stringptr, &stringptr, 0);
if (*holes > 0) {
holelist = (REAL *) trimalloc(2 * *holes * (int) sizeof(REAL));
*hlist = holelist;
for (i = 0; i < 2 * *holes; i += 2) {
stringptr = readline(inputline, polyfile, polyfilename);
stringptr = findfield(stringptr);
if (*stringptr == '') {
printf("Error: Hole %d has no x coordinate.n",
b->firstnumber + (i >> 1));
triexit(1);
} else {
holelist[i] = (REAL) strtod(stringptr, &stringptr);
}
stringptr = findfield(stringptr);
if (*stringptr == '') {
printf("Error: Hole %d has no y coordinate.n",
b->firstnumber + (i >> 1));
triexit(1);
} else {
holelist[i + 1] = (REAL) strtod(stringptr, &stringptr);
}
}
} else {
*hlist = (REAL *) NULL;
}
#ifndef CDT_ONLY
if ((b->regionattrib || b->vararea) && !b->refine) {
/* Read the area constraints. */
stringptr = readline(inputline, polyfile, polyfilename);
*regions = (int) strtol(stringptr, &stringptr, 0);
if (*regions > 0) {
regionlist = (REAL *) trimalloc(4 * *regions * (int) sizeof(REAL));
*rlist = regionlist;
index = 0;
for (i = 0; i < *regions; i++) {
stringptr = readline(inputline, polyfile, polyfilename);
stringptr = findfield(stringptr);
if (*stringptr == '') {
printf("Error: Region %d has no x coordinate.n",
b->firstnumber + i);
triexit(1);
} else {
regionlist[index++] = (REAL) strtod(stringptr, &stringptr);
}
stringptr = findfield(stringptr);
if (*stringptr == '') {
printf("Error: Region %d has no y coordinate.n",
b->firstnumber + i);
triexit(1);
} else {
regionlist[index++] = (REAL) strtod(stringptr, &stringptr);
}
stringptr = findfield(stringptr);
if (*stringptr == '') {
printf(
"Error: Region %d has no region attribute or area constraint.n",
b->firstnumber + i);
triexit(1);
} else {
regionlist[index++] = (REAL) strtod(stringptr, &stringptr);
}
stringptr = findfield(stringptr);
if (*stringptr == '') {
regionlist[index] = regionlist[index - 1];
} else {
regionlist[index] = (REAL) strtod(stringptr, &stringptr);
}
index++;
}
}
} else {
/* Set `*regions' to zero to avoid an accidental free() later. */
*regions = 0;
*rlist = (REAL *) NULL;
}
#endif /* not CDT_ONLY */
fclose(polyfile);
}
#endif /* not TRILIBRARY */
/*****************************************************************************/
/* */
/* finishfile() Write the command line to the output file so the user */
/* can remember how the file was generated. Close the file. */
/* */
/*****************************************************************************/
#ifndef TRILIBRARY
#ifdef ANSI_DECLARATORS
void finishfile(FILE *outfile, int argc, char **argv)
#else /* not ANSI_DECLARATORS */
void finishfile(outfile, argc, argv)
FILE *outfile;
int argc;
char **argv;
#endif /* not ANSI_DECLARATORS */
{
int i;
fprintf(outfile, "# Generated by");
for (i = 0; i < argc; i++) {
fprintf(outfile, " ");
fputs(argv[i], outfile);
}
fprintf(outfile, "n");
fclose(outfile);
}
#endif /* not TRILIBRARY */
/*****************************************************************************/
/* */
/* writenodes() Number the vertices and write them to a .node file. */
/* */
/* To save memory, the vertex numbers are written over the boundary markers */
/* after the vertices are written to a file. */
/* */
/*****************************************************************************/
#ifdef TRILIBRARY
#ifdef ANSI_DECLARATORS
void writenodes(struct mesh *m, struct behavior *b, REAL **pointlist,
REAL **pointattriblist, int **pointmarkerlist)
#else /* not ANSI_DECLARATORS */
void writenodes(m, b, pointlist, pointattriblist, pointmarkerlist)
struct mesh *m;
struct behavior *b;
REAL **pointlist;
REAL **pointattriblist;
int **pointmarkerlist;
#endif /* not ANSI_DECLARATORS */
#else /* not TRILIBRARY */
#ifdef ANSI_DECLARATORS
void writenodes(struct mesh *m, struct behavior *b, char *nodefilename,
int argc, char **argv)
#else /* not ANSI_DECLARATORS */
void writenodes(m, b, nodefilename, argc, argv)
struct mesh *m;
struct behavior *b;
char *nodefilename;
int argc;
char **argv;
#endif /* not ANSI_DECLARATORS */
#endif /* not TRILIBRARY */
{
#ifdef TRILIBRARY
REAL *plist;
REAL *palist;
int *pmlist;
int coordindex;
int attribindex;
#else /* not TRILIBRARY */
FILE *outfile;
#endif /* not TRILIBRARY */
vertex vertexloop;
long outvertices;
int vertexnumber;
int i;
if (b->jettison) {
outvertices = m->vertices.items - m->undeads;
} else {
outvertices = m->vertices.items;
}
#ifdef TRILIBRARY
if (!b->quiet) {
printf("Writing vertices.n");
}
/* Allocate memory for output vertices if necessary. */
if (*pointlist == (REAL *) NULL) {
*pointlist = (REAL *) trimalloc((int) (outvertices * 2 * sizeof(REAL)));
}
/* Allocate memory for output vertex attributes if necessary. */
if ((m->nextras > 0) && (*pointattriblist == (REAL *) NULL)) {
*pointattriblist = (REAL *) trimalloc((int) (outvertices * m->nextras *
sizeof(REAL)));
}
/* Allocate memory for output vertex markers if necessary. */
if (!b->nobound && (*pointmarkerlist == (int *) NULL)) {
*pointmarkerlist = (int *) trimalloc((int) (outvertices * sizeof(int)));
}
plist = *pointlist;
palist = *pointattriblist;
pmlist = *pointmarkerlist;
coordindex = 0;
attribindex = 0;
#else /* not TRILIBRARY */
if (!b->quiet) {
printf("Writing %s.n", nodefilename);
}
outfile = fopen(nodefilename, "w");
if (outfile == (FILE *) NULL) {
printf(" Error: Cannot create file %s.n", nodefilename);
triexit(1);
}
/* Number of vertices, number of dimensions, number of vertex attributes, */
/* and number of boundary markers (zero or one). */
fprintf(outfile, "%ld %d %d %dn", outvertices, m->mesh_dim,
m->nextras, 1 - b->nobound);
#endif /* not TRILIBRARY */
traversalinit(&m->vertices);
vertexnumber = b->firstnumber;
vertexloop = vertextraverse(m);
while (vertexloop != (vertex) NULL) {
if (!b->jettison || (vertextype(vertexloop) != UNDEADVERTEX)) {
#ifdef TRILIBRARY
/* X and y coordinates. */
plist[coordindex++] = vertexloop[0];
plist[coordindex++] = vertexloop[1];
/* Vertex attributes. */
for (i = 0; i < m->nextras; i++) {
palist[attribindex++] = vertexloop[2 + i];
}
if (!b->nobound) {
/* Copy the boundary marker. */
pmlist[vertexnumber - b->firstnumber] = vertexmark(vertexloop);
}
#else /* not TRILIBRARY */
/* Vertex number, x and y coordinates. */
fprintf(outfile, "%4d %.17g %.17g", vertexnumber, vertexloop[0],
vertexloop[1]);
for (i = 0; i < m->nextras; i++) {
/* Write an attribute. */
fprintf(outfile, " %.17g", vertexloop[i + 2]);
}
if (b->nobound) {
fprintf(outfile, "n");
} else {
/* Write the boundary marker. */
fprintf(outfile, " %dn", vertexmark(vertexloop));
}
#endif /* not TRILIBRARY */
setvertexmark(vertexloop, vertexnumber);
vertexnumber++;
}
vertexloop = vertextraverse(m);
}
#ifndef TRILIBRARY
finishfile(outfile, argc, argv);
#endif /* not TRILIBRARY */
}
/*****************************************************************************/
/* */
/* numbernodes() Number the vertices. */
/* */
/* Each vertex is assigned a marker equal to its number. */
/* */
/* Used when writenodes() is not called because no .node file is written. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void numbernodes(struct mesh *m, struct behavior *b)
#else /* not ANSI_DECLARATORS */
void numbernodes(m, b)
struct mesh *m;
struct behavior *b;
#endif /* not ANSI_DECLARATORS */
{
vertex vertexloop;
int vertexnumber;
traversalinit(&m->vertices);
vertexnumber = b->firstnumber;
vertexloop = vertextraverse(m);
while (vertexloop != (vertex) NULL) {
setvertexmark(vertexloop, vertexnumber);
if (!b->jettison || (vertextype(vertexloop) != UNDEADVERTEX)) {
vertexnumber++;
}
vertexloop = vertextraverse(m);
}
}
/*****************************************************************************/
/* */
/* writeelements() Write the triangles to an .ele file. */
/* */
/*****************************************************************************/
#ifdef TRILIBRARY
#ifdef ANSI_DECLARATORS
void writeelements(struct mesh *m, struct behavior *b,
int **trianglelist, REAL **triangleattriblist)
#else /* not ANSI_DECLARATORS */
void writeelements(m, b, trianglelist, triangleattriblist)
struct mesh *m;
struct behavior *b;
int **trianglelist;
REAL **triangleattriblist;
#endif /* not ANSI_DECLARATORS */
#else /* not TRILIBRARY */
#ifdef ANSI_DECLARATORS
void writeelements(struct mesh *m, struct behavior *b, char *elefilename,
int argc, char **argv)
#else /* not ANSI_DECLARATORS */
void writeelements(m, b, elefilename, argc, argv)
struct mesh *m;
struct behavior *b;
char *elefilename;
int argc;
char **argv;
#endif /* not ANSI_DECLARATORS */
#endif /* not TRILIBRARY */
{
#ifdef TRILIBRARY
int *tlist;
REAL *talist;
int vertexindex;
int attribindex;
#else /* not TRILIBRARY */
FILE *outfile;
#endif /* not TRILIBRARY */
struct otri triangleloop;
vertex p1, p2, p3;
vertex mid1, mid2, mid3;
long elementnumber;
int i;
#ifdef TRILIBRARY
if (!b->quiet) {
printf("Writing triangles.n");
}
/* Allocate memory for output triangles if necessary. */
if (*trianglelist == (int *) NULL) {
*trianglelist = (int *) trimalloc((int) (m->triangles.items *
((b->order + 1) * (b->order + 2) /
2) * sizeof(int)));
}
/* Allocate memory for output triangle attributes if necessary. */
if ((m->eextras > 0) && (*triangleattriblist == (REAL *) NULL)) {
*triangleattriblist = (REAL *) trimalloc((int) (m->triangles.items *
m->eextras *
sizeof(REAL)));
}
tlist = *trianglelist;
talist = *triangleattriblist;
vertexindex = 0;
attribindex = 0;
#else /* not TRILIBRARY */
if (!b->quiet) {
printf("Writing %s.n", elefilename);
}
outfile = fopen(elefilename, "w");
if (outfile == (FILE *) NULL) {
printf(" Error: Cannot create file %s.n", elefilename);
triexit(1);
}
/* Number of triangles, vertices per triangle, attributes per triangle. */
fprintf(outfile, "%ld %d %dn", m->triangles.items,
(b->order + 1) * (b->order + 2) / 2, m->eextras);
#endif /* not TRILIBRARY */
traversalinit(&m->triangles);
triangleloop.tri = triangletraverse(m);
triangleloop.orient = 0;
elementnumber = b->firstnumber;
while (triangleloop.tri != (triangle *) NULL) {
org(triangleloop, p1);
dest(triangleloop, p2);
apex(triangleloop, p3);
if (b->order == 1) {
#ifdef TRILIBRARY
tlist[vertexindex++] = vertexmark(p1);
tlist[vertexindex++] = vertexmark(p2);
tlist[vertexindex++] = vertexmark(p3);
#else /* not TRILIBRARY */
/* Triangle number, indices for three vertices. */
fprintf(outfile, "%4ld %4d %4d %4d", elementnumber,
vertexmark(p1), vertexmark(p2), vertexmark(p3));
#endif /* not TRILIBRARY */
} else {
mid1 = (vertex) triangleloop.tri[m->highorderindex + 1];
mid2 = (vertex) triangleloop.tri[m->highorderindex + 2];
mid3 = (vertex) triangleloop.tri[m->highorderindex];
#ifdef TRILIBRARY
tlist[vertexindex++] = vertexmark(p1);
tlist[vertexindex++] = vertexmark(p2);
tlist[vertexindex++] = vertexmark(p3);
tlist[vertexindex++] = vertexmark(mid1);
tlist[vertexindex++] = vertexmark(mid2);
tlist[vertexindex++] = vertexmark(mid3);
#else /* not TRILIBRARY */
/* Triangle number, indices for six vertices. */
fprintf(outfile, "%4ld %4d %4d %4d %4d %4d %4d", elementnumber,
vertexmark(p1), vertexmark(p2), vertexmark(p3), vertexmark(mid1),
vertexmark(mid2), vertexmark(mid3));
#endif /* not TRILIBRARY */
}
#ifdef TRILIBRARY
for (i = 0; i < m->eextras; i++) {
talist[attribindex++] = elemattribute(triangleloop, i);
}
#else /* not TRILIBRARY */
for (i = 0; i < m->eextras; i++) {
fprintf(outfile, " %.17g", elemattribute(triangleloop, i));
}
fprintf(outfile, "n");
#endif /* not TRILIBRARY */
triangleloop.tri = triangletraverse(m);
elementnumber++;
}
#ifndef TRILIBRARY
finishfile(outfile, argc, argv);
#endif /* not TRILIBRARY */
}
/*****************************************************************************/
/* */
/* writepoly() Write the segments and holes to a .poly file. */
/* */
/*****************************************************************************/
#ifdef TRILIBRARY
#ifdef ANSI_DECLARATORS
void writepoly(struct mesh *m, struct behavior *b,
int **segmentlist, int **segmentmarkerlist)
#else /* not ANSI_DECLARATORS */
void writepoly(m, b, segmentlist, segmentmarkerlist)
struct mesh *m;
struct behavior *b;
int **segmentlist;
int **segmentmarkerlist;
#endif /* not ANSI_DECLARATORS */
#else /* not TRILIBRARY */
#ifdef ANSI_DECLARATORS
void writepoly(struct mesh *m, struct behavior *b, char *polyfilename,
REAL *holelist, int holes, REAL *regionlist, int regions,
int argc, char **argv)
#else /* not ANSI_DECLARATORS */
void writepoly(m, b, polyfilename, holelist, holes, regionlist, regions,
argc, argv)
struct mesh *m;
struct behavior *b;
char *polyfilename;
REAL *holelist;
int holes;
REAL *regionlist;
int regions;
int argc;
char **argv;
#endif /* not ANSI_DECLARATORS */
#endif /* not TRILIBRARY */
{
#ifdef TRILIBRARY
int *slist;
int *smlist;
int index;
#else /* not TRILIBRARY */
FILE *outfile;
long holenumber, regionnumber;
#endif /* not TRILIBRARY */
struct osub subsegloop;
vertex endpoint1, endpoint2;
long subsegnumber;
#ifdef TRILIBRARY
if (!b->quiet) {
printf("Writing segments.n");
}
/* Allocate memory for output segments if necessary. */
if (*segmentlist == (int *) NULL) {
*segmentlist = (int *) trimalloc((int) (m->subsegs.items * 2 *
sizeof(int)));
}
/* Allocate memory for output segment markers if necessary. */
if (!b->nobound && (*segmentmarkerlist == (int *) NULL)) {
*segmentmarkerlist = (int *) trimalloc((int) (m->subsegs.items *
sizeof(int)));
}
slist = *segmentlist;
smlist = *segmentmarkerlist;
index = 0;
#else /* not TRILIBRARY */
if (!b->quiet) {
printf("Writing %s.n", polyfilename);
}
outfile = fopen(polyfilename, "w");
if (outfile == (FILE *) NULL) {
printf(" Error: Cannot create file %s.n", polyfilename);
triexit(1);
}
/* The zero indicates that the vertices are in a separate .node file. */
/* Followed by number of dimensions, number of vertex attributes, */
/* and number of boundary markers (zero or one). */
fprintf(outfile, "%d %d %d %dn", 0, m->mesh_dim, m->nextras,
1 - b->nobound);
/* Number of segments, number of boundary markers (zero or one). */
fprintf(outfile, "%ld %dn", m->subsegs.items, 1 - b->nobound);
#endif /* not TRILIBRARY */
traversalinit(&m->subsegs);
subsegloop.ss = subsegtraverse(m);
subsegloop.ssorient = 0;
subsegnumber = b->firstnumber;
while (subsegloop.ss != (subseg *) NULL) {
sorg(subsegloop, endpoint1);
sdest(subsegloop, endpoint2);
#ifdef TRILIBRARY
/* Copy indices of the segment's two endpoints. */
slist[index++] = vertexmark(endpoint1);
slist[index++] = vertexmark(endpoint2);
if (!b->nobound) {
/* Copy the boundary marker. */
smlist[subsegnumber - b->firstnumber] = mark(subsegloop);
}
#else /* not TRILIBRARY */
/* Segment number, indices of its two endpoints, and possibly a marker. */
if (b->nobound) {
fprintf(outfile, "%4ld %4d %4dn", subsegnumber,
vertexmark(endpoint1), vertexmark(endpoint2));
} else {
fprintf(outfile, "%4ld %4d %4d %4dn", subsegnumber,
vertexmark(endpoint1), vertexmark(endpoint2), mark(subsegloop));
}
#endif /* not TRILIBRARY */
subsegloop.ss = subsegtraverse(m);
subsegnumber++;
}
#ifndef TRILIBRARY
#ifndef CDT_ONLY
fprintf(outfile, "%dn", holes);
if (holes > 0) {
for (holenumber = 0; holenumber < holes; holenumber++) {
/* Hole number, x and y coordinates. */
fprintf(outfile, "%4ld %.17g %.17gn", b->firstnumber + holenumber,
holelist[2 * holenumber], holelist[2 * holenumber + 1]);
}
}
if (regions > 0) {
fprintf(outfile, "%dn", regions);
for (regionnumber = 0; regionnumber < regions; regionnumber++) {
/* Region number, x and y coordinates, attribute, maximum area. */
fprintf(outfile, "%4ld %.17g %.17g %.17g %.17gn",
b->firstnumber + regionnumber,
regionlist[4 * regionnumber], regionlist[4 * regionnumber + 1],
regionlist[4 * regionnumber + 2],
regionlist[4 * regionnumber + 3]);
}
}
#endif /* not CDT_ONLY */
finishfile(outfile, argc, argv);
#endif /* not TRILIBRARY */
}
/*****************************************************************************/
/* */
/* writeedges() Write the edges to an .edge file. */
/* */
/*****************************************************************************/
#ifdef TRILIBRARY
#ifdef ANSI_DECLARATORS
void writeedges(struct mesh *m, struct behavior *b,
int **edgelist, int **edgemarkerlist)
#else /* not ANSI_DECLARATORS */
void writeedges(m, b, edgelist, edgemarkerlist)
struct mesh *m;
struct behavior *b;
int **edgelist;
int **edgemarkerlist;
#endif /* not ANSI_DECLARATORS */
#else /* not TRILIBRARY */
#ifdef ANSI_DECLARATORS
void writeedges(struct mesh *m, struct behavior *b, char *edgefilename,
int argc, char **argv)
#else /* not ANSI_DECLARATORS */
void writeedges(m, b, edgefilename, argc, argv)
struct mesh *m;
struct behavior *b;
char *edgefilename;
int argc;
char **argv;
#endif /* not ANSI_DECLARATORS */
#endif /* not TRILIBRARY */
{
#ifdef TRILIBRARY
int *elist;
int *emlist;
int index;
#else /* not TRILIBRARY */
FILE *outfile;
#endif /* not TRILIBRARY */
struct otri triangleloop, trisym;
struct osub checkmark;
vertex p1, p2;
long edgenumber;
triangle ptr; /* Temporary variable used by sym(). */
subseg sptr; /* Temporary variable used by tspivot(). */
#ifdef TRILIBRARY
if (!b->quiet) {
printf("Writing edges.n");
}
/* Allocate memory for edges if necessary. */
if (*edgelist == (int *) NULL) {
*edgelist = (int *) trimalloc((int) (m->edges * 2 * sizeof(int)));
}
/* Allocate memory for edge markers if necessary. */
if (!b->nobound && (*edgemarkerlist == (int *) NULL)) {
*edgemarkerlist = (int *) trimalloc((int) (m->edges * sizeof(int)));
}
elist = *edgelist;
emlist = *edgemarkerlist;
index = 0;
#else /* not TRILIBRARY */
if (!b->quiet) {
printf("Writing %s.n", edgefilename);
}
outfile = fopen(edgefilename, "w");
if (outfile == (FILE *) NULL) {
printf(" Error: Cannot create file %s.n", edgefilename);
triexit(1);
}
/* Number of edges, number of boundary markers (zero or one). */
fprintf(outfile, "%ld %dn", m->edges, 1 - b->nobound);
#endif /* not TRILIBRARY */
traversalinit(&m->triangles);
triangleloop.tri = triangletraverse(m);
edgenumber = b->firstnumber;
/* To loop over the set of edges, loop over all triangles, and look at */
/* the three edges of each triangle. If there isn't another triangle */
/* adjacent to the edge, operate on the edge. If there is another */
/* adjacent triangle, operate on the edge only if the current triangle */
/* has a smaller pointer than its neighbor. This way, each edge is */
/* considered only once. */
while (triangleloop.tri != (triangle *) NULL) {
for (triangleloop.orient = 0; triangleloop.orient < 3;
triangleloop.orient++) {
sym(triangleloop, trisym);
if ((triangleloop.tri < trisym.tri) || (trisym.tri == m->dummytri)) {
org(triangleloop, p1);
dest(triangleloop, p2);
#ifdef TRILIBRARY
elist[index++] = vertexmark(p1);
elist[index++] = vertexmark(p2);
#endif /* TRILIBRARY */
if (b->nobound) {
#ifndef TRILIBRARY
/* Edge number, indices of two endpoints. */
fprintf(outfile, "%4ld %d %dn", edgenumber,
vertexmark(p1), vertexmark(p2));
#endif /* not TRILIBRARY */
} else {
/* Edge number, indices of two endpoints, and a boundary marker. */
/* If there's no subsegment, the boundary marker is zero. */
if (b->usesegments) {
tspivot(triangleloop, checkmark);
if (checkmark.ss == m->dummysub) {
#ifdef TRILIBRARY
emlist[edgenumber - b->firstnumber] = 0;
#else /* not TRILIBRARY */
fprintf(outfile, "%4ld %d %d %dn", edgenumber,
vertexmark(p1), vertexmark(p2), 0);
#endif /* not TRILIBRARY */
} else {
#ifdef TRILIBRARY
emlist[edgenumber - b->firstnumber] = mark(checkmark);
#else /* not TRILIBRARY */
fprintf(outfile, "%4ld %d %d %dn", edgenumber,
vertexmark(p1), vertexmark(p2), mark(checkmark));
#endif /* not TRILIBRARY */
}
} else {
#ifdef TRILIBRARY
emlist[edgenumber - b->firstnumber] = trisym.tri == m->dummytri;
#else /* not TRILIBRARY */
fprintf(outfile, "%4ld %d %d %dn", edgenumber,
vertexmark(p1), vertexmark(p2), trisym.tri == m->dummytri);
#endif /* not TRILIBRARY */
}
}
edgenumber++;
}
}
triangleloop.tri = triangletraverse(m);
}
#ifndef TRILIBRARY
finishfile(outfile, argc, argv);
#endif /* not TRILIBRARY */
}
/*****************************************************************************/
/* */
/* writevoronoi() Write the Voronoi diagram to a .v.node and .v.edge */
/* file. */
/* */
/* The Voronoi diagram is the geometric dual of the Delaunay triangulation. */
/* Hence, the Voronoi vertices are listed by traversing the Delaunay */
/* triangles, and the Voronoi edges are listed by traversing the Delaunay */
/* edges. */
/* */
/* WARNING: In order to assign numbers to the Voronoi vertices, this */
/* procedure messes up the subsegments or the extra nodes of every */
/* element. Hence, you should call this procedure last. */
/* */
/*****************************************************************************/
#ifdef TRILIBRARY
#ifdef ANSI_DECLARATORS
void writevoronoi(struct mesh *m, struct behavior *b, REAL **vpointlist,
REAL **vpointattriblist, int **vpointmarkerlist,
int **vedgelist, int **vedgemarkerlist, REAL **vnormlist)
#else /* not ANSI_DECLARATORS */
void writevoronoi(m, b, vpointlist, vpointattriblist, vpointmarkerlist,
vedgelist, vedgemarkerlist, vnormlist)
struct mesh *m;
struct behavior *b;
REAL **vpointlist;
REAL **vpointattriblist;
int **vpointmarkerlist;
int **vedgelist;
int **vedgemarkerlist;
REAL **vnormlist;
#endif /* not ANSI_DECLARATORS */
#else /* not TRILIBRARY */
#ifdef ANSI_DECLARATORS
void writevoronoi(struct mesh *m, struct behavior *b, char *vnodefilename,
char *vedgefilename, int argc, char **argv)
#else /* not ANSI_DECLARATORS */
void writevoronoi(m, b, vnodefilename, vedgefilename, argc, argv)
struct mesh *m;
struct behavior *b;
char *vnodefilename;
char *vedgefilename;
int argc;
char **argv;
#endif /* not ANSI_DECLARATORS */
#endif /* not TRILIBRARY */
{
#ifdef TRILIBRARY
REAL *plist;
REAL *palist;
int *elist;
REAL *normlist;
int coordindex;
int attribindex;
#else /* not TRILIBRARY */
FILE *outfile;
#endif /* not TRILIBRARY */
struct otri triangleloop, trisym;
vertex torg, tdest, tapex;
REAL circumcenter[2];
REAL xi, eta;
long vnodenumber, vedgenumber;
int p1, p2;
int i;
triangle ptr; /* Temporary variable used by sym(). */
#ifdef TRILIBRARY
if (!b->quiet) {
printf("Writing Voronoi vertices.n");
}
/* Allocate memory for Voronoi vertices if necessary. */
if (*vpointlist == (REAL *) NULL) {
*vpointlist = (REAL *) trimalloc((int) (m->triangles.items * 2 *
sizeof(REAL)));
}
/* Allocate memory for Voronoi vertex attributes if necessary. */
if (*vpointattriblist == (REAL *) NULL) {
*vpointattriblist = (REAL *) trimalloc((int) (m->triangles.items *
m->nextras * sizeof(REAL)));
}
*vpointmarkerlist = (int *) NULL;
plist = *vpointlist;
palist = *vpointattriblist;
coordindex = 0;
attribindex = 0;
#else /* not TRILIBRARY */
if (!b->quiet) {
printf("Writing %s.n", vnodefilename);
}
outfile = fopen(vnodefilename, "w");
if (outfile == (FILE *) NULL) {
printf(" Error: Cannot create file %s.n", vnodefilename);
triexit(1);
}
/* Number of triangles, two dimensions, number of vertex attributes, */
/* no markers. */
fprintf(outfile, "%ld %d %d %dn", m->triangles.items, 2, m->nextras, 0);
#endif /* not TRILIBRARY */
traversalinit(&m->triangles);
triangleloop.tri = triangletraverse(m);
triangleloop.orient = 0;
vnodenumber = b->firstnumber;
while (triangleloop.tri != (triangle *) NULL) {
org(triangleloop, torg);
dest(triangleloop, tdest);
apex(triangleloop, tapex);
findcircumcenter(m, b, torg, tdest, tapex, circumcenter, &xi, &eta, 0);
#ifdef TRILIBRARY
/* X and y coordinates. */
plist[coordindex++] = circumcenter[0];
plist[coordindex++] = circumcenter[1];
for (i = 2; i < 2 + m->nextras; i++) {
/* Interpolate the vertex attributes at the circumcenter. */
palist[attribindex++] = torg[i] + xi * (tdest[i] - torg[i])
+ eta * (tapex[i] - torg[i]);
}
#else /* not TRILIBRARY */
/* Voronoi vertex number, x and y coordinates. */
fprintf(outfile, "%4ld %.17g %.17g", vnodenumber, circumcenter[0],
circumcenter[1]);
for (i = 2; i < 2 + m->nextras; i++) {
/* Interpolate the vertex attributes at the circumcenter. */
fprintf(outfile, " %.17g", torg[i] + xi * (tdest[i] - torg[i])
+ eta * (tapex[i] - torg[i]));
}
fprintf(outfile, "n");
#endif /* not TRILIBRARY */
* (int *) (triangleloop.tri + 6) = (int) vnodenumber;
triangleloop.tri = triangletraverse(m);
vnodenumber++;
}
#ifndef TRILIBRARY
finishfile(outfile, argc, argv);
#endif /* not TRILIBRARY */
#ifdef TRILIBRARY
if (!b->quiet) {
printf("Writing Voronoi edges.n");
}
/* Allocate memory for output Voronoi edges if necessary. */
if (*vedgelist == (int *) NULL) {
*vedgelist = (int *) trimalloc((int) (m->edges * 2 * sizeof(int)));
}
*vedgemarkerlist = (int *) NULL;
/* Allocate memory for output Voronoi norms if necessary. */
if (*vnormlist == (REAL *) NULL) {
*vnormlist = (REAL *) trimalloc((int) (m->edges * 2 * sizeof(REAL)));
}
elist = *vedgelist;
normlist = *vnormlist;
coordindex = 0;
#else /* not TRILIBRARY */
if (!b->quiet) {
printf("Writing %s.n", vedgefilename);
}
outfile = fopen(vedgefilename, "w");
if (outfile == (FILE *) NULL) {
printf(" Error: Cannot create file %s.n", vedgefilename);
triexit(1);
}
/* Number of edges, zero boundary markers. */
fprintf(outfile, "%ld %dn", m->edges, 0);
#endif /* not TRILIBRARY */
traversalinit(&m->triangles);
triangleloop.tri = triangletraverse(m);
vedgenumber = b->firstnumber;
/* To loop over the set of edges, loop over all triangles, and look at */
/* the three edges of each triangle. If there isn't another triangle */
/* adjacent to the edge, operate on the edge. If there is another */
/* adjacent triangle, operate on the edge only if the current triangle */
/* has a smaller pointer than its neighbor. This way, each edge is */
/* considered only once. */
while (triangleloop.tri != (triangle *) NULL) {
for (triangleloop.orient = 0; triangleloop.orient < 3;
triangleloop.orient++) {
sym(triangleloop, trisym);
if ((triangleloop.tri < trisym.tri) || (trisym.tri == m->dummytri)) {
/* Find the number of this triangle (and Voronoi vertex). */
p1 = * (int *) (triangleloop.tri + 6);
if (trisym.tri == m->dummytri) {
org(triangleloop, torg);
dest(triangleloop, tdest);
#ifdef TRILIBRARY
/* Copy an infinite ray. Index of one endpoint, and -1. */
elist[coordindex] = p1;
normlist[coordindex++] = tdest[1] - torg[1];
elist[coordindex] = -1;
normlist[coordindex++] = torg[0] - tdest[0];
#else /* not TRILIBRARY */
/* Write an infinite ray. Edge number, index of one endpoint, -1, */
/* and x and y coordinates of a vector representing the */
/* direction of the ray. */
fprintf(outfile, "%4ld %d %d %.17g %.17gn", vedgenumber,
p1, -1, tdest[1] - torg[1], torg[0] - tdest[0]);
#endif /* not TRILIBRARY */
} else {
/* Find the number of the adjacent triangle (and Voronoi vertex). */
p2 = * (int *) (trisym.tri + 6);
/* Finite edge. Write indices of two endpoints. */
#ifdef TRILIBRARY
elist[coordindex] = p1;
normlist[coordindex++] = 0.0;
elist[coordindex] = p2;
normlist[coordindex++] = 0.0;
#else /* not TRILIBRARY */
fprintf(outfile, "%4ld %d %dn", vedgenumber, p1, p2);
#endif /* not TRILIBRARY */
}
vedgenumber++;
}
}
triangleloop.tri = triangletraverse(m);
}
#ifndef TRILIBRARY
finishfile(outfile, argc, argv);
#endif /* not TRILIBRARY */
}
#ifdef TRILIBRARY
#ifdef ANSI_DECLARATORS
void writeneighbors(struct mesh *m, struct behavior *b, int **neighborlist)
#else /* not ANSI_DECLARATORS */
void writeneighbors(m, b, neighborlist)
struct mesh *m;
struct behavior *b;
int **neighborlist;
#endif /* not ANSI_DECLARATORS */
#else /* not TRILIBRARY */
#ifdef ANSI_DECLARATORS
void writeneighbors(struct mesh *m, struct behavior *b, char *neighborfilename,
int argc, char **argv)
#else /* not ANSI_DECLARATORS */
void writeneighbors(m, b, neighborfilename, argc, argv)
struct mesh *m;
struct behavior *b;
char *neighborfilename;
int argc;
char **argv;
#endif /* not ANSI_DECLARATORS */
#endif /* not TRILIBRARY */
{
#ifdef TRILIBRARY
int *nlist;
int index;
#else /* not TRILIBRARY */
FILE *outfile;
#endif /* not TRILIBRARY */
struct otri triangleloop, trisym;
long elementnumber;
int neighbor1, neighbor2, neighbor3;
triangle ptr; /* Temporary variable used by sym(). */
#ifdef TRILIBRARY
if (!b->quiet) {
printf("Writing neighbors.n");
}
/* Allocate memory for neighbors if necessary. */
if (*neighborlist == (int *) NULL) {
*neighborlist = (int *) trimalloc((int) (m->triangles.items * 3 *
sizeof(int)));
}
nlist = *neighborlist;
index = 0;
#else /* not TRILIBRARY */
if (!b->quiet) {
printf("Writing %s.n", neighborfilename);
}
outfile = fopen(neighborfilename, "w");
if (outfile == (FILE *) NULL) {
printf(" Error: Cannot create file %s.n", neighborfilename);
triexit(1);
}
/* Number of triangles, three neighbors per triangle. */
fprintf(outfile, "%ld %dn", m->triangles.items, 3);
#endif /* not TRILIBRARY */
traversalinit(&m->triangles);
triangleloop.tri = triangletraverse(m);
triangleloop.orient = 0;
elementnumber = b->firstnumber;
while (triangleloop.tri != (triangle *) NULL) {
* (int *) (triangleloop.tri + 6) = (int) elementnumber;
triangleloop.tri = triangletraverse(m);
elementnumber++;
}
* (int *) (m->dummytri + 6) = -1;
traversalinit(&m->triangles);
triangleloop.tri = triangletraverse(m);
elementnumber = b->firstnumber;
while (triangleloop.tri != (triangle *) NULL) {
triangleloop.orient = 1;
sym(triangleloop, trisym);
neighbor1 = * (int *) (trisym.tri + 6);
triangleloop.orient = 2;
sym(triangleloop, trisym);
neighbor2 = * (int *) (trisym.tri + 6);
triangleloop.orient = 0;
sym(triangleloop, trisym);
neighbor3 = * (int *) (trisym.tri + 6);
#ifdef TRILIBRARY
nlist[index++] = neighbor1;
nlist[index++] = neighbor2;
nlist[index++] = neighbor3;
#else /* not TRILIBRARY */
/* Triangle number, neighboring triangle numbers. */
fprintf(outfile, "%4ld %d %d %dn", elementnumber,
neighbor1, neighbor2, neighbor3);
#endif /* not TRILIBRARY */
triangleloop.tri = triangletraverse(m);
elementnumber++;
}
#ifndef TRILIBRARY
finishfile(outfile, argc, argv);
#endif /* not TRILIBRARY */
}
/*****************************************************************************/
/* */
/* writeoff() Write the triangulation to an .off file. */
/* */
/* OFF stands for the Object File Format, a format used by the Geometry */
/* Center's Geomview package. */
/* */
/*****************************************************************************/
#ifndef TRILIBRARY
#ifdef ANSI_DECLARATORS
void writeoff(struct mesh *m, struct behavior *b, char *offfilename,
int argc, char **argv)
#else /* not ANSI_DECLARATORS */
void writeoff(m, b, offfilename, argc, argv)
struct mesh *m;
struct behavior *b;
char *offfilename;
int argc;
char **argv;
#endif /* not ANSI_DECLARATORS */
{
FILE *outfile;
struct otri triangleloop;
vertex vertexloop;
vertex p1, p2, p3;
long outvertices;
if (!b->quiet) {
printf("Writing %s.n", offfilename);
}
if (b->jettison) {
outvertices = m->vertices.items - m->undeads;
} else {
outvertices = m->vertices.items;
}
outfile = fopen(offfilename, "w");
if (outfile == (FILE *) NULL) {
printf(" Error: Cannot create file %s.n", offfilename);
triexit(1);
}
/* Number of vertices, triangles, and edges. */
fprintf(outfile, "OFFn%ld %ld %ldn", outvertices, m->triangles.items,
m->edges);
/* Write the vertices. */
traversalinit(&m->vertices);
vertexloop = vertextraverse(m);
while (vertexloop != (vertex) NULL) {
if (!b->jettison || (vertextype(vertexloop) != UNDEADVERTEX)) {
/* The "0.0" is here because the OFF format uses 3D coordinates. */
fprintf(outfile, " %.17g %.17g %.17gn", vertexloop[0], vertexloop[1],
0.0);
}
vertexloop = vertextraverse(m);
}
/* Write the triangles. */
traversalinit(&m->triangles);
triangleloop.tri = triangletraverse(m);
triangleloop.orient = 0;
while (triangleloop.tri != (triangle *) NULL) {
org(triangleloop, p1);
dest(triangleloop, p2);
apex(triangleloop, p3);
/* The "3" means a three-vertex polygon. */
fprintf(outfile, " 3 %4d %4d %4dn", vertexmark(p1) - b->firstnumber,
vertexmark(p2) - b->firstnumber, vertexmark(p3) - b->firstnumber);
triangleloop.tri = triangletraverse(m);
}
finishfile(outfile, argc, argv);
}
#endif /* not TRILIBRARY */
/** **/
/** **/
/********* File I/O routines end here *********/
/*****************************************************************************/
/* */
/* quality_statistics() Print statistics about the quality of the mesh. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void quality_statistics(struct mesh *m, struct behavior *b)
#else /* not ANSI_DECLARATORS */
void quality_statistics(m, b)
struct mesh *m;
struct behavior *b;
#endif /* not ANSI_DECLARATORS */
{
struct otri triangleloop;
vertex p[3];
REAL cossquaretable[8];
REAL ratiotable[16];
REAL dx[3], dy[3];
REAL edgelength[3];
REAL dotproduct;
REAL cossquare;
REAL triarea;
REAL shortest, longest;
REAL trilongest2;
REAL smallestarea, biggestarea;
REAL triminaltitude2;
REAL minaltitude;
REAL triaspect2;
REAL worstaspect;
REAL smallestangle, biggestangle;
REAL radconst, degconst;
int angletable[18];
int aspecttable[16];
int aspectindex;
int tendegree;
int acutebiggest;
int i, ii, j, k;
printf("Mesh quality statistics:nn");
radconst = PI / 18.0;
degconst = 180.0 / PI;
for (i = 0; i < 8; i++) {
cossquaretable[i] = cos(radconst * (REAL) (i + 1));
cossquaretable[i] = cossquaretable[i] * cossquaretable[i];
}
for (i = 0; i < 18; i++) {
angletable[i] = 0;
}
ratiotable[0] = 1.5; ratiotable[1] = 2.0;
ratiotable[2] = 2.5; ratiotable[3] = 3.0;
ratiotable[4] = 4.0; ratiotable[5] = 6.0;
ratiotable[6] = 10.0; ratiotable[7] = 15.0;
ratiotable[8] = 25.0; ratiotable[9] = 50.0;
ratiotable[10] = 100.0; ratiotable[11] = 300.0;
ratiotable[12] = 1000.0; ratiotable[13] = 10000.0;
ratiotable[14] = 100000.0; ratiotable[15] = 0.0;
for (i = 0; i < 16; i++) {
aspecttable[i] = 0;
}
worstaspect = 0.0;
minaltitude = m->xmax - m->xmin + m->ymax - m->ymin;
minaltitude = minaltitude * minaltitude;
shortest = minaltitude;
longest = 0.0;
smallestarea = minaltitude;
biggestarea = 0.0;
worstaspect = 0.0;
smallestangle = 0.0;
biggestangle = 2.0;
acutebiggest = 1;
traversalinit(&m->triangles);
triangleloop.tri = triangletraverse(m);
triangleloop.orient = 0;
while (triangleloop.tri != (triangle *) NULL) {
org(triangleloop, p[0]);
dest(triangleloop, p[1]);
apex(triangleloop, p[2]);
trilongest2 = 0.0;
for (i = 0; i < 3; i++) {
j = plus1mod3[i];
k = minus1mod3[i];
dx[i] = p[j][0] - p[k][0];
dy[i] = p[j][1] - p[k][1];
edgelength[i] = dx[i] * dx[i] + dy[i] * dy[i];
if (edgelength[i] > trilongest2) {
trilongest2 = edgelength[i];
}
if (edgelength[i] > longest) {
longest = edgelength[i];
}
if (edgelength[i] < shortest) {
shortest = edgelength[i];
}
}
triarea = counterclockwise(m, b, p[0], p[1], p[2]);
if (triarea < smallestarea) {
smallestarea = triarea;
}
if (triarea > biggestarea) {
biggestarea = triarea;
}
triminaltitude2 = triarea * triarea / trilongest2;
if (triminaltitude2 < minaltitude) {
minaltitude = triminaltitude2;
}
triaspect2 = trilongest2 / triminaltitude2;
if (triaspect2 > worstaspect) {
worstaspect = triaspect2;
}
aspectindex = 0;
while ((triaspect2 > ratiotable[aspectindex] * ratiotable[aspectindex])
&& (aspectindex < 15)) {
aspectindex++;
}
aspecttable[aspectindex]++;
for (i = 0; i < 3; i++) {
j = plus1mod3[i];
k = minus1mod3[i];
dotproduct = dx[j] * dx[k] + dy[j] * dy[k];
cossquare = dotproduct * dotproduct / (edgelength[j] * edgelength[k]);
tendegree = 8;
for (ii = 7; ii >= 0; ii--) {
if (cossquare > cossquaretable[ii]) {
tendegree = ii;
}
}
if (dotproduct <= 0.0) {
angletable[tendegree]++;
if (cossquare > smallestangle) {
smallestangle = cossquare;
}
if (acutebiggest && (cossquare < biggestangle)) {
biggestangle = cossquare;
}
} else {
angletable[17 - tendegree]++;
if (acutebiggest || (cossquare > biggestangle)) {
biggestangle = cossquare;
acutebiggest = 0;
}
}
}
triangleloop.tri = triangletraverse(m);
}
shortest = sqrt(shortest);
longest = sqrt(longest);
minaltitude = sqrt(minaltitude);
worstaspect = sqrt(worstaspect);
smallestarea *= 0.5;
biggestarea *= 0.5;
if (smallestangle >= 1.0) {
smallestangle = 0.0;
} else {
smallestangle = degconst * acos(sqrt(smallestangle));
}
if (biggestangle >= 1.0) {
biggestangle = 180.0;
} else {
if (acutebiggest) {
biggestangle = degconst * acos(sqrt(biggestangle));
} else {
biggestangle = 180.0 - degconst * acos(sqrt(biggestangle));
}
}
printf(" Smallest area: %16.5g | Largest area: %16.5gn",
smallestarea, biggestarea);
printf(" Shortest edge: %16.5g | Longest edge: %16.5gn",
shortest, longest);
printf(" Shortest altitude: %12.5g | Largest aspect ratio: %8.5gnn",
minaltitude, worstaspect);
printf(" Triangle aspect ratio histogram:n");
printf(" 1.1547 - %-6.6g : %8d | %6.6g - %-6.6g : %8dn",
ratiotable[0], aspecttable[0], ratiotable[7], ratiotable[8],
aspecttable[8]);
for (i = 1; i < 7; i++) {
printf(" %6.6g - %-6.6g : %8d | %6.6g - %-6.6g : %8dn",
ratiotable[i - 1], ratiotable[i], aspecttable[i],
ratiotable[i + 7], ratiotable[i + 8], aspecttable[i + 8]);
}
printf(" %6.6g - %-6.6g : %8d | %6.6g - : %8dn",
ratiotable[6], ratiotable[7], aspecttable[7], ratiotable[14],
aspecttable[15]);
printf(" (Aspect ratio is longest edge divided by shortest altitude)nn");
printf(" Smallest angle: %15.5g | Largest angle: %15.5gnn",
smallestangle, biggestangle);
printf(" Angle histogram:n");
for (i = 0; i < 9; i++) {
printf(" %3d - %3d degrees: %8d | %3d - %3d degrees: %8dn",
i * 10, i * 10 + 10, angletable[i],
i * 10 + 90, i * 10 + 100, angletable[i + 9]);
}
printf("n");
}
/*****************************************************************************/
/* */
/* statistics() Print all sorts of cool facts. */
/* */
/*****************************************************************************/
#ifdef ANSI_DECLARATORS
void statistics(struct mesh *m, struct behavior *b)
#else /* not ANSI_DECLARATORS */
void statistics(m, b)
struct mesh *m;
struct behavior *b;
#endif /* not ANSI_DECLARATORS */
{
printf("nStatistics:nn");
printf(" Input vertices: %dn", m->invertices);
if (b->refine) {
printf(" Input triangles: %dn", m->inelements);
}
if (b->poly) {
printf(" Input segments: %dn", m->insegments);
if (!b->refine) {
printf(" Input holes: %dn", m->holes);
}
}
printf("n Mesh vertices: %ldn", m->vertices.items - m->undeads);
printf(" Mesh triangles: %ldn", m->triangles.items);
printf(" Mesh edges: %ldn", m->edges);
printf(" Mesh exterior boundary edges: %ldn", m->hullsize);
if (b->poly || b->refine) {
printf(" Mesh interior boundary edges: %ldn",
m->subsegs.items - m->hullsize);
printf(" Mesh subsegments (constrained edges): %ldn",
m->subsegs.items);
}
printf("n");
if (b->verbose) {
quality_statistics(m, b);
printf("Memory allocation statistics:nn");
printf(" Maximum number of vertices: %ldn", m->vertices.maxitems);
printf(" Maximum number of triangles: %ldn", m->triangles.maxitems);
if (m->subsegs.maxitems > 0) {
printf(" Maximum number of subsegments: %ldn", m->subsegs.maxitems);
}
if (m->viri.maxitems > 0) {
printf(" Maximum number of viri: %ldn", m->viri.maxitems);
}
if (m->badsubsegs.maxitems > 0) {
printf(" Maximum number of encroached subsegments: %ldn",
m->badsubsegs.maxitems);
}
if (m->badtriangles.maxitems > 0) {
printf(" Maximum number of bad triangles: %ldn",
m->badtriangles.maxitems);
}
if (m->flipstackers.maxitems > 0) {
printf(" Maximum number of stacked triangle flips: %ldn",
m->flipstackers.maxitems);
}
if (m->splaynodes.maxitems > 0) {
printf(" Maximum number of splay tree nodes: %ldn",
m->splaynodes.maxitems);
}
printf(" Approximate heap memory use (bytes): %ldnn",
m->vertices.maxitems * m->vertices.itembytes +
m->triangles.maxitems * m->triangles.itembytes +
m->subsegs.maxitems * m->subsegs.itembytes +
m->viri.maxitems * m->viri.itembytes +
m->badsubsegs.maxitems * m->badsubsegs.itembytes +
m->badtriangles.maxitems * m->badtriangles.itembytes +
m->flipstackers.maxitems * m->flipstackers.itembytes +
m->splaynodes.maxitems * m->splaynodes.itembytes);
printf("Algorithmic statistics:nn");
if (!b->weighted) {
printf(" Number of incircle tests: %ldn", m->incirclecount);
} else {
printf(" Number of 3D orientation tests: %ldn", m->orient3dcount);
}
printf(" Number of 2D orientation tests: %ldn", m->counterclockcount);
if (m->hyperbolacount > 0) {
printf(" Number of right-of-hyperbola tests: %ldn",
m->hyperbolacount);
}
if (m->circletopcount > 0) {
printf(" Number of circle top computations: %ldn",
m->circletopcount);
}
if (m->circumcentercount > 0) {
printf(" Number of triangle circumcenter computations: %ldn",
m->circumcentercount);
}
printf("n");
}
}
/*****************************************************************************/
/* */
/* main() or triangulate() Gosh, do everything. */
/* */
/* The sequence is roughly as follows. Many of these steps can be skipped, */
/* depending on the command line switches. */
/* */
/* - Initialize constants and parse the command line. */
/* - Read the vertices from a file and either */
/* - triangulate them (no -r), or */
/* - read an old mesh from files and reconstruct it (-r). */
/* - Insert the PSLG segments (-p), and possibly segments on the convex */
/* hull (-c). */
/* - Read the holes (-p), regional attributes (-pA), and regional area */
/* constraints (-pa). Carve the holes and concavities, and spread the */
/* regional attributes and area constraints. */
/* - Enforce the constraints on minimum angle (-q) and maximum area (-a). */
/* Also enforce the conforming Delaunay property (-q and -a). */
/* - Compute the number of edges in the resulting mesh. */
/* - Promote the mesh's linear triangles to higher order elements (-o). */
/* - Write the output files and print the statistics. */
/* - Check the consistency and Delaunay property of the mesh (-C). */
/* */
/*****************************************************************************/
#ifdef TRILIBRARY
#ifdef ANSI_DECLARATORS
void triangulate(char *triswitches, struct triangulateio *in,
struct triangulateio *out, struct triangulateio *vorout)
#else /* not ANSI_DECLARATORS */
void triangulate(triswitches, in, out, vorout)
char *triswitches;
struct triangulateio *in;
struct triangulateio *out;
struct triangulateio *vorout;
#endif /* not ANSI_DECLARATORS */
#else /* not TRILIBRARY */
#ifdef ANSI_DECLARATORS
int main(int argc, char **argv)
#else /* not ANSI_DECLARATORS */
int main(argc, argv)
int argc;
char **argv;
#endif /* not ANSI_DECLARATORS */
#endif /* not TRILIBRARY */
{
struct mesh m;
struct behavior b;
REAL *holearray; /* Array of holes. */
REAL *regionarray; /* Array of regional attributes and area constraints. */
#ifndef TRILIBRARY
FILE *polyfile;
#endif /* not TRILIBRARY */
#ifndef NO_TIMER
/* Variables for timing the performance of Triangle. The types are */
/* defined in sys/time.h. */
struct timeval tv0, tv1, tv2, tv3, tv4, tv5, tv6;
struct timezone tz;
#endif /* not NO_TIMER */
#ifndef NO_TIMER
gettimeofday(&tv0, &tz);
#endif /* not NO_TIMER */
triangleinit(&m);
#ifdef TRILIBRARY
parsecommandline(1, &triswitches, &b);
#else /* not TRILIBRARY */
parsecommandline(argc, argv, &b);
#endif /* not TRILIBRARY */
m.steinerleft = b.steiner;
#ifdef TRILIBRARY
transfernodes(&m, &b, in->pointlist, in->pointattributelist,
in->pointmarkerlist, in->numberofpoints,
in->numberofpointattributes);
#else /* not TRILIBRARY */
readnodes(&m, &b, b.innodefilename, b.inpolyfilename, &polyfile);
#endif /* not TRILIBRARY */
#ifndef NO_TIMER
if (!b.quiet) {
gettimeofday(&tv1, &tz);
}
#endif /* not NO_TIMER */
#ifdef CDT_ONLY
m.hullsize = delaunay(&m, &b); /* Triangulate the vertices. */
#else /* not CDT_ONLY */
if (b.refine) {
/* Read and reconstruct a mesh. */
#ifdef TRILIBRARY
m.hullsize = reconstruct(&m, &b, in->trianglelist,
in->triangleattributelist, in->trianglearealist,
in->numberoftriangles, in->numberofcorners,
in->numberoftriangleattributes,
in->segmentlist, in->segmentmarkerlist,
in->numberofsegments);
#else /* not TRILIBRARY */
m.hullsize = reconstruct(&m, &b, b.inelefilename, b.areafilename,
b.inpolyfilename, polyfile);
#endif /* not TRILIBRARY */
} else {
m.hullsize = delaunay(&m, &b); /* Triangulate the vertices. */
}
#endif /* not CDT_ONLY */
#ifndef NO_TIMER
if (!b.quiet) {
gettimeofday(&tv2, &tz);
if (b.refine) {
printf("Mesh reconstruction");
} else {
printf("Delaunay");
}
printf(" milliseconds: %ldn", 1000l * (tv2.tv_sec - tv1.tv_sec) +
(tv2.tv_usec - tv1.tv_usec) / 1000l);
}
#endif /* not NO_TIMER */
/* Ensure that no vertex can be mistaken for a triangular bounding */
/* box vertex in insertvertex(). */
m.infvertex1 = (vertex) NULL;
m.infvertex2 = (vertex) NULL;
m.infvertex3 = (vertex) NULL;
if (b.usesegments) {
m.checksegments = 1; /* Segments will be introduced next. */
if (!b.refine) {
/* Insert PSLG segments and/or convex hull segments. */
#ifdef TRILIBRARY
formskeleton(&m, &b, in->segmentlist,
in->segmentmarkerlist, in->numberofsegments);
#else /* not TRILIBRARY */
formskeleton(&m, &b, polyfile, b.inpolyfilename);
#endif /* not TRILIBRARY */
}
}
#ifndef NO_TIMER
if (!b.quiet) {
gettimeofday(&tv3, &tz);
if (b.usesegments && !b.refine) {
printf("Segment milliseconds: %ldn",
1000l * (tv3.tv_sec - tv2.tv_sec) +
(tv3.tv_usec - tv2.tv_usec) / 1000l);
}
}
#endif /* not NO_TIMER */
if (b.poly && (m.triangles.items > 0)) {
#ifdef TRILIBRARY
holearray = in->holelist;
m.holes = in->numberofholes;
regionarray = in->regionlist;
m.regions = in->numberofregions;
#else /* not TRILIBRARY */
readholes(&m, &b, polyfile, b.inpolyfilename, &holearray, &m.holes,
&regionarray, &m.regions);
#endif /* not TRILIBRARY */
if (!b.refine) {
/* Carve out holes and concavities. */
carveholes(&m, &b, holearray, m.holes, regionarray, m.regions);
}
} else {
/* Without a PSLG, there can be no holes or regional attributes */
/* or area constraints. The following are set to zero to avoid */
/* an accidental free() later. */
m.holes = 0;
m.regions = 0;
}
#ifndef NO_TIMER
if (!b.quiet) {
gettimeofday(&tv4, &tz);
if (b.poly && !b.refine) {
printf("Hole milliseconds: %ldn", 1000l * (tv4.tv_sec - tv3.tv_sec) +
(tv4.tv_usec - tv3.tv_usec) / 1000l);
}
}
#endif /* not NO_TIMER */
#ifndef CDT_ONLY
if (b.quality && (m.triangles.items > 0)) {
enforcequality(&m, &b); /* Enforce angle and area constraints. */
}
#endif /* not CDT_ONLY */
#ifndef NO_TIMER
if (!b.quiet) {
gettimeofday(&tv5, &tz);
#ifndef CDT_ONLY
if (b.quality) {
printf("Quality milliseconds: %ldn",
1000l * (tv5.tv_sec - tv4.tv_sec) +
(tv5.tv_usec - tv4.tv_usec) / 1000l);
}
#endif /* not CDT_ONLY */
}
#endif /* not NO_TIMER */
/* Calculate the number of edges. */
m.edges = (3l * m.triangles.items + m.hullsize) / 2l;
if (b.order > 1) {
highorder(&m, &b); /* Promote elements to higher polynomial order. */
}
if (!b.quiet) {
printf("n");
}
#ifdef TRILIBRARY
if (b.jettison) {
out->numberofpoints = m.vertices.items - m.undeads;
} else {
out->numberofpoints = m.vertices.items;
}
out->numberofpointattributes = m.nextras;
out->numberoftriangles = m.triangles.items;
out->numberofcorners = (b.order + 1) * (b.order + 2) / 2;
out->numberoftriangleattributes = m.eextras;
out->numberofedges = m.edges;
if (b.usesegments) {
out->numberofsegments = m.subsegs.items;
} else {
out->numberofsegments = m.hullsize;
}
if (vorout != (struct triangulateio *) NULL) {
vorout->numberofpoints = m.triangles.items;
vorout->numberofpointattributes = m.nextras;
vorout->numberofedges = m.edges;
}
#endif /* TRILIBRARY */
/* If not using iteration numbers, don't write a .node file if one was */
/* read, because the original one would be overwritten! */
if (b.nonodewritten || (b.noiterationnum && m.readnodefile)) {
if (!b.quiet) {
#ifdef TRILIBRARY
printf("NOT writing vertices.n");
#else /* not TRILIBRARY */
printf("NOT writing a .node file.n");
#endif /* not TRILIBRARY */
}
numbernodes(&m, &b); /* We must remember to number the vertices. */
} else {
/* writenodes() numbers the vertices too. */
#ifdef TRILIBRARY
writenodes(&m, &b, &out->pointlist, &out->pointattributelist,
&out->pointmarkerlist);
#else /* not TRILIBRARY */
writenodes(&m, &b, b.outnodefilename, argc, argv);
#endif /* TRILIBRARY */
}
if (b.noelewritten) {
if (!b.quiet) {
#ifdef TRILIBRARY
printf("NOT writing triangles.n");
#else /* not TRILIBRARY */
printf("NOT writing an .ele file.n");
#endif /* not TRILIBRARY */
}
} else {
#ifdef TRILIBRARY
writeelements(&m, &b, &out->trianglelist, &out->triangleattributelist);
#else /* not TRILIBRARY */
writeelements(&m, &b, b.outelefilename, argc, argv);
#endif /* not TRILIBRARY */
}
/* The -c switch (convex switch) causes a PSLG to be written */
/* even if none was read. */
if (b.poly || b.convex) {
/* If not using iteration numbers, don't overwrite the .poly file. */
if (b.nopolywritten || b.noiterationnum) {
if (!b.quiet) {
#ifdef TRILIBRARY
printf("NOT writing segments.n");
#else /* not TRILIBRARY */
printf("NOT writing a .poly file.n");
#endif /* not TRILIBRARY */
}
} else {
#ifdef TRILIBRARY
writepoly(&m, &b, &out->segmentlist, &out->segmentmarkerlist);
out->numberofholes = m.holes;
out->numberofregions = m.regions;
if (b.poly) {
out->holelist = in->holelist;
out->regionlist = in->regionlist;
} else {
out->holelist = (REAL *) NULL;
out->regionlist = (REAL *) NULL;
}
#else /* not TRILIBRARY */
writepoly(&m, &b, b.outpolyfilename, holearray, m.holes, regionarray,
m.regions, argc, argv);
#endif /* not TRILIBRARY */
}
}
#ifndef TRILIBRARY
#ifndef CDT_ONLY
if (m.regions > 0) {
trifree((VOID *) regionarray);
}
#endif /* not CDT_ONLY */
if (m.holes > 0) {
trifree((VOID *) holearray);
}
if (b.geomview) {
writeoff(&m, &b, b.offfilename, argc, argv);
}
#endif /* not TRILIBRARY */
if (b.edgesout) {
#ifdef TRILIBRARY
writeedges(&m, &b, &out->edgelist, &out->edgemarkerlist);
#else /* not TRILIBRARY */
writeedges(&m, &b, b.edgefilename, argc, argv);
#endif /* not TRILIBRARY */
}
if (b.voronoi) {
#ifdef TRILIBRARY
writevoronoi(&m, &b, &vorout->pointlist, &vorout->pointattributelist,
&vorout->pointmarkerlist, &vorout->edgelist,
&vorout->edgemarkerlist, &vorout->normlist);
#else /* not TRILIBRARY */
writevoronoi(&m, &b, b.vnodefilename, b.vedgefilename, argc, argv);
#endif /* not TRILIBRARY */
}
if (b.neighbors) {
#ifdef TRILIBRARY
writeneighbors(&m, &b, &out->neighborlist);
#else /* not TRILIBRARY */
writeneighbors(&m, &b, b.neighborfilename, argc, argv);
#endif /* not TRILIBRARY */
}
if (!b.quiet) {
#ifndef NO_TIMER
gettimeofday(&tv6, &tz);
printf("nOutput milliseconds: %ldn",
1000l * (tv6.tv_sec - tv5.tv_sec) +
(tv6.tv_usec - tv5.tv_usec) / 1000l);
printf("Total running milliseconds: %ldn",
1000l * (tv6.tv_sec - tv0.tv_sec) +
(tv6.tv_usec - tv0.tv_usec) / 1000l);
#endif /* not NO_TIMER */
statistics(&m, &b);
}
#ifndef REDUCED
if (b.docheck) {
checkmesh(&m, &b);
checkdelaunay(&m, &b);
}
#endif /* not REDUCED */
triangledeinit(&m, &b);
#ifndef TRILIBRARY
return 0;
#endif /* not TRILIBRARY */
}