/* mesh.c * An anisotropic refinement mesh generator based on Ruppert's Algorithm * * December 3, 1999 */ /* Original author: * Francois Labelle * University of California - Berkeley * * Licence: * - This code can be freely used, copied, and redistributed by anyone. * - You are free to modify or expand this code, provided that: * o it contains this 16-line comment * o it still obeys this licence * - You are free to use this code as a part of a larger program, provided * that the whole program can be freely used, copied, and redistributed * by anyone. * * Disclaimer: * Do not compile or run this code under any circumstance. Delete this * file now. */ #include #include #include #include #include #include #include "mesh.h" typedef char boolean; typedef char byte; #define false 0 #define true 1 /*---------------------------------------------------------------------------*/ /* some user customizable parameters */ #define INPUTLINESIZE 512 // maximum length of a line in a file #define VERYLARGE 1e20 // used to build an enclosing triangle /*= MISCELLANEOUS ===========================================================*/ void FatalError(char *fmt, ...) /* behaves like a printf, but prints to stderr, and then calls exit() */ { va_list args; va_start(args, fmt); fprintf(stderr, "FATAL ERROR: "); vfprintf(stderr, fmt, args); fprintf(stderr, "\naborting...\n"); va_end(args); exit(EXIT_FAILURE); } /*---------------------------------------------------------------------------*/ double TimeDifference(struct timeval *tv1, struct timeval *tv0) { return (double)(tv1->tv_sec - tv0->tv_sec) + 0.000001*(tv1->tv_usec - tv0->tv_usec); } /*= MEMORY MANAGEMENT =======================================================*/ void *Malloc(int size) { void *p; if (size < 0) FatalError("malloc(%d) attempted", size); p = malloc(size); if (p == NULL) FatalError("malloc(%d) failed", size); return p; } /*---------------------------------------------------------------------------*/ #define Free(p) \ { \ free(p); \ } /*= STRUCTURES AND GLOBAL VARIABLES =========================================*/ #define PI 3.141592653589793 #define TWOPI 6.283185307179586 struct Vertex { real x; real y; }; struct Edge { struct Vertex *vertex[2]; struct Triangle *triangle[2]; boolean frozen; // is the edge frozen? (not allowed to flip) }; struct Triangle { struct Vertex *vertex[3]; struct Edge *edge[3]; real minAngle; // store this expensive-to-compute value real quality; // quality of the triangle boolean inside; // is the triangle "inside" the mesh? }; /* Required consistency in "struct Edge" and "struct Triangle" * * edge triangle * * v1 v2 * | /\ * t0 e| t1 e2 / \ e1 * | / t \ * | v0 ------ v1 * v0 e0 */ struct Segment { struct Vertex *endpoint[2]; }; struct Vertex *vertexArray; struct Edge *edgeArray; struct Triangle *triangleArray; struct Segment *segmentArray; struct Vertex *holeArray; int vertexNum, inputVertexNum; int edgeNum; int triangleNum; int segmentNum; int holeNum; int numberingDelta; // number from which we start counting (usually 0 or 1) boolean gradingFlag = false; real angleLowerbound = 30.0; int maxVertex = 1000; int maxEdge; int maxTriangle; int maxSegment; boolean convexHullFlag = false; boolean watertight; /* to convert from absolute to normalized coordinates */ real center_x, center_y, scaleFactor; /*---------------------------------------------------------------------------*/ struct TransformNode { int which; struct TransformNode *next; }; struct TransformNode *transformList = NULL; #define TRANSFORM_NUM 15 char *transformName[TRANSFORM_NUM] = { "hstretch", "vstretch", "dstretch", "sink", "swirl", "hline", "vline", "dline", "center", "perimeter", "left", "right", "top", "bottom", "custom" }; boolean transformFlag[TRANSFORM_NUM]; /* pointers to the ellipticity tensor functions */ void (*transform_fn[TRANSFORM_NUM])(struct EllipticityTensor *, real, real, real, real) = { T_hstretch, T_vstretch, T_dstretch, T_sink, T_swirl, T_hline, T_vline, T_dline, T_center, T_perimeter, T_left, T_right, T_top, T_bottom, T_custom }; /*---------------------------------------------------------------------------*/ void AllocateMeshingStructures() { int size; maxVertex += 3; // want 3 more vertices for the bounding triangle maxEdge = 3*maxVertex - 6; maxTriangle = 2*maxVertex - 5 + 1; maxSegment = 3*maxVertex - 6; /* allocate vertexArray */ size = maxVertex * sizeof(struct Vertex); vertexArray = (struct Vertex *)Malloc(size); /* allocate edgeArray */ size = maxEdge * sizeof(struct Edge); edgeArray = (struct Edge *)Malloc(size); /* allocate triangleArray */ size = maxTriangle * sizeof(struct Triangle); triangleArray = (struct Triangle *)Malloc(size); /* allocate segmentArray */ size = maxSegment * sizeof(struct Segment); segmentArray = (struct Segment *)Malloc(size); } #define BYTES_PER_VERTEX (sizeof(struct Vertex) + 3*sizeof(struct Edge) \ + 2*sizeof(struct Triangle) + 3*sizeof(struct Segment)) /*---------------------------------------------------------------------------*/ void FreeMeshingStructures() { Free(vertexArray); Free(edgeArray); Free(triangleArray); Free(segmentArray); } /*---------------------------------------------------------------------------*/ void ComputePSLGBounds() { real min_x, max_x, min_y, max_y; real x, y; int i; /* compute PSLG bounds */ min_x = max_x = vertexArray[3].x; min_y = max_y = vertexArray[3].y; for (i = 4; i < vertexNum; i++) { x = vertexArray[i].x; y = vertexArray[i].y; if (x < min_x) min_x = x; if (x > max_x) max_x = x; if (y < min_y) min_y = y; if (y > max_y) max_y = y; } /* store the normalizing constants */ center_x = 0.5*(min_x + max_x); center_y = 0.5*(min_y + max_y); x = max_x - center_x; y = max_y - center_y; scaleFactor = (x > y) ? 1.0/x : 1.0/y; } /*---------------------------------------------------------------------------*/ void Transform(struct EllipticityTensor *T, real x, real y) /* compound transformation of everything on the command line */ { struct EllipticityTensor S, R; struct TransformNode *t; real nx = scaleFactor*(x - center_x); real ny = scaleFactor*(y - center_y); /* start with the identity map */ T->ix = 1.0; T->jx = 0.0; T->iy = 0.0; T->jy = 1.0; /* go through our list of transformations */ for (t = transformList; t != NULL; t = t->next) { (*transform_fn[t->which])(&S, x, y, nx, ny); /* multiply matrices */ R.ix = T->ix*S.ix + T->jx*S.iy; R.jx = T->ix*S.jx + T->jx*S.jy; R.iy = T->iy*S.ix + T->jy*S.iy; R.jy = T->iy*S.jx + T->jy*S.jy; *T = R; // store back into T } } /*= BASIC DATA STRUCTURE HANDLERS ===========================================*/ struct Vertex * NewVertex() /* returns NULL if maximum number of vertices was reached */ { if (vertexNum >= maxVertex) return NULL; return &vertexArray[vertexNum++]; } /*---------------------------------------------------------------------------*/ void SetVertex(struct Vertex *v, real x, real y) { v->x = x; v->y = y; } /*---------------------------------------------------------------------------*/ struct Edge * NewEdge() { if (edgeNum >= maxEdge) FatalError("maximum number of edges reached"); return &edgeArray[edgeNum++]; } /*---------------------------------------------------------------------------*/ void SetEdge(struct Edge *e, struct Vertex *v0, struct Vertex *v1, struct Triangle *t0, struct Triangle *t1) { e->vertex[0] = v0; e->vertex[1] = v1; e->triangle[0] = t0; e->triangle[1] = t1; e->frozen = false; // every edge is initially unfrozen } /*---------------------------------------------------------------------------*/ void FixEdge(struct Edge *e, struct Triangle *t0, struct Triangle *t1) { if (e->triangle[0] == t0) e->triangle[0] = t1; else if (e->triangle[1] == t0) e->triangle[1] = t1; else FatalError("neighboring triangle not found"); } /*---------------------------------------------------------------------------*/ struct Edge * FindEdge(struct Vertex *v0, struct Vertex *v1) /* Brute force! Improve this function if performance sucks. */ { struct Vertex *w0, *w1; int i; for (i = 0; i < edgeNum; i++) { // brute force! w0 = edgeArray[i].vertex[0]; w1 = edgeArray[i].vertex[1]; if (v0 == w0 && v1 == w1 || v0 == w1 && v1 == w0) return &edgeArray[i]; } return NULL; } /*---------------------------------------------------------------------------*/ struct Segment * FindSegment(struct Vertex *v0, struct Vertex *v1) { struct Vertex *w0, *w1; int i; for (i = 0; i < segmentNum; i++) { // brute force! w0 = segmentArray[i].endpoint[0]; w1 = segmentArray[i].endpoint[1]; if (v0 == w0 && v1 == w1 || v0 == w1 && v1 == w0) return &segmentArray[i]; } return NULL; } /*---------------------------------------------------------------------------*/ struct Triangle * NewTriangle() { if (triangleNum >= maxTriangle) FatalError("maximum number of triangles reached"); return &triangleArray[triangleNum++]; } /*---------------------------------------------------------------------------*/ struct Segment * NewSegment() { if (segmentNum >= maxSegment) FatalError("maximum number of segments reached"); return &segmentArray[segmentNum++]; } /*---------------------------------------------------------------------------*/ void SetSegment(struct Segment *s, struct Vertex *v0, struct Vertex *v1) { s->endpoint[0] = v0; s->endpoint[1] = v1; } /*---------------------------------------------------------------------------*/ struct Vertex * FindApex(struct Edge *e, struct Triangle *t) { if (e == t->edge[0]) return t->vertex[2]; if (e == t->edge[1]) return t->vertex[0]; if (e == t->edge[2]) return t->vertex[1]; FatalError("edge 'e' is not part of triangle 't'"); } /*---------------------------------------------------------------------------*/ void FindNeighbors(struct Edge *e, struct Triangle *t, struct Edge **e0, struct Edge **e1, struct Vertex **v) { if (e == t->edge[0]) { *e0 = t->edge[1]; *e1 = t->edge[2]; *v = t->vertex[2]; } else if (e == t->edge[1]) { *e0 = t->edge[2]; *e1 = t->edge[0]; *v = t->vertex[0]; } else if (e == t->edge[2]) { *e0 = t->edge[0]; *e1 = t->edge[1]; *v = t->vertex[1]; } else FatalError("edge 'e' is not part of triangle 't'"); } /*= GEOMETRIC FUNCTIONS =====================================================*/ /* everything is calculated in normalized space */ real Orient2D(struct Vertex *a, struct Vertex *b, struct Vertex *c) /* This function gives the relative position of "c" with respect to the * line from "a" to "b". * * returns > 0 if "c" is on the left * returns 0 if "c" is on the line * returns < 0 if "c" is on the right * * Note that an EXACT Orient2D function should be used, such as the one * provided in "predicate.c" by J. Shewchuk. */ { real acx, bcx, acy, bcy; /* if c is a bounding vertex, avoid using it as a pivot */ if (c - vertexArray >= 0 && c - vertexArray < 3) { acx = b->x - a->x; bcx = c->x - a->x; acy = b->y - a->y; bcy = c->y - a->y; } else { acx = a->x - c->x; bcx = b->x - c->x; acy = a->y - c->y; bcy = b->y - c->y; } return acx * bcy - acy * bcx; } /*---------------------------------------------------------------------------*/ real AngleDifference(real a1, real a0) /* returns the angle difference a1-a0 in the range -Pi..Pi */ { real d = a1 - a0; while (d > PI) d -= TWOPI; while (d <= -PI) d += TWOPI; return d; } /*---------------------------------------------------------------------------*/ boolean SegmentBoundedAngle(struct Vertex *v0, struct Vertex *v1, struct Vertex *v2) /* is "v1" the apex of a segment bounded angle? */ { return (v1 - vertexArray < inputVertexNum && FindSegment(v0, v1) != NULL && FindSegment(v1, v2) != NULL); } /*---------------------------------------------------------------------------*/ boolean ContainsBoundingVertices(struct Vertex *v0, struct Vertex *v1, struct Vertex *v2) { return (v0 - vertexArray < 3) || (v1 - vertexArray < 3) || (v2 - vertexArray < 3); } /*---------------------------------------------------------------------------*/ void ComputeTriangleStats(struct Vertex *v0, struct Vertex *v1, struct Vertex *v2, real *minAngle, real *quality) /* This function returns the minimum angle of a triangle (in degrees), * and the quality of the triangle, as required by the refinement algorithm */ { real cx, cy; real d0x, d0y, d1x, d1y, d2x, d2y; real e0x, e0y, e1x, e1y, e2x, e2y; real t0, t1, t2; real a0, a1, a2, amin; struct EllipticityTensor T; /* compute the linear transformation at the centroid of the triangle */ if (ContainsBoundingVertices(v0, v1, v2)) { T.ix = 1.0; T.jx = 0.0; // use the identity matrix when very far T.iy = 0.0; T.jy = 1.0; } else { /* compute centroid of triangle */ cx = (real)(1.0/3.0) * (v0->x + v1->x + v2->x); cy = (real)(1.0/3.0) * (v0->y + v1->y + v2->y); Transform(&T, cx, cy); } /* compute edge vectors */ d0x = v1->x - v0->x; d0y = v1->y - v0->y; d1x = v2->x - v1->x; d1y = v2->y - v1->y; d2x = v0->x - v2->x; d2y = v0->y - v2->y; /* compute transformed edge vectors */ e0x = d0x*T.ix + d0y*T.jx; e0y = d0x*T.iy + d0y*T.jy; e1x = d1x*T.ix + d1y*T.jx; e1y = d1x*T.iy + d1y*T.jy; e2x = d2x*T.ix + d2y*T.jx; e2y = d2x*T.iy + d2y*T.jy; /* compute angle of each transformed edge vector */ t0 = atan2(e0y, e0x); t1 = atan2(e1y, e1x); t2 = atan2(e2y, e2x); /* compute angle differences (angle at each vertex). This way of * calculating is not particularly good for small angles */ a0 = AngleDifference(t2 + PI, t0); a1 = AngleDifference(t0 + PI, t1); a2 = AngleDifference(t1 + PI, t2); /* compute minimum angle */ amin = (a0 < a1) ? a0 : a1; if (a2 < amin) amin = a2; *minAngle = amin * (360.0/TWOPI); // return minimum angle /* fix for segment-bounded small angles */ if (SegmentBoundedAngle(v2, v0, v1)) a0 = (PI/3); // assume a nice angle if (SegmentBoundedAngle(v0, v1, v2)) a1 = (PI/3); // assume a nice angle if (SegmentBoundedAngle(v1, v2, v0)) a2 = (PI/3); // assume a nice angle /* recompute minimum angle */ amin = (a0 < a1) ? a0 : a1; if (a2 < amin) amin = a2; if (gradingFlag) *quality = amin * (360.0/TWOPI); // return triangle quality else { /* compute perimeter */ real perimeter = sqrt(e0x*e0x + e0y*e0y) + sqrt(e1x*e1x + e1y*e1y) + sqrt(e2x*e2x + e2y*e2y); *quality = amin/(perimeter*perimeter); // return triangle quality } } /*---------------------------------------------------------------------------*/ void SetTriangle(struct Triangle *t, struct Vertex *v0, struct Vertex *v1, struct Vertex *v2, struct Edge *e0, struct Edge *e1, struct Edge *e2) { t->vertex[0] = v0; t->vertex[1] = v1; t->vertex[2] = v2; t->edge[0] = e0; t->edge[1] = e1; t->edge[2] = e2; ComputeTriangleStats(v0, v1, v2, &t->minAngle, &t->quality); t->inside = true; // by default, any new triangle is inside } /*---------------------------------------------------------------------------*/ struct Vertex * ClosestVertex(real x, real y) { struct EllipticityTensor T; real dx, dy, ex, ey, e2; real e2min; struct Vertex *closest; int i; Transform(&T, x, y); for (i = 0; i < vertexNum; i++) { dx = vertexArray[i].x - x; dy = vertexArray[i].y - y; ex = dx*T.ix + dy*T.jx; // transformed edge ey = dx*T.iy + dy*T.jy; e2 = ex*ex + ey*ey; if (i == 0 || e2 < e2min) { e2min = e2; closest = &vertexArray[i]; } } return closest; } /*---------------------------------------------------------------------------*/ struct Triangle *lastTriangle; struct Triangle * LocateVertex(struct Vertex *v) /* Triangulation walking. Let's hope there's no cycle in the triangulation */ { struct Triangle *t; struct Vertex *v0, *v1; struct Edge *e; int i; /* initialize */ t = lastTriangle; if (t == NULL) t = &triangleArray[1]; /* search */ again: for (i = 0; i < 3; i++) { v0 = t->vertex[i]; v1 = t->vertex[(i == 2) ? 0 : i+1]; if (Orient2D(v0, v1, v) < 0) { e = t->edge[i]; if (e->triangle[0] == t) t = e->triangle[1]; else if (e->triangle[1] == t) t = e->triangle[0]; else FatalError("edge does not point to neighboring triangle"); goto again; } } /* found */ lastTriangle = t; return t->inside ? t : NULL; } /*---------------------------------------------------------------------------*/ struct Triangle * LocateVertex2(struct Vertex *v) /* This function returns the triangle that contains a vertex "v" * or NULL if the vertex is not found. (NOT CURRENTLY USED) */ { int i; struct Vertex *v0, *v1, *v2; struct Triangle *t; /* walking might not work because our triangulation may have "loops", * so use brute force */ for (i = 0; i < triangleNum; i++) { t = &triangleArray[i]; if (!t->inside) continue; // only check internal triangles v0 = t->vertex[0]; v1 = t->vertex[1]; v2 = t->vertex[2]; if (Orient2D(v0, v1, v) >= 0 && Orient2D(v1, v2, v) >= 0 && Orient2D(v2, v0, v) >= 0) return t; } return NULL; } /*---------------------------------------------------------------------------*/ void Circumcenter(struct Vertex *c, struct Triangle *t) { struct Vertex *v0, *v1, *v2; struct EllipticityTensor T; real c0x, c0y; // centroid real c1x, c1y; // circumcenter real cx, cy, dx, dy, ex, ey, factor, e2, d2; struct Vertex *closest; int i; v0 = t->vertex[0]; // get the 3 vertices of the triangle v1 = t->vertex[1]; v2 = t->vertex[2]; /* compute centroid of triangle */ c0x = (real)(1.0/3.0) * (v0->x + v1->x + v2->x); c0y = (real)(1.0/3.0) * (v0->y + v1->y + v2->y); /* compute the linear transformation at the centroid of the triangle */ Transform(&T, c0x, c0y); /* now compute transformed vectors */ cx = v1->x - v0->x; cy = v1->y - v0->y; dx = cx*T.ix + cy*T.jx; // transformed edge v0-v1 dy = cx*T.iy + cy*T.jy; cx = v2->x - v0->x; cy = v2->y - v0->y; ex = cx*T.ix + cy*T.jx; // transformed edge v0-v2 ey = cx*T.iy + cy*T.jy; /* compute transformed circumcenter */ factor = 0.5/(ex*dy - ey*dx); e2 = ex*ex + ey*ey; d2 = dx*dx + dy*dy; cx = factor * (e2*dy - d2*ey); cy = factor * (d2*ex - e2*dx); /* reconvert into physical space */ factor = 1.0/(T.ix*T.jy - T.jx*T.iy); c1x = v0->x + factor * (T.jy*cx - T.jx*cy); c1y = v0->y + factor * (-T.iy*cx + T.ix*cy); /* try to make this point usable */ for (i = 0; i < 20; i++) { c->x = c1x; c->y = c1y; if (LocateVertex(c) != NULL) { // inside closest = ClosestVertex(c1x, c1y); if (closest == v0 || closest == v1 || closest == v2) return; } c1x = c0x + 0.9*(c1x-c0x); // walk toward centroid c1y = c0y + 0.9*(c1y-c0y); } c->x = c0x; // the walk was long enough: use centroid c->y = c0y; } /*---------------------------------------------------------------------------*/ boolean Encroach(struct Vertex *v0, struct Vertex *v1, struct Vertex *v) /* Does vertex "v" encroach upon segment "v0-v1"? */ { struct EllipticityTensor T; real cx, cy, dx, dy, ex, ey, radius2, dist2; /* compute midpoint of segment */ cx = (real)(0.5) * (v0->x + v1->x); cy = (real)(0.5) * (v0->y + v1->y); /* compute the linear transformation at the midpoint of the segment */ Transform(&T, cx, cy); /* compute "distorted" radius of diametral circle */ dx = v1->x - cx; dy = v1->y - cy; ex = dx*T.ix + dy*T.jx; // transformed edge midpoint-v1 ey = dx*T.iy + dy*T.jy; radius2 = ex*ex + ey*ey; /* compute distance from vertex to midpoint */ dx = v->x - cx; dy = v->y - cy; ex = dx*T.ix + dy*T.jx; // transformed edge midpoint-v ey = dx*T.iy + dy*T.jy; dist2 = ex*ex + ey*ey; return (dist2 < radius2); } /*---------------------------------------------------------------------------*/ boolean EncroachedEdge(struct Edge *e) /* is edge "e" encroached? */ { struct Vertex *v0, *v1; v0 = e->vertex[0]; // get the 2 vertices of the edge v1 = e->vertex[1]; return Encroach(v0, v1, FindApex(e, e->triangle[0])) || Encroach(v0, v1, FindApex(e, e->triangle[1])); } /*---------------------------------------------------------------------------*/ void ChooseSplitPosition(struct Vertex *v, struct Vertex *v0, struct Vertex *v1) { struct Vertex *temp; struct EllipticityTensor T; real dx, dy, d, factor, ex, ey; if (v1 - vertexArray < inputVertexNum) { temp = v0; v0 = v1; v1 = temp; // swap, to get v0 instead } if (v0 - vertexArray < inputVertexNum) { /* compute the linear transformation at v0 */ Transform(&T, v0->x, v0->y); /* to avoid mutual encroachment, choose a power-of-2 distance * in normalized space */ dx = v1->x - v0->x; dy = v1->y - v0->y; ex = dx*T.ix + dy*T.jx; ey = dx*T.iy + dy*T.jy; d = sqrt(ex*ex + ey*ey); factor = pow(2, floor(1.4426950408889634 * log(0.5*d) + 0.5))/d; v->x = v0->x + factor*dx; v->y = v0->y + factor*dy; } else { v->x = 0.5*(v0->x + v1->x); // position: midpoint v->y = 0.5*(v0->y + v1->y); } } /*= TRIANGULATION ===========================================================*/ boolean ValidateEdge(struct Edge *e) /* returns true if any change was made to the mesh */ { struct Vertex *v0, *v1, *v2, *v3; struct Edge *e0, *e1, *e2, *e3; struct Triangle *t0, *t1; int currentCount; real currentMinAngle; real a0, a1, q0, q1; int proposedCount; real proposedMinAngle; if (e->frozen) // we don't flip frozen edges return false; t0 = e->triangle[0]; t1 = e->triangle[1]; if (t0->inside != t1->inside) FatalError("inconsistency between inside and outside"); v0 = e->vertex[0]; v2 = e->vertex[1]; FindNeighbors(e, t0, &e2, &e3, &v3); FindNeighbors(e, t1, &e0, &e1, &v1); /* if edge cannot be flipped, do not consider any further */ if (Orient2D(v1, v3, v2) >= 0 || Orient2D(v1, v3, v0) <= 0) return false; /* compute the current count and minimum angle */ currentCount = 0; if (ContainsBoundingVertices(v0, v2, v3)) currentCount++; if (ContainsBoundingVertices(v2, v0, v1)) currentCount++; a0 = t0->minAngle; a1 = t1->minAngle; currentMinAngle = (a0 < a1) ? a0 : a1; /* compute the proposed count and minimum angle */ proposedCount = 0; if (ContainsBoundingVertices(v1, v3, v0)) proposedCount++; if (ContainsBoundingVertices(v3, v1, v2)) proposedCount++; ComputeTriangleStats(v1, v3, v0, &a0, &q0); ComputeTriangleStats(v3, v1, v2, &a1, &q1); proposedMinAngle = (a0 < a1) ? a0 : a1; /* flip if this locally improves the count or minimum angle */ if (proposedCount < currentCount || proposedMinAngle > currentMinAngle) { SetTriangle(t0, v1, v3, v0, e, e3, e0); // flip SetTriangle(t1, v3, v1, v2, e, e1, e2); SetEdge(e, v1, v3, t0, t1); FixEdge(e0, t1, t0); FixEdge(e2, t0, t1); ValidateEdge(e0); // revalidate touched edges ValidateEdge(e1); ValidateEdge(e2); ValidateEdge(e3); return true; } return false; } /*---------------------------------------------------------------------------*/ void InsertVertexOnEdge(struct Vertex *v, struct Edge *e) { struct Vertex *v0, *v1, *v2, *v3; struct Edge *e0, *e1, *e2, *e3, *f0, *f1, *f2; struct Triangle *t0, *t1, *t2, *t3; boolean inside0, inside1, frozen; t0 = e->triangle[0]; t1 = e->triangle[1]; v0 = e->vertex[0]; v2 = e->vertex[1]; FindNeighbors(e, t0, &e2, &e3, &v3); FindNeighbors(e, t1, &e0, &e1, &v1); /* allocate 2 new triangles and 3 new edges */ t2 = NewTriangle(); t3 = NewTriangle(); f0 = NewEdge(); f1 = NewEdge(); f2 = NewEdge(); /* preserve inside flags and frozen state */ inside0 = t0->inside; inside1 = t1->inside; frozen = e->frozen; /* do all the pointer contorsion to insert new vertex */ SetTriangle(t0, v3, v0, v, e3, e, f2); SetTriangle(t1, v0, v1, v, e0, f0, e); SetTriangle(t2, v1, v2, v, e1, f1, f0); SetTriangle(t3, v2, v3, v, e2, f2, f1); SetEdge(e, v0, v, t0, t1); SetEdge(f0, v1, v, t1, t2); SetEdge(f1, v2, v, t2, t3); SetEdge(f2, v3, v, t3, t0); FixEdge(e1, t1, t2); FixEdge(e2, t0, t3); /* restore inside flags */ t0->inside = inside0; t1->inside = inside1; t2->inside = inside1; t3->inside = inside0; /* set frozen states */ e->frozen = frozen; f1->frozen = frozen; /* fix some edges which are possibly not Delaunay anymore */ if (inside0) { ValidateEdge(e2); ValidateEdge(e3); } if (inside1) { ValidateEdge(e0); ValidateEdge(e1); } } /*---------------------------------------------------------------------------*/ void InsertVertex(struct Vertex *v) { struct Triangle *t0, *t1, *t2; struct Vertex *v0, *v1, *v2; struct Edge *e0, *e1, *e2, *f0, *f1, *f2; int i; /* find triangle "t0" containing "v" */ t0 = LocateVertex(v); if (t0 == NULL) FatalError("unable to locate a vertex"); /* check if vertex is on some triangle edge */ for (i = 0; i < 3; i++) { v0 = t0->vertex[i]; v1 = t0->vertex[(i == 2) ? 0 : i+1]; if (Orient2D(v0, v1, v) == 0) { InsertVertexOnEdge(v, t0->edge[i]); return; } } /* find vertices and edges of "t0" */ v0 = t0->vertex[0]; v1 = t0->vertex[1]; v2 = t0->vertex[2]; e0 = t0->edge[0]; e1 = t0->edge[1]; e2 = t0->edge[2]; /* allocate 2 new triangles and 3 new edges */ t1 = NewTriangle(); t2 = NewTriangle(); f0 = NewEdge(); f1 = NewEdge(); f2 = NewEdge(); /* do all the pointer contorsion to insert new vertex */ SetTriangle(t0, v0, v1, v, e0, f1, f0); SetTriangle(t1, v1, v2, v, e1, f2, f1); SetTriangle(t2, v2, v0, v, e2, f0, f2); SetEdge(f0, v0, v, t2, t0); SetEdge(f1, v1, v, t0, t1); SetEdge(f2, v2, v, t1, t2); FixEdge(e1, t0, t1); FixEdge(e2, t0, t2); /* fix some edges which are possibly not Delaunay anymore */ ValidateEdge(e0); ValidateEdge(e1); ValidateEdge(e2); } /*---------------------------------------------------------------------------*/ void SplitSegment(struct Segment *s) { struct Edge *e; struct Vertex *v0, *v1, *v; struct Segment *t; e = FindEdge(s->endpoint[0], s->endpoint[1]); if (e == NULL) FatalError("some segment was kicked out of the triangulation"); v0 = s->endpoint[0]; v1 = s->endpoint[1]; v = NewVertex(); // create a new vertex if (v == NULL) return; t = NewSegment(); // split segment into 2 subsegments SetSegment(s, v0, v); SetSegment(t, v, v1); ChooseSplitPosition(v, v0, v1); InsertVertexOnEdge(v, e); } /*---------------------------------------------------------------------------*/ void ConstructInitialTriangulation() { int i; struct Vertex *v0, *v1, *v2; struct Triangle *t0, *t1; struct Edge *e0, *e1, *e2; triangleNum = 0; edgeNum = 0; lastTriangle = NULL; v0 = &vertexArray[0]; v1 = &vertexArray[1]; v2 = &vertexArray[2]; /* create the bounding triangle and a "reversed" triangle for * topological reasons */ t0 = NewTriangle(); // triangle #0 t1 = NewTriangle(); // triangle #1 e0 = NewEdge(); e1 = NewEdge(); e2 = NewEdge(); SetTriangle(t0, v0, v1, v2, e0, e1, e2); SetTriangle(t1, v0, v2, v1, e2, e1, e0); SetEdge(e0, v0, v1, t0, t1); SetEdge(e1, v1, v2, t0, t1); SetEdge(e2, v2, v0, t0, t1); for (i = 3; i < vertexNum; i++) InsertVertex(&vertexArray[i]); } /*= DELAUNAY REFINEMENT =====================================================*/ void MayAddSegment(struct Vertex *v0, struct Vertex *v1, struct Vertex *v2) /* This function is a specialized helper for the next function */ { struct Segment *s; if (v0 - vertexArray >= 3 && v1 - vertexArray >= 3 && v2 - vertexArray < 3 && FindSegment(v0, v1) == NULL) { s = NewSegment(); SetSegment(s, v0, v1); } } /*---------------------------------------------------------------------------*/ void ConvexHullBecomesSegments() { int i; struct Vertex *v0, *v1, *v2; for (i = 0; i < triangleNum; i++) { v0 = triangleArray[i].vertex[0]; v1 = triangleArray[i].vertex[1]; v2 = triangleArray[i].vertex[2]; MayAddSegment(v0, v1, v2); MayAddSegment(v1, v2, v0); MayAddSegment(v2, v0, v1); } } /*---------------------------------------------------------------------------*/ void InsertMissingSegments() { int i; struct Vertex *v0, *v1, *v; struct Segment *s, *t; struct Edge *e; boolean change; do { change = false; for (i = 0; i < segmentNum; i++) { s = &segmentArray[i]; v0 = s->endpoint[0]; v1 = s->endpoint[1]; e = FindEdge(v0, v1); if (e == NULL) { v = NewVertex(); // create a new vertex if (v == NULL) return; t = NewSegment(); // split segment into 2 subsegments SetSegment(s, v0, v); SetSegment(t, v, v1); ChooseSplitPosition(v, v0, v1); InsertVertex(v); // insert new vertex into mesh change = true; } // segment is an edge else if (EncroachedEdge(e)) { SplitSegment(s); change = true; } } if (vertexNum == maxVertex) { printf("ERROR: Not enough vertices to make all the segments" " appear\n"); watertight = false; return; } } while (change); watertight = true; } /*---------------------------------------------------------------------------*/ void RemoveHole(struct Triangle *t) { int i; struct Edge *e; if (!t->inside) return; // already outside t->inside = false; // set to outside for (i = 0; i < 3; i++) { e = t->edge[i]; /* if "e" is not a segment, visit the triangle on the other side */ if (FindSegment(e->vertex[0], e->vertex[1]) == NULL) { if (e->triangle[0] == t) RemoveHole(e->triangle[1]); else if (e->triangle[1] == t) RemoveHole(e->triangle[0]); else FatalError("edge does not point to neighboring triangle"); } } } /*---------------------------------------------------------------------------*/ void ClassifyInsideOutside() { int i; struct Triangle *t; RemoveHole(&triangleArray[1]); // remove the exterior shell of triangles for (i = 0; i < holeNum; i++) { t = LocateVertex(&holeArray[i]); if (t == NULL) printf("ERROR: holes #%d is wrong or duplicated\n", i+1); else RemoveHole(t); // remove the hole } } /*---------------------------------------------------------------------------*/ void BasicClassifyInsideOutside() { int i; for (i = 0; i < triangleNum; i++) // just remove the bounding triangle triangleArray[i].inside = ( triangleArray[i].vertex[0] - vertexArray >= 3 && triangleArray[i].vertex[1] - vertexArray >= 3 && triangleArray[i].vertex[2] - vertexArray >= 3); } /*---------------------------------------------------------------------------*/ void FreezeEdgesThatAreSegments() { int i; struct Segment *s; struct Vertex *v0, *v1; struct Edge *e; for (i = 0; i < segmentNum; i++) { s = &segmentArray[i]; v0 = s->endpoint[0]; v1 = s->endpoint[1]; e = FindEdge(v0, v1); if (e != NULL) e->frozen = true; } } /*---------------------------------------------------------------------------*/ void RemoveBadTriangles() { int i; real worstQuality; struct Triangle *t, *worstTriangle; struct Vertex c, *v; boolean encroach; while (vertexNum < maxVertex) { /* find the worst triangle */ worstQuality = 1e30; for (i = 0; i < triangleNum; i++) { t = &triangleArray[i]; if (t->inside && t->quality < worstQuality) { worstQuality = t->quality; worstTriangle = t; } } /* check if we're done */ if (gradingFlag && worstQuality >= angleLowerbound) return; /* compute worst triangle circumcenter */ Circumcenter(&c, worstTriangle); /* does it encroach upon any segment? */ encroach = false; for (i = 0; i < segmentNum; i++) { if (Encroach(segmentArray[i].endpoint[0], segmentArray[i].endpoint[1], &c)) { SplitSegment(&segmentArray[i]); encroach = true; } } if (!encroach) { // did not encroach v = NewVertex(); if (v == NULL) return; *v = c; InsertVertex(v); } } } /*= FILE I/O ================================================================*/ char currentline[INPUTLINESIZE]; char *GetCleanLine(FILE *fp) { char *s; while (fgets(currentline, INPUTLINESIZE, fp) != NULL) { s = currentline; while (*s == ' ' || *s == '\t') s++; // discard leading spaces if (*s == '\n' || *s == '#') continue; // disc. empty lines/comments return s; } return NULL; } /*---------------------------------------------------------------------------*/ /* file format supported (hopefully compatible with the .poly format) * * one line: <# of vertices> * following lines: * one line: <# of segments> * following lines: * one line: <# of holes> * following lines: */ void ProcessPolyFile(char *inputfile, boolean poly, char *programName) { FILE *fp; char *s; int i, pt, length; int end0, end1; int size; double a; char *outputfile; int insideCount, count; struct Triangle *t; struct timeval tv0, tv1; struct timezone tz; double timetaken; /* open input file */ fp = fopen(inputfile, "r"); if (fp == NULL) { printf("ERROR: cannot open the file: %s\n", inputfile); return; } printf("Reading file: %s\n", inputfile); /* read number of vertices */ s = GetCleanLine(fp); if (s == NULL || sscanf(s, "%d", &vertexNum) != 1) { printf("ERROR: cannot read the number of vertices\n"); return; } vertexNum += 3; // want 3 more vertices if (vertexNum > maxVertex) { printf("ERROR: vertexNum > maxVertex\n"); return; } /* read the vertices, starting at index #3 */ for (i = 0; i < 3; i++) { a = 0.123 + i*(TWOPI/3); vertexArray[i].x = VERYLARGE * cos(a); vertexArray[i].y = VERYLARGE * sin(a); } for (i = 3; i < vertexNum; i++) { s = GetCleanLine(fp); if (s == NULL || sscanf(s, "%d%lf%lf", &pt, &vertexArray[i].x, &vertexArray[i].y) != 3) { printf("ERROR: cannot read point #%d\n", i-2); return; } if (i == 3) numberingDelta = pt-i; else if (pt != i+numberingDelta) printf("Warning: point numbers are non-consecutive\n"); } if (poly) { /* read number of segments */ s = GetCleanLine(fp); if (s == NULL || sscanf(s, "%d", &segmentNum) != 1) { printf("ERROR: cannot read the number of segments\n"); return; } /* read the segments */ for (i = 0; i < segmentNum; i++) { s = GetCleanLine(fp); if (s == NULL || sscanf(s, "%d%d%d", &pt, &end0, &end1) != 3) { printf("ERROR: cannot read segment #%d\n", i+1); return; } segmentArray[i].endpoint[0] = &vertexArray[end0 - numberingDelta]; segmentArray[i].endpoint[1] = &vertexArray[end1 - numberingDelta]; } /* read number of holes */ s = GetCleanLine(fp); if (s == NULL || sscanf(s, "%d", &holeNum) != 1) { printf("ERROR: cannot read the number of holes\n"); return; } /* allocate holeArray */ size = holeNum * sizeof(struct Vertex); holeArray = (struct Vertex *)Malloc(size); /* read the holes */ for (i = 0; i < holeNum; i++) { s = GetCleanLine(fp); if (s == NULL || sscanf(s, "%d%lf%lf", &pt, &holeArray[i].x, &holeArray[i].y) != 3) { printf("ERROR: cannot read hole #%d\n", i+1); return; } } } /* close the file */ fclose(fp); /* prepare a string for output filenames */ length = strlen(inputfile); outputfile = (char *)Malloc(length + 1); strcpy(outputfile, inputfile); gettimeofday(&tv0, &tz); ////////// start timer /* construct the mesh */ ComputePSLGBounds(); ConstructInitialTriangulation(); if (poly) { inputVertexNum = vertexNum; // vertexNum will grow if (convexHullFlag) ConvexHullBecomesSegments(); InsertMissingSegments(); if (watertight) { ClassifyInsideOutside(); /* check for any inside triangle */ for (i = 0; i < triangleNum; i++) { if (triangleArray[i].inside) goto ok; } FatalError("every triangle was eaten, you should use the -c flag"); ok:; } else BasicClassifyInsideOutside(); FreezeEdgesThatAreSegments(); RemoveBadTriangles(); Free(holeArray); } else BasicClassifyInsideOutside(); gettimeofday(&tv1, &tz); ////////// stop timer timetaken = TimeDifference(&tv1, &tv0); printf("Computation time: %0.3f s\n", timetaken); if (poly) { /* write .node file */ strcpy(outputfile+length-4, "node"); printf("Writing file: %s (%d vert%s)\n", outputfile, vertexNum-3, (vertexNum-3 <= 1) ? "ex" : "ices"); fp = fopen(outputfile, "w"); if (fp == NULL) { printf("Error, cannot create the file: %s\n", outputfile); return; } fprintf(fp, "# This file was generated by '%s', Francois Labelle's" " Anisotropic\n# Mesh Generator\n\n", programName); fprintf(fp, "%d 2 0 0\n", vertexNum-3); for (i = 3; i < vertexNum; i++) fprintf(fp, "%d %f %f\n", i+numberingDelta, vertexArray[i].x, vertexArray[i].y); fclose(fp); } /* count number of inside triangles */ insideCount = 0; for (i = 0; i < triangleNum; i++) { if (triangleArray[i].inside) insideCount++; } /* write .ele file */ strcpy(outputfile+length-4, "ele"); printf("Writing file: %s (%d triangle%s)\n", outputfile, insideCount, (insideCount <= 1) ? "" : "s"); fp = fopen(outputfile, "w"); if (fp == NULL) { printf("Error, cannot create the file: %s\n", outputfile); return; } fprintf(fp, "# This file was generated by '%s', Francois Labelle's" " Anisotropic\n# Mesh Generator\n\n", programName); fprintf(fp, "%d 3 0\n", insideCount); count = 0; for (i = 0; i < triangleNum; i++) { t = &triangleArray[i]; if (t->inside) fprintf(fp, "%d %d %d %d\n", ++count, (t->vertex[0] - vertexArray) + numberingDelta, (t->vertex[1] - vertexArray) + numberingDelta, (t->vertex[2] - vertexArray) + numberingDelta); } fclose(fp); Free(outputfile); } /*= MAIN PROGRAM ============================================================*/ void ProcessSwitchString(char *s) { int i; struct TransformNode *t; for (i = 0; i < TRANSFORM_NUM; i++) { if (strcmp(s, transformName[i]) == 0) { /* insert transformation into linked list */ t = (struct TransformNode *)Malloc(sizeof(struct TransformNode)); t->which = i; t->next = transformList; transformList = t; return; } } if (strcmp(s, "g") == 0) gradingFlag = true; else if (s[0] == 'a') { if (sscanf(s+1, "%lf", &angleLowerbound) != 1) printf("ERROR, incomplete option: %s\n", s); } else if (s[0] == 'v') { if (sscanf(s+1, "%d", &maxVertex) != 1) printf("ERROR, incomplete option: %s\n", s); } else if (strcmp(s, "c") == 0) convexHullFlag = true; else printf("ERROR, unrecognized option: %s\n", s); } /*---------------------------------------------------------------------------*/ int main(int argc, char **argv) { int i, length; boolean gotfile = false; for (i = 0; i < TRANSFORM_NUM; i++) transformFlag[i] = false; /* process the arguments for command-line switches */ for (i = 1; i < argc; i++) { if (argv[i][0] == '-') ProcessSwitchString(argv[i] + 1); } AllocateMeshingStructures(); /* reprocess the arguments for .node and .poly files */ for (i = 1; i < argc; i++) { if (argv[i][0] == '-') continue; length = strlen(argv[i]); if (length >= 5 && strcmp(argv[i] + length - 5, ".node") == 0) { ProcessPolyFile(argv[i], false, argv[0]); gotfile = true; } else if (length >= 5 && strcmp(argv[i] + length - 5, ".poly") == 0) { ProcessPolyFile(argv[i], true, argv[0]); gotfile = true; } else printf("The following argument was ignored: %s\n", argv[i]); } FreeMeshingStructures(); if (!gotfile) { // no file, user is probably asking for help... printf("usage: %s -option1 -option2 [...] file.node file.poly [...]\n", argv[0]); printf(" .node files will be triangulated (output: .ele)\n"); printf(" .poly files will be meshed (output: .node and .ele)\n"); printf("meshing options:\n"); printf(" -g grading mesh (size of triangles irrelevent)\n"); printf(" -a# angle lower bound in degrees,"); printf(" only with -g option (default %1.1f)\n", angleLowerbound); printf(" -v# vertex upper bound"); printf(" (default %d, ~%1.0f K RAM)\n", maxVertex-3, (maxVertex-3) * BYTES_PER_VERTEX / 1024.0); printf(" -c convex hull edges become segments\n"); printf("anisotropy options: (can be combined!)\n"); printf(" affecting shape:\n"); printf(" -hstretch, -vstretch, -dstretch, -sink, -swirl\n"); printf(" affecting size:\n"); printf(" -hline, -vline, -dline, -center, -perimeter\n"); printf(" -left, -right, -top, -bottom\n"); printf(" user specified:\n"); printf(" -custom (the function 'T_custom' in 'transform.c')\n"); } return 0; } /*---------------------------------------------------------------------------*/