diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..ccc9762 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,19 @@ +{ + "files.associations": { + "__hash_table": "cpp", + "__split_buffer": "cpp", + "array": "cpp", + "bitset": "cpp", + "deque": "cpp", + "initializer_list": "cpp", + "iterator": "cpp", + "queue": "cpp", + "random": "cpp", + "stack": "cpp", + "string": "cpp", + "string_view": "cpp", + "unordered_map": "cpp", + "utility": "cpp", + "vector": "cpp" + } +} \ No newline at end of file diff --git a/delaunay.h b/delaunay.h index 160ea4a..ba4b535 100644 --- a/delaunay.h +++ b/delaunay.h @@ -1,14 +1,14 @@ /** * @defgroup DELAUNAY * - * @brief An implementation of the 2D Delaunay triangulation using the Bowyer-Watson algorithm. + * @brief An implementation of the 2D Delaunay triangulation using the flip algorithm. * * @author Yi Zhang - * @date 2021-09-12 + * @date 2021-09-19 */ -#ifndef _BW_2D_DELAUNAY_H -#define _BW_2D_DELAUNAY_H +#ifndef _FLIP_2D_DELAUNAY_H +#define _FLIP_2D_DELAUNAY_H #include "cmath" #include "vector" @@ -37,56 +37,181 @@ bool operator ==(const vertex2dc &a, const vertex2dc &b) // overload the == oper } return false; } -// End vertex definition -// Start edge definition -struct edge +void circumcircle(vertex2dc *v0, vertex2dc *v1, vertex2dc *v2, double &cx, double &cy, double &cr) // calculate the circumcircle from three points { - vertex2dc *vert[2]; // vertex of the edge + double s = 0.5 / ((v1->x - v0->x) * (v2->y - v0->y) - (v1->y - v0->y) * (v2->x - v0->x)); + double m = v1->x*v1->x - v0->x*v0->x + v1->y*v1->y - v0->y*v0->y; + double u = v2->x*v2->x - v0->x*v0->x + v2->y*v2->y - v0->y*v0->y; - edge() {vert[0] = vert[1] = nullptr;} - edge(vertex2dc *v0ptr, vertex2dc *v1ptr) {set(v0ptr, v1ptr);} - void set(vertex2dc *v0ptr, vertex2dc *v1ptr) - { - vert[0] = v0ptr; vert[1] = v1ptr; - return; - } -}; - -bool operator ==(const edge &a, const edge &b) // overload the == operator for edge type -{ - if((a.vert[0] == b.vert[0] && a.vert[1] == b.vert[1]) || - (a.vert[0] == b.vert[1] && a.vert[1] == b.vert[0])) - { - return true; - } - return false; + cx = ((v2->y - v0->y)*m + (v0->y - v1->y)*u)*s; + cy = ((v0->x - v2->x)*m + (v1->x - v0->x)*u)*s; + cr = (v0->x - cx)*(v0->x - cx) + (v0->y - cy)*(v0->y - cy); // not need to calculate the squared root here + return; } -// End edge definition +// End vertex definition // Start triangle definition struct triangle { + int id; // index of the triangle vertex2dc *vert[3]; // vertex of the triangle + triangle *neigh[3]; // neighbors of the triangle double cx, cy; // center of the triangle's circumcircle double cr; // radius of the circumcircle - triangle() {vert[0] = vert[1] = vert[2] = nullptr;} + triangle() {vert[0] = vert[1] = vert[2] = nullptr; neigh[0] = neigh[1] = neigh[2] = nullptr;} triangle(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) {set(v0ptr, v1ptr, v2ptr);} void set(vertex2dc *v0ptr, vertex2dc *v1ptr, vertex2dc *v2ptr) { vert[0] = v0ptr; vert[1] = v1ptr; vert[2] = v2ptr; - - double s = 0.5 / ((vert[1]->x - vert[0]->x) * (vert[2]->y - vert[0]->y) - (vert[1]->y - vert[0]->y) * (vert[2]->x - vert[0]->x)); - double m = vert[1]->x * vert[1]->x - vert[0]->x * vert[0]->x + vert[1]->y * vert[1]->y - vert[0]->y * vert[0]->y; - double u = vert[2]->x * vert[2]->x - vert[0]->x * vert[0]->x + vert[2]->y * vert[2]->y - vert[0]->y * vert[0]->y; - - cx = ((vert[2]->y - vert[0]->y) * m + (vert[0]->y - vert[1]->y) * u) * s; - cy = ((vert[0]->x - vert[2]->x) * m + (vert[1]->x - vert[0]->x) * u) * s; - cr = (vert[0]->x - cx) * (vert[0]->x - cx) + (vert[0]->y - cy) * (vert[0]->y - cy); // not need to sqrt() here + neigh[0] = neigh[1] = neigh[2] = nullptr; + circumcircle(vert[0], vert[1], vert[2], cx, cy, cr); return; } + + void set_neighbor(triangle *n0ptr, triangle *n1ptr, triangle *n2ptr) + { + neigh[0] = n0ptr; neigh[1] = n1ptr; neigh[2] = n2ptr; + return; + } + + bool bound_location(double inx, double iny) // Test if the location is inside the triangle + { + double l1x, l1y, l2x, l2y; + for (int i = 0; i < 3; i++) + { + l1x = vert[(i+1)%3]->x - vert[i]->x; + l1y = vert[(i+1)%3]->y - vert[i]->y; + l2x = inx - vert[i]->x; + l2y = iny - vert[i]->y; + + if ((l1x*l2y - l1y*l2x) < 0) // This condition includes points on the triangle's edge + { + return false; + } + } + return true; + } }; + +/** + * @brief Flip neighboring triangles and their neighbors + * + * original + * + * /\ + * / \ + * / \ + * / t \ + * t_id-------\ t_id (0, 1 or 2) + * \--------/ + * \ / + * \ n / + * \ / + * \/ + * n_id (0, 1 or 2) + * + * fliped + * + * /|\ + * / | \ + * / | \ + * / | \ + * t_id | \ t_id (0, 1 or 2) + * \ t | n / + * \ | / + * \ | / + * \ | / + * \|/ + * n_id (0, 1 or 2) + * + * @param t target triangle + * @param n neighboring triangle + * @param t_vid reference index of the target triangle + * @param n_vid reference index of the neighboring triangle + */ +void flip_neighboring_triangles(triangle *t, triangle *n, int t_id, int n_id) +{ + t->vert[(t_id+1)%3] = n->vert[n_id]; // flip t + circumcircle(t->vert[0], t->vert[1], t->vert[2], t->cx, t->cy, t->cr); // update circumcircle + + n->vert[(n_id+2)%3] = t->vert[(t_id+2)%3]; // flip n + circumcircle(n->vert[0], n->vert[1], n->vert[2], n->cx, n->cy, n->cr); // update circumcircle + + // set side neighbors + t->neigh[t_id] = n->neigh[(n_id+2)%3]; + n->neigh[(n_id+1)%3] = t->neigh[(t_id+1)%3]; + + // set opposite neighbors + t->neigh[(t_id+1)%3] = n; + n->neigh[(n_id+2)%3] = t; + + // set oppsite neighbors + if (t->neigh[t_id] != nullptr) + { + for (int i = 0; i < 3; i++) + { + if (t->neigh[t_id]->neigh[i] == n) + { + t->neigh[t_id]->neigh[i] = t; + break; + } + } + } + + if (n->neigh[(n_id+1)%3] != nullptr) + { + for (int i = 0; i < 3; i++) + { + if (n->neigh[(n_id+1)%3]->neigh[i] == t) + { + n->neigh[(n_id+1)%3]->neigh[i] = n; + break; + } + } + } + return; +} + +/** + * @brief Make sure that the input triangle meet the empty circumcircle condition + * + * @param t Input triangle + */ +void make_delaunay(triangle *t) +{ + double dist; + vertex2dc *n_vert; + triangle *n_tri; + for (int n = 0; n < 3; ++n) + { + if (t->neigh[n] != nullptr) // must has neighbor on this side + { + n_tri = t->neigh[n]; + for (int v = 0; v < 3; ++v) + { + n_vert = n_tri->vert[v]; + if (n_vert != t->vert[n] && n_vert != t->vert[(n+1)%3]) // find the opposite vertex + { + dist = (t->cx - n_vert->x) * (t->cx - n_vert->x) + + (t->cy - n_vert->y) * (t->cy - n_vert->y); + + if ((dist - t->cr) <= ZERO) // need to be flipped + { + flip_neighboring_triangles(t, n_tri, n, v); + // Make sure the triangles also meet the empty circumcircle condition after flipping + make_delaunay(t); + make_delaunay(n_tri); + return; // Neighborhood changed. The current loop is not valid any more. + } + break; // no need to search more + } + } + } + } + return; +} // End triangle definition /** @@ -95,7 +220,7 @@ struct triangle * @param in_verts Input vertexes. Defined by the user. * @param out_tris Output triangles. Compute by the function. */ -void triangulation(std::vector &in_verts, std::vector &out_tris) +void triangulation(std::vector &in_verts, std::vector &out_tris) { if (!out_tris.empty()) out_tris.clear(); if (in_verts.size() < 3) return; @@ -116,125 +241,155 @@ void triangulation(std::vector &in_verts, std::vector &out_ double maxi_s = std::max(xmax - xmin, ymax - ymin); // use an four times bigger rectangle to include all points vertex2dc *tmp_vert = nullptr; - std::vector assit_vert; + std::vector box_vert; tmp_vert = new vertex2dc(midx - maxi_s, midy - maxi_s); // lower left corner - assit_vert.push_back(tmp_vert); + box_vert.push_back(tmp_vert); tmp_vert = new vertex2dc(midx + maxi_s, midy - maxi_s); // lower right corner - assit_vert.push_back(tmp_vert); + box_vert.push_back(tmp_vert); tmp_vert = new vertex2dc(midx + maxi_s, midy + maxi_s); // upper right corner - assit_vert.push_back(tmp_vert); + box_vert.push_back(tmp_vert); tmp_vert = new vertex2dc(midx - maxi_s, midy + maxi_s); // upper left corner - assit_vert.push_back(tmp_vert); + box_vert.push_back(tmp_vert); - triangle *tmp_tri = nullptr; - std::vector exist_tri, cnst_tri; + triangle *old_tri = nullptr, *tmp_tri = nullptr; + triangle *cnst_tri[3]; std::vector::iterator t_iter; - tmp_tri = new triangle(assit_vert[0], assit_vert[1], assit_vert[2]); // order the vertex anti-clock wise - exist_tri.push_back(tmp_tri); + tmp_tri = new triangle(box_vert[0], box_vert[1], box_vert[2]); // order the vertex anti-clock wise + out_tris.push_back(tmp_tri); - tmp_tri = new triangle(assit_vert[0], assit_vert[2], assit_vert[3]); // order the vertex anti-clock wise - exist_tri.push_back(tmp_tri); + tmp_tri = new triangle(box_vert[0], box_vert[2], box_vert[3]); // order the vertex anti-clock wise + out_tris.push_back(tmp_tri); + + out_tris[0]->set_neighbor(nullptr, nullptr, out_tris[1]); + out_tris[1]->set_neighbor(out_tris[0], nullptr, nullptr); // loop all input vertice - bool removed; - double dist; - edge tmp_edge; - std::vector cnst_edge; - std::vector::iterator e_iter; for (int i = 0; i < in_verts.size(); ++i) { - // determine triangles that include the point and add the triangle to the cnst_tri and remove the triangle from exist_tri - // this is the part that could take a lot of time if we are working with a large amount of points. We will fix this later - cnst_tri.clear(); - for (t_iter = exist_tri.begin(); t_iter != exist_tri.end(); ) + // determine the triangle that includes the new vertex and remove it from out_tris + for (t_iter = out_tris.begin(); t_iter != out_tris.end(); ) { - tmp_tri = *t_iter; - - dist = (tmp_tri->cx - in_verts[i].x) * (tmp_tri->cx - in_verts[i].x) + - (tmp_tri->cy - in_verts[i].y) * (tmp_tri->cy - in_verts[i].y); - - if ((dist - tmp_tri->cr) <= ZERO) // Points on the circumcircle are also included + old_tri = *t_iter; + if (old_tri->bound_location(in_verts[i].x, in_verts[i].y)) { - t_iter = exist_tri.erase(t_iter); - cnst_tri.push_back(tmp_tri); + t_iter = out_tris.erase(t_iter); + break; } else t_iter++; } - // loop to remove duplicate edges - cnst_edge.clear(); - for (int c = 0; c < cnst_tri.size(); ++c) + // build three new triangles + for (int n = 0; n < 3; ++n) { - for (int e = 0; e < 3; ++e) + tmp_tri = new triangle(old_tri->vert[n], old_tri->vert[(n+1)%3], &in_verts[i]); + cnst_tri[n] = tmp_tri; + out_tris.push_back(tmp_tri); + } + + // sort neighbors + for (int n = 0; n < 3; ++n) + { + if (old_tri->neigh[n] == nullptr) { - tmp_edge.set(cnst_tri[c]->vert[e], cnst_tri[c]->vert[(e+1)%3]); - - removed = false; - for (e_iter = cnst_edge.begin(); e_iter != cnst_edge.end(); ) + cnst_tri[n]->set_neighbor(nullptr, cnst_tri[(n+1)%3], cnst_tri[(n+2)%3]); + } + else + { + cnst_tri[n]->set_neighbor(old_tri->neigh[n], cnst_tri[(n+1)%3], cnst_tri[(n+2)%3]); + for (int k = 0; k < 3; ++k) // replace neighbor for the oppositing triangle { - if (tmp_edge == *e_iter) // duplicate edge, remove from cnst_edge + if (old_tri->neigh[n]->neigh[k] == old_tri) { - e_iter = cnst_edge.erase(e_iter); - removed = true; - break; // no need to search more + old_tri->neigh[n]->neigh[k] = cnst_tri[n]; + break; } - else e_iter++; - } - - if (!removed) // not a duplicate edge, add to the cnst_edge - { - cnst_edge.push_back(tmp_edge); } } } - // construct new triangles and add to exist_tri - for (int c = 0; c < cnst_edge.size(); ++c) - { - tmp_tri = new triangle(cnst_edge[c].vert[0], cnst_edge[c].vert[1], &in_verts[i]); // order the vertex anti-clock wise - exist_tri.push_back(tmp_tri); - } + // delete the old triangle + delete old_tri; old_tri = nullptr; - // destroy memories used by cnst_edge - for (int c = 0; c < cnst_tri.size(); ++c) + // Make sure cnst_tri meet the empty circumcircle condition + for (int n = 0; n < 3; ++n) { - tmp_tri = cnst_tri[c]; - delete tmp_tri; tmp_tri = nullptr; + make_delaunay(cnst_tri[n]); } } - // remove any triangles has an assistant vertex from exist_tri - for (t_iter = exist_tri.begin(); t_iter != exist_tri.end(); ) + // remove any triangles has an box vertex from out_tris + for (t_iter = out_tris.begin(); t_iter != out_tris.end(); ) { tmp_tri = *t_iter; - if (tmp_tri->vert[0] == assit_vert[0] || tmp_tri->vert[0] == assit_vert[1] || tmp_tri->vert[0] == assit_vert[2] || tmp_tri->vert[0] == assit_vert[3] || - tmp_tri->vert[1] == assit_vert[0] || tmp_tri->vert[1] == assit_vert[1] || tmp_tri->vert[1] == assit_vert[2] || tmp_tri->vert[1] == assit_vert[3] || - tmp_tri->vert[2] == assit_vert[0] || tmp_tri->vert[2] == assit_vert[1] || tmp_tri->vert[2] == assit_vert[2] || tmp_tri->vert[2] == assit_vert[3]) + if (tmp_tri->vert[0] == box_vert[0] || tmp_tri->vert[0] == box_vert[1] || tmp_tri->vert[0] == box_vert[2] || tmp_tri->vert[0] == box_vert[3]) { + if (tmp_tri->neigh[1] != nullptr) + { + for (int k = 0; k < 3; ++k) + { + if (tmp_tri->neigh[1]->neigh[k] == tmp_tri) + { + tmp_tri->neigh[1]->neigh[k] = nullptr; + break; + } + } + } // destroy the memories located and remove from the vector - t_iter = exist_tri.erase(t_iter); + t_iter = out_tris.erase(t_iter); + delete tmp_tri; tmp_tri = nullptr; + } + else if (tmp_tri->vert[1] == box_vert[0] || tmp_tri->vert[1] == box_vert[1] || tmp_tri->vert[1] == box_vert[2] || tmp_tri->vert[1] == box_vert[3]) + { + if (tmp_tri->neigh[2] != nullptr) + { + for (int k = 0; k < 3; ++k) + { + if (tmp_tri->neigh[2]->neigh[k] == tmp_tri) + { + tmp_tri->neigh[2]->neigh[k] = nullptr; + break; + } + } + } + // destroy the memories located and remove from the vector + t_iter = out_tris.erase(t_iter); + delete tmp_tri; tmp_tri = nullptr; + } + else if (tmp_tri->vert[2] == box_vert[0] || tmp_tri->vert[2] == box_vert[1] || tmp_tri->vert[2] == box_vert[2] || tmp_tri->vert[2] == box_vert[3]) + { + if (tmp_tri->neigh[0] != nullptr) + { + for (int k = 0; k < 3; ++k) + { + if (tmp_tri->neigh[0]->neigh[k] == tmp_tri) + { + tmp_tri->neigh[0]->neigh[k] = nullptr; + break; + } + } + } + // destroy the memories located and remove from the vector + t_iter = out_tris.erase(t_iter); delete tmp_tri; tmp_tri = nullptr; } else t_iter++; } - // copy exist_tri to out_tris and destroy memories located - out_tris.resize(exist_tri.size()); - for (int i = 0; i < exist_tri.size(); ++i) + // assign triangles index + for (int i = 0; i < out_tris.size(); i++) { - out_tris[i].set(exist_tri[i]->vert[0], exist_tri[i]->vert[1], exist_tri[i]->vert[2]); - delete exist_tri[i]; exist_tri[i] = nullptr; + out_tris[i]->id = i; } - - // destroy memories located for assit_vert + + // destroy memories located for box_vert for (int i = 0; i < 4; ++i) { - delete assit_vert[i]; assit_vert[i] = nullptr; + delete box_vert[i]; box_vert[i] = nullptr; } return; } @@ -271,9 +426,9 @@ bool duplicated_vertex(const std::vector &in_verts) * * @return If the triangulation is fully delaunay */ -bool fully_delaunay(const std::vector &in_tris, const std::vector &in_verts) +bool fully_delaunay(const std::vector &in_tris, const std::vector &in_verts) { - if (in_tris.empty()) return true; + if (in_tris.empty()) return false; int count; double dist; @@ -282,10 +437,10 @@ bool fully_delaunay(const std::vector &in_tris, const std::vectorcx - in_verts[j].x) * (in_tris[i]->cx - in_verts[j].x) + + (in_tris[i]->cy - in_verts[j].y) * (in_tris[i]->cy - in_verts[j].y); - if ((dist - in_tris[i].cr) <= ZERO) + if ((dist - in_tris[i]->cr) <= ZERO) { count++; } @@ -300,4 +455,4 @@ bool fully_delaunay(const std::vector &in_tris, const std::vectorx - vert[0]->x) * (vert[2]->y - vert[0]->y) - (vert[1]->y - vert[0]->y) * (vert[2]->x - vert[0]->x)); + double m = vert[1]->x * vert[1]->x - vert[0]->x * vert[0]->x + vert[1]->y * vert[1]->y - vert[0]->y * vert[0]->y; + double u = vert[2]->x * vert[2]->x - vert[0]->x * vert[0]->x + vert[2]->y * vert[2]->y - vert[0]->y * vert[0]->y; + + cx = ((vert[2]->y - vert[0]->y) * m + (vert[0]->y - vert[1]->y) * u) * s; + cy = ((vert[0]->x - vert[2]->x) * m + (vert[1]->x - vert[0]->x) * u) * s; + cr = (vert[0]->x - cx) * (vert[0]->x - cx) + (vert[0]->y - cy) * (vert[0]->y - cy); // not need to sqrt() here + return; + } +}; +// End triangle definition + +/** + * @brief 2D Delaunay triangulation of some given points using the Bowyer-Watson algorithm. + * + * @param in_verts Input vertexes. Defined by the user. + * @param out_tris Output triangles. Compute by the function. + */ +void triangulation(std::vector &in_verts, std::vector &out_tris, int valid_size) +{ + if (!out_tris.empty()) out_tris.clear(); + if (in_verts.size() < 3) return; + + // locate the surrounding box and initiate the staring two triangles + double xmin = in_verts[0].x, xmax = in_verts[0].x; + double ymin = in_verts[0].y, ymax = in_verts[0].y; + for (int i = 0; i < in_verts.size(); ++i) + { + xmin = std::min(xmin, in_verts[i].x); + xmax = std::max(xmax, in_verts[i].x); + ymin = std::min(ymin, in_verts[i].y); + ymax = std::max(ymax, in_verts[i].y); + } + + double midx = 0.5*(xmin + xmax); + double midy = 0.5*(ymin + ymax); + double maxi_s = std::max(xmax - xmin, ymax - ymin); // use an four times bigger rectangle to include all points + + vertex2dc *tmp_vert = nullptr; + std::vector assit_vert; + + tmp_vert = new vertex2dc(midx - maxi_s, midy - maxi_s, 0); // lower left corner + assit_vert.push_back(tmp_vert); + + tmp_vert = new vertex2dc(midx + maxi_s, midy - maxi_s, 1); // lower right corner + assit_vert.push_back(tmp_vert); + + tmp_vert = new vertex2dc(midx + maxi_s, midy + maxi_s, 2); // upper right corner + assit_vert.push_back(tmp_vert); + + tmp_vert = new vertex2dc(midx - maxi_s, midy + maxi_s, 3); // upper left corner + assit_vert.push_back(tmp_vert); + + triangle *tmp_tri = nullptr; + std::vector exist_tri, cnst_tri; + std::vector::iterator t_iter; + + tmp_tri = new triangle(assit_vert[0], assit_vert[1], assit_vert[2]); // order the vertex anti-clock wise + exist_tri.push_back(tmp_tri); + + tmp_tri = new triangle(assit_vert[0], assit_vert[2], assit_vert[3]); // order the vertex anti-clock wise + exist_tri.push_back(tmp_tri); + + // loop all input vertice + bool removed; + double dist; + edge tmp_edge; + std::vector cnst_edge; + std::vector::iterator e_iter; + //for (int i = 0; i < in_verts.size(); ++i) + for (int i = 0; i < valid_size; ++i) + { + // determine triangles that include the point and add the triangle to the cnst_tri and remove the triangle from exist_tri + // this is the part that could take a lot of time if we are working with a large amount of points. We will fix this later + cnst_tri.clear(); + for (t_iter = exist_tri.begin(); t_iter != exist_tri.end(); ) + { + tmp_tri = *t_iter; + + dist = (tmp_tri->cx - in_verts[i].x) * (tmp_tri->cx - in_verts[i].x) + + (tmp_tri->cy - in_verts[i].y) * (tmp_tri->cy - in_verts[i].y); + + if ((dist - tmp_tri->cr) <= ZERO) // Points on the circumcircle are also included + { + t_iter = exist_tri.erase(t_iter); + cnst_tri.push_back(tmp_tri); + } + else t_iter++; + } + + // loop to remove duplicate edges + cnst_edge.clear(); + for (int c = 0; c < cnst_tri.size(); ++c) + { + for (int e = 0; e < 3; ++e) + { + tmp_edge.set(cnst_tri[c]->vert[e], cnst_tri[c]->vert[(e+1)%3]); + + removed = false; + for (e_iter = cnst_edge.begin(); e_iter != cnst_edge.end(); ) + { + if (tmp_edge == *e_iter) // duplicate edge, remove from cnst_edge + { + e_iter = cnst_edge.erase(e_iter); + removed = true; + break; // no need to search more + } + else e_iter++; + } + + if (!removed) // not a duplicate edge, add to the cnst_edge + { + cnst_edge.push_back(tmp_edge); + } + } + } + + // construct new triangles and add to exist_tri + for (int c = 0; c < cnst_edge.size(); ++c) + { + tmp_tri = new triangle(cnst_edge[c].vert[0], cnst_edge[c].vert[1], &in_verts[i]); // order the vertex anti-clock wise + exist_tri.push_back(tmp_tri); + } + + // destroy memories used by cnst_edge + for (int c = 0; c < cnst_tri.size(); ++c) + { + tmp_tri = cnst_tri[c]; + delete tmp_tri; tmp_tri = nullptr; + } + } + + // Start testing code + std::ofstream outfile("test_backup.msh"); + outfile << "$MeshFormat" << std::endl << "2.2 0 8" << std::endl << "$EndMeshFormat "<x << " " << assit_vert[i]->y << " 0.0" << std::endl; + } + for (int i = 0; i < valid_size; i++) + { + outfile << i + 5 << " " << std::setprecision(16) + << in_verts[i].x << " " << in_verts[i].y << " 0.0" << std::endl; + } + outfile<<"$EndNodes"<vert[j]->id + 1; + } + outfile << std::endl; + } + outfile << "$EndElements"<< std::endl; + outfile.close(); + // End testing code + + // remove any triangles has an assistant vertex from exist_tri + for (t_iter = exist_tri.begin(); t_iter != exist_tri.end(); ) + { + tmp_tri = *t_iter; + if (tmp_tri->vert[0] == assit_vert[0] || tmp_tri->vert[0] == assit_vert[1] || tmp_tri->vert[0] == assit_vert[2] || tmp_tri->vert[0] == assit_vert[3] || + tmp_tri->vert[1] == assit_vert[0] || tmp_tri->vert[1] == assit_vert[1] || tmp_tri->vert[1] == assit_vert[2] || tmp_tri->vert[1] == assit_vert[3] || + tmp_tri->vert[2] == assit_vert[0] || tmp_tri->vert[2] == assit_vert[1] || tmp_tri->vert[2] == assit_vert[2] || tmp_tri->vert[2] == assit_vert[3]) + { + // destroy the memories located and remove from the vector + t_iter = exist_tri.erase(t_iter); + delete tmp_tri; tmp_tri = nullptr; + } + else t_iter++; + } + + // copy exist_tri to out_tris and destroy memories located + out_tris.resize(exist_tri.size()); + for (int i = 0; i < exist_tri.size(); ++i) + { + out_tris[i].set(exist_tri[i]->vert[0], exist_tri[i]->vert[1], exist_tri[i]->vert[2]); + delete exist_tri[i]; exist_tri[i] = nullptr; + } + + // destroy memories located for assit_vert + for (int i = 0; i < 4; ++i) + { + delete assit_vert[i]; assit_vert[i] = nullptr; + } + return; +} + +/** + * @brief Check for duplicated vertex + * + * @param[in] in_verts Input vertexes + * + * @return If there is duplicated vertex + */ +bool duplicated_vertex(const std::vector &in_verts) +{ + if (in_verts.empty()) return false; + + for (int i = 0; i < in_verts.size()-1; ++i) + { + for (int j = i+1; j < in_verts.size(); ++j) + { + if (in_verts[i] == in_verts[j] && in_verts[i].id != in_verts[j].id) + { + return true; + } + } + } + return false; +} + +/** + * @brief Check to see if the triangulation is fully delaunay + * + * @param[in] in_tris Input triangles + * @param[in] in_verts Input vertexes + * + * @return If the triangulation is fully delaunay + */ +bool fully_delaunay(const std::vector &in_tris, const std::vector &in_verts) +{ + if (in_tris.empty()) return true; + + int count; + double dist; + for (int i = 0; i < in_tris.size(); ++i) + { + count = 0; + for (int j = 0; j < in_verts.size(); ++j) + { + dist = (in_tris[i].cx - in_verts[j].x) * (in_tris[i].cx - in_verts[j].x) + + (in_tris[i].cy - in_verts[j].y) * (in_tris[i].cy - in_verts[j].y); + + if ((dist - in_tris[i].cr) <= ZERO) + { + count++; + } + } + + if (count > 3) + { + return false; + } + } + + return true; +} + +#endif // _BW_2D_DELAUNAY_H \ No newline at end of file diff --git a/demo.cpp b/demo.cpp index 4d68961..a349490 100644 --- a/demo.cpp +++ b/demo.cpp @@ -5,6 +5,7 @@ int main(int argc, char const *argv[]) { + std::vector points(21); points[0].set(-0.8, -0.8, 0); points[1].set(0.4, -1.2, 1); @@ -28,13 +29,38 @@ int main(int argc, char const *argv[]) points[19].set(3.5, 1.8, 19); points[20].set(3.6, 3.1, 20); + /* + std::vector points(21); + points[0].set(-0.8, -0.8, 4); + points[1].set(0.4, -1.2, 5); + points[2].set(1.2, -0.9, 6); + points[3].set(1.6, 0.1, 7); + points[4].set(2.5, 0.5, 8); + points[5].set(4.1, 0.7, 9); + points[6].set(5.7, 1.8, 10); + points[7].set(5.1, 3.4, 11); + points[8].set(2.5, 4.4, 12); + points[9].set(1.2, 3.7, 13); + points[10].set(-1.2, 3.9, 14); + points[11].set(-3.2, 5.1, 15); + points[12].set(-4.3, 2.9, 16); + points[13].set(-3.1, 0.7, 17); + points[14].set(-1.3, 0.6, 18); + points[15].set(-2.1, 2.9, 19); + points[16].set(0.6, 1.2, 20); + points[17].set(0.1, 2.4, 21); + points[18].set(2.4, 2.8, 22); + points[19].set(3.5, 1.8, 23); + points[20].set(3.6, 3.1, 24); + */ + if (duplicated_vertex(points)) { - std::cerr << "Duplicated vertexes detected.\n"; + std::cerr << "Duplicated vertice detected.\n"; return 0; } - std::vector elements; + std::vector elements; triangulation(points, elements); if (fully_delaunay(elements, points)) @@ -59,13 +85,35 @@ int main(int argc, char const *argv[]) outfile << i + 1 << " 2 0"; for (int j = 0; j < 3; j++) { - outfile << " " << elements[i].vert[j]->id + 1; + outfile << " " << elements[i]->vert[j]->id + 1; } outfile << std::endl; } outfile << "$EndElements"<< std::endl; outfile.close(); + + // write a neighbor file + outfile.open("demo.neigh"); + outfile << elements.size() << std::endl; + for (int i = 0; i < elements.size(); i++) + { + outfile << i + 1; + for (int j = 0; j < 3; j++) + { + if (elements[i]->neigh[j] != nullptr) + { + outfile << " " << elements[i]->neigh[j]->id + 1; + } + else outfile << " -1"; + } + outfile << std::endl; + } + outfile.close(); + + // destroy allocated memories + for (int i = 0; i < elements.size(); i++) + { + delete elements[i]; + } return 0; } - - diff --git a/demo.msh b/demo.msh index fb98d1d..2ce06d2 100644 --- a/demo.msh +++ b/demo.msh @@ -27,34 +27,34 @@ $Nodes $EndNodes $Elements 30 -1 2 0 2 3 4 -2 2 0 4 3 5 -3 2 0 5 3 6 -4 2 0 10 9 11 -5 2 0 11 9 12 -6 2 0 14 1 15 -7 2 0 11 12 16 -8 2 0 12 13 16 -9 2 0 13 14 16 -10 2 0 14 15 16 +1 2 0 12 13 16 +2 2 0 17 5 19 +3 2 0 3 4 2 +4 2 0 17 2 4 +5 2 0 4 3 5 +6 2 0 7 8 20 +7 2 0 5 3 6 +8 2 0 18 16 15 +9 2 0 10 9 11 +10 2 0 11 9 12 11 2 0 1 2 17 -12 2 0 2 4 17 -13 2 0 4 5 17 -14 2 0 15 1 17 -15 2 0 10 11 18 -16 2 0 11 16 18 -17 2 0 16 15 18 -18 2 0 15 17 18 -19 2 0 9 10 19 -20 2 0 17 5 19 -21 2 0 10 18 19 -22 2 0 18 17 19 -23 2 0 6 7 20 -24 2 0 7 8 20 -25 2 0 5 6 20 -26 2 0 19 5 20 -27 2 0 8 9 21 -28 2 0 9 19 21 -29 2 0 19 20 21 -30 2 0 20 8 21 +12 2 0 17 18 15 +13 2 0 14 1 15 +14 2 0 13 14 16 +15 2 0 14 15 16 +16 2 0 11 12 16 +17 2 0 15 1 17 +18 2 0 4 5 17 +19 2 0 10 11 18 +20 2 0 11 16 18 +21 2 0 17 19 18 +22 2 0 10 18 19 +23 2 0 5 6 20 +24 2 0 9 10 19 +25 2 0 19 5 20 +26 2 0 6 7 20 +27 2 0 8 21 20 +28 2 0 19 20 21 +29 2 0 8 9 21 +30 2 0 9 19 21 $EndElements diff --git a/demo.neigh b/demo.neigh new file mode 100644 index 0000000..4eb1cef --- /dev/null +++ b/demo.neigh @@ -0,0 +1,31 @@ +30 +1 -1 14 16 +2 18 25 21 +3 5 4 -1 +4 11 3 18 +5 3 7 18 +6 -1 27 26 +7 5 -1 23 +8 20 15 12 +9 24 10 19 +10 9 -1 16 +11 -1 4 17 +12 21 8 17 +13 -1 17 15 +14 -1 15 1 +15 13 8 14 +16 10 1 20 +17 13 11 12 +18 5 2 4 +19 9 20 22 +20 16 8 19 +21 2 22 12 +22 19 21 24 +23 7 26 25 +24 9 22 30 +25 2 23 28 +26 -1 6 23 +27 29 28 6 +28 25 27 30 +29 -1 30 27 +30 24 28 29 diff --git a/demo_backup.cpp b/demo_backup.cpp new file mode 100644 index 0000000..5758dd6 --- /dev/null +++ b/demo_backup.cpp @@ -0,0 +1,96 @@ +#include "delaunay_backup.h" +#include "iostream" +#include "fstream" +#include "iomanip" + +int main(int argc, char const *argv[]) +{ + /* + std::vector points(21); + points[0].set(-0.8, -0.8, 0); + points[1].set(0.4, -1.2, 1); + points[2].set(1.2, -0.9, 2); + points[3].set(1.6, 0.1, 3); + points[4].set(2.5, 0.5, 4); + points[5].set(4.1, 0.7, 5); + points[6].set(5.7, 1.8, 6); + points[7].set(5.1, 3.4, 7); + points[8].set(2.5, 4.4, 8); + points[9].set(1.2, 3.7, 9); + points[10].set(-1.2, 3.9, 10); + points[11].set(-3.2, 5.1, 11); + points[12].set(-4.3, 2.9, 12); + points[13].set(-3.1, 0.7, 13); + points[14].set(-1.3, 0.6, 14); + points[15].set(-2.1, 2.9, 15); + points[16].set(0.6, 1.2, 16); + points[17].set(0.1, 2.4, 17); + points[18].set(2.4, 2.8, 18); + points[19].set(3.5, 1.8, 19); + points[20].set(3.6, 3.1, 20); + */ + + std::vector points(21); + points[0].set(-0.8, -0.8, 4); + points[1].set(0.4, -1.2, 5); + points[2].set(1.2, -0.9, 6); + points[3].set(1.6, 0.1, 7); + points[4].set(2.5, 0.5, 8); + points[5].set(4.1, 0.7, 9); + points[6].set(5.7, 1.8, 10); + points[7].set(5.1, 3.4, 11); + points[8].set(2.5, 4.4, 12); + points[9].set(1.2, 3.7, 13); + points[10].set(-1.2, 3.9, 14); + points[11].set(-3.2, 5.1, 15); + points[12].set(-4.3, 2.9, 16); + points[13].set(-3.1, 0.7, 17); + points[14].set(-1.3, 0.6, 18); + points[15].set(-2.1, 2.9, 19); + points[16].set(0.6, 1.2, 20); + points[17].set(0.1, 2.4, 21); + points[18].set(2.4, 2.8, 22); + points[19].set(3.5, 1.8, 23); + points[20].set(3.6, 3.1, 24); + + if (duplicated_vertex(points)) + { + std::cerr << "Duplicated vertexes detected.\n"; + return 0; + } + + std::vector elements; + triangulation(points, elements, atoi(argv[1])); + + if (fully_delaunay(elements, points)) + { + std::clog << "The triangulation is fully delaunay.\n"; + } + else std::clog << "The triangulation is not fully delaunay\n"; + + // Write a Gmsh's .msh file + std::ofstream outfile("demo.msh"); + outfile << "$MeshFormat" << std::endl << "2.2 0 8" << std::endl << "$EndMeshFormat "<id + 1; + } + outfile << std::endl; + } + outfile << "$EndElements"<< std::endl; + outfile.close(); + return 0; +} + + diff --git a/test.msh b/test.msh new file mode 100644 index 0000000..2270361 --- /dev/null +++ b/test.msh @@ -0,0 +1,78 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +25 +1 -9.300000000000001 -8.050000000000001 0.0 +2 10.7 -8.050000000000001 0.0 +3 10.7 11.95 0.0 +4 -9.300000000000001 11.95 0.0 +5 -0.8 -0.8 0.0 +6 0.4 -1.2 0.0 +7 1.2 -0.9 0.0 +8 1.6 0.1 0.0 +9 2.5 0.5 0.0 +10 4.1 0.7 0.0 +11 5.7 1.8 0.0 +12 5.1 3.4 0.0 +13 2.5 4.4 0.0 +14 1.2 3.7 0.0 +15 -1.2 3.9 0.0 +16 -3.2 5.1 0.0 +17 -4.3 2.9 0.0 +18 -3.1 0.7 0.0 +19 -1.3 0.6 0.0 +20 -2.1 2.9 0.0 +21 0.6 1.2 0.0 +22 0.1 2.4 0.0 +23 2.4 2.8 0.0 +24 3.5 1.8 0.0 +25 3.6 3.1 0.0 +$EndNodes +$Elements +44 +1 2 0 1 2 6 +2 2 0 7 8 6 +3 2 0 5 1 6 +4 2 0 6 2 7 +5 2 0 5 6 8 +6 2 0 7 2 10 +7 2 0 8 7 9 +8 2 0 11 12 10 +9 2 0 9 7 10 +10 2 0 2 3 11 +11 2 0 10 2 11 +12 2 0 3 4 13 +13 2 0 11 3 12 +14 2 0 12 3 13 +15 2 0 8 9 23 +16 2 0 13 4 16 +17 2 0 14 13 15 +18 2 0 17 18 16 +19 2 0 15 13 16 +20 2 0 4 1 17 +21 2 0 16 4 17 +22 2 0 1 5 18 +23 2 0 5 19 18 +24 2 0 17 1 18 +25 2 0 20 18 19 +26 2 0 5 8 21 +27 2 0 15 16 20 +28 2 0 16 18 20 +29 2 0 19 21 20 +30 2 0 19 5 21 +31 2 0 8 14 21 +32 2 0 22 20 21 +33 2 0 21 14 22 +34 2 0 14 15 22 +35 2 0 15 20 22 +36 2 0 9 10 24 +37 2 0 13 14 23 +38 2 0 14 8 23 +39 2 0 23 9 24 +40 2 0 10 12 24 +41 2 0 25 23 24 +42 2 0 13 23 25 +43 2 0 24 12 25 +44 2 0 12 13 25 +$EndElements diff --git a/test_backup.msh b/test_backup.msh new file mode 100644 index 0000000..c1d3196 --- /dev/null +++ b/test_backup.msh @@ -0,0 +1,45 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +14 +1 -9.300000000000001 -8.050000000000001 0.0 +2 10.7 -8.050000000000001 0.0 +3 10.7 11.95 0.0 +4 -9.300000000000001 11.95 0.0 +5 -0.8 -0.8 0.0 +6 0.4 -1.2 0.0 +7 1.2 -0.9 0.0 +8 1.6 0.1 0.0 +9 2.5 0.5 0.0 +10 4.1 0.7 0.0 +11 5.7 1.8 0.0 +12 5.1 3.4 0.0 +13 2.5 4.4 0.0 +14 1.2 3.7 0.0 +$EndNodes +$Elements +22 +1 2 0 4 1 5 +2 2 0 1 2 6 +3 2 0 5 1 6 +4 2 0 6 2 7 +5 2 0 5 6 8 +6 2 0 6 7 8 +7 2 0 8 7 9 +8 2 0 7 2 10 +9 2 0 9 7 10 +10 2 0 2 3 11 +11 2 0 10 2 11 +12 2 0 11 3 12 +13 2 0 10 11 12 +14 2 0 3 4 13 +15 2 0 12 3 13 +16 2 0 9 10 13 +17 2 0 10 12 13 +18 2 0 4 5 14 +19 2 0 13 4 14 +20 2 0 5 8 14 +21 2 0 8 9 14 +22 2 0 9 13 14 +$EndElements