search

meshpro / meshplex

:spider_web: Compute interesting points, areas, and volumes in simplex meshes of any dimension.
104 stars 22 forks source link

Usage Example and MeshTetra.show() Borked (And We Have No Idea How to Use voropy) #22

Closed leycec closed 5 years ago

leycec commented 7 years ago

The Usage example in README.md appears to be fundamentally broken on newer voropy versions – and should instead resemble:

mesh, point_data, cell_data, field_data = voropy.read('pacman.msh')

print(mesh.node_coords)
print(mesh.get_control_volumes())

mesh.show()

Likewise, the MeshTetra.show() method appears to be equally broken. For example, given a three-dimensional tetrahedral mesh exported by netgen to Gmsh format:

>>> mesh, point_data, cell_data, field_data = voropy.read("poppy.msh")
>>> mesh.show()
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-57-cbb52bec63cc> in <module>()
----> 1 mesh.show()

/usr/local/lib/python3.5/dist-packages/voropy/mesh_tetra.py in show(self)
    462                     [cc[..., 1], face_cc[..., 1]],
    463                     [cc[..., 2], face_cc[..., 2]],
--> 464                     'b-'
    465                     )
    466         return

/usr/lib/python3/dist-packages/mpl_toolkits/mplot3d/axes3d.py in plot(self, xs, ys, *args, **kwargs)
   1539             zs = np.ones(len(xs)) * zs
   1540 
-> 1541         lines = Axes.plot(self, xs, ys, *args[argsi:], **kwargs)
   1542         for line in lines:
   1543             art3d.line_2d_to_3d(line, zs=zs, zdir=zdir)

/usr/lib/python3/dist-packages/matplotlib/__init__.py in inner(ax, *args, **kwargs)
   1812                     warnings.warn(msg % (label_namer, func.__name__),
   1813                                   RuntimeWarning, stacklevel=2)
-> 1814             return func(ax, *args, **kwargs)
   1815         pre_doc = inner.__doc__
   1816         if pre_doc is None:

/usr/lib/python3/dist-packages/matplotlib/axes/_axes.py in plot(self, *args, **kwargs)
   1422             kwargs['color'] = c
   1423 
-> 1424         for line in self._get_lines(*args, **kwargs):
   1425             self.add_line(line)
   1426             lines.append(line)

/usr/lib/python3/dist-packages/matplotlib/axes/_base.py in _grab_next_args(self, *args, **kwargs)
    384                 return
    385             if len(remaining) <= 3:
--> 386                 for seg in self._plot_args(remaining, kwargs):
    387                     yield seg
    388                 return

/usr/lib/python3/dist-packages/matplotlib/axes/_base.py in _plot_args(self, tup, kwargs)
    357 
    358         if len(tup) == 2:
--> 359             x = _check_1d(tup[0])
    360             y = _check_1d(tup[-1])
    361         else:

