Add support for tetrahedral connectivity and use scipy to generate it #6303
Conversation
nilsdeppe
force-pushed
the
tetrahedral_connectivity
branch
from
9c6c9eb
to
7a21ee1
Compare
| @@ -219,10 +223,16 @@ def _xmf_grid( | |||
| ) | |||
| else: | |||
| # Cover volume | |||
| topology_type = {3: "Hexahedron", 2: "Quadrilateral"}[topo_dim] | |||
| connectivity_name = "connectivity" | |||
| if use_tetrahedral_connectivity: | |||
There was a problem hiding this comment.
In Python it's easier to read if you move the code above into an else clause (scopes are more relaxed in Python so the variable definitions are available outside the else)
| default=False, | ||
| help=( | ||
| "Use a tetrahedral connectivity called tetrahedral_connectivity in " | ||
| "the HDF5 file" |
There was a problem hiding this comment.
Can you add here and/or in other docs what effect this will have, why you might want this, etc? Just looking for some context to tell the user why this exists and how to use it. Can also just point to the docs of GenerateTetrahedralConnectivity.py
| @@ -72,6 +72,46 @@ def test_generate_xdmf(self): | |||
| ), | |||
| ) | |||
|
|
|||
| def test_tetrahedral_generate_xdmf(self): | |||
There was a problem hiding this comment.
Is this test the same as above? Then you can move the two tests into a for loop
| can be used to fill gaps between finite-difference or Gauss elements. | ||
|
|
||
| If this algorithm is too slow, one possible improvement is to apply | ||
| qhull to each block of the domain and then connect the blocks to each |
There was a problem hiding this comment.
Note that this could also be done in parallel for each block or H5 file? I think parallelization over H5 files with multiprocessing.Pool would be the easiest way to speed this up, since it probably isn't IO bound.
There was a problem hiding this comment.
Yes, but then you wouldn't connect between blocks, which is one of the things I'm looking for this to do... However, we could still parallelize over blocks and then just connect the blocks afterward.
| stride: Optional. View only every stride'th time step. | ||
| coordinates: Optional. Name of coordinates dataset. Default: | ||
| "InertialCoordinates". | ||
| delete_tetrahedral_connectivity: Optional. Delete the dataset and return. |
There was a problem hiding this comment.
What's the purpose of this option? If the idea is to overwrite the dataset then for consistency with other tools I suggest you rename the option to force.
| "This may take a few minutes and may even fail depending on " | ||
| "the grid structure." | ||
| ) | ||
| delaunay = spatial.Delaunay(coords) |
There was a problem hiding this comment.
This won't connect points across H5 files, right? Does that lead to visualization artifacts? Delaunay supports iteratively adding points so you could triangulate across files, though I'm not sure where to store this global triangulation.
There was a problem hiding this comment.
That's right, so you need to first combine HDF5 files. I think it is doable to generate across HDF5 files, but considering the cost scales as N^2, I don't expect it to be practically useful to want to do that a lot... The bookkeeping is tricky because you need to get all the right indices in the end.
There was a problem hiding this comment.
Ah ok maybe mention in the docs that you should run combine-h5 first to avoid gaps in the rendering? And/or make this script take only a single H5 file rather than a list
|
|
||
|
|
||
| @click.command( | ||
| name="generate-tetrahedral-connectivity", |
There was a problem hiding this comment.
(optional) Maybe just name triangulate?
There was a problem hiding this comment.
I don't feel too strongly, I used this so it matches up with GenerateXdmf use-tetrahedral-connectivity but I agree it's very long
| type=float, | ||
| help=( | ||
| "The time at which to stop visualizing. The stop-time value is " | ||
| "not included." |
There was a problem hiding this comment.
Actually I think it is included, no? This might be a documentation bug in GenerateXdmf.py as well.
| if not subfile_name.endswith(".vol"): | ||
| subfile_name += ".vol" | ||
|
|
||
| for h5file, filename in h5files: |
There was a problem hiding this comment.
You can relatively easily show a progress bar here so the script doesn't appear "frozen" without any indication that it is doing anything. Here's how:
import rich.progress
progress = rich.progress.Progress(
rich.progress.TextColumn("[progress.description]{task.description}"),
rich.progress.BarColumn(),
rich.progress.MofNCompleteColumn(),
rich.progress.TimeRemainingColumn(),
disable=(len(h5files) == 1),
)
task_id = progress.add_task("Processing files", total=len(h5files))
with progress:
for h5file, filename in progress.track(h5files, task_id=task_id):
# ...
progress.update(task_id, advance=1) # don't remember if this is neededDocs are here: https://rich.readthedocs.io/en/stable/progress.html
| coordinates: str = "InertialCoordinates", | ||
| delete_tetrahedral_connectivity: bool = False, | ||
| ): | ||
| """Generate tetrahedral connectivity using scipy.spatial.Delaunay |
There was a problem hiding this comment.
Question: how does this handle shared grid points on LGL element boundaries? The tetrahedra connecting these points would be degenerate.
There was a problem hiding this comment.
I haven't read the paper that describes qhull, but as far as I can tell it works fine and doesn't produce any double connectivity.
nilsdeppe
force-pushed
the
tetrahedral_connectivity
branch
from
7a21ee1
to
459d56e
Compare
Proposed changes
This provides an easy way to fill gaps in for domains with Gauss points. There are several optimizations possible to make it run faster if needed.
Upgrade instructions
Code review checklist
make docto generate the documentation locally intoBUILD_DIR/docs/html.Then open
index.html.code review guide.
bugfixornew featureif appropriate.Further comments