[WIP] Polyregions #113
[WIP] Polyregions #113
Conversation
…eate proper regions.data from setrun.py
|
To reduce the added overhead we should look into implementing the following:
|
|
Also, this would be a good candidate for a clawpack enhancement proposal... |
|
the bounding box is already calculated and checked before the ray casting as is. Not sure about the ray-casting...looks correct to me and seems to work based on preliminary tests. However, the refinement doesn't seem to work. ray-casting deserves a look from others. |
|
I tried it on one example where it did seem to work, but I will check it On Sat, Nov 15, 2014 at 11:18 AM, David George [email protected]
|
|
this seems to work for my tests. python data class is implemented as well. regions can be specified as rectangles in the old way. or as polyregions : [minlevel,maxlevel, t1, tend,[(x1,y1),(x2,y2),...,(xN,yN)]] not sure about Travis CI stuff. |
|
The new version gives different refined patches in some of the unit tests -- not sure why yet. |
|
It's an issue of rounding when a region is aligned with the grid -- one more row of grid cells might be flagged in one version than the other, which can lead to changes in the gauge solution as the solution evolves. It would be nice if results agreed when regions are rectangles. Zoomed view of acoustics_2d_radial results at time t=0. The region is -0.2 < y < 0.2 in this case.
|
|
it should/could work if we are careful to always compute the locations using the exact same arithmetic formula. Marsha On Nov 21, 2014, at 9:05 PM, Randall J. LeVeque [email protected] wrote:
|
|
The original version is actually adding one more row of cells to the refined region than it should, while the new version is getting it right -- in this particular case. However -- I realized there's a major problem with the new approach of only checking whether the cell center lies in the polyregion via ray casting. In GeoClaw we often specify a region for refinement that is a small fraction of a cell on Level 1, and that probably doesn't include the cell center. The new method would not refine the cell that contains this region, whereas it should. We really need to test if the intersection of the grid cell with the polyregion is non-empty rather than just testing if the cell center lies in it. If we assume the polyregion is convex then I think this can be done fairly easily using the Separating Axis Theorem (checking if all vertices of the polygon lie outside the cell, which is easy, and if not testing if all vertices of the cell lie outside the polygon by checking planes defined by each side of the polygon via inner products). But if we allow non-convex polyregions then it seems harder -- requires splitting it up into a set of convex polygons. |
|
yes, good point, hadn't thought of that. rethinking this as I write:
seems like that would get all cases (?) but would be fairly expensive for non intersections. |
|
I think that might work, but if there's no segment-side intersection then I Computing segment intersections still takes a bit of work but wouldn't On Mon, Nov 24, 2014 at 10:59 AM, David George [email protected]
|
|
good points...(I think it managed to send an email before I thought through and edited my comment.) Also, we could first add a check to see if the polygon bounding box and cell are disjoint so that the expensive checks are only implemented for cells "near" the polygon (?). (Currently it just checks the bounding box and cell center.) |
|
This is starting to sound more and more like it is going to be really expensive to do right. Is it worthwhile trying to come up with another strategy before you spend too much time implementing this? |
|
Kyle and I had discussed some possible preprocessing that would make it more efficient, for example, preprocessing the polyregions into a list of rectangles that are easier to check. Marsha On Nov 25, 2014, at 11:16 AM, Kyle Mandli [email protected] wrote:
|
|
good idea, that might be the way to go, creating a set of rectangles...the refinement becomes a set of rectangles eventually anyway. On the flip side, with the polygon route it's in fairly limited cases where multiple checks would be implemented...the cell has to intersect the bounding box of the polygon...the check for which is the same as previously for rectangles...and if a cell is interior to the polygon then a single ray casting check gets it...which might be similar to checking many rectangles, depending on the number of vertices etc (?). Either way, if we take the polygon route we should probably flag polygons that happen to be rectangles so that the implementation is no less efficient than previously, if a user simply selects rectangles. there is of course a trade-off between checking regions and having a tightly defined area of allowed refinement...so maybe preprocessing the polygons allows more flexibility for optimizing the trade-off (?) |
|
If the rectangles are processing using the same “grid generation” mechanism that makes refined rectangle patches, there is an efficiency field that returns the number of “flagged points” over the total number of points in the rectangle. Something similar could be used here, in a preprocessing step. Then all the rectangles would be saved, and the same thing as currently done would be repeated for regridding. Marsha On Nov 25, 2014, at 1:07 PM, David George [email protected] wrote:
|
|
It would be best to keep the polyregions as (not necessarily convex) polygons if possible since I envision the user creating a polygon that follows some complicated coastline or valley, for example, where a high level of refinement is required after some time in order for the fgmax routines to work well (which monitors the maximum inundation depth over the full computation, for example). If the region is relatively easy to cover with a set of rectangles then the current code is sufficient, particularly with the Leaflet Widget tools @cekees is helping develop -- The attached screen shot shows a complicated polygon, but you can also easily select a bunch of rectangles. Pre-processing the regions should certainly be done as much as possible to speed up whatever intersection routines we use.
|
|
Are you allowing arbitrary polygonal regions? If so it might be useful to triangulate any polygon and test against bounding boxes of the sub-triangles. This should give you pretty tight bounding box tests. Also since you know that you are intersecting a shape against a mesh can you use the connectivity information of your mesh to avoid having to test in far away regions? Basically find the intersection of one polygon with a cell and then walk around to neighboring cells testing if they intersect the polygon? |
|
If you flagged every cell intersected by the polygonal region edges like @smoe1 suggests, is it then easy to flag every cell within the region (now you have a boundary of cells rather than segments)? Here's a link to a paper that is supposed to give a fast and robust algorithm for computing the cells intersected by the path. @mfarthin has used it for general unstructured meshes, and I believe it was efficient. Maybe it's even easier to trace paths for clawpack grids. http://www.sciencedirect.com/science/article/pii/S002199910700126X |
|
That said, it seems like @dlgeorge is close to adding the additional checks to flag the general polygon/cell interesections in a reasonably efficient way. |
|
I've been thinking about this some more and I don't see any way to make this sufficiently fast for practical geoclaw situations. Even if we require the polygon to be convex, it's easy to come up with cases where no corner of a cell is in the polygon and no vertex of the polygon is in the cell and yet there's a non-zero intersection (e.g. when the polygon is a rotation of the rectangular cell) and so we'd have to check for intersection of edges along with the various cases where one is entirely inside the other. To have any hope of making this efficient we'd first have to do it on a patchwise basis: First check if a grid patch lies entirely inside or entirely outside the polygon before bothering to check individual cells. Maybe using information from neighboring cells (or patches) can speed things up, but would be even more complex to implement. I'm in favor of leaving this as a future project and instead have focused on making the leafletwidget tool work well for specifying rectangles. I'll submit a PR in the apps directory for this. |
|
Maybe this would then be best accomplished as a preprocessing step then? Say have one of those leaflets allow for an arbitrary polygon but then construct the rectangles that would properly contain the polygon. |
|
Need to re-evaluate this later. |
don't merge:
polyregions are implemented....but doesn't seem to work yet. Note also that python data has not been altered yet. regions.data should be altered. See regions module for form of regions.data.