From 44e46d131de06ceb92cc8fa887101ecf187aafa7 Mon Sep 17 00:00:00 2001 From: Keith Roberts Date: Tue, 9 Feb 2021 07:04:46 -0300 Subject: [PATCH 1/3] Faster unique edges (#182) * faster unique edges Co-authored-by: Keith Roberts --- README.md | 2 + SeismicMesh/generation/cpp/delaunay_class.cpp | 101 +++++++++++----- .../generation/cpp/delaunay_class3.cpp | 113 ++++++++++++++++++ benchmarks/run_BP2004.py | 38 +++--- benchmarks/run_EAGE.py | 37 ++++-- benchmarks/run_ball.py | 31 +++-- benchmarks/run_disk.py | 34 +++--- 7 files changed, 274 insertions(+), 82 deletions(-) diff --git a/README.md b/README.md index 0fa159a6..ff955c83 100644 --- a/README.md +++ b/README.md @@ -626,6 +626,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Mesh improvement now solves Lapl. smoothing as a fixed-point problem using AMG solver. - User can now mesh user-defined sizing functions in parallel (not from :class:SizeFunction) - Ability to specify data type `dtype` of floating point number inside binary files. +### Improved +- Faster unique edge calculation. ## [3.3.0] -2021-01-08 ### Added diff --git a/SeismicMesh/generation/cpp/delaunay_class.cpp b/SeismicMesh/generation/cpp/delaunay_class.cpp index b24ed612..491f5d57 100644 --- a/SeismicMesh/generation/cpp/delaunay_class.cpp +++ b/SeismicMesh/generation/cpp/delaunay_class.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include @@ -26,7 +27,9 @@ using DT = CGAL::Delaunay_triangulation_2; using Point = K::Point_2; using Vertex_handle = DT::Vertex_handle; +using Face_handle = DT::Face_handle; using Vi = DT::Finite_vertices_iterator; +using Fi = DT::Finite_faces_iterator; template class TypedInputIterator { public: @@ -61,6 +64,23 @@ template class TypedInputIterator { py::iterator py_iter_; }; + +class Timer { +public: + Timer() { clock_gettime(CLOCK_REALTIME, &beg_); } + + double elapsed() { + clock_gettime(CLOCK_REALTIME, &end_); + return end_.tv_sec - beg_.tv_sec + + (end_.tv_nsec - beg_.tv_nsec) / 1000000000.; + } + + void reset() { clock_gettime(CLOCK_REALTIME, &beg_); } + +private: + timespec beg_, end_; +}; + PYBIND11_MODULE(delaunay_class, m) { py::class_(m, "Point") .def(py::init(), py::arg("x"), py::arg("y")) @@ -137,37 +157,58 @@ PYBIND11_MODULE(delaunay_class, m) { return dt; }) - // get all the finite edges of the triangulation + // Put the face handles of these faces in a std::set. + // Then for each face f in the list look at the three edges (f,i) for i = 0,1,2. + // For each edge get the opposite face, and check if it is in the set. + // If it is in the set report the unique edge by comparing the &* of the two face handles. If it is not report the edge e. .def("get_edges", - [](DT &dt){ - std::vector edges; - - for(DT::Edge_iterator ei=dt.finite_edges_begin();ei!=dt.finite_edges_end(); ei++){ - // Get a vertex from the edge - DT::Face& f = *(ei->first); - int i = ei->second; - Vertex_handle vs = f.vertex(f.cw(i)); - Vertex_handle vt = f.vertex(f.ccw(i)); - - edges.push_back(vs->info()); - edges.push_back(vt->info()); - //std::cout << vs->info() << vt->info() << std::endl; - } - ssize_t soint = sizeof(int); - ssize_t num_edges = edges.size() / 2; - ssize_t ndim = 2; - std::vector shape = {num_edges, 2}; - std::vector strides = {soint * 2, soint}; - - // return 2-D NumPy array - return py::array(py::buffer_info( - edges.data(), /* data as contiguous array */ - sizeof(int), /* size of one scalar */ - py::format_descriptor::format(), /* data type */ - 2, /* number of dimensions */ - shape, /* shape of the matrix */ - strides /* strides for each axis */ - )); + [](const DT &dt, const std::vector &faces_to_get){ + + // User supplies face_to_get + std::unordered_set face_handles; + face_handles.reserve(faces_to_get.size()); + int ix = 0; + for (Fi fi = dt.finite_faces_begin(); fi != dt.finite_faces_end(); fi++) { + if(faces_to_get[ix] == 1){ + face_handles.insert(fi); + } + ix += 1; + } + std::vector edges; + edges.reserve(3*faces_to_get.size()); + // If it is in the set report the unique edge by comparing the &* of the two face handle + int n = 0; + for (const auto& face : face_handles){ + // for each face look at its three edges + for (std::size_t i = 0; i < 3; ++i){ + Face_handle nei_face = face->neighbor(i); + if (face > nei_face && face_handles.count(nei_face)) continue; // report interior edges once + + Vertex_handle vs = face->vertex(face->cw(i)); + Vertex_handle vt = face->vertex(face->ccw(i)); + + edges.push_back(vs->info()); + edges.push_back(vt->info()); + n += 2; + } + } + edges.resize(n); + + ssize_t soint = sizeof(int); + ssize_t num_edges = edges.size() / 2; + ssize_t ndim = 2; + std::vector shape = {num_edges, 2}; + std::vector strides = {soint * 2, soint}; + + // return 2-D NumPy array + return py::array(py::buffer_info( + edges.data(), /* data as contiguous array */ + sizeof(int), /* size of one scalar */ + py::format_descriptor::format(), /* data type */ + 2, /* number of dimensions */ + shape, /* shape of the matrix */ + strides /* strides for each axis */ + )); }) .def("number_of_vertices", &DT::number_of_vertices) diff --git a/SeismicMesh/generation/cpp/delaunay_class3.cpp b/SeismicMesh/generation/cpp/delaunay_class3.cpp index 0ea1b3ae..106a8a47 100644 --- a/SeismicMesh/generation/cpp/delaunay_class3.cpp +++ b/SeismicMesh/generation/cpp/delaunay_class3.cpp @@ -25,7 +25,32 @@ using DT = CGAL::Delaunay_triangulation_3; using Point = K::Point_3; using Vertex_handle = DT::Vertex_handle; +using Cell_handle = DT::Cell_handle; +using Edge = DT::Edge; using Vi = DT::Finite_vertices_iterator; +using Cc = DT::Cell_circulator; +using Ci = DT::Finite_cells_iterator; + +template + std::vector vectorSortIntArr(std::vector> v) { + std::sort(v.begin(), v.end()); + // double t = tmr.elapsed(); + // tmr.reset(); + auto iter = std::unique(v.begin(), v.end()); + // t = tmr.elapsed(); + // std::cout << t << std::endl; + + size_t len = iter - v.begin(); + std::vector outvec; + outvec.reserve(len * 2); + for (auto i = v.begin(); i != iter; ++i) { + outvec.push_back(i->at(0)); + outvec.push_back(i->at(1)); + } + return outvec; + } + + template class TypedInputIterator { public: @@ -144,6 +169,94 @@ PYBIND11_MODULE(delaunay_class3, m) { .def("number_of_vertices", [](DT &dt) { return dt.number_of_vertices(); }) + .def("get_edges",[](const DT &dt, const std::vector &cells_to_get){ + + // User supplies cells_to_get + std::unordered_set cell_handles; + cell_handles.reserve(cells_to_get.size()); + int ix = 0; + for (Ci ci = dt.finite_cells_begin(); ci != dt.finite_cells_end(); ci++) { + if(cells_to_get[ix] == 1){ + cell_handles.insert(ci); + } + ix += 1; + } + std::vector edges; + edges.reserve(6*cells_to_get.size()); + // this creates each edge of the tetrahedron + int indices[6][2]; + // first edge + indices[0][0]=0; + indices[0][1]=1; + // second edge + indices[1][0]=1; + indices[1][1]=2; + // third edge + indices[2][0]=2; + indices[2][1]=0; + // fourth edge + indices[3][0]=0; + indices[3][1]=3; + // fifth edge + indices[4][0]=1; + indices[4][1]=3; + // sixth edge + indices[5][0]=2; + indices[5][1]=3; + // If it is in the set report the unique edge by comparing the &* of the two face handle + int n = 0; + for (const auto& cell : cell_handles){ + // for each cell look at its six edges + for (std::size_t i = 0; i < 6; ++i){ + // build edge structure + int indx0 = indices[i][0]; + int indx1 = indices[i][1]; + auto edge = Edge(cell, indx0, indx1); + // circulate around the edge + Cc nei_cells = dt.incident_cells(edge, cell); + Cc done = nei_cells; + bool is_smallest = true; + do + { + Cell_handle nei_cell = nei_cells; + if (cell > nei_cell && cell_handles.count(nei_cell)) { + ++nei_cells; // advance to next neighoring cell + is_smallest = false; + break; // report interior edges once + } + ++nei_cells; + } + while(nei_cells != done); + + if (is_smallest) + { + Vertex_handle vs = cell->vertex(indx0); + Vertex_handle vt = cell->vertex(indx1); + + edges.push_back(vs->info()); + edges.push_back(vt->info()); + n += 2; + } + } + } + edges.resize(n); + ssize_t soint = sizeof(int); + ssize_t num_edges = edges.size() / 2; + ssize_t ndim = 2; + std::vector shape = {num_edges, 2}; + std::vector strides = {soint * 2, soint}; + + // return 2-D NumPy array + return py::array(py::buffer_info( + edges.data(), /* data as contiguous array */ + sizeof(int), /* size of one scalar */ + py::format_descriptor::format(), /* data type */ + 2, /* number of dimensions */ + shape, /* shape of the matrix */ + strides /* strides for each axis */ + )); + }) + .def("number_of_cells", [](DT &dt) { int count = 0; diff --git a/benchmarks/run_BP2004.py b/benchmarks/run_BP2004.py index f7c7c7ac..e8771c5e 100644 --- a/benchmarks/run_BP2004.py +++ b/benchmarks/run_BP2004.py @@ -1,14 +1,16 @@ # run BP2004 benchmark and plot a figure import matplotlib.pyplot as plt +import matplotlib import numpy from benchmark_BP2004 import run_gmsh, run_SeismicMesh, _build_sizing +matplotlib.use("TkAgg") -colors1 = ["ko-", "ro-"] -colors2 = ["ko--", "ro--"] -labels = ["gmsh", "SeismicMesh"] +colors1 = ["C7o-", "C3o-"] +colors2 = ["C7o--", "C3o--"] +labels = ["Gmsh", "SeismicMesh"] -plt.rcParams.update({"font.size": 18}) +# plt.rcParams.update({"font.size": 18}) entries = [] # minimize mesh size @@ -39,20 +41,28 @@ entries.append(h) plt.subplot(1, 2, 1) -plt.title("Number of cells vs. mesh generation time") +plt.title("# of cells vs. mesh generation time [s]") plt.legend() -plt.xlabel("Number of cells") -plt.ylabel("Elapsed time (s)") +plt.xlabel("# of cells") +# plt.ylabel("Elapsed time (s)") plt.ticklabel_format(style="sci", axis="x", scilimits=(0, 0)) -plt.grid() +axes = plt.gca() +axes.yaxis.grid() +axes.set_frame_on(False) plt.subplot(1, 2, 2) -plt.title("Number of cells vs. cell quality") -plt.xlabel("Number of cells") -plt.ylabel("Cell quality") -plt.xlabel("Number of cells") +plt.title("# of cells vs. cell quality") +plt.xlabel("# of cells") +# plt.ylabel("Cell quality") +plt.xlabel("# of cells") plt.ylim(ymax=1, ymin=0) plt.ticklabel_format(style="sci", axis="x", scilimits=(0, 0)) -plt.grid() +axes = plt.gca() +axes.yaxis.grid() +axes.set_frame_on(False) -plt.show() + +# plt.show() + +plt.tight_layout(pad=3.0) +plt.savefig("out2.svg", transparent=True, bbox_inches="tight", pad_inches=0) diff --git a/benchmarks/run_EAGE.py b/benchmarks/run_EAGE.py index 6f112f88..0609a2d8 100644 --- a/benchmarks/run_EAGE.py +++ b/benchmarks/run_EAGE.py @@ -1,14 +1,16 @@ # run EAGE benchmark and plot a figure +import matplotlib import matplotlib.pyplot as plt import numpy from benchmark_EAGE import run_gmsh, run_SeismicMesh, _build_sizing +matplotlib.use("TkAgg") -colors1 = ["ko-", "ro-", "bo-"] -colors2 = ["ko--", "ro--", "bo--"] -labels = ["gmsh", "SeismicMesh", "cgal"] +colors1 = ["C7o-", "C3o-", "C0o-"] +colors2 = ["C7o--", "C3o--", "C0o--"] +labels = ["Gmsh", "SeismicMesh", "CGAL"] -plt.rcParams.update({"font.size": 18}) +# plt.rcParams.update({"font.size": 18}) entries = [] # minimal mesh size @@ -44,18 +46,27 @@ entries.append(h) plt.subplot(1, 2, 1) -plt.title("Number of cells vs. mesh generation time") +plt.title("# of cells vs. mesh generation time [s]") plt.ticklabel_format(style="sci", axis="x", scilimits=(0, 0)) plt.legend() -plt.xlabel("Number of cells") -plt.ylabel("Elapsed time (s)") -plt.grid() +plt.xlabel("# of cells") +plt.ylim(ymax=60, ymin=0) + +axes = plt.gca() +axes.yaxis.grid() +axes.set_frame_on(False) + plt.subplot(1, 2, 2) -plt.title("Number of cells. vs cell quality") +plt.title("# of cells. vs cell quality") plt.ticklabel_format(style="sci", axis="x", scilimits=(0, 0)) -plt.xlabel("Number of cells") -plt.ylabel("Cell quality") -plt.grid() +plt.xlabel("# of cells") +plt.ylim(ymax=1, ymin=0) + +axes = plt.gca() +axes.yaxis.grid() +axes.set_frame_on(False) + -plt.show() +plt.tight_layout(pad=3.0) +plt.savefig("out4.svg", transparent=True, bbox_inches="tight", pad_inches=0) diff --git a/benchmarks/run_ball.py b/benchmarks/run_ball.py index 186d9b13..8aa1018a 100644 --- a/benchmarks/run_ball.py +++ b/benchmarks/run_ball.py @@ -1,13 +1,16 @@ # run ball benchmark and plot a figure +import matplotlib import matplotlib.pyplot as plt import numpy from benchmark_ball import run_gmsh, run_SeismicMesh, run_cgal -plt.rcParams.update({"font.size": 18}) -colors1 = ["ko-", "ro-", "bo-"] -colors2 = ["ko--", "ro--", "bo--"] -labels = ["gmsh", "SeismicMesh", "cgal"] +matplotlib.use("TkAgg") +# plt.rcParams.update({"font.size": 18}) + +colors1 = ["C7o-", "C3o-", "C0o-"] +colors2 = ["C7o--", "C3o--", "C0o--"] +labels = ["Gmsh", "SeismicMesh", "CGAL"] entries = [] # minimal mesh size @@ -36,18 +39,24 @@ entries.append(h) plt.subplot(1, 2, 1) -plt.title("Number of cells vs. mesh generation time") +plt.title("# of cells vs. mesh generation time [s]") plt.legend() -plt.xlabel("Number of cells") -plt.ylabel("Elapsed time (s)") +plt.xlabel("# of cells") plt.ticklabel_format(style="sci", axis="x", scilimits=(0, 0)) -plt.grid() +axes = plt.gca() +axes.yaxis.grid() +plt.ylim(ymax=19, ymin=0) +axes.set_frame_on(False) plt.subplot(1, 2, 2) plt.title("Number of cells. vs cell quality") plt.ticklabel_format(style="sci", axis="x", scilimits=(0, 0)) plt.xlabel("Number of cells") -plt.ylabel("Cell quality") -plt.grid() +plt.ylim(ymax=1, ymin=0) +axes = plt.gca() +axes.yaxis.grid() +axes.set_frame_on(False) + -plt.show() +plt.tight_layout(pad=3.0) +plt.savefig("out3.svg", transparent=True, bbox_inches="tight", pad_inches=0) diff --git a/benchmarks/run_disk.py b/benchmarks/run_disk.py index c561c695..c8e6019a 100644 --- a/benchmarks/run_disk.py +++ b/benchmarks/run_disk.py @@ -1,15 +1,18 @@ # run disk benchmark and plot a figure +import matplotlib import matplotlib.pyplot as plt + import numpy from benchmark_disk import run_gmsh, run_SeismicMesh, run_cgal +matplotlib.use("TkAgg") -plt.rcParams.update({"font.size": 19}) +# plt.rcParams.update({"font.size": 14}) -colors1 = ["ko-", "ro-", "bo-"] -colors2 = ["ko--", "ro--", "bo--"] -labels = ["gmsh", "SeismicMesh", "cgal"] +colors1 = ["C7o-", "C3o-", "C0o-"] +colors2 = ["C7o--", "C3o--", "C0o--"] +labels = ["Gmsh", "SeismicMesh", "CGAL"] entries = [] # minimize mesh size @@ -38,20 +41,23 @@ entries.append(h) plt.subplot(1, 2, 1) -plt.title("Number of cells vs. mesh generation time") +plt.title("# of cells vs. mesh generation time [s]") plt.legend() -plt.xlabel("Number of cells") +plt.xlabel("# of cells") plt.ticklabel_format(style="sci", axis="x", scilimits=(0, 0)) -plt.ylabel("Elapsed time (s)") -plt.grid() +axes = plt.gca() +axes.yaxis.grid() +axes.set_frame_on(False) plt.subplot(1, 2, 2) -plt.title("Number of cells. vs cell quality") -plt.xlabel("Number of cells") +plt.title("# of cells. vs cell quality") +plt.xlabel("# of cells") plt.ticklabel_format(style="sci", axis="x", scilimits=(0, 0)) -plt.ylim(ymax=1, ymin=0) -plt.ylabel("Cell quality") +plt.ylim(ymax=1.1, ymin=0) +axes = plt.gca() +axes.yaxis.grid() +axes.set_frame_on(False) -plt.grid() -plt.show() +plt.tight_layout(pad=3.0) +plt.savefig("out.svg", transparent=True, bbox_inches="tight", pad_inches=0) From 8751bad903b1a81dd4ee83eb7efec542bfa369f7 Mon Sep 17 00:00:00 2001 From: Keith Roberts Date: Fri, 12 Feb 2021 13:01:37 -0300 Subject: [PATCH 2/3] Adding another example (#191) * adding example specifying 2d boundaries Co-authored-by: Keith Roberts --- README.md | 81 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) diff --git a/README.md b/README.md index ff955c83..0c302380 100644 --- a/README.md +++ b/README.md @@ -41,6 +41,7 @@ Table of contents * [Intersection](#intersection) * [Difference](#difference) * [Immersion](#immersion) + * [Boundaries](#boundaries) * [Periodic](#periodic) * [Parallelism](#parallelism) * [Performance comparison](#performance) @@ -277,6 +278,7 @@ if comm.rank == 0: **The user can still specify their own signed distance functions and sizing functions to `generate_mesh` (in serial or parallel) just like the original DistMesh algorithm but now with quality bounds in 3D. Try the codes below!** + Cylinder -------- @@ -555,6 +557,84 @@ meshio.write_points_cells( ) ``` +Boundaries +----------- + +Boundary conditions can also be prescribed and written to `gmsh` compatible files using `mehsio`. In the following example, we immerse a disk into the connectivity and then prescribe boundary conditions around the circle and each wall of the domain for later usage inside a finite element solver. + +Screen Shot 2021-02-12 at 12 04 03 PM + + +```python +import numpy as np +import meshio +import SeismicMesh as sm + +bbox = (0.0, 10.0, 0.0, 1.0) +channel = sm.Rectangle(bbox) +suspension = sm.Disk([0.5, 0.5], 0.25) + +hmin = 0.10 +fh = lambda p: 0.05 * np.abs(suspension.eval(p)) + hmin +points, cells = sm.generate_mesh( + domain=channel, + edge_length=fh, + h0=hmin, + subdomains=[suspension], + max_iter=1000, + ) +# This gets the edges of the mesh in a winding order (clockwise or counterclockwise). +ordered_bnde = sm.geometry.get_winded_boundary_edges(cells) +# We use the midpoint of the edge to determine its boundary label +mdpt = points[ordered_bnde].sum(1) / 2 +infl = ordered_bnde[mdpt[:, 0] < 1e-6, :] # x=0.0 +outfl = ordered_bnde[mdpt[:, 0] > 9.9 + 1e-6, :] # x=10.0 +walls = ordered_bnde[ + (mdpt[:, 1] < 1e-6) | (mdpt[:, 1] > 0.99 + 1e-6), : +] # y=0.0 or y=1.0 +cells_prune = cells[suspension.eval(sm.geometry.get_centroids(points, cells)) < 0] +circle = sm.geometry.get_winded_boundary_edges(cells_prune) + +# Write to gmsh22 format with boundary conditions for the walls and disk/circle. +meshio.write_points_cells( + "example.msh", + points, + cells=[ + ("triangle", cells), + ("line", np.array(infl)), + ("line", np.array(outfl)), + ("line", np.array(walls)), + ("line", np.array(circle)), + ], + field_data={ + "InFlow": np.array([11, 1]), + "OutFlow": np.array([12, 1]), + "Walls": np.array([13, 1]), + "Circle": np.array([14, 1]), + }, + cell_data={ + "gmsh:physical": [ + np.repeat(3, len(cells)), + np.repeat(11, len(infl)), + np.repeat(12, len(outfl)), + np.repeat(13, len(walls)), + np.repeat(14, len(circle)), + ], + "gmsh:geometrical": [ + np.repeat(1, len(cells)), + np.repeat(1, len(infl)), + np.repeat(1, len(outfl)), + np.repeat(1, len(walls)), + np.repeat(1, len(circle)), + ], + }, + file_format="gmsh22", + binary=False, +) +``` + + + Periodic ------------- Periodic torus @@ -626,6 +706,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Mesh improvement now solves Lapl. smoothing as a fixed-point problem using AMG solver. - User can now mesh user-defined sizing functions in parallel (not from :class:SizeFunction) - Ability to specify data type `dtype` of floating point number inside binary files. +- Example how to specify and write boundary conditions. ### Improved - Faster unique edge calculation. From de58fb6e5c015c2d065d736d5fa6f7e2caf15bef Mon Sep 17 00:00:00 2001 From: Keith Roberts Date: Sat, 13 Feb 2021 16:04:11 -0300 Subject: [PATCH 3/3] typo in sizing function --- SeismicMesh/sizing/mesh_size_function.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SeismicMesh/sizing/mesh_size_function.py b/SeismicMesh/sizing/mesh_size_function.py index 9ee22616..04c244f6 100644 --- a/SeismicMesh/sizing/mesh_size_function.py +++ b/SeismicMesh/sizing/mesh_size_function.py @@ -166,7 +166,7 @@ def get_sizing_function_from_segy(filename, bbox, comm=None, **kwargs): "byte_order", "axes_order", "axes_order_sort", - "dytpe", + "dtype", }: pass else: