1
$\begingroup$

I have a simple square mesh, and a curve (discretised by another mesh) inside it. Here a picture worths thousand words. What I want to achieve is to find, for every cell $K$ of the circular (discretised) grid, I want to get the intersection points with the cells T of other grid. Here they are marked with a red dot.

enter image description here

What is the correct algorithmic to follow? Of course, if two consecutive support points for the circular grid are in different cells, then I can compute the intersection, but in general it may be that a cell $T$ is cut and no support points of the other mesh are in there, see again the picture above.

I started with FEnics, but I think in this case deal.II is the way to go. So, any explanation in the deal.II "lingo" is more than welcome.

Improve this question
$\endgroup$
11
  • $\begingroup$ Is it possible that one edge of your curve crosses multiple edges of the background mesh? $\endgroup$
    – Daniel Shapero
    Commented Jan 25, 2022 at 18:19
  • $\begingroup$ @DanielShapero yes, that is possible, and it is indeed what happens in the top right cell! $\endgroup$
    – FEGirl
    Commented Jan 25, 2022 at 18:35
  • $\begingroup$ Ah of course, I should have looked more carefully! In any case, if you're using deal.II I think a good place to start would be to find the cell containing an initial point of the curve and work by breadth-first search from there. The functions you'll want are probably compute_point_locations and get_neighbors. $\endgroup$
    – Daniel Shapero
    Commented Jan 25, 2022 at 19:48
  • $\begingroup$ When I had to deal with non-conforming meshes, I used Clipper to determine the cut areas needed for the domain integrals. Until now, I couldn't find software that can clip arbitrary polyhedra, so going to 3D could be difficult. $\endgroup$
    – Zoltan Csati
    Commented Jan 25, 2022 at 22:33
  • $\begingroup$ @DanielShapero I started with those functions indeed, but I don't know what you mean with "breadth-first search" in this case. $\endgroup$
    – FEGirl
    Commented Jan 25, 2022 at 22:51

2 Answers 2

Reset to default
2
$\begingroup$

The other answers and comments have good suggestions already. In practice, you will probably find that on coarse meshes most cells have no intersection with the curve and, assuming that your curve is contained in a relatively small part of the domain, you can easily exclude the case that a cell is intersected by the curve using bounding boxes. To this end, you split the curve into a moderate number of bounding boxes each of which covers one line segment either in its entirety, or a smaller piece of a line segment, and intersect them with the bounding box of a cell. If the intersection is empty, you know the curve doesn't go through the cell.

Comparing lots of bounding boxes of cells against lots of bounding boxes of line segments is expensive. You can make this much more efficient if you use an rtree of bounding boxes. deal.II has classes both for bounding boxes and for rtrees.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ Thanks @WolfgangBangerth. I understood the idea, and I like it a lot as indeed I'd be using deal.II capabilities to the full. What I am missing is how to compute those coordinates. From your answer, it seems that I only can say which cells are intersected. But that could've been achieved also by calling compute_point_locations() repeatedly. $\endgroup$
    – FEGirl
    Commented Jan 26, 2022 at 22:59
  • 1
    $\begingroup$ @bobinthebox I don't know about the details of how deal.II library works, but generally the process for using spatial data structures such as rtrees is you use the rtree to filter the number of elements you need to perform a complex intersection check on. So you would first find what cells are intersected, then manually (or using another library call if present) to check if your curve intersects any elements which are in the filtered cells. $\endgroup$
    – helloworld922
    Commented Jan 27, 2022 at 1:03
  • $\begingroup$ @bobinthebox compute_point_locations() is expensive on its own. You want to make sure that you only look at cells for which there is a reasonable chance that a point is inside. That's why GridTools::Cache exists. $\endgroup$
    – Wolfgang Bangerth
    Commented Jan 27, 2022 at 14:52
1
$\begingroup$

When I faced this problem in the past, in the context of domain decomposition across mismatched surface tessellations, I ended up basically using the Sutherland-Hodgman algorithm (polygon-to-polygon clipping) to clip every polygon against every other, combined with a broad-phase culling based on BVH/BSH over the interface elements. It can definitely be applied to this problem, though it is possibly overkill if you only want intersection points (and don't care about topology/polygonization of those points). One of the tricky bits with this is that the resulting clipped region can basically be an arbitrary polygon, so you might find that you need further algorithmic machinery (ie triangulation/meshing) in order to do any FE-like things with the result (like interpolation/resampling/bases).

Sutherland-Hodgman Algorithm

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.