/usr/lib/python3/dist-packages/matplotlib/cbook.py in _check_1d(x)
   2205     '''
   2206     if not hasattr(x, 'shape') or len(x.shape) < 1:
-> 2207         return np.atleast_1d(x)
   2208     else:
   2209         try:

/usr/lib/python3/dist-packages/numpy/core/shape_base.py in atleast_1d(*arys)
     48     res = []
     49     for ary in arys:
---> 50         ary = asanyarray(ary)
     51         if len(ary.shape) == 0:
     52             result = ary.reshape(1)

/usr/lib/python3/dist-packages/numpy/core/numeric.py in asanyarray(a, dtype, order)
    531 
    532     """
--> 533     return array(a, dtype, copy=False, order=order, subok=True)
    534 
    535 def ascontiguousarray(a, dtype=None):

ValueError: setting an array element with a sequence.

These are the least of our concerns, however. My wife and I were contemplating integrating voropy into our downstream finite volume multiphysics biology simulator for generation of three-dimensional Voronoi diagrams, but frankly have no idea how to actually use voropy for that purpose.

voropy's documentation is sparse at best (and where documented is frequently obsolete, which is almost worse), the codebase itself is largely uncommented, the usage examples mostly pertain to two-dimensional finite volume discretizations requiring the equally undocumented pyfvm, and there doesn't actually appear to be any support for three-dimensional Voronoi tesselations yet. ...not that we could grep, anyway.

voropy currently appears to wrap three-dimensional meshio return values with instances of the slightly higher-level MeshTetra class, which only appears to encapsulate various geometric properties of the input mesh rather than generate a three-dimensional Voronoi tesselation. Maybe? Or maybe we have no idea what we're doing here.

Ignoring our obvious interest in three-dimensional Voronoi tesselations, we'd settle for merely being able to access any (ideally, all) of the following MeshTetra properties:

Documentation for existing MeshTetra properties – particularly the idx_hierarchy and node_face_cells Numpy arrays, which we suspect would probably be of interest to us if we actually understood them – would also be awesome.

Our working assumption is that voropy is currently optimized for two-dimensional meshes and Voronoi tesselations. If so, we'll probably just unroll our own application-specific Voronoi implementation and revisit voropy in a few years... after the dust has settled, fully works, and is well-documented.

In short, it's been a frustrating Sunday. Hey, why are we working on Sunday, anyway!? </sigh>

nschloe commented 7 years ago

Ah, I didn't anyone except myself was using voropy to some extend.

You're right: The examples in the readme are broken. I'll fix that tomorrow and will also include those snippets in the nightly tests so this doesn't happen again.

the codebase itself is largely uncommented,

That's a bit harsh. There's even beautiful ascii art! :) https://github.com/nschloe/voropy/blob/master/voropy/mesh_tri.py

Face normals. Face surface areas. Neighbouring faces. Neighbouring tetrahedrons.

I'll look at it this week and see if this info could be made easily accessible.

Documentation for existing MeshTetra properties – particularly the idx_hierarchy and node_face_cells Numpy arrays, which we suspect would probably be of interest to us if we actually understood them – would also be awesome.

Right, right. I'll make sure to add a few lines of documentation this week.

Early spoiler: Input meshes typically consist of a bunch of points with their locations, and simplices given by their corner points (sometimes call the element connectivity table). In particular there are no edges and faces. One can of course compute those with the given information, but this is costly and almost never required. Instead, faces are treated local to a cell, edges local a face. This constitutes the idx_hierarchy array:

cells -> faces -> edges -> nodes

In 3D, with tets, It has shape [n, 4, 3, 2], so to get the 1st node of the 3rd edge of the 2nd face of cell number 377, you go idx_hierarchy[377, 2, 3, 1].

pietakio commented 7 years ago

Thank you for the rapid reply! And thanks for providing your work with voropy -- it will be super useful. By the way, the ASCII art is impressive, and we got fairly far just looking into the code (which you're right, is pretty well commented).

The explanation you give for the idx_hierarchy will help a lot. Thank you.

Looking at older documentation (0.1.6) it seemed like Voropy's MeshTetra used to support derivation of face normals...but I can't find that in the module or base presently. Having access to face normal vectors would be very helpful.

Also, for MeshTetra's show() method, if you're fixing it (I wasn't able to get it working), you might want to try matplotlib's Poly3DCollection, which may plot the information faster than the repeated plot command.

Finally, are we correct in the interpretation that the MeshTetra isn't actually calculating a 3D Voronoi tessellation, but rather, is providing useful information for constructing finite volume method calculations on a tetrahedral mesh (e.g. the idx heirarchy, tetrahedral volumes, centroids, and hopefully surface areas, cell neighbours, and normal vectors for the faces?).

Thank you for your help with this!

nschloe commented 7 years ago

Having access to face normal vectors would be very helpful.

It is one defining "feature" of voropy that it starts off from a simplex mesh, and creates the dual voronoi mesh from that. This sets it apart from many other voronoi codes which take a point set with some notion of domain boundary as input.

From the presence of the dual mesh it follows that some entities are readily available, for example the face normals of the Voronoi formulation: They are simply the edges in the simplex mesh.

Also, for MeshTetra's show() method, if you're fixing it (I wasn't able to get it working), you might want to try matplotlib's Poly3DCollection, which may plot the information faster than the repeated plot command.

Thanks for the hint! The Voronoi mesh representation might actually be one of the hardest things to get in order. Even if all the intersection points can be computed, I'm not sure how to display Voronoi regions in 3D without messing up the picture.

Finally, are we correct in the interpretation that the MeshTetra isn't actually calculating a 3D Voronoi tessellation, but rather, is providing useful information for constructing finite volume method calculations on a tetrahedral mesh (e.g. the idx heirarchy, tetrahedral volumes, centroids, and hopefully surface areas, cell neighbours, and normal vectors for the faces?).

Well – what constitures a Voronoi tesselation? Given a simplex mesh, voropy mainly computes the size of the Voronoi regions associated at each point and the area of the faces between two adjacent boxes (and does so quite cleverly). For geometrical purposes, the circumcenters of the boxes and faces are of interest (which are also computed).

pietakio commented 7 years ago

Thank you for confirming that voropy is indeed calculating the 3D Voronoi dual from the simplex mesh. That's fantastic. Voropy truly has a lot of capability. I think I'm starting to get it now and see that all the information I'm asking for is indeed cleverly calculated and there. BUT I still don't see how I can make use of all that cleverly calculated Voronoi mesh information...

Now, what I'd like is to understand how the Voronoi dual data is packaged with respect to the Voronoi tesselation, not the original tetrahedral mesh.

For instance, you mention that the edges of the original graph are the normals to the faces of the Voronoi cell...OK, great, but what are the indices to the edges of the original mesh that are associated with each cell in the Voronoi dual?

From what you said, it seems that the original tetrahedral mesh is taken to be the 3D Delaunay triangulation of the Voronoi dual, so voropy is going back to reconstruct the Voronoi dual from that. Therefore, each node of the original mesh has a Voronoi dual cell built around it (analogous to 2D, which is easy to visualize with the pacman.msh example). We could further say that the index of the Voronoi cell (not the index of the original mesh tetrahedra) is the integer set 'q'. And I think the number of Voronoi cells (length of 'q') corresponds to the total number of nodes in the original mesh, right?

So what I'm looking to understand is given the Voronoi dual, how do I access:

Thanks for your help with understanding how to use voropy!

PS

Incidentally, these are some of the data structures I understand and do not understand after looking at voropy code. If I don't need to understand the following to achieve the above-mentioned goals relating to the Voronoi dual, please don't worry about explaining them -- I'm just mentioning this so you're aware of my knowledge gaps.

I'm reading in the mesh using the following command, and performing the method calls to fill in what seems to be missing information: 

mesh, point_data, cell_data, field_data = voropy.read("/home/pietakio/Documents/ForceFields/poppy.msh")

mesh.mark_boundary()
mesh.create_face_edge_relationships()
mesh._compute_cell_circumcenters()
mesh.get_control_volumes()

Now, this is what I understand the various resulting data fields of the mesh object to represent:

# The following I think I understand:
mesh._circumcenters  # the center of the Voronoi cell? (why is this variable private??)
mesh._control_volumes  # the volume of the Voronoi cell? (also, why is this private??)
mesh.cell_volumes   # the volume of the original mesh tetrahedra? 
mesh.circumcenter_face_distances # distance between Voronoi cell centre and each face of the Voronoi cell?
mesh.node_coords # the coordinates of the original mesh points
mesh.idx_hierarchy # a fantastic way to structure the tet, face, edge, node information of the original mesh
mesh.edge_coords   # we already have coords of the original mesh nodes and the idx_hierarchy, so these are...redundant?   
mesh.edges           # indices to edges of the original mesh tetrahedra? 
mesh.faces           # but reason for 'edges' and 'nodes' indexing not understood...? 

# The following I don't understand at all: 

mesh.ei_dot_ei        # what is this? 
mesh.ei_dot_ej      # what is this? 
mesh.cells   # nodes index points of the original mesh, what is "opposing vertex"...? 
mesh.node_face_cells #??? I don't get what this is packaging...or why...???
mesh.local_idx    # ????
mesh.local_idx_inv    # ????

Thanks!

nschloe commented 7 years ago

For instance, you mention that the edges of the original graph are the normals to the faces of the Voronoi cell...OK, great, but what are the indices to the edges of the original mesh that are associated with each cell in the Voronoi dual?

Edges and faces are expensive to compute which is why voropy doesn't do that (by default). You almost always only need local edges and faces; idx_hierarchy is your friend.

From what you said, it seems that the original tetrahedral mesh is taken to be the 3D Delaunay triangulation of the Voronoi dual, so voropy is going back to reconstruct the Voronoi dual from that. Therefore, each node of the original mesh has a Voronoi dual cell built around it (analogous to 2D, which is easy to visualize with the pacman.msh example). We could further say that the index of the Voronoi cell (not the index of the original mesh tetrahedra) is the integer set 'q'. And I think the number of Voronoi cells (length of 'q') corresponds to the total number of nodes in the original mesh, right?

Correct.

for an individual Voronoi cell 'q', the centroid (which I believe is the private variable _circumcentre? )

A circumcenter and a centroid are different things. The circumcenter is, well, the center of the circumcircle. This only makes sense for simplices, in our case the triangles/tetrahedra of the original mesh. A centroid is the center of mass of a domain, and voropy can compute the centroid of the Voronoi cells -- check out get_control_volume_centroids(). (The voronoi cells are called control volumes in the code, clearly the code comes from a finite volume context. Perhaps I should rename those.)

for a particular Voronoi cell 'q', the set of face circumcenter points associated with each face of the Voronoi cell and packaged as a set of values for each 'q' (I think I see these calculated in the MeshTet.show() method, but not as an official stored variable?).

In 3D, the faces of the Voronoi cells are polygons and might not have a circumcenter in the sense of a circle that goes through all polygon corners. Perhaps you mean something different?

for a particular Voronoi cell 'q', the volume of that Voronoi cell (which believe is the private variable _control_volumes?)

get_control_volumes()

for a particular Voronoi cell of index 'q', the surface area of each of the faces associated with that Voronoi cell.

For some background: It turns out that instead of computing the face area, the ratio face_area / edge_length can be computed much more easily. Instead of calling it face_area, the code calls is covolume -- perhaps there's a better name for it. Anyways, get_ce_ratio() gives you the ratio of face area and edge length. (Incidentally, this is the only value you need for most FVM computations. -- Perhaps for you too?)

Again, to avoid global summation, the ratios are not computed for each edge, but for each edge in a face in a cell. The return value of get_ce_ratio() will have the same shape as idx_hierarchy.

for a particular Voronoi cell of index 'q' the normal vectors to each face (it makes sense that these are edges of the original mesh, but how do we package them to correspond to the 'qth' Voronoi cell?) for a particular Voronoi cell of index 'q', the set of neighboring Voronoi cell indices for each 'q'.

Those two are the same, and they are hard to get. create_edges is perhaps what you're looking for, but really I don't remember any actual data being associated with an edge. Everything is an edge in a face in a cell.

nschloe commented 7 years ago

mesh._circumcenters # the center of the Voronoi cell? (why is this variable private??)

Nope, circumcenter of a polygon doesn't make sense. These are the circumcenters of the simplices of the dual mesh. (The circumcenters are mostly the points where the Voronoi cells have corners.)

mesh._control_volumes # the volume of the Voronoi cell? (also, why is this private???)

Indeed. I could perhaps make it public, but that means computing it every time, which is not always necessary. Check get_control_volumes().

mesh.cell_volumes # the volume of the original mesh tetrahedra?

Right.

mesh.circumcenter_face_distances # distance between Voronoi cell centre and each face of the Voronoi cell?

No. The distance between the simplex circumcenter and the faces of the simplex.

mesh.node_coords # the coordinates of the original mesh points

Yes.

mesh.idx_hierarchy # a fantastic way to structure the tet, face, edge, node information of the original mesh

Yes.

mesh.edge_coords # we already have coords of the original mesh nodes and the idx_hierarchy, so these are...redundant?

Indeed, though they should probably be called half_edge_coords (half-edge in the sense of an edge in a face in a cell).

mesh.edges # indices to edges of the original mesh tetrahedra?

It's a dictionary, but yeah. That's a variable that I wouldn't rely on though, since it requires singling out the edges. Again, almost always the half-edges are enough.

mesh.faces # but reason for 'edges' and 'nodes' indexing not understood...?

Again, a dictionary.

mesh.ei_dot_ei # what is this? mesh.ei_dot_ej # what is this?

It turns out that many values can be computed if the dot-products of the simplex edges are known. Luckily, they can be computed very fast, and the resuts are stores in the above variables. ei_dot_ei are the dot products of all simplex edges onto themselves, ei_dot_ej is the dot product of all simplex edges onto all other simplex edges in the same simplex.

mesh.cells # nodes index points of the original mesh, what is "opposing vertex"...?

opposing vertex is outdated.

mesh.node_face_cells #??? I don't get what this is packaging...or why...???

=idx_hierarchy... I think.

nschloe commented 7 years ago

Altogether I have to say that I haven't looked at tet meshes for too long, and they are in need of updating. If you want to play it safe, you start off with triangular meshes.

pietakio commented 7 years ago

Thanks for all your answers!

So final question for the day: what is/was your ultimate intended use scenario was for MeshTetra? Why did you create it?

Altogether I have to say that I haven't looked at tet meshes for too long, and they are in need of updating. If you want to play it safe, you start off with triangular meshes.

Ha, ha, that's hilarious :) No problems.

So I already have a 2D voronoi FVM solver system going, but am interested in moving to 3D, and the specific interest is in setting up FVM on a nicely bounded 3D Voronoi mesh with clearly defined Voronoi cells, centres, faces, neighbours, etc.

Tetgen makes Voronoi meshes with nice output, but they're prone to critical errors, like funny inserted points that mess up the Voronoi cell faces, and missing Voronoi cell faces, so I can't use them :(

Other 3D Voronoi mesh builders (Scipy's Qhull implementation, Tess) don't seem to clip to the bounding surface of the 3D object, so I'm out of luck there too as clipping the 3D Voronoi mesh and reoptimizing is one of the hardest parts.

Thus, the hope for voropy...

But that gets us back to what we were saying originally: even if voropy doesn't really output a 3D Voronoi mesh for FVM use, but calculates tetrahedron volumes, etc, and makes it easy to set up FVM on the tet mesh, then it's still of interest.

Anyway, thanks for your help with this.

pietakio commented 7 years ago

By the way, I read through your recent comments and get what's going on more now. It might work for my application, we'll see. Anyway, thanks!

nschloe commented 7 years ago

So I already have a 2D voronoi FVM solver system going, but am interested in moving to 3D, and the specific interest is in setting up FVM on a nicely bounded 3D Voronoi mesh with clearly defined Voronoi cells, centres, faces, neighbours, etc.

Sounds good, I'm trying to go somwhere similar with pyfvm. I'll have to go do that for 3D such that I can clean up mesh_tet along the way. Your interest will gives me some incentive to tackle it next. (It's been on my work list for a while.)

Other 3D Voronoi mesh builders (Scipy's Qhull implementation, Tess) don't seem to clip to the bounding surface of the 3D object, so I'm out of luck there too as clipping the 3D Voronoi mesh and reoptimizing is one of the hardest parts.

Always starting of from a simplex mesh (which is the way they come out of most mesh generators), I was surprised to see there were no tools that make a simplex mesh into a voronoi with, specifically respecting the boundaries. That's why I created voropy.

But that gets us back to what we were saying originally: even if voropy doesn't really output a 3D Voronoi mesh for FVM use, but calculates tetrahedron volumes, etc, and makes it easy to set up FVM on the tet mesh, then it's still of interest.

Yeah, I wouldn't hope for the the visualization. For too large meshes, one graphics card wouldn't be good enough anyways, and I'm personally more interested in the solution of FVM systems. (Those can easily be visualized by an specialized tool like ParaView once the solution has been written to a file.)

Let's say we keep the report open, I dig around pyfvm/voropy until it all looks better, and I'll get back to you then.

pietakio commented 7 years ago

Sounds good, I'm trying to go somwhere similar with pyfvm. I'll have to go do that for 3D such that I can clean up mesh_tet along the way. Your interest will gives me some incentive to tackle it next. (It's been on my work list for a while.)

Cool! Thanks!

Yeah, I wouldn't hope for the the visualization. For too large meshes, one graphics card wouldn't be good enough anyways, and I'm personally more interested in the solution of FVM systems. (Those can easily be visualized by an specialized tool like ParaView once the solution has been written to a file.)

I kind of disagree on this point. If one's graphics card is good enough to visualize the tet mesh, then one can visualize the Voronoi dual -- it's not that much more expensive. For example, I've got Tetgen's Voronoi, Tess' Voronoi, and Scipy's Qhull Voronoi's visualized no problem using matplotlib. It is true that trying to visualize the Voronoi graph as repeated matplotlib plot() calls will chug forever with no results (I tried the MeshTet.show() call with the broken plot sequence removed and it didn't return a plot for over several minutes so I killed it). BUT if you have the coordinates of each Voronoi cell's face vertices organized into a data structure, I find you can plot them very quickly as a Poly3DCollection (in addition to visualizing it in ParaView if that's desired instead).

For my application I absolutely need to be able to visualize the Voronoi cells as 3D geometric entities, as relevant data is placed on the centre of each face of the Voronoi cell, as well as at the centre of each Voronoi cell. What I'm doing isn't exactly a conventional FVM, and having the ability to work with as few as 16 to 250 cells in 3D would actually be enough.

So, just for crazy people like me, even if it's stupidly expensive to calculate, the essential feature on my voropy wish list is an optional method that allows the user to calculate the geometric features specific to each reconstructed Voronoi cell of the mesh, for plotting, but also for non-conventional discrete volume applications like mine. That bullet-point list I made before is basically what I'd like to be able to calculate in voropy, and to repeat the list with some revisions after our discussions:

But if the above is too much work, or not on your application list, or not possible to reconstruct given the nature of the problem space, no worries -- it's just what I'd cross my fingers and hope desperately for :)

Let's say we keep the report open, I dig around pyfvm/voropy until it all looks better, and I'll get back to you then.

Thank you! I really appreciate you looking into this, and look forwards to future updates.

nschloe commented 5 years ago

Since this issue was opened, meshplex (the artist formerly known as voropy) has developed quite a bit and some of the issues are resolved, including the 3D representation. Tentatively closing. Feel free to reopen or open a new issue.