Skip to content
This repository has been archived by the owner on Nov 13, 2023. It is now read-only.

Commit

Permalink
Merge branch 'master' of https://github.com/krober10nd/SeismicMesh
Browse files Browse the repository at this point in the history
  • Loading branch information
Keith Roberts committed Feb 14, 2021
2 parents d0b123a + de58fb6 commit e9103a5
Show file tree
Hide file tree
Showing 8 changed files with 356 additions and 83 deletions.
83 changes: 83 additions & 0 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,7 @@ Table of contents
* [Intersection](#intersection)
* [Difference](#difference)
* [Immersion](#immersion)
* [Boundaries](#boundaries)
* [Periodic](#periodic)
* [Parallelism](#parallelism)
* [Performance comparison](#performance)
Expand Down Expand Up @@ -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
--------

Expand Down Expand Up @@ -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.

<img width="1221" alt="Screen Shot 2021-02-12 at 12 04 03 PM" src="https://user-images.githubusercontent.com/18619644/107784877-b1902500-6d2a-11eb-98f3-e01c1175f498.png">


```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
-------------
<img alt="Periodic torus" src="https://user-images.githubusercontent.com/18619644/101163708-bfcb1200-3612-11eb-9c6d-4f664a754d01.png" width="30%">
Expand Down Expand Up @@ -626,6 +706,9 @@ 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.

## [3.3.0] -2021-01-08
### Added
Expand Down
101 changes: 71 additions & 30 deletions SeismicMesh/generation/cpp/delaunay_class.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
#include <iterator>
#include <string>
#include <vector>
#include <ctime>

#include <boost/lexical_cast.hpp>

Expand All @@ -26,7 +27,9 @@ using DT = CGAL::Delaunay_triangulation_2<K, Tds>;

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 <typename T> class TypedInputIterator {
public:
Expand Down Expand Up @@ -61,6 +64,23 @@ template <typename T> 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_<Point>(m, "Point")
.def(py::init<int, int>(), py::arg("x"), py::arg("y"))
Expand Down Expand Up @@ -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<Face_handle>.
// 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<int> 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<ssize_t> shape = {num_edges, 2};
std::vector<ssize_t> 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<int>::format(), /* data type */
2, /* number of dimensions */
shape, /* shape of the matrix */
strides /* strides for each axis */
));
[](const DT &dt, const std::vector<int> &faces_to_get){

// User supplies face_to_get
std::unordered_set<Face_handle> 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<int> 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<ssize_t> shape = {num_edges, 2};
std::vector<ssize_t> 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<int>::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)
Expand Down
113 changes: 113 additions & 0 deletions SeismicMesh/generation/cpp/delaunay_class3.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,32 @@ using DT = CGAL::Delaunay_triangulation_3<K, Tds>;

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 <typename T>
std::vector<T> vectorSortIntArr(std::vector<std::array<T, 2>> 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<T> 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 <typename T> class TypedInputIterator {
public:
Expand Down Expand Up @@ -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<int> &cells_to_get){

// User supplies cells_to_get
std::unordered_set<Cell_handle> 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<int> 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<ssize_t> shape = {num_edges, 2};
std::vector<ssize_t> 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<int>::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;
Expand Down
2 changes: 1 addition & 1 deletion SeismicMesh/sizing/mesh_size_function.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
Loading

0 comments on commit e9103a5

Please sign in to comment.